电磁场有限元Matlab解法

《电磁场有限元Matlab解法》由会员分享，可在线阅读，更多相关《电磁场有限元Matlab解法（12页珍藏版）》请在金锄头文库上搜索。

1、电磁场有限元Matlab数值解法摘要有限元法(FEM)在电磁场的问题中十分有效，能帮我们解决很多电磁场的问 题。本文主要对电磁场的有限元方(FEM)法进行学习和研究.利用Matlab进行 编程仿真.对特定问题进行求解.模拟仿真出在一定边界条件下的电势分布 图。该方法十分方便有效.能直观的得到我们想要的电势分布图。0引言有限元的原理在数学上首先由K Courant于1943年提出，早期在力学中用 于结构分析。五十年代初期.由于工程分析的需要.有限元法在复杂的航空结构 分析中最先的到应用.而有限元法(Finite Element Method)这一名称在20世纪 60年代由R. W. Clough

2、的首先提出。三十多年来.以变分原理为基础建立起来的 有限元法.因其理论依据的普遍性.不仅被广泛地应用于各种结构工程，而且作 为一种声誉很高的数值分析方法己被普遍推广并成功地用來解决其他工程领域 中的问题.例如热传导、流体力学、空气动力学、机械零件强度分析等。1968年 有限元法开始用于求解电磁场的问题有限元方法在电磁场中的应用至今己有近 50年的历史。随着科学技术的发展.计算机性能得到提高.算法的不断被优化.有限元理 论及技术在电磁应用方面己愈加成熟并取得许多研究成果。有限元法作为一种强 有力的工程分析方法被广泛地应用于各种研究领域有限元法同样是用于各类电 磁场、电磁波工程问题定量分析与优化设

3、计的最主要的数值方法，并且无一例外 地是构成各种先进、有效的计算软件包的基础。目前，有限元分析己成为计算机 辅助设计的一个重要组成部分。1. 有限元简介有限元法的基本思想是将结构离散化用有限个容易分析的单元来表示复杂 的对象.单元之间通过有限个节点相互连接.然后根据边界条件综合求解。由于 单元的数目是有限的.节点的数目也是有限的.所以称为有限元法(FEM, Finite Element Method)。有限元法是以微分方程为基础的数值方法。最初主要是利用变分原理将微分 方程变为等价的变分方程.经改进的Ritz法则将微分方程的求解变为代数方程的求解问题。由于并非任意的微分方程都能找到与其等价的变

4、分方程，使得上 述形式的有限元法的应用受到了一定的限制。用加权余量的伽辽金法直接从微分 方程出发构造有限元方程，就可以突破变分原理的限制。有限元最大的特点是：先通过各种适当的形式将求解域划分成有限个单元. 再在每个单元中构造分域基函数.利用Ritz法或伽辽金法构造代数形式的有限 元方程。有限元法的最大优点是其离散单元的灵活性。相对而言.有限元法可以更精 确地模拟各种复杂的几何结构，并通过选取样点的疏密情况來适应场分布的不同 情况.既能保证计算精度的要求.又不增加过大的计算量。另一优点是所形成的 有限尤方程组的系数矩阵是稀疏的对称阵.这非常有利于代数方程组的求解。与 其他数值计算方法相比较，有限

5、元法在适应场域边界几何形状以及媒质物理性质 变异情况复杂的问题求解上有突出的优点即方法应用不受边界形状和媒质性质 的限制.而且不同媒质分界面上的边界条件是自动满足的，第二三类边界条件不 必作单独处理。此外离散点配置比较随意.并且取决于有限单元剖分密度和笊元 插值两数的选取.可以获得令人满意的数值计算精度。有限元法还可以方便地编 写通用计算程序，使之构成模块化的子程序集合，适应计算功能延拓的需要.从 而构成各种高效能的计算软件包.有限元法的发展与应用前景令人瞩目.2. 有限元问题分析所需求解的电场区域如图2-1所示。该图描述的为屏蔽微带传输线结构。假 设边界上的即屏蔽导体上电势10,内导体即条带

6、上电势U=l内部上部分介丿页 相对介电常数电荷密度P=0下部分介质相对介电常数e=9.电荷密度p =0。由于其几何结构是对称的.我们取其一半进行研究.如图2-2所示。其中电 势在左边边界上由于对称分布在法向方向导数为0o图2-1求解区域图2-2等效求解区域设求解区域长宽均为1而条带处于1/4位置且长度为1/4。电势满足如下的方程：-（刃） = 0 on Q 严 0Don TDri-（rV） = 0 on TN把求解域按如图2-3所示方法进行剖分并对节点和单元进行编号.21222324一25345图2-3网格划分每个区域的电势可以用0（匕刃二“ +加亠刃 来表示。把三角形的三个顶点带 入可得疔二

7、加解上述三个方程得到的abc带入 八（1）以+肚+尸 中得:1对方程（刃0） = 0乘权函数wi并在整个区域积分得一1严V-（rV（p）t/2=0利用矢量恒等式叫（“）（忙BUV叫和高斯定理 J （wV（p）t/Q. = ；（vVjV（pd得J四叫VMQ=切乳（凶甲）书带入边界条件鬧 = 0切TN得使用伽辽金方法取权函数 叱二M心1,2，N带入上式我们得以L NN, TN】dQ =eVM ZN： dQJl=l uN 工&阴#心1,2,N可以写为小，可描述为1KIW = WK = eVN VN dQ L二刃汕竝 其中，;Ju 11.11Kij可表示为各个单元贡献的总和：clrsl每个单元的&；=

