%亚音速-超音速,利用U1,U2,U3守恒形式,激波捕捉_
%A,gamma,dt,dx
clear,clc
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(1,NUM);
% 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+1,NUM);V_ar=zeros(step+1,NUM);T_ar=zeros(step+1,NUM);
Ma_ar=zeros(step+1,NUM);flux_m=zeros(step+1,NUM);P_ar=zeros(step+1,NUM);
% drou_t=zeros(1,NUM);dV_t=zeros(1,NUM);dT_t=zeros(1,NUM);
% drou_t2=zeros(1,NUM);dV_t2=zeros(1,NUM);dT_t2=zeros(1,NUM);
% drou_t3=zeros(1,NUM);dV_t3=zeros(1,NUM);dT_t3=zeros(1,NUM);
% rou_e=zeros(1,NUM);V_e=zeros(1,NUM);T_e=zeros(1,NUM);
rou_r=zeros(1,NUM);V_r=zeros(1,NUM);T_r=zeros(1,NUM);
dt_ar=zeros(1,step);
rou=zeros(1,NUM);T=zeros(1,NUM);
dU1_t=zeros(1,NUM);dU2_t=zeros(1,NUM);dU3_t=zeros(1,NUM);
U1_e=zeros(1,NUM);U2_e=zeros(1,NUM);U3_e=zeros(1,NUM);
dU1_t2=zeros(1,NUM);dU2_t2=zeros(1,NUM);dU3_t2=zeros(1,NUM);
dU1_t3=zeros(1,NUM);dU2_t3=zeros(1,NUM);dU3_t3=zeros(1,NUM);
U1_r=zeros(1,NUM);U2_r=zeros(1,NUM);U3_r=zeros(1,NUM);
U1_ar=zeros(step+1,NUM);U2_ar=zeros(step+1,NUM);U3_ar=zeros(step+1,NUM);
S1_e=zeros(1,NUM);S2_e=zeros(1,NUM);S3_e=zeros(1,NUM);
S1_r=zeros(1,NUM);S2_r=zeros(1,NUM);S3_r=zeros(1,NUM);
%初始化
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);
V_e(1)=U2_e(1)/U1_e(1);
U3_e(1)=U1_e(1)*(1/(gamma-1)+(gamma/2)*(V_e(1).^2));
U1_e(NUM)=2*U1_e(NUM-1)-U1_e(NUM-2);
U2_e(NUM)=2*U2_e(NUM-1)-U2_e(NUM-2);
V_e_N=U2_e(NUM)/U1_e(NUM);
U3_e(NUM)=(p_e*A(NUM))/(gamma-1)+(gamma/2)*U2_e(NUM)*V_e_N;
%原始变量
rou_e=U1_e./A;
V_e=U2_e./U1_e;
T_e=(gamma-1)*(U3_e./U1_e-(gamma/2)*(V_e.^2));
P_e=rou_e.*T_e;
%F变量
F1_e=U2_e;
F2_e=(U2_e.^2)./U1_e+(gamma-1)/gamma*(U3_e-(gamma/2)*(U2_e.^2)./U1_e);
F3_e=gamma*U2_e.*U3_e./U1_e-(gamma*(gamma-1)/2)*((U2_e.^3)./(U1_e.^2));
J2_e=(gamma-1)/gamma*(U3_e-gamma/2*(U2_e.^2)./U1_e).*(dlnA_dx);
%%%%%%%%%%%%
%预测步结束
%修正步
for i=2:1:(NUM-1)
dU1_t2(i)=-(F1_e(i)-F1_e(i-1))/dx;
dU2_t2(i)=-(F2_e(i)-F2_e(i-1))/dx+J2_e(i);
dU3_t2(i)=-(F3_e(i)-F3_e(i-1))/dx;
dU1_t3(i)=0.5*(dU1_t(i)+dU1_t2(i));
dU2_t3(i)=0.5*(dU2_t(i)+dU2_t2(i));
dU3_t3(i)=0.5*(dU3_t(i)+dU3_t2(i));
S1_r(i)=Cx*((abs(P_e(i+1)-2*P_e(i)+P_e(i-1)))/(P_e(i+1)+2*P_e(i)+P_e(i-1)))*(U1_e(i+1)-2*U1_e(i)+U1_e(i-1));
S2_r(i)=Cx*((abs(P_e(i+1)-2*P_e(i)+P_e(i-1)))/(P_e(i+1)+2*P_e(i)+P_e(i-1)))*(U2_e(i+1)-2*U2_e(i)+U2_e(i-1));
S3_r(i)=Cx*((abs(P_e(i+1)-2*P_e(i)+P_e(i-1)))/(P_e(i+1)+2*P_e(i)+P_e(i-1)))*(U3_e(i+1)-2*U3_e(i)+U3_e(i-1));
U1_r(i)=U1(i)+dU1_t3(i)*dt+S1_r(i);
U2_r(i)=U2(i)+dU2_t3(i)*dt+S2_r(i);
U3_r(i)=U3(i)+dU3_t3(i)*dt+S3_r(i);
end
%BC
U1_r(1)=A(1);
U2_r(1)=2*U2_r(2)-U2_r(3);
V_r(1)=U2_r(1)/U1_r(1);
U3_r(1)=U1_r(1)*(1/(gamma-1)+(gamma/2)*(V_r(1).^2));
U1_r(NUM)=2*U1_r(NUM-1)-U1_r(NUM-2);
U2_r(NUM)=2*U2_r(NUM-1)-U2_r(NUM-2);
V_r_N=U2_r(NUM)/U1_r(NUM);
U3_r(NUM)=(p_e*A(NUM))/(gamma-1)+(gamma/2)*U2_r(NUM)*V_r_N;
%原始变量
rou_r=U1_r./A;
V_r=U2_r./U1_r;
T_r=(gamma-1)*(U3_r./U1_r-(gamma/2)*(V_r.^2));
%F变量
% F1_r=U2_r;
% F2_r=(U2_r.^2)/U1_r+(gamma-1)/gamma.*(U3_r-(gamma/2)*(U2_r.^2)./U1_r);
% F3_r=gamma*U2_r.*U3_r./U1_r-(gamma*(gamma-1)/2)*((U2_r.^3)/(U1_r.^2));
% J2_r=(gamma-1)/gamma*(U3_r-gamma/2*(U2_r.^2)./U1_r).*(dlnA_dx);
%%%%%%%%%%
U1_ar(k+1,:)=U1_r;
U2_ar(k+1,:)=U2_r;
U3_ar(k+1,:)=U3_r;
rou_ar(k+1,:)=rou_r;
V_ar(k+1,:)=V_r;
T_ar(k+1,:)=T_r;
Ma_ar(k+1,:)=V_r./sqrt(T_r);
flux_m(k+1,:)=rou_r.*V_r.*A;
P_ar(k+1,:)=rou_r.*T_r;
dt_ar(k)=dt;
Time=Time+dt;
end %模拟的时间步数
% %理论解的计算
% x_theory=0:0.01:3;
% A_theory=1+2.2*(x_theory-1.5).^2;
% num_x_th=size(x_theory,2);
% Ma_theory=zeros(1,num_x_th);
% %求解理论上的Ma解
% for i=1:1:num_x_th
% F=@(ma)(1/ma^2)*((2/(gamma+1))*(1+(gamma-1)/2*ma^2))^((gamma+1)/(gamma-1))-A_theory(i)^2;
% if i<=(num_x_th-1)/2+1
% Ma_theory(i)=fzero(F,[1e-8,1]);
% else
% Ma_theory(i)=fzero(F,[1,10]);
% end
% end
%
% P_theory=(1+((gamma-1)/2)*(Ma_theory(end,:).^2)).^(-gamma/(gamma-1));
% rou_theory=(1+((gamma-1)/2)*(Ma_theory(end,:).^2)).^(-1/(gamma-1));
% T_theory=(1+((gamma-1)/2)*(Ma_theory(end,:).^2)).^(-1);
%
% %模拟解
P_simu=P_ar(end,:);
rou_simu=rou_ar(end,:);
T_simu=T_ar(end,:);
Ma_simu=Ma_ar(end,:);
没有合适的资源?快使用搜索试试~ 我知道了~
Anderson CFD 代码 matlab
共39个文件
m:31个
asv:6个
jpg:2个
需积分: 50 94 下载量 165 浏览量
2018-12-07
10:10:25
上传
评论 9
收藏 100KB RAR 举报
温馨提示
anderson 计算流体力学书中的matlab代码,共同学习了。
资源推荐
资源详情
资源评论
收起资源包目录
CFD code in Matlab for anderson.rar (39个子文件)
Chapter9_Version_3
Visualize_1.m 855B
FTAOXY.m 896B
MACORMACK.m 5KB
THERMC.m 64B
BC.m 2KB
FTAOYY.m 510B
untitled.jpg 55KB
DYNVIS.m 79B
FTAOXX.m 525B
CONVER.m 251B
MACORMACK.asv 4KB
plate_super_shock.m 4KB
FQY.m 338B
TSTEP.m 457B
plate_super_shock.asv 4KB
FQX.m 342B
FBC.m 1KB
Chapter7
visualize_3.m 1KB
Untitled.m 388B
visualize_1.m 1KB
one_D_supersonic_all_1.m 3KB
visualize_4.m 1KB
one_D_shockwave_4.m 6KB
one_D_subsonic_all_2.m 4KB
one_D_supersonic_all_3.m 5KB
Untitled3.m 5KB
visualize_2.m 1KB
Chapter8
Pressure_revised.jpg 32KB
Pressure_revised_2.m 5KB
visualize_2.asv 349B
visualize_1.m 508B
plot_look.m 262B
Thomas.m 395B
Crank_Nicolson_1.m 808B
Pressure_revised_2.asv 4KB
Crank-Nicolson_1.asv 370B
Crank_Nicolson_1.asv 637B
visualize_2.m 710B
Untitled2.m 86B
共 39 条
- 1
资源评论
rl900911
- 粉丝: 2
- 资源: 5
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功