有限元法基础11热传导与热应力

1、 11 11 传热分析与热应力传热分析与热应力 本章参考书本章参考书 1.孔祥谦孔祥谦. 有限单元法在传热学中的应用有限单元法在传热学中的应用. 第三版，科学出版社，第三版，科学出版社， 1998 2.王勖成王勖成. 有限单元法有限单元法. 清华大学出版社，清华大学出版社，2003. 第第12章章 3.曾攀曾攀. 有限元分析及应用有限元分析及应用. 清华大学出版社，清华大学出版社，2004 . 第第8章章 1 11. 11. 传热分析与热应力传热分析与热应力 l传热是广泛存在的自然现象，只要有温度差存在，就会有热量的传传热是广泛存在的自然现象，只要有温度差存在，就会有热量的传 递，只要有热量的

2、输入和输出，就会引起温度的变化。递，只要有热量的输入和输出，就会引起温度的变化。 l传热分析的对象是固体、液体和气体，其应用包括热量交换、化学传热分析的对象是固体、液体和气体，其应用包括热量交换、化学 反应、材料相变、能量转换等。反应、材料相变、能量转换等。 l温度的变化和不均匀分布，引起结构出现应力变化，称为热应力。温度的变化和不均匀分布，引起结构出现应力变化，称为热应力。 l当以应力分析为目的时，为确定温度场，需要对固体进行传热计算，当以应力分析为目的时，为确定温度场，需要对固体进行传热计算， 以便确定相关的热应力。以便确定相关的热应力。 2 11 11 传热分析与热应力传热分析与热应力

3、11.1 传热问题的基本方程传热问题的基本方程 l固体热传导的现象固体热传导的现象 3 11 11 传热分析与热应力传热分析与热应力 l术语和单位术语和单位 在国际标准单位制中，传热分析的术语和单位在国际标准单位制中，传热分析的术语和单位 c比热容比热容J/(kgK)T温度（温度（K 或或 C） Q热流（热流（W/m2） h对流换热系数对流换热系数 W/(m2K) t时间（时间（s） k导热系数导热系数W/(mK)质量密度（质量密度（kg/m3） q单位体积热生成率单位体积热生成率 (W/m3) Stefan-Boltzman常数常数 =5.6710-8 W/(m2K4) T /(K /s)T

4、t 4 11 11 传热分析与热应力传热分析与热应力 l控制方程控制方程 对于微元对于微元dxdydz，生成的热量为，生成的热量为 微元体内的净流出热流量为微元体内的净流出热流量为 由于热量的储存使内能增加，即由于热量的储存使内能增加，即 由能量守恒定律，在微元内有由能量守恒定律，在微元内有 d d dqx y z ()d d d y xz Q QQ x y z xyz d d dcx y z T d d dd d d()d d d y xz Q QQ cx y z Tq x y zx y z xyz () y xz Q QQ cTq xyz 5 11 11 传热分析与热应力传热分析与热应力

5、l由于方程不封闭，需补充条件。各向同性均匀材料的由于方程不封闭，需补充条件。各向同性均匀材料的Fourier定律定律 于是，有于是，有 l第一类边界条件：边界上已知温度第一类边界条件：边界上已知温度T，即即 l第二类边界条件：边界上已知热流密度，即第二类边界条件：边界上已知热流密度，即 l第三类边界条件：已知与物体相接触的流体介质的温度和换热系数第三类边界条件：已知与物体相接触的流体介质的温度和换热系数 为已知，即为已知，即 gradkT Q 2 c TkTq ( , , , ) w TTTf x y z t 或 ( , , , ) nwn TT QkQQkg x y z t nn 或 ()

6、nf T Qkh TT n 6 11 11 传热分析与热应力传热分析与热应力 11.2 变分原理与有限元变分原理与有限元 l 瞬态热传导问题变分泛函为瞬态热传导问题变分泛函为 l 稳态问题，温度不随时间变化稳态问题，温度不随时间变化 l泛函的变分取驻值，可得控制方程和第二类和第三类边界条件泛函的变分取驻值，可得控制方程和第二类和第三类边界条件 第一类边界条件应强制满足，称为本质边界条件；第一类边界条件应强制满足，称为本质边界条件； 第二、第三类边界条件是自然边界条件。第二、第三类边界条件是自然边界条件。 23 2 11 ()dd() d 222 wf SS k TqTc TTQ T Ah TT

