基于 PD-SPH 方法的细长体破冰及入水过程中冰损伤机理研究

《International Journal of Impact Engineering》：Study on the ice damage mechanism during the ice-breaking and water entry of a slender body based on the PD-SPH methods

　　

摘要 随着全球变暖加剧以及极地地区战略地位的提升，细长体高速跨介质破冰已成为极地资源开发和军事防御的关键研究方向。该过程涉及冰、水和固体介质之间的强非线性相互作用以及瞬态冰层破裂，呈现出复杂的力学机制。因此，全面理解穿透过程中极地冰层的损伤机制对于设计跨介质载具至关重要，这包括理解内部裂纹萌生、扩展网络演化以及宏观破碎模式。本研究基于耦合近场动力学（PD）与光滑粒子流体动力学（SPH）的数值方法，辅以物理模型试验，对这些损伤机制进行了系统研究。首先，利用改进的分裂霍普金森压杆对具有不同角度预裂纹的冰试件进行动态压缩试验，分析了缺陷角度对应力响应和裂纹特征的影响。随后，基于非普通状态基近场动力学（NOSB-PD）理论，结合Drucker-Prager屈服准则和最大主应力断裂准则，建立了冰的弹塑性损伤本构模型。之后，构建了一个双向PD-SPH流体-结构相互作用数值模型，以模拟整个破冰和入水过程。该模型引入了一种基于声频黎曼求解器的弱可压缩SPH公式，通过原位虚拟粒子映射和改进的单边黎曼边界条件实现高效界面耦合。最后，将该耦合框架应用于细长体穿透平冰的大规模数值模拟，系统研究了径向和环向裂纹的形成机制及其对冲击条件的参数依赖性。结果表明，当冲击速度从50 m/s增加到200 m/s时，失效模式从纯径向断裂转变为密集的"蛛网状"网络。同时，冰层厚度增加增强了局部弯曲抗力，阻止了环向裂纹的萌生和径向扩展。此外，弹头形状显著影响损伤机制。锥形弹头通过点接触引发径向裂纹，半球形弹头导致整体弯曲形成复合网络，而平顶弹头则引发严重的瞬态冲剪失效，伴随密集的径向和远端环向裂纹。

引言 随着极地地区战略意义的日益提升，传统的破冰方法已不能满足深度探索和隐蔽侦察的需求。依赖巨大瞬时动能的细长体高速跨介质破冰入水已成为一种创新且高效的方法。然而，该过程涉及"冰-水-弹体"三相介质在极短时间内的强烈非线性耦合，伴随严重的冰层断裂、大规模水体变形和空化演化。鉴于传统物理实验面临的挑战和高成本，开发高保真数值模拟方法来揭示多相耦合下的动态失效机制已成为该领域的重要科学难题。

冰作为一种天然准脆性多晶材料[1]，其宏观失效模式对多种因素高度敏感，包括应变率[2]、环境温度[3]、盐度、内部微观结构[4]和缺陷[5]。弹体穿透作为一种有效的破冰方法，受到了学术界的广泛关注。现有的基于弹体的破冰研究主要集中在水下无冰环境中，系统揭示了穿透参数对冰失效模式的影响[[6],[7],[8],[9],[10]]。Zhang等人[6,7]通过实验和基于键的近场动力学（BBPD）方法阐明了冲击下的动态断裂和损伤演化机制。Galadima等人将速率相关损伤本构模型引入非普通状态基近场动力学（NOSB-PD）[8]。该模型考虑了冰的应变率敏感性，进一步揭示了阻尼比和冲击速度如何影响其失效特征。然而，无冰环境中的力和失效机制与真实的冰水耦合环境存在显著差异。针对破冰入水问题，Wang等人利用LS-DYNA研究了弹头形状和水体支撑对冰盖冲击碎片的影响[11]。Liu等人通过将BBPD与更新拉格朗日粒子流体动力学方法耦合，模拟了刚性球体冲击漂浮冰板时的裂纹扩展[12]。然而，该模型未能模拟冰层破碎后球体的入水过程，且其裂纹萌生顺序与实验结果不符。Ren和Zhao提出了一种新型模型，耦合BBPD、流体体积法（VOF）和浸没边界法，模拟了刚性球体破冰入水的整个过程[13]。然而，他们的冰板表现出非真实的、类似塑性的大局部变形。最近，Zhang等人采用耦合欧拉-拉格朗日（CEL）方法研究了细长体破冰入水过程中的动态特性，确定了初始冲击、持续穿透、入水拍击和稳定水动力运动四个特征阶段[14]。

