稳态热传导问题有限元法

1、6.稳态热传导问题的有限元法本章的内容如下：6.1热传导方程与换热边界6.2稳态温度场分析的一般有限元列式6.3三角形单元的有限元列式6.4温度场分析举例6.1热传导方程与换热边界在分析工程问题时， 经常要了解工件内部的温度分布情况，例如发动机的工作温度、 金属工件在热处理过程中的温度变化、流体温度分布等。物体内部的温度分布取决于物体内部的热量交换，以及物体与外部介质之间的热量交换，一般认为是与时间相关的。物体内部的热交换采用以下的热传导方程（Fourier方程）来描述，TTTTcxyzQ（ 6-1）t x x y y z z式中 为密度，kg/m3； c为比热容，J/（kg K） ； x,

2、y, z为导热系数， w m k ； T 为温度，C； t为时间，s； Q为内热源密度，w/m3。对于各向同性材料，不同方向上的导热系数相同，热传导方程可写为以下形式，T2T2T2Tc 22 Q（6-2）txyz除了热传导方程，计算物体内部的温度分布，还需要指定初始条件和边界条件。初始条件是指物体最初的温度分布情况，T t 0 T0 x, y,z（6-3）边界条件是指物体外表面与周围环境的热交换情况。在传热学中一般把边界条件分为三类。1）给定物体边界上的温度，称为第一类边界条件。物体表面上的温度或温度函数为已知，Ts Ts或 Ts Ts（x,y,z,t）（6-4）2）给定物体边界上的热量输入或

3、输出，称为第二类边界条件。 已知物体表面上热流密度，(xT nx xTynyyTznz)zsqsTTT 、或(xnxynyznz)sqs(x, y, z,t)（6-5）xyz3) 给定对流换热条件，称为第三类边界条件。物体与其相接触的流体介质之间的对流换热系数和介质的温度为已知。Tx nxxTynyyz 丄 nzh(Tf Ts)z(6-6)其中h为换热系数,W/(m 2 K);Ts是物体表面的温度；Tf是介质温度。如果边界上的换热条件不随时间变化，物体内部的热源也不随时间变化，在经过一定时间的热交换后，物体内各点温度也将不随时间变化，即丄0 t的变化，而是指温度分布稳定后的状态， 到最后的稳定

4、温度场。 随时间变化的瞬态 三维问题的稳态热传导方程为，这类问题称为稳态(Steady state)热传导问题。稳态热传导问题并不是温度场不随时间 我们不关心物体内部的温度场如何从初始状态过渡热传导方程就退化为稳态热传导方程,(Transient)(6-7)对于各向同性的材料,222TTT222xyz可以得到以下的方程，称为(6-8)Poisson 方程,Laplace 方程,考虑物体不包含内热源的情况，各向同性材料中的温度场满足(6-9)2t2t2t02220x y z在分析稳态热传导问题时，不需要考虑物体的初始温度分布对最后的稳定温度场的影 响，因此不必考虑温度场的初始条件，而只需考虑换热

5、边界条件。计算稳态温度场实际上是求解偏微分方程的边值问题。温度场是标量场，将物体离散成有限单元后，每个单元结点上 只有一个温度未知数，比弹性力学问题要简单。进行温度场计算时有限单元的形函数与弹性 力学问题计算时的完全一致，单元内部的温度分布用单元的形函数，由单元结点上的温度来确定。由于实际工程问题中的换热边界条件比较复杂，在许多场合下也很难进行测量，如何定义正确的换热边界条件是温度场计算的一个难点。6.2稳态温度场分析的一般有限元列式在前面我们已经介绍了有限元方法可以用来分析场问题，稳态温度场计算是一个典型的场问题。我们可以采用虚功方程建立弹性力学问题分析的有限元格式，推导出的单元刚度矩阵有明

6、确的力学含义。在这里，介绍如何用加权余量法(Weighted Residual Method )建立稳态温度场分析的有限元列式。微分方程的边值问题，可以一般地表示为未知函数u满足微分方程组，未知函数u还满足边界条件，B(u)B1（u）B2（u）0（在边界上）(6-11)A(u)A(u)A2(u)0（在域 内）(6-10)如果未知函数 u是上述边值问题的精确解，则在域中的任一点上u都满足微分方程（6-10），在边界的任一点上都满足边界条件（6-11 ）。对于复杂的工程问题，这样的精确解往往很难找到，需要设法寻找近似解。所选取的近似解是一族带有待定参数的已知函数， 般表示为n(6-12)u u N

