有限元法的基础理论

一、里兹法与迦辽金法（摘自电磁场有限元方法 金建铭） 1. 里兹法

里兹法是一种变分方法，其中边值问题用变分表达式（也称泛函）表示，泛函的极小值对应于给定边界条件下的控制微分方程。通过求泛函相对于其变量的极小值可得到近似解。 2. 伽辽金法

伽辽金法属于残数加权方法类型，它通过对微分方程的残数求加权的方法得到方程的解。

若u

是方程的近似解，将u 代入方程可得到非零的残数： r Lu

f =- u

的最佳近似应能使残数r 在Ω内所有点上有最小值。残数加权方法要求： 0i i R rd ωΩ

=Ω=?

这里i R 表示残数的加权积分，i ω是所选的加权函数。

在伽辽金法中，加权函数与近似解展开中所用的函数相同。通常，这样可得到最精确的

解。

二、有限元方法

里兹法和伽辽金法中，在整个解域内找出能表示或至少近似表示问题真实解的试探函数是非常重要的。然而对于许多问题，这个步骤是十分困难的，对二维和三维问题尤其如此。为此，我们可将整个区域划分成小子域，并应用定义在每个子域上的试探函数。因为子域是小区域，因而在每一子域内函数的变化不大，所以定义在子域上的试探函数通常比较简单。这正是有限元法的基本思想。应用里兹法的过程通常称为里兹有限元法或变分有限元法，而应用伽辽金方法的过程通常称为伽辽金有限元方法。

有限元法与经典里兹法和伽辽金法的不同之处是在试探函数的公式上。在经典里兹法和伽辽金法中，试探函数由定义在全域上的一组基函数组成。这种组合必须能够（至少近似）表示真实解，也必须满足适当的边界条件。在有限元法中，试探函数是由定义在组成全域的子域上的一组基函数构成。因为子域很小，所以定义在子域上的基函数能够十分简单。 三、关于形函数（摘自有限元法在电磁计算中的应用 张榴晨）

对于一个待求的微分方程，用一组线性独立的尝试函数i ψ和待定系数i C 来表示方程的近似解，并用加权余数法（迦辽金法）来求解这些待定系数。求解待定系数的代数方程组为：

1

[]1,2,,n

i j i j i d C q d j n ψψψΩ

Ω

=??Ω=Ω

=∑?

?

这里j ψ为所选择的加权函数，应用迦辽金法时，所选取的加权函数即为尝试函数。 有限元中应用的尝试函数代表了单元上近似解的一种插值关系，它决定了近似解在单元上的形状。因此尝试函数在有限元法中又称为形函数。对于一维有限元来说，形函数为一个直线段；对一维高阶有限元来说，形函数为一个曲线段；对二维一阶有限元来说，形函数为一个平面；对二维高阶有限元来说，形函数为一个曲面；三维有限元来说，形函数为多维平面或曲面。选择形函数时可以使一个任意元上的函数只与该元所对应的节点势函数值有关，而与其它各点的值无关。 1. 一维有限元

对于一维有限元来说，形函数分段线性。该形函数i ψ在节点i 上的值为1，并在与节点

i 相邻的两个单元上线性减小，直到在相邻节点1i -和1i +上分别减小为零。其形函数的形

式为：

i i i x ψαβ=+

根据形函数的性质即可得到两个线性方程组，解这两个线性方程组即得到系数i α和i β，进而得到i ψ。应用同样的方法可以得到每一个节点上的形函数表达式。 2. 二维有限元

以三角形单元为例，单元的顶点分别为,,i j k 。每个顶点都对应于一个形函数。考虑一阶有限元法，即用线性插值的方法表示一个有限元上的势函数分布。这样，三角形作为一个单元所对应的三个形函数都是分段线性的函数，并可以看作由几个平面组成。对应于节点i 的形函数i ψ在节点i 上取单位值（即为1），并由此在单元e 上直线下降，直到在其它另外两个节点j 和k 上降为零，在单元之外的区域则一直保持零值。于是这三个节点值就决定了该形函数的形状。形函数i ψ可由一个线性函数表达为：

i i i i x y ψαβγ=++

四、势函数分布

