光滑粒子流体动力学(smoothed particle hydrodynamics, SPH)方法^[1-2]是近年来兴起并快速发展的一种无网格Lagrange型粒子方法。该方法中的计算域通过粒子实现离散，由于其不依赖网格，且SPH粒子本身的Lagrange特性，因此可以轻松实现物质运动的跟踪和结构/材料变形的描述，克服了传统有限元方法(finite element method, FEM)处理大变形问题时存在的固有缺陷，广泛应用于爆炸冲击动力学等问题研究。杨刚等^[3]采用SPH方法对不同药型罩线性聚能射流及后效侵彻过程进行模拟，充分展示出SPH方法模拟大变形问题的优越性。然而采用SPH法计算时，需要通过搜索粒子影响域内相邻粒子确定每个粒子的物理参量，当粒子数较大时，计算耗时巨大，导致计算效率远远低于FEM。刘天生等^[4]分别采用FEM和SPH法对大变形冲击动力学问题开展研究，发现在同一模拟中FEM的计算耗时约为SPH法的1/5。

为了充分发挥SPH和FEM的各自优势，通过少量SPH粒子得出冲击载荷，而主体结构的响应采用FEM进行计算，SPH耦合FEM(SPH-FEM)算法研究由此展开。Johnson等^[5]和Attaway等^[6]较早应用该方法进行侵彻贯穿数值模拟，采用SPH粒子计算大变形区域，既消除了FEM中的网格畸变问题，又提高了计算精度，同时结合FEM处理小变形区域，大大提高了计算效率。王吉等^[7]在此基础上，对二维碰撞问题进行分析，初始结构采用FEM建模，计算过程中大变形区域的FEM网格单元自动转换成SPH粒子，该耦合算法取得了不错的计算结果。卞梁等^[8]基于王吉的思想，引入金属和陶瓷的本构模型，成功模拟了陶瓷复合靶板的侵彻过程。初文华等^[9]、朱东俊等^[10]基于SPH-FEM耦合算法，计算分析了简单结构在爆炸分离过程中冲击环境及结构响应特性。张志春等^[11]提出了一种新型耦合SPH-FEM的接触转换算法，对弹体冲击钢板进行了三维数值模拟，计算结果与实验结果吻合良好。姜忠涛等^[12]采用SPH-FEM耦合算法模拟了水下近场接触爆炸，实现了从装药爆轰到结构物变形失效的全过程模拟。然而，虽然人们在SPH-FEM耦合算法开发上取得了一系列成果，但是在算法应用方面仍然处于较简单的模型计算分析阶段，对于大型具有复杂力学特性的复合材料结构的计算分析鲜有探索。

迄今为止，复合材料已经广泛运用于航空航天、船舶等领域。复合材料蒙皮结构既满足了飞行器结构的强度和刚度要求，又有效地减轻了结构质量^[10]。对于复合材料蒙皮结构飞行器的爆炸分离研究，现今主要依靠试验研究，然而试验研究的花费巨大，试验周期长，具有一定的危险性，相关报道甚少，相应的数值研究依然欠缺^[11-12]。本研究中, 基于SPH理论，在初文华等^[9]、朱东俊等^[10]工作的基础上，采用SPH-FEM耦合算法建立复合材料蒙皮结构飞行器的爆炸分离模型，对其爆炸分离过程进行计算模拟，分析并总结复合材料蒙皮结构飞行器的相关冲击动力学特性，旨在为SPH-FEM耦合算法的工程实际应用以及航天飞行器结构设计与防护提供参考。

1.   计算方法

1.1   SPH算法模型

SPH方法的控制方程为:

