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

拟一维喷管流动的数值解法（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); %入口边界值