一种用于考虑准脆性断裂的水弹性流固耦合的EISPH-TLSPH-DM耦合方法

《Applied Mathematical Modelling》：A coupled EISPH-TLSPH-DM method for hydroelastic fluid-structure interaction with quasi-brittle fracture

　　

•开发了一种基于连续介质的统一SPH框架，用于水弹流固耦合（FSI）分析。
•在该统一框架中捕获了大变形和准脆性断裂行为。
•梯度增强的非局部损伤模型消除了裂纹尖端处的导数奇异性。
•在强耦合流固耦合条件下，应变局部化与裂纹演化得到了有效解析。
•准静态与动态流固耦合的预测结果与实验和仿真结果吻合。

引言

涉及动态结构断裂的流固耦合（FSI）精确模拟在许多工程应用中具有关键意义。然而，高度非线性流场和严重材料不连续性的存在往往带来重大的数值挑战。作为一种成熟的无网格方法，光滑粒子流体动力学（SPH）将连续介质离散化为粒子，并通过邻近粒子间的核函数相互作用来求解质量、动量和能量守恒的控制偏微分方程（PDEs）[1]。由于质量和自由度的共置，SPH提供了一种相对简单的方式处理固液相变、结构损伤与断裂以及强耦合流固耦合问题。

根据不同的压力公式化方案来封闭控制方程，已开发了多种基于SPH的流体求解器。在弱可压缩SPH（WCSPH）中，假设流体具有轻微可压缩性，通过状态方程计算流体压力，这会导致由于人工声速而引起的较小时间步限制[2]。相比之下，隐式不可压缩SPH（IISPH）采用一阶投影法，其中粘性力和体力显式求解，而不可压缩条件通过迭代求解压力泊松方程（PPE）来施加[3]。该方法对时间步的CFL约束较宽松，避免了WCSPH中常见的压力噪声问题。显式不可压缩SPH（EISPH）通过使用流体压力的显式Jacobi型求解方法提高了计算效率，消除了全局矩阵组装和迭代求解器，代价是精度的轻微损失[[4], [5], [6], [7]]。

在基于SPH的流固耦合模拟中，流体域和结构域之间采用两种耦合策略：混合方法和统一方法。在混合FSI方法中，流体域使用基于SPH的求解器来发挥其在流体建模方面的优势，而结构域则使用独立的非SPH求解器处理。相比之下，统一方法采用完全基于SPH的描述来描述耦合的流体-结构系统。更新拉格朗日SPH（ULSPH）[8]相对于不断演化的当前构型评估粒子相互作用，全拉格朗日SPH（TLSPH）[9,10]则相对于初始参考构型进行评估，二者被广泛用于结构动力学，并在水弹流固耦合问题中展现了显著的精度和鲁棒性[11]。

尽管SPH具有巨大潜力，但在SPH中实现精确的断裂建模并非易事，因为裂纹表面的出现引入了导数奇异性，破坏了核函数的一致性。这导致导数评估中误差放大，裂纹面附近粒子分布畸变，往往导致变形梯度张量的雅可比行列式为负值。因此，混合FSI方法通过采用替代性的非连续介质结构表征，绕开了这些病态微分问题。该方法中使用的显著结构求解器包括有限粒子法（FPM）、离散元法（DEM）、格点粒子法（LPM）和近场动力学（PD）。这些方法的共同特征是结构响应通过局部粒子相互作用而非基于连续介质的微分方程来描述。

基于向量力学的FPM中[12]，粒子运动通过求和满足牛顿第二定律的局部力矢量来推进，确保内节点力与外力之间始终处于动态平衡[13,14]。该公式不需要组装全局刚度矩阵或迭代求解过程，使得引入额外粒子来表示断裂过程中的材料分离或碎裂变得简便[15]。

在DEM中[16,17]，固体的力学行为通过邻近粒子之间的轴向类弹簧连接来建模。弹簧刚度通过离散粒子系统与其相应连续表征之间的能量等效来确定。材料分离自然地通过键断裂规则描述，一旦满足失效准则，力学连接即被移除。LPM[18]被认为是DEM的扩展，允许使用任意泊松比。除了最近邻成对局部键伸长外，它还纳入了非局部体积势能，缓解了与初始粒子堆积构型相关的裂纹路径偏置。弹性模型、屈服面或粒子堆积的任何变化都需要相应地重新公式化键力。为提高分辨率收敛性，最近的LPM公式采用了稳健的断裂度量，包括与材料断裂韧性对比的逐粒子应力强度因子[19]，或由定义为粒子尺寸、材料刚度和断裂能量的演化定律控制的损伤相关变量[20]。

