电磁场有限元分析.ppt

1、工程电磁场数值分析（4）（电磁场有限元法）,华中科技大学电机与控制工程系 陈德智 Email: Tel:Office: Room 108, 电机楼 2010.12,第4章 电磁场有限元法（Finite Element Method, FEM）,有限元法可以基于变分原理导出，也可以基于加权余量法导出，本章以加权余量法作为有限元法的基础，以静电场问题的求解为例介绍有限元法的基本原理与实施步骤。并介绍有限元法中的一些特殊问题。,第4章 电磁场有限元法（FEM）,有限元基本原理与实施步骤：1D FEM 有限元基本原理与实施步骤：2D FEM 有限元方程组的求解 二维有限元工程

2、应用 三维有限元原理与工程应用 矢量有限元,加权余量法回顾： 对算子方程 用 作为该方程的近似解（试探解）： 代入方程得余量：,1. 有限元法基本原理与实施步骤：一维问题,在有限元法中，基函数一般用 表示。采用Galerkin方案，取权函数与基函数相同。使与余量正交化：,设L为线性算子，代入 ，得,或,记,得代数方程组：,加权余量法回顾（续）,利用有限元法求解一维边值问题： （1）单元剖分 如图5个单元，6个节点 （2）选取基函数,（3）方程离散 （计算系数阵 K 和右端项 b） 基函数 Ni 只是一阶可导 的，不能严格满足微分方 程，称为“弱解”。,（3）方程离散,第一项在 xj 处为0，在

3、 xi 处的值 被来自 (i-1) 单元的贡献抵消，故只剩下第二项。,由于基函数 Ni 局域支撑，显见只有 不为0。 使用分步积分：,（3）方程离散,故,类似，当 j = i 时,右端项：,总体方程,强加边界条件：u1 = 0, u6 = 0,（4）求解方程,思考：（1）有限元的解跟有限差分法的解有何根本不同？ （2）有限元的系数阵总是对称的吗？,与有限差分法（FDM）相比，有限差分法是对点的离散，得到一系列离散点上的解；而有限元（FEM）是对区域的离散（单元），尽管所求的是节点上的自由度，但它的解在场域中每一个点上都有定义。 所以，即是有限元节点上的解是精确的，有限元的整个解仍然是近似的。好

4、的数据处理技术可以从该近似解中提取更精确的分析结果。 线性单元中，如果所求的自由度是电位j，单元中的电场 E是场量；节点上的 E 取邻近单元的平均。,一些补充说明：关于有限元的解,计算系数阵是有限元分析的主要工作量。所涉及到的积分，如果不是解析可积的，通常要用到数值积分。其中最常用的数值积分方法是Gauss数值积分。,一些补充说明： 高斯数值积分,先将积分区间变换到-1,1上；按照固定的积分点计算若干函数值 P(xi), 以固定权值 wi 累加即可。具(2n+1)阶精度。,n=4 x(1)= 0.861136311594053d0 x(2)= 0.339981043584856d0 w(1)=

5、 0.347854845137454d0 w(2)= 0.652145154862546d0 n=5 x(1)= 0.906179845938664d0 x(2)= 0.538469310105683d0 x(3)= 0.0d0 w(1)= 0.236926885056189d0 w(2)= 0.478628670499366d0 w(3)= 0.568888888888889d0 n=6 x(1)= 0.932469514203152d0 x(2)= 0.661209386466265d0 x(3)= 0.238619186083197d0 w(1)= 0.171324492379170d0

6、 w(2)= 0.360761573048139d0 w(3)= 0.467913934572691d0,n=16 x(1) = 0.9894003948d0 x(2) = 0.9445750231d0 x(3) = 0.8656312024d0 x(4) = 0.7554044084d0 x(5) = 0.6178762444d0 x(6) = 0.4580167777d0 x(7) = 0.2816035508d0 x(8) = 0.0950125098d0 w(1) = 0.0271524594d0 w(2) = 0.0622535239d0 w(3) = 0.0951585117d0

7、w(4) = 0.1246289713d0 w(5) = 0.1495959888d0 w(6) = 0.1691565194d0 w(7) = 0.1826034150d0 w(8) = 0.1894506105d0,一些Gauss积分点和权值：（关于x=0对称，只给出一半）,为提高有限元分析精度，有两种方法： 其一：增加节点，细化网格称为h方法。 其二：增加有限元的阶数称为p方法。,一些补充说明： 线性单元与高阶单元,一些补充说明： 二阶单元,一些补充说明： 三阶单元,h方法和p方法的求解精度,By Jianming Jin. The Finite Element Method in El

