大作业参考样例

1、工程电磁场数值计算报告 小组编号：01 小组成员： 上海交通大学电气工程系2013年12月小组成员信息第 01 小组组员姓名学 号班 级目录 TOC o 1-4 h z u HYPERLINK l _Toc344072861 1.实验任务 PAGEREF _Toc344072861 h 3 HYPERLINK l _Toc344072862 2.实验内容 PAGEREF _Toc344072862 h 4 HYPERLINK l _Toc344072863 2.1有限元法 PAGEREF _Toc344072863 h 4 HYPERLINK l _Toc344072864 2.1.1概述 P

2、AGEREF _Toc344072864 h 4 HYPERLINK l _Toc344072865 2.1.2应用步骤 PAGEREF _Toc344072865 h 4 HYPERLINK l _Toc344072866 2.2实验过程 PAGEREF _Toc344072866 h 4 HYPERLINK l _Toc344072867 2.2.1问题分析 PAGEREF _Toc344072867 h 4 HYPERLINK l _Toc344072868 2.2.2前处理 PAGEREF _Toc344072868 h 5 HYPERLINK l _Toc344072869 2.2.

3、2.1 圆形域划分 PAGEREF _Toc344072869 h 5 HYPERLINK l _Toc344072870 2.2.2.2 节点和单元编号 PAGEREF _Toc344072870 h 6 HYPERLINK l _Toc344072871 2.2.3一些推导 PAGEREF _Toc344072871 h 6 HYPERLINK l _Toc344072872 2.2.3.1 概述 PAGEREF _Toc344072872 h 6 HYPERLINK l _Toc344072873 2.2.3.2 单元内系数矩阵C(e) PAGEREF _Toc344072873 h 6

4、 HYPERLINK l _Toc344072874 2.2.3.4 代数方程推导与求解 PAGEREF _Toc344072874 h 8 HYPERLINK l _Toc344072875 2.2.4边界条件 PAGEREF _Toc344072875 h 9 HYPERLINK l _Toc344072876 2.2.5后处理 PAGEREF _Toc344072876 h 9 HYPERLINK l _Toc344072877 2.2.5.1 概述 PAGEREF _Toc344072877 h 9 HYPERLINK l _Toc344072878 2.2.5.2 磁力线 PAGER

5、EF _Toc344072878 h 9 HYPERLINK l _Toc344072879 2.2.5.3 磁感应强度B PAGEREF _Toc344072879 h 10 HYPERLINK l _Toc344072880 2.2.6 程序介绍 PAGEREF _Toc344072880 h 10 HYPERLINK l _Toc344072881 2.2.6.1 流程图 PAGEREF _Toc344072881 h 10 HYPERLINK l _Toc344072882 2.2.6.2 函数简介 PAGEREF _Toc344072882 h 11 HYPERLINK l _Toc

6、344072883 2.3实验结果 PAGEREF _Toc344072883 h 13 HYPERLINK l _Toc344072884 2.3.1磁感应强度分布 PAGEREF _Toc344072884 h 13 HYPERLINK l _Toc344072885 2.3.2B-曲线 PAGEREF _Toc344072885 h 14 HYPERLINK l _Toc344072886 2.3.3空腔中心B PAGEREF _Toc344072886 h 14 HYPERLINK l _Toc344072887 2.4结果讨论 PAGEREF _Toc344072887 h 17 H

7、YPERLINK l _Toc344072888 3.心得体会 PAGEREF _Toc344072888 h 18 HYPERLINK l _Toc344072889 4.程序附录 PAGEREF _Toc344072889 h 19 HYPERLINK l _Toc344072890 5.参考文献 PAGEREF _Toc344072890 h 27实验任务一置于均匀外磁场H0中的长直圆柱形空腔导磁体如图1.1所示。其结构尺寸1=0.1m, 2=0.2m；圆柱形导磁体的磁导率=0r=10000，外磁场B0=0H0=1T。（1）以标量磁位m为待求量，求该圆柱形空腔磁屏蔽问题的磁场分布；（2）