式中：ρ、u、e、V分别为粒子密度、速度、比内能和体积；h为光滑长度；g为重力加速度；r为粒子位置矢量；p_JWL为JWL状态方程函数；p为压力；W_ij为SPH近似时的核函数，本文中选用三次样条函数作为核函数；$\nabla$_iW_ij为W_ij关于粒子i的偏导数，其中j为粒子i的相邻粒子；α为人工黏性系数，本研究中α取0.5；c为爆轰波的传播速度，ρ₀为炸药初始密度，本研究中采用TNT炸药，即c=6 930 m/s，ρ₀=1 630 kg/m³。式(1)中前3个方程分别为质量、动量和能量守恒方程。通过粒子近似，计算ρ、u和e对时间t的全导数，在时域进行预测校正积分，实现炸药爆轰过程模拟。通过Jones-Wilkins-Lee(JWL)状态方程^[1-3]建立ρ、e和p之间的显式关系，将3个守恒方程解构，分别求积分。此外，粒子体积V通过粒子质量与密度的比值进行更新。

1.2   SPH-FEM耦合算法

在复合材料蒙皮结构飞行器的爆炸分离过程中，爆炸冲击载荷主要由若干个爆炸螺栓的解锁分离效应产生。爆炸螺栓在航天飞行器分离前后起着关键作用：分离前，需要保证两个分离部件可靠连接；分离后，要求顺利分离，且不影响续航飞行器舱段的正常任务执行^[13]。

对于爆炸螺栓的解锁分离过程，初文华等^[9-10]做了相关数值研究工作，结果显示，其采用的SPH-FEM耦合算法能够较有效地模拟爆炸分离过程中螺栓的冲击特性及简单分离结构的响应规律。为此，本研究中在基于初文华等^[9-10]耦合算法基本原理的基础上，引入复合材料的损伤模型及三维复杂结构的数值建模，对具有复杂结构的复合材料蒙皮结构飞行器的爆炸分离过程开展进一步研究。

计算结构中爆炸螺栓的分布如图 1(a)所示。与简单结构相比，三维复杂结构的计算量大大增加，在模拟螺栓爆炸分离过程中，如果对整个结构端框铺设背景网格，即将6个螺栓的计算置于同一计算域中，则不可避免地将螺栓爆炸区域外的大面积非计算区域引入，使计算效率大幅下降。因此，在SPH算法的实现过程中引入复合链表搜索法^[6]，结合计算模型的结构特点，对关键的螺栓爆炸区域展开相互独立的区域搜索计算，不同螺栓的粒子仅在各自爆炸区域内的背景网格进行链表搜索，不存在其他螺栓粒子搜索计算的影响与非计算区域的引入，从而实现在保证计算精度的前提下提高耦合算法的计算效率，为SPH-FEM耦合算法的实际工程应用提供了可能性。

图  1  复合链表搜索原理

Figure  1.  Search principle of compound linked list

下载: 全尺寸图片 幻灯片

1.3   复合材料损伤模型

在常用的复合材料层合板的单层板失效准则中，一般简单地认为在应力满足一定条件下单层板随即发生损伤破坏，在材料失效准则中未考虑失效模式或损伤机理的影响^[13]。考虑到层合板在实际高频冲击作用下伴随着混合损伤失效模式，忽略层间应力在层合板损伤过程中的影响，采用二维Hashin准则对层合板在高频爆炸冲击过程中的损伤进行判定^[14]，具体表达式如下。

纤维拉伸失效：

纤维压缩断裂：

基体拉伸或剪切失效：

基体压缩或剪切失效：

式中：σ₁、σ₂、τ₁₂分别为复合材料层合板单层板的纵向平面应力、横向平面应力及剪应力；S₁₂为单层板1-2方向的剪切强度；X_t、X_c分别为单层板的纵向拉伸和压缩强度；Y_t和Y_c分别为单层板的横向拉伸和压缩强度。

在爆炸冲击过程中，如果层合板中的材料满足式(2)~(5)中的任一失效准则，则需对材料刚度进行折减。现行的折减方式主要有两种：一种是直接对材料参数进行折减，然而折减系数需要通过大量试验确定，不易取得；另一种是基于能量对参数进行折减，从而有效地降低甚至消除有限元计算中软化段的网格依赖性，其中断裂能参数可以通过相应文献计算得到。本模型中采用ABAQUS软件自带的基于能量的折减模式^[15]。在材料损伤演化阶段，引入一个简单的基于能量的线弹性模型，满足初始损伤阀值后，材料积累的损伤能量达到破坏临界阀值(断裂能)时，材料将遵循线性软化定义的应力-应变关系进行折减。在ABAQUS的计算定义中，只要设置对应的4个临界断裂能参数即可。

