电磁场数值计算作业报告

《电磁场数值计算》—有限元法报告第一题（一）问题描述及数学物理模型的建立有一矩形场区，尺寸为6×9，如图1所示，在区域内部J=0，的左边界A=0,右边界A=100，上下边界满足：0An∂=∂，媒质均为的磁导率为μ，利用有限元法求A 的分布首先，2222000(0,0)0;1000x x a y y b A Ax a y b x y A A A A yy ====⎧∂∂+=<<<<⎪∂∂⎪⎪==⎨⎪∂∂⎪==⎪∂∂⎩图1 求解场域及网格划分如果利用有限元来求解，那么对应的泛函为（注：24,l l 分别为上边界和下边界）：2422221()221min(0;0)2S l l S A A AF A JA dxdy dl x y n A A A dxdy J x y n μ+⎡⎤⎛⎫∂∂∂⎛⎫=+--⎢⎥ ⎪ ⎪∂∂∂⎝⎭⎢⎥⎝⎭⎣⎦⎡⎤⎛⎫∂∂∂⎛⎫=+===⎢⎥⎪ ⎪∂∂∂⎝⎭⎢⎥⎝⎭⎣⎦⎰⎰⎰⎰⎰第一类边界条件为： 00;100x x aAA====（二）场域的离散化及单元节点的编号处理网格的划分如图1所示，先求出每个单元的刚度矩阵，然后再进行整体合成，得到系数矩阵，具体的做法可参考教材，注意，这里左边界和右边界属于第一类边界，因此，利用下式进行计算：11I 112II=-⎧⎨=⎩K A R K dA d 这里，因为J=0,所以1R =0，11I 12=-K A K d节点编号如图1所示，这里，由于网格的划分和节点的编号都不是很规则，但节点数目和网格数目不是很多，所以，在程序中，通过穷举的办法，找到相应单元的（i ，j ，m ）的值，以及相应的节点坐标。

各单元的（i,j,m）（按逆时针方向）如下：1(3,7,1) 2(3,8,7) 3(9,8,3) 4(5,9,3) 5(6,5,3) 6(6,3,4) 7(4,3,1) 8(4,1,2) 9(4,2,10) 10(11,4,10) 11(12,4,11) 12(12,6,4) 节点i的坐标为：if i=1,3,5 then x(i)=3, y(i)=1.5(i-1)if i=2,4,6 then x(i)=6, y(i)=1.5iif i=7,8,9 then x(i)=0, y(i)=3(i-7)if i=10,11,12 then x(i)=9, y(i)=3(i-10)边界节点及其位函数的值为：7（0） 8（0） 9（0）10（100） 11（100） 12（100）然后建立有限元方程组，进行求解，方程组的建立过程参考教材（三）利用数值模拟计算按照有限元法计算的相关步骤，用FORTRAN 90语言编写的程序清单如下：程序运行结果如下：理论结果与实际结果的比较本题对应的偏微分方程及相应的定界问题为：2222000(0,0)0;1000x x a y y b A Ax a y b x y A A A A yy ====⎧∂∂+=<<<<⎪∂∂⎪⎪==⎨⎪∂∂⎪==⎪∂∂⎩数学上不难求出下面偏微分方程的解为：100(,)A x y x a=可以看出，理论结果与数值计算的结果几乎是相同，这是因为数值计算的方法给予有限元的离散，而每个单元都是用线性插值函数来近似的，因为这里A 与坐标x,y 本来就是线性的关系，因此，用线性插值函数来近似真实函数可以得到0误差，这时，误差的主要来源是解方程组的迭代误差。

上述程序中迭代误差取为0.0000000001（1010-）因此，本题得到的精确度极高。

但是，一般的问题A 与坐标x,y 是线性的关系的情况很少，这种情况下，由于线性插值函数近似代替真实函数的截断误差便不可被忽略了，而且一般来说，在两个单元形状类似的情况下，单元面积越小，截断误差也越小。

下面的几个例子将会说明这一点。

第二题2.1（一）问题描述及数学物理模型的建立 如图2,矩形场域ABCD 内的电荷密度为0，AB=CD=a;AC=BD=0.6a;在边界AC,CD,BD 上，电位为0，在顶端BD 上电位分布为f(x)。

分别求出()100f x =（本题的情形） 以及()100sinxfx aπ=时对应的电位分布。

网格可以用边长为1的正方形网格，对于 本题，将计算结果与有限差分法进行对比图2此问题对应的偏微分方程及相应的定解问题为：22220,00.60,00.5,00.50.5,00((0,),(0,0.6))0()0x y a y x a y a x a x a y ax a y a x y f x x ϕϕϕϕϕϕ=<<=<<=<<=<<∂∂⎧+=∈∈⎪∂∂⎪⎪==⎨⎪∂⎪==⎪∂⎩这里：()100f x =对应于本题的情况，而()100sin xf x aπ=则对应于下题的情况。

此题对应的泛函及等价的变分问题为：()()()22/20,00.500.5,000.5,0.6()d d min(0;0)20;100D x y a x a y x a y a F x y x y n εϕϕϕϕρϕϕϕ=≤≤<≤=<≤=⎧⎡⎤⎛⎫∂∂∂⎛⎫⎪=+===⎢⎥ ⎪ ⎪⎪∂∂∂⎝⎭⎢⎥⎝⎭⎨⎣⎦⎪===⎪⎩⎰⎰ 由于题目没有给我们划分网格，也没有给我们定好节点，这时，要自己划分网格，为了容易编写程序，减少用户输入网格和节点的信息这一步，采用网格自动剖分的办法，网格划分如图3所示(网格单元及节点的编号规则下面介绍)：图3 网格的划分（注：上题未给出标号的节点及单元可按照图中的规律仿照进行编号。