8、取=0.15m处，由m的离散解，求B的离散解，并画出B-关系曲线（此处为角度）；（3）计算B=0。图1.1 匀强磁场下圆柱形空腔导磁体截面示意图实验内容有限元法概述有限元法是一种高效能、常用的计算方法。有限元法在早期是以变分原理为基础发展起来的，所以它广泛地应用于以拉普拉斯方程和泊松方程所描述的各类物理场中（这类场与泛函的极值问题有着紧密的联系）。自从1969年以来，某些学者在流体力学中应用加权余数法中的伽辽金法或最小二乘法等同样获得了有限元方程，因而有限元法可应用于以任何微分方程所描述的各类物理场中，而不再要求这类物理场和泛函的极值问题有所联系。有限元法的基本思想是：由解给定的泊松方程化为求

9、解泛函的极值问题。应用步骤有限元分析方法大致可以分为以下三个步骤：（1）剖分：将待解区域进行分割，离散成有限个单元的集合。单元的形状原则上是任意的。二维问题一般采用三角形单元或矩形单元，三维空间可采用四面体或多面体等。每个单元的顶点称为节点；（2）单元分析：进行分片插值，即将分割单元中任意点的未知函数用该分割单元中形状函数及离散网格点上的函数值展开，即建立一个线性插值函数；（3）求解近似变分方程：每个单元的场函数是只包含有限个待定节点参量的简单场函数，这些单元场函数的集合就能近似代表整个连续体的场函数。根据能量方程或加权残量方程可建立有限个待定参量的代数方程组，求解此离散方程组就得到有限元法的

10、数值解。实验过程问题分析如图1.1所示，该场域具有明显的对称性，包括沿着z轴方向的平移不变性，和在x-y平面内关于x,y轴的轴对称性。基于场域这种空间位置的对称性，可以将该场域简化成一个1/4圆域。为了分析方便，可以考虑位于第一象限内的场域。由于整个空间为无源区域，可以使用标量磁位求解，标量磁位m满足拉普拉斯方程2m=0，可以借助于有限元法求解。前处理2.2.2.1 圆形域划分将该1/4圆以等间隔划分成层，然后再在每层上以等差数列的形式选取节点数，即第0 层有1个节点，第1 层有2个节点，第2 层有3个节点然后将所有节点与相邻节点连接成三角形即可。实现圆形域划分的程序详见附录中函数div2tr

11、i 和2.2.6 程序介绍部分。划分完成的区域如下图所示：图2.2.2.1 1/4圆形区域划分示意图图2.2.2.1中，ratio=4，d=0.02，红色边界表示=0.1m，绿色边界表示=0.2m。2.2.2.2 节点和单元编号单元和节点编码自最小半径的圆周（含半径为零的圆心）起始，等间隔逐个圆周按逆时针方向规范地连续编排，直至最大半径处的圆周。注意到周向剖分节点数构成公差d=1的等差数列；周边单元数构成公差为d=2的等差数列。利用这个特征，可以很容易对单元和节点进行编号。关于单元和节点的编号可以参见附录程序中函数L2G 函数和2.2.6 程序介绍部分。2.2.3一些推导2.2.3.1 概述根

12、据图2.2.2.1 圆域划分所得的每个单元，下面将推导每个划分单元满足的方程，得到每个单元的系数矩阵C(e)的具体表达式，进而得出整个场域的系数矩阵C的表达式，并将推导得到便于进行数值求解的一组代数方程组。2.2.3.2 单元内系数矩阵C(e)对于图2.2.2.1 中任意三角形单元e，如果已知每个三角形的三个顶点的磁势值，记为e1，e2，e3，那么可以通过插值的方法来求三角形区域内部任意一点的磁势值，记为e(x,y)：ex,y=a+bx+cy (2.1)考虑到单元e的三个顶点的磁势值1，2，3可以用(2.1)表达，于是12,3=1x1y11x2y21x3y3abc (2.2)求出系数a，b，c

