电磁场数值计算4-西安交通大学电气工程学院

1、第四章 二维电磁场有限元分析二维电磁场：平行平面场、轴对称场的基本概念§4-1 电磁场的微分方程及泛函常用的泊松方程、拉普拉斯方程和亥姆霍兹方程是方程 （4-1）的几种特殊形式，它适用于非线性、非均匀、各向异性材料的静、稳态和简谐电磁场问题。在扩散方程中，在波动方程中。边界条件 电场 （4-2a） 磁场（） （4-2b） E线与对称面平行，或B线与对称面垂直 E线与对称面垂直，或B线与对称面平行式（4-1）对应的泛函求极值问题为 （4-3）（1） 若存在第一类边界条件，需要专门处理。（2） 若存在媒质分界面，变分中将自动满足（3） 不论是实数还是复数，式（4-3）均成立。如果这些参数

2、只是实数则可用另一泛函 （4-4）若为电场问题： （4-5）若为磁场问题：或 （4-6）若为非线性问题，§4-2 有限元分析4.2.1 离散化（前处理）将面区域离散为许多小面单元：三角形（三节点、六节点）、四边形（四节点、八节点）1、 离散的基本要求* 单元之间既没有间隔，又没有重叠，单元通过它们的顶点相连；* 尽量避免较小内角的单元产生，可以证明，有限元解的误差反比与最小内角的正弦。最理想的三角形是等边三角形；* 一个单元的3个节点不能同时在一个边界上，若64123是第一类边界，则单元的剖分法是错误的，这将导致24边等位；* 单元越小，数值解结果越好，但未知量多会增加计算机内存和计

3、算量，会带来计算误差。因此，解变化剧烈的地方剖分可以较密些。2、 编码要求为了压缩内存和简化程序，单元和节点编码应该注意：* 单元用一组整数编码（）；* 节点编码用另一组整数，它反映两个信息：局部编码，表示节点在相应单元中的位置；全局编码，表示该节点在全域中的位置。两个功能用1个二维数组来表示，如，其中，i =1,2,3反映三角单元3个顶点局部位置，是单元编码，的值是节点的全局编码。例：写出4单元6节点的节点编码注：局部编码要逆时针进行。 * 节点坐标数组：* 每个单元的材料参数及源密度4.2.2 单元插值及基函数以三角形3节点单元为例，采用线性插值函数逼近e单元内的未知函数 （4-7）为待定

4、常数，e单元上3个节点的未知数分别为 （4-8）从中解出，再代入式（4-7），得到 （4-9）式中， （4-10）其中， （4-11）单元面积 由此可见，基函数是坐标的线性插值函数，只与三角形的形状及节点分布有关，故又称之为形状函数。基函数具有性质：（1） 定义域上任意点p的坐标代入基函数，有体现方程具有稀疏矩阵的特性。（2） 若p点是单元节点，则，基函数为表明： 场点在节点i上，； 当场点在节点i的对边上时，。所以，一个单元边上的值与它的相对节点处的值无关，它仅由该边两端点处的值决定。这就保证了单元两侧解的连续性。（实际上，是面积坐标一般而言，离P点最近的节点的电位值所占比例最大。与式（4-

5、10）对照可知。当P点落在j、m连线上，所以，。当P点落在i节点上，）这样，每个三角形单元的由该单元上各顶点的值的线性插值去逼近。另外，由于两个相邻单元公共边和公共节点上的值相等，就有可能保持线性插值函数的连续性。若为电场，根据式（4-5），电场强度的数值解为 （4-12）显然，单元内的电场强度是一常矢量，为了提高精度，常采用高阶插值。若为磁场，根据式（4-6），磁场强度的数值解为 （4-13）同理，单元内的磁感应场强度是一常矢量。4.2.3 建立有限元方程不论是基于变分原理的有限元法，还是伽辽金有限元法，最终的有限元方程是一样的其中，系数矩阵中的元素及右端列向量中的元素分别为 （4-14）

6、（4-15）设单元中媒质参数不变，为（若不是常数，可取相应的单元的平均值），通过三角形上的数值积分公式 （4-16）其中，a,b,c分别为1，1，0。若上式积分等于根据基函数的性质，即积分单元区域中，只有与节点i，或节点j关联的单元参加积分，其余为零，可以得到式（4-14）、式（4-15）的解析表达式 （4-17） （4-18）矩阵形式为（电场：，且） （4-19） （4-20）以下图为例，K为66阶方阵 与节点1，1关联的单元为（1），与节点1，2关联的单元为（1），与节点1，3关联的单元为0，与节点1，4关联的单元为（1），与节点1，5关联的单元为0，与节点1，6关联的单元为0， 与节点2

