维超声速流动的数值解普朗特迈耶稀疏波.docx
- 文档编号:24309215
- 上传时间:2023-05-26
- 格式:DOCX
- 页数:15
- 大小:241.83KB
维超声速流动的数值解普朗特迈耶稀疏波.docx
《维超声速流动的数值解普朗特迈耶稀疏波.docx》由会员分享,可在线阅读,更多相关《维超声速流动的数值解普朗特迈耶稀疏波.docx(15页珍藏版)》请在冰豆网上搜索。
维超声速流动的数值解普朗特迈耶稀疏波
二维超声速流动的数值解——普朗特迈耶稀疏波
主程序数值计算
%二维超声速流动的数值解——普朗特迈耶稀疏波
clc;
clear;
%常量设置
CFL=;
Cy=;
r=;
R=287;
%初始物理条件
Ma0=2;
P0=*10A5;
M0=;
T0=;
%流场几何外形描述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;
fork=1:
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;
M(:
1)=M0;
P(:
1)=P0;
T(:
1)=T0;
Ma(:
1)=Ma0;
F1=M(:
1).*U(:
1);
F2=M(:
1).*U(:
1).A2+P(:
1);
F3=M(:
1).*U(:
1).*V(:
1);
F4=r/(r-1)*P(:
1).*U(:
1)+M(:
1).*U(:
1).*(U(:
1).A2+V(:
1).A2)/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);
%麦考马克方法空间推进计算fork=1:
i
G1=M(:
k).*F3./F1;
G2=F3;
G3=M(:
k).*(F3./F1).A2+F2-F142./M(:
k);
G4=r/(r-1)*(F2-F192./M(:
k)).*F3./F1+M(:
k)/2.*F3./F1.*((F1./M(:
k)).A2+(F3./F1)42);
%预估偏导数,前向差分
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;
%计算空间空间步长计算
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,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*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:
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)42./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).A3;
M_1(1:
j-1)=(-B(1:
j-1)+sqrt(B(1:
j-1).A2-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).A2./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)).A2+F2_1(1:
j-1)-F1_1(1:
j-1).A2./M_1(1:
j-1);
G4_1(1:
j-1)=r/(r-1)*(F2_1(1:
j-1)-F1_1(1:
j-1).A2./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)).A2+(F3_1(1:
j-1)./F1_1(1:
j-1)).A2);
%修正偏导数计算
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).*(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)=*(F1_eps_1(2:
j-1)+F1_eps(2:
j-1));
F2_eps_av(2:
j-1)=*(F2_eps_1(2:
j-1)+F2_eps(2:
j-1));
F3_eps_av(2:
j-1)=*(F3_eps_1(2:
j-1)+F3_eps(2:
j-1));
F4_eps_av(2:
j-1)=*(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(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)=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:
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);
%下边界单侧差分
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=F342./F1/2-F4;
B=r/(r-1)*F1.*F2;
C=-(r+1)/(r-1)/2*F193;
M(:
k+1)=(-B+sqrt(B.A2-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;
Ma(:
k+1)=sqrt(U(:
k+1).A2+V(:
k+1).A2)./sqrt(r*R*T(:
k+1));
%下边界处物理量修正
Ma_cal=sqrt(U(1,k+1)A2+V(1,k+1)A2)/sqrt(r*R*T(1,k+1));f_cal=sqrt((r+1)/(r-1))*atan(sqrt((r-1)*(Ma_calA2-1)/(r+1)))-atan(sqrt(Ma_calA2-1));
ifk fi=atan(abs(V(1,k+1))/U(1,k+1)); else psai=atan(abs(V(1,k+1))/U(1,k+1)); fi=sita-psai; end f_act=f_cal+fi; %二分法求马赫数 xm=[15];n=0; whileabs(xm (1)-xm(3))>=&&n<100 fm=sqrt((r+1)/(r-1))*atan(sqrt((r-1)/(r+1)*(xm.A2-1)))-atan(sqrt(xm.A2-1))-f_act;iffm (2)*fm(3)<0 xm=[xm (2)(xm (2)+xm(3))/2xm(3)]; elseiffm (2)*fm(3)>0xm=[xm (1)(xm (1)+xm (2))/2xm (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_calA2)/(1+(r-1)/2*Ma(1,k+1)A2))A(r/(r-1));T(1,k+1)=T(1,k+1)*((1+(r-1)/2*Ma_calA2)/(1+(r-1)/2*Ma(1,k+1)A2)); M(1,k+1)=P(1,k+1)/R/T(1,k+1); ifk V(1,k+1)=0; else V(1,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)A2+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)A2+V(1,k+1)A2)/2;k end 计算结果后处理 closeall; %横向速度 point=5*floor(1: 20: 130); figure fori=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([660760-1540]) %axisoff end %纵向速度 point=5*floor(1: 20: 130); figure fori=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([-25050-1540]) %axisoff end %马赫数 point=5*floor(1: 20: 130); figure fori=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([13-1540]) %axisoff end %密度 point=5*floor(1: 20: 130); figure fori=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([-1540]) %axisoff end %压力 point=5*floor(1: 20: 130); figure fori=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([*10A52*10A5-1540]) %axisoff end %温度 point=5*floor(1: 20: 130); figure fori=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([200300-1540]) %axisoff End 计算结果(粘性+角度5度) 0 ■ ■ ■ - ■- - 1-5'■ - 1□- - 5■ - 0• &1: ? £'Ve rwu X 30 M : d O r P: -E1 D JO 2Q 2D W: 1 呼ill」」 HQ ■-» Wlkii'irU hl他旷曰 F*V|k>.||l! >」 trihJi 讯化屮di »-1k-41»-IO 伽TAMg映)TAh)Wo •-R X.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 超声速 流动 数值 解普朗特迈耶 稀疏