clear
clc
tic
%-----Code?By?WangNianhua?2013.10.17
%?亞音速-超音速等熵噴管，非守恒形式方程MacCormack方法求解，教材P195-222
%---------------第一步，噴管形狀參數和初始條件--------
L=3;N=31;dx=L/（N-1）;C=0.5;
x=0:dx:L;?
A=1+2.2*（x-1.5）.^2;
subplot（221）plot（xA‘-‘）;
hold?on;
plot（x-A‘-‘）;
plot（xzeros（1N）‘.-‘）
title（‘噴管外形（網格點分布如圖）‘）
hold?off;
Rou=1-0.3146*x;
T=1-0.2314*x;
V=（0.1+1.09*x）.*sqrt（T）;
for?n=1:1400
????Y1=Rou（16）;
%-------------第二步，預估步------------------------
for?i=2:30
????P_rou_t（i）=-V（i）*（（Rou（i+1）-Rou（i））/dx）-Rou（i）*（（V（i+1）-V（i））/dx）-Rou（i）*V（i）*（（log（A（i+1））-log（A（i）））/dx）;
????P_v_t（i）=-V（i）*（（V（i+1）-V（i））/dx）-（（T（i+1）-T（i））/dx+（（Rou（i+1）-Rou（i））/dx）*T（i）/Rou（i））*1/1.4;
????P_T_t（i）=-V（i）*（（T（i+1）-T（i））/dx）-0.4*T（i）*（（（V（i+1）-V（i））/dx）+V（i）*（（log（A（i+1））-log（A（i）））/dx））;
end
%-------------第三步，求Δt求rou?V?T的預估值-----
dt=C*（dx./（V（2:30）+sqrt（T（2:30））））;
dt=min（dt）;
Rou1（2:30）=Rou（2:30）+P_rou_t（2:30）.*dt;
V1（2:30）=V（2:30）+P_v_t（2:30）.*dt;
T1（2:30）=T（2:30）+P_T_t（2:30）.*dt;
V1（1）=V（1）;T1（1）=T（1）;Rou1（1）=Rou（1）;
%--------------第四步，校正步--------------
for?i=2:30
????P_rou_t_2（i）=-V1（i）*（（Rou1（i）-Rou1（i-1））/dx）-Rou1（i）*（（V1（i）-V1（i-1））/dx）-Rou1（i）*V1（i）