📄 hw-nozzle.m
字号:
else
Uh_B(:,i,j)=U_t_T(:,i,j);
Uh_T(:,i,j)=U_t_B(:,i+1,j);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function dt=tstep(Uh_L,Uh_R,Uh_B,Uh_T,para)
cfl=para.cfl;
gam=para.gam;
BTX=para.BTPX;
BTY=para.BTPY;
LRX=para.LRPX;
LRY=para.LRPY;
old_dt=1;
for i=2:para.N-2
for j=2:para.M-2
dx=LRX(i,j)-LRX(i,j-1);
D=sqrt(Uh_R(1,i,j-1)/Uh_L(1,i,j-1));
u=(Uh_L(2,i,j-1)/Uh_L(1,i,j-1)+D*Uh_R(2,i,j-1)/Uh_R(1,i,j-1))/(1+D);
v=(Uh_L(3,i,j-1)/Uh_L(1,i,j-1)+D*Uh_R(3,i,j-1)/Uh_R(1,i,j-1))/(1+D);
E=(Uh_L(4,i,j-1)/Uh_L(1,i,j-1)+D*Uh_R(4,i,j-1)/Uh_R(1,i,j-1))/(1+D);
a=((gam-1)*gam*(E-0.5*(u^2+v^2)))^0.5;
lamta_x(i,j)=sqrt(u^2+v^2)+a;
dt_x=dx*cfl/lamta_x(i,j);
dy=abs(BTY(i,j)-BTY(i-1,j));
D=sqrt(Uh_T(1,i-1,j)/Uh_B(1,i-1,j));
u=(Uh_B(2,i-1,j)/Uh_B(1,i-1,j)+D*Uh_T(2,i-1,j)/Uh_T(1,i-1,j))/(1+D);
v=(Uh_B(3,i-1,j)/Uh_B(1,i-1,j)+D*Uh_T(3,i-1,j)/Uh_T(1,i-1,j))/(1+D);
E=(Uh_B(4,i-1,j)/Uh_B(1,i-1,j)+D*Uh_T(4,i-1,j)/Uh_T(1,i-1,j))/(1+D);
a=((gam-1)*gam*(E-0.5*(u^2+v^2)))^0.5;
lamta_y(i,j)=sqrt(u^2+v^2)+a;
dt_y=dy*cfl/lamta_y(i,j);
dt=dt_x;
if dt_x>dt_y;
dt=dt_y;
end
if old_dt>dt
old_dt=dt;
end
end
end
dt=old_dt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [FhX,GhX,FhY,GhY]=GetFh_MUSCL(U,Uh_L,Uh_R,Uh_B,Uh_T,para)
M=para.M;
N=para.N;
gam=para.gam;
%%%
for i=1:N
for j=1:M-1
Fh_L(:,i,j)=[ Uh_L(2,i,j);
(gam-1)*Uh_L(4,i,j)+0.5*(3-gam)*(Uh_L(2,i,j))^2/Uh_L(1,i,j)-0.5*(gam-1)*(Uh_L(3,i,j))^2/Uh_L(1,i,j);
Uh_L(2,i,j)*Uh_L(3,i,j)/Uh_L(1,i,j);
gam*Uh_L(2,i,j)*Uh_L(4,i,j)/Uh_L(1,i,j)-0.5*(gam-1)*((Uh_L(2,i,j))^3+Uh_L(2,i,j)*(Uh_L(3,i,j))^2)/(Uh_L(1,i,j)^2 )
];
Fh_R(:,i,j)=[ Uh_R(2,i,j);
(gam-1)*Uh_R(4,i,j)+0.5*(3-gam)*Uh_R(2,i,j)^2/Uh_R(1,i,j)-0.5*(gam-1)*Uh_R(3,i,j)^2/Uh_R(1,i,j);
Uh_R(2,i,j)*Uh_R(3,i,j)/Uh_R(1,i,j);
gam*Uh_R(2,i,j)*Uh_R(4,i,j)/Uh_R(1,i,j)-0.5*(gam-1)*((Uh_R(2,i,j))^3+Uh_R(2,i,j)*(Uh_R(3,i,j))^2)/(Uh_R(1,i,j)^2 )
];
D=sqrt(abs(Uh_R(1,i,j)/Uh_L(1,i,j)));
u=(Uh_L(2,i,j)/Uh_L(1,i,j)+D*Uh_R(2,i,j)/Uh_R(1,i,j))/(1+D);
v=(Uh_L(3,i,j)/Uh_L(1,i,j)+D*Uh_R(3,i,j)/Uh_R(1,i,j))/(1+D);
E=(Uh_L(4,i,j)/Uh_L(1,i,j)+D*Uh_R(4,i,j)/Uh_R(1,i,j))/(1+D);
a=((gam-1)*gam*(E-0.5*(u^2+v^2)))^0.5;
kx=1;
ky=0;
V=u^2+v^2;
qk=kx*u+ky*v;
TkF=[ -(gam-1)/(a^2) 0 -1/2/a 1/2/a
-(gam-1)*u/(a^2) -ky kx/2-u/2/a kx/2+u/2/a
-(gam-1)*v/(a^2) kx ky/2-v/2/a ky/2+v/2/a
-(gam-1)*V/2/(a^2) kx*v-ky*u (a*(kx*u+ky*v)-V/2-(a^2)/(gam-1))/2/a (a*(kx*u+ky*v)+V/2+(a^2)/(gam-1))/2/a
];
TkF_1=[V/2-(a^2)/(gam-1) -u -v 1
-kx*v+ky*u -ky kx 0
-kx*u-ky*v-(gam-1)*V/2/a kx+(gam-1)*u/a ky+(gam-1)*v/a -(gam-1)/a
-kx*u-ky*v+(gam-1)*V/2/a kx-(gam-1)*u/a ky-(gam-1)*v/a (gam-1)/a
];
lamta1=getlamta(qk,0.15*(abs(qk)+a));
lamta2=getlamta(qk,0.15*(abs(qk)+a));
lamta3=getlamta(qk-a,0.15*(abs(qk)+a));
lamta4=getlamta(qk+a,0.15*(abs(qk)+a));
Diag=[lamta1 0 0 0;
0 lamta2 0 0;
0 0 lamta3 0
0 0 0 lamta4];
TF=TkF*abs(Diag)*TkF_1;
FhX(:,i,j)=0.5*(Fh_R(:,i,j)+Fh_L(:,i,j))-0.5*TF*(Uh_R(:,i,j)-Uh_L(:,i,j));
end
end
%%%
for i=1:N
for j=1:M-1
Gh_L(:,i,j)=[ Uh_L(3,i,j);
Uh_L(2,i,j)*Uh_L(3,i,j)/Uh_L(1,i,j);
(gam-1)*Uh_L(4,i,j)+0.5*(3-gam)*(Uh_L(3,i,j))^2/Uh_L(1,i,j)-0.5*(gam-1)*(Uh_L(2,i,j))^2/Uh_L(1,i,j);
gam*Uh_L(3,i,j)*Uh_L(4,i,j)/Uh_L(1,i,j)-0.5*(gam-1)*((Uh_L(3,i,j))^3+Uh_L(3,i,j)*(Uh_L(2,i,j))^2)/(Uh_L(1,i,j)^2 )
];
Gh_R(:,i,j)=[ Uh_R(3,i,j);
Uh_R(2,i,j)*Uh_R(3,i,j)/Uh_R(1,i,j);
(gam-1)*Uh_R(4,i,j)+0.5*(3-gam)*(Uh_R(3,i,j))^2/Uh_R(1,i,j)-0.5*(gam-1)*(Uh_R(2,i,j))^2/Uh_R(1,i,j);
gam*Uh_R(3,i,j)*Uh_R(4,i,j)/Uh_R(1,i,j)-0.5*(gam-1)*((Uh_R(3,i,j))^3+Uh_R(3,i,j)*(Uh_R(2,i,j))^2)/(Uh_R(1,i,j)^2 )
];
D=sqrt(abs(Uh_R(1,i,j)/Uh_L(1,i,j)));
u=(Uh_L(2,i,j)/Uh_L(1,i,j)+D*Uh_R(2,i,j)/Uh_R(1,i,j))/(1+D);
v=(Uh_L(3,i,j)/Uh_L(1,i,j)+D*Uh_R(3,i,j)/Uh_R(1,i,j))/(1+D);
E=(Uh_L(4,i,j)/Uh_L(1,i,j)+D*Uh_R(4,i,j)/Uh_R(1,i,j))/(1+D);
a=((gam-1)*gam*(E-0.5*(u^2+v^2)))^0.5;
kx=0;
ky=1;
V=u^2+v^2;
qk=kx*u+ky*v;
TkG=[ -(gam-1)/(a^2) 0 -1/2/a 1/2/a
-(gam-1)*u/(a^2) -ky kx/2-u/2/a kx/2+u/2/a
-(gam-1)*v/(a^2) kx ky/2-v/2/a ky/2+v/2/a
-(gam-1)*V/2/(a^2) kx*v-ky*u (a*(kx*u+ky*v)-V/2-(a^2)/(gam-1))/2/a (a*(kx*u+ky*v)+V/2+(a^2)/(gam-1))/2/a
];
TkG_1=[V/2-(a^2)/(gam-1) -u -v 1
-kx*v+ky*u -ky kx 0
-kx*u-ky*v-(gam-1)*V/2/a kx+(gam-1)*u/a ky+(gam-1)*v/a -(gam-1)/a
-kx*u-ky*v+(gam-1)*V/2/a kx-(gam-1)*u/a ky-(gam-1)*v/a (gam-1)/a
];
lamta1=getlamta(qk,0.15*(abs(qk)+a));
lamta2=getlamta(qk,0.15*(abs(qk)+a));
lamta3=getlamta(qk-a,0.15*(abs(qk)+a));
lamta4=getlamta(qk+a,0.15*(abs(qk)+a));
Diag=[lamta1 0 0 0;
0 lamta2 0 0;
0 0 lamta3 0
0 0 0 lamta4];
TG=TkG*abs(Diag)*TkG_1;
GhX(:,i,j)=0.5*(Gh_R(:,i,j)+Gh_L(:,i,j))-0.5*TG*(Uh_R(:,i,j)-Uh_L(:,i,j));
end
end
%%%
for i=1:N-1
for j=1:M
Fh_B(:,i,j)=[ Uh_B(2,i,j);
(gam-1)*Uh_B(4,i,j)+0.5*(3-gam)*(Uh_B(2,i,j))^2/Uh_B(1,i,j)-0.5*(gam-1)*(Uh_B(3,i,j))^2/Uh_B(1,i,j);
Uh_B(2,i,j)*Uh_B(3,i,j)/Uh_B(1,i,j);
gam*Uh_B(2,i,j)*Uh_B(4,i,j)/Uh_B(1,i,j)-0.5*(gam-1)*((Uh_B(2,i,j))^3+Uh_B(2,i,j)*(Uh_B(3,i,j))^2)/(Uh_B(1,i,j)^2 )
];
Fh_T(:,i,j)=[ Uh_T(2,i,j);
(gam-1)*Uh_T(4,i,j)+0.5*(3-gam)*Uh_T(2,i,j)^2/Uh_T(1,i,j)-0.5*(gam-1)*Uh_T(3,i,j)^2/Uh_T(1,i,j);
Uh_T(2,i,j)*Uh_T(3,i,j)/Uh_T(1,i,j);
gam*Uh_T(2,i,j)*Uh_T(4,i,j)/Uh_T(1,i,j)-0.5*(gam-1)*((Uh_T(2,i,j))^3+Uh_T(2,i,j)*(Uh_T(3,i,j))^2)/(Uh_T(1,i,j)^2 )
];
D=sqrt(abs(Uh_T(1,i,j)/Uh_B(1,i,j)));
u=(Uh_B(2,i,j)/Uh_B(1,i,j)+D*Uh_T(2,i,j)/Uh_T(1,i,j))/(1+D);
v=(Uh_B(3,i,j)/Uh_B(1,i,j)+D*Uh_T(3,i,j)/Uh_T(1,i,j))/(1+D);
E=(Uh_B(4,i,j)/Uh_B(1,i,j)+D*Uh_T(4,i,j)/Uh_T(1,i,j))/(1+D);
a=((gam-1)*gam*(E-0.5*(u^2+v^2)))^0.5;
kx=1;
ky=0;
V=u^2+v^2;
qk=kx*u+ky*v;
TkF=[ -(gam-1)/(a^2) 0 -1/2/a 1/2/a
-(gam-1)*u/(a^2) -ky kx/2-u/2/a kx/2+u/2/a
-(gam-1)*v/(a^2) kx ky/2-v/2/a ky/2+v/2/a
-(gam-1)*V/2/(a^2) kx*v-ky*u (a*(kx*u+ky*v)-V/2-(a^2)/(gam-1))/2/a (a*(kx*u+ky*v)+V/2+(a^2)/(gam-1))/2/a
];
TkF_1=[V/2-(a^2)/(gam-1) -u -v 1
-kx*v+ky*u -ky kx 0
-kx*u-ky*v-(gam-1)*V/2/a kx+(gam-1)*u/a ky+(gam-1)*v/a -(gam-1)/a
-kx*u-ky*v+(gam-1)*V/2/a kx-(gam-1)*u/a ky-(gam-1)*v/a (gam-1)/a
];
lamta1=getlamta(qk,0.15*(abs(qk)+a));
lamta2=getlamta(qk,0.15*(abs(qk)+a));
lamta3=getlamta(qk-a,0.15*(abs(qk)+a));
lamta4=getlamta(qk+a,0.15*(abs(qk)+a));
Diag=[lamta1 0 0 0;
0 lamta2 0 0;
0 0 lamta3 0
0 0 0 lamta4];
TF=TkF*abs(Diag)*TkF_1;
FhY(:,i,j)=0.5*(Fh_T(:,i,j)+Fh_B(:,i,j))-0.5*TF*(Uh_T(:,i,j)-Uh_B(:,i,j));
end
end
%%%
for i=1:N-1
for j=1:M
Gh_B(:,i,j)=[ Uh_B(3,i,j);
Uh_B(2,i,j)*Uh_B(3,i,j)/Uh_B(1,i,j);
(gam-1)*Uh_B(4,i,j)+0.5*(3-gam)*(Uh_B(3,i,j))^2/Uh_B(1,i,j)-0.5*(gam-1)*(Uh_B(2,i,j))^2/Uh_B(1,i,j);
gam*Uh_B(3,i,j)*Uh_B(4,i,j)/Uh_B(1,i,j)-0.5*(gam-1)*((Uh_B(3,i,j))^3+Uh_B(3,i,j)*(Uh_B(2,i,j))^2)/(Uh_B(1,i,j)^2 )
];
Gh_T(:,i,j)=[ Uh_T(3,i,j);
Uh_T(2,i,j)*Uh_T(3,i,j)/Uh_T(1,i,j);
(gam-1)*Uh_T(4,i,j)+0.5*(3-gam)*(Uh_T(3,i,j))^2/Uh_T(1,i,j)-0.5*(gam-1)*(Uh_T(2,i,j))^2/Uh_T(1,i,j);
gam*Uh_T(3,i,j)*Uh_T(4,i,j)/Uh_T(1,i,j)-0.5*(gam-1)*((Uh_T(3,i,j))^3+Uh_T(3,i,j)*(Uh_T(2,i,j))^2)/(Uh_T(1,i,j)^2 )
];
D=sqrt(abs(Uh_T(1,i,j)/Uh_B(1,i,j)));
u=(Uh_B(2,i,j)/Uh_B(1,i,j)+D*Uh_T(2,i,j)/Uh_T(1,i,j))/(1+D);
v=(Uh_B(3,i,j)/Uh_B(1,i,j)+D*Uh_T(3,i,j)/Uh_T(1,i,j))/(1+D);
E=(Uh_B(4,i,j)/Uh_B(1,i,j)+D*Uh_T(4,i,j)/Uh_T(1,i,j))/(1+D);
a=((gam-1)*gam*(E-0.5*(u^2+v^2)))^0.5;
kx=0;
ky=1;
V=u^2+v^2;
qk=kx*u+ky*v;
TkG=[ -(gam-1)/(a^2) 0 -1/2/a 1/2/a
-(gam-1)*u/(a^2) -ky kx/2-u/2/a kx/2+u/2/a
-(gam-1)*v/(a^2) kx ky/2-v/2/a ky/2+v/2/a
-(gam-1)*V/2/(a^2) kx*v-ky*u (a*(kx*u+ky*v)-V/2-(a^2)/(gam-1))/2/a (a*(kx*u+ky*v)+V/2+(a^2)/(gam-1))/2/a
];
TkG_1=[V/2-(a^2)/(gam-1) -u -v 1
-kx*v+ky*u -ky kx 0
-kx*u-ky*v-(gam-1)*V/2/a kx+(gam-1)*u/a ky+(gam-1)*v/a -(gam-1)/a
-kx*u-ky*v+(gam-1)*V/2/a kx-(gam-1)*u/a ky-(gam-1)*v/a (gam-1)/a
];
lamta1=getlamta(qk,0.15*(abs(qk)+a));
lamta2=getlamta(qk,0.15*(abs(qk)+a));
lamta3=getlamta(qk-a,0.15*(abs(qk)+a));
lamta4=getlamta(qk+a,0.15*(abs(qk)+a));
Diag=[lamta1 0 0 0;
0 lamta2 0 0;
0 0 lamta3 0
0 0 0 lamta4];
TG=TkG*abs(Diag)*TkG_1;
GhY(:,i,j)=0.5*(Gh_T(:,i,j)+Gh_B(:,i,j))-0.5*TG*(Uh_T(:,i,j)-Uh_B(:,i,j));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function lamta=getlamta(temp,delta)
if abs(temp)>=delta
lamta=abs(temp);
else
lamta=(temp^2+delta^2)/2/delta;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function run=stopcalculate(t,sumt,step,sumstep,constant,error)
run=error>constant;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -