求解导热问题的有限单元法

1、第四章 求解导热问题的有限单元法第4.1节 概述第4.2节 泛函变分原理第4.3节 有限单元法第4.1节 概述粗略地讲：有限元法是获得微分方程近似解的一种方法，是一种适合计算机来求解的数值计算方法。（元素特性方程和总体合成方程的建立可以采用直接法，变分法，加权余数法和能量平衡法等四种方法之一，所以粗略地说有限元法是获得微分方程近似解的一种方法也有道理）比较严格的定义：有限单元法是求解泛函变分问题的一种近似方法。那么这两种说法有什么联系，或者说是共同之处呢？变分和微分是对未知函数的不同描述，同一连续介质问题往往都可以找到微分和变分的等价表达方式。变分和微分几乎是同时发展起来的两个数学分支，其目的

2、是相同的，都是求解未知函数，但是方法上有很大差别。在已知边界条件的情况下，求微分方程的精确解析虽然已有完整的理论，但是真正能解出的只有极少数的几种简单情况，因为在很多情况下，微分方程并不存在初等函数解析解。（对于各种各样的映射，初等函数的表达能力实在太有限了，初等函数包括：冥函数、指数函数、对数、三角函数，以及它们的四则运算等。）由于寻求微分方程的初等函数解析解有困难，所以我们在前一章讲述了微分方程的近似解法，即差分法。泛函变分原理虽然也可以用解析法（即积分）求得未知函数，但是因为有很多被积函数根本无法找到初等原函数，也就不能积分，尤其是对于二维和三维问题，解析法更加困难。所以我们也要寻求泛函

3、变分的近似解法。泛函变分的近似解法包括里兹法和有限元法（里兹法是有限元法的前身），这两种方法的原理完全相同，即：构造一个近似的初等函数，用近似的初等函数去逼近未知函数。因为任何未知函数都可以找到它的近似初等函数（如：包含待定系数的多项式或三角函数），所以从根本上克服了解析法（无法找到初等原函数）的局限性牺牲极小的理论计算精度，却换回了对大量复杂二维和三维工程问题的适用性。微分方程的近似解法：差分法泛函变分的近似解法：里兹法，有限元法第4.2节 泛函变分原理一、泛函的概念（借助讲解）二、变分的概念 借助普通函数微分的概念，用类比法讲解三、泛函的极值条件 借助普通函数的极值条件，用类比法讲解四、里

4、兹法（补充内容，但是很重要）泛函变分的近似解法一、泛函的概念通过教材4.1.2的运算实例引入一个泛函表达方式（作图讲解）泛函的概念：函数的函数泛函与普通函数的区别就在于：函数的自变量是数；而泛函的自变量则是函数，泛函的定义域由具有一定条件的一组函数组成。泛函是一个函数集到一个数集的映射；普通函数则是一个数集到另一个数集的映射。泛函的表达式：J=J(y)=Jy(x) J=J(T)=JT(x，y)泛函的一般式：从物理意义上讲，暂时你也可以把泛函理解成熵，自由能等。对于泛函的具体数值我们并不是特别关心，而更关注它何时取得极值，即取什么样的自变量函数，泛函有极值。二、变分的概念普通函数 y(x)泛函

5、Jy(x)自变量的增量自变函数的增量其中：是一个有限小的量， 自变量的微分自变函数的变分 函数的增量泛函的增量函数的微分泛函的变分函数微分运算规则：有函数u(x)和v(x)，则：有函数y(x)，则 泛函变分运算规则：有函数u(x)和v(x)，则：有函数y(x)，则 函数的极值条件泛函的极值条件 三、泛函的极值条件（欧拉方程）泛函的极值条件等价于欧拉方程给出了泛函极值条件与微分方程的关系！利用欧拉方程解教材4.1.2的例题。（注意这是求解变分问题的解析法）在变分问题中，使泛函J(y)为极小值的条件，除外，还应有（二阶变分）四、里兹法（泛函变分的近似解法）（变分原理在求解微分方程中的应用）连续介质