8、ectromagnetics, 2nd Ed., 2002,作业：要独立完成，凡雷同者没分！,编写有限元程序，计算一维边值问题。改变剖分单元数目，观察解的精度变化。（建议也同时做一个有限差分法的程序，比较二者的精度差别）,以二维静电场泊松方程的求解为例。,2. 有限元法基本原理与实施步骤：二维问题,目标：依据加权余量法，利用分域基，建立离散的代数方程组，即确定系数Kij 和bi。,场域离散 二维问题常使用三角形单元离散，便于处理复杂的场域形状，容易实现。,单元：互不重叠，覆盖全部场域；每个单元内介质是 单一、均匀的。 节点：网格的交点，待求变量的设置点。,该步骤需要记录的信息： 节点编号、节点

9、坐标 节点属性(激励源、是否边界等) 单元编号 单元节点编号 单元介质,基函数,有限元采用分片逼近的思想，类似于一维情况下使用折线逼近一条任意曲线。 使用分域基Ni，基函数的个数等于节点的个数；每个基函数Ni的作用区域是与该节点i相关联的所有单元。,三角形单元内的基函数 设三角形三个顶点处待求函数值分别为u1, u2, u3。如果单元足够小，可以采用线性近似，将单元内任意p点的u(x,y)表示为,代入三个顶点的坐标和函数值，可以解出a、b、c。得到,单元节点的编号按逆时针方向排列！,其中，,记住我们的任务 寻找基函数,对比,可得,基函数Ni常被称为插值函数或者形状函数，具有以下性质： （1）是

10、插值的； （2） （3）在相邻单元的公共边界上， Ni是连续的，从而通过Ni构造的逼近函数也是连续的。,在积分 中，对于确定的 i，j的有效取值为i 本身以及与节点i相联的周围节点，积分的有效区域为以i、j 为公共节点的所有三角形单元 ，在这些单元中Ni、Nj才有交叠。,计算系数阵,这些积分可以分单元进行。例如对右图所示的局部编码，K01、K00以及b0的计算公式为：,计算系数阵,以下把单元e的贡献记为,这样，就有,每个 或 的计算都在具体的单元内单独考虑（称为单元分析）。,单元分析：计算单元内积分对系数阵和右端项元素的贡献。,系数阵元素：,当L为拉普拉斯算子时，由于Ni在单元内是(x, y)

11、的线性函数，经Laplace算子作用后值为0。 但是，在相邻单元的边界上， Ni是连续但是不光滑的，因此对积分的贡献主要来自边界。 为考虑单元边界的影响，需要借助于格林公式：,故 ，,格林公式：,因：,写成一般形式，若一个三角形三个顶点编号为i, j, m（逆时针顺序），则,从而,再看边界部分：,（1）在节点 i 的对边Gjm上，Ni0，故积分贡献为0；,结论：单元边界对积分的贡献为0。所以单元e对系数阵元素的贡献为：,（2）在节点 i 的邻边Gij上，由于计算Kij时需要把具有公共邻边的单元的积分累加，此二单元的Ni是连续的；对于单一均匀媒质，要求相邻单元满足 ，故积分的贡献相互抵消。,由于

12、单元很小，做单元分析时通常可以取 f (e) 为常数值（可以认为等于三个顶点上的平均值）。因此,右端项元素：,公式：,上述以节点为序的分析过程对于有限元原理的说明是易于理解的。而在实际编程中，更有效率的是以单元为序，逐个计算单元系数阵K(e)，然后合成整体系数阵K。单元系数阵K(e)定义为 设 i, j, m 是节点的整体编号，元素Kij在整体矩阵中的实际位置是第i行、j列；因此 必须合成到整体矩阵的第i行、j列元素上。,单元矩阵：,整体矩阵合成：,通过上述过程，对于一个“正常”的内部节点就建立起了一个代数方程。“非正常”的节点包括：媒质交界面衔接条件和场域边界条件。,对于静电场问题，媒质分界

