1、工程电磁场数值分析（有限元法）华中科技大学电机与控制工程系陈德智2007.12第4章 电磁场有限元法（Finite Element Method,FEM）有限元法可以基于变分原理导出，也可以基于加权余量法导出，本章以加权余量法作为有限元法的基础，以静电场问题的求解为例介绍有限元法的基本原理与实施步骤。并介绍有限元法中的一些特殊问题。第4章 电磁场有限元法（FEM）有限元的基本原理与实施步骤有限元方程组的求解前处理与后处理技术渐近边界条件矢量有限元法求解运动导体涡流问题的迎风有限元法 在有限元法中，基函数一般用 表示。采用Galerkin方案，取权函数与基函数相同。使与余量正交化：加权余量法回顾
2、：对算子方程用 作为该方程的近似解（试探解）：代入方程得余量：1.有限元法的基本原理与实施步骤()L uf1niiiu u()RL uf(,)()d0iiN RN L uf(1,2,)in,1,2,iNin设L为线性算子，代入 ，得11()d()d0nnijjijjjjN LNfNL Nf 1niiiuN或(1,2,)in1()ddnjijijN L NN f 记diibN f,()di jijKN L N得代数方程组：Kb加权余量法回顾（续）(,)()d0iiN RN L uf 场域离散以二维静电场泊松方程的求解为例。二维问题常使用三角形单元离散，便于处理复杂的场域形状，容易实现。单元：互不
3、重叠，覆盖全部场域；每个单元内介质是 单一、均匀的。节点：网格的交点，待求变量的设置点。需要记录信息：节点编号、节点坐标节点属性(激励源、是否边界等)单元编号单元节点编号单元介质目标：建立节点变量之间满足的代数方程组，即确定系数Kij 和bi。依据的原理是加权余量法使用的基函数为分域基。基函数 有限元采用分片逼近的思想，跟使用折线逼近一条任意曲线的做法相同。使用分域基Ni，基函数的个数等于节点的个数；每个基函数Ni的作用区域是与该节点i相关联的所有单元。diibN f,()di jijKN L N在积分 中，对于确定的 i，j的有效取值为i本身以及与节点i相联的周围节点，积分的有效区域为以i、
4、j为公共节点的所有三角形单元，在这些单元中Ni、Nj才有交叠。()dijijKN L N这些积分可以分单元进行。例如对右图所示的局部编码，K01、K00以及b0的计算公式为：1234560000()dKN L N 160101()dKN L N 12345600dbN f 以下把单元e的贡献记为()()()()deeeeijijKNL N()()()deeeeiibNf这样，就有(1)(2)(3)(4)(5)(6)00000000000000KKKKKKK(1)(6)010101KKK(1)(2)(3)(4)(5)(6)0000000bbbbbbb 每个 或 的计算都在具体的单元内单独考虑（称
5、为单元分析）。()eijK()eib三角形单元内的基函数设三角形三个顶点处待求函数值分别为u1,u2,u3。如果单元足够小，可以采用线性近似，将单元内任意p点的u(x,y)表示为(,)u x yabxcy312123(,)(,)(,)(,)x yx yx yu x yuuu代入三个顶点的坐标和函数值，可以解出a、b、c。得到312123(,)(,)(,)(,)x yx yx yu x yuuu12312311112xxxyyy 1232311112xxxyyy 2131311112xxxyyy 3121211112xxxyyy 单元节点的编号按逆时针方向排列！其中，112233(,)u x y
6、NNN记住我们的任务寻找基函数对比312123(,)(,)(,)(,)x yx yx yu x yuuu(1,2,3)i 可得(,)iix yN基函数Ni常被称为插值函数或者形状函数，具有以下性质：（1）是插值的；（2）（3）在相邻单元的公共边界上，Ni是连续的，从而通过Ni构造的逼近函数也是连续的。1 ()(,)0 ()ijjijN xyij单元分析：计算单元内积分对系数阵和右端项元素的贡献。()()()()deeeeijijKNL N系数阵元素：当L为拉普拉斯算子时，由于Ni在单元内是(x,y)的线性函数，经Laplace算子作用后值为0。但是，在相邻单元的边界上，Ni是连续但是不光滑的，
7、因此对积分的贡献主要来自边界。为考虑单元边界的影响，需要借助于格林公式：故 ，()2()d dd ddeeejeijijijiNKNNx yNNx yNn 2()d dVSVS格林公式：1233223321()()()2x yx yyy xxxy 23321()()()2yyxxNij31132()()()2yyxxNij23313213122()()()()()()4yyyyxxxxNN(,)iix yN因：写成一般形式，若一个三角形三个顶点编号为i,j,m（逆时针顺序），则13231323122()()()()()()4yyyyxxxxNN d d()()()()4eijimjmimjmN
8、Nx yyyyyxxxx从而2()()()()()()4ijimjmimjmNNyyyyxxxx dejiNNn再看边界部分：（1）在节点 i 的对边jm上，Ni0，故积分贡献为0；()()()()()d d4eimjmimjmeijijyyyyxxxxKNNx y 结论：单元边界对积分的贡献为0。所以单元e为系数阵元素的贡献为：（2）在节点 i 的邻边ij上，由于计算Kij时需要把具有公共邻边的单元的积分累加，此二单元的Ni是连续的；对于单一均匀媒质，要求相邻单元满足 ，故积分的贡献相互抵消。12()()/eeijijnn 由于单元很小，做单元分析时通常可以取 f(e)为常数值（可以认为等于
9、三个顶点上的平均值）。因此()()()deeeeiibNf右端项元素：123!()()()d2(2)!elmnl m nNNNlmn 公式：()()()()d3eeeeeiibfNf 通过上述过程，对于一个“正常”的内部节点就建立起了一个代数方程。“非正常”的节点包括：媒质交界面衔接条件和场域边界条件，稍后再讨论。上述以节点为序的分析过程对于有限元原理的说明是易于理解的。而在实际编程中，更有效率的是以单元为序，逐个计算单元系数阵K(e)，然后合成整体系数阵K。单元系数阵K(e)定义为设i,j,k是节点的整体编号，元素Kij在整体矩阵中的实际位置是第i行、j列；因此 必须合成到整体矩阵的第i行、
10、j列元素上。()()()()()()()()()()eeeiiijikeeeejijjjkeeekikjkkKKKKKKKKKK()eijK对于静电场问题，媒质分界面衔接条件为媒质交界面衔接条件121212nndejiNNn第一个条件是自动满足的（Why？），无须格外处理。20/r 0/r 1212rrnn对于第二个条件，前面计算单元边界上积分 时，默认两边 的法向导数相等，使内边界上的积分结果抵消。因此只要把泊松方程写成 或满足的条件将是 ,从而也无需另行处理。由于有限元方法能够自动满足媒质交界面条件，因此有限元法特别适合于处理多层复杂媒质问题。这是其它方法无可比拟的。媒质交界面衔接条件第一
11、类边界条件（强加边界条件）第一类边界节点是指边界上函数值 已知。因此处理方法是，合成整体系数阵之后，将该节点所在行的主元素置1，其它元素均置零，同时将右端项中对应元素设为已知函数值。iif要保持对称性；有更简便的做法第二类边界条件（自然边界条件）第二类边界节点是指边界上函数法向导数 已知。对于内部单元，相邻单元边界的积分 相互抵消。但是对于场域边界，如果给定第二类边界条件不为0，则积分结果要计入右端项中。但是若给定的是齐次第二类边界条件，则积分结果为0，无需另行处理，非常方便。/jjng dejiNNn有限元方法的推导过程虽然看起来有些复杂，但是最终结果是非常简单而且优美的。因为边界条件的处理
12、和媒质交界面条件的处理都非常方便，使得有限元方法在处理复杂媒质问题和复杂场域问题时得心应手，获得了广泛的应用，称为最重要的数值分析手段。有人用“功盖四方”来形容有限元，实不为过。中国人在有限元的发明中有自己独特的贡献。作业：（1）研究方向为数值计算的同学：编写一个二维静电场有限元程序，计算右图所示问题，或其它自己找一个问题。（2）研究方向非数值计算的同学：简要叙述有限元的原理，试分析计算精度可能跟哪些因素有关；并归纳一下，有限元法与有限差分法有那些相同点和不同点？稀疏矩阵技术ICCG法2.有限元方程组的求解建模自动剖分技术误差估计，h方法与p方法可视化问题：等位线与电力线电场力的计算电容、电感与电阻3.有限元的前处理与后处理技术场域的封闭渐近边界条件4.渐近边界条件节点元存在的问题矢量有限元5.矢量有限元方法速度效应产生的问题迎风法6.运动导体的涡流问题（迎风有限元）
