第五章导热问题的数值解ppt课件.ppt

《第五章导热问题的数值解ppt课件.ppt》由会员分享，可在线阅读，更多相关《第五章导热问题的数值解ppt课件.ppt（40页珍藏版）》请在三一办公上搜索。

1、第五章 导热问题的数值解,数值求解方法是以离散数学为基础，以计算机为工具的一种求解方法。与各种分析求解的方法相比，它在应用方面表现出很大的适用性，对于处理诸如非线性、复杂几何形状、复杂边界条件的问题以及藕合的偏微分方程组等都能较好地解决。尤其是随着计算机技术的迅速发展，数值计算的精度和速度都大大提高。因此，稍微复杂一些的导热问题现在主要依靠数值法求解。目前用于求解偏微分方程的数值方法主要有有限差分法(FDM)、有限元法(FEM) 1,2和边界元法(BEM)。有限元法起源于固体力学和结构分析。与有限差分法相比，有限元法在整个区域内单元的划分比较随意，更易于处理不规则几何形状的问题。有限差分法涉及

2、的数学基础及其表达式比较简单，本章将主要介绍有限差分法及其在导热中的应用。,5-1 导数的有限差分近似表达式,有限差分的数学基础是用差商代替微商，即用有限差分代替导数。若f(x)是连续函数，则它的导数定义为 (5-1-1) 在这里，df /dx称为微商(导数)，f/x 称为有限差商。微商是有限差商当x趋于零时的极限。在x没有达到零以前，f/x只是df/dx的近似，而两者的差值就是用差商代替微商的偏差。 利用泰勒级数可以得到一个函数的导数的有限差分近似表达式，如图5-1所示 朝前差分 (5-l-7 )朝后差分 (5-1-8),5-1 导数的有限差分近似表达式,图5-1 函数的微商与差商,5-1

3、导数的有限差分近似表达式,式中的O(x)是用差商代替微商的误差，即舍去泰勒级数的高阶项引起的误差，称为截断误差。这里O 表示数量级，O(x)表示该截断误差是与x 同量级的小量。一阶导数的中心差分表达式 (5-1-9) 注意中心差分格式的截断误差是(x)2量级的。二阶导数的中心差分格式 (5-l-10)前面介绍的一阶导数的差分格式都只用到两个点的函数值，为了得到更高精度的差分格式，可以利用更多点的函数值。 二阶导数的朝前差分： (5-1-14) 将以上方法推广，可以得到更高精度的各种有限差分表达式。,这里以二维稳态导热为例讨论有限差分法的应用。为了进行数值求解，首先要把求解的区域离散化。有限差分

4、法要求来用正交的网格，在直角坐标系中就是矩形网格，如图5-2 所示。网线的交点称为节点，在区域内的节点称为内节点，落在边界上的节点为边界节点。网线之间的距离称为步长。x 方向的步长记为x，y 方向的步长记为y。为了便于计算，需要对节点编号。若节点P的坐标为(x, y)，xix，yiy，则用(i, j )表示节点的位置，节点的温度t(x，y)则记为ti,j。数值求解的第二步是建立节点方程，即对每一个节点写出一个代数方程。建立节点方程的方法之一是从微分形式出发，也就是在微分方程或边界条件中用差商近似微商，可以得到以节点温度为未知量的代数方程。,5-2 稳态导热的数值分析,5-2 稳态导热的数值分析

5、,图5-2 差分区域的离散化、节点和步长,有内热源的二维稳态问题，在常物性条件下可由泊松方程描述： (5-2-1)区域内的所有点，包括内节点(i, j)都应满足以上的方程。把节点(i, j)处的二阶偏导数用对应的差商来近似，有 (5-2-2) (5-2 -3) 把以上两式舍去截断误差并代入式(5-2-1)，就得到内节点的差分方程 (5-2-4),5-2 稳态导热的数值分析,如果采用正方形的网格，即xy ，且无内热源(qV0) ，则式(5-2-4 ) 简化为 ( 5-2-5 ) 注意到该差分方程的截断误差是O(x)2 + (y)2。利用边界条件也可导得边界节点的差分方程。若边界节点(1, j)和

6、(n, j)分别满足绝热边界条件和第三类边界条件： (5-2-6) (5-2-7) 把式(5-l-7 )的朝前差分格式代入式(5-2-6)，可得边界节点(1，j)的差分方程 (5-2-8),5-2 稳态导热的数值分析,注意到，该节点方程的截断误差为x量级，与内节点差分方程的截断误差相比，其精度低于一个量级。为了使各节点方程的精度能够均衡，可以利用“虚节点”的概念对边界节点进行适当的处理。 假设在边界节点(l, j)的外面还有个节点(0, j )，且 。注意该处的边界还是维持绝热，而节点(1, j)可以按内节点处理，得到 (5-2-9)很显然，该节点方程的截断误差是(x)2量级的。对于第三类边界

