1、計算傳熱學程序報告題目：一維非穩態導熱問題的數值解姓名：學號：學院：能源與動力工程學院專業：工程熱物理日期：2014年5月25日一維非穩態導熱問題數值解求解下列熱傳導問題:62T'xr=O (0x) a tT(x) = 0T(0) = l,T(L) = OL = .a = 11方程離散化對方程進行控制體積分得到:X非穩態項：選取T隨X階梯式變化，有(廠+A _ TJdx = (%7 - TlP )x擴散項：選取一階導數隨時間做顯示變化，有(¥)(#) J力=(#);(¥)!J進一步取T隨X呈分段線性變化，有竺 _TE-TPx e (x)e(OT _TP-TW整理可以

2、得到總的離散方程為:THT； T：-2T；+T：rx22計算空間和時間步長取空間步長為:h=LN網格FOUriCr數為:Fo=(小于o5時穩定)時間步長為：a3. 建立溫度矩陣與邊界條件T=OnCS(N+1,M+1)T(J)=Ti(初始條件溫度都為0)T(l,：)=TO0力界條件x=()處溫度為1)T(N+1,:)=TU(邊界條件X=L處溫度為0)4. 差分法求解溫度由離散方程可得到：T=F.(TlE-TrIPV)-TtP轉化為相應的溫度矩陣形式：Tm.k + l) = FQ T(m + l,k) + T(m-ik)-2T(m,k) + T(myk)5. 輸入界面考慮到方程的變量，采用inpu

3、tdig函數設置5個輸入變量，對這5個變量設置了默認值，如圖1所示。在計算中可以改變不同的數值，得到不同的結果，特別注 意穩定條件的臨界值是().5。根據設置的默認值，得到的計算結果如圖2所示。圖1 matbb變量輸入界面圖2默認值的計算結果6結果分析根據上面的分析，給出了程序的輸入界面，以及默認值狀態下的數值解。可以通過改變不同的輸入值，得到需要的分析結杲，總結出了下面4點結論：(I)取FR).48,得到一維非穩態導熱結果如下圖所示圖2 Fll=O.48時一維非穩態導熱從圖中可以看出，對于長度L=I的細桿，初始時刻尸0時溫度為0,邊界條 件x=0時，T=I,邊界條件x=l時，T=OO隨著時間

4、的增加，溫度從x=()通過導熱 的形式傳遞到x=l,不同時刻不同位置桿的溫度都不同，并且隨著時間的增加, 桿的溫度也逐漸增加。(2)取Fn=O.48,可以得到不同位置的溫度響應曲線，如下圖所示圖中紅色曲線代表x=()l位置的溫度瞬態響應，黑色曲線代表x=0,2位置的溫度瞬態響應，藍色曲線代表x=()4位置的溫度瞬態響應。從圖中可以看出，隨著X的增加，曲線與X軸的交點值越大，溫度開始傳遞到該位置的所需的時間越 長。隨著X的增加，溫度響應曲線的變化速率越慢，最終的達到的溫度也越低。(3)取FlI=O.25,得到不同位置的溫度響應曲線如下圖所示圖中三條曲線分別是x=0.1, x=0.2, x=().

5、4位置的溫度瞬態響應。與圖3的 Fo=O.48迸行對比，兩種情況下的珂值不同，值越大表明熱擴散系數Q的值越 大。從圖中可以看出熱擴散系數對于導熱的影響，F(I=O.25時，與F=0.48相比較, 各位置開始響應時所需的時間較長，而且各位置響應曲線的變化速率較小，最終 的達到的溫度也較低，說明了熱擴散系數越小，熱傳導越慢，傳遞效率越低。(4)取Fll-0.51 ,得到非穩定的數值解如圖所示O 1002003004005006007008009001000時間t圖6 Fo=0.51時不同X位直處的混度響應從圖中可以看出，對于顯示格式的離散方程，并不是所有的F。值都能得到有意義的解，必須要求F<

6、;0.5時才能得到穩定的數值解，當Fo>0.5時，會出現物理上不真實的解。附件:（matlab程序）function hcat_conducri（）n0 %一維齊次熱傳導方程%設亶輸入界面OPtionS彳空間桿長LT空間點數Z ,時間點數MT擴散系數/穩定條件的值Fo(臨界值0.5); topics維非穩態導熱：標題欄顯示IinCS= 1;%輸入行為1行dcf=vooooo7i7o.48i% 默認值輸入f=inputdlg()ptions,topic,lines,dc0i%y 入框設直L=CVal(fl);%設直輸入值N=CVal(F2);MF 郵3);JI=CV 狄 1(F4);Fo=

7、CVaKf5);%FO的值必須小于0.5,小于0.5波動%計算空間步長與時間步長h=L/N;%空間步長xl=O:h:L;x=;n=Fo*h2ai%時間步長tm=n*M;%傳導總時間tl=O：n：Im;t=t;%計算初始條件與邊界條件Ti=x.*O;%初始條件To=l+t*0j%=0的邊界條件Tc=t.*O;%x=L的邊界條件%逹立溫度矩陣TT=OnCS(N+1,M+1);T(:,I)=Ti;%第一列為初始條件T(1,:)=To;%第一行為X=O邊界條件T(N+l,:)=Tc;%>后一行為X=L邊界條件%利用差分法求解溫度矩陣Tfor k=l:Mm=2;WhiIC m<=N;T(mR

8、+l)=Fo%T(m+ljk)+T(m")2*T(mrk)+T(myk);m=m+l;CndCIId%將時間空間的一維坐標轉化為二維坐標PiX=mcshgrid(tl ,);%根據溫度矩陣T繪圖SUbPIot(2,2,1);mcsh(xy,T)；%三維圖繪制(1,-1,1);%調整視圖角度 山IUC非穩態導熱)；圖像若稱 MlbUlC長度) 軸合稱 ylabclf時間tt)%y軸君稱 ZIabCIC溫度Tt)；%Z軸各稱 SUbPlot(2,2,2);A=T(II取矩陣第11列的值 PlOt(A>*);%二維曲線繪制 kgend(Q=(Ilr);%顯示函數名稱 UdCCX=0.1瞬態響應*)； XlabCIf時間 t*);ybbul(t溫度 T); ais(0 1000 0 1)；%坐標軸敵值范圍 SUbPIOt(2,2t3);B=T取矩陣第21列PIot(B,*k); ICgend(IB=O.2t);UtlCCX=0.2瞬態響應)； Xlabel(I時間 tf);Vlabclf 溫度 T)； ais(O IoOOO 1);SUbPIot(2,2,4);C=T(41,:);%取矩陣第41列 Plot(A