工程电磁场数值分析有限元法

工程电磁场数值分析有限元法第一页，共八十二页，编辑于2023年，星期六第4章电磁场有限元法

（FiniteElementMethod,FEM）有限元法可以基于变分原理导出，也可以基于加权余量法导出，本章以加权余量法作为有限元法的基础，以静电场问题的求解为例介绍有限元法的基本原理与实施步骤。并介绍有限元法中的一些特殊问题。第二页，共八十二页，编辑于2023年，星期六第4章电磁场有限元法（FEM）有限元基本原理与实施步骤：1DFEM有限元基本原理与实施步骤：2DFEM有限元方程组的求解二维有限元工程应用三维有限元原理与工程应用矢量有限元第三页，共八十二页，编辑于2023年，星期六加权余量法回顾： 对算子方程用作为该方程的近似解（试探解）：代入方程得余量：1.有限元法基本原理与实施步骤：一维问题在有限元法中，基函数一般用表示。采用Galerkin方案，取权函数与基函数相同。使与余量正交化：第四页，共八十二页，编辑于2023年，星期六设L为线性算子，代入，得或记得代数方程组：加权余量法回顾（续）第五页，共八十二页，编辑于2023年，星期六利用有限元法求解一维边值问题：（1）单元剖分 如图5个单元，6个节点（2）选取基函数

第六页，共八十二页，编辑于2023年，星期六（3）方程离散（计算系数阵[K]和右端项[b]）

基函数Ni只是一阶可导 的，不能严格满足微分方 程，称为“弱解”。第七页，共八十二页，编辑于2023年，星期六（3）方程离散第一项在xj处为0，在xi处的值被来自(i-1)单元的贡献抵消，故只剩下第二项。由于基函数Ni局域支撑，显见只有不为0。使用分步积分：第八页，共八十二页，编辑于2023年，星期六（3）方程离散故类似，当j=i时右端项：第九页，共八十二页，编辑于2023年，星期六总体方程强加边界条件：u1=0,u6=0第十页，共八十二页，编辑于2023年，星期六（4）求解方程思考：（1）有限元的解跟有限差分法的解有何根本不同？ （2）有限元的系数阵总是对称的吗？x0000.20.03610.03600.40.06280.06250.60.07100.07080.80.05250.05231.000第十一页，共八十二页，编辑于2023年，星期六与有限差分法（FDM）相比，有限差分法是对点的离散，得到一系列离散点上的解；而有限元（FEM）是对区域的离散（单元），尽管所求的是节点上的自由度，但它的解在场域中每一个点上都有定义。所以，即是有限元节点上的解是精确的，有限元的整个解仍然是近似的。好的数据处理技术可以从该近似解中提取更精确的分析结果。线性单元中，如果所求的自由度是电位j，单元中的电场E是场量；节点上的E取邻近单元的平均。一些补充说明：

关于有限元的解第十二页，共八十二页，编辑于2023年，星期六计算系数阵是有限元分析的主要工作量。所涉及到的积分，如果不是解析可积的，通常要用到数值积分。其中最常用的数值积分方法是Gauss数值积分。一些补充说明：

高斯数值积分先将积分区间变换到[-1,1]上；按照固定的积分点计算若干函数值P(xi),以固定权值wi累加即可。具(2n+1)阶精度。第十三页，共八十二页，编辑于2023年，星期六n=4x(1)=0.861136311594053d0x(2)=0.339981043584856d0w(1)=0.347854845137454d0w(2)=0.652145154862546d0n=5x(1)=0.906179845938664d0x(2)=0.538469310105683d0x(3)=0.0d0w(1)=0.236926885056189d0w(2)=0.478628670499366d0w(3)=0.568888888888889d0n=6x(1)=

0.932469514203152d0x(2)=

0.661209386466265d0x(3)=

0.238619186083197d0w(1)=0.171324492379170d0w(2)=0.360761573048139d0w(3)=0.467913934572691d0n=16x(1)=0.9894003948d0x(2)=0.9445750231d0x(3)=0.8656312024d0x(4)=0.7554044084d0x(5)=0.6178762444d0x(6)=0.4580167777d0x(7)=0.2816035508d0x(8)=0.0950125098d0w(1)=0.0271524594d0w(2)=0.0622535239d0w(3)=0.0951585117d0w(4)=0.1246289713d0w(5)=0.1495959888d0w(6)=0.1691565194d0w(7)=0.1826034150d0w(8)=0.1894506105d0一些Gauss积分点和权值：

（关于x=0对称，只给出一半）第十四页，共八十二页，编辑于2023年，星期六为提高有限元分析精度，有两种方法： 其一：增加节点，细化网格——称为h方法。 其二：增加有限元的阶数——称为p方法。一些补充说明：

线性单元与高阶单元第十五页，共八十二页，编辑于2023年，星期六一些补充说明：