7、，1关联的单元为（1），与节点2，2关联的单元为（1）、（2）、（3），与节点2，3关联的单元为（3），与节点2，4关联的单元为（1）、（2），与节点2，5关联的单元为（2）、（3），与节点2，6关联的单元为0，以此类推，最终有 - D - 系数矩阵K的4个特点：（1） K是对称阵， （2） 对角线元素占优，即，且为正定阵；（3） K是稀疏矩阵，第i行的n个元素中，只有与节点i关联的节点上是非零元素，其余均为零，所以存储量小；（4） K中非零元素的位置呈现带状，带宽为D=1+H , H指三角形单元节点编码最大差值的绝对值，如图中单元或中节点的最大差值为H = 4-1=3或H = 5-2=3，带

8、宽D=4。不同的求解方法有不同的存储技术，大体分为：（1） 等带宽存储：下三角带宽内的各元素依次（按行或按列）存储；（2） 变带宽存储：下三角带宽内从每一行第一个非零元素开始依次存储；（3） 非零元存储：下三角带宽内按行只存储非零元素。前两者用于直接解法，后者用于迭代解法，如ICCG法。4.2.4 边界条件的处理1、第一类边界条件处理设节点3，5，6位于第一类边界上，其值分别为，首先强加做法一：令，这样将导致系数矩阵不对称。做法二：令,令，并以此类推，方程为显然，方程可以降阶，且依然对称2、第三类边界条件处理第二类边界条件是第三类边界条件的特殊情况，故只讨论后者。若为齐次边界条件，不需要增加任

9、何项。若为非齐次边界条件，需要增加一项 （4-21）设边界离散后的单元边有个，上式近似表示为每一线段积分的叠加 （4-22）每一线段上的未知函数可以近似为 （4-23）式中，其中是沿线段从节点1到节点2的归一化距离，即节点1处，节点2处，两个节点之间按线性规律变化，是线段长度。将式（4-22）代入式（4-21），并对求导，得 （4-24）写成矩阵形式 （4-25）式中系数阵是对称的 （4-26） (4-27)如果在线段上是常数，分别用表示，则可用解析解表示 （4-28） （4-29）式中，将式（4-27）与式（4-17）合并，式（4-28）与式（4-18）合并，得到最终方程 （4-30）

10、67;4-3 轴对称场电场与磁场的轴对称问题处理方法略有不同，分别介绍。4.3.1 电场中的轴对称问题设Z轴为对称轴，泊松方程在圆柱坐标下的表达式为等式两边同乘以半径，便可得到 （4-31）做下面替换 （4-32）便可得到与平行平面场中泊松方程形式完全相同的方程。也就是说，前面建立的有限元方程可以直接应用于该问题的求解。轴对称情况下的开域问题可用点源公式设置人工边界 （如 ， ） (4-33)式中，。设的源仅限于附近，当时，的渐进式为 （4-34）与式（4-14）、式（4-15）及式（4-17）、式（4-18）比较有是三角形的3个顶点的半径。等中的替代。电场强度(三角形3节点) （4-35）（

11、第七次课）4.3.2 磁场中的轴对称问题在磁场中，磁矢位满足的是矢量泊松方程，轴对称场中，电流只有方向的分量，那么，考虑非线性磁介质，泊松方程为 （4-36）变形后 （4-37）做下面替换 （4-38）电场中的边界条件与平行平面场的表达形式相同，不再重述。下面讨论作为变量的边界条件是否相应于二维场讨论的一般形式1、第一类边界条件可以证明，轴对称场中，磁感应线方程为 ，因此第一类边界条件 （4-39）2、 第二类边界条件若已知边界上切向磁密（或，即已知面电流K），第二类边界条件是非齐次的，可以写为 （4-40）由于 其中 所以得到（是roz平面边界上的法向分量） 由于 又因 ， 所以 所以，边界

