擬一維噴管流動的數(shù)值解法(MATLAB)代碼

1、-作者xxxx-日期xxxx擬一維噴管流動的數(shù)值解法(MATLAB)代碼【精品文檔】擬一維噴管流動的數(shù)值解法（MATLAB）代碼數(shù)值計算代碼%擬一維噴管流動的數(shù)值解%亞聲速-超聲速，非守恒形式function main()clear;clc;r=1.4; %絕熱指數(shù)N=1001; %時間步長i=31; %網(wǎng)格數(shù)目L=3; %噴管長度C=0.5; %柯朗數(shù)dx=L/(i-1); %空間步長dt(N)=0; %時間步長x=linspace(0,L,i); %網(wǎng)格點橫坐標A=1+2.2*(x-1.5).2; %噴管面積%賦值M(N,i)=0;T(N,i)=0;V(N,i)=0;%初始條件M(1,:)

2、=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 %預估偏導數(shù) M_t(1:i-1)=-V(k,1:i-1).*(M(k,2:i)-M(k,1:i-1)/dx-M(k,1:i-1).*(V(k,2:i)-V(k,1:i-1)/dx-M(k,1:i-1).*V(k,1:i-1).*log(A(2:i)./A(1:i-1)/dx; V_t(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

3、)/dx+T(k,1:i-1)./M(k,1:i-1).*(M(k,2:i)-M(k,1:i-1)/dx); T_t(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)./A(1:i-1)/dx); %求取內(nèi)部網(wǎng)格點處最小時間步長 t=C*dx./(V(k,2:i-1)+sqrt(T(k,2:i-1); dt(k)=min(t); %預估值 M1(1:i-1)=M(k,1:i-1)+M_t(1:i-1)*dt(k); V1(1:

4、i-1)=V(k,1:i-1)+V_t(1:i-1)*dt(k); T1(1:i-1)=T(k,1:i-1)+T_t(1:i-1)*dt(k); %校正偏導數(shù) M_t_1(2:i-1)=-V1(2:i-1).*(M1(2:i-1)-M1(1:i-2)./dx-M1(2:i-1).*(V1(2:i-1)-V1(1:i-2)./dx-M1(2:i-1).*V1(2:i-1).*log(A(2:i-1)./A(1:i-2)./dx; V_t_1(2:i-1)=-V1(2:i-1).*(V1(2:i-1)-V1(1:i-2)./dx-1/r.*(T1(2:i-1)-T1(1:i-2)./dx+T1(2

5、:i-1)./M1(2:i-1).*(M1(2:i-1)-M1(1:i-2)./dx); T_t_1(2:i-1)=-V1(2:i-1).*(T1(2:i-1)-T1(1:i-2)./dx-(r-1).*T1(2:i-1).*(V1(2:i-1)-V1(1:i-2)./dx+V1(2:i-1).*log(A(2:i-1)./A(1:i-2)./dx); %偏導數(shù)平均值 M_t_av(2:i-1)=0.5*(M_t(2:i-1)+M_t_1(2:i-1); V_t_av(2:i-1)=0.5*(V_t(2:i-1)+V_t_1(2:i-1); T_t_av(2:i-1)=0.5*(T_t(2:i

6、-1)+T_t_1(2:i-1); %內(nèi)部網(wǎng)格點修正值 M(k+1,2:i-1)=M(k,2:i-1)+M_t_av(2:i-1)*dt(k); V(k+1,2:i-1)=V(k,2:i-1)+V_t_av(2:i-1)*dt(k); T(k+1,2:i-1)=T(k,2:i-1)+T_t_av(2:i-1)*dt(k); %出口邊界值 M(k+1,i)=2*M(k+1,i-1)-M(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); %入口邊界值 M(k+1,1)=1; V(k+1,1)=

7、2*V(k+1,2)-V(k+1,3); T(k+1,1)=1;endend圖形處理代碼close all;a=1 51 101 151 201 701;plot(x,M(a(1),:).*A(:).*V(a(1),:),r-)hold onplot(x,M(a(2),:).*A(:).*V(a(2),:),m.-)plot(x,M(a(3),:).*A(:).*V(a(3),:),g*-)plot(x,M(a(4),:).*A(:).*V(a(4),:),co-)plot(x,M(a(5),:).*A(:).*V(a(5),:),bh-)plot(x,M(a(6),:).*A(:).*V(a

8、(6),:),k-)axis equallegend(0dt,50dt,100dt,150dt,200dt,700dt)figureax,h1,h2=plotyy(x,M(a(6),:),x,sqrt(M(a(6),:).(1-r)-1)*5),plot);set(h1,linestyle,-,marker,o,color,b); set(h2,linestyle,-,marker,*,color,r);set(ax,xtick,0:0.3:3)set(ax(1),ytick,0:0.1:1)set(ax(2),ytick,0:0.4:4)set(h1,linewidth,1.5)set(h2,linewidth,1.5)set(get(ax(1),xlabel),string,無量綱軸向距離