7、 iai Nai 1其中ai为待定系数，Ni为已知函数，被称为试探函数。试探函数要取自完全的函数序列， 是线性独立的。由于试探函数是完全的函数序列，任一函数都可以用这个序列来表示。采用这种形式的近似解不能精确地满足微分方程和边界条件，所产生的误差就称为余量。微分方程（6-10）的余量为，R A（Na）（ 6-13）边界条件（6-11）的余量为，R B（Na）（ 6-14）选择一族已知的函数， 使余量的加权积分为零， 强迫近似解所产生的余量在某种平均意 义上等于零，TT Wj RdWj Rd 0（ 6-15）Wj和Wj称为权函数，通过公式（6-15）可以选择待定的参数ai。这种采用使余量的加权积

8、分为零来求得微分方程近似解的方法称为加权余量法。对权函数的不同选择就得到了不同的加权余量法，常用的方法包括配点法、子域法、最小二乘法、 力矩法和伽辽金法（Galerkin method ）。在很多情况下，采用Galerkin法得到的方程组的系数矩阵是对称的，在这里也采用Galerkin法建立稳态温度场分析的一般有限元列式。在Galerkin法中，直接采用试探函数序列作为权函数，取WjN j，WjNj。下面用求解二阶常微分方程为例，说明Galerkin法（参见，王勖成编著“有限元法基本原理和数值方法”的节）。d2udx2例，求解二阶常微分方程(0 x 1)边界条件：当x 0时，u 0 ;当x 1

9、时，u 0。 取两项近似解：Ni x(1 x)2N2 x (1 x)u N1a1 N2a2a1x(1 x)a2x (1 x)W1N1 ,W2N2由公式(6-15)可以得到两个加权积分方程,x(12x)x 印(2 x x ) a2(2 6x23、x x )dx 0x2(12x)x a1( 2 x x )23a2(2 6x x x )dx 0积分后可以得到一个二元一次方程组，解得,a1 0.1924,a20.1707近似解为，x(1x)( 0.1924 0.1707x)sin x该方程的精确解为，uxsi n1近似解与精确解的结果比较见表6-1，表6-1近似解与精确解比较x=0.25x=0.5x=

10、0.75sin x uxsin 10.044010.069750.06006x(1x)(0.19240.1707x)0.044080.069440.06008假定单元的形函数为，N N12 心单元结点的温度为，Te T1 T2. TnT单元内部的温度分布为,T NTe.维问以二维问题为例，说明用Galerkin法建立稳态温度场的一般有限元格式的过程。题的稳态热传导方程为,(6-16a)第一类换热边界为Tss第二类换热边界条件为,(6-16b)TTx n xy n y q sxy第三类边界条件为，(6-16c)T nx xTnyh(TfTs)y(6-16d)eTTW1 -(x) (y)Qd0xx

11、 yy由分部积分得，T、W1,T、T、(w1x)：(x-)W1( x )xxxxxxT、W1 ,T、T、(w1y)1( y)W1( y)yyyyyy在一个单元内的加权积分公式为,(6-17)应用Green定理，一个单元内的加权积分公式写为,x(-(yT nx x(6-18)Tn y)d0y采用Galerkin方法，选择权函数为,W!Ni将单元内的温度分布函数和换热边界条件代入(6-18)式，单元的加权积分公式为,6 旦(xxN)旦(x yy3)Ted yNiQde2Mqsd(6-19)MhNTedMhTfd 0换热边界条件代入后,在(6-19)式内相应出现了第二类换热边界项e3 Niqsd ，

12、第三3类换热边界项Nih NTed3e3 NihTfd ，但没有出现与第一类换热边界对应的3项。这是因为，采用 M作为权函数，第一类换热边界被自动满足。写成矩阵形式有,e (型)T(Xe tNTQd(6-20)型)(叫x ye t2NTqsd)Ted yTehN NT d3 NThTfd0n个结点的温度Ti。按有限元格式将公式（6-20 ）是n个联立的线性方程组，可以确定（6-20 ）表示为，KeTe Pe(6-21)其中矩阵Ke为单元的导热矩阵或称为温度刚度矩阵, 称为单元的温度载荷向量或热载荷向量（ 热矩阵Ke和温度载荷向量KjNi N-j)dye3 hNiNjd(6-22)MqsdNih