12、条件为 或 （4-41）显然，这是第二类边界条件。若为齐次边界条件3、 分界面衔接条件分界面连续边界条件 ，用磁矢位表示为 ， （4-42）这样，以为变量的所有边界条件和连续性衔接条件均与的条件相对应，有限元方程自动满足分界面衔接条件。单元磁感应强度 （4-43）实际单元磁感应强度 （4-44）磁感应线就是等线 （4-45）节点的磁矢位 （4-46）§4-4 非线性方程组解法磁场多为非线性问题，材料参数是场量的函数，而，因此，有限元方程是非线性方程。需要采用专门的解法。4.4.1 非线性材料特性的数值逼近方法基本磁化曲线（分析恒定磁场）、磁滞回线（分析磁滞损耗）、等效磁化曲线（分析时

13、变场）都是非线性曲线。它表明，如果激励源是正弦电流，那么，通过有限元法计算出来的磁矢位是正弦函数，由此取旋度得到的磁感应强度也是正弦函数，但通过得到的磁场强度一定不是正弦的。如果激励源是正弦电压，那么，由于电感是非常数、非线性的梯度电感，因此电流一定不是正弦的，由电流计算得到的磁矢位、磁感应强度都不是正弦的。在有限元法计算中，需要将B-H曲线用函数形式表达，常用插值法和拟合法，每一种方法又有多种表达形式：插值法：分段线性插值、拉格朗日插值、高次多项式插值、样条插值等；拟合法：最小二乘法、分段线性拟合、相应表面法等。以分段插值为例：先将B-H曲线分成m个子区间，每个子区间两个端点分别为及，用线性

14、插值，在此区间上任一点的为 （4-47）这种插值方法使每一子区间段上的插值函数只依赖于本段上的节点值，与其它节点无关，优点是计算稳定。但由于总体光滑度不高，导致计算精度不高。拉格朗日插值是二阶插值函数，计算精度高于线性插值。样条插值要求节点的一阶导数连续，光滑度高，计算精度也高。4.4.2 线性化迭代算法将非线性方程组分解为一系列线性方程组的迭代格式，然后逐次逼近非线性方程组的解。迭代步骤如下是否k=0设的初值，赋材料参数，可取B-H曲线线性段斜率作为初值建立方程解“线性”代数方程，即用解线性代数方程的方法解之求解查插值曲线，得到，结束检验A是否满足计算精度k=k+1线性化迭代法的特点是：方法

15、简单，有良好的收敛性，但收敛速度较慢，尤其当磁场进入饱和区时，收敛很慢，因此，当方程阶数不是很高时用此方法比较合适。否则用迭代法。4.4.3 牛顿拉夫逊（Newton-Raphson）算法该方法也属于线性化方法，但迭代格式不同，适用于非线性方程组求解。尤其对饱和磁场问题收敛性能好，比线性化迭代法更为优越。1、一个未知量的N-R法先从一个未知量的非线性方程的解法来说明N-R法。设 （4-48）该式的解是图4-1(a)中焦点对应的，设A的近似初解为，用泰勒级图4-1牛顿拉夫逊（Newton-Raphson）算法示意图数在点展开 （4-49）忽略第三项以后各项， （4-50） 是在处的斜率，称为雅可

16、比系数，是常数。移项后有 （4-51）当该式等于零时，得到解。若不为零，求得的解作为第一个近似解，上式方程为以此类推，有通式 （4-52）当满足 时，是最终的近似解。这是一种变步长的方法，利用梯度确定步长,如图4-1(b)所示。这种方法也有不收敛的情况。当是非单调的，可能出现在相邻两次迭代中符号相反的情况，如图4-1(c)所示，不收敛。形成雅可比矩阵很费时，提出一种修正的N-R法。 2、修正的N-R法每次迭代都用，一般的迭代公式为 （4-53）物理意义如图4-1(d)所示，但迭代次数增多。3、对n个未知函数的N-R法对有限元方程 采用修正的N-R法 非线性方程组的迭代格式为 （4-54）其中，

17、是阶方阵。一个三角形单元中 （4-55）而 其中，下面逐一推导。由于中有，而是A的函数，所以 由于，取分段线性插值函数则有 因为，前面已给出单元磁感应强度的离散式其模值为 因此，最终同理可得以及用通式表示 （4-56）三角形单元的右端项矩阵为 （4-57）式中， 其中，那么 （4-59）最后，总雅可比矩阵J由单元雅可比矩阵组合而成，只需将下标相同的相应项叠加即可 （4-60）雅可比矩阵J与刚度矩阵相似，具有对称性和稀疏性，因此可以采用与K一样的存储技术。牛顿拉夫逊法的特点是收敛快，远优于线性化迭代法。缺点是：（1）要有较好的初值，否则可能迭代不收敛。如图4- 3 所示。（2）计算量大，每解一次

