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

1、第六章第六章 三维问题有限元分析三维问题有限元分析讲讲 授：陈得良授：陈得良 TelTel：1332731100813327311008QQQQ：416501065416501065EmailEmail：deliang_deliang_1四四 教学基本内容教学基本内容第五章第五章 三维问题有限元分析三维问题有限元分析 第一节第一节 三维应力状态三维应力状态 第二节第二节 4节点四面体单元节点四面体单元 第三节第三节 8节点六面体等参单元节点六面体等参单元 第四节第四节 20节点等参单元节点等参单元 第五节第五节 ansys空间问题实例空间问题实例 第六节第六节 空间轴对称问题有限元法空间轴对称

2、问题有限元法 第七节第七节 Ansys轴对称旋转问题实例轴对称旋转问题实例2工程实际中的很多问题工程实际中的很多问题难于简化难于简化为平面问题，如受任为平面问题，如受任意空间载荷作用的任意形状几何体，受对称于轴线载意空间载荷作用的任意形状几何体，受对称于轴线载荷作用的回转体荷作用的回转体,这类问题经典弹性力学往往无能为力。这类问题经典弹性力学往往无能为力。在在FEM中，空间问题只要求中，空间问题只要求0阶连续，因此构造单元阶连续，因此构造单元方便方便空间问题简介空间问题简介3空间问题的主要困难：空间问题的主要困难：（1）离散化不直观；）离散化不直观；（网格自动生成）（网格自动生成）（2）分割的

3、单元数量多，未知量的数目剧增。）分割的单元数量多，未知量的数目剧增。 （对某些问题简化）（对某些问题简化） （轴对称问题）（轴对称问题）空间分析的优点空间分析的优点 精确精确46.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由弹性力学知，应变与位移间的几何关系是 xuxxvyuxyyvyywzvyzzwzzuxwzx

4、 (6-1)三维弹性体的应变分量，用矩阵表示为 wvuxzyzxyzyxzxyzxyzyx000000000(6-2)6弹性体受力作用，内部任意一点的应力状态也是三维的，用列向量表示为 Tzxyzxyzyx在线弹性范围内，应力与应变间的物理关系矩阵表达式为 D 对于各向同性弹性体，在三维应力状态下，弹性矩阵 的形式为(6-3) D 1221000001221000012210001111112111称对ED(6-4)7(1)空间问题常用单元：四面体单元、长方体单元、直边六面体单元、曲边六面体单元 、轴对称单元。o4结点四面体单元：结点四面体单元：是空间问题最简单的单元，也是常应变、常应力单元，

5、可以类似平面问题三结点三角形单元进行分析。o8结点长方体单元结点长方体单元：可以类似平面四结点矩形单元进行分析。o8结点直边六面体单元结点直边六面体单元：可以类似平面四结点任意四边形等参元分析 。o20结点曲边六面体单元：结点曲边六面体单元：等参单元,可以类似平面八结点曲边四边形等参元进行分析 。o轴对称单元轴对称单元：一平面单元绕一对称轴旋转形成的空间问题。只需在rz平面划分网格，就像平面问题xy平面中的网格一样，这样这类空间问题可以得到简化。 （环向位移等于零）（2）结点位移3个分量。（3）基本方程比平面问题多。3个平衡方程，6个几何方程，6个物理方程。86.2 四节点四面体单元四节点四面

6、体单元图6-2表示任一简单四面体单元，其中四个结点编号设为 i、j、m、n (或1、2、3、4)。单元变形时，各结点沿x、y、z方向上的位移，以列向量表示为图6-2四面体单元8 Tnnnmmmjjjiiiewvuwvuwvuwvu 单元变形时，单元内各点也有沿x、y、z方向的位移u、v、w，一般应为坐标x、y、z的函数。对于这种简单的四面体单元，其内部位移可假设为坐标的线性函数，为满足变形协调条件，取为zyxwzyxvzyxu121110987654321(6-5)式(6-5)含有12个待定系数a，可由单元的12项结点位移决定.将4个结点的坐标值代入式(6-5)的u式中。 i、j、m、n共4个

7、结点，分别有nnnnmmmmjjjjiiiizyxuzyxuzyxuzyxu4321432143214321(6-6)1 1 单元形函数单元形函数9nnmmjjiiuNuNuNuNu 其中 zdycxbaVNiiiii61式中，V为四面体的体积，且有nnnmmmjjjiiizyxzyxzyxzyxV111161nnmmjjinnmmjjinnmmjjinnnmmmjjjiyxyxyxdzxzxzxczyzyzybzyxzyxzyxa111111111(6-7)由式(6-6)求出 ，再代回式(6-5) 中，整理后得1234,a a a a10 为使四面体的体积V不为负值，在右手坐标系中，使右手旋

8、转按着由i- j- m的转向转动时，且法向n方向前进。用求位移u的同样方法，可求得nmjiiinnmmjjiivNvNvNvNvNv,nmjiiinnmmjjiiwNwNwNwNwNw,将位移的3个线性方程形成的线性方程组用矩阵表示为 eNwvu(6-8) ININININNnmji式中(6-9)112 单元刚度矩阵单元刚度矩阵 将式(6-8)代入几何方程式(6-2)，经过微分运算，可得单元内应变为 enmjieBBBBB(6-10) 式中 iiiiiiiiiiiiiiiiiiibdcdbcdcbVxNzNyNzNxNyNzNyNxNB00000000061000000000(6-11) 简单

9、四面体单元内，各点的应变都是一样的，这是一种常应变单元（是三维单元中精度最低的单元）。这一点与平面问题的简单三角形单元相似，由于单元内位移都假定为线性变化的，因而由位移一阶导数组成的应变也为常量。 12 同样，用虚功原理建立结点力和结点位移间的关系式，从而得出简单四面体单元的刚度矩阵。 eVedvBDBdxdydzBDBkTT(6-12) eeVBDBkT(6-13)按结点分块表示，此单元刚度矩阵可表示为 nnnmnjnimnmmmjmijnjmjjjiinimijiiekkkkkkkkkkkkkkkkk(6-14)13(r=i、j、m、n， S=i、j、m、n ) (6-15)式中11A12

10、212A 弹性体三维（空间）问题的原始平衡方程组，即 KF eneekK1其中其中任一子矩阵为 esrrsVBDBkTVE211361srsrsrsrsrsrsrdbAbdAcbAbcAddccAbb21212srsrsrsrsrsrsrsrsrsrsrsrsrsrccbbAdddcAcdAcdAdcAddbbAccbdAdbAbcAcbA2212122121143 整体结构载荷列向量整体结构载荷列向量整体结构的结点载荷列向量 11eenneeeeepsCeeFFFFFF(6-16)式中 单元上集中力等效结点载荷列向量；单元上表面力等效结点载荷列向量；单元上体积力等效结点载荷列向量；单元结点载

11、荷列向量。 epF eSF eF eCF等效结点力公式为 TepFNF TeeSSSFNpds TeevVFNpdV TxyzFFFF式中 156.3 8节点六面体等参单元节点六面体等参单元8 (x5,y5,z5)1234 (x4,y4,z4)5 (x5,y5,z5)67xzy 23(1,-1,1)48(1,1,-1)657如同二维等参单元一样，三维等参单元的有关公式的建立也是采用局部自然如同二维等参单元一样，三维等参单元的有关公式的建立也是采用局部自然坐标（曲面坐标）坐标（曲面坐标）。可参考的母体单元则为一正六面体，图可参考的母体单元则为一正六面体，图6-36-3则表示了任则表示了任意六面体