13、，并代入(2.1)得e=1xy12s(x2y3-x3y2)(x3y1-x1y3)(x1y2-x2y1)(y2-y3)(y3-y1)(y1-y2)(x3-x2)(x1-x3)(x2-x1)12,3或者e=i=13i(x,y)i (2.3)令a1=x2y3-x3y2，b1=y2-y3，c1=x3-x2，其余系数a2，b2，c2，a3等可以依次按1，2，3的顺序置换而得，则：1=12sa1+b1x+c1y,2=12sa2+b2x+c2y, (2.4)3=12sa3+b3x+c3y,S为三角形单元e的面积：S=12b2c3-b3c2 (2.5)场域拉普拉斯方程2m=0，可以从下式中导出：Wme=2H2

14、dS=2m2dS=2(x)2+(y)2dS (2.6)其中x=12s(b11+b22+b33)(x)2dS=4sb11+b22+b332 =4s123b1b1b1b2b1b3b2b1b2b2b2b3b3b1b3b2b3b3123同理(y)2dS=4sc11+c22+c332 =4s123c1c1c1c2c1c3c2c1c2c2c2c3c3c1c3c2c3c3123因此Wme=2etCee (2.7)其中e=123Ce=C11(e)C12(e)C13(e)C21(e)C22(e)C23(e)C31(e)C32(e)C33(e)Cij(e)=14s(bibj+cicj)2.2.3.3 整个场域系数

15、矩阵C将所有三角形单元的单位长度磁场能量累加可得：Wm=e=1NWme=02tC (2.8)注意：考虑到整体场域有三部分构成，而且这三个部分的本构参数不同，所以由三角形单元e的系数矩阵Ce组合成整个区域的系数矩阵C时，需要做如下处理：即将单元内的系数矩阵Ce以rCe代替，然后再合成整体的系数矩阵C。如此处理以后便可以使用(2.8)来计算这个题目中整个场域的单位长度总磁场能量Wm。2.2.3.4 代数方程推导与求解根据多元函数 取极值的必要条件知，Wmk=0 (2.9)于是0=i=1niCik (2.10)其中，n为总的节点个数，k=1，2，3，n。考虑到节点比较多，所以采用迭代法求解(2.10

16、)。为此将所有的节点分成两类，一类节点的磁势值已知，称作“强制节点”；另一类节点的磁势值位置，称作“自由节点”。因此，迭代过程只是针对那些“自由节点”进行，并取初值为0，如式(2.11)所示：k=-1Ckkj=1,ikniCki (2.11)2.2.4边界条件本题目包含有三个区域，即铁磁介质内部的空腔（区域1），圆环状铁磁介质（区域2），外部空间（区域3），因而包含有三个交界面，即区域 1和2之间（边界1），区域 2 和3 之间（边界2），还有无穷远的边界（边界3）。由于有限元法的数学本质是变分原理，所以在相应泛函取得极值的时候，自然能够满足由于介质的本构参数不同而产生的自然边界，所以边界1和

17、2 自然成立，只需要考虑边界 3。边界 3对应的强制边界条件为：m|=3=H0 x=H03cos (2.12)式(2.12)中3应足够远，从而可以逼近均匀磁场。此外，由于考虑了对称性，标量磁位场还应满足以下两个条件：（1）x轴与B线相重合，mn|y=0=0。这类边界条件自然成立，不需单独考虑。（2）y轴为对称线，可设m|x=0=const=0，为强制边界条件。总之，所有节点可以分为两类：自由节点和强制节点。强制节点为：y轴上的所有结点，磁势为0；最外层结点，磁势值等于其x坐标乘以匀强磁场的磁场强度值H0。2.2.5后处理2.2.5.1 概述后处理部分的功能是：通过有限元法求得的所有节点磁势m的

18、值计算磁感应强度B值，再进行可视化处理，画出磁感应强度的分布情况。2.2.5.2 磁力线若使用矢量磁位A求解，则A等于定值的轨迹即为B线。当使用标量磁位m求解时，磁场强度H=-m，进而可求出每个单元的磁感应强度B。利用栅格分割函数meshgrid，可以在平面上取得一个矩形区域的点阵。判断点阵中每个点所在单元的编号，求得该单元的B值也就是该点的B值。通过矢量场绘制函数quiver，可以以箭头的形式绘出点阵中各点磁感应强度B的大小和方向。2.2.5.3 磁感应强度B通过迭代法求解出了所有节点的磁势值之后，利用插值公式(2.3)和(2.4)，可以求出场域内每个节点的磁势值。然后利用公式H=-m=-x

