拟一维喷管流动的数值解——亚声速-超声速等熵喷管流动的守恒型CFD解法（MacCormack方法）

一、Matlab代码片

%亚声速-超声速等熵喷管流动守恒形CFD解法 MacCormack方法
clear; %清理内存变量
clc; %清理工作窗中的所有显示内容
r=1.4; %比热比
L=3; %喷管长度
i=31; %网格数目
C=0.5; %科朗数
x=linspace(0,L,i); %网格点横坐标
N=1401; %时间步
%预分配内存
rho=zeros(1401,31); V=zeros(1401,31); T=zeros(1401,31);
U1=zeros(1401,31); U2=zeros(1401,31); U3=zeros(1401,31);
dx=L/(i-1); %空间步长
dt(1:N)=0; %时间步长
A=1+2.2*(x-1.5).*(x-1.5); %喷管面积
%初始条件
rho(1,:)=[ones(1,6),1.0-0.366.*(x(7:16)-0.5),0.634-0.3879.*(x(17:31)-1.5)];
T(1,:)=[ones(1,6),1.0-0.167.*(x(7:16)-0.5),0.833-0.3507.*(x(17:31)-1.5)]; 
U1(1,:)=rho(1,:).*A(1:31);
U2(1,:)=0.59;
V(1,:)=U2(1,:)./U1(1,:);
U3(1,:)=U1(1,:).*((T(1,:)./(r-1))+r/2.*V(1,:).^2);
%MacCormack方法迭代
for k=1:N-1
    %确定最小时间步长
    t=C*dx./(V(k,:)+sqrt(T(k,:)));
    dt(k)=min(t);
    %预估步通量
    F1(1:i)=U2(k,1:i);
    F2(1:i)=U2(k,1:i).^2./U1(k,1:i)+(r-1)/r.*...
    (U3(k,1:i)-r/2.*U2(k,1:i).^2./U1(k,1:i));
    F3(1:i)=r.*U2(k,1:i).*U3(k,1:i)./U1(k,1:i)...
    -r*(r-1)/2.*U2(k,1:i).^3./U1(k,1:i).^2;
    J2(1:i-1)=1/r.*rho(k,1:i-1).*T(k,1:i-1).*(A(2:i)-A(1:i-1))/0.1; %源项
    %预估偏导数
    U1t(1:i-1)=-(F1(2:i)-F1(1:i-1))/dx;
    U2t(1:i-1)=-(F2(2:i)-F2(1:i-1))/dx+J2(1:i-1);
    U3t(1:i-1)=-(F3(2:i)-F3(1:i-1))/dx;
    %预估值
    U1_(1:i-1)=U1(k,1:i-1)+U1t(1:i-1).*dt(k);
    U2_(1:i-1)=U2(k,1:i-1)+U2t(1:i-1).*dt(k);
    U3_(1:i-1)=U3(k,1:i-1)+U3t(1:i-1).*dt(k);
    rho_(1:i-1)=U1_(1:i-1)./A(1:i-1);
    T_(1:i-1)=(r-1).*( U3_(1:i-1)./U1_(1:i-1)-r/2.*U2_(1:i-1).^2./U1_(1:i-1).^2);
    F1_(1:i-1)=U2_(1:i-1);
    F2_(1:i-1)=U2_(1:i-1).^2./U1_(1:i-1)+(r-1)/r.*...
    (U3_(1:i-1)-r/2.*U2_(1:i-1).^2./U1_(1:i-1));
    F3_(1:i-1)=r.*U2_(1:i-1).*U3_(1:i-1)./U1_(1:i-1)...
    -r*(r-1)/2.*U2_(1:i-1).^3./U1_(1:i-1).^2;
    J2_(2:i-1)=1/r.*rho_(2:i-1).*T_(2:i-1).*(A(2:i-1)-A(1:i-2))/0.1; %注意
    %校正步 校正偏导数
    U1t_(2:i-1)=-(F1_(2:i-1)-F1_(1:i-2))/dx;
    U2t_(2:i-1)=-(F2_(2:i-1)-F2_(1:i-2))/dx+J2_(2:i-1);
    U3t_(2:i-1)=-(F3_(2:i-1)-F3_(1:i-2))/dx;
    %时间导数平均值
    U1tav(2:i-1)=0.5.*(U1t(2:i-1)+U1t_(2:i-1));
    U2tav(2:i-1)=0.5.*(U2t(2:i-1)+U2t_(2:i-1));
    U3tav(2:i-1)=0.5.*(U3t(2:i-1)+U3t_(2:i-1));
    %最终校正值
    U1(k+1,2:i-1)=U1(k,2:i-1)+U1tav(2:i-1).*dt(k);
    U2(k+1,2:i-1)=U2(k,2:i-1)+U2tav(2:i-1).*dt(k);
    U3(k+1,2:i-1)=U3(k,2:i-1)+U3tav(2:i-1).*dt(k);
    rho(k+1,2:i-1)=U1(k+1,2:i-1)./A(2:i-1);
    V(k+1,2:i-1)=U2(k+1,2:i-1)./U1(k+1,2:i-1);
    T(k+1,2:i-1)=(r-1).*(U3(k+1,2:i-1)./U1(k+1,2:i-1)-r/2.*V(k+1,2:i-1).^2);
    %入流边界值
    rho(k+1,1)=1;
    T(k+1,1)=1;
    U1(k+1,1)=A(1);
    U2(k+1,1)=2.*U2(k+1,2)-U2(k+1,3);
    V(k+1,1)=U2(k+1,1)./U1(k+1,1);
    U3(k+1,1)=U1(k+1,1).*((T(k+1,1)./(r-1))+r/2.*V(k+1,1).^2);
    %出流边界值
    U1(k+1,i)=2.*U1(k+1,i-1)-U1(k+1,i-2);
    U2(k+1,i)=2.*U2(k+1,i-1)-U2(k+1,i-2);
    U3(k+1,i)=2.*U3(k+1,i-1)-U3(k+1,i-2);
    %由于不需要返回原变量的出流边界值，故这里省去了原变量出流边界值的计算
end
%绘图 无量纲质量流量在六个不同时刻的瞬时分布
figure;
plot(x,U2(1,:),'k--');
hold on;
plot(x,U2(101,:),'b-');
hold on;
plot(x,U2(201,:),'y-');
hold on;
plot(x,U2(701,:),'m-');
axis([0,3,0.35,0.7]);
title('不同时刻质量流量的变化');
ylabel('无量纲质量流量'),xlabel('无量纲轴向距离');
text(2.6,0.58,'700\Deltat');
text(1.8,0.56,'200\Deltat'); 
text(1.5,0.66,'100\Deltat');
text(1.5,0.6,'0\Deltat');
    

二、计算结果展示

无量纲质量流量在四个不同时刻的瞬时分布