二维超声速流动的数值解——普朗特迈耶稀疏波.docx

二维超声速流动的数值解——普朗特迈耶稀疏波.docx

《二维超声速流动的数值解——普朗特迈耶稀疏波.docx》由会员分享，可在线阅读，更多相关《二维超声速流动的数值解——普朗特迈耶稀疏波.docx（9页珍藏版）》请在冰点文库上搜索。

二维超声速流动的数值解——普朗特迈耶稀疏波

主程序数值计算

%二维超声速流动的数值解——普朗特迈耶稀疏波

clc;

clear;

%常量设置

CFL=0.5;

Cy=0.2;

r=1.4;

R=287;

%初始物理条件

Ma0=2;

P0=1.01*10^5;

M0=1.23;

T0=286.1;

%流场几何外形描述

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）.^2+P（:

1）;

F3=M（:

1）.*U（:

1）.*V（:

1）;

F4=r/（r-1）*P（:

1）.*U（:

1）+M（:

1）.*U（:

1）.*（U（:

1）.^2+V（:

1）.^2）/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）.^2+F2-F1.^2./M（:

k）;

G4=r/（r-1）*（F2-F1.^2./M（:

k））.*F3./F1+M（:

k）/2.*F3./F1.*（（F1./M（:

k））.^2+（F3./F1）.^2）;

%预估偏导数，前向差分

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）.^2./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）.^3;

M_1（1:

j-1）=（-B（1:

j-1）+sqrt（B（1:

j-1）.^2-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）.^2./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））.^2+F2_1（1:

j-1）-F1_1（1:

j-1）.^2./M_1（1:

j-1）;

G4_1（1:

j-1）=r/（r-1）*（F2_1（1:

j-1）-F1_1（1:

j-1）.^2./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））.^2+（F3_1（1:

j-1）./F1_1（1:

j-1））.^2）;

%修正偏导数计算

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）=0.5*（F1_eps_1（2:

j-1）+F1_eps（2:

j-1））;

F2_eps_av（2:

j-1）=0.5*（F2_eps_1（2:

j-1）+F2_eps（2:

j-1））;

F3_eps_av（2:

j-1）=0.5*（F3_eps_1（2:

j-1）+F3_eps（2:

j-1））;

F4_eps_av（2:

j-1）=0.5*（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=F3.^2./F1/2-F4;

B=r/（r-1）*F1.*F2;

C=-（r+1）/（r-1）/2*F1.^3;

M（:

k+1）=（-B+sqrt（B.^2-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）.^2+V（:

k+1）.^2）./sqrt（r*R*T（:

k+1））;

%下边界处物理量修正

Ma_cal=sqrt（U（1,k+1）^2+V（1,k+1）^2）/sqrt（r*R*T（1,k+1））;

f_cal=sqrt（（r+1）/（r-1））*atan（sqrt（（r-1）*（Ma_cal^2-1）/（r+1）））-atan（sqrt（Ma_cal^2-1））;

ifk

fi=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=[12.55];

n=0;

whileabs（xm

（1）-xm（3））>=0.00001&&n<100

fm=sqrt（（r+1）/（r-1））*atan（sqrt（（r-1）/（r+1）*（xm.^2-1）））-atan（sqrt（xm.^2-1））-f_act;

iffm

（2）*fm（3）<0

xm=[xm

（2）（xm

（2）+xm（3））/2xm（3）];

elseiffm

（2）*fm（3）>0

xm=[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_cal^2）/（1+（r-1）/2*Ma（1,k+1）^2））^（r/（r-1））;

T（1,k+1）=T（1,k+1）*（（1+（r-1）/2*Ma_cal^2）/（1+（r-1）/2*Ma（1,k+1）^2））;

M（1,k+1）=P（1,k+1）/R/T（1,k+1）;

ifk

V（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）^2+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）^2+V（1,k+1）^2）/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（[0.51.5-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（[0.1*10^52*10^5-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

计算结果（粘性0.2+角度5度）