第五章 导热问题的数值方法

1、5 热传导问题的数值方法5.1一维稳态导热一维稳态导热在直角坐标系下的控制方程可表示为： (5-1)式中k为导热系数，T是温度，s是单位容积的热产生率。 首先选定控制体和网格，如图5.1所示，并对方程（5-1）在所选定的控制体进行积分，即得： (5-2)图5.1 控制体和网格然后进行离散化。如果用分线段性分布来计算方程（5-2）中的微商，那么最终的方程为: (5-3)假设源项s在任一控制体中之值可以表示为温度的线性函数，即，则导出的离散化方程为： (5-4)式中 (5-5)式(5-4)就是一维稳态导热方程的离散形式，系数aE和aW分别代表了节点P与E间及W与P间导热阻力的倒数，它们的大小反映了

2、节点W和E处的温度对P点的影响程度。式中的ke和kw是控制容积中的e和w界面上的当量导热系数。进行计算时，物理参数值存储在节点的位置上。为了确定ke和kw，还需规定由节点上的物理量来计算相应界面上的量的方法。常用的方法由两种，即算术平均法与调和平均法。1、算术平均法假定k与x呈线性关系，由P与E点的导数系数确定的公式为： （5-6）2、调和平均法利用传热学的基本公式可以导出确定界面上当量导热系数的调和平均公式。控制容积中P和E的导热系数不相等，但界面上热流密度应该连续，则由Fourier定律可得： （5-7）而 则 （5-8）这就是确定界面上当量导热系数的调和平均公式，它反映了串联过程热阻的迭

3、加原则。3、两种方法比较设，网格均分时， （5-9）即P和E两点间的导热阻力为，表明此时P和E间的热阻主要是由导热系数大的物体所决定，这显然不符合传热学的基本原理。实际上，此时控制体E构成了热阻的主要部分。P和E间的热阻可表示为： （5-10）从中可以看出与调和平均一致。令 ， 则 （5-11）因此，总体上看，调和平均要比算术平均更好一些。5.2 边界条件与源项的处理式（5-4）导出的离散方程只适用于内部节点。为了对某个特定的导热问题进行计算，还应加上边界条件。传热的边界条件有三类，即（1）给定边界温度；（2）给出边界热流量；（3）通过换热系数以及周围流体温度给定边界热流量。 如果是第一类边界

4、条件，则就比较简单。但如果是第二或第三类边界，则需要对边界条件进行处理。在有限差分范围内，有两种处理方法，即对边界点补充代数方程和附加源项法。首先介绍边界点补充代数方程的方法。先讨论区域离散外节点法的情形。在边界给定热流量，如图5-2所示，即给定第二类边界条件，可表示为： （5-12）上式可化为： （5-13）即 (5-14)图5-2 边界节点控制体上式的截断误差为一阶，而内节点上如采用中心差分，则截断误差为二阶。一般希望内节点与边界点离散方程截断误差等级保持一致，否则会影响计算结果的准确度。为得出具有二阶截断误差的公式，可采用虚拟点法。即在右边界外虚设M+1点，这样节点M就可视为内节点，其一

5、阶导数就可采用中心差分，即 （5-15）为消去，由一维稳态含内热源的控制方程（5-1）可得到在 M点的离散形式，即 （5-16）从以上两式消去可得 (5-17)其中，是节点M所代表的控制体的厚度。对于第三类边界条件，以代入式(5-14)和(5-17)可得相应于一阶与二阶截断误差的节点离散方程： （5-18） （5-19）用控制容积平衡法来推导。对图5-2所示边界节点的控制体作能量平衡，即可得 （5-20）解出即得(5-17)。为了使第二和第三类边界条件的离散方程具有统一的形式，把边界上的热流密度表示成以下形式，其中为边界温度。a等于给定的热流密度(进入为+)，b=0相当于第二类边界条件；a =

6、，b=相当于第三类边界条件。，0 （5-21），0 （5-22）当区域离散化采用第二种方法（即内节点法）时，边界节点可以看成是第一种区域离散法中当边界条件所代表的控制容积厚度趋近于0时的极限，式（5-22）变为，0 （5-23）其中为边界节点与第一个内节点之间的距离。上式虽然在形式上与区域离散外节点方法具有一阶截断误差的公式一样，但它却是区域离散内节点方法中具有二阶截断误差的公式。离散化过程中源项的处理问题：如果源项是常数，则在离散方程的建立过程中不带来任何困难。当源项是新求解未知量的函数时，源项处理显得十分重要，有时甚至是数值求解成败的关键所在。目前应用较广泛的一种处理方法是把源项局部线性化

