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

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 示出了在局部坐标系中三角 形壳体单元刚度矩阵用平面应力和平面弯曲刚度矩阵的构成方法。
1. 划分单元，选定整体坐标系 oxyz ，定出节点在整体坐标系中 的坐标值。
2. 对于各个单元利用节点坐标值，建立一个局部坐标系 ox' y ' z ' 例如三角形单元123，可以选取节点1为局部坐标系的原点，并且以
1-2边为 x ' 轴的正方向，如图8-3所示。于是，x ' 方向的单位e1求得
是
单元e中任意结点i的平衡方程，在两个坐标系中分别为
R'i k 'ij 'i ， Ri kij i
n
n
(g)
式中 kij 是刚度矩阵 k 的子矩阵。而对于局部坐标和整体坐标之 间的变换公式是

j 1
j 1
'
j j
Ri R'i
n
(h)
把(h)式代入(g)式得
R'i kij 'i
j 1
将公式(g)中的第一式左乘矩阵 ，并且同上式进行比较，可
以得到
k 'ij kij
由于 t 是正交阵，容易证明 也是正交阵，即 T 1 。这样 就得到关于矩阵 kij 的转换公式
系中的单元刚度矩阵。如果将单元刚度矩阵 k 和 k '对应于单元节 点划分为n×n个子矩阵，每个子矩阵都是6×6的，于是 k ' 的子矩 阵有如下形式
4.建立局部坐标系中的单元刚度矩阵 k ' ，从而求出整体坐标
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

e


T 1

T 2

T T n

(d)
而所对应的单元节点力（包括等效节点力）列阵是
R' F ' Q'
或
e


T R'1
R'T 2

T T Rn
T T R' nBiblioteka (e)Re


T R1
T R2


(f)
式中n = 3是对应于三角形单元；n = 4对应于四边形单元。本节以下 的n所指的意义均是如此，不再重复说明。
的附加条件。
7．计算应力。首先是按照公式 'i T i 求出局部坐标系
p 中的结点位移，再按第二章中所给出的公式计算应力 xp 、 yp 和 xy
b ；通过第七章所给出的公式计算 M x 、M xy 和 M xy 进而求得应力 x
、 y 和 xy 。于是，壳体应力可以由简单的叠加求得；即
一. 单元几何形状的确定 在图8-6中所示的壳单元，象空间等参数单元一样引进一个自然 坐标系 o 。命 , 为壳体中面上的曲线坐标；对应于 1 的表 面称为顶面（或上表面），对应于 1 的表面称为底面（或下表 面）。在单元的中面上选取八个点称为结点，过各结点i(i=1,2,…,8) 作中面的法线，交顶面和底面的点称为结点i的对点。结点i相对应 的对点，它的整体坐标值分别记作
式中 是总的结点位移列阵。特别值得注意，在局部坐标系中单 元刚度矩阵[K]对于三角形单元它的第6、12、及18行和列全是零元 素，对于四边形单元它的第6、12、18及24行和列全是零元素，其 原因是转角 ' zi ，并不包含在平面应力的单元结点位移列阵中。当 所有在一个结点相连接的单元共面时，壳体结构的刚度矩阵将是奇 异的。避免这个奇异性的一个办法是引入关于壳体法线的转动为零
通常选 x ' 轴和x轴均沿柱面母线方向。如图8-4中所示，由矢量12确 定单位矢量e1,再由矢量14确定单位矢量e2，于是e3 = e1 ×e2。
3.对于各个单元，确立在局部坐标系 ox' y ' z ' 中的结点载荷列阵
R'i 。壳体载荷可以分解成二组：一组作用在平面内，另一组垂直 于平面。为此，在计算各个单元的结点载荷列阵 R' i （包括等效结

p p p k '11 k '13 k '12
p k '11
p k '12
b k '11
b k '12 p k ' 22
p k '13
p p k ' 22 k ' 23 k ' 21
p p p k ' 31 k ' 32 k ' 33 p k ' 21
p
b k '13
p k ' 23

k k'
ij ij
T
(8-8)
5.集和单元刚度矩阵及等效结点力。线作简单求和
k
ij e 1
ne
R
i e 1
ne
然后将它们放入整体刚度矩阵[K]和等效结点荷载列阵 R 的相应位 置上去。
6.修改整体刚度矩阵，然后求解平衡方程
K R
i 'i
式中 而
Fi F 'i
(8-6)
t e1
t 0 0 t
(a)
e2
e3
(b)
于是，壳体单元e在局部坐标下的结点位移列阵是
'
或
e


T '1

'T 2

T T 'n

(c)
xi yi zi 顶
xi yi zi 底
图8-6 八结点四十个自由度 的一般壳体单元
于是，中面上的结点i的整体坐标值是
x x xi 1 i i yi yi yi 2 z z z i i i 顶 底
(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），即
(8-11)
式中形函数 N i ( , ) 由公式（5-1）表示。 这样，我们就可以通过八对点的整体坐标值，按照（8-11）式
近似地确定了单元的形状。当 1 时，分别确定上下表面各点；
当 0 时，确定中面各点；而单元的侧表面是由中面法线（或近 似的中面法线）所构成。
b k '11
b b k '13 k '12
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 三角形壳体单元刚度矩阵用平面应力和平板弯曲刚度矩阵的构成方法
F 'i U 'i
V 'i W 'i
M 'xi

在整体坐标系中对各特征量的计算，我们引进
'i u'i
F 'i U 'i
v'i
w'i ' xi ' yi ' zi
M 'xi M 'yi

T
V 'i W 'i
M 'zi

(8-4)
T
显然，在上式中 M 'zi 实际上总是等于零的。
定精度的解答。另外，在薄壳理论中都是用中面位移来表示中面转
动。正如在第七章中所述，这将要求在单元交界面上有横向位移及 其一阶导数的连续性，于是增加了选择位移模式的困难。如果考虑 横向剪切变形的影响就可以认为中面转动是独立变量而不依赖于位 移的一阶导数。因此，只要利用单元交界面上位移函数的连续性就 可以了，并不要求其一阶导数的连续性。 现在我们来论述一个考虑横向剪切影响的曲面单元，称为八结 点四十个自由度的一般壳体单元，如图8-6所示。
以求出结构在整体坐标下的
节点载荷列阵。
显然，平面单元在局部坐标系中，结点i有五个广义位移：即
'i u'i
v' i
w'i ' xi ' yi T ，其中前两个对应于平面应力问题，
T M 'yi 。为了经坐标变换后不影响
后三个对应于平板弯曲问题。类似地，所对应的结点力列阵
z y
(b)
xi h i y V3i i 2 z i
V3i
i ( xi , y i , z i )
o
x
于是，单元内任意点的坐标值可以通过形函数 N i ( , ) 的插值得到，