matlaB程序有限元法解泊松方程.doc

基于matlaB编程的有限元法 一、待求问题： 泛定方程： 边界条件：以（0,-1），（0,1），（1,0）为顶点的三角形区域边界上 编程思路及方法 给节点和三角形编号，并设定节点坐标 采用上述剖分方法，节点位置也比较规则。然后利用循环从区域内部（非边界）的节点开始编号，格式为NN(i,j)=n1，i，j分别表示节点所在列数与行数，并将节点坐标存入相应矩阵X(n1)，Y(n1)。由于区域上下两部分形状不同因此，分两个循环分别编号赋值，然后再对边界节点编号赋值。 然后再每个单元的节点进行局部编号，由于求解区域和剖分单元的特殊性，分别对内部节点对应左上角正方形的两个三角形单元，上左，左上，下斜边界节点要对应三个单元，上左，左上，左下，右顶点的左下、左上，右上边界的左上，分别编号以保证覆盖整个区域。 求解 首先一次获得每个单元节点的整体编号，然后根据其坐标求出每个三角形单元的面积。利用有限元方法的原理，分别求出系数矩阵和右端项，并且由于边界条件特殊，边界上，因此做积分时只需对场域单元积分而不必对边界单元积分。求的两个矩阵后很容易得到节点电位向量，即泊松方程的解。 画 2、f=1，x方向75份，y方向150份时，解函数平面图和曲面图 对比：当f=1时，界函数平面图 3、输出文本文件由于节点多较大，列在本文最末 四、结果简析 由于三角形区域分布的是正电荷，因此必定电位最高点在区域中部，且沿x轴对称，三角形边界电位最低等于零。当f=1时，发现电位最高点向x轴负方向移动了，这是由于此时电荷在三角形区域上均匀分布，而f=x时，x越大面电荷密度越大，附近相应电位越高，所得图像与实际情况相符。 五、MatlaB源程序 1、Finite_element_tri.m文件 function Finite_element_tri(Imax) % 用有限元法求解三角形形区域上的Possion方程,其中a为1，f=x Jmax=2*Imax; % 其中Imax Jmax分别表示x轴和y轴方向的网格数，其中Jmax等于Imax的两倍 % 定义一些全局变量 global ndm nel na % ndm 总节点数 % nel 基元数 % na 表示活动节点数 V=0; J=0;X0=1/Imax;Y0=X0;%V=0为边界条件 domain_tri % 调用函数画求解区域 [X,Y,NN,NE]=setelm_tri(Imax,Jmax); % 给节点和三角形元素编号，并设定节点坐标 % 以下求解有限元方程的求系数矩阵 T=zeros(ndm,ndm); for n=1:nel n1=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;%三角形面积 for k=1:3 if n1=na|n2=na T(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(n1,n2); T(n1,n1)=T(n1,n1)+((Y(n2)-Y(n3))^2+(X(n3)-X(n2))^2)/(4*s);%V=0则边界积分为零，非零时积分编程类似，再加边界积分。 end k=n1;n1=n2;n2=n3;n3=k; % 轮换坐标将值赋入3阶主子矩阵中 end end M=T(1:na,1:na); % 求有限元方程的右端项 f=X; G=zeros(na,1); for n=1:nel n1=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; for k=1:3 if n1=na G(n1)=G(n1)+(2*f(n1)+f(n2)+f(n3))*s/12;%f在单元上为线性差值时场域单元的积分公式 end n4=n1; n1=n2; n2=n3; n3=n4; % 轮换坐标标 end end F=M\G; % 求解方程得结果 NNV=zeros(Imax+1