资源描述：
《二维超声速流动的数值解——普朗特迈耶稀疏波.doc》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、二维超声速流动的数值解——普朗特迈耶稀疏波主程序数值计算%二维超声速流动的数值解——普朗特迈耶稀疏波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
2、=H-ys;h(1:point1)=H;y(j,i)=0;fork=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
3、,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_ep
4、s=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
5、=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_
6、av=zeros(j-1,1);F3_eps_av=zeros(j-1,1);F4_eps_av=zeros(j-1,1);%麦考马克方法空间推进计算fork=1:iG1=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
7、).*(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
8、-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(sit
