穩(wěn)態(tài)熱傳導(dǎo)問(wèn)題有限元法

1、6.穩(wěn)態(tài)熱傳導(dǎo)問(wèn)題的有限元法本章的內(nèi)容如下：6.1熱傳導(dǎo)方程與換熱邊界6.2穩(wěn)態(tài)溫度場(chǎng)分析的一般有限元列式6.3三角形單元的有限元列式6.4溫度場(chǎng)分析舉例6.1熱傳導(dǎo)方程與換熱邊界在分析工程問(wèn)題時(shí)， 經(jīng)常要了解工件內(nèi)部的溫度分布情況，例如發(fā)動(dòng)機(jī)的工作溫度、 金屬工件在熱處理過(guò)程中的溫度變化、流體溫度分布等。物體內(nèi)部的溫度分布取決于物體內(nèi)部的熱量交換，以及物體與外部介質(zhì)之間的熱量交換，一般認(rèn)為是與時(shí)間相關(guān)的。物體內(nèi)部的熱交換采用以下的熱傳導(dǎo)方程（Fourier方程）來(lái)描述，TTTTcxyzQ（ 6-1）t x x y y z z式中 為密度，kg/m3； c為比熱容，J/（kg K） ； x,

2、y, z為導(dǎo)熱系數(shù)， w m k ； T 為溫度，C； t為時(shí)間，s； Q為內(nèi)熱源密度，w/m3。對(duì)于各向同性材料，不同方向上的導(dǎo)熱系數(shù)相同，熱傳導(dǎo)方程可寫(xiě)為以下形式，T2T2T2Tc 22 Q（6-2）txyz除了熱傳導(dǎo)方程，計(jì)算物體內(nèi)部的溫度分布，還需要指定初始條件和邊界條件。初始條件是指物體最初的溫度分布情況，T t 0 T0 x, y,z（6-3）邊界條件是指物體外表面與周?chē)h(huán)境的熱交換情況。在傳熱學(xué)中一般把邊界條件分為三類。1）給定物體邊界上的溫度，稱為第一類邊界條件。物體表面上的溫度或溫度函數(shù)為已知，Ts Ts或 Ts Ts（x,y,z,t）（6-4）2）給定物體邊界上的熱量輸入或

3、輸出，稱為第二類邊界條件。 已知物體表面上熱流密度，(xT nx xTynyyTznz)zsqsTTT 、或(xnxynyznz)sqs(x, y, z,t)（6-5）xyz3) 給定對(duì)流換熱條件，稱為第三類邊界條件。物體與其相接觸的流體介質(zhì)之間的對(duì)流換熱系數(shù)和介質(zhì)的溫度為已知。Tx nxxTynyyz 丄 nzh(Tf Ts)z(6-6)其中h為換熱系數(shù),W/(m 2 K);Ts是物體表面的溫度；Tf是介質(zhì)溫度。如果邊界上的換熱條件不隨時(shí)間變化，物體內(nèi)部的熱源也不隨時(shí)間變化，在經(jīng)過(guò)一定時(shí)間的熱交換后，物體內(nèi)各點(diǎn)溫度也將不隨時(shí)間變化，即丄0 t的變化，而是指溫度分布穩(wěn)定后的狀態(tài)， 到最后的穩(wěn)定

4、溫度場(chǎng)。 隨時(shí)間變化的瞬態(tài) 三維問(wèn)題的穩(wěn)態(tài)熱傳導(dǎo)方程為，這類問(wèn)題稱為穩(wěn)態(tài)(Steady state)熱傳導(dǎo)問(wèn)題。穩(wěn)態(tài)熱傳導(dǎo)問(wèn)題并不是溫度場(chǎng)不隨時(shí)間 我們不關(guān)心物體內(nèi)部的溫度場(chǎng)如何從初始狀態(tài)過(guò)渡熱傳導(dǎo)方程就退化為穩(wěn)態(tài)熱傳導(dǎo)方程,(Transient)(6-7)對(duì)于各向同性的材料,222TTT222xyz可以得到以下的方程，稱為(6-8)Poisson 方程,Laplace 方程,考慮物體不包含內(nèi)熱源的情況，各向同性材料中的溫度場(chǎng)滿足(6-9)2t2t2t02220x y z在分析穩(wěn)態(tài)熱傳導(dǎo)問(wèn)題時(shí)，不需要考慮物體的初始溫度分布對(duì)最后的穩(wěn)定溫度場(chǎng)的影 響，因此不必考慮溫度場(chǎng)的初始條件，而只需考慮換熱

