二维三角形线性插值有限元法的原理与应用

二维三角形单元的有限元程序编写及应用葛青宇 1182201150 研电1803班一、二维有限元理论介绍1、电磁场边值问题电磁场边值问题一般满足不同的微分方程和相应的边界条件。

泊松方程是电磁场位函数满足方程中的一大类，如下所示：2a u f -∇=(1)10uu Γ=(2)2,3u abu c nΓ∂=+∂(3)其中：（1）式称为泊松方程，u 为电磁场的位函数，a 为表征空间介质材料的常数，f 表示 产生电磁场位函数的源；（2）式表示第一类边界条件， u 0为边界上已知的位函数；(3)式为第三类边界条件，b 和 c 为已知的边界条件系数。

2、二维伽辽金加权余量方程以插值方式构造近似函数，u 的近似解表示为1nn n n n u M u ==∑(4)其中，基函数为12,, ,n n M M M ⋅⋅⋅，相应待定常数为 12,, ,n n u u u ⋅⋅⋅ 根据基函数和权函数与单元形状函数的关系，在有限元法中一个单元内基函数和权函数的值等于相应的形状函数的值，表示为在单元e 中，,e i n i M N =，,e j m j M N =。

采用伽辽金加权余量形式，设一组与基函数相同的权函数12,, ,n n M M M ⋅⋅⋅，将权函数代入伽辽金加权余量方程，得一组方程2()d d m m M a u M f ΩΩ-∇Ω=Ω⎰⎰(5)其中，m =1，2，…，n n , 对应节点的总体编号。

对（4）式应用格林定理，并将第二、三类边界条件代入得d ()d mm m a Mu M bu c d M f •ΩΓΩ∇∇Ω-+Γ=Ω⎰⎰⎰(6)其中m 同上。