二阶单元第十六页，共八十二页，编辑于2023年，星期六一些补充说明：

三阶单元第十七页，共八十二页，编辑于2023年，星期六h方法和p方法的求解精度ByJianmingJin.TheFiniteElementMethodinElectromagnetics,2ndEd.,2002第十八页，共八十二页，编辑于2023年，星期六作业：要独立完成，凡雷同者没分！！ 编写有限元程序，计算一维边值问题。改变剖分单元数目，观察解的精度变化。（建议也同时做一个有限差分法的程序，比较二者的精度差别）第十九页，共八十二页，编辑于2023年，星期六以二维静电场泊松方程的求解为例。2.有限元法基本原理与实施步骤：

二维问题

目标：依据加权余量法，利用分域基，建立离散的代数方程组，即确定系数{Kij}

和{bi}。第二十页，共八十二页，编辑于2023年，星期六场域离散 二维问题常使用三角形单元离散，便于处理复杂的场域形状，容易实现。

单元：互不重叠，覆盖全部场域；每个单元内介质是 单一、均匀的。

节点：网格的交点，待求变量的设置点。 该步骤需要记录的信息： 节点编号、节点坐标 节点属性(激励源、是否边界等) 单元编号 单元节点编号 单元介质 第二十一页，共八十二页，编辑于2023年，星期六基函数

有限元采用分片逼近的思想，类似于一维情况下使用折线逼近一条任意曲线。使用分域基Ni，基函数的个数等于节点的个数；每个基函数Ni的作用区域是与该节点i相关联的所有单元。 第二十二页，共八十二页，编辑于2023年，星期六三角形单元内的基函数 设三角形三个顶点处待求函数值分别为u1,u2,u3。如果单元足够小，可以采用线性近似，将单元内任意p点的u(x,y)表示为 代入三个顶点的坐标和函数值，可以解出a、b、c。得到第二十三页，共八十二页，编辑于2023年，星期六

单元节点的编号按逆时针方向排列！其中，第二十四页，共八十二页，编辑于2023年，星期六记住我们的任务—寻找基函数对比可得基函数Ni常被称为插值函数或者形状函数，具有以下性质：（1）是插值的；（2）（3）在相邻单元的公共边界上，Ni是连续的，从而通过Ni构造的逼近函数也是连续的。第二十五页，共八十二页，编辑于2023年，星期六在积分中，对于确定的i，j的有效取值为i本身以及与节点i相联的周围节点，积分的有效区域为以i、j为公共节点的所有三角形单元，在这些单元中Ni、Nj才有交叠。计算系数阵

第二十六页，共八十二页，编辑于2023年，星期六

这些积分可以分单元进行。例如对右图所示的局部编码，K01、K00以及b0的计算公式为：

计算系数阵

第二十七页，共八十二页，编辑于2023年，星期六以下把单元e的贡献记为这样，就有

每个或的计算都在具体的单元内单独考虑（称为单元分析）。第二十八页，共八十二页，编辑于2023年，星期六单元分析：计算单元内积分对系数阵和右端项元素的贡献。

系数阵元素：当L为拉普拉斯算子时，由于Ni在单元内是(x,y)的线性函数，经Laplace算子作用后值为0。但是，在相邻单元的边界上，Ni是连续但是不光滑的，因此对积分的贡献主要来自边界。为考虑单元边界的影响，需要借助于格林公式：第二十九页，共八十二页，编辑于2023年，星期六故，格林公式：因：第三十页，共八十二页，编辑于2023年，星期六 写成一般形式，若一个三角形三个顶点编号为i,j,m（逆时针顺序），则 从而第三十一页，共八十二页，编辑于2023年，星期六再看边界部分：

（1）在节点i的对边Gjm上，Ni＝0，故积分贡献为0；

结论：单元边界对积分的贡献为0。所以单元e对系数阵元素的贡献为：

（2）在节点i的邻边Gij上，由于计算Kij时需要把具有公共邻边的单元的积分累加，此二单元的Ni是连续的；对于单一均匀媒质，要求相邻单元满足 ，故积分的贡献相互抵消。第三十二页，共八十二页，编辑于2023年，星期六

由于单元很小，做单元分析时通常可以取f(e)为常数值（可以认为等于三个顶点上的平均值）。因此

右端项元素：公式：第三十三页，共八十二页，编辑于2023年，星期六上述以节点为序的分析过程对于有限元原理的说明是易于理解的。而在实际编程中，更有效率的是以单元为序，逐个计算单元系数阵[K(e)]，然后合成整体系数阵[K]。单元系数阵[K(e)]定义为

设i,j,m是节点的整体编号，元素Kij在整体矩阵中的实际位置是第i行、j列；因此必须合成到整体矩阵的第i行、j列元素上。

