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

b k '11
b b k '12 k '13
k 'b 21
p k ' 31
k 'b 22
p k ' 32
p k ' 33
k 'b 23
k 'b 21
k 'b 22
k 'b 32
k 'b 23
k 'b 33
k 'b 31
k 'b 31
k 'b 32
k 'b 33
图8-5 三角形壳体单元刚度矩阵用平面应力和平板弯曲刚度矩阵的构成方法
容易看出，把以上结点位移和结点力变换到整体坐标中后，他 的结点位移和结点力列阵具有如下形式
T u v w i i i i xiyizi T F U M M i iV iW i M xi yi zi
(8-5)
上式右端的前三项分别表示位移和力，后三项分别表示转角和力矩， 它们都是有明显物理意义的矢量。因此，（8-4）式和（8-5）式之 间的坐标变换公式是
n
(h)
把(h)式代入(g)式得
R ki j 'i 'i
j 1
将公式(g)中的第一式左乘矩阵 ，并且同上式进行比较，可
以得到
k ' k ij ij
1 T 由于 t 是正交阵，容易证明 也是正交阵，即 。这样

(f)
式中n = 3是对应于三角形单元；n = 4对应于四边形单元。本节以下 的n所指的意义均是如此，不再重复说明。
系中的单元刚度矩阵。如果将单元刚度矩阵 k 和 k ' 对应于单元节 点划分为n×n个子矩阵，每个子矩阵都是6×6的，于是 k ' 的子矩 阵有如下形式
4.建立局部坐标系中的单元刚度矩阵 k ' ，从而求出整体坐标

(c)
或

e
e
T T T T 1 2 n

(d)
而所对应的单元节点力（包括等效节点力）列阵是

或
T TT T R ' F ' Q ' R ' R ' R ' 1 2 n

(e)

e
T T T T R R R 1R 2 n
是
12 e1 12
(8-1)
z
z'
y'
3
x'
y
1
2
o
x
图8-3 三角形单元局部坐标系
2 2 2 式中 12 是矢量12的长 l x x y y z z 12 2 1 2 1 2 1
度。取单元的外法线方向作为 z ' 轴的正方向，于是它的单位矢量
壳体实质上是从平板演变而来的，它的中面是一个曲面。在分
析壳中应力时，虽然平板的基本假定同样有效，但是壳体的变形有
着很大程度的不同，它除了弯曲变形外还存在中面变形。因而，壳 中内力包括有弯曲内力和中面内力。 应用有限单元法分析壳体结构时，广泛地采用了平面单元和曲 面的单元这两类壳体单元。本章首先介绍平面单元，它是平面应力 问题和平板弯曲问题的组合；这种单元虽然简单，但是相当有效。 然后讨论一个考虑横向剪切影响的曲面单元，称为八结点40个自由 度的一般壳单元，可以适用于厚壳和薄壳。
现在把平面单元的计算步骤归纳如下
1. 划分单元，选定整体坐标系 oxyz ，定出节点在整体坐标系中 的坐标值。
2. 对于各个单元利用节点坐标值，建立一个局部坐标系 ox' y' z' 例如三角形单元123，可以选取节点1为局部坐标系的原点，并且以
1-2边为 x ' 轴的正方向，如图8-3所示。于是，x ' 方向的单位e1求得
(8-10)
式中l3i、m3i和n3i是结点i处中面法线方向对于整体坐标轴oxyz的方向
余弦，而hi是结点i处的壳体厚度，即
2 2 2 (a) h x x y y z z i i 顶 i 顶 i 顶 i 底 i 底 i 底




结点i处法线上任意点的整体坐标值，可以通过矢量相加得到 （图8-7），即
p 0 0 0 k 'rs 0 0 0 0 0 k'rs 0 0 k 'b rs 0 0 0 0 0 0 0


0 0 0 0 0 0
r, s 1 ,2 , , n
(8-7)
式中 k' p 和 k'b 分别是平面应力问题和平面弯曲问题的相应子 rs rs 矩阵，它们是2×2和3×3矩阵。图8-5 示出了在局部坐标系中三角 形壳体单元刚度矩阵用平面应力和平面弯曲刚度矩阵的构成方法。
在整体坐标系中对各特征量的计算，我们引进


T ' u ' ' ' ' ' i i v i w i ' xi yi zi
Байду номын сангаас

T F ' U ' ' ' ' M ' M ' i iV iW iM xi yi zi

(8-4)
显然，在上式中 M 'zi 实际上总是等于零的。
显然，结点i处的中面法线方向可以由下列单位矢量所确定
(8-9)
x x l3 i 1 i i V m y y 3 i 3 i i i h i n z z i i i 3 顶 底
体，如图8-2所示。
图8-1 任意壳体作为平面三角形单元的集合
图8-2 圆柱壳作为平面矩形单元的集合
壳体平面单元的应力状态是由平面应力和弯曲应力的叠加而成
的，因此在构造壳体平面单元时，只要将第二章和第七章所讨论的
相应单元进行简单的组合就可以了。同样，前述二章所导出的刚度 矩阵可作为建立壳体平面单元刚度矩阵的基础。
就得到关于矩阵 k ij 的转换公式