7、条件的节点(n, j ) ，出于同样的考虑，也可以假设有一个虚节点(n+1, j )。这样，一方面边界节点(n, j )可以按内节点处理，得到 (5-2-10),5-2 稳态导热的数值分析,另一方面该节点还应满足边界条件式(5-2-14)。采用截断误差为O(x)2的中心差分格式的式(5-l-9)代入边界条件，得 (5-2-11) 联立式(5-2-10)、(5-2-11)消去虚节点的温度，整理得 (5-2-12)其他边界节点，如两个边界相交处的节点，也可以按同样的思路加以处理。,5-2 稳态导热的数值分析,建立节点方程的另一个方法是“元体热平衡法”。这一方法不是从偏微分方程的边值问题出发考虑问题

8、，而是以积分形式对一些有限的元体从能量守恒关系建立方程，因此方程的物理概念更加清晰。两者处理问题的步骤基本相同，对大多数情况也得到形式相同的差分方程。但是，从微分方程出发的方法要求温度场在节点处的函数值和一阶、二阶导数值连续，而元体热平衡法则没有这样的要求。因此，在处理位于边界上、复合介质的结合面上、有接触热阻处等温度分布不光滑、不连续的节点时，元体热平衡法就更为方便和灵活。用元体热平衡法建立节点方程，第一步仍是把区域离散化，需要把整个区域划分为有限个小单元体。在每个元体内近似地取一点的温度代表整个元体的温度，该点也称为节点。从原则上说，节点位置的选取可以是任意的，但从建立节点方程的方便考虑，

9、内部单元的节点总是取在它的中心，边界单元的节点则取在边界上。,5-2 稳态导热的数值分析,把能量守恒关系应用于每个元体，在稳态导热的情况下，从所有相邻的元体导入的热量和该元体本身的发热量之代数和应为零。这样，对于图5-3 所示的内部元体，有以下的关系： ( 5-2-13 )在根据傅里叶定律计算相邻两个元体的导热量时，假设两个节点间的温度是线性分布的。这实际上也是采用了以一阶差商近似偏导数的概念。例如： ， (5-2-14) 代入式(5-2-13)并整理，得内部元体的节点方程为与从微分方程出发得到的节点方程完全相同。,5-2 稳态导热的数值分析,5-2 稳态导热的数值分析,下面再以图5-4中位于

10、角上的单元(1, l)为例，说明怎样用元体热平衡法处理外边界从属于两种不同的边界条件时的混合边界元体。参见图5-4 ，边界元体(l, l)的一个边界处于对流换热，另一个边界受到给定热流qw的作用。用元体热平衡法可以写出 (5-2-15)如果x y ，且qV0上式可改写为 (5-2-16)各种不同组合的混合边界元体以及处于复合介质交界处的元体都可按同样的思路列出节点方程。如果微分方程和边界条件都是线性的，则得到的内节点和边界节点的差分方程都将是线性代数方程。由全部节点的差分方程构成一个线性代数方程组，其中未知数的个数与方程的个数相等。求解这一线性代数方程组就可以得到节点处的温度。,5-2 稳态导

11、热的数值分析,5-2 稳态导热的数值分析,下面介绍求解线性代数方程组的方法。从物理意义上说，物体具有稳态温度分布的条件是，单位时间里在全部边界上流出的热量应等于物体内部发出的热量；或者，在没有内热源的情况下，在全部边界上流出的总热量应等于零。如果在全部边界上给定热流边界条件，则稳态导热问题可能无解(不满足上述条件时)，或温度场的解不确定(满足上述条件时)。所以，全部给定第一类边界条件的稳态导热的提法是不充分的。线性代数方程组的数值解法可分为直接法和迭代法两类。直接法是指在没有舍入误差的条件下经过有限次的运算即可得到方程组的精确解的方法。高斯消元法和系数矩阵求逆的方法就是常用的直接法。对于阶数不

