清华大学：《有限元分析》课程教学资源（电子讲义）第六章 稳态热传导问题的有限元法

6.稳态热传导问题的有限元法 本章的内容如下: 6热传导方程与换热边界 62稳态温度场分析的一般有限元列式 6.3三角形单元的有限元列式 64温度场分析举例 61热传导方程与换热边界 在分析工程问题时,经常要了解工件内部的温度分布情况,例如发动机的工作温度、金 属工件在热处理过程中的温度变化、流体温度分布等。物体内部的温度分布取决于物体内部 的热量交换,以及物体与外部介质之间的热量交换,一般认为是与时间相关的。物体内部的 热交换采用以下的热传导方程( Fourier方程)来描述, T pc ot ox r ox)ayl"yay)az 式中P为密度,kgm;c为比热容,JkgK):,,为导热系数,w/(mk):T 为温度,℃;t为时间,s;Q为内热源密度,w/m3 对于各向同性材料,不同方向上的导热系数相同,热传导方程可写为以下形式 o=dT, aT 除了热传导方程,计算物体内部的温度分布,还需要指定初始条件和边界条件。初始条 件是指物体最初的温度分布情况, T=0=To(xy, z) (63) 边界条件是指物体外表面与周围环境的热交换情况。在传热学中一般把边界条件分为三 1)给定物体边界上的温度,称为第一类边界条件。 物体表面上的温度或温度函数为已知, T 或7,=7,(x,y1) (6-4) 2)给定物体边界上的热量输入或输出,称为第二类边界条件。 已知物体表面上热流密度, h+2O7 n+2 =4s 2n),=9,(xy=D

6. 稳态热传导问题的有限元法 本章的内容如下： 6.1 热传导方程与换热边界 6.2 稳态温度场分析的一般有限元列式 6.3 三角形单元的有限元列式 6.4 温度场分析举例 6.1 热传导方程与换热边界 在分析工程问题时，经常要了解工件内部的温度分布情况，例如发动机的工作温度、金 属工件在热处理过程中的温度变化、流体温度分布等。物体内部的温度分布取决于物体内部 的热量交换，以及物体与外部介质之间的热量交换，一般认为是与时间相关的。物体内部的 热交换采用以下的热传导方程（Fourier 方程）来描述， Q z T y z T x y T t x T c  +          +             +          =     x  y  z （6-1） 式中  为密度，kg/m3； c 为比热容， J/(kg  K) ；  x  y  z , , 为导热系数， w (m k) ；T 为温度，℃；t 为时间，s； Q 为内热源密度，w/m3。 对于各向同性材料，不同方向上的导热系数相同，热传导方程可写为以下形式， Q z T y T x T t T c 2 2 2 2 2 2 +   +   +   =       （6-2） 除了热传导方程，计算物体内部的温度分布，还需要指定初始条件和边界条件。初始条 件是指物体最初的温度分布情况， T T (x, y, z) t=0 = 0 （6-3） 边界条件是指物体外表面与周围环境的热交换情况。在传热学中一般把边界条件分为三 类。 1）给定物体边界上的温度，称为第一类边界条件。 物体表面上的温度或温度函数为已知， s s T = T 或 T T (x, y,z,t) s s = （6-4） 2）给定物体边界上的热量输入或输出，称为第二类边界条件。 已知物体表面上热流密度， s s x x y y z nz q z T n y T n x T =   +   +   (   ) 或 ( n ) q (x, y,z,t) z T n y T n x T s s x x y y z z =   +   +      （6-5）

