傳熱學課程上機實習

傳熱學計算機實習指導(dǎo)書本指導(dǎo)書是為配合本科生傳熱學課中計算機應(yīng)用方面的教學而編寫的。應(yīng)用計算機解決工程實際問題，是現(xiàn)代工程技術(shù)人員所必備的技能。在傳熱學課程中引入計算機實習的目的，是使學生初步掌握用計算機求解傳熱問題的技能，從而提高學生應(yīng)用計算機解決工程實際問題的能力。大量的傳熱問題能夠用計算機求解。研究如何用計算機求解傳熱問題的專門知識數(shù)值傳熱學（或稱計算傳熱學）已經(jīng)發(fā)展成了傳熱學的一個分支學科。傳熱學課中所涉及的只是數(shù)值傳熱學的初步知識。因此，本次計算機實習也僅僅是作為數(shù)值傳熱學的入門。本指導(dǎo)書給出了三個練習題及相應(yīng)的算法。這三個練習題分別涉及了一維穩(wěn)態(tài)導(dǎo)熱、二維穩(wěn)態(tài)導(dǎo)熱和一維非穩(wěn)態(tài)導(dǎo)熱。要求學生在掌握問題的數(shù)值計算方法的基礎(chǔ)上，獨立編寫計算機程序并用所編的程序計算出這三個練習題的數(shù)值結(jié)果。1練習題一：一維穩(wěn)態(tài)導(dǎo)熱的數(shù)值計算1物理問題圖1示出了一個等截面直肋，處于溫度tK80°C的流體中。肋表面與流休之間的對流換熱系數(shù)為h=45W/m2.C，肋基處溫度tw=300C，肋端絕熱。肋片由鋁合金制成，其導(dǎo)熱系數(shù)為入=110W/mC，肋片厚度為0=0.01m，高度為H=0.1m。試計算肋內(nèi)的溫度分布及肋的總換熱量。i-1i AlJH-l_J ?-圖1等截面直肋傳熱示意圖2數(shù)學描述及其解析解t—t引入無量綱過余溫度。件為：= -引入無量綱過余溫度。件為：(1-1)(1-2)(1-3)迎一m2。=0dx(1-1)(1-2)(1-3)x=0,0=18。x=H,—=08x

■hp其中m=1(其中符號含義與教科書楊世銘陶文銓編著《傳熱學》相同，以下可。人At^=▼t^=▼^^ch\m(x-H)]8w8~ch(mH)-(1-4)卜匝—t)h(mH)按式(1-4)計算得到的在肋內(nèi)各點的溫度由表1給出。表1等截面直肋內(nèi)各點的溫度坐標(X)00.010.020.030.040.050.060.070.080.090.1溫度(T)300.00286.56274.84264.76256.13248.97243.19238.1235.58233.70233.091.3數(shù)值離散1.3.1 區(qū)域離散在對方程(1-1)?(1-3)進行數(shù)值離散之前,應(yīng)首先進行計算區(qū)域的離散。計算區(qū)域的離散如圖1所示，總節(jié)點數(shù)取N。1.3.2 微分方程的離散由于方程(1-1)在計算區(qū)域內(nèi)部處處成立，因而對圖1所示的各離散點亦成立。對任一節(jié)點i有：—m26=0(d26)—m26=0[dx2J

i用。在節(jié)點i的二階差分代替。在節(jié)點i的二階導(dǎo)數(shù)，得:氣1-老+11-m26Ax2 i整理上式成迭代形式:(1-5)(i=2,3,…,N-1)(1-5)1.3.3 邊界條件離散上面得到的離散方程式(1-5),對所有內(nèi)部節(jié)點都成立，因此每個內(nèi)部節(jié)點都可得出一個類似的方程。事實上，式(1-5)表達的是一個代數(shù)方程組。但這個方程組的個數(shù)少于未知數(shù) 6i(i=1,2,……,N)的個數(shù)。因此，還需要根據(jù)邊界條件補充進兩個方程后代數(shù)方程組才封閉。左邊界(x=0)為第一類邊界條件，溫度為已知，因此可以根據(jù)式(1-2)直接補充一個方程為：I"'8M廣1右邊界為第三類邊界條件，由圖1中邊界節(jié)點N的向后差分來代替式(1-3)中的導(dǎo)數(shù)，得：6-6N N-1=0Ax

將此式整理為迭代形式，得:1.3.4 最終的離散格式6廣6=1Z=27打"s+"1)SW) (1-6)6=6N N-11.3.5 代數(shù)方程組的求解及其程序應(yīng)首先假定一個溫度場的初始分布，即給出代數(shù)方程組有各種求解方法，較為有效而簡便的方法是高斯-賽德爾迭代方法。式(1-6)已給出了代數(shù)方程組的迭代形式。在實際計算中各節(jié)點的溫度初值：應(yīng)首先假定一個溫度場的初始分布，即給出60,6060,60,...,601 2N(1-6)中進行迭代計算，：

