在现代军事斗争领域，高速动能弹对各种地下结构（如：地下储藏库、地下隧道和地下军事指挥所等）具有很大的威胁，弹体侵彻问题一直是热点研究课题。岩体材料作为地下空间结构的天然防护材料，对抵抗杀伤武器有良好的防护能力，在各种防护工程中具有广泛的应用（沈俊等，2018）。另一方面，冲击破岩技术广泛应用于矿山开采与石油开采领域。近年来，为了提高硬岩中的钻井速度，一项新的钻井技术——粒子冲击钻井技术逐渐发展起来。该技术通过球形钢粒子对岩体的高速侵彻冲击实现破岩，有效提高了硬岩中的钻井效率，具有广阔的应用前景（况雨春等，2012；赵健等，2013）。花岗岩作为一种常见的岩石材料，研究其在弹体高速撞击下的力学损伤响应以及成坑规律对于防护工程设计和钻井破岩技术的发展具有重要意义。

弹体和靶体目标之间存在一种复杂的相互作用，涉及到大应变、非线性问题，侵彻深度的预测是研究的重点，而根据大量试验数据总结出的经验公式是过去预测弹体侵彻深度的主要方法。在试验研究方面，早在20世纪50年代，美军就对几类航弹撞击花岗岩和砂岩时的过程进行了系统的试验研究（Livingston et al.，1951）。Forrestal et al.（1994）进行了卵形头弹对混凝土靶体的一系列正侵彻试验，得到一系列弹体撞击速度与侵彻深度之间的关系，总结出了关于混凝土靶体侵彻深度的经验公式，该公式中侵彻深度是关于无量纲经验常数S的函数，S的值取决于混凝土靶体的无侧限抗压强度。Frew et al.（2000）将以上经验公式拓展到石灰岩材料靶体，计算结果与试验数据对比吻合度良好。张德志等（2005，2006）利用压缩气体炮开展了一系列卵形钢弹侵彻高强度花岗岩靶体的试验研究，提出了关于侵彻深度的经验公式。然而，这类利用试验数据拟合得到的经验公式（Livingston et al.，1951；Forrestal et al.，1994；Frew et al.，2000；张德志等，2005；Young，1967，1997；Berard，1975，1977，1978）均是基于各自的试验数据总结出来的，大多没有广泛的适用性，在实际应用中难以保证高准确性。在理论研究方面，空腔膨胀理论应用最广泛。研究者们针对不同的侵彻介质（如土壤、混凝土和岩体等），运用空腔膨胀理论进行了大量的研究（Forrestal et al，1992；王明洋等，2003，2005；邓勇军等，2020）。

随着现代数值计算能力的发展，数值计算方法成为解决岩体冲击问题的重要方法。光滑粒子流体动力学（SPH）方法作为无网格方法中应用最广泛的计算方法，既可以跟踪物质的运动，又不存在网格畸变严重而终止计算的问题，很适合解决高速冲击等问题。但SPH方法也存在着一定的缺陷，如：边界条件难以设置以及计算资源消耗过大等。基于此，本文采用显式动力学分析软件LS-DYNA，应用SPH-有限元（FEM）耦合方法，采用经SHPB试验验证的白麻花岗岩材料参数，开展了白麻花岗岩在刚性弹体高速侵彻条件下的损伤响应及侵彻深度预测研究。

本研究中的岩石材料为产自湖北随州的白麻花岗岩，经取芯加工后制成直径和高度均约为50 mm的圆柱体试样，如图1所示。经室内试验得到白麻花岗岩的相关物理静力学参数主要为：密度ρ=2 607 kg/m³，单轴抗压强度f_c=89.4 MPa，单轴抗拉强度f_t=4.56 MPa，弹性模量E=24.1 GPa，泊松比λ=0.22。