3)给定对流换热条件,称为第三类边界条件 物体与其相接触的流体介质之间的对流换热系数和介质的温度为已知 ,+,0n,+2n=h(T-T,) (6-6) 其中h为换热系数,W/(m2K):T是物体表面的温度;T是介质温度。 如果边界上的换热条件不随时间变化,物体内部的热源也不随时间变化,在经过一定时 间的热交换后,物体内各点温度也将不随时间变化,即 T 这类问题称为稳态( Steady state)热传导问题。稳态热传导问题并不是温度场不随时间 的变化,而是指温度分布稳定后的状态,我们不关心物体内部的温度场如何从初始状态过渡 到最后的稳定温度场。随时间变化的瞬态( Transient)热传导方程就退化为稳态热传导方程, 三维问题的稳态热传导方程为, (a),(,( +Q=0 对于各向同性的材料,可以得到以下的方程,称为 Poisson方程 aT...o (6-8) 考虑物体不包含内热源的情况,各向同性材料中的温度场满足 Laplace方程 OT OT aT (6-9) 在分析稳态热传导问题时,不需要考虑物体的初始温度分布对最后的稳定温度场的影 响,因此不必考虑温度场的初始条件,而只需考虑换热边界条件。计算稳态温度场实际上是 求解偏微分方程的边值问题。温度场是标量场,将物体离散成有限单元后,每个单元结点上 只有一个温度未知数,比弹性力学问题要简单。进行温度场计算时有限单元的形函数与弹性 力学问题计算时的完全一致,单元内部的温度分布用单元的形函数,由单元结点上的温度来 确定。由于实际工程问题中的换热边界条件比较复杂,在许多场合下也很难进行测量,如何 定义正确的换热边界条件是温度场计算的一个难点。 62稳态温度场分析的一般有限元列式 在前面我们已经介绍了有限元方法可以用来分析场问题,稳态温度场计算是一个典型的 场问题。我们可以采用虚功方程建立弹性力学问题分析的有限元格式,推导出的单元刚度矩 阵有明确的力学含义。在这里,介绍如何用加权余量法( Weighted Residual Method)建立稳 态温度场分析的有限元列式。 微分方程的边值问题,可以一般地表示为未知函数u满足微分方程组

3）给定对流换热条件，称为第三类边界条件。 物体与其相接触的流体介质之间的对流换热系数和介质的温度为已知。 ( ) x x y y z nz h Tf Ts z T n y T n x T = −   +   +      （6-6） 其中 h 为换热系数，W/(m2 K)； Ts 是物体表面的温度； Tf 是介质温度。 如果边界上的换热条件不随时间变化，物体内部的热源也不随时间变化，在经过一定时 间的热交换后，物体内各点温度也将不随时间变化，即 = 0   t T 这类问题称为稳态（Steady state）热传导问题。稳态热传导问题并不是温度场不随时间 的变化，而是指温度分布稳定后的状态，我们不关心物体内部的温度场如何从初始状态过渡 到最后的稳定温度场。随时间变化的瞬态（Transient）热传导方程就退化为稳态热传导方程， 三维问题的稳态热传导方程为， Q 0 z T y z T x y T x  + =          +             +           x  y  z （6-7） 对于各向同性的材料，可以得到以下的方程，称为 Poisson 方程， 0 z T y T x T 2 2 2 2 2 2 + =   +   +    Q （6-8） 考虑物体不包含内热源的情况，各向同性材料中的温度场满足 Laplace 方程， 0 z T y T x T 2 2 2 2 2 2 =   +   +   （6-9） 在分析稳态热传导问题时，不需要考虑物体的初始温度分布对最后的稳定温度场的影 响，因此不必考虑温度场的初始条件，而只需考虑换热边界条件。计算稳态温度场实际上是 求解偏微分方程的边值问题。温度场是标量场，将物体离散成有限单元后，每个单元结点上 只有一个温度未知数，比弹性力学问题要简单。进行温度场计算时有限单元的形函数与弹性 力学问题计算时的完全一致，单元内部的温度分布用单元的形函数，由单元结点上的温度来 确定。由于实际工程问题中的换热边界条件比较复杂，在许多场合下也很难进行测量，如何 定义正确的换热边界条件是温度场计算的一个难点。 6.2 稳态温度场分析的一般有限元列式 在前面我们已经介绍了有限元方法可以用来分析场问题，稳态温度场计算是一个典型的 场问题。我们可以采用虚功方程建立弹性力学问题分析的有限元格式，推导出的单元刚度矩 阵有明确的力学含义。在这里，介绍如何用加权余量法（Weighted Residual Method）建立稳 态温度场分析的有限元列式。 微分方程的边值问题，可以一般地表示为未知函数 u 满足微分方程组

A1(a) A(u)={A2(u)}=0 (在域Ω内) (6-10) 未知函数u还满足边界条件, B,(u B()={B2()}=0 (在边界r上) (6-11) 如果未知函数u是上述边值问题的精确解,则在域中的任一点上u都满足微分方程 (6-10),在边界的任一点上都满足边界条件(6-11)。对于复杂的工程问题,这样的精确解 往往很难找到,需要设法寻找近似解。所选取的近似解是一族带有待定参数的已知函数, 般表示为 l≈l=Nan=Na (6-12) 其中a,为待定系数,N为已知函数,被称为试探函数。试探函数要取自完全的函数序列 是线性独立的。由于试探函数是完全的函数序列,任一函数都可以用这个序列来表示。 采用这种形式的近似解不能精确地满足微分方程和边界条件,所产生的误差就称为余 微分方程(6-10)的余量为, R=A(Na) (6-13) 边界条件(6-11)的余量为, R=B(Na 选择一族已知的函数,使余量的加权积分为零,强迫近似解所产生的余量在某种平均意 义上等于零 W RdQ2+ w RdT=0 (6-15) W和W称为权函数,通过公式(615)可以选择待定的参数a1 这种采用使余量的加权积分为零来求得微分方程近似解的方法称为加权余量法。对权函 数的不同选择就得到了不同的加权余量法,常用的方法包括配点法、子域法、最小二乘法、 力矩法和伽辽金法( Galerkin method)。在很多情况下,采用 Galerkin法得到的方程组的系 数矩阵是对称的,在这里也采用 Galerkin法建立稳态温度场分析的一般有限元列式。在 Galerkin法中,直接采用试探函数序列作为权函数,取H=N,W=-Ny 下面用求解二阶常微分方程为例,说明 Galerkin法(参见,王勖成编著“有限元法基本 原理和数值方法”的1.23节) 例,求解二阶常微分方程 d 2u +l+x=0 ≤1)

0 ... ( ) ( ) ( ) 2 1 =           = A u A u A u （在域  内） （6-10） 未知函数 u 还满足边界条件， 0 .... ( ) ( ) ( ) 2 1 =           = B u B u B u （在边界  上） （6-11） 如果未知函数 u 是上述边值问题的精确解，则在域中的任一点上 u 都满足微分方程 （6-10），在边界的任一点上都满足边界条件（6-11）。对于复杂的工程问题，这样的精确解 往往很难找到，需要设法寻找近似解。所选取的近似解是一族带有待定参数的已知函数，一 般表示为  =  = Na = i i n i u u N a 1 （6-12） 其中 i a 为待定系数， Ni 为已知函数，被称为试探函数。试探函数要取自完全的函数序列， 是线性独立的。由于试探函数是完全的函数序列，任一函数都可以用这个序列来表示。 采用这种形式的近似解不能精确地满足微分方程和边界条件，所产生的误差就称为余 量。 微分方程（6-10）的余量为， R = A(Na) （6-13） 边界条件（6-11）的余量为， R = B(Na) （6-14） 选择一族已知的函数，使余量的加权积分为零，强迫近似解所产生的余量在某种平均意 义上等于零， +  = 0   d d T j T Wj R W R （6-15） Wj和Wj 称为权函数，通过公式（6-15）可以选择待定的参数 i a 。 这种采用使余量的加权积分为零来求得微分方程近似解的方法称为加权余量法。对权函 数的不同选择就得到了不同的加权余量法，常用的方法包括配点法、子域法、最小二乘法、 力矩法和伽辽金法（Galerkin method）。在很多情况下，采用 Galerkin 法得到的方程组的系 数矩阵是对称的，在这里也采用 Galerkin 法建立稳态温度场分析的一般有限元列式。在 Galerkin 法中，直接采用试探函数序列作为权函数，取 Wj = N j ，Wj = −Nj 。 下面用求解二阶常微分方程为例，说明 Galerkin 法（参见，王勖成编著“有限元法基本 原理和数值方法”的 1.2.3 节）。 例，求解二阶常微分方程 0 (0 1) 2 2 + u + x =  x  dx d u

边界条件:当x=0时,=0:当x=1时,u=0。 取两项近似解: N2=x2(1-x) l=Na1+N2a2=a1x(1-x)+a2x2(1-x) W2=N2 由公式(6-15)可以得到两个加权积分方程, x(1-x)x+a1(-2+x-x2)+a2(2-6x+x2-x3)ar=0 x(-xx+a(-2+x-x)+a(2-6x+x2-x)=0 积分后可以得到一个二元一次方程组,解得, a1=0.1924,a2=0.1707 近似解为,=x(1-x)(0.1924+0.1707x) 该方程的精确解为,u sin x 近似解与精确解的结果比较见表6-1, 表6-1近似解与精确解比较 x=0.75 x 0.04401 0.06975 0.06006 =x(1-x)0.1924+0.1707x)0048 0.06944 0.06008 假定单元的形函数为, [N]=[M1N2….Nn 单元结点的温度为, {T}°=[T1T2…Tn 单元内部的温度分布为, T=INITI 以二维问题为例,说明用 Galerkin法建立稳态温度场的一般有限元格式的过程。二维问 题的稳态热传导方程为 (“丿()+0=0 (6-16a 第一类换热边界为

边界条件：当 x = 0 时， u = 0 ；当 x =1 时， u = 0。 取两项近似解： (1 ) 1 N = x − x (1 ) 2 2 N = x − x (1 ) (1 ) ~ 2 1 1 2 2 1 2 u = N a + N a = a x − x + a x − x W1 = N1， W2 = N2 由公式（6-15）可以得到两个加权积分方程， (1 )[ ( 2 ) (2 6 )] 0 2 3 2 2 1 1 0 − + − + − + − + − =  x x x a x x a x x x dx (1 )[ ( 2 ) (2 6 )] 0 2 3 2 2 1 2 1 0 − + − + − + − + − =  x x x a x x a x x x dx 积分后可以得到一个二元一次方程组，解得， a1 = 0.1924, a2 = 0.1707 近似解为， (1 )(0.1924 0.1707 ) ~ u = x − x + x 该方程的精确解为， x x u = − sin 1 sin 近似解与精确解的结果比较见表 6-1， 表 6-1 近似解与精确解比较 x=0.25 x=0.5 x=0.75 x x u = − sin 1 sin 0.04401 0.06975 0.06006 (1 )(0.1924 0.1707 ) ~ u = x − x + x 0.04408 0.06944 0.06008 假定单元的形函数为， [ ] [ ... ] N = N1 N2 Nn 单元结点的温度为， T n e {T} [T T ... T ] = 1 2 单元内部的温度分布为， e T = [N]{T} 以二维问题为例，说明用 Galerkin 法建立稳态温度场的一般有限元格式的过程。二维问 题的稳态热传导方程为， Q 0 y T x y T x + =              +           x  y （6-16a） 第一类换热边界为

T=T 第二类换热边界条件为 +1 n= 第三类边界条件为, A-n,+a,a-n,=h(t-t, (6-16d) 在一个单元内的加权积分公式为, )+(,)+Q1s2=0 (6-17) 由分部积分得, aT 1)=(4a)+"ax(ax 应用 Green定理,一个单元内的加权积分公式写为 Q ayay f w a, an,+a, a n, dr=0 采用 Galerkin方法,选择权函数为, 将单元内的温度分布函数和换热边界条件代入(6-18)式,单元的加权积分公式为, aN2o[N]、,aN/oN ){T}°c2 N, OdQ2-.N, q,dr (6-19) S, N mnity'dr-N hT dr=0 换热边界条件代入后,在(6-19)式内相应出现了第二类换热边界项-Nqd,第三 类换热边界项NAMT-与1NTdm,但没有出现与第一类换热边界对应的 项。这是因为,采用N作为权函数,第一类换热边界被自动满足。写成矩阵形式有

s s T = T (6-16b) 第二类换热边界条件为， x x y ny qs y T n x T =   +     (6-16c) 第三类边界条件为， ( ) x x y ny h Tf Ts y T n x T = −   +     （6-16d） 在一个单元内的加权积分公式为， ) ] 0 ~ ) ( ~ [ ( 1 +  =     +      Q d y T x y T x w x y e   （6-17） 由分部积分得， ) ~ ) ( ~ ) ( ~ ( 1 1 1 x T x w x T x w x T w x x x x     +     =        ) ~ ) ( ~ ) ( ~ ( 1 1 1 y T y w y T y w y T w y y y y     +     =        应用 Green 定理，一个单元内的加权积分公式写为， ) 0 ~ ~ ( ) ] ~ ) ( ~ [ ( 1 1 1 1  =   +   + −      +     −     n d y T n x T w w Q d y T y w x T x w x x y y e x y e     （6-18） 采用 Galerkin 方法，选择权函数为， w1 = Ni 将单元内的温度分布函数和换热边界条件代入(6-18)式，单元的加权积分公式为， [ ]{ } 0 )]{ } [ ] ) ( [ ] [ ( 3 3 2 +  −  = −  −       +               N h N T d N hT d N Qd N q d T d y N y N x N x N i f e e i e e i i s e e y i x i e   （6-19） 换热边界条件代入后，在(6-19)式内相应出现了第二类换热边界项−   Niqsd e 3 ，第三 类换热边界项  −    N h N T d NihTf d e e i e 3 3 [ ]{ } ，但没有出现与第一类换热边界对应的 项。这是因为，采用 Ni 作为权函数，第一类换热边界被自动满足。写成矩阵形式有