计算传热学第4讲扩散方程的数值解

1、计算传热学计算传热学 第第4讲讲主要内容主要内容主要目的主要目的阅读与作业阅读与作业l阅读要求：陶文铨数值传热学第4章l作业：P124 题41；P125 题47l完成课外作业第一题和第二题4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解(1) 0)()(1SdxdTxAdxdxA4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE图图 1 1 一维问题空间区域的离散化4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解在通常情况下，SS（T）线性化： SScSpT （2）其中，按负斜率源项原则， SpSp(T*

2、) 0 （3）4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解将方程（1）两边通乘A（x），并对x从w到e积分：x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE图图 1 1 一维问题空间区域的离散化 0)()(xSAdxdTxAdxd(4) 0)()(ewewdxxSAdxdxdTxAdxd4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解weewdxdTxAdxdTxAdxdxdTxAdxd)()()(wweedxdTAdxdTA)()(5) )()()()(wWPwePEexTTAxTTAx(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE

3、图图 1 1 一维问题空间区域的离散化4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解待求变量T在控制容积P上为常数整个控制容积的A(x)为常数，且等于P点的值。ewpcewdxxATSSdxxSA)()()(6) )()(PPpPcxTASxAS4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解l将（5）和（6）代入方程（4），)()(WPwPEeTTxATTxA(7) 0)()(PPpPcxTASxASl整理后得到，(8) PEEWWPPbTaTaTal其中，4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解(8) PEEWWPPbTaTaTal其中，(8a) wWxA

4、a(8b) eExAa(8c) )(xASaaaPpEWPl下标：大写字母表示在节点处取值，小写字母表示在相应的控制面处取值(8d) )(xASbPcP4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解：对源项积分时采用线性分布(8a) )(2)()(PPwwwwWASxxxxAa(8c) )()()()()()(22)(eeewwwPpEWPxxxxxxxASaaa(8b) )(2)()(PPeeeeEASxxxxAa4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解l

5、把跃变界面作为边界把跃变界面作为边界可以考虑接触热阻rc (Wm2)/K满足流的唯一性原则，iiidxdTdxdTqciiirqTTiTiT4.1 一维稳态导热问题的数值解一维稳态导热问题的数值解直角坐标左边界，第二类边界条件qBx=0(8) 0BxqdxdT注意：qB的正方向与x轴的正方向一致！ee边界条件的处理边界条件的处理l边界上出现半个控制容积qBx=0BxqxTTdxdT11210)( 123(x)1(x)2l边界节点的差分方程可以用下述方法推出：(9) )( 1120 xTTdxdTx边界条件的处理边界条件的处理l整理后得到：(10) )(1121BqxTTl特点：l只有边界条件的

6、处理边界条件的处理l在研究边界节点所代表的控制容积（元体）的能量平衡流入CV的能量内热源发出的热量流出CV的能量1121)( xTTqe123eeqBx=0(x)1(x)21(x)1Sq1流入CV的能量通过边界流入的热量qB通过控制面流入的热量q1内热源发出的热量 (x)1S流出CV的能量0边界条件的处理边界条件的处理l代入能量守恒关系，0)(21)(1112SxxTTqeBl整理后得到，(11) )(21)(1121SxqxTTBel特点l灵活，便于处理各种复杂的边界条件l二阶精度，与内部节点精度等级匹配（）l推导过程较繁边界条件的处理边界条件的处理在直角坐标的条件下，方程（1）变为， 0S

7、dxdTdxd123eeqBx=0(x)1(x)21假定在节点12之间的导热系数为常数，且恒等于e，则有，对于点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推导繁琐边界条件的处理边界条件的处理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边界节点的控制容积或它所代表的求解区域？v边界节点的控制容积0v于是，按元体能量平衡法或其它二阶精度的方法，令与源项对应的项等于0，得到，(19) )(1121BqxTT(15) )(21)(11121SxqxTTBee边界条件的处理边界条件的处理尽管它与一阶Taylor级数展开法的结果形式上相同，但它却是二阶精度的！请大家证明这一结论。l采用划分网格时，即使在均匀网络的

10、前提下，第1个也的。qBx=0(x)2( x)3123(x)1(x)2l从图中可以清楚地看出这一点l即使 (x)2= (x)3 ( x)1 ( x)2l所以要对第一个内部节点给予特别注意。边界条件的处理边界条件的处理l例如，对于直角坐标系，对于节点2，qBx=0(x)2( x)3123(x)1(x)2(20) 2312bTaTaTaEWP )(11xxawW )()(22xASaaapEWP )(2xxaeeE22)()(xASbcP注意：边界条件的处理边界条件的处理l或写成，(21) )()( )()()()(23211122211xSTxTxTxSxxcepe边界条件的处理边界条件的处理l

11、附加源项法（Additional source term method)以内节点法为例由方程（19）解出边界节点上的待求变量T1，(22) )(1121BqxTT代入与第1个近边界节点的差分方程（21），边界条件的处理边界条件的处理l代入与第1个近边界节点的差分方程（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对于，也可以做类似的处理，但是这时，l请大家证明，(31) )(2,xTSfadc(32) )(2,xSadp边界条件的处理边界条件的处理边界节点消去法不仅能用于内节点网格，也能用于外节点网格计算附加源项：Sc,ad，Sp,ad把附加源项计入该控制容积中的源项中令与边界节点对应的系数（aW）等于0是

13、传热问题数值计算之一的基础地位l尽可能采用划分网格广泛应用提高收敛速度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）系数阵A的特殊性，通常称之为三对角方程（Tri-diagonal equation)l三对角方程可以采用非常高效的追赶法（TDMA法）求解l基于矩阵分解l属于必须掌握的内容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 计算机实现：算例计算机实现：算例(36) 0SdxdTdxd(37) 21 10 33xbTxaT(38) ),(xTSS (36) 1 ,200 xxTTT4.1.8 计算机实现：算例计算机实现：算例内节点法：先划分控制容积，在确定节点均匀网格：x=x将整个求解区域划分为（N-2）个控制容积，N个节点（包括2个边界节点）(39) 11iiiiiiidTcTbTa1)-( 4,., , 3 , 2Ni 4.1.8 计算机实现：算例计算机实现：算例l其中，(40) eiwixcxa(40)

16、)(xScabipiii(40) )(xSdici(40) 2 21111iiiieiiiiw4.1.8 计算机实现：算例计算机实现：算例：采用内节点法划分网格时，近边界节点与其它内部节点不尽相同，所以必须单独考虑。123NN-1(x)w(x)e(x)ev(x)w= xv w= W= 1v所以，当i2时4.1.8 计算机实现：算例计算机实现：算例l所以，当i2时， 2 23232111xxcxxaewxScabp2111)( )(0121TaxSdc(41) 13121dTcTb4.1.8 计算机实现：算例计算机实现：算例l当i3，4，。，N-2时，xcxaeiwi11 xScabipiii)

