稳态热传导问题的有限元法

6. 稳态热传导问题的有限元法本章的内容如下：6.1热传导方程与换热边界6.2稳态温度场分析的一般有限元列式 6.3三角形单元的有限元列式 6.4温度场分析举例6.1热传导方程与换热边界在分析工程问题时，经常要了解工件内部的温度分布情况，例如发动机的工作温度、金属工件在热处理过程中的温度变化、流体温度分布等。

物体内部的温度分布取决于物体内部的热量交换，以及物体与外部介质之间的热量交换，一般认为是与时间相关的。

物体内部的热交换采用以下的热传导方程（Fourier 方程）来描述，Q z T z y T y x T x t T c+⎪⎭⎫⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂z y x λλλρ （6-1）式中ρ为密度，kg/m 3； c 为比热容，K)J/(kg ⋅；z y x λλλ,,为导热系数，()k m w ⋅；T为温度，℃；t 为时间，s ；Q 为内热源密度，w/m 3。

对于各向同性材料，不同方向上的导热系数相同，热传导方程可写为以下形式，Q zT y T x T t T c 222222+∂∂+∂∂+∂∂=∂∂λλλρ（6-2）除了热传导方程，计算物体内部的温度分布，还需要指定初始条件和边界条件。

初始条件是指物体最初的温度分布情况，() z y,x,T T 00t ==（6-3）边界条件是指物体外表面与周围环境的热交换情况。

在传热学中一般把边界条件分为三类。

1） 给定物体边界上的温度，称为第一类边界条件。

物体表面上的温度或温度函数为已知，s s T T =或),,,(t z y x T T s s =（6-4）2） 给定物体边界上的热量输入或输出，称为第二类边界条件。

已知物体表面上热流密度，s sz z y y x xq n zT n y T n x T =∂∂+∂∂+∂∂)(λλλ或),,,()(t z y x q n zT n y T n x T s sz z y y x x=∂∂+∂∂+∂∂λλλ（6-5）3） 给定对流换热条件，称为第三类边界条件。

物体与其相接触的流体介质之间的对流换热系数和介质的温度为已知。

)(s f z z y y x xT T h n zT n y T n x T -=∂∂+∂∂+∂∂λλλ （6-6）其中h 为换热系数，W/(m 2K)；s T 是物体表面的温度；f T 是介质温度。

如果边界上的换热条件不随时间变化，物体内部的热源也不随时间变化，在经过一定时间的热交换后，物体内各点温度也将不随时间变化，即0=∂∂tT这类问题称为稳态（Steady state ）热传导问题。

稳态热传导问题并不是温度场不随时间的变化，而是指温度分布稳定后的状态，我们不关心物体内部的温度场如何从初始状态过渡到最后的稳定温度场。

随时间变化的瞬态（Transient ）热传导方程就退化为稳态热传导方程，三维问题的稳态热传导方程为，0Q z T z y T y x T x =+⎪⎭⎫⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂z y x λλλ （6-7）对于各向同性的材料，可以得到以下的方程，称为Poisson 方程，0zT y T x T 222222=+∂∂+∂∂+∂∂λQ（6-8）考虑物体不包含内热源的情况，各向同性材料中的温度场满足Laplace 方程，0z Ty T x T 222222=∂∂+∂∂+∂∂ （6-9）在分析稳态热传导问题时，不需要考虑物体的初始温度分布对最后的稳定温度场的影响，因此不必考虑温度场的初始条件，而只需考虑换热边界条件。

计算稳态温度场实际上是求解偏微分方程的边值问题。

温度场是标量场，将物体离散成有限单元后，每个单元结点上只有一个温度未知数，比弹性力学问题要简单。

进行温度场计算时有限单元的形函数与弹性力学问题计算时的完全一致，单元内部的温度分布用单元的形函数，由单元结点上的温度来确定。

由于实际工程问题中的换热边界条件比较复杂，在许多场合下也很难进行测量，如何定义正确的换热边界条件是温度场计算的一个难点。

6.2稳态温度场分析的一般有限元列式在前面我们已经介绍了有限元方法可以用来分析场问题，稳态温度场计算是一个典型的场问题。

我们可以采用虚功方程建立弹性力学问题分析的有限元格式，推导出的单元刚度矩阵有明确的力学含义。

在这里，介绍如何用加权余量法（Weighted Residual Method ）建立稳态温度场分析的有限元列式。

微分方程的边值问题，可以一般地表示为未知函数u 满足微分方程组，0...)()()(21=⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧=u A u A u A（在域Ω内）（6-10）未知函数u 还满足边界条件，0....)()()(21=⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧=u B u B u B（在边界Γ上）（6-11）如果未知函数u 是上述边值问题的精确解，则在域中的任一点上u 都满足微分方程（6-10），在边界的任一点上都满足边界条件（6-11）。

对于复杂的工程问题，这样的精确解往往很难找到，需要设法寻找近似解。

所选取的近似解是一族带有待定参数的已知函数，一般表示为Na ==≈∑=i i ni a N u u 1（6-12）其中i a 为待定系数，i N 为已知函数，被称为试探函数。

试探函数要取自完全的函数序列，是线性独立的。

由于试探函数是完全的函数序列，任一函数都可以用这个序列来表示。

采用这种形式的近似解不能精确地满足微分方程和边界条件，所产生的误差就称为余量。

微分方程（6-10）的余量为， )(Na A R =（6-13）边界条件（6-11）的余量为，B(Na)R =（6-14）选择一族已知的函数，使余量的加权积分为零，强迫近似解所产生的余量在某种平均意义上等于零， 0=Γ+Ω⎰⎰ΓΩd d Tj T j R W R W（6-15）j j W W 和称为权函数，通过公式（6-15）可以选择待定的参数i a 。

这种采用使余量的加权积分为零来求得微分方程近似解的方法称为加权余量法。

对权函数的不同选择就得到了不同的加权余量法，常用的方法包括配点法、子域法、最小二乘法、力矩法和伽辽金法（Galerkin method ）。

在很多情况下，采用Galerkin 法得到的方程组的系数矩阵是对称的，在这里也采用Galerkin 法建立稳态温度场分析的一般有限元列式。

在Galerkin 法中，直接采用试探函数序列作为权函数，取j j N W =，j j N W -=。

下面用求解二阶常微分方程为例，说明Galerkin 法（参见，王勖成编著“有限元法基本原理和数值方法”的1.2.3节）。

例，求解二阶常微分方程)10(022≤≤=++x x u dxud边界条件：当0=x 时，0=u ；当1=x 时，0=u 。

取两项近似解： )1(1x x N -= )1(22x x N -=)1()1(~2212211x x a x x a a N a N u -+-=+= 11N W =， 22N W =由公式（6-15）可以得到两个加权积分方程，0)]62()2()[1(3222110=-+-+-+-+-⎰dx x x x a x x a x x x0)]62()2()[1(3222121=-+-+-+-+-⎰dx x x x a x x a x x x积分后可以得到一个二元一次方程组，解得，1707.0,1924.021==a a近似解为，)1707.01924.0)(1(~x x x u+-= 该方程的精确解为，x xu -=1sin sin近似解与精确解的结果比较见表6-1，假定单元的形函数为，]...[][21n N N N N =单元结点的温度为，T n e T T T T ]...[}{21=单元内部的温度分布为，e T N T }]{[=以二维问题为例，说明用Galerkin 法建立稳态温度场的一般有限元格式的过程。

二维问题的稳态热传导方程为，0Q y T y x T x =+⎪⎪⎭⎫⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂y x λλ （6-16a ）第一类换热边界为s sT T=(6-16b)第二类换热边界条件为，s y y x xq n yTn x T =∂∂+∂∂λλ (6-16c)第三类边界条件为，)(s f y y x xT T h n yT n x T -=∂∂+∂∂λλ （6-16d ）在一个单元内的加权积分公式为，0])~()~([1=Ω+∂∂∂∂+∂∂∂∂⎰Ωd Q yT y x T x w y x eλλ（6-17）由分部积分得，)~()~()~(111x Tx w x T x w x T w x x x x ∂∂∂∂+∂∂∂∂=∂∂∂∂λλλ)~()~()~(111yTy w y T y w y T w y y y y ∂∂∂∂+∂∂∂∂=∂∂∂∂λλλ 应用Green 定理，一个单元内的加权积分公式写为，)~~(])~()~([1111=Γ∂∂+∂∂+Ω-∂∂∂∂+∂∂∂∂-⎰⎰ΓΩd n y T n x T w d Q w yT y w x T x w y y x x ey x eλλλλ （6-18）采用Galerkin 方法，选择权函数为，i N w =1将单元内的温度分布函数和换热边界条件代入(6-18)式，单元的加权积分公式为，}]{[})]{][()][([332=Γ-Γ+Γ-Ω-Ω∂∂∂∂+∂∂∂∂⎰⎰⎰⎰⎰ΓΓΓΩΩd hT N d T N h N d q N d Q N d T yN y N x N x N f i eei ee s i i ee y i x i eλλ （6-19）换热边界条件代入后，在(6-19)式内相应出现了第二类换热边界项Γ-⎰Γd q N s i e3，第三类换热边界项Γ-Γ⎰⎰ΓΓd hT N d T N h N f i eei e33}]{[，但没有出现与第一类换热边界对应的项。

这是因为，采用i N 作为权函数，第一类换热边界被自动满足。

写成矩阵形式有，][}]{[][][][})]{][()][()][()][[(332=Γ-Γ+Γ-Ω-Ω∂∂∂∂+∂∂∂∂⎰⎰⎰⎰⎰ΓΓΓΩΩd hT N d T N N h d q N d Q N d T yN y N x N x N f T eeTee s T T ee y T x T eλλ （6-20）公式（6-20）是n 个联立的线性方程组，可以确定n 个结点的温度i T 。

按有限元格式将（6-20）表示为，e e e P T K }{}{][=（6-21）其中矩阵[K]e为单元的导热矩阵或称为温度刚度矩阵，{T}e为单元的结点温度向量，{P}e称为单元的温度载荷向量或热载荷向量（Thermal load vector ）。