Matlab 程式設計 《計算流體力學基礎及應用（約翰D安德森）》 亞聲速-超聲速等熵噴管流動CFD解法 拉瓦爾噴管 守恒形式方程解法

《計算流體力學基礎及應用（約翰D安德森）》 亞聲速-超聲速等熵噴管流動CFD解法，拉瓦爾噴管 守恒形式方程解法 Matlab 程式設計

• 問題之 亞聲速一超聲速等熵噴管流動的CFD解法
• • 初始化參數
  • 疊代過程
  • 繪圖
  • 結果

最近學習經典的計算流體力學入門書籍《計算流體力學基礎及應用（約翰D安德森）》，跟着作者給出的過程寫了計算程式，程式寫的比較初級。

問題之 亞聲速一超聲速等熵噴管流動的CFD解法

問題的提法見《計算流體力學基礎及應用（約翰D安德森）》中的7.3節，下面直接給出寫的程式。

初始化參數

%亞聲速-超聲速等熵噴管流動CFD解法，拉瓦爾噴管 守恒形式方程解法 
%MacCormack方法 預估-校正 具有二階精度
%初始參數（無量綱）
clc;clear;clf
x=0:0.1:3;
gm=1.4;%比熱比
C=0.5;%科朗數
N=length(x);%節點數量
A=1+2.2.*(x-1.5).^2;%截面函數
rho=[ones(1,6),1.0-0.366.*(x(7:16)-0.5),0.634-0.3879.*(x(17:31)-1.5)];%密度
T=[ones(1,6),1.0-0.167.*(x(7:16)-0.5),0.833-0.3507.*(x(17:31)-1.5)];%溫度
U1=rho.*A;%控制方程因變量
U2=repmat(0.59,1,31);%控制方程因變量
V=U2./U1;%速度
U3=U1.*(T./(gm-1)+gm/2.*V.^2);%控制方程因變量
U1_t=zeros(1,N);%因變量1一階偏導
U1_mid=U1;%預估步因變量1
U1_tt=zeros(1,N);%校正步因變量1
U1_av=zeros(1,N);%校正因變量1
U2_t=zeros(1,N);%因變量2一階偏導
U2_mid=U2;%預估步因變量2
U2_tt=zeros(1,N);%校正步因變量2
U2_av=zeros(1,N);%校正因變量1
U3_t=zeros(1,N);%因變量3一階偏導
U3_mid=U3;%預估步因變量2
U3_tt=zeros(1,N);%校正步因變量3
U3_av=zeros(1,N);%校正因變量1
F1=U2;%通量1
F1_mid=U2;%預估步通量1初始化
F2=U2.^2./U1+(gm-1)/gm.*(U3-gm/2.*U2.^2./U1);%通量2
F2_mid=U2.^2./U1+(gm-1)/gm.*(U3-gm/2.*U2.^2./U1);%預估步通量2初始化
F3=gm.*U2.*U3./U1-gm*(gm-1)/2.*U2.^3./U1.^2;%通量3
F3_mid=gm.*U2.*U3./U1-gm*(gm-1)/2.*U2.^3./U1.^2;%預估步通量3初始化
J2=zeros(1,N);%源項2
J2_mid=zeros(1,N);%預估步源項2初始化
rho_mid=zeros(1,N);%密度預估值
T_mid=zeros(1,N);%溫度預估值
t=zeros(1,N);%局部步長
k=0;%疊代次數初始化
k0=2000;%指定疊代次數
           

疊代過程

%疊代過程
while 1
    k=k+1;
% 計算預估步 向前差分
for i=1:N
    t(i)=C.*0.1./(V(i)+sqrt(T(i)));%計算滿足穩定性條件的各個局部步長
