第五章--导热问题的数值解

第五章导热问题的数值解数值求解方法是以离散数学为基础，以计算机为工具的一种求解方法。与各种分析求解的方法相比，它在应用方面表现出很大的适用性，对于处理诸如非线性、复杂几何形状、复杂边界条件的问题以及藕合的偏微分方程组等都能较好地解决。尤其是随着计算机技术的迅速发展，数值计算的精度和速度都大大提高。因此，稍微复杂一些的导热问题现在主要依靠数值法求解。目前用于求解偏微分方程的数值方法主要有有限差分法(FDM)、有限元法(FEM)[1,2]和边界元法(BEM)。有限元法起源于固体力学和结构分析。与有限差分法相比，有限元法在整个区域内单元的划分比较随意，更易于处理不规则几何形状的问题。有限差分法涉及的数学基础及其表达式比较简单，本章将主要介绍有限差分法及其在导热中的应用。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)00()()limlimxxdtfxxfxfdxxx()()()dffxxfxxOxdxx()()()dffxfxxOxdxx5-1导数的有限差分近似表达式图5-1函数的微商与差商5-1导数的有限差分近似表达式式中的O(△x)是用差商代替微商的误差，即舍去泰勒级数的高阶项引起的误差，称为截断误差。这里O表示数量级，O(△x)表示该截断误差是与△x同量级的小量。一阶导数的中心差分表达式(5-1-9)注意中心差分格式的截断误差是(△x)2量级的。二阶导数的中心差分格式(5-l-10)前面介绍的一阶导数的差分格式都只用到两个点的函数值，为了得到更高精度的差分格式，可以利用更多点的函数值。二阶导数的朝前差分：(5-1-14)将以上方法推广，可以得到更高精度的各种有限差分表达式。2()()()2dffxxfxxOxdxx2222()2()()()()dffxxfxfxxOxdxx22222()5()4(2)(3)()()dffxfxxfxxfxxOxdxx这里以二维稳态导热为例讨论有限差分法的应用。为了进行数值求解，首先要把求解的区域离散化。有限差分法要求来用正交的网格，在直角坐标系中就是矩形网格，如图5-2所示。网线的交点称为节点，在区域内的节点称为内节点，落在边界上的节点为边界节点。网线之间的距离称为步长。x方向的步长记为△x，y方向的步长记为△y。为了便于计算，需要对节点编号。若节点P的坐标为(x,y)，x＝i△x，y＝i△y，则用(i,j)表示节点的位置，节点的温度t(x，y)则记为ti,j。数值求解的第二步是建立节点方程，即对每一个节点写出一个代数方程。建立节点方程的方法之一是从微分形式出发，也就是在微分方程或边界条件中用差商近似微商，可以得到以节点温度为未知量的代数方程。5-2稳态导热的数值分析5-2稳态导热的数值分析图5-2差分区域的离散化、节点和步长有内热源的二维稳态问题，在常物性条件下可由泊松方程描述：(5-2-1)区域内的所有点，包括内节点(i,j)都应满足以上的方程。把节点(i,j)处的二阶偏导数用对应的差商来近似，有(5-2-2)(5-2-3)把以上两式舍去截断误差并代入式(5-2-1)，就得到内节点的差分方程(5-2-4)5-2稳态导热的数值分析22220Vqttxy21,,1,222,2()()ijijijijtttdtOxdxx2,1,,1222,2()()ijijijijtttdtOydyy1,,1,,1,,1,,22220()()ijijijijijijVijttttttqxy如果采用正方形的网格，即△x＝△y，且无内热源(qV＝0)，则式(5-2-4)简化为(5-2-5)注意到该差分方程的截断误差是O[(△x)2+(△y)2]。利用边界条件也可导得边界节点的差分方程。若边界节点(1,j)和(n,j)分别满足绝热边界条件和第三类边界条件：(5-2-6)(5-2-7)把式(5-l-7)的朝前差分格式代入式(5-2-6)，可得边界节点(1，j)的差分方程(5-2-8)5-2稳态导热的数值分析,1,1,,1,114ijijijijijttttt0,0txx,()0ftxLhttx1,2,0jjtt注意到，该节点方程的截断误差为△x量级，与内节点差分方程的截断误差相比，其精度低于一个量级。为了使各节点方程的精度能够均衡，可以利用“虚节点”的概念对边界节点进行适当的处理。假设在边界节点(l,j)的外面还有个节点(0,j)，且。注意该处的边界还是维持绝热，而节点(1,j)可以按内节点处理，得到(5-2-9)很显然，该节点方程的截断误差是(△x)2量级的。对于第三类边界条件的节点(n,j)，出于同样的考虑，也可以假设有一个虚节点(n+1,j)。这样，一方面边界节点(n,j)可以按内节点处理，得到(5-2-10)5-2稳态导热的数值分析0,2,jjtt1,2,1,1,1124jjjijtttt,1,1,,1,114njnjnjnjnjttttt另一方面该节点还应满足边界条件式(5-2-14)。采用截断误差为O[(△x)2]的中心差分格式的式(5-l-9)代入边界条件，得(5-2-11)联立式(5-2-10)、(5-2-11)消去虚节点的温度，整理得(5-2-12)其他边界节点，如两个边界相交处的节点，也可以按同样的思路加以处理。5-2稳态导热的数值分析1,1,,02njnjnjftthttx,1,,1,12224njnjnjnjfhxhxttttt建立节点方程的另一个方法是“元体热平衡法”。这一方法不是从偏微分方程的边值问题出发考虑问题，而是以积分形式对一些有限的元体从能量守恒关系建立方程，因此方程的物理概念更加清晰。两者处理问题的步骤基本相同，对大多数情况也得到形式相同的差分方程。但是，从微分方程出发的方法要求温度场在节点处的函数值和一阶、二阶导数值连续，而元体热平衡法则没有这样的要求。因此，在处理位于边界上、复合介质的结合面上、有接触热阻处等温度分布不光滑、不连续的节点时，元体热平衡法就更为方便和灵活。用元体热平衡法建立节点方程，第一步仍是把区域离散化，需要把整个区域划分为有限个小单元体。在每个元体内近似地取一点的温度代表整个元体的温度，该点也称为节点。从原则上说，节点位置的选取可以是任意的，但从建立节点方程的方便考虑，内部单元的节点总是取在它的中心，边界单元的节点则取在边界上。5-2稳态导热的数值分析把能量守恒关系应用于每个元体，在稳态导热的情况下，从所有相邻的元体导入的热量和该元体本身的发热量之代数和应为零。这样，对于图5-3所示的内部元体，有以下的关系：(5-2-13)在根据傅里叶定律计算相邻两个元体的导热量时，假设两个节点间的温度是线性分布的。这实际上也是采用了以一阶差商近似偏导数的概念。例如：，(5-2-14)代入式(5-2-13)并整理，得内部元体的节点方程为与从微分方程出发得到的节点方程完全相同。5-2稳态导热的数值分析0EWNSVQQQQqxy1,,ijijEttQyx1,,ijijWttQyx1,,1,,1,,1,,22220()()ijijijijijijVijttttttqxy5-2稳态导热的数值分析下面再以图5-4中位于角上的单元(1,l)为例，说明怎样用元体热平衡法处理外边界从属于两种不同的边界条件时的混合边界元体。参见图5-4，边界元体(l,l)的一个边界处于对流换热，另一个边界受到给定热流qw的作用。用元体热平衡法可以写出(5-2-15)如果△x＝△y，且qV＝0上式可改写为(5-2-16)各种不同组合的混合边界元体以及处于复合介质交界处的元体都可按同样的思路列出节点方程。如果微分方程和边界条件都是线性的，则得到的内节点和边界节点的差分方程都将是线性代数方程。由全部节点的差分方程构成一个线性代数方程组，其中未知数的个数与方程的个数相等。求解这一线性代数方程组就可以得到节点处的温度。5-2稳态导热的数值分析2,11,11,21,11,1()022224fwVttttyxyxxyhttqqxy1,12,11,22wffqyhyhyttttt5-2稳态导热的数值分析下面介绍求解线性代数方程组的方法。从物理意义上说，物体具有稳态温度分布的条件是，单位时间里在全部边界上流出的热量应等于物体内部发出的热量；或者，在没有内热源的情况下，在全部边界上流出的总热量应等于零。如果在全部边界上给定热流边界条件，则稳态导热问题可能无解(不满足上述条件时)，或温度场的解不确定(满足上述条件时)。所以，全部给定第一类边界条件的稳态导热的提法是不充分的。线性代数方程组的数值解法可分为直接法和迭代法两类。直接法是指在没有舍入误差的条件下经过有限次的运算即可得到方程组的精确解的方法。高斯消元法和系数矩阵求逆的方法就是常用的直接法。对于阶数不是很高的方程组，采用直接法是很有效的。所谓迭代法，是把求解方程组的问题化为构造一个无限序列，来逐步逼近所求的精确解，因而在有限步的迭代中将只能得到方程组的近似解。注意到用有限差分法求解稳态导热问题时，不管总节点数的多少，每个节点方程中未知数的个数对二维问题不超过5个，对三维问题不超过7个。因此，产生的代数方程组的系数知阵中包含有大量的零元素，这样的系数矩阵称为稀疏矩阵。对于稳态导热的有限差分法得到的线性方程组，用迭代法求解时收敛较快，计算程序也比较简单，因此被广泛采用。5-2稳态导热的数值分析5-2稳态导热的数值分析用迭代法求解线性代数方程组的基本思路如下。有线性代数方程组(5-2-17)可以把其中任一方程改写为(5-2-18)11112211211222221122nnnnnnnnnnatatatbatatatbatatatbLLLLLLL11(),1,2,,niiijjjiijitbatinaL5-2稳态导热的数值分析任意给定各节点的温度值作为解的初次近似值。把这些近似值代入式(5-2-18)的右端，可计算得到解的第一次近似值。把第一次近似值再次代入式(5-2-18)的右端，可得到解的第二次近似值。一般地说，在得到解的第k次近似值后，可由式(5-2-18)得到第k＋1次近似值：(5-2-19)在实际计算中，当然只能迭代有限的次数，而取迭代的结果作为方程组的近似解。通常以相邻的两次迭代值之间的差值不大于预先设定的小量ε作为结束迭代的判据，即(0)1,2,,itinL(1)1,2,,itinL(1)()11(),1,2,,nkhiiijjjiijitbatinaL()1,2,,kitinL()(1),(1,2,,)kkiittinL5-2稳态导热的数值分析以上介绍的迭代过程称为简单迭代，或同步迭代。它在计算k＋1次迭代值时全部采用各节点的第k次近似值。实际上，在进行第k＋l次迭代的过程中，有一部分先计算的节点的