电磁场有限元Matlab解法

电磁场有限元Matlab数值解法摘要有限元法(FEM)在电磁场的问题中十分有效，能帮我们解决很多电磁场的问题。本文主要对电磁场的有限元方(FEM)法进行学习和研究，利用Matlab进行编程仿真，对特定问题进行求解，模拟仿真出在一定边界条件下的电势分布图。该方法十分方便有效，能直观的得到我们想要的电势分布图。0引言有限元的原理在数学上首先由R.Courant于1943年提出，早期在力学中用于结构分析。五十年代初期，由于工程分析的需要，有限元法在复杂的航空结构分析中最先的到应用，而有限元法(FiniteElementMethod)这一名称在20世纪60年代由R.W.Clough的首先提出。三十多年来，以变分原理为基础建立起来的有限元法，因其理论依据的普遍性，不仅被广泛地应用于各种结构工程，而且作为一种声誉很高的数值分析方法已被普遍推广并成功地用来解决其他工程领域中的问题，例如热传导、流体力学、空气动力学、机械零件强度分析等。1968年有限元法开始用于求解电磁场的问题，有限元方法在电磁场中的应用至今已有近50年的历史。随着科学技术的发展，计算机性能得到提高，算法的不断被优化，有限元理论及技术在电磁应用方面已愈加成熟并取得许多研究成果。有限元法作为一种强有力的工程分析方法被广泛地应用于各种研究领域，有限元法同样是用于各类电磁场、电磁波工程问题定量分析与优化设计的最主要的数值方法，并且无一例外地是构成各种先进、有效的计算软件包的基础。目前，有限元分析已成为计算机辅助设计的一个重要组成部分。1.有限元简介有限元法的基本思想是将结构离散化，用有限个容易分析的单元来表示复杂的对象，单元之间通过有限个节点相互连接，然后根据边界条件综合求解。由于单元的数目是有限的，节点的数目也是有限的，所以称为有限元法(FEM，FiniteElementMethod)。有限元法是以微分方程为基础的数值方法。最初主要是利用变分原理将微分方程变为等价的变分方程，经改进的Ritz法，则将微分方程的求解变为代数方程的求解问题。由于并非任意的微分方程都能找到与其等价的变分方程，使得上述形式的有限元法的应用受到了一定的限制。用加权余量的伽辽金法直接从微分方程出发构造有限元方程，就可以突破变分原理的限制。有限元最大的特点是：先通过各种适当的形式将求解域划分成有限个单元，再在每个单元中构造分域基函数，利用Ritz法或伽辽金法构造代数形式的有限元方程。有限元法的最大优点是其离散单元的灵活性。相对而言，有限元法可以更精确地模拟各种复杂的几何结构，并通过选取样点的疏密情况来适应场分布的不同情况，既能保证计算精度的要求，又不增加过大的计算量。另一优点是所形成的有限元方程组的系数矩阵是稀疏的对称阵，这非常有利于代数方程组的求解。与其他数值计算方法相比较，有限元法在适应场域边界几何形状以及媒质物理性质变异情况复杂的问题求解上有突出的优点，即方法应用不受边界形状和媒质性质的限制，而且不同媒质分界面上的边界条件是自动满足的，第二三类边界条件不必作单独处理。此外离散点配置比较随意，并且取决于有限单元剖分密度和单元插值函数的选取，可以获得令人满意的数值计算精度。有限元法还可以方便地编写通用计算程序，使之构成模块化的子程序集合，适应计算功能延拓的需要，从而构成各种高效能的计算软件包。有限元法的发展与应用前景令人瞩目。2.有限元问题分析所需求解的电场区域如图2-1所示。该图描述的为屏蔽微带传输线结构。假设边界上的即屏蔽导体上电势U=0,内导体即条带上电势U=1V，内部上部分介质相对介电常数ε=1，电荷密度ρ=0，下部分介质相对介电常数ε=9，电荷密度ρ=0。由于其几何结构是对称的，我们取其一半进行研究，如图2-2所示。其中电势在左边边界上由于对称分布在法向方向导数为0。图2-1求解区域图2-2等效求解区域设求解区域长宽均为1，而条带处于1/4位置且长度为1/4。电势满足如下的方程：0=DonD()0ˆonnNon把求解域按如图2-3所示方法进行剖分并对节点和单元进行编号。图2-3网格划分每个区域的电势可以用来表示。把三角形的三个顶点带入可得解上述三个方程得到的a，b，c带入中得：为插值函数其中为三角形面积。对各个单元的电位值进行求和，则整个求解区域的电位值可表示为对方程乘权函数wi并在整个区域积分得0利用矢量恒等式和高斯定理得带入边界条件得使用伽辽金方法取权函数带入上式我们得可以写为，可描述为其中，Kij可表示为各个单元贡献的总和：每个单元的，ρ=0。以上为理论上的分析。然后利用Matlab进行编程仿真。仿真程序见程序附录。3.仿真结果及分析仿真结果如图所示，分别为：4-1平面电位分布图；4-2电位分布三维立体图；4-3电位等势线图。图4-1平面电位分布图图4-2电位分布三维立体图图4-3电位等势线图从上面几幅图中我们可以直观的看到电势在所给的求解区域中的分布，生动形象。且边界处的条件吻合得很好。例如上下边和右边电势都为0，中心处电势为1，左边界处可以看到等势线垂直于边界，则电场法向分量为零，也与边界条件很符合。通过有限元的方法，我们用Matlab编程仿真出了静电场的电场分布图。5.总结通过此次有限元法的研究和学习，成功地仿真出电场的分布图，收获很大。在刚开始接触有限元法时，完全没有思路，对Matlab软件也一窍不通。后来通过大量的文献查阅和学习，逐渐理解了有限元法的基本思路和编程方法。在仿真的过程中，也遇到过各种问题，挑战极大，后来通过反复看文献理解有限元的思想，学习Matlab的编程方法才慢慢有了一定的体会，最终才得到了较为正确的结果，也学到了对其他电场区域仿真的基本思路。这次的仿真，不仅仅是一次完成作业，更是一次挑战，不仅学到了课本上的知识，还学到了Matlab的编程技巧，收获远不止一次完成任务这么简单。虽然仿真结果基本吻合，但是此次仿真还有很多不足之处。首先是只仿真出了电势的分布图，而没有电场分布图。然后是程序算法不够简单，当剖分的网格数较多时计算速度很慢。还有就是网格的划分算法不够好，应该尽量划分为正三角形，来减小误差，以后还要多学习这方面的知识。最后，由于时间比较紧，理论推到中的公式又特别复杂，所以没来得及把公式亲自编辑一遍，而是在书上进行截图，还请老师见谅！下次大作业一定做得更好！6.参考文献[1]JinJ.TheoryandComputationofElectromagneticFields[M]//Theoryandcomputationofelectromagneticfields.Wiley:,2010.[2]佚名.电磁场数值计算[M].高等教育出版社,1996.[3]何红雨.电磁场数值计算法与MATLAB实现[J].2004.[4]JinJ.电磁场有限元方法[M].西安电子科技大学出版社,1998.[5]王长清.现代计算电磁学基础[J].2005.[6]陈涌频,孟敏,方宙奇.电磁场数值方法[M].科学出版社,2016.程序附录functionFEM(Imax)%输入参数纵向最大网格数的1/4，最好在20~30，能满足精度且计算时间不会过长globalndmnelne%全局常量ndm总节点数，nel基元数，ne表示边界上的节点数Jmax=4*Imax;%Jmax表示横向纵向网格数ndm=(Jmax+1)*(Jmax+1);%ndm表示总点数V1=1;V2=0;%输入边界电位值eps1=9;eps2=1;rho1=0;rho2=0;%输入上下两部分的ε和ρ值dx=1/Jmax;dy=1/Jmax;%横向纵向求网格间距X=1:ndm;Y=1:ndm;NN=zeros(Jmax+1,Jmax+1);%节点编号%************************网格剖分**********************************n1=0;forj=1:Jmax+1fori=1:Jmax+1n1=n1+1;NN(i,j)=n1;%X=i列，Y=j行处节点X(n1)=(i-1)*dx;Y(n1)=(j-1)*dy;%求节点横纵坐标endendNE=zeros(3,2*ndm);%描述各单元局部节点编号与总体编号对应的矩阵n1=0;forj=1:Jmaxfori=2:Jmax+1n1=n1+1;NE(1,n1)=NN(i,j);%各单元局部节点编号与总体编号对应NE(2,n1)=NN(i-1,j+1);NE(3,n1)=NN(i-1,j);n1=n1+1;NE(1,n1)=NN(i,j);NE(2,n1)=NN(i,j+1);NE(3,n1)=NN(i-1,j+1);endend%每个小正方形分为两个三角形网格，三角形区域上下两部分节点坐标分别求nel=n1;%总网格数%******************定义各个单元的常量和矩阵************************K=zeros(ndm,ndm);%定义K矩阵Ke=zeros(3,3);%单元Ke矩阵s=0.5/(Jmax*Jmax);%单元面积b=zeros(ndm,1);%b矩阵be=1:3;%单元be矩阵eps=1:nel;rho=1:nel;%定义ε和ρ数组forn=1:2*Jmax*Imax%定义上下两部分的ε和ρ值，,两部分的ε分别为9和1，ρ都为0eps(n)=eps1;rho(n)=rho1;endforn=2*Jmax*Imax+1:neleps(n)=eps2;rho(n)=rho2;end%****************计算系统的[K][b]矩阵*************************forn=1:nelfori=1:3n1=NE(1,n);n2=NE(2,n);n3=NE(3,n);%给每个单元的点进行编号bn(1)=Y(n2)-Y(n3);bn(2)=Y(n3)-Y(n1);bn(3)=Y(n1)-Y(n2);cn(1)=X(n3)-X(n2);cn(2)=X(n1)-X(n3);cn(3)=X(n2)-X(n1);forj=1:3Ke(i,j)=eps(n)*(bn(i)*bn(j)+cn(i)*cn(j))/(4*s);be(i)=s*rho(n)/3;%计算每个单元的Ke和be矩阵endendfori=1:3forj=1:3K(NE(i,n),NE(j,n))=K(NE(i,n),NE(j,n))+Ke(i,j);b(NE(i,n))=b(NE(i,n))+be(i);%把Ke和be分别相加求总矩阵endendend%*****************边值问题的处理，得到最终电位分布*****************n1=0;nd=1:3*Jmax+Imax+2;%定义描述边界节点编号与总体编号对应的数组p=1:3*Jmax+Imax+2;%定义边值数组fori=1:Jmax+1;%定义各边的边值n1=n1+1;nd(n1)=NN(i,1);p(n1)=V2;endforj=2:Jmax+1;n1=n1+1;nd(n1)=NN(Jmax+1,j);p(n1)=V2;endfori=1+Jmax:-1:1n1=n1+1;nd(n1)=NN(i,Jmax+1);p(n1)=V2;endfori=1:Imax+1n1=n1+1;nd(n1)=NN(i,Imax+1);p(n1)=V1;end%对边界上的节点进行赋值和单独编号ne=n1;%边界上的节点总数fori=1:neb(nd(i))=p(i);K(nd(i),nd(i))=1;forj=1:ndmifj~=nd(i)b(j)=b(j)-K(j,nd(i))*p(i);K(j,nd(i))=0;K(nd(i),j)=0;endendend%把边值带入进行处理计算新的K和b矩阵U=K\b;%根据K和b求出各点的电位值%****************************电势分布图***************************NF=zeros(Jmax+1,Jmax