1、擬一維噴管流動的數(shù)值解法(matlab)代碼擬一維噴管流動的數(shù)值解法(matlab)代碼 編輯整理：尊敬的讀者朋友們：這里是精品文檔編輯中心，本文檔內容是由我和我的同事精心編輯整理后發(fā)布的，發(fā)布之前我們對文中內容進行仔細校對，但是難免會有疏漏的地方，但是任然希望（擬一維噴管流動的數(shù)值解法(matlab)代碼）的內容能夠給您的工作和學習帶來便利。同時也真誠的希望收到您的建議和反饋，這將是我們進步的源泉，前進的動力。本文可編輯可修改，如果覺得對您有幫助請收藏以便隨時查閱，最后祝您生活愉快 業(yè)績進步，以下為擬一維噴管流動的數(shù)值解法(matlab)代碼的全部內容。擬一維噴管流動的數(shù)值解法（matlab

2、）代碼數(shù)值計算代碼%擬一維噴管流動的數(shù)值解%亞聲速-超聲速，非守恒形式function main（)clear;clc;r=1.4； 絕熱指數(shù)n=1001； %時間步長i=31; 網格數(shù)目l=3; 噴管長度c=0。5; 柯朗數(shù)dx=l/(i1）; 空間步長dt(n)=0; 時間步長x=linspace(0,l，i）; %網格點橫坐標a=1+2.2*（x1。5).2; 噴管面積%賦值m（n，i)=0；t(n,i）=0;v（n，i）=0；%初始條件m(1，：)=10.3146x；t（1，：）=1-0。2314*x；v(1，：)=(0.1+1。09x).*（10.2314x).0。5；按時間步長推進

3、for k=1：n1 預估偏導數(shù) m_t(1：i1)=-v(k，1:i1).（m(k,2:i)m（k，1：i1)）/dx-m(k,1：i1）。(v(k,2：i）v(k，1：i-1）/dx-m(k，1:i1)。v(k，1:i1)。*log(a(2:i)。/a(1:i-1）/dx； v_t（1：i1)=v(k，1：i-1）.（v（k，2：i）v(k,1:i-1）/dx1/r。*（t(k,2：i)t(k,1：i-1)）/dx+t(k,1:i1)./m（k,1：i1）.*(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：i1)）

4、/dx(r-1）.t（k,1:i-1）.(v(k，2:i）-v（k,1:i1）)/dx+v(k,1：i-1)。log（a(2:i)./a（1:i-1)/dx）; 求取內部網格點處最小時間步長 t=cdx./(v(k,2:i1)+sqrt(t（k,2：i-1)）； dt（k)=min(t）; 預估值 m1(1：i-1)=m(k,1：i-1）+m_t(1：i1）dt(k)； v1(1：i1）=v（k,1:i-1)+v_t（1：i1)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:

5、i1)m1（1:i-2)。/dxm1（2:i1)。*(v1（2：i1)-v1（1:i-2)。/dxm1（2:i-1）。*v1(2：i1)。*log（a(2:i-1)./a(1：i2）)./dx； v_t_1（2:i1)=-v1(2:i1)。（v1（2：i-1)v1（1：i2）。/dx1/r。*(t1(2：i1）t1（1：i-2)）。/dx+t1(2:i-1)./m1(2：i1).*(m1(2：i-1)m1（1:i2)）./dx）； t_t_1（2:i1）=-v1（2：i-1).*(t1（2:i-1)-t1(1：i2)./dx-(r1)。t1(2:i1）。(v1（2:i-1)-v1（1：i-2)

6、./dx+v1（2:i1)。*log（a（2：i1）。/a（1：i2）)。/dx)； %偏導數(shù)平均值 m_t_av（2:i-1）=0.5(m_t（2：i-1）+m_t_1（2:i1)）; v_t_av（2:i1)=0.5*(v_t（2：i1)+v_t_1（2:i1)； t_t_av(2：i1)=0.5（t_t（2:i1)+t_t_1(2：i1）)； 內部網格點修正值 m(k+1，2:i1)=m(k,2：i1）+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：i1）=t(k，2：i1）+t_t_a

7、v（2：i-1）*dt（k)； %出口邊界值 m（k+1，i）=2*m(k+1，i1)-m(k+1，i-2）; v(k+1，i）=2v(k+1,i-1）v(k+1,i-2）； t(k+1,i）=2*t（k+1,i1)t（k+1,i2)； 入口邊界值 m(k+1,1)=1; v(k+1，1）=2v(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

8、。-)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（6)，：),k-)axis equallegend(0dt，50dt,100dt，150dt,200dt,700dt）figureax,h1,h2=plotyy（x，m(a(6）,：),x,sqrt（(m(a(6），:）.(1r)-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），str