第八章一般壳体问题的有限元法

第八章一般壳体问题的有限元法第一页，共四十八页，2022年，8月28日

将壳体曲面划分为有限个单元，它们都是曲面单元。但是在单元细分时，用平面单元组成的一个单向或双向折板来近似壳体的几何形状将会得到良好的结果。通常对于任意形状的壳体，采用三角形单元比较方便，如图8-1所示。如果在壳体上容易找到同一平面上的四个点，可以采用平面四边形单元。例如具有正交边界的柱面壳体，如图8-2所示。图8-1任意壳体作为平面三角形单元的集合图8-2圆柱壳作为平面矩形单元的集合§8-1平面壳体单元第二页，共四十八页，2022年，8月28日壳体平面单元的应力状态是由平面应力和弯曲应力的叠加而成的，因此在构造壳体平面单元时，只要将第二章和第七章所讨论的相应单元进行简单的组合就可以了。同样，前述二章所导出的刚度矩阵可作为建立壳体平面单元刚度矩阵的基础。现在把平面单元的计算步骤归纳如下

1.划分单元，选定整体坐标系oxyz，定出节点在整体坐标系中的坐标值。

2.对于各个单元利用节点坐标值，建立一个局部坐标系例如三角形单元123，可以选取节点1为局部坐标系的原点，并且以1-2边为轴的正方向，如图8-3所示。于是，方向的单位e1求得是(8-1)第三页，共四十八页，2022年，8月28日式中是矢量12的长度。取单元的外法线方向作为轴的正方向，于是它的单位矢量图8-3三角形单元局部坐标系第四页，共四十八页，2022年，8月28日容易看出，矢量12和13的矢性积的模等于三角形面积Δ的一倍，即|12×13|=2Δ。最后，按右手定则可以决定y轴的正方向,它的单位矢量e2是e2=e3×e1(8-3)利用上述方法确定的局部坐标系，三角形单元123是在平面内，它的三个角点的局部坐标值是很容易确定的。对于柱面上的矩形单元，局部坐标的原点选在矩形的形心，通常选轴和x轴均沿柱面母线方向。如图8-4中所示，由矢量12确定单位矢量e1,再由矢量14确定单位矢量e2，于是e3=e1×e2。第五页，共四十八页，2022年，8月28日

3.对于各个单元，确立在局部坐标系中的结点载荷列阵。壳体载荷可以分解成二组：一组作用在平面内，另一组垂直于平面。为此，在计算各个单元的结点载荷列阵（包括等效结点力）可以直接引用第二章和第七章中所叙述的载荷计算的相应公式。各个单元的结点载荷列阵图8-4矩形单元局部坐标求得后，建立变换矩阵公式，从而把转换到整体坐标系中去求出在整体坐标下的单元节点载荷列阵，然后经各单元的简单叠加可以求出结构在整体坐标下的节点载荷列阵。第六页，共四十八页，2022年，8月28日显然，平面单元在局部坐标系中，结点i有五个广义位移：即

，其中前两个对应于平面应力问题，后三个对应于平板弯曲问题。类似地，所对应的结点力列阵

显然，在上式中实际上总是等于零的。。为了经坐标变换后不影响在整体坐标系中对各特征量的计算，我们引进(8-4)第七页，共四十八页，2022年，8月28日容易看出，把以上结点位移和结点力变换到整体坐标中后，他的结点位移和结点力列阵具有如下形式上式右端的前三项分别表示位移和力，后三项分别表示转角和力矩，它们都是有明显物理意义的矢量。因此，（8-4）式和（8-5）式之间的坐标变换公式是

式中而(b)

(8-5)(8-6)(a)第八页，共四十八页，2022年，8月28日于是，壳体单元e在局部坐标下的结点位移列阵是或而所对应的单元节点力（包括等效节点力）列阵是或式中n=3是对应于三角形单元；n=4对应于四边形单元。本节以下的n所指的意义均是如此，不再重复说明。(c)(d)

(e)

(f)第九页，共四十八页，2022年，8月28日

