维超声速流动的数值解普朗特迈耶稀疏波.docx
《维超声速流动的数值解普朗特迈耶稀疏波.docx》由会员分享，可在线阅读，更多相关《维超声速流动的数值解普朗特迈耶稀疏波.docx（15页珍藏版）》请在冰点文库上搜索。
维超声速流动的数值解普朗特迈耶稀疏波
二维超声速流动的数值解——普朗特迈耶稀疏波
主程序数值计算
%二维超声速流动的数值解——普朗特迈耶稀疏波
clc;
clear;
%常量设置
CFL=;
Cy=;
r=;
R=287;
%初始物理条件
Ma0=2;
P0=*10A5;
M0=;
T0=;
%流场几何外形描述sita=5*pi/180;
H=40;
L=65;
%网格节点描述
i=6501;
j=4001;
x=linspace（0,L,i）;
point1=length（find（x<=10））;
ys=-（x-10）*tan（sita）;
ys（1:
point1）=0;
h=H-ys;
h（1:
point1）=H;
y（j,i）=0;
fork=1:
i
y（:
k）=linspace（ys（k）,H,j）';
end
%贴体网格生成
eps=x;
eta=（y-repmat（ys,j,1））./（repmat（H+abs（ys）,j,1））;
eta_x（j,i）=0;
eta_x（:
point1+1:
end）=（1-eta（:
point1+1:
end））./repmat（h（point1+1:
end）,j,1）*tan（sita）;deta=1/（j-1）;
%初始条件-入口条件
U（j,i）=0;
V（j,i）=0;
M（j,i）=0;
P（j,i）=0;
T（j,i）=0;
Ma（j,i）=0;
U（:
1）=Ma0*sqrt（r*R*T0）;
V（:
1）=0;
M（:
1）=M0;
P（:
1）=P0;
T（:
1）=T0;
Ma（:
1）=Ma0;
F1=M（:
1）.*U（:
1）;
F2=M（:
1）.*U（:
1）.A2+P（:
1）;
F3=M（:
1）.*U（:
1）.*V（:
1）;
F4=r/（r-1）*P（:
1）.*U（:
1）+M（:
1）.*U（:
1）.*（U（:
1）.A2+V（:
1）.A2）/2;
%计算初始化话工作，其实也可以不需要的
F1_eps=zeros（j-1,1）;
F2_eps=zeros（j-1,1）;
F3_eps=zeros（j-1,1）;
F4_eps=zeros（j-1,1）;
SF1=zeros（j-1,1）;
SF2=zeros（j-1,1）;
SF3=zeros（j-1,1）;
SF4=zeros（j-1,1）;
F1_1=zeros（j-1,1）;
F2_1=zeros（j-1,1）;
F3_1=zeros（j-1,1）;
F4_1=zeros（j-1,1）;
SF1_1=zeros（j,1）;
SF2_1=zeros（j,1）;
SF3_1=zeros（j,1）;
SF4_1=zeros（j,1）;
A=zeros（j,1）;
B=zeros（j,1）;
C=zeros（j,1）;
M_1=zeros（j,1）;
P_1=zeros（j,1）;
G1_1=zeros（j-1,1）;
G2_1=zeros（j-1,1）;
G3_1=zeros（j-1,1）;
G4_1=zeros（j-1,1）;
F1_eps_1=zeros（j-1,1）;
F2_eps_1=zeros（j-1,1）;
F3_eps_1=zeros（j-1,1）;
F4_eps_1=zeros（j-1,1）;
F1_eps_av=zeros（j-1,1）;
F2_eps_av=zeros（j-1,1）;
F3_eps_av=zeros（j-1,1）;
F4_eps_av=zeros（j-1,1）;
%麦考马克方法空间推进计算fork=1:
i
G1=M（:
k）.*F3./F1;
G2=F3;
G3=M（:
k）.*（F3./F1）.A2+F2-F142./M（:
k）;
G4=r/（r-1）*（F2-F192./M（:
k））.*F3./F1+M（:
k）/2.*F3./F1.*（（F1./M（:
k））.A2+（F3./F1）42）;
%预估偏导数，前向差分
F1_eps（1:
j-1）=eta_x（1:
j-1,k）.*（F1（1:
j-1）-F1（2:
j））/deta+1/h（k）*（G1（1:
j-1）-G1（2:
j））/deta;
F2_eps（1:
j-1）=eta_x（1:
j-1,k）.*（F2（1:
j-1）-F2（2:
j））/deta+1/h（k）*（G2（1:
j-1）-G2（2:
j））/deta;
F3_eps（1:
j-1）=eta_x（1:
j-1,k）.*（F3（1:
j-1）-F3（2:
j））/deta+1/h（k）*（G3（1:
j-1）-G3（2:
j））/deta;
F4_eps（1:
j-1）=eta_x（1:
j-1,k）.*（F4（1:
j-1）-F4（2:
j））/deta+1/h（k）*（G4（1:
j-1）-G4（2:
j））/deta;
%计算空间空间步长计算
miu=asin（1./Ma（1:
j-1,k））;
deps1=abs（tan（sita+miu））;
deps2=abs（tan（sita-miu））;
deps=CFL*h（k）/（j-1）/max（[deps1;deps2]）;
%人工粘性预估值
SF1（2:
j-1）=Cy*abs（P（3:
j,k）-2*P（2:
j-1,k）+P（1:
j-2,k））./（P（3:
j,k）+2*P（2:
j-1,k）+P（1:
j-2,k））.*（F1（3:
j）-2*F1（2:
j-1）+F1（1:
j-2））;
SF2（2:
j-1）=Cy*abs（P（3:
j,k）-2*P（2:
j-1,k）+P（1:
j-2,k））./（P（3:
j,k）+2*P（2:
j-1,k）+P（1:
j-2,k））.*（F2（3:
j）-2*F2（2:
j-1）+F2（1:
j-2））;
SF3（2:
j-1）=Cy*abs（P（3:
j,k）-2*P（2:
j-1,k）+P（1:
j-2,k））./（P（3:
j,k）+2*P（2:
j-1,k）+P（1:
j-2,k））.*（F3（3:
j）-2*F3（2:
j-1）+F3（1:
j-2））;
SF4（2:
j-1）=Cy*abs（P（3:
j,k）-2*P（2:
j-1,k）+P（1:
j-2,k））./（P（3:
j,k）+2*P（2:
j-1,k）+P（1:
j-2,k））.*（F4（3:
j）-2*F4（2:
j-1）+F4（1:
j-2））;SF1
（1）=2*SF1
（2）-SF1（3）;
SF2
（1）=2*SF2
（2）-SF2（3）;
SF3
（1）=2*SF3
（2）-SF3（3）;
SF4
（1）=2*SF4
（2）-SF4（3）;
%预估值
F1_1（1:
j-1）=F1（1:
j-1）+F1_eps（1:
j-1）*deps+SF1（1:
j-1）;
F2_1（1:
j-1）=F2（1:
j-1）+F2_eps（1:
j-1）*deps+SF2（1:
j-1）;
F3_1（1:
j-1）=F3（1:
j-1）+F3_eps（1:
j-1）*deps+SF3（1:
j-1）;
F4_1（1:
j-1）=F4（1:
j-1）+F4_eps（1:
j-1）*deps+SF4（1:
j-1）;
A（1:
j-1）=F3_1（1:
j-1）42./F1_1（1:
j-1）/2-F4_1（1:
j-1）;
B（1:
j-1）=r/（r-1）*F1_1（1:
j-1）.*F2_1（1:
j-1）;
C（1:
j-1）=-（r+1）/（r-1）/2*F1_1（1:
j-1）.A3;
M_1（1:
j-1）=（-B（1:
j-1）+sqrt（B（1:
j-1）.A2-4*A（1:
j-1）.*C（1:
j-1）））./A（1:
j-1）/2;
P_1（1:
j-1）=F2_1（1:
j-1）-F1_1（1:
j-1）.A2./M_1（1:
j-1）;
G1_1（1:
j-1）=M_1（1:
j-1）.*F3_1（1:
j-1）./F1_1（1:
j-1）;
G2_1（1:
j-1）=F3_1（1:
j-1）;
G3_1（1:
j-1）=M_1（1:
j-1）.*（F3_1（1:
j-1）./F1_1（1:
j-1））.A2+F2_1（1:
j-1）-F1_1（1:
j-1）.A2./M_1（1:
j-1）;
G4_1（1:
j-1）=r/（r-1）*（F2_1（1:
j-1）-F1_1（1:
j-1）.A2./M_1（1:
j-1））.*F3_1（1:
j-1）./F1_1（1:
j-1）+...
M_1（1:
j-1）/2.*F3_1（1:
j-1）./F1_1（1:
j-1）.*（（F1_1（1:
j-1）./M_1（1:
j-1））.A2+（F3_1（1:
j-1）./F1_1（1:
j-1））.A2）;
%修正偏导数计算
F1_eps_1（2:
j-1）=eta_x（2:
j-1,k）.*（F1_1（1:
j-2）-F1_1（2:
j-1））/deta+1/h（k）*（G1_1（1:
j-2）-G1_1（2:
j-1））/deta;
F2_eps_1（2:
j-1）=eta_x（2:
j-1,k）.*（F2_1（1:
j-2）-F2_1（2:
j-1））/deta+1/h（k）*（G2_1（1:
j-2）-G2_1（2:
j-1））/deta;
F3_eps_1（2:
j-1）=eta_x（2:
j-1,k）.*（F3_1（1:
j-2）-F3_1（2:
j-1））/deta+1/h（k）*（G3_1（1:
j-2）-G3_1（2:
j-1））/deta;
F4_eps_1（2:
j-1）=eta_x（2:
j-1,k）.*（F4_1（1:
j-2）-F4_1（2:
j-1））/deta+1/h（k）*（G4_1（1:
j-2）-G4_1（2:
j-1））/deta;
F1_eps_av（2:
j-1）=*（F1_eps_1（2:
j-1）+F1_eps（2:
j-1））;
F2_eps_av（2:
j-1）=*（F2_eps_1（2:
j-1）+F2_eps（2:
j-1））;
F3_eps_av（2:
j-1）=*（F3_eps_1（2:
j-1）+F3_eps（2:
j-1））;
F4_eps_av（2:
j-1）=*（F4_eps_1（2:
j-1）+F4_eps（2:
j-1））;
%人工粘性修正值
SF1_1（2:
j-2）=Cy*abs（P_1（3:
j-1）-2*P_1（2:
j-2）+P_1（1:
j-3））./（P_1（3:
j-1）+2*P_1（2:
j-2）+P_1（1:
j-3））.*（F1_1（3:
j-1）-2*F1_1（2:
j-2）+F1_1（1:
j-3））;
SF2_1（2:
j-2）=Cy*abs（P_1（3:
j-1）-2*P_1（2:
j-2）+P_1（1:
j-3））./（P_1（3:
j-1）+2*P_1（2:
j-2）+P_1（1:
j-3））.*（F2_1（3:
j-1）-2*F2_1（2:
j-2）+F2_1（1:
j-3））;
SF3_1（2:
j-2）=Cy*abs（P_1（3:
j-1）-2*P_1（2:
j-2）+P_1（1:
j-3））./（P_1（3:
j-1）+2*P_1（2:
j-2）+P_1（1:
j-3））.*（F3_1（3:
j-1）-2*F3_1（2:
j-2）+F3_1（1:
j-3））;
SF4_1（2:
j-2）=Cy*abs（P_1（3:
j-1）-2*P_1（2:
j-2）+P_1（1:
j-3））./（P_1（3:
j-1）+2*P_1（2:
j-2）+P_1（1:
j-3））.*（F4_1（3:
j-1）-2*F4_1（2:
j-2）+F4_1（1:
j-3））;
SF1_1（j-1）=2*SF1_1（j-2）-SF1_1（j-3）;
SF2_1（j-1）=2*SF2_1（j-2）-SF2_1（j-3）;
SF3_1（j-1）=2*SF3_1（j-2）-SF3_1（j-3）;
SF4_1（j-1）=2*SF4_1（j-2）-SF4_1（j-3）;
%修正值
F1（2:
j-1）=F1（2:
j-1）+F1_eps_av（2:
j-1）*deps+SF1_1（2:
j-1）;
F2（2:
j-1）=F2（2:
j-1）+F2_eps_av（2:
j-1）*deps+SF2_1（2:
j-1）;
F3（2:
j-1）=F3（2:
j-1）+F3_eps_av（2:
j-1）*deps+SF3_1（2:
j-1）;
F4（2:
j-1）=F4（2:
j-1）+F4_eps_av（2:
j-1）*deps+SF4_1（2:
j-1）;
%上边界单侧差分
F1（j）=2*F1（j-1）-F1（j-2）;
F2（j）=2*F2（j-1）-F2（j-2）;
F3（j）=2*F3（j-1）-F3（j-2）;
F4（j）=2*F4（j-1）-F4（j-2）;
%下边界单侧差分
F1
（1）=2*F1
（2）-F1（3）;
F2
（1）=2*F2
（2）-F2（3）;
F3
（1）=2*F3
（2）-F3（3）;
F4
（1）=2*F4
（2）-F4（3）;
%所有网格点处物理量计算
A=F342./F1/2-F4;
B=r/（r-1）*F1.*F2;
C=-（r+1）/（r-1）/2*F193;
M（:
k+1）=（-B+sqrt（B.A2-4*A.*C））./A/2;
U（:
k+1）=F1./M（:
k+1）;
V（:
k+1）=F3./F1;
P（:
k+1）=F2-F1.*U（:
k+1）;
T（:
k+1）=P（:
k+1）./M（:
k+1）/R;
Ma（:
k+1）=sqrt（U（:
k+1）.A2+V（:
k+1）.A2）./sqrt（r*R*T（:
k+1））;
%下边界处物理量修正
Ma_cal=sqrt（U（1,k+1）A2+V（1,k+1）A2）/sqrt（r*R*T（1,k+1））;f_cal=sqrt（（r+1）/（r-1））*atan（sqrt（（r-1）*（Ma_calA2-1）/（r+1）））-atan（sqrt（Ma_calA2-1））;
ifkfi=atan（abs（V（1,k+1））/U（1,k+1））;
else
psai=atan（abs（V（1,k+1））/U（1,k+1））;
fi=sita-psai;
end
f_act=f_cal+fi;
%二分法求马赫数
xm=[15];n=0;
whileabs（xm
（1）-xm（3））>=&&n<100
fm=sqrt（（r+1）/（r-1））*atan（sqrt（（r-1）/（r+1）*（xm.A2-1）））-atan（sqrt（xm.A2-1））-f_act;iffm
（2）*fm（3）<0
xm=[xm
（2）（xm
（2）+xm（3））/2xm（3）];
elseiffm
（2）*fm（3）>0xm=[xm
（1）（xm
（1）+xm
（2））/2xm
（2）];
else
xm=[xm
（2）xm
（2）xm
（2）];
end
n=n+1;
end
Ma（1,k+1）=xm
（2）;
%下边界处物理量修正值
P（1,k+1）=P（1,k+1）*（（1+（r-1）/2*Ma_calA2）/（1+（r-1）/2*Ma（1,k+1）A2））A（r/（r-1））;T（1,k+1）=T（1,k+1）*（（1+（r-1）/2*Ma_calA2）/（1+（r-1）/2*Ma（1,k+1）A2））;
M（1,k+1）=P（1,k+1）/R/T（1,k+1）;
ifkV（1,k+1）=0;
else
V（1,k+1）=-U（1,k+1）*tan（sita）;
end
F1
（1）=M（1,k+1）*U（1,k+1）;F2
（1）=M（1,k+1）*U（1,k+1）A2+P（1,k+1）;
F3
（1）=M（1,k+1）*U（1,k+1）*V（1,k+1）;
F4
（1）=r/（r-1）*P（1,k+1）*U（1,k+1）+M（1,k+1）*U（1,k+1）*（U（1,k+1）A2+V（1,k+1）A2）/2;k
end
计算结果后处理
closeall;
%横向速度
point=5*floor（1:
20:
130）;
figure
fori=1:
length（point）
subplot（1,length（point）,i）plot（U（:
point（i））,y（:
point（i））,'-','linewidth',2）ylabel（'y/cm'）
xlabel（'u/（m/s）'）
title（strcat（'x=',num2str（point（i））））
axis（[660760-1540]）
%axisoff
end
%纵向速度
point=5*floor（1:
20:
130）;
figure
fori=1:
length（point）
subplot（1,length（point）,i）plot（V（:
point（i））,y（:
point（i））,'-','linewidth',2）ylabel（'y/cm'）
xlabel（'v/（m/s）'）title（strcat（'x=',num2str（point（i））））axis（[-25050-1540]）
%axisoff
end
%马赫数
point=5*floor（1:
20:
130）;
figure
fori=1:
length（point）
subplot（1,length（point）,i）plot（Ma（:
point（i））,y（:
point（i））,'-','linewidth',2）ylabel（'y/cm'）
xlabel（'Ma'）title（strcat（'x=',num2str（point（i））））axis（[13-1540]）
%axisoff
end
%密度
point=5*floor（1:
20:
130）;
figure
fori=1:
length（point）
subplot（1,length（point）,i）plot（M（:
point（i））,y（:
point（i））,'-','linewidth',2）ylabel（'y/cm'）
xlabel（'M/（kg/m3）'）title（strcat（'x=',num2str（point（i））））
axis（[-1540]）
%axisoff
end
%压力
point=5*floor（1:
20:
130）;
figure
fori=1:
length（point）
subplot（1,length（point）,i）plot（P（:
point（i））,y（:
point（i））,'-','linewidth',2）ylabel（'y/cm'）
xlabel（'P/（Pa）'）title（strcat（'x=',num2str（point（i））））
axis（[*10A52*10A5-1540]）
%axisoff
end
%温度
point=5*floor（1:
20:
130）;
figure
fori=1:
length（point）
subplot（1,length（point）,i）
plot（T（:
point（i））,y（:
point（i））,'-','linewidth',2）ylabel（'y/cm'）
xlabel（'T/（k）'）
title（strcat（'x=',num2str（point（i））））axis（[200300-1540]）
%axisoff
End
计算结果（粘性+角度5度）
0
■
■
■
-
■-
-
1-5'■
-
1□-
-
5■
-
0•
&1:
?
£'Ve
rwu
X
30
M
:
d
O
r
P：
-E1
D
JO
2Q
2D
W：
1
呼ill」」
HQ
■-»
Wlkii'irU
hl他旷曰
F*V|k>.||l!
>」
trihJi
讯化屮di
»-1k-41»-IO伽TAMg映）TAh）Wo
•-R
X.
升级会员 