有限元板壳单元

有限元板壳单元有限元板壳单元 板壳结构在工程上应用十分广泛。在设板壳结构在工程上应用十分广泛。在设计分析中采用板壳单元进行结构分析，可以计分析中采用板壳单元进行结构分析，可以得到足够的精度和良好的效果。得到足够的精度和良好的效果。第八章第八章 关于板壳单元关于板壳单元 第八章第八章 关于板壳单元关于板壳单元 8.1 板壳结构板壳结构 板壳结构板壳结构是指板的厚度t与其它两个方向的尺寸相比小得多。板壳结构的板可以是平板也可以是单曲面或双曲面板，同时可以承受任意方向上的载荷，也就是既有作用在平面内的载荷，又作用有垂直于平面的载荷。一般板壳结构处于三维应力状态。结构是否为板壳问题，需要确定厚度与其它方位尺寸的比值，如果 1/80t1/10可以归结为板（薄壳）问题，若介于1/10 1/5 之间属于厚壳问题，若大于 1/5 则不属于板壳结构问题。板壳单元的力学模型取为结构单元的中性面，即以各中性面来代表为不同厚度的板或壳单元的组合体，以此来模拟结构体。在工程有限单元法的软件设计中，常常将板壳结构划分成薄板、厚板以及壳单元。如图8-1所示平板，取其中性面为坐标面，z轴垂直于中性面。其中 t 为板厚。当板受有垂直于板中性面的外力时，板的中性面将发生弯扭变形，从而变成一个曲面。板变形的同时，在板的横截面上将存在内力弯矩和扭矩。第八章第八章 关于板壳单元关于板壳单元 图8-1 平板弯曲 对于薄板弯曲问题采用如下假设：a.板的法线没有伸缩；b.板的法线在板变形后仍垂直于中性面；8.2 薄板基础理论知识薄板基础理论知识 图8-2所示为板的一个微元体。为方便计，取 和 的方向的宽度均为1。在垂直于 轴的横截面上的正应力与z 坐标成正比，并可合成为一个力偶，从而构成该横截面上的弯矩（单位宽度上的弯矩）。同理，弯矩，和合成扭矩和。由于剪应力互=。内力列向量为合成等，因此图8-2 薄板微元体内力与应力示意图 c.板内各点没有平行于中性面的位移；d.垂直于板面挤压应力可以不计。=（8-11）是对中性面力矩的合成（见图8-2），即引用记号（8-12）则（8-13）式中 弹性薄板的应力应变转换矩阵，它等于平面应力问题中的 与的乘积。根据与 之间的关系，不难由（8-13）和（8-10）式求出板上下表面 的应力（8-15）（8-14）综上所述，薄板的中性面挠度w 是基本的未知量。由w即可计算出位移、应变、应力及内力。8.3 3结点三角形薄板单元结点三角形薄板单元8.3.1坐标变换坐标变换 图8-3为一个任意形状的3结点三角形板单元，结点编号 1、2、3 按右手法则排序。图8-3（a）为单元直角坐标系 ,图8-3（b）为单元自然坐标系 。（a）单元直角坐标系 （b）单元自然坐标系图8-3 3结点三角形板单元坐标系 单元坐标变换（8-16）为面积坐标。式中（8-17）面积坐标 具有插值函数的性质，即（8-18）8.3.2 位移向量位移向量根据薄板理论，薄板结点位移如图8-4所示。图8-4 薄板结点位移示意图单元任一结点位移列向量为（8-19）单元结点位移列向量（8-20）单元内任意点的位移w用结点位移插值表示如下（8-21）其中，、和 为插值函数，是 的行阵（8-22）插值函数具体形式如下（8-23）其中，。8.3.3 应变位移转换矩阵应变位移转换矩阵 为了建立单元刚度矩阵，需要建立位移应变转换矩阵，即建立 与单元结点位移 的关系式。将式（8-21）代入式（8-5），可得（8-24）式中（8-25）矩阵是插值函数 的二阶导数。是 的函数，它们对x和y的偏导数按复合函数求导法则（8-26）类似地有（8-27）对式（8-26）和式（8-2）二阶求导（8-28）式中 为二阶微分算子。（8-29）由式（8-23）可得（8-30）同样，是与单元结点坐标有关的数，见第3章。8.3.4 单元刚度矩阵单元刚度矩阵由虚功原理得到薄板的单元刚度矩阵（8-31）一般采用哈默值积分来计算式（8-31）比较方便。（8-32）8.4 厚板基础理论知识厚板基础理论知识 厚板理论假设如下：a板的挠度w微小；b板中性面法线在变形后仍保持直线,但不再垂直变形后的中曲面；c垂直于中性面的应力可以忽略。由此确定了板的独立位移分量为在薄板理论中，因不考虑横向剪切变形，即因此，与薄板理论类似，板的曲率和扭率为（8-33）式（8-33）与薄板的区别在于，这里还要考虑由于剪切而产生的应变（8-34）图8-6 厚板微元体内力与应力示意图厚板的应力应变关系如下（8-35）对于各向同性材料有（8-35）8.5 4结点四边形板单元结点四边形板单元8.5.1 坐标变换坐标变换 图8-7 所示为任意四边形板单元，结点编号按逆时针排序。图8-7（a）为直角坐标系，图8-7（b）为自然坐标系。实际单元与基本单元的对应关系为（8-39）(a)直角坐标系与实际单元 (b)自然坐标系与基本单元图8-7 四结点四边形厚板单元 8.5.2 单元位移场与应变位移转换矩阵单元位移场与应变位移转换矩阵单元任意结点位移有三个独立分量 单元结点位移列向量（8-40）同样，用表示单元的几何坐标变换的结点插值函数来表示单元内任意点的位移变换关系。（8-41）将位移表达式（8-40）代入应变式（8-33）和式（8-34）可得（8-42）式中（8-43）（8-44）8.5.3 单元刚度矩阵单元刚度矩阵单元刚度矩阵由弯曲和剪切两部分构成（8-45）式中 由于弯曲引起的刚度；由于剪切引起的刚度。8.6 壳单元壳单元8.6.1 板壳结构物理特性板壳结构物理特性 壳体的中性面是一个曲面，壳单元受力状态及应力状态见图8-8。(a)壳单元受力状态示意图 (b)壳单元应力状态示意图图8-8 壳单元 在作结构分析时，一般采用平面单元（板）或者曲面单元处理。平面单元是平面应力单元和平面弯曲单元的组合体，它依赖于平板理论。在几何上以平板代替壳体，结构模拟是一种近似。但是，这种单元简单，只要结构离散化分合理，完全可以满足工程上的要求。曲面单元能够更好地模拟真实结构，相应得到的计算结果会更有效。但是，曲面壳体的变形与平板变形有所区别。壳体的中性面变形不能忽略，在壳体中的内力包括弯曲内力和中性面内力。对于曲面单元，现常采用考虑横向剪切变形的超参数曲面壳单元。有关内容本节不作介绍，需要时读者自己可以查阅相关资料。平面壳单元的物理特性如下：平面壳单元可以视为平面应力单元与板弯曲单元的平面应力单元与板弯曲单元的组合体。组合体。平面应力单元（亦称膜单元）仅仅能够承受作用于平面内的载荷，不能够承受其它载荷。假设z方向上的位移w=0，每一结点仅存在沿x轴和y轴的位移（a）平面单元模拟 （b）平面弯曲单元模拟图8-9 壳单元模拟。板弯曲单元仅仅承受弯曲载荷，此类单元只有沿坐标z方向的位移，见图8-10。图中阴影部分为中性 面，、分别为板壳结构的顶面、中性面及底面应 力。第八章第八章 关于板壳单元关于板壳单元 图8-10 四结点板弯曲单元示意图平面壳单元结点的位移列向量为结点力列向量为 单元内应力分量可以简单地将相应的平面应力单元和板弯曲单元的应力分量进行叠加。（8-46）式中、平面应力单元（膜单元）应力分量；、壳单元应力分量；板弯曲单元应力分量。8.6.2 三结点平面壳单元三结点平面壳单元 三结点平面壳单元的整体坐标与局部坐标一般如图8-11（a）所示。三结点三角形壳单元结点位移及内力列向量如图8-11（b）所示。(a)整体坐标系与局部坐标系 (b)位移列向量示意图图8-11 三结点三角形平面壳单元在局部坐标系下在整体坐标系下 三角形平面壳单元的3个结点i，j，m共有15个自由度，位移函数可采用如下形式。（8-47）（8-48）三结点三角形平面壳单元的位移函数可以看成是平面应力单元位移函数式(8-47)与板弯曲单元的位移函数式(8-48)的组合。按照有限单元法求解步骤，不难得出在局部坐标下的单元刚度矩阵。经坐标变换后得出整体坐标下的单元刚度矩阵8.7 ANSYS板壳单元计算示例板壳单元计算示例 开口的圆弧为 180，半 径 为 0.05m，CA长 为0.1m，壳 体 厚 0.001m，=90，弹 性 模 量E=210GPa，=0.3。约束：边 CD和 边 DE全约束。承受载荷：B点作用集中载荷 F10N，方向水平向右。8.7.1 问题描述问题描述 根据平衡条件建立原始平衡方程同样，原始平衡方程经过约束出理后即可求解。8.7.2 ANSYS求解操作过程求解操作过程 （1）选择单元类型 运 行 PreprocessorElement TypeAdd/Edit/Delete，弹出Element Types对话框如图 8-6所示，然后单击 Add，弹出Library of Element Types窗口，如图 8-7所示，选择Shell Elastic 4node 63，选择完毕单击OK按钮。图8-5 壳体示意图 图8-6 单元类型对话框 第八章第八章 关于板壳单元关于板壳单元 （2）设置材料属性 运 行 PreprocessorMaterial PropsMaterial Models，弹出图 8-8所 示对话框，双击 Isotropic，弹出如图 8-9所示对话框，在EX选项栏中设置数值2.1e11，在 PRXY选项栏中 设 置 数 值 0.3。设置完毕单击OK按钮。图8-7 单元类型库对话框 图8-8 选择材料属性对话框 第八章第八章 关于板壳单元关于板壳单元 第八章第八章 关于板壳单元关于板壳单元 （3）设置壳厚参数 运 行 PreprocessorReal ConstantsAdd/Edit/Delete，弹出如图 8-10所 示对话框，设置 TI(k)选项栏为0.001。设置完毕单击 OK按钮完成设置。图8-9 设置材料属性对话框 图8-10 设置实常数对话框 （4）建立模型 生成圆环，选择PreprocessorModelingCreateLinesArcsBy Cent&Radius，弹出如图 8-11所 示对 话 框，填 写 数 据 0,0,0，单 击Apply按钮，接着对话框提示再选取一点，填入 0.05作为半径，单击OK按钮，弹出如图 8-12所示对话框。图8-11依圆心和半径生成弧线对话框 图8-12圆弧旋转角度设置对话框 第八章第八章 关于板壳单元关于板壳单元 在ARC选项栏中填写旋转角度180，生成圆弧。建立关键点（0.05,0,-0.1），（-0.05,0,-0.1），（0,0.05,-0.1）。建立完关键点后沿关键点生成线，将图形连接如图8-13所示图形。运行 ModelingOperateBooleansAddLines，利用布尔运算将线相加：将 CD，DE相加为一条线，另外将圆弧的两段线加为一条线。生生成成面面：运行PreprocessorModelingCreateAreasArbitraryBy lines，拾取相交后的线可得到如图 8-14所示的图形。图8-13 壳体线框模型 图8-14 壳体模型 第八章第八章 关于板壳单元关于板壳单元 （5）划分网格 设置网格大小：运行 MeshingSize Cntrls ManualSizeAreasAll Areas，弹出如图 8-15所示对话框，设置Element edge length为0.005。设置完毕点击 OK。运行MeshingMeshAreasFree，拾取壳体。划分网格完毕，如图8-16所示。图8-15 设置网格尺寸对话框 图8-16 壳体有限元模型 第八章第八章 关于板壳单元关于板壳单元 （6）施加约束 选择菜单SolutionDefine LoadsApply StructureDisplacementOn Lines，拾取线段CD和DE，选择All DOF在这两条线段上施加全约束。（7）施加载荷 执行命令SolutionDefine LoadsApply StructureForce/MomentOn Keypoints，拾取B点,施加载荷FX10N，如图8-17所示。图8-17 壳体约束与载荷 （8）求解 选择SolutionSolveCurrent LS，开始计算，计算结束会弹出计算完毕对话框，单击 C l o s e关闭对话框，计算完毕。第八章第八章 关于板壳单元关于板壳单元 （9）后处理 运 行 General PostprocPlot ResultsContour PlotNodal Solu，弹出如图 8-18所示对话框。运行 DOF SolutionDisplacement vector sum和Stressvon Mises stress，分别显示壳体位移和应力云图。结果如图 8-19和图8-20所示。图8-18 云图显示对话框 第八章第八章 关于板壳单元关于板壳单元 从图8-19壳体的位移云图可知，最大位移发生在B点处，且最大位移值为0.000707m。从图8-20壳体的应力云图可知，最大应力发生在B点处，且最大应力值为50.2MPa。第八章第八章 关于板壳单元关于板壳单元 图8-19 位移变形云图 图8-20 等效应力云图 8.7.3 结论结论 结束语结束语谢谢大家聆听！谢谢大家聆听！40