有限元计算注意事项

1、精选优质文档-倾情为你奉上1、计算原理任何有限元模拟的第一步都是利用一个有限单元的集合离散结构的实际几何形状,每个单元（element）代表这个实际结构的一个离散部分.这些单元通过公用的节点（node）来连接.节点和单元的集合成为网格（mesh）.在一个特定网格中的单元数目称为网格密度（mesh density）.在ansys计算过程中,程序以每个节点的每个自由度建立平衡方程,以节点的位移作为未知量,利用矩阵求解节点的位移.一旦节点位移求出,整个结构的应力和应变都很容易计算出来.这种计算的过程和方法,数学上称之为隐式方法.从上叙述来看,整个计算过程中就是求解n个n元一次方程组（n表示节点数量）

2、,当计算模型复杂而且庞大时,隐式求解方法的计算量还是相当大的.与之相对应的,显式求解方法.显式求解方法是通过动态方法从一个增量步前推到下一个增量步得到的.具体显式求解方法和隐式求解方法例子如下：（1） 隐式求解（2） 显式求解隐式求解中,计算的精度完全控制于计算步数,在一般的计算软件中（flac、abaqus）,软件均是利用不平衡力来控制计算步数（当不平衡力<10-5时,停止计算）.不平衡力A+B.A表示施加在节点上的集中力；B表示：在n步数下,根据第n步计算出来的应力,求出节点的内力. Flac软件中 B,以上公式是根据虚功原理推倒而得到.具体推倒过程见flac 原理.2、Ansys计

3、算注意事项：计算单位、参数、荷载、标准值、设计值,计算过程中系数的加入.（1） beam单元对于beam单元.Ansys软件中我们常用的有两种梁单元：beam188和beam4.这两种单元均是三维的梁单元,每个节点都具有6个自由度（ux、uy、uz、mx、my、mz）,并且单元坐标系x轴是i点指向j点.Beam188单元是基于铁木辛哥理论的梁,beam4单元是我们常用的经典结构力学梁.（铁木辛哥理论考虑了梁的剪切变形,而我们常用的经典结构力学梁只考虑了弯矩对结构的变形影响）所以说,beam188可以更精确的计算梁单元,因此我们结构计算中,一般都采用beam188单元.当然还有beam189单元

4、,189单元属于三维二次的梁单元（beam188属于三维一次梁单元）,精度比beam188更加高.定义beam188单元,一般采用如下形式：!定义单元/prep7 !进入前处理et,1,beam188 !定义单元188号标号为1!定义材料属性mp,ex,1,2.55e7 !定义弹性模量(kn/m2)mp,nuxy,1,0.167 !定义泊松比mp,dens,1,2.5 !定义密度（KN/N*KG/M3）nummrg,all ！合并重合节点numcmp,all ！压缩编号!定义梁截面SECTYPE,1,BEAM,RECT,A1,0 ! 1表示梁编号 ; RECT表示是矩形梁（还有其他t型等等,具

5、体见ansys帮助）; A1 表示梁的名称 ; 0表示薄壁梁单元网格划分精细程度（05）.SECDATA,1,3,4,12 !1表示梁b ; 3表示梁h ; 4和12定义对应宽长等分份数.SECOFFSET,CENT !cent质心 ; shrc剪切中心 ; origin原始中心 ; user用户定义;!注意：当梁单元和壳单元一起使用时,可以设置梁单元的偏心,使梁的一面和壳的一面共面.（secoffset,user,offsety,offsetz）,如下图：!划分网格LSEL,S,1 !选中编号为1的线.LATT,1,1,1,1 !mp,r,et,方向点,SECTYPE截面号.LESIZE,A