2.   数值模型有效性验证

在实际爆炸分离试验中，飞行器舱段采用如图 2(a)所示的结构，主要由主体舱段壳体和后端分离结构组成，其中主体舱段壳体分为常规铝合金内壳和复合材料蒙皮外壳，分离结构由分离前、后板组成。分离前、后板通过6个爆炸螺栓固定；在爆炸分离阶段，6个螺栓同时起爆，在螺栓爆炸分离冲击激励作用下，前、后舱段实现安全平稳脱离。

图  2  复合材料蒙皮结构飞行器数值模型

Figure  2.  Numerical model of composite skin aircraft structure

下载: 全尺寸图片 幻灯片

在建立有限元模型时，需要对实际飞行器结构进行细致的分析和必要的简化处理。爆炸分离时螺栓发生爆炸大变形，与之直接相连的分离结构属于关键结构，结合SPH-FEM耦合算法原理，将分离结构采用实体单元划分，忽略一些装配小孔。然而如果整个飞行器模型都采用实体单元建模，则单元数量巨大，且不可避免地产生大量扭曲单元，不仅降低数值结果的精确性，还可能导致计算不收敛。考虑到主体舱段壳体属于薄壁结构，计及剪切变形的影响，本研究中采用带沙漏控制的壳单元(S4R)加以描述，不仅可以提高计算效率，还能保证计算结果的准确度。最终建立的飞行器结构有限元模型见图 2(b)。参照具体试验条件，主体舱段的蒙皮外壳采用航天领域常用的T300复合材料进行模拟^[16]，相应的层合板铺层形式为[0°/90°]2S，单层板厚度d₀=1.75 mm，共8层，层合板铺层形式见图 2(c)。复合材料层合板的材料参数见表 1，其中：ρ为密度，E₁₁和E₂₂为纵向和横向弹性模量，G₁₂为层板平面的剪切模量，ν₁₂为主泊松比，G_f^t、G_f^c、G_m^t、G_m^c分别为材料损伤阶段中纤维拉伸、纤维压缩、基体拉伸、基体压缩各损伤模式下的临界断裂能参数。

表  1  复合材料层合板材料参数

Table  1.  Material parameters of composite laminate

ρ/(g·cm-3)	E11/GPa	E22/GPa	G12/GPa	ν12	Gft/(kN·m-1)	Gfc/(kN·m-1)
1.472	146.8	11.4	6.1	0.30	89.83	78.27
Yt/MPa	Yc/MPa	S12/MPa	Xt/MPa	Xc/MPa	Gmt/(kN·m-1)	Gmc/(kN·m-1)
66.5	268.2	58.7	1730.0	1379.0	0.23	0.76

下载: 导出CSV 

| 显示表格

在爆炸分离过程中爆炸螺栓为主要研究结构，数值模拟时有必要对常用的爆炸螺栓进行适当简化，建立简化的爆炸螺栓三维SPH模型，如图 3(a)所示。简化模型中，爆炸螺栓包括防护盒、螺栓头、炸药以及螺母4部分。根据飞行器爆炸分离结构，结合实际装配爆炸螺栓工艺要求，可以准确得出分离结构中6个螺栓孔对应的爆炸螺栓空间位置。依据各螺栓的标准空间位置，建立三维SPH爆炸螺栓模型，如图 3(b)所示。

图  3  爆炸螺栓三维SPH模型

Figure  3.  Three-dimensional SPH model of explosive bolt

下载: 全尺寸图片 幻灯片

在爆炸冲击模型中，确保爆炸激励源稳定可靠是十分重要的前提。参考实际爆炸分离试验条件，控制6个螺栓同时起爆，对飞行器舱段的爆炸分离展开数值模拟研究，在此对爆炸激励源附近的后端框选取测点进行数值模型的有效性和准确性验证。为了方便描述，将飞行器舱段的后端框进行象限划分，形成Ⅰ、Ⅱ、Ⅲ、Ⅳ 4个象限(见图 1(b))。

