用于光滑/非光滑壳的稳定化新型协同转动4 节点四边形壳单元

李忠学，胡万波

(浙江大学土木工程系，杭州 310058)

在单元中引入协同转动法，可采用线性单元理论求解几何非线性问题，而且协同转动框架不依赖于特定的某个单元，因此可以利用现有的、性能较好的线性单元作为协同转动单元公式的“核”，建立新的单元公式，使一些原本只能用于解决小位移与小转角问题的单元公式可用于解决大位移、大转角问题[1−4]。

在进行壳单元有限元分析时，往往会遇到闭锁现象[5−6]。闭锁现象的实质是单元的假设位移模式无法准确地描述实际位移模式，导致单元刚度偏硬，位移的计算值偏小[7]。在壳单元有限元计算中，剪切闭锁和膜闭锁的影响较为严重[8]。对于以受弯为主的单元，当其膜应力很小或接近于零时，采用数值积分将会导致膜闭锁问题，使单元的膜刚度被夸大而使求解失真。对于薄壁梁元或薄板单元、薄壳单元，在采用数值积分计算单元刚度矩阵时，常会由于剪切闭锁而导致单元剪切刚度被夸大，从而使得到的结构变形偏小[9]。

解决闭锁问题常采用2 类方法：一类是通过改变积分方案来消除闭锁，如缩减积分法；另一类是改变应变能表达式，如各种假定应变法[10−14]。缩减积分法有简单、精度高、计算效率高等优点，但由于使用的积分点少于高斯积分所需的积分点数，在求解过程中，单元刚度矩阵的秩有时会出现缺陷，导致无法捕捉到某些变形模式，从而导致求解失败。

自Zienkiewicz 等[15]首次提出缩减积分法以来，为了控制因缩减积分引起的零能模态，使单 元 更 稳 定，Kosloff[16]、Taylor[17]、Flanagan[18]、Liu[19−20]、王丽等[21]相继提出了各种改进方法。MacNeal[22]基于释放多余约束的等参原理，发展了QUAD4 单元。Hughes 等[23]在MacNeal[22]启发下，基于缩减积分法，横向变形采用9 节点双二次幂形函数进行插值，转动变量采用4 节点双线性形函数插值，构造出一种没有零能模态的新型4 节点四边形双线性等参单元。Belytschko 和Bindeman[24]基于稳定性刚度与不完全积分刚度两者最大特征值成比例的思路，提出了稳定性参数的取值方法，选用合理的参数来控制零能模态，又不至于引起闭锁现象。Belytschko 和Tsay[25]通过增广B 矩阵并设置稳定性参数实现稳定化过程，提出了扰动零能模态控制法。然而，扰动法也无法完全消除潜在的零能模态，且计算结果容易受稳定性参数的影响。Belytschko 和Leviathan[26]基于Hu-Washizu 变分原理，在4 节点四边形单元中，利用混合张量插值法来消除剪切闭锁，通过广义应变来控制零能模态，提出了采用物理稳定化零能模态控制法的QPH 单元。李忠学等[27]采用假定应变法构造出稳定应变来消除零能模态，发展了一种适用于光滑壳的新型协同转动9 节点四边形单元。

本文发展了一种新型协同转动4 节点四边形壳单元，采用矢量型转动变量描述旋转，这些变量在非线性增量求解过程中可以采用简单的加法进行更新，在整体坐标系和局部坐标系下都能得到对称的单元切线刚度矩阵，可提高计算效率和节省存储空间。为消除单元中可能出现的闭锁现象，本文采用了单点积分方案计算内力矢量和单元切线刚度矩阵，并引入Belytschko 等[26]的物理稳定化零能模态控制法，以避免因切线刚度矩阵缺秩而产生的零能模态。新型单元在4 种典型算例中均表现出良好的性能。

1 单元的新型协同转动法描述

单元协同转动框架如图1 所示。初始构形下，四边形两条对角线的方向矢量v130和v240的定义如下：

图1 协同转动框架Fig. 1 Description of the co-rotational framework

2 局部坐标系下的单元公式

将应变 ε分解为3 部分，即：膜应变 εm、弯曲应变zlχ、出平面的剪切应变 γ。各应变按下式计算：

3 稳定化的单点积分法

图2 剪切应变插值关联点分布图Fig. 2 Interpolation points for shear strain

通过混合张量插值应变法[28]得到自然坐标系下的单元剪切应变：

4 整体坐标系下的单元公式

5 算例分析

将引入物理稳定化零能模态控制法的4 节点四边形壳单元简记为PO4Q。

5.1 径向受拉的开口圆柱壳

如图3 所示的开口圆柱壳，两端自由，几何尺寸为：长度L=10.35 m ，半径R=4.953 m，厚度a=0.094 m ，材料的弹性模量E=10.5 MPa，泊松比μ=0.3125。圆柱壳在中部受一对大小相等、方向相反的径向集中拉力F=40 kN作用。利用结构和荷载的对称性，取壳的1/8 结构进行分析。