6、问题经常有着不同的但是等价的表达公式微分表达公式和变分表达公式，从例4.1.2中我们可以看到，泛函的变分计算可用微分方程的求解来代替，反之，微分方程的求解也可用泛函的变分计算来代替。求解变分精确解的过程中需要进行各种积分运算，而许多情况下被积函数根本无法找到与相应的初等函数形式的原函数，这说明通过求原函数来计算积分有它的局限性，甚至于可以说这种形式的变分运算根本无法体现出它的运算较微分解方程有什么优越性。变分法的优越性体现在：我们可以找到一种适用于求得以变分形式表达问题的近似解的简便方法，这种方法叫里玆法，是有限元法的前身 。例：用里兹法解微分方程： 边界条件解：构造泛函，在满足边界条件情况下

7、，该泛函的极值条件与微分方程同解。利用欧拉方程可以很容易的证明两者同解（，）令近似函数（或称为试探函数）其中为待定系数，因为近似函数必须满足边界条件，所以我们构造了这样一个函数。将近似函数代入泛函，则此时泛函J已经“变质”了，它不再是函数的函数（泛函）。而实际上是关于未知数的多元函数J()普通多元函数取极值的条件： k=1,2,3,n从这n个代数方程中显然可以求得等n个未知数这里，我们令 0得：，所以：检验：近似解：由解析法得精确解：述评：i）一般而言近似函数（试探函数）的项数愈多，达到的精度愈高。ii）这个方法只是从一族假定中给我们最好的解，因此，非常明显，近似解的精确度和试探函数的选择有关

8、。iii）我们要求试探函数应定义在整个求解区域上，而且它们至少要满足一些边界条件，通常是要满足全部边界条件。iv）通常试探函数是由次数连续增大的多项式构成，但在有些情况下，采用其它类型的函数可能是有好处的。第4.3节 有限单元法一、里兹法的不足（有限元法与里兹法的异同点）相同点：两者实质是相同，数学基础都是泛函变分，求解方法都是以初等函数（多项式）去近似未知函数，利用普通多元函数的极值条件来求解多项式中的待定系数。不同点（里兹法的缺点）：）里兹法将近似函数定义在整个定义域上，而构造近似函数时，又要求它满足所有的边界条件，所以说定义域的边界只能是简单的多边形或多面体（它只能应用在形状相当简单的求

9、解区域上），而不能是复杂形状的边界，故适用范围有限；有限元素法将近似函数定义在一个单元内（给出未知函数的分片近似函数），虽然对单元同样存在几何上的限制，但是由于形状简单的单元可以被集合起来表示非常复杂的几何形状，因此，有限单元法比里兹法的用途要广泛得多。有限单元法的近似函数的自变量（即待定系数）是节点上的场变量，应变量当然是单元内区域上各点的场变量。）对里兹法，一般来说，试探函数的项数越多精度越高；对有限元法，一般可选用线性（一次）多项式试探函数，通过缩小单元格尺寸提高计算精度。）一般情况下，里玆法要求解出试探函数表达式；而有限元法则一般只需求得离散节点上的未知函数值。二、有限元法的解题思路（

10、即一般步骤）1有限元方法的一般步骤（摘录自其它教材）：（1） 连续介质的离散化把连续介质划分成很多元素（单元）。即使在一个求解区域中亦可应用不同形状的元素。（2） 指定节点选择插值函数场变量可能是一个标量，或是一个更高阶的张量。通常选择多项式作为场变量的插值函数，因为它们易于积分和微分，多项式次数的选择取决于每个元素上指定的节点数目，每个节点上未知数的性质和数目以及加在节点上和元素边界上的某些连续性的要求。节点上场变量的大小以及它的导数的大小很可能是未知的。（3） 求出元素特性应用直接法，变分法，加权余数法和能量平衡法等四种方法之一，建立表示各个元素特性的矩阵方程。（4） 集合元素特性以求得系

