计算传热学第讲扩散方程的数值解

计算传热学第讲扩散方程的数值解第1页，共83页，2023年，2月20日，星期四主要内容一维稳态问题的数值解一维非稳态问题多维非稳态问题的离散化差分方程的求解第2页，共83页，2023年，2月20日，星期四主要目的掌握用数值方法求解传热问题的整体步骤数值方法的计算机实现边界条件的处理第3页，共83页，2023年，2月20日，星期四阅读与作业阅读要求：陶文铨《数值传热学》第4章作业：P124题4－1；P125题4－7完成课外作业第一题和第二题第4页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解控制方程：其中，A(x)是面积函数。定义如下：直角坐标系：A(x)＝1（无限大平板导热问题）柱坐标系：A(x)＝x（极坐标系中的一维问题，无限长圆筒壁导热问题）球坐标系：A(x)＝x2（通过球壁的导热）变截面肋片：A(x)第5页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解4.1.1求解区域的离散化x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE图1一维问题空间区域的离散化第6页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解4.1.2源项的线性化在通常情况下，S＝S（T）线性化：

S＝Sc＋SpT

（2）其中，按负斜率源项原则，

Sp＝Sp(T*)0（3）第7页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解4.1.3控制方程的离散化将方程（1）两边通乘A（x），并对x从w到e积分：x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE图1一维问题空间区域的离散化第8页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE图1一维问题空间区域的离散化第9页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解在上面的积分过程中，我们假定：待求变量T在控制容积P上为常数整个控制容积的A(x)为常数，且等于P点的值。第10页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解将（5）和（6）代入方程（4），整理后得到，其中，第11页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解其中，下标：大写字母表示在节点处取值，小写字母表示在相应的控制面处取值第12页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解可能的改进方案：对源项积分时采用线性分布第13页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解4.1.4交界面参数的计算线性插值法（算术平均）调和平均法待求变量插值Kirchhoff变换法第14页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解4.1.5跃变界面的处理把跃变界面作为控制面调和平均法把跃变界面作为节点算术平均法Kirchhoff变换法待求变量插值法把跃变界面放置在其它位置所有方法都适用把跃变界面作为边界第15页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解把跃变界面作为边界可以考虑接触热阻rc

(W·m2)/K满足流的唯一性原则，第16页，共83页，2023年，2月20日，星期四4.1一维稳态导热问题的数值解4.1.6边界条件的处理直角坐标左边界，第二类边界条件qBx=0注意：qB的正方向与x轴的正方向一致！！第17页，共83页，2023年，2月20日，星期四ee边界条件的处理网格是用外节点法划分的边界上出现半个控制容积qBx=0123(x)1(x)2边界节点的差分方程可以用下述方法推出：一阶精度的Taylor级数展开法第18页，共83页，2023年，2月20日，星期四边界条件的处理整理后得到：特点：最简单的处理方法只有一阶精度与控制方程的精度不匹配第19页，共83页，2023年，2月20日，星期四边界条件的处理元体能量平衡法：在研究边界节点所代表的控制容积（元体）的能量平衡流入CV的能量＋内热源发出的热量＝流出CV的能量123eeqBx=0(x)1(x)21½(x)1Sq1流入CV的能量通过边界流入的热量qB通过控制面流入的热量q1内热源发出的热量＝½(x)1S流出CV的能量＝0第20页，共83页，2023年，2月20日，星期四边界条件的处理代入能量守恒关系，整理后得到，特点灵活，便于处理各种复杂的边界条件二阶精度，与内部节点精度等级匹配（请证明！！）推导过程较繁第21页，共83页，2023年，2月20日，星期四边界条件的处理控制方程法在直角坐标的条件下，方程（1）变为，123eeqBx=0(x)1(x)21假定在节点1～2之间的导热系数为常数，且恒等于e，则有，对于点e，第22页，共83页，2023年，2月20日，星期四边界条件的处理另一方面，参照附图，123eeqBx=0(x)1(x)21第23页，共83页，2023年，2月20日，星期四边界条件的处理所以，将之代入式（8）第24页，共83页，2023年，2月20日，星期四边界条件的处理整理后得到，特点二阶精度不具有一般性推导繁琐第25页，共83页，2023年，2月20日，星期四边界条件的处理二阶精度的Taylor级数展开法123eeqBx=0(x)1(x)21按二阶精度的差商公式第26页，共83页，2023年，2月20日，星期四边界条件的处理代入式（16），整理后得到，代入方程（8），整理后得到，第27页，共83页，2023年，2月20日，星期四边界条件的处理二阶精度具有一般性增加计算工作量一般很少采用第28页，共83页，2023年，2月20日，星期四边界条件的处理求解区域是用内节点法离散化的qBx=0(x)2(x)3123(x)1(x)2边界节点的控制容积或它所代表的求解区域？边界节点的控制容积＝0于是，按元体能量平衡法或其它二阶精度的方法，令与源项对应的项等于0，得到，第29页，共83页，2023年，2月20日，星期四边界条件的处理说明：尽管它与一阶Taylor级数展开法的结果形式上相同，但它却是二阶精度的！请大家证明这一结论。采用内节点法划分网格时，即使在均匀网络的前提下，第1个近边界节点也不是等步长的。qBx=0(x)2(x)3123(x)1(x)2从图中可以清楚地看出这一点即使(x)2=(x)3