在解决复杂裂纹萌生和扩展问题时，传统的基于网格的连续介质方法，如有限元法（FEM）[[15],[16],[17],[18]]，经常因网格畸变而面临挑战[19]。相比之下，相场模型（PFM）已成为模拟复杂断裂过程的有效方法，通过标量相场演化来模拟裂纹萌生、扩展、汇合和分叉，无需显式裂纹追踪[20]。同时，无网格离散元方法（DEM）[[21],[22],[23],[24]]在连续变形阶段通常存在计算精度不足的问题。与之相比，Silling[25]提出的近场动力学（PD）通过积分方程重建控制方程，能够自发捕获裂纹演化。Dorduncu等人[26]对PD理论及其各种公式——包括基于键的（BB）、普通状态基（OSB）和非普通状态基（NOSB）PD——以及非局部算子及其计算机实现进行了全面综述。除了标准PD公式外，非局部算子方法与显式相场模型的结合也已开发用于动态脆性断裂模拟，提供了替代的非局部离散化策略[27]。早期的BBPD已广泛应用于船舶/海洋结构-冰相互作用[28,29]、水下爆炸[30]、螺旋桨-冰铣削[31]和气泡诱导破冰[32]等问题。然而，BBPD受到固定泊孙比的限制和理想脆性假设的约束[33]。虽然最近的努力已将塑性修正引入BBPD用于延性材料的疲劳分析[34]，但对于冰等准脆性材料的复杂弹塑性损伤建模，状态基PD公式仍然更为适用。PD也已扩展到功能梯度材料和双材料结构，其中近场动力学微分算子可以自然地适应空间性质变化和材料不连续性[35,36]。为了克服这些局限性并开发更复杂的冰材料冲击损伤本构模型，Silling进一步提出了状态基近场动力学（SBPD）[37,38]。SBPD分为普通状态基近场动力学（OSBPD）和非普通状态基近场动力学（NOSB-PD）。OSBPD克服了固定泊松比的约束，已被证明在模拟冰结构相互作用方面具有高度适用性[39]。同时，通过引入非局部变形梯度，NOSB-PD可以直接引入经典连续介质力学本构模型[40]。在NOSB-PD基础上，Song等人应用Drucker-Prager准则分析了球-冰板冲击，但由于忽略了动态力学参数，未能产生明显的裂纹[41]。Galadima等人通过引入速率相关损伤模型成功模拟了裂纹，尽管产生的裂纹形态相对无序[8]。目前，NOSB-PD在精确控制裂纹形态和准确反映冰的动态损伤演化方面仍存在一定的局限性。

对于涉及大自由面变形的流体动力学问题，光滑粒子流体动力学（SPH）[42,43]作为一种拉格朗日无网格方法，具有显著优势。与传统的基于网格的计算流体动力学方法相比，SPH将连续体分解为粒子。这些粒子携带质量、速度和压力等物理性质。SPH通过核函数插值计算场变量及其导数。其拉格朗日方法避免了由网格畸变和纠缠引起的数值不稳定性。这在模拟大变形自由面流动方面具有明显优势[[44],[45],[46]]。根据流体可压缩性的处理策略，SPH方法分为不可压缩SPH（ISPH）和弱可压缩SPH（WCSPH）。ISPH严格遵循流体的不可压缩性约束，通过求解压力泊松方程（PPE）获得压力场。这提供了光滑且准确的压力场，但需要更多的计算资源。相比之下，WCSPH通过状态方程（EOS）将压力与密度关联，计算过程更简单，并能自动跟踪流体界面[46]。然而，SPH核函数的近似离散特性在模拟高雷诺数或无黏性流体运动时可能导致非物理的数值振荡。为了解决这个问题，研究人员提出了人工阻尼、密度控制和黎曼求解器集成等修正策略。在动量方程中添加人工粘性项是一种广泛使用的平滑方法[[46],[47],[48]]。然而，参数设置过高可能导致动能异常下降和结果不准确。Antuono和Marrone通过在连续性方程中引入密度耗散项提出了δ-SPH模型[49,50]。随后，Sun等人[51]和Zhang等人[52]将此方法应用于精确复制楔形和圆柱形结构入水和出水过程中压力分布、非线性爬升、射流形成和流体分离等现象。密度过滤通过定期重置粒子密度来消除波动，但如果过于频繁地应用，可能会违反质量守恒并引入额外的耗散。为了消除对经验人工参数的依赖，Inutsuka提出了一种基于黎曼求解器的SPH方法，通过求解界面状态隐式引入必要的耗散[53]。在此基础上，Zhang等人于2017年提出了一种低耗散黎曼-SPH方法[54]。该方法使用自适应限幅器动态去除过量的能量损失，在经典测试场景中展现出优异的机械能守恒和光滑压力场特性。与此同时，SPH固体边界处理技术也在不断改进。从排斥力模型[46]和动态边界条件[55]发展到低耗散黎曼框架内的优化策略[54]，大大提高了边界计算精度。总之，SPH及其改进算法已成为模拟强非线性流体-结构相互作用的可信工具。

