Anderson计算流体力学入门第7章matlab程序

Anderson计算流体力学入门第7章,四个代码程序，分别对应四个visualize画图程序。参数可以任意修改。matlab编程。

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

%Agammadtdx
clearclc
NUM=61;
C=0.5;%Courant number
dx=3/（NUM-1）;
step=1400;%时间模拟的步数
Time=0;
gamma=1.4;
p_e=0.6784; %出口对总压的比值，此程序模拟激波
x=0:dx:3;
Cx=0.2;%人工粘性常数 0.01~0.3


A=1+2.2*（x-1.5）.^2;
dlnA_dx=4.4*（x-1.5）./A;%对A求导
% 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是更配合稳态解的数值。

%
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   %模拟的时间步数

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时间间隔
dt_list=C*dx./（T.^0.5+V）;
dt=min（dt_list）;
%
%预测步
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）;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       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