在工程应用中，爆炸冲击环境评估的主要参数是冲击加速度及其相应的冲击谱。相比于非规则性随机振荡的冲击加速度响应，频域分析更适于复杂冲击环境下的结构研究。对于航天器分级、飞行器舱段间爆炸分离等复杂的瞬态振荡冲击环境，采用冲击响应谱描述，以便观察和分析模拟结果与试验值的差异。选择爆炸冲击试验中常用的最大谱加以分析，图谱的横坐标为固有频率，纵坐标为单自由度系统的响应峰值^[10]。

通过选取与实际爆炸分离试验相同的采样频率，计算得到相应的最大冲击响应谱。取特征长度为飞行器舱段长度，即L=3.10 m，特征速度取为声速v=340 m/s，对相应的冲击响应谱进行无量纲化处理，定义无量纲化的加速度和固有频率分别为a和ω，最终得到试验结果与数值结果的无量纲化最大冲击响应谱对比情况。图 4和图 5分别给出了第Ⅱ象限和第Ⅳ象限内螺栓孔附近测点的冲击响应谱对比结果。为了便于分析比较，相应的无量纲化冲击响应谱峰值对比情况列于表 2。

图  4  第Ⅱ象限内螺栓孔附近区域测点的无量纲化冲击响应谱

Figure  4.  Dimensionless shock response spectra of measuring point in quadrant Ⅱ

下载: 全尺寸图片 幻灯片

图  5  第Ⅳ象限内螺栓孔附近区域测点的无量纲化冲击响应谱

Figure  5.  Dimensionless shock response spectra of measuring point in quadrant Ⅳ

下载: 全尺寸图片 幻灯片

表  2  无量纲化冲击响应谱峰值对比

Table  2.  Comparison of dimensionless shock response spectrum peak

象限	方向	无量纲化冲击响应谱峰值
		试验值	计算值	相对误差/%
Ⅱ	轴向 径向	0.082 0.115	0.098 0.142	19.51 23.49
Ⅳ	轴向 径向	0.174 0.127	0.199 0.144	14.37 13.38

下载: 导出CSV 

| 显示表格

结合图 4、图 5和表 2可知，在螺栓爆炸冲击激励下，由爆炸分离模型计算得到的分离结构中螺栓孔附近冲击响应与实际试验结果确定的冲击环境在趋势上一致，数值计算结果与试验结果的冲击谱值具有相同量级，虽然数值计算结果与试验结果存在一定差异，但在允许误差范围内，相对误差不超过30%。由此验证了基于SPH-FEM耦合算法建立的飞行器结构爆炸分离数值模型的有效性和准确性，该模型可用于研究飞行器结构在螺栓爆炸分离下的动态冲击响应和爆炸冲击环境预示。

3.   爆炸分离过程冲击环境分析

6个爆炸螺栓在爆炸分离过程中的冲击压力分布如图 6所示。可以看出，由于各螺栓的尺寸和装药情况相同，因此螺栓的爆炸冲击过程基本一致。当螺栓内的装药引爆后，爆炸冲击波从起爆药心处开始向螺栓两端与周围结构迅速传播，距离药心最近的螺栓头和螺母发生剧烈的爆炸膨胀变形；当t=50 μs时，螺栓内的冲击波逐渐衰减，螺母也产生较大的变形破坏；随着炸药继续爆轰以及螺栓结构的持续膨胀变形，当t=100 μs时，装药附近的螺栓头和螺母结构发生大变形并断裂破坏，螺栓彻底失去连接作用。

图  6  螺栓爆炸分离过程中的冲击压力分布

Figure  6.  Pressure profiles in explosion process of explosive bolt

下载: 全尺寸图片 幻灯片