）这里，先讲一下节点和单元的编号规则： ①首先，有三条边上的节点属于第一类边界条件的节点，因此，这些点的编号要靠后，先对电位未知的节点进行编号，然后再对电位已知的节点编号。

②对于内部电位未知的节点，单元和节点的编号采用与课本相同的处理办法，按照一定的顺序和规律进行编号。

设长和宽方向每行分别有y ,x N N 个结点。

对于内部的及节点，节点编号h （i,j ）与（i ，j ）的关系为：()y y y y (2)(2)12;21(1)(2)(1)(,)(1)(2)1(1;1)(1)(2)1(1;)x y xx x xx y i N j i N j N N N i j h i j N N N j j i N N N j i j N ⎧--+-≤≤≤≤-⎪--+=⎪=⎨--++->=⎪⎪--++->=⎩ 先确定各个节点的编号信息，设∆x 与∆y 分别是x 轴和y 轴方向直角三角形的直角边长。

则节点（i ，j ）的坐标为：(,)(1)(,)(1)x i j i xy i j j y=-∆⎧⎨=-∆⎩ 再确定节点的统一编号。

这里，做一个约定，用tp 代表网格的类型，tp=0代表a 类三角形单元，tp=1代表了b 类三角形单元。

如图4所示：图4 两类三角形单元因此，如果网格位于第i 行，j 列，那么三个节点相应的位置就应该为： 对于a 类三角形单元:I(i+1,j+1) J(i,j+1) M(i,j) 对于b 类三角形单元:I(i,j) J(i+1,j) M(i+1,j+1) 统一起来，可以写为：I(i+1-tp,j+1-tp) J(i+tp,j+1-tp) M(i+tp,j+tp)因此，对于任一个三角形单元，如果已知单元的编号，只要知道该三角形单元位于哪一行哪一列，三角形单元是什么样的（a 类还是b 类），就可以知道该三角形单元三个节点I,J,M 的统一编号。

下面讨论如何根据三角形单元的编号，确定该单元所在的行列和三角形单元的类型的信息。

对于编号规则的三角形单元（这部分三角形单元的三个顶点的节点电位都是未知的），编号E 与单元最左下方的节点所在行列(i,j)的关系为： 对于a 类三角形单元：JMI J(,)2(2)(3)1(2,31;2,32)y x y E i j i N j i N j N =--+-=-=-对于b 类三角形单元：(,)(23)(3)1(2,31;2,32)y x y E i j i N j i N j N =--+-=-=-对于这些三角形单元，编号的最大值为: 0max 2(2)(3)x y E N N =--（一）当0max E E ≤时，根据上面计算E(i,j)的公式可知：a 类三角形单元与b 类三角形单元的区别在于3y N -前面是奇数还是偶数，因为2y j N ≤-，因此，24y j N -≤-，而2(2)(3)2(,)1(23)(3)2y y i N j E i j i N j --+-⎧⎪-=⎨--+-⎪⎩可知，①当(,)13y E i j N --的整数部分为偶数时，单元为a 类三角形单元，即tp=0，反之，则为b 类三角形单元，tp=1.②2j -为((,)1)(3)y E i j N --的余数，用[x]表示对x 进行取整运算，mod(a,b)表示整数a 除以整数b 的余数。

那么tp 与j 的计算如下：(,)1mod ,23mod((,)1,3)2y y E i j tp N j E i j N ⎧⎛⎫⎡⎤-⎪= ⎪⎢⎥⎪ ⎪-⎨⎢⎥⎣⎦⎝⎭⎪=--+⎪⎩ 知道了这两个数之后，i 也不难计算，i 的计算公式如下：1(,)1423y E i j i tp N ⎧⎫⎡⎤-⎪⎪=+-⎢⎥⎨⎬-⎢⎥⎪⎪⎣⎦⎩⎭上面讨论的是当 0max E E ≤的情形（二）当0max E E >时，处理起来比较简单，只是相应的情形比较多，例如，按照上面的标号规则，如果0max 0max 2(1)x E E E N <≤+-,在上图中相当于：（32<E ≤42）那么j=1,此时，E 与i 的关系为：0max 0max(0)(11)1(1)x x E iif tp E i N E N i if tp +=⎧=≤≤-⎨+-+=⎩采用与上面类似的办法，可以求出i 和tp 的表达式，下面仅给出E 取值不同情况下(i,j,tp)的计算公式：①0max 0max 2(1)x E E E N <≤+-时[]()()0max 0max mod 11,2mod 1,111x x tp E E N i E E N j ⎧=---⎪⎪=---+⎨⎪=⎪⎩②0max 0max 2(1)2(1)2(2)x x y E N E E N N +-<≤+-+-时：()()0max 0max mod 212,21mod 1,12x y y tp E E N N i j E E N ⎧⎡⎤=--+-⎣⎦⎪⎪=⎨⎪=---+⎪⎩③当0max 0max 2(1)2(2)2(1)2(2)2(2)x y x y x E N N E E N N N +-+-<≤+-+-+- 时，即0max 2(1)2(2)2(1)(1)x y x y E N N E N N +-+-<≤--时：()()0max 0max mod 2252,2mod 225,221x y x x y x y tp E E N N N i E E N N N j N ⎧⎡⎤=---+-⎣⎦⎪⎪=---+-+⎨⎪=-⎪⎩以上给出了如何根据网格的编号计算三个顶点的位置，并给出三角形单元的类型，知道了这些就可以计算三角形单元三个节点的整体编号，进一步计算系数矩阵，求解有限元方程组。