6、LL,0.2,1 !0.2是单元大小,1是确认细分规则.LMESH,ALL !用beam单元离散模型,形成网格.!对于划分网格,空间的beam单元,由于需要确定b、h的方向,ansys软件利用方向点来控制b、h的方向.方向点的编号最好定义的很大,如果定义太小,会影响后面的加载.具体方向点如何控制见上面的latt命令和ansys帮助.自己试两下就知道怎么用了.AllsFINISH!加载加约束/SOLUACEL,9.8 !重力加速度.注意方向,数值和整体坐标相反,比如重力指向z轴负向,则为正值.SFCUM,ALL,ADD !设置单元荷载是叠加还是替代,只对加在单元和节点上的荷载有效,对于加在面、线

7、上的荷载,都只有替代作用（对同一个面,第二次加的荷载替代第一次加的荷载）!对于beam单元,只能根据sfbeam命令增加均布荷载 等大小的均布荷载.Lsel,s,1ESLL,S,1sfbeam,ALL,1,PRES,-161.5 !1表示作用在beam单元的面上(如下图,面表示beam单元的轴向,面表示单元侧面,面表示beam单元顶面),-161.5表示均布荷载大小,正负号可以控制作用力的方向. 梯形均布荷载Sfbeam命令是对每个单元进行加载.如果一根梁承受10100的梯形均布荷载,而且这根梁被分成了10个beam单元,这样施加荷载就非常困难.因此我将这种加载过程写成命令流,让软件自动进行加

8、载.命令流如下：LSEL,S,1 !选中要加载的那根梁（线）ESLL,S,1 !选中属于这根梁（线）的beam单元*GET,Nelem,ELEM,COUNT, , , , !获得当前所选单元个数,赋予参数Nelem*GET,Ne,ELEM,NUM,MIN, , , , !获得当前所选单元最小编号,赋予参数Ne*DO,I,1,Nelem !循环加载,循环次数单元个数ESEL,S,NeNSLE,S,1*GET,Nnode,NODE,COUNT, , , , !获得当前所选节点个数,赋予参数Nnode*GET,Nn,NODE,NUM,MIN, , , , !获得当前所选节点最小编号,赋予参数NnNN

9、1X=NX(NN) !将nn节点的x坐标赋予NN1X(NX表示x坐标,NY表示y坐标)NN=NDNEXT(NN) !NN=当前所选节点的下一个编号NN2X=NX(NN) !将nn节点的x坐标赋予NN2Xsfbeam,ALL,1,PRES,-1630.76/3.23*(3.23-NN1X),-1630.76/3.23*(3.23-NN2X)!以上荷载公式应根据实际情况进行调整LSEL,S,1ESLL,S,1NE=ELnext(NE) !NE=当前所选单元的下一个编号 *ENDDO!对于此命令流,根据不同的实际情况,ABC部分需要修改,其他不需要修改.!后处理 (XY平面) （大拇指指向y，就是m

10、y）etable,ImY,smisc,2 !显示弯距etable,JmY,smisc,15PLLS,IMY,JMYETABLE,IFX,SMISC,1 !显示轴力ETABLE,JFX,SMISC,14PLLS,IFX,JFXETABLE,IFY,SMISC,5 !显示剪力ETABLE,JFY,SMISC,18PLLS,IFY,JFY!注意：beam单元的结果输出都是以单元坐标系输出的,且拉为正、压为负.前面我们已经知道,单元坐标系x轴就是i点指向j点,其他坐标可以根据整体坐标系推出.详细内容见ansys 帮助.（2） Shell单元对于shell单元应用的范围,ansys软件并没有强制规定,只

11、是从字面上区分了薄壳和厚壳.我以前看过一本电子教案仿真在线,里面说一般规定壳体的主尺寸是厚度的10倍左右,都是可以用壳体来模拟的.一般高度与跨度之比（非与单元尺寸比较）<1/15,可以当作薄壳处理,>1/15 & <1/10,可以当作厚壳来处理.shell63是薄壳单元,他包含弯曲和薄膜效应,但是忽略横向剪切变形；shell43,shell143,shell181,shell91,shell93和shell99,都属于厚壳单元,不仅有弯曲、薄膜效应,他也包含了横向剪切效应横向剪力所引起的变形和弯曲变形属同一量级,中厚板应该考虑.横向剪切被表示为整个厚度上的常剪切应变.

