基于matlab的空间三维桁架结构受力分析通用程序设计Word文档格式.docx

1、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.）,男,空军第一航空学院一系发动机教

2、研室助教,现为北京航空航天大学能源与动力工程学 院在读硕士研究生,主要研究方向为结构振动与应力. 一 一, ?36?喀什师范学院第26卷 和）,因此,其每个单元有6个自由度.则对于一 个有;r/个结点的桁架结构而言,其整体刚度矩阵 K是33;r/的.在得到整体刚度矩阵K后,就可 以得到如下的方程组:K【,=F,（2） 式中,【,是结构结点位移矢量,F是结构结点力矢 量.对方程应用高斯消去法就可以得出未知的位移 和支反力,就可以用下式求出每个单元（杆件）的结 点力:,= （3） （3）式中,是单元结点力（标量）,lI是61的单 元结点位移矢量.将单元结点力除以横截面积A 就可以得到单元应力. 2

3、.2刚度矩阵组集方法 对每个单元分析得到单元刚度矩阵后,由于其 是66阶的矩阵,要组集成3n37/的总刚矩阵, 必须首先对其进行拓展,按单元刚度矩阵的结点编 号叠加,放入总刚矩阵相应的位置.单元刚度矩阵 的这个变换可以起到两个作用【:（1）将单元刚度矩阵K扩大到与结构（总）刚 度矩阵同阶,以便进行矩阵相加;（2）将单元刚度矩阵中的各个子块按照单元 结点的实际编码安放在扩大的矩阵中,可以反映出 该单元对结构（总）刚度矩阵K的贡献. 2.3约束条件的处理方法 在依据单元刚度矩阵组集得到的总刚度矩阵 K具有奇异性,即任意给定结构的结点位移所得 到的结构结点力总体上是满足力和力矩的平衡的, 因此,不能

4、直接对式（2）求解得出结构的位移【, 来,这反映了结构可能发生任意的刚体位移.为了 消除K的奇异性,则需要对结构施加约束条件,限 定结构的刚体位移,才能解出结点位移.这样做的 另一个好处是还可以使总刚矩阵的阶次降低,便于 程序求解. 为保证结构的结点位移的顺序性不被破坏,这 里采用置1赋0法处理约束条件.即在刚度矩阵 K中,与零结点位移相对应的行列中,将主对角元 素改为1,其他元素改为0;在载荷列阵中,将与零 结点位移相对应的元素改为0. 2Matlab通用程序设计 2.1通用程序组成及设计思想 该程序是基于上述对空间三维桁架结构受力 分析的有限元计算格式,将各个需要完成的模块编 制成相应的子

5、程序,再对整体求解过程编制主程序 实现通用求解功能.各子程序分别完成数据输入, 求解单元刚度矩阵,组集总刚度矩阵,约束条件处 理和位移,应力求解等任务.其程序流程如图2所 示. 图2程序流程图 2.2各个子程序说明 2.2.1M,N=SpaeeTrussdata（filename） 该函数用于获得（从键盘输入）所有关于求解 问题的相关数据,并将这些数据按要求分别各全局 变量,以便在子程序和主程序中调用.源程序如下: funcfion ndema,nnodes=SpaeeTrussdata（filename） globalTITLEN0DEFIXXGYGZGELEMDEFFORCEA EFKDI

6、SPMNn 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）

7、=input（GlobalXnodalcoordinate YG（i）=input（GlobalYnodalcoordinate ZG（i）=input（GlobalZnodalcoordinate FoRCE（i,1）=input（ForceinXdirectionFoRCE（i,2）=input（ForceinYdirectionFoRCE（i,3）=input（ForceinZdirectionend fori=1:ndems;第3期李罡:基于matlab的空间三维桁架结构受力分析通用程序设计?37?ELEMDEF（i,1）=input（LeftNodeELEMDEF（i,2）=inpu

8、t（RightNode E（i）=input（YoungsmoduluSA（i）=input（CrossSectionArea2.2.2ke=SpaceTrussElementStiffness（jde. ment） 该函数首先根据各个桁架元的杨式模量E, 横截面积A和单元左右结点的三个方向的坐标, 用前面给出的有限元计算公式求得单元刚度矩阵 （66）中的四个子块中的一个子块（33）,再利用 单元刚度矩阵对称性的特点,将四个子块按对应的 节点编号放人已经定义过的单元刚度矩阵中.源程 序如下:function ke=SpaceTrussElementStiffness（jdement） glob

