第四章导热问题的数值求解

传热学：第四章导热问题的数值求解1第四章导热问题的数值求解随着计算机的普及应用和性能的不断改善，以及相关的数值计算方法的发展和应用程序的开发，传热学数值计算方法作为数值求解传热问题的有效工具也得到了相应的发展，利用计算机求解传热学问题愈来愈受到人们的普遍重视，而且在计算复杂传热问题中显示出它的优越性，因而成为传热学的一个重要的分支。数值传热的相关内容也很自然地成为工程类学生学习传热学课程的不可缺少的部分。为了使学生能简要地掌握传热学数值计算的基本方法，在这里我们以导热问题为例对传热学数值计算方法做一个简单的介绍。4-1导热问题数值解概述在第二章和第三章中我们对较为简单的导热问题，如一维、二维简单几何形状和边界条件的稳态导热和非稳态导热、以及一些特殊导热问题，象通过肋片的导热和忽略内阻的集总导热系统，进行了分析求解。然而对于一些更为复杂的导热问题，如复杂的几何形状和边界条件以及物性变化较大的情况，分析求解往往很复杂或者根本不可能。此时求解问题的唯一途径是利用数值分析的办法获得数值解。数值求解通常是对微分方程直接进行数值积分或者把微分方程转化为一组代数方程组再求解。这里要介绍的是后一种方法。如何实现从微分方程到代数方程的转化又可以采用不同的数学方法，如有限差分法、有限元法和边界元法等。作为一本入门的教材，这里仅向读者简要地介绍用有限差分析方法从微分方程确立代数方程的处理过程。有限差分法的基本思想是把原来在时间和空间坐标中连续变化的物理量（如温度、压力、速度和热流等），用有限个离散点上的数值集合来近似表示。有限差分的数学基础是用差商代替微商（导数），而几何意义是用函数在某区域内的平均变化率代替函数的真实变化率。在图4-1中可以看出有限差分表示的温度场与真实温度场的区别。图中用T0、T1、T2…表示连续的温度场T；Δx为步长，它将区域的x方向划分为有限个数的区域，Δx0、Δx1、Δx2…，它们可以相等，也可以不相等。当Δx相等时，T1处的真实变化率a可以用平均变化率b、c或d来表示，其中b、c和d分别表示三种不同差分格式下的温度随时间的变化率，即：b为向后差分格式xxxTxTdxdT)()(111；Δx10T3Δx2Δx3T2T1T0Δx0bdcaxT图4-1温度场的有限差分表示传热学：第四章导热问题的数值求解2c为向前差分格式xxTxTdxdT)()(111＋；d为中心差分格式xxxTxxTdxdT)()(111。这种差分格式也可以推广到高阶微商的情形。对于二阶微商的差分格式可以在一阶差分格式的基础上得出：2111212)()()(2)(xxxTxTxxTdxTd。采用这样的处理之后，反映温度场随时间、空间连续变化的微分方程就可以用反映离散点间温度线性变化规律的代数方程来表示。当利用相应的数学办法求解这些代数方程组之后，我们就能获得离散点上的温度值。这些温度值就可以近似表示温度场的连续的温度分布。从上面的分析不难看出，当我们要对导热问题进行数值求解时一定要采取三个大的步骤，即研究区域的离散化；离散点（节点）差分方程的建立；节点方程（代数方程）的求解。下面我们将导热问题的数值求解进行较为详细的讨论。4-2研究区域温度场的离散化导热问题的温度场是假设为时间和空间的连续函数，当进行数值求解时首先要做的事情是在所研究的时间和空间区域内把时间和空间分割成为有限大小的小区域，尤如地球被人为地划分为不同的地域且冠以不同的名称，时间被年、月、日和时、分、秒分割。如果在所分割的每一个时间间隔和空间区域内均用同一个温度值来表示，那么原来连续变化的温度场就被一个离散的阶跃变化的温度分布所代替。这就是连续变化的温度场离散化处理的基本思路。这里我们以一个矩形长柱体的非稳态导热过程为例来讨论区域离散化问题。如果不考虑矩形长柱体长度方向上的温度变化，那么它是一个二维非稳态导热问题，图4-2表示了长柱体矩形截面上区域离散化的情况。图中可见，对于给定的空间区域，在x方向上的步长为Δx，在y方向上的步长为Δy，用它们作为空间尺度可以将矩形区域划分成纵横交错的网格，交点称为节点。然后以节点为中心，在两个节点的中心处划分界限，定出节点的控制面积，对于三维情况则为控制体积SWENPSWENPSWENPτxy图4-2矩形长柱体截面区域离散化K-1时刻K时刻K+1时刻ΔxΔy传热学：第四章导热问题的数值求解3或控制容积，因而常在一般意义上称之为节点的控制体。控制体的形状是随着坐标系的不同而改变的，这里的控制体是一个个的矩形面积。网格的步长在每一个方向上可以均匀划分，也可以不均匀的划分。因此，选用不同的步长和不同的划分方法，可以将同一区域划分出不同大小、不同数目的控制区域，以及不同数目的节点数。显然，随着步长的不断减小，节点数目的不断增加，由节点温度表示的离散的温度场就会更加接近连续的温度场，但计算工作量也会随之增加。在时间方向上离散化的步长常用Δτ来表示，Δτ的选取也是可大可小的，也可以随时间的进程而变化。显然，无限小的时间步长Δτ亦会使得离散温度变化接近连续的温度改变，但随之而来的是相应的计算工作量的增加。4-3温度场节点方程的建立为了得出所研究区域的节点温度，必须建立相应的节点方程。建立节点方程可以采用不同的方法，为了更好地理解节点方程的物理意义和掌握节点方程的建立方法，我们采用控制体热平衡法来建立节点方程。下面我们以实例的形式介绍不同节点的节点方程的建立过程。1控制体的内节点方程控制体热平衡法建立节点方程的过程是将能量守恒方程应用于控制体，建立该节点与周围节点之间的能量平衡关系式，再利用傅立叶的导热定律，最后获得控制体节点温度与周围节点温度之间的关系式。考察图4-2中的节点P及其控制体，由能量平衡关系应有EQQQQQVNSEW＋＋＋＋，4-1式中，QW、QE、QS和QN分别为邻近节点W、E、S和N通过传导方式传给节点P的热流量；QV为单位时间控制体内热源的发热量；ΔΕ为控制体单位时间内热能的增加量。由导热傅立叶定律，在线性温度分布的假设下，时刻K周围节点传给节点P的热流量分别为：，1.)(;1.)(;1.)(;1.)(xTTyQxTTyQyTTxQyTTxQKPKNNKPKSSKPKEEKPKWW传热学：第四章导热问题的数值求解4以及控制体的发热流量1yxqQVV,（qV为内热源强度，即单位时间单位体积的内热源发热量。）控制体单位时间的内能增加量为1111yxTTcEyxTTcEKPKPKPKP－＋或，前者为时间上的向前差分，而后者为时间上向后差分。以上关系式中温度T的上标为所在时刻，下标为所在空间位置。将以上关系式一并代入方程4-1中，且假设Δx=Δy，经整理可以得出二维非稳态导热问题的内节点的两种差分格式的差分方程，即显示差分格式cqTxaTTTTxaTVKPKNKSKEKWKP)41()(221＋4-2和隐示差分格式])([11122cqTTTTTxaxaTVKPKNKSKEKWKP＋。4-3比较上面两种差分格式可以看出，显示差分格式最突出的优点是节点温度表达式的右边只涉及K时刻的节点温度值，那么只要知道这一时刻周围节点的温度值就可以求出该节点的下一时刻的温度值；而隐示差分格式却相反，温度表达式的两端都是同一K时刻的节点温度值，这就意味着必须同时计算同一时刻所有节点的温度值，即必须联立求解K时刻所有节点的差分方程组，增大计算工作量是显而易见的。虽然显示差分格式计算比较方便，但它却存在着一个缺点，即计算式中aΔτ/Δx2值必须满足一定的条件才不至于引起数值计算出现不收敛的问题，这在数值计算中称为差分格式的不稳定性。这里差分方程稳定性的条件是方程4-2中的变量T前面的系数必须大于或等于零，分析一下差分方程中的各项系数，在TPK前的系数应为0412xa－，改写成为4142xa。4-4此式称为显示差分格式的稳定性判据，从中看出时间步长和空间步长是相互制约的。为了获得较为精确的节点温度值，空间步长Δx的选择不能太小，按照稳定性判据的要求势必会使时间步长Δτ也要相应地不能太大，因而必须在增加节点数目的同时增多时间间隔，从而使计算工作量加大。与显示差分格式相反，由于隐示差分格式的节点方程中没有会使方程系数成为负值的系数项，因而不存在方程求解的不稳定性的问题。也就是说，对于隐示差分格式，无论传热学：第四章导热问题的数值求解5aΔτ/Δx2中的Δx和Δτ取什么样的数值，均不会出现数值计算结果的不收敛问题，因而是无条件稳定的。这样就使得我们能在满足一定精度的情况下尽可能地加大时间步长或空间步长，亦可以在计算过程中随意改变步长，而不必担心会造成计算结果的不收敛。2控制体的边界节点方程在数值计算中所研究的区域的边界条件是通过边界节点的节点方程来反映的，因而边界节点的差分方程的建立十分重要。这里同样采用边界节点的控制体热平衡来确立边界节点的差分方程。下面将具体讨论一些典型的边界节点的节点方程。对流换热边界条件下的节点方程。如图4-3所示的边界节点P，其控制体的热平衡关系式为EQQQQQCVNSW，4-5式中各项可以写成：（隐示格式）。（显示格式）；1212;1)(;12;12;12;111yxTTcEyxTTcEyTTQyxqQxyTTQxyTTQyxTTQKPKPKPKPKPKcVVKPKNNKPKSSKPKWW将上述关系式一并代入4-5式中，经整理可得出在二维非稳态导热过程中处于平直边界上的边界节点在对流换热边界条件下的节点差分方程（当Δx=Δy），即：对于显示差分格式有cqTxaxxaTxaxTTTxaTVKPKKNKSKWKP)241(2)2(22221＋；QWSNWPQNQSQCT∞ΔyΔxα图4-3对流换热边界节点示意图传热学：第四章导热问题的数值求解64-6对于隐示差分格式有]2)2([21112222cqTTxaxTTTxaxaxxaTVKPKKNKSKWKP＋。4-7分析公式4-6，可以得到其稳定性判据为xxa24142＋，不难看出条件要比内节点的情形更加苛刻一些。绝热边界条件下的边界节点方程。如图4-4所示的绝热边界节点P，其控制体的热平衡关系式为EQQQQVNSE，4-8式中各项可写成：（隐示格式）。（显示格式）；1212;12;12;12;111yxTTcEyxTTcEyxqQxyTTQxyTTQyxTTQKPKPKPKPVVKPKNNKPKSSKPKEE将上述关系式代入4-8式中，可以整理得到绝热边界条件下的节点方程（当Δx=Δy），即：对于显示差分格式有cqTxaTTTxaTVKPKNKSKEKP)41()2(221＋；4-9对于隐示差分格式有图4-4热绝缘边界节点示意图QESNEPQNQSΔyΔx绝热边界传热学：第四章导热问题的数值求解7])2([11122cqTTTTxaxaTVKPKNKSKEKP＋。4-10对于其它边界节点也可以采用上述的控制体热平衡的办法建立对应的节点差分方程式。从以上的节点方程中可以发现两个特别的系数项，即//2xxa和,它们分别是控制体（元体）的傅立叶数和毕欧数，其物理意义是前者表征控制体的导热性能与热储蓄性能的对比关系，反映控制体温度随时间变化的动态特性；而后者则体现控制体和环境间的换热性能与其导热性能的对比关系。另外，还要一个常数项)/(cqV,它则反映了内热源产生的热流引起的节点温度升高值。下