三维问题有限元分析(包括轴对称问题).ppt

第六章三维问题有限元分析,讲授：陈得良TelQ：416501065Email：deliang_chen,1,四教学基本内容,第五章三维问题有限元分析第一节三维应力状态第二节4节点四面体单元第三节8节点六面体等参单元第四节20节点等参单元第五节ansys空间问题实例第六节空间轴对称问题有限元法第七节Ansys轴对称旋转问题实例,2,工程实际中的很多问题难于简化为平面问题，如受任意空间载荷作用的任意形状几何体，受对称于轴线载荷作用的回转体,这类问题经典弹性力学往往无能为力。在FEM中，空间问题只要求0阶连续，因此构造单元方便,空间问题简介,3,空间问题的主要困难：（1）离散化不直观；（网格自动生成）（2）分割的单元数量多，未知量的数目剧增。（对某些问题简化）（轴对称问题）,空间分析的优点精确,4,6.1三维应力状态,工程结构一般都是空间的弹性体。受力作用后，其内部各点将沿x、y、z坐标轴方向产生位移，是三维空间问题，其应力状态如图6-1所示。,图6-1空间结构应力状态,各点沿x、y、z方向的位移以u、v、w表示，这些位移为各点坐标的函数，即：,u=u(x、y、z)v=v(x、y、z)w=w(x、y、z),5,由弹性力学知，应变与位移间的几何关系是,(6-1),三维弹性体的应变分量，用矩阵表示为,(6-2),6,弹性体受力作用，内部任意一点的应力状态也是三维的，用列向量表示为,在线弹性范围内，应力与应变间的物理关系矩阵表达式为,对于各向同性弹性体，在三维应力状态下，弹性矩阵的形式为,(6-3),(6-4),7,(1)空间问题常用单元：四面体单元、长方体单元、直边六面体单元、曲边六面体单元、轴对称单元。,4结点四面体单元：是空间问题最简单的单元，也是常应变、常应力单元，可以类似平面问题三结点三角形单元进行分析。8结点长方体单元：可以类似平面四结点矩形单元进行分析。8结点直边六面体单元：可以类似平面四结点任意四边形等参元分析。20结点曲边六面体单元：等参单元,可以类似平面八结点曲边四边形等参元进行分析。轴对称单元：一平面单元绕一对称轴旋转形成的空间问题。只需在rz平面划分网格，就像平面问题xy平面中的网格一样，这样这类空间问题可以得到简化。（环向位移等于零）,（2）结点位移3个分量。（3）基本方程比平面问题多。3个平衡方程，6个几何方程，6个物理方程。,8,6.2四节点四面体单元,图6-2表示任一简单四面体单元，其中四个结点编号设为i、j、m、n(或1、2、3、4)。单元变形时，各结点沿x、y、z方向上的位移，以列向量表示为,图6-2四面体单元,8,单元变形时，单元内各点也有沿x、y、z方向的位移u、v、w，一般应为坐标x、y、z的函数。对于这种简单的四面体单元，其内部位移可假设为坐标的线性函数，为满足变形协调条件，取为,(6-5),式(6-5)含有12个待定系数a，可由单元的12项结点位移决定.将4个结点的坐标值代入式(6-5)的u式中。i、j、m、n共4个结点，分别有,(6-6),1单元形函数,9,其中,式中，V为四面体的体积，且有,(6-7),由式(6-6)求出，再代回式(6-5)中，整理后得,10,为使四面体的体积V不为负值，在右手坐标系中，使右手旋转按着由i-j-m的转向转动时，且法向n方向前进。用求位移u的同样方法，可求得,将位移的3个线性方程形成的线性方程组用矩阵表示为,(6-8),式中,(6-9),11,2单元刚度矩阵,将式(6-8)代入几何方程式(6-2)，经过微分运算，可得单元内应变为,(6-10),式中,(6-11),简单四面体单元内，各点的应变都是一样的，这是一种常应变单元（是三维单元中精度最低的单元）。这一点与平面问题的简单三角形单元相似，由于单元内位移都假定为线性变化的，因而由位移一阶导数组成的应变也为常量。,12,同样，用虚功原理建立结点力和结点位移间的关系式，从而得出简单四面体单元的刚度矩阵。,(6-12),(6-13),按结点分块表示，此单元刚度矩阵可表示为,(6-14),13,(r=i、j、m、n，S=i、j、m、n)(6-15),式中,弹性体三维（空间）问题的原始平衡方程组，即,其中,其中任一子矩阵为,14,3整体结构载荷列向量,整体结构的结点载荷列向量,(6-16),式中,单元上集中力等效结点载荷列向量；单元上表面力等效结点载荷列向量；单元上体积力等效结点载荷列向量；单元结点载荷列向量。,等效结点力公式为,式中,15,6.38节点六面体等参单元,8(x5,y5,z5),7,如同二维等参单元一样，三维等参单元的有关公式的建立也是采用局部自然坐标（曲面坐标）。可参考的母体单元则为一正六面体，图6-3则表示了任意六面体单元与母体单元、局部的三维自然坐标与整体的直角坐标系的几何关系,母体单元,任意六面体单元（对面不全平行）,图6-38结点三维等参单元,16,用形函数表示的位移插值形式的位移模式，可以直接利用拉格朗日插值公式，得到单元位移函数为,根据等参单元的定义。自然坐标与整体直角坐标之间的关系可以写为,其中形函数为,例,其中为节点坐标值（角点）,17,由节点位移求单元应变时，他要求形函数在整体坐标下的导数，但形函数是建立在局部坐标下的，这就需要将局部坐标中的表达式转换到整体坐标系中，如同平面等参单元一样，需要通过雅克比矩阵来实现，由偏导法则,同理可得,写成矩阵,18,求单元刚度矩阵，尚需对积分的单元体积进行积分变换,反之,有了上式很容易得到单元的应变应力矩阵,19,6.420结点等参元,为适应三维结构的曲面边界，可以采用曲面六面体单元。正方体基本单元内任一点与实际曲面单元内的点一一对应，结点也一一对应。这里，实际单元边界线中间的结点9、10、20，都“映射”成为正方体的棱边中点。,8结点单元是线性单元，其位移模式是三维线性的，在8结点单元的基础上每边增加一个中点作为节点就构成了20节点单元。此时六面体单元每条边上有3个节点，他们既可以是直线的，也可以是曲线的，因此每个面也可以是平面的，也可以是曲面的。,1形状函数,20,(a)直角坐标系与实际单元(b)自然坐标系与基本单元,图6-320结点三维等参单元,位移函数和几何坐标的变换式应取为相同的参数，其坐标变换关系可表示为,(6-17),则单元的位移函数可写成,21,(6-18),在自然坐标系（局部坐标系）中，各结点的形状函数可写成如下形式，对于8个顶角结点（i1，2，8）,式中xi、yi、zi结点i的坐标；ui、vi、wi结点i沿x、y、z方向的位移；Ni对应于i结点的形状函数。,22,对于的边上点（i17，18，19，20）,(6-19),对于的边上点（i9，11，13，15）,对于的边上点（i10，12，14，16）,23,2单元刚度矩阵,三维变形状态下，一点的应变与位移的几何关系为,(6-20),24,(6-22),为便于以下计算，弹性矩阵D可分块写为,(6-23),令,则，,为6060的方阵，可按结点写为子块形式,25,式中第i行j列的子矩阵为,(6-24),将(6-20)、(6-22)分块式代入(6-23)，其被积函数可写为,(6-25),式中,与式(5-9)相似，按坐标变换式(6-17)，应有,26,同样可有,(6-26),三维六面体的雅可比矩阵为,(6-27),同理可采用三维高斯求积公式计算单元刚度矩阵。即,27,式中，L，M，N为沿、方向的积分点数目，而积分点坐标及权重可由高斯积分表查得。,对于20节点的三维单元，通常可取积分点数目m=3(即333），查表可得对应的积分点坐标和权重为,28,参照方法，可以很方便的得到另外两个方向的积分点坐标值和权重。,29,6.5ANSYS空间问题示例,1问题描述,如图6-4所示，一个圆柱实体。柱高0.2m，圆柱横截面直径为0.1m。约束方式：底面全约束。承受载荷：A点承受Z方向集中载荷Fz=5000N和Y方向集中载荷Fy=-5000N；B点承受X方向集中载荷Fx=5000N；C点承受Z方向集中载荷Fz=-5000N；D点承受X方向集中载荷Fx=-5000N。弹性模量为EX=210GP，泊松比=0.3,2ANSYS求解操作过程,30,（1）选择单元类型运行PreprocessorElementTypeAdd/Edit/Delete，弹出ElementTypes对话框，如图6-5所示。然后单击Add，弹出LibraryofElementTypes窗口，如图6-6所示，选择SOLID45单元，单击OK。,图6-5单元类型对话框,图6-6单元类型库对话框,31,（2）设置材料属性运行PreprocessorMaterialPropsMaterialModels，弹出如图6-7所示对话框。双击Isotropic，弹出LinearIsotropicPropertiesforMaterialNumber1对话框，如图6-8所示，在EX选项栏中设置数值2.1e11，在PRXY选项栏中设置数值0.3。设置完毕单击OK按钮。,图6-7选择材料属性对话框,图6-8设置材料属性对话框,32,（3）建立模型运行PreprocessorModelingCreateAreasRectangleBy2Corners，弹出如图6-9所示对话框，在WPX选项栏中填写0，在WPY选项栏中填写0，在Width选项栏中填写0.05，在Height选项栏中填写0.2，点击OK。生成如图6-10所示图形。,图6-9两点建立矩形对话框,图6-10生成的长方形面,33,将长方形旋转成柱体，运行PreprocessorModelingOperateExtrudeAreasAboutAxis，弹出如图6-11所示拾取框。选择图7中长方形后弹出单击OK，再选择长方形左上角和左下角结点后，单击OK.。弹出如图6-12所示对话框。在ARC选项栏中填入旋转角度360度，设置完毕单击OK按钮，生成如图6-13所示圆柱体。,图6-11拾取对称轴对话框,图6-12设置绕轴旋转参数对话框,图6-13圆柱模型,34,运行MeshingSizeCntrlsManualSizeGlobalSize弹出如图6-14所示对话框，设置SIZE选项栏中的数据为0.01。运行MeshingMeshVolumesFree自由划分网格后得到如图6-15所示图形。,图6-14设置网格尺寸对话框,图6-15圆柱有限元模型,（5）施加约束运行SolutionDefineLoadsApplyDisplacementOnAreas，拾取圆柱的底面，施加全约束。,35,（6）施加载荷显示图形的关键点，运行PlotCtrlsNumbering弹出如图6-16所示对话框，激活KPNumbers后面的选框，使它变成on形式。选择菜单SolutionDefineLoadsApply,StructureForce/MomentOnKeypoints，载荷分别如下：8点承受Z方向集中载荷Fz=5000N和Y方向集中载荷Fy=-5000N；10点承受X方向集中载荷Fx=5000N；3点承受Z方向集中载荷Fz=-5000N；6点承受X方向集中载荷Fx=-5000N。施加载荷，图形如图6-17所示。,图6-16编号显示设置对话框,36,图6-17圆柱实体示意图,（7）求解选择SolutionSolveCurrentLS，开始计算，计算结束会弹出计算完毕对话框，单击Close。关闭对话框计算完毕。,（8）后处理运行GeneralPostprocPlotResultsContourPlotNodalSolu。弹出如图6-18所示对话框，运行DOFSolutionDisplacementvectorsum和StressvonMisesstress，分别显示圆柱体的位移和应力云图。,37,图6-18云图显示对话框,结果显示如图6-19和图6-20所示。,38,图6-19位移云图图6-20应力云图,39,6.6空间轴对称问题的有限元法,对空间轴对称问题，常采用圆柱坐标系。r表示径向坐标，z表示轴向坐标，任一对称面为rz面。在有限元分析时，可采用轴对称的环形单元进行。环形单元可以是任何平面单元。,某一平面图形绕平面上某一轴旋转形成的回转体称为轴对称物体，此平面称为子午面。在动力机械，特别是叶轮机械中，有很多零件都具有轴对称特性，比如轮盘、旋转轴、承力环等。对于直齿圆柱齿轮，由于齿的存在，严格地说它并非轴对称物体。如果忽略齿的部分(将齿用外载荷表示)，则所得到的齿根以内的旋转体部分为轴对称物体。,轴对称物体的变形及应力分布不一定是轴对称的，只有当其约束和载荷都对称于旋转轴时，轴对称物体的变形和应力分布才是轴对称的。轴对称物体+轴对称约束+轴对称载荷=轴对称系统对轴对称系统的应力分析=轴对称物体,40,1）几何形状关于轴线对称；2）作用于其上的载荷关于轴线对称。3）约束条件关于轴线对称。因过z轴的任一子午面都是对称面，其上任一点p只在该平面上发生位移，即弹性体内任一点的位移、应力与应变只与坐标r、z有关，与无关。从而，轴对称问题可转化为二维问题，但因与平面问题有区别，常称为二维半问题。,柱坐标系,41,注意：应变虽然与无关，但是周向应变，周向应力，由径向位移引起，因为径向位移会导致周长的改变。,1、基本方程,位移分量,应力分量,应变分量,42,虚功方程,应变分量,轴对称问题的弹性矩阵：,43,2、轴对称问题的离散化,对于轴对称问题，利用其轴对称特性，在对其进行网格划分时可知取任意通过Z轴的截面进行，类似平面问题的网格形式。,本节以三角形单元为例。,1、位移模式,轴对称问题的环向位移恒等于零，径向r位移与轴向z位移不等于零。对于图示情形，依照平面问题的三角形单元分析，取位移模式为,代入结点位移后，可解出a1-a6，再代入上式，得,xr,yz,44,其中形函数：,；,单元中位移,根据弹性力学理论，空间轴对称问题的几何方程为,2、单元中应变,45,将u,w表达式代入上式，整理后,46,式中,其中,B矩阵中含有变量r，z，因此它不是常数矩阵，即轴对称问题的三角形环形单元不是常应变单元。,47,3、单元中应力,根据弹性力学理论，空间轴对称问题的应力-应变关系为,弹性矩阵:,48,单元中任意一点的应力：,4、单元刚度矩阵,由于被积函数与无关，故在三角形截面的环单元的积分可简化为在三角形截面上的积分。故有：,49,单元刚度矩阵的积分参照图示分区，按下式采用数值积分的方法进行,50,当单元较小时，可把各个单元中的r，z近似看作常数，并且分别等于各单元形心的坐标，即,这样，就可把各个单元近似地当做常应变单元,51,单元刚度矩阵k的分块形式,其中的近似子矩阵为,52,5、等效结点荷载,类似平面问题。对于作用于三角形环单元上的体积力、表面力的等效结点力为：,体力,53,面力：,（1）均布表面力设单元ij边上作用均布表面力，其集度为,l,当ri=rj时，静力等效原则,54,(2)三角形分布表面力沿单元ij边作用了三角形分布的表面力，表面力在i点集度为,当ri=rj时，静力等效原则。2/3集中在i点，1/3集中在j点。,55,56,圆筒直径0.4m，高度0.6m，壁厚0.005m；材料Q235，弹性模量E=2.1e11Pa，泊松比=0.3；约束：圆筒的下部在轴线方向固定，其它方向自由；载荷：顶部环线上承受轴向线压力P-200000N/m。,图6-4圆筒示意图,图6-5单元类型对话框,1问题描述,6.7ANSYS轴对称旋转单元计算示例,57,（1）选择单元类型运行PreprocessorElementTypeAdd/Edit/Delete，弹出ElementTypes对话框单击Add，弹出LibraryofElementTypes对话框，如图7-6所示，选择SHELL51单元。,2ANSYS求解操作过程,图7-6单元类型库对话框,图7-7选择材料属性对话框,58,（2）设置材料属性运行PreprocessorMaterialPropsMaterialModels，弹出DefineMaterialModelBehavior对话框，如图7-7所示。双击Isotropic选项，弹出LinearIsotropicPropertiesforMaterialNumber1对话框，如图7-8所示。,图7-8设置材料属性对话框,59,（3）定义单元实常数选择MainMenuPreprocessorRealConstantsAdd/Edit/Delete，弹出如图7-9所示对话框，单击Add按钮弹出ElementTypeforRealConstants对话框，如图7-10所示，选择Type1SHELL51，单击OK，弹出RealConstantSetNumber1,forSHELL51对话框，如图7-11所示，在TK(I)项输入0.005，单击OK。,图7-9实常数对话框图7-10选择要设置实常数的单元类型,60,图7-11设置SHELL51实常数对话框,（4）建立模型首先生成关键点，运行主菜单PreprocessorModelingCreateKeypointsInActiveCS，弹出如图7-12所示对话框。创建关键点1（0.2,0,0），2（0.2,0.6,0）。生成圆筒母线：运行MainMenuPreprocessorModelingCreateLinesLinesStraightLine，弹出拾取关键点对话框，拾取关键点1、2，单击OK。,61,图7-12创建关键点对话框,（5）设置单元属性运行MainMenuPreprocessorMeshingMeshTool，弹出MeshTool对话框，在ElementAttributes下拉列表中选择Lines，然后单击其后的Set按钮弹出拾取线对话框，单击PickAll，弹出分配线单元属性对话框，将MAT,TEAL,TYPE依次设置为1，1，1，单击OK。,（6）划分网格,62,单击MeshTool中Lines后的Set按钮，弹出拾取线对话框，单击PickAll弹出控制线单元尺寸对话框，将NDIV设置为10，单击OK。在MeshTool对话框中的Mesh下拉列表中选择Lines单击Mesh，弹出拾取线对话框，单击PickAll，划