基于Matlab语言按平面三角形单元划分结构有限元程序设计Word格式.docx
《基于Matlab语言按平面三角形单元划分结构有限元程序设计Word格式.docx》由会员分享,可在线阅读,更多相关《基于Matlab语言按平面三角形单元划分结构有限元程序设计Word格式.docx(21页珍藏版)》请在冰豆网上搜索。
求解式①所示的有限元线性方程,得到节点的位移。
在该步骤中,若有限元模型的节点越多,则线性方程的数量就越多,随之有限元分析的计算量也将越大。
(7)求解单元应力及应变根据求出的节点位移求解单元的应力和应变。
(8)结果处理与显示进入有限元分析的后处理部分,对计算出来的结果进行加工处理,并以各种形式将计算结果显示出。
2.Matlab简介
在用有限元法进行结构分析时,将会遇到大量的数值计算,因而在实用上是一定要借助于计算机和有限元程序,才能完成这些复杂而繁重的数值计算工作。
而Matlab是当今国际科学界最具影响力和活力的软件。
它起源于矩阵运算,并已经发展成一种高度集成的计算机语言。
它提供了强大的科学计算,灵活的程序设计流程,高质量的图形可视化与界面设计,便捷的与其他程序和语言接口的功能。
Matlab在各国高校与研究单位起着重大的作用。
“工欲善其事,必先利其器”。
如果有一种十分有效的工具能解决在教学与研究中遇到的问题,那么Matlab语言正是这样的一种工具。
它可以将使用者从繁琐、无谓的底层编程中解放出来,把有限的宝贵时间更多地花在解决问题中,这样无疑会提高工作效率。
目前,Matlab已经成为国际上最流行的科学与工程计算的软件工具,现在的Matlab已经不仅仅是一个“矩阵实验室”了,它已经成为了一种具有广泛应用前景的全新的计算机高级编程语言了,有人称它为“第四代”计算机语言,它在国内外高校和研究部门正扮演着重要的角色。
Matlab语言的功能也越来越强大,不断适应新的要求提出新的解决方法。
可以预见,在科学运算、自动控制与科学绘图领域Matlab语言将长期保持其独一无二的地位。
为此,本例采用Matlab语言编程,以利用其快捷强大的矩阵数值计算功能。
二、问题描述
一矩形薄板,一边固定,承受150kN集中荷载,结构简图及按平面三角形单元划分的有限元模型图如下所示。
材料参数:
弹性模量
;
泊松比:
薄板厚度
。
在本例中,所取结构模型及数据主要用于程序设计理论分析,与工程实际无关。
三、参数输入:
单元个数NELEM=4
节点个数NNODE=6受约束边界点数NVFIX=2
节点荷载个数NFORCE=1
弹性模量YOUNG=2e8
泊松比POISS=0.2
厚度THICK=0.002
单元节点编码数组LNODS=
节点坐标数组COORD=
节点力数组FORCE=[40-150]
约束信息数组FIXED=
以上数值数据为程序运行前输入的初始数据,存为“471220580.txt”文本格式,同时必须放在Matlab工作目录下,路径不对程序不能自动读取指定初始文件,运行出错。
初始数据文本文件输入格式如下图:
四、Matlab语言程序源代码:
1.程序中变量说明
NNODE单元节点数
NPION总结点数
NELEM单元数
NVFIX受约束边界点数
FIXED约束信息数组
NFORCE节点力数
FORCE节点力数组
COORD结构节点坐标数组
LNODS单元定义数组
YOUNG弹性模量
POISS泊松比
THICK厚度
B单元应变矩阵(3*6)
D单元弹性矩阵(3*3)
S单元应力矩阵(3*6)
A单元面积
ESTIF单元刚度矩阵
ASTIF总体刚度矩阵
ASLOD总体荷载向量
ASDISP节点位移向量
ELEDISP单元节点位移向量
STRESS单元应力
FP1数据文件指针
2.Matlab语言程序代码
%******************************************************************************
%初始化及数据调用
formatshorte%设定输出类型
clear%清除内存变量
FP1=fopen('
471220580.txt'
'
rt'
);
%打开输入数据文件,读入控制数据
NELEM=fscanf(FP1,'
%d'
1),%单元个数(单元编码总数)
NPION=fscanf(FP1,'
1),%结点个数(结点编码总数)
NVFIX=fscanf(FP1,'
1)%受约束边界点数
NFORCE=fscanf(FP1,'
1),%结点荷载个数
YOUNG=fscanf(FP1,'
%e'
1),%弹性模量
POISS=fscanf(FP1,'
%f'
1),%泊松比
THICK=fscanf(FP1,'
1)%厚度
LNODS=fscanf(FP1,'
[3,NELEM])'
%单元定义数组(单元结点号)
%相应为单元结点号(编码)、按逆时针顺序输入
COORD=fscanf(FP1,'
[2,NPION])'
%结点坐标数组
%坐标:
x,y坐标(共NPOIN组)
FORCE=fscanf(FP1,'
[3,NFORCE])'
%结点力数组(受力结点编号,x方向,y方向)
FIXED=fscanf(FP1,'
[3,NVFIX])'
%约束信息(约束点,x约束,y约束)
%有约束为1,无约束为0
%*****************************************************************************
%生成单元刚度矩阵并组成总体刚度矩阵
ASTIF=zeros(2*NPION,2*NPION);
%生成特定大小总体刚度矩阵并置0
fori=1:
NELEM
%生成弹性矩阵D
D=[1POISS0;
POISS10;
00(1-POISS)/2]*YOUNG/(1-POISS^2)
%计算当前单元的面积A
A=det([1COORD(LNODS(i,1),1)COORD(LNODS(i,1),2);
1COORD(LNODS(i,2),1)COORD(LNODS(i,2),2);
1COORD(LNODS(i,3),1)COORD(LNODS(i,3),2)])/2
%生成应变矩阵B
forj=0:
2
b(j+1)=COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);
c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1);
end
B=[b
(1)0b
(2)0b(3)0;
...
0c
(1)0c
(2)0c(3);
c
(1)b
(1)c
(2)b
(2)c(3)b(3)]/(2*A);
B1(:
:
i)=B;
%求应力矩阵S=D*B
S=D*B;
ESTIF=B'
*S*THICK*A;
%求解单元刚度矩阵
a=LNODS(i,:
%临时向量,用来记录当前单元的节点编号
forj=1:
3
fork=1:
ASTIF((a(j)*2-1):
a(j)*2,(a(k)*2-1):
a(k)*2)…
=ASTIF((a(j)*2-1):
a(k)*2)+ESTIF(j*2-1:
j*2,k*2-1:
k*2);
%根据节点编号对应关系将单元刚度分块叠加到总刚
%度矩阵中
%将约束信息加入总体刚度矩阵(对角元素改一法)
fori=1:
NVFIX
ifFIXED(i,2)==1
ASTIF(:
(FIXED(i,1)*2-1))=0;
%一列为零
ASTIF((FIXED(i,1)*2-1),:
)=0;
%一行为零
ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1;
%对角元素为1
ifFIXED(i,3)==1
ASTIF(:
FIXED(i,1)*2)=0;
%一列为零
ASTIF(FIXED(i,1)*2,:
%一行为零
ASTIF(FIXED(i,1)*2,FIXED(i,1)*2)=1;
%生成荷载向量
ASLOD(1:
2*NPION)=0;
%总体荷载向量置零
NFORCE
ASLOD((FORCE(i,1)*2-1):
FORCE(i,1)*2)=FORCE(i,2:
3);
%求解内力
ASDISP=ASTIF\ASLOD'
%计算节点位移向量
ELEDISP(1:
6)=0;
%当前单元节点位移向量
ELEDISP(j*2-1:
j*2)=ASDISP(LNODS(i,j)*2-1:
LNODS(i,j)*2);
%取出当前单元的节点位移向量
i
STRESS=D*B1(:
:
i)*ELEDISP'
%求内力
fclose(FP1);
%关闭数据文件
五、程序运行结果
NELEM=
4
NPION=
6
NVFIX=
2
NFORCE=
1
YOUNG=
200000000
POISS=
2.0000e-001
THICK=
2.0000e-003
LNODS=
126
234
245
256
COORD=
00
10
20
21
11
01
FORCE=
40-150
FIXED=
111
611
D=
2.0833e+0084.1667e+0070
4.1667e+0072.0833e+0080
008.3333e+007
A=
5.0000e-001
ASDISP=
0
-8.0607e-004
-1.5848e-003
-1.0281e-003
-4.4727e-003
1.1937e-003
-4.6947e-003
8.6670e-004
-1.8880e-003
(说明:
以上12个ASDISP输出结果数据表示节点位移向量,即各节点分别在X,Y方向上产生的位移。
)
i=
STRESS=
-1.6793e+005
-3.3586e+004
-1.3207e+005
-5.5503e+004
3
5.5503e+004
-4.9526e+004
-9.4497e+004
1.6793e+005
-2.7040e+004
-1.7932e+004
以上四组STRESS输出结果数据表示单元应力SX,SY,SXY,i为单元号。
六、ANSYS建模比较分析
利用ANSYS有限元分析软件,完全按照前面Matlab程序设计的各项条件参数以及单元划分方式建立ANSYS模型,其求解结果与以上程序计算结果比较,验证程序是否正确。
ANSYS模型节点单元如下图所示:
ANSYS模型求解变形图如下所示:
ANSYS求解节点位移结果列表显示如下:
(单位:
m)
PRINTUNODALSOLUTIONPERNODE
*****POST1NODALDEGREEOFFREEDOMLISTING*****
THEFOLLOWINGDEGREEOFFREEDOMRESULTSAREINTHEGLOBALCOORDINATESYSTEM
NODEUXUYUZUSUM
10.00000.00000.00000.0000
2-0.80607E-03-0.15848E-020.00000.17780E-02
3-0.10281E-02-0.44727E-020.00000.45893E-02
40.11937E-02-0.46947E-020.00000.48441E-02
50.86670E-03-0.18880E-020.00000.20774E-02
60.00000.00000.00000.0000
MAXIMUMABSOLUTEVALUES
NODE4404
VALUE0.11937E-02-0.46947E-020.00000.48441E-02
ANSYS求解单元应力结果列表显示如下:
KN/m2)
PRINTSELEMENTSOLUTIONPERELEMENT
*****POST1ELEMENTNODALSTRESSLISTING*****
THEFOLLOWINGX,Y,ZVALUESAREINGLOBALCOORDINATES
ELEMENT=1PLANE182
NODESXSYSZSXYSYZ
1-0.16793E+06-33586.0.0000-0.13207E+060.0000
2-0.16793E+06-33586.0.0000-0.13207E+060.0000
6-0.16793E+06-33586.0.0000-0.13207E+060.0000
ELEMENT=2PLANE182
2-55503.-55503.0.0000-55503.0.0000
3-55503.-55503.0.0000-55503.0.0000
4-55503.-55503.0.0000-55503.0.0000
ELEMENT=3PLANE182
255503.-49526.0.0000-94497.0.0000
455503.-49526.0.0000-94497.0.0000
555503.-49526.0.0000-94497.0.0000
ELEMENT=4PLANE182
20.16793E+06-27040.0.0000-17932.0.0000
50.16793E+06-27040.0.0000-17932.0.0000
60.16793E+06-27040.0.0000-17932.0.0000
结论
通过比较Matlab语言设计程序运行结果和ANSYS建模分析结果可知,两种方式算出的结果完全一致,说程序设计正确。
所以本程序适用于按三角形单元划分的平面结构有限元分析。
formatshorte
clear
LinearTriangleElementofThinplate.txt'
1)
[3,NELEM])
[2,NPION])
[3,NFORCE])
[3,NVFIX])
%生成特定大小总体刚度矩阵并置0
D=YOUNG/(1-POISS^2)*[1POISS0;
POISS10;
00(1-POISS)/2]
A=det([1COORD(LNODS(i,1),1)COORD(LNODS(i,1),2);
1COORD(LNODS(i,2),1)COORD(LNODS(i,2),2);
1COORD(LNODS(i,3),1)COORD(LNODS(i,3),2)])/2
forj=0:
b(j+1)=COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);
c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1);
end
B=[b
(1)0b
(2)0b(3)0;
0c
(1)0c
(2)0c(3);
c
(1)b
(1)c
(2)b
(2)c(3)b(3)]/(2*A);
a=LNODS(i,:
forj=1:
fork=1:
ASTIF((a(j)*2-1):
=ASTIF((a(j)*2-1):
%将约束信息加入总体刚度矩阵
ifFIXED(i,2)==1
ASTIF(:
ASTIF((FIXED(i,1)*2-1),:
ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1;
ifFIXED(i,3)==1
ASTIF(:
ASTIF(FIXED(i,1)*2,:
ASTIF(FIXED(i,1)*2,FIXED(i,1)*2)=1;
ASLOD(1:
ASLOD((FORCE(i,1)*2-1):
3)
%求解内力
ASDISP=ASTIF\ASLOD'
ELEDISP(1:
ELEDISP(j*2-1:
i
STRESS=D*B1(:
i)*ELE