壳单元和实体单元模拟爆炸荷载作用钢板动力响应的差异性比较

朱俊杰，卢永刚，冯星奎，陶俊林

(1. 西南科技大学 土木工程与建筑学院，四川 绵阳 621010；2. 中国工程物理研究院 总体工程研究所，四川 绵阳 621000；3. 火箭军驻中国工程物理研究院 型号办事处，四川 绵阳 621900)

壳单元和实体单元模拟爆炸荷载作用钢板动力响应的差异性比较

朱俊杰1，卢永刚2，冯星奎3，陶俊林1

(1. 西南科技大学 土木工程与建筑学院，四川 绵阳 621010；2. 中国工程物理研究院 总体工程研究所，四川 绵阳 621000；3. 火箭军驻中国工程物理研究院 型号办事处，四川 绵阳 621900)

以舰船舱室钢板为研究对象，利用动力分析有限元程序 Ls-dyan，分别采用壳单元和实体单元 2 种建模方式对爆炸荷载作用下钢板的动力响应进行数值模拟研究。对比钢板的应力分布云图、中线上各点最大位移曲线图，发现 2 种方法计算结果与试验值大致吻合，采用壳单元建模在保证计算结果准确性的前提下可以大大提高计算效率。如何正确地选择单元进行分析是数值模拟很重要的一个环节，通过对壳单元和实体单元计算结果的初步对比，为舰船舱室数值模拟研究如何选取单元提供参考。

船用钢板；爆炸荷载；数值模拟；壳单元

0 引 言

舰艇是由纵横构架和壳板结构组成，其中壳板主要承受水平荷载。舰船舱室结构承受的主要载荷之一是战斗部爆炸产生的冲击波载荷[1]，钢板在受到具有传播速度快、峰值荷载大、作用时间短等特征的爆炸冲击波载荷作用时，其受力特性和动态响应显著异于静载情况下，且十分复杂，故而研究钢板在爆炸荷载作用下的动力响应很有必要，是研究舱室结构的基础。目前国内对舰船舱室方面的研究大多集中在舱室结构的仿真分析上，朱建方等[2]对舰艇舱室在爆炸荷载下的动态响应过程进行了数值模拟，得到了舰艇舱室在爆炸冲击波作用下的响应特性和破坏机理。鄢顺伟等[3]仿真计算了不同壁厚舱室结构在反舰导弹战斗部内爆作用下的毁伤效应。但是，对船用钢板在爆炸荷载作用下的响应研究较少。本文采用数值模拟的方法对爆炸载荷下四边固支的钢板进行研究。在利用数值分析研究时，单元类型的选取对计算效率有很大的影响。钢板属于板结构，实体单元与壳单元均可以对其进行模拟，但是 2 种建模方法得到结果的差异值得探讨。通过分析两者的差别，可以更好地为钢板使用何种单元在保证其结果准确性的基础上又能提高计算效率提供参考。为此，本文以文献中试验为基础[4–5]，在验证数值计算方法可靠的基础上，分别采用壳单元与实体单元对四边固支的方形板在爆炸荷载作用下的动力响应进行数值模拟研究，考察使用不同的建模单元类型对计算结果的影响[6]，并对结果的差异性进行了比较与分析。

1 Shell163 单元介绍

Shell163 是一个具有 4 节点单元，有弯曲和膜特征，可施加平面和法向载荷。单元在每个节点上有 12个自由度，分别为在节点 x，y 和 z 方向的平动、加速度、速度和绕 x，y 和 z 轴的转动。Shell163 单元一共有 12 种算法，这里选第 10 种算法，对于翘曲情况可以得到正确结果。和实体单元一样，积分点的个数直接影响着 CPU 的时间。典型地，对于弹性材料沿厚度方向需要 2 个积分点，而对于塑性材料则需要 3 个或更多的积分点[7]。

2 数值计算模型

2.1 算法选取

分析主要涉及到流体流动、固体大变形以及流体与固体相互作用等问题，Ls-dyna 程序中提供了 Lagrange，Euler 和 ALE（任意拉格郎日欧拉算法）3 种不同的算法，其中 ALE 算法兼具 Lagrange 算法和 Euler算法的优点，可以描述流体与固体的相互耦合作用；同时在求解过程中网格还可以根据定义的参数适当调整位置，克服网格严重畸变的问题[8]。因此，本文中选取 Ls-dyna 程序提供的 ALE 算法建立流体与固体间的联系，实现流固耦合动态分析，模拟方形钢板在爆炸荷载作用下的响应和破坏。

2.2 材料本构与状态方程

结构材料在爆炸冲击载荷作用下会产生巨大的变形，与受静载荷作用的情况相比，动态屈服应力比静态屈服应力有明显提高，材料力学性能表现相对也更为复杂。固支方板采用低碳钢，钢板选用与应变率相关的各向同性塑性随动硬化模型，材料参数见表 1。