结构爆炸冲击响应情况如图 7所示。螺栓内装药起爆后，爆炸冲击波从药心处向外传播，通过螺栓传递到与之相连的分离结构；当t=0.01 ms时，从螺栓孔传出的应力波开始传至飞行器外壳上；当t=0.1 ms时，应力波基本传遍整个后部分离结构，飞行器后段的应力响应开始逐渐增强，应力波持续向舱段前端传播；当t=5 ms时，整个飞行器舱段都已经受到爆炸冲击激励作用，分离结构的内应力维持在较高幅值，整个舱段结构内的Mises应力响应开始衰减；当t=100 ms时，爆炸冲击激励不断衰减，整体飞行器舱段内的应力响应也随之减小并最终趋于稳定。在整个爆炸分离过程中，螺栓孔附近区域始终维持着较大的应力响应。如图 8所示，在螺栓爆炸冲击激励作用下，最终分离结构实现顺利平稳分离。

图  7  复合材料蒙皮结构飞行器的应力分布

Figure  7.  Stress distribution of composite skin aircraft structure

下载: 全尺寸图片 幻灯片

图  8  分离结构实现爆炸分离

Figure  8.  Explosive separation of separated structure

下载: 全尺寸图片 幻灯片

图 9(a)给出了爆炸分离后板中典型区域的冲击加速度时程曲线。为了考察在高频爆炸冲击下复合材料蒙皮结构的响应情况，选取蒙皮结构前、中、后段典型区域进行冲击响应分析，图 9(b)~(d)给出了蒙皮结构上各段区域典型冲击加速度时程曲线，其中a_a为轴向加速度。

图  9  飞行器结构典型区域冲击加速度时程曲线

Figure  9.  Acceleration time histories of typical areas for aircraft structure

下载: 全尺寸图片 幻灯片

由图 9可知，螺栓的爆炸冲击激励作用具有高频瞬态特点，结构在极短时间内受到剧烈的冲击振荡。在短时间的高频冲击后，结构中的冲击效应迅速衰减，进入低频振动衰减阶段。在爆炸分离中的高频冲击阶段，结构上的冲击加速度呈瞬态高幅振荡，冲击脉宽基本维持在20 ms左右，与实际试验结果以及国内外相关文献描述一致^{[8, 17]}，进一步验证了本耦合算法与所建立模型的有效性和准确性。

在爆炸分离过程中，冲击波率先在分离结构面内传播，之后传入复合材料蒙皮结构，并由后段继续向前传播，从冲击响应时间上也可以看出其在各段的传播过程。分离结构最先受到螺栓爆炸冲击激励作用，从图 9(a)可以看到分离结构后板上的冲击加速度幅值瞬时达到峰值；爆炸冲击波传入复合材料蒙皮结构后，不断向前传播并伴随着能量衰减，从图 9(b)~(d)可以清晰地看出冲击加速度的下降趋势，与试验中实际爆炸冲击响应过程一致。

结合图 7和图 9可以看出，在爆炸冲击的极短作用时间内，飞行器后段区域的冲击响应最剧烈，当冲击波传递到前段区域时，结构上的冲击响应已经衰减得较缓和。在爆炸冲击过程中，分离结构后板中的冲击加速度可达10⁴ m/s²量级；随着冲击波的继续传播与衰减，到复合材料蒙皮结构上的加速度下降到10³ m/s²量级，衰减程度较大。实际上，在爆炸分离过程中，分离后板会与飞行器舱段实现脱离，它与蒙皮外壳的连接部位属于结构接合部。爆炸冲击波通过结构接合部时，爆炸冲击能量与幅值都会出现较大的衰减。结构接合部的连接形式不同时，爆炸冲击的衰减程度也不同。结构接合部的设计特点主要包括以下因素：爆炸螺栓的分布、连接形式和数量，以及分离后板的设计特点。在设计分离后板时主要考虑相关装置盒的安装分布。如图 2所示，模型选用6个爆炸螺栓对称安装连接，在特定的装置盒安装及配重分布情况下，通过分析结构接合部连接情况及相应的冲击响应情况可知，模型中的爆炸冲击幅值衰减约70%。因此，在航天工程的实际爆炸分离中，往往通过改变与优化接合部的设计形式和冲击传递路径，最大限度地降低爆炸冲击引起的冲击响应。

4.   模型计算效率对比分析

