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

一、Matlab代码片

%全亚声速等熵喷管流动 非守恒型麦考马克方法数值求解
clear; %清理内存变量
clc; %清理工作窗中的所有显示内容
r=1.4; %比热比
L=3; %喷管长度
i=31; %网格数目
C=0.5; %科朗数
x=linspace(0,L,i); %网格点横坐标
N=5001; %时间步
dx=L/(i-1); %空间步长
dt(1:N)=0; %时间步长
Pe=0.93; %喷管两端压力比
A(1:31)=0; %预分配内存
for j=1:31
if(x(j)<1.5)
    A(j)=1+2.2*(x(j)-1.5).*(x(j)-1.5); %喷管面积
else
    A(j)=1+0.2223*(x(j)-1.5).*(x(j)-1.5);%喷管面积
end
end
rho(N,i)=0; %流场密度赋值
T(N,i)=0; %流场温度赋值
V(N,i)=0; %流场速度赋值
%预分配内存
rhotav2(1:1400,1:30)=0;  Vtav2(1:1400,1:30)=0;
%初始条件
rho(1,:)=1-0.023*x; %流场密度初值
T(1,:)=1-0.009333*x; %流场温度初值
V(1,:)=0.05+0.11*x; %流场速度初值
%按时间步长推进
for k=1:N-1
    %计算预估步偏导数
    rhot(1:i-1)=-V(k,1:i-1).*(rho(k,2:i)-rho(k,1:i-1))/dx-rho(k,1:i-1)...
    .*(V(k,2:i)-V(k,1:i-1))/dx-rho(k,1:i-1).*V(k,1:i-1)...
    .*(log(A(2:i))-log(A(1:i-1)))/dx;
    Vt(1:i-1)=-V(k,1:i-1).*(V(k,2:i)-V(k,1:i-1))/dx-1/r.*((T(k,2:i)...
    -T(k,1:i-1))/dx+T(k,1:i-1)./rho(k,1:i-1).*(rho(k,2:i)-rho(k,1:i-1))/dx);
    Tt(1:i-1)=-V(k,1:i-1).*(T(k,2:i)-T(k,1:i-1))/dx-(r-1).*T(k,1:i-1)...
    .*((V(k,2:i)-V(k,1:i-1))/dx+V(k,1:i-1).*(log(A(2:i))-log(A(1:i-1)))/dx);
    %确定最小时间步长
    t=C*dx./(V(k,2:i-1)+sqrt(T(k,2:i-1)));
    dt(k)=min(t);
    %计算预估值
    rho_(1:i-1)=rho(k,1:i-1)+rhot(1:i-1).*dt(k);
    V_(1:i-1)=V(k,1:i-1)+Vt(1:i-1).*dt(k);
    T_(1:i-1)=T(k,1:i-1)+Tt(1:i-1).*dt(k);
    %校正偏导数
    rhot_(2:i-1)=-V_(2:i-1).*(rho_(2:i-1)-rho_(1:i-2))/dx-rho_(2:i-1)...
    .*(V_(2:i-1)-V_(1:i-2))/dx-rho_(2:i-1).*V_(2:i-1).*(log(A(2:i-1))-log(A(1:i-2)))/dx;
    Vt_(2:i-1)=-V_(2:i-1).*(V_(2:i-1)-V_(1:i-2))/dx-1/r.*((T_(2:i-1)...
    -T_(1:i-2))/dx+T_(2:i-1)./rho_(2:i-1).*(rho_(2:i-1)-rho_(1:i-2))/dx);
    Tt_(2:i-1)=-V_(2:i-1).*(T_(2:i-1)-T_(1:i-2))/dx-(r-1).*T_(2:i-1)...
    .*((V_(2:i-1)-V_(1:i-2))/dx+V_(2:i-1).*(log(A(2:i-1))-log(A(1:i-2)))/dx);
    %时间导数平均值
    rhotav(2:i-1)=0.5*(rhot(2:i-1)+rhot_(2:i-1)); rhotav2(k,2:i-1)=abs(rhotav(2:i-1));
    Vtav(2:i-1)=0.5*(Vt(2:i-1)+Vt_(2:i-1));Vtav2(k,2:i-1)=abs(Vtav(2:i-1));
    Ttav(2:i-1)=0.5*(Tt(2:i-1)+Tt_(2:i-1));
    %流场变量校正值
    rho(k+1,2:i-1)=rho(k,2:i-1)+rhotav(2:i-1)*dt(k);
    V(k+1,2:i-1)=V(k,2:i-1)+Vtav(2:i-1)*dt(k);
    T(k+1,2:i-1)=T(k,2:i-1)+Ttav(2:i-1)*dt(k);
    %入流边界值
    V(k+1,1)=2*V(k+1,2)-V(k+1,3);
    rho(k+1,1)=1;
    T(k+1,1)=1;
    %出流边界值
    rho(k+1,i)=2*rho(k+1,i-1)-rho(k+1,i-2);
    V(k+1,i)=2*V(k+1,i-1)-V(k+1,i-2);
    T(k+1,i)= Pe./rho(k+1,i);
    %马赫数
    Ma=V(1:k+1,1:i)./(sqrt(T(1:k+1,1:i)));
    %流场压强
    P=rho(1:k+1,1:i).*T(1:k+1,1:i);