表 1 钢板材料参数Tab. 1 Steel plate material parameters

炸药在空气中爆炸时，会在极短的时间内把化学能全部转变为热能，导致爆轰产物向周围迅速膨胀而做功，形成空气冲击波。文中采用 *HIGH_EXPLOSIVE_BURN 模型作为炸药的材料模型，采用 Jones-Wilkins-Lee（JWL）状态方程来描述，即

式中：A ，B，R1，R2和 ω 为 JWL 状态方程参数；E为炸药的内能；V 为当前相对体积。

对空气采用 *MAT_NULL 材料模型以及线性多项式 *EOS_LINEAR_POLYNOMI 状态方程加以描述，即

式中：C0～C6为常量；e 为初始体积内能；μ = 1/V0–1；V0为相对体积；C1～C4均为 0。

2.3 有限元模型

1） 空气网格的验证

在数值模拟中，网格的划分至关重要，网格的大小不仅关系着数值模拟计算的精度，也关系着数值模拟计算所用的时间。一般而言，网格越小则计算的结果越接近于真实情况，但随着网格的减少计算时间也急剧增大，所以选择合理的网格尺寸很有必要。网格验证几何模型尺寸为空气 120 mm × 120 mm × 60 mm，炸药为 27 mm × 27 mm × 27 mm。利用对称性取 1/4 模型进行构建，模型如图 2 所示。

表 2 炸药与状态方程参数Tab. 2 Explosive and equation of state parameters

表 3 多项式状态方程参数Tab. 3 Polynomial equation of state parameters

空气与炸药均采用 8 节点实体单元，TNT 当量为 20 g，炸药中心位置起爆。在对称面上施加对称约束，四周（除对称面）施加固定边界条件。为了模拟出无限空气域，空气模型边界施加透射边界。工况如表 4 所示。

表 4 空气网格验证工况表Tab. 4 Air grid authentication working table

上述工况在距离炸药中心 6 cm 和 8 cm 位置处的超压如图 3～图 4 所示。

由图 5～图6 可以看出，随着网格的减少，压力的峰值越来越大，当网格取 0.3 cm 和 0.2 cm 时压力峰值的差值约为 5%，考虑到计算精度和时间的因素空气网格取 0.3 cm。

2） 炸药网格的验证

根据上面网格确定炸药网格的验证工况如表 5 所示。

表 5 炸药网格验证工况表Tab. 5 Explosive grid authentication working table

由图 5～图 6 可知，随着炸药网格的减少，压力峰值有所增加，当网格取为 0.2 cm 和 0.1 cm 时的压力峰值差值约为 2%，炸药网格选取 0.2 cm 较为合适。

3） 2 种单元计算模型

采用 Ansys/Ls-dyna 建立有限元模型，对于固支钢板在爆炸载荷的作用过程进行了数值模拟研究。固支钢板尺寸为 63 mm × 63 mm，厚度为 1.6 mm，建立 1/4模型。方板尺寸如图 7 所示。

数值模型由炸药、空气、钢板 3 部分组成，采用壳单元和实体单元分别建立板模型，2 种模拟方式除了方板的单元类型不同，其他参数均一致。空气单元尺寸为 3 mm，炸药单元尺寸为 2 mm，板单元为 0.4 mm，有限元模型如图 8～图9 所示。

2.4 失效准则

采用 *MAT_ADD_EROSION 对钢板材料失效准则进行定义，塑性失效应变 εmax= 0.30，考虑到在冲击荷载作用下薄壁材料容易出现断裂、绝热剪切等破坏情况，在计算模型中增加最大剪应变失效准则，剪切失效应变 τmax= 0.20[9]。利用单元删除技术来模拟钢板的破坏情况，即当单元变形引起的塑性应变或剪切应变超过临界值时，就认为该单元破坏。

3 壳单元与实体单元模拟结果比较

3.1 方形钢板的动态响应

1）t = 0.15 ms 时刻，方形钢板 Von Mises 应力云图对比。

从图 10 可知，无论是采用壳单元还是实体单元，方形板中的应力分布大致相同，最大值均出现在炸药最上方中心位置处。同一时刻，实体单元中最大应力值为 503.6 MPa，壳单元中最大应力值为 504.7 MPa。

从图 11 可知，方形钢板的变形大致呈现为中心处的外鼓，呈“金字塔”状，开始先发生中心处的翘曲，随着时间的增加，翘曲越来越严重，而后慢慢下降，2 种方法模拟的钢板变形发展趋势与试验结果变形趋势相同。

2）钢板中线上各点处最大位移值对比

图 12给出了在 20 g 炸药量产生的爆炸荷载作用下钢板中线上各点的最大位移曲线。由图中可以看出，2种方法得到的各点处最大位移曲线基本一致，与试验值较为吻合，各点处最大位移差值均在 10% 以内。最大的位移发生在钢板中心位置处，2 种模拟方法得到的结果与试验值差值均在 5% 以内，其他位置处的差值较小。

