一、Matlab代碼片

%亞聲速-超聲速等熵噴管流動守恒形CFD解法 MacCormack方法 clear; %清理內存變量 clc; %清理工作窗中的所有顯示內容 r=1.4; %比熱比 L=3; %噴管長度 i=31; %網格數目 C=0.5; %科朗數 x=linspace(0,L,i); %網格點橫坐標 N=1401; %時間步 %預分配內存 rho=zeros(1401,31); V=zeros(1401,31); T=zeros(1401,31); U1=zeros(1401,31); U2=zeros(1401,31); U3=zeros(1401,31); dx=L/(i-1); %空間步長 dt(1:N)=0; %時間步長 A=1+2.2*(x-1.5).*(x-1.5); %噴管面積 %初始條件 rho(1,:)=[ones(1,6),1.0-0.366.*(x(7:16)-0.5),0.634-0.3879.*(x(17:31)-1.5)]; T(1,:)=[ones(1,6),1.0-0.167.*(x(7:16)-0.5),0.833-0.3507.*(x(17:31)-1.5)]; U1(1,:)=rho(1,:).*A(1:31); U2(1,:)=0.59; V(1,:)=U2(1,:)./U1(1,:); U3(1,:)=U1(1,:).*((T(1,:)./(r-1))+r/2.*V(1,:).^2); %MacCormack方法迭代 for k=1:N-1%確定最小時間步長t=C*dx./(V(k,:)+sqrt(T(k,:)));dt(k)=min(t);%預估步通量F1(1:i)=U2(k,1:i);F2(1:i)=U2(k,1:i).^2./U1(k,1:i)+(r-1)/r.*...(U3(k,1:i)-r/2.*U2(k,1:i).^2./U1(k,1:i));F3(1:i)=r.*U2(k,1:i).*U3(k,1:i)./U1(k,1:i)...-r*(r-1)/2.*U2(k,1:i).^3./U1(k,1:i).^2;J2(1:i-1)=1/r.*rho(k,1:i-1).*T(k,1:i-1).*(A(2:i)-A(1:i-1))/0.1; %源項%預估偏導數U1t(1:i-1)=-(F1(2:i)-F1(1:i-1))/dx;U2t(1:i-1)=-(F2(2:i)-F2(1:i-1))/dx+J2(1:i-1);U3t(1:i-1)=-(F3(2:i)-F3(1:i-1))/dx;%預估值U1_(1:i-1)=U1(k,1:i-1)+U1t(1:i-1).*dt(k);U2_(1:i-1)=U2(k,1:i-1)+U2t(1:i-1).*dt(k);U3_(1:i-1)=U3(k,1:i-1)+U3t(1:i-1).*dt(k);rho_(1:i-1)=U1_(1:i-1)./A(1:i-1);T_(1:i-1)=(r-1).*( U3_(1:i-1)./U1_(1:i-1)-r/2.*U2_(1:i-1).^2./U1_(1:i-1).^2);F1_(1:i-1)=U2_(1:i-1);F2_(1:i-1)=U2_(1:i-1).^2./U1_(1:i-1)+(r-1)/r.*...(U3_(1:i-1)-r/2.*U2_(1:i-1).^2./U1_(1:i-1));F3_(1:i-1)=r.*U2_(1:i-1).*U3_(1:i-1)./U1_(1:i-1)...-r*(r-1)/2.*U2_(1:i-1).^3./U1_(1:i-1).^2;J2_(2:i-1)=1/r.*rho_(2:i-1).*T_(2:i-1).*(A(2:i-1)-A(1:i-2))/0.1; %注意%校正步 校正偏導數U1t_(2:i-1)=-(F1_(2:i-1)-F1_(1:i-2))/dx;U2t_(2:i-1)=-(F2_(2:i-1)-F2_(1:i-2))/dx+J2_(2:i-1);U3t_(2:i-1)=-(F3_(2:i-1)-F3_(1:i-2))/dx;%時間導數平均值U1tav(2:i-1)=0.5.*(U1t(2:i-1)+U1t_(2:i-1));U2tav(2:i-1)=0.5.*(U2t(2:i-1)+U2t_(2:i-1));U3tav(2:i-1)=0.5.*(U3t(2:i-1)+U3t_(2:i-1));%最終校正值U1(k+1,2:i-1)=U1(k,2:i-1)+U1tav(2:i-1).*dt(k);U2(k+1,2:i-1)=U2(k,2:i-1)+U2tav(2:i-1).*dt(k);U3(k+1,2:i-1)=U3(k,2:i-1)+U3tav(2:i-1).*dt(k);rho(k+1,2:i-1)=U1(k+1,2:i-1)./A(2:i-1);V(k+1,2:i-1)=U2(k+1,2:i-1)./U1(k+1,2:i-1);T(k+1,2:i-1)=(r-1).*(U3(k+1,2:i-1)./U1(k+1,2:i-1)-r/2.*V(k+1,2:i-1).^2);%入流邊界值rho(k+1,1)=1;T(k+1,1)=1;U1(k+1,1)=A(1);U2(k+1,1)=2.*U2(k+1,2)-U2(k+1,3);V(k+1,1)=U2(k+1,1)./U1(k+1,1);U3(k+1,1)=U1(k+1,1).*((T(k+1,1)./(r-1))+r/2.*V(k+1,1).^2);%出流邊界值U1(k+1,i)=2.*U1(k+1,i-1)-U1(k+1,i-2);U2(k+1,i)=2.*U2(k+1,i-1)-U2(k+1,i-2);U3(k+1,i)=2.*U3(k+1,i-1)-U3(k+1,i-2);%由于不需要返回原變量的出流邊界值，故這里省去了原變量出流邊界值的計算 end %繪圖 無量綱質量流量在六個不同時刻的瞬時分布 figure; plot(x,U2(1,:),'k--'); hold on; plot(x,U2(101,:),'b-'); hold on; plot(x,U2(201,:),'y-'); hold on; plot(x,U2(701,:),'m-'); axis([0,3,0.35,0.7]); title('不同時刻質量流量的變化'); ylabel('無量綱質量流量'),xlabel('無量綱軸向距離'); text(2.6,0.58,'700\Deltat'); text(1.8,0.56,'200\Deltat'); text(1.5,0.66,'100\Deltat'); text(1.5,0.6,'0\Deltat');

二、計算結果展示

無量綱質量流量在四個不同時刻的瞬時分布

總結

以上是生活随笔為你收集整理的拟一维喷管流动的数值解——亚声速-超声速等熵喷管流动的守恒型CFD解法（MacCormack方法）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。