基于有限元算法的编程实践
学号：2011043010031 姓名：廖校毅
电子科技大学物理电子学院
【摘要】本文就有限元算法在电磁场分析中的应用展开研究与实践，从电磁场的Maxwell方程出发,根据电磁场的边值问题及变分公式建立了有限元方程组。具体在实践中，将这些知识运用到C++语言和Matlab 中，并将这两种语言有机结合，编程并实现二维FEM。程序最后通过计算含两种介质的电位槽电位分布来验证其正确性。
关键词: 有限元变分法C++ Matlab
The Programming Practice Based on The Finite Element Algorithm Student ID：2011043010031 Name：Liao Xiaoyi University of Electronic Science and technology &School of Physical
Electronics
Abstract In this paper, we take the application of finite element method in electromagnetic field analysis into research and practice. Starting from Maxwell equations of electromagnetic field, the electromagnetic field boundary value problems and variational formula established the finite element equations. In specific practice, this knowledge will be applied to the C++ language and Matlab, and the organic combination of two languages, programming and implementation of two-dimensional FEM. Finally, through the program to verify the validity of the calculation of potential distribution in channel potential containing two kinds of medium.
Key words The finite element method The variational method C++ Matlab
一、前言
在数学中，有限元法（FEM，Finite Element Method）是一种为求得偏微分方程边值问题近似解的数值技术。它通过变分方法[1]，使得误差函数达到最小值并产生稳定解。类比于连接多段微小直线逼近圆
的思想，有限元法包含了一切可能的方法，这些方法将许多被称为有限元的小区域上
的简单方程联系起来，并用其去估计更大区域上的复杂方程。它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的(较简单的）近似解，然
后推导求解这个域总的满足条件(如结构的
平衡条件），从而得到问题的解。这个解不
是准确解，而是近似解，因为实际问题被较简单的问题所代替。由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。因其理论依据的普遍性，作为一种声誉很高的数值分析方法，,已被普遍推广并成功地用来解决各种结构工程、热传导、渗流、流体力学、空气动力学、土壤力学、机械零件强度分析、电磁工程问题等，并且大量可靠的基于有限元算法的软件如
ANSYS 被开发和使用。
有限元法有以下特点：一，离散化过程保持了明显的物理意义。因为变分原理描述了支配物理现象的物理学中的最小作用原理（如力学中的最小势能原理、静电学中的汤姆逊定理等）,基于问题固有的物理特性而予以离散化处理，列出计算公式，是保证方法的正确性、数值解的存在与稳定性等前提要素。二，优异的解题能力。与其他数值方法相比较，有限元法在适应场域边界几何形状及媒质物理性质变异情况的复杂问题求解上，有突出优点：不受几何形状和媒质分布的复杂程度限制；不同媒质分界面上的边界条件是自动满足的；不必单独处理第二、三类边界条件；离散点配置比较随意，通过控制有限单元剖分密度和单元插值函数的选取，可以充分保证所需的数值计算精度。
二、 算法原理 2.1 变分基本原理：
对于一般的问题，可以给出下列对应于一个
自变量x 的简单形式的泛函
2
1
[](,,')x x J y F x y y dx =?
函数y 稍有变化时:
2
1
[](,,'')x x J y y F x y y y y dx δδδ+=++?
2
1
[][]
[(,,)(,,)]d x x J J y y J y F x y y y y F x y y x δδδ?=+-'''=++-?
上式用多元函数的泰勒公式加以展开：
22
222
22231{[][()2d 2
()]}F
F F F y y y y y y
J x
F
y y y y y y J J J J
δδδδδδδδδδ???'++'????=??''+++''
???=+++≈
?
称为泛函的一阶变分。泛函同一般函数的极值特性相似，泛函的取极值的条件是泛函的一阶变分为零，泛函值为拐点的条件是二阶变分为零。 所以，根据一阶变分为零，有下式成立：
22
1
1-0
'[](,,')0
'x x x x F d F y dx y J y F x y y dx F
e h y ???
??=? ?????
?=??????=?????
???
2.2 有限元算法基本原理[2]：
2.2.1 有限元网格划分的基本原则
鉴于实际应用中三角形的网格剖分方式比较常见，下面以三角形网格剖分为例，
介绍有限元网格划分基本原则: 1.不同媒质的分界线，不容许跨越分界线的三角元
2.三角元的边逼近边界
3.三角形要求不要太尖或太钝
4.顶点编号相差不能太悬殊，对多区域的编号，按区域连续编号，
5.三角形节点编
号按逆时针顺序编号
J δ
图 2 网格划分示意图
2.2.1 有限元分片差值
取任意一个网格如图，在有限单元（三角形）上进行分片线性插值，基函数为:
123(,)e x y x y ?ααα=++
123123123i i i
j j j m m m
x y x y x y ?ααα?ααα?ααα?=++?
=++??
=++? 由线性代数知识从而有：
()()()123/2/2/2i i j j m m i i j j m m i i j j m m a a a b b b c c c α???α???α????=++?
??=++???
=++??? 其中 ()1
2
i j m m j i j m i m j
i j j i a x y x y b y y c x x b c b c =-=-=-?=- 把 代入三角元上的线性插值函数,
[],,1
(,)()2()() =(,)
e i i i i j j j j m m m m e s s i j m
x y a b x c y a b x c y a b x c y N x y ?????=
++?
++++++∑1
(,)() 2 (,,)
e s s s s N x y a b x c y s i j m =++?
=其中：
{}{}(,), , i e e e e i
j m j e e m x y N N N N ???????
????==????????
下面进行单元的泛函分析
[][]22
2e
e e e e D J J dxdy x y ε??????
??????≈=+?? ? ???????????
??
()
2/2e
i i j j m m b b b x ?α????==++??
()224e
e i i j j m m D dxdy b b b x ?ε
ε??????=++ ?????
??
{}{}{}{}1 ,, 4 i i i j i m
i T
i j m j i j j j m j e e e m i m j m m m bb bb bb b b b b b b K b b b b b b ?ε???????????????=?????????
??
??＝
{}1 4 i i i j i m
j i j j j m e
m i m j m m bb bb bb where K b b b b b b b b b b b b ε???
?=?????
??
图 2 插值示意图 123,,ααα
同理：{}{}{}2
2e e T
e
e e D dxdy K y ?ε?????= ???
??? {}2 4 i i i j i m j i j j j m e
m i m j m m c c c c c c where K c c c c c c c c c c c c ε??
?
?=?????
??
单元上总的泛函:
[]2
2
2e
e
e
e D J dxdy
x y ε?????
??
????=+?? ? ?????????????
{}{}{}{}{}{}{}{}{}121122
12T T
e e e e e e T
e e e K K K ??????=
+=
{}{}{}12 e e e
e e e ii ij im e e e
ji jj jm e e e mi mj mm K K K K K K K K K K K K =+??????=???????
?其中()
4e e rs sr r s r s where K K b b c c ε
==
+?
总体合成：改写扩充单位阵 到所有单位，即把扩充部分添零，为方便总体矩阵的处理。
[][]{}{}{}
{}{}{}111212
N
N
T e e e e T
J J K K ??????==?
?≈= ???
=
∑∑
()1
,1,2,
N
e
ij ij e where K K i j N ===∑
变分问题被离散化的多元二次函数的极值问题：
[](){}{}{}12,1
1,,,2
1min 2T
N N
ij i j i j J J K K ????????=≈=
==∑
根据函数极值理论:
0 (1,2,,)i
J
i N ??==?
所以有下式成立：
,1
0 (1,2,,)N ij j i j K j N ?===∑ {}{}{}0K ?=
三、 仿真实现[4][5]
3.1 仿真实例
采用有限元法仿真并求解如图的含有介质体的静场电位分布。其中a,b 自定。
图 3长直接地金属槽的横截面
e K
3.3 网格划分情况
情况一：总共划分191个点，共323个网格。
图 5 网格和节点划分情况
情况二:网格加密，共719个点，1352
个网格
图 6 加密后的网格和节点
情况三:介质高度改为7，网格加密，共727个点，1360个网格
3.4 仿真结果
由于计算出的结果是离散无规律的点
的电位，不能直接用matlab画图,下面展示
的是这些点的电压情况，越高的点电压越高。后面用插值拟合的方法，通过matlab画二
维电位图颜色越亮表示电位越高。
这是情况一的结果，网格不够密：图7 三维示意图---点图图9 二维电势图
情况二：网格加密后的情况
图8 插值计算后的三维示意图
图11加密计算后的三维示意图
图12加密后的二维电势图
情况三：介质高度增加为7的情况下图13 增加介质高度且加密后的三维点图图14增加介质高度且加密后的三维图
图15增加介质高度且加密后的二维图
四、分析总结
由上面的电位图可以看出电位分布从
右到左呈现出一种递减的趋势,且从分界面
开始，真空（空气视为真空电子环境）和介质里面的斜率不一样，介质里面斜率更加低，这与常识符合，说明以上算法和程序有效，即：通过有限元算法可以有效仿真实现电磁场的计算和模拟。对比网格加密前后的图像，可以明显看出图像的平滑性显著加强，边缘变化呈现出一种缓慢变化而非突变，这与实际物理过程像吻合。说明通过网格加密可以增加计算的精度，但是，程序运行时间显著增加，由原来的3s增加到10s。对比介质高度增加后的情况可以看出，介质高度增加以后，真空里面部分的斜率更大，即电势变化更剧烈。
图10加密后的三维示意图---点图
五、引文
[1] 老大中变分法基础(第2版) 国防工业出版社
[2] 许学军有限元方法的数学理论科学出版社
[3] 孙志中吴宏伟等计算方法与实习东南大学出版社
[4] 徐士良葛兵计算机软件基础清华大学出版社
[5] 苏彦华张宏林Matlab7.0从入门到精通人民邮电出版社
surf(x,y,z)