計(jì)算傳熱學(xué)第4講擴(kuò)散方程的數(shù)值解

1、計(jì)算傳熱學(xué)計(jì)算傳熱學(xué) 第第4講講主要內(nèi)容主要內(nèi)容主要目的主要目的閱讀與作業(yè)閱讀與作業(yè)l閱讀要求：陶文銓數(shù)值傳熱學(xué)第4章l作業(yè)：P124 題41；P125 題47l完成課外作業(yè)第一題和第二題4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解(1) 0)()(1SdxdTxAdxdxA4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖圖 1 1 一維問題空間區(qū)域的離散化4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解在通常情況下，SS（T）線性化： SScSpT （2）其中，按負(fù)斜率源項(xiàng)原則， SpSp(T*

2、) 0 （3）4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解將方程（1）兩邊通乘A（x），并對x從w到e積分：x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖圖 1 1 一維問題空間區(qū)域的離散化 0)()(xSAdxdTxAdxd(4) 0)()(ewewdxxSAdxdxdTxAdxd4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解weewdxdTxAdxdTxAdxdxdTxAdxd)()()(wweedxdTAdxdTA)()(5) )()()()(wWPwePEexTTAxTTAx(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE

3、圖圖 1 1 一維問題空間區(qū)域的離散化4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解待求變量T在控制容積P上為常數(shù)整個(gè)控制容積的A(x)為常數(shù)，且等于P點(diǎn)的值。ewpcewdxxATSSdxxSA)()()(6) )()(PPpPcxTASxAS4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解l將（5）和（6）代入方程（4），)()(WPwPEeTTxATTxA(7) 0)()(PPpPcxTASxASl整理后得到，(8) PEEWWPPbTaTaTal其中，4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解(8) PEEWWPPbTaTaTal其中，(8a) wWxA

4、a(8b) eExAa(8c) )(xASaaaPpEWPl下標(biāo)：大寫字母表示在節(jié)點(diǎn)處取值，小寫字母表示在相應(yīng)的控制面處取值(8d) )(xASbPcP4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解：對源項(xiàng)積分時(shí)采用線性分布(8a) )(2)()(PPwwwwWASxxxxAa(8c) )()()()()()(22)(eeewwwPpEWPxxxxxxxASaaa(8b) )(2)()(PPeeeeEASxxxxAa4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解l

5、把躍變界面作為邊界把躍變界面作為邊界可以考慮接觸熱阻rc (Wm2)/K滿足流的唯一性原則，iiidxdTdxdTqciiirqTTiTiT4.1 一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解直角坐標(biāo)左邊界，第二類邊界條件qBx=0(8) 0BxqdxdT注意：qB的正方向與x軸的正方向一致！ee邊界條件的處理邊界條件的處理l邊界上出現(xiàn)半個(gè)控制容積qBx=0BxqxTTdxdT11210)( 123(x)1(x)2l邊界節(jié)點(diǎn)的差分方程可以用下述方法推出：(9) )( 1120 xTTdxdTx邊界條件的處理邊界條件的處理l整理后得到：(10) )(1121BqxTTl特點(diǎn)：l只有邊界條件的

6、處理邊界條件的處理l在研究邊界節(jié)點(diǎn)所代表的控制容積（元體）的能量平衡流入CV的能量內(nèi)熱源發(fā)出的熱量流出CV的能量1121)( xTTqe123eeqBx=0(x)1(x)21(x)1Sq1流入CV的能量通過邊界流入的熱量qB通過控制面流入的熱量q1內(nèi)熱源發(fā)出的熱量 (x)1S流出CV的能量0邊界條件的處理邊界條件的處理l代入能量守恒關(guān)系，0)(21)(1112SxxTTqeBl整理后得到，(11) )(21)(1121SxqxTTBel特點(diǎn)l靈活，便于處理各種復(fù)雜的邊界條件l二階精度，與內(nèi)部節(jié)點(diǎn)精度等級匹配（）l推導(dǎo)過程較繁邊界條件的處理邊界條件的處理在直角坐標(biāo)的條件下，方程（1）變?yōu)椋?0S