鉴于SPH在模拟大变形自由面流体方面的有效性和PD在模拟不连续介质断裂方面的优势，开发无网格PD-SPH耦合流体-结构相互作用模型具有重要意义。Li等人通过定义界面区实现两种算法之间的相互作用，模拟了爆炸冲击下的土壤碎裂[[56],[57],[58],[59]]。2021年，Zhou等人尝试了光滑近场动力学（SPD）方法[60]，在控制方程层面实现了两种方法的深度整合。SPD在SPH粒子的支撑域内引入PD的"键"连接机制，旨在结合SPH的流体建模能力与PD的断裂捕获优势。从2020年到2022年，Sun提出了一种基于移动虚拟粒子的PD-SPH耦合模型[61]。该模型在流体-固体边界粒子之间建立"伙伴关系"，利用虚拟粒子进行机械数据传输，并引入排斥力保持界面几何稳定。最近的研究在此基础上进一步推进了该方法的工程应用，并深入研究了复杂流体-结构相互作用中的物理失效机制。Liang等人创建了一个综合的三维数值冰池，模拟了垂直圆柱形结构连续破冰的典型场景[62]。针对高速水射流破冰技术，Shi等人通过数值模拟精确再现了水射流冲击海冰的整个动态过程，从能量和应力波传播角度揭示了裂纹萌生和扩展控制机制[63]。尽管PD-SPH已被应用于多个冰工程问题，但尚无学者采用该方法研究弹体穿透破冰入水过程中复杂的三相耦合场景。此外，现有用于冰工程的PD-SPH框架很少集成完整的冰弹塑性速率相关损伤本构模型，在处理高速强冲击流体-结构耦合问题时仍存在数值稳定性和界面精度方面的局限性。

本研究结合物理实验和数值模拟，揭示了细长体高速破冰入水过程中的冰损伤机制，其贡献明确区分为方法适配和工程应用两个方面。在方法方面，建立并专门适配了一个整合NOSB-PD和SPH的耦合数值框架，用于冰-水-固体三相冲击断裂场景。具体来说，将Drucker-Prager屈服准则、速率相关抗压强度关系和最大主应力键断裂准则纳入NOSB-PD框架，以表征冰的全过程弹塑性动态损伤行为。

首先进行了分裂式霍普金森压杆（SHPB）动态冲击试验，并利用NOSB-PD方法建立了动态冰冲击数值模型。通过将数值模拟结果与试验数据进行对比，验证了所开发的基于NOSB-PD的冰损伤模型在捕捉冰动态断裂特性方面的有效性。随后，将SPH公式与经过验证的NOSB-PD框架耦合，建立了耦合PD-SPH框架，用于研究高速穿透过程中冰-水-弹体的相互作用。对于流体域，采用了一种基于黎曼的低耗散弱可压缩SPH方案，并开发了一种改进的虚拟粒子映射算法，以实现精确的双向流固耦合。同时引入了接触惩罚算法、数值阻尼和沙漏控制，以确保强冲击载荷下数值计算的稳定性。

在应用方面，通过刚体球体破冰入水的物理试验，进一步验证了该耦合模型的可靠性。在此基础上，本研究进一步将PD-SPH耦合框架扩展到高速尖体破冰及入水的大规模数值应用，这代表了一种利用PD-SPH方法很少被研究过的复杂三相耦合场景。大规模模拟被用作参数化数值研究，以探究径向和周向裂缝的形成机制，并评估冲击速度、冰层厚度和弹体头部形状对冰层宏观断裂响应的影响。

**截面片段：NOSB-PD应用于冰**

如图1所示，PD中的连续体被离散为具有质量密度$\rho_i$和增量体积$dV_i$的材料点$X_i$。这些点之间的非局部相互作用通过"键"在半径为$\delta$的近场动力学视界$H_i$内发生。初始构型中的形状张量$K_i$定义为：

$$K_i = \int_{H_i} \bar{\omega}\langle\xi_{ij}\rangle \xi_{ij} \otimes \xi_{ij} \, dV_j$$

其中$\overline{\langle\xi_{ij}\rangle} = \xi_{ij}$是表示初始键几何向量的位置矢量状态，$\bar{\omega}\langle\xi_{ij}\rangle$是描述键强度的影响函数，其Validation of the PD-SPH coupled model译为"PD-SPH耦合模型的验证"，作者从之前开展的SHPB冲击试验中选取了一个代表性算例进行数值复现[64]，以验证所提出的NOSB-PD框架和冰塑性-损伤本构模型。

图4展示了含有45°预存缺陷的冰试件在动态SHPB加载下的整个损伤演化过程。上行呈现了物理试验的高速摄影图像，下行则展示了相应的NOSB-PD模拟结果。

**结论**

本研究采用耦合的PD-SPH数值框架和物理模型试验，系统研究了尖体高速破冰及入水过程中冰的动态损伤机制。基于NOSB-PD理论，结合Drucker-Prager屈服准则和最大主应力断裂准则，建立了冰的弹塑性损伤本构模型。此外，利用低耗散方法建立了双向PD-SPH流固耦合模型。

**CRediT作者贡献声明**

Yiheng Zhang：撰写—审阅与编辑，监督，项目管理，经费获取，概念构思。Ziyang Liu：撰写—初稿，方法学，研究。Chen Chi：撰写—初稿，方法学，研究。Ying Tang：方法学。Yiquan Yu：研究。Xiaobin Li：监督。

武汉理工大学 船舶冲击与防护湖北省重点学科实验室，中国武汉，430063