11、统方程组（总体合成）将表示元素性态的矩阵方程组加以合并，形成表示整个求解区域或系统性态的矩阵方程组系统的矩阵方程组和一个单独元素的方程组具有相同的形式，只是系统方程组包括更多的项，因为它们包括所有的节点。总体集合成的基础是：对于每个共有某个互连节点的元素而言，在这个节点上的场变量值是相同的在准备求解系统方程组之前，必须引入边界条件对系统方程组加以修正（5） 求解系统方程组求解系统方程组，可求得场变量的未知节点值。（6） 按照需求要进行附加计算。利用系统方程组的解来计算其它的重要参数。有限元素法的应用范围： 平衡问题（或不依赖于时间的稳态问题）。 定常状态的特征值问题，如固体或流体振动的固有频率

12、。 传播问题（或非稳态问题）。将时间量纲加入前两类问题而产生的新问题。2有限元法的解题思路（我的总结） 构造泛函JDTD（x,y）= 单元体泛函JeTe(x,y)与整个定义域泛函形式相同,使泛函的极值条件与需要求解的微分方程组等价（如傅立叶导热微分方程和边界条件构成的微分方程组）。 网格剖分：如果将求解区域D（以二维平面为例）划分为E个单元（如三角单元）和n个节点。（注意单元和节点编号规则） 单元分析：在每个单元内部应用温度近似函数（是待定系数），通过解方程可以将待定系数转换为节点温度Tie,Tje ,Tme（待定系数的转换）。这样处理以后，温度近似函数写作：Te=NieTie+NjeTje+

13、NmeTme, 一般称之为温度插值函数，其中Nie,Nje,Nme称为形函数。将温度插值函数Te=NieTie+NjeTje+NmeTme代入单元体泛函JeTe(x,y)，单元体泛函JeTe(x,y)实际上变质为普通多元函数Je(Tei,Tje,Tme)。即：由关于未知函数Te(x,y)的泛函JeTe(x,y)转化成关于待定系数的普通多元函数Je(Tei,Tje,Tme) 然后求解普通多元函数Je(Tei,Tje,Tme)对节点温度Tie,Tje ,Tme的偏导数,。（因为里兹法只有一个单元，而有限元素法包括很多单元，所以它们均不等于零，只有总体合成后才等于零。） 对全部E个单元都进行类似的运

14、算 总体合成：为使泛函JeTe(x,y)取得极值，要求其中k=1,2,3,n上式包含n个代数方程，解之可得n个节点的温度。三、温度场泛函（与傅里叶导热微分方程等价的泛函表达式）一般来说，泛函可根据物理学原理（如熵的最大化，能量最低原量等）或者由微分方程推导而来，具体的推导过程已超出了课程的范围。对于一个二维稳态温度场，函数族所构成的泛函:式中D表示温度场所对应的积分区域。对于一个二维非稳态温度场，函数族所构成的泛函，理论上:但实际上，目前非稳态温度场的处理并不采用上式，而是先把t（或）当作常数（即先把非稳态温度场作为稳态温度场处理），然后将用一价差商代替。泛函J的形式与原微分方程及其边界条件密

15、切相关，原微分方程和边界条件决定泛函J的形式。（原微分方程决定了泛函在区域内部的积分式，边界条件决定了泛函在区域边界上的积分式）1 一般形式傅立叶方程所描述的温度场的泛函微分方程的一般形式其中：表示单位体积内热源的发热强度；kx，ky，kz表示不同坐标方向上的导热系数，且导热系数随温度不同亦即随坐标位置（x,y,z）不同而变化。边界条件的一般形式可表示为其中：h为综合换热系数； T是环境温度；q为单位边界表面上的热流强度； 分别为边界曲面S外法线的方向余弦。 与上述二式对应的泛函 2 无内热源三维稳态温度场,但时，泛函由于q=0,故 3 无内热源三维稳态温度场的泛函，且 (并且是一常数) 1)

