📄 hw-nozzle.m
字号:
U(4,jindex,iindex)=P_in/(gam-1)+rho_in*Vt^2/2;
end
end
U_in = U(:,1,1);
U_out = U(:,1,M);
for jindex=1:N
for iindex=3:M-2
U(:,jindex,iindex) = U_in*((M - 2 - iindex)/(M-4)) + U_out*((iindex - 2)/(M-4));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [U,para]=getpara(U,Data)
gam=Data.gam;
cfl=Data.cfl;
X=Data.CyPX;
Y=Data.CyPY;
M=Data.Nx+1;
N=Data.Ny+1;
for iindex=1:N % N 为行数，M为列数 第2列为进口 M－1 列为出口
for jindex=1:M
rho0(iindex,jindex)=U(1,iindex,jindex); %初始时刻的流场参数值
u0(iindex,jindex)=U(2,iindex,jindex)/rho0(iindex,jindex);
v0(iindex,jindex)=U(3,iindex,jindex)/rho0(iindex,jindex);
E0(iindex,jindex)=U(4,iindex,jindex)/rho0(iindex,jindex);
p0(iindex,jindex)=(gam-1)*rho0(iindex,jindex)*(E0(iindex,jindex)-0.5*((u0(iindex,jindex))^2+(v0(iindex,jindex))^2));
a0(iindex,jindex)=sqrt(gam*p0(iindex,jindex)/rho0(iindex,jindex));
M0(iindex,jindex)=sqrt((u0(iindex,jindex))^2+(v0(iindex,jindex))^2)/a0(iindex,jindex);
end
end
para=struct('gam',gam,'cfl',cfl,'rho',rho0,'p',p0,'u',u0,'v',v0,'a',a0,'N',N,'M',M,'LRPX',Data.CyPX,'LRPY',Data.CyPY,'BTPX',Data.CxPX,'BTPY',Data.CxPY,'X',Data.UPX,'Y',Data.UPY,'PX',Data.PX,'PY',Data.PY,'E0',E0,'M0',M0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [U,dt]=Update(U,para)
k=5;
RK_ceof(1)= 0.1067; 龙格－裤它 格式 各项时间多部系数 5 步
RK_ceof(2)= 0.1979;
RK_ceof(3)= 0.3232;
RK_ceof(4)= 0.5201;
RK_ceof(5)= 1.0000;
while k
[Uh_L,Uh_R,Uh_B,Uh_T]=GetU_L_R_B_T_Vanleer(U,para);
dt=tstep(Uh_L,Uh_R,Uh_B,Uh_T,para);
[FhX,GhX,FhY,GhY]=GetFh_MUSCL(U,Uh_L,Uh_R,Uh_B,Uh_T,para);
U=Runge_Kutta(RK_ceof(6-k),U,FhX,GhX,FhY,GhY,dt,para);
k=k-1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function U=Runge_Kutta(afa,U,FhX,GhX,FhY,GhY,dt,para)
gam=para.gam;
M=para.M;
N=para.N;
X=para.X;
Y=para.Y;
PX=para.PX;
PY=para.PY;
for i=2:N-1
for j=2:M-1
dx12=PX(i-1,j-1)-PX(i,j-1);
dy12=PY(i-1,j-1)-PY(i,j-1);
dx34=PX(i,j)-PX(i-1,j);
dy34=PY(i,j)-PY(i-1,j);
dx23=PX(i-1,j)-PX(i-1,j-1);
dy23=PY(i-1,j)-PY(i-1,j-1);
dx41=PX(i,j-1)-PX(i,j);
dy41=PY(i,j-1)-PY(i,j);
S=0.5*abs((PX(i,j)-PX(i-1,j-1))*(PY(i-1,j)-PY(i,j-1))-(PY(i,j)-PY(i-1,j-1))*(PX(i-1,j)-PX(i,j-1)));
U(:,i,j)=U(:,i,j)-(dt/S)*afa*( (FhX(:,i,j-1)*dy12 - GhX(:,i,j-1)*dx12) + (FhX(:,i,j)*dy34 - GhX(:,i,j)*dx34) + (FhY(:,i-1,j)*dy23 - GhY(:,i-1,j)*dx23) + (FhY(:,i,j)*dy41 - GhY(:,i,j)*dx41) );
end
end
%在这里进行边界条件的设置
%进口,出口边界
P_in=130000;%入口静压Pa
P_out=101325;%出口反压
Vt=200.0;%入口平均速度m
afa=0;%入口角
for iindex=2:N-1
for jindex=2:M-1
rho0(iindex,jindex)=U(1,iindex,jindex);
u0(iindex,jindex)=U(2,iindex,jindex)/rho0(iindex,jindex);
v0(iindex,jindex)=U(3,iindex,jindex)/rho0(iindex,jindex);
E0(iindex,jindex)=U(4,iindex,jindex)/rho0(iindex,jindex);
p0(iindex,jindex)=(gam-1)*rho0(iindex,jindex)*(E0(iindex,jindex)-0.5*((u0(iindex,jindex))^2+(v0(iindex,jindex))^2));
a0(iindex,jindex)=sqrt(gam*p0(iindex,jindex)/rho0(iindex,jindex));
M0(iindex,jindex)=sqrt((u0(iindex,jindex))^2+(v0(iindex,jindex))^2)/a0(iindex,jindex);
end
%进口处理
rho0(iindex,1)=1.579;
u0(iindex,1)=Vt*cos(afa);
v0(iindex,1)=Vt*sin(afa);
E0(iindex,1)=2*E0(iindex,2)-E0(iindex,3);
p0(iindex,1)=(gam-1)*rho0(iindex,1)*(E0(iindex,1)-0.5*((u0(iindex,1))^2+(v0(iindex,1))^2));
a0(iindex,1)=sqrt(gam*p0(iindex,1)/rho0(iindex,1));
M0(iindex,1)=sqrt((u0(iindex,1))^2+(v0(iindex,1))^2)/a0(iindex,1);
U(1,iindex,1)=rho0(iindex,1);
U(2,iindex,1)=rho0(iindex,1)*u0(iindex,1);
U(3,iindex,1)=rho0(iindex,1)*v0(iindex,1);
U(4,iindex,1)=p0(iindex,1)/(gam-1)+rho0(iindex,1)*((u0(iindex,1))^2+(v0(iindex,1))^2)/2;
rho0(iindex,2)=(rho0(iindex,1)+rho0(iindex,3))/2;
u0(iindex,2)=(u0(iindex,1)+u0(iindex,3))/2;
v0(iindex,2)=(v0(iindex,1)+v0(iindex,3))/2;
U(1,iindex,2)=rho0(iindex,2);
U(2,iindex,2)=rho0(iindex,2)*u0(iindex,2);
U(3,iindex,2)=rho0(iindex,2)*v0(iindex,2);
U(4,iindex,2)=p0(iindex,2)/(gam-1)+rho0(iindex,2)*((u0(iindex,2))^2+(v0(iindex,2))^2)/2;
%出口处理
u0(iindex,M)=2*u0(iindex,M-1)-u0(iindex,M-2);
v0(iindex,M)=2*v0(iindex,M-1)-v0(iindex,M-2);
rho0(iindex,M)=2*rho0(iindex,M-1)-rho0(iindex,M-2);
p0(iindex,M)=P_out;
E0(iindex,M)=p0(iindex,M)/(gam-1)/rho0(iindex,M)+(u0(iindex,M))^2+(v0(iindex,M))^2;
a0(iindex,M)=sqrt(gam*p0(iindex,M)/rho0(iindex,M));
M0(iindex,M)=sqrt((u0(iindex,M))^2+(v0(iindex,M))^2)/a0(iindex,M);
U(1,iindex,M)=rho0(iindex,M);
U(2,iindex,M)=rho0(iindex,M)*u0(iindex,M);
U(3,iindex,M)=rho0(iindex,M)*v0(iindex,M);
U(4,iindex,M)=p0(iindex,M)/(gam-1)+rho0(iindex,M)*((u0(iindex,M))^2+(v0(iindex,M))^2)/2;
p0(iindex,M-1)=(p0(iindex,M)+p0(iindex,M-2))/2;
U(1,iindex,M-1)=rho0(iindex,M-1);
U(2,iindex,M-1)=rho0(iindex,M-1)*u0(iindex,M-1);
U(3,iindex,M-1)=rho0(iindex,M-1)*v0(iindex,M-1);
U(4,iindex,M-1)=p0(iindex,M-1)/(gam-1)+rho0(iindex,M-1)*((u0(iindex,M-1))^2+(v0(iindex,M-1))^2)/2;
end
%壁面边界
for iindex=2:M-1
temp_rho0=U(1,2,iindex);
temp_u0=U(2,2,iindex)/temp_rho0;
temp_v0=U(3,2,iindex)/temp_rho0;
temp_E0=U(4,2,iindex)/temp_rho0;
temp_p0=(gam-1)*temp_rho0*(temp_E0-0.5*(temp_u0^2+temp_v0^2));
%法向速度置为0
V = sqrt(temp_u0^2+temp_v0^2);
tg_a1=temp_v0/temp_u0;
a1=atan(tg_a1);
tg_a2=(Y(N-1,iindex)-Y(N-1,iindex-1))/(X(N-1,iindex)-X(N-1,iindex-1));
a2=atan(tg_a2);
Vt=V*cos(a1-a2);
temp_u0=Vt*cos(a2);
temp_v0=-Vt*sin(a2);
temp_rho02=U(1,3,iindex);
temp_u02=U(2,3,iindex)/temp_rho02;
temp_v02=U(3,3,iindex)/temp_rho02;
temp_E02=U(4,3,iindex)/temp_rho02;
temp_p02=(gam-1)*temp_rho02*(temp_E02-0.5*(temp_u02^2+temp_v02^2));
temp_p0=temp_p02;
temp_E0=temp_E02;
%再根据新值来确定U
U(1,2,iindex) = temp_rho0;
U(2,2,iindex) = temp_rho0*temp_u0;
U(3,2,iindex) = temp_rho0*temp_v0;
U(4,2,iindex) = temp_rho0*temp_E0;
temp_rho0=U(1,N-1,iindex);
temp_u0=U(2,N-1,iindex)/temp_rho0;
temp_v0=U(3,N-1,iindex)/temp_rho0;
temp_E0=U(4,N-1,iindex)/temp_rho0;
temp_p0=(gam-1)*temp_rho0*(temp_E0-0.5*(temp_u0^2+temp_v0^2));
%法向速度置为0
V = sqrt(temp_u0^2+temp_v0^2);
tg_a1=temp_v0/temp_u0;
a1=atan(tg_a1);
tg_a2=(Y(N-1,iindex)-Y(N-1,iindex-1))/(X(N-1,iindex)-X(N-1,iindex-1));
a2=atan(tg_a2);
Vt=V*cos(a1-a2);
temp_u0=Vt*cos(a2);
temp_v0=Vt*sin(a2);
temp_rho02=U(1,N-2,iindex);
temp_u02=U(2,N-2,iindex)/temp_rho02;
temp_v02=U(3,N-2,iindex)/temp_rho02;
temp_E02=U(4,N-2,iindex)/temp_rho02;
temp_p02=(gam-1)*temp_rho02*(temp_E02-0.5*(temp_u02^2+temp_v02^2));
temp_p0=temp_p02;
temp_E0=temp_E02;
%再根据新值来确定U
U(1,N-1,iindex) = temp_rho0;
U(2,N-1,iindex) = temp_rho0*temp_u0;
U(3,N-1,iindex) = temp_rho0*temp_v0;
U(4,N-1,iindex) = temp_rho0*temp_E0;
U(:,1,iindex) = 2*U(:,2,iindex) - U(:,3,iindex);
U(:,N,iindex) = 2*U(:,N-1,iindex) - U(:,N-2,iindex);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Uh_L,Uh_R,Uh_B,Uh_T]=GetU_L_R_B_T_Vanleer(U,para); % 构造 MUSCL 格式
gam=para.gam;
M=para.M;
N=para.N;
X=para.X;
Y=para.Y;
temp_e=0.00000001;
temp_g=1/3;
%单方向分步求解
for i=1:N
for j=2:M-1
for k=1:4
deltaUp_X(k)=U(k,i,j+1)-U(k,i,j);
deltaDown_X(k)=U(k,i,j)-U(k,i,j-1);
q_X(k)=2* deltaUp_X(k) * deltaDown_X(k)/(temp_e + deltaUp_X(k)^2 + deltaDown_X(k)^2);
S_L(k)=0.5*q_X(k)*((1-q_X(k)*temp_g)*deltaUp_X(k)+(1+q_X(k)*temp_g)*deltaDown_X(k));
S_R(k)=0.5*q_X(k)*((1+q_X(k)*temp_g)*deltaUp_X(k)+(1-q_X(k)*temp_g)*deltaDown_X(k));
U_t_L(k,i,j)=U(k,i,j)+S_L(k)/2;
U_t_R(k,i,j)=U(k,i,j)-S_R(k)/2;
end
end
end
for i=1:N
for j=1:M-1
if(j==1)
Uh_L(:,i,j)=U(:,i,j);
Uh_R(:,i,j)=U_t_L(:,i,j+1);
elseif(j==M-1)
Uh_L(:,i,j)=U_t_R(:,i,j);
Uh_R(:,i,j)=U(:,i,j+1);
else
Uh_L(:,i,j)=U_t_R(:,i,j);
Uh_R(:,i,j)=U_t_L(:,i,j+1);
end
end
end
for j=1:M
for i=2:N-1
for k=1:4
deltaUp_Y(k)=U(k,i+1,j)-U(k,i,j);
deltaDown_Y(k)=U(k,i,j)-U(k,i-1,j);
q_Y(k)=2* deltaUp_Y(k) * deltaDown_Y(k)/(temp_e + deltaUp_Y(k)^2 + deltaDown_Y(k)^2);
S_B(k)=0.5*q_Y(k)*((1-q_Y(k)*temp_g)*deltaUp_Y(k)+(1+q_Y(k)*temp_g)*deltaDown_Y(k));
S_T(k)=0.5*q_Y(k)*((1+q_Y(k)*temp_g)*deltaUp_Y(k)+(1-q_Y(k)*temp_g)*deltaDown_Y(k));
U_t_B(k,i,j)=U(k,i,j)+S_B(k)/2;
U_t_T(k,i,j)=U(k,i,j)-S_T(k)/2;
end
end
end
for j=1:M
for i=1:N-1
if(i==1)
Uh_B(:,i,j)=U(:,i,j);
Uh_T(:,i,j)=U_t_B(:,i+1,j);
elseif(i==N-1)
Uh_B(:,i,j)=U_t_T(:,i,j);
Uh_T(:,i,j)=U(:,i+1,j);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -