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

1、第六章 三维问题有限元分析,讲 授：陈得良 TelQQ：416501065 Email：deliang_,1,四 教学基本内容,第五章 三维问题有限元分析 第一节 三维应力状态 第二节 4节点四面体单元 第三节 8节点六面体等参单元 第四节 20节点等参单元 第五节 ansys空间问题实例 第六节 空间轴对称问题有限元法 第七节 Ansys轴对称旋转问题实例,2,工程实际中的很多问题难于简化为平面问题，如受任意空间载荷作用的任意形状几何体，受对称于轴线载荷作用的回转体,这类问题经典弹性力学往往无能为力。在FEM中，空间问题只要求0阶连续，因此构造单元方便,空间问题简介

2、,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

3、),6,弹性体受力作用，内部任意一点的应力状态也是三维的，用列向量表示为,在线弹性范围内，应力与应变间的物理关系矩阵表达式为,对于各向同性弹性体，在三维应力状态下，弹性矩阵 的形式为,(6-3),(6-4),7,(1)空间问题常用单元：四面体单元、长方体单元、直边六面体单元、曲边六面体单元 、轴对称单元。,4结点四面体单元：是空间问题最简单的单元，也是常应变、常应力单元，可以类似平面问题三结点三角形单元进行分析。 8结点长方体单元：可以类似平面四结点矩形单元进行分析。 8结点直边六面体单元：可以类似平面四结点任意四边形等参元分析 。 20结点曲边六面体单元：等参单元,可以类似平面八结点曲边四边

4、形等参元进行分析 。 轴对称单元：一平面单元绕一对称轴旋转形成的空间问题。只需在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的函数。对于这种简单的四面体单