16、 若同时具备第三类边界条件，即（对流和辐射换热边界条件），为边界外法线方向上的温度梯度 其中:是常数2）若具备第二类边界条件（边界上热流已知）,则此时 3）若具备第一类边界条件（表面温度已知） 则, 因为在这种边界条件下,边界上温度是固定的，边界上的积分为常数，其变分为零，因而在泛函J中没有与有关的项。4无内热源二维稳态温度场，且 (并且是一常数) 1）若具备第三类边界条件：（改）教材p117写作 两者差一常数,完全等价。 2）若具备第二类边界条件： 3）若具备第一类边界条件： 四、网格剖分网络剖分亦称为时空离散化，就是将时间和空间分割成若干有限的小单元。（一）时间离散化对于稳态导热问题，只涉

17、及空间的剖分。而对于非稳态导热问题，目前对这类方程的泛函变分问题尚未很好解决，通常的处理方法是在空间域内用有限单元网格划分，而在时间域内则用有限差分网格划分（即：有限元法处理空间变量，有限差分法处理时间变量）。具体而言，有两种处理方式（两种处理方法得到的结果完全相同）：(1) 令时间变量暂时固定，即先考虑在某一瞬间对泛函变分，然后再考虑t的变化，把用差分展开为（亦可采用其它差分格式）。(2) 先把用差分展开为（亦可采用其它差分格式），然后进行泛函变分运算，而在变分运算过程中，把与均作为常量处理（，即作为变量）。（二）空间离散化有限元法（Finite Element Method ，简称FEM）

18、，它的最大特点是单元形状和疏密程度可以任意变化，因而对具有复杂几何形状和条件的物体极为适用。1网格剖分的规则从理论上说，有限元法可以采用多种形状的多边形或多面体网格单元，但是通常以三角形和四边形单元用得最多。平面温度场计算一般取三角形网络，三角形的形状和疏密程度可任意变化，对复杂形状的物体尤其能显示其优点。网格剖分的规则：（结合课本P115图4.7讲解规则的应用）（1）一个单元中只能包含一种材料。（2）不能把一个三角形的顶点取在另一个三角形的边上。（3）把求解区域划分为内部单元和边界单元，规定边界单元只能有一条边落在边界上（注意：对曲线边界，实际是两个节点落在边界上）。求解区域的边界为曲线时，

19、剖分时用直线代替，并取为三角形单元的一边。（4）三角形的三条边的长度不宜相过大，图为计算精度受单元最长边长与最短边长之比值的控制。2单元和节点编号规则为了减少失误，提高计算效率，便于计算机实现。单元和节点编号有一定规则：（1）单元编号规则（结合课本P115图4.7讲解规则的应用）单元分为两大类：内部单元，边界单元单元编号规则：先编内部单元后编边界单元，一、二、三类边界单元依次编排。（原因是单元体泛函由依序由简单到复杂，有助于简化总体合成后矩阵方程中的系数矩阵，从而提高算法效率。其原因仍然有待进一步深入思考） 记作：（2）节点编号规则 每一个节点分别有两个序号：局部序号，全局序号节点全局序号编号

20、规则：（结合课本P136图4.11讲解规则的应用）依序编排，力求使邻近节点的编号尽可能接近，尤其是同一单元中三个节点的编号相差不宜太大。（原因：合理地编排节点全局序号，可以使得总体合成后矩阵方程中的系数矩阵有规则地分布在主对角线附近相对狭小的宽度内，构成所谓带状矩阵。这种对称正定带状矩阵对利用消去法求解极为有利，有助于大大提高算法效率。）记作：1，2，3，4节点局部序号编号规则：()内部单元：i，j，m按逆时针方向依次编排，起始位置可以是任意的；()边界单元：对于三角形单元，节点i与边界相对（即内节点），i，j，m逆时针方向依次编排；对于其他多边形单元，保证j，m落在边界上。记作：i，j，m，