8、 rVN；叭VNjdQ =（曙曙 +娜于）空T p =0o以上为理论上的分析。然后利用Matlab进行编程仿真。仿真程序见程序附录。3. 仿真结果及分析仿真结果如图所示，分别为：4-1平面电位分布图；4-2电位分布三维立体图；4-3电位等势线图。102030405060708090100100908070605040302010图4-1平面电位分布图108060.402000 4-2电位分布三维立体因10.90.80.70.60.50 40 30.20.1000.10203040.50.60.70.80.914-3电位等势线图从上面几幅图中我们可以直观的看到电势在所给的求解区域中的分布生动 形

9、象。且边界处的条件吻合得很好。例如上下边和右边电势都为0中心处电势 为1左边界处可以看到等势线垂直于边界.则电场法向分量为零，也与边界条 件很符合。通过有限元的方法.我们用Matlab编程仿真出了静电场的电场分布 图。5 总结通过此次有限元法的研究和学习，成功地仿真出电场的分布图.收获很大。 在刚开始接触有限元法时.完全没有思路.对Matlab软件也一窍不通。后来通 过大量的文献査阅和学习，逐渐理解了有限元法的基本思路和编程方法。在仿真 的过程中.也遇到过各种问题.挑战极大.后来通过反复看文献理解有限元的思 想.学习Matlab的编程方法才慢慢有了一定的体会.最终才得到了较为正确的 结果，也学

10、到了对其他电场区域仿真的基本思路。这次的仿真.不仅仅是一次完 成作业.更是一次挑战，不仅学到了课本上的知识，还学到了 Mat lab的编程技 巧.收获远不止一次完成任务这么简单。虽然仿真结果基本吻合.但是此次仿真还有很多不足之处。首先是只仿真出了电势的分布图.而没有电场分布图。然后是程序算法不够简单.当剖分的网格 数较多时计算速度很慢。还有就是网格的划分算法不够好，应该尽量划分为正三 角形，來减小误差，以后还要多学习这方面的知识。最后，由于时间比较紧.理 论推到中的公式又特别复杂.所以没来得及把公式亲自编辑一遍.而是在书上进 行截图.还请老师见谅！下次大作业一定做得更好！6.参考文献1 Jin

11、 J Theoiy and Confutation of Electromagnetic FieldsM/ Theory and confutation of electromagnetic fields Wiley , 20102 佚名.电磁场数值计算MJ.高等教育出版社,1996阿红雨.电磁场数值计算法与MATLAB实现J 2004.4 Jin J.电磁场有限元方法M西安电子科技大学出版社,19985 王长淸.现代计算电磁学基础卩.2005陈涌频，孟敏,方宙奇电磁场数值方法M科学出版|, 2016.程序附录fimction FEM (Im ax)%输入参数纵向最大网格数的1/4最好在20-

12、30.能满足精度且计算时间不会过长global ndm nel ne%全局常量总节点数.nel基兀数，i】e表不边界上的节点数Jmax=4*Imax;%Jmax表示横向纵向网格数n din=(Jm ax+1 )*(Jm a x+1);%ndm表示总点数Vl=l; V2=0;%输入边界电位值epsl =9; eps2=l; rhol =0;rho2=0; %输入上下两部分的和p值dx=l/Jmax; dy=l/Jniax;%横向纵向求网格间距X=l:ndin;Y=l:ndin;NN=zeios(Jmax+l Jmax+1);%节 点编号|xx| 格咅 ij 分111=0;forj=l:Jmax+

13、lfor i=l:Jmax+l111=111+1;NN(ij)=nl;X(nl)=(i-l)*dx; Y(nl)=(j-l)*dy; endendNE=zeros(3 2 *nclm); nl=0;forj=l:Jmaxfor i=2:Jmax+l*%X=i列，Y=j行处节点%求节点横纵坐标%描述各单元局部节点编号与总体编号对应的矩阵%各单元局部节点编号与总体编号对应nl=nl+l;NE(luil)=NN(ij); NE(2jil)=NN(i-lj+l); NE(3ail)=NN(i-lj); nl=nl+l;NE(14il)=NN(ij);NE(2jil)=NN(ij+l);NE(34il)

14、=NN(i-lj+l); endend%每个小正方形分为两个三角形网格.三 角形区域上下两部分节点坐标分别求nel=nl;%总网格数*K=zeros(n dm .ndin);%定义K矩阵Ke=zeros(33)；%单元Ke矩阵s=0.5/(Jmax*Jmax);%单元面积b=zeros(ndm,l);%b矩阵be=l:3; eps=l:nel;%单元be矩阵rho=l :nel;%定义和P数组for n=l :2* Jmax*Iinax%定义上下两部分的和p值，两部分的分别 为9和1. p都为0eps(n)=epsl; rlio(n)=rhol; endfor n=2* Jinax*Iinax+l :nel*水*木*水水*水*eps(n)=eps2; rho(n)=iiio2; end%*计算系统的Kb矩阵for n=l:nelfor i=l:3nl=NE