单元矩阵：第三十四页，共八十二页，编辑于2023年，星期六

整体矩阵合成：第三十五页，共八十二页，编辑于2023年，星期六 通过上述过程，对于一个“正常”的内部节点就建立起了一个代数方程。“非正常”的节点包括：媒质交界面衔接条件和场域边界条件。第三十六页，共八十二页，编辑于2023年，星期六 对于静电场问题，媒质分界面衔接条件为媒质交界面衔接条件

第一个条件是自动满足的（Why？），无须格外处理。

对于第二个条件，前面计算单元边界上积分时，默认两边u

的法向导数相等，使内边界上的积分结果抵消。因此只要把泊松方程写成或 满足的条件将是,从而也无需另行处理。第三十七页，共八十二页，编辑于2023年，星期六 由于有限元方法能够自动满足媒质交界面条件，因此有限元法特别适合于处理多层复杂媒质问题。这是其它方法无可比拟的。媒质交界面衔接条件第三十八页，共八十二页，编辑于2023年，星期六第一类边界条件（强加边界条件）第一类边界节点是指边界上函数值已知。因此处理方法是，合成整体系数阵之后，将该节点所在行的主元素置1，其它元素均置零，同时将右端项中对应元素设为已知函数值。要保持对称性；有更简便的做法第三十九页，共八十二页，编辑于2023年，星期六第二类边界条件（自然边界条件）第二类边界节点是指边界上函数法向导数已知。对于内部单元，相邻单元边界的积分相互抵消。但是对于场域边界，如果给定第二类边界条件不为0，则积分结果要计入右端项中。但是若给定的是齐次第二类边界条件，则积分结果为0，无需另行处理，非常方便。——这是ANSYS中自动满足的边界条件。第四十页，共八十二页，编辑于2023年，星期六有限元方法的推导过程虽然看起来有些复杂，但是最终结果是非常简单而且优美的。因为边界条件的处理和媒质交界面条件的处理都非常方便，使得有限元方法在处理复杂媒质问题和复杂场域问题时得心应手，获得了广泛的应用，成为最重要的数值分析手段，广泛应用于各个领域。有人用“功盖四方”来形容有限元，实不为过。中国人在有限元的发明中有自己独特的贡献。第四十一页，共八十二页，编辑于2023年，星期六作业：(2)对于研究方向为数值计算的同学： 编写一个二维静电场有限元程序，计算右图所示问题，或其它自己找一个问题。（1）推导三角形单元的2次和3次插值函数。第四十二页，共八十二页，编辑于2023年，星期六3.有限元方程组的求解代数方程组求解方法概述

所有的数值方法最终都归结为求解一个代数方程组： 系数阵A也称系统矩阵或刚度矩阵。不同离散方法得到的系统矩阵具有不同的特点，方程组的解法也就不同。

基于微分方程（如FEM、FDM等）得到的系统矩阵是稀疏的，有时还是对称的； 而基于积分方程得到的系统矩阵则是稠密的，如BEM、模拟电荷法等。第四十三页，共八十二页，编辑于2023年，星期六 代数方程组的求解是数值计算（计算数学）研究的核心内容。求解代数方程组的方法归纳起来有两类：直接法和迭代法。 3.有限元方程组的求解直接法：直接法都是基于高斯消去法，经过确定次数的运算，理论上可以得到方程组的精确解。适用于小型、稠密方程组的计算。

迭代法：是一种间接方法，从某个预定的初值出发，按照一定的迭代步骤，逐渐逼近方程组的真解。得到一个满足给定精度要求的近似解。适用于大型、稀疏方程组的计算。第四十四页，共八十二页，编辑于2023年，星期六3.有限元方程组的求解直接法（LU分解算法）第四十五页，共八十二页，编辑于2023年，星期六LU分解算法：回带：消元：

计算量：需要的乘除法次数：O(n3)

稳定性：选主元第四十六页，共八十二页，编辑于2023年，星期六迭代法迭代法的基本思想：（等价方程组）从一组猜测的初值开始迭代直至不再变化为止，即为方程组的解（收敛）。好的迭代法应该：对初值不敏感；收敛速度快。第四十七页，共八十二页，编辑于2023年，星期六例如：高斯—赛德尔迭代（有限差分法常用）第四十八页，共八十二页，编辑于2023年，星期六高斯—赛德尔迭代实际运算过程：。。。。。这就是通常的高斯塞德尔迭代格式，矩阵中的零元素不参与运算，矩阵甚至可以不出现。大大减少了内存需求量和计算量。第四十九页，共八十二页，编辑于2023年，星期六共轭梯度法（ConjugateGradientMethod,CG法） 共轭梯度法在原理上可以通过n步迭代得到方程的准确解，因而也称为半直接法或半迭代法。把迭代法表示为更一般的形式：

