一、Matlab代碼片

%全亞聲速等熵噴管流動 非守恒型麥考馬克方法數值求解 clear; %清理內存變量 clc; %清理工作窗中的所有顯示內容 r=1.4; %比熱比 L=3; %噴管長度 i=31; %網格數目 C=0.5; %科朗數 x=linspace(0,L,i); %網格點橫坐標 N=5001; %時間步 dx=L/(i-1); %空間步長 dt(1:N)=0; %時間步長 Pe=0.93; %噴管兩端壓力比 A(1:31)=0; %預分配內存 for j=1:31 if(x(j)<1.5)A(j)=1+2.2*(x(j)-1.5).*(x(j)-1.5); %噴管面積 elseA(j)=1+0.2223*(x(j)-1.5).*(x(j)-1.5);%噴管面積 end end 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.023*x; %流場密度初值 T(1,:)=1-0.009333*x; %流場溫度初值 V(1,:)=0.05+0.11*x; %流場速度初值 %按時間步長推進 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,2:i-1)+sqrt(T(k,2:i-1)));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)= Pe./rho(k+1,i);%馬赫數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:5001,rho(1:5001,16),'b-'); ylabel('\rho/\rho_0'),xlabel('Timestep'); title('無量綱密度變化'); subplot(2,2,2),plot(1:5001,T(1:5001,16),'r-'); ylabel('T/T_0'),xlabel('Timestep');title('無量綱溫度變化'); subplot(2,2,3),plot(1:5001,P(1:5001,16),'k-'); ylabel('P/P_0'),xlabel('Timestep');title('無量綱總壓變化'); subplot(2,2,4),plot(1:5001,Ma(1:5001,16),'m-'); ylabel('Ma'),xlabel('Timestep');title('馬赫數變化'); %繪圖 噴管喉道處無量綱密度和速度時間導數的變化 figure; plot(1:5000,rhotav2(1:5000,16),'k-'); ylabel('殘差'),xlabel('Timestep'); title('噴管喉道處無量綱密度和速度時間導數的變化'); hold on; plot(1:5000,Vtav2(1:5000,16),'k--'); %繪圖 無量綱質量流量在三個不同時刻的瞬時分布 figure; plot(x,rho(1,1:i).*V(1,1:i).*A(1:i),'k--'); hold on; plot(x,rho(501,1:i).*V(501,1:i).*A(1:i),'r-'); hold on; plot(x,rho(5001,1:i).*V(5001,1:i).*A(1:i),'b-'); title('不同刻質量流量的變化'); ylabel('無量綱質量流量'),xlabel('無量綱軸向距離'); text(3.2,0.47,'5000\Deltat'); text(3.2,0.40,'500\Deltat'); text(2.0,0.25,'0\Deltat'); %繪圖 不同時刻壓力分布的變化 figure; plot(x,rho(1,1:i).*T(1,1:i),'k--'); hold on; plot(x,rho(1001,1:i).*T(1001,1:i),'m-'); hold on; plot(x,rho(501,1:i).*T(501,1:i),'r-'); hold on; plot(x,rho(5001,1:i).*T(5001,1:i),'b-'); hold on; title('不同時刻壓力分布的變化'); ylabel('P/P_0'),xlabel('x/L'); text(1.5,0.96,'0\Deltat'); text(2.0,0.91,'500\Deltat'); text(2.2,0.89,'1000\Deltat'); text(2.3,0.867,'5000\Deltat');

二、計算結果展示

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

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

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

4.無量綱壓力在四個不同時刻的瞬時分布

總結

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

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