在OpenSees中加入新型几何非线性三角形壳单元

OpenSees软件因其开源的特性，与众多商业有限元软件相比具有较强的扩展性和灵活性，在工程结构的动力灾变模拟中受到了越来越多的关注。在OpenSees中，常用的壳单元有MITC4和DKGQ/NLDKGQ单元两种。在OpenSees中集成的MITC4单元是一种三角形/四边形通用单元，但其目前无法考虑几何非线性行为，同时在计算中可能存在闭锁现象。此外MITC4单元只有5个自由度，缺少面外转角自由度，使得该单元和梁柱等单元连接时可能出现问题，给建模和计算带来了不便。DKGQ单元是一种四边形平板壳单元，NLDKGQ是其拓展形式，能够模拟几何非线性，可以应用于大变形问题。但该单元在模拟结构屈曲时容易出现计算不收敛的问题。针对OpenSees在壳单元方面功能的不足，有必要提出一种可以应用于几何非线性问题的三角形壳单元。

清华大学土木工程系研究生张书豪、田源等在之前工作的基础上，基于广义协调理论的三角形平面膜单元GT9以及三角形薄板弯曲单元DKT来构造新型三角形壳元DKGT。平面膜元GT9采用单元间协调条件，通过构造节点刚体转角，使其位移场具有更高阶次，提高了三角形单元的计算精度。同时使GT9的每个节点拥有三个自由度。DKT单元是基于Kirchhoff薄板弯曲理论构造，可以有效避免剪切闭锁的发生。 

图1 三角形平板壳单元 DKGT单元的基本构成

算例1  Scordelis-Lo屋盖算例

图2所示是一个受竖向均布荷载作用的单曲壳体，均布荷载g=90。壳体两端用刚体隔板支撑，考察其自由边中点的竖向位移，本文以文献（MacNeal et al, 1985）中提供的理论解作为计算结果的比较值。利用结构的对称性，取图中所示壳体的四分之一进行有限元计算，并采用不同的网格划分，得到结果如表1所示。同时表1还给出了DKT-CST-15RB单元（Nicholas et al, 1986）和OLSON单元（Olson et al, 1979）在相同网格划分下的计算结果，结果证明，DKGT单元具有较好的收敛性和计算精度。

图2 Scordelis-Lo屋盖问题

表1   Scordelis-Lo屋盖问题的自由边中点位移

参考文献

Macneal R H, Harder R L. 1985. A proposed standard set of problems to test finite element accuracy. Finite Elements in Analysis & Design, 1(1):3-20.

Nicholas C, Henryk S, Ted B. 1986. Improvements in 3‐node triangular shell elements. International Journal for Numerical Methods in Engineering, 23(9):1643-1667.

Olson M D, Bearden T W. 1979. A simple flat triangular shell element revisited. International Journal for Numerical Methods in Engineering, 14(1):51-68.

算例2 扭曲梁算例

扭曲悬臂梁问题是由MacNeal和Harder（1985）提出，被广泛采用作为三角形壳单元的经典验证算例。如图3所示是扭曲悬臂梁的示意图，从固定端到自由端梁的转角为90度，采用2x12的网格对梁进行离散，为梁的自由端添加两种不同荷载工况，（1）P=1, Q=0（2）P=0, Q=1。表2是单元DKGT沿荷载方向的位移计算结果，由数据可知DKGT单元具有较好的精度。

图3 扭曲悬臂梁示意图

表2 DKGT扭曲悬臂梁问题计算结果对比

参考文献

Macneal R H, Harder R L. 1985. A proposed standard set of problems to test finite element accuracy. Finite Elements in Analysis & Design, 1(1):3-20.

算例3 工字钢屈曲算例

如图4所示是一个工字梁。工字梁算例采用弹性材料，取 E=206GPa，泊松比v=0.3。工字梁两端简支，在腹板的每个节点上添加一个面外方向的压力以模拟初始缺陷，荷载大小取为p=0.5，然后在一端施加压力来模拟工字梁的屈曲。将工字梁的腹板采用4x20的网格划分，翼缘采用2x20网格划分，并分别采用能模拟几何非线性的四边形壳单元NLDKGQ和三角形壳单元NLDKGT进行计算。图5为梁端中点沿压力方向的位移-荷载曲线。图6为工字钢腹板变形情况。

图4 工字梁尺寸和网格划分示意图

图5 工字梁端中点沿压力方向位移-荷载曲线

（a）工字梁腹板中线变形图

（b）工字钢变形图（变形10倍放大）

图6工字钢腹板变形图

根据图5可知，四边形单元NLDKGQ的位移曲线在刚刚过了拐点就终止了，这是因为在NLKDGQ单元进行计算时，迭代次数超过了最大限制，即迭代不收敛，这可能是屈曲问题引起的四边形单元翘曲导致的。而NLDKGT单元则在拐点之后可以稳定的计算，说明三角形壳单元NLDKGT对翘曲问题有较好的改善。从图6可知，在NLDKGQ曲线终点时，两种单元计算得到的腹板变形形状相同，这证明了NLDKGT单元拥有较好的计算精度，同时三角形壳单元NLDKGT相比于四边形壳单元对屈曲问题有更好的适应性。

上述算例可以在以下网址下载

http://pan.baidu.com/s/1mit4FBq

欢迎大家下载试用，并及时和我们交流发现的问题。相应单元的源代码会在相关测试完成后发布。

需要说明的是：在有限元过去半个多世纪的发展历程中，诸多优秀的研究者提出了很多非常好的单元类型。我们的工作只是这些工作中的沧海一粟，我们只是希望给大家提供一个可以用的单元，抛砖引玉，共同推动结构弹塑性分析的发展。

相关研究：