2有限差分法及热传导数值计算.ppt

1、第二章 有限差分法及热传导的数值计算,本章要点：1. 着重掌握导热问题数值解法的基本思想 2. 掌握节点离散方程的建立及求解 本章难点：离散方程的建立（有限差分方程）,本章主要内容： 第一节 导热问题数值求解的基本思想 第二节 内节点离散方程的建立方法 第三节 边界节点离散方程的建立及代数方程的求解 第四节 非稳态导热问题的数值解法,求解导热问题实际上就是对导热微分方程在定解条件下的积分求解，从而获得分析解。近100年来,对大量几何形状及边界条件比较简单的问题获得了分析解,但对于工程技术中遇到的许多几何形状或边界条件复杂的导热问题,由于数学上的困难目前还无法得出其分析解.随着计算机技术的迅速发

2、展，并得到日益广泛的应用.对物理问题进行离散求解的数值方法发展得十分迅速，这些数值解法主要有以下几种：,（1）有限差分法 （2）有限元方法 （3）边界元方法, 数值解法的实质 对物理问题进行数值解法的基本思路可以概括为：把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场等，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。该方法称为数值解法。 这些离散点上被求物理量值的集合称为该物理量的数值解。,2.1 导热问题数值解法的基本思想离散化 理论解 在规定的边界条件下积分，有很大局限性； 数值解 借助计算机，前景广阔。 1

3、.有限差分法原理（连续的问题 离散的问题） 以有限差分 无限微分 无限划分 实质 达到精度 以差分代数方程 微分方程 计算机帮助 （当离散点足够多时可以满足要求）,建立控制方程及定解条件,确定节点（区域离散化）,建立节点物理量的代数方程,设立温度场的迭代初值,求解代数方程,是否收敛,解的分析,改进初场,是,否,下面先对稳态导热问题中位于计算区域内部的节点(简称内节点)介绍其离散方程的建立方法，而位于边界上的节点及非稳态导热中的非稳态项的离散将在以后讨论。 为讨论方便，把如图中的节点(m，n)及其邻点取出并放大，如图所示。,图4-3 内节点离散方程的建立,2.2 内节点离散方程的建立方法,（b）

4、,基本概念：控制单元、网格划分、节点、边界、步长等,二维矩形域内稳态、常物性的导热问题,下面以一个二维导热问题为例进行分析(有限差分法)： 把一个二维物体在X及Y方向上分别以 及 距离分割成矩形网格。则其中节点（m，n)的坐标为：X=m , Y=n , 其余节点类推。（举例） 三种基本差分格式：以节点(m，n)为例 （1）向前差分： （2）向后差分： （3）中心差分：,对无内热源、稳态、二阶导热微分方程，有： 用中心差分格式 因为： 所以：,最终得： 如果取正方形网格，即取 ，则上式为： tm+1,n+tm-1,n+tm,n+1+tm,n-1-4tm,n=0 上式说明：在导热系数为常量时，热量

5、的转移可用温度差来表达； 在稳态下，流向任何节点的热量的总和必须为零。 对于每个节点写出上式，然后联立求解方程组，即可求解。 （如边界温度已知，可逐步递推求解）,泰勒级数展开法,根据泰勒级数展开式，用节点(m,n)的温度tm,n 来表示节点(m+1,n)而温度tm+1,n 用节点(m,n)的温度tm,n来表示节点(m-1,n)的 温度tm-1,n,将上两式相加可得,将上式改写成 的表达式，有,同样可得：,根据导热问题的控制方程 ( 导热微分方程 ),若 x=y 则有,得,如图所示 边界节点 (m,n) 只能代表半个元体，若边界上有向该元体传递的热流密度为q ，据能量守恒定律对该元体有：,1.

