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

1、6.穩(wěn)態(tài)熱傳導(dǎo)問題的有限元法本章的內(nèi)容如下：6熱傳導(dǎo)方程與換熱邊界6.2穩(wěn)態(tài)溫度場分析的一般有限元列式6.3三角形單元的冇限元列式6.4溫度場分析舉例6熱傳導(dǎo)方程與換熱邊界在分析工程問題時，經(jīng)常要了解工件內(nèi)部的溫度分布情況，例如發(fā)動機的工作溫度、金 屬工件在熱處理過程屮的溫度變化、流體溫度分布等。物體內(nèi)部的溫度分布収決于物體內(nèi)部 的熱量交換，以及物體與外部介質(zhì)z間的熱量交換，一般認為是與時間相關(guān)的。物體內(nèi)部的 熱交換采用以下的熱傳導(dǎo)方程（fourier方程）來描述，dt dx i_8_+ 內(nèi)1dy)d+ 一8z2魯卜。i qz丿(6-1)式中。為密度，kg/m c為比熱容，j/（kgk）；,

2、av, 為導(dǎo)熱系數(shù)，w/（mk）； t為溫度，°c；（為時間，s； 為內(nèi)熱源密度，w/naa2t “ a2t對于各向同性材料，不同方向上的導(dǎo)熱系數(shù)相同，熱傳導(dǎo)方程可寫為以下形式,(6-2)顯:7卑+2卑+2卑+6t 5x2 sy2 az2除了熱傳導(dǎo)方程，計算物體內(nèi)部的溫度分布，還需要指定初始條件和邊界條件。初始條 件是指物體最初的溫度分布情況，tt=0 =t0(x,y,z)(6-3)邊界條件是指物體外表而與周圍環(huán)境的熱交換情況。在傳熱學(xué)中一般把邊界條件分為三 類。1）給定物體邊界上的溫度，稱為第一類邊界條件。物體表血上的溫度或溫度函數(shù)為已知，或 t =t,（x,y,z,/）（6-4）

3、s2）給定物體邊界上的熱量輸入或輸出，稱為第二類邊界條件。 已知物體表面上熱流密度，3) 給定對流換熱條件，稱為第三類邊界條件。物體與其相接觸的流體介質(zhì)z間的對流換熱系數(shù)和介質(zhì)的溫度為已知。dt , dt(6-6)dtdz 、其中h為換熱系數(shù)，w/(m2k)； 7；是物體表面的溫度；7是介質(zhì)溫度。如果邊界上的換熱條件不隨時間變化，物體內(nèi)部的熱源也不隨時間變化，在經(jīng)過一定時 間的熱交換后，物體內(nèi)各點溫度也將不隨時間變化，即竺=0dt這類問題稱為穩(wěn)態(tài)(steady state)熱傳導(dǎo)問題。穩(wěn)態(tài)熱傳導(dǎo)問題并不是溫度場不隨時間 的變化，而是指溫度分布穩(wěn)定后的狀態(tài)，我們不關(guān)心物體內(nèi)部的溫度場如何從初始狀

4、態(tài)過渡 到最后的穩(wěn)定溫度場。隨時間變化的瞬態(tài)(transient)熱傳導(dǎo)方程就退化為穩(wěn)態(tài)熱傳導(dǎo)方程, 三維問題的穩(wěn)態(tài)熱傳導(dǎo)方程為，a5xdx)乜丿聖、dz丿(6-7)對于各向同性的材料，可以得到以卜-的方程，稱為poisson方程,卑*卑*啤+2=odx2 dy2 qz a(6-8)考慮物體不包含內(nèi)熱源的情況，各向同性材料中的溫度場滿足laplace方程,(6-9)a2t a2t a2t在分析穩(wěn)態(tài)熱傳導(dǎo)問題吋，不需要考慮物體的初始溫度分布對最后的穩(wěn)定溫度場的影 響，因此不必考慮溫度場的初始條件，而只需考慮換熱邊界條件。計算穩(wěn)態(tài)溫度場實際上是 求解偏微分方程的邊值問題。溫度場是標(biāo)量場，將物體離散