end
delta_t=min(t);%取最小步長
for i=2:N-1
    J2(i)=1/gm.*rho(i).*T(i).*(A(i+1)-A(i))/0.1;%預估步源項2
    U1_t(i)=-(F1(i+1)-F1(i))/0.1;%計算因變量1一階偏導
    U2_t(i)=-(F2(i+1)-F2(i))/0.1+J2(i);%計算因變量2一階偏導
    U3_t(i)=-(F3(i+1)-F3(i))/0.1;%計算因變量3一階偏導;
    U1_mid(i)=U1(i)+U1_t(i).*delta_t;%因變量1預估
    U2_mid(i)=U2(i)+U2_t(i).*delta_t;%因變量2預估
    U3_mid(i)=U3(i)+U3_t(i).*delta_t;%因變量3預估
    rho_mid(i)=U1_mid(i)./A(i);%密度預估
    T_mid(i)=(gm-1).*(U3_mid(i)./U1_mid(i)-gm/2.*U2_mid(i).^2./U1_mid(i).^2);%溫度預估
    F1_mid(i)=U2_mid(i);%預估步通量1
    F2_mid(i)=U2_mid(i).^2/U1_mid(i)+(gm-1)/gm.*(U3_mid(i)- ...
        gm/2.*U2_mid(i).^2./U1_mid(i));%預估步通量2
    F3_mid(i)=gm.*U2_mid(i).*U3_mid(i)./U1_mid(i)-gm*(gm-1)/ ...
        2.*U2_mid(i).^3./U1_mid(i).^2;%預估步通量3
end
U2_mid(1)=2.*U2_mid(2)-U2_mid(3);%入口點1處通過點2 3的插值獲得
rho_mid(1)=U1_mid(1)./A(1);
V_mid(1)=U2_mid(1)./U1_mid(1);
U3_mid(1)=U1_mid(1)*(1/(gm-1)+gm/2*V_mid(1)^2);
F1_mid(1)=U2_mid(1);
F2_mid(1)=U2_mid(1).^2./U1_mid(1)+(gm-1)/gm.*(U3_mid(1)-gm/2.*U2_mid(1).^2./U1_mid(1));
F3_mid(1)=gm.*U2_mid(1).*U3_mid(1)./U1_mid(1)-gm*(gm-1)/2.*U2_mid(1).^3./U1_mid(1).^2;
%計算校正步 向後差分
for i=2:N-1
    J2_mid(i)=1/gm.*rho_mid(i).*T_mid(i).*(A(i)-A(i-1))/0.1;%校正步源項2
    U1_tt(i)=-(F1_mid(i)-F1_mid(i-1))/0.1;%計算因變量1一階偏導
    U2_tt(i)=-(F2_mid(i)-F2_mid(i-1))/0.1+J2_mid(i);%計算因變量2一階偏導
    U3_tt(i)=-(F3_mid(i)-F3_mid(i-1))/0.1;%計算因變量3一階偏導;
    U1_av(i)=0.5.*(U1_t(i)+U1_tt(i));%變量1校正
    U2_av(i)=0.5.*(U2_t(i)+U2_tt(i));%變量2校正
    U3_av(i)=0.5.*(U3_t(i)+U3_tt(i));%變量3校正
end
U1=U1+U1_av.*delta_t;%下一時間步
U1(N)=2.*U1(N-1)-U1(N-2);%出口處點N通過插值獲得
U2=U2+U2_av.*delta_t;%下一時間步
U2(1)=2.*U2(2)-U2(3);%入口點1處通過點2 3的插值獲得
U2(N)=2.*U2(N-1)-U2(N-2);%出口處點N通過插值獲得
rho=U1./A;
V=U2./U1;
U3=U3+U3_av.*delta_t;%下一時間步
U3(1)=U1(1)*(1/(gm-1)+gm/2*V(1)^2);
U3(N)=2.*U3(N-1)-U3(N-2);%出口處點N通過插值獲得
T=(gm-1).*(U3./U1-gm/2.*V.^2);
p=rho.*T;
Ma=V./sqrt(T);
F1=U2;%通量1
F2=U2.^2./U1+(gm-1)/gm.*(U3-gm/2.*U2.^2./U1);%通量2
F3=gm.*U2.*U3./U1-gm*(gm-1)/2.*U2.^3./U1.^2;%通量3
%喉道參數擷取
rhot(k)=rho(16);
Vt(k)=V(16);
Tt(k)=T(16);
Mat(k)=Ma(16);
pt(k)=p(16);
re_rho(k)=abs(rho(16));
%繪制不同時間步品質流量變化
if k==1
    xtick=1:N;
    rhoAV=rho(xtick).*A(xtick).*V(xtick);
    plot(xtick/10,rhoAV,'-b')
    hold on
