Anderson CFD 代码 matlab

anderson 计算流体力学书中的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）;

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

     文件        370  2012-12-10 23:47  Chapter8\Crank-Nicolson_1.asv

     文件        637  2012-12-11 00:23  Chapter8\Crank_Nicolson_1.asv

     文件        808  2012-12-12 20:33  Chapter8\Crank_Nicolson_1.m

     文件        262  2012-12-12 20:40  Chapter8\plot_look.m

     文件      32436  2012-12-17 08:40  Chapter8\Pressure_revised.jpg

     文件       3680  2012-12-12 12:47  Chapter8\Pressure_revised_2.asv

     文件       4874  2012-12-12 18:45  Chapter8\Pressure_revised_2.m

     文件        395  2012-12-10 22:26  Chapter8\Thomas.m

     文件         86  2012-12-12 14:09  Chapter8\Untitled2.m

     文件        508  2012-12-12 20:38  Chapter8\visualize_1.m

     文件        349  2012-12-12 14:28  Chapter8\visualize_2.asv

     文件        710  2012-12-12 18:46  Chapter8\visualize_2.m

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

     文件       4694  2012-12-15 20:17  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

     文件       2210  2012-12-16 18:25  Chapter9_Version_3\BC.m

     文件        251  2012-12-16 15:20  Chapter9_Version_3\CONVER.m

     文件         79  2012-12-15 23:51  Chapter9_Version_3\DYNVIS.m

     文件       1211  2012-12-17 09:22  Chapter9_Version_3\FBC.m

     文件        342  2012-12-16 13:50  Chapter9_Version_3\FQX.m

     文件        338  2012-12-16 13:52  Chapter9_Version_3\FQY.m

     文件        525  2012-12-16 13:34  Chapter9_Version_3\FTAOXX.m

     文件        896  2012-12-16 20:16  Chapter9_Version_3\FTAOXY.m

............此处省略15个文件信息