k k '
ij ij
T
(8-8)
5.集和单元刚度矩阵及等效结点力。线作简单求和
k
ij e 1
ne
R
i e 1
ne
然后将它们放入整体刚度矩阵[K]和等效结点荷载列阵 R 的相应位 置上去。
6.修改整体刚度矩阵，然后求解平衡方程
xi yi zi 顶
xi yi zi 底
图8-6 八结点四十个自由度 的一般壳体单元
于是，中面上的结点i的整体坐标值是
x x xi 1 i i yi yi yi 2z z z i i i 顶 底
以求出结构在整体坐标下的
节点载荷列阵。
'xi ' u ' ' ' i i v i w i
显然，平面单元在局部坐标系中，结点i有五个广义位移：即
T 'yi ，其中前两个对应于平面应力问题，
后三个对应于平板弯曲问题。类似地，所对应的结点力列阵
T 。为了经坐标变换后不影响 F ' U ' V ' W ' M ' M ' i i i i xi yi
定精度的解答。另外，在薄壳理论中都是用中面位移来表示中面转
动。正如在第七章中所述，这将要求在单元交界面上有横向位移及 其一阶导数的连续性，于是增加了选择位移模式的困难。如果考虑 横向剪切变形的影响就可以认为中面转动是独立变量而不依赖于位 移的一阶导数。因此，只要利用单元交界面上位移函数的连续性就 可以了，并不要求其一阶导数的连续性。 现在我们来论述一个考虑横向剪切影响的曲面单元，称为八结 点四十个自由度的一般壳体单元，如图8-6所示。
12 13 e3 12 13
容易看出，矢量12和13的矢性积的模等于三角形面积Δ的一倍，即 |12×13|=2Δ。最后，按右手定则可以决定y轴的正方向,它的单位 矢量e2是 e2 = e3 ×e1 (8-3)
利用上述方法确定的局部坐标系，三角形单元123是在 x ' y ' 平面内 ，它的三个角点的局部坐标值是很容易确定的。 对于柱面上的矩形单元，局部坐标的原点 o ' 选在矩形的形心，
F F ' ' i i i i
式中
(8-6)
而
t e 1e 2e 3
t 0 0 t
(a) (b)
于是，壳体单元e在局部坐标下的结点位移列阵是
'
e
T T T T ' ' 1 ' 2 n
点力）可以直接引用第二章和第七章中所叙述的载荷计算的相应公
式。 各个单元的结点载荷列阵 R'i 求得后，建立变换矩阵
z y
1
y ' ( )
公式，从而把 R'i 转换到整
z'
4
o'
3
x ' ( )
体坐标系中去求出在整体坐
标下的单元节点载荷列阵， 然后经各单元的简单叠加可
2
x o
图8-4 矩形单元局部坐标
一. 单元几何形状的确定 在图8-6中所示的壳单元，象空间等参数单元一样引进一个自然 坐标系 o 。命 , 为壳体中面上的曲线坐标；对应于 1 的表 面称为顶面（或上表面），对应于 1 的表面称为底面（或下表 面）。在单元的中面上选取八个点称为结点，过各结点i(i=1,2,…,8) 作中面的法线，交顶面和底面的点称为结点i的对点。结点i相对应 的对点，它的整体坐标值分别记作
的附加条件。
T ' 7．计算应力。首先是按照公式 i i 求出局部坐标系
中的结点位移，再按第二章中所给出的公式计算应力 xp 、 yp 和 ；通过第七章所给出的公式计算 M x 、M 、
b y
p xy
b x
xy
和M
xy
进而求得应力
和
b xy
。于是，壳体应力可以由简单的叠加求得；即
单元e中任意结点i的平衡方程，在两个坐标系中分别为
R'i k'ij 'i ， Ri kij i
n
n
(g)
式中 k ij
是刚度矩阵 k 的子矩阵。而对于局部坐标和整体坐标之
j 1
j1
间的变换公式是
'
j j
R R ' i i
p b p b p b xy xy 。 y y 和 xy x x x ， y
对于一个壳体结构，如果采用上节所述的平面单元，将会引起
几何上的离散误差。人们希望采用曲面单元来描述壳体的真正几何 形状，使之用不太多的单元来替代复杂形状的壳体，并得到具有一
将壳体曲面划分为有限个单元，它们都是曲面单元。但是在单 元细分时，用平面单元组成的一个单向或双向折板来近似壳体的几 何形状将会得到良好的结果。通常对于任意形状的壳体，采用三角 形单元比较方便，如图8-1所示。如果在壳体上容易找到同一平面上 的四个点，可以采用平面四边形单元。例如具有正交边界的柱面壳
K R
式中 是总的结点位移列阵。特别值得注意，在局部坐标系中单 元刚度矩阵[K]对于三角形单元它的第6、12、及18行和列全是零元 素，对于四边形单元它的第6、12、18及24行和列全是零元素，其 原因是转角 ' zi ，并不包含在平面应力的单元结点位移列阵中。当 所有在一个结点相连接的单元共面时，壳体结构的刚度矩阵将是奇 异的。避免这个奇异性的一个办法是引入关于壳体法线的转动为零

p p p k '11 k '13 k '12
p p p k ' 22 k ' 21 k ' 23
p p p k ' 31 k ' 32 k ' 33
p k '11
p k '12
p k '13
b k '12
b k '11
p k ' 21 p k ' 22
b k '13
p k ' 23
通常选 x ' 轴和x轴均沿柱面母线方向。如图8-4中所示，由矢量12确 定单位矢量e1,再由矢量14确定单位矢量e2，于是e3 = e1 ×e2。
3.对于各个单元，确立在局部坐标系 ox' y' z' 中的结点载荷列阵
R'i 。壳体载荷可以分解成二组：一组作用在平面内，另一组垂直 于平面。为此，在计算各个单元的结点载荷列阵 R'i （包括等效结