7、，即把源项写成，其中为常数，为S随T而变化的曲线在P点的斜率。源项的线性化处理需要注意以下几点：1、当源项为未知量函数时，线性化的处理比假定源项为常数更为合理。因为如果S=f(T)，那么把各个控制体中的S作为常数处理就是将上一次迭代计算得到的T*计算S，也就是说源项的计算相对于当前T的计算有一个滞后，而如果按线性化处理，就不存在这个问题。2、线性化处理是建立线性代数方程所必须的。因为如果采用更高阶的多项式，则所得到的离散化方程就不可能成为线性代数方程。3、为了保证代表方程迭代求解的收敛，要求0。因为，（下标nb表示邻点，为控制体的体积），线性代数迭代求解收敛的一个充分条件是对角占优，这就要求0

8、。4、由代数方程迭代求解的公式为 (5-24)的大小影响到迭代过程中温度的变化速度，的绝对值越大，所对应系统的惯性越大，相邻两次迭代之间的变化越小，因而收敛速度下降，但有利于克服迭代过程的发散。例题：设有一导热型方程，边界条件为x=0，T=0；x=1，。试将该区域三等分，分别用外节点法和内节点法求解该问题。解：（1）采用区域离散方法外节点法时，网格划分如图5-3(a)所示。内节点采用中心差分，即。右端点采用一阶截差时，离散方程为： 图5-3 两种离散方法的网格对节点： => 对： => 边界条件： (x=0) (x=1) => 右端点采用二阶截差时,，即 => 此问题精

9、确解为： 边界条件采用不同离散表达式，对计算结果会产生影响，由表5-1可知，右端点采用二阶截差时的结果远比采用一阶截差时的结果准确。 表5-1 采用外节点法时的数值解格式精确解0.22000.46480.7616一阶截差0.24770.52290.8563二阶截差0.21680.48670.7497（2）采用离散方法内节点法求解此问题。网格选取如图5-3(b)，则对于节点2，3，4和5的离散化方程分别为：数值解与精确解的比较见表5-2。表5-2 采用内节点法时的数值解格式精确解0.10850.33770.60480.7616数值解0.10340.33720.60350.77025.3 一维非稳

10、态导热一维非稳态导热控制方程可写为: (5-25)为了建立其离散化方程，在t, t+t时间间隔内对控制体P作积分，为简便起见，设与时间无关，则可得到:(5-26)为了完成上式的积分，需对上式右端项中T如何随时间而变的方式作出选择。常用的有显式、隐式和Crank-Nicolson三种，可用下式来表示： （5-27）式中f是在0与1之间的加权因子。上述积分式最后可化为: （5-28）整理上式可得: （5-29）其中，式(5-29)就是一维非稳态导热的离散化方程，取f=0，1，以及1/2可依次得到显式、隐式以及C-N格式。在直角坐标系中，当网格均匀时，无内热源、常物性导热问题有三种格式，它们分别为：

11、1、显式： (5-30)2、隐式： （5-31）3、C-N格式： （5-32）采用von Neumann分析方法可以证明，对于源项不随时间而变的问题，对于式(5-29)，当0.5f1时，绝对稳定；而当0f<1时，稳定的条件则为，其中。当采用显式时，由式（5-30）可知，TP并不与TE和TW等其它未知数有关。任何f0的格式都将是隐式，TP必定和未知数TE和TW有关，并且必须求解一组联立方程。相比之下，显式就比较方便，但这一点将被它局限性所抵消。因为式(5-30)中的系数可能为负。事实上，为了使这个系数为正，必须使时间步长小到足以使大于。对于均匀导热，均匀网格，可以表达为：。从中可以看出，当

12、为了提高计算精度而减小时，只能采用更小的。至此，可以总结出保持离散化方程计算收敛和不失真的四个基本规则:1、控制容积界面上的连续性原则。在两个控制容积的离散化方程中，通过界面的流量必须用相同的表达式来表示。如对于扩散项，由于取了分段线性分布假设，因而保证了在控制容积界面上热量和热通量的连续性。如果如图5-4所示的二次曲线来计算界面上的热通量（-kdT/dx），则在公共界面上会造成梯度的不连续性。图5-4 由二次曲线分布引发的热通量不连续示意图1-右边斜率；2-左边斜率2、系数为正原则。某个网格节点处的因变量值只是通过对流和扩散的过程才受到相邻节点上值的影响。如对于一维导热问题，如果TE的增加必

