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

二维超声速流动的数值解——普朗特迈耶稀疏波 主程序数值计算 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(xlt;=10)); ys=-(x-10)*tan(sita); ys(1:point1)=0; h=H-ys; h(1:point1)=H; y(j,i)=0; for k=1:i y(:,k)=linspace(ys(k),H,j)#39;; 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); %麦考马克方法空间推进计算 for k=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.