有限元的作用在于求解分布场的势函数在每个节点上的近似值，而势函数在单元的其它位置的值，可用插值的原理来表示。如果采用线性插值的方法来表示分布势函数，则称为一阶有限元法，如果有限元法采用高阶插值法表示分布势函数，则称为高阶有限元法。

任意一个单元上的势函数分布由这个单元上的节点势函数值及相应的形函数表示，对于一个一维单元有：

1

1e i i i i ??ψ?ψ++=+ 整个区域的势函数分布则由每一个单元的势函数分布相加得到。 五、关于对称性的利用

利用对称性可以减少节点和单元的数目，从而节省用于建立模型和计算近似解的时间。从另外一个角度来讲，如果维持节点和单元的数目不变，则利用对称性可以对这四分之一区域作更为详细地划分，即单元的尺寸可大为减小，从而提高近似解的精度。 case 1（a ）：一个两极电容器的静电场分布问题。假设两个极板都接在0.5v 电源端，极板间的距离为2，极板间充有密度为ρ的自由电荷。这样，电容器的激励和几何形状都关于y 轴对称。由于这种关系的存在，我们只要求解电场在整个区域的一半的分布即可，而另一半的电场分布则能从对称关系而得到。从边界条件的观点出发，这种对称的结构导致电力线垂直穿过y 轴，使电势在该对称轴上沿x 方向的变化率为零。

描述这一问题的微分方程和边界条件为：

21

(0,1)

0.50

x x x n

?ρ

?

?==?=-∈=?=?

见原文page62 case 1（b ）：当考虑极板端部的边缘效应时（考虑一个方向的边缘效应），这个问题则为一个二维边值问题。假设电势分布在与极板垂直平等的截面上分布相同，那么只需计算一个截面的电势分布即可。在y 轴两侧，几何结构及施加电压都对称，因此y 轴为对称轴，即沿y 轴的电势对法线方向（x 方向）的变化率为零（是否可以理解为等势面与y 轴垂直而与x 轴平行），这样y 轴便构成一个齐次诺伊曼边界条件。从另外一个角度看，该问题的几何结构x 轴对称，而施加电压则x 轴反对称，因此沿x 轴的电势应取两极板施加电压的中值，即零值。因此x 轴构成一个狄利克莱边界条件。于是计算区域被减为原区域的四分之一。

为了把边缘效应考虑进去，计算区域还应考虑电容器周围空间，并应将无限远处设为零电位参考点。但从实际意义来讲，假设电容器外某一定距离的空间处为零电位参考点即可满足实际需要。设沿x 轴方向的2m 处电位为零。

描述这一问题的微分方程和边界条件为：

12

3

4

20

1000

n

???

?

?ΓΓΓΓ?====?=?

case 2：一个同轴传输线，两个同芯长方形导体间充满了线性介质，其介电常数为ε。假设两导体间加有直流电压10V ，导体间贮有密度为ρ的自由电荷，传输线的长度远远大于其

截面的长和宽，那么可以认为电场在传输线各个截面上的分布都相同，因此只需求解电场在某个截面的分布，从而检查绝缘材料的工况。该问题的几何形状、介质及激励都x 轴和y 轴对称，因此只需求解整个截面的四分之一即可。在对称轴的两侧，电势对于该轴线的法向

变化率为零。从电力线的观点出发，也可以说电力线垂直对称轴线。

描述这一问题的微分方程和边界条件为：

13

2

4

2,()010

q

x y q n

n

ρ?ε

??

??ΓΓΓΓ?=-∈Ω=

==??=

=??

见原文page68。

case 3：一个简单的变压器的静电场分布。变压器结构及激励关于x 轴对称，因此沿x 轴有齐次的诺依曼条件成立；变压器的几何结构y 轴对称，且激励y 轴反对称，因此沿y 轴存在着磁势为零的狄利克莱条件。这样，只需选取变压器的四分之一就能完全解出磁场分布。

这一静磁场的描述方程为帕松方程及拉普拉斯方程，由于只考虑二维结构，磁势只在沿z 轴的方向不为零，而沿其他方向的分量均为零。对于这样一个二维问题的矢量磁势A 实际上被简化为一个标量（或为沿一个方向的矢量）z A ：

201

00

z r

A J

n