7、 T A 0 T t 23 2 11 ()dd() d 222 wf SS k TqTQ T Ah TT T A 7 11 11 传热分析与热应力传热分析与热应力 l有限元法有限元法 将物体离散为将物体离散为n个单元体，即个单元体，即 ，将单元内的温度场用节，将单元内的温度场用节 点上的温度插值，有点上的温度插值，有 代入泛函，泛函称为代入泛函，泛函称为 e n e T NT 3 23 1 ( ) ( ) 2 ( ) ( )d , d d ,R= d d ee eeee eTeeTe n TT S T TTT wf SS TkhTTc TR kkNNhh NNA ccNNq NdQNAhTNA

8、 8 11 11 传热分析与热应力传热分析与热应力 l将单元矩阵组合后，得到物体的总体矩阵，对泛函变分取极值，即将单元矩阵组合后，得到物体的总体矩阵，对泛函变分取极值，即 得有限元方程组得有限元方程组 KT称为热传导矩阵，称为热传导矩阵，C为热容矩阵，为热容矩阵，RT为等效节点温度载荷列阵。为等效节点温度载荷列阵。 l对稳态问题，上述方程是线性的，热传导方程是对称正定的，求解对稳态问题，上述方程是线性的，热传导方程是对称正定的，求解 类方法似于结构分析。类方法似于结构分析。 l对瞬态问题，需采用有限差分法，将对瞬态问题，需采用有限差分法，将 离散，或采用显式时间积离散，或采用显式时间积 分，如

9、中心差分，或采用隐式时间积分，如分，如中心差分，或采用隐式时间积分，如Newmark差分。差分。 /0T TT C TKTR T 9 11 11 传热分析与热应力传热分析与热应力 l 对线性瞬态问题，还可采用结构分析中的模态迭加法对线性瞬态问题，还可采用结构分析中的模态迭加法 首先求特征问题首先求特征问题 每个特征向量每个特征向量T i相对于相对于C正则化，即正则化，即 令令 是模态矩阵，它的每一列是正则化的特征向量是模态矩阵，它的每一列是正则化的特征向量T i，于是，于是 可将节点温度表示为广义温度可将节点温度表示为广义温度Z的关系的关系 将其代入有限元方程，并左乘将其代入有限元方程，并左乘

10、 得到得到n个解耦的方程组个解耦的方程组 积分上述方程组后，得积分上述方程组后，得Z(t)，由此可得到节点，由此可得到节点T(t)。 ( ) 0 T KCT ) 1 T ii TCT , TT T CIK TZ T , T iiiiiiT ZZPPR 10 11 11 传热分析与热应力传热分析与热应力 11.3热辐射热辐射 考虑两个无限大的平行平面，由于无限大，不用考虑边界效应。设考虑两个无限大的平行平面，由于无限大，不用考虑边界效应。设 每个平面都有均匀温度，平面每个平面都有均匀温度，平面1的温度为的温度为T1，平面，平面2的温度为的温度为T2，平面，平面 都是理想的黑体，因此每个平面都是理

11、想的吸收体和辐射体，平面表都是理想的黑体，因此每个平面都是理想的吸收体和辐射体，平面表 面的热流量为面的热流量为 是是Stefan-Boltzman常数。由于实际的辐射面并非理想黑体，也不常数。由于实际的辐射面并非理想黑体，也不 是无穷大的平面，也不一定平行，因此把面积为是无穷大的平面，也不一定平行，因此把面积为A1和温度和温度T1表面所接表面所接 收的热流量表示为收的热流量表示为 包含了各种因素引起的辐射折减，包括视图因子和辐射率等。包含了各种因素引起的辐射折减，包括视图因子和辐射率等。 2 4444 11212 ()()QTTQTT 2 44 11 ()QTT 11 11 11 传热分析与

