第6章热传导问题有限元.pdf

187第六章第六章热传导问题有限元热传导问题有限元§§6- 1 引言引言一、温度应力和温度场在变温条件下工作的结构和部件，通常都存在温度应力，有的是稳定的温度应力，有的是随时间变化的瞬态温度应力 这些应力在结构应力中经常占有相当的比重， 甚至成为设计结构或部件的控制应力要计算这些应力首先要确定结构或构件工作所在的稳态或瞬态的温度场由于结构的形状以及变温条件的复杂性，依靠传统的解析方法要精确地确定温度场往往是不可能的，有限单元法是解决上述问题的方便和有效的工具稳态或瞬态温度场问题，即稳态或瞬态热传导问题，在空间域的离散与前面几章中讨论的弹性力学问题类似，采用0C型的插值函数，前面弹性力学问题中讨论过的单元和相应的位移模式在这里都可以使用主要的不同在于场变量，在弹性力学问题中场变量是位移，是向量场；在热传导问题中场变量是温度，是标量场对于瞬态温度场，除了空间域的离散外还有时间域的离散二、热传导微分方程及边界条件（1）传导定义：传导是在没有任何材料质量纯运动的情况下，热通过材料的传递沿x方向上传导的热流速率为xTAkqx（6- 1- 1）其中xk是x方向上材料的导热系数；A是垂直于x方向热流通过的面积；T是温度。

2）对流定义：对流是固体与周围物体之间进行热能传递的过程对流的热流速率可表示为 TThAq（6- 1- 2）其中h是热传导系数；A是热流通过物体表面的面积；T是物体表面的温度；T是环境介质的温度3）辐射定义：辐射热传导是在服从电磁学定律的两个表面之间的热能交换过程辐射热流速率由下述关系确定 TTAq （6- 1- 3）其中是斯特芬－波尔兹曼（Stefen- Baltzmann）常数；是表面的放射率；A是热流通过物体表面的面积；T是物体表面的绝对温度；T是环境介质的绝对温度4）固体中产生的能量当其他形式的能量如化学能、核能或电能转换为热能时，在固体中就会产生能量生成热的速率由下述方程确定VqEg（6- 1- 4）188其中q 是热流的强度（单位时间单位体积生成热的速率） ；V是物体的体积5）固体中储存的能量当固体中的温度增加时，热能将会储存在固体中描述这种现象的方程为tTcVEs（6- 1- 5）其中sE是固体中能量储存的速率；是材料的密度；c是材料的比热；V是物体的体积；T是物体的温度；t是时间参数分析物体微元，由能量平衡可得三维物体中热传导的控制微分方程：tTcqzTkzyTkyxTkxzyx           （6- 1- 6）上式是控制正交各向异性体中的热传导微分方程。

如果假设x，y和z方向的热传导率相同，即kkkkzyx，则上式可写成tT kq zT yT xT  1222222（6- 1- 7）其中，常数ck称为放热系数如果物体中没有热源，方程（6- 1- 7）可进一步简化为傅里叶方程tT zT yT xT  1222222 （6- 1- 8）如果物体处于稳定状态（有热源） ，则方程（6- 1- 7）可简化为泊松方程0222222  kq zT yT xT（6- 1- 9）如果物体处于没有任何热源的稳定状态，则方程（6- 1- 7）可简化为拉普拉斯方程0222222  zT yT xT（6- 1- 10）由于微分方程（6- 1- 6）或（6- 1- 7）是二阶的，所以需要规定两个边界条件可能的边界条件是在1上：0,,,TtzyxT（6- 1- 11a）在2上：0qlzTklyTklxTkzzyyxx（6- 1- 11b）在3上：TThlzTklyTklxTkzzyyxx （6- 1- 11c）其中，q是边界上的热流；h是对流热传导系数；zyxlll,,是垂直于边界向外的方向余弦；1是温度值规定为 tT0的边界；2是热流量规定为q的边界；3是对流热损耗规定为TTh的边界。

此外，由于微分方程（6- 1- 6）或（6- 1- 7）在时间t内是一阶的，因而要求一个初始条件，通常所189用的初始条件为在V内：zyxTtzyxT,,0,,,0（6- 1- 11d）其中，V表示固体的区域（或体积） ；0T表示在时间为零时所规定的温度分布求固体内温度分布的问题， 就是在满足边界条件方程（6- 1- 11a）至 （6- 1- 11c）及初始条件 （6- 1- 11d）的情况下，解方程（6- 1- 6）或（6- 1- 7） 为了求解上述热传导问题的数值解，可以将三维热传导问题用下述等价的变分形式来描述，然后导出求解此类问题的有限元方程求固体内的温度分布tzyxT,,,，该分布应使下述泛函取极小值             VzyxdVT tTcqzTkyTkxTk~222 221（6- 1- 12）并应满足边界条件方程（6- 1- 11a）至（6- 1- 11c）及初始条件（6- 1- 11d） 可以证明方程（6- 1- 6）是对应于泛函（6- 1- 12）的欧拉－拉格郎日方程。

