Matlab 编程 《计算流体力学基础及应用（约翰D安德森）》 全亚声速等熵喷管流动CFD解法 拉瓦尔喷管 非守恒形式方程解法

问题之 全亚声速等熵喷管流动的CFD解法

问题的提法见《计算流体力学基础及应用（约翰D安德森）》中的7.4节，下面直接给出写的程序。

初始化参数

%全亚声速的喷管流动 非守恒形式方程解法
%MacCormack方法 预估-校正 具有二阶精度
%初始参数（无量纲）
clc;clear;clf
x1=0:0.1:1.5;
x2=1.6:0.1:3;
x=[x1,x2];
N=length(x);
A1=1+2.2*(x1-1.5).^2;%截面函数
A2=1+0.2223*(x2-1.5).^2;
A=[A1,A2];
rho=1-0.023.*x;%密度
T=1-0.009333.*x;%温度
V=0.05+0.11.*x;%速度
p_e=0.93;%出口压力
gm=1.4;%比热比
C=0.5;%科朗数
rho_t=zeros(1,N);%一阶密度偏导数
rho_mid=1-0.023.*x;%预估步密度
rho_tt=zeros(1,N);%校正步密度
rho_av=zeros(1,N);%校正步密度
V_t=zeros(1,N);%一阶速度偏导数
V_mid=0.05+0.11.*x;%预估步速度
V_tt=zeros(1,N);%校正步速度
V_av=zeros(1,N);%校正步速度
T_t=zeros(1,N);%一阶温度偏导数
T_mid=1-0.009333.*x;%预估步温度
T_tt=zeros(1,N);%校正步温度
T_av=zeros(1,N);%校正步温度
t=zeros(1,N);%局部步长
k=0;%迭代次数
k0=5000;%指定最大迭代次数

迭代过程

%迭代过程
while 1
    k=k+1;
% 计算预估步 向前差分
for i=1:N
    t(i)=C.*0.1./(V(i)+sqrt(T(i)));%计算满足稳定性条件的各个局部步长
end
delta_t=min(t);%取最小步长
for i=2:N-1
    rho_t(i)=-rho(i).*(V(i+1)-V(i))/0.1-V(i).*(rho(i+1)-rho(i))/0.1- ...
        rho(i).*V(i).*(log(A(i+1))-log(A(i)))/0.1;%计算密度的一阶偏导
    V_t(i)=-V(i).*(V(i+1)-V(i))/0.1-1/gm.*((T(i+1)-T(i))/0.1+T(i)./rho(i) ...
        .*(rho(i+1)-rho(i))/0.1);%计算速度的一阶偏导
    T_t(i)=-V(i).*(T(i+1)-T(i))/0.1-(gm-1).*T(i).*((V(i+1)-V(i))/0.1+V(i) ...
        .*(log(A(i+1))-log(A(i)))/0.1);%计算温度的一阶偏导
    rho_mid(i)=rho(i)+rho_t(i).*delta_t;%密度预估
    V_mid(i)=V(i)+V_t(i).*delta_t;%速度预估
    T_mid(i)=T(i)+T_t(i).*delta_t;%温度预估
end
V_mid(1)=2.*V_mid(2)-V_mid(3);%入口点1处通过点2 3的插值获得
%计算校正步 向后差分
for i=2:N-1
    rho_tt(i)=-rho_mid(i).*(V_mid(i)-V_mid(i-1))/0.1-V_mid(i).*(rho_mid(i) ...
        -rho_mid(i-1))/0.1-rho_mid(i).*V_mid(i).*(log(A(i))-log(A(i-1)))/0.1;%计算密度的一阶偏导
    V_tt(i)=-V_mid(i).*(V_mid(i)-V_mid(i-1))/0.1-1/gm.*((T_mid(i)-T_mid(i-1)) ...
        /0.1+T_mid(i)./rho_mid(i).*(rho_mid(i)-rho_mid(i-1))/0.1);%计算速度的一阶偏导
    T_tt(i)=-V_mid(i).*(T_mid(i)-T_mid(i-1))/0.1-(gm-1).*T_mid(i).* ...
        ((V_mid(i)-V_mid(i-1))/0.1+V_mid(i).*(log(A(i))-log(A(i-1)))/0.1);%计算温度的一阶偏导
    rho_av(i)=1/2.*(rho_t(i)+rho_tt(i));%密度校正
    V_av(i)=1/2.*(V_t(i)+V_tt(i));%速度校正
    T_av(i)=1/2.*(T_t(i)+T_tt(i));%温度校正
