第四章_导热问题数值方法.doc

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间导热阻力的倒数，它们的大小反映了节点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、两种方法比较设，网格均分时， （5-9）即P和E两点间的导热阻力为，表明此时P和E间的热阻主要是由导热系数大的物体所决定，这显然不符合传热学的基本原理。

实际上，此时控制体E构成了热阻的主要部分P和E间的热阻可表示为： （5-10）从中可以看出与调和平均一致令 ， 则 （5-11）因此，总体上看，调和平均要比算术平均更好一些5.2 边界条件与源项的处理式（5-4）导出的离散方程只适用于内部节点为了对某个特定的导热问题进行计算，还应加上边界条件传热的边界条件有三类，即（1）给定边界温度；（2）给出边界热流量；（3）通过换热系数以及周围流体温度给定边界热流量 如果是第一类边界条件，则就比较简单但如果是第二或第三类边界，则需要对边界条件进行处理在有限差分范围内，有两种处理方法，即对边界点补充代数方程和附加源项法首先介绍边界点补充代数方程的方法先讨论区域离散外节点法的情形在边界给定热流量，如图5-2所示，即给定第二类边界条件，可表示为： （5-12）上式可化为： （5-13）即 (5-14)图5-2 边界节点控制体上式的截断误差为一阶，而内节点上如采用中心差分，则截断误差为二阶。

一般希望内节点与边界点离散方程截断误差等级保持一致，否则会影响计算结果的准确度为得出具有二阶截断误差的公式，可采用虚拟点法即在右边界外虚设M+1点，这样节点M就可视为内节点，其一阶导数就可采用中心差分，即 （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 =，b=相当于第三类边界条件0 （5-21），0 （5-22）当区域离散化采用第二种方法（即内节点法）时，边界节点可以看成是第一种区域离散法中当边界条件所代表的控制容积厚度趋近于0时的极限，式（5-22）变为，0 （5-23）其中为边界节点与第一个内节点之间的距离上式虽然在形式上与区域离散外节点方法具有一阶截断误差的公式一样，但它却是区域离散内节点方法中具有二阶截断误差的公式离散化过程中源项的处理问题：如果源项是常数，则在离散方程的建立过程中不带来任何困难当源项是新求解未知量的函数时，源项处理显得十分重要，有时甚至是数值求解成败的关键所在目前应用较广泛的一种处理方法是把源项局部线性化，即把源项写成，其中为常数，为S随T而变化的曲线在P点的斜率源项的线性化处理需要注意以下几点：1、当源项为未知量函数时，线性化的处理比假定源项为常数更为合理。

因为如果S=f(T)，那么把各个控制体中的S作为常数处理就是将上一次迭代计算得到的T*计算S，也就是说源项的计算相对于当前T的计算有一个滞后，而如果按线性化处理，就不存在这个问题2、线性化处理是建立线性代数方程所必须的因为如果采用更高阶的多项式，则所得到的离散化方程就不可能成为线性代数方程3、为了保证代表方程迭代求解的收敛，要求≤0因为，（下标nb表示邻点，为控制体的体积），线性代数迭代求解收敛的一个充分条件是对角占优，这就要求≤04、由代数方程迭代求解的公式为 (5-24)的大小影响到迭代过程中温度的变化速度，的绝对值越大，所对应系统的惯性越大，相邻两次迭代之间的变化越小，因而收敛速度下降，但有利于克服迭代过程的发散例题：设有一导热型方程，，边界条件为x=0，T=0；x=1，试将该区域三等分，分别用外节点法和内节点法求解该问题解：（1）采用区域离散方法外节点法时，网格划分如图5-3(a)所示内节点采用中心差分，即右端点采用一阶截差时，离散方程为： 图5-3 两种离散方法的网格对节点： => 对： => 边界条件： (x=0) (x=1) => 右端点采用二阶截差时,，即 => 此问题精确解为： 边界条件采用不同离散表达式，对计算结果会产生影响，由表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 一维非稳态导热一维非稳态导热控制方程可写为: (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格式在直角坐标系中，当网格均匀时，无内热源、常物性导热问题有三种格式，它们分别为：1、显式： (5-30)2、隐式： （5-31）3、C-N格式： （5-32）采用von Neumann分析方法可以证明，对于源项不随时间而变的问题，对于式(5-29)，当0.5≤f≤1时，绝对稳定；而当0≤f<1时，稳定的条件则为，其中当采用显式时，由式（5-30）可知，TP并不与TE和TW等其它未知数有关任何f≠0的格式都将是隐式，TP必定和未知数TE和TW有关，并且必须求解一组联立方程相比之下，显式就比较方便，但这一点将被它局限性所抵消因为式(5-30)中的系数可能为负。

事实上，为了使这个系数为正，必须使时间步长小到足以使大于对于均匀导热，均匀网格，可以表达为：≤从中可以看出，当为了提高计算精度而减小时，只能采用更小的至此，可以总结出保持离散化方程计算收敛和不失真的四个基本规则:1、控制容积界面上的连续性原则在两个控制容积的离散化方程中，通过界面的流量必须用相同的表达式来表示如对于扩散项，由于取了分段线性分布假设，因而保证了在控制容积界面上热量和热通量的连续性如果如图5-4所示的二次曲线来计算界面上的热通量（-kdT/dx），则在公共界面上会造成梯度的不连续。