SPH-FEM耦合算法在爆炸螺栓解锁分离过程中的应用

初文华，张阿漫，明付仁，杨文山

(哈尔滨工程大学 船舶工程学院，哈尔滨 150001)

爆炸螺栓的分离解锁是一个相当复杂的过程，目前尚无十分理想的理论计算方法。由于高能炸药的爆轰及冲击波对结构的破坏作用将导致结构大变形等问题，传统的有限元方法在模拟爆炸螺栓的解锁分离过程时将遇到极大阻碍，甚至是计算崩溃[1]。相比于传统的基于网格的有限元法，光滑粒子流体动力学(SPH)方法[1]则在这方面具有很大的优势。它将系统的状态用一系列包含着独自的材料性质且按照守恒控制方程规律运动的粒子来描述，同时具有无网格特性和拉格朗日性质，因此能够克服在计算中与大变形相关的困难。但是光滑粒子流体动力学(SPH)方法在每个时间步内都要对问题域内的所有粒子进行搜索，对于三维问题来说，这将消耗大量的计算时间，其计算效率相比于有限元方法大大降低。为了解决上述问题，综合了FEM方法计算效率高和SPH方法计算精度高的优点的耦合SPH-FEM计算方法被提出。

关于 SPH 与 FEM 的耦合算法，Johnson[2]和 Attaway 等[3]较早开展了相关研究，此后国内外学者[4－8]在其工作的基础上分别对SPH-FEM耦合算法进行了不同程度的改进，且逐渐将其应用到高速冲击、侵彻等问题的数值模拟。然而对于该方法在爆炸冲击问题方面的应用，相关研究开展的则相对较少。因此，本文基于SPH-FEM耦合算法，采用复合链表搜索方式，在充分提高计算效率的同时，计及爆炸螺栓与其连接结构间的相互作用，建立爆炸螺栓解锁分离过程三维数值模型，模拟分析爆炸螺栓的分离冲击特性及结构在分离冲击激励作用下的动响应特性。

1 计算方法

1.1 控制方程

在爆炸冲击的强烈作用下，固体材料的力学特性将会产生类似流体性质的变化，此时运动方程和高压状态方程成为控制材料力学行为特性的关键描述。本文采用计及材料强度的流体动力学控制方程[1]计算在爆炸螺栓冲击激励作用下结构的动力学特性，基于SPH方法离散后的控制方程如式(1)所示。

式中:p，ρ，e，vα，xα，σαβ，t分别代表压力、密度、内能、速度分量、空间坐标、总应力张量和时间;为光滑核函数及其导数，本文采用分段三次样条函数作为光滑核函数;同时，式中还考虑了人工粘度∏ij和人工热量Hi的影响。

1.2 复合链表搜索

在SPH方法的实现过程中，邻近粒子搜索很大程度上决定了计算的效率。在光滑长度为空间常量的情况下，链表搜索法[1]是提高搜索效率的一种有效的方法。它通过将胞元看做簿记式的记录形式来节省大量的计算时间。将所有的粒子都分布在胞元内，并且使用链表定义，在NNPS过程中只要使用任一组粒子进行搜索即可，从而在很大程度上节省计算时间。

在实现链表算法时，要在问题域上铺设一临时网格。同时将每个粒子都分布在网格单元内，并通过简单的存储规则将每个网格内的所有粒子连接起来。若每个单元内的平均粒子数量足够小，链表搜索法将在极大程度上提高计算效率。然而，当问题域的跨度较大，由此导致的临时网格区域铺设过大时，对于三维问题来说，搜索的网格总数将会大大增加，从而导致搜索效率在一定程度上降低。本文所计算的模型中，每个螺栓的尺寸并不是很大，然而两个螺栓之间的距离却是螺栓尺寸的几倍，如果将两个螺栓放在同一个问题域中，采用链表搜索法不但不会提高搜索效率，反而会因过多没有粒子存在的网格参与搜索循环而使得计算时间大大增加。

因此，基于模型距离跨度大的特点，本文采用复合链表搜索方法来提高搜索效率。如图1所示，在每个螺栓所在的区域分别铺设一局部临时网格，由于两个螺栓距离很远，一个螺栓的粒子不会存在于另一个螺栓的支持域内。由此，每个局部临时网格内的粒子仅参与该螺栓所在区域的粒子搜索，而不会对另一个螺栓区域的粒子搜索产生影响。此时，两个螺栓之间的区域由于没有粒子的存在则不需要铺设网格，网格总数大大降低，同时每个螺栓所在区域分别采用各自的链表搜索系统，则每个搜索区域的网格数仅为整个问题域网格总数的一半，因此搜索效率在极大程度上得到了提高。

图1 复合链表搜索方式示意图Fig.1 Compound linked list search mode

1.3 耦合SPH-FEM技术