12、热应力传热分析与热应力 l由于辐射面是有限的、非平行的，用视图因子表示由于辐射面是有限的、非平行的，用视图因子表示 l对于两个无限大的平行面为对于两个无限大的平行面为1，对于两个相互看不见的平面是，对于两个相互看不见的平面是0 12 11 11 传热分析与热应力传热分析与热应力 l与面积为与面积为A1交换辐射能的表面有多少个，就有多少个式子。如果交换辐射能的表面有多少个，就有多少个式子。如果A1 不是很大，可认为不是很大，可认为Q1在在A1上是个常数，因此上是个常数，因此 l与对流边界比较，该式与之相同。因此在有限元方程中将与对流边界比较，该式与之相同。因此在有限元方程中将h项中换项中换 成成

13、 就可用有限元分析辐射问题。就可用有限元分析辐射问题。 l需注意，需注意， 是与温度有关的，故辐射问题是高度非线性问题。是与温度有关的，故辐射问题是高度非线性问题。 l当材料常数是温度的函数，问题也是非线性的。当材料常数是温度的函数，问题也是非线性的。 l非线性方程可采用非线性方程可采用Newton-Raphson法求解，但有一些特殊的适合法求解，但有一些特殊的适合 传热问题的处理方法。传热问题的处理方法。 l为了避免辐射的强非线性，实际问题的处理，有时将为了避免辐射的强非线性，实际问题的处理，有时将 也处理为也处理为 常数，如在热锻时的传热分析。常数，如在热锻时的传热分析。 22 12112

14、21 (),()() radrad QhTThTTTT rad h rad h rad h 13 11 11 传热分析与热应力传热分析与热应力 11.4伴有相变的导热问题伴有相变的导热问题 l特点：控制方程是非稳态导热方程，特点：控制方程是非稳态导热方程， 区域内存在一个随时间移动的两相界区域内存在一个随时间移动的两相界 面，在界面上放出或吸收潜热。面，在界面上放出或吸收潜热。 1891年年J. Stefan关于地极冰层厚度关于地极冰层厚度 的研究首次讨论这一课题。的研究首次讨论这一课题。 l当越过相变区间时，热流密度不连当越过相变区间时，热流密度不连 续，在数学上是一个强非线性问题，续，在数

15、学上是一个强非线性问题， 计算发生困难。计算发生困难。 14 11 11 传热分析与热应力传热分析与热应力 l相变界面的两边各自满足非稳态导热控制方程，一般为了简单略相变界面的两边各自满足非稳态导热控制方程，一般为了简单略 去液相区的自然对流或强制对流等作用。去液相区的自然对流或强制对流等作用。 l在相变界面在相变界面S(t)上，满足温度连续条件上，满足温度连续条件 能量守恒条件能量守恒条件 设设L J/kg为物质的相变潜热，则为物质的相变潜热，则 l由于相变界面的移动，给数值方法带来困难。由于相变界面的移动，给数值方法带来困难。 , SL LSLLSS TT QQQkQk xx 相变潜热 (

16、 ( ), )( ( ), ) SLm T S t tT S t tT 2 d ( ) W/m d S t L t 相变潜热项 15 11 11 传热分析与热应力传热分析与热应力 l焓法模型焓法模型 采用焓（采用焓（HcT）和温度同时作为）和温度同时作为 待求函数。由于相变界面上温度随待求函数。由于相变界面上温度随 时间的变化曲线是间断的，但焓随时间的变化曲线是间断的，但焓随 时间的变化曲线是连续的，因此用时间的变化曲线是连续的，因此用 数值方法求解焓分布时不需跟踪两数值方法求解焓分布时不需跟踪两 相界面，从而使液相区和固相区统相界面，从而使液相区和固相区统 一处理称为可能，焓场解出后，温一处

17、理称为可能，焓场解出后，温 度场可容易得到。度场可容易得到。 16 11 11 传热分析与热应力传热分析与热应力 l焓表示的传热方程为焓表示的传热方程为 其中温度与焓的关系为其中温度与焓的关系为 () h k T t / ()/ SSm mSmSm LSm h chc T TTc Thc TL hLchc TL 17 11 11 传热分析与热应力传热分析与热应力 l有限元法中的处理有限元法中的处理 （1）相变潜热作为附加比热）相变潜热作为附加比热 用加权余量法推导有限元方程，有用加权余量法推导有限元方程，有 使用使用T 的插值的形函数作为权函数，可得的插值的形函数作为权函数，可得Galerki

