资源描述：
《matlaB程序的有限元法解泊松方程.doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、基于matlaB编程的有限元法一、待求问题：泛定方程：边界条件：以（0,-1），（0,1），（1,0）为顶点的三角形区域边界上二、编程思路及方法1、给节点和三角形单元编号，并设定节点坐标画出以（0,-1），（0,1），（1,0）为顶点的三角形区域figure1由于积分区域规则，故采用特殊剖分单元，将区域沿水平竖直方向分等份，此时所有单元都是等腰直角三角形，剖分单元个数由自己输入，但竖直方向份数（用Jmax表示）必须是水平方向份数（Imax）的两倍，所以用户只需输入水平方向的份数Imax。采用上述剖分方法，节点位置也比
2、较规则。然后利用循环从区域内部（非边界）的节点开始编号，格式为NN(i,j)=n1，i，j分别表示节点所在列数与行数，并将节点坐标存入相应矩阵X(n1)，Y(n1)。由于区域上下两部分形状不同因此，分两个循环分别编号赋值，然后再对边界节点编号赋值。然后再每个单元的节点进行局部编号，由于求解区域和剖分单元的特殊性，分别对内部节点对应左上角正方形的两个三角形单元，上左，左上，下斜边界节点要对应三个单元，上左，左上，左下，右顶点的左下、左上，右上边界的左上，分别编号以保证覆盖整个区域。1、求解泊松方程首先一次获得每个单元节
3、点的整体编号，然后根据其坐标求出每个三角形单元的面积。利用有限元方法的原理，分别求出系数矩阵和右端项，并且由于边界条件特殊，边界上，因此做积分时只需对场域单元积分而不必对边界单元积分。求的两个矩阵后很容易得到节点电位向量，即泊松方程的解。2、画解函数的平面图和曲面图由节点单位向量得到，j行i列节点的电位，然后调用绘图函数imagesc(NNV)与surf(X1,Y1,NNV')分别得到解函数的平面图figure2和曲面图figure3。3、将结果输出为文本文件输出节点编号，坐标，电位值三、计算结果1、积分区域：2、f
4、=1，x方向75份，y方向150份时，解函数平面图和曲面图对比：当f=1时，界函数平面图3、输出文本文件由于节点多较大，列在本文最末四、结果简析由于三角形区域分布的是正电荷，因此必定电位最高点在区域中部，且沿x轴对称，三角形边界电位最低等于零。当f=1时，发现电位最高点向x轴负方向移动了，这是由于此时电荷在三角形区域上均匀分布，而f=x时，x越大面电荷密度越大，附近相应电位越高，所得图像与实际情况相符。五、MatlaB源程序1、Finite_element_tri.m文件functionFinite_element_
5、tri(Imax)%用有限元法求解三角形形区域上的Possion方程,其中a为1，f=xJmax=2*Imax;%其中ImaxJmax分别表示x轴和y轴方向的网格数，其中Jmax等于Imax的两倍%定义一些全局变量globalndmnelna%ndm总节点数%nel基元数%na表示活动节点数V=0;J=0;X0=1/Imax;Y0=X0;%V=0为边界条件domain_tri%调用函数画求解区域[X,Y,NN,NE]=setelm_tri(Imax,Jmax);%给节点和三角形元素编号，并设定节点坐标%以下求解有限元
6、方程的求系数矩阵T=zeros(ndm,ndm);forn=1:neln1=NE(1,n);n2=NE(2,n);n3=NE(3,n);%整体编号s=abs((X(n2)-X(n1))*(Y(n3)-Y(n1))-(X(n3)-X(n1))*(Y(n2)-Y(n1)))/2;%三角形面积fork=1:3ifn1<=na
7、n2<=naT(n1,n2)=T(n1,n2)+((Y(n2)-Y(n3))*(Y(n3)-Y(n1))+(X(n3)-X(n2))*(X(n1)-X(n3)))/(4*s);T(n2,n1)=T(n
8、1,n2);T(n1,n1)=T(n1,n1)+((Y(n2)-Y(n3))^2+(X(n3)-X(n2))^2)/(4*s);%V=0则边界积分为零，非零时积分编程类似，再加边界积分。endk=n1;n1=n2;n2=n3;n3=k;%轮换坐标将值赋入3阶主子矩阵中endendM=T(1:na,1:na);%求有限元方程的右端项f=X;%场源函数G=zeros(na,1);forn=1:neln1=NE(1,n);n2=NE(2,n);n3=NE(3,n);s=abs((X(n2)-X(n1))*(Y(n3)-Y(
9、n1))-(X(n3)-X(n1))*(Y(n2)-Y(n1)))/2;fork=1:3ifn1<=naG(n1)=G(n1)+(2*f(n1)+f(n2)+f(n3))*s/12;%f在单元上为线性差值时场域单元的积分公式endn4=n1;n1=n2;n2=n3;n3=n4;%轮换坐标标endendF=MG;%求解方程得结果NNV=zero