在本文中，光滑粒子流体动力学(SPH)方法与有限元(FEM)方法耦合的基本思想是在初始时刻爆炸螺栓及其连接结构分别采用SPH粒子和有限元网格进行建模，在计算过程中，使交界面上的粒子与有限元边界保持接触但不允许相互侵入。

如图2所示，当接触面处的光滑粒子侵入有限元界面时，按下式对该边界粒子进行位置和速度的调整:

图2 SPH-FEM耦合计算方法原理图Fig.2 Schematic diagram of coupling SPH-FEM method

与此同时，为了实现耦合计算，接触面处的有限元界面上的单元压力由各单元中心附近的光滑粒子压力值取平均值得到:

其中:Pb_FEM为接触面处有限元单元上的压力，Pjb_SPH为该单元中心附近的光滑粒子压力值;K为该有限元单元中心点“支持域”内的光滑粒子总数，本文中该“支持域”尺寸定义为与光滑粒子的支持域尺寸相同。

2 计算模型

本文针对爆炸螺栓基本结构及工作原理，采用耦合SPH-FEM算法建立爆炸螺栓爆炸解锁分离数值模型。爆炸螺栓及其内部装药采用SPH粒子进行离散，而其连接结构则采用有限元网格单元进行模拟计算。

2.1 材料模型

本文的计算模型中，爆炸螺栓及其连接结构均采用钢材料进行模拟。在采用SPH方法对爆炸螺栓进行建模时，将Jaumann应力率应用到材料本构模型的基本方程中，即在式(1)中，剪切应力ταβi与应变率εαβi存在如下关系:

高校应高度重视政府会计制度落实，高校主要领导要高度重视，确保这项工作能够落到实处。高校要成立由单位主要领导牵头，财务、资产等相关业务部门组成的领导小组，规划方案、明确目标、落实责任，加强部门协同，按责任分工。

式中:Rβγ为扭转率应变率 εαβ=

在应力更新的过程中，采用Von Mises屈服准则来判断材料的屈服状态，材料的屈服应力Y采用Johnsoncook材料屈服模型计算，其表达式为:

若SPH方法的计算过程中，粒子i的等效应力σeff超过材料的屈服强度Y，则该粒子进入塑性屈服状态，其剪切应力ταβ按下式进行修正:

2.2 爆炸螺栓三维SPH 模型

本文基于常用爆炸螺栓的基本结构，建立简化的爆炸螺栓三维SPH数值模型如图3所示，对于爆炸螺栓基本结构中的螺栓头、螺母、装药及用以减少爆炸螺栓冲击环境的防护盒，该SPH模型中均进行了模拟。

在SPH-FEM耦合计算方法中，一个关键问题就是SPH粒子与有限元网格在接触面上的耦合，而解决这一问题的首要因素就是保证二者初始时刻在接触面处吻合。本文在自主编写计算程序进行爆炸螺栓三维SPH建模的过程中，严格按照有限元模型中两个螺栓孔的尺寸和位置，分别建立了上下两个爆炸螺栓的数值模型。

图3 爆炸螺栓三维简化SPH模型Fig.3 Three-dimensional simplified SPH model of explosion bolt

2.3 分离结构有限元模型

在数值建模的过程中，将爆炸螺栓的连接结构简化为如图4(a)所示的对称半圆框结构。图中的圆孔部分代表四个相同尺寸爆炸螺栓的安装位置。

图4 分离结构及爆炸螺栓位置示意图Fig.4 Separating structure and position of the four explosion bolts

3 计算结果

3.1 爆炸螺栓爆炸冲击特性分析

为分析爆炸螺栓解锁分离过程冲击激励的特点，图5给出了爆炸螺栓解锁分离过程中典型时刻的压力云图。

从图中可以看到，两爆炸螺栓从各自的装药中心同时起爆后，冲击波从炸药中心向螺栓头、螺母等周围结构中传播，并随时间迅速衰减;炸药周围的结构逐渐膨胀变形，直至最终产生破坏断裂。由于模型尺寸较小，且为三维模型，为了更清楚地看到冲击压力在结构中的分布及其随时间的变化，图中给出了放大后的爆炸螺栓结构压力分布图，并将图例中压力最大值的显示范围调小。此外还可以看到，由于上下两螺栓结构尺寸完全相同，装药量也完全相同，两爆炸螺栓在解锁分离过程中的冲击特性也几乎完全相同。

图5 爆炸螺栓解锁分离过程冲击压力分布图Fig.5 Impact pressure at different moments in the unlocking process of explosion bolt

图6 冲击激励作用下爆炸螺栓连接结构不同时刻应力分布图Fig.6 Stress distribution of the structure connected with explosion bolts under the impact excitation at different moments

3.2 爆炸螺栓冲击激励作用下的结构响应分析