5、邊界條件。計(jì)算穩(wěn)態(tài)溫度場(chǎng)實(shí)際上是求解偏微分方程的邊值問(wèn)題。溫度場(chǎng)是標(biāo)量場(chǎng)，將物體離散成有限單元后，每個(gè)單元結(jié)點(diǎn)上 只有一個(gè)溫度未知數(shù)，比彈性力學(xué)問(wèn)題要簡(jiǎn)單。進(jìn)行溫度場(chǎng)計(jì)算時(shí)有限單元的形函數(shù)與彈性 力學(xué)問(wèn)題計(jì)算時(shí)的完全一致，單元內(nèi)部的溫度分布用單元的形函數(shù)，由單元結(jié)點(diǎn)上的溫度來(lái)確定。由于實(shí)際工程問(wèn)題中的換熱邊界條件比較復(fù)雜，在許多場(chǎng)合下也很難進(jìn)行測(cè)量，如何定義正確的換熱邊界條件是溫度場(chǎng)計(jì)算的一個(gè)難點(diǎn)。6.2穩(wěn)態(tài)溫度場(chǎng)分析的一般有限元列式在前面我們已經(jīng)介紹了有限元方法可以用來(lái)分析場(chǎng)問(wèn)題，穩(wěn)態(tài)溫度場(chǎng)計(jì)算是一個(gè)典型的場(chǎng)問(wèn)題。我們可以采用虛功方程建立彈性力學(xué)問(wèn)題分析的有限元格式，推導(dǎo)出的單元?jiǎng)偠染仃囉忻?/p>

6、確的力學(xué)含義。在這里，介紹如何用加權(quán)余量法(Weighted Residual Method )建立穩(wěn)態(tài)溫度場(chǎng)分析的有限元列式。微分方程的邊值問(wèn)題，可以一般地表示為未知函數(shù)u滿足微分方程組，未知函數(shù)u還滿足邊界條件，B(u)B1（u）B2（u）0（在邊界上）(6-11)A(u)A(u)A2(u)0（在域 內(nèi)）(6-10)如果未知函數(shù) u是上述邊值問(wèn)題的精確解，則在域中的任一點(diǎn)上u都滿足微分方程（6-10），在邊界的任一點(diǎn)上都滿足邊界條件（6-11 ）。對(duì)于復(fù)雜的工程問(wèn)題，這樣的精確解往往很難找到，需要設(shè)法尋找近似解。所選取的近似解是一族帶有待定參數(shù)的已知函數(shù)， 般表示為n(6-12)u u N

7、 iai Nai 1其中ai為待定系數(shù)，Ni為已知函數(shù)，被稱為試探函數(shù)。試探函數(shù)要取自完全的函數(shù)序列， 是線性獨(dú)立的。由于試探函數(shù)是完全的函數(shù)序列，任一函數(shù)都可以用這個(gè)序列來(lái)表示。采用這種形式的近似解不能精確地滿足微分方程和邊界條件，所產(chǎn)生的誤差就稱為余量。微分方程（6-10）的余量為，R A（Na）（ 6-13）邊界條件（6-11）的余量為，R B（Na）（ 6-14）選擇一族已知的函數(shù)， 使余量的加權(quán)積分為零， 強(qiáng)迫近似解所產(chǎn)生的余量在某種平均意 義上等于零，TT Wj RdWj Rd 0（ 6-15）Wj和Wj稱為權(quán)函數(shù)，通過(guò)公式（6-15）可以選擇待定的參數(shù)ai。這種采用使余量的加權(quán)積

8、分為零來(lái)求得微分方程近似解的方法稱為加權(quán)余量法。對(duì)權(quán)函數(shù)的不同選擇就得到了不同的加權(quán)余量法，常用的方法包括配點(diǎn)法、子域法、最小二乘法、 力矩法和伽遼金法（Galerkin method ）。在很多情況下，采用Galerkin法得到的方程組的系數(shù)矩陣是對(duì)稱的，在這里也采用Galerkin法建立穩(wěn)態(tài)溫度場(chǎng)分析的一般有限元列式。在Galerkin法中，直接采用試探函數(shù)序列作為權(quán)函數(shù)，取WjN j，WjNj。下面用求解二階常微分方程為例，說(shuō)明Galerkin法（參見(jiàn)，王勖成編著“有限元法基本原理和數(shù)值方法”的節(jié)）。d2udx2例，求解二階常微分方程(0 x 1)邊界條件：當(dāng)x 0時(shí)，u 0 ;當(dāng)x 1

