电磁场有限元Matlab解法

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