一般在假定温度分布时，不难满足边界条件（6- 1- 11a） ，而满足边界条件（6- 1- 11b）和（6- 1- 11c）会有困难为了克服这种困难，可以把边界条件（6- 1- 11b）和（6- 1- 11c）的积分形式引入泛函式（6- 1- 12）中，组成新的泛函方程（6- 1- 11b）和（6- 1- 11c）的积分如下322 2321dSTThqTdSSS  因此，组成的新泛函为322~2222321221dSTThqTdSdVT tTcqzTkyTkxTkSSVzyx               （6- 1- 13）§§6- 2有限元方程的推导有限元方程的推导步骤1：将空间域离散成E个有限单元体，每个单元有p个节点步骤2：在每个单元内，各点的温度T可以近似地用单元的节点温度插值表示为 piee iiTNTzyxNT1)(,,（6- 2- 1）步骤3：由于近似场函数是构造在单元中的，因此泛函（6- 1- 13）式可改写为对单元积分的总和。

32 )( 2)()( ~)(2)(2)(2)(1)( 2)( 321221dSTThdSqTdVT tTcqzTkyTkxTkeeeSSeeVeeezeyexEee                     （6- 2- 2）为了使泛函取极值，利用必要条件190 EeieiTT1)( 0，（Mi,, 2 , 1）其中：M是节点温度未知数的总个数将泛函式（6- 2- 2）代入上式，并由式（6- 2- 1） ，对于每个单元可得   0)()()( 3)()( 2)()( 1)()( eeeeeee ee PTKTKTKT（6- 2- 3）其中：)( 1eK表示各单元对热传导矩阵的贡献，)( 2eK表示热交换边界条件对热传导矩阵的修正，)( 3eK表示非稳态导致的附加项，称为单元的热容量矩阵，)(eP表示温度载荷列阵它们的元素可由以下各式给出      )()( 1eVji zji yji xe ijdVzNzNkyNyNkxNxNkK（6- 2- 3a） )( 33)( 2eSjie ijdSNhNK（6- 2- 3b）  )()( 3eVjie ijdVNcNK（6- 2- 3c）    )()( 2)( 332)(eeeVSSiie idShTdSqNdVNqP（6- 2- 3d）式（6- 2- 3d）中的三项分别为热源、给定热流（第二边界条件）和热交换（第三边界条件）引起的温度载荷。

步骤 4：根据式（6- 2- 2）组集各单元，可以得到  PTKTK.3（6- 2- 4）其中  EeeKK1)( 33； EeeeKKK1)( 2)( 1； EeePP1)(式中T是系统所有节点温度未知数向量；P是系统所有节点的温度载荷向量步骤 5：式（6- 2- 4）就是热传导问题的有限元方程，当引入在1上的边界条件式（6- 1- 11a）和初始条件式（6- 1- 11d）后，就可求解方程式（6- 2- 4） 以下以最简单而十分有用的 3 节点三角形单元为例， 推导二维稳态热传导和瞬态热传导的有限元列式§§6- 3稳态二维热传导稳态二维热传导根据有限元部分的§2- 1 节的第（2- 1- 2a）式，3 节点有限元的插值函数为ycxbaANiiii21（mji,,）对于任一单元ijm，可将插值函数求导代入式（6- 2- 3a） ，得到热传导矩阵元素191jiy jixe ijccAkbbAkK44)( 1（6- 3- 1）单元热传导矩阵为    mmmjjjmijiii ymmmjjjmijiii xeccsymccccccccccAkbbsymbbbbbbbbbbAkK44)( 1（6- 3- 2）对于具有第三类边界条件的单元，如rsp单元，除按式（6- 3- 2）计算单元热传导矩阵外，还应计算由于第三类边界条件引起的对热传导矩阵的修正。

修正项可将插值函数代入式（6- 2- 3b）得到hLdlNhNKKhLdlNhNKKlrre rre sslsre rse sr3161)( 2)( 2)( 2)( 2（6- 3- 3）式中L是对流边界sr 的边长如单元只有sr边为对流换热边界，则单元对热传导矩阵的修正是   0000210126)( 2hLKe（6- 3- 4）单元的温度载荷可由式（6- 2- 3d）求得，假设单元厚度为 1        1113.0.3210.)( 1)(AVTeAqdA LLL qdVNqPe（6- 3- 4a）如果ij边位于对流边界2上，则有033LN，dstdsdS2，从而     0112021 )( 2ijssseqsdsLLqPji（6- 3- 4b）与此类似，如果ij边位于流交换边界3上，可得     0112021 )( 3ijssseshTdsLL hTPji（6- 3- 4c）所以，二维结构稳态热传导的有限元求解方程为 PTK（6- 3- 5）其中热传导矩阵K由各单元的子矩阵式（6- 3- 2）和（6- 3- 4）叠加而成；温度载荷向量P由各单元的温度载荷式（6- 3- 4a） 、 （6- 3- 4b）和（6- 3- 4c）叠加而成。

192§§6- 4瞬态二维热传导瞬态二维热传导在热传导中，时间相关或非稳态问题是很普遍的这种瞬态热传导问题的控制方程由方程（6- 1- 6）给定，有关的边界条件和初始条件由式（6- 1- 11a）和（6- 1- 11d）给定全部参数zy。