(x)1也不等于

(x)2所以要对第一个内部节点给予特别注意。第30页，共83页，2023年，2月20日，星期四边界条件的处理例如，对于直角坐标系，对于节点2，qBx=0(x)2(x)3123(x)1(x)2注意：CV2的左控制面w与节点W（节点1）重合，即与左边界重合！控制面e！第31页，共83页，2023年，2月20日，星期四边界条件的处理或写成，第32页，共83页，2023年，2月20日，星期四边界条件的处理附加源项法（Additionalsourcetermmethod)以内节点法为例由方程（19）解出边界节点上的待求变量T1，代入与第1个近边界节点的差分方程（21），第33页，共83页，2023年，2月20日，星期四边界条件的处理代入与第1个近边界节点的差分方程（21），整理后得到，第34页，共83页，2023年，2月20日，星期四边界条件的处理整理后得到，或者写成，其中，第35页，共83页，2023年，2月20日，星期四边界条件的处理其中Additionalsourceterm!第36页，共83页，2023年，2月20日，星期四边界条件的处理对于第3类边界条件，也可以做类似的处理，但是这时，请大家证明，第37页，共83页，2023年，2月20日，星期四边界条件的处理附加源项法的实质边界节点消去法不仅能用于内节点网格，也能用于外节点网格实施方法：计算附加源项：Sc,ad，Sp,ad把附加源项计入该控制容积中的源项中令与边界节点对应的系数（aW）等于0第38页，共83页，2023年，2月20日，星期四特别提示边界条件的处理是传热问题数值计算最重要的环节之一元体能量平衡法的基础地位尽可能采用外节点法划分网格边界节点消去法广泛应用提高收敛速度尽可能采用均匀网格边界节点的离散化方程在形式上与内部节点的相同第39页，共83页，2023年，2月20日，星期四4.1.7差分方程的求解将前面得到的差分方程（8）改写为，或者简单的写成矩阵的形式，其中，第40页，共83页，2023年，2月20日，星期四4.1.7差分方程的求解第41页，共83页，2023年，2月20日，星期四4.1.7差分方程的求解与方程（33）对比，知，由方程（33）系数阵[A]的特殊性，通常称之为三对角方程（Tri-diagonalequation)三对角方程可以采用非常高效的追赶法（TDMA法）求解基于矩阵分解属于必须掌握的内容第42页，共83页，2023年，2月20日，星期四4.1.7差分方程的求解TDMA法Fortran源程序第43页，共83页，2023年，2月20日，星期四4.1.8计算机实现：算例求解下面的一维稳态导热问题：第44页，共83页，2023年，2月20日，星期四4.1.8计算机实现：算例求解区域的离散化：内节点法：先划分控制容积，在确定节点均匀网格：x=x将整个求解区域划分为（N-2）个控制容积，N个节点（包括2个边界节点）内部节点的差分方程第45页，共83页，2023年，2月20日，星期四4.1.8计算机实现：算例其中，第46页，共83页，2023年，2月20日，星期四4.1.8计算机实现：算例注意：采用内节点法划分网格时，近边界节点与其它内部节点不尽相同，所以必须单独考虑。123NN-1(x)w(x)e(x)e当i＝2时(x)w=½xw=W=1所以，当i＝2时第47页，共83页，2023年，2月20日，星期四4.1.8计算机实现：算例所以，当i＝2时，第48页，共83页，2023年，2月20日，星期四4.1.8计算机实现：算例当i＝3，4，。。。，N-2时，第49页，共83页，2023年，2月20日，星期四4.1.8计算机实现：算例同样，当i＝N-1时(x)e=½xe=E=N所以，当i＝N-1时123NN-1(x)w(x)e(x)e第50页，共83页，2023年，2月20日，星期四4.1.8计算机实现：算例所以，当i＝N-1时，第51页，共83页，2023年，2月20日，星期四4.1.8计算机实现：算例最后得到由（N－2）个方程构成的方程组为求解上面的方程，即可得到（N－2）个未知数，即，T2,T3,T4,…….,TN-1。第52页，共83页，2023年，2月20日，星期四4.1.8计算机实现：算例注意：上面的方程组是非线性的，必须用迭代法求解求解方法：假定一个温度分布：Ti，i＝1，2，3，。。。，N计算i，i＝1，2，3，。。。，N计算a,b,c,d用TDMA法求解方程组，得到新的温度分布：Ti’计算：Max{abs(Ti-Ti’),i＝1,2,3,……,N}判断：abs(Ti-Ti’)是否小于（精度要求）如果不能满足精度要求，令Ti＝Ti’，重复上面的计算满足精度要求：计算结束第53页，共83页，2023年，2月20日，星期四4.1.8计算机实现：算例希望大家用计算机完成上面的计算，并与下面的分析解结果比较：第54页，共83页，2023年，2月20日，星期四特别提示计算机实现的基础地位关键：掌握循环变量的使用基础：对算法清晰透彻的把握保障：细心细心再细心第55页，共83页，2023年，2月20日，星期四4.2一维非稳态导热问题与稳态问题的区别：增加了时间项数学上：常微分方程变为偏微分方程数值方法上：与稳态问题没有本质的区别控制方程:第56页，共83页，2023年，2月20日，星期四4.2.1非稳态项的处理时间坐标的离散化tk=(k-1)t,k=1,2,3,……(46)非稳态项的的离散化至少有三种方案：方案1、二阶中心差分（Centraldifference)第57页，共83页，2023年，2月20日，星期四4.2.1非稳态项的处理非稳态项的的离散化方案2、向前差分（Forwarddifference)方案3、向后差分（Backwarddifference)第58页，共83页，2023年，2月20日，星期四4.2.2控制方程的离散化参照图示的节点组，对于任意一个CV，x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE图1一维问题空间区域的离散化第59页，共83页，2023年，2月20日，星期四4.2.2控制方程的离散化将方程（49）代入方程（48）第60页，共83页，2023年，2月20日，星期四4.2.2控制方程的离散化假定节点间按线性分布，则x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE图1一维问题空间区域的离散化而按式（47c），第61页，共83页，2023年，2月20日，星期四4.2.2控制方程的离散化整理后得到，将式（51）代入方程（50），第62页，共83页，2023年，2月20日，星期四4.2.2控制方程的离散化其中，第63页，共83页，2023年，2月20日，星期四4.2.2控制方程的离散化将之与稳态问题对比发现：非稳态项的存在并没有改变差分方程的基本形式发生变化的只是：aP中多出了a0P项bP中多出了a0PTk-1P项求解方法与稳态问题一样每前进一个时间步长都要求解一个方程组第64页，共83页，2023年，2月20日，星期四4.2.3边界条件的处理处理方法与稳态问题类似元体能量平衡法：考虑CV内能随时间的变化考虑图示的左边界，按能量守恒，显然有123eeqBx=0(x)1(x)21½(x)1Sq1其中Et是CV单位时间内能的变化，显然，第65页，共83页，2023年，2月20日，星期四边界条件的处理而，将q1和Et代入方程（54），有整理后得到，第66页，共83页，2023年，2月20日，星期四边界条件的处理它与稳态问题相比，多出了内能变化项第67页，共83页，2023年，2月20日，星期四4.3.4求解过程开始，t＝0t=t+t计算有关系数形成方程组，进行求解输出结果t<tcal是STOP否第68页，共83页，2023年，2月20日，星期四特别提示非稳态问题的差分格式有显式格式（Explicitformulations）时间步进法：Time-marching不需要求解方程组程序简单，对计算机的内存要求低稳定性差隐式格式（Implicitformulations）每推进一个时间步长，都需要求解一个方程组程序复杂，要求计算机有较大的内存稳定性好第69页，共83页，2023年，2月20日，星期四4.3多维非稳态导热问题的离散化4.3.1控制方程以直角坐标系中的二维导热为例第70页，共83页，2023年，2月20日，星期四4.3多维非稳态导热问题的离散化4.3.2求解区域的离散化SNWEPsnwexy(y)n(y)s(y)-n(y)+s(x)w(x)e(x)-e(x)+w第71页，共83页，2023年，2月20日，星期四4.3多维非稳态导热问题的离散化4.3.3控制方程的离散化基本思路与一维问题完全相同对方程（59）在控制容积CV上积分，第72页，共83页，2023年，2月20日，星期四4.3.3控制方程的离散化其中SNWEPsnwexy(y)n(y)s(y)-n(y)+s(x)w(x)e(x)-e(x)+w第73页，共83页，2023年，2月20日，星期四4.3.3控制方程的离散化节点间按线性分布，则代入（63）得到，第74页，共83页，2023年，2月20日，星期四4.3.3控制方程的离散化同样SNWEPsnwexy(