1.2.1 二维静电场有限元法初步

本节以一个二维静电场问题为例，首先给出问题描述及微分方程边值问题的数学表达，其次介绍采用有限元方法求解时的弱形式推导过程，采用线性基函数时的单元分析及整体合成过程，边界条件的处理方法，代数方程组的求解方法，最后介绍如何进行有限元结果的精度验证。

1.问题描述

如图1-3a所示二维区域（为方便起见，未按真实尺寸绘制），其中圆形电极（半径为0.2m）区域给定电位φ₁=1V，正方形区域（2m×2m）的四条边界为薄片导体且给定电位φ₀=0V。电极之间充满了均匀的线性电介质，相对介电常数ε_r=5，同时假设空间中不存在自由电荷。由于问题的几何模型、材料参数及边界条件关于x，y轴都对称，静电场及电位分布同样有相应的对称性，故可将该问题简化为图1-3b所示的1/4区域计算模型（甚至是1/8模型），从而通过缩小区域并施加对称边界条件来降低计算量，达到加快求解速度的目的，这在三维问题计算中也非常常见。

对于图1-3b给出的1/4区域计算模型，施加以下的边界条件。其中在对称边界S₁₁、S₁₂上，由于电场强度只存在切向分量，故法向分量为0

边界S₂₁、S₂₂上，给定电位φ=0；边界S₂₀上，给定电位φ=1。

由于均匀电介质中不存在自由电荷，所以此时直角坐标系下的静电场控制方程可以表示为拉普拉斯方程，即

2.微分方程定解问题

综上，该静电场边值问题可以归纳为

边界S₁₁与S₁₂上

边界S₂₁与S₂₂上

边界S₂₀上

3.有限元弱形式

有限元计算时，首先需要将上述偏微分方程定解问题转化为与之等价的积分形式的变分形式（variational form）或者弱形式（weak form）。在对弱形式进行离散求解时，需要对计算区域进行网格划分，其目的是将连续的计算区域离散为一定数量的网格单元（对于二维几何模型，单元可以有不同的几何形状，如三角形、四边形及多边形等），并在单元上定义自由度（即通常所说的未知数，或者有限元基函数下的展开系数），自由度与基函数的线性组合构成了试探函数空间，即有限元空间。有限元方法即在有限维试探空间中求解微分方程定解问题的近似数值解，该数值解满足有限元弱形式并与一组自由度唯一对应。对于弱形式，根据有限元网格进行空间离散并计算积分，可得到与有限元网格对应的所有自由度的代数方程组。求解方程组之后，即可得到所有单元自由度的具体取值，进而每个单元任意位置的电位函数值均可通过有限元基函数及展开系数计算得到。本节以一阶线性三角形网格单元为例说明有限元方法的具体求解过程。将图1-3b中的计算区域划分为n个三角形单元，得到图1-4所示的网格划分示意图。

推导变分方程时，可以采用加权余量法，将微分方程两边同乘以权函数（或叫检验函数）并进行分部积分[7]，最后取权函数与有限元空间基函数相同，即得到有限元变分方程或弱形式。将式（1-41）两边同乘以任意的W_i并在整个计算区域积分，应用分部积分公式可将原微分方程中关于φ的二阶偏导数转化为以下两个一阶导数乘积的形式：

式中，W_i为权函数；R_i表示当以W_i做权函数时微分方程的积分残量或余量。注意：若有φ满足微分方程边值问题，则其自然满足式（1-45）；反之未必。故式（1-45）是在弱意义[微分方程式（1-41）的加权积分为0]下建立的方程。但可以证明若式（1-45）的弱解同时满足充分的光滑性条件，则该弱解也满足强意义[逐点满足微分方程式（1-41）]下的微分方程定解问题[7]。为区域边界微元处的外法向量；Г为计算区域Ω的边界。注意：在边界S₁₁、S₁₂上有，限制权函数在边界S₂₁、S₂₂及S₂₀上取值为零，即W_i=0，故式（1-45）中的边界积分项为0。令方程式（1-45）对选取的有限维函数空间中任意的权函数或检验函数成立，则可以建立相应的有限维代数方程。

4.单元分析及整体合成

整个计算域Ω上的积分式（1-45）可表示为各网格单元积分加和的形式，即

任意单元e上的积分对加权残量的贡献为

所有计算区域中任意空间位置的电位可以通过展开系数（或节点电位值）的插值函数获得，即

式中，n为有限元网格节点的总数；N_j为给定的已知函数，称为节点j的基函数[7]；φ_j是待定展开系数即自由度。假设φ在每个三角单元内线性变化并满足整体C0连续，则线性有限元在基函数下，展开系数φ_j的意义为节点j上的电位函数值。

