第4章-热传导问题的数值解法 (1)

1、第四章热传导问题的数值解法第四章热传导问题的数值解法4-14-1导热问题数值求解基本思想导热问题数值求解基本思想4-24-2内节点离散方程的建立内节点离散方程的建立 4-3 4-3 边界节点离散方程的建立及代数边界节点离散方程的建立及代数 方程的求解方程的求解4-4 4-4 非稳态导热问题的数值解法非稳态导热问题的数值解法1 1 、重点内容：、重点内容： 掌握导热问题数值解法的基本思路；掌握导热问题数值解法的基本思路； 利用热平衡法和泰勒级数展开法建立利用热平衡法和泰勒级数展开法建立节点的离散方程。节点的离散方程。2 2 、掌握内容：、掌握内容：数值解法的实质。数值解法的实质。 求解导热问题的

2、三种基本方法求解导热问题的三种基本方法：(1)(1)实验法实验法; (2); (2)理论分析法；理论分析法；(3)(3)数值计算法数值计算法三种方法的特点三种方法的特点实验法实验法: : 是传热学的基本研究方法。是传热学的基本研究方法。 a a 适应性不好；适应性不好； b b 费用昂贵费用昂贵分析法分析法: : a a 能获得所研究问题的精确解，可以为实验和数能获得所研究问题的精确解，可以为实验和数值计算提供比较依据；值计算提供比较依据；b b 局限性很大，对复杂的局限性很大，对复杂的问题无法求解；问题无法求解；c c 分析解具有普遍性，各种情况分析解具有普遍性，各种情况的影响清晰可见的影响

3、清晰可见 数值计算法数值计算法 有效解决复杂问题的方法；是具有一定精度的近有效解决复杂问题的方法；是具有一定精度的近似方法。在很大程度上弥补了分析法的缺点，适似方法。在很大程度上弥补了分析法的缺点，适应性强，特别对于复杂问题更显其优越性；与实应性强，特别对于复杂问题更显其优越性；与实验法相比成本低。验法相比成本低。数值解法：数值解法： 有限差分法（有限差分法（finite-differencefinite-difference） 有限元法（有限元法（finite-lementfinite-lement） 边界元法（边界元法（boundary-elementboundary-element） 分

4、子动力学模拟（分子动力学模拟（MDMD） 格子格子BoltzmannBoltzmann方法模拟（方法模拟（LBMLBM）分析解法与数值解法的异同点：分析解法与数值解法的异同点： 相同点：相同点：根本目的是相同的，即确定：根本目的是相同的，即确定： t=f(x,y,z) t=f(x,y,z) ； 热流量。热流量。 不同点：不同点：数值解法求解的是区域或时间空间坐标数值解法求解的是区域或时间空间坐标系中离散点的温度分布代替连续的温度场；分析解系中离散点的温度分布代替连续的温度场；分析解法求解的是连续的温度场的分布特征，而不是分散法求解的是连续的温度场的分布特征，而不是分散点的数值。点的数值。 对物

5、理问题进行数值解法的基本思路可以概括对物理问题进行数值解法的基本思路可以概括为：把原来在时间、空间坐标系中连续的物理量的为：把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场等，用有限个离散点上的场，如导热物体的温度场等，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理关于这些值的代数方程，来获得离散点上被求物理量的值，该方法称为量的值，该方法称为数值解法数值解法。 这些离散点上被求物理量值的集合称为该物理这些离散点上被求物理量值的集合称为该物理量的量的数值解数值解。 4-1 4-

6、1 导热问题数值求解的基本思想导热问题数值求解的基本思想4.1.1 4.1.1 基本思想基本思想建立控制方程及定解条件建立控制方程及定解条件确定节点（区域离散化）确定节点（区域离散化）建立节点物理量的代数方程建立节点物理量的代数方程设立温度场的迭代初值设立温度场的迭代初值求解代数方程求解代数方程是否收敛是否收敛解的分析解的分析改进初场改进初场是是否否4.1.2 4.1.2 物理问题的数值求解过程物理问题的数值求解过程0tyf3thf2thf1thx二维矩形域内稳态无内热二维矩形域内稳态无内热源，常物性的导热问题源，常物性的导热问题2 2 例题条件例题条件（a）（1 1）建立控制方程及定解条件）