7、dxdTdxd123eeqBx=0(x)1(x)21假定在節(jié)點(diǎn)12之間的導(dǎo)熱系數(shù)為常數(shù)，且恒等于e，則有，對于點(diǎn)e，(12) 022eeSdxTd邊界條件的處理邊界條件的處理另一方面，參照附圖，123eeqBx=0(x)1(x)21)(2)(211220 xOxdxTddxdTdxdTeex)(2)(211xOxSdxdTee(13) )(2)()(211112xOxSxTTe邊界條件的處理邊界條件的處理l所以，(14) 2)()(11120 xSxTTdxdTexl將之代入式（8）(8) 0BxqdxdT010 xxdxdTdxdTBeqxSxTT2)()(11121邊界條件的處理邊界條件的

8、處理(15) )(21)(11121SxqxTTBeel二階精度l不具有一般性l推導(dǎo)繁瑣邊界條件的處理邊界條件的處理123eeqBx=0(x)1(x)21(16) )()(21122220 xOxdxTddxdTdxdTx 按二階精度的差商公式 ) 1()1 ()(12212322TLTLTLxdxdTxxx21123222)()()1 (12xxTLTLTLdxTdxxx邊界條件的處理邊界條件的處理l代入式（16），整理后得到， (17) )2() 1()1 ()(1122320TLLTLTLxdxdTxxxxxl代入方程（8），BxxxxqTLLTLTLx)2() 1()1 ()(1223

9、21(8) 0BxqdxdTl整理后得到，(18) )1 ()() 1()2(123221BxxxxqLxTTLTLL邊界條件的處理邊界條件的處理邊界條件的處理邊界條件的處理qBx=0(x)2( x)3123(x)1(x)2v邊界節(jié)點(diǎn)的控制容積或它所代表的求解區(qū)域？v邊界節(jié)點(diǎn)的控制容積0v于是，按元體能量平衡法或其它二階精度的方法，令與源項(xiàng)對應(yīng)的項(xiàng)等于0，得到，(19) )(1121BqxTT(15) )(21)(11121SxqxTTBee邊界條件的處理邊界條件的處理盡管它與一階Taylor級數(shù)展開法的結(jié)果形式上相同，但它卻是二階精度的！請大家證明這一結(jié)論。l采用劃分網(wǎng)格時(shí)，即使在均勻網(wǎng)絡(luò)的

10、前提下，第1個(gè)也的。qBx=0(x)2( x)3123(x)1(x)2l從圖中可以清楚地看出這一點(diǎn)l即使 (x)2= (x)3 ( x)1 ( x)2l所以要對第一個(gè)內(nèi)部節(jié)點(diǎn)給予特別注意。邊界條件的處理邊界條件的處理l例如，對于直角坐標(biāo)系，對于節(jié)點(diǎn)2，qBx=0(x)2( x)3123(x)1(x)2(20) 2312bTaTaTaEWP )(11xxawW )()(22xASaaapEWP )(2xxaeeE22)()(xASbcP注意：邊界條件的處理邊界條件的處理l或?qū)懗桑?21) )()( )()()()(23211122211xSTxTxTxSxxcepe邊界條件的處理邊界條件的處理l

11、附加源項(xiàng)法（Additional source term method)以內(nèi)節(jié)點(diǎn)法為例由方程（19）解出邊界節(jié)點(diǎn)上的待求變量T1，(22) )(1121BqxTT代入與第1個(gè)近邊界節(jié)點(diǎn)的差分方程（21），邊界條件的處理邊界條件的處理l代入與第1個(gè)近邊界節(jié)點(diǎn)的差分方程（21）， (23) )()( )()()()()(2321121122211xSTxqxTxTxSxxceBpel整理后得到， 邊界條件的處理邊界條件的處理l整理后得到， (24) )()()()(232222BcepeqxSTxTxSxl或者寫成， (25) )(232xSTaTacEPl其中， 邊界條件的處理邊界條件的處理l其

12、中(25) )()( )()(,2222adccBcceEpePSSxqSSxaxSxa(26) )(2,xqSBadcAdditional source term!邊界條件的處理邊界條件的處理(28) TSSSpc(27) )(1TTqfB(29) ,adcccSSS(30) ,adpppSSSl對于，也可以做類似的處理，但是這時(shí)，l請大家證明，(31) )(2,xTSfadc(32) )(2,xSadp邊界條件的處理邊界條件的處理邊界節(jié)點(diǎn)消去法不僅能用于內(nèi)節(jié)點(diǎn)網(wǎng)格，也能用于外節(jié)點(diǎn)網(wǎng)格計(jì)算附加源項(xiàng)：Sc,ad，Sp,ad把附加源項(xiàng)計(jì)入該控制容積中的源項(xiàng)中令與邊界節(jié)點(diǎn)對應(yīng)的系數(shù)（aW）等于0是

13、傳熱問題數(shù)值計(jì)算之一的基礎(chǔ)地位l盡可能采用劃分網(wǎng)格廣泛應(yīng)用提高收斂速度l盡可能采用的離散化方程在形式上與的相同4.1.7 差分方程的求解差分方程的求解l將前面得到的差分方程（8）改寫為，(8) PEEWWPPbTaTaTa(33) )(PEEPPWWbTaTaTal或者簡單的寫成矩陣的形式，(34) DTA其中，(34a) ,.,1321TnnTTTTTTT(34b) ,.,1321TnnTDDDDDD4.1.7 差分方程的求解差分方程的求解(34c) . 11133322211nnnnnbacbacbacbacbA4.1.7 差分方程的求解差分方程的求解l與方程（33）對比，知，(35) ,

14、 (35) ,PEPWbDacabaal由方程（33）系數(shù)陣A的特殊性，通常稱之為三對角方程（Tri-diagonal equation)l三對角方程可以采用非常高效的追趕法（TDMA法）求解l基于矩陣分解l屬于必須掌握的內(nèi)容4.1.7 差分方程的求解差分方程的求解lTDMA法Fortran源程序 SUBROUTINE TDMA(A,B,C,D,N) REAL A(N),B(N),C(N),D(N) DO 5 I=1,N-1 F=A(I+1)/B(I) B(I+1)=B(I+1)-F*C(I) 5 D(I+1)=D(I+1)-F*D(I) D(N)=D(N)/B(N) DO 10 I=N-1,