5、成有限單元后，每個單元結(jié)點上 只有-個溫度未知數(shù)，比彈性力學(xué)問題要簡單。進行溫度場計算時有限單元的形函數(shù)與彈性 力學(xué)問題計算時的完全一致，單元內(nèi)部的溫度分布用單元的形函數(shù)，由單元結(jié)點上的溫度來 確定。由于實際工程問題屮的換熱邊界條件比較復(fù)雜，在許多場合下也很難進行測最，如何 定義正確的換熱邊界條件是溫度場計算的-個難點。6.2穩(wěn)態(tài)溫度場分析的一般冇限元列式在前面我們已經(jīng)介紹了有限元方法可以用來分析場問題，穩(wěn)態(tài)溫度場計算是一個典型的 場問題。我們可以采用虛功方程建立彈性力學(xué)問題分析的有限元格式，推導(dǎo)出的單元剛度矩 陣有明確的力學(xué)含義。在這里，介紹如何用加權(quán)余.量法(weighted residu

6、al method)建立穩(wěn) 態(tài)溫度場分析的有限元列式。微分方程的邊值問題，可以一燉地表示為未知函數(shù)u滿足微分方程組，a（w）a（w） = < a2（u） > = 0（在域。內(nèi)）（6-10） 未知函數(shù)u還滿足邊界條件，bgb（w） = < b2（m） > = 0（在邊界上）（6-11） 如果未知函數(shù)u是上述邊值問題的精確解，則在域屮的任一點上u都滿足微分方程 （6-10）,在邊界的任一點上都滿足邊界條件（6-11）。對于復(fù)雜的工程問題，這樣的精確解 往往很難找到，需耍設(shè)法尋找近似解。所選取的近似解是一族帶有待定參數(shù)的已知函數(shù)，一 般表示為« u 帀=工 niai

7、 = nai=(6-12)其小勺為待定系數(shù)，n,為己知函數(shù)，被稱為試探函數(shù)。試探函數(shù)要取占完全的函數(shù)序列,是線性獨立的。由于試探函數(shù)是完全的函數(shù)序列，任一函數(shù)都對以用這個序列來表示。 采用這種形式的近似解不能精確地滿足微分方程和邊界條件，所產(chǎn)生的謀差就稱為余fio微分方程（6-10）的余量為，r = a(na)(6-13)邊界條件（611）的余量為，r = b(na)(6-14)選擇一族已知的函數(shù)，使余量的加權(quán)積分為零，強迫近似解所產(chǎn)生的余量在某種平均意 義上等于零，£ w：rdq+ wjr/r = 0（6-15）w/和兀稱為權(quán)函數(shù)，通過公式（6-15）可以選擇待定的參數(shù)勺。這種采用

8、使余量的加權(quán)積分為零來求得微分方程近似解的方法稱為加權(quán)余量法。對權(quán)函 數(shù)的不同選擇就得到了不同的加權(quán)余量法，常用的方法包括配點法、了域法、最小二乘法、 力矩法和伽遼金法（galcrkin method）o在很多情況下，采用galcrkin法得到的方程組的系 數(shù)矩陣是對稱的，在這里也采用galerkin法建立穩(wěn)態(tài)溫度場分析的一般有限元列式。在 galerkin法屮，直接采用試探函數(shù)序列作為權(quán)函數(shù)，取= n., wj = -n . 0下而用求解二階常微分方程為例，說hj galerkin法（參見，王勖成編苦“有限元法基本 原理和數(shù)值方法”的123節(jié)）。例，求解二階常微分方程j 2牛 + u + x