9、alTITLENODEFIXXGYGZGELEMDEFFORCEA 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=xrxl:ly=yryl;lz=zrzl;le=sqrt（Ix-2+ly2+）; nxx=lxAe） nxy=ly/le;I1）【z=lzAe;ske=A（jdement）*E（idement）*nxx2nxx*nxynxx* n

10、xz;nxx*nxynxy2nxynxz;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.3K=SpaceTrussAssemble（jdement,ke,k） 此函数的功能是将扩展后的单元刚度矩阵依 据其节点编号顺序放人总刚度矩阵中,叠加生成总 刚度矩阵.源程序如下:K=SpaceTrussAssemble（jdement,ke,K） globalTITLE

11、NODEFIXXGYGZGELEMDEFFORCEA EFKDISPMNI正:K=K+ke:2.2.4K,F=SpaceTrussconstrain（jdof,K,F） 该函数用置1赋0法对总刚度矩阵及力向量 进行处理,使总刚度矩阵变为非奇异,为求解位移 向量作准备.源程序如下:K,F=SpaceTrussconstrain（jdof,K,F） globalTITLENODEFIXXGYGZGELEMDEFFORCEA EFKDISPMNRESS K（idof,:）=0;K（:,jdof）=0;K（jdof,ido1）=1;F（jdof）:0;2.2.5DISP=SpaceTrusssolve

12、（K,F） 该函数用于求解单元结点位移向量,并将该向 量转化为N3维的矩阵.源程序如下: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.6STRESS=SpaceTrussElementStress （jelement） 该函数根据求出的位移向量最后解出单元的 应力.源程序如下:STRF=SpaceTlementStress（

13、jdement） l=ELEMDEF（jdement,1）; r=ELEMDEF（jdement,2）; xl=XG（1）;zl=ZG（1）;】【r=XG（r）;zr=ZG（r）;Ix=xrxl;ly=yryl;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主程序编制 根据上述各个子程序,利用空间三维桁架结构 的有限元计算格式,编制出主程序,完成求解任务.

14、 源程序如下:dear dc m眦e=input（输人文件名n,8） M,N=SpaceTmsscta（filename） K=os（3*N,3*N）;forjdement=1:M jelement ke=SpaceTrussElementSfiffness（jdement） K=SpaceTble（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）:=1 jdof=3*（i一1）+j;K,F=SpaceTnm

15、seonslwain（jdd,K,F） DISP=os（N,3）; DISP=SpaceTrusssolve（K,F） RESS=踟（1,M） forjeerent=1:M STRESS=SpaeeTrtmsElementStress（jelement） TITLE TITLE XGXG YGYG ZGZG FORCE FORcE NODEFIX NODEFIX ELEMDEFELEN伍）EF d/sp（AA （EE 赋1置0后K=K DISPDISP STRESSRESS 3结论 本程序经多个算例验证,都得到正确结果.应 用此程序,即可对任意条件的空间三维桁架结构进 行快速的力学分析,判断结

16、构的受力情况,为进一 步的工程设计奠定基础.因此,具有较强的工程实 用价值. 参考文献:1LoganD.AFirstCourseinFiniteElementMethodM. SecondEdition.PWSPublishingCompany,1993. 2SennettR.MatrixAnalysisofStructureM.Prentice Hau.1994. 3王瑁成.有限单元法M.北京:清华大学出版社, 2003.7. 4张志涌,徐彦琴.等,MATLAB教程基于6.x版本 M.北京:北京航空航天大学出版社,2001.4. Designof3DSpaceTrussSolverBased

17、onMATLAB LIGang （Aeronautics&A.tronauties,BdjingUniversity,Beijing100083,China） Abstract:Thefiniteelementmethodof3Dspacetrussforceanalysisisintroduced.Theprogramsettlesthe problemshoWtoassemblethestiffnessmatrixandhowtoaispotheconstrainbasedonmatlablanguage.It isallpurposeprogram.anditcancalculateallkindofforceanalysisproblemsabout3Dspacetruss.Thepro gramCAInbeveryusefulonprojects. Keywords:3Dtruss;Finiteelementmethod;Matlab;Anpurposeprogramdesign