12、是很高的方程组，采用直接法是很有效的。所谓迭代法，是把求解方程组的问题化为构造一个无限序列，来逐步逼近所求的精确解，因而在有限步的迭代中将只能得到方程组的近似解。注意到用有限差分法求解稳态导热问题时，不管总节点数的多少，每个节点方程中未知数的个数对二维问题不超过5个，对三维问题不超过7个。因此，产生的代数方程组的系数知阵中包含有大量的零元素，这样的系数矩阵称为稀疏矩阵。对于稳态导热的有限差分法得到的线性方程组，用迭代法求解时收敛较快，计算程序也比较简单，因此被广泛采用。,5-2 稳态导热的数值分析,5-2 稳态导热的数值分析,用迭代法求解线性代数方程组的基本思路如下。有线性代数方程组 (5-2

13、-17) 可以把其中任一方程改写为 (5-2-18),5-2 稳态导热的数值分析,任意给定各节点的温度值 作为解的初次近似值。把这些近似值代入式(5-2-18 )的右端，可计算得到解的第一次近似 值 。把第一次近似值再次代入式(5-2-18)的右端，可得到解的第二次近似值。一般地说，在得到解的第k 次近似值后，可由式(5-2- 18 )得到第k1次近似值： (5-2-19) 在实际计算中，当然只能迭代有限的次数，而取迭代的结果 作为方程组的近似解。通常以相邻的两次迭代值之间的差值不大于预先设定的小量作为结束迭代的判据，即,5-2 稳态导热的数值分析,以上介绍的迭代过程称为简单迭代，或同步迭代。

14、它在计算k1 次迭代值时全部采用各节点的第k 次近似值。实际上，在进行第kl 次迭代的过程中，有一部分先计算的节点的温度已经得到了第k十1 次近似值。显然它们优于第k次近似值，应优先于第k次近似值而被采用。因此，可对式(5-2-19)进行修改，得 (5-2-21) 这种迭代过程称为异步迭代，或高斯一赛德尔迭代。与同步迭代相比，它的收敛速度加快。此外，由于异步迭代不需要用两套工作单元存放节点温度的老值和新值，而可以只用一套工作单元，因此可以节省计算机的存储单元。,5-2 稳态导热的数值分析,当节点的个数很多时，通常收敛速度较慢，需要的迭代次数很大。为此提出了各种改进收敛速度的措施。其中之一是在高

15、斯一赛德尔迭代基础上的超松弛迭代法。在每一轮迭代中，首先按高斯一赛德尔迭代计算得到节点温度的新值 ，然后再用一个适当选取的松弛因子来改善这一结果，即 (5-2-22)当1 时，就相当于普通的高斯赛德尔迭代；0l 称为欠松弛，不利于收敛；l 称为超松弛。对于特定的问题，可以找出一个收敛速度最快的的值，称为最佳松弛因子*。最佳松弛因子的范围为l *2 。在处理节点个数很多的问题时，另一个加快收敛的方法是先采用大网格进行粗略的计算。用计算结果对小网格的节点温度进行插值，作为迭代的初值。由于选取的迭代初值比较接近精确解，迭代过程能大大加快。,5-3 非稳态导热的数值分析,首先，求解区域的离散化不仅涉及

16、空间坐标的离散化，而且需要将时间离散化。例如，对于一维的温度响应t = t(x, )，取适当的空间步长x 和时间步长，数值分析就是把一个求连续函数t = t(x, )的问题转化为求特定的节点上在特定的时间间隔上 的温度值 离散问题。以一维非稳态导热为例，在直角坐标系中，区域内的温度满足以下偏微分方程： (5-3-1 ) 如果采用朝前差分，p时刻对时间的偏导数近似为 ( 5-3-2 ),5-3 非稳态导热的数值分析,对坐标x的二阶偏导数用p 时刻的中心差分近似为 ( 5-3-3 )在微分方程中用有限差分代替导数，经整理可以得到内节点的朝前差分形式的节点方程 ( 5-3-4 ) 其截断误差为 。因

17、为在p+1 时刻的节点温度 可以由p 时刻的已知温度 、 和 根据式(5-3-4)直接求得，所以这样的差分格式称为显式格式。,5-3 非稳态导热的数值分析,如果p + l 时刻对时间的偏导数采用朝后差分近似为 ( 5-3-5)对坐标x 的二阶偏导数用p + 1时刻的中心差分近似为 ( 5-3-6)在微分方程中用有限差分代替导数，经整理可以得到内节点的朝后差分形式的节点方程 (5-3-7),5-3 非稳态导热的数值分析,其中， ，其截断误差也是 。现在，在一个节点方程中有3 个未知的节点温度。因此，为了确定p + l 时刻各节点的温度，需要求解一组联立代数方程。这样的差分格式称为隐式格式。如果对