21、五、单元分析1. 温度插值函数温度插值函数实际上就是里兹法试探函数（近似函数）的变形。插值函数一般由不同次幂的多项式组成，但也有采用三角函数或指数函数组成的乘积表示，但最常用的多项式插值函数。有限元插值函数分为两大类，一类只要求插值多项式本身在插值点取已知值，称为拉格朗日(Lagrange)多项式插值；另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标，有对称和不对称等。常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状，一维看作长度比，二维看作面积比，三维看作体积比。在二维有限元中，

22、三角形单元应用的最早，近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元，常采用的插值函数为有Lagrange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。1.1 单元内部的温度插值函数设三角形单元e上的温度是x，y的线性函数，即 其中：为待定系数将节点坐标及温度代入得 写成矩阵形式： 利用矩阵求逆的方法，可以解出待定系数 式中：为三角形单元的面积，其值为 （为了使面积不为负值，故要求节点的局部序号i,j,m必须按逆时针方向排列） 将式代入式，得 上式通常可以简写成式中 称为形状因子，形状函数或简称形函数。1.2 单元边

23、界上的温度插值函数4式适用于整个单元区域，对于单元边界当然也是适用的，但是，根据线性插值的概念，既然单元边界jm上两节点的温度分别为和，那么直线jm上任意一点的温度应在和之间呈线性变化，而与无关，这样在边界上就可以构造一个更加简单的插值函数式中g为一参变量，；g=0对应于节点j，g=1对应于节点m。1.2 边界弧长（积分变量）s与参变量g的关系显然jm的边长是一个已知数而曲线积分中的边界弧长变量s与间的关系可利用g联系起来 2. 单元变分运算（即普通多元函数偏导数计算）2.1 二维稳态温度场单元变分计算（实际上可理解为偏导计算，因为用温度插值函数代泛函以后，泛函实际上已经成为一个普通多元函数）

24、2.1.1第三类边界单元变分计算 （我不成熟的理解：泛函表达式的第二部分原来也是对整个单元的面积分，但是被积函数只有在边界上才有非零值，所以转换成了边界的线积分）泛函表达式的第二部分,对整个求解区域而言，边界一般是封闭的，记作“”；然而对某个单元而言，边界一般是不封闭的，记作“”。将三角形单元内部的插值公式以及边界上的插值公式代入泛函，泛函实际上变成了普通多元函数求，可以采用两种方法：(1) 先将插值公式代入式，然后再求， 单元内部插值函数： 单元边界插值函数： 积分变量的转换： 接下去，重点看下列推导过程的后半部分。 同理： 写成矩阵形式 其中(2) 先推导变分公式，再将插值函数代入令，则同

25、理可求得，（后续推导略）2.1.2 内部单元变分计算内部单元的泛函 将三角形单元内部的插值公式以及边界上的插值公式代入泛函，泛函实际上变成了普通多元函数求，可以采用两种方法：(1) 先将插值公式代入式，然后再求， 单元内部插值函数： 同理 同理： 写成矩阵形式其中从形式上说，内部单元的变分结果与边界单元基本相同；从本质上说，仅需令边界单元变分结果中的，即可得内部单元的变分结果。(2) 先推导变分公式，再将插值函数代入令，则同理可求得，（后续推导略）2.2 二维非稳态温度场单元变分计算，无内热源二维非稳态温度场泛函变分问题尚未很好解决，目前采用的处理方法有两种：先用有限元法处理空间域，后用有限差分法处理时间域先用有限差分法处理时间域，后用有限元法处理空间域我们采用前一种处理方法。二维非稳态温度场泛函（先固定时间变量，即先把看作常数）：2.2.1第三类边界单元变分计算无内热源二维非稳态温度场第三类边界单元泛函无内热源二维稳态温度场第三类边界单元泛函两者相比较，仅有项是新的，并且我们在整个变分运算（实际上是普通函数积分运算）过程中，把作为常数处理。 因此可以证明 得到：同理可得： 写成矩阵形式 注意：此处的并非插值函数中的形状函