7、建立控制方程及定解条件 控制方程（即控制方程（即导热微分方程导热微分方程） 22220ttxy0tyf3thf2thf1thx二维矩形域内无内热源、稳态、常物性的导热问二维矩形域内无内热源、稳态、常物性的导热问题采用数值解法的步骤：题采用数值解法的步骤：（2 2）区域离散化（确立节点）区域离散化（确立节点） 用一系列与坐标轴平行的网格线把求解区域用一系列与坐标轴平行的网格线把求解区域划分成若干个子区域，用网格线的交点作为划分成若干个子区域，用网格线的交点作为需要确定温度值的空间位置，称为需要确定温度值的空间位置，称为节点节点 ( ( 结结点点 ) ) ，节点的位置用该节点在两个方向上的，节点的

8、位置用该节点在两个方向上的标号标号 m m ， n n 表示。表示。 相邻两节点间的距离相邻两节点间的距离称称步长步长。xyxynm(m,n)MN（b）xyxynm(m,n)MN基本概念：网格线、节点、界面线、步长、基本概念：网格线、节点、界面线、步长、控制容积控制容积二维矩形二维矩形域内稳态域内稳态无内热源，无内热源，常物性的常物性的导热问题导热问题 （3 3）建立节点物理量的代数方程（离散方程）建立节点物理量的代数方程（离散方程） 节点上物理量的代数方程称离散方程。节点上物理量的代数方程称离散方程。 首先划分各节点的类型；首先划分各节点的类型； 其次，建立节点离散方程；其次，建立节点离散方

9、程； 最后，代数方程组的形成。最后，代数方程组的形成。 对节点对节点 (m,n) (m,n) 的代数方程，当的代数方程，当 x=x=y y 时，时，有：有： ,1,1,1,11()4m nmnmnm nm nttttt（4 4） 设立迭代初场设立迭代初场 代数方程组的求解方法有直接解法与迭代数方程组的求解方法有直接解法与迭代解法，传热问题的有限差分法中主要采用代解法，传热问题的有限差分法中主要采用迭代法。采用迭代法求解时，需对被求的温迭代法。采用迭代法求解时，需对被求的温度场预先设定一个解，这个解称为度场预先设定一个解，这个解称为初场初场，并，并在求解过程中不断改进。在求解过程中不断改进。 （

10、5 5）求解代数方程组）求解代数方程组 本例中除本例中除 m=1 m=1 的左边界上的左边界上各节点的温度已知外，其余各节点的温度已知外，其余 (M-1)N (M-1)N 个节点均需建立离散个节点均需建立离散方程，共有方程，共有 (M-1)N (M-1)N 个方程，个方程，则构成一个封闭的代数方程则构成一个封闭的代数方程组。组。xyxynm(m,n)MN 求解时遇到的问题：求解时遇到的问题： 线性；线性； 非线性；非线性； 收敛性等。收敛性等。 2 2 ）非线性代数方程组：）非线性代数方程组：代数方程一经建立，其中代数方程一经建立，其中各项系数在整个求解过程中不断更新。各项系数在整个求解过程中

11、不断更新。 3 3 ）是否收敛判断：）是否收敛判断：是指用迭代法求解代数方程是是指用迭代法求解代数方程是否收敛，即本次迭代计算所得之解与上一次迭代计否收敛，即本次迭代计算所得之解与上一次迭代计算所得之解的偏差是否小于允许值。算所得之解的偏差是否小于允许值。 1 1 ）线性代数方程组：）线性代数方程组：代数方程一经建立，其中代数方程一经建立，其中各项系数在整个求解过程中不再变化；各项系数在整个求解过程中不再变化；（6 6） 解的分析解的分析 通过求解代数方程，获得物体中的温度分布，通过求解代数方程，获得物体中的温度分布，根据温度场应进一步计算通过的热流量，热根据温度场应进一步计算通过的热流量，热