15、1,-1 10 D(I)=(D(I)-C(I)*D(I+1)/B(I) RETURN END 4.1.8 計(jì)算機(jī)實(shí)現(xiàn)：算例計(jì)算機(jī)實(shí)現(xiàn)：算例(36) 0SdxdTdxd(37) 21 10 33xbTxaT(38) ),(xTSS (36) 1 ,200 xxTTT4.1.8 計(jì)算機(jī)實(shí)現(xiàn)：算例計(jì)算機(jī)實(shí)現(xiàn)：算例內(nèi)節(jié)點(diǎn)法：先劃分控制容積，在確定節(jié)點(diǎn)均勻網(wǎng)格：x=x將整個(gè)求解區(qū)域劃分為（N-2）個(gè)控制容積，N個(gè)節(jié)點(diǎn)（包括2個(gè)邊界節(jié)點(diǎn)）(39) 11iiiiiiidTcTbTa1)-( 4,., , 3 , 2Ni 4.1.8 計(jì)算機(jī)實(shí)現(xiàn)：算例計(jì)算機(jī)實(shí)現(xiàn)：算例l其中，(40) eiwixcxa(40)

16、)(xScabipiii(40) )(xSdici(40) 2 21111iiiieiiiiw4.1.8 計(jì)算機(jī)實(shí)現(xiàn)：算例計(jì)算機(jī)實(shí)現(xiàn)：算例：采用內(nèi)節(jié)點(diǎn)法劃分網(wǎng)格時(shí)，近邊界節(jié)點(diǎn)與其它內(nèi)部節(jié)點(diǎn)不盡相同，所以必須單獨(dú)考慮。123NN-1(x)w(x)e(x)ev(x)w= xv w= W= 1v所以，當(dāng)i2時(shí)4.1.8 計(jì)算機(jī)實(shí)現(xiàn)：算例計(jì)算機(jī)實(shí)現(xiàn)：算例l所以，當(dāng)i2時(shí)， 2 23232111xxcxxaewxScabp2111)( )(0121TaxSdc(41) 13121dTcTb4.1.8 計(jì)算機(jī)實(shí)現(xiàn)：算例計(jì)算機(jī)實(shí)現(xiàn)：算例l當(dāng)i3，4，。，N-2時(shí)，xcxaeiwi11 xScabipiii)

