拟一维喷管流动的数值解法（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_...
