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

1、武器传热学武器传热学主讲：张领科主讲：张领科1 、重点内容：、重点内容： 掌握导热问题数值解法的基本思路；掌握导热问题数值解法的基本思路； 利用热平衡法和泰勒展开法建立节点利用热平衡法和泰勒展开法建立节点的离散方程。的离散方程。2 、掌握内容：、掌握内容：数值解法的实质。数值解法的实质。 3 、重点掌握内容：、重点掌握内容：非稳态导热问题的两非稳态导热问题的两种差分格式及其稳定性。种差分格式及其稳定性。 第四章第四章 导热问题的数值解法导热问题的数值解法传热学研究手段传热学研究手段- -数值解法数值解法 对物理问题进行数值解法的基本思路可对物理问题进行数值解法的基本思路可以概括为：以概括为：把

2、原来在时间、空间坐标系中连把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场等，续的物理量的场，如导热物体的温度场等，用有限个离散点上的值的集合来代替，通过用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。数方程，来获得离散点上被求物理量的值。该方法称为数值解法。该方法称为数值解法。 这些离散点上被求物理量值的集合称为这些离散点上被求物理量值的集合称为该物理量的数值解。该物理量的数值解。 4.1.1 4.1.1 基本思想基本思想建立控制方程及定解条件建立控制方程及定解条件确定节点（

3、区域离散化）确定节点（区域离散化）建立节点物理量的代数方程建立节点物理量的代数方程(组组)设立温度场的初值、边值设立温度场的初值、边值求解代数方程（组）求解代数方程（组）是否收敛是否收敛解的分析解的分析改进初场改进初场/边边界界/步长步长是是否否4.1.2 4.1.2 导热问题数值求解基本步骤导热问题数值求解基本步骤图图4-1 物理问题数值求解过程物理问题数值求解过程 4.1.2 4.1.2 导热问题数值求解基本步骤导热问题数值求解基本步骤图图4-2a 二维矩形域导热问题二维矩形域导热问题 以图以图4-2a4-2a所示二维矩形域内无内热源、稳态、常物性所示二维矩形域内无内热源、稳态、常物性的导

4、热问题对数值解法的六个步骤作的导热问题对数值解法的六个步骤作进一步说明：进一步说明： （1 1）建立控制方程及定解条件）建立控制方程及定解条件 针对图示的导热问题，它的控制方程（即针对图示的导热问题，它的控制方程（即导热微分方程导热微分方程）为：）为： 22220ttxy描述物理问题的微分方程常称控制方程。描述物理问题的微分方程常称控制方程。由图可知，四个边界条件分别为第一类和第三类边界条件。由图可知，四个边界条件分别为第一类和第三类边界条件。4.1.2 4.1.2 导热问题数值求解基本步骤导热问题数值求解基本步骤（2 2）区域离散化（确立节点）区域离散化（确立节点） 4.1.2 4.1.2

5、导热问题数值求解基本步骤导热问题数值求解基本步骤二维矩二维矩形域内形域内稳态无稳态无内热源，内热源，常物性常物性的导热的导热问题问题图图4-2b 导热问题数值求解网格划分导热问题数值求解网格划分4.1.2 4.1.2 导热问题数值求解基本步骤导热问题数值求解基本步骤 （3 3）建立节点物理量的代数方程（离散方程）建立节点物理量的代数方程（离散方程） 节点上物理量的代数方程称离散方程节点上物理量的代数方程称离散方程(discretization(discretization equation) equation)。其过程如下：。其过程如下： l 首先划分各节点的类型；首先划分各节点的类型； l

6、其次，建立节点离散方程；其次，建立节点离散方程； l 最后，代数方程组的形成。最后，代数方程组的形成。 对节点对节点 (m,n(m,n) ) 的代数方程，当的代数方程，当 x=x=y y 时，有：时，有： ,1,1,1,11()4m nmnmnm nm nttttt4.1.2 4.1.2 导热问题数值求解基本步骤导热问题数值求解基本步骤（4 4） 设立迭代初场设立迭代初场 代数方程组的求解方法有代数方程组的求解方法有直接解法直接解法与与迭代迭代解法解法，传热问题的，传热问题的有限差分法有限差分法中主要采用迭代中主要采用迭代法。采用迭代法求解时，法。采用迭代法求解时，需对被求的温度场预需对被求的

