1、第九章 导热问题的数值解法1.重点内容： 掌握导热问题数值解法的基本思路； 利用热平衡法和泰勒级数展开法建立节点的离散方程。 2.掌握内容：数值解法的实质。3.了解内容：了解非稳态导热问题的两种差分格式及其稳定性。 由前述可知，求解导热问题实际上就是对导热微分方程在定解条件下的积分求解，从而获得分析解。但是，对于工程中几何形状及定解条件比较复杂的导热问题，从数学上目前无法得出其分析解。随着计算机技术的迅速发展，对物理问题进行离散求解的数值方法发展得十分迅速，并得到广泛应用，并形成为传热学的一个分支计算传热学（数值传热学），这些数值解法主要有以下几种： （1）有限差分法；（2）有限元方法；（3）
2、边界元方法。数值解法能解决的问题原则上是一切导热问题，特别是分析解方法无法解决的问题。如：几何形状、边界条件复杂、物性不均、多维导热问题。 分析解法与数值解法的异同点：相同点：根本目的是相同的，即确定 ； 。x,yzftx,yzgQ不同点：数值解法求解的是区域或时间空间坐标系中离散点的温度分布代替连续的温度场；分析解法求解的是连续的温度场的分布特征，而不是分散点的数值。 4-1 导热问题数值求解的基本思想及内节点离散方程的建立一.数值解法的基本概念1.实质：对物理问题进行数值解法的基本思路可以概括为：把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场等，用有限个离散点上的值的集合来
3、代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。该方法称为数值解法。 这些离散点上被求物理量值的集合称为该物理量的数值解。2.基本思路：数值解法的求解过程可用框图 4-1 表示。 由此可见：（1）物理模型简化成数学模型是基础； （2）建立节点离散方程是关键； （3）一般情况微分方程中，某一变量在某一坐标方向所需边界条件的个数等于该变量在该坐标方向最高阶导数的阶数。二.数值求解的步骤如图 4-2（a），二维矩形域内无内热源、稳态、常物性的导热问题采用数值解法的步骤如下：（1）建立控制方程及定解条件控制方程：是指描写物理问题的微分方程。针对图示的导热问题，它的
4、控制方程（即导热微分方程）为：（a）0y2tx边界条件： 时， ；0t时，HxfHxtytht，2当 时，0yfyttt010，当 时，Wfytxtht，3（2）区域离散化（确立节点）用一系列与坐标轴平行的网格线把求解区域划分成若干个子区域，用网格线的交点作为需要确定温度值的空间位置，称为节点(结点)，节点的位置用该节点在两个方向上的标号 ， 表示。mn相邻两节点间的距离称步长，计为 、 。每个节点都可以看成是以它xy为中心的一个小区域的代表，把节点代表的小区域称为元体（又叫控制容积），如图 4-2(b)。（3）建立节点物理量的代数方程（离散方程）节点上物理量的代数方程称离散方程。其过程如下：
5、 首先划分各节点的类型； 其次，建立节点离散方程；最后，代数方程组的形成。对节点（ ， ）的代数方程，当 时，有： mnxy1,1,4nmnmttt（b）（4）设立迭代初场代数方程组的求解方法有直接解法与迭代解法，传热问题的有限差分法中主要采用迭代法。采用迭代法求解时，需对被求的温度场预先设定一个解，这个解称为初场，并在求解过程中不断改进。（5）求解代数方程组如图 4-2（b），除 的左边界上各节点的温度已知外，其余 个1mNM1节点均需建立离散方程，共有 个方程，则构成一个封闭的代数方程组。NM求解时遇到的问题： 线性； 非线性； 收敛性等。线性代数方程组：代数方程一经建立，其中各项系数在整
6、个求解过程中不再变化；非线性代数方程组：代数方程一经建立，其中各项系数在整个求解过程中不断更新。是否收敛判断：是指用迭代法求解代数方程是否收敛，即本次迭代计算所得之解与上一次迭代计算所得之解的偏差是否小于允许值。关于变物性(物性为温度的函数)导热问题，建立的离散方程，四个邻点温度的系数不是常数，而是温度的函数。在迭代计算时，这些系数应不断更新，这是非线性问题。（6）解的分析通过求解代数方程，获得物体中的温度分布，根据温度场应进一步计算通过的热流量，热应力及热变形等。因此，对于数值分析计算所得的温度场及其它物理量应作详细分析，以获得定性或定量上的结论。三、稳态导热中位于计算区域内部的节点离散方程
7、的建立方法1.基本概念（1）内节点：位于计算区域内部的节点，称内节点。（2）差分格式：差商中的差分可以用向前、向后、中心差分表示的格式称差分格式。 2.基本方法方法： 泰勒级数展开法； 热平衡法。以下分述之。（1）泰勒级数展开法如图 4-3 所示，以节点（ ， ）处的二阶偏导数为例，对节点（ ，mn 1m） 及（ ， ）分别写出函数 对（ ， ）点的泰勒级数展开式：n1mnt对（ ， ）： nmnmnmnmn xtxtxtxttt ,4,3,2,1 26（c）对（ ， ）： nmnmnmnmnm xtxtxtxttt ,4,3,2,1 26（d）（a）+（b）得： nmnmnmnm xtxtt
8、t ,4,2,1, 12（e）变形为 的表示式得： nmxt,222,1,2 0xxttt nmn （f）上式是用三个离散点上的值计算二阶导数的严格表达式，其中： 称截断误nmxt,220x差，误差量级为 ，即表示未明确写出的级2x数余项中 的最低阶数为 2。在数值计算时，用三个相邻节点上的值近似表示二阶导数的表达式即可，则相应的略去 。于是得：20x（4-2,1,2xttxt nmnm1a）同理： （4-2,1,2yttyt nmnmn1b）根据导热问题的控制方程(导热微分方程) 得：0y2tx(4-2)22,1,1, ttxtt nmnmnmn 若 = ,则有：xy1,1,4nnn ttt
9、(2)热平衡法：其本质是傅里叶导热定律和能量守恒定律的体现。对每个元体，可用傅里叶导热定律写出其能量守恒的表达式。如图 4-3 所示，元体在垂直纸面方向取单位长度，通过元体界面(w,e,n,s) 所传导的热流量可以对有关的两个节点根据傅里叶定律写出：从节点（ ， ）通过界面 W 传导到节点（ ， ）的热流量为： 1mnmnxtynmw,1（g）同理：通过界面 e,n,s 传导给节点（ ， ）的热流量：xtynme,1（h）ytnmw,1,（i）ytxnmw,1,（j）对元体（ ， ）,根据能量守恒定律可知：mn（4-0snew3）其中规定：导入元体（ ， ）的热流量为正；导出元体（ ， ）的热
10、流量mmn为负。 将式（g）、（h）、（i）、（j）代入式（4-3），当 时即得式(b)。yx说明： 上述分析与推导是在笛卡儿坐标系中进行的； 热平衡法概念清晰，过程简捷； 热平衡法与 22 建立微分方程的思路与过程一致，但不同的是前者是有限大小的元体，后者是微元体。4-2 边界节点离散方程的建立及代数方程的求解对于第一类边界条件的导热问题，所有内节点的离散方程组成一个封闭的代数方程组，即可求解； 第二类或第三类边界条件的导热问题，所有内节点的离散方程组成的代数方程组是不封闭的，因未知边界温度，因而应对位于该边界上的节点补充相应的代数方程，才能使方程组封闭，以便求解。一、用热平衡法导出典型边界
11、点上的离散方程在下面的讨论中，先把第二类边界条件及第三类边界条件合并起来考虑，并以 代表边界上已知的热流密度值或热流密度表达式，用热平衡方法导出三wq类典型边界节点的离散方程，然后针对 的三种不同情况使导得的离散方程进wq一步具体化，为使结果更具一般性，假设物体具有内热源 (不必均匀分布)。1.位于平直边界上的节点如图 4-4 所示有阴影线的区域，边界节点 只能代表半个元体，设边界nm,上有向该元体传递的热流密度为 ，据能量守恒定律对该元体有：wq（4-0222 ,1,1,1 wnmnmnmnm yqxytxytxty 4a）若 时，则： （4-yx wnmnmnmn xqxttt 2241,
12、21, 4b）2.外部角点如图 4-5 所示，二维墙角计算区域中，节点 AE 均为外部角点，其特点是每个节点仅代表 1/4 个以 为边长的元体。假设边界上有向该元体传递的yx、热流密度为 ，则据能量守恒定律得其热平衡式为：wq（4-02422 ,1,1 wnmnmnm qyxyxytxty 5a）若 时，则： （4-yx wnmnmnm xxtt 221,1, 5b）3.内部角点：图 4-5 中的 F 点为内部角点，代表了 3/4 个元体，同理得：（4-02432,1,1,1, wnmnm nmnn qyxyxtty tttxt 6a）若 时，则：yx（4- wnmnmnmnnm xqxttt
13、t 23261 ,2,1,1, 6b）4.讨论有关 的三种情况：wq（1）若是绝热边界则 ，即令上式 即可。0w0w（2）若时 q则以给定的 值代入上述方程，注意：流入元体， 取正，流出元体，w wq取负。wq（3）若属对流边界则 ，将此表达式代入式（4-4）（4-6），并将此项中nmfwth,与等号前的 合并。对于 的情形，有：nmt, , yx平直边界：fnmnmnmn txhxttttxh 222 ,21,1,1, （4-7）对外角点：（4-fnmnmnm txhxttxh 212,21,1, 8）对内角点：（4-fnmnmnmnnm txhxttttxh 23232 ,2,1,1,1,
14、 9）其中无量纲数 是以网格步长 为特征长度的毕渥数，即为 ，是在xxBi对流边界条件的离散过程中引入的。二、代数方程的求解方法 1.直接解法：通过有限次运算获得精确解的方法，如：矩阵求逆、高斯消元法等。这一方法的缺点是计算所需的计算机内存较大，当代数方程的数目较多时使用不便。2.迭代法：先对要计算的场作出假设（设定初场），在迭代计算中不断予以改进，直到计算前的假定值与计算结果相差小于允许值为止的方法，称迭代计算收敛。目前应用较多的是：（1）高斯赛德尔迭代法：每次迭代计算，均使用节点温度的最新值。（2）用雅可比迭代法：每次迭代计算，均用上一次迭代计算出的值。设有一个三元方程组，记为：32312
15、21btatatt（a）其中 及 是已知的系数（设均不为零）及常数。,2;3,1, jij ,ib采用高斯赛德尔迭代法的步骤：（1）将三元方程变形为迭式方程：23133223111tabattabat （b）（2）假设一组解（迭代初场），记为： ，并代入迭代方程求03201tt、得第一次解 ，每次计算都用 的最新值代入。例如当由式（b）中的1321tt、第三式计算 时代入的是 之值。312t及（3）以计算所得之值作为初场，重复上述计算，直到相邻两次迭代值之差小于允许值，则称迭代收敛，计算终止。 三、判断迭代收敛的准则判断迭代是否收敛的准则一般有以下三种：（4-1maxkiit10a）（4-ki
16、t1max10b）（4-kiitmax110c）其中上角标 及 表示迭代次数， 为第 次迭代计算所得的计算区域中的k1ktax最大值。若计算区域中有接近于零的 时，采用式（4-10c）比较合适。说明：（1）对于一个代数方程组，若选用的迭代方式不合适，有可能导致发散，即称迭代过程发散；（2）对于常物性导热问题组成的差分方程组，迭代公式的选择应使一个迭代变量的系数总是大于或等于该式中其他变量系数绝对值的代数和，此时，结果一定收敛。这一条件数学上称主对角线占优（对角占优），即：， ，1312a 123a 132a（3）采用热平衡法导出差分方程时，若每一个方程都选用导出该方程中心节点的温度作为迭代变量
17、，则上述条件必满足，迭代一定收敛。4-3 非稳态导热问题的数值解法由前可知：非稳态导热和稳态导热二者微分方程的区别在于控制方程中多了一个非稳态项，其中扩散项的离散方法与稳态导热相同。本节重点讨论：（1）非稳态项离散的方法；（2）扩散项离散时所取时间层的不同对计算带来的影响。一、一维非稳态导热时间空间区域的离散化1.基本概念如图 4-8 所示， 为空间坐标，x为时间坐标。（1）时间步长 ：指从一个时间层到下一个时间层的时间间隔 。（2）节点（n，i）表示空间网格线与时间网格线的交点，即表示了时间空间区域中一个节点的位置，相应的记为 。int2.非稳态项的离散非稳态项的离散有三种不同的格式，向前差
18、分、向后差分、中心差分。（1）向前差分将函数 在节点（n，i+1）对点（n，i）作泰勒展开，则有：t n,in,iini ttt 21（a）于是有：ttinin,i O1（b）其中 表示余项中 的最低阶为一次。O由式（b）可得在点（n,i）处一阶导数的向前差分表示式：（4-ttinin,i111）（2）向后差分将 在节点（n,i-1）对点 (n,i) 作泰勒展开，可得 的向后差分表示t n,it式： ttinin,i1（4-12）（3）中心差分的向前差分与向后差分之和，即得 的中心差分表达式：n,it n,it（4-ttiin,i2113）二、一维非稳态导热微分方程的离散方法1.泰勒级数展开法
19、（1）显式差分格式一维非稳态导热微分方程中的扩散项离散与稳态导热微分方程中的方法相同，对一维非稳态导热微分方程中的扩散项采用中心差分；非稳态项采用向前差分。则有：2111xttat iniiniin（414a）变形得： ininiin txatxat 2121（414b）求解非稳态导热微分方程，是从已知的初始温度分布出发，根据边界条件依次求得以后各个时间层上的温度值。由此可见，只要 i 时层上各节点的温度已知，那么 i+1 时层上各节点的温度即可算出，且不需设立方程组求解。此关系式即为显式差分格式。优点：计算工作量小；缺点：受时间及空间步长的限制。（2）隐式差分格式对一维非稳态导热微分方程中的
20、扩散项在(i+1)时层上采用中心差分，非稳态项在节点（n,i+1）处对节点（n,i）采用向前差分，得：(4-15)2111xttat iniiniin 式中已知的是 i 时层上的值 ，而未知量有 3 个，因此不能直接由上式立即算int出 ，而必须求解（i+1）时层上的一个联立方程组，才能算出(i+1)时层各1int节点的温度，此种差分格式称隐式差分格式。优点：不受时间及空间的步长影响；缺点：计算工作量大。综上可知： 非稳态导热微分方程中，扩散项采用中心差分，非稳态项采用向前差分得到显式差分格式； 非稳态导热微分方程中，扩散项采用中心差分，非稳态项采用向后差分得到隐式差分格式。2.热平衡法（1）
21、优点：不受网格是否均匀限制；不受物性是否为常数限制。（2）求解方法以一维非稳态导热边界节点为例，应用热平衡法建立节点离散方程。如图4-9 所示，一无限大平板，右侧面受周围流体的冷却，表面传热系数为 ，此时h边界节点 N 代表了宽为 的元体。2x根据傅立叶定律，在 i 时层上，从节点(N-1)传导给节点 N 的热流量，即从(N-1)传给元体 单位面积的热流量为：xtqiNiiN1,1 （a）根据牛顿冷却公式，平板右侧被冷却时，在 i 时层上其单位面积的热流量为：Nfcthq（b）在 i 时层上元体热力学能的增量为：iNijNtxctcq1, 2（c）据能量守恒定律可知：在 i 时层通过导热和对流
22、进入元体的能量应等于元体热力学能的变化量，则有：（4-iNiiNfiNi txcthxt 11 216a）经整理得：（4- fiNiNi txchtxaxcht 2211216b）其中 是以 为特征长度的傅里叶数，称2xa网格傅里叶数，记为： 。Fo一项可作如下变化：xchBiFoxhcxh 2称为网格毕渥数，记为 。xhBi于是式（4-16b）可改写为：（4- fiNiNi tBiFotFoit 221116c）说明：对多维非稳态导热问题应用热平衡法来建立离散方程的原则与过程与之类似。至此，可以把第三类边界条件下、厚度为 的无限大平板的数值计算问题2作一归纳。由于问题的对称性，只要求解一半厚
23、度即可，其数学描写见式（3-11）到（3-14）。此处不再重复。设将计算区域等分为（N1）等份（N 个节点），节点 1 为绝热的对称面，节点 N 为对流边界，则与微分形式的数学描写相对应的离散形式为：1,2211ntFotFot ininii （417）Ntn,01（418）（4- fiNiNi tBiFotFoBit 2211 19）iit12（420）其中式（420）是绝热边界的一种离散方式，在确定 之值时需要用到 。1it it1根据对称性该值等于 。这样从已知的初始分布 出发，利用式（417）至it2 0（4-19）可以依次求得第二时层、第三时层直到 I 时层上的温度值（见图48）。至
24、于空间步长 及时间步长 的选取，原则上步长越小，计算结果x越接近于精确解，但所需的计算机内存及计算时间则大大增加。此外， 与之间的关系还受到显式格式稳定性的影响。x三、讨论一维导热问题显式差分格式稳定性限制的物理意义从离散方程的结构分析，对于一维导热显式格式的内节点方程，点 n 上i+1 时刻的温度是在该点 i 时刻温度的基础上计及了左右两邻点温度的影响后得出的。若两邻点的影响保持不变，则合理的情况是： 越高，则 越高；int1int越低，则 越低。int1int在上式中，满足这种合理性是有条件的，即上式中 前的系数必大于等于int零，即： （4-21xaFo21）否则，将出现不合理情况。若
25、，则表明节点（n,i）在 i 时刻的01越高，经 时段后， 越低，这种节点温度随时间的跳跃式变化是不符合intint物理规律的，所以称该方程具有不稳定性。对于一维导热显示格式的对流边界节点方程，得出合理解的条件是：（4-021FoBi22a）即： （4-i2122b）由此可见：（1）对流边界节点要得到的合理的解，其限制条件比内节点更为严格，所以，当由边界条件及内节点的稳定性条件得出的 不同时，应选Fo较小的 来确定允许采用的时间步长 。Fo（2）对于第一、二类边界条件，其限制条件只有内节点的限制条件。（3）内边界节点差分方程的稳定性条件不同，但在数值计算时，二节点又必须选择相同的 、 。因此，
26、在选择了 后，则 的选择就要受到稳定条件xx的限制，不能任意选择，而必须按两节点的稳定性条件分别计算 ，取其中较小的 作为时间步长，方能满足二者稳定性要求。四、数值解法的求解步骤1.首先选择空间坐标间隔 ，即距离步长。对二维问题一般使 ；xyx2.对显式格式差分方程，根据方程的稳定性条件选择允许的最大时间步长 ； 在稳定性条件允许范围内， 越大，计算工作量越小，但精度较差；对一维问题，一般取 ，即可满足工程精度要求；对于隐式差分方程，214Fo、 可任意选取，不必进行稳定性条件校核；x3.按题意给定的初始温度分布，确定各节点上的温度初值 ；0nt4.根椐建立的差分方程组，求 时刻各节点的温度
27、；15.再由 为初值，重复计算得出 ，如此反复，最后得到 i 时刻的 。1nt 2nt int例题：4-6 极坐标中常物性、无内热源的非稳态导热方程为：221tttarr试利用本题附图中的符号，列出节点 的,ij差分方程式。解：利用热平衡法，对 节点列出能量平,ij衡式，有 inE 1 1, ,.k kijij ijijjt tcVcr ,1, ,11,1,22k kijij ijijinj jkkijijijijt trrrrttr 整理，得： 1, ,1, ,121,1,22k k kj jijij ijij ijijj jkkijijijijjjrrt ttattrr 与由控制方程得到的结果相同。