18、n法的有限元方程。法的有限元方程。 式中式中fs是固相率是固相率,为无因次量为无因次量, , 液相液相0, 固相为固相为1。 2 ()0 S l fT W kTqcLd Tt ()()0 lllS lll WWWfTTTTT kWdSkqWcLWd nxxyyzztt 01 S f 18 11 11 传热分析与热应力传热分析与热应力 （2）相变潜热作为源项）相变潜热作为源项 在前面的公式中，相变潜热项在前面的公式中，相变潜热项 与内热源项与内热源项 作相同作相同 处理，得有限元方程处理，得有限元方程 是与相变项有关的等效载荷。是与相变项有关的等效载荷。 l由于相变得复杂性，在实际数值模拟中，尤

19、其在大型结构计算时，由于相变得复杂性，在实际数值模拟中，尤其在大型结构计算时， 不考虑相变的两相问题，相变的影响只考虑潜热，把它作为一个常不考虑相变的两相问题，相变的影响只考虑潜热，把它作为一个常 热源处理。热源处理。 S l f LW t l qW 0 eeeeee KTNTp p e p 19 11 11 传热分析与热应力传热分析与热应力 11.5 热应力（热应力（Thermal stress） 大多数情况下，传热问题所确定的温度场将直接影响物体的热应力大多数情况下，传热问题所确定的温度场将直接影响物体的热应力, 而热应力对温度场的耦合影响不大，因而可将物体的热问题看成是单而热应力对温度场

20、的耦合影响不大，因而可将物体的热问题看成是单 向耦合过程。向耦合过程。 （一）热应力问题中的物理方程（一）热应力问题中的物理方程 热膨胀系数（热膨胀系数（Thermal expansion coefficent）, 为为Lam系数系数 1 () ijijklijij ijklijkliljkikjlijkl ST SC , 20 11 11 传热分析与热应力传热分析与热应力 l将热应变看作初始应变，上式写为将热应变看作初始应变，上式写为 其中初应变为其中初应变为 ，矩阵形式，矩阵形式 （二）虚功原理（二）虚功原理 弹性力学问题的虚功原理为内力虚功等于外力虚功，即弹性力学问题的虚功原理为内力虚功

21、等于外力虚功，即 将物理方程代入，得将物理方程代入，得 0 ijij T 0 () ijijklijij C 0 ,0,0,0TTTT d(dd )0 t ijijiiii S b ut u A 0 d(ddCd )0 t ijklklijiiiiijklklij S Cb ut u A 21 11 11 传热分析与热应力传热分析与热应力 （三）有限元列式（三）有限元列式 设单元节点列阵为设单元节点列阵为 假设单元内位移由节点位移表示的插值函数为假设单元内位移由节点位移表示的插值函数为 应变可表示为应变可表示为 虚位移与虚应变为虚位移与虚应变为 代入虚功原理，得代入虚功原理，得 其中其中 组装

22、到总体矩阵后，由于组装到总体矩阵后，由于 的任意性，得到有限元方程的任意性，得到有限元方程 111 , eT nnn u v wu v wq e uNq ee DNqBq , ee uNqBq ()()0 e Tee Te m qKqqQ 0 d ,dd eeee t eTeTTT S d KB CBQN bN tB C Kq = Q q 22 11 11 传热分析与热应力传热分析与热应力 （三）求解热应力的方法（三）求解热应力的方法 在有限元分析程序中解热应力问题有两种方法，即直接法和间接法。在有限元分析程序中解热应力问题有两种方法，即直接法和间接法。 l直接法直接法 直接将传热分析和热应力

23、耦合起来分析的方法。在求解时，直接将传直接将传热分析和热应力耦合起来分析的方法。在求解时，直接将传 热边界条件、力学边界条件施加在有限元模型上，以节点温度和位移作热边界条件、力学边界条件施加在有限元模型上，以节点温度和位移作 为未知变量求解。为未知变量求解。 在有限元商业软件中，有多场耦合单元，如在有限元商业软件中，有多场耦合单元，如ANSYS中中Solid5和和 Solid98有热、电、磁、压电和结构的耦合场单元，每节点有有热、电、磁、压电和结构的耦合场单元，每节点有6DOF， 即温度、电势、磁场强度势函数，即温度、电势、磁场强度势函数，3个方向位移。个方向位移。 23 11 11 传热分析