3.2 两种方法计算时间对比

2 种模型的炸药和空气域单元数目相同，而使用壳单元建模使得钢板单元数目减少 50%。由图 13 可以看出使用壳单元建模计算时间降低约 50%，显著提高计算效率。尤其对于单元数目比较多的模型，壳单元的优势更为突出。

无论是从钢板的动力响应历程还是钢板的变形特征层面上的对比验证，其结构均表明：本次数值模拟的建模技术和参数选取合理、可靠，其数值模拟的结果可以舰船舱室内爆选取单元提供依据。

4 结 语

对方形钢板在爆炸荷载作用下的响应，分别采用壳单元和实体单元进行模拟。对比 2 种方法中板的应力分布和中线处各点的最大位移值，得到的结果大致相同，差值在 10% 以内。结果表明使用实体单元建模和使用壳单元建模对于等厚度方板的计算结果影响不大。而采用壳单元建模相比采用实体单元建模可以显著提高计算效率。

数值分析方法在爆炸荷载作用下舰船舱室的动力响应研究中有着广泛的应用，但如何正确地选择单元进行分析是数值模拟中很重要的一个环节。本文通过对壳单元和实体单元模拟方形钢板计算结果的初步对比，为舰船舱室模型中如何选择单元提供参考。

[1]陈攀, 刘志忠. 舱室内爆冲击波载荷特性及影响因素[J]. 舰船科学技术, 2016, 38(2): 43–48. CHEN Pan, LIU Zhi-zhong. Research on loading of explosive and influencing factors inside closed cabin [J]. Ship Science and Technology, 2016, 38(2): 43–48.

[2]朱建方, 王伟力, 曾亮. 舰艇舱室内爆毁伤的建模与仿真分析[J]. 系统仿真学报, 2009, 21(22): 7066–7068.

[3]杜志鹏, 李晓彬, 夏利娟, 等. 反舰导弹攻击舰船舷侧防护结构过程数值仿真[J]. 哈尔滨工程大学学报, 2006. 27(4): 484–487.

[4]LONGDON G S, CHUNG Ki Y S, NURICK G N. Experimental and numercial studiesion the response of quadrangular stiffened plates. Part I: subjected to uniform blast load [J]. In-ternational Journal of Impact Engineering. 2005, (31): 55–83.

[5]LONGDON G S, CHUNG Ki Y S, NURICK G N. Experimental and numercial studiesion the response of quadrangular stiffened plates. Part II: localist blast loading [J]. International Journal of Impact Engineering. 2005, (31): 85–11.

[6]徐新光. 梁单元和实体单元模拟框架结构的差异性比较[J].四川建材, 2014, 40(5): 46–48.

[7]ANSYS工程结构数值分析[M]. 北京. 人民交通出版社, 2012: 11–13.

[8]ANSYS/LSDYNA在爆炸与冲击领域内的工程应用[M]. 北京:中国建筑工业出版社, 2011.

[9]李帆. 船舰室内爆炸破坏的数值模拟和试验研究[D]. 北京: 国防科技大学, 2012.

Comparision of differentiation with the steel plates' dynamic response subjected to blast loads using shell element and solid element modeling methods

ZHU Jun-jie1, LU Yong-gang2, FENG Xing-kui3, TAO Jun-lin1
(1. School of Civil Engineering and Architecture, Southwest University of Science and Technology, Mianyang 621010, China; 2. Institute of Structural Mechanics,China Academy of Engineering Physics, Mianyang 621900, China; 3. The Office of Rocket Force, China Academy of Eugineering Physicies, Mianyang 621900, China)

Take the steel plates of warship cabins as research objects, numerical simulation study on the response of the quadrangular steel plates subjected to blast loads is presented by dynamic analysis finite element procedure Ls-dyan3d, using modeling methods of shell element and solid element respectively. Comparing pressure distribution cloud picture and each point of midcourt line miximum displacement curve graph, it can be found that transformation trend of steel plates is similar for different thickness while each point of midcourt line miximum displacement is discriminating. However results from two methods are alike for equal thickness. Meanwhile, it can improve calculative effciency. How to choose which kind of element types correctly to analysis is an important link. Through the preliminary contrast study of calculative results of shell element and solid element, it can provide reference for choosing element type in warship cabins numerical simulation.

marine steel plate；blast loads；numerical simulation；shell element

TJ03

A

1672–7619(2017)03–0018–05

10.3404/j.issn.1672–7619.2017.03.004

2016–09–18；

2016–10–27

国家自然科学基金资助项目(11502258)

朱俊杰(1991–)，女，硕士研究生，研究方向为结构冲击与爆炸。

陶俊林(1972–)，男，博士，教授，E-mail:junlintao@126.com。