拟一维喷管流动的数值解法(MATLAB)代码

拟一维喷管流动的数值解法（MATLAB）代码
数值计算代码
%拟一维喷管流动的数值解
%亚声速-超声速，非守恒形式
function main()
clear;
clc;
r=1.4; %绝热指数
N=1001; %时间步长
i=31; %网格数目
L=3; %喷管长度
C=0.5; %柯朗数
dx=L/(i-1); %空间步长
dt(N)=0; %时间步长
x=linspace(0,L,i); %网格点横坐标
A=1+2.2*(x-1.5).^2; %喷管面积
%赋值
M(N,i)=0;
T(N,i)=0;
V(N,i)=0;
%初始条件
M(1,:)=1-0.3146*x;
T(1,:)=1-0.2314*x;
V(1,:)=(0.1+1.09*x).*(1-0.2314*x).^0.5;
%按时间步长推进
for k=1:N-1
%预估偏导数
M_t(1:i-1)=-V(k,1:i-1).*(M(k,2:i)-M(k,1:i-1))/dx-M(k,1:i-1).*(V(k,2:i)-V(k,1:i-1))/dx-M(k,1:i-1).*V(k,1:i-1).*log(A(2:i)./A(1:i-1))/dx;
V_t(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)./M(k,1:i-1).*(M(k,2:i)-M(k,1:i-1))/dx);
T_t(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)./A(1:i-1))/dx);
%求取内部网格点处最小时间步长
t=C*dx./(V(k,2:i-1)+sqrt(T(k,2:i-1)));
dt(k)=min(t);
%预估值
M1(1:i-1)=M(k,1:i-1)+M_t(1:i-1)*dt(k);
V1(1:i-1)=V(k,1:i-1)+V_t(1:i-1)*dt(k);
T1(1:i-1)=T(k,1:i-1)+T_t(1:i-1)*dt(k);
%校正偏导数
M_t_1(2:i-1)=-V1(2:i-1).*(M1(2:i-1)-M1(1:i-2))./dx-M1(2:i-1).*(V1(2:i-1)-V1(1:i-2))./dx-M1(2:i-1).*V1(2:i-1).*log(A(2:i-1)./A(1:i-2))./dx;
V_t_1(2:i-1)=-V1(2:i-1).*(V1(2:i-1)-V1(1:i-2))./dx-1/r.*((T1(2:i-1)-T1(1:i-2))./dx+T1(2:i-1)./M1(2:i-1).*(M1(2:i-1)-M1(1:i-2))./dx);
T_t_1(2:i-1)=-V1(2:i-1).*(T1(2:i-1)-T1(1:i-2))./dx-(r-1).*T1(2:i-1).*((V1(2:i-1)-V1(1:i-2))./dx+V1(2:i-1).*log(A(2:i-1)./A(1:i-2))./dx);
%偏导数平均值
M_t_av(2:i-1)=0.5*(M_t(2:i-1)+M_t_1(2:i-1));
V_t_av(2:i-1)=0.5*(V_t(2:i-1)+V_t_1(2:i-1));
T_t_av(2:i-1)=0.5*(T_t(2:i-1)+T_t_1(2:i-1));
%内部网格点修正值
M(k+1,2:i-1)=M(k,2:i-1)+M_t_av(2:i-1)*dt(k);
V(k+1,2:i-1)=V(k,2:i-1)+V_t_av(2:i-1)*dt(k);
T(k+1,2:i-1)=T(k,2:i-1)+T_t_av(2:i-1)*dt(k);
%出口边界值
M(k+1,i)=2*M(k+1,i-1)-M(k+1,i-2);
V(k+1,i)=2*V(k+1,i-1)-V(k+1,i-2);
T(k+1,i)=2*T(k+1,i-1)-T(k+1,i-2);
%入口边界值
M(k+1,1)=1;
V(k+1,1)=2*V(k+1,2)-V(k+1,3);
T(k+1,1)=1;
end
end
图形处理代码
close all;
a=[1 51 101 151 201 701];
plot(x,M(a(1),:).*A(:)'.*V(a(1),:),'r--')
hold on
plot(x,M(a(2),:).*A(:)'.*V(a(2),:),'m.-')
plot(x,M(a(3),:).*A(:)'.*V(a(3),:),'g*-')
plot(x,M(a(4),:).*A(:)'.*V(a(4),:),'co-')
plot(x,M(a(5),:).*A(:)'.*V(a(5),:),'bh-')
plot(x,M(a(6),:).*A(:)'.*V(a(6),:),'k<-')
axis equal
legend('0dt','50dt','100dt','150dt','200dt','700dt')
figure
[ax,h1,h2]=plotyy(x,M(a(6),:),x,sqrt((M(a(6),:).^(1-r)-1)*5),'plot');
set(h1,'linestyle','-','marker','o','color','b');
set(h2,'linestyle','-','marker','*','color','r');
set(ax,'xtick',0:0.3:3)
set(ax(1),'ytick',0:0.1:1)
set(ax(2),'ytick',0:0.4:4)
set(h1,'linewidth',1.5)
set(h2,'linewidth',1.5)
set(get(ax(1),'xlabel'),'string','无量纲轴向距离(x)') set(get(ax(1),'ylabel'),'string','无量纲密度(M)') set(get(ax(2),'ylabel'),'string','马赫数(Ma)')
title('无量纲密度和马赫数定常分布数值解') legend('无量纲密度数值解','马赫数数值解')。