則K+1次迭代的計算式為:將這些初值代入方程組6f,62,...,6k為已知，6k+1=6叩+1=irmkr',+1K+6-1K+16K+1=6K+1(60=6=1)直至收斂。假如第K步迭代已完成，即(i=2,3,…,N-1)(1-7)根據(jù)式(1-7)編寫程序的工1作由學生自行完成。計算結(jié)果可與解析解比較。2練習題二：二維穩(wěn)態(tài)導(dǎo)熱的數(shù)值計算1物理問題圖2示出了一矩形區(qū)域，其邊長L=W=1，假設(shè)區(qū)域內(nèi)無內(nèi)熱源，導(dǎo)熱系數(shù)為常數(shù)，三個邊溫度為T]=0，一個邊溫度為T2=1，求該矩形區(qū)域內(nèi)的溫度分布。L副二^形區(qū)睇態(tài)導(dǎo).熱

L副二^形區(qū)睇態(tài)導(dǎo).熱2.2數(shù)學描述對上述問題的微分方程及其邊界條件為:ST ST八+ =0dST ST八+ =0dx2 dy2x=0,T=T=0ix=1,T=T=0iy=0,T=T1=0y=1,T=T=1作為參考，以下給出該問題的解析解： 2——1-T2題1-(-1)n.(n兀一工 sin——兀n\Ln=11(皿、sh-ykL7s^竺-W'kL7(2-1)(2-2)(2-3)表2列出了由式(2-3)計算得到的在平面區(qū)域內(nèi)各不同位置的溫度值。表2平面區(qū)域內(nèi)溫度分布坐標溫度T坐標溫度Txyxy0.250.250.0680.250.750.4320.500.250.0950.500.750.5400.750.250.0680.750.750.4321.000.250.4x10-61.000.753x10-60.250.500.1820.251.000.9990.500.500.250.501.001.0240.750.500.1820.751.000.9991.000.502x10-61.001.0065x10-63數(shù)值離散區(qū)域離散區(qū)域離散如圖2所示，x方向總節(jié)點數(shù)為N,y方向總節(jié)點數(shù)為M,區(qū)域內(nèi)任一節(jié)點用i,j表示。方程的離散對于圖2中所有的內(nèi)部節(jié)點方程(2-1)都適用，因此可寫為：‘仞t) (d21)八〔相7+〔程7i,j i,j用i,j節(jié)點的二階中心差分代替上式中的二階導(dǎo)數(shù)，得：

T*j—T*j—2T*+TAx2山+T，LT，+T，頊=0Ay2上式整理成迭代形式：T= ——^T+T)+ ——^T+T)ij2'Ax2+Ay2i+i,jtj2'Ax2+Ay2/j+1ij-1(i=2,3,…,N-1),(j=2,3,…,M-1)補充四個邊界上的第一類邊界條件得:(j=1,2,…,M)(j=1,2,…,M)(i=1,2,…(j=1,2,…,M)(j=1,2,…,M)(i=1,2,…,N)(i=1,2,…,N)T\=TT"=TT5i,M2.4 計算程序計算程序由學生自行完成。3練習題三：一維非穩(wěn)態(tài)導(dǎo)熱的數(shù)值計算非穩(wěn)態(tài)導(dǎo)熱問題由于有時間變量，其數(shù)值計算出現(xiàn)了一些新的特點。在非穩(wěn)態(tài)導(dǎo)熱微分方程中，與時間因素相關(guān)的非穩(wěn)態(tài)項是溫度對時間的一階導(dǎo)數(shù)，這給差分離散帶來了新的特點。由于這個特點，可以采用不同的方法構(gòu)造差分方程，從而得到幾種不同的差分格式，即所謂的顯式、隱式和半隱式。我們?nèi)匀粡囊粋€具體問題出發(fā)來研究非穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值計算。3.1問題一塊無限大平板（如圖3所示），其一半厚度為L=0.1m，初始溫度To=1000°C，突然將其圖3平板一細曜態(tài)皇熱插入溫度L=20C的流體介質(zhì)中。平板的導(dǎo)熱系數(shù)入=34.89W/mC，密度p=7800kg/m3，比熱c=0.712x10；J/kgC，平板與介質(zhì)的對流換熱系數(shù)為h=233W/m2C，求平板內(nèi)各點的溫度分布。2 數(shù)學描述由于平板換熱關(guān)于中心線是對稱的，僅對平板一半?yún)^(qū)域進行計算即可。坐標x的原點選在平板中心線上，因而一半?yún)^(qū)域的非穩(wěn)態(tài)導(dǎo)熱的數(shù)學描述為：TOC\o"1-5"\h\z6TST (3-1)—=a Stdx2\o"CurrentDocument"T=0,T=T (3-2)0stx=0,絲=0 (3-3)dx\o"CurrentDocument"x=L,-xST.=h(T-T) (3-4)該數(shù)學模型的解析解為：\o"CurrentDocument"T=T+T-T般——2^^^——cos[^x'-隊氣 (3-5)8 0 8 1|li+sin日cos日knL)\o"CurrentDocument"aT hL其中F0=~l2H為方程ctg^=旦/B.的根，B.—。表3給出了在平板表面(x=L)處由式(3-5)計算得到的不同時刻的溫度值。表3 平板表面各不同時刻溫度值。時間(S)12345678910溫度(°C)981.84974.47968.88964.20960.11956.14953.08949.97947.07944.343數(shù)值離散3.3.1 計算區(qū)域的離散一維非穩(wěn)態(tài)導(dǎo)熱指的是空間坐標是一維的。若考慮時間坐標，則所謂的一維非穩(wěn)態(tài)導(dǎo)熱實圖4X-T平面區(qū)域的離散

際上是二維問題(見圖4),即：有時間坐標T和空間坐標x兩個變量。但要注意，時間坐標是單向的，就是說，前一時刻的狀態(tài)會對后一時刻的狀態(tài)有影響，但后一時刻的狀態(tài)卻影響不到前一時刻，圖4示出了以x和T為坐標的計算區(qū)域的離散，時間從！=0開始，經(jīng)過一個個時層增加到K時層和K+1時層。3.3.2 微分方程的離散對于i節(jié)點，在K和K+1時刻可將微分方程(3-1瀉成下面式子:dT)K+1 (a2T=aSt) ^3x2dT)K+1 (a2T=aSt) ^3x2)ii將式(3-6)?(3-8)的左端溫度對時間的偏導(dǎo)數(shù)進行差分離散為:3T]k Tk+1-tk I —i iSt) Ati觀察式(3-8)和(3-9)，(3-6)(3-7)(3-8)(3-9)這兩個式子的右端差分式完全相同，但在兩個式子中卻有不同含義。對式 (dTV (3-8)，右瑞項相對i點在K時刻的導(dǎo)數(shù)—是向前差分。而在式(3-9)中，右端項是I點在Idt)i(dT\K+1K+1時刻的導(dǎo)數(shù)— 的向后差分。將式(3-8)和(3-9)分別代入式(3-6)和(3-7)，并將式(3-6)和ldT)i(3-7)右端關(guān)于x的二階導(dǎo)數(shù)用相應(yīng)的差分代替，則可得到下列顯式和隱式兩種不同的差分格式:顯式：Tk+1=fTK+(1-2f)Tk+fTKi i+1 i i-1(K=0,1,2, ,i=2,3,…,N-1)(3-10)全隱式:Tk+1=―1—^fTK+1+fTK+1+Tk)