称为步长，p称为搜索方向，用r表示残差。第五十页，共八十二页，编辑于2023年，星期六预优共轭梯度法（PreconditionedConjugateGradientMethod,PCG法）当系数阵的特征值较均匀地分布在一个很长的区间上时，称矩阵具有很大的条件数；此时共轭梯度法的收敛速度可能很慢。解决的办法是选取合适的非奇异矩阵C进行处理：矩阵称为预处理矩阵。第五十一页，共八十二页，编辑于2023年，星期六第五十二页，共八十二页，编辑于2023年，星期六 预优矩阵M应具有如下特性：稀疏性与A相近；矩阵的特征值分布集中；形如的方程组容易求解；易于寻找。目前公认有效的方法是对系数阵A做不完全Cholesky分解，以M=LDLT

为预优矩阵。这种方法称为不完全分解预优共轭梯度法（IncompleteCholeskydecompositionpreconditionedConjugateGradientMethod，ICCG法）。ICCG法第五十三页，共八十二页，编辑于2023年，星期六稀疏矩阵技术

没有稀疏矩阵技术就没有有限元的成功。带状矩阵技术

通过合适的节点编码，可使系统矩阵的非零元素集中于主对角线附近的带形区域内。在使用LU分解法求解方程组的过程中，带形区域以外的0元素无需计算。

缺点：节点编码优化；带形区域内仍然有大量零元素第五十四页，共八十二页，编辑于2023年，星期六非零元素存储技术

只存储矩阵的非零元素。一般用一个一维数组存储矩阵的非零元素，用另一个索引数组存储这些非零元素在原矩阵中的行和列值。例如：第五十五页，共八十二页，编辑于2023年，星期六非零元素存储技术

在迭代法中，系统矩阵参与的运算只有矩阵左乘以某个相量。如果采用上面的存储方法，则实现q=Ap的算法如下： CComputethematrix–vectorproduct DOk=1,not i=IA(k) j=JA(k) q(i)=q(i)+A(k)*p(j) ENDDO第五十六页，共八十二页，编辑于2023年，星期六非零元素存储技术

下列的二维存储方案是我的发明：直接将矩阵挤扁。数组A有n行m列，m是A各行非零元素个数的最大值。对角元素放在该行的第一个位置上，且不管是否为0，都要存贮。用一个整型数组JA存放各元素的列号。JA中第一列可以用来记录每行的非零元素个数。第五十七页，共八十二页，编辑于2023年，星期六本节更多的参考文献：金建铭.电磁场有限元方法，西安电子科技大学出版社，1998徐树方，矩阵计算的理论与方法，北京大学出版社，1995杨绍祺，谈根林，稀疏矩阵——算法及其程序实现,高等教育出版社,1985刘万勋,刘长学,华伯浩等，大型稀疏线性方程组的解法.国防工业出版社,1981RP梯华森,稀疏矩阵，科学出版社，1981第五十八页，共八十二页，编辑于2023年，星期六有限元分析精度的影响因素静电场问题静磁场问题涡流问题波的传输与散射4.二维有限元的工程应用第五十九页，共八十二页，编辑于2023年，星期六有限元分析精度的影响因素（1）数学模型对工程问题的近似；（2）材料电磁参数的不确定性；（3）数学模型的有限元近似：场域拟合精度——单元大小、未知数个数与局部 场 的变化情况；边界拟合精度——曲线边界；系数阵计算过程中数值积分精度；方程求解精度、数字误差；计算结果的数据处理；好的处理技术可以提高分析 精度。第六十页，共八十二页，编辑于2023年，星期六h方法、p方法；高次单元单元阶数与计算精度的关系

ByJin。第六十一页，共八十二页，编辑于2023年，星期六曲边三角形单元：更好地拟合曲面边界第六十二页，共八十二页，编辑于2023年，星期六数值积分——三角元的高斯数值积分 单元比较小的时候，单元内的函数可以近似认为是常数，通常可以获得满意的精度。当单元内函数变化比较快，或者采用曲边三角形单元、高次单元时时，都会用到数值积分。

高斯-勒让德数值积分公式：——L1、L2、L3

是位置的面积坐标。——权值wi如下页表格。第六十三页，共八十二页，编辑于2023年，星期六三角形单元高斯勒让德数值积分点和权值。第六十四页，共八十二页，编辑于2023年，星期六三角形单元高斯勒让德数值积分点和权值。第六十五页，共八十二页，编辑于2023年，星期六二维边值问题的通用形式在媒质交界面G12上，第六十六页，共八十二页，编辑于2023年，星期六

单元矩阵元素计算公式：

单元右端项计算公式：二维边值问题一般形式在媒质交界面G12上的条件自动满足；第一类边界条件需要强加；第三类边界条件的计算参看金建铭《电磁场有限元方法》。第六十七页，共八十二页，编辑于2023年，星期六