μμ??

?=-=?=?

四、高阶有限元

一阶有限元的形函数是一个线性函数，应用一阶有限元形成的系数矩阵为稀疏矩阵，其求解简单但精度较低。解决精度问题的方法之一为采用高阶有限元法。

对于一维问题用曲线来逼近，即采用高阶插值的方法通常可以得到更为精确的结果。例如采用二阶插值：

2

(1,2,3)i i i i a b x c x i ψ=++=

二阶有限元由三个点构成，其中的两个点为构成该有限元区域的端点，而第三个点通常选在这两点中间，这三个点分别对应三个二阶形函数，也就是说在一个单元中，每个节点都对应一个二阶形函数。

一般来说高阶有限元允许我们我用尺寸较大的单元来描述具体工程问题，特别是几何形状较复杂的问题。从另外一个角度来说，如果保持单元的尺寸不变，相对于一阶有限元来说，高阶有限元能够提高近似解的精确度。但高阶有限元也增加了计算量，因此更耗时。

Page106给出了一个结构分析中所用的几种三维有限元法的比较，其结果很能说明有限元形状和疏密对计算精度的影响。

五、单元与插值函数（有限单元法基本原理和数值方法 王勖成） 一个函数在域内其本身连续，它的一阶导数具有有限个不连续点但在域内可积，这样的函数称之为具有0C 连续性的函数。类似地，如果微分算子A 出现的最高阶导数是n 阶，则要求函数u 必须有连续的1n -阶导数，即函数应具有1n C -阶连续性。

关于单元插值函数的形式，有限单元法中几乎全部采用不同阶次幂函数的多项式。这是因为它们具有便于运算和易于滑坡路收敛性要求的优点。如果采用幂函数多项式作为单元的插值函数，对于只满足0C 连续性的单元（称0C 型单元），单元内的未知场函数的线性变化能够仅用角（或端）结点的参数表示。对于它的二次变化则必须在角（或端）结点之间的边界上适当配置一个边内结点（二次单元）。它的三次变化则必须在每个边界上配置二个边内结点（三次单元）。配置边内结点的另一原因是常常要求单元的边界是曲线的，沿边界配置适当的边内结点可以构成二次或更高次多项式来描述它们。 1. 一维单元

对于n 个结点的一维单元，i N 可以采用1n -次Lagrange 插值多项式(1)()n i l x -，即令：

(1)

12111,1211()()()()()()()()()()()()(1,2,,)

n

j n i i n i i

j j i

i j

i i i i i i i n x x x x x x x x x x x x N x l x x x x x x x x x x x x x i n --+=≠-+------==∏

=

------= 构造插值函数时，一般用无量纲表示，一维单元中的无量纲为长度坐标，更一般化的可称为自然坐标。

为构造形式的Lagrange 单元方便，还可把插值多项式写为：

(1)

1,()()()()

n

j n i i

j j i

j i f N x l x f ξξ-=≠==∏

其中()j j f ξξξ=-表示任一点ξ至点j ξ的距离，也是j 点坐标j ξξ=表示成方程形式

()0j j f ξξξ=-=的左端项。显然，()0j j f ξ=。(1)n i i N l -=的展开式中包含了除()i f ξ外

的所有()(1,2,,1,1,,)j f j i i n ξ=-+ 的因子，从而保证了()0()i j N j i ξ=≠这一要求。()j i i j f ξξξ=-是点I 的坐标代入()j f ξ的结果，

这一因子引入(1)n i l -的分母是为了保证满足()1i i N ξ=的要求。

2. 二维单元

1) 三角形单元

I. 一次单元

对于3结点三角形单元，引入面积坐标

/(1,2,3)i i L A A

i ==

则单元插值函数可表示为：

(1,2,3)i i

N L i ==

II. 二次单元

二次单元有六个结点，各结点的面积坐标分别标注在图中。(page98)需要构造的插值函数表示成： ()2

1,2,3()1

1,2,3()()()

i j i i j j i i i f L L L N x f L L L ==∏

其中的()

1,2,3()i j f L L L ，()

1,2,3()i j i i i f L L L 赋与了和三角形单元相对应的几何意义。 ()

1,2,3()i j f L L L 是通过除结点i 以外所有结点的二根直线方程()