17、(111 )(11xSdici 2 21111iiiieiiiiw(42) 111111iiiiiiidTcTbTa4.1.8 計(jì)算機(jī)實(shí)現(xiàn)：算例計(jì)算機(jī)實(shí)現(xiàn)：算例v(x)e= xv e= E= Nv所以，當(dāng)i N-1時(shí)123NN-1(x)w(x)e(x)e4.1.8 計(jì)算機(jī)實(shí)現(xiàn)：算例計(jì)算機(jī)實(shí)現(xiàn)：算例 2 2 212122xxcxxaNeNNNNNwNxScabNpNNN1222)( 1)(T )(212NNNcNcxSd(44) 21222NNNNNdTbTa4.1.8 計(jì)算機(jī)實(shí)現(xiàn)：算例計(jì)算機(jī)實(shí)現(xiàn)：算例l最后得到由（N2）個(gè)方程構(gòu)成的方程組為(41) 13121dTcTb(42) 111111i

18、iiiiiidTcTbTa2)- 5,., , 4 , 3(Ni (44) 21222NNNNNdTbTal求解上面的方程，即可得到（N2）個(gè)未知數(shù)，即，T2, T3, T4,., TN-1。4.1.8 計(jì)算機(jī)實(shí)現(xiàn)：算例計(jì)算機(jī)實(shí)現(xiàn)：算例l注意：上面的方程組的，必須用求解l求解方法：v假定一個(gè)溫度分布：Ti，i1，2，3，。，Nv計(jì)算i ，i1，2，3，。，Nv計(jì)算a, b, c, dv用TDMA法求解方程組，得到新的溫度分布： Tiv計(jì)算：Maxabs(Ti -Ti), i1,2,3,Nv判斷： abs(Ti -Ti)是否小于（精度要求）v如果不能滿足精度要求，令Ti Ti，v滿足精度要求：4

19、.1.8 計(jì)算機(jī)實(shí)現(xiàn)：算例計(jì)算機(jī)實(shí)現(xiàn)：算例v希望大家用計(jì)算機(jī)完成上面的計(jì)算，并與下面的分析解結(jié)果比較：10 2)1 (110 )1 (4040404xbaxTbxbaxTbTT特別提示特別提示：掌握循環(huán)變量的使用：對算法清晰透徹的把握：細(xì)心再細(xì)心4.2 一維非穩(wěn)態(tài)導(dǎo)熱問題一維非穩(wěn)態(tài)導(dǎo)熱問題(45) )()(1tTcSxTxAxxA4.2.1 非穩(wěn)態(tài)項(xiàng)的處理非穩(wěn)態(tài)項(xiàng)的處理tk=(k-1)t, k=1,2,3, (46)v至少有三種方案：(47a) )O( 2211ttTTtTkPkPkP4.2.1 非穩(wěn)態(tài)項(xiàng)的處理非穩(wěn)態(tài)項(xiàng)的處理(47b) )O( 1ttTTtTkPkPkP(47c) )O( 1t

20、tTTtTkPkPkP4.2.2 控制方程的離散化控制方程的離散化v參照圖示的節(jié)點(diǎn)組，對于任意一個(gè)CV，x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖圖 1 1 一維問題空間區(qū)域的離散化(48) )()(dxtTcdxxASdxxTxAxkewewkkew(49a) )(kwkekewxTAxTAdxxTxAx4.2.2 控制方程的離散化控制方程的離散化v將方程（49）代入方程（48）(49b) )()( )()(kPkPpkPcewkpcewkxTASxASdxATSSdxxAS(49c) xtTcdxtTckPkew(50) )( )( xtTcxTASxASxTA

21、xTAkPkPkPpkPckwke4.2.2 控制方程的離散化控制方程的離散化l假定節(jié)點(diǎn)間按線性分布，則(51a) )(kPkEekeTTxAxTAx(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖圖 1 1 一維問題空間區(qū)域的離散化(51b) )(kWkPwkwTTxAxTAl而按式（47c），(51c) )O( 1ttTTtTkPkPkP4.2.2 控制方程的離散化控制方程的離散化v整理后得到，(52) )( )( )( )(1xtTTcxTASxASTTxATTxAkPkPkPkPpkPckWkPwkPkEe(53) PkEEkWWkPPbTaTaTav將式（51）代

22、入方程（50），4.2.2 控制方程的離散化控制方程的離散化l其中，(53a) wWxAa(53b) eExAa(53c) )(0 xASaaaaPpEWPP(53e) )(10 xASTabPckPPP(53d) )(0txcAaPP4.2.2 控制方程的離散化控制方程的離散化l將之發(fā)現(xiàn)：的存在并差分方程的：vaP中多出了a0P項(xiàng)vbP中多出了a0PTk-1P項(xiàng)v每前進(jìn)一個(gè)時(shí)間步長都要求解一個(gè)方程組4.2.3 邊界條件的處理邊界條件的處理123eeqBx=0(x)1(x)21(x)1Sq1(54) )(2111tBESxqq0)(21)(1112SxxTTqeBv其中Et是CV單位時(shí)間內(nèi)能的

23、變化，顯然，(55) )(211111tTTcxEkkt邊界條件的處理邊界條件的處理v而，(56) )(1121xTTqkkev將q1和 Et代入方程（54），有(57) )(21)(21)(11111112tTTcxSxxTTqkkkkeBv整理后得到，邊界條件的處理邊界條件的處理v它與穩(wěn)態(tài)問題相比，多出了內(nèi)能變化項(xiàng)(58) )(21)( )(2)(2111111212121SxqxTtxcTTtxcBekekke4.3.4求解過程求解過程開始，開始，t0t=t+ t計(jì)算有關(guān)系數(shù)計(jì)算有關(guān)系數(shù)形成方程組，進(jìn)行求解形成方程組，進(jìn)行求解輸出結(jié)果輸出結(jié)果ttcal是STOP否特別提示特別提示v顯式格

24、式（顯式格式（Explicit formulations）v時(shí)間步進(jìn)法：Time-marchingv不需要求解方程組v程序簡單，對計(jì)算機(jī)的內(nèi)存要求低v穩(wěn)定性差v每推進(jìn)一個(gè)時(shí)間步長，都需要求解一個(gè)方程組v程序復(fù)雜，要求計(jì)算機(jī)有較大的內(nèi)存v穩(wěn)定性好4.3 多維非穩(wěn)態(tài)導(dǎo)熱問題的離散化多維非穩(wěn)態(tài)導(dǎo)熱問題的離散化v以直角坐標(biāo)系中的二維導(dǎo)熱為例(59) tTcSyTyxTx(60) )(T(61) TSSSpc4.3 多維非穩(wěn)態(tài)導(dǎo)熱問題的離散化多維非穩(wěn)態(tài)導(dǎo)熱問題的離散化SNWEPsnwexy(y)n(y)s(y)-n(y)+s(x)w(x)e(x)-e(x)+w4.3 多維非穩(wěn)態(tài)導(dǎo)熱問題的離散化多維非穩(wěn)態(tài)

25、導(dǎo)熱問題的離散化基本思路與一維問題完全相同對方程（59）在控制容積CV上積分，(62) CVCVCVCVdxdytTcdxdySdxdyyTydxdyxTx4.3.3 控制方程的離散化控制方程的離散化l其中SNWEPsnwexy(y)n(y)s(y)-n(y)+s(x)w(x)e(x)-e(x)+wdydxxTxdxdyxTxnsewCV dyxTxTnswenswedyxTxT(63) yxTxTwe4.3.3 控制方程的離散化控制方程的離散化l節(jié)點(diǎn)間按線性分布，則)(PEeeTTxxT)(WPwwTTxxTl代入（63）得到，yTTxTTxdxdyxTxWPwPEeCV)()(64) )()()()(WwwPeewwEeeTxyTxyxyTxy4.3.3 控制方程的離散化控制方程的離散化l同樣SNWEPsnwexy(y)n(y)s(y)-n(y)+s(x)w(x)e(x)-e(x)+wdxdyyTydxdyyTyewnsCV dxyTyTewsnewsn