一、Matlab代碼片

%亞聲速-超聲速等熵噴管流動 非守恒型麥考馬克方法數值求解 clear; %清理內存變量 clc; %清理工作窗中的所有顯示內容 r=1.4; %比熱比 L=3; %噴管長度 i=31; %網格數目 C=0.5; %科朗數 x=linspace(0,L,i); %網格點橫坐標 N=1401; %時間步 dx=L/(i-1); %空間步長 dt(1:N)=0; %時間步長 A=1+2.2*(x-1.5).*(x-1.5); %噴管面積 rho(N,i)=0; %流場密度賦值 T(N,i)=0; %流場溫度賦值 V(N,i)=0; %流場速度賦值 %預分配內存 rhotav2(1:1400,1:30)=0; Vtav2(1:1400,1:30)=0; %初始條件 rho(1,:)=1-0.3146*x; %流場密度初值 T(1,:)=1-0.2314*x; %流場溫度初值 V(1,:)=(0.1+1.09*x).*(1-0.2314*x).^0.5; %流場速度初值 %按時間步長推進 for k=1:N-1%計算預估步偏導數rhot(1:i-1)=-V(k,1:i-1).*(rho(k,2:i)-rho(k,1:i-1))/dx-rho(k,1:i-1)....*(V(k,2:i)-V(k,1:i-1))/dx-rho(k,1:i-1).*V(k,1:i-1).*(log(A(2:i))-log(A(1:i-1)))/dx;Vt(1:i-1)=-V(k,1:i-1).*(V(k,2:i)-V(k,1:i-1))/dx-1/r.*((T(k,2:i)...-T(k,1:i-1))/dx+T(k,1:i-1)./rho(k,1:i-1).*(rho(k,2:i)-rho(k,1:i-1))/dx);Tt(1:i-1)=-V(k,1:i-1).*(T(k,2:i)-T(k,1:i-1))/dx-(r-1).*T(k,1:i-1)....*((V(k,2:i)-V(k,1:i-1))/dx+V(k,1:i-1).*(log(A(2:i))-log(A(1:i-1)))/dx);%確定最小時間步長t=C*dx./(V(k,:)+sqrt(T(k,:)));dt(k)=min(t);%計算預估值rho_(1:i-1)=rho(k,1:i-1)+rhot(1:i-1).*dt(k);V_(1:i-1)=V(k,1:i-1)+Vt(1:i-1).*dt(k);T_(1:i-1)=T(k,1:i-1)+Tt(1:i-1).*dt(k);%校正偏導數rhot_(2:i-1)=-V_(2:i-1).*(rho_(2:i-1)-rho_(1:i-2))/dx-rho_(2:i-1)....*(V_(2:i-1)-V_(1:i-2))/dx-rho_(2:i-1).*V_(2:i-1).*(log(A(2:i-1))-log(A(1:i-2)))/dx;Vt_(2:i-1)=-V_(2:i-1).*(V_(2:i-1)-V_(1:i-2))/dx-1/r.*((T_(2:i-1)...-T_(1:i-2))/dx+T_(2:i-1)./rho_(2:i-1).*(rho_(2:i-1)-rho_(1:i-2))/dx);Tt_(2:i-1)=-V_(2:i-1).*(T_(2:i-1)-T_(1:i-2))/dx-(r-1).*T_(2:i-1)....*((V_(2:i-1)-V_(1:i-2))/dx+V_(2:i-1).*(log(A(2:i-1))-log(A(1:i-2)))/dx);%時間導數平均值rhotav(2:i-1)=0.5*(rhot(2:i-1)+rhot_(2:i-1)); rhotav2(k,2:i-1)=abs(rhotav(2:i-1)); Vtav(2:i-1)=0.5*(Vt(2:i-1)+Vt_(2:i-1));Vtav2(k,2:i-1)=abs(Vtav(2:i-1));Ttav(2:i-1)=0.5*(Tt(2:i-1)+Tt_(2:i-1));%流場變量校正值rho(k+1,2:i-1)=rho(k,2:i-1)+rhotav(2:i-1)*dt(k);V(k+1,2:i-1)=V(k,2:i-1)+Vtav(2:i-1)*dt(k);T(k+1,2:i-1)=T(k,2:i-1)+Ttav(2:i-1)*dt(k);%入流邊界值V(k+1,1)=2*V(k+1,2)-V(k+1,3);rho(k+1,1)=1;T(k+1,1)=1;%出流邊界值rho(k+1,i)=2*rho(k+1,i-1)-rho(k+1,i-2);V(k+1,i)=2*V(k+1,i-1)-V(k+1,i-2);T(k+1,i)=2*T(k+1,i-1)-T(k+1,i-2);%馬赫數Ma=V(1:k+1,1:i)./(sqrt(T(1:k+1,1:i)));%流場壓強P=rho(1:k+1,1:i).*T(1:k+1,1:i); end %繪圖 噴管喉道處密度、溫度、壓力和馬赫數的變化 figure; subplot(2,2,1),plot(1:1401,rho(1:1401,16),'b-'); ylabel('\rho/\rho_0'),xlabel('Timestep'); title('無量綱密度變化'); subplot(2,2,2),plot(1:1401,T(1:1401,16),'r-'); ylabel('T/T_0'),xlabel('Timestep');title('無量綱溫度變化'); subplot(2,2,3),plot(1:1401,P(1:1401,16),'k-'); ylabel('P/P_0'),xlabel('Timestep');title('無量綱總壓變化'); subplot(2,2,4),plot(1:1401,Ma(1:1401,16),'m-'); ylabel('Ma'),xlabel('Timestep');title('馬赫數變化'); %繪圖 噴管喉道處無量綱密度和速度時間導數的變化 figure; plot(1:1400,rhotav2(1:1400,16),'k-'); ylabel('殘差'),xlabel('Timestep'); title('噴管喉道處無量綱密度和速度時間導數的變化'); hold on; plot(1:1400,Vtav2(1:1400,16),'k--'); %繪圖 無量綱質量流量在六個不同時刻的瞬時分布 figure; plot(x,rho(1,1:i).*V(1,1:i).*A(1:i),'k--'); hold on; plot(x,rho(51,1:i).*V(51,1:i).*A(1:i),'r-'); hold on; plot(x,rho(101,1:i).*V(101,1:i).*A(1:i),'b-'); hold on; plot(x,rho(151,1:i).*V(151,1:i).*A(1:i),'g-'); hold on; plot(x,rho(201,1:i).*V(201,1:i).*A(1:i),'y-'); hold on; plot(x,rho(701,1:i).*V(701,1:i).*A(1:i),'m-'); title('不同時刻質量流量的變化'); ylabel('無量綱質量流量'),xlabel('無量綱軸向距離'); text(2.4,0.63,'700\Deltat'); text(2.4,0.75,'150\Deltat'); text(2.4,0.53,'200\Deltat'); text(2.4,0.9,'50\Deltat'); text(2.4,1.05,'100\Deltat'); text(2.4,1.3,'0\Deltat');

二、計算結果展示

1.噴管喉道處密度、溫度、壓力和馬赫數的變化

2.噴管喉道處無量綱密度和速度時間導數的變化

3.無量綱質量流量在六個不同時刻的瞬時分布

總結

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

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