end
%绘图 喷管喉道处密度、温度、压力和马赫数的变化
figure;
subplot(2,2,1),plot(1:5001,rho(1:5001,16),'b-'); 
ylabel('\rho/\rho_0'),xlabel('Timestep'); title('无量纲密度变化');
subplot(2,2,2),plot(1:5001,T(1:5001,16),'r-');   
ylabel('T/T_0'),xlabel('Timestep');title('无量纲温度变化');
subplot(2,2,3),plot(1:5001,P(1:5001,16),'k-');   
ylabel('P/P_0'),xlabel('Timestep');title('无量纲总压变化');
subplot(2,2,4),plot(1:5001,Ma(1:5001,16),'m-');  
ylabel('Ma'),xlabel('Timestep');title('马赫数变化');
%绘图 喷管喉道处无量纲密度和速度时间导数的变化
figure;
plot(1:5000,rhotav2(1:5000,16),'k-');
ylabel('残差'),xlabel('Timestep'); 
title('喷管喉道处无量纲密度和速度时间导数的变化');
hold on;
plot(1:5000,Vtav2(1:5000,16),'k--');
%绘图 无量纲质量流量在三个不同时刻的瞬时分布
figure;
plot(x,rho(1,1:i).*V(1,1:i).*A(1:i),'k--');
hold on;
plot(x,rho(501,1:i).*V(501,1:i).*A(1:i),'r-');
hold on;
plot(x,rho(5001,1:i).*V(5001,1:i).*A(1:i),'b-');
title('不同刻质量流量的变化');
ylabel('无量纲质量流量'),xlabel('无量纲轴向距离');
text(3.2,0.47,'5000\Deltat');
text(3.2,0.40,'500\Deltat');
text(2.0,0.25,'0\Deltat'); 
%绘图 不同时刻压力分布的变化
figure;
plot(x,rho(1,1:i).*T(1,1:i),'k--');
hold on;
plot(x,rho(1001,1:i).*T(1001,1:i),'m-');
hold on;
plot(x,rho(501,1:i).*T(501,1:i),'r-');
hold on;
plot(x,rho(5001,1:i).*T(5001,1:i),'b-');
hold on;
title('不同时刻压力分布的变化');
ylabel('P/P_0'),xlabel('x/L');
text(1.5,0.96,'0\Deltat');
text(2.0,0.91,'500\Deltat');
text(2.2,0.89,'1000\Deltat');
text(2.3,0.867,'5000\Deltat'); 

二、计算结果展示

1.喷管喉道处密度、温度、压力和马赫数的变化

2.喷管喉道处无量纲密度和速度时间导数的变化

3.无量纲质量流量在三个不同时刻的瞬时分布

4.无量纲压力在四个不同时刻的瞬时分布