19、ex-yey可以求得所有单元内部的磁场强度H的值。对r值不同的区域，通过B=H=0rH可以求出磁感应强度B的值。对于三角形顶点处的磁感应强度B，需要单独处理。设对于节点 i，设与其关联的单元分别为ei1，ei2，ei3eim，那么计算节点的磁感应强度可以按照下式进行：Bi=j=1mBeijSeijj=1mSeij (2.13)上式中是Beij与节点i相邻的单元eij内部的磁感应强度值。需要注意的是，由于插值公式是关于x和y是一阶的，所以每个三角形单元内部的磁感应强度值是一个定值。2.2.6 程序介绍2.2.6.1 流程图如图2.2.6.1所示，输入参数 ratio，d，r，然后使用函数cmat

20、r计算出场域的系数矩阵 C，然后作为函数iteration的参数，通过迭代求解出所有节点的磁势值，存入行向量v中。以v为参数传递给函数calc_b ，计算所有单元内部磁感应强度B的幅值B及其x方向和y方向的分量Bx和By。以Bx和By为参数传递给函数draw_B可以画出磁场分布情况，以B为参数传递给函数draw可以画出=0.15m处的B-曲线。图2.2.6.1 程序流程图2.2.6.2 函数简介(1)matr=cmatr(ratio,d,ur)功能：求取系数矩阵的函数。参数说明：ratio表示无穷远边界距离0.2m边界相对外径0.2m的倍率，这里取 5；ur为磁导体的相对磁导率，根据题目，其值

21、为1000；d 为划分有限元时的步长，取d=0.01m；matr 为场域的系数矩阵。(2)v=iteration(ni,ratio,d,matr)功能：求取磁势的函数。参数说明：ni为迭代次数，取 10000 ；ratio 含义同上，取值为5；d 为划分步长，取值 d=0.01m；matr为函数 cmatr 计算得到的场域的系数矩阵；v 为所有节点的磁势，以行向量表示。另一迭代函数iteration2和函数iteration几乎相同，区别在于：iteration2函数无初始化条件，可以进行迭代后的再迭代。(3)bx,by,b=calc_b(ratio,v,d)功能：求取磁感应强度的函数。Rat

22、io含义同上，取值为5；v为函数iteration迭代得到的所有节点的磁势，以行向量表示；d为划分步长，取值0.01m 。(4)n=L2G(ele,loc)功能：节点单元内编号向全局编号转换。参数说明：ele为单元编号；loc为节点在单元内的编号，取值 1，2，3；n为节点全局编号。(5)x,y=calc_coord(g,d)功能：由节点全局编号求取节点x,y坐标。参数说明：g为节点全局编号，d为划分单元时的划分步长，x,y 为节点坐标。(6)draw(b,d)功能：画出=0.15m处的B-关系曲线，并显示曲线的图像。参数说明：b为行向量，存储所有单元内磁感应强度值；d为划分步长。(7)div

23、2tri(ratio,d)功能：将圆域分割成三角形单元，并以图像显示出来。参数说明：ratio含义同前，d为划分步长。(8) draw_B(step,d,bx,by)功能：绘制0.6m*0.6m区域内的磁场分布情况。参数说明：step为栅格相邻两列间的距离，d为划分步长，bx和by为各节点磁势的x分量和y分量值。(9) area=triangle_area(n,d)功能：返回指定序号三角形单元的面积。参数说明：n为单元的全局编号，d为步长。(10) n=which_triangle(x,y,d)功能：判断某一点位于哪一块三角形区域内。参数说明：x，y为该点坐标，d为步长，返回值n为三角形单元的全局编号。实验结果磁感应强度分布通过draw_B函数可以绘制出圆柱形空腔磁屏蔽问题的磁场分布图，如图2.3.1所示。绘图区域为0.6m*0.6m的正方形。虽然图 2.3.1中只绘制了0.6m*0.6m区