近场动力学（PD）是一种非局部粒子公式，用积分相互作用取代空间导数，使其适用于涉及不连续性的断裂问题。基于键的近场动力学（BBPD）采用成对键力，但限制了可实现泊松比的范围[21]。普通态基近场动力学（OSPD）在非局部框架内允许具有任意泊松比的一般弹性本构行为[22]，而非普通态基近场动力学（NOSPD）进一步消除了键方向约束，实现了与基于连续介质的弹塑性本构模型的完全对应[23,24]。由于上述优势，这些混合方法在基于SPH的流固耦合框架中获得了广泛应用，包括WCSPH-FPM [15]、WCSPH-DEM [25]、WCSPH-LPM [19,20]、WCSPH-BBPD [26,27]、WCSPH-OSPD [[28], [29], [30]]、WCSPH-NOSPD [31,32]以及EISPH-OSPD [7]。

尽管混合流固耦合方法已取得相当大的成功，但统一的完全基于SPH的方法具有若干实际优势。首先，流体和结构求解器共享许多算法特征，如PDE离散化过程和最近邻粒子搜索，便于简单实现并行计算。其次，守恒一致的界面模型更容易构建，便于在子域之间传递质量、动量、能量和化学物种。这一特征对多物理场问题尤为有利。第三，许多实际工程模型，特别是涉及弹塑性或弹塑性-蠕变断裂行为的模型，都是在有限元（FE）或有限体积（FV）框架内开发的。为了使这些成熟的模型适用于SPH框架，采用连续介质一致的SPH断裂公式可能更为理想。

据作者所知，目前提出过少数包含水弹脆性断裂的全SPH流固耦合框架，主要包括伪弹簧法（SPH-PS）和相场法（SPH-PF）。SPH-PS [[33], [34], [35]]将基于核的粒子相互作用视为虚拟弹簧连接，一旦满足预设准则便移除虚拟连接来建模宏观断裂，方法概念简单且易于实现。然而，纯粹的应力型断裂准则往往导致收敛性较差，因为裂纹尖端处预测应力随粒子间距减小而增大。此外，在没有核函数不一致处理的情况下，新形成的裂纹面会经历粒子分布畸变和噪声应力场。在SPH-PF [36,37]中，通过求解由离散断裂问题Lagrange泛函推导出的双曲PDE来获得时间演化相场参数。它用完全损伤区域与未损伤区域之间的平滑过渡来替代离散裂纹面。平滑裂纹近似的宽度由用户定义的尺度参数控制，这是实现SPH-PF框架中稳定断裂计算的关键因素。

基于裂纹均匀分布在单元（或粒子）上的假设，损伤力学（DM）通过由损伤相关变量确定的刚度降低来表示结构退化。通过引入基于黏聚力模型的损伤演化定律[38]，FE-DM方法能够实现稳健的失效模拟，减少与单元形状、尺寸和排列相关的裂纹路径偏置[[39], [40], [41], [42]]。在FE-DM中建模宏观断裂需要对自由度进行额外富集以表征裂纹张开，伴随复杂的编码算法。然而，在SPH中，可以通过简单地排除完全损伤粒子之间的相互作用来建模宏观断裂，自然地捕获从裂纹萌生、扩展、分支到碎裂的整个断裂过程。此外，如TLSPH-DM方法所示[43]，将梯度增强非局部性引入应变度量，有效缓解了裂纹面附近的核函数不一致问题，实现了更稳定的损伤演化和应力场。尽管其数学基础与相场不同，但在数值正则化方面起着类似的作用。凭借其相对简化的算法，TLSPH-DM达到了与复杂断裂模型如FE-PF [44]、FE-DM [[39], [40], [41], [42]]、LPM [45,46]和SPH-PF [36]相当的精度。