4.建立局部坐标系中的单元刚度矩阵，从而求出整体坐标系中的单元刚度矩阵。如果将单元刚度矩阵和对应于单元节点划分为n×n个子矩阵，每个子矩阵都是6×6的，于是的子矩阵有如下形式式中和分别是平面应力问题和平面弯曲问题的相应子矩阵，它们是2×2和3×3矩阵。图8-5示出了在局部坐标系中三角形壳体单元刚度矩阵用平面应力和平面弯曲刚度矩阵的构成方法。(8-7)第十页，共四十八页，2022年，8月28日图8-5三角形壳体单元刚度矩阵用平面应力和平板弯曲刚度矩阵的构成方法第十一页，共四十八页，2022年，8月28日单元e中任意结点i的平衡方程，在两个坐标系中分别为，式中是刚度矩阵的子矩阵。而对于局部坐标和整体坐标之间的变换公式是把(h)式代入(g)式得(g)(h)第十二页，共四十八页，2022年，8月28日将公式(g)中的第一式左乘矩阵，并且同上式进行比较，可以得到由于是正交阵，容易证明也是正交阵，即。这样就得到关于矩阵的转换公式(8-8)第十三页，共四十八页，2022年，8月28日

5.集和单元刚度矩阵及等效结点力。线作简单求和然后将它们放入整体刚度矩阵[K]和等效结点荷载列阵的相应位置上去。第十四页，共四十八页，2022年，8月28日

6.修改整体刚度矩阵，然后求解平衡方程式中是总的结点位移列阵。特别值得注意，在局部坐标系中单元刚度矩阵[K]对于三角形单元它的第6、12、及18行和列全是零元素，对于四边形单元它的第6、12、18及24行和列全是零元素，其原因是转角，并不包含在平面应力的单元结点位移列阵中。当所有在一个结点相连接的单元共面时，壳体结构的刚度矩阵将是奇异的。避免这个奇异性的一个办法是引入关于壳体法线的转动为零的附加条件。第十五页，共四十八页，2022年，8月28日

7．计算应力。首先是按照公式求出局部坐标系中的结点位移，再按第二章中所给出的公式计算应力、和；通过第七章所给出的公式计算、和进而求得应力和。、和。于是，壳体应力可以由简单的叠加求得；即，第十六页，共四十八页，2022年，8月28日

对于一个壳体结构，如果采用上节所述的平面单元，将会引起几何上的离散误差。人们希望采用曲面单元来描述壳体的真正几何形状，使之用不太多的单元来替代复杂形状的壳体，并得到具有一定精度的解答。另外，在薄壳理论中都是用中面位移来表示中面转动。正如在第七章中所述，这将要求在单元交界面上有横向位移及其一阶导数的连续性，于是增加了选择位移模式的困难。如果考虑横向剪切变形的影响就可以认为中面转动是独立变量而不依赖于位移的一阶导数。因此，只要利用单元交界面上位移函数的连续性就可以了，并不要求其一阶导数的连续性。

现在我们来论述一个考虑横向剪切影响的曲面单元，称为八结点四十个自由度的一般壳体单元，如图8-6所示。8-2考虑横向剪切变形影响的壳体单元第十七页，共四十八页，2022年，8月28日

在图8-6中所示的壳单元，象空间等参数单元一样引进一个自然坐标系。命为壳体中面上的曲线坐标；对应于的表面称为顶面（或上表面），对应于的表面称为底面（或下表面）。在单元的中面上选取八个点称为结点，过各结点i(i=1,2,…,8)作中面的法线，交顶面和底面的点称为结点i的对点。结点i相对应的对点，它的整体坐标值分别记作图8-6八结点四十个自由度的一般壳体单元一.单元几何形状的确定第十八页，共四十八页，2022年，8月28日于是，中面上的结点i的整体坐标值是显然，结点i处的中面法线方向可以由下列单位矢量所确定(8-10)(8-9)第十九页，共四十八页，2022年，8月28日式中l3i、m3i和n3i是结点i处中面法线方向对于整体坐标轴oxyz的方向余弦，而hi是结点i处的壳体厚度，即

(a)结点i处法线上任意点的整体坐标值，可以通过矢量相加得到（图8-7），即(b)第二十页，共四十八页，2022年，8月28日于是，单元内任意点的坐标值可以通过形函数的插值得到，即式中形函数由公式（5-1）表示。这样，我们就可以通过八对点的整体坐标值，按照（8-11）式近似地确定了单元的形状。当时，分别确定上下表面各点；当时，确定中面各点；而单元的侧表面是由中面法线（或近似的中面法线）所构成。(8-11)第二十一页，共四十八页，2022年，8月28日

假设中面法线变形后仍为直线，但是不再是变形后的中面法线。也就是说，中面法线有绕两不同轴的转动，这两不同轴分别垂直于法线。设V1i和V2i分别是这两轴上的单位矢量，显然它们具有一定的任意性。为了唯一地确定起见，不妨设式中(c)