18、方程都要重新形成雅可比矩阵，因为每一次的及都变了。为了加快运算速度，出现了许多的改进牛顿拉夫逊法，如，形成一次J，以后每次迭代都用相同的斜率J，但K随磁导率变化。这样做的代价是增加迭代次数，但节省了每次计算的CPU，因此，总的CPU会节省。可以采用两种方法相结合的方法：开始使用线性化方法，进入饱和区后再使用牛顿拉夫逊法。（第八次课）4.4.4 预处理共轭梯度法（ICCG法）ICCG法是解无约束优化问题中的一种基本优化方法。那么，泛函极值问题 （4-61）就可以作为无约束优化问题，未知函数是优化变量，是目标函数。ICCG法的特点：（a）只需存储非零元素，节约内存；（b）降低系数矩阵的条件数（条件

19、数的大小主要取决于主对角线元素数值之间的差值，差值越大，条件数越大，计算误差也越大）；（c）解大型稀疏矩阵方程时，CPU时间比牛顿拉夫逊法的时间少。ICCG法是在梯度法的基础上发展起来的，因此，先介绍梯度法的基本步骤：（1）第k次迭代点用表示，设初值优化变量 （2）从出发，选择搜索方向，在该方向上寻找使目标函数下降的新点（3）在方向上找出新点，其中，称为方向上的步长因子（4）如果处的目标函数满足 或 迭代结束，为优化解。由于优化解是使目标函数取极小值的解，而是极值解的必要条件，故可将其作为迭代终止的判据。上述基本步骤可见，迭代过程的两个关键步骤是 搜索方向； 步长因子。不同的和将构成不同的计算

20、方法。1、最速下降法该方法是将搜索方向定为目标函数下降的方向，即泛函在迭代点处的负梯度方向，即下降的方向 （4-62）步长因子选择 （4-63）可以证明，与正交 这意味着两次搜索方向相互垂直，搜索过程呈现“锯齿”状。 图4-2 最速下降法物理含义其物理含义如图4-2所示。函数在几何上是一空间曲面，它与面相切的点即是它的零极小值点，如图4-2(a)所示。对于空间曲面，如果用一系列平行于的平面相截之，可以得到一族平面曲线，将他们投影于面，得到如图4-2(a)所示的曲线族，称为曲面的等高线。从初值点出发，沿着使I值下降的方向逐步地下降I值，直到降到它的零极小值。该方法的优点： 方法简单、直观，计算存

21、储量小； 开始收敛速度快，对任意初值均可收敛。缺点：越接近极值点收敛越慢，对于条件数大的方程情况尤其严重。通常与其他方法组合使用，如与牛顿拉夫逊法组合使用，先用最速下降法，后用牛顿拉夫逊法。2、共轭梯度法（CG）该方法的第一步用最速下降法，即用目标函数在处的负梯度作为搜索方向，其余各步均用共轭方向作为搜索方向。而共轭方向是用泛函的梯度来构造的。共轭方向的定义：设K为阶正定对称矩阵，P为阶矩阵，若满足 （4-64）则称和为对刚度矩阵K共轭，且和为线性独立。上式表明：第j次的搜索方向经过K矩阵变化后才与正交，这是与最速下降法区别之处。基本步骤为：（1）设初值，计算，令第一次搜索方向（实际上，一开始

22、余量不为零，这是最速下降法）（2）如果，则为解，程序停止。否则按（3）进行（3）按下面计算第k次步长因子 k =0，1，2，计算新点 第k+1次梯度 计算系数 第k+1次搜索方向 转入步骤（2）。该方法优点：只需存储非零元素，节约内存；从理论上讲，在无舍入误差的前提下，对N个变量迭代N次后必然收敛。缺点：计算量大，且对条件数大的方程难以收敛。3、预处理共轭梯度法（ICCG）ICCG法是将刚度矩阵做不完全分解，得到乔列斯基的阵作为预处理阵（L为K的下三角阵，D为K的对角阵），做这样处理的系数阵可以降低条件数，将处理后的系数阵放入CG法的迭代公式中去，令替代CG法中的K，就形成了ICCG法。4.4