HJC材料模型考虑了高静水压力、高应变率和材料的损伤效应，能够应用于岩石在高应变率下的大变形问题研究，能够较好地描述岩石在高速侵彻下的力学行为响应。HJC本构模型主要包含3个部分，即屈服面方程、损伤定义和状态方程，具体的模型描述可参见文献Holmquist et al.（1993）。在室内测试得到的白麻花岗岩物理力学基本参数的基础上，结合SHPB试验测得的岩石动态压缩强度，通过理论计算和参数研究确定了岩石的HJC材料模型参数，列于表1。

表2为利用MTS和SHPB试验系统测得的白麻花岗岩在不同应变率下的压缩强度，试样强度的提高不仅受应变率的影响，而且受静水压力的影响，采用文献（Holmquist et al.，1993）的方法，消除静水压力的影响，计算得到了HJC本构材料中的应变率效应参数C=0.0036。

根据动态空腔膨胀理论，在一定的速度范围内，侵彻弹体可以视为刚体。本研究中，假设高强钢弹体材料为刚体材料，质量密度ρ₀=7 850 kg，弹性模量E=206 GPa，泊松比λ=0.3。

为了验证上述本构材料模型在高应变率下岩石材料动力响应预测中的有效性，将该模型用于模拟2种应变率下的花岗岩试件的SHPB试验，将模拟所得的岩体应力—应变规律与实际室内试验中的岩体真实应力—应变规律进行比较。

入射杆和透射杆的长度均为2 000 mm，直径为50 mm。2种杆件均采用Cr40钢制作，并采用弹性材料模型进行模拟。钢体的主要材料参数：杨氏模量为240 GPa，质量密度为7 810 kg/m³，泊松比为0.285。对尺寸为50 mm×50 mm（直径×长度）的岩样进行建模。数值模拟过程中使用表1中的HJC模型参数来模拟岩石材料。图2为SHPB试验模型的局部示意图，杆件和岩样均采用三维实体单元建模。通过使用关键字*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE，将试件与入射杆和透射杆之间的接触设置为光滑接触。同时由于数值模拟的重点是岩石的响应行为，因此将实际室内试验中记录在入射杆上的应力脉冲作为入射杆的应力边界施加在入射杆端，而不是对纺锤形冲头进行建模，由此可以合理简化数值模拟。

图3为室内SHPB试验中典型加载条件下入射杆和投射杆上应变片所记录到的入射波、反射波和透射波的应变信号以及对应的试样应力平衡图。由图3可以看出，当达到应力峰值时，试样已经达到了应力平衡状态，说明试验数据有效。

将室内试验入射杆上记录到的入射应力波作为边界条件施加到图2所示的数值模型入射杆端，得到白麻花岗岩试件应力—应变规律。图4所示为花岗岩试件在140 s^-1和155 s^-1应变率下的真实应力—应变曲线和数值模拟的岩石应力—应变曲线对比。结果表明，采用上述HJC模型常数的数值模拟对岩样的行为预测较合理。如图4所示，在应变率为140 s^-1的条件下，应力—应变曲线的峰值变化趋势基本一致，SHPB试验中测得的真实抗压强度为137 MPa，而数值模拟中白麻花岗岩试件在同样应变率下的动态抗压强度为141 MPa，比实际的岩石强度稍高。这种差别可部分归因于钢杆与岩样之间的接触摩擦力，这在数值模拟中一般不予考虑。由室内试验和数值模拟的应力—应变曲线对比可以看出，采用上述HJC材料模型能够较合理地预测白麻花岗岩高应变率下的力学行为。

侵彻过程伴随着网格的大变形及高应变率，传统的侵彻过程数值模拟通常采用基于纯网格的有限元方法（FEM），如Lagrange算法。通常来说，Lagrange算法通过关键字添加设置材料侵蚀失效准则来删除畸变网格以保证计算结果的收敛。但由于很难得到一个通用的材料失效准则，该方法的数值模拟结果的可信度难以保证。为了解决这类材料的大变形问题，作为一种基于拉格朗日运动方程的无网格方法，光滑粒子流体动力学（Smoothed Particle Hydrodynamics，SPH）方法应运而生，且应用较为广泛。

实现SPH方法的2个关键步骤是核估计与粒子近似（Liu et al.，2005）。第一步是核估计，粒子场函数f（x）的核估计表达式为