(8-12)二.位移模式第二十二页，共四十八页，2022年，8月28日于是有而另一轴的单位矢量可以由V3i和V1i的矢量积求得注意，在（8-12）式中如果i×V3i=0，则可用j×V3i来代替，而(j)(d)(e)

(8-13)第二十三页，共四十八页，2022年，8月28日设结点i处的中面法线V3i绕V1i和V2i两轴的转角分别为和，如图8-8所示。于是，转动矢量可以写成图8-8结点的自由度(8-14)第二十四页，共四十八页，2022年，8月28日现在我们来计算结点i处法线上任意点的位移值。这可以通过结点i在整体坐标中的位移ui、vi、wi以及法线绕i点转动而产生的位移相加而得到；关于转动引起的位移可以根据运动学中的公式计算。于是，在结点i处法线上各点在整体坐标系中的位移是显然。于是，上式就可以写成(h)(g)第二十五页，共四十八页，2022年，8月28日这样，利用插值方法可以得到单元内任意结点的位移列阵是式中(8-16)(8-15)第二十六页，共四十八页，2022年，8月28日可以将（8-15）式写成标准形式式中[I]是三阶单位阵。(8-17)(i)第二十七页，共四十八页，2022年，8月28日

在整体坐标系中，利用几何方程和（8-17）式可以将应变列阵写成标准形式是（8-18）三.应变计算第二十八页，共四十八页，2022年，8月28日式中而(j)(8-19)第二十九页，共四十八页，2022年，8月28日因此，为了计算应变需要算出六个偏导数、、和、前三个偏导数可以按照公式（5-54）进行计算，后三个则按定义分别是[J]-1的第三列的三个元素，即由（5-58）式得到(8-20)第三十页，共四十八页，2022年，8月28日出现在所提到的各式中的S、T和V分别按公式（5-55）所定义。将（8-11）式代入（5-55）式得(8-21)第三十一页，共四十八页，2022年，8月28日

我们还假设中面法线上的线段不伸长也不缩短。如果在单元中面各点处建立局部的直角坐标系，而轴是中面法线方向，于是有；这里表示方向的位移分量。建立局部坐标系的关键是中面上各点的法线方向余弦矢量V3，一个合理的简单做法是由单元各结点处中面法线的方向余弦矢量V3i通过内插法构造而成，即(8-22)四.坐标变换第三十二页，共四十八页，2022年，8月28日中面法线的方向余弦矢量构成后，可以按照如下方法作出轴及轴的方向余弦矢量V1和V2于是，在坐标系中的矢量都可以通过下列变换矩阵变换为整体坐标系oxyz中的矢量。(8-23)

(k)第三十三页，共四十八页，2022年，8月28日在局部坐标系中，它的应变列阵是根据应变张量的定义，立刻可以写出应变的坐标形式式中(8-26)(8-24)(8-25)第三十四页，共四十八页，2022年，8月28日

在壳体中还假设，于是在局部坐标中有如下的应力应变关系式中以及(8-28)(8-27)(1)五.弹性矩阵的变换第三十五页，共四十八页，2022年，8月28日在整体坐标系oxyz中，应力应变关系虽然仍可写成如下形式式中。但是弹性矩阵[D]不再像以前那样具有简单的形式，而需要用坐标将变换成[D]。由于这两组应力表示同一个应力状态，因此在任意虚位移上，两个坐标系中单位体积内的应力所做的虚功必须相等，即是(8-30)(8-29)第三十六页，共四十八页，2022年，8月28日将（8-27）式代入上式并利用（8-25）式，由于的任意性得和（8-29）式相比较，得到这就是坐标变换情况下弹性矩阵的变换公式。可以证明这一形式具有普遍意义。(m)(8-31)第三十七页，共四十八页，2022年，8月28日

习惯下，我们要计算的壳体应力是对应于局部坐标系的。利用（8-25）式和（8-18）式，公式（8-27）可以写成其中(8-33)(8-32)根据（8-19）式和（8-26）式，很容易地算出矩阵的显式。为了使公式书写不致太长，我们给出它的分块形式，即(8-34)六.应力计算第三十八页，共四十八页，2022年，8月28日式中(n)(o)(p)(q)第三十九页，共四十八页，2022年，8月28日而(j=1,2,3)(r)于是(8-35)第四十页，共四十八页，2022年，8月28日式中(s)第四十一页，共四十八页，2022年，8月28日在进行强度校核时，只要计算壳体上下表面上的应力，即(8-36)第四十二页，共四十八页，2022年，8月28日容易算出(8-37)(8-38)第四十三页，共四十八页，2022年，8月28日

若将单元刚度矩阵写成如（5-61）式的形式，于是每个5