热问题的数值计算

1、热流问题的数值计算热流问题的数值计算numerical simulations of thermal & fluid problems第四章第四章 扩散方程的数值解法及应用扩散方程的数值解法及应用主讲主讲 李炎锋李炎锋2008年7月 北京4.1 一维导热问题 边界条件的处理，源项的线性化，代数方程的求解方法。 4.1.1 一维导热问题的通用控制方程 （1） 0)()(1sdxdtxadxdxa4.1.2 用控制容积积分法导出通用方程的 离散形式 假定源项s在任一控制容积内中的值可以表示温度的线性函数： （2） sc为常数，sp为s = f(t) 的曲线在p点的斜率，且规定恒为负值。pp

2、ctsss对方程（1）在控制容积内作积分，得 经整理，得 (3) 0)()(ewewdxxsadxdxdtxadxd0)()()(xatssdxdtadxdtapppcwwweeexasxatxatxasxaxatpcwwwweeeeppwwweeep)()()()(令 则式（3）变为：系数ae, aw分别代表了节点p, e之间以及w, p之间导热阻力的倒数，反映了节点e，w处的温度对p点温度的影响程度，具有影响系数的意义。eeeexaa)(wwwwxaa)(xasaaappwepxasbpcbtatatawweepp4.1.3 界面上当量导热系数的确定方法1.算术平均法（arithmethi

3、c mean） 设在节点p、e之间，与x成线性关系，则e可由p和e确定为：)()()()(eeeeepexxxx2.调和平均（harmonic mean） 利用界面上热流连续的原理，由fourier 定律： 又 eepepeeeeepepeexxttxttxttq)()()()()epeeettqx 因此 3.两种方法的比较 若 ，按算术平均法，当网格均分时， 表明p，e间的导热热阻主要由导热系数大的物体所决定，不符合传热学的基本原理。peeeeexxx)()()(ep22pepe按调和平均法：p, e间的热阻与调和平均方法的结果一致。调和平均法在多数场合都优于算术平均法。由此，系数ae, a

4、w可以表示为： eepeeexxx)()()(eepeeexxaa)()(wwpwewxxaa)()(4.导热系数发生阶跃性变化时，阶跃面的两种处理方式 把物性阶跃面作为控制容积的分界面； 把物性的阶跃面设置成一个节点的位置，计算阶跃面上的热流密度时，精度高于第一种方法。4.1.4 一维非稳态导热方程及其离散化 一维非稳态导热的通用控制方程为：在t, t+t时间间隔内对控制体p作积分，得 sxtxaxxattc)()(1tttppcptttwwpwwepeeetpttpppxdttssadtxttaxttattxac)()()()()()()(常用的t随 时间变化的型线可以表示为： f为0与1

5、之间的加权因子。将t表达式代入（2），得：ttfftttffttdttttttt)1 ()1 (0 xatssftssfxttaxttafxttaxttaftttxacpppcppcwwpwwepeeewwpwwepeeepppp)(1 ()()()()()()1 ()()()()()()(000000化简整理后，得 其中 xasxasfafafattfftatfftatapcppweppwwweeepp)1 ()1 ()1 ()1 ()1 (0000eeeexaa)(wwwwxaa)(txacappp)(0 xafsafafaapppwep0 取f = 0, 1,1/2可依次得显式，隐式和c

6、-n格式。在直角坐标系中 当网格划分均匀时，对无内热源、常物性的 导热问题三种格式可表示为：显式： 隐式： 000022epwpptttttatx022epwpptttttatxc-n格式： 通过von neumann分析法可证明，对源项不随时间而变的问题，当 时，格式（5）是绝对稳定的；当时，稳定的条件为 12/1 f)21 (212fxta222220000 xtttxtttatttwpewpepp4.1.5 稳定的格式导致有物理意义的解 数学上认为是稳定的初值问题的格式，未必能保证在所有的时间步长下均获得具有物理意义的解。如c-n格式。例如：例如：一块无限大平板，初始温度t0 0，被 置

7、于温度为0oc，对流换热系数为无限大的流 中，试在平板内取3个节点，两个在边界上，用式（5）来确定板的中点温度随时间变化的情况。由式（5），可得 为网格fourier数。001 2(1)(1)(1)1 2opewpppof faf af attaff2oatfx 从物理概念上，tp与前一时层的温度之比应永远为正值，并随时间的增长而趋近于0。在不同的f下， 随 的变化如图所示：图4-3。 只有全隐格式才能满足要求，任何f1的格式， 当 大于一定值后都会出现物理 上不真实的解。 0/ppttofof4.3 多维非稳态导热方程的全隐格式 4.2.1 三种正交坐标系下的全隐离散方程 1.直角坐标系：二

8、维非稳态导热方程为：在时间间隔t, t+t内对控制容积p作积分，假定控制容积的界面上的热流密度是均匀的，采用全隐格式，可得：sytyxtxttc)()(直角坐标的网格系统直角坐标的网格系统非稳态项的积分：扩散项： 0() ()ttneppptswtcdxdydtcttx yt txyttytttyxttxttdydxdtytydxdydtxtxsspsnpnnwwpwepeetttewnstttnsew )()()()()()()()()()(源项：整理，得:tyxtssdxdydtsppcnsewttt )(btatatatatassnnwweepp其中界面上的当量导热系数按调和平均法计算。

9、 eeexya/)(wwwxya/)(nnnyxa/)(sssyxa/)(yxsaaaaaappsnwep0tyxcapp)(000ppctayxsb2.圆柱轴对称坐标系 在圆柱轴对称坐标系下，非稳态导热问题的控制方程为：整理，得 srtrrrxtxttc)(1)(btatatatatassnnwweepp其中eepexrra/)(wwpwxrra/)(nnnnyxra/)(ssssyxra/)(vsaaaaaappsnwep0tvcapp)(000ppctavsbxrrrvsn)(5 . 03、极坐标系极坐标系下的控制方程：离散方程为：strrrtrrrttc)(1)(1btatatatat

10、assnnwweepp其中 都是相邻两节点间导热热阻的倒数， 具有热惯性意义， 热惯性越大，上一时层的温度对下一时层的影响越大。eeeerra/)(wwwwrra/)(nnnnrra/)(ssssrra/)(vsaaaaaappsnwep0tvcapp)(000ppctavsbrrrvsn)(5 . 00pansweaaaa,4.3 源项及边界条件的处理 4.3.1 非常数源项的线性化处理 源项是指不能包括到控制方程的非稳态项、对流项与扩散项中的所有其它各项之和。 源项的局部线性化：假定在未知量微小的变动范围内，源项s可以表示成为该未知量的线性函数，在控制容积p内，可表示为： ppctsss关

11、于源项线性化的说明1. 当源项为未知量的函数时，线性化处理比假定源项为常数更合理。因为如果s = f(t)，则把各控制容积内的s作为常数处理就是以上一次迭代计算所得的t*来计算s，这样源项对于t永远有一个滞后； 2. 线性化处理是建立线性代数方程所必须的。如果采用二阶或高阶的 多项式，则所形成的离散方程就不是代数方程；3.为保证代数方程迭代求解的收敛，要求 。0ps因为离散方程都可表示为 线性代数方程迭代求解收敛的一个充分必要条件是对角占优，即： 要求 btatanbnbppvsaapnbpnbpaa0ps4. 由代数方程迭代求解的公式 sp绝对值的大小影响到迭代过程中温度的变化速度，sp的绝

12、对值越大，系统的惯性越大，相邻两次迭代之间tp的变化越小，收敛速度下降，但有利于克服迭代过程的发散；sp的绝对值小，可使变化率加快，但容易引起发散。vsabtatpnbnbnbp例1： ，取，取s sc c = 3 = 3，s sp p = -5 = -5。例2： ，可取，可取， 例3：所以所以 ts53ts95*59cst0ps3ps*125tsc224tstttttttttdtdsss*2*2*4)24()(4()24()(*242cst*4tst 当计算区域的边界为第二、三类边界条件，可采用补充节点法和附加源项法 1.补充节点法： 对于区域离散法a，以无限大平板的第二类边界条件为例，图4

13、-9边界条件的表达式： bxqdxdt4.3.2 边界条件的处理用差分表达式表示为： （2） 为得到二阶截差的公式，采用虚拟点法，在右边设一个虚拟点m1+1，m1即可视为内节点，其一阶导数的中心差分为： （3） bmmxqtt111bmmqxtt21111为消去tm1+1，由一维、稳态、含内热源的控制方程可得在m1点的离散形式为： (4)由式（3）和（4）消去tm1+1得 (5) 其中02211111sxtttmmmxqsxxttbmm)(1112/xx对于第三类边界条件，用 代入式（2）、（5），得相应于一阶与二阶截差的节点离散方程为： （6） （7）)(1mfbtthq)1/()(111x

14、htxhttfmm)1/()()(111xhtxhsxxttfmm用控制容积平衡法，对边界节点m1的控制容积作能量平衡，得 求解上式，得 （8）控制容积平衡法在边界节点离散方程的建立中得到广泛的应用。 0111xsxttqmmbxqsxxttbmm)(111二阶精度二阶精度 当采用区域离散法b时，边界节点可看作是第一种区域离散法中当边界节点所代表的控制容积趋近于零时的极限，对右端点，由式（5）和（7）得第二类边界条件： （9） 第三类边界条件： （10） x为边界节点距第一个内节点的距离。bmmxqtt111)1/()(111xhtxhttfmm例例1 1：设有一导热型方程 ，边界条件 为x=

15、0，t=0；x=1, dt/dx = 1。试将该区域三等分，分别用区域离散法a及b求解。022tdxtd （a） (b) 采用区域离散方法a时，网格划分如图 a所示。内点上采用中心差分右端采用一阶截差时，离散方程为：t1t2t3t4t1t2t3t4t5091223tt0912234ttt3/134 tt右端点采用二阶截差时，第三式应为：这一问题的精确解为： 3/2291234 tt)(12xxeeeet若采用区域离散法b，网格划分如图b，离散方程为： 028932tt0199199432ttt02818289543ttt6/154 tt例例2 2：对于不规则边界的处理不规则边界如图4-11所示

16、，试用区域离散法a进行区域离散并列出内节点(i, j)及边界节点的离散方程。解：解：设边界与网格线交于节点1，2及3。内节点记为（i, j）。对于常数无内热源问题，可以用taylor展开法导出具有一阶截差的内节点（i, j）的离散方程： (1) 0)1()1()1()1(2,522,14ybbtbbttxaatatatjiji不规则边界的处理令 ，则有 (2) 边界节点2的离散方程可以从能量守恒观点导出： (3) jitbaaatatbtbbt,2142522)11() 1(1) 1() 1(xy /02,2321jicv 对流换热量为： (4) 从节点(i, j)到节点2的导热量为： (5)

17、 )()(1(2122222hfcvtthxcba)(21)(212,2,2,ttbaxbttxaxjijiji 为计算节点1与2及2与3之间的导热量，近似地取导热面积为 为 ，得： (6)2222122221212)(2battbxbattxbxbyb2121 (7) 将式（5）-（7）代入（2），整理得 22232312cttbfjitcbaxhtbatcbtbabtcbaxhbacbbab)1(11)1(1122222,3221222222222222222.附加源项法 对与边界相邻接的控制容积p可写出方程式： (1) 其中 ，为了在tp的代数方程中不出现未知的边界温度，需利用已知的条件

18、消去tw,为此，作如下的变换：btatatatatassnnwweeppbwwxya/)(附加源项法附加源项法 (2) (3)式（2）可化为： (4)bttatatatataapwwssnneepwp)()(yqxttyttabwpwbpww)()()(byqtatatatabssnneepp 对于第二类边界条件，qb已知，可把它与b组成一个新项： 附加源项法的实质：把 作为与边界相邻的控制容积的附加常数源项，同时令aw = 0.yxssyxyxyqsbyqadccbcb)()(,yxsaaaaaapsnewpp0yxyqb /对于第三类边界条件，qb可表示为：由fourier定律，得 (7)

19、于是， (8) )(wfbtthqwpwbbxttq)()(bwpfbwpwwfbxhttxtthttq/)(/1/)()/1)(将（8）代入（4）并归并同类项，得 式中a为所研究控制容积在边界上的传热面积，v为控制容积的体积。vxhvatstatatatxhaabwfcssnneepbwp/)(/1 )/)(/1( 对第三类边界条件： 同时令aw = 0可实现使未知边界温度不进入离散方程。bwfadcxhtvas/)(/1,bwadpxhvas/)(/11,附加源项法的实施步骤：1.计算与边界相邻的内节点的控制容积的附加源 项sc,ad及sp,ad，并将它们分别加入到该控制容积 原有的sc，

20、sp中；2.令该边界上节点的导热系数b = 0, 以使 ; 3.按常规方法建立内节点的离散方程，并在内节 点的范围内求解代数方程组； 4.获得收敛解后按fourier导热定律或newton冷却公式解出未知边界的温度。0wa5.两种处理方法的比较 附加源项法比补充节点法更为简单、有效： 有利于用统一模式来处理三类边界条件； 可以缩小计算区域； 采用补充节点法时，如果把求解代数方程的区域也限在内节点，然后通过边界节点方程不断更新边界节点上的值，并以此作为下一步迭代计算的边界条件，则附加源项法的计算时间可以比这种方法节约一个数量级。 三对角阵算法三对角阵算法 求解一维导热问题求解一维导热问题交替方向

21、隐式方法交替方向隐式方法 求解多维隐式格式求解多维隐式格式 4.4.1 三对角阵算法 (1)btatatawweepp4.4 求解离散方程的三对角阵算法及交替方向隐式算法1、thomas算法 将离散方程写作： (2) 假定共由m1个节点，当i=1时，ci = 0；i = m1时，bi = 0，首尾两个节点方程中仅由两个未知数。 thomas算法分消元和回代两步。iiiiiiidtctbta11消元的目的是把（2）变为如下形式： (3) 以（3） ci + (2), 得 111iiiiqtpt111111iiiiiiiiiiiiiqctpcdtctbtcta整理，得 与（3）相比，得 1111i

22、iiiiiiiiiiipcaqcdtpcabt1iiiiipcabp11iiiiiiipcaqcdq 要计算pi, qi，需知道pi-1, qi-1, 最终需知道p1及q1。p1及q1可由左端点的离散方 程来确定： 其中 所以1012111dtctbta001tc111/ abp 111/ adq 当消元到最后一行时，有： ,而 所以将tm1回代即可得出ti。 该算法又称为三对角阵算法（tridiagonal matrix method ，tdma算法）11111mmmmqtpt0111mmtp11mmqt4.4.2 交替方向隐式算法 adi原理：原理：在某一方向用隐式求解而使y，z方向成为显式，然后再类似地在y, z方向上各自用隐式求