5、元，其内部位移可假设为坐标的线性函数，为满足变形协调条件，取为,(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、,将式(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),式中,弹性体三维（空间）问题的原始平衡方程组，即,

7、其中,其中任一子矩阵为,14,3 整体结构载荷列向量,整体结构的结点载荷列向量,(6-16),式中,单元上集中力等效结点载荷列向量； 单元上表面力等效结点载荷列向量； 单元上体积力等效结点载荷列向量； 单元结点载荷列向量。,等效结点力公式为,式中,15,6.3 8节点六面体等参单元,8 (x5,y5,z5),7,如同二维等参单元一样，三维等参单元的有关公式的建立也是采用局部自然坐标（曲面坐标）。可参考的母体单元则为一正六面体，图6-3则表示了任意六面体单元与母体单元、局部的三维自然坐标与整体的直角坐标系的几何关系,母体单元,任意六面体单元 （对面不全平行）,图6-3 8结点三维等参单元,16,

8、用形函数表示的位移插值形式的位移模式，可以直接利用拉格朗日插值公式，得到单元位移函数为,根据等参单元的定义。自然坐标与整体直角坐标之间的关系可以写为,其中形函数为,例,其中 为节点坐标值（角点）,17,由节点位移求单元应变时，他要求形函数在整体坐标下的导数，但形函数是建立在局部坐标下的，这就需要将局部坐标中的表达式转换到整体坐标系中，如同平面等参单元一样，需要通过雅克比矩阵来实现，由偏导法则,同理可得,写成矩阵,18,求单元刚度矩阵，尚需对积分的单元体积进行积分变换,反之,有了上式很容易得到单元的应变应力矩阵,19,6.4 20结点等参元,为适应三维结构的曲面边界，可以采用曲面六面体单元。正方

9、体基本单元内任一点与实际曲面单元内的点一一对应，结点也一一对应。这里，实际单元边界线中间的结点9、10、20，都“映射”成为正方体的棱边中点。,8结点单元是线性单元，其位移模式是三维线性的，在8结点单元的基础上每边增加一个中点作为节点就构成了20节点单元。此时六面体单元每条边上有3个节点，他们既可以是直线的，也可以是曲线的，因此每个面也可以是平面的，也可以是曲面的。,1 形状函数,20,(a)直角坐标系与实际单元 (b) 自然坐标系与基本单元,图6-3 20结点三维等参单元,位移函数和几何坐标的变换式应取为相同的参数，其坐标变换关系可表示为,(6-17),则单元的位移函数可写成,21,(6-1

10、8),在自然坐标系（局部坐标系）中，各结点的形状函数可写成如下形式， 对于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

11、列的子矩阵为,(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(即3 3 3），查表可得对应的积分点坐标和权重为,28,参照方法，可以很方便的得到另外两个方向的积分点坐标值和权重。,29,6.5 ANSYS空间问题示例,1 问题

12、描述,如图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,2 ANSYS求解操作过程,30,（1）选择单元类型 运行PreprocessorElement TypeAdd/Edit/Delete，弹出Element Types对话框，如图6-5所示。然后单击 Add，弹出Library of Elem

13、ent Types窗口，如图6-6所 示，选择SOLID45单元，单击OK。,图6-5 单元类型对话框,图6-6 单元类型库对话框,31,（2）设置材料属性 运行PreprocessorMaterial PropsMaterial Models，弹出如图6-7所示对话框。双击Isotropic，弹出Linear Isotropic Properties for Material Number1对话框，如图6-8所示，在EX选项栏中设置数值2.1e11，在PRXY选项 栏中设置数值0.3。设置完毕单击OK按钮。,图6-7 选择材料属性对话框,图6-8 设置材料属性对话框,32,（3）建立模型 运

14、行PreprocessorModelingCreateAreasRectangleBy 2 Corners，弹出如图6-9所示对话框，在WP X选项栏中填写0，在WP Y选项栏中填写0，在Width选项栏中填写0.05，在Height选项栏中填写0.2，点击OK。生成如图6-10所示图形。,图6-9 两点建立矩形对话框,图6-10 生成的长方形面,33,将长方形旋转成柱体，运行 PreprocessorModeling Operate ExtrudeAreasAbout Axis，弹出如图6-11所示拾取框。选择图7中长方形后弹出单击OK，再选择长方形左上角和左下角结点后，单击OK.。弹出如图

15、6-12所示对话框。在ARC选项栏中填入旋转角度360度，设置完毕单击OK按钮，生成如图6-13所示 圆柱体。,图6-11 拾取对称轴对话框,图6-12 设置绕轴旋转参数对话框,图6-13 圆柱模型,34,运行MeshingSize CntrlsManualSizeGlobalSize弹出如图6-14所示对话框，设置SIZE选项栏中的数据为0.01。运行MeshingMeshVolumesFree自由划分网格后得到如图6-15所示图形。,图6-14 设置网格尺寸对话框,图6-15 圆柱有限元模型,（5）施加约束 运行SolutionDefine LoadsApplyDisplacementOn

16、 Areas，拾取圆柱的底面，施加全约束。,35,（6）施加载荷 显示图形的关键点，运行PlotCtrlsNumbering弹出如图6-16所示对话框，激活KP Numbers后面的选框，使它变成on形式。选择菜单SolutionDefine LoadsApply,Structure Force/Moment On Keypoints，载荷分别如下：8点承受Z方向集中载荷Fz=5000N和Y方向集中载荷Fy=-5000N；10点承受X方向集中载荷Fx=5000N；3点承受Z方向集中载荷Fz=-5000N；6点承受X方向集中载荷Fx=-5000N。施加载荷，图形如图6-17 所示。,图6-16

17、编号显示设置对话框,36,图6-17圆柱实体示意图,（7）求解 选择Solution SolveCurrent LS，开始计算，计算结束会弹出计算完毕对话框，单击Close。关闭对话框计算完毕。,（8）后处理 运行 General PostprocPlot ResultsContour Plot Nodal Solu。 弹出如图6-18所示对话框，运行DOF SolutionDisplacement vector sum和Stressvon Mises stress，分别显示圆柱体的位移和应力云图。,37,图6-18 云图显示对话框,结果显示如图6-19和图6-20所示。,38,图6-19 位

18、移云图 图6-20 应力云图,39,6.6空间轴对称问题的有限元法,对空间轴对称问题，常采用圆柱坐标系。r表示径向坐标，z表示轴向坐标，任一对称面为rz面。在有限元分析时，可采用轴对称的环形单元进行。环形单元 可以是任何平面单元。,某一平面图形绕平面上某一轴旋转形成的回转体称为轴对称物体，此平面称为子午面。在动力机械，特别是叶轮机械中，有很多零件都具有轴对称特性，比如轮盘、旋转轴、承力环等。 对于直齿圆柱齿轮，由于齿的存在，严格地说它并非轴对称物体。如果忽略齿的部分(将齿用外载荷表示)，则所得到的齿根以内的旋转体部分为轴对称物体。,轴对称物体的变形及应力分布不一定是轴对称的，只有当其约束和载荷

19、都对称于旋转轴时，轴对称物体的变形和应力分布才是轴对称的。 轴对称物体+轴对称约束+轴对称载荷=轴对称系统 对轴对称系统的应力分析=轴对称物体,40,1）几何形状关于轴线对称； 2）作用于其上的载荷关于轴线对称。 3）约束条件关于轴线对称。因过z轴的任一子午面都是对称面，其 上任一点p只在该平面上发生位移，即弹性体内任一点的位移、应力与应变只与坐标r、z有关，与 无关。从而，轴对称问题可转化为二维问题，但因与平面问题有区别，常称为二维半问题。,柱坐标系,41,注意：应变 虽然与 无关，但是周向应变 ，周向应力 ，由径向位移 引起，因为径向位移会导致周长的改变。,1、基本方程,位移分量,应力分量

20、,应变分量,42,虚功方程,应变分量,轴对称问题的弹性矩阵：,43,2、轴对称问题的离散化,对于轴对称问题，利用其轴对称特性，在对其进行网格划分时可知取任意通过Z轴的截面进行，类似平面问题的网格形式。,本节以三角形单元为例。,1、位移模式,轴对称问题的环向位移恒等于零，径向r位移与轴向z位移不等于零。对于图示情形，依照平面问题的三角形单元分析，取位移模式为,代入结点位移后，可解出a1-a6，再代入上式，得,xr, yz,44,其中形函数：,；,单元中位移,根据弹性力学理论，空间轴对称问题的几何方程为,2、单元中应变,45,将u,w表达式代入上式，整理后,46,式中,其中,B矩阵中含有变量r，z

21、，因此它不是常数矩阵，即轴对称问题的三角形环形单元不是常应变单元。,47,3、单元中应力,根据弹性力学理论，空间轴对称问题的应力-应变关系为,弹性矩阵:,48,单元中任意一点的应力：,4、单元刚度矩阵,由于被积函数与无关，故在三角形截面的环单元的积分可简化为在三角形截面上的积分。故有：,49,单元刚度矩阵的积分参照图示分区，按下式采用数值积分的方法进行,50,当单元较小时，可把各个单元中的r，z 近似看作常数，并且分别等于各单元形心的坐标，即,这样，就可把各个单元近似地当做常应变单元,51,单元刚度矩阵k的分块形式,其中的近似子矩阵为,52,5、等效结点荷载,类似平面问题。对于作用于三角形环单

22、元上的体积力、表面力的等效结点力为：,体力,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

23、.7 ANSYS轴对称旋转单元计算示例,57,（1）选择单元类型 运行PreprocessorElement TypeAdd/Edit/Delete，弹出Element Types对话框单击Add，弹出Library of Element Types对话框，如图7-6所示，选择SHELL51单元。,2 ANSYS求解操作过程,图7-6 单元类型库对话框,图7-7 选择材料属性对话框,58,（2）设置材料属性 运行PreprocessorMaterial PropsMaterial Models，弹出Define Material Model Behavior对话框，如图7-7所示。双击Isot

24、ropic选项，弹出Linear Isotropic Properties for Material Number1对话框，如图7-8所示。,图7-8 设置材料属性对话框,59,（3）定义单元实常数 选择Main MenuPreprocessor Real Constants Add/Edit/Delete，弹出如图7-9所示对话框，单击Add按钮弹出Element Type for Real Constants对话框，如图7-10所示，选择Type 1 SHELL51，单击OK，弹出Real Constant Set Number 1,for SHELL51对话框，如图7-11所示，在TK(

25、I)项输入0.005，单击OK。,图7-9 实常数对话框图 7-10选择要设置实常数的 单元类型,60,图7-11设置SHELL51实常数对话框,（4）建立模型 首先生成关键点，运行主菜单Preprocessor ModelingCreateKeypointsIn Active CS，弹出如图 7-12所示对话框。 创建关键点1（0.2,0,0），2（0.2,0.6,0）。生成圆筒母线：运行Main MenuPreprocessorModeling CreateLinesLinesStraight Line，弹出拾取关键点对 话框，拾取关键点1、2，单击OK。,61,图7-12创建关键点对话框

26、,（5）设置单元属性 运行Main MenuPreprocessorMeshingMesh Tool，弹出Mesh Tool对话框，在Element Attributes下拉列表中选择Lines，然后单击其后的Set按钮弹出拾取线对话框，单击Pick All，弹出分配线单元属性对话框，将 MAT,TEAL,TYPE依次设置为1，1，1，单击OK。,（6）划分网格,62,单击Mesh Tool中Lines后的Set按钮，弹出拾取线对话框，单击Pick All弹出控制线单元尺寸对话框，将NDIV设置为10，单击OK。在Mesh Tool对话框中的Mesh下拉列表中选择Lines单击Mesh，弹出拾取线对话框，单击Pick All，划分网格完毕。 运行Plot CtrlsSty