%亞音速-超音速，利用U1，U2，U3守恒形式，激波捕捉_

%Agammadtdx
clearclc
NUM=61;
C=0.5;%Courant?number
dx=3/（NUM-1）;
step=1400;%時(shí)間模擬的步數(shù)
Time=0;
gamma=1.4;
p_e=0.6784;?%出口對總壓的比值，此程序模擬激波
x=0:dx:3;
Cx=0.2;%人工粘性常數(shù)?0.01~0.3


A=1+2.2*（x-1.5）.^2;
dlnA_dx=4.4*（x-1.5）./A;%對A求導(dǎo)
%?dlnA_dx=zeros（1NUM）;
%?for?i=2:1:NUM
%?dlnA_dx（i）=（log（A（i））-log（A（i-1）））/dx;
%?end
%?dlnA_dx（1）=dlnA_dx（2）;


rou_ar=zeros（step+1NUM）;V_ar=zeros（step+1NUM）;T_ar=zeros（step+1NUM）;
Ma_ar=zeros（step+1NUM）;flux_m=zeros（step+1NUM）;P_ar=zeros（step+1NUM）;
%?drou_t=zeros（1NUM）;dV_t=zeros（1NUM）;dT_t=zeros（1NUM）;
%?drou_t2=zeros（1NUM）;dV_t2=zeros（1NUM）;dT_t2=zeros（1NUM）;
%?drou_t3=zeros（1NUM）;dV_t3=zeros（1NUM）;dT_t3=zeros（1NUM）;
%?rou_e=zeros（1NUM）;V_e=zeros（1NUM）;T_e=zeros（1NUM）;
rou_r=zeros（1NUM）;V_r=zeros（1NUM）;T_r=zeros（1NUM）;
dt_ar=zeros（1step）;
rou=zeros（1NUM）;T=zeros（1NUM）;
dU1_t=zeros（1NUM）;dU2_t=zeros（1NUM）;dU3_t=zeros（1NUM）;
U1_e=zeros（1NUM）;U2_e=zeros（1NUM）;U3_e=zeros（1NUM）;
dU1_t2=zeros（1NUM）;dU2_t2=zeros（1NUM）;dU3_t2=zeros（1NUM）;
dU1_t3=zeros（1NUM）;dU2_t3=zeros（1NUM）;dU3_t3=zeros（1NUM）;
U1_r=zeros（1NUM）;U2_r=zeros（1NUM）;U3_r=zeros（1NUM）;
U1_ar=zeros（step+1NUM）;U2_ar=zeros（step+1NUM）;U3_ar=zeros（step+1NUM）;
S1_e=zeros（1NUM）;S2_e=zeros（1NUM）;S3_e=zeros（1NUM）;
S1_r=zeros（1NUM）;S2_r=zeros（1NUM）;S3_r=zeros（1NUM）;

%初始化
for?i=1:1:NUM
????if?x（i）>=0?&&?x（i）<=0.5
???????rou（i）=1;
???????T（i）=1;
????elseif?x（i）>0.5?&&?x（i）<=1.5
????????rou（i）=1-0.366*（x（i）-0.5）;
????????T（i）=1-0.167*（x（i）-0.5）;
????elseif?x（i）>1.5?&&?x（i）<=2.1
????????rou（i）=0.634-0.702*（x（i）-1.5）;
????????T（i）=0.833-0.4908*（x（i）-1.5）;
????elseif?x（i）>2.1?&&?x（i）<=3.0
????????rou（i）=0.5892+0.10228*（x（i）-2.1）;
????????T（i）=0.93968+0.0622*（x（i）-2.1）;
????end
end
V=0.59./（rou.*A）;??%0.59是更配合穩(wěn)態(tài)解的數(shù)值。

%
rou_ar（1:）=rou;
V_ar（1:）=V;
T_ar（1:）=T;
P_ar（1:）=rou.*T;
Ma_ar（1:）=V./sqrt（T）;
flux_m（1:）=rou.*V.*A;

for?k=1:1:step???%模擬的時(shí)間步數(shù)

rou=rou_ar（k:）;
V=V_ar（k:）;
T=T_ar（k:）;
P=P_ar（k:）;

U1=rou.*A;
U2=rou.*V.*A;
U3=rou.*A.*（（T/（gamma-1））+（gamma/2）*V.^2）;
F1=U2;
F2=（U2.^2）./U1+（gamma-1）/gamma*（U3-（gamma/2）*（U2.^2）./U1）;
F3=gamma*U2.*U3./U1-（gamma*（gamma-1）/2）*（（U2.^3）./（U1.^2））;
J2=（gamma-1）/gamma*（U3-gamma/2*（U2.^2）./U1）.*（dlnA_dx）;


%確定本次模擬的dt時(shí)間間隔
dt_list=C*dx./（T.^0.5+V）;
dt=min（dt_list）;
%
%預(yù)測步
for?i=2:1:（NUM-1）

????dU1_t（i）=-（F1（i+1）-F1（i））/dx;
????dU2_t（i）=-（F2（i+1）-F2（i））/dx+J2（i）;
????dU3_t（i）=-（F3（i+1）-F3（i））/dx;
????
????S1_e（i）=Cx*（（abs（P（i+1）-2*P（i）+P（i-1）））/（P（i+1）+2*P（i）+P（i-1）））*（U1（i+1）-2*U1（i）+U1（i-1））;
????S2_e（i）=Cx*（（abs（P（i+1）-2*P（i）+P（i-1）））/（P（i+1）+2*P（i）+P（i-1）））*（U2（i+1）-2*U2（i）+U2（i-1））;
????S3_e（i）=Cx*（（abs（P（i+1）-2*P（i）+P（i-1）））/（P（i+1）+2*P（i）+P（i-1）））*（U3（i+1）-2*U3（i）+U3（i-1））;
????
????U1_e（i）=U1（i）+dU1_t（i）*dt+S1_e（i）;
????U2_e（i）=U2（i）+dU2_t（i）*dt+S2_e（i）;
????U3_e（i）=U3（i）+dU3_t（i）*dt+S3_e（i）;
????
end

%BC
U1_e（1）=A（1）;
U2_e（1）=2*U2_e（2）-U2_e（3）;

?屬性????????????大小?????日期????時(shí)間???名稱
-----------?---------??----------?-----??----

?????文件???????6109??2012-12-09?18:40??Chapter7\one_D_shockwave_4.m

?????文件???????3934??2012-12-08?18:12??Chapter7\one_D_subsonic_all_2.m

?????文件???????3402??2012-12-09?00:46??Chapter7\one_D_supersonic_all_1.m

?????文件???????5056??2012-12-09?15:50??Chapter7\one_D_supersonic_all_3.m

?????文件????????388??2012-12-09?15:43??Chapter7\Untitled.m

?????文件???????5033??2012-12-09?00:23??Chapter7\Untitled3.m

?????文件???????1208??2012-12-08?11:19??Chapter7\visualize_1.m

?????文件???????1450??2012-12-08?18:04??Chapter7\visualize_2.m

?????文件???????1213??2012-12-09?00:28??Chapter7\visualize_3.m

?????文件???????1218??2012-12-09?18:43??Chapter7\visualize_4.m

?????目錄??????????0??2012-12-10?23:15??Chapter7

-----------?---------??----------?-----??----

????????????????29011????????????????????11