7、温度场预先设定一个解，这个解称为初场先设定一个解，这个解称为初场(initial (initial field)field)，并在求解过程中不断改进。，并在求解过程中不断改进。 4.1.2 4.1.2 导热问题数值求解基本步骤导热问题数值求解基本步骤（5 5） 求解代数方程组求解代数方程组 求解时遇到的问题：求解时遇到的问题： 线性；线性； 非线性；非线性； 收敛性收敛性等。等。 如图如图4-2b4-2b，除，除 m=1 m=1 的左边界上各节点的温的左边界上各节点的温度已知外，其余度已知外，其余 (M-1)N (M-1)N 个节点均需建立离个节点均需建立离散方程，共有散方程，共有 (M-1)

8、N (M-1)N 个方程，则构成一个个方程，则构成一个封闭的代数方程组。封闭的代数方程组。 1 1 ）线性代数方程组：）线性代数方程组：代数方程一经建立，代数方程一经建立，其中各项系数在整个求解过程中不再变化；其中各项系数在整个求解过程中不再变化； 4.1.2 4.1.2 导热问题数值求解基本步骤导热问题数值求解基本步骤2 2 ）非线性代数方程组：）非线性代数方程组：代数方程一经建立，代数方程一经建立，其中各项系数在整个求解过程中不断更新。其中各项系数在整个求解过程中不断更新。 ( )(1)( )(1)( )( )(1)( )maxmaxmaxmaxkkiikkiikikkiikttttttt

9、t3 3 ）是否收敛判断：）是否收敛判断：是指用是指用迭代法求解代数方程是否收敛，迭代法求解代数方程是否收敛，即本次迭代计算所得之解与上即本次迭代计算所得之解与上一次迭代计算所得之解的偏差一次迭代计算所得之解的偏差是否小于允许值。是否小于允许值。 4.1.2 4.1.2 导热问题数值求解基本步骤导热问题数值求解基本步骤（6） 解的分析解的分析 通过求解代数方程，获得物体中的温度通过求解代数方程，获得物体中的温度分布，根据温度场应进一步计算通过的热流分布，根据温度场应进一步计算通过的热流量，热应力及热变形等。因此，对于数值分量，热应力及热变形等。因此，对于数值分析计算所得的温度场及其它物理量应作

10、详细析计算所得的温度场及其它物理量应作详细分析，以获得定性或定量上的结论。分析，以获得定性或定量上的结论。 4.1.2 4.1.2 导热问题数值求解基本步骤导热问题数值求解基本步骤 以上六个步骤中，控制方程和定解条件前面已以上六个步骤中，控制方程和定解条件前面已经详细叙述；对于规则的区域，网格的划分与节点经详细叙述；对于规则的区域，网格的划分与节点的生成比较容易；对数值解的分析是解决实际问题的生成比较容易；对数值解的分析是解决实际问题过程中的重要一步，但不涉及数值求解技术；过程中的重要一步，但不涉及数值求解技术； 本章的重点放在如何建立离散方程组，以及本章的重点放在如何建立离散方程组，以及如何

11、求解离散方程组。对于非稳态导热，除了空间如何求解离散方程组。对于非稳态导热，除了空间区域的离散，还要离散时间坐标。区域的离散，还要离散时间坐标。4.1.2 4.1.2 导热问题数值求解基本步骤导热问题数值求解基本步骤(1) Taylor(1) Taylor（泰勒）级数展开法（泰勒）级数展开法；(2) (2) 多项式拟合法；多项式拟合法；(3) (3) 控制容积积分法；控制容积积分法；(4) (4) 控制容积平衡法控制容积平衡法( (也称为热平衡法也称为热平衡法) ) 4.2 4.2 内节点离散方程的建立方法内节点离散方程的建立方法接下来我们介绍稳态导热问题中位于计算区域内接下来我们介绍稳态导热

