电磁场数值计算作业报告

《电磁场数值计算》—有限元法报告

第一题

（一）问题描述及数学物理模型的建立

有一矩形场区，尺寸为6×9，如图1所示，在区域内部J=0，的左边界A=0,右边界A=100，上下边界满足：

0A

n

∂=∂，媒质均为的磁导率为μ，利用有限元法求A 的分布

首先，2222000(0,0)0;100

0x x a y y b A A

x a y b x y A A A A y

y ====⎧∂∂+=<<<<⎪∂∂⎪⎪

==⎨⎪

∂∂⎪==⎪∂∂⎩

图1 求解场域及网格划分

如果利用有限元来求解，那么对应的泛函为（注：24,l l 分别为上边界和下边界）：

2422

22

1()221min(0;0)2S l l S A A A

F A JA dxdy dl x y n A A A dxdy J x y n μ+⎡⎤⎛⎫∂∂∂⎛⎫=+--⎢⎥ ⎪ ⎪∂∂∂⎝⎭⎢⎥⎝⎭⎣⎦⎡⎤⎛⎫∂∂∂⎛⎫=+===⎢⎥

⎪ ⎪∂∂∂⎝⎭⎢⎥⎝⎭⎣⎦

⎰⎰⎰⎰⎰

第一类边界条件为： 0

0;100

x x a

A

A

====

（二）场域的离散化及单元节点的编号处理

网格的划分如图1所示，先求出每个单元的刚度矩阵，然后再进行整体合成，得到系数矩阵，具体的做法可参考教材，注意，这里左边界和右边界属于第一类边界，因此，利用下式进行计算：

11I 112II

=-⎧⎨

=⎩K A R K d

A 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.5i

if 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 A

x a y b x y A A A A y

y ====⎧∂∂+=<<<<⎪∂∂⎪⎪

==⎨⎪

∂∂⎪==⎪∂∂⎩

数学上不难求出下面偏微分方程的解为：

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 =（本题的情形） 以及()100sin

x

f

x 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 a

x a y a x y f x x ϕϕ

ϕϕϕϕ=<<=<<=<<=<<∂∂⎧+=∈∈⎪∂∂⎪⎪

==⎨⎪

∂⎪==⎪∂⎩

这里：()100f x =对应于本题的情况，而()100sin x

f 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 个结点。