end
rho=rho+rho_av.*delta_t;%下一时间步密度
rho(N)=2.*rho(N-1)-rho(N-2);%出口处点N密度通过插值获得
V=V+V_av.*delta_t;%下一时间步速度
V(1)=2.*V(2)-V(3);%入口点1处通过点2 3的插值获得
V(N)=2.*V(N-1)-V(N-2);%出口处点N速度通过插值获得
T=T+T_av.*delta_t;%下一时间步温度
T(N)=p_e./rho(N);%出口处点N温度通过条件式获得
Ma=V./sqrt(T);
%喉道参数获取
rhot(k)=rho(16);
Vt(k)=V(16);
Tt(k)=T(16);
Mat(k)=Ma(16);
pt(k)=rho(16).*T(16);
re_rho(k)=abs(rho_av(16));
%绘制不同时间步质量流量变化
if k==1
    xtick=1:N;
    rhoAV=rho(xtick).*A(xtick).*V(xtick);
    plot(xtick/10,rhoAV,'-b')
    hold on
    pt1=rho.*T;
    rho1=rho;
    T1=T;
    Ma1=Ma;
elseif k==500
    xtick=1:N;
    rhoAV=rho(xtick).*A(xtick).*V(xtick);
    plot(xtick/10,rhoAV,'-r')
    hold on
    pt500=rho.*T;
    rho500=rho;
    T500=T;
    Ma500=Ma;    
elseif k==5000
    xtick=1:N;
    rhoAV=rho(xtick).*A(xtick).*V(xtick);
    plot(xtick/10,rhoAV,'-g')
    xlabel('x/L'),ylabel('\rhoAV/\rho_0A_1a_0')
    title('不同时刻质量流量的变化')
    text(2,rhoAV(N)+0.02,'5000\Deltat')
    text(2,rhoAV(N)-0.05,'500\Deltat')
    text(2,0.31,'1\Deltat')
    axis([0,3.5,0,0.6])
    hold on
    pt5000=rho.*T;
    rho5000=rho;
    T5000=T;
    Ma5000=Ma;
else
    
end
%迭代次数
if k==k0
    break;
else
    
end 
end

绘图

%绘图
figure
x0=1:1:k;%时间步
%绘制密度-时间步图像
rho_change=rhot(x0);
subplot(2,2,1),plot(x0,rho_change)
ylabel('\rho/\rho_0'),xlabel('时间步数')
title('无量纲密度变化')
%绘制温度-时间步图像
T_change=Tt(x0);
subplot(2,2,2),plot(x0,T_change)
ylabel('T/T_0'),xlabel('时间步数')
title('无量纲温度变化')
%绘制压力-时间步图像
p_change=pt(x0);
subplot(2,2,3),plot(x0,p_change)
ylabel('p/p_0'),xlabel('时间步数')
title('无量纲总压变化')
%绘制马赫数-时间步图像
Ma_change=Mat(x0);
subplot(2,2,4),plot(x0,Ma_change)
ylabel('Ma'),xlabel('时间步数')
title('马赫数变化')
fprintf('迭代最终结果:\n密度rho=%6.4e\n温度T=%6.4e\n速度V=%6.4e\n总压p=%6.4e\n马赫数Ma=%6.4e\n' ...
    ,rhot(k),Tt(k),pt(k),Vt(k),Mat(k))
%绘制残差-时间步图像
figure 
re_rho_change=re_rho(x0);
semilogy(x0,re_rho_change)
ylabel('残差'),xlabel('时间步数')
title('残差随时间步变化')
x00=1:1:N;
figure
p0=pt1;
p1=pt500;
p2=pt5000;
subplot(2,2,1)
plot(x00/10,p0)
hold on
plot(x00/10,p1)
plot(x00/10,p2)
hold off
text(2.4,0.86,'5000\Deltat')
text(1.6,0.89,'500\Deltat')
text(1.8,0.95,'1\Deltat')
ylabel('p/p_0'),xlabel('x/L')
title('喷管内压力分布')
Ma0=Ma1;
Ma1=Ma500;
Ma2=Ma5000;
subplot(2,2,2)
plot(x00/10,Ma0)
hold on
plot(x00/10,Ma1)
plot(x00/10,Ma2)
hold off
text(2.2,0.52,'5000\Deltat')
text(1.8,0.4,'500\Deltat')
text(1.2,0.15,'1\Deltat')
ylabel('Ma'),xlabel('x/L')
title('喷管内马赫数分布')
T0=T1;
T1=T500;
T2=T5000;
subplot(2,2,3)
plot(x00/10,T0)
hold on
plot(x00/10,T1)
plot(x00/10,T2)
hold off
text(2.3,0.96,'5000\Deltat')
text(1.5,0.97,'500\Deltat')
text(1.3,1,'1\Deltat')
ylabel('T/T_0'),xlabel('x/L')
title('喷管内温度分布')
rho0=rho1;
rho1=rho500;
rho2=rho5000;
subplot(2,2,4)
plot(x00/10,rho0)
hold on
plot(x00/10,rho1)
plot(x00/10,rho2)
hold off
text(2.2,0.91,'5000\Deltat')
text(1.7,0.93,'500\Deltat')
text(1.5,0.97,'1\Deltat')
ylabel('\rho/\rho_0'),xlabel('x/L')
title('喷管内密度分布')

结果

迭代最终结果:
密度rho=8.6164e-01
温度T=9.4219e-01
速度V=8.1182e-01
总压p=5.4207e-01
马赫数Ma=5.5845e-01

Figure 1

Figure 2

Figure 3

Figure 4