12、这种一阶近似只适用于中等厚度壳体. 线形分析时,如果不包含横向剪切应变,使用63,163单元;如果横向剪切变形重要,则遵守以下原则:均匀材料,使用43,93,143单元,复合材料使用91,99,181.我们土木工程中,一般利用shell43计算.!定义单元/prep7 !进入前处理et,1,shell43 !定义单元43号标号为1!定义材料属性mp,ex,1,2.55e7 !定义弹性模量(kn/m2)mp,nuxy,1,0.167 !定义泊松比mp,dens,1,2.5 !定义密度（KN/N*KG/M3）!定义墙体厚度!等厚度板R,1,2 !1表示编号,2表示厚度（m）R,2,3Asel,s,

13、1 !选中1号面Aatt,1,1,1 !mp,real,typeESIZE,0.2 !定义单元大小为0.2左右MSHAPE,0,2D !规定划分单元形状,0表示四边形（1表示三角形）,2d表示划分面（3d表示划分体）MSHKEY,2 !指定是自由划分还是映射划分,2表示：尽量用映射划分,不符合要求就自动使用自由划分,具体参见ansys帮助的eshkey命令.AMESH,ALL !划分面单元.！注意：在网格剖面方面，最好全部用四边形，而且形状尽量规则、均匀！因为将来后处理内力提取的时候，提取出来的力和单元的大小有直接的关系。后处理部分将详细讲解。!渐变厚度的板对于渐变的板,首先选中要划分的面,然

14、后进行网格划分.但是在网格划分的时候,赋予的厚度编号开始选用100（前面没有定义的就可以）,然后利用以下命令进行厚度赋值.ASEL,S,1ESLAMXNODE = NDINQR(0,14)*DIM,THICK,MXNODE*DO,NODE,1,MXNODE*IF,NSEL(NODE),EQ,1,THENTHICK(node) = 3-(NZ(NODE)-7.5)/25.5*1.5 !厚度和点坐标系的关系,NZ(NODE)表示编号为node的点的Z坐标值.*ENDIF*ENDDONODE = $ MXNODE =RTHICK,THICK(1),1,2,3,4!注意:根据不同的实际情况,以上命令流

15、只需要修改ABC部分内容.!划分好网格后,有一个很重要的事情：重新调整单元坐标系,因为好多单元的单元坐标系是乱七八糟的,但是将来结果文件输出值的方向是单元坐标系的.因此必须把所有单元的单元坐标系调整的和整体坐标系一样(只是方向一样,不一定x、y、z对应,因为壳体的z轴永远是壳体的法向).命令如下：Asel,s,1 !选中1号面Esla,s,1 !选中属于1号面的所有单元Local,11,0,1,0,0 !11局部坐标系编号;0表示笛卡尔坐标系;1表示x值Emodif,all,esys,11 !修改单元坐标系为局部坐标系11注意:修改单元坐标系的时候,必须先定义一个局部坐标系.因为emodif命

16、令只能利用局部坐标系修改单元坐标系.!将壳体法向反方向：有时候需要调整壳体法向的时候,可以利用如下命令：Asel,s,1AREVERSE,ALL !将面的方向反过来,ansys会提示你：是否也要将属于面的单元法向也反过来?AllsFINISH!加载加约束/SOLUACEL,9.8 !施加重力加速度,详细解释见上.等大小的均布荷载 对壳体施加等大小的均布荷载比较简单,以下命令均可以实现,只要你选中相应的面、点或者单元.Sfa,ALL,1,PRES,100 Sf,ALL,PRES,100 !个人不推荐此命令,因为此命令对施加在实体面上的节点有限制.Sfe,ALL,1,PRES,0,100 !主力推

