有限元C++编程实践

基于有限元算法的编程实践学号：2011043010031姓名：廖校毅电子科技大学物理电子学院【摘要】本文就有限元算法在电磁场分析中的应用展开研究与实践，从电磁场的Maxwell方程出发,根据电磁场的边值问题及变分公式建立了有限元方程组。

具体在实践中，将这些知识运用到C++语言和Matlab中，并将这两种语言有机结合，编程并实现二维FEM。

程序最后通过计算含两种介质的电位槽电位分布来验证其正确性。

关键词: 有限元变分法C++ MatlabThe Programming Practice Based on The Finite Element Algorithm Student ID：2011043010031Name：Liao Xiaoyi University of Electronic Science and technology &School of PhysicalElectronicsAbstract 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 的简单形式的泛函21[](,,')x x J y F x y y dx =⎰函数y 稍有变化时:21[](,,'')x x J y y F x y y y y dx δδδ+=++⎰21[][][(,,)(,,)]d x x J J y y J y F x y y y y F x y y xδδδ∆=+-'''=++-⎰上式用多元函数的泰勒公式加以展开：2222222231{[][()2d 2()]}FF F F y y y y y yJ xFy y y y y y J J J Jδδδδδδδδδδ∂∂∂'++'∂∂∂∆=∂∂''+++''∂∂∂=+++≈⎰称为泛函的一阶变分。

泛函同一般函数的极值特性相似，泛函的取极值的条件是泛函的一阶变分为零，泛函值为拐点的条件是二阶变分为零。

所以，根据一阶变分为零，有下式成立：2211-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.三角形节点编号按逆时针顺序编号图 2 网格划分示意图J δ2.2.1 有限元分片差值取任意一个网格如图，在有限单元（三角形）上进行分片线性插值，基函数为:123(,)e x y x y ϕααα=++123123123i i ij j j m m mx 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 αϕϕϕαϕϕϕαϕϕϕ⎧=++∆⎪⎪=++∆⎨⎪=++∆⎪⎩ 其中()12i j m m ji j m i m ji 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 mx 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 ϕϕϕϕϕ⎡⎤⎢⎥⎡⎤==⎢⎥⎣⎦⎢⎥⎣⎦下面进行单元的泛函分析[][]222ee e e e D J J dxdy x y εϕϕϕϕ⎡⎤⎛⎫⎛⎫∂∂≈=+⎢⎥ ⎪ ⎪∂∂⎢⎥⎝⎭⎝⎭⎣⎦⎰⎰()2/2ei i j j m m b b b xϕαϕϕϕ∂==++∆∂()224ee 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 b b b b b b where K b b b b b b b b b b b b ε⎧⎫⎪⎪=⎨⎬∆⎪⎪⎩⎭ 同理：{}{}{}22e e Te e e D dxdy K y ϕεϕϕ⎛⎫∂= ⎪∂⎝⎭⎰⎰ {}2 4 i i i j i m j i j j j m em 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 ε⎧⎫⎪⎪=⎨⎬∆⎪⎪⎩⎭单元上总的泛函:[]222ee e eD J dxdyx y εϕϕϕ⎡⎤⎛⎫⎛⎫∂∂=+⎢⎥ ⎪⎪∂∂⎢⎥⎝⎭⎝⎭⎣⎦⎰⎰{}{}{}{}{}{}{}{}{}12112212T Te e e e e e Te e eK K K ϕϕϕϕϕϕ=+=图 2 插值示意图 123,,ααα{}{}{}12 e e ee e e ii ij im e e eji 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 ε==+∆总体合成：改写扩充单位阵 到所有单位，即把扩充部分添零，为方便总体矩阵的处理。

[][]{}{}{}{}{}{}111212NNT e e e e TJ J K K ϕϕϕϕϕϕ==⎛⎫≈= ⎪⎝⎭=∑∑()1,1,2,Neij ij e where K K i j N ===∑变分问题被离散化的多元二次函数的极值问题：[](){}{}{}12,11,,,21min 2TN Nij i j i j J J K K ϕϕϕϕϕϕϕϕ=≈===∑ 根据函数极值理论:0 (1,2,,)iJi N ϕ∂==∂ 所以有下式成立：,10 (1,2,,)Nijji j K j N ϕ===∑{}{}{}0K ϕ=三、 仿真实现[4][5]3.1 仿真实例采用有限元法仿真并求解如图的含有介质体的静场电位分布。