18、空间变量x的偏导数的有限差分取显式格式式(5-3-3 )和隐式格式式(5-3-6 )的算术平均值来表示，则得到另一种差分格式，称为克兰克-尼科尔森(Crank -Nicolson )格式，或六点差分格式。内节点的差分方程变为 (5-3-8) 或利用以上定义的符号f，上式可整理为 (5-3-9 ),5-3 非稳态导热的数值分析,在这一差分格式中，对时间的差分相当于在 时刻的中心差分，因此这一差分格式的截断误差是 ，精度高于朝前差分和朝后差分格式。对于边界节点，同样地利用元体热平衡法建立节点方程更加方便。注意，与稳态导热问题不同，在热平衡关系中增加了一项因元体温度改变而引起的内能的变化。例如，对于

19、图5-5 中给出的一维问题的边界节点n，满足对流边界条件。 (5-3-10) 如果温度对时间的变化率采用朝前差分计算，则对于单位面积的边界元体的热平衡方程可写作 (5-3-11),5-3 非稳态导热的数值分析,图5-5 时空区域的离散化,5-3 非稳态导热的数值分析,整理得 (5-3-12) 其中， ， ，这一节点方程是显式格式。如果采用朝后差分格式或克兰克-尼科尔森格式，则可导得同一个边界元体的节点方程分别为 (5-3-13) (5-3-14)以上两个节点方程都属隐式格式。其他边界条件下的边界节点，二维问题的边界节点、二维混合边界节点都可以按同样的思路导得不同差分格式的节点方程。,5-3 非

20、稳态导热的数值分析,5-3-1 差分格式的稳定性用有限差分法求解非稳态导热问题时存在差分格式稳定性的问题。用有限差分法求导热问题的数值解时，存在两类误差。其一是用有限差分近似导数时带来的截断误差。减小差分步长可以减小截断误差。此外，数值计算只能进行到有限位有效数字，在每一步计算中都会带来舍入误差。因此，计算得到的代数方程组的解通常不可能是它的精确解，而差分方程组的精确解一般地也不等同于微分方程的精确解。对于非稳态导热问题(较一般地说，对于求解抛物线型的“扩散方程”) ，由于问题本身的特性，在数值求解时存在解的稳定性问题。数值求解非稳态导热问题的过程是随时间逐步推进的。即：先由初始温度分布算出时

21、刻的温度分布，再算出2时刻的温度分布，然后渐次推进。这就有一个以前各时刻计算的温度产生误差对后面的计算有什么影响的问题。如果误差趋向于不断积累和放大，就会引起解的不稳定性。因此，保证解的稳定性在实际计算中是十分重要的。,5-3 非稳态导热的数值分析,这里仅介绍一种称为傅里叶级数法的比较简单的方法。这种方法的不足之处是仅适用于线性齐次的差分方程，因此不能给出边界条件对稳定性的影响。以一维问题为例，如果在时刻计算得到的温度分布为 (5-3-15) 其中 (x, )为精确的温度分布，以(x, )为温度分布的误差。当然，若 = 0，(x, 0)可以看作初始温度分布的误差，由前面章节的分离变量法的分析中

22、知，假设函数在0，区间上逐段可积，就可以在该区间上用特征函数展开成广义傅里叶级数，级数的每一项都满足一组特定的边界条件。作为一个例子，假定傅里叶级数有以下的形式： (5-3-16),5-3 非稳态导热的数值分析,由于微分方程是线性的，级数中每一项对其后各时刻解的影响可以叠加，因此在这里只需研究级数中的任一项。以下去掉求和符号和下标m，且常数A 也可不加考虑。实际上，这就是可以简化地考虑一个简谐波形式的分布误差随计算过程的传递。进行离散化处理，则节点i, i+1, i-1 处的误差分别为 (5-3-17) 经过一次步进运算，到时刻，误差分布变为 。记 (5-3-18) 称为误差的放大倍数。如果某

23、种差分格式在运算过程中能始终保证 ，则误差分布在计算过程中不会被放大，这种差分格式就是稳定的；反之，若 ，则差分格式就是不稳定的。,5-3 非稳态导热的数值分析,由于我们讨论的差分方程是线性齐次的，而计算的温度分布看成是精确的温度分布和分布误差的叠加，因而和两者应各自满足差分方程。对于内节点的朝前差分格式，即式(5-3-4) ，有 (5-3-19) 把的节点值式(5-3-17 )代入，可得 (5-3-20),5-3 非稳态导热的数值分析,根据以上讨论的稳定性条件，为了使显式差分方程(5-3-4 )稳定，需有 (5-3-22) 这一稳定性条件意味着在选定空间步长x以后，时间步长不能随意选取，必须