1,2,3()0i j f L L L =的左端项。如当1i =时，(1)

j

f 分别是通过结点4,6的直线方程(1)

11,2,31()1/20f L L L L =-=和通过

结点2,5,3的直线方程(1)

21,2,31()0f L L L L ==的左端项。()

1,2,3()i j i i i f L L L 中的1,2,3i i i L L L 是结点i 的面积坐标。所以得到：

11

1111/2(21)1/21

L L N L L -=

=-

类似的方法可以得到其它的形函数。

III. 三次单元

为保证二维域三次多项式的完备性，三次单元应有10个结点，按与二次单元相同的步骤按画线法构造它的插值函数。Page98

2) 矩形单元

如果所研究问题的总体域是矩形的，采用矩形单元将比三角形单元更有效。

I. Lagrange 矩形单元

构造任意的Lagrange 矩形单元插值函数的一个简便而系统的方法是利用二个坐标方向适当方次Lagrange 多项式的乘积。在第J 行和第I 列的结点上的插值函数是：

()()()()r p i IJ I J N N l l ξη==

虽然构造Lagrange 插值函数很容易，但是这一类型的单元存在一定缺点，主要是出现了随插值函数方次增高而增加的内结点，从而增加了单元的自由度数，而这些自由度的增加通常并不能提高单元的精度。

II. Serendipity 单元

通常我们希望结点仅配置在单元的边界上，并且在实际应用中常希望同一单元的不同边界有不同数目的结点，这样可以实现不同阶次单元之间的过渡，从而可能在求解的不同区域采用不同精度的单元。假定开始只有四处角结点，对应这些结点的插值函数可利用双一次Lagrange 多项式构造，即：

001

(1)(1)(1,2,3)4

i N i ξη=

++=

其中，00,i i ξξξηηη==

如果增加边内结点，则与它对应的插值函数可以按划线法构造或直接表示成ξ（或η）方向二次和η（或ξ）方向一次Lagrange 多项式的乘积。例如增加结点5（位于结点1,2之间），则：

251

(1)(1)2

N ξη=

--

需要指出的是，对于现在5个结点的情况，5N 满足55(1,2,,5)j j N j δ== 要求，而

(1,2,3,4)i N i = 不再满足55(1,2)i i N j δ==

的要求了。

（是否理解为插值函数在该结点上为1而在其它结点上为零，在1,2结点间增加了结点5后若仍然采用线性插值函数则对于1N 将不满足插值函数的特性，因为1N 在结点5处不为零，为了保证在结点5处1N 为零则必须使插值函数曲线在该点过零。如图page103图3.15所示，这时的插值函数是一条曲线。）为满足要求，1N ，2N 需要修正为：

11512N N N =- ，22512

N N N =-

类似地，可以讨论增加边内结点6，7，8的情况，

11582256

33674478

2256227811

11,

222211

11,

22221

1

(1)(1),

(1)(1)

221

1

(1)(1),

(1)(1)

22

N N N N N N N N N N N N N N N N N N N N ξηηξξηηξ=--=--=--=--=--=-+=-+=-- 若5，6，7，8结点中任一个不存在，则对应的插值函数为0。

对于高次的Serendipity 单元，可用同样的方法构造它的插值函数。例如对于p 次单元的边内结点，它的插值函数可表示成ξ（或η）方向p 次和和η（或ξ）方向一次Lagrange 多项式的乘积，而对于角结点则其插值函数可表示成一双线性函数和用适当分数分别乘以相邻两个边界上的各个边界结点插值函数之和，以保证它在边内结点上的值也为0。因为角结点插值函数中双线性函数之附加项可随相邻边界上的边内结点的增减而变化，所以边内结点数是灵活的，这对于构造变结点数的过渡单元是很适合的。 3. 三维单元

1) 四面体单元

四面体单元与二维情况的三角形单元类似，插值函数是在三维坐标内的各次完全多项式。在各个面上的结点配置与同次的二维三角形单元相同，函数是相应二维的完全多项式。

根据三维四面体单元的几何特点，引进的自然坐标是体积坐标。 2) Serendipity 单元 3) Lagrange 单元

以上内容详见原文page106。 六、等参单元

