一. 单元几何形状的确定 在图8-6中所示的壳单元，象空间等参数单元一样引进一个自然 坐标系 oξηζ 。命 ξ ,η 为壳体中面上的曲线坐标；对应于 ζ = 1 的表 面称为顶面（或上表面），对应于 ζ = −1 的表面称为底面（或下表 面）。在单元的中面上选取八个点称为结点，过各结点i(i=1,2,…,8) 作中面的法线，交顶面和底面的点称为结点i的对点。结点i相对应 的对点，它的整体坐标值分别记作
[ ]
[k ] = [λ ] [k ' ][λ ]
ij ij
T
(8-8)
5.集和单元刚度矩阵及等效结点力。线作简单求和
∑ [k ]
ne ij e =1
∑ [R ]
i e =1
ne
然后将它们放入整体刚度矩阵[K]和等效结点荷载列阵 {R} 的相应位 置上去。
6.修改整体刚度矩阵，然后求解平衡方程
[ K ] {δ } = {R}
图8-1 任意壳体作为平面三角形单元的集合
图8-2 圆柱壳作为平面矩形单元的集合
壳体平面单元的应力状态是由平面应力和弯曲应力的叠加而成 的，因此在构造壳体平面单元时，只要将第二章和第七章所讨论的 相应单元进行简单的组合就可以了。同样，前述二章所导出的刚度 矩阵可作为建立壳体平面单元刚度矩阵的基础。 现在把平面单元的计算步骤归纳如下 1. 划分单元，选定整体坐标系 oxyz ，定出节点在整体坐标系中 的坐标值。 2. 对于各个单元利用节点坐标值，建立一个局部坐标系 ox' y ' z ' 例如三角形单元123，可以选取节点1为局部坐标系的原点，并且以 1-2边为 x ' 轴的正方向，如图8-3所示。于是，x ' 方向的单位e1求得 是
{δ i } = [λ ]{δ 'i }
式中 而
{Fi } = [λ ]{F 'i }
(8-6)
t 0 [λ ] = 0 t
(a)
[ t ] = [ e1
e2
e3 ]
(b)
于是，壳体单元e在局部坐标下的结点位移列阵是
{δ '}
或
e
= δ
[
T '1
δ
'T 2
⋯ δ
T T 'n
]
(c) (d)
[ ] [ ]
p p p k '11 k '12 k '13
p k '11
p k '12
b k '11
b k '12
p k ' 22
p k '13
p p k ' 21 k ' 22 k ' 23
p p p k ' 31 k ' 32 k ' 33 p k ' 21
p
b k '13
p k ' 23
b k '11
显然，结点i处的中面法线方向可以由下列单位矢量所确定
(8-9)
xi l 3i xi 1 V3i = m3i = y i − y i n3i hi z i 顶 z i 底
{δ i } = [ ui {Fi } = [ U i
vi
wi
θ xi θ yi θ zi ]T
M θxi M θyi M θzi
Vi Wi
]
(8-5)
T
上式右端的前三项分别表示位移和力，后三项分别表示转角和力矩， 它们都是有明显物理意义的矢量。因此，（8-4）式和（8-5）式之 8-4 8-5 间的坐标变换公式是
y' (η )
4 3
公式，从而把 {R' i } 转换到整 体坐标系中去求出在整体坐 标下的单元节点载荷列阵， 然后经各单元的简单叠加可 以求出结构在整体坐标下的 节点载荷列阵。
z y
z'
o'
1 2
x' (ξ )
x o
图8-4 矩形单元局部坐标
显然，平面单元在局部坐标系中，结点i有五个广义位移：即
12 e1 = 12
(8-1)
z
y'
3
z'
x'
y
1
2
o
x
图8-3 三角形单元局部坐标系
式中 12 = l12 =
(x2 − x1 )2 + ( y2 − y1 )2 + (z2 − z1 )2
12 × 13 e3 = 12 × 13
是矢量12的长
度。取单元的外法线方向作为 z ' 轴的正方向，于是它的单位矢量
{δ 'i } = [ u 'i
{F 'i } = [ U 'i
v 'i
w'i θ ' xi
M 'θxi
θ ' yi θ ' zi ]
M 'θyi
T
V 'i W 'i
M 'θzi
]
(8-4)
T
显然，在上式中 M 'θzi 实际上总是等于零的。
容易看出，把以上结点位移和结点力变换到整体坐标中后，他 的结点位移和结点力列阵具有如下形式
7．计算应力。首先是按照公式 {δ ' i } = [λ ] T {δ i } 求出局部坐标系
p 中的结点位移，再按第二章中所给出的公式计算应力 σ xp 、σ yp 和 τ xy
b ；通过第七章所给出的公式计算 M x 、M xy 和 M xy 进而求得应力 σ x
、σ y 和 τ xy 。于是，壳体应力可以由简单的叠加求得；即
{δ ' i } = [ u ' i
v' i
w' i θ ' xi θ ' yi ] T ，其中前两个对应于平面应力问题，
T M 'θyi 。为了经坐标变换后不影响
后三个对应于平板弯曲问题。类似地，所对应的结点力列阵
{F 'i } = [U 'i
V 'i
W 'i
M 'θxi
]
在整体坐标系中对各特征量的计算，我们引进
xi yi zi 顶
xi yi zi 底
图8-6 八结点四十个自由度 的一般壳体单元
于是，中面上的结点i的整体坐标值是
xi xi xi 1 yi = yi + yi zi 2 zi 顶 zi 底
壳体实质上是从平板演变而来的，它的中面是一个曲面。在分 析壳中应力时，虽然平板的基本假定同样有效，但是壳体的变形有 着很大程度的不同，它除了弯曲变形外还存在中面变形。因而，壳 中内力包括有弯曲内力和中面内力。 应用有限单元法分析壳体结构时，广泛地采用了平面单元和曲 面的单元这两类壳体单元。本章首先介绍平面单元，它是平面应力 问题和平板弯曲问题的组合；这种单元虽然简单，但是相当有效。 然后讨论一个考虑横向剪切影响的曲面单元，称为八结点40个自由 度的一般壳单元，可以适用于厚壳和薄壳。
b b k '12 k '13
k' 21 k 'b
p k ' 31
k' 22 k 'b
p k ' 32
p k ' 33
k 23 k'' b
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 三角形壳体单元刚度矩阵用平面应力和平板弯曲刚度矩阵的构成方法
将壳体曲面划分为有限个单元，它们都是曲面单元。但是在单 元细分时，用平面单元组成的一个单向或双向折板来近似壳体的几 何形状将会得到良好的结果。通常对于任意形状的壳体，采用三角 形单元比较方便，如图8-1所示。如果在壳体上容易找到同一平面上 的四个点，可以采用平面四边形单元。例如具有正交边界的柱面壳 体，如图8-2所示。
容易看出，矢量12和13的矢性积的模等于三角形面积∆的一倍，即 |12 13 12×13 12 13|=2∆。最后，按右手定则可以决定y轴的正方向,它的单位 矢量e2是 e2 = e3 ×e1 ，它的三个角点的局部坐标值是很容易确定的。 对于柱面上的矩形单元，局部坐标的原点 o '选在矩形的形心， 通常选 x ' 轴和x轴均沿柱面母线方向。如图8-4中所示，由矢量12确 定单位矢量e1,再由矢量14确定单位矢量e2，于是e3 = e1 ×e2。 (8-3) 利用上述方法确定的局部坐标系，三角形单元123是在 x' y ' 平面内
(8-10)
式中l3i、m3i和n3i是结点i处中面法线方向对于整体坐标轴oxyz的方向 余弦，而hi是结点i处的壳体厚度，即
hi =
(xi顶 − xi底 ) + ( yi顶 − yi底 ) + (z i顶 − z i底 )
2 2
2
(a)
结点i处法线上任意点的整体坐标值，可以通过矢量相加得到 （图8-7），即
3.对于各个单元，确立在局部坐标系 ox' y ' z ' 中的结点载荷列阵
{R'i } 。壳体载荷可以分解成二组：一组作用在平面内，另一组垂直 于平面。为此，在计算各个单元的结点载荷列阵 {R' i }（包括等效结
点力）可以直接引用第二章和第七章中所叙述的载荷计算的相应公 式。 各个单元的结点载荷列阵 {R'i } 求得后，建立变换矩阵