13、然导致TP的增加，则要求TE的系数aE和TP的系数aP具有相同的符号。换言之，中心系数与相邻系数的符号必须相同。3、源项的负斜率线性化原则。由式（5-5）可知，当SP为正值时，中心节点的系数仍有可能为负值。只有当SP不为正值时，才能完全避免这种情况的出现。因此，当源项线性化处理时，。从计算方法上讲，保持SP负值，使之不致产生计算过程中的不稳定以及导致产生物理意义上的不真实解，这是至关重要的。4、相邻节点系数加和原则。基本微分方程往往只包含有因变量的导数项，于是，如果T代表因变量，则函数T与T+C（C为任一常数）两者均满足微分方程，微分方程所具有的这个特性必定要反映到与之相对应的离散化方程中。因

14、此，当TP以及所有的Tnb都增加同一常数时，离散化方程应该仍然适用，这就要求=。但当源项与T有关时，特别进行线性化处理时，相邻节点系数不服从此原则，这只是一个特例，如SP为零，仍满足此原则。需要指出的是在数学上认为是稳定的初值问题的差分格式，未必能保证在所有的时间步长下均获得具有物理意义的解。C-N格式就属于这种情况。当时，方程中的系数变成。对于均匀导热系数和均匀网格问题，可以得出该系数为。从中可以看出只要时间步长不是足够小，该系数就有可能变为负，从而有可能出现物理上不真实的解。如果要求方程中的系数永远不为负的话，只有f为1才能确保这个要求。因此，全隐身格式(f=1)能满足格式简单且在物理上不

15、失真的要求。正是因为此，我们通常采用全隐式格式。但必须承认，对小的时间步长，全隐式格式不如C-N格式精确。全隐式离散化方程可表示为: （5-33）式中, , , , 全隐式格式的主要原理是新值在整个时间步长占控制地位。5.4二维非稳态导热问题的全隐式格式在直角坐标系下二维非稳态导热方程为: （5-34）在时间间隔t, t+t内，对图5-4所示的控制容积P做积分，假设在控制容积界面上热流密度均匀，采用全隐式格式，则有：1、时间项： （5-35）2、扩散项: （5-36） 图5-5 二维直角坐标系下的网格体系3、源项: （5-37）将（5-35）、（5-36）和（5-37）代入（5-34）并整理可

16、得: （5-38）其中: , , , , , 为控制体的容积，界面上当量导热系数按调解和平均计算。轴对称的圆柱坐标中，非稳态导热问题的控制方程为: （5-39）取一弧度为中心角所包含的范围作为研究对象，网格体系见图5-6所示，采用类似的推导方法可得形式上与上式完全相同的离散方程式，即 （5-40）, , , , , , ,图5-6 二维柱坐标系下的网格体系在极坐标中, 非稳态导热问题的控制方程为： （5-41） （5-42）, , , , , ,图5-7 极坐标系下的网格体系5.5三维问题的离散化方程在z方向上再加两个相邻近节点T和B(顶和底),即可完成三维情况的离散化，离散化方程的形式与二维

17、的类似，即 （5-43）式中, , , , , , , 5.6多维导热问题边界条件的处理对于第一类边界条件的问题，由于边界温度为已知值，这样就不必补充边界节点方程，求解的区域仅限于内部节点。对于二维、三维问题的第二、三类边界条件问题，如果也能做到这样，则不仅能节省计算时间，而且有助于编程。目前比较成功的处理方法就是附加源项法。所谓附加源项法就是把由第二类或第三类边界条件所规定的进入或导出计算区域的热量作为与边界相邻的控制容积的当量源项。以直角坐标中的情形为例，如图5-8所示。与边界相邻的控制容积中的节点为P。此控制容积的离散化方程可表达为： （5-44）图5-8 二维边界控制体为了在的代数方程中不出现未知的边界温度，就需要利用已知的边界条件把消去，为此将上式作如下变换：由于其中为边界节点的导热系数，为进入该控制容积的热流密度，以进入为正。这样关于P点的方程可化为： （5-45）对于第二类边界条件，为已知，故可把它与b组成一个新项：同时， 这就是说，对第二类边界条件，如果把作为与边界相邻的控制容积的附加常数源项，记为，同时令，则所得的离散方程既