图3 径向受拉圆柱壳Fig. 3 Pull-out of an open cylindrical shell

计算荷载作用点的挠度时，分别采用8×16、12×24 的单元网格，得到结构的荷载-位移曲线与Jiang 和 Chernuka[30]采用8×14 单元网格的4 节点ANS 退化壳单元以及Campello 等[31]采用8×16×2单元网格的6 节点三角形曲壳单元的计算结果进行比较，如图4 所示。可见，当荷载达到 20 kN左右时，壳体出现了轻微的跳跃屈曲现象，该单元成功跟踪到了这一过程。得到的临界荷载和荷载-位移曲线与其他文献[30−31]中的结果吻合较好。

图4 径向受拉圆柱壳在加载点处的荷载-位移曲线Fig. 4 Load-displacement curve of cylindrical shell

5.2 悬臂梯形板

图5 所示的悬臂梯形板，一端固支，另一端受到平面内均布剪切荷载。结构几何尺寸如下：长L=48 mm ，高H1=16 mm、H2=44 mm。材 料的弹性模量E=1.0 MPa ，泊松比μ=1/3，厚度a=1.0 mm 。均布荷载p=P/H1，其中，P=λPref，Pref=1 N。

图5 悬臂梯形板Fig. 5 A cantilever trapezoidal plate

图6 给出了分别采用8×8、16×16 和32×32 个PO4Q 单元和ANSYS[32]有限元软件采用8×8、16×16 和32×32 个shell181 单元的计算结果。表1给出了荷载参数为1 时采用16×16、32×32 个PO4Q单元和shell181 单元的计算结果。可见，PO4Q 单元略微偏硬，但采用32×32 个PO4Q 单元与shell181单元的计算结果的相对误差仅5.40%，PO4Q 单元的收敛性和计算精度基本令人满意。

图6 悬臂板在点A 处的荷载-位移曲线Fig. 6 Load-displacement curves at tip A of cantilever plate

表1 悬臂板在点A 处的位移(λ=1)Table 1 Displacement at tip A of cantilever plate (λ=1)

5.3 直角形折板

图7 所示的直角形折板，一端完全约束，另一端自由，并施加竖向均布荷载。结构几何尺寸为：长L=6 m ，宽B=3 m ，厚度a=0.03 m。材料的弹性模量E=30 MPa ，泊松比μ=0.0。自由端受均布荷载p=P/B，其中，P=λPref，Pref=3 N。

图7 直角形折板Fig. 7 A right-angle cantilever subjected to distributed force

为检验单元的收敛性，分别采用3×6×2 和6×12×2 的PO4Q 单元网格对该折板进行分析，自由端的荷载-位移曲线在图8 中给出。作为对比，还给出了采用ANSYS[32]软件中的shell181、shell281以及Li 等[33]采用QUAD9 单元的计算结果。可见，采用PO4Q 单元网格得到的结果与ANSYS[32]和Li 等[33]的计算结果吻合较好。

图8 直角形折板在自由端处的荷载-位移曲线Fig. 8 Load-displacement curves of right-angle cantilever

5.4 镰刀形壳

如图9 所示的镰刀形壳，结构左端固支，在右侧自由端的中点A 处作用一平面内集中荷载。结构的几何尺寸为：L=10 m，R=5 m ，宽B=1 m，厚度a=0.01 m 。材料的弹性模量E=30 MPa，泊松比μ=0.3 。自由端作用的集中荷载值为：P=λPref，Pref=0.01 N。

图9 镰刀形壳Fig. 9 A cantilever sickle shell subjected to a lateral tip force

分别采用(15+15)×3 和(20+20)×4 的单元网格计算该镰刀形壳自由边中点的位移。其结果和Chróścielewski 等[34]采用基于合应力的半混合壳单元SEMe9、基于位移-转动的壳单元CAMe9 以及标准退化壳单元SELe9 得到的结果吻合较好(见图10)。

图10 自由边中点的荷载-位移曲线Fig. 10 Load-displacement curves of cantilever sickle shell

该镰刀形壳在不同荷载作用下的变形如图11所示。

图11 不同荷载作用下镰刀形壳变形图Fig. 11 Deformed shape of cantilever sickle shell under different levels of tip force

6 结论

本文发展了一种适用于解决光滑和非光滑板壳大位移和大转角弹性变形问题的新型协同转动4 节点四边形壳单元。单元中采用矢量型转动变量描述旋转，采用单点积分法消除闭锁现象，引入物理稳定化方法消除单点积分引起的零能模态。通过对2 个光滑壳和2 个非光滑壳的非线性分析，并将结果与参考文献中的单元或ANSYS 软件中的单元的计算结果进行对比，表明该单元已有效地减轻了闭锁现象，并消除了零能模态的不利影响，单元的可靠性、计算效率与计算精度是令人满意的。