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

二维超声速流动的数值解——普朗特迈耶稀疏波主程序数值计算%二维超声速流动的数值解——普朗特迈耶稀疏波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;for k=1:iy(:,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);%麦考马克方法空间推进计算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.*((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; %计算空间空间步长计算mi