12、应力及热变形等。应力及热变形等。因此，对于数值分析计算所得的温度场及其因此，对于数值分析计算所得的温度场及其它物理量应作详细分析，以获得定性或定量它物理量应作详细分析，以获得定性或定量上的结论。上的结论。 4.2 4.2 内节点离散方程的建立方法内节点离散方程的建立方法(1) Taylor(1) Taylor（泰勒）级数展开法；（泰勒）级数展开法；(2) (2) 控制容积平衡法控制容积平衡法( (热平衡法热平衡法) )4.2.1 4.2.1 泰勒级数展开法泰勒级数展开法( )2( )( )( )()( )1!2!nnfxfxfxf xxf xxxxn 2233441,234,2624mnmnm

13、nmnmntxtxtxtttxxxxx用节点用节点(m,n)(m,n)的温度的温度t tm,nm,n来表示节点来表示节点(m-1,n)(m-1,n)的的温度温度t tm-1,nm-1,n根据泰勒级数展开式，用节点根据泰勒级数展开式，用节点( (m,nm,n) )的温度的温度t tm,nm,n来表示节点来表示节点( (m+1,nm+1,n) )的温度的温度t tm+1,nm+1,n2233441,234,2624mnmnmnmnmntxtxtxtttxxxxx将上两式相加可得将上两式相加可得24421,1,24,212mnmnm nm ntxttttxxx22,mntx将上式改写成将上式改写成

