壳单元的公式推导

s
（1-6）
（1-7） （1-8）
2、 厚板弯曲单元
位移函数为
u
=
⎧θ ⎪⎨θ
x y
⎫ ⎪ ⎬
=
[
N1I
N2I
N3I
N4I ]ae = Nae
⎪⎩ w ⎪⎭
其中 I 为三阶单位矩阵；单元节点位移向量为
⎧a1 ⎫
ae
=
⎪⎪⎨⎪aa23
⎪⎪⎬ , ⎪
⎪⎩a4 ⎪⎭
a
i=
⎧θ ⎪⎨θ
xi yi
⎫ ⎪ ⎬
⎪⎩ wi ⎪⎭
=
⎡⎢⎢kk12bb11 ⎢⎢⎢⎣kk34bb11
k1b2 k2b2 k3b2 k4b2
k1b3 k2b3 k3b3 k4b3
k1b4 k2b4 k3b4 k4b4
⎤ ⎥ ⎥ ⎥ ⎥ ⎥⎦
⎧⎪⎪⎨⎪aaa132bbb ⎪⎩a4b
⎫ ⎪⎪ ⎬ ⎪ ⎪⎭
=
K
a eb eb
（3-2）
{ } { } 其中 Fib = Fzi
[ ] [ ] ⎡⎣K e ⎤⎦node = R T ⎡⎣K e ⎤⎦ plane R
（5-1）
[R]是面内位移和节点位移的变换矩阵，即[u] = [R][u]
plane
node
（5-2）
[R]T 是节点力及力矩和面内力及力矩的变换矩阵，即[F ] = [R]T [F ]
node
plane
（5-3）
−
Ni
Bsi
=
⎢ ⎢
⎣⎢
0
0⎥
⎦
0 −Ni
∂Ni ⎤
∂x ∂Ni
⎥ ⎥ ⎥
,
∂y ⎦⎥
(i = 1, 2,3, 4)
根据复合函数求导法则，可得
（2-4）
⎧ ∂Ni
⎪⎪ ∂x
⎨ ⎪
∂NNi
⎪⎪ ∂ξ
⎨ ⎪
∂Ni
⎫ ⎪⎬⎪ , ⎪
⎩⎪ ∂y ⎭⎪ ⎩⎪ ∂η ⎭⎪
⎧ ∂x
x, z y, z z, z
⎤ ⎥ ⎥ ⎥ ⎥ ⎥⎦
⎧ ⎪⎪ ⎨ ⎪ ⎩⎪
x y z
⎫ ⎪⎪ ⎬ ⎪ ⎭⎪
=
φ
⎧ ⎪⎪ ⎨ ⎪ ⎩⎪
x⎫
y ⎪⎬⎪
z
⎪ ⎪⎭
其中φ 为局部坐标系 xyz 与整体坐标系 x y z 之间的方向余弦矩阵。
⎡λ 0 0 0⎤
节点位移变换公式为 ae
=
T ae
=
⎢ ⎢
0
⎢0
⎢ ⎣
⎪F3
⎪⎩F4
⎫ ⎪⎪ ⎬ ⎪ ⎪⎭
=
⎡⎢⎢kk1211 ⎢⎢k31 ⎢⎣k41
k12 k22 k32 k42
k13 k23 k33 k43
k14 k24 k34 k44
⎤ ⎥ ⎥ ⎥ ⎥ ⎥⎦
⎧⎪⎪⎨aa12 ⎪a3 ⎪⎩a4
⎫ ⎪⎪ ⎬ ⎪ ⎪⎭
=
K
eae
（3-4）
⎡⎢kijp ⎢
0 0 0 0⎤ 0 0 0 0⎥⎥
其中单元刚度矩阵的字块为 kij
=
⎢ ⎢
0
⎢
0 0
kibj
0⎥ ⎥
0⎥
⎢ ⎢
0
0
0⎥⎥
⎢⎣ 0 0 0 0 0 0⎥⎦
4、 单元局部坐标系 单元局部坐标系的建立如图 1 所示，取 G1G3 连线与 G2G4 连线的交点为坐标原点，x
轴沿 G1G3 连线与 G2G4 连线所形成的角平分线方向，正方向由 G1 指向 G2；y 轴垂直于 x 轴，位于单元平面内，正方向由 G1 指向 G4；z 轴垂直于 xoy 面，正方向按 G1、G2、G3、 G4，并符合右手法则。
K
e s
=
Ae BsT Ds Bs dA
（2-7）
⎡
⎤
⎢1 ν 0 ⎥
( ) 其中 Db
=
t3 12
D1
= 12
Et 3 1−ν 2
⎢⎢ν ⎢
1
⎥ 0 ⎥, 1−ν ⎥
Ds
= tD2
=
Gt k
⎡1 ⎢⎣0
0⎤ 1⎥⎦
⎢0 0
⎥
⎣
2⎦
∫ 等效节点载荷 Pe = N T pdA Ae
（2-8）
3、 壳元的单元刚度矩阵 在单元局部坐标系下，与面内变形有关的情况同平面应力问题，单元基本方程可写为
⎧ ∂Ni
⎪⎪ ∂x
⎨ ⎪
∂Ni
⎫ ⎪⎪ ⎬ ⎪
=
J
−1
⎧ ∂Ni
⎪⎪ ∂ξ
⎨ ⎪
∂Ni
⎫ ⎪⎬⎪ , ⎪
⎩⎪ ∂y ⎭⎪ ⎩⎪ ∂η ⎭⎪
⎧ ∂x
J
=
⎪⎪ ⎨ ⎪
∂ξ ∂x
⎩⎪ ∂η
∂y ⎫ ⎧
∑ ∂ξ
∂x
⎪⎪ ⎬ ⎪
=
⎪⎪ ⎨ ⎪
∑ ∂η ⎭⎪ ⎩⎪
∂Ni ∂ξ
xi
∂Ni ∂η
xi
∑ ∂Ni
∂ξ
0
0
⎥ ⎥
⎢
其中
D
=
E 1−ν
2
⎢0 ⎢ ⎢⎢0
⎢
⎢⎢⎣0
0 1−ν 2
00
00
0 1−ν 2k
0
⎥
0 0
⎥ ⎥ ⎥
=
⎡ ⎢ ⎣
D1 0
⎥
1
−ν
⎥ ⎥
2k ⎥⎦
0⎤
D2
⎥ ⎦
k 为考虑剪应力分布不均匀影响的系数，可取 1.20。
∫ ∫ 单元刚度矩阵
K
e bending
=
Kbe
+
K
e s
，
Kbe
=
Ae BbT Db Bb dA,
∂θx − ∂θ y
⎫ ⎪ ⎪ ⎬⎪⎪ , ⎪ ⎪ ⎪
⎪⎩ ∂y ∂x ⎪⎭
γ
=
⎧γ ⎨⎩γ
xz xz
⎫ ⎬ ⎭
=
⎧ ⎪⎪
∂w ∂x
−
θ
x
⎫ ⎪⎪
⎨ ⎪ ⎩⎪
∂w ∂y
−
θ
y
⎬ ⎪ ⎪⎭
（2-3）
⎧εx ⎫
代入几何方程可得 ε
=
⎪⎪⎪⎨γεxyy
⎪ ⎪⎪ ⎬
⎪⎪γ
xz
⎪ ⎪
=
⎧zκ ⎫
⎨ ⎩
γ
⎬, ⎭
Ke membrane
=
BT DBdΩ =
Ωe
BT DBtdA =
Ae
1 1 BT DB J tdξ dη
−1 −1
等效节点载荷
∫ ∫ ∫ ∫ 体力
Pfe =
N T fdΩ =
Ωe
N T ftdA =
Ae
1 1 NT f
−1 −1
J tdξ dη
∫ ∫ 面力
Ppe =
N T pdΓ =
Γe
N T ptds
⎪⎩γ yz ⎪⎭
κ = Bbae ,
γ = Bsae
[ ] [ ] 其中 Bb = Bb1 Bb2 Bb3 Bb4 , Bs = Bs1 Bs2 Bs3 Bs4
⎡ ⎢
−
∂Ni
⎢ ∂x
⎢ Bbi = ⎢ 0
⎢
⎢ ⎢ ⎣
−
∂Ni ∂y
0
− ∂Ni ∂y
− ∂Ni ∂x
⎤
0⎥ ⎥ ⎥
0⎥ , ⎥ ⎥
⎡ ⎢
理。其中一种方法就是在节点的单元刚度矩阵对应位置给以任意的虚拟刚度系数 Kθzθz ，这
样局部坐标系中θz 方向的平衡方程为 Kθzθzθz = 0 。经过坐标变化，整体坐标系中的该节点
平衡方程将满足有唯一解条件。由于θz 并不影响应力，而且与其他节点平衡方程无关，故
Kθzθz 可赋予任何非零值。
7、 壳单元的通用性 根据壳单元的厚度与面内尺寸的比例关系，壳可以分为薄壳和厚壳。Quad4 单元以厚壳
(i = 1, 2,3, 4)
（2-1） （2-2）
⎧εx ⎫
对厚板 ε
=
⎪⎨⎪⎪γεxyy ⎪⎪γ xz
⎪ ⎪⎪ ⎬ ⎪ ⎪
=
⎧ zκ
⎨ ⎩
γ
⎫ ⎬, ⎭
⎪⎩γ yz ⎪⎭
κ
=
⎧ ⎪ ⎨
κ κ
x y
⎩⎪κ xy
⎫ ⎪ ⎬ ⎪ ⎭
=
⎧ ⎪ ⎪ ⎪⎪ ⎨ ⎪ ⎪ ⎪−
− ∂θx ∂x
− ∂θ y ∂y
( ) 即 ΔFz1
=
−ΔFz 4
=
h l14
F14 − F41
（5-4）
对所有的边都进行变换后，即可得到[R]T 。
6、 单元刚度矩阵的奇异性问题
在 Quad4 单元中认为壳体法线的旋转自由度θz = 0 。对于平的壳体单元就意味着每个
节点的刚度矩阵只耦合 5 个自由度，这将导致求解时刚度矩阵出现奇异性问题。除非法向旋 转连接着其他单元（如 Cbar）或相临单元不位于同一平面内，必须对单元刚度矩阵进行处
0
λ 0 0
0 λ 0
0
⎥ ⎥
,
0⎥
λ
⎥ ⎦
λ
=
⎡φ ⎢⎣0
0⎤ φ ⎥⎦
单元节点力变换公式为 F e = T F e
故整体坐标系下单元刚度矩阵为 K e = T T K eT
（4-1）
（4-2） （4-3） （4-4）
5、 单元 G1、G2、G3、G4 节点不共面问题 Quad4 单元并不要求四个节点共面，当采用平板壳单元模拟曲面时，四个节点可能不在
同一平面上，如图 2 所示。为了定义单元 xy 平面，首先定义 2h 为 G1G3 连线与 G2G4 连线 的距离，然后 G1、G3 下移 h，G2、G4 上移 h，则新形成的四个点将构成单元 xy 平面。