对于图1-4中所示的任一网格单元e，其节点编号按逆时针顺序依次记为k、m、n。假设与编号节点相对应的电位值分别为φ_k、φ_m、φ_n，则单元任意位置的电位φ（x,y）可通过以下的线性插值函数表示：

式中，a、b和c为待定系数，由节点电位插值条件确定。将节点坐标及对应的电位值代入（1-49）有

解上述方程可以得到系数a、b和c的表达式

式中

p_k=x_my_n-y_mx_n,p_m=x_ny_k-y_nx_k,p_n=x_ky_m-y_kx_m

q_k=y_m-y_n,q_m=y_n-y_k,q_n=y_k-y_m

r_k=x_n-x_m,r_m=x_k-x_n,r_n=x_m-x_k

这里，Δ为三角形单元的有向面积

代入式（1-48）中可得

通常简写为

在应用加权余量法进行数值计算时，需要选择适当的基函数和权函数。若基函数N_i与权函数W_i一致，这种方法称为伽辽金（Galerkin）有限元法，那么式（1-47）中的各导数项为

最后式（1-47）中的残量表达式为

这里，又可称为关于单元e第i个节点的单元余量。式（1-57）写作矩阵形式为

式中，矩阵K_e称作单元刚度矩阵，类似地可以得到其他所有单元的刚度矩阵。

各个单元刚度矩阵中的元素按其所在行列对应节点的全局编号，可以叠加到总体刚度矩阵K中相应的位置，从而形成整体刚度矩阵[7]。具体实现方法举例说明如下：首先将总体矩阵K（维数为n×n）中的所有元素置零，此时

接下来，将1号三角单元（假设三个网格节点按照逆时针顺序的全局编号为k=1，m=2，n=4）的单元刚度矩阵为

将1号三角单元刚度矩阵K₁中各元素按其下标地址累加到总体矩阵K的相应位置，此时有

其他单元的刚度矩阵也类似叠加到总体刚度矩阵中，就可以组装生成整体系数矩阵K（也叫全局刚度矩阵）。注意，由于拉普拉斯方程右端项（即激励项）为0，故形成的有限元线性方程组右端荷载向量为0。

5.边界条件的处理及求解

在形成整体系数矩阵K后，还不能直接求解线性代数方程组，因为尚未处理第一类边界条件。对应第一类边界条件的节点，因为其电位为已知，所以可以移除这些节点对应的线性方程。把包含第一类边界条件的系数项作为已知量移到方程的右端之后，可以得到剩余自由度的方程[7]。另外也可以保留第一类边界条件对应的自由度，将这些自由度编号（设第i个节点位于第一类边界，相应自由度取值为φ_i）对应的系数矩阵的对角线位置设置为一个充分大的实数tanv，同时第i个右端项分量修改为φ_i×tanv，则既保持了原有矩阵K的对称性，同时又方便地解决了第一类边界条件的施加。最终形成的有限元线性代数方程组为

K φ=F

通过对以上线性方程组求解，可获得各节点的电位值。

6.数值计算结果的精度验证

为了验证有限元计算结果的正确性和精度，下面与解析解进行比较。对于该二维静电场模型，其解析解表达式可以参考本章参考文献[16]

式中，r表示圆形电极半径，a=b=1，r=0.2，k=0.08290，（x，y）为求解区域任一点的坐标。如果公式推导及编程正确，则有限元方法的计算精度将随着网格单元数量的增加而提高。当然，网格越密，计算量越大，计算时间越长。本节通过给出四种不同密度的有限元网格来说明网格节点数量对计算精度的影响，所采用的四套网格如图1-5所示。通过计算每套网格下数值解与精确解之间的L2误差，可以发现加密网格对计算精度的提升（见表1-1），并且达到了与理论分析相当的收敛阶[7]。

7.高阶有限元基函数

采用有限元方法求解出标量电位φ之后，实际应用中还需要计算其梯度以得到电场强度。由于在线性有限元空间中求解得到的φ仅具有分片线性逼近精度，因此在任意单元对φ求一阶空间导数将得到分片常数，对的逼近精度存在不足。为了提升对的逼近精度，可以对φ采用二次有限元基函数（三角单元上的六个二次有限元基函数如图1-6所示[10]）进行展开，求出φ之后，空间求导将得到分片线性分布的。另外注意采用拉格朗日（Lagrange）有限元基函数，所得到的解仅满足C0连续性，网格节点处的值还需要额外的光滑处理方法进行计算。