图6 给出了爆炸螺栓分离冲击激励作用下不同时刻连接结构的MISES应力分布图。从图6中可以看出，当爆炸螺栓内装药起爆后，冲击波通过接触面向连接结构中传播;随着时间增加，冲击波在传播过程中逐渐衰减，结构中的MISES应力响应也逐渐减小，最终趋于稳定;而爆炸螺栓连接的前后两部分结构也在爆炸分离冲击激励的作用下最终分离。

从图6中还可以看到，应力较大的位置始终出现在上下两螺栓孔附近区域。图7进一步给出了冲击激励作用下连接结构上下两螺栓孔附近位置处的MISES应力时程曲线。

从图7中可以看到，上下两螺栓孔附近测点的MISES应力响应峰值均达近240 MPa，因此在爆炸螺栓解锁分离过程的前期，强烈的冲击激励极有可能对爆炸螺栓附近连接结构可能产生较高的破坏效应。

为进一步分析爆炸螺栓解锁分离过程中连接结构在冲击激励作用下的响应特性，图8给出了连接结构在上下两螺栓孔附近区域的三方向冲击加速度时程曲线。

在爆炸冲击激励作用下，应力波在结构中以三维空间波的方式传播，结构中的加速度响应在极短时间内迅速突跃到最大值，然后随着时间增加幅值逐渐减小为零，整个持续时间不过几十毫秒，然而加速度脉宽仅为零点几毫秒，可见应力波作用期间波形变化极快。此外，即使是同一点，加速度峰值在三个方向的大小也是不同的。

图7 上下两螺栓孔附近MISES应力时程曲线Fig.7 MISES stress-time history curves of points near bolt holes

工程应用中，主要把加速度及其计算谱作为爆炸冲击环境评估的主要参数，设计规范和试验规范都是以冲击响应谱为参数给出的。最大谱是爆炸冲击试验中最常用的一种，它是一组理想单自由度振子对基础运动的最大响应随振子频率ω变化的图谱。图9给出了根据图8中六组冲击加速度响应时程曲线得到的六组最大冲击响应谱。

图8 上下两螺栓三个方向的加速度响应时程曲线Fig.8 Acceleration-time history curves of the two explosion bolts in the three direction

图9 上下两螺栓三个方向的加速度响应谱Fig.9 Acceleration response spectrums of the two explosion bolts in the three direction

4 结论

本文针对爆炸螺栓解锁分离过程的工作原理及载荷特点，采用光滑粒子流体动力学(SPH)方法与有限元(FEM)方法相结合的技术手段，对爆炸螺栓的解锁分离过程进行了数值模拟研究，得出了以下主要结论:

(1)两爆炸螺栓分别从各自的装药中心同时起爆后，冲击波分别向两爆炸螺栓结构中传播，并随时间迅速衰减;同时，炸药周围的结构逐渐膨胀变形直至最终产生破坏断裂;

(2)在爆炸螺栓分离冲击激励作用下，冲击波从螺栓孔附近开始向结构四周传播，并随时间增加逐渐衰减，结构中的MISES应力响应也随之逐渐减小，最终趋于稳定;爆炸螺栓连接的前后两部分结构在爆炸分离冲击激励的作用下最终分离;

(3)在爆炸冲击激励作用下，应力波在结构中以三维空间波的方式传播，结构中的加速度响应在极短时间内迅速突跃到最大值，然后随着时间增加幅值逐渐减小为零，整个持续时间不过几十毫秒。

[1]韩 旭，杨 刚，强洪夫译.光滑粒子流体动力学——一种无网格粒子法[M].长沙:湖南大学出版社，2005.

[2]Johnson G R，Stryk R A，Beissel S R，et al.An algorithm to automatically convert distorted finite elements into meshless particle during dynamics deformation[J].Int J Impact Eng.，2002，27(10):997－1013.

[3]Attaway S，Heinstein M，Swegle J.Coupling of smoothed particle hydrodynamics with the finite element method [J].Nucl.Eng.Des.，1994，150(2 －3):199 －205.

[4]王 吉，王肖钧，卞 梁.光滑粒子法与有限元的耦合算法及其在冲击动力学中的应用[J].爆炸与冲击，2007，27(6):522－528.

[5]卞 梁，王肖钧，章 杰.SPH/FEM耦合算法在陶瓷复合靶抗侵彻数值模拟中的应用[J].高压物理学报，2010，24(3):161－167.

[6]Jun-xiang X U，LIU X L.Analysis of structural response under blast loads using the coupled SPH-FEM approach[J].Journal of Zhejiang University，2008，9(9):1184 －1192.

[7]Zang M Y，Chen H，Lei Z.Simulation on high velocity impact process of windshield by SPH/FEM coupling method [J].2010 WASE International Conference on Information Engineering，2010:381 －384.

[8]Fourey G，Oger G，Touze D Le，et al.Violent fluid-structure Interaction simulations using a coupled SPH/FEM method[J].Materials Science and Engineering，2010，10:1－8.