9、時(shí)，u 0。 取兩項(xiàng)近似解：Ni x(1 x)2N2 x (1 x)u N1a1 N2a2a1x(1 x)a2x (1 x)W1N1 ,W2N2由公式(6-15)可以得到兩個(gè)加權(quán)積分方程,x(12x)x 印(2 x x ) a2(2 6x23、x x )dx 0x2(12x)x a1( 2 x x )23a2(2 6x x x )dx 0積分后可以得到一個(gè)二元一次方程組，解得,a1 0.1924,a20.1707近似解為，x(1x)( 0.1924 0.1707x)sin x該方程的精確解為，uxsi n1近似解與精確解的結(jié)果比較見(jiàn)表6-1，表6-1近似解與精確解比較x=0.25x=0.5x=

10、0.75sin x uxsin 10.044010.069750.06006x(1x)(0.19240.1707x)0.044080.069440.06008假定單元的形函數(shù)為，N N12 心單元結(jié)點(diǎn)的溫度為，Te T1 T2. TnT單元內(nèi)部的溫度分布為,T NTe.維問(wèn)以二維問(wèn)題為例，說(shuō)明用Galerkin法建立穩(wěn)態(tài)溫度場(chǎng)的一般有限元格式的過(guò)程。題的穩(wěn)態(tài)熱傳導(dǎo)方程為,(6-16a)第一類換熱邊界為T(mén)ss第二類換熱邊界條件為,(6-16b)TTx n xy n y q sxy第三類邊界條件為，(6-16c)T nx xTnyh(TfTs)y(6-16d)eTTW1 -(x) (y)Qd0xx

11、 yy由分部積分得，T、W1,T、T、(w1x)：(x-)W1( x )xxxxxxT、W1 ,T、T、(w1y)1( y)W1( y)yyyyyy在一個(gè)單元內(nèi)的加權(quán)積分公式為,(6-17)應(yīng)用Green定理，一個(gè)單元內(nèi)的加權(quán)積分公式寫(xiě)為,x(-(yT nx x(6-18)Tn y)d0y采用Galerkin方法，選擇權(quán)函數(shù)為,W!Ni將單元內(nèi)的溫度分布函數(shù)和換熱邊界條件代入(6-18)式，單元的加權(quán)積分公式為,6 旦(xxN)旦(x yy3)Ted yNiQde2Mqsd(6-19)MhNTedMhTfd 0換熱邊界條件代入后,在(6-19)式內(nèi)相應(yīng)出現(xiàn)了第二類換熱邊界項(xiàng)e3 Niqsd ，

12、第三3類換熱邊界項(xiàng)Nih NTed3e3 NihTfd ，但沒(méi)有出現(xiàn)與第一類換熱邊界對(duì)應(yīng)的3項(xiàng)。這是因?yàn)?，采?M作為權(quán)函數(shù)，第一類換熱邊界被自動(dòng)滿足。寫(xiě)成矩陣形式有,e (型)T(Xe tNTQd(6-20)型)(叫x ye t2NTqsd)Ted yTehN NT d3 NThTfd0n個(gè)結(jié)點(diǎn)的溫度Ti。按有限元格式將公式（6-20 ）是n個(gè)聯(lián)立的線性方程組，可以確定（6-20 ）表示為，KeTe Pe(6-21)其中矩陣Ke為單元的導(dǎo)熱矩陣或稱為溫度剛度矩陣, 稱為單元的溫度載荷向量或熱載荷向量（ 熱矩陣Ke和溫度載荷向量KjNi N-j)dye3 hNiNjd(6-22)MqsdNih