为进一步分析本数值算法的计算效率优势，选取螺栓孔处的连接区域(见图 10)进行分析，为了便于SPH粒子离散建模，根据实际模型情况，定义简易局部模型1和模型2。连接区域结构的厚度为35.5 mm，在有限元结构中在厚度方向划分6层。基于模型1和模型2的结构计算区域，进行SPH离散，联立SPH爆炸螺栓模型，得到对应的纯SPH算法模型。对于模型1中的结构计算区域，使用有限元划分，得到96个网格单元，使用30 520个SPH粒子离散；对于模型2中的结构计算区域，也使用有限元划分，得到432个网格单元，使用238 280个SPH粒子离散。因为SPH法采用粒子作用域搜索算法，所以其粒子建模对一致化及均匀化的要求较高。对比模型1和模型2，可以发现：随着结构计算区域的增大，建模离散的SPH粒子数量大幅增加。无论是SPH算法还是有限元法，其计算效率都与计算粒子或单元个数直接相关。对比实际计算效率，若采用内存为6 GB的个人计算机进行计算，则采用纯SPH算法时模型1和模型2的计算时间分别为1.8和2.6 h，而采用SPH-FEM耦合算法时模型1和模型2的计算时间分别为2.3和2.5 h。

图  10  分离结构

Figure  10.  Separation structure

下载: 全尺寸图片 幻灯片

由此可知：当结构计算区域较小时，本研究中的SPH-FEM耦合算法在计算效率上的优势并未体现，纯SPH算法的计算时间更少；但是，随着结构计算区域的增大，SPH粒子数较FEM明显增加，相应的计算时间也大幅增加，SPH-FEM耦合算法的计算效率优势逐步显现。结构计算区域对实际计算耗时有重要影响，结构复杂程度则直接影响建模难易程度及耗时情况。在实际工程中，结构往往呈大型化、复杂化、系统化，因此建模的可实现性及便捷化也是一个重要的考虑因素。以图 10所示的结构为例，若将结构的计算区域再扩大，即包含连接螺栓和装置盒等结构，则采用SPH建模并不容易实现；对于飞行器等大型复杂结构而言，全部采用SPH粒子建模几乎是不可能实现的。所以对于大型复杂结构模型，采用SPH-FEM耦合算法不仅可以解决单纯SPH建模的局限性问题，而且有效提高模型分析的计算效率。

5.   结论

基于SPH理论，在总结前人工作的基础上，采用SPH-FEM耦合算法对具有复杂结构的复合材料蒙皮结构飞行器的爆炸分离过程进行三维数值模拟，计算分析其螺栓爆炸冲击环境及飞行器结构的冲击动力学特性，初步实现了SPH-FEM耦合算法的工程领域应用，并得到如下结论。

(1) 采用少量的SPH粒子便可成功预测炸药爆轰载荷，而结构的主体响应采用FEM进行计算，最终SPH-FEM的计算效率主要由FEM的计算效率决定。

(2) 将复合链表搜索法运用到SPH-FEM耦合算法的计算过程中，在一定程度上提高了三维复杂结构的计算及分析效率，由此实现了飞行器舱段结构的螺栓爆炸解锁分离的成功模拟；通过对比数值模拟结果与实际爆炸分离试验结果，验证了SPH-FEM耦合算法及相应数值模型的有效性和准确性。

(3) 螺栓爆炸冲击激励通过耦合界面传递到飞行器结构中，螺栓孔附近区域始终维持着较大的应力响应，结构中的冲击加速度均在20 ms内基本完成主要衰减。

(4) 在螺栓爆炸分离过程中，飞行器后段区域的响应较剧烈，冲击波通过分离结构后板与复合材料蒙皮结构之间的结构接合部时发生较大程度的衰减。在6个爆炸螺栓对称安装连接及后板中装置盒的特定分布模型中，爆炸冲击幅值衰减可达70%，之后继续向前段传播，离螺栓爆炸激励源越远的结构受到的冲击响应越小。

(5) 对于大型复杂结构模型，SPH-FEM耦合算法不仅解决了单纯SPH建模的局限性问题，还有效提高了计算效率。

研究结果可为航天飞行器结构设计与防护提供一定的参考。