基于matlab的空间三维桁架结构受力分析通用程序设计Word文档格式.docx
《基于matlab的空间三维桁架结构受力分析通用程序设计Word文档格式.docx》由会员分享，可在线阅读，更多相关《基于matlab的空间三维桁架结构受力分析通用程序设计Word文档格式.docx（11页珍藏版）》请在冰豆网上搜索。
1空间三维桁架结构受力分析的有限元法1.1空间三维桁架结构的有限元计算格式空间桁架元是既有局部坐标,又有总体坐标的三维有限元,单元类型用线性函数描述.空间桁架元的系数有弹性模量E,横截面积A和长度L.每个空间桁架单元有2个结点,并且在从总体坐标系x,y,Z轴到局部坐标系z轴下的倾斜角分别为,和,如图1所示.假设=cos,=
cos0.,则单元刚度矩阵给定如式
（1）所示[1--2]:
一一一
ll一一—
一如一一:
一{l蕾{l,{l
:
扎
（1）
图1空间桁架元
空间桁架元的每个结点有三个自由度（,收稿日期:
2oo5.O3.O2
作者简介:
李罡（1978.）,男,空军第一航空学院一系发动机教研室助教,现为北京航
空航天大学能源与动力工程学
院在读硕士研究生,主要研究方向为结构振动与应力.一
一,
?
36?
喀什师范学院第26卷
和）,因此,其每个单元有6个自由度.则对于一个有;
r/个结点的桁架结构而言,其整体刚度矩阵K是3×
3;
r/的.在得到整体刚度矩阵K后,就可以得到如下的方程组:
[K]{【,}={F},
（2）
式中,【,是结构结点位移矢量,F是结构结点力矢量.对方程应用高斯消去法就可以得出未知的位移和支反力,就可以用下式求出每个单元（杆件）的结点力:
=[]
（3）
（3）式中,,是单元结点力（标量）,lI是6×
1的单元结点位移矢量.将单元结点力除以横截面积A就可以得到单元应力.
2.2刚度矩阵组集方法
对每个单元分析得到单元刚度矩阵后,由于其是6×
6阶的矩阵,要组集成3n×
37"
/的总刚矩阵,必须首先对其进行拓展,按单元刚度矩阵的结点编号叠加,放入总刚矩阵相应的位置.单元刚度矩阵的这个变换可以起到两个作用【]:
（1）将单元刚度矩阵K扩大到与结构（总）刚度矩阵同阶,以便进行矩阵相加;
（2）将单元刚度矩阵中的各个子块按照单元结点的实际编码安放在扩大的矩阵中,可以反映出该单元对结构（总）刚度矩阵K的"
贡献"
.2.3约束条件的处理方法
在依据单元刚度矩阵组集得到的总刚度矩阵
K具有奇异性,即任意给定结构的结点位移所得到的结构结点力总体上是满足力和力矩的平衡的,因此,不能直接对式
（2）求解得出结构的位移【,来,这反映了结构可能发生任意的刚体位移.为了消除K的奇异性,则需要对结构施加约束条件,限定结构的刚体位移,才能解出结点位移.这样做的另一个好处是还可以使总刚矩阵的阶次降低,便于程序求解.
为保证结构的结点位移的顺序性不被破坏,这里采用置1赋0法处理约束条件[.即在刚度矩阵K中,与零结点位移相对应的行列中,将主对角元素改为1,其他元素改为0;
在载荷列阵中,将与零结点位移相对应的元素改为0.
2Matlab通用程序设计
2.1通用程序组成及设计思想
该程序是基于上述对空间三维桁架结构受力分析的有限元计算格式,将各个需要完成的模块编制成相应的子程序,再对整体求解过程编制主程序实现通用求解功能.各子程序分别完成数据输入,求解单元刚度矩阵,组集总刚度矩阵,约束条件处理和位移,应力求解等任务.其程序流程如图2所示.
图2程序流程图
2.2各个子程序说明
2.2.1[M,N]=SpaeeTrussdata（filename）
该函数用于获得（从键盘输入）所有关于求解问题的相关数据,并将这些数据按要求分别各全局变量,以便在子程序和主程序中调用.源程序如下:
funcfion
[ndema,nnodes]=SpaeeTrussdata（filename）
globalTITLEN0DEFIXXGYGZGELEMDEFFORCEAEFKDISPMNn
nderm=input（'
PleaseentertheNumberofElements'
）;
nnodes=input（'
PleaseentertheNumberofNodes'
fori=1:
nnodes
disp（'
PleaseenterthefollowingforNode'
disp（i）;
NODEFIX（i,1）=input（'
FixedinXdirection?
（0=free,1=fixed）'
N0DEFIX（i,2）=input（'
FixedinYdirection?
（0=flee,1=fixed）'
N0DEFIX（i,3）=input（'
FixedinZdirection?
XG（i）=input（'
GlobalXnodalcoordinate'
YG（i）=input（'
GlobalYnodalcoordinate'
ZG（i）=input（'
GlobalZnodalcoordinate'
FoRCE（i,1）=input（'
ForceinX—direction'
FoRCE（i,2）=input（'
ForceinY—direction'
FoRCE（i,3）=input（'
ForceinZ—direction'
end
fori=1:
ndems;
第3期李罡:
基于matlab的空间三维桁架结构受力分析通用程序设计?
37?
ELEMDEF（i,1）=input（'
LeftNode'
ELEMDEF（i,2）=input（'
RightNode'
E（i）=input（'
YoungsmoduluS'
A（i）=input（'
CrossSectionArea'
2.2.2[ke]=SpaceTrussElementStiffness（jde.
ment）
该函数首先根据各个桁架元的杨式模量E,横截面积A和单元左右结点的三个方向的坐标,用前面给出的有限元计算公式求得单元刚度矩阵（6×
6）中的四个子块中的一个子块（3×
3）,再利用单元刚度矩阵对称性的特点,将四个子块按对应的节点编号放人已经定义过的单元刚度矩阵中.源程序如下:
function
[ke]=SpaceTrussElementStiffness（jdement）
globalTITLENODEFIXXGYGZGELEMDEFFORCEA
EFKDISPMNSr】SS
ske=zeros（3.3）;
ke=zeros（3*N,3*N）;
l:
ELEMDEF（jdement,1）;
r:
ELEMDEF（jdement,2）;
xl=XG（I）;
yl=YG（I）;
zl=ZG（I）;
xr=XG（r）:
yr=YG（r）;
zl-=ZG（r）;
Ix=xr—xl:
ly=yr——yl;
lz=zr——zl;
le=sqrt（Ix-2+ly^2+）;
nxx=lxAe）
nxy=ly/le;
I1）
【z=lzAe;
ske=A（jdement）*E（idement）*[nxx'
2nxx*nxynxx*nxz;
nxx*nxynxy'
2nxynxz;
nxx*nxznxynxzI吁2]/
le;
ke（[（3*l一2）:
（3*1）],[（3*l一2）:
（3*1）]）=ske;
ke（[（3*r一2）:
（3*r）],[（3*r一2）:
（3*r）]）=ske;
（3*1）],[（3*r一2）:
（3*r）]）=一ske;
（3*r）],[（3*l一2）:
（3*1）]）=一ske;
2.2.3[K]=SpaceTrussAssemble（jdement,ke,k）此函数的功能是将扩展后的单元刚度矩阵依
据其节点编号顺序放人总刚度矩阵中,叠加生成总
刚度矩阵.源程序如下:
[K]=SpaceTrussAssemble（jdement,ke,K）globalTITLENODEFIXXGYGZGELEMDEFFORCEAEFKDISPMNI正:
K=K+ke:
2.2.4[K,F]=SpaceTrussconstrain（jdof,K,F）该函数用置1赋0法对总刚度矩阵及力向量
进行处理,使总刚度矩阵变为非奇异,为求解位移
向量作准备.源程序如下:
[K,F]=SpaceTrussconstrain（jdof,K,F）globalTITLENODEFIXXGYGZGELEMDEFFORCEAEFKDISPMNRESS
K（idof,:
）=0;
K（:
jdof）=0;
K（jdof,ido1）=1;
F（jdof）:
0;
2.2.5[DISP]=SpaceTrusssolve（K,F）
该函数用于求解单元结点位移向量,并将该向量转化为N×
3维的矩阵.源程序如下:
[DISP]=SpaceTrusssolve（K,F）globalTITLEN0DEFIXXGYGZGELEMDEFF0RCEA
EFKDISPMNRESS
di=my（K）*F:
fori=l:
N
DISP（i,1）=di（3*i一2）;
DISP（i,2）=di（3*i一1）;
DISP（i,3）=di（3*i）;
2.2.6[STRESS]=SpaceTrussElementStress
（jelement）
该函数根据求出的位移向量最后解出单元的应力.源程序如下:
[STRF~]=SpaceT~lementStress（jdement）
l=ELEMDEF（jdement,1）;
r=ELEMDEF（jdement,2）;
xl=XG
（1）;
zl=ZG
（1）;
】【r=XG（r）;
zr=ZG（r）;
Ix=xr—xl;
ly=yr—yl;
lz=一zl:
nxx--~lxAe;
Iu【y=lyAe;
38?
c=[（DISP（r,1）一DISP（I,1））;
（DISP（r,2）一DISP（1,2））;
（DISP（r,3）一DISP（I,3））];
[STRESS（t,ielement）]=E（ielement）nxxnxym【z]*c/le;
2.3主程序编制
根据上述各个子程序,利用空间三维桁架结构
的有限元计算格式,编制出主程序,完成求解任务.源程序如下:
dear
dc
m眦e=input（'
输人文件名\n'
'
8'
）
[M,N]=SpaceTmssc~ta（filename）K=~os（3*N,3*N）;
forjdement=1:
M
jelement
[ke]=SpaceTrussElementSfiffness（jdement）
[K]=SpaceT~ble（jdement,ke,K）end
F=z翻（3*N,1）;
F（3*i一2）=FORCE（i,1）;
F（3*i一1）:
FoRCE（i,2）;
F（3*i）=FORCE（i,3）;
end
forj=1:
3
ifNODEFIX（i,j）:
=1jdof=3*（i一1）+j;
[K,F]=SpaceTnmseonslwain（jdd,K,F）
DISP=~os（N,3）;
DISP=SpaceTrusssolve（K,F）
RESS=踟（1,M）
forje[er~ent=1:
M[STRESS]=SpaeeTrtmsElementStress（jelement）
TITLE'
TITLE
XG'
XG
YG'
YG
ZG'
ZG
FORCE'
FORcE
NODEFIX'
NODEFIX
ELEMDEF'
ELEN伍）EF
d/sp（'
A'
A
（'
E'
E
赋1置0后K='
K
DISP'
DISP
STRESS'
RESS
3结论
本程序经多个算例验证,都得到正确结果.应
用此程序,即可对任意条件的空间三维桁架结构进
行快速的力学分析,判断结构的受力情况,为进一
步的工程设计奠定基础.因此,具有较强的工程实
用价值.
参考文献:
[1]LoganD.AFirstCourseinFiniteElementMethod[M].
SecondEdition.PWSPublishingCompany,1993.
[2]SennettR.MatrixAnalysisofStructure[M].Prentice
Hau.1994.
[3]王瑁成.有限单元法[M].北京:
清华大学出版社,
2003.7.
[4]张志涌,徐彦琴.等,MATLAB教程——基于6.x版本
[M].北京:
北京航空航天大学出版社,2001.4.
Designof3DSpaceTrussSolverBasedonMATLAB
LIGang
（Aeronautics&
amp;
A..tronauties,BdjingUniversity,Beijing100083,China）
Abstract:
Thefiniteelementmethodof3Dspacetrussforceanalysisisintroduced.Theprogramsettlesthe
problemshoWtoassemblethestiffnessmatrixandhowtoaispo~theconstrainbasedonmatlablanguage.It
isall—
purposeprogram.anditcancalculateallkindofforceanalysisproblemsabout3Dspacetruss.Thepro—
gramCAInbeveryusefulonprojects.
Keywords:
3Dtruss;
Finiteelementmethod;
Matlab;
An—purposeprogramdesign
升级会员 