elseif k==500
    xtick=1:N;
    rhoAV=rho(xtick).*A(xtick).*V(xtick);
    plot(xtick/10,rhoAV,'-r')
    hold on
elseif k==1000
    xtick=1:N;
    rhoAV=rho(xtick).*A(xtick).*V(xtick);
    plot(xtick/10,rhoAV,'-g')
    xlabel('x/L'),ylabel('\rhoAV/\rho_0A_1a_0')
    title('不同時刻品質流量的變化')
    text(2,rhoAV(N)-0.0025,'1000\Deltat')
    text(2,rhoAV(N)+0.0025,'500\Deltat')
    text(2,0.602,'1\Deltat')
    hold on
else
    
end
%疊代次數
if k==k0
    break;
else
    
end
end
           

繪圖

figure
x0=1:1:k;%時間步
%繪制密度-時間步圖像
rho_change=rhot(x0);
subplot(2,2,1),plot(x0,rho_change)
ylabel('\rho/\rho_0'),xlabel('時間步數')
title('無量綱密度變化')
%繪制溫度-時間步圖像
T_change=Tt(x0);
subplot(2,2,2),plot(x0,T_change)
ylabel('T/T_0'),xlabel('時間步數')
title('無量綱溫度變化')
%繪制壓力-時間步圖像
p_change=pt(x0);
subplot(2,2,3),plot(x0,p_change)
ylabel('p/p_0'),xlabel('時間步數')
title('無量綱總壓變化')
%繪制馬赫數-時間步圖像
Ma_change=Mat(x0);
subplot(2,2,4),plot(x0,Ma_change)
ylabel('Ma'),xlabel('時間步數')
title('馬赫數變化')
fprintf('疊代最終結果:\n密度rho=%6.4e\n溫度T=%6.4e\n速度V=%6.4e\n總壓p=%6.4e\n馬赫數Ma=%6.4e\n' ...
    ,rhot(k),Tt(k),pt(k),Vt(k),Mat(k))
figure 
re_rho_change=re_rho(x0);
semilogy(x0,re_rho_change)
ylabel('殘差'),xlabel('時間步數')
title('殘差随時間步變化')
x00=1:1:N;
figure
p0=p(x00);
subplot(2,2,1)
plot(x00/10,p0)
ylabel('p/p_0'),xlabel('x/L')
title('噴管内壓力分布')
Ma0=Ma(x00);
subplot(2,2,2)
plot(x00/10,Ma0)
ylabel('Ma'),xlabel('x/L')
title('噴管内馬赫數分布')
T0=T(x00);
subplot(2,2,3)
plot(x00/10,T0)
ylabel('T/T_0'),xlabel('x/L')
title('噴管内溫度分布')
rho0=rho(x00);
subplot(2,2,4)
plot(x00/10,rho0)
ylabel('\rho/\rho_0'),xlabel('x/L')
title('噴管内密度分布')
           

結果

疊代最終結果:

密度rho=6.4917e-01

溫度T=8.3937e-01

速度V=5.4489e-01

總壓p=9.0153e-01

馬赫數Ma=9.8401e-01

Figure 1

Figure 2

Figure 3

Figure 4