全波数值计算方法之有限元法（FEM）

有限元法是一种将偏微分方程转化为线性代数方程组，进而求解边值问题的数值方法，最早由Courant于1973年提出，用来求解势论中的变分问题，自此以后该方法得到了极大的发展，并被广泛应用于结构力学分析以及其他领域问题的求解。由于有限元法不仅能适应各种复杂结构，而且计算精度高，因此成为了处理微波工程和电磁学问题的一种通用方法。

本文从有限元法的一般原理出发，推导出矢量场的边值问题，从而建立有限元公式。

1.1 有限元法的一般原理

使用加权残差法或变分法可以建立有限元公式。加权残差法是从边值问题的偏微分方程出发，而变分法则是从边值问题的变分形式出发。本小节采用加权残差法来建立有限元公式，

考虑如下偏微分方程：

式中， L为微分算子， 为待求未知解或者称作自由度，f 为激励函数。为了寻找 的解，首先使用如下一组基函数将其展开：

式中， ( =1,2,…,N )为基函数，其线性组合可表示未知解； 为相应的未知展开系数。加权残差法确定 的思想是：将式（2-2）代入式（2-1）中，然后乘以加权函数并在整个求解区域 中进行积分，得到式：

给定一组加权函数，式(2-3)就定义了一个代数方程组，在满足边界条件的要求下求解该代数方程组即可得到 。加权函数 一般等于 ，以此构建有限元公式的过程称为伽辽金法。式(2-3)变成

或

式中，

对于自共轭问题，则有

因此有 =，即式对应的线性系统的系数矩阵是对称矩阵。

在构建有限元公式过程中，最关键的一步是找到一组可以用来展开未知解的基函数，但对于不规则形状的二维三维问题，这一步骤极其困难。因此有限元法的基本思想是将求解区域 划分为许多子域，称为有限单元(有限元)，然后使用简单的基函数来近似单元内的未知解。

在构建有限元公式过程中，最关键的一步是找到一组可以用来展开未知解的基函数，但对于不规则形状的二维三维问题，这一步骤极其困难。因此有限元法的基本思想是将求解区域划分为许多子域，称为有限单元(有限元)，然后使用简单的基函数来近似单元内的未知解。

1.2 矢量场的边值问题以及有限元公式建立

1.2.1 边值问题

在介电常数为 、磁导率为 的区域 中，需要求解由电流密度 产生的电场强度 ，求解二维或三维区域 。求解服从给定边界条件的Maxwell方程组：

消去式(2-9)中和式(2-10)中的 ，可得到关于 的矢量波动方程为

式中，和分别为相对磁导率和相对介电常数；和 分别为自由空间中的波数和本征阻抗。

处理两种不同的电场边界条件-理想导体表面的齐次Dirichlet条件和阻抗表面的混合边界条件，将边界条件假设为：

式中，P是边界上的切向电场，是上的归一化表面阻抗， 表示边界上的边界源。

式(2-13)至式(2-15)所描述的边值问题通常很复杂，尤其是当求解区域 不规则以及相对介电常数 非均匀时，很难得到封闭形式的解析解，而数值算法中的有限元法因具有能够处理任意形状边界和非均匀媒质的能力成为了唯一选择。

1.2.2 有限元公式建立

采用加权残差法建立有限元公式，给式(2-13)乘以一个合适的加权函数，并在求解区域 中进行积分，从而得到原边值问题的弱式表达式，即

对式(2-16)应用矢量恒等式

后，其中

再应用高斯定理

然后，应用式(2-15)的边界条件，得到式(2-13)的弱式表达式为

式(2-20)作为有限元的基本方程，被用于有限元计算区域划分出的每个子区域。

1.3 矩阵填充与求解

将区域划成小的有限单元后，在划分的每一个小单元内，使用一组离散值插值可得到电场强度。在选择给定单元每一条棱边上电场的切向分量后，使用一组矢量基函数对其他位置的进行插值。下面以四面体单元中的场插值为例，采用式(2-21)进行插值：

式中，表示单元e中连接节点l和k的棱边上的电场切向分量，表示相应的插值函数或基函数。将四面体单元中与节点l 和k 相关的线性标量基函数分别表示为 和 ，则式(2-21)中的矢量基函数可以写成

式中，为连接节点l 和k 的带符号的边长。式(2-22)定义的基函数为矢量函数，相应的单元称为矢量元或棱边元。

由于每一个单元中的电场E 都可以用该单元中棱边上的切向电场分量进行插值，因此整个区域 中的电场可以表示为

式中，为除了上的棱边的所有棱边总数，为第j条棱边上的切向分量，为相应的矢量基函数。此外，为上的棱边总数，和 表示这些棱边上的切向电场和相应的基函数。

将式(2-23)代入式(2-20)中，使用矢量基函数作为加权函数，由于上有 ，故式公式参考此处在上积分为0，因此可以得到

式中，

式(2-24)可以写成紧凑形式

求解上式可得到。由于式(2-25)中单元之间的相互作用是局部的，因此是一个稀疏且对称的矩阵，它可以由稀疏矩阵求解算法高效求解。求出后，由式(2-27)可以求出中每一处的场值。

求解该矩阵方程一般有两种方法：迭代求解和直接求解。前者包括共轭梯度法和广义最小残差法等，迭代求解法是一种逼近精确解的近似方法，该方法因程序设计简单和消耗低计算机内存，常常被用于大型稀疏矩阵求解。而后者有高斯消元法、LU分解法、LDLT分解法等，这些方法通过有限步运算便可以得到精确解。