13、面衔接条件为,媒质交界面衔接条件,第一个条件是自动满足的（Why？），无须格外处理。,对于第二个条件，前面计算单元边界上积分 时，默认两边 u 的法向导数相等，使内边界上的积分结果抵消。因此只要把泊松方程写成 或 满足的条件将是 , 从而也无需另行处理。,由于有限元方法能够自动满足媒质交界面条件，因此有限元法特别适合于处理多层复杂媒质问题。这是其它方法无可比拟的。,媒质交界面衔接条件,第一类边界条件（强加边界条件）,第一类边界节点是指边界上函数值 已知。因此处理方法是，合成整体系数阵之后，将该节点所在行的主元素置1，其它元素均置零，同时将右端项中对应元素设为已知函数值。,要保持对称性；有更简便

14、的做法,第二类边界条件（自然边界条件）,第二类边界节点是指边界上函数法向导数 已知。对于内部单元，相邻单元边界的积分 相互抵消。但是对于场域边界，如果给定第二类边界条件不为0，则积分结果要计入右端项中。但是若给定的是齐次第二类边界条件，则积分结果为0，无需另行处理，非常方便。,这是ANSYS中自动满足的边界条件。,有限元方法的推导过程虽然看起来有些复杂，但是最终结果是非常简单而且优美的。因为边界条件的处理和媒质交界面条件的处理都非常方便，使得有限元方法在处理复杂媒质问题和复杂场域问题时得心应手，获得了广泛的应用，成为最重要的数值分析手段，广泛应用于各个领域。有人用“功盖四方”来形容有限元，实不

15、为过。 中国人在有限元的发明中有自己独特的贡献。,作业：,(2) 对于研究方向为数值计算的同学： 编写一个二维静电场有限元程序，计算右图所示问题，或其它自己找一个问题。,（1）推导三角形单元的2次和3次插值函数。,3. 有限元方程组的求解,代数方程组求解方法概述 所有的数值方法最终都归结为求解一个代数方程组：,系数阵 A也称系统矩阵或刚度矩阵。不同离散方法得到的系统矩阵具有不同的特点，方程组的解法也就不同。 基于微分方程（如FEM、FDM等）得到的系统矩阵是稀疏的，有时还是对称的； 而基于积分方程得到的系统矩阵则是稠密的，如BEM、模拟电荷法等。,代数方程组的求解是数值计算（计算数学）研究的核

16、心内容。求解代数方程组的方法归纳起来有两类：直接法和迭代法。,3. 有限元方程组的求解,直接法：直接法都是基于高斯消去法，经过确定次数的运算，理论上可以得到方程组的精确解。适用于小型、稠密方程组的计算。,迭代法：是一种间接方法，从某个预定的初值出发，按照一定的迭代步骤，逐渐逼近方程组的真解。得到一个满足给定精度要求的近似解。适用于大型、稀疏方程组的计算。,3. 有限元方程组的求解,直接法（LU分解算法）,LU分解算法：,回带：,消元：,计算量：需要的乘除法次数： O(n3) 稳定性：选主元,迭代法,迭代法的基本思想：,（等价方程组）,从一组猜测的初值 开始迭代,直至 不再变化为止，即为方程组的

17、解（收敛）。,好的迭代法应该：对初值不敏感；收敛速度快。,例如：高斯赛德尔迭代（有限差分法常用）,高斯赛德尔迭代,实际运算过程：,。,这就是通常的高斯塞德尔迭代格式，矩阵中的零元素不参与运算，矩阵甚至可以不出现。大大减少了内存需求量和计算量。,共轭梯度法（Conjugate Gradient Method, CG法） 共轭梯度法在原理上可以通过n步迭代得到方程的准确解，因而也称为半直接法或半迭代法。,把迭代法表示为更一般的形式：,称为步长，p 称为搜索方向，用r 表示残差。,预优共轭梯度法（Preconditioned Conjugate Gradient Method, PCG法）,当系数阵

18、的特征值较均匀地分布在一个很长的区间上时，称矩阵具有很大的条件数；此时共轭梯度法的收敛速度可能很慢。解决的办法是选取合适的非奇异矩阵C进行处理：,矩阵 称为预处理矩阵。,预优矩阵M应具有如下特性：稀疏性与A相近；矩阵 的特征值分布集中；形如 的方程组容易求解；易于寻找。目前公认有效的方法是对系数阵A做不完全Cholesky分解，以M = LDLT 为预优矩阵。这种方法称为不完全分解预优共轭梯度法（Incomplete Cholesky decomposition preconditioned Conjugate Gradient Method，ICCG法）。,ICCG法,稀疏矩阵技术 没有稀疏

