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

1、2T2 y(6-2)除了熱傳導(dǎo)方程,計(jì)算物體內(nèi)部的溫度分布,還需要指定初始條件和邊界條件。初始條件是指物體最初的溫度分布情況,T t 0To x, y,z(6-3)邊界條件是指物體外表面與周圍環(huán)境的熱交換情況。 類。1)給定物體邊界上的溫度，稱為第一類邊界條件。物體表面上的溫度或溫度函數(shù)為已知，在傳熱學(xué)中一般把邊界條件分為三TsTs或 Ts Ts(x,y,z,t)(6-4)2)給定物體邊界上的熱量輸入或輸出，稱為第二類邊界條件。 已知物體表面上熱流密度，(xT nx xT ynyyTznz)zsqs或(xT nx xTynyyznz)zsqs(x, y, z,t)(6-5)6.穩(wěn)態(tài)熱傳導(dǎo)問題的

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

3、z為導(dǎo)熱系數(shù)，w/ m k ； T為溫度，C； t為時(shí)間，s; Q為內(nèi)熱源密度，w/m3。對(duì)于各向同性材料，不同方向上的導(dǎo)熱系數(shù)相同，熱傳導(dǎo)方程可寫為以下形式,2T2x3)給定對(duì)流換熱條件，稱為第三類邊界條件。物體與其相接觸的流體介質(zhì)之間的對(duì)流換熱系數(shù)和介質(zhì)的溫度為已知。Tx nxxTy nyyz nzh(Tf Ts)z(6-6)其中h為換熱系數(shù),W/(m 2 K);Ts是物體表面的溫度；Tf是介質(zhì)溫度。如果邊界上的換熱條件不隨時(shí)間變化，物體內(nèi)部的熱源也不隨時(shí)間變化，在經(jīng)過一定時(shí)間的熱交換后，物體內(nèi)各點(diǎn)溫度也將不隨時(shí)間變化，即二0t這類問題稱為穩(wěn)態(tài)(Steady state熱傳導(dǎo)問題。穩(wěn)態(tài)熱傳

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

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

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

7、為nu uNiaiNa(6-12)i 1其中ai為待定系數(shù)，Ni為已知函數(shù)，被稱為試探函數(shù)。試探函數(shù)要取自完全的函數(shù)序列，是線性獨(dú)立的。由于試探函數(shù)是完全的函數(shù)序列，任一函數(shù)都可以用這個(gè)序列來表示。采用這種形式的近似解不能精確地滿足微分方程和邊界條件，所產(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)四和皿 稱為權(quán)函數(shù)，通過公式(6-15)可以選擇待定的參數(shù) a-這種采用使余

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

9、dx邊界條件：當(dāng)x 0時(shí)，u 0;當(dāng)x 1時(shí)，u 0。取兩項(xiàng)近似解：Ni x(1 x)2N2 x (1 x)2 一 、uN1a1N2a2a1x(1x) a2x(1 x)W1N1 , W2N2由公式(6-15)可以得到兩個(gè)加權(quán)積分方程，1 223、.° x(1 x)x a1 ( 2 x x ) a2(2 6x x x )dx 01 2223、.° x (1 x)x a1( 2 x x ) a2(2 6x x x )dx 0積分后可以得到一個(gè)二元一次方程組，解得，a10.1924, a20.1707近似解為， x(1 x)(0.1924 0.1707x)sin x該方程的精確解

10、為，u 叫三xsin1近似解與精確解的結(jié)果比較見表6-1,表6-1近似解與精確解比較x=0.25x=0.5x=0.75sin x uxsin 10.044010.069750.06006 x(1 x)(0.1924 0.1707x)0.044080.069440.06008假定單元的形函數(shù)為，N N14Nn單元結(jié)點(diǎn)的溫度為，Te T1 T2. TnT單元內(nèi)部的溫度分布為，T NTe以二維問題為例，說明用Galerkin法建立穩(wěn)態(tài)溫度場的一般有限元格式的過程。二維問題的穩(wěn)態(tài)熱傳導(dǎo)方程為，TT 二八一 x一 yQ 0(6-16a)x x y y第一類換熱邊界為Ts(6-16b)第二類換熱邊界條件為

11、,Tx nxXTy-ny qs y(6-16c)第三類邊界條件為,Tx nxXy ny h(Tf Ts) y(6-16d)在一個(gè)單元內(nèi)的加權(quán)積分公式為,eWi( xXT一)xT 一)Qd(6-17)由分部積分得,一(w1XT)XWi (XT)Xwi 一(XT)X一(w1yT) yWi(y T) ywi 一(y T) y應(yīng)用Green定理，一個(gè)單元內(nèi)的加權(quán)積分公式寫為,e w1(采用Wl x-) X T -nx xWlyTy -nyyT)wQd y)d 0Galerkin方法，選擇權(quán)函數(shù)為,Ni將單元內(nèi)的溫度分布函數(shù)和換熱邊界條件代入(6-18)(6-18)式，單元的加權(quán)積分公式為,e g( X

12、Xyg)/%yNiQdiqsd(6-19)NihNTedNihTfd 0換熱邊界條件代入后,在(6-19)式內(nèi)相應(yīng)出現(xiàn)了第二類換熱邊界項(xiàng)e& Mqsd3類換熱邊界項(xiàng)NihNTed3eNihTfd ，但沒有出現(xiàn)與第一類換熱邊界對(duì)應(yīng)的3項(xiàng)。這是因?yàn)?，采?Ni作為權(quán)函數(shù)，第一類換熱邊界被自動(dòng)滿足。寫成矩陣形式有,e14(Xe t 一NTQd回)4(y NTqsdyN)Ted y(6-20), _ _Te3 hN NT d3NThTfd0公式(6-20)是n個(gè)聯(lián)立的線性方程組，可以確定n個(gè)結(jié)點(diǎn)的溫度Ti o按有限元格式將(6-20)表不為,KeTe Pe(6-21)T e為單元的結(jié)點(diǎn)溫度向量