作为Kim和Kim [43]工作的后续研究，本文将TLSPH-DM扩展至涉及高度动态和非线性自由表面流动的水弹流固耦合问题。为此，TLSPH-DM与EISPH耦合，利用后者在计算流体粒子运动方面的高效性和稳定性。虽然EISPH和TLSPH-DM各自都是成熟的求解器，但与基于粒子键断裂、弹簧或格点连接来表示失效的FSI断裂方法（如SPH-PD、SPH-LPM、SPH-FPM、SPH-PS）相比，连续介质损伤SPH在双向不可压缩自由表面载荷下的耦合应用尚待探索。本研究的主要贡献如下：(i)一个物理一致的流固界面模型，包含双侧非穿透处理和损伤依赖的外部加速度缩放，配合具有时间自适应性的显式交错时间积分器，实现了准静态和动态流固耦合载荷下近脆性断裂中流体压力、结构响应和损伤演化的稳定预测；(ii)非局部梯度增强的等效应变公式有效抑制了非物理的沙漏模式，确保了水弹流固耦合条件下脆性断裂模拟的收敛性；(iii)将耦合EISPH-TLSPH-DM方法与六个基准问题进行全面验证，与上述断裂专用流固耦合求解器结果吻合良好。

本文其余部分组织如下：第2节介绍EISPH公式，包括控制方程、一阶投影法和拟拉格朗日方法。第3节描述TLSPH-DM框架，涵盖控制方程和线弹性断裂力学的本构模型。第4节介绍流固耦合策略，包括流固界面模型和时间积分器。第5节呈现耦合EISPH-TLSPH-DM模型的验证与性能评估。最后，第6节总结本研究的主要发现和贡献。附录A提供了有限应变线弹性-断裂公式的验证结果。所有模拟均使用自行开发的SPH求解器SOPHIA完成，该求解器使用CUDA C实现，并完全在图形处理器（GPU）上并行化[[47], [48], [49]]。

控制方程

对于流体求解器，假设为等温不可压缩牛顿流体。控制方程由质量和动量守恒定律组成，表示为

?·u = 0,
du/dt = -(1/ρ)?P + ν?2u + g,

其中u表示流速，ρ为流体密度，P为流体压力，ν为流体运动粘度，g = -g·ez为重力加速度，g = 9.81 m/s2。

一阶投影法

采用基于Helmholtz-Leray分解[52]的一阶投影法[50,51]来推进流体运动。

控制方程

全拉格朗日SPH（TLSPH）与连续介质损伤力学（DM）耦合，用于模拟线弹性断裂力学（以下简称TLSPH-DM）。质量和动量守恒定律表示为：

ρ? = det(F)ρ,
du/dt = (1/ρ?)??·P + g + a^COL + a^JST + a^FS + a^R,

其中上标0表示在初始（参考）构型中评估的量；F为变形梯度张量；P为第一Piola-Kirchhoff应力张量；??·(·)为在初始构型中定义的散度算子。

流固耦合策略

第2节给出用于流体求解器的显式不可压缩SPH，第3节给出用于结构求解器的含损伤力学的全拉格朗日SPH，分别呈现了统一全SPH流固耦合框架的流体和结构求解器。两个求解器之间的耦合通过流固界面处的粒子级动量交换来实现，详细内容见以下小节。

模型验证

为验证所提出的EISPH-TLSPH-DM方法，本研究考虑了六个基准问题。前三个案例涉及纯水弹响应，后三个案例涉及具有准脆性断裂的水弹行为。除第5.4节采用二维设置以降低GPU内存占用外，所有案例均执行三维模拟。

结论

本研究提出了耦合的EISPH-TLSPH-DM方法，用于三维涉及准脆性断裂的水弹流固耦合分析。EISPH基于投影法，其中流体压力通过PPE的显式求解获得（第2节）。通过在PPE中避免矩阵组装和迭代求解器，该方法实现了计算效率的提升。对于结构求解器，一种无张拉不稳定性TLSPH公式与连续介质损伤力学相结合……

CRediT作者贡献声明

金宇（Jin-Woo Kim）：原始草稿撰写、可视化、验证、方法论、调查研究、正式分析、数据整理、概念化。

禹熙相（Hee Sang Yoo）：方法论、软件、可视化、审稿与编辑。

金泳洙（Eung Soo Kim）：审稿与编辑、监督、软件、资源、项目管理、方法论、资金获取。

罗敏（Min Luo）：审稿与编辑、项目管理、资金获取。

金宇 | 禹熙相 | 金泳洙 | 罗敏