19、矩阵技术就没有有限元的成功。,带状矩阵技术 通过合适的节点编码，可使系统矩阵的非零元素集中于主对角线附近的带形区域内。在使用LU分解法求解方程组的过程中，带形区域以外的0元素无需计算。,缺点：节点编码优化；带形区域内仍然有大量零元素,非零元素存储技术 只存储矩阵的非零元素。一般用一个一维数组存储矩阵的非零元素，用另一个索引数组存储这些非零元素在原矩阵中的行和列值。例如：,非零元素存储技术 在迭代法中，系统矩阵参与的运算只有矩阵左乘以某个相量。如果采用上面的存储方法，则实现 q = A p 的算法如下： C Compute the matrixvector product DO k=1, not

20、 i=IA(k) j=JA(k) q(i)=q(i)+A(k)*p(j) ENDDO,非零元素存储技术 下列的二维存储方案是我的发明：直接将矩阵挤扁。,数组A有n行m列，m是A各行非零元素个数的最大值。对角元素放在该行的第一个位置上，且不管是否为0，都要存贮。用一个整型数组JA存放各元素的列号。JA中第一列可以用来记录每行的非零元素个数。,本节更多的参考文献： 金建铭. 电磁场有限元方法，西安电子科技大学出版社，1998 徐树方，矩阵计算的理论与方法，北京大学出版社，1995 杨绍祺，谈根林，稀疏矩阵算法及其程序实 现 ,高等教育出版社, 1985 刘万勋, 刘长学, 华伯浩等，大型稀疏线性方

21、程组的解法. 国防工业出版社, 1981 R P 梯华森, 稀疏矩阵，科学出版社，1981,有限元分析精度的影响因素 静电场问题 静磁场问题 涡流问题 波的传输与散射,4. 二维有限元的工程应用,有限元分析精度的影响因素,（1）数学模型对工程问题的近似； （2）材料电磁参数的不确定性； （3）数学模型的有限元近似： 场域拟合精度单元大小、未知数个数与局部场的变化情况； 边界拟合精度曲线边界； 系数阵计算过程中数值积分精度； 方程求解精度、数字误差； 计算结果的数据处理；好的处理技术可以提高分析精度。,h方法、p方法；高次单元,单元阶数与计算精度的关系By Jin。,曲边三角形单元：更好地拟合曲

22、面边界,数值积分 三角元的高斯数值积分 单元比较小的时候，单元内的函数可以近似认为是常数，通常可以获得满意的精度。当单元内函数变化比较快，或者采用曲边三角形单元、高次单元时时，都会用到数值积分。 高斯-勒让德数值积分公式：, L1、L2、L3 是位置的面积坐标。 权值 wi 如下页表格。,三角形单元高斯勒让德数值积分点和权值。,三角形单元高斯勒让德数值积分点和权值。,二维边值问题的通用形式,在媒质交界面G12 上，,单元矩阵元素计算公式：,单元右端项计算公式：,二维边值问题一般形式,在媒质交界面G12 上的条件自动满足； 第一类边界条件需要强加； 第三类边界条件的计算参看金建铭电磁场有限元方法

23、。,第三类边界条件中，如果g = 0则成为第二类边界条件；如果再有 q = 0 则自动满足，无需任何处理。,二维静电场：,无限远边界条件,静电场分析,电力线平行与垂直边界条件,注意与下式的区别：,轴对称静电场：,无限远边界条件,静电场分析,如果不将 r 看作系数，有 系数阵将不再对称。,二维静磁场：,无限远边界条件,静磁场分析,磁力线垂直与平行边界条件,轴对称静磁场：,静磁场分析,在平行平面场中，等 A 线平行于磁力线； 在轴对称静磁场中，等 rA 线平行于磁力线。,注意：与轴对称静电 场的表达式非常不同,如果直接以 A 为求解量， 系数阵将不再对称。,故：,轴对称静磁场问题无限远边界条件：,

24、小电流环产生的磁位：( ),一阶渐近边界条件，适合于电流回路，球心位于回路中心的球面。,其中，,二维静磁场分析（包括轴对称静磁场分析）一般不使用磁标位，因为正如所看到的，磁矢位 A 只有一个分量，而且激励电流和边界条件都很容易处理。 使用磁标位时必须单独处理电流，有时候还需要规定磁障碍面，对有限元来说很不方便。 但是三维的情况将非常不同，此为后话。,二维涡流问题：,涡流场分析,相量的含义； 参数是位置的函数； 直流与交流的区别； 透入深度与单元尺寸；,电场与涡流的计算； 电压激励线圈的处理； 载流块导体的处理；,轴对称涡流问题：,涡流场分析,在不引起混淆的情况下，可以省略表示相量的点。 线圈电

