第八章一般壳体问题的有限元法课件.ppt

壳体实质上是从平板演变而来的，它的中面是一个曲面。在分析壳中应力时，虽然平板的基本假定同样有效，但是壳体的变形有着很大程度的不同，它除了弯曲变形外还存在中面变形。因而，壳中内力包括有弯曲内力和中面内力。应用有限单元法分析壳体结构时，广泛地采用了平面单元和曲面的单元这两类壳体单元。本章首先介绍平面单元，它是平面应力问题和平板弯曲问题的组合；这种单元虽然简单，但是相当有效。然后讨论一个考虑横向剪切影响的曲面单元，称为八结点40个自由度的一般壳单元，可以适用于厚壳和薄壳。将壳体曲面划分为有限个单元，它们都是曲面单元。但是在单元细分时，用平面单元组成的一个单向或双向折板来近似壳体的几何形状将会得到良好的结果。通常对于任意形状的壳体，采用三角形单元比较方便，如图8-1所示。如果在壳体上容易找到同一平面上的四个点，可以采用平面四边形单元。例如具有正交边界的柱面壳体，如图8-2所示。图8-1 任意壳体作为平面三角形单元的集合 图8-2 圆柱壳作为平面矩形单元的集合 壳体平面单元的应力状态是由平面应力和弯曲应力的叠加而成的，因此在构造壳体平面单元时，只要将第二章和第七章所讨论的相应单元进行简单的组合就可以了。同样，前述二章所导出的刚度矩阵可作为建立壳体平面单元刚度矩阵的基础。现在把平面单元的计算步骤归纳如下 1.划分单元，选定整体坐标系 oxyz，定出节点在整体坐标系中的坐标值。2.对于各个单元利用节点坐标值，建立一个局部坐标系例如三角形单元123，可以选取节点1为局部坐标系的原点，并且以1-2边为 轴的正方向，如图8-3所示。于是，方向的单位e1求得是(8-1)式中 是矢量12的长度。取单元的外法线方向作为 轴的正方向，于是它的单位矢量图8-3 三角形单元局部坐标系容易看出，矢量12和13的矢性积的模等于三角形面积的一倍，即|12121313|=2。最后，按右手定则可以决定y轴的正方向,它的单位矢量e2是e2=e3 e1 (8-3)利用上述方法确定的局部坐标系，三角形单元123是在 平面内，它的三个角点的局部坐标值是很容易确定的。对于柱面上的矩形单元，局部坐标的原点 选在矩形的形心，通常选 轴和x轴均沿柱面母线方向。如图8-4中所示，由矢量12确定单位矢量e1,再由矢量14确定单位矢量e2，于是e3=e1 e2。3.对于各个单元，确立在局部坐标系 中的结点载荷列阵 。壳体载荷可以分解成二组：一组作用在平面内，另一组垂直于平面。为此，在计算各个单元的结点载荷列阵 （包括等效结点力）可以直接引用第二章和第七章中所叙述的载荷计算的相应公式。各个单元的结点载荷列阵图8-4 矩形单元局部坐标 求得后，建立变换矩阵公式，从而把 转换到整体坐标系中去求出在整体坐标下的单元节点载荷列阵，然后经各单元的简单叠加可以求出结构在整体坐标下的节点载荷列阵。显然，平面单元在局部坐标系中，结点i有五个广义位移：即 ，其中前两个对应于平面应力问题，后三个对应于平板弯曲问题。类似地，所对应的结点力列阵 显然，在上式中 实际上总是等于零的。为了经坐标变换后不影响在整体坐标系中对各特征量的计算，我们引进(8-4)容易看出，把以上结点位移和结点力变换到整体坐标中后，他的结点位移和结点力列阵具有如下形式上式右端的前三项分别表示位移和力，后三项分别表示转角和力矩，它们都是有明显物理意义的矢量。因此，（8-4）式和（8-5）式之间的坐标变换公式是 式中而(b)(8-5)(8-6)(a)于是，壳体单元e在局部坐标下的结点位移列阵是或而所对应的单元节点力（包括等效节点力）列阵是或式中n=3是对应于三角形单元；n=4对应于四边形单元。本节以下的n所指的意义均是如此，不再重复说明。(c)(d)(e)(f)4.建立局部坐标系中的单元刚度矩阵 ，从而求出整体坐标系中的单元刚度矩阵。如果将单元刚度矩阵 和 对应于单元节点划分为nn个子矩阵，每个子矩阵都是66的，于是 的子矩阵有如下形式式中 和 分别是平面应力问题和平面弯曲问题的相应子矩阵，它们是22和33矩阵。图8-5 示出了在局部坐标系中三角形壳体单元刚度矩阵用平面应力和平面弯曲刚度矩阵的构成方法。(8-7)图8-5 三角形壳体单元刚度矩阵用平面应力和平板弯曲刚度矩阵的构成方法 单元e中任意结点i的平衡方程，在两个坐标系中分别为，式中 是刚度矩阵 的子矩阵。而对于局部坐标和整体坐标之间的变换公式是把(h)式代入(g)式得(g)(h)将公式(g)中的第一式左乘矩阵 ，并且同上式进行比较，可以得到由于 是正交阵，容易证明 也是正交阵，即 。这样就得到关于矩阵 的转换公式(8-8)5.集和单元刚度矩阵及等效结点力。线作简单求和然后将它们放入整体刚度矩阵K和等效结点荷载列阵 的相应位置上去。6.修改整体刚度矩阵，然后求解平衡方程式中 是总的结点位移列阵。特别值得注意，在局部坐标系中单元刚度矩阵K对于三角形单元它的第6、12、及18行和列全是零元素，对于四边形单元它的第6、12、18及24行和列全是零元素，其原因是转角 ，并不包含在平面应力的单元结点位移列阵中。当所有在一个结点相连接的单元共面时，壳体结构的刚度矩阵将是奇异的。避免这个奇异性的一个办法是引入关于壳体法线的转动为零的附加条件。7计算应力。首先是按照公式 求出局部坐标系中的结点位移，再按第二章中所给出的公式计算应力 、和；通过第七章所给出的公式计算 、和 进而求得应力和。、和 。于是，壳体应力可以由简单的叠加求得；即，对于一个壳体结构，如果采用上节所述的平面单元，将会引起几何上的离散误差。人们希望采用曲面单元来描述壳体的真正几何形状，使之用不太多的单元来替代复杂形状的壳体，并得到具有一定精度的解答。另外，在薄壳理论中都是用中面位移来表示中面转动。正如在第七章中所述，这将要求在单元交界面上有横向位移及其一阶导数的连续性，于是增加了选择位移模式的困难。如果考虑横向剪切变形的影响就可以认为中面转动是独立变量而不依赖于位移的一阶导数。因此，只要利用单元交界面上位移函数的连续性就可以了，并不要求其一阶导数的连续性。现在我们来论述一个考虑横向剪切影响的曲面单元，称为八结点四十个自由度的一般壳体单元，如图8-6所示。在图8-6中所示的壳单元，象空间等参数单元一样引进一个自然坐标系 。命 为壳体中面上的曲线坐标；对应于 的表面称为顶面（或上表面），对应于 的表面称为底面（或下表面）。在单元的中面上选取八个点称为结点，过各结点i(i=1,2,8)作中面的法线，交顶面和底面的点称为结点i的对点。结点i相对应的对点，它的整体坐标值分别记作图8-6 八结点四十个自由度 的一般壳体单元一.单元几何形状的确定于是，中面上的结点i的整体坐标值是显然，结点i处的中面法线方向可以由下列单位矢量所确定(8-10)(8-9)式中l3i、m3i和n3i是结点i处中面法线方向对于整体坐标轴oxyz的方向余弦，而hi是结点i处的壳体厚度，即(a)结点i处法线上任意点的整体坐标值，可以通过矢量相加得到（图8-7），即(b)于是，单元内任意点的坐标值可以通过形函数 的插值得到，即式中形函数 由公式（5-1）表示。这样，我们就可以通过八对点的整体坐标值，按照（8-11）式近似地确定了单元的形状。当 时，分别确定上下表面各点；当 时，确定中面各点；而单元的侧表面是由中面法线（或近似的中面法线）所构成。(8-11)假设中面法线变形后仍为直线，但是不再是变形后的中面法线。也就是说，中面法线有绕两不同轴的转动，这两不同轴分别垂直于法线。设V1i和V2i分别是这两轴上的单位矢量，显然它们具有一定的任意性。为了唯一地确定起见，不妨设式中(c)(8-12)二.位移模式于是有而另一轴的单位矢量可以由V3i和V1i的矢量积求得注意，在（8-12）式中如果i iV V3i=0，则可用j jV V3i来代替，而(j)(d)(e)(8-13)设结点i处的中面法线V3i绕V1i和V2i两轴的转角分别为 和 ，如图8-8所示。于是，转动矢量可以写成图8-8 结点的自由度(8-14)现在我们来计算结点i处法线上任意点的位移值。这可以通过结点i在整体坐标中的位移ui、vi、wi以及法线绕i点转动而产生的位移相加而得到；关于转动引起的位移可以根据运动学中的公式计算。于是，在结点i处法线上各点在整体坐标系中的位移是显然 。于是，上式就可以写成(h)(g)这样，利用插值方法可以得到单元内任意结点的位移列阵是式中(8-16)(8-15)可以将（8-15）式写成标准形式式中I是三阶单位阵。(8-17)(i)在整体坐标系中，利用几何方程和（8-17）式可以将应变列阵写成标准形式是（8-18）三.应变计算式中而(j)(8-19)因此，为了计算应变需要算出六个偏导数 、和 、前三个偏导数可以按照公式（5-54）进行计算，后三个则按定义分别是J-1的第三列的三个元素，即由（5-58）式得到(8-20)出现在所提到的各式中的S S、T T和V分别按公式（5-55）所定义。将（8-11）式代入（5-55）式得(8-21)我们还假设中面法线上的线段不伸长也不缩短。如果在单元中面各点处建立局部的直角坐标系 ，而 轴是中面法线方向，于是有 ；这里 表示 方向的位移分量。建立局部坐标系 的关键是中面上各点的法线方向余弦矢量V3，一个合理的简单做法是由单元各结点处中面法线的方向余弦矢量V3i通过内插法构造而成，即(8-22)四.坐标变换中面法线的方向余弦矢量构成后，可以按照如下方法作出 轴及轴的方向余弦矢量V1和V2于是，在 坐标系中的矢量都可以通过下列变换矩阵变换为整体坐标系oxyz中的矢量。(8-23)(k)在局部坐标系中，它的应变列阵是根据应变张量的定义，立刻可以写出应变的坐标形式式中 (8-26)(8-24)(8-25)在壳体中还假设 ，于是在局部坐标中有如下的应力应变关系式中以及(8-28)(8-27)(1)五.弹性矩阵的变换 在整体坐标系oxyz中，应力应变关系虽然仍可写成如下形式式中 。但是弹性矩阵D不再像以前那样具有简单的形式，而需要用坐标将 变换成D。由于这两组应力表示同一个应力状态，因此在任意虚位移上，两个坐标系中单位体积内的应力所做的虚功必须相等，即是(8-30)(8-29)将（8-27）式代入上式并利用（8-25）式，由于 的任意性得和（8-29）式相比较，得到这就是坐标变换情况下弹性矩阵的变换公式。可以证明这一形式具有普遍意义。(m)(8-31)习惯下，我们要计算的壳体应力是对应于局部坐标系 的。利用（8-25）式和（8-18）式，公式（8-27）可以写成其中(8-33)(8-32)根据（8-19）式和（8-26）式，很容易地算出矩阵 的显式。为了使公式书写不致太长，我们给出它的分块形式，即(8-34)六.应力计算式中(n)(o)(p)(q)而 (j=1,2,3)(r)于是(8-35)式中(s)在进行强度校核时，只要计算壳体上下表面 上的应力，即(8-36)容易算出 (8-37)(8-38)若将单元刚度矩阵写成如（5-61）式的形式，于是每个55的子矩阵可以按照下式进行计算把（8-31）式代入上式，于是有(8-40)(8-39)七.单元刚度矩阵式中利用（8-33）、（8-34）和（8-35）式可以得到若利用(n)(q)，不难将上式各子矩阵的乘积之和表示为显式。(8-41)(8-42)假设对应于结点i的五个自由度 的广义力是 1当单元内作用的体积力为 时，等效结点力的计算公式是式中Ni由(i)式给出。(8-43)八.等效结点力的计算 2当单元某面上作用有分布表面力 时，等效结点力是式中曲面积分是在单元的某一面上进行，其上作用有分布力。特别是在 的表面作用有法向荷载q的情形，可以按照下式计算(8-45)(8-44)在确定单元形状时，我们使用了八对点的整体坐标，每双对点都在一条中面法线上。在实际计算时，这一要求希望能尽可能地得到满足，因为关于变形的一系列推导，是与这个要求密切有关的。经过实际计算与一些理论分析证明，当使用222的高斯积分则计算刚度矩阵（8-40）式时，上述壳单元不仅适用于厚壳，而且适用于薄壳。作为一般规律，应力都在积分点上进行计算，然后外推到结点上去。九.其他