12、单元与母体单元、局部的三维自然坐标与整体的直角坐标系的几何意六面体单元与母体单元、局部的三维自然坐标与整体的直角坐标系的几何关系关系母体单元母体单元任意六面体单元任意六面体单元（对面不全平行）（对面不全平行）图6-3 8结点三维等参单元16 用形函数表示的位移插值形式的位移模式，可以直接利用拉格朗日插值公式，得到单元位移函数为888111, , ,iiiiiiiiiuN uuN uuN u 根据等参单元的定义。自然坐标与整体直角坐标之间的关系可以写为81000000iiiiiiizyxNNNzyx 其中形函数为1(1)(1)(1)8iiiiN )1)(1)(1(818N例 其中 为节点坐标值（

13、角点） (1 8)iiii， ，17 由节点位移求单元应变时，他要求形函数在整体坐标下的导数，但形函数是建立在局部坐标下的，这就需要将局部坐标中的表达式转换到整体坐标系中，如同平面等参单元一样，需要通过雅克比矩阵来实现，由偏导法则iiiiNNNNxyzxyz同理可得,iiNN写成矩阵iiiiiiiiiNxyzNNxxNNNxyzyyNxyzNNzzJ181iiiiiiNNxNNyNNzJ求单元刚度矩阵，尚需对积分的单元体积进行积分变换反之有了上式很容易得到单元的应变应力矩阵 eBeDB111111eTTkB DBdB DB J d d d 19 6.4 20 6.4 20结点等参元结点等参元

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

15、为 iiiiizyxNzyx201(6-17) 则单元的位移函数可写成21iiiiiwvuNwvu201(6-18)在自然坐标系（局部坐标系）中，各结点的形状函数可写成如下形式， 对于8个顶角结点（ i1，2，8）式中 xi、yi、zi结点i的坐标； ui、vi、wi结点i沿x、y、z方向的位移； Ni对应于i结点的形状函数。)2)(1)(1)(1 (81iiiiiiiN22)1)(1)(1 (412iiiN0i对于 的边上点（i17，18，19，20）)1)(1)(1 (412iiiN(6-19)对于 的边上点（i9，11，13，15）0i)1)(1)(1 (412iiiN0i对于 的边上点

16、（i10，12，14，16）232 单元刚度矩阵单元刚度矩阵三维变形状态下，一点的应变与位移的几何关系为 ezxyzxyzyxBwvuxzyzxyzyx000000000(6-20)24 eeVVedVBBBDBBBdVBDBk2021T2021T(6-22) 为便于以下计算，弹性矩阵D可分块写为 2100DDD(6-23) 令 ,)21)(1 (E12EG 则 ，GGGD2221GGGD0000002 202020220122022211201211kkkkkkkkkke为6060的方阵，可按结点写为子块形式 ek25式中第i行j列的子矩阵为 eVjieijdVBDBkT33(6-24)将将

17、(6-20)、(6-22)分块式代入分块式代入(6-23)，其被积函数可写为，其被积函数可写为 jjiijiSTDDSTBDB21TT00zzzyzxyzyyyxxzxyxxHHHHHHHHH)()2(zNzNyNyNGxNxNGHjijijixxxNyNGyNxNHjijixyyNxNGxNyNHjijiyz(6-25)式中与式(5-9)相似，按坐标变换式(6-17)，应有26 iiiiiiNNNJzNyNxN1同样可有(6-26)dVJ ddd三维六面体的雅可比矩阵为 iiiiiiiiiiiiiiiiiizNyNxNzNyNxNzNyNxNzyxzyxzyxJ(6-27)同理可采用三维高斯

18、求积公式计算单元刚度矩阵。即 27 eVjTieijdVBDBkdddJHHHHHHHHHzzxzzxyzyyyxxzxyxx 111111kjiLiMjNkzzxzzxyzyyyxxzxyxxkjiHHHHHHHHHwww 式中，L，M，N为沿 、 、 方向的积分点数目，而积分点坐标 及权重 可由高斯积分表查得。kji、kjiwww、 对于20节点的三维单元，通常可取积分点数目m=3(即3 3 3），查表可得对应的积分点坐标和权重为28d1121310.774 596 669 241 384, 0.555 555 555 555 5560.000 000 000 000 000, 0.888

19、 888 888 888 8880.774 596 669 241 384, 0.555 55www5 555 555 556 参照方法，可以很方便的得到另外两个方向的积分点坐标值和权重。296.5 ANSYS空间问题示例空间问题示例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

20、 2 ANSYS求解操作过程求解操作过程30 （1）选择单元类型选择单元类型 运行PreprocessorElement TypeAdd/Edit/Delete，弹出Element Types对话框，如图6-5所示。然后单击Add，弹出Library of Element Types窗口，如图6-6所示，选择SOLID45单元，单击OK。图6-5 单元类型对话框 图6-6 单元类型库对话框31 （2）设置材料属性设置材料属性 运行PreprocessorMaterial PropsMaterial Models，弹出如图6-7所示对话框。双击Isotropic，弹出Linear Isotrop

21、ic Properties for Material Number1对话框，如图6-8所示，在EX选项栏中设置数值2.1e11，在PRXY选项栏中设置数值0.3。设置完毕单击OK按钮。图6-7 选择材料属性对话框 图6-8 设置材料属性对话框 32 （3）建立模型建立模型 运行PreprocessorModelingCreateAreasRectangleBy 2 Corners，弹出如图6-9所示对话框，在WP X选项栏中填写0，在WP Y选项栏中填写0，在Width选项栏中填写0.05，在Height选项栏中填写0.2，点击OK。生成如图6-10所示图形。图6-9 两点建立矩形对话框 图6

22、-10 生成的长方形面 33 将长方形旋转成柱体，运行PreprocessorModeling OperateExtrudeAreasAbout Axis，弹出如图6-11所示拾取框。选择图7中长方形后弹出单击OK，再选择长方形左上角和左下角结点后，单击OK.。弹出如图6-12所示对话框。在ARC选项栏中填入旋转角度360度，设置完毕单击 O K 按 钮 ， 生 成 如 图 6 - 1 3 所 示圆柱体。 图6-11 拾取对称轴对话框 图6-12 设置绕轴旋转参数对话框图6-13 圆柱模型 34 运行MeshingSize CntrlsManualSizeGlobalSize弹出如图6-14所

23、示对话框，设置SIZE选项栏中的数据为0.01。运行MeshingMeshVolumesFree自由划分网格后得到如图6-15所示图形。图6-14 设置网格尺寸对话框 图6-15 圆柱有限元模型 （5）施加约束施加约束 运行SolutionDefine LoadsApplyDisplacementOn Areas，拾取圆柱的底面，施加全约束。 35 （6）施加载荷施加载荷 显示图形的关键点，运行PlotCtrlsNumbering弹出如图6-16所示对话框，激活KP Numbers后面的选框，使它变成on形式。选择菜单SolutionDefine LoadsApplyStructure For

24、ce/Moment On Keypoints，载荷分别如下：8点承受Z方向集中载荷Fz=5000N和Y方向集中载荷Fy=-5000N；10点承受X方向集中载荷Fx=5000N；3点承受Z方向集中载荷Fz=-5000N；6点承受X方向集中载荷Fx=-5000N。施加载 荷 ， 图 形 如 图 6 - 1 7所示。图6-16 编号显示设置对话框36图6-17圆柱实体示意图 （7）求解求解 选 择 S o l u t i o n S o l v e C u r r e n t LS，开始计算，计算结束会弹出计算完毕对话框，单击Close。关闭对话框计算完毕。 （8）后处理后处理 运行 General

25、 PostprocPlot ResultsContour PlotNodal Solu。 弹出如图6-18所示对话框，运行DOF SolutionDisplacement vector sum和Stressvon Misesstress，分别显示圆柱体的位移和应力云图。 37图6-18 云图显示对话框 结果显示如图6-19和图6-20所示。 38 图6-19 位移云图 图6-20 应力云图396.6空间轴对称问题的有限元法 对空间轴对称问题，常采用圆柱坐标系。r表示径向坐标，z表示轴向坐标，任一对称面为rz面。在有限元分析时，可采用轴对称的环形单元进行。环形单元 可以是任何平面单元。某一平面图

26、形绕平面上某一轴旋转形成的回转体称为轴对称物体，此某一平面图形绕平面上某一轴旋转形成的回转体称为轴对称物体，此平面称为子午面。在动力机械，特别是叶轮机械中，有很多零件都具平面称为子午面。在动力机械，特别是叶轮机械中，有很多零件都具有轴对称特性，比如轮盘、旋转轴、承力环等。有轴对称特性，比如轮盘、旋转轴、承力环等。对于直齿圆柱齿轮，由于齿的存在，严格地说它并非轴对称物体。如对于直齿圆柱齿轮，由于齿的存在，严格地说它并非轴对称物体。如果忽略齿的部分果忽略齿的部分( (将齿用外载荷表示将齿用外载荷表示) )，则所得到的齿根以内的旋转体，则所得到的齿根以内的旋转体部分为轴对称物体。部分为轴对称物体。轴

27、对称物体的变形及应力分布不一定是轴对称的，只有当其约束和载轴对称物体的变形及应力分布不一定是轴对称的，只有当其约束和载荷都对称于旋转轴时，轴对称物体的变形和应力分布才是轴对称的。荷都对称于旋转轴时，轴对称物体的变形和应力分布才是轴对称的。轴对称物体轴对称物体+轴对称约束轴对称约束+轴对称载荷轴对称载荷=轴对称系统轴对称系统对轴对称系统的应力分析对轴对称系统的应力分析=轴对称物体轴对称物体401）几何形状关于轴线对称；2）作用于其上的载荷关于轴线对称。3）约束条件关于轴线对称。因过z轴的任一子午面都是对称面，其上任一点p只在该平面上发生位移，即弹性体内任一点的位移、应力与应变只 与坐标r、z有关

28、，与 无关。从而，轴对称问题可转化为二维问题，但因与平面问题有区别，常称为二维半问题。zrxp( , , )rz柱坐标系柱坐标系41注意注意：应变应变 虽然与虽然与 无关，但是周向应变无关，但是周向应变 ，周向应力，周向应力 ，由径向位移，由径向位移 引起，因为径引起，因为径向位移会导致周长的改变。向位移会导致周长的改变。 =0 Tru wu Trzrz 1、基本方程、基本方程位移分量位移分量应力分量应力分量 Trzrz 应变分量应变分量00ru42虚功方程虚功方程 = TrzrzTrrruuuwwrrzzr 2*02 2 TTdFrdrdz则 应变分量应变分量轴对称问题的弹性矩阵：轴对称问题

29、的弹性矩阵： 1 21012(1)100(1) 1(1)(1 2 )11 22(1)ED对称432、轴对称问题的离散化、轴对称问题的离散化 对于轴对称问题，利用其轴对称特性，在对其进行网格划分时可知取任意通过Z轴的截面进行，类似平面问题的网格形式。本节以三角形单元为例。1、位移模式o轴对称问题的环向位移恒等环向位移恒等于零于零，径向r位移与轴向z位移不等于零。对于图示情形，依照平面问题的三角形单元分析，取位移模式为zrwzru654321代入结点位移后，可解出a1-a6，再代入上式，得 mmjjiimmjjiiwNwNwNwuNuNuNuxr,yz44o 其中形函数：),)(21mjizcrb

30、aANiiiiemjieNNNNfjmimjijmmjirrczzbzrzra ； ),(mji单元中位移根据弹性力学理论，空间轴对称问题的几何方程为2、单元中应变45rwzuzwrururzzr将u,w表达式代入上式，整理后emjieBBBB46式中),(00021mjibccfbABiiiiii),(mjirzcbrafiiii其中 B矩阵中含有变量r，z，因此它不是常数矩阵，即轴对称问题的三角形环形单元不是常应变单元。473、单元中应力根据弹性力学理论，空间轴对称问题的应力-应变关系为 Drzzr121000111111)21)(1 ()1 (称对ED弹性矩阵:48单元中任意一点的应力：

31、eeSBD4、单元刚度矩阵VTVTdzrdrdBDBdVBDBk 由于被积函数与无关，故在三角形截面的环单元的积分可简化为在三角形截面上的积分。故有： ATrdrdzBDBk249)()(312112),(),(),(),(rzrzbczrzabAAdzzrgdrdzzrgdrdrdzzrgdrdzzrgG单元刚度矩阵的积分参照图示分区，按下式采用数值积分的方法进行50 当单元较小时，可把各个单元中的r，z 近似看作常数，并且分别等于各单元形心的坐标，即)(31)(31mjimjizzzzzrrrrrrzcbraffiiiii这样，就可把各个单元近似地当做常应变单元 2BDBArkT:),(成

32、为mjirzcbrafiiii获得。代替中用在zrzrB,51 单元刚度矩阵k的分块形式mmmjmijmjjjiimijiikkkkkkkkkksrsrsrssrsrrrssrsrrsrsrbbAcccbAfbcAbcAfbcAccAbfbAffbbAAr22121213)()()(2ArBDBksrrs2其中的近似子矩阵为)21)(1 ()1 ()1 (2211321EAAA525、等效结点荷载类似平面问题。对于作用于三角形环单元上的体积力、表面力的等效结点力为：rdrdzpNPAVTeV2体力AieVziVrieVirdrdzNPPP02),()3(60mjirrAi)(31mjimmjj

33、iirrrrLrLrLrr53o面力：面力：lSTeSrdlpNP2（1）均布表面力 设单元ij边上作用均布表面力，其集度为zrSppP zrjieSipprrlP)2(3zrjieSjpprrlP)2(3zrSppP l当 ri=rj 时，静力等效原则54(2)三角形分布表面力 沿单元ij边作用了三角形分布的表面力，表面力在i点集度为 zrsppp zrjieSipprrlP)3(6zrjieSjpprrlP)(6当 ri=rj 时，静力等效原则。2/3集中在i点，1/3集中在j点。5556 圆筒直径0.4m，高度0.6m，壁厚0.005m；材料Q235，弹性模量E=2.1e11Pa，泊松比

34、=0.3；约束：圆筒的下部在轴线方向固定，其它方向自由；载荷：顶部环线上承受轴向线压力P-200000N/m。 图6-4 圆筒示意图 图6-5单元类型对话框 1 问题描述问题描述 6.7 ANSYS轴对称旋转单元计算示例轴对称旋转单元计算示例 57（1）选择单元类型）选择单元类型 运行PreprocessorElement TypeAdd/Edit/Delete，弹出Element Types对话框单击Add，弹出Library of Element Types对话框，如图7-6所示，选择SHELL51单元。2 ANSYS求解操作过程求解操作过程 图7-6 单元类型库对话框 图7-7 选择材料

35、属性对话框 58 （2）设置材料属性 运行PreprocessorMaterial PropsMaterial Models，弹出Define Material Model Behavior对话框，如图7-7所示。双击Isotropic选项，弹出Linear Isotropic Properties for Material Number1对话框，如图7-8所示。 图7-8 设置材料属性对话框 59 （3）定义单元实常数 选择Main MenuPreprocessor Real Constants Add/Edit/Delete，弹出如图7-9所示对话框，单击Add按钮弹出Element Ty

36、pe for Real Constants对话框，如图7-10所示，选择Type 1 SHELL51，单击OK，弹出Real Constant Set Number 1,for SHELL51对话框，如图7-11所示，在TK(I)项输入0.005，单击OK。 图7-9 实常数对话框图 7-10选择要设置实常数的 单元类型 60图7-11设置SHELL51实常数对话框 （4）建立模型 首先生成关键点，运行主菜单PreprocessorModelingCreateKeypointsIn Active CS，弹出如图7-12所示对话框。 创建关键点1（0.2,0,0），2（0.2,0.6,0）。生成圆筒母线：运行Main MenuPreprocessorModelingCreateLinesLinesStraight Line，弹出拾取关键点对话框，拾取关键点1、2，单击OK。 61 图7-12创建关键点对话框 （5）设置单元