6、边界节点离散方程的建立：,(1) 平直边界上的节点,2.3 边界节点离散方程的建立及代数方程的求解,(2) 外部角点,如图所示，二维墙角计算区域中，该节点外角点仅代表 1/4 个以 为边长的元体。假设边界上有向该元体传递的热流密度为 ，则据能量守恒定律得其热平衡式为：,(3) 内部角点,如图所示内部角点代表了 3/4 个元体，在同样的假设条件下有,讨论关于边界热流密度的三种情况：,（1）绝热边界,即令上式 即可。,（2） 值不为零,流入元体， 取正，流出元体， 取负使用上述公式,（3）对流边界,此时 ，将此表达式代入上述方程，并将此项中的 与等号前的 合并。对于 的情形有,（a）平直边界,（b

7、）外部角点,（c）内部角点,2. 代数方程的求解方法,2） 迭代法：先对要计算的场作出假设（设定初场），在迭代计算中不断予以改进，直到计算前的假定值与计算结果相差小于允许值为止的方法，称迭代计算收敛。,1） 直接解法：通过有限次运算获得精确解的方法，如：矩阵求解，高斯消元法。,迭代法目前应用较多的是：,1 ）高斯赛德尔迭代法：每次迭代计算，均是使用节点温度的最新值。 2 ）用雅可比迭代法：每次迭代计算，均用上一次迭代计算出的值。,设有一三元方程组：,其中 （ i=1,2,3 ； j=1,2,3 ）及 是已知的系数（均不为零）及常数。,采用高斯赛德尔迭代法的步骤： （1）将三元方程变形为迭式方程

8、：,（2）假设一组解（迭代初场），记为： 并代入迭代方程求得第一 次解 每次计算均用最新值代入。,（3）以新的初场 重复计算，直到相邻两次迭代值之差小于允许值，则称迭代收敛，计算终止。,判断迭代是否收敛的准则：,当有接近于零的t 时，第三个较好,迭代过程结束,说明： 1 ）对于一个代数方程组，若选用的迭代方式不合适，有可能导致发散，即称迭代过程发散；,2 ）对于常物性导热问题，组成的差分方程组，迭代公式的选择应使一个迭代变量的系数总是大于或等于该式中其他变量系数绝对值的代数和，此时，结果一定收敛。,3 ）采用热平衡法导出差分方程时，若每一个方程都选用导出该方程中心节点的温度作为迭代变量，则上述

9、条件必满足，迭代一定收敛。,这一条件数学上称主对角线占优（对角占优）；,当计算区域中出现曲线边界或倾斜的边界时，常常用阶梯形的折线来模拟真实边界，然后再用上述方法建立起边界节点的离散方程。例如，如要用数值方法确定如图4-6a所示二维区域的形状因子，显然，根据对称性我们只要考虑四分之一的计算区域即可。图4-6a中的内圆边界可以来用图4-6b所示的阶梯形的折线边界来近似。只要网格取得足够密，这种近似处理方法仍能获得相当准确的结果。处理不规则边界的更好的方法要用到坐标变换，这里不做介绍。,图4-6 不规则区域的处理,2.4 非稳态导热问题的数值解法 非稳态导热与稳态导热的主要差别在于控制方程中多了一

10、个非稳态项，而其中扩散项的离散方法与稳态导热是一样的。因此，本节讨论重点将放在非稳态项的离散以及扩散项离散时所取时间层的不同对计算带来的影响上。 1.泰勒展开法 首先以一维非稳态导热为例讨论时间空间区域的离散化。如图4-8所示，x为空间坐标，我们将计算区域划分为(N-1)等份，得到N个空间节点；为时间坐标，我们将时间坐标上的计算区域划分为(I-1)等份，得到I个时间节点。从一个时间层到下一个时间层的间隔称为时间步长。空间网格线与时间网格线的交点，如(n，i)，代表了时间空间区域中的一个节点的位置，相应的温度记为tn(i)。 非稳态项 的离散有三种不同的格式。如果将函数在节点(n，i+1)对点(

11、n，i)作泰勒展开，可有,于是有,由式(b)可得在点(n，i)处一阶导数的一种差分表示式 , 的向前差分: 类似地，将t在点(n，i-1)对点(n，i)作泰勒展开，可得 的向后差分的表达式： 如果将t在点(n，i+1)及(n，i-1)处的展开式相加，则可得一阶导数的中心差分的表达式： 在非稳态导热问题的数值计算中，非稳态项的上述三种差分格式都有人采用，本书主要采用向前差分的格式，但也简单介绍了向后差分的格式。,至此，对于形如式(3-10)所示的一维非稳态导热方程，如扩散项取中心差分，非稳态项取向前差分，则有 此式可进一步改写为 求解非稳态导热方程就是从已知的初始温度分布出发，根据边界条件依次求