对于一个给定问题的求解域，预期用较少的单元即可获得需要精度的解答，但是用较少的形状规则的单元离散几何形状比较复杂的求解域常会遇到困难，因此需要寻找适当的方法将规则形状的单元转化为其边界为曲线或曲面的相应单元。在有限元法中最普遍采用的变换方法是等参变换，即单元几何形状的变换和单元的场函数采用相同数目的结点参数及相同的插值函数进行变换。采用等参变换的单元称为等参单元。 1. 等参变换

为将局部（自然）坐标中几何形状规则的单元转换成总体（笛卡尔）坐标中几何形状扭曲的单元以满足对一般形状求解域进行离散化的需要，要建立一个坐标变换。

最方便的方法是将这种坐标变换也表示成插值函数的形式：

1

1

1

m

m

m

i i

i i

i i

i i i x N x y N y z N z ==='''===∑∑∑

其中m 是用以进行坐标变换的单元结点数，,,i i i x y z 是这些结点在总体（笛卡尔）坐标内的坐标值，N '为形函数，实际上它也是用局部（自然）坐标表示的插值函数。通过上式

建立起两个坐标系间的变换，从而将自然坐标内的形状规则的单元（母单元）变换为笛卡尔坐标内的形状扭曲的单元（子单元）。

坐标变换关系式和函数的插值表示式在形式上是相同的。如果坐标变换和函数插值采用相同的结点并采用相同的插值函数则称这种变换为等参变换。 2. 导数间的变换

函数i N 对ξ的偏导数可表示成为：

i i i i N N N N x y z

x y z ξξξξ

???????=++

??????? 对于其它两个坐标(,)η?，可写出类似的表达式，将它们集合成矩阵形式则有：

i i i i i i i i i N x y

z N N x x N N N x y z J y y N x y z N N z z ξξξξηηηη?

?

??????????????

????????????????????????

??????????????

==??????????????????????

??????

??????????????????????????

???? 式中J 为Jacobi 矩阵，可记作(,,)

(,,)

x y z ξη???，根据上述的插值关系，J 可以显式地表示为

自然坐标的函数：

'

'

'''

'

12'

'''''1

2'''1(,,)(,,)m m

m

i i i m i i i i i i i i i m m

m

i i i m i

i i i i i i

i i m m

m

i i i i i i i i i i i i N N N N N N x y z N N N N N N x y z J x y z N N N N x y z ξ

ξ

ξξ

ξξξη?η

ηηηηη?

?

?=========?????????????????????????????≡==??????????????????????????∑∑∑∑∑∑∑∑∑ 1

11222''

'2m

m

m m x y z x y z x y z N N ?

?

?????

???????????

????

?????????????

??

???????

,,d d d ξη?在笛卡儿坐标系内形成的体积微元是：

()dV d d d ξη?=??

而：

x y z

d d i d j d k x y z d d i d j d k x y z d d i d j d k

ξξξξξξξηηηηηηη??????????=

++??????=++??????=++???

其中，,,i j k 是笛卡儿坐标,,x y z 方向的单位向量。将,,d d d ξη?代入dV 可得：

x y z x y z

dV d d d J d d d x y z ξξξξη?ξη?ηηη???

?????????=

=????????? 3. 等参变换的条件

两个坐标间一对一变换的条件是Jacobi 行列式J 不为0，等参变换作为一种坐标变换也必须服从这个条件。

如果0J =则表明笛卡尔坐标中体积微元（或面积微元）为0，即在自然坐标中的体积微元d d d ξη?（或面积微元d d ξη）对应笛卡尔坐标中的一个点，这种变换显然不是一一对应的。

笛卡尔坐标中的面积微元可直接表示成：

sin(,)dA d d d d d d ξηξηξη=?=

而： dA J d d ξη=

于是有：sin(,)

d d d d J d d ξηξηξη

=

可见，只在以下三种情况之一成立，即

0,0,sin(,)0d d d d ξηξη===或或就会出现0J =。因此，为了保证变换的一一

对应，应防止因任意的二个结点退化为一个结点而导致,,d d d ξη?中的任意一个为0，还应防止因单元过分歪曲而导致的,,d d d ξη?中的任何二个发生共线的情况。