有限元板壳单元

在板理论中经常用内力，即弯矩和剪力来表示， 它们与应力之间的关系为：
Mx = ∫ Qx = ∫
hБайду номын сангаас/2
− h /2
σ x zdz M y = ∫
Qy = ∫
h /2
− h /2
σ y zdz M xy = ∫
h /2
− h /2
τ xy zdz
h /2
− h /2
τ xz dz
h /2
− h /2
根据Mindlin板理论的假设，板内任意一点的位移由3个 广义位移
w,ψ x ,ψ y
确定，为了与有限元的节点位移
T
相对应，采用的位移列阵为：
δ = {w1 θ x1 θ y1 w2 θ x 2 θ y 2 w3 θ x 3 θ y 3 w4 θ x 4 θ y 4 }
e
θ x = ψ y , θ y = −ψ
把位移模式代入得到：
⎡ zBb ⎤ e ε = ⎢ ⎥δ ⎣ Bs ⎦ 式中，Bb = [ Bb1 Bb 2 " Bbn ] , Bs = [ Bs1 Bs 2 " Bsn ⎡ 0 ⎢0 ⎢ ⎢ ∂N i Bbi = ⎢0 − ∂y ⎢ ⎢ ∂N i ⎢0 − ∂x ⎣ ∂N i ⎤ ⎥ ∂x ⎥ ⎥ 0 ⎥ ⎥ ∂N i ⎥ ⎥ ∂y ⎦ ⎡ ∂N i ⎢ ∂y Bsi = ⎢ ⎢ ∂N i ⎢ ⎣ ∂x
2 ⎡ ⎤ a 2 a22 = Hb ⎢ 2 (1 − μ ) ξ 0 ( 3 + 5η0 ) + 5 2 ( 3 + ξ 0 )( 3 + η0 ) ⎥ b ⎣ ⎦
a23 = −15 H μ ab (ξi + ξ j )(ηi + η j ) ⎡⎛ ⎤ b2 ⎞ b2 a31 = 3Ha ⎢⎜ 2 + 3μ + 5 2 ⎟ ξ jη0 + 15 2 ξ j + 5μξiη0 ⎥ a ⎠ a ⎣⎝ ⎦ a32 = −15 H μ ab (ξi + ξ j )(ηi + η j )
7.1.1 Kirchhoff板理论 由于忽略横向剪切变形，即 w' x = ψ x 和w' y = ψ y ， 因此板内所有的力学量都能用挠度w表示：
⎧M x ⎫ ⎡1 μ ⎪ ⎪ ⎢ ⎨M y ⎬ = − D ⎢ μ 1 ⎪ ⎪ ⎢0 0 M ⎣ ⎩ xy ⎭
3
⎤ ⎧ w' xx ⎫ ⎪ ⎥⎪ ⎥ ⎨ w' yy ⎬ ⎪ ⎪ ⎥ − 1 μ 2 ( ) ⎦ ⎩2w' xy ⎭ 0 0
2 ⎡ ⎤ b 2 a33 = Ha ⎢ 2 (1 − μ )η0 ( 3 + 5ξ 0 ) + 5 2 ( 3 + ξ 0 )( 3 + η0 ) ⎥ a ⎣ ⎦ D , ξ 0 = ξiξ j ,η0 = ηiη j 式中 H = 60ab
7.3 基于Mindlin板理论的四边形单元
基于Kirchhoff 薄板理论的薄板矩形单元忽略了 剪切变形的影响。由于Kirchhoff 板理论要求挠 度的导数连续，给构造协调单元带来了不少麻 烦。为此，采用考虑剪切变形的Mindlin 板理论 来克服。这种方法比较简单，精度较好，并且 能利用等参变换，得到任意四边形甚至曲边四 边形单元，因而实用价值较高。
由节点位移条件可求得待定系数a1至a12：再代入 位移模式整理后得到：
w = ∑ ( N i wi + N xiθ xi + N yiθ yi ) = N δ e
i =1 4
其中形函数： N i = (1 + ξiξ )(1 + ηiη ) ( 2 + ξiξ + ηiη − ξ 2 − η 2 ) / 8 N xi = −bηi (1 + ξiξ )(1 + ηiη ) (1 − η 2 ) / 8 N yi = aξi (1 + ξiξ )(1 + ηiη ) (1 − ξ 2 ) / 8
式中，抗弯刚度D为 Eh D= 12 (1 − μ 2 )
7.1.2 Mindlin板理论
根据Mindlin板理论的假设，中面法线在变性后不 再垂直于中面，因此必须采用3个位移分量
w、ψ x 和ψ y 来描述板内的变形，即
⎧ε x ⎫ ⎧ψ x ' x ⎫ ⎪ ⎪ ⎪ ⎪ ⎨ε y ⎬ = − z ⎨ψ y ' y ⎬ ⎪ ⎪ ⎪ ⎪ ⎩γ xy ⎭ ⎩ψ x ' y +ψ y ' x ⎭
τ yz dz
式中M x、M y 是弯矩，M xy 是扭矩，Qx、Qy 是剪力
对于线弹性材料，板内的应力应变关系为：
⎧σ x ⎫ ⎪ ⎪ E ⎨σ y ⎬ = 2 − μ 1 ⎪ ⎪ ⎩τ xy ⎭
⎡1 μ ⎢ ⎢μ 1 ⎢0 0 ⎣
⎧ε x ⎫ ⎤ ⎧ε x ⎫ ⎪ ⎪ ⎥⎪ ⎪ 0 ⎥ ⎨ε y ⎬ = D ⎨ε y ⎬ ⎪ ⎪ ⎪ ⎪ (1 − μ ) 2⎥ γ ⎦ ⎩ xy ⎭ ⎩γ xy ⎭ 0
可以得到转角：
∂w ∂w 1 = = (a3 + a3ξ + 2a6η + a8ξ 2 θx = ∂y b∂η b + 2a9ξη + 3a10η 2 + a11ξ 3 + 3a12ξη 2 )
∂w ∂w 1 θy = = = − (a2 + 2a4ξ + a5η + 3a7ξ 2 ∂x a∂ξ a + 2a8ξη + a9η 2 + 3a11ξ 2η + a12η 3 )
（3）单元应变力的表达
由物理方程有：
σ = [ D] ⋅ ε = S ⋅ δ
e
e
⎡1 E ⎢ μ [ D] = 2 ⎢ 1− μ ⎢ ⎣0
⎤ ⎥ 1 ⎥ 0 1 − μ 2⎥ ⎦ 0 0
μ
（4）单元刚度矩阵
单元刚度矩阵的分块子矩阵的公式：
3 1 1 h 2 T K ij = ∫∫∫ z Bi DBi dxdydz = ∫ ∫ BiT DBi abd ξ dη (i, j = 1, 2,3, 4) 12 −1 −1
x
（1）单元位移场的表达
w = ∑ ( N i wi + N xiθ xi + N yiθ yi ) = N ⋅ δ e
i =1
4
其中形函数： N = [ N1 N 2 " N N
]
（2）单元应变场的表达
Mindlin 板理论考虑了横向剪切变形，因此应变有5 个分量，即
z ∂θ y ∂x ⎫ ⎧ε x ⎫ ⎧ ⎪ ⎪ ⎪ ⎪ − z ∂θ x ∂y ⎪ ⎪ε y ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ε = ⎨γ xy ⎬ = ⎨ z ( ∂θ y ∂y − ∂θ x ∂x ) ⎬ ⎪ ⎪ ⎪ ⎪ ∂w ∂y − θ x ⎪γ yz ⎪ ⎪ ⎪ ⎪γ ⎪ ⎪ ⎪ θ w x ∂ ∂ + y ⎩ zx ⎭ ⎩ ⎭
7.2 矩形薄板单元
每一个节点有三个位移分量，即挠度w， x θ 和y θ 。 则单元位移列阵和节点力列阵分别为：
δ = {w1 θ x1 θ y1 w2 θ x 2 θ y 2 w3 θ x 3 θ y 3 w4 θ x 4 θ y 4 }
e
e
T
P =⎡ ⎣ f z1 M θ x1 M θ y1 f z 2 M θ x 2 M θ y 2 f z 3 M θ x 3 M θ y 3 f z 4 M θ x 4 M θ y 4 ⎤ ⎦
把应变矩阵B和弹性矩阵D代入并运算得到：
⎡ a11 a12 a13 ⎤ ⎢ ⎥ K ij = ⎢ a21 a22 a23 ⎥ ⎢ ⎥ a a a ⎣ 31 32 33 ⎦
式中9个元素为：
⎡ ⎛ b2 ⎤ a2 ⎞ ⎛ b2 a2 ⎞ a11 = 3H ⎢15 ⎜ 2 ξ 0 + 2 η0 ⎟ + ⎜14 − 4 μ + 5 2 + 5 2 ⎟ ξ 0η0 ⎥ b a b ⎠ ⎠ ⎝ ⎣ ⎝a ⎦ ⎡⎛ ⎤ a2 ⎞ a2 a12 = −3Hb ⎢⎜ 2 + 3μ + 5 2 ⎟ ξ 0ηi + 15 2 ηi + 5μξ 0η j ⎥ b ⎠ b ⎣⎝ ⎦ ⎡⎛ ⎤ b2 ⎞ b2 a13 = 3Ha ⎢⎜ 2 + 3μ + 5 2 ⎟ ξiη0 + 15 2 ξi + 5μξ jη0 ⎥ a ⎠ a ⎣⎝ ⎦ ⎡⎛ ⎤ a2 ⎞ a2 a21 = −3Hb ⎢⎜ 2 + 3μ + 5 2 ⎟ ξ 0η j + 15 2 η j + 5μξ 0ηi ⎥ b ⎠ b ⎣⎝ ⎦
T
（1）单元位移场的表达
与平面问题矩形单元类似，引入自然坐标系（ξ,η）。由于每个节 点有3个位移分量，所以选取含有12个参数的多项式作为位移模 式：
w = a1 + a2ξ + a3η + a4ξ 2 + a5ξη + a6η 2 + a7ξ 3 + a8ξ 2η + a9ξη 2 + a10η 3 + a11ξ 3η + a12ξη 3
u = − zψ x
v = − zψ y ε x = − zψ x ' x ε y = − zψ y ' y
γ xy = − z (ψ x ' y +ψ y ' x ) γ yz = w' y −ψ y γ xz = w' x −ψ x
这里w是板的横向挠度，假设它沿板的厚度方向不变，即εz =0 。上式是Mindlin 板理论的基本假定。如果假定变形后的 法线仍然是变形后中面的法线，即w,x =ψ x和w,y =ψ y ，则 式中的两个横向剪切应变γyz 和γxz 为零，这就退化为 Kirchhoff 板理论。当板足够薄时，用Kirchhoff 板理论能得到 符合实际的结果。
Qx = κ Gh ( w' x −ψ x )
Qy = κ Gh ( w' y −ψ y )
应用有限元法求解板弯曲问题时，用一些 离散的板单元代替原来连续的结构。每个 结点有三个广义位移分量，即挠度w，绕x 轴的转角xθ 和绕y 轴的转角y θ 。挠度w的 正方向跟z 轴一致，转角则以按右手螺旋法 则标出的矢量沿坐标轴正向为正。
则内力与位移的关系为：
⎧M x ⎫ ⎡1 μ ⎪ ⎪ ⎢ ⎨M y ⎬ = − D ⎢ μ 1 ⎪ ⎪ ⎢0 0 M ⎣ ⎩ xy ⎭ ⎫ ⎤ ⎧ψ x ' x ⎪ ⎥⎪ ⎬ ⎥ ⎨ψ y ' y ⎥ ⎪ψ +ψ ⎪ (1 − μ ) 2⎦ y'x ⎭ ⎩ x' y 0 0
另外为了修正横向剪应力沿板厚均匀分布导致的误差， 引入了所谓剪切修正因子κ 来修正剪力，即
（4）单元刚度矩阵
]
− Ni 0
⎤ 0 ⎥ ⎥ (i = 1, 2,...n) ⎥ Ni ⎥ ⎦
（3）单元应力场的表达
相应的应力也有5 个分量，它们与应变的关系是：
⎧σ x ⎫ ⎧ε x ⎫ ⎪ ⎪ ⎪ ⎪ ⎪σ y ⎪ ⎪ε y ⎪ ⎤⎪ ⎪ ⎪ ⎪ ⎡ Db σ = ⎨τ xy ⎬ = ⎢ ⎨γ xy ⎬ ⎥ Ds ⎦ ⎪ ⎪ ⎪ ⎪ ⎣ τ γ yz ⎪ ⎪ ⎪ yz ⎪ ⎪τ ⎪ ⎪γ ⎪ ⎩ zx ⎭ ⎩ zx ⎭ ⎡ Db ⎤ ⎡ zBb ⎤ e =⎢ ⋅δ ⎥ ⎢ ⎥ Ds ⎦ ⎣ Bs ⎦ ⎣
（2）单元应变场的表达
由弹性力学几何方程有：
⎧ε x ⎫ ⎧ w' xx ⎫ ⎪ ⎪ ⎪ ⎪ e ⎨ε y ⎬ = − z ⎨ w' yy ⎬ = z [ B1 B2 B3 B4 ] δ ⎪ ⎪ ⎪ ⎪ γ w 2 ⎩ xy ⎭ ⎩ ' xy ⎭ 式中 ⎧ N i 'ξξ / a 2 ⎫ ⎧ N i ' xx ⎫ ⎧bN i 'ξξ / a ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ 1 ⎪ ⎪ 2 ⎪ Bi = − ⎨ N i ' yy ⎬ = − ⎨ N i 'ηη / b ⎬ = − ⎨aN i 'ηη / b ⎬ (i = 1, 2,3, 4) ab ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ N ab N 2 / 2 N 2 i 'ξη ⎭ ⎪ i 'ξη ⎪ ⎩ ⎩ i ' xy ⎭ ⎩ ⎭
第七章
板壳单元
对于小挠度弹性薄板弯曲问题，板的变形完 全由垂直于板面的挠度w确定，在一般情况 下，取w和它的一阶、二阶导数为参数进行数 值计算。当前用来离散薄板的单元多用四边 形或三角形单元，相邻之间有弯矩传递，所 以将节点看成刚性的。
7.1 弹性板的弯曲
设板的中面在xy 平面上，即z = 0表示板的中面，在 板理论中，一般假设板的中面是一中性面，也就是在 没有面内力时，中面上的三个应变 ε x =ε y =γ xy = 0 。另一个基本假设即为所谓的直法线假定：变形前 垂直于中面的法线变形后仍然保持直线，但是不一定 仍然垂直于变形后的中面。这条直线有绕y 和x 轴的 转角分别为ψ x 和ψ y 。则距离中面距离为z 的任意 点的位移和应变分别是