12、问题中位于计算区域内部的节点的离散方程的建立方法。部的节点的离散方程的建立方法。二维矩二维矩形域内形域内稳态无稳态无内热源，内热源，常物性常物性的导热的导热问题问题图图4-2b 导热问题数值求解网格划分导热问题数值求解网格划分4.1.2 4.1.2 导热问题数值求解基本步骤导热问题数值求解基本步骤根据泰勒级数展开式，用节点根据泰勒级数展开式，用节点(m,n(m,n) )的温度的温度t tm,nm,n来表示节点来表示节点(m+1,n)(m+1,n)而温度而温度t tm+1,nm+1,n用节点用节点(m,n(m,n) )的温度的温度t tm,nm,n来表示节点来表示节点(m-1,n)(m-1,n)

13、的的温度温度t tm-1,nm-1,n2233441,234,2624mnmnmnmnmntxtxtxtttxxxxx223341,234,42624mnmnmnmnmntxtxtxtttxxxxx4.2.1 4.2.1 泰勒基数展开法泰勒基数展开法将上两式相加可得将上两式相加可得24421,1,24,212mnmnm nm ntxttttxxx22,mntx将上式改写成将上式改写成 的表达式，有的表达式，有)(222, 1, 1,22xoxtttxtnmnmnmnm)(2221,1,22yoytttytnmnmnmnm同样可得：同样可得：表示未明确写出的表示未明确写出的级数余项中的级数余项中

14、的X X的最低阶数为的最低阶数为2 24.2.1 4.2.1 泰勒基数展开法泰勒基数展开法 根据导热问题的控制方程根据导热问题的控制方程 ( ( 导热微分方程导热微分方程 ) )1,1,1,122220mnmnmnmnmnmnttttttxy若若 x=x=y y 则有则有 ,1,1,1,11()4m nmnmnm nm nttttt22220ttxy得得4.2.1 4.2.1 泰勒基数展开法泰勒基数展开法 在传热学问题的控制方程中，主要遇到的是一阶与二阶导在传热学问题的控制方程中，主要遇到的是一阶与二阶导数，在均分网格中，一阶二阶常见的离散表达式数，在均分网格中，一阶二阶常见的离散表达式( (

15、差分表示式差分表示式) )如下表所示。如下表所示。4.2.1 4.2.1 泰勒基数展开法泰勒基数展开法基本思想：基本思想：是傅里叶导热定律和能量守恒定律的体是傅里叶导热定律和能量守恒定律的体现。对每个元体，可用傅里叶导热定律写出其能量现。对每个元体，可用傅里叶导热定律写出其能量守恒的表达式。守恒的表达式。4.2.2 4.2.2 热平衡法热平衡法如图所示，如图所示， 从节点从节点 (m-1,n) (m-1,n) 通过界面通过界面 w w 传导到传导到节点节点 (m,n(m,n) ) 的热流量：的热流量： 1,mnm nwttyx 同理：通过界面同理：通过界面 e,n,se,n,s 传导给节点（传

16、导给节点（ m,nm,n ）的热流量）的热流量4.2.2 4.2.2 热平衡法热平衡法1,mnm nettyx ,1,m nm nnttyx ,1,m nm nsttyx 对元体对元体 (m,n(m,n). ). 根据能量守恒定律可知：根据能量守恒定律可知： 0ewns 其中，其中，规定：规定：导入元体（导入元体（ m,nm,n ）的热流量）的热流量为正；导出元体（为正；导出元体（ m,nm,n ）的热流量为负。）的热流量为负。 4.2.2 4.2.2 热平衡法热平衡法若若 x=x=y y 则有则有 ,1,1,1,11()4m nmnmnm nm nttttt说明：说明： 上述分析与推导是在笛

17、卡儿坐标系中进行的；上述分析与推导是在笛卡儿坐标系中进行的； 热平衡法概念清晰，过程简捷；热平衡法概念清晰，过程简捷； 热平衡法与建立微分方程的思路与过程一致，热平衡法与建立微分方程的思路与过程一致，但不同的是前者是有限大小的元体，后者是微元但不同的是前者是有限大小的元体，后者是微元体。体。 4.2.2 4.2.2 热平衡法热平衡法 对于第一类边界条件的热传导问题，处理比较简单，因对于第一类边界条件的热传导问题，处理比较简单，因为已知边界的温度，可将其以数值的形式加入到内节点的离为已知边界的温度，可将其以数值的形式加入到内节点的离散方程中，组成封闭的代数方程组，直接求解。散方程中，组成封闭的代

18、数方程组，直接求解。 而对于第二类或第三类边界条件的导热问题，所有内节而对于第二类或第三类边界条件的导热问题，所有内节点的离散方程组成的代数方程组是不封闭的，因未知边界温点的离散方程组成的代数方程组是不封闭的，因未知边界温度，因而应对位于该边界上的节点补充相应的代数方程，才度，因而应对位于该边界上的节点补充相应的代数方程，才能使方程组封闭，以便求解。能使方程组封闭，以便求解。 为了求解方便，这里我们将第二类边界条件及第三类边为了求解方便，这里我们将第二类边界条件及第三类边界条件合并起来考虑，用界条件合并起来考虑，用q qw w表示边界上的热流密度或热流密度表示边界上的热流密度或热流密度表达式。

19、表达式。为使结果更具一般性，假设物体具有内热源为使结果更具一般性，假设物体具有内热源 ( ( 不不必均匀分布必均匀分布 ) ) 。4.3.1 4.3.1 边界节点离散方程的建立边界节点离散方程的建立如图所示如图所示 边界节点边界节点 (m,n(m,n) ) 只能代表半个元只能代表半个元体，若边界上有向该元体，若边界上有向该元体传递的热流密度，据体传递的热流密度，据能量守恒定律对该元体能量守恒定律对该元体有：有： (1) (1) 平直边界上的节点平直边界上的节点1,1,1,0222mnm nm nm nm nm nm nwttttttxxxyyyqxyy yx2,1,1,12124m nwm n

20、mnm nm nx xqtttt4.3.1 4.3.1 边界节点离散方程的建立边界节点离散方程的建立(2) (2) 外部角点外部角点2,1,12122m nwm nmnm nx xqttt1,1,02242mnm nm nm nm nwttttx yyxxyqxy yx如图所示，二维墙角计算区域中，该节点外角点仅如图所示，二维墙角计算区域中，该节点外角点仅代表代表 1/4 1/4 个以个以 为边长的元体。假设边为边长的元体。假设边界上有向该元体传递的热流密度为界上有向该元体传递的热流密度为 ，则据能量，则据能量守恒定律得其热平衡式为：守恒定律得其热平衡式为： xy、wq(3) (3) 内部角点

21、内部角点22,1,1,11,213(22)62wm nmnm nm nmnx qxttttt1,1,1,1,230242mnm nm nm nm nm nmnm nm nwttttttxyxxyyttyx yxyqx yx如图所示内部角点代表了如图所示内部角点代表了 3/4 3/4 个元体，在同样的个元体，在同样的假设条件下有假设条件下有讨论关于边界热流密度的三种情况：讨论关于边界热流密度的三种情况： （1 1）绝热边界）绝热边界即令上式即令上式 即可。即可。 0wq （2 2） 值不为零值不为零wq流入元体，流入元体， 取正，流出元体，取正，流出元体， 取负使取负使用上述公式用上述公式 wq

22、wq（3 3）对流边界）对流边界此时此时 ，将此表达式代入上述方程，将此表达式代入上述方程，并将此项中的并将此项中的 与等号前的与等号前的 合并。对合并。对于于 的情形有的情形有)(,nmfwtthq,m nt,m ntxy （a a）平直边界）平直边界（b b）外部角点）外部角点（c c）内部角点）内部角点2,1,1,12222m nm nmnm nm nfh xxh xttttt2,1,12212m nm nmnm nfh xxh xtttt2,1,11,1322322m nm nmnm nmnm nfh xxh xtttttt2 代数方程的求解方法代数方程的求解方法 2 2） 迭代法：迭

23、代法：先对要计算的场作出假设（设先对要计算的场作出假设（设定初场），在迭代计算中不断予以改进，直定初场），在迭代计算中不断予以改进，直到计算前的假定值与计算结果相差小于允许到计算前的假定值与计算结果相差小于允许值为止的方法，称迭代计算收敛。值为止的方法，称迭代计算收敛。1 1） 直接解法：直接解法：通过有限次运算获得精确通过有限次运算获得精确解的方法，如：矩阵求解，高斯消元法。解的方法，如：矩阵求解，高斯消元法。 2 迭代法目前应用较多的是：迭代法目前应用较多的是： 1 1 ）高斯）高斯赛德尔迭代法：赛德尔迭代法：每次迭代计算，每次迭代计算，均是使用节点温度的最新值。均是使用节点温度的最新值。

24、 2 2 ）用雅可比迭代法：）用雅可比迭代法：每次迭代计算，均用每次迭代计算，均用上一次迭代计算出的值。上一次迭代计算出的值。 设有一三元方程组设有一三元方程组： 11 112 213 3121 122 223 3231 132 233 33a ta ta tba ta ta tba ta ta tb其中其中 （ i=1,2,3 i=1,2,3 ； j=1,2,3 j=1,2,3 ）及）及 是已知的系数（均不为零）及常数。是已知的系数（均不为零）及常数。, i jaib采用高斯采用高斯赛德尔迭代法的步骤：赛德尔迭代法的步骤： （1 1）将三元方程变形为迭式方程：）将三元方程变形为迭式方程： 1

25、112 213 3112221 123 3223331 132 2331()1()1()tba ta tatba ta tatba ta ta（2 2）假设一组解（迭代初场），记为）假设一组解（迭代初场），记为： 并代入迭代方程求得第一并代入迭代方程求得第一 次解次解 每次计算均用最新值每次计算均用最新值代入。代入。 (0)(0)(0)123ttt、(1)(1)(1)123ttt、 、（3 3）以新的初场）以新的初场 重复计算，直到相邻两重复计算，直到相邻两次迭代值之差小于允许值，则称迭代收敛，次迭代值之差小于允许值，则称迭代收敛，计算终止。计算终止。 判断迭代是否收敛的准则：判断迭代是否收敛

26、的准则：)(max)() 1()()() 1()() 1(maxmaxmaxkkikikikikikikittttttttk k及及k+1k+1表示迭代次数；表示迭代次数；第第k k次迭代得到的最大值次迭代得到的最大值(k)maxt当有接近于零的当有接近于零的t t 时，第三个较好时，第三个较好36 1010 允许的偏差；相对偏差 值一般取说明：说明： 1 1 ）对于一个代数方程组，若选用的迭代方）对于一个代数方程组，若选用的迭代方式不合适，有可能导致发散，即称式不合适，有可能导致发散，即称迭代过程迭代过程发散发散； 2 2 ）对于常物性导热问题，组成的差分方程）对于常物性导热问题，组成的差分

27、方程组，迭代公式的选择应使一个迭代变量的系组，迭代公式的选择应使一个迭代变量的系数总是大于或等于该式中其他变量系数绝对数总是大于或等于该式中其他变量系数绝对值的代数和，此时，结果一定收敛。值的代数和，此时，结果一定收敛。 3 3 ）采用热平衡法导出差分方程时，若每）采用热平衡法导出差分方程时，若每一个方程都选用导出该方程中心节点的温一个方程都选用导出该方程中心节点的温度作为迭代变量，则上述条件必满足，迭度作为迭代变量，则上述条件必满足，迭代一定收敛。代一定收敛。 121321233132112233111aaaaaaaaa，这一这一条件条件数学上称主对角线占优（对角占数学上称主对角线占优（对角

28、占优）；优）； 与稳态导热相比，非稳态导热多了一个非稳态项。与稳态导热相比，非稳态导热多了一个非稳态项。4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法22220ttxy二维稳态导热：二维稳态导热：一维非稳态导热：一维非稳态导热：22ttax4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法管装药管装药圆柱七孔药圆柱七孔药底排药柱底排药柱固体发射药常用的点火模型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法1 1、固定边界、火药不透明且表面接触温度保持不变、固定边界、火药不透明且表面接触温度保持不变 2222()(0)()()(

29、0)cccccccciibgTTTTcaatytycTTtTTTyTTTy 或或固体发射药常用的点火模型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法1 1、固定边界、火药不透明且表面接触温度保持不变、固定边界、火药不透明且表面接触温度保持不变 固体发射药常用的点火模型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法2 2、固定边界、不透明、表面接受辐射热通量、固定边界、不透明、表面接受辐射热通量2222()(0)()(0)cccccccciicTTTTcaatytycTTtTTTxTqxy 或固体发射药常用的点火模

30、型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法2 2、固定边界、不透明、表面接受辐射热通量、固定边界、不透明、表面接受辐射热通量固体发射药常用的点火模型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法3 3、热气体对流加热简化模型、热气体对流加热简化模型 22220()(0)()()()cccccccciicgsyTTTTcaatytycTTtTTTyTh TThconsty 或固体发射药常用的点火模型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法3 3、热气体对流加热简化模

31、型、热气体对流加热简化模型 固体发射药常用的点火模型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法1 1、壁面固定温度加热模型、壁面固定温度加热模型22exp(/)(0)()(0)ccccciisiTTcBERTtyTTtTTTyTTTy 固体发射药常用的点火模型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法1 1、壁面固定温度加热模型、壁面固定温度加热模型固体发射药常用的点火模型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法2 2、壁面有辐射热通量的点火、壁面有辐射热通量的

32、点火220exp(/)(0)()ccccciireyTTcBERTtyTTtTTTyTqy 固体发射药常用的点火模型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法2 2、壁面有辐射热通量的点火、壁面有辐射热通量的点火0.000000.000050.000100.000150.00020300350400450500550600650700T/Kx/m t1.6ms t3.2ms t6.3ms 8.4ms t8.6ms固体发射药常用的点火模型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法3 3、两相流点火模型、两相流

33、点火模型12122112exp()()( , ) ()()|cccccccccccfffffffrrffttirTTTTnEccrAQryrtyyyyRTTTTncrdyrtyyyTy tTTh TTy初始条件11,11()0 ( , )( , )crd tfcfr tr tcfTTTTyyT r tTr t（边界条件）（连接条件）固体发射药常用的点火模型固体发射药常用的点火模型4.44.4非稳态导热问题的非稳态导热问题的 数值解法数值解法3 3、两相流点火模型、两相流点火模型0.000.020.040.060.080.100.120.140.1630035040045050055060065

34、0700T/Ky/mm t=0.30ms t=0.94ms t=4.07ms t=6.81ms t=8.72ms0123456789300350400450500550600650700T/Kt/ms y=0.0um y=2.0um y=4.0um y=6.0um4.4.1 4.4.1 时间时间- -空间区域的离散化空间区域的离散化一维非稳态导热时间一维非稳态导热时间空间区域离散空间区域离散4.4.1 4.4.1 时间时间- -空间区域的离散化空间区域的离散化将函数将函数t在节点（在节点（n,i+1）对点）对点（n,i）作泰勒展开：作泰勒展开：223344(1)( )234,2624iinnn in in in itttttt因此：因此：(1)( ),()iinnn itttO(1)( ),iinnn ittt向前差