发布时间 : 星期一 文章matlaB程序的有限元法解泊松方程更新完毕开始阅读176f49d649649b6648d747f4
基于matlaB编程的有限元法
一、待求问题:
泛定方程:???=x
边界条件:以(0,-1),(0,1),(1,0)为顶点的三角形区域边界上?=0
2二、编程思路及方法
1、给节点和三角形单元编号,并设定节点坐标
画出以(0,-1),(0,1),(1,0)为顶点的三角形区域figure1
由于积分区域规则,故采用特殊剖分单元,将区域沿水平竖直方向分等份,此时所有单元都是等腰直角三角形,剖分单元个数由自己输入,但竖直方向份数(用Jmax表示)必须是水平方向份数(Imax)的两倍,所以用户只需输入水平方向的份数Imax。
采用上述剖分方法,节点位置也比较规则。然后利用循环从区域内部(非边界)的节点开始编号,格式为NN(i,j)=n1,i,j分别表示节点所在列数与行数,并将节点坐标存入相应矩阵X(n1),Y(n1)。由于区域上下两部分形状不同因此,
分两个循环分别编号赋值,然后再对边界节点编号赋值。
然后再每个单元的节点进行局部编号,由于求解区域和剖分单元的特殊性,分别对内部节点对应左上角正方形的两个三角形单元,上左,左上,下斜边界节点要对应三个单元,上左,左上,左下,右顶点的左下、左上,右上边界的左上,分别编号以保证覆盖整个区域。
2、求解泊松方程
首先一次获得每个单元节点的整体编号,然后根据其坐标求出每个三角形单元的面积。利用有限元方法的原理,分别求出系数矩阵和右端项,并且由于边界条件特殊,边界上?=0,因此做积分时只需对场域单元积分而不必对边界单元积分。求的两个矩阵后很容易得到节点电位向量,即泊松方程的解。
3、画解函数的平面图和曲面图
由节点单位向量得到,j行i列节点的电位,然后调用绘图函数imagesc(NNV)与surf(X1,Y1,NNV')分别得到解函数的平面图figure2和曲面图figure3。
4、将结果输出为文本文件
输出节点编号,坐标,电位值
三、计算结果
1、积分区域:
10.80.60.40.20-0.2-0.4-0.6-0.8-100.10.20.30.40.50.60.70.80.91
2、f=1,x方向75份,y方向150份时,解函数平面图和曲面图
10203040500.030.0250.020.0150.01600.0057020406080100120140
0.030.040.0250.030.020.020.0150.010.01021.510.5000.40.20.60.810.005
对比:当f=1时,界函数平面图
0.071020300.04400.03500.02600.0170204060801001201400.060.05
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 % 调用函数画求解区域