24、与热应力传热分析与热应力 l间接法间接法 热应力问题是一个单向耦合问题，在多数情况下，温度变化会产生热应力问题是一个单向耦合问题，在多数情况下，温度变化会产生 热应变，但热应变部不引起温度变化。这样，可以将热应力问题分成热应变，但热应变部不引起温度变化。这样，可以将热应力问题分成 两个过程来计算分析，即传热分析和热应力计算。两个过程来计算分析，即传热分析和热应力计算。 首先通过传热分析获得结构的温度场分布，然后在已知温度分布的首先通过传热分析获得结构的温度场分布，然后在已知温度分布的 情况下求热应力。情况下求热应力。 优点优点：1）求解规模比直接法小，有更高的计算效率；求解规模比直接法小，有更

25、高的计算效率； 2）在瞬态问题时，可先求解温度场，然后在关心的时间点上）在瞬态问题时，可先求解温度场，然后在关心的时间点上 求热应力，可节省大量的存储空间，即使是非线性问题也求热应力，可节省大量的存储空间，即使是非线性问题也 可这样处理。可这样处理。 24 11 11 传热分析与热应力传热分析与热应力 l例例 稳态温度分析稳态温度分析 一圆柱型容器上，垂直接有一小管道，一圆柱型容器上，垂直接有一小管道， 容器内装有容器内装有450F的液体，小管内有的液体，小管内有100F 的流体。假设容器足够长，远端的温度的流体。假设容器足够长，远端的温度 为为450F。 25 11 11 传热分析与热应力传

26、热分析与热应力 l例例 稳态温度长分析（续）稳态温度长分析（续） 网格图网格图 温度分布温度分布26 11 11 传热分析与热应力传热分析与热应力 l例例 稳态温度常分析（续）稳态温度常分析（续） 27 11 11 传热分析与热应力传热分析与热应力 l例例 稳态温度长分析（续）稳态温度长分析（续） 等效应力分布图等效应力分布图 28 11 11 传热分析与热应力传热分析与热应力 l例例 液固相变分析液固相变分析 有厚度为有厚度为a得液体，初始温度为得液体，初始温度为T0, 在其上表面突加温度在其上表面突加温度TsT0， ，分析液体 分析液体 的温度变化及相变情况，假设液体的的温度变化及相变情况

27、，假设液体的 其他边界是绝热的。其他边界是绝热的。 29 11 11 传热分析与热应力传热分析与热应力 l例例 液固相变分析（续）液固相变分析（续） 相变潜热（相变潜热（latent heat）的影响以焓的快速变化来反映。）的影响以焓的快速变化来反映。 焓（焓（enthalpy）与温度的变化曲线）与温度的变化曲线 30 l例例 液固相变分析（续）液固相变分析（续） 经计算分析，在经计算分析，在789s797s间液体开始完全固化，即节点间液体开始完全固化，即节点2开始降开始降 温至温至1 1 。 在在501s时温度沿厚度的分布时温度沿厚度的分布 各节点温度随时间的变化各节点温度随时间的变化 11

28、 11 传热分析与热应力传热分析与热应力 oC oC 31 11 11 传热分析与热应力传热分析与热应力 l例例 热处理分析热处理分析 将工件放在一定介质中加热到适宜的温度、并在此温度中保持一定将工件放在一定介质中加热到适宜的温度、并在此温度中保持一定 时间后、又以不同速度冷却的一种工艺方法。热处理不改变工件的形时间后、又以不同速度冷却的一种工艺方法。热处理不改变工件的形 状和整体化学成分，而是通过改变工件内部的显微组织，或改变工件状和整体化学成分，而是通过改变工件内部的显微组织，或改变工件 的表面的化学成分，赋予或改变工件的使用性能。的表面的化学成分，赋予或改变工件的使用性能。 特点：特点：