12、得以后各个时间层上的温度值，式(4-14b)是对平板中各内点进行这种计算的公式。由该式可见，一旦i层上各节点的温度已知，可立即算出(i+1)时层上各内点的温度，而不必求解联立方程，因而式(4-14)所代表的计算格式称为显式差分格式。显式的优点是计算工作量小，缺点是对时间步长及空间步长有一定的限制，否则会出现不合理的结果。,如果把式(4-14a)中的扩散项也用(i+1)时层上的值来表示，则有 式中已知的是i时层的值tn(i)，而未知量有3个，因此不能直接由上式立即算出tn(i+1)之值，而必须求解(i+1)时层的一个联立方程才能得出(i十1)时层各节点的温度，因而式(4-15)称为隐式差分格式。

13、从时空坐标系中的节点(n，i+1)来看，式(4-15)的左端是非稳态项的一种向后差分。隐式格式的缺点是计算工作量大，但它对步长没有限制，不会出现解的振荡现象。 以上是将一维非稳态导热方程中的两个导数项用相应的差分表示式代替而建立差分方程的，这种方法称为泰勒展开法。,2.热平衡法 这种方法不受网格是否均分及物性是否为常数等限制，是更为一般的方法。,图4-9示出了一无限大平板的右面部，其右侧面受到周围流体的冷却，表面传热系数为h。此时边界节点N代表宽度为x/2的元体(图中有阴影线的部分)。对该元体应用能量守恒定律可得,经整理得,一维非稳态导热,式中 是以x特征长度的傅里叶数，称为网格傅里叶数， 一

14、项可作如下变化： 式中Fo及Bi分别为网格傅里叶数及网格毕渥数。于是式(4-16b)又可改写为 至此，我们可以把第三类边界条件下、厚度为2的无限大平板的数值计算问题作一归纳。由于问题的对称性，只要求解一半厚度即可。设将计算区域等分为N-1等份(N个节点，见图4-10)，节点1为绝热的对称面，节点N为对流边界，则与微分形式的数学描写相对应的离散形式为,其中式(4-20)是绝热边界的一种离散方式，在确定t1(i+1)之值时需要用到t-1(i) 。根据对称性该值等于t2(i)。这样，从已知的初始分布t0出发，利用式(4-17)及(4-19)可以依次求得第二时层、第三时层直到 i 时层上的温度值(见图

15、4-8)。至于空间步长x及时间步长的选取，原则上步长越小，计算结果越接近于精确解，但所需的计算机内存及计算时间则大大增加。此外，x及之间的关系还受到显式格式稳定性的影响。下面我们从离散方程的结构来分析，说明稳定性限制的物理意义。 式(4-17)的物理意义是很明确的。该式表明，点n上i+l时刻的温度是在该点i时刻温度的基础上计及了左右两邻点温度的影响后得出的。假如两邻点的影响保持不变，合理的情况是：i时刻点n的温度越高，则其相继时刻的温度也较高；反之，i时刻点n的温度越低，则其相继时刻的温度也较低。在差分方程中要满足这种合理性是有条件的，即式(4-17)中tn(i)前的系数必须大于或等于零。如用判别式表示，则为必须保证,否则将会出现十分不合理的情况。 式(4-21)是从一维问题显