壳单元的公式推导

( ) 即 ΔFz1
=
−ΔFz 4
=
h l14
F14 − F41
（5-4）
对所有的边都进行变换后，即可得到[R]T 。
6、 单元刚度矩阵的奇异性问题
在 Quad4 单元中认为壳体法线的旋转自由度θz = 0 。对于平的壳体单元就意味着每个
节点的刚度矩阵只耦合 5 个自由度，这将导致求解时刚度矩阵出现奇异性问题。除非法向旋 转连接着其他单元（如 Cbar）或相临单元不位于同一平面内，必须对单元刚度矩阵进行处
[ ] { } 在节点 i 处
Fi
T plane
=
Fxi
Fyi
Fzi
M xi
M yi
0 ，变化到节点后将产生附加的力
矩 ΔM xi = ±hFyi 和 ΔM yi = ±hFxi ，正负取决于节点在 xy 平面的上方还是下方。
由于单元的膜拉伸刚度远大于弯曲刚度，因此弯矩引起的弯曲变形远大于拉力引起的 拉伸变形，所以常用拉伸力取代弯矩。由图 2 可见，在 G1 和 G4 上的修正力是一个垂直力，
中增加节点角位移θz ，即单元平面内的转角或法线绕 z 轴的转角。即
{ } { } Fi p = Fxip
Fypi
T
M zi ,
aip = uip
vip
θT zi
（3-3）
这样单元刚度矩阵被扩大为 24 阶方阵，其中与θz 对应的行和列皆为零元素。
所以，平板壳单元的基本方程为
F
e
=
⎧ ⎪⎪ ⎨
F1 F2
其中单元刚度矩阵的字块为 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，并符合右手法则。
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,
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 ⎪⎭
M xi
T
M yi ,
aib = wi
θ xi
θ yi
T
，单元刚度矩阵同
Ke bending
。
将上述两部分单元基本方程组合起来就可以得到壳元的基本方程和刚度矩阵。对于平板型矩 形壳元，每个节点有 5 个自由度，其刚度矩阵为 20 阶方阵。然而在整体坐标系中，壳体单 元每个节点有 6 个自由度，为了便于从局部坐标系到整体坐标系的转换，在单元局部坐标系
∂θx − ∂θ y
⎫ ⎪ ⎪ ⎬⎪⎪ , ⎪ ⎪ ⎪
⎪⎩ ∂y ∂x ⎪⎭
γ
=
⎧γ ⎨⎩γ
xz xz
⎫ ⎬ ⎭
=
⎧ ⎪⎪
∂w ∂x
−
θ
x
⎫ ⎪⎪
⎨ ⎪ ⎩⎪
∂w ∂y
−
θ
y
⎬ ⎪ ⎪⎭
（2-3）
⎧εx ⎫
代入几何方程可得 ε
=
⎪⎪⎪⎨γεxyy
⎪ ⎪⎪ ⎬
⎪⎪γ
xz
⎪ ⎪
=
⎧zκ ⎫
⎨ ⎩
γ
⎬, ⎭
理论为出发点，所以在壳厚度逐渐变薄时将会出现剪切闭锁现象，主要原因就是在厚度很薄 时，转角为挠度的导数而不再独立，采用分别插值就会出现问题。
为了避免剪切闭锁现象，可以采用减缩积分方案。但是，这种方法并不总是有效，因为 减缩积分可能造成零能模式，即采用非刚体位移模式时系统的变形能为零（实际是不可能 的）。因此，对于 Quad4 单元的剪切应变场需要进行特别处理。
yi
⎫ ⎪⎪
∑ ∂Ni
∂η
⎬
yi
⎪ ⎪⎭
(i = 1, 2,3, 4) （1-4）
代入本构方程可得σ = Dε = DBae = Sae = [S1 S2 S3 ]S4 ae
（1-5）
⎧
⎫
⎪1 ν 0 ⎪
其中
D
=
E 1−ν
2
⎪⎨ν ⎪
⎪0
1 0
⎪ 0⎬
1−ν
⎪ ⎪
⎩
2⎭
∫ ∫ ∫ ∫ 单元刚度矩阵
=
1 4
(1+ ξ )(1+η)
，N4
=
1 4
(1− ξ )(1+η)
（0-1）
1、 平面应力单元
位移函数为
u
=
⎧u ⎨⎩v
⎫ ⎬ ⎭
=
[
N1I
N2I
N3I
N4I ] ae = Nae
其中 I 为二阶单位矩阵；单元节点位移向量为
⎧a1 ⎫
ae
=
⎪⎪⎨⎪aa23
⎪⎪ ⎬ ⎪
,
⎪⎩a4 ⎪⎭
a
i=
⎨⎧⎩uvii
⎫ ⎬ ⎭
(i = 1, 2,3, 4)
（1-1） （1-2）
[ ] 代入几何方程可得 ε = B1 B2 B3 B4 ae = Bae
（1-3）
⎧ ⎪
∂N
i
⎪ ∂x
⎪ 其中 Bi = ⎨ 0
⎪
⎪ ⎪ ⎩
∂Ni ∂y
⎫
0⎪
⎪
∂Ni ∂y
⎪ ⎬ ⎪
∂Ni ∂x
⎪ ⎪ ⎭
(i = 1, 2,3, 4)
根据复合函数求导法则，可得
⎪⎩γ 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⎥ , ⎥ ⎥
⎡ ⎢
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
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
∂ξ
−
Ni
Bsi
=
⎢ ⎢
⎣⎢
0
0⎥
⎦
0 −Ni
∂Ni ⎤
∂x ∂Ni
⎥ ⎥ ⎥
,
∂y ⎦⎥
(i = 1, 2,3, 4)
根据复合函数求导法则，可得
（2-4）
⎧ ∂Ni
⎪⎪ ∂x
⎨ ⎪
∂Ni
⎫ ⎪⎪ ⎬ ⎪
=
J
−1
⎧ ∂Ni
⎪⎪ ∂ξ
⎨ ⎪
∂Ni
⎫ ⎪⎬⎪ , ⎪
⎩⎪ ∂y ⎭⎪ ⎩⎪ ∂η ⎭⎪
⎧ ∂x
[ ] [ ] ⎡⎣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）
=
K
a ep ep
（3-1）
{ } { } 其中 Fi p = Fxip
Fypi
T
,
aip = uip
vip
T
，单元刚度矩阵同
Ke membrane
。
在单元局部坐标系下，与弯曲变形有关的情况同厚板弯曲问题，单元基本方程可写为
F eb
=
⎧ ⎪⎪ ⎨ ⎪
F1b F2b F3b
⎫ ⎪⎪ ⎬ ⎪
⎪⎩F4b ⎪⎭
J
=
⎪⎪ ⎨ ⎪
∂ξ ∂x
⎩⎪ ∂η
∂y ⎫ ⎧
∑ ∂ξ
∂x
⎪⎪ ⎬ ⎪
=
⎪⎪ ⎨ ⎪
∑ ∂η ⎭⎪ ⎩⎪
∂Ni ∂ξ
xi
∂Ni ∂η
xi
∑ ∂Ni
∂ξ
yi
⎫ ⎪⎪
∑ ∂Ni
∂η
⎬
yi
⎪ ⎭⎪
(i = 1, 2,3, 4) （2-5）
代入本构方程可得σ = Dε
（2-6）
⎡1 ν 0 0 0 ⎤
⎢⎢ν 1 0
=
⎡⎢⎢kk12bb11 ⎢⎢⎢⎣kk34bb11
k1b2 k2b2 k3b2 k4b2
k1b3 k2b3 k3b3 k4b3
k1b4 k2b4 k3b4 k4b4
⎤ ⎥ ⎥ ⎥ ⎥ ⎥⎦
⎧⎪⎪⎨⎪aaa132bbb ⎪⎩a4b
⎫ ⎪⎪ ⎬ ⎪ ⎪⎭
=
K
a eb eb
（3-2）
{ } { } 其中 Fib = Fzi
F ep
=
⎧⎪⎪⎨⎪FFF132ppp
⎫ ⎪⎪ ⎬ ⎪
⎪⎩F4p ⎪⎭
=
⎡⎢⎢kk12pp11 ⎢⎢⎢⎣kk34pp11
k1p2 k2p2 k3p2 k4p2
k1p3 k2p3 k3p3 k4p3
k1p4 k2p4 k3p4 k4p4
⎤ ⎥ ⎥ ⎥ ⎥ ⎥⎦
⎧⎪⎪⎨⎪aaa132ppp ⎪⎩a4p
⎫ ⎪⎪ ⎬ ⎪ ⎪⎭
同一平面上，如图 2 所示。为了定义单元 xy 平面，首先定义 2h 为 G1G3 连线与 G2G4 连线 的距离，然后 G1、G3 下移 h，G2、G4 上移 h，则新形成的四个点将构成单元 xy 平面。
图 2 单元四个节点不共面
在有翘曲的情况下，单元将承受不平衡的力和力矩，必须 xy 面变换到实际的几何节点 位置上，以计算刚度。节点的单元刚度阵与平面内的单元刚度阵变换关系为：
x, z y, z z, z
⎤ ⎥ ⎥ ⎥ ⎥ ⎥⎦
⎧ ⎪⎪ ⎨ ⎪ ⎩⎪
x y z
⎫ ⎪⎪ ⎬ ⎪ ⎭⎪
=
φ
⎧ ⎪⎪ ⎨ ⎪ ⎩⎪
x⎫
y ⎪⎬⎪
z
⎪ ⎪⎭
其中φ 为局部坐标系 xyz 与整体坐标系 x y z 之间的方向余弦矩阵。
⎡λ 0 0 0⎤
节点位移变换公式为 ae
=
T ae
=
⎢ ⎢
0
⎢0
⎢ ⎣
(i = 1, 2,3, 4)
（2-1） （2-2）
⎧εx ⎫
对厚板 ε
=
⎪⎨⎪⎪γεxyy ⎪⎪γ xz
⎪ ⎪⎪ ⎬ ⎪ ⎪
=
⎧ zκ
⎨ ⎩
γ
⎫ ⎬, ⎭
⎪⎩γ yz ⎪⎭
κ
=
⎧ ⎪ ⎨
κ κ
x y
⎩⎪κ xy
⎫ ⎪ ⎬ ⎪ ⎭
=
⎧ ⎪ ⎪ ⎪⎪ ⎨ ⎪ ⎪ ⎪−
− ∂θx ∂x
− ∂θ y ∂y
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 单元并不要求四个节点共面，当采用平板壳单元模拟曲面时，四个节点可能不在
⎪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⎥⎥
理。其中一种方法就是在节点的单元刚度矩阵对应位置给以任意的虚拟刚度系数 Kθzθz ，这
样局部坐标系中θz 方向的平衡方程为 Kθzθzθz = 0 。经过坐标变化，整体坐标系中的该节点
平衡方程将满足有唯一解条件。由于θz 并不影响应力，而且与其他节点平衡方程无关，故
Kθzθz 可赋予任何非零值。
7、 壳单元的通用性 根据壳单元的厚度与面内尺寸的比例关系，壳可以分为薄壳和厚壳。Quad4 单元以厚壳
壳单元的公式推导
Nastran 中的 Quad4 单元是由平面应力单元和厚板弯曲单元组成，单元包括 4 个节点，每个
节点有 5 个自由度（ u, v, w,θx ,θ y ），没有绕单元法线的旋转自由度θz 。
单元形函数：
N1
=
1 4
(1− ξ )(1−η) ，N2
=
1 4
(1+ ξ )(1−η)
，N3
图 1 单元局部坐标系 局部坐标与整体坐标之间的关系为
( ) ( ) ( ) ⎧x⎫
⎡cos x, x ⎢
⎪
( ) ( ) ( ) ⎨ ( ) ( ) ( ) ⎪⎩
y z
⎪ ⎬ ⎪⎭
=
⎢⎢cos ⎢⎣⎢ cos
y, x z, x
cos x, y cos y, y cos z, y
cos cos cos