13、Tf d2如果某个单元完全处于物体的内部,NiQd(6-23)Kj(NiNj(xx xNiyyT e为单元的结点温度向量，PeThermal load vector）。对于某个特定单元，单元导P e的元素分别为，NiQd在整个物体上的加权积分方程是单元积分方程的和,宀T(X凹）NTQde t2NTqsd(6-24)3 hNTNTed3 NThTfd0e根据单元结点的局部编号与整体编号的关系，直接求和得到整体刚度矩阵，整体方程组 为,KTP6.3三角形单元的有限元列式(心亠)图6-1三角形单元回顾第三章的内容可以发现，与计算弹性力学平面问题时所采用的方法一样,场问题计算中所采用的三角形单元可以使

14、用相同的形函数，壬(ai bi x Ci y)2 A±(aj bjX Cjy) 2A(am bmX Cmy)2A二维温度NiNjNmaiajXj ymXmyjbiyjymCiXmXjXmYiXiYmbjymyiCjXiXmXiyjXj yibmyiyjCmXjXiXiyiXjyjTXmym2A单元内的温度分布用结点上的温度值表示为,T Ni NjNm TjTm(6-25)在三角形单兀上,采用Galerki n法可得，A NTL (AXx T)X(y _ yT)QdA 0 y(6-26)Ti-)XNt t假定单元内的导热系数为常数,NT TXdAX X14Abi bj bi bmbjb

15、mTi Tj TmdA4Ab2bibjbibmbi bjb2bjbmbibmbjbmbmTiTjTm(6-27)Nt-)y(N y-)yNt tyNT TydA y14ACCj CicmCjCmTiTjTmdA4A CiCjCi CmGC2CjCj CmCi CmC j Cm2CmTiTjTm(6-28)单元的刚度矩阵为,Ke/4bi2bibjbi bj2CiCiCjCi Cmbibjbj2bjbmyGCj2CjCj Cm c；bibmbjbmb；4Ci CmCj CmXNi1T Z-QaQANTQdAAN j dA13Nm1由Green公式可得A (NTXT)(NT y TdAAXXyy-

16、(NT xT nNT y-Tn y)dSsXy显然，单元的导热矩阵是对称的。如果单元的内部热源为常数，由内部热源产生的温度载荷项为,方便起见，把换热边界统一表示为第三类换热边界，(6-29)(6-30)A (NT X 丄)一(NT y 丄dAAxx yy( 6-30)TTTe：hN (Tf Ts)dS ； hN TfdS hN NT dS如果在单元边上存在热交换，各条边上的边界换热条件在单元刚度矩阵中生成的附加项为，hlj210Ke6120(6-31)000000Kehl jm021(6-32)6012201Kehl mi000(6-33)6102由边界换热条件生成的温度载荷向量为,1hTfl

17、jPe21(6-34)00PehT f l jm21(6-35)11ehT f l miP20(6-36)16.4温度场分析举例正方形截面的烟囱如图 6-2所示，烟囱由混凝土建造，边长为60cm，通道的边长为20cm,混凝土的导热系数为 k 1.4W /(m K)。假定烟囱内表面的温度为 100°C,烟囱外表面暴露在空气中，空气的温度为30C,换热系数为h20W/(m2 K)。计算烟囱截面内的稳态温度场。(参见，Finite Element Method Theory and Application with ANSYS, p279)60cjl图6-2烟囱截面图6-3有限元模型G冷0863胡勻 11州"拥图6-4稳态温度分布Ml92.617B5.3r?4ICld图6-5热流量分布稳态温度场分布与物体的初始状态无关，那么是否与材料的导热系数相关？我们把烟囱的模型做些修改，假定烟囱壁由两层材料构成。内层材料为混凝土，外表面的截面尺寸为 30cm 30cm，烟囱通道 的尺寸不 变，仍为20cm 20cm。外 层材料的 导热 系数为 k 0.1W/(m K)，外部表面的截面尺寸不变，内部表面的截面尺寸为30cm 30cm。换热边界条件不变，双层烟囱的有限元模型如图6-6所示。图6-6双层烟囱的有限元模型JTlIiAL gCLEMH 3TTP-J.T