23、.5 磁性参数的提取从有限元方程中得到磁矢位后，可以通过下式得到磁感应强度平行平面场 （4-65）轴对称场 （4-66）求出其模值，通过查磁化曲线，或通过插值函数得到，单元的磁阻率，若为轴对称场，§4-5 等参元有限元法在有限元法中，形状函数有两个作用：（1）用于未知函数的展开，即已知单元节点上的位函数值，采用形状函数便可计算单元内任一点的位函数值。称之为形状函数1。（2）用于局部坐标与整体坐标之间的变换，即将整体坐标系中形状复杂的单元通过形状函数（与（1）中的形状函数的概念不同）变换到局部坐标系中的标准单元，称之为形状函数2。这样，单元刚度矩阵的各元素可以在局部坐标系下，用前面论述

24、的方法建立。当形状函数1用于作用（1）时，根据计算精度的要求，取低阶或高阶插值的形状函数。当形状函数2用于作用（2）时，整体坐标中单元的节点数可以与局部坐标中的单元节点数不同，或者说两者的插值函数的阶数可以不同：形状函数1的阶数与形状函数2的阶数相同时，称为等参元有限元法，本节介绍内容；形状函数1的阶数小于形状函数2的阶数时，称为亚参元有限元法；形状函数1的阶数大于形状函数2的阶数时，称为超参元有限元法。4.5.1坐标系整体坐标系：表示场域中所有节点的位置；局部坐标系：表示单元节点在单元中的位置，在单元之外无任何意义。局部坐标的表示方法是用面积比来表示的，故又称为面积坐标系。以三角形单元为例，

25、 P点是三 图4-3 面积坐标定义 角形单元中任一点，定义面积坐标： （6-67）其中，是三角形单元面积，如图4-3所示。面积坐标特点：（1）面积坐标均小于等于1。（2）平行于棱边的任一直线，如直线上的任意节点都具有相同的值，因为三角形的底和高相同；同理，对直线上的任意节点都具有相同的值。如图4-4所示。 图4-4 面积坐标系 （3）将图4-1中的点拉到节点1，即可得到节点1的面积坐标，显然此时和都等于零，而与相等，同理对节点2和节点3，因此有节点1的面积坐标 （4-68）节点2的面积坐标 （4-69）节点3的面积坐标 （4-70）三角形重心（中线的交点）的面积坐标 （4-71） 可见，任一节

26、点对应的面积坐标有3个，但其中只有两个是独立的，因为 （4-72） 所以只要用任意两个即可，即任一节点对应的面积坐标有2个。建立实际单元局部坐标系如图4-5所示，指定的直线为轴，显然，当的那些直线恰是等于定值的直线。指定的直线为轴，显然，当的那些直线恰是等于定值的直线。因此有 图4-5 实际单元局部坐标系 因此，在坐标系中，只需用两个量即可确定单元中各点的位置。将实际单元局部坐标系转换成正交标准单元局部坐标系，如图4-5所示。两坐标之间的关系为： 在标准单元局部坐标系中，直角边边长均为 1。 对于平面上的曲边三角形，如图4-7所示。 图4-6 标准- 坐标系图4-8所示6节点曲边三角形的标准单

27、元坐标系，如图4-9所示。 图4-7 曲边三角形 图4-8 六节点曲边三角形局部坐标系图4-9 标准- 坐标系 图4-10 四边形局部坐标系图4-11 标准- 坐标系 4.5.2 整体坐标与局部坐标之间的转换关系在图4-2中，的面积为用前面定义过的符号可以看出：将其代入面积坐标可以得到局部坐标（归一化后的坐标均小于1） （4-73）用矩阵形式表示为 （4-74）反之，也可得到用面积坐标表示整体坐标的关系： （4-75）比较式（4-71）（4-73）、与前面讲的形状函数可知所以式（4-73）可以写成 （4-74）可以证明，整体坐标与局部坐标变换后，单元顶点的坐标经变换后保持不变；单元之间不会重叠，也不会分离；分界处保持连续。4.5.3 等参元有限元法根据前面所述已知 由于是的函数，而，反推：由于写成矩阵形式J成为雅可比矩阵，它表示了坐标之间的转换关系。其中各元素等于 由此可得其中 可以证明，。这样最后得到有限与方程 等参元法的优点：经过坐标变换后，单元在局部坐标中非常“标准”（等腰直角三角形、正方形），容易写出形状函数的解析表达式。对于边界无法用折线或用折线误差大时，可以用这种方法来提高计算精度。缺点：需要求雅可比矩阵的逆阵