《2有限差分法及热传导数值计算》由会员分享，可在线阅读，更多相关《2有限差分法及热传导数值计算（41页珍藏版）》请在金锄头文库上搜索。
1、第二章 有限差分法及热传导的数值计 算本章要点：本章要点：1. 1. 着重掌握导热问题数值解法的基本思想着重掌握导热问题数值解法的基本思想2. 2. 掌握掌握节点离散方程的建立及求解节点离散方程的建立及求解本章难点：本章难点：离散方程的建立（离散方程的建立（有限差分方程有限差分方程）本章主要内容：本章主要内容：第一节第一节 导热问题数值求解的基本思想导热问题数值求解的基本思想第二节第二节 内节点离散方程的建立方法内节点离散方程的建立方法第三节第三节 边界节点边界节点离散方程的建立及代数方程的求解离散方程的建立及代数方程的求解第四节第四节 非非稳态导热稳态导热问题的数值解法问题的数值解法求解导热
2、问题实际上就是对导热微分方程在定解条件下的积分求解，从而获得分析解。近100年来,对大量几何 形状及边界条件比较简单的问题获得了分析解,但对于工 程技术中遇到的许多几何形状或边界条件复杂的导热问题 ,由于数学上的困难目前还无法得出其分析解.随着计算机 技术的迅速发展，并得到日益广泛的应用.对物理问题进 行离散求解的数值方法发展得十分迅速，这些数值解法主 要有以下几种： （1）有限差分法 （2）有限元方法 （3）边界元方法 数值解法的实质 对物理问题进行数值解法的基本思路可 以概括为：把原来在时间、空间坐标系中连续的物 理量的场，如导热物体的温度场等，用有限个离散 点上的值的集合来代替，通过求解
3、按一定方法建立 起来的关于这些值的代数方程，来获得离散点上被 求物理量的值。该方法称为数值解法。 这些离散点上被求物理量值的集合称为 该物理量的数值解。 2.1 2.1 导热问题数值解法的基本思想导热问题数值解法的基本思想离散化离散化理论解 在规定的边界条件下积分，有很大局限性；数值解 借助计算机，前景广阔。 1.有限差分法原理（连续的问题 离散的问题）以有限差分 无限微分 无限划分 实质 达到精度以差分代数方程 微分方程 计算机帮助（当离散点足够多时可以满足要求）建立控制方程及定解条件确定节点（区域离散化）建立节点物理量的代数方程设立温度场的迭代初值求解代数方程是否收敛解的分析改进初场是否物
4、理问题的数值求解过程下面先对稳态导热问题中位于计算区域内部的节点(简称 内节点)介绍其离散方程的建立方法，而位于边界上的节点及 非稳态导热中的非稳态项的离散将在以后讨论。为讨论方便，把如图中的节点(m，n)及其邻点取出并放 大，如图所示。图4-3 内节点离散方程的建立2.2 2.2 内节点离散方程的建立方法内节点离散方程的建立方法（b）xynm(m,n)MN基本概念：控制单元、网格划分、节点、 边界、步长等二维矩 形域内 稳态、 常物性 的导热 问题下面以一个二维导热问题为例进行分析(有限差分法)： 把一个二维物体在X及Y方向上分别以 及 距离分割成矩 形网格。则其中节点（m，n)的坐标为：X
5、=m , Y=n ,其余节点类推。（举例）三种基本差分格式：以节点(m，n)为例 （1）向前差分：（2）向后差分：（3）中心差分：对无内热源、稳态、二阶导热微分方程，有：用中心差分格式因为：所以：最终得：如果取正方形网格，即取 ，则上式为：tm+1,n+tm-1,n+tm,n+1+tm,n-1-4tm,n=0上式说明：在导热系数为常量时，热量的转移可用温度差来表达；在稳态下，流向任何节点的热量的总和必须为零。对于每个节点写出上式，然后联立求解方程组，即可求解。（如边界温度已知，可逐步递推求解）泰勒级数展开法根据泰勒级数展开式，用节点(m,n)的温度tm,n来表示节点(m+1,n)而温度tm+1
6、,n用节点(m,n)的温度tm,n来表示节点(m-1,n)的温度tm-1,n将上两式相加可得将上式改写成 的表达式，有同样可得：表示未明确写出的 级数余项中的X 的最低阶数为2根据导热问题的控制方程 ( 导热微分方程 )若 x=y 则有 得xy如图所示 边界节点 (m,n) 只能代表半个元体，若边界上有 向该元体传递的热流密度为q ，据能量守恒定律对该元体有： 1. 边界节点离散方程的建立：(1) 平直边界上的节点2.3 2.3 边界节点边界节点离散方程的建立及代数方程的求解离散方程的建立及代数方程的求解傅里叶定律(2) 外部角点如图所示，二维墙角计算区域中，该节点外角点仅代表 1/4 个以
7、为边长的元体。假设边界上有向该元体传递 的热流密度为 ，则据能量守恒定律得其热平衡式为： xyqw(3) 内部角点如图所示内部角点代表了 3/4 个元体，在同样的假设条 件下有xyqw讨论关于边界热流密度的三种情况： （1）绝热边界 即令上式 即可。 （2） 值不为零流入元体， 取正，流出元体， 取负使 用上述公式 （3）对流边界此时 ，将此表达式代入上述方程， 并将此项中的 与等号前的 合并。 对于 的情形有（a）平直边界（b）外部角点（c）内部角点2. 代数方程的求解方法 2） 迭代法：先对要计算的场作出假设（设 定初场），在迭代计算中不断予以改进，直 到计算前的假定值与计算结果相差小于允
8、许 值为止的方法，称迭代计算收敛。1） 直接解法：通过有限次运算获得精确 解的方法，如：矩阵求解，高斯消元法。 迭代法目前应用较多的是： 1 ）高斯赛德尔迭代法：每次迭代计算，均是使用节点温度的最新值。 2 ）用雅可比迭代法：每次迭代计算，均用上一次迭代计算出的值。 设有一三元方程组： 其中 （ i=1,2,3 ； j=1,2,3 ）及 是已知的系数（均不为零）及常数。采用高斯赛德尔迭代法的步骤： （1）将三元方程变形为迭式方程： （2）假设一组解（迭代初场），记为： 并代入迭代方程求得第一 次解 每次计算均用最新值 代入。 （3）以新的初场 重复计算，直到相邻两 次迭代值之差小于允许值，则称
9、迭代收敛 ，计算终止。 判断迭代是否收敛的准则：k及k+1表示迭代次数；第k次迭代得到的最大值当有接近于零的t 时，第三个较好迭代过程结束说明： 1 ）对于一个代数方程组，若选用的迭代方 式不合适，有可能导致发散，即称迭代过程 发散； 2 ）对于常物性导热问题，组成的差分方程 组，迭代公式的选择应使一个迭代变量的系 数总是大于或等于该式中其他变量系数绝对 值的代数和，此时，结果一定收敛。 3 ）采用热平衡法导出差分方程时，若每一 个方程都选用导出该方程中心节点的温度作 为迭代变量，则上述条件必满足，迭代一定 收敛。 这一条件数学上称主对角线占优（对角占优 ）； 当计算区域中出现曲线边界或倾斜的
10、边界时，常常用阶 梯形的折线来模拟真实边界，然后再用上述方法建立起边界 节点的离散方程。例如，如要用数值方法确定如图4-6a所示 二维区域的形状因子，显然，根据对称性我们只要考虑四分 之一的计算区域即可。图4-6a中的内圆边界可以来用图4-6b 所示的阶梯形的折线边界来近似。只要网格取得足够密，这 种近似处理方法仍能获得相当准确的结果。处理不规则边界 的更好的方法要用到坐标变换，这里不做介绍。图4-6 不规则区域的处理2.4 非稳态导热问题的数值解法非稳态导热与稳态导热的主要差别在于控制方程中多了 一个非稳态项，而其中扩散项的离散方法与稳态导热是一样的 。因此，本节讨论重点将放在非稳态项的离散
11、以及扩散项离散 时所取时间层的不同对计算带来的影响上。 1.泰勒展开法首先以一维非稳态导热为例讨论时间空间区域的离散 化。如图4-8所示，x为空间坐标，我们将计算区域划分为(N- 1)等份，得到N个空间节点；为时间坐标，我们将时间坐标上 的计算区域划分为(I-1)等份，得到I个时间节点。从一个时间 层到下一个时间层的间隔称为时间步长。空间网格线与时 间网格线的交点，如(n，i)，代表了时间空间区域中的一个 节点的位置，相应的温度记为tn(i)。非稳态项 的离散有三种不同的格式。如果将函数在节 点(n，i+1)对点(n，i)作泰勒展开，可有于是有由式(b)可得在点(n，i)处一阶导数的一种差分表
12、示式 , 的向前差分: 类似地，将t在点(n，i-1)对点(n，i)作泰勒展开，可得 的向后差分的表达式： 如果将t在点(n，i+1)及(n，i-1)处的展开式相加，则可得 一阶导数的中心差分的表达式：在非稳态导热问题的数值计算中，非稳态项的上述三种 差分格式都有人采用，本书主要采用向前差分的格式，但也简 单介绍了向后差分的格式。至此，对于形如式(3-10)所示的一维非稳态导热方程，如 扩散项取中心差分，非稳态项取向前差分，则有此式可进一步改写为求解非稳态导热方程就是从已知的初始温度分布出发， 根据边界条件依次求得以后各个时间层上的温度值，式(4- 14b)是对平板中各内点进行这种计算的公式。
13、由该式可见，一 旦i层上各节点的温度已知，可立即算出(i+1)时层上各内点的 温度，而不必求解联立方程，因而式(4-14)所代表的计算格式 称为显式差分格式。显式的优点是计算工作量小，缺点是对时 间步长及空间步长有一定的限制，否则会出现不合理的结果。如果把式(4-14a)中的扩散项也用(i+1)时层上的值来表示 ，则有式中已知的是i时层的值tn(i)，而未知量有3个，因此不能 直接由上式立即算出tn(i+1)之值，而必须求解(i+1)时层的一个联 立方程才能得出(i十1)时层各节点的温度，因而式(4-15)称为 隐式差分格式。从时空坐标系中的节点(n，i+1)来看，式(4- 15)的左端是非稳
14、态项的一种向后差分。隐式格式的缺点是计 算工作量大，但它对步长没有限制，不会出现解的振荡现象。以上是将一维非稳态导热方程中的两个导数项用相应的差 分表示式代替而建立差分方程的，这种方法称为泰勒展开法。2.热平衡法 这种方法不受网格是否均分及物性是否为常数等限制，是 更为一般的方法。图4-9示出了一无限大平板的右面部， 其右侧面受到周围流体的冷却，表面传 热系数为h。此时边界节点N代表宽度 为x/2的元体(图中有阴影线的部分)。 对该元体应用能量守恒定律可得经整理得一维非稳态导热式中 是以x特征长度的傅里叶数，称为网格傅里叶数，一项可作如下变化：式中Fo及Bi分别为网格傅里叶数及网格毕渥数。于是
15、式(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 时层上的温度值(见图4-8)。至于空间步长x及时 间步长的选取，原则上步长越小，计算结果越接近于精确解 ，但所需的计算机内存及计算时间则大大增加。此外，x及 之间的关系还受到显式格式稳定性的影响。下面我们从离散方 程的结构来分析，说明稳定性限制的物理意义。式(4-17)的物理意义是很明