13、Tf d2如果某個(gè)單元完全處于物體的內(nèi)部,NiQd(6-23)Kj(NiNj(xx xNiyyT e為單元的結(jié)點(diǎn)溫度向量，PeThermal load vector）。對(duì)于某個(gè)特定單元，單元導(dǎo)P e的元素分別為，NiQd在整個(gè)物體上的加權(quán)積分方程是單元積分方程的和,宀T(X凹）NTQde t2NTqsd(6-24)3 hNTNTed3 NThTfd0e根據(jù)單元結(jié)點(diǎn)的局部編號(hào)與整體編號(hào)的關(guān)系，直接求和得到整體剛度矩陣，整體方程組 為,KTP6.3三角形單元的有限元列式(心亠)圖6-1三角形單元回顧第三章的內(nèi)容可以發(fā)現(xiàn)，與計(jì)算彈性力學(xué)平面問(wèn)題時(shí)所采用的方法一樣,場(chǎng)問(wèn)題計(jì)算中所采用的三角形單元可以使

14、用相同的形函數(shù)，壬(ai bi x Ci y)2 A±(aj bjX Cjy) 2A(am bmX Cmy)2A二維溫度NiNjNmaiajXj ymXmyjbiyjymCiXmXjXmYiXiYmbjymyiCjXiXmXiyjXj yibmyiyjCmXjXiXiyiXjyjTXmym2A單元內(nèi)的溫度分布用結(jié)點(diǎn)上的溫度值表示為,T Ni NjNm TjTm(6-25)在三角形單兀上,采用Galerki n法可得，A NTL (AXx T)X(y _ yT)QdA 0 y(6-26)Ti-)XNt t假定單元內(nèi)的導(dǎo)熱系數(shù)為常數(shù),NT TXdAX X14Abi bj bi bmbjb

15、mTi Tj TmdA4Ab2bibjbibmbi bjb2bjbmbibmbjbmbmTiTjTm(6-27)Nt-)y(N y-)yNt tyNT TydA y14ACCj CicmCjCmTiTjTmdA4A CiCjCi CmGC2CjCj CmCi CmC j Cm2CmTiTjTm(6-28)單元的剛度矩陣為,Ke/4bi2bibjbi bj2CiCiCjCi Cmbibjbj2bjbmyGCj2CjCj Cm c；bibmbjbmb；4Ci CmCj CmXNi1T Z-QaQANTQdAAN j dA13Nm1由Green公式可得A (NTXT)(NT y TdAAXXyy-

16、(NT xT nNT y-Tn y)dSsXy顯然，單元的導(dǎo)熱矩陣是對(duì)稱的。如果單元的內(nèi)部熱源為常數(shù)，由內(nèi)部熱源產(chǎn)生的溫度載荷項(xiàng)為,方便起見(jiàn)，把換熱邊界統(tǒng)一表示為第三類換熱邊界，(6-29)(6-30)A (NT X 丄)一(NT y 丄dAAxx yy( 6-30)TTTe：hN (Tf Ts)dS ； hN TfdS hN NT dS如果在單元邊上存在熱交換，各條邊上的邊界換熱條件在單元?jiǎng)偠染仃囍猩傻母郊禹?xiàng)為，hlj210Ke6120(6-31)000000Kehl jm021(6-32)6012201Kehl mi000(6-33)6102由邊界換熱條件生成的溫度載荷向量為,1hTfl

17、jPe21(6-34)00PehT f l jm21(6-35)11ehT f l miP20(6-36)16.4溫度場(chǎng)分析舉例正方形截面的煙囪如圖 6-2所示，煙囪由混凝土建造，邊長(zhǎng)為60cm，通道的邊長(zhǎng)為20cm,混凝土的導(dǎo)熱系數(shù)為 k 1.4W /(m K)。假定煙囪內(nèi)表面的溫度為 100°C,煙囪外表面暴露在空氣中，空氣的溫度為30C,換熱系數(shù)為h20W/(m2 K)。計(jì)算煙囪截面內(nèi)的穩(wěn)態(tài)溫度場(chǎng)。(參見(jiàn)，F(xiàn)inite Element Method Theory and Application with ANSYS, p279)60cjl圖6-2煙囪截面圖6-3有限元模型G冷0863胡勻 11州"擁圖6-4穩(wěn)態(tài)溫度分布Ml92.617B5.3r?4ICld圖6-5熱流量分布穩(wěn)態(tài)溫度場(chǎng)分布與物體的初始狀態(tài)無(wú)關(guān)，那么是否與材料的導(dǎo)熱系數(shù)相關(guān)？我們把煙囪的模型做些修改，假定煙囪壁由兩層材料構(gòu)成。內(nèi)層材料為混凝土，外表面的截面尺寸為 30cm 30cm，煙囪通道 的尺寸不 變，仍為20cm 20cm。外 層材料的 導(dǎo)熱 系數(shù)為 k 0.1W/(m K)，外部表面的截面尺寸不變，內(nèi)部表面的截面尺寸為30cm 30cm。換熱邊界條件不變，雙層煙囪的有限元模型如圖6-6所示。圖6-6雙層煙囪的有限元模型JTlIiAL gCLEMH 3TTP-J.T