13、，Pe其中矩陣Ke為單元的導(dǎo)熱矩陣或稱為溫度剛度矩陣,稱為單元的溫度載荷向量或熱載荷向量( Thermal load vector)。對(duì)于某個(gè)特定單元，單元導(dǎo)熱矩陣Ke和溫度載荷向量P e的元素分別為,KijNi NjNiNj j)d yhNiN jd(6-22)MqsdNihTf dNiQd(6-23)如果某個(gè)單元完全處于物體的內(nèi)部,KijNi Nj(XX XNiyyNj)d yPiNiQd在整個(gè)物體上的加權(quán)積分方程是單元積分方程的和,四)與T(yEdTNTQdNTqsd(6-24)3 hNTNTed3NThTfd0根據(jù)單元結(jié)點(diǎn)的局部編號(hào)與整體編號(hào)的關(guān)系，直接求和得到整體剛度矩陣，整體方程組

14、 為,KT P6.3三角形單元的有限元列式(號(hào)圖6-1三角形單元回顧第三章的內(nèi)容可以發(fā)現(xiàn)， 與計(jì)算彈性力學(xué)平面問題時(shí)所采用的方法一樣，二維溫度場問題計(jì)算中所采用的三角形單元可以使用相同的形函數(shù)，Ni2A(aiNj2A(ajNm2A(ambi x Ci y)bjX Cj y)bmx Cmy)ai Xj ymXmyjbiyj ym Ci XmXjajXmyiXiymbjYm N cjXiXmam Xiyj Xj yibmyyjCmXjXi1Xiyi1XjyjT1Xmy mT 2A單元內(nèi)的溫度分布用結(jié)點(diǎn)上的溫度值表示為，TiT NiNj Nm Tj(6-25)Tm在三角形單元上，采用 Galerki

15、n法可得，A NT ( x) ( y) QdA 0(6-26)x x y yNT7(xT I(NT)Xx NT TX假定單元內(nèi)的導(dǎo)熱系數(shù)為常數(shù),<1dA14A2bibj bi bmbjbmTi Tj TmdA4Ah2bibjbibmbi bjb2bjbmbib bjb bmTi Tj Tm(6-27)NTT) y一(NyT) yNT T yNTyc©4A jGcm單元的剛度矩陣為,Kebi2bibjbibmdA yG52 cjcjcmbibjb2bjbm顯然，單元的導(dǎo)熱矩陣是對(duì)稱的。如果單元的內(nèi)部熱源為常數(shù),A NTQdA14A2Gcj cicmCjcm Ti Tj TmdAc

16、icm沁2cmbi bjbjbm廿由Green公式可得INTT) X(NT xT一nx xTi Tj Tm2ci(6-28)GCjcicmG62 cjcjcmcicm由內(nèi)部熱源產(chǎn)生的溫度載荷項(xiàng)為,Ni Nj NmdA1QA .131(6-29)一(NTyTdA y(6-30)NT y ny)dSy方便起見，把換熱邊界統(tǒng)一表示為第三類換熱邊界,A (NT x) (NT y-dA(6-30)A xx yy_ _T _ _ _ T _ _ _Te _os h N (Tf Ts)dS os hN TfdS os hN NT dS如果在單元邊上存在熱交換，各條邊上的邊界換熱條件在單元?jiǎng)偠染仃囍猩傻母郊?/p>

17、項(xiàng)為,2 1 0 e hljKe -12 060 0 00 0 0 hl.Kej- 0 2 160 1 22 0 1Ke 皿 0 0 0 61 0 2(6-31)(6-32)(6-33)由邊界換熱條件生成的溫度載荷向量為,P1 hTflj120(6-34)P0 hTfljm / 121(6-35)P1 hTflmi021(6-36)6.4溫度場分析舉例正方形截面的煙囪如圖 6-2所示，煙囪由混凝土建造，邊長為60cm,通道的邊長為20cm, 混凝土的導(dǎo)熱系數(shù)為 k 1.4W/(m K)。假定煙囪內(nèi)表面的溫度為100C,煙囪外表面暴露在空氣中，空氣的溫度為 30C,換熱系數(shù)為h 20W/(m2 K)。計(jì)算煙囪截面內(nèi)的穩(wěn)態(tài)溫度場。(參見，F(xiàn)inite Element Method Theory and Application with ANSYS, p279)圖6-2煙囪截面圖6-3有限元模型圖6-5熱流量分布穩(wěn)態(tài)溫度場分布與物體的初始狀態(tài)無關(guān)，那么是否與材料的導(dǎo)熱系數(shù)相關(guān)？我們把煙囪的模型做些修改，假定煙囪壁由兩層材料構(gòu)成。內(nèi)層材料為混凝土，外表面的截面尺寸為30cm 30cm ,煙囪通道 的尺寸不 變，仍為20cm 20cm