29、1）有热传导）有热传导 2）有相变）有相变 3）热力耦合）热力耦合 32 11 11 传热分析与热应力传热分析与热应力 l例例 热处理分析热处理分析(续续) 淬火分析的网格和淬火分析的网格和 温度分布温度分布 33 2 2 传热分析与热应力传热分析与热应力 34 11 11 传热分析与热应力传热分析与热应力 l热处理分析（续）热处理分析（续） 典型的导热系数曲线典型的导热系数曲线 35 11 11 传热分析与热应力传热分析与热应力 l热处理分析（续）热处理分析（续） 典型的质量密度曲线典型的质量密度曲线 36 11 11 传热分析与热应力传热分析与热应力 l热处理分析（续）热处理分析（续） 典

30、型的比热容曲线典型的比热容曲线 37 11 11 传热分析与热应力传热分析与热应力 l热处理分析（续）热处理分析（续） 典型的相变焓曲线典型的相变焓曲线 38 11 11传热分析与热应力传热分析与热应力 l热处理分析（续）热处理分析（续） 弹性模量曲线弹性模量曲线 39 11 11 传热分析与热应力传热分析与热应力 l热处理分析（续）热处理分析（续） 屈服应力曲线屈服应力曲线 40 11 11 传热分析与热应力传热分析与热应力 l热处理分析（续）热处理分析（续） 淬火实验设备及淬火介质的换热系数曲线淬火实验设备及淬火介质的换热系数曲线 41 11 11 传热分析与热应力传热分析与热应力 l热处

31、理分析（续）几种钢的金相对比热处理分析（续）几种钢的金相对比 42 2 2 传热分析与热应力传热分析与热应力 l热处理分析（续）热处理分析（续） 43 11 11 传热分析与热应力传热分析与热应力 l热处理分析（续）热处理分析（续） 44 11 11 传热分析与热应力传热分析与热应力 l热处理分析（续）热处理分析（续） 45 11 11 传热分析与热应力传热分析与热应力 l焊接残余应力通过改变加辅助热源改变残余应力分布焊接残余应力通过改变加辅助热源改变残余应力分布 46 11 11 传热分析与热应力传热分析与热应力 l焊接残余应力通过改变加辅助热源改变残余应力分布焊接残余应力通过改变加辅助热源

32、改变残余应力分布 47 11 11 传热分析与热应力传热分析与热应力 T T形焊接形焊接 48 11 11 传热分析与热应力传热分析与热应力 l横梁的焊接残余应力焊接现场横梁的焊接残余应力焊接现场 49 11 11 传热分析与热应力传热分析与热应力 l横梁的焊接残余应力有限元模型横梁的焊接残余应力有限元模型 50 11 11 传热分析与热应力传热分析与热应力 l横梁的焊接残余应力温度场分布横梁的焊接残余应力温度场分布 51 11 11 传热分析与热应力传热分析与热应力 l横梁的焊接残余应力沿焊缝方向的残余应力分布横梁的焊接残余应力沿焊缝方向的残余应力分布 52 11 11 传热分析与热应力传热

33、分析与热应力 l横梁的焊接残余应力垂直焊缝方向的残余应力分布横梁的焊接残余应力垂直焊缝方向的残余应力分布 53 11 11 传热分析与热应力传热分析与热应力 本章参考书本章参考书 1.孔祥谦孔祥谦. 有限单元法在传热学中的应用有限单元法在传热学中的应用. 第三版，科学出版社，第三版，科学出版社， 1998 2.王勖成王勖成. 有限单元法有限单元法. 清华大学出版社，清华大学出版社，2003. 第第12章章 3.曾攀曾攀. 有限元分析及应用有限元分析及应用. 清华大学出版社，清华大学出版社，2004 . 第第8章章 54 11 11 传热分析与热应力传热分析与热应力 l术语和单位术语和单位 在国际标准单位制中，传热分析的术语和单位在国际标准单位制中，传热分析的术语和单位 c比热容比热容J/(kgK)T温度（温度（K 或或 C） Q热流（热流（W/m2） h对流换热系数对流换热系数 W/(m2K) t时间（时间（s） k导热系数导热系数W/(mK)质量密度（质量密度（kg/m3） q单位体积热生成率单位体积热生成率 (W/m3) Stefan-Boltzman常数常数 =5.6710-8 W/(m2K4) T /(K /s)Tt 55 11 11 传热分析与热应力传热分析与热应力 11.4伴有相变的导热问题伴有相变