17、荐此命令.1表示作用壳体的哪个面（如下图所示 .）;0表示确实存在的压力（SURF153 or SURF154单元选 用2）;100表示压力(kn/m2).梯形荷载SFGRAD,PRES,0,Z,5,-10*1.2 !为下面的加载形成梯度。pres表示施加压力；0表示 坐标系代号(笛卡尔坐标系）；z表示梯度沿z方向；5 表示：在0坐标系下z5的位置(下面命令的300就加 在这个位置）；-10*1.2表示施加的梯度，负号表示随 着z的增大，压力在减小.SFe,ALL,2,PRES,0,300 !见上sfe命令解释.!后处理（如下图）!壳体单元输出的是每个单元在单元坐标系下，单位长度的值。(这句话

18、是ansys帮助里翻译而来的）有点含糊哈，大家凑合理解。我举个列子给你们听听。如下例题： 计算图 选取单元图求基础和板之间的剪力？用屁股想都想的出来，精确解是100kn。如此简单的问题，我们选用shell63做计算（当然我们也可以选用平面单元，这里只是为了试验shell63求解结果的含义）对这个模型，我们分两种情况计算，两种情况唯一不同之处就是网格大小，第一种情况网格大小1*1,第二种情况网格大小0.5*0.5。取出靠近地面的最近的一排单元（如上图），利用单元table输出每个单元的剪力值（单元坐标系下），然后相加。第一种情况，相加后A99.8 kn第二种情况，相加后A198.9 kn由上可得

19、:剪力值A*单元大小这里可以看出，划分单元尽量规则、均匀的重要性，特别是将来你要提取内力的位置，在划分网格的时候就要做好准备！而且计算前单元坐标系要统一也非常重要！本人编写了一个提取内力的通用程序，各位可以参考参考。只要选中你要提取内力的那排单元（如同上图的那排），然后输入以下命令即可：JLzonghe=0 !剪力ZLZONGHEx=0 !X方向轴力ZLZONGHEy=0 !Y方向轴力mzonghex=0 !X方向弯矩mzonghey=0 !Y方向弯矩etab,JL,smisc,3ETAB,ZL2,SMISC,2ETAB,ZL1,SMISC,1etab,m1,smisc,4etab,m2,sm

20、isc,5*GET,MAXN,ELEM,0,NUM,MAX !所选单元的最大编号给maxn*GET,MINN,ELEM,0,NUM,MIN*DO,I,MINN,MAXN *IF,ESEL(I),EQ,1,THEN !判断单元号I是否属于所选范围 *get,jianli,elem,i,etab,JL *get,ZHOULIx,elem,i,etab,ZL1 *get,ZHOULIy,elem,i,etab,ZL2 *get,wanjux,elem,i,etab,m1 *get,wanjuy,elem,i,etab,m2 JLzonghe=JLzonghe+jianli ZLzonghex=ZLz

21、onghex+ZHOULIx ZLzonghey=ZLzonghey+ZHOULIy mzonghex=mzonghex+wanjux mzonghey=mzonghey+wanjuy *ENDIF*ENDDO注意：（1）以上命令求解出来的 Jlzonghe 等等还没有乘以单元大小。（2）以上求解出来的X方向、Y方向的轴力和弯矩都是单元坐标系下的值（具体方向如下图,SX、SY就是单元坐标系），单元坐标系的X、Y方向可以通过视图显示出来，命令如下：/psymb,esys,1蓝线表示Z轴，白线表示X轴，青线表示Y轴（3） Solid单元实体单元其实是ansys最简单的单元，但是也是ansys最难操

22、作的单元。实体单元从建模到形成计算书，70的时间用在划分网格上，20的时间用在后处理上，还有10就是一些小杂事了。前期操作方面，有两种方法可以完成：（1）利用ansys建模同时划分网格。我有一个很牛的同事，做窝壳方面很有研究，他曾经做过一个工程的窝壳，写了整整2万行的命令流，那窝壳真是漂亮，不得不让人佩服。不过这种方法需要非常非常熟悉ansys，不是一般人可以的哦（汗颜！我为什么不行，55555）。（2）利用三维软件建模(个人推荐solidwork,和ansys的兼容性非常ok)，然后导入ansys。这种方法比较通用，但是一般solidwork导入的都是一个整体块，需要人工不停的用wp（工作面

23、）切割，切到可以映射划分为止。这个切割技术还是比较难的，需要自己多去试试。一定要自己动手去切割、去划分，只要认真做两个例题，你就能自学成才了。上面废话一堆，实体单元的前处理我就不说了，重点说一下我们土建专业比较关注的后处理：内力（剪力、轴力、弯矩）提取。由于solid单元最后输出的都是应力，无法在工程中利用（当然也可以根据规范提到的应力配筋，但个人感觉不太稳妥，我曾经比较过：同一截面，应力配筋和提取内力的结果，内力配筋结果都大于应力配筋），但是我们可以利用这些应力积分得到我们所需要的内力。积分的原理就不说了，哪个兄弟要是忘了就去再看看结构力学书吧。首先将工作平面的xy面移到你要计算内力的位置，

24、然后你要判断：（1） 你准备计算内力的截面属于整体坐标系的xy平面SUCR,SUZ2,CPLANE,3 !移动工作平面、创建面SUZ2 SUMAP,MYSX,S,X$SUMAP,MYSY,S,Y !映射X和Y方向应力 SUMAP,MYSZ,S,Z$SUMAP,MYSXY,S,XY !映射Z和XY方向应力 SUMAP,MYSYZ,S,YZ$SUMAP,MYSXZ,S,XZ !映射YZ和XZ方向应力 SUEVAL,XFORCE,MYSXZ,INTG !求截面上FXSUEVAL,YFORCE,MYSYZ,INTG !求截面上FYSUEVAL,ZFORCE,MYSZ,INTG !求截面上FZSUEVA

25、L,MYA,DA,SUM !求截面面积并赋给MYA变量 SUEVAL,MYYA,GCY,INTG !求关于X轴的面积矩并赋给变量MYYA MYYA=MYYA/MYA !得到面积重心到X轴的距离=面积矩/面积 SUEVAL,MYXA,GCX,INTG !求关于Y轴的面积矩并赋给变量MYXA MYXA=MYXA/MYA !得到面积重心到Y轴的距离=面积矩/面积 SUCALC,SZGCY,MYSZ,MULT,GCY !计算MYSZ×GCY，并赋给面项SZGCY SUEVAL,MX1,SZGCY,INTG !对面项SZGCY在面上积分得到MX1 SUCALC,SZGCX,MYSZ,MULT,

26、GCX !计算MYSZ×GCZ，并赋给面项SZGCX SUEVAL,MY1,SZGCX,INTG !对面项SZGCX在面上积分得到MY1 !上述弯矩基于总体直角坐标系原点而言的，应对面积重心取矩，将内力简化到面积重心上 MX1=MX1-ZFORCE*MYYA !MY1=MY1-ZFORCE*MYXA !（2） 你准备计算内力的截面属于整体坐标系的yz平面SUCR,SUZ2,CPLANE,3 !移动工作平面、创建面SUZ2 SUMAP,MYSX,S,X$SUMAP,MYSY,S,Y !映射X和Y方向应力 SUMAP,MYSZ,S,Z$SUMAP,MYSXY,S,XY !映射Z和XY方向

27、应力 SUMAP,MYSYZ,S,YZ$SUMAP,MYSXZ,S,XZ !映射YZ和XZ方向应力 SUEVAL,XFORCE,MYSx,INTG !求截面上FXSUEVAL,zFORCE,MYSXZ,INTG !求截面上FzSUEVAL,yFORCE,MYSXY,INTG !求截面上FySUEVAL,MYA,DA,SUM !求截面面积并赋给MYA变量 SUEVAL,MYZA,GCZ,INTG !求关于X轴的面积矩并赋给变量MYYA MYZA=MYZA/MYA !得到面积重心到Y轴的距离=面积矩/面积 SUEVAL,MYYA,GCY,INTG !求关于Z轴的面积矩并赋给变量MYYA MYYA=

28、MYYA/MYA !得到面积重心到Z轴的距离=面积矩/面积 SUCALC,SZGCZ,MYSX,MULT,GCZ !计算MYSX×GCZ，并赋给面项SZGCZ SUEVAL,MY1,SZGCZ,INTG !对面项SZGCZ在面上积分得到MY1 SUCALC,SZGCY,MYSX,MULT,GCY !计算MYSX×GCY，并赋给面项SZGCY SUEVAL,MZ1,SZGCY,INTG !对面项SZGCY在面上积分得到MZ1 !上述弯矩基于总体直角坐标系原点而言的，应对面积重心取矩，将内力简化到面积重心上 MY1=MY1-XFORCE*MYZA !MZ1=MZ1-XFORCE

29、*MYYA !（3） 你准备计算内力的截面属于整体坐标系的xz平面SUCR,SUZ2,CPLANE,3 !移动工作平面、创建面SUZ2 SUMAP,MYSX,S,X$SUMAP,MYSY,S,Y !映射X和Y方向应力 SUMAP,MYSZ,S,Z$SUMAP,MYSXY,S,XY !映射Z和XY方向应力 SUMAP,MYSYZ,S,YZ$SUMAP,MYSXZ,S,XZ !映射YZ和XZ方向应力 SUEVAL,XFORCE,MYSxy,INTG !求截面上FXSUEVAL,zFORCE,MYSYZ,INTG !求截面上FzSUEVAL,yFORCE,MYSy,INTG !求截面上FySUEVA

30、L,MYA,DA,SUM !求截面面积并赋给MYA变量 SUEVAL,MYZA,GCZ,INTG !求关于X轴的面积矩并赋给变量MYYA MYZA=MYZA/MYA !得到面积重心到X轴的距离=面积矩/面积 SUEVAL,MYXA,GCX,INTG !求关于Z轴的面积矩并赋给变量MYXA MYXA=MYXA/MYA !得到面积重心到Z轴的距离=面积矩/面积 SUCALC,SZGCZ,MYSy,MULT,GCZ !计算MYSY×GCZ，并赋给面项SZGCZ SUEVAL,MX1,SZGCZ,INTG !对面项SZGCZ在面上积分得到MX1 SUCALC,SZGCX,MYSy,MULT,

31、GCX !计算MYSY×GCX，并赋给面项SZGCX SUEVAL,MZ1,SZGCX,INTG !对面项SZGCX在面上积分得到MY1 !上述弯矩基于总体直角坐标系原点而言的，应对面积重心取矩，将内力简化到面积重心上 MX1=MX1-yFORCE*MYZA !MZ1=MZ1-yFORCE*MYXA !注意：（1） XFORCE、YFORCE、ZFORCE、MX1、MY1、MZ1为最后的计算结果，而且是在整体坐标系下的结果。（2） 对于xy平面，ZFORCE为该断面的轴力；XFORCE、YFORCE分别表示x、y方向的剪力；MX1、MY1分别表示右手准则，大拇指分别指向整体坐标系下的

32、x、y轴。其他yz、xz平面类似。（3） 这里强调一下求得内力的方向问题。我们求解内力的时候，其实就是对各应力进行面积积分，所得内力正负值的意义和各应力一样。在ansys中各应力的正负的意义是和有限元定义正负号的方法是一致的，和单纯的整体坐标系的正负号是不一样的。举例如下： 比如你求一个属于xy平面的截面的内力，求出剪力100kn。那么对于这个剪力方向是什么呢？根据有限元正负号定义原则！可参照有限元原理母校教材哈，徐之伦写的经典的。实在不想看，就q我，我qq：。（4） 以前我一个同事发现了一个命令，很简单的一个命令输入，然后就可以输出以上结果。但是对于此命令得到结果，本人不太明白他输出的意思，