25、压与阻抗的计算。,另一种形式：,涡流场分析中的无限远边界条件,如果所有的电流（包括源电流和涡流）都位于有限区域内，有限元模型的外围区域为空气，那么涡流问题的无限远边界条件与恒定磁场的无限远边界条件一致。 如果导体区域一直延伸到无限远处，有限元模型的外围区域包含导体，普通的无限远边界条件是无效的，需要重新推导。,轴对称涡流场中线圈电压与阻抗的计算 （1） （2）,二维波导问题：,波的传输与散射,二维波散射问题：,波的传输与散射,散射场一阶吸收边界条件（ABC）,总场等于入射场加散射场。,二维有限元的工程问题,非线性问题 电磁力和力矩的计算 永磁问题 功率与阻抗 ANSYS中的二维电磁场分析,非线

26、性问题的求解,设求解静磁场问题,如果磁导率 是常数，使用有限元法可将上述问题化为线性代数方程组求解：,但如果磁导率 是磁感应强度B的函数，从而也是磁矢位A的函数。当A还没获得之前， 是未知的，从而矩阵 K 也未知。 非线性问题需要用迭代法求解。,非线性问题的求解,使用有限元法求解非线性问题的思路是：,k=0, 猜测一个,计算系数阵,求解有限元方程,得到各点的B(k+1),k=k+1， 估计新的,以上是线性迭代的流程，收敛缓慢。更好的方法是牛顿-拉普逊方法。,非线性问题的求解,使用有限元法求解非线性问题的思路是：,k=0, 猜测一个,计算系数阵,求解有限元方程,得到各点的B(k+1),k=k+1

27、， 估计新的,以上是线性迭代的流程，收敛缓慢。更好的方法是牛顿-拉夫逊方法。,求解非线性方程 ： 从一个初值 x0 出发，求得近似方程 的解 。 依次迭代，逐渐逼近问题的真解 x*。,牛顿拉夫逊法（Newton-Raphson method ）,左端记为 ，是一个向量函数； 如果K 不是 u 的函数，则 导数矩阵 ，称为雅克比矩阵，记做 J。 方程组的解一次即可获得 如果K 是 u 的函数，则雅克比矩阵的元素为 迭代求解方程组,非线性方程组 的牛顿拉夫逊算法,牛顿-拉夫逊法对初值很敏感。有一些改进的算法，但是目前没有从根本上解决问题的路子。 使用ANSYS计算铁磁材料磁场时，如果磁化曲线和初始

28、点不合适，也会碰到类似问题。,牛顿拉夫逊法不收敛的例子,B-H曲线的表示,为了计算非线性有限元矩阵，必须输入铁材料的磁化曲线。 磁化曲线可以通过若干点插值获得（如ANSYS）；也可以采用多项式或者其他方式（如指数形式）的拟合表示。 磁化曲线通常不可能很准 确，计算时如果需要可以 根据实际情况作适当的调 整，以保障非线性迭代的 收敛和计算的精度。,正弦涡流问题B-H曲线的等效处理,严格说来，含有非线性材料的场不可能是正弦电磁场。 当只关心力、力矩、功率、损耗等随时间平均意义的量，而不怎么关心波形的时候，正弦近似有很好的精度。 否则只能按照瞬态场计算（计算量很大）。 正弦近似B-H曲线的等效处理方

29、法： 令等效磁导率得到的磁能密度与瞬时形式磁能在一个周期内的平均值相等： 得到,He为磁场有效值,可用数值积分获得。,作业：计算等效磁化曲线,H (Am-1) B(T) 0 0 90.0000000 0.50000000 270.000000 1.00000000 318.250000 1.10000000 384.500000 1.20000000 479.500000 1.30000000 608.562000 1.38750000 755.437000 1.45000000 939.185000 1.50000000 1188.93000 1.54500000 1407.93000 1.57500000 2077.31000 1.62750000 3117.93000 1.67375000 3969.37000 1.70225000,H (A m-1) B(T) 4843.66000 1.72750000 6081.34000 1.75825000 8581.09000 1.80875000 11066.4000 1.850