将不含未知函数u 的项移至方程右端，得d d d d mm m m a Mu bM u M f cM •ΩΓΩΓ∇∇Ω-Γ=Ω+Γ⎰⎰⎰⎰(7)将由基函数线性组合构成的近似函数代入方程，得 11[()]d [()d ()d ()d nnn n m n nmn nmm n n a M M u bM M u Mf cM •==ΩΓΩΓ∇∇Ω-Γ=Ω+Γ∑∑⎰⎰⎰⎰ (8)将近似函数的梯度运算表示出来，得11{[()]}d [()]d ()d ()d nnn n mnn m n n n n m m a M Mu bM M u M f cM •==ΩΓΩΓ∇∇Ω-Γ=Ω+Γ∑∑⎰⎰⎰⎰ (9)权函数与求和无关，可以将其移到求和号之内，得11{[()]}d [()]d ()d ()d nnn n mn n m n n n n m m a MM u b M M u M f cM •==ΩΓΩΓ∇∇Ω-Γ=Ω+Γ∑∑⎰⎰⎰⎰ (10)将未知节点变量从积分中提出得：11()d ()d ()d ()d nnn n mn n m n n m m n n a MM u bM M u M f cM •==ΩΓΩΓ∇∇Ω-Γ=Ω+Γ∑∑⎰⎰⎰⎰ (11)可以看出代数方程如下所示：,11,22,,+n n m m m n n m n n m K u K u K u K u R ++++=(12)对比式（10），可得： ,()d ()d m n m n m n K a M M bM M •ΩΓ=∇∇Ω-Γ⎰⎰(13)()d ()d m m m R M f cM ΩΓ=Ω+Γ⎰⎰(14)将（12）、（13）式化作单元上积分之和得形式，得：,11()d ()d eebbn n m n m n m n e e EE K a M M bM M •===∇∇Ω-Γ∑∑⎰⎰(15)11()d ()d e ebbn n m m m e e EE R M f cM ===Ω+Γ∑∑⎰⎰(16)n e 和n eb 分别是场域单元数和边界单元数，E 和E b 分别是单个场域单元积分区域和单个边界单元积分区域。

将基函数和权函数更换为单元形状函数，在单元上将基函数和权函数用双层下标标记：,,,,,11()d ()d e ebe i e j e i e j bn n m n m n m n e e EE K a M M bM M •===∇∇Ω-Γ∑∑⎰⎰(17),,,,,11()d ()d e ebe i e j e i e j bn n m n m n m n e e EE K a M M bM M •===∇∇Ω-Γ∑∑⎰⎰(18)给定m 和n 的值，根据单元节点局部编号与整体编号的对应关系，寻找单元上的i 和j 。

单元形状函数如下所示：,11()d ()d eebbn n m n i j i j e e EE K a N N bN N •===∇∇Ω-Γ∑∑⎰⎰(19)11()d ()d e ebbn n m i i e e EE R N f cN ===Ω+Γ∑∑⎰⎰(20)得到得方程组以矩阵形式表示为： =Ku R (21)3、二维三角形有限元矩阵计算在二维三角形单元中（1）场域单元系数矩阵K e,i,j ：[],,12122121212122121()d ()()(2)()()()()4e i j i j Ej j i i i i j j i i j j i i j j K a N N y y a y y x x x x a y y y y x x x x •++++++++++++++++=∇∇Ω-⎡⎤--⎢⎥-⎣⎦=∆∆=--+--∆⎰ (22)其中，Δ为三角形单元面积。

单元节点编号下标采用1到3循环计数。

（19）式中当计算出的下标大于 3 时， 取计算值减 3 为实际下标。

（2）边界单元系数矩阵K eb,i,j 为()1,,1d d eeb i j i j e i j K bN N N N ξΓ-=-Γ=-⎰(23)令l =116161eb bl ⎡⎤=-⎢⎥⎣⎦K(24)（3）场域单元右端项R e,i(),d e i i ER N f =Ω⎰(25)若源在单元上为常数，根据三角形单元解析积分公式：1131e f ⎡⎤∆⎢⎥=⎢⎥⎢⎥⎣⎦R (26)否则根据源在单元节点上的函数值进行插值，这样场域右端项为：,,,21112112112e i e j e k n e n n f f f ⎡⎤⎡⎤⎢⎥∆⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦R(27)（4）边界右端项R eb,i ,可以表示为(),d beb i i E R cN =Γ⎰(28)若c 在单元上为常数，则112eb cl ⎡⎤=⎢⎥⎣⎦R (29)若在单元上，c 为线性插值函数，则,,21126eb i sb j n eb n c l c ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦R (30)二、二维三角形单元有限元程序编写及验证1、程序思路图 1 二维三角形单元有限元法流程图2、程序验证 （1）算例1采用PDE 工具箱给出的一个例子，对程序进行验证。

两个圆形的金属导体被放置在一张浸过盐水的吸墨纸上，吸墨纸起着平面薄导体的作用。

这个问题的物理模型由拉普拉斯方程组成。

求解场域以及边界条件如图 2 模型示意图图 2与式(31)所示。

图 2 模型示意图22222=01.2,0.60.6=00.6, 1.2 1.21.2,0.60.60.6, 1.2 1.21,(0.6)0.091,(0.6)0.09u n u x y y x x y y x u x y u x y σσ⎧⎪∇⎪⎪⎧⎪⎪=--≤≤⎪⎪⎪⎪=-≤≤⎨⎨⎪⎪=-≤≤⎪⎪=--≤≤⎪⎪⎩⎪=++=⎪⎪=--+=⎩∂∂， (31)这里令 σ=1。

通过PDE 工具箱进行网格剖分，共得到512个单元，305个节点，如图 3所示。

图 3 网格剖分结果使用PDE工具箱求解此问题，得到电压等势图及电流场，如图4所示。

图 4 PDE工具箱求得的电位及电流密度导出剖分所得单元、节点信息。

使用自编程序进行求解，得到结果如下。

图 5 自编程序求得的电位及电流密度比较图 4 PDE 工具箱求得的电位及电流密度图 4与图 5可知，自编程序得到了PDE 工具箱给出的结果。

图上些许的差异，特别是等势线与矢量场的不同，是由于自己绘图方式与PDE 工具箱绘图方式不同所导致。

进一步比较得到的电位值，最大相对误差为4.9982e-10，可说明自编程序的正确性。

（2）算例2在单位圆盘内求解泊松方程：022222=1,110u y u x y x ερρ⎧∇⎪⎪⎨⎪⎪=+=⎩=+ (32)网格剖分如图 6所示。

图 6 单位圆盘上的网格剖分通过PDE工具箱求解，得到电位分布如图7所示。

图7 PDE工具箱得到的单位圆盘上的电位分布导入网格后，用自编程序求解，得到结果如图8所示。

图8 自编程序得到的单位圆盘上的电位分布、等势线与电场分布从图7图8可以看到，对于电荷密度在求解域变化的情况，自编程序与PDE工具箱的结果出现了偏差，最大相对误差达到了0.9966。

三、结语本文首先叙述二维伽辽金加权余量方程的基本原理、以及基于该方程的有限元矩阵方程。

在ansys软件剖分的基础上，根据伽辽金理论和有限元理论，编写matlab程序，计算边界和场域单元系数矩阵，并对剖分所得各节点电位进行求解，作出等位线图。

在编写过程中，遇到了一些问题。

特别是课本115页第6项中的式（1），它在书中放在了6个公式的最前面，但在程序中，式（1）的放在式（5）的后面，也就是最后，否则会是的单元系数矩阵亏秩，造成无法求解。

周遭同学大多结合Ansys软件来剖分网格与验证对比，但由于Ansys软件难于使用，我决定采用熟悉的PDE工具箱来替代Ansys。

在求解后的结果图绘制过程中，为了尽可能的美观，查阅了matlab对图形对象控制的文档资料，了解了一些深层次的图像属性，收获颇丰。

本程序参考了白婉欣师姐的代码，借鉴了她对intersec函数的使用，同时结合课本给出的公式，自主完成了本程序。

虽然本程序无法媲美PDE工具箱等成熟有限元计算软件，但通过此次编程实践，我了解了有限元的原理以及编程实现，对我的认知与编程能力有了不小的提升。

参考文献[1] 王泽忠. 简明电磁场数值计算[M]. 机械工业出版社, 2011.[2] 白婉欣. 二维三角形单元的有限元法程序编写及应用.附录二维三角形单元有限元matlab程序%% My FEM Programme%% input mesh data% 由pde工具箱导入得到，三个矩阵:p,e,t,帮助条目：Mesh Data as [p,e,t] Triples 给出的解释% 翻译一下，就是% p:[横坐标；纵坐标]% e:[点1编号；点2编号；...(第5行)几何边界编号(确定边界条件类型);...]%边界单元% t:[点1编号；点2编号；点3编号;子域编号]%三角形单元，针对线性插值% 边界条件% bc:[类型；参数值1;参数值2]1表示第1类条件，g,q;2表示第2类条件，h,r% load data1 %包含example1网格信息p e t,方程系数f a,边界条件bc以及求解结果u% load data2 %example2%% 泊松问题% a=1;%本题目中σ% f=0;%%下列参数在边界条件处理时从变量中读取% b=0;% c=0;%% 预定义K=zeros(length(p(1,:)));R=zeros(length(p(1,:)),1);%非边界单元for i=1:length(t(1,:))%提取三个节点node(1)=t(1,i);node(2)=t(2,i);node(3)=t(3,i);x(1)=p(1,node(1));y(1)=p(2,node(1));x(2)=p(1,node(2));y(2)=p(2,node(2));x(3)=p(1,node(3));y(3)=p(2,node(3));%先算三角形面积delta=0.5*((x(1)-x(3))*(y(2)-y(3))-(x(2)-x(3))*(y(1)-y(3)));%计算系数矩阵fun1=@(x) (x+1)-3*(ceil((x+1)/3)-1);%定义一个匿名函数，用于计算节点编号下标fun2=@(x) (x+2)-3*(ceil((x+2)/3)-1);for ii=1:3for jj=1:3yProduct=(y(fun1(ii))-y(fun2(ii)))*(y(fun1(jj))-y(fun2(jj)));%纵坐标差值的乘积xProduct=(x(fun2(ii))-x(fun1(ii)))*(x(fun2(jj))-x(fun1(jj)));%横坐标差值的乘积K(node(ii),node(jj))=K(node(ii),node(jj))+a/(4*delta)*(yProduct+xProduct);endend%场域单元有段向量RA=ones(3,3)+eye(3);if length(f)>1F=f([node(1),node(2),node(3)]);elseF=f*ones(1,3);endB=R([node(1),node(2),node(3)]).'+delta/12*A*F.';R(node(1))=B(1);R(node(2))=B(2);R(node(3))=B(3);% R(node(1))=R(node(1))+delta*f/3;% R(node(2))=R(node(2))+delta*f/3;% R(node(3))=R(node(3))+delta*f/3;end%边界单元for i=1:length(e(1,:))node(1)=e(1,i);node(2)=e(2,i);boundary=e(5,i);%所在的边界编号if bc(1,boundary)==2%第二、第三类边界条件b=-bc(2,boundary);%将pdeTool的公式转化为书上的参数c=bc(3,boundary);%单元系数矩阵x(1)=p(1,node(1));y(1)=p(2,node(1));x(2)=p(1,node(2));y(2)=p(2,node(2));%长度l=sqrt( (x(1)-x(2))^2+(y(2)-y(1))^2);K(node(1),node(1))=K(node(1),node(1))+b*l/3;K(node(1),node(2))=K(node(1),node(2))+b*l/6;K(node(2),node(1))=K(node(1),node(2))+b*l/6;K(node(2),node(2))=K(node(2),node(2))+b*l/3;% 右端向量RR(node(1))=R(node(1))+c*l/2;R(node(2))=R(node(2))+c*l/2;else%第一类边界条件u0=bc(3,boundary)/bc(2,boundary);%求得给定的位函数值R=R-K(:,node(1))*u0;%（2）R(node(1))=u0;%（3）K(node(1),:)=0;%（4）K(:,node(1))=0;%（5）K(node(1),node(1))=1;%(1)%注意课本上115页这一语句在语句（2）的前面...%这样会使K矩阵对角线元素为零，使其不满秩，无法求得正确的结果%猜测边界上任何一点都会作为边界元的开头，不过为保险起见，先加入末端点%经过验证，确实不用下面这些% R=R-K(:,node(2))*u0;% R(node(2))=u0;% K(node(2),:)=0;% K(:,node(2))=0;% K(node(2),node(2))=1;endend%% 求解u1=K\R;%% 最大误差error=max(abs(u-u1)./u);%% 绘图figure(1)x=linspace(min(p(1,:)),max(p(1,:)),20);y=linspace(min(p(2,:)),max(p(2,:)),20);[xx,yy]=meshgrid(x,y);uxy=tri2grid(p,t,u1,x,y);[dx, dy]=gradient(uxy);% 绘制彩色三角形单元[X,Y]=meshgrid(p(1,:),p(2,:));trisurf(t(1:3,:)',p(1,:),p(2,:),u1,'EdgeColor','none');view(2);axis equalcolormap('cool')hold on% 绘制等势线与电场线contourf(xx,yy,uxy,'ContourZLevel', max(u1)+1);quiver3(xx,yy,(max(u1)+1)*ones(size(xx)),-dx,-dy,zeros(size(dx)),'r'); hold off。