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=0.00001elsexm=xm(2) xm(2) xm(2);endn=n+1;endMa(1,k+1)=xm(2);%下边界处物理量修正值P(1,k+
2、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 kpoint1V(1,k+1)=0;elseV(1,k+1)=-U(1,k+1)*tan(sita);endF1(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)
3、;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),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 4
4、0)% 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(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,l
5、ength(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);figurefor 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)
6、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(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 度)