i1+2f i+1 i-1 i(3-11)(K=0,1,2, i=2,3,…,N-1)以上兩式中的fp。從式(3-10)可見，其右端只涉及K時刻的溫度，當從K=0(即T=0時刻)開始計算時，在K=0時等號右端都是已知值，因而直接可計算出K=1時刻各點的溫度。由K=1時刻的各點的溫度值，又可以直接利用式(3-10)計算K=2時刻的各點的溫度，這樣一個時層一個時層的往下推，各時層的溫度都能用式(3-10)直接計算出來，不要求解代數(shù)方程組。而對于式(3-11)等號右端包含了與等號左端同一時刻但不同節(jié)點的溫度，因而必須通過求解代數(shù)方程組才能求得這些節(jié)點的溫度值。3.3.3邊界條件的離散對于式(3-3)和(3-4)所給出的邊界條件，可以直接用差分代替微分，也可以用元體平衡法給出相應(yīng)的邊界條件，亦有顯式和隱式之分。通常，當內(nèi)部節(jié)點采用顯式時，邊界節(jié)點也用顯式離散；內(nèi)部節(jié)點用隱式時，邊界節(jié)點亦用隱式。邊界節(jié)點的差分格式是顯示還是隱式，取決于如何與內(nèi)部節(jié)點的差分方程組合。用K+1時刻相應(yīng)節(jié)點的差分，代替式(3-3)和(3-4)中的微分，可得到邊界節(jié)點的差分方程：TK+1=TK+1=TK+1

1TK+1N21hAx.—+1(—Tk+1+(3-12)3.3.4 最終的離散格式顯式：初始值： T=T (i=1,2,3,…,N) (3-13)TK+1=^fTk+fTK+G-2fM) (i=2,3,…,N-1) (3-14)i i+1 i-1 iTK+1'=TK+1Tk+1= —fTk+1+竺T) (3-15)N hAx1"N-1 入勺\其中K=0,1,2,??…隱式