一維非穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解

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

2、得到總的離散方程為:THT； T：-2T；+T：rx22計(jì)算空間和時(shí)間步長(zhǎng)取空間步長(zhǎng)為:h=LN網(wǎng)格FOUriCr數(shù)為:Fo=(小于o5時(shí)穩(wěn)定)時(shí)間步長(zhǎng)為：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轉(zhuǎn)化為相應(yīng)的溫度矩陣形式：Tm.k + l) = FQ T(m + l,k) + T(m-ik)-2T(m,k) + T(myk)5. 輸入界面考慮到方程的變量，采用inpu

3、tdig函數(shù)設(shè)置5個(gè)輸入變量，對(duì)這5個(gè)變量設(shè)置了默認(rèn)值，如圖1所示。在計(jì)算中可以改變不同的數(shù)值，得到不同的結(jié)果，特別注 意穩(wěn)定條件的臨界值是().5。根據(jù)設(shè)置的默認(rèn)值，得到的計(jì)算結(jié)果如圖2所示。圖1 matbb變量輸入界面圖2默認(rèn)值的計(jì)算結(jié)果6結(jié)果分析根據(jù)上面的分析，給出了程序的輸入界面，以及默認(rèn)值狀態(tài)下的數(shù)值解?？梢酝ㄟ^改變不同的輸入值，得到需要的分析結(jié)杲，總結(jié)出了下面4點(diǎn)結(jié)論：(I)取FR).48,得到一維非穩(wěn)態(tài)導(dǎo)熱結(jié)果如下圖所示圖2 Fll=O.48時(shí)一維非穩(wěn)態(tài)導(dǎo)熱從圖中可以看出，對(duì)于長(zhǎng)度L=I的細(xì)桿，初始時(shí)刻尸0時(shí)溫度為0,邊界條 件x=0時(shí)，T=I,邊界條件x=l時(shí)，T=OO隨著時(shí)間

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

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

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

7、CVaKf5);%FO的值必須小于0.5,小于0.5波動(dòng)%計(jì)算空間步長(zhǎng)與時(shí)間步長(zhǎng)h=L/N;%空間步長(zhǎng)xl=O:h:L;x=;n=Fo*h2ai%時(shí)間步長(zhǎng)tm=n*M;%傳導(dǎo)總時(shí)間tl=O：n：Im;t=t;%計(jì)算初始條件與邊界條件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%將時(shí)間空間的一維坐標(biāo)轉(zhuǎn)化為二維坐標(biāo)PiX=mcshgrid(tl ,);%根據(jù)溫度矩陣T繪圖SUbPIot(2,2,1);mcsh(xy,T)；%三維圖繪制(1,-1,1);%調(diào)整視圖角度 山IUC非穩(wěn)態(tài)導(dǎo)熱)；圖像若稱 MlbUlC長(zhǎng)度) 軸合稱 ylabclf時(shí)間tt)%y軸君稱 ZIabCIC溫度Tt)；%Z軸各稱 SUbPlot(2,2,2);A=T(II取矩陣第11列的值 PlOt(A>*);%二維曲線繪制 kgend(Q=(Ilr);%顯示函數(shù)名稱 UdCCX=0.1瞬態(tài)響應(yīng)*)； XlabCIf時(shí)間 t*);ybbul(t溫度 T); ais(0 1000 0 1)；%坐標(biāo)軸敵值范圍 SUbPIOt(2,2t3);B=T取矩陣第21列PIot(B,*k); ICgend(IB=O.2t);UtlCCX=0.2瞬態(tài)響應(yīng))； Xlabel(I時(shí)間 tf);Vlabclf 溫度 T)； ais(O IoOOO 1);SUbPIot(2,2,4);C=T(41,:);%取矩陣第41列 Plot(A