17、(111 )(11xSdici 2 21111iiiieiiiiw(42) 111111iiiiiiidTcTbTa4.1.8 计算机实现：算例计算机实现：算例v(x)e= xv e= E= Nv所以，当i N-1时123NN-1(x)w(x)e(x)e4.1.8 计算机实现：算例计算机实现：算例 2 2 212122xxcxxaNeNNNNNwNxScabNpNNN1222)( 1)(T )(212NNNcNcxSd(44) 21222NNNNNdTbTa4.1.8 计算机实现：算例计算机实现：算例l最后得到由（N2）个方程构成的方程组为(41) 13121dTcTb(42) 111111i

18、iiiiiidTcTbTa2)- 5,., , 4 , 3(Ni (44) 21222NNNNNdTbTal求解上面的方程，即可得到（N2）个未知数，即，T2, T3, T4,., TN-1。4.1.8 计算机实现：算例计算机实现：算例l注意：上面的方程组的，必须用求解l求解方法：v假定一个温度分布：Ti，i1，2，3，。，Nv计算i ，i1，2，3，。，Nv计算a, b, c, dv用TDMA法求解方程组，得到新的温度分布： Tiv计算：Maxabs(Ti -Ti), i1,2,3,Nv判断： abs(Ti -Ti)是否小于（精度要求）v如果不能满足精度要求，令Ti Ti，v满足精度要求：4

19、.1.8 计算机实现：算例计算机实现：算例v希望大家用计算机完成上面的计算，并与下面的分析解结果比较：10 2)1 (110 )1 (4040404xbaxTbxbaxTbTT特别提示特别提示：掌握循环变量的使用：对算法清晰透彻的把握：细心再细心4.2 一维非稳态导热问题一维非稳态导热问题(45) )()(1tTcSxTxAxxA4.2.1 非稳态项的处理非稳态项的处理tk=(k-1)t, k=1,2,3, (46)v至少有三种方案：(47a) )O( 2211ttTTtTkPkPkP4.2.1 非稳态项的处理非稳态项的处理(47b) )O( 1ttTTtTkPkPkP(47c) )O( 1t

20、tTTtTkPkPkP4.2.2 控制方程的离散化控制方程的离散化v参照图示的节点组，对于任意一个CV，x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE图图 1 1 一维问题空间区域的离散化(48) )()(dxtTcdxxASdxxTxAxkewewkkew(49a) )(kwkekewxTAxTAdxxTxAx4.2.2 控制方程的离散化控制方程的离散化v将方程（49）代入方程（48）(49b) )()( )()(kPkPpkPcewkpcewkxTASxASdxATSSdxxAS(49c) xtTcdxtTckPkew(50) )( )( xtTcxTASxASxTA

21、xTAkPkPkPpkPckwke4.2.2 控制方程的离散化控制方程的离散化l假定节点间按线性分布，则(51a) )(kPkEekeTTxAxTAx(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE图图 1 1 一维问题空间区域的离散化(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将之发现：的存在并差分方程的：vaP中多出了a0P项vbP中多出了a0PTk-1P项v每前进一个时间步长都要求解一个方程组4.2.3 边界条件的处理边界条件的处理123eeqBx=0(x)1(x)21(x)1Sq1(54) )(2111tBESxqq0)(21)(1112SxxTTqeBv其中Et是CV单位时间内能的

23、变化，显然，(55) )(211111tTTcxEkkt边界条件的处理边界条件的处理v而，(56) )(1121xTTqkkev将q1和 Et代入方程（54），有(57) )(21)(21)(11111112tTTcxSxxTTqkkkkeBv整理后得到，边界条件的处理边界条件的处理v它与稳态问题相比，多出了内能变化项(58) )(21)( )(2)(2111111212121SxqxTtxcTTtxcBekekke4.3.4求解过程求解过程开始，开始，t0t=t+ t计算有关系数计算有关系数形成方程组，进行求解形成方程组，进行求解输出结果输出结果ttcal是STOP否特别提示特别提示v显式格

24、式（显式格式（Explicit formulations）v时间步进法：Time-marchingv不需要求解方程组v程序简单，对计算机的内存要求低v稳定性差v每推进一个时间步长，都需要求解一个方程组v程序复杂，要求计算机有较大的内存v稳定性好4.3 多维非稳态导热问题的离散化多维非稳态导热问题的离散化v以直角坐标系中的二维导热为例(59) tTcSyTyxTx(60) )(T(61) TSSSpc4.3 多维非稳态导热问题的离散化多维非稳态导热问题的离散化SNWEPsnwexy(y)n(y)s(y)-n(y)+s(x)w(x)e(x)-e(x)+w4.3 多维非稳态导热问题的离散化多维非稳态

25、导热问题的离散化基本思路与一维问题完全相同对方程（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节点间按线性分布，则)(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