1、6. 稳态热传导问题的有限元法本章的内容如下：6.1 热传导方程与换热边界6.2 稳态温度场分析的一般有限元列式6.3 三角形单元的有限元列式6.4 温度场分析举例6.1 热传导方程与换热边界在分析工程问题时，经常要了解工件内部的温度分布情况，例如发动机的工作温度、金属工件在热处理过程中的温度变化、流体温度分布等。物体内部的温度分布取决于物体内部的热量交换，以及物体与外部介质之间的热量交换，一般认为是与时间相关的。物体内部的热交换采用以下的热传导方程（Fourier 方程）来描述，（6-1）QzTyxTtc 式中 为密度，kg/m 3； 为比热容， ； 为导热系数，cK)J/(kg zyx,；
2、T 为温度，；t 为时间，s ； 为内热源密度，w/m 3。kmw对于各向同性材料，不同方向上的导热系数相同，热传导方程可写为以下形式，（6-2）QzTyxtc222除了热传导方程，计算物体内部的温度分布，还需要指定初始条件和边界条件。初始条件是指物体最初的温度分布情况，（6-3） z,T0t边界条件是指物体外表面与周围环境的热交换情况。在传热学中一般把边界条件分为三类。1）给定物体边界上的温度，称为第一类边界条件。物体表面上的温度或温度函数为已知， sT或 （6-4）),(tzyxs2）给定物体边界上的热量输入或输出，称为第二类边界条件。已知物体表面上热流密度， szzyx qnTnT)(或
3、 （6-5）),()( tzyxzsyx3）给定对流换热条件，称为第三类边界条件。物体与其相接触的流体介质之间的对流换热系数和介质的温度为已知。（6-6）)(sfzzyx ThnTn其中 h 为换热系数，W/(m 2 K)； 是物体表面的温度； 是介质温度。s f如果边界上的换热条件不随时间变化，物体内部的热源也不随时间变化，在经过一定时间的热交换后，物体内各点温度也将不随时间变化，即 0tT这类问题称为稳态（Steady state）热传导问题。稳态热传导问题并不是温度场不随时间的变化，而是指温度分布稳定后的状态，我们不关心物体内部的温度场如何从初始状态过渡到最后的稳定温度场。随时间变化的瞬
4、态（Transient）热传导方程就退化为稳态热传导方程，三维问题的稳态热传导方程为，（6-7）0QzTyxT对于各向同性的材料，可以得到以下的方程，称为 Poisson 方程，（6-8）0zy22Q考虑物体不包含内热源的情况，各向同性材料中的温度场满足 Laplace 方程，（6-9）zTx22在分析稳态热传导问题时，不需要考虑物体的初始温度分布对最后的稳定温度场的影响，因此不必考虑温度场的初始条件，而只需考虑换热边界条件。计算稳态温度场实际上是求解偏微分方程的边值问题。温度场是标量场，将物体离散成有限单元后，每个单元结点上只有一个温度未知数，比弹性力学问题要简单。进行温度场计算时有限单元的
5、形函数与弹性力学问题计算时的完全一致，单元内部的温度分布用单元的形函数，由单元结点上的温度来确定。由于实际工程问题中的换热边界条件比较复杂，在许多场合下也很难进行测量，如何定义正确的换热边界条件是温度场计算的一个难点。6.2 稳态温度场分析的一般有限元列式在前面我们已经介绍了有限元方法可以用来分析场问题，稳态温度场计算是一个典型的场问题。我们可以采用虚功方程建立弹性力学问题分析的有限元格式，推导出的单元刚度矩阵有明确的力学含义。在这里，介绍如何用加权余量法（Weighted Residual Method）建立稳态温度场分析的有限元列式。微分方程的边值问题，可以一般地表示为未知函数 u 满足微
6、分方程组，（在域 内） （6-10）0.)()(21uA未知函数 u 还满足边界条件，（在边界 上） （6-11）0.)()(21B如果未知函数 u 是上述边值问题的精确解，则在域中的任一点上 u 都满足微分方程（6-10） ，在边界的任一点上都满足边界条件（ 6-11） 。对于复杂的工程问题，这样的精确解往往很难找到，需要设法寻找近似解。所选取的近似解是一族带有待定参数的已知函数，一般表示为（6-12）Nainiu1其中 为待定系数， 为已知函数，被称为试探函数。试探函数要取自完全的函数序列，iai是线性独立的。由于试探函数是完全的函数序列，任一函数都可以用这个序列来表示。采用这种形式的近似
7、解不能精确地满足微分方程和边界条件，所产生的误差就称为余量。微分方程（6-10）的余量为，（6-13）)(NaAR边界条件（6-11）的余量为，（6-14）B()选择一族已知的函数，使余量的加权积分为零，强迫近似解所产生的余量在某种平均意义上等于零，（6-15）0ddTjTj RW称为权函数，通过公式（6-15 ）可以选择待定的参数 。jj和 ia这种采用使余量的加权积分为零来求得微分方程近似解的方法称为加权余量法。对权函数的不同选择就得到了不同的加权余量法，常用的方法包括配点法、子域法、最小二乘法、力矩法和伽辽金法（Galerkin method） 。在很多情况下，采用 Galerkin 法
8、得到的方程组的系数矩阵是对称的，在这里也采用 Galerkin 法建立稳态温度场分析的一般有限元列式。在 Galerkin 法中，直接采用试探函数序列作为权函数，取 ， 。jjNWjj下面用求解二阶常微分方程为例，说明 Galerkin 法（参见，王勖成编著“有限元法基本原理和数值方法”的 1.2.3 节） 。例，求解二阶常微分方程 )10(2 xxud边界条件：当 时， ；当 时， 。0xu1x0u取两项近似解： )1(N2x)1()1(221 xaxau，NW由公式（6-15）可以得到两个加权积分方程， 0)62()2()1( 3210 dxxaxax2积分后可以得到一个二元一次方程组，解
9、得， 170.,194.02aa近似解为， )7.94)(xxu该方程的精确解为， 1sin近似解与精确解的结果比较见表 6-1，表 6-1 近似解与精确解比较x=0.25 x=0.5 x=0.75xu1sin0.04401 0.06975 0.06006)170.924.)(x0.04408 0.06944 0.06008假定单元的形函数为， .21nNN单元结点的温度为， TneT.21单元内部的温度分布为， eN以二维问题为例，说明用 Galerkin 法建立稳态温度场的一般有限元格式的过程。二维问题的稳态热传导方程为，（6-16a）0QyTx第一类换热边界为(6-16b)sT第二类换热
10、边界条件为，(6-16c)syxqn第三类边界条件为，（6-16d）)(sfyx ThT在一个单元内的加权积分公式为，（6-17）0)()(1 dQyxwe 由分部积分得， )()()( 111 xTwxTxT)()()( 111 yyyw应用 Green 定理，一个单元内的加权积分公式写为，（6-18）0)()1 111 dnyTxwQwe yxe采用 Galerkin 方法，选择权函数为，iN1将单元内的温度分布函数和换热边界条件代入(6-18)式，单元的加权积分公式为，（6-19）0)()(332 dhTNdThNqQyxfieeiesii eiie 换热边界条件代入后，在(6-19)式
11、内相应出现了第二类换热边界项 ，第dqNsie3三类换热边界项 ，但没有出现与第一类换热边界对应 fieeie 33的项。这是因为，采用 作为权函数，第一类换热边界被自动满足。写成矩阵形式有，iN（6-20）0)()()332 dhTNdTNhqQyxfeee s eTTe 公式（6-20）是 n 个联立的线性方程组，可以确定 n 个结点的温度 。按有限元格式iT将（6-20 ）表示为，（6-21）eePTK其中矩阵K e 为单元的导热矩阵或称为温度刚度矩阵，T e 为单元的结点温度向量，P e称为单元的温度载荷向量或热载荷向量（Thermal load vector） 。对于某个特定单元，单
12、元导热矩阵K e 和温度载荷向量P e 的元素分别为，（6-22） dNhdyNxKjiejiyjixij 3)(（6-23）QhTdqNPiefiesiei 32如果某个单元完全处于物体的内部， dyNxKjiyjixeij )(dQNPiei在整个物体上的加权积分方程是单元积分方程的和，（6-24）0)()()332 dhTNdTNhqdyxfeee sTeTe 根据单元结点的局部编号与整体编号的关系，直接求和得到整体刚度矩阵，整体方程组为， PK6.3 三角形单元的有限元列式图 6-1 三角形单元回顾第三章的内容可以发现，与计算弹性力学平面问题时所采用的方法一样，二维温度场问题计算中所采
13、用的三角形单元可以使用相同的形函数， )(21ycxbaANiiii jjjj )(ycxmm jmijijji xcbyxa ijimjiimj yijjiijji xcyxT1mjjiiyxA2单元内的温度分布用结点上的温度值表示为，（6-25）mjiji TNT在三角形单元上，采用 Galerkin 法可得，（6-26）0)()( dAQyxTAxTNxTNxTxNT )()( 假定单元内的导热系数为常数，（6-27）mjimji jjijiix mjijimjAxTxA Tbb dATbd2224 41yTNyNyTyNT )()( （6-28）mjimji jjijiix mjiji
14、mjAxTyA Tcc dATccd2224 41单元的刚度矩阵为， 222244 mjmi jjijiiymmi jjjiiixe ccbbKj 显然，单元的导热矩阵是对称的。如果单元的内部热源为常数，由内部热源产生的温度载荷项为，（6-29） 13 AQdNdAQNmjiTA由 Green 公式可得（6-30）dSnyTnxTsyxA )( ()方便起见，把换热边界统一表示为第三类换热边界，（6-30）dSTNhdSTNhdSTNhAyyxx esfssfsA )( 如果在单元边上存在热交换，各条边上的边界换热条件在单元刚度矩阵中生成的附加项为，（6-31）0216ijehlK（6-32）
15、21jmel（6-33）06miehlK由边界换热条件生成的温度载荷向量为，（6-34）012ijfelhTP（6-35）1jmfel（6-36）02mifelhTP6.4 温度场分析举例正方形截面的烟囱如图 6-2 所示，烟囱由混凝土建造，边长为 60cm，通道的边长为20cm，混凝土的导热系数为 。假定烟囱内表面的温度为 100，烟囱外)/(4.1KmWk表面暴露在空气中，空气的温度为 30，换热系数为 。计算烟囱截面)/(202KmWh内的稳态温度场。 （参见，Finite Element Method Theory and Application with ANSYS, p279）图 6-2 烟囱截面 图 6-3 有限元模型图 6-4 稳态温度分布图 6-5 热流量分布
