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

2、y, z為導熱系數， w m k ； T 為溫度，C； t為時間，s； Q為內熱源密度，w/m3。對于各向同性材料，不同方向上的導熱系數相同，熱傳導方程可寫為以下形式，T2T2T2Tc 22 Q（6-2）txyz除了熱傳導方程，計算物體內部的溫度分布，還需要指定初始條件和邊界條件。初始條件是指物體最初的溫度分布情況，T t 0 T0 x, y,z（6-3）邊界條件是指物體外表面與周圍環境的熱交換情況。在傳熱學中一般把邊界條件分為三類。1）給定物體邊界上的溫度，稱為第一類邊界條件。物體表面上的溫度或溫度函數為已知，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) 給定對流換熱條件，稱為第三類邊界條件。物體與其相接觸的流體介質之間的對流換熱系數和介質的溫度為已知。Tx nxxTynyyz 丄 nzh(Tf Ts)z(6-6)其中h為換熱系數,W/(m 2 K);Ts是物體表面的溫度；Tf是介質溫度。如果邊界上的換熱條件不隨時間變化，物體內部的熱源也不隨時間變化，在經過一定時間的熱交換后，物體內各點溫度也將不隨時間變化，即丄0 t的變化，而是指溫度分布穩定后的狀態， 到最后的穩定

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

5、邊界條件。計算穩態溫度場實際上是求解偏微分方程的邊值問題。溫度場是標量場，將物體離散成有限單元后，每個單元結點上 只有一個溫度未知數，比彈性力學問題要簡單。進行溫度場計算時有限單元的形函數與彈性 力學問題計算時的完全一致，單元內部的溫度分布用單元的形函數，由單元結點上的溫度來確定。由于實際工程問題中的換熱邊界條件比較復雜，在許多場合下也很難進行測量，如何定義正確的換熱邊界條件是溫度場計算的一個難點。6.2穩態溫度場分析的一般有限元列式在前面我們已經介紹了有限元方法可以用來分析場問題，穩態溫度場計算是一個典型的場問題。我們可以采用虛功方程建立彈性力學問題分析的有限元格式，推導出的單元剛度矩陣有明

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

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

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

9、時，u 0。 取兩項近似解：Ni x(1 x)2N2 x (1 x)u N1a1 N2a2a1x(1 x)a2x (1 x)W1N1 ,W2N2由公式(6-15)可以得到兩個加權積分方程,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積分后可以得到一個二元一次方程組，解得,a1 0.1924,a20.1707近似解為，x(1x)( 0.1924 0.1707x)sin x該方程的精確解為，uxsi n1近似解與精確解的結果比較見表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假定單元的形函數為，N N12 心單元結點的溫度為，Te T1 T2. TnT單元內部的溫度分布為,T NTe.維問以二維問題為例，說明用Galerkin法建立穩態溫度場的一般有限元格式的過程。題的穩態熱傳導方程為,(6-16a)第一類換熱邊界為Tss第二類換熱邊界條件為,(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在一個單元內的加權積分公式為,(6-17)應用Green定理，一個單元內的加權積分公式寫為,x(-(yT nx x(6-18)Tn y)d0y采用Galerkin方法，選擇權函數為,W!Ni將單元內的溫度分布函數和換熱邊界條件代入(6-18)式，單元的加權積分公式為,6 旦(xxN)旦(x yy3)Ted yNiQde2Mqsd(6-19)MhNTedMhTfd 0換熱邊界條件代入后,在(6-19)式內相應出現了第二類換熱邊界項e3 Niqsd ，

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

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

14、用相同的形函數，壬(ai bi x Ci y)2 A±(aj bjX Cjy) 2A(am bmX Cmy)2A二維溫度NiNjNmaiajXj ymXmyjbiyjymCiXmXjXmYiXiYmbjymyiCjXiXmXiyjXj yibmyiyjCmXjXiXiyiXjyjTXmym2A單元內的溫度分布用結點上的溫度值表示為,T Ni NjNm TjTm(6-25)在三角形單兀上,采用Galerki n法可得，A NTL (AXx T)X(y _ yT)QdA 0 y(6-26)Ti-)XNt t假定單元內的導熱系數為常數,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顯然，單元的導熱矩陣是對稱的。如果單元的內部熱源為常數，由內部熱源產生的溫度載荷項為,方便起見，把換熱邊界統一表示為第三類換熱邊界，(6-29)(6-30)A (NT X 丄)一(NT y 丄dAAxx yy( 6-30)TTTe：hN (Tf Ts)dS ； hN TfdS hN NT dS如果在單元邊上存在熱交換，各條邊上的邊界換熱條件在單元剛度矩陣中生成的附加項為，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溫度場分析舉例正方形截面的煙囪如圖 6-2所示，煙囪由混凝土建造，邊長為60cm，通道的邊長為20cm,混凝土的導熱系數為 k 1.4W /(m K)。假定煙囪內表面的溫度為 100°C,煙囪外表面暴露在空氣中，空氣的溫度為30C,換熱系數為h20W/(m2 K)。計算煙囪截面內的穩態溫度場。(參見，Finite Element Method Theory and Application with ANSYS, p279)60cjl圖6-2煙囪截面圖6-3有限元模型G冷0863胡勻 11州"擁圖6-4穩態溫度分布Ml92.617B5.3r?4ICld圖6-5熱流量分布穩態溫度場分布與物體的初始狀態無關，那么是否與材料的導熱系數相關？我們把煙囪的模型做些修改，假定煙囪壁由兩層材料構成。內層材料為混凝土，外表面的截面尺寸為 30cm 30cm，煙囪通道 的尺寸不 變，仍為20cm 20cm。外 層材料的 導熱 系數為 k 0.1W/(m K)，外部表面的截面尺寸不變，內部表面的截面尺寸為30cm 30cm。換熱邊界條件不變，雙層煙囪的有限元模型如圖6-6所示。圖6-6雙層煙囪的有限元模型JTlIiAL gCLEMH 3TTP-J.T