14、的表达式，有的表达式，有)(222, 1, 1,22xoxtttxtnmnmnmnm)(2221,1,22yoytttytnmnmnmnm同样可得：同样可得：表示未明确写出的表示未明确写出的级数余项中的级数余项中的X X的最低阶数为的最低阶数为2 2 根据导热问题的控制方程根据导热问题的控制方程 ( ( 导热微分方程导热微分方程 ) )1,1,1,122220mnmnmnmnmnmnttttttxy若若 x=x=y y 则有则有 ,1,1,1,11()4m nmnmnm nm nttttt22220ttxy得得一阶一阶基本思想：基本思想：对每个有限大小的控制容积应用能量守对每个有限大小的控制容

15、积应用能量守恒，从而获得温度场的代数方程组，它从基本物理恒，从而获得温度场的代数方程组，它从基本物理现象和基本定律出发，不必事先建立控制方程，依现象和基本定律出发，不必事先建立控制方程，依据能量守恒和据能量守恒和FourierFourier导热定律即可。导热定律即可。能量守恒：能量守恒：流入控制体的总热流量控制体内热源生成热流入控制体的总热流量控制体内热源生成热 流出控制体的总热流量控制体内能的增量流出控制体的总热流量控制体内能的增量ovi4.2.2 4.2.2 控制容积平衡法控制容积平衡法( (热平衡法热平衡法) )voi)(ovi 从所有方向流入控制体的净热流量从所有方向流入控制体的净热流

16、量 控制体内热源生成热控制体内能的增量控制体内热源生成热控制体内能的增量注意：上面的公式对内部节点和边界节点均适用注意：上面的公式对内部节点和边界节点均适用稳态、无内热源时：稳态、无内热源时：从所有方向流入控制体的总热流量从所有方向流入控制体的总热流量0 01,mnm nettyx ,1,m nm nnttxy ,1,m nm nsttxy 从节点通过界面传导到节点从节点通过界面传导到节点 (m,n) (m,n) 的热流量：的热流量：1,mnm nwttyx 对元体对元体 (m,n). (m,n). 根据能量守恒定律可知：根据能量守恒定律可知： 0ewns 1,mnm nwttyx 1,mnm

17、 nettyx ,1,m nm nnttxy ,1,m nm nsttxy +=0稳态、无内热源时：稳态、无内热源时：从所有方向流入控制体的总热流量从所有方向流入控制体的总热流量0 01,1,1,12222 0mnm nmnm nm nm nttttttxy化简得说明：说明： 上述分析与推导在笛卡儿坐标系中进行的；上述分析与推导在笛卡儿坐标系中进行的； 热平衡法概念清晰，过程简捷；热平衡法概念清晰，过程简捷； 热平衡法与建立微分方程的思路与过程一致，热平衡法与建立微分方程的思路与过程一致，但不同的是前者是有限大小的元体，后者是微但不同的是前者是有限大小的元体，后者是微元体。元体。 4.3 4.

18、3 边界节点离散方程的建立边界节点离散方程的建立及代数方程的求解及代数方程的求解 对于对于第一类边界条件第一类边界条件的热传导问题，处理比较简单，因的热传导问题，处理比较简单，因为已知边界的温度，可将其以数值的形式加入到内节点的离为已知边界的温度，可将其以数值的形式加入到内节点的离散方程中，组成封闭的代数方程组，直接求解。散方程中，组成封闭的代数方程组，直接求解。 对于对于第二类或第三类边界条件第二类或第三类边界条件的导热问题，所有内节点的的导热问题，所有内节点的离散方程组成的代数方程组是不封闭的，因未知边界温度，离散方程组成的代数方程组是不封闭的，因未知边界温度，因此应对边界上的节点补充相应

19、的代数方程，才能使方程组因此应对边界上的节点补充相应的代数方程，才能使方程组封闭，以便求解。封闭，以便求解。 为了求解方便，将第二类边界条件及第三类边界条件合并为了求解方便，将第二类边界条件及第三类边界条件合并起来考虑，用起来考虑，用q qww表示边界上的热流密度或热流密度表达式。表示边界上的热流密度或热流密度表达式。为使结果更具一般性，假设物体具有内热源为使结果更具一般性，假设物体具有内热源( ( 不必均匀分不必均匀分布布 ) ) 。xyqw边界节点边界节点 (m,n) (m,n) 只代表半个元体，若边界上有向该元只代表半个元体，若边界上有向该元体传递的热流密度为体传递的热流密度为q qww

20、 ，据能量守恒定律：，据能量守恒定律： 4.3.1 4.3.1 边界节点离散方程的建立边界节点离散方程的建立( (1) 1) 平直边界上的节点平直边界上的节点1,1,1,2022mnmnmnmnmnmnmnwttttxyxyttxxyyqy yx2,1,1,12124m nwm nmnm nm nx xqtttt(2) (2) 外部角点外部角点2,1,12122m nwm nmnm nx xqttt1,1,22042mnm nm nm nm nwttttyxxyx yxyq yx如图所示，二维墙角计算区域中，该节点外角点仅如图所示，二维墙角计算区域中，该节点外角点仅代表代表 1/4 1/4 个

21、以个以 为边长的元体。假设边界上有为边长的元体。假设边界上有向该元体传递的热流密度为向该元体传递的热流密度为 ，则据能量守恒定律，则据能量守恒定律得其热平衡式为：得其热平衡式为： x、 ywq(3) (3) 内部角点内部角点22,1,1,11,213(22)62wm nmnm nm nmnx qxttttt1,1,1,1,230242mnm nm nm nm nm nmnm nm nwttttttxyxxyyttyx yxyqx yx内部角点代表了内部角点代表了 3/4 3/4 个元体，在同样的假设条件下个元体，在同样的假设条件下xyqw讨论关于边界热流密度的三种情况：讨论关于边界热流密度的三

22、种情况： （1 1）绝热边界）绝热边界即令上式即令上式 即可。即可。 0wq （2 2） 值不为零值不为零wq（3 3）对流边界）对流边界此时此时 ，将此表达式代入上述方程，并，将此表达式代入上述方程，并将此项中的将此项中的 与等号前的与等号前的 合并。合并。对于对于 的情形有：的情形有：)(,nmfwtthq,m nt,m ntxy wq流入元体，流入元体， 取正，流出元体，取正，流出元体， 取负取负wq（a a）平直边界）平直边界（b b）外部角点）外部角点（c c）内部角点）内部角点2,1,1,12222m nm nmnm nm nfh xxh xttttt2,1,12212m nm n

23、mnm nfh xxh xtttt2,1,11,1322322m nm nmnm nmnm nfh xxh xtttttt4.3.2 4.3.2 处理不规则区域的阶梯型逼进法处理不规则区域的阶梯型逼进法当计算区域出现曲线边界或倾斜边界时，常常采当计算区域出现曲线边界或倾斜边界时，常常采用用阶梯形的折线阶梯形的折线来模拟真实边界，然后用上述方来模拟真实边界，然后用上述方法建立边界节点的离散方程。法建立边界节点的离散方程。4.3.3 4.3.3 代数方程的求解方法代数方程的求解方法 2 2）迭代法：）迭代法：先对要计算的场作出假设（设先对要计算的场作出假设（设定初场），在迭代计算中不断予以改进，直

24、定初场），在迭代计算中不断予以改进，直到计算前的假定值与计算结果相差小于允许到计算前的假定值与计算结果相差小于允许值为止的方法，称迭代计算收敛。值为止的方法，称迭代计算收敛。1 1）直接解法：）直接解法：通过有限次运算获得精确解通过有限次运算获得精确解的方法，如：矩阵求解，高斯消元法。的方法，如：矩阵求解，高斯消元法。 2 2 迭代法目前应用较多的是：迭代法目前应用较多的是： 1 1 ）雅可比迭代法（简单迭代）：）雅可比迭代法（简单迭代）：每次迭代每次迭代计算，均用上一次迭代计算出的值。计算，均用上一次迭代计算出的值。 2 2 ）高斯）高斯赛德尔迭代法：赛德尔迭代法：每次迭代计算，每次迭代计算

25、，均是使用节点温度的最新值。均是使用节点温度的最新值。 在计算后面的节点温度时应按下式（采用最新值）在计算后面的节点温度时应按下式（采用最新值）例如：根据第例如：根据第 k k 次迭代的数值次迭代的数值(k)n(k)2(k)1.ttt、可以求得节点温度：可以求得节点温度：)(1)(1)(212)(111) 1(1.kknnkkkbtatatat)()() 1(11) 1(22) 1(11) 1()(3)(3) 1(232) 1(131) 1(3)(2)(2)(222) 1(121) 1(2.knknnnknnnknknknkknnkkkkknnkkkbtatatatatbtatatatbtat

26、atat设有一三元方程组设有一三元方程组： 11112 213 3121122 223 3231132 233 33a ta ta tba ta ta tba ta ta tb其中其中 （ i=1,2,3 i=1,2,3 ； j=1,2,3 j=1,2,3 ）及）及 是已是已知的系数（均不为零）及常数。知的系数（均不为零）及常数。, i jaib采用高斯采用高斯赛德尔迭代法的步骤：赛德尔迭代法的步骤： （1）将三元方程变形为迭式方程：）将三元方程变形为迭式方程： 1112 213 3112221 123 3223331 132 2331()1()1()tba ta tatba ta tatba

27、 ta ta （2 2）假设一组解（迭代初场），记为）假设一组解（迭代初场），记为： 并代入迭代方程求得第一并代入迭代方程求得第一 次次 解解 每次计算均用最新值代入。每次计算均用最新值代入。 (0)(0)(0)123ttt、(1)(1)(1)123ttt、 、（3 3）以新的初场重复计算，直到相邻两次）以新的初场重复计算，直到相邻两次迭代值之差小于允许值，则称迭代收敛，迭代值之差小于允许值，则称迭代收敛，计算终止。计算终止。 判断迭代是否收敛的准则：判断迭代是否收敛的准则：)(max)() 1()()() 1()() 1(maxmaxmaxkkikikikikikikittttttttk及及

28、k+1表示迭代次数；表示迭代次数；第第k次迭代得到的最大值次迭代得到的最大值(k)maxt当有接近于零的当有接近于零的t 时，第三个较好时，第三个较好36 1010 允许的偏差；相对偏差 值一般取迭代能否收敛的判据迭代能否收敛的判据 1 1 ）对于一个代数方程组，若选用的迭代方式不合）对于一个代数方程组，若选用的迭代方式不合适，有可能导致发散，即称适，有可能导致发散，即称迭代过程发散迭代过程发散； 2 2 ）对于常物性导热问题，组成的差分方程组，迭）对于常物性导热问题，组成的差分方程组，迭代公式的选择应使一个迭代变量的系数总是大于或代公式的选择应使一个迭代变量的系数总是大于或等于该式中其他变量

29、系数绝对值的代数和，此时，等于该式中其他变量系数绝对值的代数和，此时，结果一定收敛。结果一定收敛。 3 3 ）采用）采用热平衡法热平衡法导出差分方程时，若每一个方导出差分方程时，若每一个方程都选用导出该方程中心节点的温度作为迭代变程都选用导出该方程中心节点的温度作为迭代变量，则上述条件必满足，迭代一定收敛。量，则上述条件必满足，迭代一定收敛。 121321233132112233111aaaaaaaaa，这一这一条件条件数学上称主对角线占优（对角占优）；数学上称主对角线占优（对角占优）； 4.4 4.4 非稳态导热问题的数值解法非稳态导热问题的数值解法4.4.1 4.4.1 时间时间- -空间

30、区域的离散化空间区域的离散化1 1、基本概念、基本概念如图4-8所示，x为空间坐标，为时间坐标。1）时间步长时间步长 ：指从一个时间层到下一个时间层的间隔 。2）节点节点（n, i）表示空间网格线与时间网格线的交点，即表示了时间空间区域中一个节点的位置，相应的记为： 。2 2、非稳态项的离散、非稳态项的离散非稳态项的离散有三种不同的格式格式：1）向前差分2）向后差分3）中心差分 )(int1）向前差分2）向后差分3）中心差分 (1)( ),iinnn ittt )1()(,inininttt2)1()1(,inininttt4.4.2 4.4.2 一维平板非稳态导热的显示格式一维平板非稳态导热

31、的显示格式泰勒级数展开法1）一维非稳态导热微分方程中的扩散项离散与稳态导热微分方程中的方法相同，则对一维非稳态导热微分方程中 的扩散项中心差分； 非稳态项向前差分（1）非稳态项： 采用向前差分为：（2）稳态项： 采用中心差分则为：)()1(,inininttt( )( )( )211,222iiinnnn ittttxx2)(1)()(1)() 1(2xtttattininininin则有：可改写为：显示差分与隐式差分格式显示差分与隐式差分格式求解非稳态导热微分方程，是从已知的初始温度分布出发，根据边界条件依次求得以后各个时间层上的温度值。显示差分格式显示差分格式定义定义：就是指若已知i时层上

32、各节点的温度值，根据该差分格式即可算出（i+1）时层上各内点的温度，而不必求解联立方程。即 是前一时刻（i）n节点及相邻两节点温度的显函数。优点优点：计算工作量小；缺点缺点：受时间及空间步长的限制。)(2)(1)(12) 1()21 ()(inininintxattxat4.4.3 4.4.3 非稳态导热方程的隐式格式非稳态导热方程的隐式格式隐式差分格式隐式差分格式对一维非稳态导热微分方程 中的扩散项在(i+1)时层上采用中心差分，非稳态项将t在节点(n,i+1)处对节点(n,i)采用向前差分，得定义定义：就是指已知i 时层上各节点的温度值 ，根据差分格式不能直接算出(i+1)时层上各节点的温度，而必须求解(i+1)时层上的一个联立方程组，才能算出(i+1)时层各节点的温度，此种差分格式称隐式差分格式。优点优点：不受时间及空间的步长影响；缺点缺点：计算工作量大。22xtat2) 1(1) 1() 1(1)() 1(2xtttattin