式中：角括号代表核近似运算；W(x-x',h)为光滑核函数；h为定义核函数W影响区域的光滑长度；x为粒子的位置。在3次样条函数基础上提出的B-样条函数（Liu et al.，2005）是研究者最常用的核函数，其形式为

式中：R=|x-x'|/h，为粒子间距。由式（2）可知，核函数的影响域为2h，超过2h范围的粒子不再对所求粒子产生影响。第二步是粒子近似，函数f（x）以及其导数∇f（x）的核估计可以转换为问题域Ω中的有限个粒子加权求和的形式。对于给定的粒子i，应用粒子近似法可得到粒子i处的函数及其导数近似式如下：

式中：N为粒子i周围支持域内的粒子总数；ρ为密度。

连续介质条件下，流体动力学基本方程的质量守恒、动量守恒和能量守恒方程的SPH离散表达式如下：

式中：σ为应力；v为粒子速度；α和β为空间坐标方向；∏ij为人工黏度，用来减少非物理数值震荡；e为内能。

弹体网格采用拉格朗日法计算，为了减少计算量并节省计算资源，假设弹体为刚体。弹体直径为20 mm；长径比为6；弹头形状系数（CRH），即弹头圆弧曲率半径和截面直径的比值为3。SPH方法存在一些局限性，例如计算效率低、边界条件难以设置等，如果将整个靶体设置成SPH粒子，计算会非常耗时。基于此，本研究中只将弹体接触的靶体中心区域采用SPH粒子建模，周围岩体采用拉格朗日网格建模。以SPH粒子作为从节点，以对应的拉格朗日网格作为主面，通过关键字*CONTACT_TIED_NODES_TO_SURFACE_OFFSET将SPH粒子与朗格朗日网格耦合，即采用罚函数算法处理接触界面。值得一提的是，本数值模型中接触界面上的SPH粒子与拉格朗日网格保持一一对应。白麻花岗岩靶体总体尺寸为50 cm×50 cm×15 cm，中心的20 cm×20 cm×15 cm区域设置为SPH粒子。靶体四周设置为固定边界，其余2个面为自由面，对称面上施加对称位移约束。同时，弹体与靶体粒子之间采用点面侵蚀接触。整体网格模型的1/2模型如图5所示。

与爆炸类似，中高速冲击引起的地冲击在一定范围内也可以采用式（6）描述（王明洋等，2018）：

式中：P_max为应力波峰值；d为冲击点到测点的距离；n为衰减指数。衰减系数与岩石压力状态相关。哈努卡耶夫（1980）总结了衰减指数与侧压力系数α之间的简单关系，表达式为

式中：正号对应高应力流体状态；负号对应低应力固体区；在高应力流体区，α值约为1；在低应力固体区，α的取值范围为0.2~0.6。

为了验证拉格朗日网格与SPH粒子之间接触的有效性，在靶体上表面距离侵彻中心一定距离的位置，间隔20 mm布置一系列测点，具体测点位置如图6所示，O点为侵彻中心点。

将撞击速度为250 m/s的情况下各测点的压力峰值以幂函数拟合，地冲击衰减规律如图7所示，得到：

式中：P_max为峰值压力（MPa）；d为距离（cm）。

图7展示了初速度为250 m/s的情况下靶体历史压力云图及各测点处的峰值压力衰减规律。由图7可以看出，应力由SPH粒子传导至拉格朗日有限元网格，应力峰值连续且衰减规律符合幂函数，所取测点区域为低应力固体区，拟合得到的衰减系数1.74符合文献（哈努卡耶夫，1980）中总结的衰减系数。这说明本研究中SPH-FEM耦合方法的SPH粒子与拉格朗日网格的接触设置有效，数值模拟结果具有较高可信度。

当正侵彻初速度为250 m/s 时，白麻花岗岩靶体的历史损伤云图如图8所示。其中，损伤值介于0~1之间，0表示岩体完全无损伤，1表示岩体完全破坏。由图8可以看出，SPH-FEM耦合算法可以很好地模拟出侵彻粒子飞溅的效果，损伤云图直观地展现了侵彻损伤的区域以及刚体弹正侵彻形成的前弹坑，靶体损伤区域较为集中。当侵彻初速度为500 m/s时，弹体穿透了白麻花岗岩靶板，如图9所示。由于撞击速度更大，在靶体前自由面上形成的损伤面积相比于初速度为250 m/s时更为明显，同时靶体前弹坑和后弹坑都能较好地展现。

张德志等（2005）开展了一系列2种直径的高强钢弹对花岗岩正侵彻的试验研究，弹体初速度为300~900 m/s，得到了相应的侵彻深度。基于Berard（1977）提出的无量纲侵彻公式，提出了如下经验公式：

式中：ρ为靶体密度；H为侵彻深度； B和V₀为待定系数；M为射弹质量；A为弹体截面面积；Y为岩体无侧限抗压强度；RQD为岩石质量指标。由式（9）可以看出，在相同侵彻条件下，侵彻深度与侵彻速度呈线性正相关，而与岩体无侧限抗压强度呈负相关。

本研究中，进行了一系列不同初速度的弹体侵彻数值模拟，将所得到的弹体侵彻深度与式（9）拟合，取RQD=90，可以确定待定系数的值：B=0.238，V₀=13.85 m/s。图10为数值模拟结果与经验公式的比较，可以看出，数值模拟结果与式（9）吻合度较好。该经验公式可以为强度相近的岩体侵彻深度预测提供参考，但其适用性还需要更多的现场侵彻试验加以验证。

为了研究弹头形状对弹体侵彻白麻花岗岩靶体的影响，设置一系列平头弹以不同的初速度侵彻岩石靶体，其中，岩体材料、弹体材料参数、弹体总质量与之前试验保持一致。得到一系列侵彻深度如图11所示。由图11可以看出，相同初动能的平头弹比卵形弹侵彻深度要小。初速度为50，100，150，200，250，300 m/s的平头弹侵彻深度分别为卵形弹侵彻深度的16.7%、27.8%、35.1%、32.1%、36.1%和40.5%，在低速撞击情况下平头弹与卵形弹的侵彻性能差别更为明显。相比于平头弹，卵形弹的侵彻过程是一种穿刺式的过程，而平头弹的侵彻过程是一种挤凿式的过程，这造成了二者侵彻性能的差别。图12为初速度为250 m/s的卵形弹与平头弹对靶体的最终侵彻损伤云图对比。由图12可以看出，相比卵形弹，平头弹在靠近入射自由面处的损伤区域更大，即产生的沙漏形状弹坑更大，但侵彻深度小得多，总体侵彻损伤区域更小。

基于经SHPB试验验证的HJC材料参数，建立了一系列刚性弹侵彻白麻花岗岩靶体模型，采用显式动力学分析软件LS-DYNA结合SPH-FEM耦合方法分析了白麻花岗岩靶体在中高速侵彻冲击下的损伤过程，得到以下主要结论：

（1）通过比较分析室内SHPB试验与数值模拟中得到的试样应力—应变规律，验证了HJC本构材料模型能够合理地描述白麻花岗岩在高应变率下的力学行为。

（2）SPH-FEM耦合方法能够有效地模拟出靶体在中高速侵彻作用下的力学损伤响应，应力在SPH粒子与有限元网格接触界面上传导连续，靶体在高速撞击下的地冲击规律与爆炸类似，符合经验公式预测结果。

（3）由数值模拟得到的侵彻深度数据，建立了白麻花岗岩侵彻深度预测的经验公式，将该公式与前人建立的经验公式进行了对比，结果表明该公式结果与数值模拟结果一致性良好，可用于相近强度的岩石侵彻深度预测。

（4）在不同侵彻初速度情况下，平头弹的侵彻深度以及对靶体的损伤区域均比相同动能下的卵形弹小，且在低速撞击情况下差别更为明显，卵形弹对于靶体的侵彻性能优于平头弹。