matlaB程序的有限元法解泊松方程.docx
- 文档编号:13245521
- 上传时间:2023-06-12
- 格式:DOCX
- 页数:229
- 大小:167.85KB
matlaB程序的有限元法解泊松方程.docx
《matlaB程序的有限元法解泊松方程.docx》由会员分享,可在线阅读,更多相关《matlaB程序的有限元法解泊松方程.docx(229页珍藏版)》请在冰点文库上搜索。
matlaB程序的有限元法解泊松方程
基于matlaB编程的有限元法
一、待求问题:
泛定方程:
边界条件:
以(0,-1),(0,1),(1,0)为顶点的三角形区域边界上
二、编程思路及方法
1、给节点和三角形单元编号,并设定节点坐标
画出以(0,-1),(0,1),(1,0)为顶点的三角形区域figure1
由于积分区域规则,故采用特殊剖分单元,将区域沿水平竖直方向分等份,此时所有单元都是等腰直角三角形,剖分单元个数由自己输入,但竖直方向份数(用Jmax表示)必须是水平方向份数(Imax)的两倍,所以用户只需输入水平方向的份数Imax。
采用上述剖分方法,节点位置也比较规则。
然后利用循环从区域内部(非边界)的节点开始编号,格式为NN(i,j)=n1,i,j分别表示节点所在列数与行数,并将节点坐标存入相应矩阵X(n1),Y(n1)。
由于区域上下两部分形状不同因此,分两个循环分别编号赋值,然后再对边界节点编号赋值。
然后再每个单元的节点进行局部编号,由于求解区域和剖分单元的特殊性,分别对内部节点对应左上角正方形的两个三角形单元,上左,左上,下斜边界节点要对应三个单元,上左,左上,左下,右顶点的左下、左上,右上边界的左上,分别编号以保证覆盖整个区域。
2、求解泊松方程
首先一次获得每个单元节点的整体编号,然后根据其坐标求出每个三角形单元的面积。
利用有限元方法的原理,分别求出系数矩阵和右端项,并且由于边界条件特殊,边界上
,因此做积分时只需对场域单元积分而不必对边界单元积分。
求的两个矩阵后很容易得到节点电位向量,即泊松方程的解。
3、画解函数的平面图和曲面图
由节点单位向量得到,j行i列节点的电位,然后调用绘图函数imagesc(NNV)与surf(X1,Y1,NNV')分别得到解函数的平面图figure2和曲面图figure3。
4、将结果输出为文本文件
输出节点编号,坐标,电位值
三、计算结果
1、积分区域:
2、f=1,x方向75份,y方向150份时,解函数平面图和曲面图
对比:
当f=1时,界函数平面图
3、输出文本文件由于节点多较大,列在本文最末
四、结果简析
由于三角形区域分布的是正电荷,因此必定电位最高点在区域中部,且沿x轴对称,三角形边界电位最低等于零。
当f=1时,发现电位最高点向x轴负方向移动了,这是由于此时电荷在三角形区域上均匀分布,而f=x时,x越大面电荷密度越大,附近相应电位越高,所得图像与实际情况相符。
五、MatlaB源程序
1、Finite_element_tri.m文件
functionFinite_element_tri(Imax)
%用有限元法求解三角形形区域上的Possion方程,其中a为1,f=x
Jmax=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);%给节点和三角形元素编号,并设定节点坐标
%以下求解有限元方程的求系数矩阵
T=zeros(ndm,ndm);
forn=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;%三角形面积
fork=1:
3
ifn1<=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);
forn=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;
fork=1:
3
ifn1<=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,Jmax+1);
fi=zeros(ndm,1);
fi(1:
na)=F(1:
na);
fi(na+1:
ndm)=V;
forj=0:
Jmax
fori=0:
Imax
n=NN(i+1,j+1);
ifn<=0
n=na+1;
end
NNV(i+1,j+1)=fi(n);
end
end
figure
(2)
imagesc(NNV);%画解函数的平面图
X1=zeros(1,Imax+1);
Y1=zeros(1,Jmax+1);
fori=1:
Imax+1
X1(i)=(i-1)*X0;
end
fori=1:
Jmax+1
Y1(i)=(i-1)*Y0;
end
figure(3)
surf(X1,Y1,NNV');%画解函数的曲面图
%以下是结果的输出
fid=fopen('Finite_element_tri.txt','w');
fprintf(fid,'\n*********有限元法求解三角形区域上Possion方程的结果**********\n\n');
L=[1:
ndm]';
fprintf(fid,'\n\n节点编号坐标分量x坐标分量yu(x,y)的值\n\n');
fori=1:
ndm
fprintf(fid,'%8d%14.5f%14.5f%14.5f\n',L(i),X(i),Y(i),fi(i));
end
fclose(fid);
End
2、domain_tri.m文件
functiondomain_tri
%画求解区域
xy=[01;0-1;10];
A=zeros(3,3);
A(1,1)=2;A(1,2)=-1;A(1,3)=-1;
A(2,2)=2;A(2,1)=-1;A(2,3)=-1;
A(3,3)=2;A(3,2)=-1;A(3,1)=-1;
A=sparse(A);
figure
(1);
gplot(A,xy);
3、setelm_tri.m文件
function[X,Y,NN,NE]=setelm_tri(Imax,Jmax)
%给节点和三角形单元编号,并设定节点坐标
%定义一些全局变量
globalndmnelna
%I1I2J1J2ImaxJmax分别描述网线纵向和横向数目的变量
%XY表示节点坐标
%NN描述节点编号
%NE描述各单元局部节点编号与总体编号对应的矩阵
%ndm总节点数
%nel单元数
%na表示不含边界的节点数
nlm=Imax*Jmax;
dx=1/Imax;
dy=1/Jmax;
X=nlm:
1;
Y=nlm:
1;
NN=zeros(Imax+1,Jmax+1);
n1=0;
forj=3:
Jmax/2
fori=2:
j-1
n1=n1+1;
NN(i,j)=n1;%X=i列,Y=j行处节点
X(n1)=(i-1)*dx;
Y(n1)=-1+(j-1)*dy;
end
end
k=Jmax/2+1;
forj=Jmax/2+1:
Jmax-1%三角形区域上下两部分节点坐标分别求
k=k-1;
fori=2:
k
n1=n1+1;
NN(i,j)=n1;
X(n1)=(i-1)*dx;
Y(n1)=1+(j-Jmax-1)*dy;
end
end
na=n1;%不含边界节点数
forj=Jmax+1:
-1:
Jmax/2+1%降序
n1=n1+1;
NN(1,j)=n1;
X(n1)=0;
Y(n1)=1+(j-Jmax-1)*dy;
end
forj=Jmax/2:
-1:
1
n1=n1+1;
NN(1,j)=n1;
X(n1)=0;
Y(n1)=-1+(j-1)*dy;
end%
fori=2:
Imax+1
n1=n1+1;
NN(i,i)=n1;
X(n1)=(i-1)*dx;
Y(n1)=-1+(i-1)*dy;
end
K=0;
fori=Imax:
-1:
2
K=K+2;
n1=n1+1;
NN(i,i+K)=n1;
X(n1)=(i-1)*dx;
Y(n1)=1+(i+K-Jmax-1)*dy;
end
%以上四个循环为对边界节点进行编号
ndm=n1;
NE=zeros(3,2*ndm);n1=0;
forj=3:
Jmax/2
fori=2:
j-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
k=Jmax/2+1;
forj=Jmax/2+1:
Jmax-1
k=k-1;
fori=2:
k
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%内部节点对应左上角正方形的两个三角形单元,上左,左上
fori=2:
Imax
n1=n1+1;
NE(1,n1)=NN(i,i);
NE(2,n1)=NN(i-1,i);
NE(3,n1)=NN(i-1,i-1);
n1=n1+1;
NE(1,n1)=NN(i,i);
NE(2,n1)=NN(i-1,i+1);
NE(3,n1)=NN(i-1,i);
n1=n1+1;
NE(1,n1)=NN(i,i);
NE(2,n1)=NN(i,i+1);
NE(3,n1)=NN(i-1,i+1);
end%下斜边界节点要对应三个单元,上左,左上,左下
n1=n1+1;
NE(1,n1)=NN(Imax+1,Imax+1);
NE(2,n1)=NN(Imax,Imax+1);
NE(3,n1)=NN(Imax,Imax);%右顶点的左下
n1=n1+1;
NE(1,n1)=NN(Imax+1,Imax+1);
NE(2,n1)=NN(Imax,Imax+2);
NE(3,n1)=NN(Imax,Imax+1);%右顶点的左上
K=0;
fori=Imax:
-1:
2
K=K+2;
n1=n1+1;
NE(1,n1)=NN(i,i+K);
NE(2,n1)=NN(i-1,i+K+1);
NE(3,n1)=NN(i-1,i+K);
end%右上边界的左上
nel=n1;%此时n1值为总的单元个数
六、输出文本文件
*********有限元法求解三角形区域上Possion方程的结果**********
节点编号坐标分量x坐标分量yu(x,y)的值
10.01333-0.986670.00000
20.01333-0.980000.00001
30.02667-0.980000.00001
40.01333-0.973330.00002
50.02667-0.973330.00003
60.04000-0.973330.00002
70.01333-0.966670.00002
80.02667-0.966670.00004
90.04000-0.966670.00005
100.05333-0.966670.00003
110.01333-0.960000.00003
120.02667-0.960000.00006
130.04000-0.960000.00007
140.05333-0.960000.00007
150.06667-0.960000.00005
160.01333-0.953330.00004
170.02667-0.953330.00008
180.04000-0.953330.00010
190.05333-0.953330.00011
200.06667-0.953330.00010
210.08000-0.953330.00007
220.01333-0.946670.00005
230.02667-0.946670.00010
240.04000-0.946670.00014
250.05333-0.946670.00016
260.06667-0.946670.00016
270.08000-0.946670.00013
280.09333-0.946670.00008
290.01333-0.940000.00006
300.02667-0.940000.00012
310.04000-0.940000.00017
320.05333-0.940000.00020
330.06667-0.940000.00022
340.08000-0.940000.00021
350.09333-0.940000.00017
360.10667-0.940000.00010
370.01333-0.933330.00007
380.02667-0.933330.00015
390.04000-0.933330.00021
400.05333-0.933330.00025
410.06667-0.933330.00028
420.08000-0.933330.00028
430.09333-0.933330.00026
440.10667-0.933330.00021
450.12000-0.933330.00012
460.01333-0.926670.00009
470.02667-0.926670.00017
480.04000-0.926670.00024
490.05333-0.926670.00030
500.06667-0.926670.00034
510.08000-0.926670.00036
520.09333-0.926670.00036
530.10667-0.926670.00032
540.12000-0.926670.00025
550.13333-0.926670.00015
560.01333-0.920000.00010
570.02667-0.920000.00020
580.04000-0.920000.00028
590.05333-0.920000.00036
600.06667-0.920000.00041
610.08000-0.920000.00045
620.09333-0.920000.00046
630.10667-0.920000.00044
640.12000-0.920000.00038
650.13333-0.920000.00030
660.14667-0.920000.00017
670.01333-0.913330.00011
680.02667-0.913330.00022
690.04000-0.913330.00032
700.05333-0.913330.00041
710.06667-0.913330.00048
720.08000-0.913330.00053
730.09333-0.913330.00056
740.10667-0.913330.00056
750.12000-0.913330.00052
760.13333-0.913330.00045
770.14667-0.913330.00034
780.16000-0.913330.00019
790.01333-0.906670.00013
800.02667-0.906670.00025
810.04000-0.906670.00037
820.05333-0.906670.00047
830.06667-0.906670.00055
840.08000-0.906670.00062
850.09333-0.906670.00066
860.10667-0.906670.00068
870.12000-0.906670.00066
880.13333-0.906670.00061
890.14667-0.906670.00052
900.16000-0.906670.00039
910.17333-0.906670.00022
920.01333-0.900000.00014
930.02667-0.900000.00028
940.04000-0.900000.00041
950.05333-0.900000.00053
960.06667-0.900000.00063
970.08000-0.900000.00071
980.09333-0.900000.00077
990.10667-0.900000.00080
1000.12000-0.900000.00080
1010.13333-0.900000.00077
1020.14667-0.900000.00070
1030.16000-0.900000.00059
1040.17333-0.900000.00044
1050.18667-0.900000.00024
1060.01333-0.893330.00016
1070.02667-0.893330.00031
1080.04000-0.893330.00045
1090.05333-0.893330.00059
1100.06667-0.893330.00071
1110.08000-0.893330.00080
1120.09333-0.893330.00088
1130.10667-0.893330.00093
1140.12000-0.893330.00095
1150.13333-0.893330.00093
1160.14667-0.893330.00089
1170.16000-0.893330.00080
1180.17333-0.893330.00067
1190.18667-0.893330.00049
1200.20000-0.893330.00027
1210.01333-0.886670.00017
1220.02667-0.886670.00034
1230.04000-0.886670.00050
1240.05333-0.886670.00065
1250.06667-0.886670.00078
1260.08000-0.886670.00090
1270.09333-0.886670.00099
1280.10667-0.886670.00106
1290.12000-0.886670.00110
1300.13333-0.886670.00110
1310.14667-0.886670.00107
1320.16000-0.886670.00101
1330.17333-0.886670.00090
1340.18667-0.886670.00074
1350.20000-0.886670.00055
1360.21333-0.886670.00030
1370.01333-0.880000.00019
1380.02667-0.880000.00037
1390.04000-0.880000.00055
1400.05333-0.880000.00071
1410.06667-0.880000.00086
1420.08000-0.880000.00100
1430.09333-0.880000.00111
1440.10667-0.880000.00119
1450.12000-0.880000.00125
1460.13333-0.880000.00127
1470.14667-0.880000.00126
1480.16000-0.880000.00122
1490.17333-0.880000.00113
1500.18667-0.880000.001
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlaB 程序 有限元 法解泊松 方程
![提示](https://static.bingdoc.com/images/bang_tan.gif)