24、满足式(5-3-22)规定的条件，否则差分格式将不稳定。同样地，对于隐式的朝前差分格式式(5-3-7)可以导得其误差放大倍数 (5-3-23) 对于任何f 0，很明显 的稳定性判据总是成立的。因此隐式差分格式式(5-3-7 )是无条件稳定的。,5-3 非稳态导热的数值分析,对于克兰克-尼科尔森差分格式，可以导得 (5-3-24)可知，对于任何f 0 , 的稳定性判据总是成立的。因此克兰克-尼科尔森差分格式也是无条件稳定的。以上是用数学的方法，根据温度分布误差在步进运算的过程中是否被放大为判据来确定某一差分格式的稳定性。另一方面，从非稳态导热的物理概念上也可以得到一些有用的结论。还是以一维问题为

25、例，对内节点的显式格式有,5-3 非稳态导热的数值分析,该差分方程表明，i节点在p+l 时刻的温度可以用i-1、i、i+1三个节点在p 时刻的温度的加权平均来表示。要使这一计算符合非稳态导 热过程的物理意义，必须使三个节点温度 , , 的系数(即权重)都不小于零。在这一具体情况下，即要求1-2f 0，亦即f 1/2。如果不满足这一条件，即的系数小于零，那就意味着i 节点在p时刻的温度越高，在p + l 时刻它的温度就反而越低，进而在p + 2 时刻它的温度又会变得更高。这样就会引起节点温度的不稳定振荡，显然不符合导热的规律。同样地，对于边界节点的显式格式差分方程，如 (5-3-25),5-3

26、非稳态导热的数值分析,5-3-2 一维差分方程组的求解 追赶法在用隐式格式，包括朝后差分和六点差分格式求解非稳态导热问题时，所导得的代数方程组在每一次步进计算中都必须联立求解。但是对于一维非稳态导热问题，可以采用一种称为“追赶法”的更为简捷的直接求解方法。从以上推导的隐式格式的节点方程可以看到，每一个内节点的方程中只包含有3个未知温度，边界节点的方程中最多只有2个未知温度。在方程组的系数矩阵中只有主对角线及相邻的两条对角线上有非零元素，方程组可以整理成式以下的形式： (5-3-26),5-3 非稳态导热的数值分析,这样的系数矩阵称为“三对角线矩阵”，可以用一种比较简单的直接法，即追赶法求解。整

27、个求解分为两段，第一部分是“顺追赶”求系数。对于方程组中的第一个方程，可以改写为 (5-3-27) 其中 ， (5-3-28) 把式(5-3-27)代入方程组的第二个方程可以消去其中的t1，经整理得 (5-3-29),5-3 非稳态导热的数值分析,或改写为 (5-3-30) 其中 ， (5-3-31) 再把式(5-3-3)代入方程组的第三个方程可以消去其中的t2，并依此类推。一般地，可以得到 ，i = 2, 3， , n-1 (5-3-32),5-3 非稳态导热的数值分析,其中“顺追赶”求系数的公式为 ， (5-3-33 ) 最后，把tn-1 的表达式代入方程组的最后一个方程，就可以得到最后一

28、个节点的温度tn： (5-3-34) 现在，把式(5-3-34)代入式(5-3-32) ，可得tn-1。利用式(5-3-32)重复这一过程，可确定各节点的温度ti, ( i = n-2, n-3， ，2, 1 )。这一过程称为“逆追赶”。,5-3 非稳态导热的数值分析,5-3-3 多维空间的非稳态导热对于无内热源的二维非稳态导热，其导热微分方程为 (5-3-35) 在微分方程中用有限差分代替导数，可以直接写出内节点的差分方程。例如，若取x y ，且温度对时间的偏导数采用朝前差分，则得到显式格式 (5-3-36)其中 。这种差分格式是有条件地稳定的，该方程的稳定性条件是f l/ 4 。,5-3 非稳态导热的数值分析,如果温度对时间的偏导数采用朝后差分，则得到隐式格式 (5-3-37) 不论 取什么值，以上差分格式是无条件稳定的。如果温度对坐标的差分取p 时刻和p + l 时刻的值的平均值，则得到十点差分格式，相当于一维问题中的六点差分格式。其内节点的差分方程为 (5-3-38) 十点差分格式也是无条件稳定的。多维问题边界节点的差分方程可以用元体热平衡法导得，较为简捷。二维问题隐式格式的每个节点方程可以包含5个未知数，不能用追赶法求解，需要用消元法或迭代法求解。,