%function charpt7_3 %安德森《计算流体力学及其应用》第七章拟一维喷管流动的数值解 7.3节 亚声速—超声速等熵喷管流动的CFD解法 。
clear;
clc;
L=1;%喷管长L=1m。
N=30;%x等分段数;喷管形状A(x)=1+2.2*(x-1.5)^2,0<=x<=3,即deltax=0.1。
x=linspace(0,3,N+1)./L;%x无量纲化。
%第一步:确定初值条件,t=0时刻。
Rho=1-0.314.*x;%密度初值。
T=1-0.2314.*x;%温度初值。
V=(0.1+1.09.*x).*sqrt(T);%速度初值。
A0=1;%喉道处面积。
A=(1+2.2.*(x-1.5).^2)./A0;%喷管沿x方向的面积无量纲化。
%输出喷管形状和初值条件表,即书中表7—1。
Biao7_1=zeros(N+1,5);%表中数据箱子。
Biao7_1(:,1)=x';
Biao7_1(:,2)=A';
Biao7_1(:,3)=Rho';
Biao7_1(:,4)=V';
Biao7_1(:,5)=T';
save '表7-1 喷管形状和初始条件.txt' Biao7_1 -ascii;%将表7-1 喷管形状和初始条件数据以十进制保存到txt文档,并命名为“表7-1 喷管形状和处世条件.txt”。
deltax=3/N;%网格尺寸。
Bushu=1400;%总的时间步数。
%喉道处密度随时间步变化值箱子。
Rhox=1:1:Bushu;%行向量
Rhoy=zeros(1,Bushu);
P0=zeros(1,Bushu);
T0=zeros(1,Bushu);
Ma0=zeros(1,Bushu);
gamma=1.4;%比热比,标准状态下的空气,gamma=1.4。
C=0.5;%柯朗数。
%第二步,预估步t=0时导数值箱子。
Rho1=zeros(1,N+1);
T1=zeros(1,N+1);
V1=zeros(1,N+1);
%第三步,预估值t=deltat时预估值箱子。
Rho2=zeros(1,N+1);
T2=zeros(1,N+1);
V2=zeros(1,N+1);
%第四步,校正步到导数值箱子。
Rho3=zeros(1,N+1);
T3=zeros(1,N+1);
V3=zeros(1,N+1);
%第五步,平均时间导数值箱子。
Rho34=zeros(1,N+1);
T4=zeros(1,N+1);
V4=zeros(1,N+1);
i=1;%第i个时间步。
while i<=Bushu%步数推进条件
%开始第二步,t=0时导数值。
for j=1:1:N
Rho1(j)=-V(j)*(Rho(j+1)-Rho(j))/deltax-Rho(j)*(V(j+1)-V(j))/deltax-Rho(j)*V(j)*(log(A(j+1))-log(A(j)))/deltax;
V1(j)=-V(j)*(V(j+1)-V(j))/deltax-1/gamma*((T(j+1)-T(j))/deltax+T(j)/Rho(j)*(Rho(j+1)-Rho(j))/deltax);
T1(j)=-V(j)*(T(j+1)-T(j))/deltax-(gamma-1)*T(j)*((V(j+1)-V(j))/deltax+V(j)*(log(A(j+1))-log(A(j)))/deltax);
end
%开始第三步,t=deltat时预估值。
deltatt=C*(deltax./(V+sqrt(T)));
deltat=min(deltatt(2:N));
Rho2=Rho+Rho1.*deltat;
V2=V+V1.*deltat;
T2=T+T1.*deltat;
%第四步,校正步。
for j=2:1:N+1
Rho3(j)=-V2(j)*(Rho2(j)-Rho2(j-1))/deltax-Rho2(j)*(V2(j)-V2(j-1))/deltax-Rho2(j)*V2(j)*(log(A(j))-log(A(j-1)))/deltax;
V3(j)=-V2(j)*(V2(j)-V2(j-1))/deltax-1/gamma*((T2(j)-T2(j-1))/deltax+T2(j)/Rho2(j)*(Rho2(j)-Rho2(j-1))/deltax);
T3(j)=-V2(j)*(T2(j)-T2(j-1))/deltax-(gamma-1)*T2(j)*((V2(j)-V2(j-1))/deltax+V2(j)*(log(A(j))-log(A(j-1)))/deltax);
end
%第五步,平均导数值。
Rho4=0.5*(Rho1+Rho3);
V4=0.5*(V1+V3);
T4=0.5*(T1+T3);
%第六步。
Rho=Rho+Rho4*deltat;
V=V+V4*deltat;
T=T+T4*deltat;
P=Rho.*T;
Ma=V./sqrt(T);
Rho(1)=1;
V(1)=2*V(2)-V(3);
T(1)=1;
Rho(N+1)=2*Rho(N)-Rho(N-1);
V(N+1)=2*V(N)-V(N-1);
T(N+1)=2*T(N)-T(N-1);
%喷管喉道处密度的随时间步的变化,即喷管中间处。
Rhoy(i)=Rho(N/2+1);
T0(i)=T(N/2+1);
P0(i)=T(N/2+1);
Ma0(i)=Ma(N/2+1);
i=i+1;
end
%输出1400时间步后的流场变量表,即书中表7—2。
Biao7_2=zeros(N+1,5);%表中数据箱子。
Biao7_2(:,1)=x';
Biao7_2(:,2)=A';
Biao7_2(:,3)=Rho';
Biao7_2(:,4)=V';
Biao7_2(:,5)=T';
Biao7_2(:,6)=P';
Biao7_2(:,7)=Ma';
save '表7-2 第一个时间步后的流场变量.txt' Biao7_2 -ascii;%将表7-2 第一个时间步后的流场变量数据以十进制保存到txt文档,并命名为“表7-2 第一个时间步后的流场变量.txt”。
figure(1);
subplot 221
plot(Rhox,Rhoy,'r-');%绘制喉道处流场密度随时间步的变化图。
title('喉道处流场密度随时间步的变化');%图形标题。
xlabel('时间步数');%x坐标名称。
ylabel('密度(无量纲)');%y坐标名称。
legend('变化曲线');%添加图例。
subplot 222
plot(Rhox,T0,'B-');%绘制喉道处流场密度随时间步的变化图。
title('喉道处流场温度随时间步的变化');%图形标题。
xlabel('时间步数');%x坐标名称。
ylabel('温度(无量纲)');%y坐标名称。
legend('变化曲线');%添加图例。
subplot 223
plot(Rhox,P0,'y-');%绘制喉道处流场密度随时间步的变化图。
title('喉道处流场压强随时间步的变化');%图形标题。
xlabel('时间步数');%x坐标名称。
ylabel('压强(无量纲)');%y坐标名称。
legend('变化曲线');%添加图例。
subplot 224
plot(Rhox,Ma0,'g-');%绘制喉道处流场密度随时间步的变化图。
title('喉道处流场Ma随时间步的变化');%图形标题。
xlabel('时间步数');%x坐标名称。
ylabel('Ma(无量纲)');%y坐标名称。
legend('变化曲线');%添加图例。
%axis([0,Bushu,0.4,2]);%设置坐标范围axis([xmin,xmax,ymin,ymax,zmin,zmax])。%将y轴向下:set(gca,'ydir','reverse')。
grid on;