第一章 有限元概貌与发展
有限元方法是近似求解数理边值问题的一种数值技术。这种方法大约有60年的历史。它首先在本世纪40年代被提出，在50年代开始用于飞机设计。后来，该方法得到了发展并被非常广泛地用于结构分析问题中。目前，作为广泛应用于工程和数学问题的一种通用方法，有限元已相当著名。
有限元法应用于电磁场中，最先是用结点上的插值基函数来表征该结点上的矢量电场或磁场分量的，称为结点有限元。但是，在使用结点有限元进行电磁仿真时，会有几个严重的问题。首先，非物理的或所谓伪解可能会出现。其次，在材料界面和导体表面强加边界条件很不方便。再次，处理导体和介质边缘及角也很困难，这是由与这些结构相关的场的奇异性造成的。在这些问题中，最后一个问题比其它两个问题更严重，因为它缺少通用的处理方法。即使对前两个问题，目前的处理状况也不能完全令人满意。因此，有必要探讨其它的可能性或其它方法，而不仅仅是改进，从而将电磁场有限元分析引入一个新的时代。
幸运的是，一种崭新的方法已经被发现。这种方法使用所谓矢量基或矢量元，它将自由度（未知量）赋予棱边而不是单元结点。因为这个原因，它也叫棱边元（edge element ）。虽然Whitney 早在35年前就描述过这些类型的单元，但它们在电磁学中的应用及其重要性直到前几年才被认识到。在80年代初，Nedelec 讨论了四面体和矩形块棱边元的构造。Bossavit 和Verite 将四面体棱边元应用于三维涡流问题。Hano 独立地导出了矩形棱边元，并用于介质加载波导的分析。Mur 和de Hoop 考虑了非均匀媒质中的电磁场问题。Van Welij 和Kameari 应用六面体棱边元进一步考虑了棱边元在涡流计算中的应用。Barton 和Cendes 将四面体棱边元应用于三维磁场计算，同时，Crowley 提出了一种更复杂的单元类型，即所谓的协变（covariant ）投影单元，它允许单元带有弯曲的棱边。在所有这些工作中，已经证明：棱边元没有前面提到的所有缺点。因为这些缺点困惑了研究者多年，可以想象，棱边元的重要性很快就被认识到了，因此，在过去的几年中，开展过大量的研究也获得了很多成功的应用。
棱边元分析时谐电磁场问题通常是基于矢量波动方程的，但这也常常引发一个问题:作为电场强加条件的[]0E ε?=i ，低频时在波动方程中不起作用。这一点使得所产生的有限元矩阵严重病态。尽管这一缺陷并不影响基于场形式波动方程的边棱元方法在高频电磁分析中的应用，但是在使用自适应网格剖分产生的一些小单元中，会局部的逼近静态的情况，同样会产生病态的矩阵的。这种病态的矩阵不利于直接法的求解，并且也大大的降低了迭代法的收敛速度。近年来发展了一种矢量标量位有限元，将边棱元与结点元相结合，能在很宽的频带内直接强加电场散度条件，同时有效的改善有限元矩阵的性态，尤其适合于结合迭代法求解。
为了降低有限元矩阵未知量的数目，不少学者对高阶有限元也作过了大量的研究工作。它的主要思想就是利用高阶的基函数对未知的场获得更精确的逼近，或者说在较稀疏的网格上获得一般的线性插值基函数在稠密网格上同样的精度。但这种高阶的有限元通常也会产生严重病态的矩阵，不利于快速求解，一直是个待解决的问题。
第二章 有限元方法入门
在本章中，我们首先回顾求解边值问题的两种经典方法，它们都包含着有限元方法的本质。然后，应用一个筒单的例子，介绍有限元方法。最后，我们不针对任何特定问题来描述该方法的基本步骤。
2．1 边值问题的经典方法
在本节中，我们首先定义边值问题，然后讨论边值问题求解的两种经典方法，一是里兹
(Ritz)变分方法，另一种是伽辽金(Galerkin)方法，它们构成了现代有限元方法的基础。 2．1．1 边值问题
边值问题出现在物理系统的数学模型中，它们的求解一直是数学物理中的研究主题。典型的边值问题可用区域内的控制微分方程和包围区域边界上的边界条件来定义。微分方程可表示
L f Φ= （2.1） 式中，L 是微分算符，f 是激励或强加函数，Φ是未知量。在电磁学中，控制微分方程包括简单的泊松方程以及复杂的标量波动方程，甚至也有更复杂的矢量波动方程。边界条件有简单的狄利克雷(Dirichlet)条件和诺曼(Neumann)条件，也有复杂的阻抗和辐射边界条件，甚至还有更复杂的高阶条件。
当然，我们希望尽可能用解析方法求解边值问题。然而，只有少数情况可得到解析解。在电磁学中，可解析求解的问题包括无限大平行板间的静电势问题，波在矩形、圆柱和椭圆波导中的传输问题，矩形、圆柱和球形腔内的腔体谐振问题，以及无限平板、劈、圆柱和球对波的散射问题等等。许多实际重要的工程问题都没有解析解。为了克服这种困难，人们已发展了各种近似方法，其中应用最广泛的是里兹方法和伽辽金方法。 2．1． 2 里兹方法
里兹方法也称为瑞利一里兹(Rayleigh-Ritz)方法，它是一种变分方法，其中边值问题 用变分表达式(也称为泛函)表示，泛函的极小值对应于给定边界条件下的控制微分方程。通过求泛函相对于其变量的极小值，可得到近似解。为了说明这种过程，首先让我们定义内积，用尖括号表示为
,d ??
ΦΨ=
ΦΨ?∫
(2．2)
式中星号表示复共轭。在这种内积定义下，如果有
,,L L ΦΨ=ΦΨ (2. 3) 则(2．1)式中的算符是自伴的，如果有 L 00
,0
L >Φ≠?ΦΦ=?
=Φ=? (2. 4) 则(2．1)式中的算符是正定的。可以证明：如果(2．1)式中算符既自伴又正定，那么，
L L (2．1)式的解可通过求下式泛函对Φ
的极小值得到： ()111,,222
F L f f Φ
=ΦΦ?Φ?Φ , (2. 5) 式中，表示试探函数。 一旦确定了泛函，即可用下述步骤来求解。为简单起见，我们假设间距是实数值的，并假定(2．5)式中的Φ
Φ
可近似展开为 (2．6) {}{}{}{}1N
T T j j
j c v c v v c =Φ===∑ 式中,叫是定义在全域上的展开函数，是待定的展开系数。
{·}仍表示列向量，上标T j v j c
表示向量的转置。将(2．6)式代入到(2．5)式中，我们得到 {}{}{}
{}{}
{}12
T T
T
F c v L v d c c v fd ??
=
??∫?∫ (2．
7) 为了求极小，我们令其对的偏导数为零，从而得到线性代数方程蛆： ()
F Φ i
c
{}{}{}{}()11122
120
1,2,3,,T
T i i
i N
j i j j i i j F v L v d c c v Lv d v fd c c v Lv v Lv d v fd i N
??
??=i
??=?+???=+???
==∫∫∑∫∫ ?
∫
(2.8)
可将其写成下列矩阵方程
． []{}{}S c b = (2．
9) []S 的元素为
(1
2
ij i j j i S v Lv v Lv ?=
+∫)d ? (2．10) []b 的元素为
i i b v fd ?
=?∫ (2.11)
显然，[]S 是对称矩阵。引用算符的自伴性质，可写成 L ij S ij i
j S v Lv d ?=
?∫
(2．12)
求解矩阵方程(2．9)式，即可得到(2．1)式的近似解。
2．1．3伽辽金方法
伽辽金方法属于残数加权方法类型，正如其名称所指，它通过对微分方程的残数求加权
方法来得到方程的解。假设Φ
是(2. 1)式的近似解，用 Φ 代替(2．1)式中的中，得到非零的残数
0r L f =Φ
?≠ (2．13) Φ
的最佳近似应能使残数r 在内所有点上有最小值。残数加权方法要求 ? 0i i R w rd ?
=
?=∫
(2．14)
这里i R 表示残数加权积分，是所选择的加权函数。
i w 在伽辽金方法中，加权函数与近似解展开中所用的函数相同。通常，这样可得到最精确
的解，因而，这是建立有限元方程的常用途径。为了更清楚地说明这种方法，我们假设方程
的解可用(2．6)式中的方法表示，那么，加权函数选为
1,2,3,,i i
w v i N == （2．15）
因此，(2．14)式变成
({}{})1,2,3,,T i i i R v L v c v f d i N ?
=??
=∫ (2．16)
这里同样得到了(2．9)式给出的矩阵方程组。在此，除非算符是自伴算符，否则，矩
阵 []不一定是对称的。在算符为自伴算符的情况下，伽辽金方法与里兹方法得到相同的方程组。
L S L 顺便指出，除了选择展开函数作为加权外，也能选择其它函数作为加权。这样可得到不同的公式。下面作简要讨论。 1．点配置法
在电磁学界，这种方法也称为点匹配法。狄拉克(Dirac)δ函数被选为加权函数(在点，
；在其它各点，)，结果，(2．14)式变成
i i w =∞0i w =[{}{}]0T i i R L v c f =?在点处＝ （2.17）
显然，这等价于在给定点上满足(2．1)式。通常，选择的匹配点数等于未知量个数。
2．子域配置法
在这种方法中，令加权函数在给定子域内等于1，在其它地方为零。因而有
({}{})0i
T i R v c f d ?=?∫?= (2．18)
式中，表示第i 个子域。通常，所选择的子域数也等于未知量个数。 i ? 3. 最小二乘法 定义新的误差项
21
2I r d ?
=
?∫ (2．19) 最小二乘法就是求上式的极小值。极小值是对近似解的未知系数求的，这相当于
({}{})0T i i
I
Lv L v c f d c ??=??∫?= (2．20) 在这种情况下，加权函数显然为。
i Lv 2．2 —个简单的例子
简单讨论里兹和伽辽金方法后，我们将它们应用于简单的边值问题来说明求解过程， 然后，再用这个例子介绍有限元方法。希望用这种方式介绍有限元方法将更加自然而不易引起混淆。
2．2．1 问题的描述
这里考虑的特定问题是求两无限大平行板间的静电势Φ。一块板位于处，;另一块板位于0x =0V Φ=x =1m 处，Φ＝lV。两块平行板间充满了介电常数为ε(F／m)的媒质，其间的电荷密度是变化的，()(1)x x ρε=?+C／。此问题可用泊松方程(1．16)式作数学描述，简化成二阶微分方程
3
m
22
10d 1x x dx Φ
=+<< (2.21)
其边界条件为
0x =Φ
= (2.22)
1
1x =Φ
= (2.23)
此问题的精确解为
32111
()623
x x x Φ=
++x (2.24) 对(2．21)式积分两次，然后再应用(2.22)式和(2．23)式定出积分常数，就可很容易地
得到上式的解。然而，许多实际问题没有这样简单的解。在许多情况下，只能得到近似解。由于这种原因，我们假设不知道这个精确解，而用里兹方法和伽辽金方法来求解。 2．2．2 用里兹方法求解
如前所述，在里兹方法中，我们用泛函给出问题的公式，泛函的极小值对应于给定边界
条件下的微分方程。对于(2．21)式至(2．23)式中定义的问题，可以证明泛函由下式给出：
211001()(1)2d F dx x dx ??Φdx Φ=++????
∫∫ Φ （2.25） 根据这种方法，我们首先用多项式展开Φ
(注意多项式只是许多试探函数中的一种)： 21234
c c x c x c x Φ=+++ 34c (2.26) 式中， (i 1，2，3，4)是待定的未知常数。应用边界条件(2．22)式和(2．23)式，我们得到。因而，(2.26)式退化为
i c =1230,1c c c ==??2334
()()()x x c x x c x x Φ=+?+? (2.27) 将(2.27)式代入(2．25)式，井进行积分得到
22433443211231562604F c c c c c c =++??+4
3
(2.28)
它对和的导数分别为 3c 4c
343113
2F c c c ?=+??1
4 (2.29)
3441422
56F c c c ?=+??3
0 (2.30)
令上两式分别等于零，我们得到两个线性方程，它们的解给出312c =
和41
6
c =，因而得到了与(2．24)式中相同的解。在这种情形下，因为所用的试探函数具有表达精确解的能力，所以，
我们能得到精确解。在许多情况下，这是不可能的，因而只能得到近似解。
2. 2．3 用伽辽金方法求解
让我们用伽辽金方法再考虑这个问题。(2. 21)式的残数加权方程为
21
20(1)i d w x dx dx
Φ
0??=∫ (2.31)
用(2．26)式中的相同方式展开Φ，并令 2
1w x x =?和3
2w x x =?，它们是与和有关的展开函数，可以得到：
3c 4
c 34111
0324c c +
?= (2.32)
341423
02560
c c +?= (2.33)
求解上面两式，再次得到312c =和41
6
c =。这与精确解相同。事实上，只要试探函数组成被
研究问题的完备基，我们总能得到精确解，即使我们应用其它残数加权方法(例如点配置
法，或于域配置法、或最小二乘法)，情况也是如此。当试探函数不能组成完备基时，我们 将得到近似解。在这种情形下，用不同的残数加权方法将得到不同的近似解。
2．2．4 用子域展开函擞求解——有限元方法
从上面描述的步骤中可以看出：在里兹方法和伽辽金方法中，在整个解城内找出能表示或至少近似表示问题真实解的试探函数，这是非常重要的步骤。然而对于许多问题，这个步骤是十分困难的，对二维和三维问题尤其如此。为了克服这种困难，我们可将整个区域划分成小子域，并应用定义在每个子域上的试探函数。因为子域是小区域，因而在每一子域内，函数()x Φ的变化不大，所以，定义在子域上的这些试探函数通常比较简单。为了说明这个过程，让我们重新考虑(2．21)～ (2．23)式所定义的问题。
首先，我们将整个解区域(0，1)分成三个子域，并用12(,)x x 、23(,)x x 和34(,)x x 表示，
其中1x 和4x 是端点，即和10x =41x = (见图2. 1)。尽管不必等分子域，但在此我们采用等分方法，因而另外二个点选为2x =1／3和3x =2／3。假设()x Φ在每个子域上线性变化，
11
11()i i i i i i i
i x x x x x x x x +?++?+?x ?Φ=Φ+Φ (2．34) 其中，1i i x x x +≤≤，i =l，2，3。这里的i Φ是待定的未知常数。仔细观察发现：当i x x =时，
代表i Φ()x Φ的值。从边界条件(2．22)式和(2．23)式，我们得到10Φ=和。然后，
我们以用里兹变分或伽辽金方法确定剩下的两个常数41Φ=2Φ和3Φ。
10x =2x 3x 41x = 图2.1 分成三个子域的解区间
首先，让我们用里兹方法，即对泛函求极小值。为此，将(2．34)式代入到(2．25)式中，得到
F ()1123
111111112i i i i x x i i
i i i i x x i i i i i i i x x x x F dx x x x x x x x +++++=++??????Φ?Φ????=++Φ+Φ??????????????
∑∫∫1dx +? （2.35） 求积分后得到
()2
3
11111111212112333i i i i i i i i i i i i i F x x x x x x x x +++++=+??
??Φ?Φ??????=?+Φ+++Φ++?????????????????
∑13(2.36)
将i x 、和的值代入到上式中后，得到 1Φ4Φ 2
2
2323234224333992F =Φ+Φ?ΦΦ+
Φ?Φ+97
(2.37) 为了求的极小值，取对和F F 2Φ3Φ的偏导数，并令它们分别为零。我们得到
232
4
639F 0?=Φ?Φ+=?Φ (2.38)
233
22
369F 0?=?Φ+Φ?=?Φ (2.39)
从上面两式解得
214
40
81
81
Φ=Φ=
3 (2．40) 上面的结果也可用伽辽金方法获得。在这种方法中，我们将(2.39)式中的展开函数选 为(2．36)式中使用的加权函数。这种选择给出：
i w 111
111i i i i i i i i i i i
x x x x x x x w x x x x x x x ???+++??≤≤???=?
??≤≤???当时
当时
（2.41）
其中2，3。然后，将(2．34)式和(2．41)式代入到(2．31)式中，得到和的线
性方程组，在这样做之前，我们注意到：由(2．34)式表示的i =2Φ3ΦΦ
只有一次微分。因此，我们必须降低(2．34)式中对微分的阶数。这可通过分部积分来实现，即将导数转移至加权函Φ
数
1
1
1
1
1
1
22i i i i i x x x i i i
x x x dw d d d w dx w dx dx dx dx
++?????ΦΦ=?????
∫
∫
i dx +Φ
(2.42)
因为在1i x ?和1i x +处的为零，所以，(2.26)式可重新写成
i w
1
11
1
(1)i i i i x x i i x x dw d dx x w dx dx dx
++??Φ
0++∫
∫ = (2.43)
将(2．34)式和(2．41)式代入上式，并进行积分运算(为方便起见，在替代前将(2．39)式中的下标i 改成j)，得到
234
639Φ?Φ+
=0 (2.44) 2322
369?Φ+Φ?=0 (2.45)
这与(2.38)式和(2.39)式中给出的方程相同，当然，其解也与(2.40)式中的解相同。
一旦得到i x 处的解，其它点的解可从(2.34)式的线性插值求得。图2．2给出了这个例子的结果，作为比较，图中也给出了精确解。可以看出:尽管在i x 处得到了精确解(不会始终
如此)，但在其它点存在着小偏差。当子域数增加时，这种偏差会变得更小。
图2．2 数值结果与精确解的比较
上面刚刚描述的求解过程，正是有限元方法的求解过程。应用里兹方法的过程通常称为里兹有限元方法，或变分有限元方法，而应用伽辽金方法的过程通常称为伽辽金有限元方法。 从上面的例子可以看出：有限元方法与经典里兹方法和伽辽金方法的不同之处是在试探函数的公式上。在经典里兹方法和伽辽金方法中，试探函数由定义在全域上的一组基函数组
成。这种组合必须能够(至少近似)表示真实解，也必须满足适当的边界条件。在有限元方法中，试探函数是由定义在组成全域的子域上的一组基函数构成。如前所述，因为子域是小的，所以，定义在子域上的基函数能够十分简单。
在此，也许会问：正如上面例子所示，既然用全域基函数求解问题比用于域基函牧求解容易得多，那么，为什么我们还不厌其烦地介绍有限元方法呢?此问题的答案有两部分。第一，在上面的例子中，我们处理的是一维问题，对几乎所有的一维问题，确实总可以找到所需的试探函数。然而，当处理二维或三维问题时，在全域上寻找所需的试探函数；是非常困难的，对不规则边界问题尤其如此。第二，在上面的例子中，我们用一支笔、一张纸靠手算来求解，对非常筒单的问题才可如此，对实际感兴趣的复杂问题，必须应用计算机，即必须用计算机语言写出描述求解过程的程序，井在计算机上运行程序来求出问题的解。正如以后将要看到的那样，有限元方法非常适合于编程计算。用这种方法，可以写出能处理任意边界和许多不同问题的通用程序。
从上述讨论可以看出，正是应用了子域基函数的思想才使得处理复杂边值问题成为可 能；正是由于计算机技术的出现才使得这种方法成为可能。正是这种原因，有限元方法通 常是在绝大多数物理和工程问题方面(包括结构分析，流体力学，振动，热传递和电磁学) 应用最广泛的一种计算机辅助分析方法。
2．3 有限元方法的基本步骤
在2．2节中简要介绍有限元方法后，我们现在用比较抽象的形式定义这种方法，然后 描述该方法的基本步骤。正如本章开始时所述，有限元方法是求边值问题的数值过程。正如在例子中已见到的那样，该方法的原理是用许多子域来代表整个连续区域。在子域中，未知函数用带有未知系数的简单插值函数来表示。因此，无限个自由度的原边值问题被转化成了有限个自由度的问题，或者换句话说，整个系统的解用有限数目的未知系数近似。然后，用里兹变分或伽辽金方法得到一组代数方程(即方程组)。最后，通过求解方程组得到边值问题的解。所以，边值问题的有限元分析应包括下列基本步骤： (1)区域的离散或划分； (2)插值函数的选择； (3)方程组的建立； (4)方程组的求解。
下面我们较详细地讨论每一步骤。 2．3．1 区域离散
假设区域为?。在任何有限元分析中，区域离散是第一步，也许是最重要的一步，因为 区域离散的方式将影响计算机内存需求、计算时间和数值结果的精确度。在这一步骤中， 全域口被分成许多小区域，用表示，这里M 表示于域总数。这些子域 (1,2,3,,e
e ?= )M 通常被称为单元。对于实际上是直线或曲线的一维区域，单元通常是短直线段，它们连接 起来组成原来的线域(见图2．3())。对于二维区域，单元遭常是小三角形或矩形(见图2．3 a (b))。当然，矩形单元最适合离散矩形区域，而三角形单元可用来离散不规则区域。在三维求解中，区域可划分成四面体、三棱柱或矩形块（见图2. 3(c)），其中四面体是最简单、最适合离散任意体积区域的单元。我们指出：线性单元、三角形单元及四面体单元是用直线段和平面块建立曲线或面模型的基本一维、二维、三维单元。在图2. 4中，我们给出了表示二维和三维区域有限元离散化的两个例子。
图2. 3 基本离散单元
(a) 一维; (b)二维; (c)三维
图2.4 有限元离散化例子
(a) 用三角形单元的二维情形; (b) 用四面体单元的三维情形
在多数有限元解中,问题是用与单元有关的结点上的未知函数Φ表达的。例如，线性元有两个结点，每端点一个，线性三角形单元有三个结点，各顶点上一个；而线性四质体单元有四个结点，各角点上一个。在此有必要描述一下这些结点。一个结点的完整描述应包括它的坐标值、局部编码和全局编码。结点的局部编码表示它在单元中的位置，而全局编码表示它在整个系统中的位置。标明坐标值是一项直截了当的工作，而结点和单元的编码需要一定的策略。
我们将要看到：有限元公式通常得到带状矩阵，其带宽由一个单元中两个结点的全局编码之差的最大值决定。因此，如果用带状矩阵求解方法来解最终的矩阵方程，那么，通过适当的结点编码使得带宽最小，就可大大减少计算机的存储量和运算时间。然而，当带宽不必最小时，编码方案可以任取，通常选择使编程简单的方法。
因为区域离散过程完全可以与其它步骤分开，所以，通常将它当作一项预处理工作：许多较完善的有限元软件包具有将任意形状线，面和体离散成相应单元的能力，同时也能提供最优全局编码。
2．3．2 插值函数的选择
有限元分析的第二步是选择能近似表达一个单元中未知解的插值函数。通常，插值函数可选择为一阶(线性)，二阶(--次)、或高阶多项式。尽管高阶多项式的精度较高，但通常得到的公式也比较复杂。因此，简单且基本的线性插值仍被广泛采用。一旦选定了多项式的阶数，我们就能导出一个单元中未知解的表达式。以单元为例，得到下列形式：
e (2.46) 1
{}{}{}{}n
e e e e T e e T j j
j N N N =Φ=Φ=Φ=Φ∑ e 式中，n 是单元中的结点数；是单元中j 结点的e j ΦΦ值，是插值函数，通常也称为展开函数或基函数。的最高阶被称为单元的阶。例如，若是线性函数，则单元是线性单元。函数的重要特征是：它们只有在单元e 内才不为零，而在单元e 外均为零。 e
j N e
j N e
j N e e
j N 2．3．3方程组公式的建立
第三步(也是有限元分析中的主要步骤)是导出方程组的公式。类似于前几节，里兹变分和伽辽金方法均可用于这种目的。让我们首先考虑里兹变分公式。 1．里兹方法的求解公式
让我们再考虑(2．1)式中定义的问题。为简单起见，假设问题是实数值的。(2．5)式给出的泛函F 可表示为
1
()()M
e e e F F =Φ
=Φ∑ (2.47) 式中M 是组成全域的单元数，而
1()2e e e e e e e F L d ?
?Φ=ΦΦ??Φ?∫∫ f d (2.48) 将(2.46)式代入(2.48)式，得到 1{}{}{}{}{}{}2
e e e
e T
e e T e e T e F N L N d
f ??=
Φ?Φ?Φ?∫∫N d (2.49)
可将其写成矩阵形式
1{}[]{}{}{}2
e
e T e
e e T F K =
ΦΦ?Φe b ? (2.50) 其中[是n ×n 矩阵，{是n ×l 的列向量，它们的元素为
]K }b e ij e e e i i K N LN d ?
=∫ (2.51)
和
e
e e i i b fN d ?=?∫
(2.52)
我们注意到：因为是自伴算符，所以，单元矩阵[是对称的。将(2.50)式代入(2.47)式，得到
L ]e
K 11()({}[]{}{}{})2
M
e T e e e T e e F K =Φ
=ΦΦ?Φ∑ b (2.53) 通过求和运算，井采用全局结点编码，上式可写成
1
{}[]{}{}{}2
T F K =ΦΦ?ΦT b (2.54)
式中[是N ×N 对称矩阵，N 是未知量或结点总数，{]K }Φ是N 1×的未知向量，其元素是未知
展开系数，{是N 的已知向量。应用驻点条件，即令对}b 1×F i Φ的偏导数为零：
()1
101,2,3,2N
ij ji j i j i F K K b i N =?=+Φ?==?Φ∑ (2.55)
即可得到方程组。因为[是对称的，]K ij ji K K =，因此，(2．60)式变成
1
01,2,3,N
ij j i j i F
K b i =?=Φ?==?Φ∑ N b (2.56) 或写成矩阵形式:
[]{}{}K Φ= (2.57)
与上式等价但推导方式略微不同的推导过程是首先取对e F e
i Φ的导数：
{}{}1,2,3,,e e e i
e e T e e
i e F N L N d fN d i n ???=?Φ??=?Φ∫∫ (2.58)
上式可写成矩阵形式
[]{}{e e e e
e F K ?????=Φ????Φ????
}b (2.59) 式中
12
,,n
T
e
e
e
e e e F F F F e e
??
??????=?????Φ?Φ?Φ?Φ??????
为了获得方程组，必须首先求出{
}F
??Φ
，这里
12
,,N T
F F F
F ????????=?
????Φ?Φ?Φ?Φ??????
既然只有与结点i 直接相连的单元才对{
}
i
F ??Φ有贡献，那么，应用局部和全局结点编码的
关系，将每一单元的{
}e
e
F
??Φ扩展成N 1
×的列向量，然后再将它们相加，得到
1e M e e F F =?
???????=????Φ?Φ???????
∑ (2.60) 这里上划线表示该向量是扩展的或增广的。然后用驻点条件，得到方程组
11
([]{}{}){0}e M M e e e
e e e F F K b ==????????==Φ??????Φ?Φ??????∑∑=] (2.61)
这里，在求和号内的所有向量和矩阵都是增广的。更具体地说，应用局部和全局结点编码间的关系，通过添零，将[叫扩展成了N e
K ×N 的矩阵[e
K ]。同样，{}e
Φ和{}e
b 也都增广到N 的列向量。后几章中将详细讨论所有这些运算。结果，(2.61)式也可写成(2.57)式。
显然，这种不同的推导过程不如前面推导漂亮。然而，正如我们马上就会看到的那样，因为它和伽辽金公式有相似性，所以，我们还是介绍了这种推导过程。 1× 2．伽辽金方法的求解公式
上述方程组也可通过伽辽金方法导出。对于(2．1)式，第个单元的残数加权为 e ()1,2,3,e
e e e i i ,R N L
f d i n ?
=
Φ??=∫
(2.62)
将(2．51)式代入(2．67)式后得到 {}{}1,2,3,,e
e e
e e T e e i i i R N L N d
f N d i n ?
?
=?Φ??
=∫
∫ (2.63)
也可写成矩阵形式
{}[]{}{}e e e e R K =Φ?b T
(2.64)
这里12{}[,,,]e
e
e
e n R R R R = ，矩阵元素和分别与(2．56)式和(2．57)式中的相同。因为这里的算符不必是自伴的，所以，单元矩阵[]也不必是对称的。既然与一结点有关的展开函数(因而加权函数)遍及所有和该结点直接相连的单元，那么，与结点有关的残数加权ij K e
i b L e
K i i R ，是对所有直接和结点i 相连的单元求和。所以，利用局部和全局编码的关系，我 们可以扩展(2.64)式，然后将它对每一单元求和，得到 1
1
{}{}([]{}{})M
M
e
e
e e e e R R K
b ===
=Φ?∑∑ (2.65)
式中，12{}[,,,]T
n R R R R = 。我们再次指出：(2.65)式中求和号内的所有向量和矩阵都按前面描述的方式扩展了。令(2.65)式等于零，可得到方程组
1
([]{}{}){0}M
e
e e e K
b =Φ?=∑ (2.66)
它也可写成(2.57)式的形式。
在方程组(2.57)式能求得定解之前，必须应用所需的边界条件。有两类边界条件经常出现：一是狄利克雷边界条件，像(2．22)式和(2．23)式那样，它给出了边界处的中值；另一类是齐次诺曼边界条件，它要求边界处中的法向导数为零。第一类边界是必要边界条件，它必须显式地强加在计算中，第二类边界条件通常在求解过程中隐含地自动满足。正是这种原因，第二类边界条件通常称为自然边界条件。我们将在后面章节中详细讨论这两类边界条件和一般边界条件。
可以看出，在这一步骤中，我们实际上有三个子步骤。首先，应用两种方法中的任一种，写出单元方程(2.50)式或(2.64)式。其次，我们将单元方程对所有单元求和，得到方程组，这个过程叫做组合。最后，我们应用边界条件来得到方程组的最终形式。注意：在用计算机实现此计算的过程中，这三个子步骤通常不是分立的，相反，它们相互交织在一起，单元矩阵的生成和边界条件的强加通常发生在组合过程中。 2．3．4 方程组的求解
方程组的求解是有限元分析的最后一步。最终的方程组是下列两种形式之一： []{}{}K b Φ= (2.67) 或者
[]{}[]{}A B λΦ=Φ (2.68)
方程(2.67)式是确定型的，它是从非齐次微分方程或非齐次边界条件或从它们两者兼有的问题中导出的。在电磁学中，确定性方程组通常与散射、辐射以及其它存在源或激励的确定性问题有关。而方程(2.68)式是本征值型的,它是从齐次控制微分方程和齐次边界条件导出的。在电磁学中,本征值方程组通常与诸如波导中传输和腔体中谐振等无源问题有关。在这种情形下，已知向量{为零，矩阵[可写成[]}b ]K []A B λ?的形式，这里λ表示未知的本征值。一旦解出{的方程组，就能计算出所需要的参数，例如电容、电感、输入阻抗和散射 或辐射图等，并能用曲线、图形或彩色图片等形式表示结果，这些形式更有意义并易于解释。这个最后步骤也称为后处理过程，可完全同其它步骤分开处理。 }Φ
第三章 一维有限元分析
在第二章中，我们介绍了有限元方法。本章我们将首先按照该方法的基本步骤来处理，一般一维边值问题。我们指出：有限元方法没有广泛应用于一维电磁问题中，这是因为只有少数这类问题的模型复杂到需要做数值解。然而，因为其简单性,一维问题是用来说明公式建立过程的最好途径。 3. 1 边值问题
考虑用下列微分方程定义的边值问题： (0,d d )f x dx dx αβΦ???
+Φ=∈????
L (3．1)
这里是未知函数，Φα和β是与区域物理性质有关的已知参数，f 是已知源或激励函数。标准一维拉普拉斯(Laplce)方程、泊松(Poisson)方程和亥姆霍兹(Helmholtz)方程是(3. 1)式的特殊形式。
Φ的边界条件由下面两式给出
x p =Φ
= (3．2)
以及
x L
d q dx αγ=Φ??+Φ=???? （3. 3） 式中，p 、γ和q 是巳知参数。在上面两式中，通常称(3．2)式为第一类边界条件，即狄利克雷(Diriehlet)条件；而(3．3)式称为第三类边界条件。第二类边界条件，条件，是(3．3)式在0γ=条件下的特殊形式。在区域两端点处取两类不同边界条件的理由 仅仅是为了说明它们的不同处理方法，并保持分析的普遍性。
在区域内某点上，例如在d x 处（0d x L <<），如果α不连续或突变，则需要满足连续性条件，即
d d x x x x =+=?Φ
=Φ
（3. 4）
式中表示0d x x =+x 从右侧趋于d x ，而0d x x =?表示x 沿左侧趋于d x 。
3．2 变分公式
为了应用变分方法列出有限元方程组，我们首先需要建立所需的变分原理。对上节定 义的问题，可以证明：通过求解下面定义的等效变分问题，能得到该问题的解：
()00
x F p δ=?Φ=??Φ=??
（3. 6）
式中
()2
2
00122L L x L d F dx f dx dx γαβ=??Φ????Φ=+Φ?Φ+Φ?Φ????2q ??????????
∫∫ （3. 7） (3. 6)式的含义是：在给定狄利克雷边界条件下，求泛函()F Φ的驻点。我们指出：为了书
写方便，我们已经用Φ代替以前泛函中的Φ
，在本书的其余部分也将这样做。但是，需要记住:每当出现在泛函或残数加权方程中时，它是试探函数的近似解，而不是精确解。
Φ为了证明上面描述的变分公式，我们取()F Φ对Φ的第一变分：
()()0
L
L x L d d F dx q dx dx δδαβδγδ=?ΦΦ?
????Φ=+ΦΦ+Φ?Φ?Φ????????????????
∫
∫f dx δ （3. 8）目前，假设α在全域中连续，对上式右边第一项作分部积分后，(3．8)式可写成
()()0
00
x L
L
x L x L d d d F dx dx dx dx q f dx
δαβδαγδδδ===?Φ?Φ????
Φ=?+ΦΦ+Φ????????????
+Φ?Φ?Φ????∫
∫ (3．9) 皆为在处，Φ值固定，所以，0x =0
0x δ=Φ
=。因此，
()00
L
x L L d d d F dx dx dx dx
f dx
δαβδαγδ=q δ?Φ??Φ?????Φ=?
+ΦΦ++Φ?Φ???????
??????????Φ∫∫ (3.10) 强加驻点条件0F δ=，我们得到
0L
x L d d d f dx q dx dx dx
αβδαγδ=?Φ??Φ?
??
???+Φ?
Φ++Φ?Φ??????????????∫
=?? （3．11）
因为δΦ是任意变分，(3．11)式中的积分和边界项必须全部为零。因此，
0d d f dx dx αβΦ??
?
+Φ?=????
(3．12) 0d q dx
α
γΦ
+Φ?= (3．13) 将(3．12)式和(3．13)式与(3．1)式和(3．3)式进行比较，我们得出结论：变分问题(3．6)式的解确实也是(3．1)～(3．3)式所定义的边值问题的解。方程(3．12)式称为(3．7)式给出泛函的欧拉(Euler)方程。正如我们所看到的那样，条件(3．13)式在求泛函的极大或极小的过程中自动满足，所以，此条件被称为自然边界条件。另一方面，狄利克雷边界条件 (3. 2)式必须被强加上，因此，它被称为必要边界条件。
()F Φ在我们的证明中，已假设α为连续函数。现在考虑。在d x 处不连续的情况。在这种情形下，我们不能直接在全域上进行分部积分。然而，我们能够将区域()0,L 分成两个子区域，一个是()0,d x 另一个是(),d x L 。然后，在每一个子域上进行分部积分，得到(3.11)式，其左边有以下两个附加项:
0d x x d dx αδ=+Φ??Φ???? 0
d x x d dx αδ=?Φ???Φ???? (3.14) 既然(3．4)式要求中在d x 处连续，因此，δΦ也连续，因而(3．14)式可写成
00d d x x x x d d dx dx αα=+=???ΦΦ??????
δ?Φ?????????????
? (3. 15) 另外，既然δΦ是任意变分，当强加驻点条件时，可得到除(3．12)式和(3．13)式以外
的另一自然条件
00
0d d x x x x d d dx dx αα=+=?ΦΦ?????=???????? (3．16) 这在求泛函极值的过程中隐含地满足。这个条件即为连续性条件(3．5)式。而另一连续
世条件(3．4)式是必须强加的必要条件。对于这种情形，严格说来，变分表达式(3. 6)式中应该补上(3．4)式，但是，因为(3. 4)式通常在Φ的有限元展开中已被满足，所以，通常在 (3．6)式中忽略(3．4)式。
3．3有限元分析
在这一节中，我们按照第二章中列出的有限元方法的基本步骤写出问题的解。为简单起见，选用线性插值基函数。
3．3．1 离散化和插值
按照有限元方法，第一步是将解域()划分成短线段的小子域。设 (＝1，2，3,
，0,L e
l e M )表示第个单元的长度，M 表示单元总数。进一步假设e i x (i =1，2，3， ，N)表
示第i 个结点的位置，10,N x x ==L 。如图3．1所示，单元和结点均按从左至右的顺序编码。结果，无需引进局部编码系统。但是，为了使公式推导与二维和三维情形保持一致，在公式建立过程中也采用局部编码系统。在下面，上标e 表示有关量的下标为局部编码，而所有其它量的下标是全局编码。所以，局部和全局系统的关系为。
1e e x x = 和 2e
e 1x x += (3.18)
这里＝l，2，3，…，M。
e 有限元方法的第二步是选择插值函数，为简单起见，我们应用线性插值。因此，在第e 个单元内，()x Φ可以近似为
()e e e x a b x Φ=+ (3．19)
式中，和是待定常数。对于线性单元，每一单元均有两个结点:一个位于e
a e
b 1e
x ，另一个位于2e
x 。定出（3.19）式在这两个结点处的值，得到
112
2
e
e e e e e
a b x a b x
Φ=+Φ=+e
e
单元 1 2 M-1 M
x
结点 1 2 3 N-2 N-1 N x=0 x=L
（）
a
e
1 2
（）
b
图 3. 1 划分为线性单元的一维区域
（a ） 单元和全局结点编码； （b ） 局部结点编码的线性单元
式中，和分别表示1e
Φ2e
Φ()x Φ在1e
x 和2e
x 处的值。然后解出和，在将它们带回（3.19）式，得到
e
a e
b （3.20）
2
1
()()e
e
j j x N
x =Φ=
Φ∑e j 式中，和表示插值函数或基函数，记为
1e N 2e
N 21
e
e
e x x N l ?= 和 12e e
e x x N l ?= （3.21）
这里。显然，2e e l x x =?1e ()e e
j i ij N x δ=。当i j =时，1ij δ=；当i j ≠时，0ij δ=，我们用图3. 2表示这种关系。
图3. 2 一维线性插值函
3.3.2 用里兹方法建立公式
有限元方法的第三步骤是建立方程组公式。为此，首先要利用里兹方法或伽辽金方法来建
立单元方程。让我们首先考虑里兹公式。
1. 单元方程推导
为简单起见，让我们暂时将边界条件（3. 3）式限于齐次诺曼条件0q γ==。我们以后 会考虑一般边界条件。在这种限制条件下，(3．1)式的泛函可写成 ()(1
M
e
e
e F F =)Φ=Φ∑ (3.22)
式中
()()2211
22
12e e e e e x x e e e e x x d F d dx αβ????ΦΦ=+Φ?Φ??????????
∫∫x fdx (3．23) 将(3．20)式代入上式，然后取对e F e
i Φ的导数，得到
2211
21e e e e e e e x x j e e e
i j i j e x x j i dN dN F N N dx N fdx dx dx αβ=???=Φ+??????Φ??
∑∫∫e i (3．24) 上式也可写成矩阵形式
{}{}e e e e
e i F K ?????=Φ?????
?Φ??
b (3．25) 在上式中，{
}12,T
e
e e ?Φ
=ΦΦ???，并且
2
1
e
e ij e e x j e
i i j x dN dN K dx dx αβ??=
+????
∫
e e N N dx ??dx ?? (3．26) (3．27)
2
1
e e x e
e i
i x b N f =
∫
我们指出：是对称的，并且如果在每一单元内和是常数或可用常数近似表示，那么，
可用解析法求其矩阵元素。结果是
e K ?? 11
22
3e
e
e
e e
e l K K l
αβ==+ (3．28)
12
21
6e
e
e
e e
e l K K l
αβ==?+ (3．29)
对于，我们同样得到
e
i b 1
2
2
e
e e e
l b b f == (3.30)
这里，e α、e β和e
f 表示α、β和f 在第个单元内的值。 e 2．组合成方程组
将(3．25)式给出的单元方程对所有单元求和，然后强加驻点条件，得到方程组：
{}{}(
{}11
0e M M e e e
e e e F F K b ==??????????==Φ???????
?Φ?Φ??????∑∑)
= (3．31) 式中，通过添零将求和号后的矩阵和向量进行了扩展。正如在第二章中所看到的那样，这
个过程叫做组合。为了清楚地说明组合，让我们考虑图3．3所示的三个单元和四个结点的
单元 1 2 3
x
结点 1 2 3 4
图3. 3 表示组合过程的例子
情形。应用局部和全局结点编码的关系，我们能够将e K ????扩展成4×4矩阵，并将{}
e Φ扩展成4×1列向量。例如，第一单元的e K ????
和{}
e
Φ能够扩展成 ()
()()
()()1111
12111212200000
0000
00K K K K K ??
??????=??????????, (){}
()()1111200??Φ??
??ΦΦ=????????
（3.32）它们的和变成
()()
{}
()()()()()()()()1111111122
11111121122200K K K K K ??Φ+Φ????Φ+Φ??Φ=????
??
????
(3.33)
同样可得到
()(){}
()()()()()()()()22222211112222
2221122200K K K K K ??
??Φ+Φ????Φ=????
Φ+Φ??????
（3.34） 以及
()()
{}
()()()()()()()()3333331111223333211
22200K K K K K ??
??????Φ=????
Φ+Φ????
Φ+Φ?? (3.35)
第1讲 抛物问题有限元方法 1、椭圆问题有限元方法 考虑椭圆问题边值问题： （1） ()? ??Ω?∈=Ω ∈=?-x u x x f u ,0, 问题（1）的变分形式：求()Ω∈1 0H u 使满足 （2） ()()()Ω∈?=1 , ,,H v v f v u a ()v u a ,的性质，广义解的正则性结果。 区域Ω的剖分，矩形剖分，三角剖分，剖分规则，正则剖分条件，拟一致剖分条件。 剖分区域上分片k 次多项式构成的有限元空间()Ω?1 0H S h 。 h S 的逼近性质，逆性质： ∞≤≤≤≤≤≤-+-+p k k m u Ch u I u p k m k p m h 1,1,0,,11, h h p m h q n p n l m q l h S v l m q p v Ch v ∈?≤∞≤≤≤---,,,1,), 0(max , 这里，h h S u I ∈为u 的插值逼近。 问题（2）的有限元近似：求h h S u ∈使满足 （3） ()()h h h h h S v v f v u a ∈?=, ,, （3）的解唯一存在，且满足f M u h ≤1 。 （3）的解()()∑== N i i i h x u x u 1φ所满足的矩阵方程（离散方程组）形式： ()()N j f u a j N i i j i Λ,2,1, ,,1 ==∑=φφφ （4） f u K ? ?= 刚度矩阵()() N N j i a K ?=φφ,的由单元刚度矩阵组装而成。 -1H 模误差分析：由（2）－（3）可得
（5） h h h h S v v u u a ∈?=-,0),( 由（5）可首先得到 ()()11 21 ,,u I u u u M u I u u u a u u u u a u u r h h h h h h h --≤--=--≤- 则得到 （6） 1,1 11 ≥≤-≤-+k u Ch u I u C u u k k h h 2L －模误差分析 设21 0H H w I ∈ 满足 h h u u C w w in u u w -≤=Ω-=?-Ω ?2, 0,, 用h u u -与此方程做内积，由（5）式和插值逼近性质得到 ()()w u u A u u w A u u h h h ,,2 -=-=- ()h h h h h h h u u u u Ch w u u Ch w I w u u C w I w u u A --≤-≤--≤--=1 2 1 11, 再利用-1H 模误差估计结果，得到 （7） 1,1 11 ≥≤-≤-++k u Ch u u Ch u u k k h h 最优阶误差估计和超收敛估计概念。 当()t u u =与时间相关时（为抛物问题准备），由（5）式得 （8） h h h t h S v v u u a ∈?=-,0),)(( 利用（7），类似分析可得 （9） ()()1,1 11 ≥≤-+-++k u Ch u u h u u k t k t h t h 2、抛物问题半离散有限元方法 考虑抛物型方程初边值问题： （10） ()()()()()()()()??? ? ???Ω∈=Ω??∈=Ω?∈=?-??x x u x u T x t x u T x t x t f u t u ,,0,0,,0,0,, ,0 （10）的变分形式：求 ()())()0(, ],0(:01 0x u u H T t u =Ω→ 使满足 （11） ()()()()Ω∈?=+1 , ,,,H v v f v u a v u t
第一章 有限元概貌与发展 有限元方法是近似求解数理边值问题的一种数值技术。这种方法大约有60年的历史。它首先在本世纪40年代被提出，在50年代开始用于飞机设计。后来，该方法得到了发展并被非常广泛地用于结构分析问题中。目前，作为广泛应用于工程和数学问题的一种通用方法，有限元已相当著名。 有限元法应用于电磁场中，最先是用结点上的插值基函数来表征该结点上的矢量电场或磁场分量的，称为结点有限元。但是，在使用结点有限元进行电磁仿真时，会有几个严重的问题。首先，非物理的或所谓伪解可能会出现。其次，在材料界面和导体表面强加边界条件很不方便。再次，处理导体和介质边缘及角也很困难，这是由与这些结构相关的场的奇异性造成的。在这些问题中，最后一个问题比其它两个问题更严重，因为它缺少通用的处理方法。即使对前两个问题，目前的处理状况也不能完全令人满意。因此，有必要探讨其它的可能性或其它方法，而不仅仅是改进，从而将电磁场有限元分析引入一个新的时代。 幸运的是，一种崭新的方法已经被发现。这种方法使用所谓矢量基或矢量元，它将自由度（未知量）赋予棱边而不是单元结点。因为这个原因，它也叫棱边元（edge element ）。虽然Whitney 早在35年前就描述过这些类型的单元，但它们在电磁学中的应用及其重要性直到前几年才被认识到。在80年代初，Nedelec 讨论了四面体和矩形块棱边元的构造。Bossavit 和Verite 将四面体棱边元应用于三维涡流问题。Hano 独立地导出了矩形棱边元，并用于介质加载波导的分析。Mur 和de Hoop 考虑了非均匀媒质中的电磁场问题。Van Welij 和Kameari 应用六面体棱边元进一步考虑了棱边元在涡流计算中的应用。Barton 和Cendes 将四面体棱边元应用于三维磁场计算，同时，Crowley 提出了一种更复杂的单元类型，即所谓的协变（covariant ）投影单元，它允许单元带有弯曲的棱边。在所有这些工作中，已经证明：棱边元没有前面提到的所有缺点。因为这些缺点困惑了研究者多年，可以想象，棱边元的重要性很快就被认识到了，因此，在过去的几年中，开展过大量的研究也获得了很多成功的应用。 棱边元分析时谐电磁场问题通常是基于矢量波动方程的，但这也常常引发一个问题:作为电场强加条件的[]0E ε?=i ，低频时在波动方程中不起作用。这一点使得所产生的有限元矩阵严重病态。尽管这一缺陷并不影响基于场形式波动方程的边棱元方法在高频电磁分析中的应用，但是在使用自适应网格剖分产生的一些小单元中，会局部的逼近静态的情况，同样会产生病态的矩阵的。这种病态的矩阵不利于直接法的求解，并且也大大的降低了迭代法的收敛速度。近年来发展了一种矢量标量位有限元，将边棱元与结点元相结合，能在很宽的频带内直接强加电场散度条件，同时有效的改善有限元矩阵的性态，尤其适合于结合迭代法求解。 为了降低有限元矩阵未知量的数目，不少学者对高阶有限元也作过了大量的研究工作。它的主要思想就是利用高阶的基函数对未知的场获得更精确的逼近，或者说在较稀疏的网格上获得一般的线性插值基函数在稠密网格上同样的精度。但这种高阶的有限元通常也会产生严重病态的矩阵，不利于快速求解，一直是个待解决的问题。 第二章 有限元方法入门 在本章中，我们首先回顾求解边值问题的两种经典方法，它们都包含着有限元方法的本质。然后，应用一个筒单的例子，介绍有限元方法。最后，我们不针对任何特定问题来描述该方法的基本步骤。 2．1 边值问题的经典方法 在本节中，我们首先定义边值问题，然后讨论边值问题求解的两种经典方法，一是里兹