1、二维超声速流动的数值解普朗特迈耶稀疏波主程序数值计算%二维超声速流动的数值解普朗特迈耶稀疏波clc;clear;%常量设置CFL=0.5;Cy=0.2;r=1.4;R=287;%初始物理条件Ma0=2;P0=1.01*105;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
2、:i y(:,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;
3、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=zer
4、os(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(
5、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
6、-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)-F
7、3(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; %计算空间空间步长计算 miu=asin(1./Ma(1:j-1,k); deps1=abs(tan(sita+miu); deps2=abs(tan(sita-miu); deps=CFL*h(k)/(j-1)/max(deps1;deps2); %人工粘性预估值 SF1(2:j-1)=Cy*abs(P(3:j,k)-2*P(2:j-1
8、,k)+P(1:j-2,k)./(P(3:j,k)+2*P(2:j-1,k)+P(1:j-2,k).*(F1(3:j)-2*F1(2:j-1)+F1(1:j-2); SF2(2:j-1)=Cy*abs(P(3:j,k)-2*P(2:j-1,k)+P(1:j-2,k)./(P(3:j,k)+2*P(2:j-1,k)+P(1:j-2,k).*(F2(3:j)-2*F2(2:j-1)+F2(1:j-2); SF3(2:j-1)=Cy*abs(P(3:j,k)-2*P(2:j-1,k)+P(1:j-2,k)./(P(3:j,k)+2*P(2:j-1,k)+P(1:j-2,k).*(F3(3:j)-2*
9、F3(2:j-1)+F3(1:j-2); SF4(2:j-1)=Cy*abs(P(3:j,k)-2*P(2:j-1,k)+P(1:j-2,k)./(P(3:j,k)+2*P(2:j-1,k)+P(1:j-2,k).*(F4(3:j)-2*F4(2:j-1)+F4(1:j-2); SF1(1)=2*SF1(2)-SF1(3); SF2(1)=2*SF2(2)-SF2(3); SF3(1)=2*SF3(2)-SF3(3); SF4(1)=2*SF4(2)-SF4(3); %预估值 F1_1(1:j-1)=F1(1:j-1)+F1_eps(1:j-1)*deps+SF1(1:j-1); F2_1(1
10、:j-1)=F2(1:j-1)+F2_eps(1:j-1)*deps+SF2(1:j-1); F3_1(1:j-1)=F3(1:j-1)+F3_eps(1:j-1)*deps+SF3(1:j-1); F4_1(1:j-1)=F4(1:j-1)+F4_eps(1:j-1)*deps+SF4(1:j-1); A(1:j-1)=F3_1(1:j-1).2./F1_1(1:j-1)/2-F4_1(1:j-1); B(1:j-1)=r/(r-1)*F1_1(1:j-1).*F2_1(1:j-1); C(1:j-1)=-(r+1)/(r-1)/2*F1_1(1:j-1).3; M_1(1:j-1)=(-B
11、(1:j-1)+sqrt(B(1:j-1).2-4*A(1:j-1).*C(1:j-1)./A(1:j-1)/2; P_1(1:j-1)=F2_1(1:j-1)-F1_1(1:j-1).2./M_1(1:j-1); G1_1(1:j-1)=M_1(1:j-1).*F3_1(1:j-1)./F1_1(1:j-1); G2_1(1:j-1)=F3_1(1:j-1); G3_1(1:j-1)=M_1(1:j-1).*(F3_1(1:j-1)./F1_1(1:j-1).2+F2_1(1:j-1)-F1_1(1:j-1).2./M_1(1:j-1); G4_1(1:j-1)=r/(r-1)*(F2_1(
12、1:j-1)-F1_1(1:j-1).2./M_1(1:j-1).*F3_1(1:j-1)./F1_1(1:j-1)+. M_1(1:j-1)/2.*F3_1(1:j-1)./F1_1(1:j-1).*(F1_1(1:j-1)./M_1(1:j-1).2+(F3_1(1:j-1)./F1_1(1:j-1).2); %修正偏导数计算 F1_eps_1(2:j-1)=eta_x(2:j-1,k).*(F1_1(1:j-2)-F1_1(2:j-1)/deta+1/h(k)*(G1_1(1:j-2)-G1_1(2:j-1)/deta; F2_eps_1(2:j-1)=eta_x(2:j-1,k).*(
13、F2_1(1:j-2)-F2_1(2:j-1)/deta+1/h(k)*(G2_1(1:j-2)-G2_1(2:j-1)/deta; F3_eps_1(2:j-1)=eta_x(2:j-1,k).*(F3_1(1:j-2)-F3_1(2:j-1)/deta+1/h(k)*(G3_1(1:j-2)-G3_1(2:j-1)/deta; F4_eps_1(2:j-1)=eta_x(2:j-1,k).*(F4_1(1:j-2)-F4_1(2:j-1)/deta+1/h(k)*(G4_1(1:j-2)-G4_1(2:j-1)/deta; F1_eps_av(2:j-1)=0.5*(F1_eps_1(2:
14、j-1)+F1_eps(2:j-1); F2_eps_av(2:j-1)=0.5*(F2_eps_1(2:j-1)+F2_eps(2:j-1); F3_eps_av(2:j-1)=0.5*(F3_eps_1(2:j-1)+F3_eps(2:j-1); F4_eps_av(2:j-1)=0.5*(F4_eps_1(2:j-1)+F4_eps(2:j-1); %人工粘性修正值 SF1_1(2:j-2)=Cy*abs(P_1(3:j-1)-2*P_1(2:j-2)+P_1(1:j-3)./(P_1(3:j-1)+2*P_1(2:j-2)+P_1(1:j-3).*(F1_1(3:j-1)-2*F1_1
15、(2:j-2)+F1_1(1:j-3); SF2_1(2:j-2)=Cy*abs(P_1(3:j-1)-2*P_1(2:j-2)+P_1(1:j-3)./(P_1(3:j-1)+2*P_1(2:j-2)+P_1(1:j-3).*(F2_1(3:j-1)-2*F2_1(2:j-2)+F2_1(1:j-3); SF3_1(2:j-2)=Cy*abs(P_1(3:j-1)-2*P_1(2:j-2)+P_1(1:j-3)./(P_1(3:j-1)+2*P_1(2:j-2)+P_1(1:j-3).*(F3_1(3:j-1)-2*F3_1(2:j-2)+F3_1(1:j-3); SF4_1(2:j-2)=
16、Cy*abs(P_1(3:j-1)-2*P_1(2:j-2)+P_1(1:j-3)./(P_1(3:j-1)+2*P_1(2:j-2)+P_1(1:j-3).*(F4_1(3:j-1)-2*F4_1(2:j-2)+F4_1(1:j-3); SF1_1(j-1)=2*SF1_1(j-2)-SF1_1(j-3); SF2_1(j-1)=2*SF2_1(j-2)-SF2_1(j-3); SF3_1(j-1)=2*SF3_1(j-2)-SF3_1(j-3); SF4_1(j-1)=2*SF4_1(j-2)-SF4_1(j-3); %修正值 F1(2:j-1)=F1(2:j-1)+F1_eps_av(2
17、:j-1)*deps+SF1_1(2:j-1); F2(2:j-1)=F2(2:j-1)+F2_eps_av(2:j-1)*deps+SF2_1(2:j-1); F3(2:j-1)=F3(2:j-1)+F3_eps_av(2:j-1)*deps+SF3_1(2:j-1); F4(2:j-1)=F4(2:j-1)+F4_eps_av(2:j-1)*deps+SF4_1(2:j-1); %上边界单侧差分 F1(j)=2*F1(j-1)-F1(j-2); F2(j)=2*F2(j-1)-F2(j-2); F3(j)=2*F3(j-1)-F3(j-2); F4(j)=2*F4(j-1)-F4(j-2)
18、; %下边界单侧差分 F1(1)=2*F1(2)-F1(3); F2(1)=2*F2(2)-F2(3); F3(1)=2*F3(2)-F3(3); F4(1)=2*F4(2)-F4(3); %所有网格点处物理量计算 A=F3.2./F1/2-F4; B=r/(r-1)*F1.*F2; C=-(r+1)/(r-1)/2*F1.3; M(:,k+1)=(-B+sqrt(B.2-4*A.*C)./A/2; U(:,k+1)=F1./M(:,k+1); V(:,k+1)=F3./F1; P(:,k+1)=F2-F1.*U(:,k+1); T(:,k+1)=P(:,k+1)./M(:,k+1)/R; M
19、a(:,k+1)=sqrt(U(:,k+1).2+V(:,k+1).2)./sqrt(r*R*T(:,k+1); %下边界处物理量修正 Ma_cal=sqrt(U(1,k+1)2+V(1,k+1)2)/sqrt(r*R*T(1,k+1); f_cal=sqrt(r+1)/(r-1)*atan(sqrt(r-1)*(Ma_cal2-1)/(r+1)-atan(sqrt(Ma_cal2-1); if k=0.00001&n100 fm=sqrt(r+1)/(r-1)*atan(sqrt(r-1)/(r+1)*(xm.2-1)-atan(sqrt(xm.2-1)-f_act; if fm(2)*fm
20、(3)0 xm=xm(1) (xm(1)+xm(2)/2 xm(2); else xm=xm(2) xm(2) xm(2); end n=n+1; end Ma(1,k+1)=xm(2); %下边界处物理量修正值 P(1,k+1)=P(1,k+1)*(1+(r-1)/2*Ma_cal2)/(1+(r-1)/2*Ma(1,k+1)2)(r/(r-1); T(1,k+1)=T(1,k+1)*(1+(r-1)/2*Ma_cal2)/(1+(r-1)/2*Ma(1,k+1)2); M(1,k+1)=P(1,k+1)/R/T(1,k+1); if kpoint1 V(1,k+1)=0; else V(1
21、,k+1)=-U(1,k+1)*tan(sita); end F1(1)=M(1,k+1)*U(1,k+1); F2(1)=M(1,k+1)*U(1,k+1)2+P(1,k+1); F3(1)=M(1,k+1)*U(1,k+1)*V(1,k+1); F4(1)=r/(r-1)*P(1,k+1)*U(1,k+1)+M(1,k+1)*U(1,k+1)*(U(1,k+1)2+V(1,k+1)2)/2; kend计算结果后处理close all;%横向速度point=5*floor(1:20:130);figurefor i=1:length(point)subplot(1,length(point)
22、,i) plot(U(:,point(i),y(:,point(i),-,linewidth,2)ylabel(y/cm)xlabel(u/(m/s)title(strcat(x=,num2str(point(i)axis(660 760 -15 40)% axis offend%纵向速度point=5*floor(1:20:130);figurefor i=1:length(point)subplot(1,length(point),i) plot(V(:,point(i),y(:,point(i),-,linewidth,2)ylabel(y/cm)xlabel(v/(m/s)title(
23、strcat(x=,num2str(point(i)axis(-250 50 -15 40)% axis offend%马赫数point=5*floor(1:20:130);figurefor i=1:length(point)subplot(1,length(point),i) plot(Ma(:,point(i),y(:,point(i),-,linewidth,2)ylabel(y/cm)xlabel(Ma)title(strcat(x=,num2str(point(i)axis(1 3 -15 40)% axis offend% 密度point=5*floor(1:20:130);fi
24、gurefor i=1:length(point)subplot(1,length(point),i) plot(M(:,point(i),y(:,point(i),-,linewidth,2)ylabel(y/cm)xlabel(M/(kg/m3)title(strcat(x=,num2str(point(i)axis(0.5 1.5 -15 40)% axis offend% 压力point=5*floor(1:20:130);figurefor i=1:length(point)subplot(1,length(point),i) plot(P(:,point(i),y(:,point(
25、i),-,linewidth,2)ylabel(y/cm)xlabel(P/(Pa)title(strcat(x=,num2str(point(i)axis(0.1*105 2*105 -15 40)% axis offend%温度point=5*floor(1:20:130);figurefor i=1:length(point)subplot(1,length(point),i) plot(T(:,point(i),y(:,point(i),-,linewidth,2)ylabel(y/cm)xlabel(T/(k)title(strcat(x=,num2str(point(i)axis(200 300 -15 40)% axis offEnd计算结果（粘性0.2+角度5度）