9、 = 0（0 < x < 1）dx2邊界條件：當(dāng)x = 0時，w=0；當(dāng)兀=1時，w=0o 取兩項近似解：n、= x(l - x)n2 =x2(1-x)u = n仇 + n2a2 = ax(y 一工)+ a2x (1 一 x)州=m, w2 = n2由公式(6-15)可以得到兩個加權(quán)積分方程，x( - x)x + ax (-2 + x-x2) + 6z2(2-6x + x2 - x3)dx = 0(1 兀)兀 + % (_2 x (2 6x + 兀2 _ )dx = 0積分后可以得到一個二元一次方程組，解得，% =0.1924,a 2 = 0.1707近似解為,u=x(l-x)(0

10、.1924 + 0707x)該方程的粕確解為，u=-xsinl近似解與楮確解的結(jié)果比較見表6-1,表6-1近似解與精確解比較x=025x=0.5x=0.75sinxu =xsin 10.044010.069750.06006m = x(l - x)(0.1924 + 0.1707x)0.044080.069440.06008假定單元的形函數(shù)為，n = n,皿nj單元結(jié)點的溫度為，te = t2 . tnt單元內(nèi)部的溫度分布為，r = 7vte以二維問題為例，說明用galerkin法建立穩(wěn)態(tài)溫度場的一般有限元格式的過程。二維問 題的穩(wěn)態(tài)熱傳導(dǎo)方程為，(6-16a)聖、dx)第一類換熱邊界為(6-

11、16c)第二類換熱邊界條件為,oxdy第三類邊界條件為,(6-16d)在一個單元內(nèi)的加權(quán)積分公式為,；叩？仇浮)+ £仇普) + odg = o1 ox ex dy dy(6-17)由分部積分得,別譏甞)=譽“乎+時仇乎知認寥)=譽久甞)+時0乎dydy dy dy dy dy應(yīng)用green定理，一個單元內(nèi)的加權(quán)積分公式寫為, r r5yvj of svvj cti嘰孰+4軌)dr = 0采jij galerkin方法，選擇權(quán)函數(shù)為，(6-18)將單元內(nèi)的溫度分布函數(shù)和換熱邊界條件代入(6-18)式,£ >(4+讐久警川門掘4 dxoy dy-£ nqd&#

12、163;l-2niqsdv+ l tv""""- nihtfdr = ody dy單元的加權(quán)積分公式為,(6-19)換熱邊界條件代入后，在(6式內(nèi)相應(yīng)岀現(xiàn)了第二類換熱邊界項-j ndt,第三類換熱邊界項匸n/n7t"-匚ngdr ,但沒有出現(xiàn)與笫一類換熱邊界對應(yīng)的項。這是因為，釆用tv,作為權(quán)函數(shù)，第一類換熱邊界被自動滿足。寫成矩陣形式有,(叨羅4 oxox oycy(6-20)-£ n! qdsl-n1 cisdr+燦77- nrhtfdv = q公式(6-20)是n個聯(lián)立的線性方程組，可以確定n個結(jié)點的溫度7；。按有限元格式將(6

13、-20)表示為,kyte=pe(6-21)其屮炬陣kf為單元的導(dǎo)熱矩陣或稱為溫度剛度矩陣，t)e為單元的結(jié)點溫度向量，p° 稱為單元的溫度載荷向量或熱載荷向量(thermal load vector)o對于某個特定單元，單元導(dǎo) 熱矩陣kf和溫度載荷向量pf的元素分別為，(6-22)的 cnjc：n(，nj 十hnjnjdtx dx dx y dy dyz nihtfdt+ £ n©dv(6-23)如果某個單元完全處于物體的內(nèi)部,k廠(人如殂+人殂殂jx dx dx )dy dyp, = £ n © 在在整個物體上的加權(quán)積分方程是單元積分方程的和

14、,ze讐)l7tdq(6-24)-e £【ng工帥15聲+ee工 匚燦ntntecir-y £ nthtfcir = 0根據(jù)單元結(jié)點的局部編號與整休編號的關(guān)系，血接求和得到整體剛度矩陣，整體方程紐 為，kt = p6.3三角形單元的有限元列式冋顧第三章的內(nèi)容可以發(fā)現(xiàn)，與計算彈性力學(xué)平而問題時所采用的方法一樣，二維溫度 場問題計算中所采用的三角形單元可以使用相同的形函數(shù)，ni =方（也+勺兀+沙nj =云（勺+巧兀+勺刃v” 二丄加 + x + 5）2aai = xjymxmyj bj =yi-ym cj =xtn-xjaj =xmyi-xiymbj = ym - cj =無

15、-兀”am =兀兒一勺幾 bm =yi-yj cm =xj-xi1 x： y, 11 xj yj =t j心兒_t =2a單元內(nèi)的溫度分布用結(jié)點上的溫度值表示為,、t=lnt nj nm t（6-25）t川丿在三角形單元上，采用galerkin法可得，(6-26)可吟 h (血 m) + da = 0 ox ox oy oy燈？(人冬ox ox oxoxdx dx假定單元內(nèi)的導(dǎo)熱系數(shù)為常數(shù),4a2(6-27)4a®bjbib,nb®bjb幾sbt,町軟島)(心浮dy dy dydydy dyj血皿竺da 山 dy dygm >dajtrri jc?(6-28)a4a單

16、元的剛度矩陣為，b-k b”jb；bbmc.2wjac ；55cjj£顯然，單元的導(dǎo)熱矩陣是對稱的。如果單元的內(nèi)部熱源為常數(shù)，由內(nèi)部熱源產(chǎn)生的溫度載荷項為,ntqda = q由green公式可得f °“入口7 0t5t dt（"人=）+（"人必ja dxdx dydy(6-29)(6-30)=j （仏軌+可九軌）ds 方便起見，把換熱邊界統(tǒng)一表示為第三類換熱邊界，夕(ns浮)+ g(n玄知必oxox dyoy(6-30)燦 nv (t-ts)ds = 燦 nytjds-?燦 nrw7tds如果在單元邊上存在熱交換，各條邊上的邊界換熱條件在單元剛度矩陣屮生

17、成的附加項為,ke21_0(6-31)2hlk“ =亠 0012(6-32)(6-33)山邊界換熱條件牛成的溫度載荷向量為,py(6-34)pe-/(6-35)py(6-36)6.4溫度場分析舉例止方形截而的煙囪如圖6-2所示，煙囪由混凝土建造，邊長為60cm,通道的邊長為20cm, 混凝土的導(dǎo)熱系數(shù)為r=1.4w/(mk)。假定煙囪內(nèi)表面的溫度為100°c,煙囪外表面暴鉤 在空氣中，空氣的溫度為30°c,換熱系數(shù)為h = 20w /(/h2-)o計算煙囪截面內(nèi)的穩(wěn)態(tài)溫 度場。(參見，finite element method theory and application

18、with ansys, p279 )20cm60cm圖62煙囪截面modal solutionstep"丄sub .丄time 丄temp(avg)b3y3-03hm *34.1623hx 丄0034.182圖6-3有限元模型46.808ore is zooz23.33.s541.49556.12170.74885.374100圖6-4穩(wěn)態(tài)溫度分布ansys圖6-5熱流量分布穩(wěn)態(tài)溫度場分布與物體的初始狀態(tài)無關(guān),那么是否與材料的導(dǎo)熱系數(shù)相關(guān)？我們把煙囪 的模型做些修改，假定煙囪壁由兩層材料構(gòu)成。內(nèi)層材料為混凝土，外表面的截面尺寸為 30cmx30cm ,煙囪通道的尺寸不變,仍為20cmx20cm o外層材料的導(dǎo)熱系數(shù)為 k = 0aw /（m-a：）,外部表面的截血尺寸不變，內(nèi)部表面的截面尺寸為30cm x 30cm o換熱邊界條件不變，雙層煙囪的冇限元模型如圖6-6所示。圖6-6雙層煙囪的有限元模型1hodal solution3tep=1sub =1圖67雙層煙囪的溫度分布dec 19 200200:11:0737.97853.4846