含颗粒流动与柔性结构相互作用的数值建模：一种耦合SPH-DEM-FEM方法

《Computers and Geotechnics》：Numerical modeling of particle?laden flows interacting with flexible structures: A coupled SPH-DEM-FEM approach

　　

摘要
颗粒载载多相流与柔性结构的相互作用在自然界和工程应用中普遍存在，例如泥石流防护网、柔性输送软管和土工布淤泥拦挡帘。然而，由于复杂的相间动量交换、严重的结构变形以及大量存在的不规则非球形颗粒的复杂形状，准确模拟这些动态过程仍然具有挑战性。为了解决这些挑战，本研究提出了一种完全解析的耦合框架，将光滑粒子流体动力学（SPH）、离散元方法（DEM）和有限元方法（FEM）集成在一起，以模拟颗粒多相流与柔性结构的相互作用。该框架通过完全解析的SPH-DEM方法处理流体-颗粒相互作用，其中非球形颗粒由DEM中的多面体模型描述。颗粒-结构相互作用通过双向DEM-FEM耦合处理，柔性结构由四面体FEM单元离散化。流体-结构相互作用使用双向SPH-FEM耦合方法求解，其中引入了基于重心坐标映射和虚拟粒子的动态边界更新技术，以准确反映柔性结构的变形。在SPH-DEM-FEM框架内，虚拟粒子描述非球形颗粒和固体结构的边界。通过实时更新其信息，跟踪移动边界并将相互作用力传递回FEM网格节点。为验证该框架，对两种溃坝冲击柔性板的工况进行了实验和模拟：一种为纯流体，另一种为含有非球形颗粒的颗粒载载多相系统。随后，对数值和实验结果进行了定性比较（例如，多相流的流动模式）和定量验证（例如，波前位移和柔性结构的动态位移）。比较结果表明，对于测试的溃坝冲击配置，所提出的SPH-DEM-FEM耦合模型在重现颗粒载载自由表面流与可变形柔性板相互作用方面提供了令人满意的准确性和可靠性。

1. 引言
颗粒载载多相流与柔性结构的相互作用在自然界和工程系统中广泛存在，例如泥石流防护网（Ishiki等，2025年；Vicari等，2022年）、柔性输送软管（Pa?doussis，2008年；Sun等，2024年）以及土工布淤泥拦挡帘或柔性网板（Cui和Bi，2026年；Radermacher等，2016年；Zhang等，2025年）。这些过程涉及流体、离散颗粒和可变形结构之间的强耦合作用。流体相将动量传递给颗粒和结构，而颗粒则经历平移、旋转、碰撞和摩擦接触。同时，柔性结构在流体和颗粒载荷下发生变形，并反过来改变局部流场和颗粒运动。非球形颗粒的存在进一步增加了复杂性，因为颗粒形状会影响接触行为和流体-颗粒相互作用。实验研究为这些问题提供了有价值的观测和基准数据（Antoci等，2007年；Song等，2018年），但同时也很难同时测量内部流场、颗粒尺度运动、接触力和瞬态结构载荷，特别是在涉及大结构变形的快速演化的多相系统中。因此，数值模拟已成为研究此类复杂耦合问题的关键替代和补充手段（Felippa等，2001年；Griffith和Patankar，2020年；Hou等，2012年）。

在颗粒载载流-结构相互作用问题中，通常采用不同的数值方法来模拟流体相、连续结构和离散颗粒。对于流体相，通常使用计算流体动力学（CFD）方法，可大致分为基于网格的方法，如有限体积法（FVM）（Patankar和Spalding，1983年），以及无网格粒子方法，以光滑粒子流体动力学（SPH）为代表（Gingold和Monaghan，1977年；Lucy，1977年）。对于连续固体结构，有多种数值方法可用，包括边界元方法、有限差分法和有限元方法（FEM）。其中，FEM因其对复杂几何形状的适应性以及纳入非线性材料行为的能力而被广泛采用。对于离散固体相，开发了基于粒子的方法，如无单元Galerkin法、近场动力学和离散元方法（DEM）。DEM显式跟踪单个颗粒的平移和旋转运动并解析其接触和碰撞相互作用，使其被广泛用于颗粒动态过程（Debnath等，2024年；Guo和Curtis，2015年；Ma等，2025c；Van Der Haven等，2024年；Zhu等，2008年）。传统DEM模型通常假设球形颗粒（Cleary，2010年；Zhu等，2007年；Zhu等，2008年）。然而，自然和工业系统中遇到的许多颗粒是非球形的，忽略颗粒形状会显著影响微观接触行为和宏观颗粒流特性，以及多相系统中的周围流体运动（Lu等，2015年；Ma等，2022年；Zhong等，2016年）。因此，开发了各种非球形粒子模型，包括多球模型（Abbaspour-Fard，2004年；Kruggel-Emden等，2008年；Liu等，2021b）、超椭球模型（Liu和Zhao，2020年；You和Zhao，2018年）和多面体模型（Eliá?，2014年；Xie等，2020年）。

在与结构相互作用的颗粒载载流中，完整的耦合问题涉及三个密切相关的相互作用：流体-颗粒相互作用、流体-结构相互作用和颗粒-结构相互作用。这些相互作用通常通过不同的二元耦合策略进行研究。对于流体-颗粒相互作用，通常使用CFD-DEM框架，其中流体相由基于网格的CFD方法求解，颗粒相由DEM描述。这些方法已被应用于沉积物输送和洪水沉积问题（Nie和Wang，2024年）、密集流化床、湿法研磨以及其他化学或矿物加工系统（Golshan等，2020年；Ma等，2025b；Zhao等，2024年），以及泥石流冲击和减灾问题（Bi等，2024年；Li等，2020年）。然而，当涉及大界面变形、移动边界或与可变形结构的强耦合时，基于网格的流体求解器可能需要额外的界面跟踪、网格运动或重新网格化处理。在这种背景下，SPH-DEM（He等，2018年；Robinson等，2014年；Sun等，2023年；Xu和Dong，2021年）耦合提供了一种有用的基于粒子的替代方案，因为流体和固体颗粒均以拉格朗日方式描述。这种共同的拉格朗日描述便于处理大自由表面变形和移动颗粒边界。

对于颗粒-结构相互作用问题，通常采用耦合的DEM-FEM方法来描述离散颗粒与连续结构之间的接触。这些方法通常可分为三类：用于脆性材料断裂的联合有限-离散元方法、耦合局部离散域和全域连续域的多尺度DEM-FEM方法，以及用于离散颗粒与连续结构相互作用的耦合DEM-FEM方法。本研究属于第三类，其中颗粒-结构接触力传递到FEM节点或单元，同时更新的结构配置反馈给DEM颗粒。早期的DEM-FEM研究主要关注球形颗粒与刚性或简单表示的结构边界之间的相互作用。最近，非球形粒子模型已被纳入DEM-FEM框架，以更好地捕捉颗粒形状效应和真实颗粒材料的接触行为（Liu等，2020b；Ma等，2025a；Qiao等，2022年；Sun等，2023年）。与此同时，颗粒与可变形或柔性结构（包括柔性屏障、可变形组件和柔性输送系统）之间的相互作用也引起了越来越多的关注（Liu等，2020a；Ma等，2025a；Michael等，2015年）。

流体-结构相互作用通常使用CFD-FEM耦合来建模，流体和结构求解器通过其界面交换运动学和力信息。对于涉及大变形柔性结构的问题，上述基于粒子的SPH方法也可以与FEM耦合，形成SPH-FEM框架，该框架无需流体网格运动或重新网格化即可适应不断演化的流体-结构界面。这种SPH-FEM方法已被应用于涉及移动或可变形结构边界的各种流体-结构相互作用问题（Caleyron等，2013年；Fuchs等，2021年；Li等，2023年；Yang等，2012年；Yao等，2023年）。最近的研究通过多分辨率处理、边界穿透防止技术和SPH粒子与结构表面网格之间的直接耦合进一步推进了三维SPH-FEM耦合（Guan等，2026年；Li等，2025年；Shi等，2026年）。然而，对于使用边界或虚拟粒子表示结构表面的耦合方案，在大结构变形下保持几何一致的移动边界并准确地将流体力传递到FEM节点仍然是一个重要的挑战。

如上所述，流体-颗粒-结构系统中涉及的三种二元相互作用，即流体-颗粒、流体-结构和颗粒-结构相互作用，已被广泛研究，相应的耦合方法相对成熟。这些发展为流体-颗粒-结构相互作用问题的集成模拟提供了重要的方法论基础。在这些二元耦合策略的基础上，几项研究已开始开发更全面的耦合框架，包括CFD-DEM-FEM、SPH-DEM-FEM和其他混合方法。然而，大多数现有应用主要关注刚性结构（Fang等，2022年；Jin等，2023年；Luo等，2022年；Ng等，2024年）。虽然几项最近的研究将这些框架扩展到柔性或可变形结构（Liu等，2021a；Luo等，2026年；Xiong等，2025年），但仍存在重要的方法论局限性。在这些研究中，通常在结构表示、颗粒形状或多相描述方面引入简化。例如，在结构建模方面，一些研究通过在DEM框架中将结构组件离散化为由弹簧-阻尼段连接的离散质量节点网络来简化柔性边界（Leonardi等，2016年；Li和Zhao，2018年；Xiong等，2025年）。然而，这种方法难以捕捉连续柔性结构的实际力学行为，也无法准确评估其内部应力分布。在多相流系统中，类似的简化也存在。例如，颗粒相通常由理想化的球形颗粒表示（Liu等，2021a），而忽略非球形形状的影响，或者将整个多相混合物直接简化为等效的单相流体（Du等，2026a）。此外，虽然已经建立了相对全面的CFD-DEM-FEM框架（Dai等，2021年；Luo等，2026年），但如上所述，使用传统的欧拉网格在涉及移动或高度可变形结构边界时可能需要额外的界面跟踪、网格运动或边界重建。这些附加处理可能会增加耦合算法的复杂性，并使维持鲁棒的流体-结构界面变得更加困难。

为了解决这些局限性，本研究开发了一种耦合的SPH-DEM-FEM框架，用于模拟颗粒载载流与柔性结构的相互作用。该框架提供了一种统一的处理方法，同时解析流体运动、非球形颗粒动力学以及柔性结构的连续变形和内部应力。在该框架内，SPH-DEM和DEM-FEM耦合方案基于作者之前的研究实现（Du等，2026b；Ma等，2025a；Shen等，2024年）。在此基础上，进一步开发了SPH-FEM界面处理以考虑不断变形的流体-结构边界。在此处理中，虚拟粒子通过重心坐标映射绑定到FEM网格的外部三角形面上，根据瞬时面几何进行动态更新，并用于将流体力传递到相应的FEM节点。所提出框架的可靠性通过与两个溃坝冲击案例进行实验评估来验证，一个涉及纯流体条件下的柔性板，另一个涉及颗粒载载条件下的柔性板。这些案例通过对流演变、前锋传播和柔性板变形的比较来评估模型。数值结果与实验观察结果吻合良好，证明了所提出框架的可靠性。

2. 数学方法

2.1. SPH方法

在本研究中，采用弱可压缩SPH（WCSPH）方法模拟流体相。在WCSPH框架内，流体相的运动由拉格朗日形式的Navier-Stokes方程控制。质量守恒方程和动量守恒方程表示如下：

（1）dρ/dt = -ρ?·u

（2）du/dt = -(1/ρ)?p + (1/ρ)?·τ + g

其中 ρ、u、p、τ 和 g 分别表示流体的密度、速度向量、压力、剪切应力张量和重力加速度。通过引入SPH方法的核近似和粒子离散化技术，上述连续控制方程被转换为以下离散粒子相互作用方程：

（3）dρ?/dt = ρ?∑?(m?/ρ?)(u?-u?)·??W?? + Γ

（4）du?/dt = -∑?m?(p?+p?)/(ρ?ρ?)??W?? + ∑?m?(τ?/ρ?2+τ?/ρ?2)·??W?? + g?

其中下标 i 和 j 分别表示参考粒子及其支撑域内的邻近粒子。m? 为粒子的质量。

Wij = W(xi-xj, h) 是光滑核函数，其中 h 表示光滑长度。Γ 是为减少密度和压力场中的高频数值振荡而引入的人工密度耗散项。其表达式及相关辅助项分别给出在公式(8)、(9)、(10)、(11)中。本研究将 h 设为初始 SPH 粒子间距 df 的 1.2 倍（即 h = 1.2df）。本研究采用经典三次样条函数作为光滑核函数，其数学表达式如下：

(5) W(q, h) = αD(2 ? 3q2 + q3), 0 ≤ q ≤ 1
(1/6)(2 ? q)3, 1 < q ≤ 2
0, q > 2

其中 q 为 r 与 h 的比值，r 为流体粒子间距离，三维问题中的归一化常数 αD 为 3/(2πh3)。

在 WCSPH 中，采用 Tait 状态方程闭合上述控制方程组。本研究采用了简洁且被广泛采用的状态方程形式（Antuono 等，2012），表达式如下：

(6) p = (c?2/2)(ρ ? ρ?)

其中 ρ? 为流体参考密度，c? 表示人工声速。根据弱可压缩假设（Monaghan，1994），为将模拟过程中流体密度波动限制在 1.0% 以内，人工声速 c? 通常由以下条件确定：

(7) c? ≥ max(10Umax, 10pmax/ρ?)

其中 Umax 和 pmax 分别表示流体的预期最大速度和最大压力。

离散连续性方程（公式(3)）中出现的人工密度耗散项 Γ，参照 Molteni 和 Colagrossi（2009）的定义为：

(8) Γ = 2δhc? ∑? (m?/ρ?) ψ?? (r?? · ??W??) / (|r??|2 + η2)

其中 δ 为无量纲耗散系数，建议取 0.1。r?? = r? ? r? 表示粒子 i 与粒子 j 之间的相对空间位置矢量。η = 0.01h 是为消除粒子恰好在同一位置时奇异性而引入的小参数，ψ?? 是与粒子密度及其梯度相关的项。为消除计算密度归一化梯度的需求，本研究采用了 Fourtakas 等人（2019）提出的 ψ?? 修正形式：

(9) ψ?? = ρ??? ? ρ???

其中上标 T 和 H 分别表示总分量和静水分量。根据公式(6)，静水密度差可由静水压力差确定。

(10) ρ??? = p??? / c?2

(11) p??? = ρ?g · z??

其中 z?? 表示粒子 i 与粒子 j 之间的垂直距离。ρ??? 和 p??? 分别表示粒子 i 与 j 之间的静水密度差和静水压力差。

本研究采用预测-校正算法（Gomez-Gesteira 等，2012；Monaghan，1989）作为时间推进算法。时间步长受扩散、声学、对流以及 Courant-Friedrichs-Lewy（CFL）条件的约束（Colagrossi 等，2016）：

(12) Δt_d = 0.125h2/νmax, Δt_a = CFL · h/amax, Δt_c = CFL · h/c?
Δt ≤ min(Δt_d, Δt_a, Δt_c)

其中 νmax 表示流体最大运动粘度，amax 表示流体和虚粒子的最大加速度，对于预测-校正格式，CFL 取 0.25（Cheng 等，2019）。

在 SPH 方法中，边界处理是不可缺少的组成部分，其主要发挥两个关键作用：防止流体粒子穿透固体壁面，以及通过施加物理边界条件来提高近壁流场的模拟精度。本研究采用 Adami 等人（2012）提出的固定虚粒子方法来处理固体边界，并结合 Du 等人（2026b）开发的改进虚粒子生成方法，生成贴合物体且均匀分布的壁面虚粒子。在模拟过程中，虚粒子相对于局部边界的相对构型保持不变，其全局坐标根据固体边界的运动或变形动态更新（关于柔性结构上虚粒子的具体位置更新方法将在第 3.1 节中详述）。

本研究对所有固体壁面（包括 DEM 中的固体粒子）施加无滑移边界条件。基于力的平衡原理和无滑移条件，通过相邻流体粒子的插值和外推获得虚粒子的压力和速度，密度则通过状态方程求解。随后，携带这些更新物理属性的虚粒子参与边界附近流体粒子守恒方程（即公式(3)、(4)）的计算。此外，为有效防止复杂流固相互作用场景中的流体穿透，本研究引入了一种改进的边界排斥力模型（Du 等，2026b），当流体粒子过于靠近虚粒子时，对其施加额外的法向排斥力。关于虚粒子物理属性计算和排斥力的详细控制方程，请参阅上述文献（Du 等，2025；Du 等，2026b）。

2.2. DEM 模型

本研究采用 DEM 模拟离散固体粒子的动态行为。在 DEM 中，每个粒子的平动和转动由牛顿第二定律控制，表达式如下：

(13) m·dv/dt = ∑F_c + F_f + mg

(14) I·dω/dt = ∑T_c + T_f

其中 m、I、v 和 ω 分别表示固体粒子的质量、转动惯量、平动速度和角速度。F_c 和 T_c 分别为粒子之间或粒子与壁面之间的接触力和接触力矩，F_f 和 T_f 分别表示施加在粒子上的颗粒-流体相互作用力和力矩。需要注意的是，本文引入了多面体模型来近似非球形粒子。本研究所用多面体粒子的接触检测算法和接触力模型已在我们的先前工作中详细描述（Liu 等，2021b；Xie 等，2020）。

2.3. FEM 模型

FEM 的核心概念是将连续域离散化为有限数量的互连单元，通过求解节点位移来确定整体结构响应。在结构动力学分析中，有限元系统的离散控制方程通常表示为：

(15) M·? + C·? + K·x = F_ext

其中 M、C 和 K 分别表示结构的全局质量矩阵、阻尼矩阵和刚度矩阵。?、? 和 x 分别表示节点加速度、速度和位移向量。F_ext 表示施加在结构上的外部载荷向量，包括本研究中流体相和固体粒子施加的相互作用力。本研究采用知名的各向同性线性弹性模型和胡克弹性定律来确定 FEM 框架中的上述位移、变形和应力。关于 FEM 实施的进一步细节已在文献中充分确立（Michael 等，2015），此处从略。

3. 耦合策略

本节描述集成 SPH-DEM-FEM 框架中使用的耦合策略。首先介绍 SPH-FEM 耦合方法，因为它为流固相互作用提供了动态边界处理和力传递策略。在该模块中，边界虚粒子通过重心坐标映射和局部正交坐标系进行更新。随后概述全解析 SPH-DEM 模块以描述流体-颗粒相互作用，接着介绍涉及非球形粒子和可变形结构的颗粒-结构相互作用的 DEM-FEM 耦合模块。最后，展示集成框架的整体求解流程和数据处理过程。

3.1. SPH-FEM 耦合算法

为解决流体与柔性结构之间的相互作用，本研究开发了一种 SPH-FEM 耦合框架。如图 1 所示，固体结构被划分为 FEM 所需的四面体网格以模拟其结构变形。同时，利用 Du 等人（2026b）提出的虚粒子生成方法，在柔性体表面布置贴合物体且均匀分布的边界虚粒子以处理流固相互作用。基于这些预处理工作，提出了一种使用局部正交坐标系和重心坐标映射的全解析耦合方案。该方案将边界虚粒子绑定到柔性体外表面三角面上，根据网格的瞬态几何信息动态更新其位置和法向量，同时促进流体与结构之间的相间力传递。值得注意的是，在模拟中，这些离散 SPH 虚粒子在结构壁内生成并布置为两层，距结构边界分别为 0.5df 和 1.5df，但为简洁和清晰起见，图 1 中仅绘制了外壁上的单层。该处理的主要步骤总结在图 2 中。

图 1. 柔性结构的网格划分与边界离散化示意图。
图 2. SPH-FEM 耦合的整体流程：(a) 虚粒子生成和映射到 FEM 表面的预处理阶段，(b) 虚粒子更新和流体力传递的时间步进阶段。

3.1.1. 边界虚粒子的几何映射和局部坐标绑定

在耦合模拟开始前，每个边界虚粒子与 FEM 表面的一个外三角面相关联。此处，外三角面是指从四面体 FEM 网格中提取的表面面。虚粒子的初始法向方向，记为 n_b，由描述初始结构边界的水平集函数的梯度获得。该法向方向用于将虚粒子投影到 FEM 表面上并确定其关联面。

对于每个虚粒子，首先在其邻域内搜索候选外三角面。然后将虚粒子沿其初始边界法向方向 n_b 投影到这些候选面上。当投影点位于三角面内时，认为该投影有效，这可以通过投影点的重心坐标进行检查。如果存在多个有效候选，则选择投影距离最小的面作为关联面。一旦建立这种关联，面的 ID 和相应的映射信息将被存储并在本次模拟中保持固定。只要 FEM 表面保持其连通性且未发生面翻转、自接触、断裂或拓扑变化，这种固定关联就仍然有效。对于涉及此类事件的结构，则需要动态重新绑定。

对于选定的具有顶点 v? = (x?, y?, z?)、v? = (x?, y?, z?) 和 v? = (x?, y?, z?) 的三角面，投影点 p' = (x_p, y_p, z_p) 使用重心坐标表示，如图 3 所示：

(16) p' = λ?v? + λ?v? + λ?v?

系数 λ?、λ? 和 λ? 由面积比获得：

(17) A = ?|(v? ? v?) × (v? ? v?)|
λ? = |(v? ? p') × (v? ? p')| / (2A)
λ? = |(v? ? p') × (v? ? p')| / (2A)
λ? = 1 ? λ? ? λ?

其中 A 为三角面的总面积。重心坐标被存储用于虚粒子更新和力传递。

对于位于平滑表面区域的虚粒子，其初始虚粒子法向 n_b 通常与关联面的法向量平行。在这种情况下，只需存储虚粒子与投影点之间的法向偏移距离，如图 3(a) 所示。然而，对于位于边缘或角落附近的虚粒子，n_b 可能不与关联面的法向量平行，如图 3(b) 所示。在这种情况下，沿面法向直接更新虚粒子可能会改变其相对于物理边界的初始相对位置。因此，引入局部坐标系以保留完整的相对偏移矢量。具体而言，计算虚粒子初始全局坐标 r_b0 与投影点 p'0 之间的全局相对偏移矢量 d_global = r_b0 ? p'0。然后将该矢量转换到关联面的局部坐标系中：

(18) d_local = E? · d_global

其中 E = [e_x, e_y, e_z] 是在初始关联面上构建的局部正交基矩阵。这里，e_z 与关联面的法向方向对齐，而 e_x 和 e_y 位于面平面内，通过对面的边矢量进行正交化确定。通过存储 d_local，即使初始水平集法向与关联面法向不平行，虚粒子也能维持其与变形 FEM 表面的初始相对几何关系。

3.1.2. 虚粒子坐标更新及力向 FEM 节点的传递

在耦合模拟过程中，柔性结构的变形由外部 FEM 求解器求解。

在每个耦合步骤中，有限元（FEM）表面单元的当前节点坐标和速度通过数据交换接口直接获取。对于与边界虚拟粒子关联的三角形单元，其三个节点的当前坐标分别记为v1(t)、v2(t)和v3(t)。利用初始绑定阶段存储的重心坐标，当前投影点可重建为：

(19) p′(t) = λ1v1(t) + λ2v2(t) + λ3v3(t)

如果虚拟粒子由法向偏移量描述，其当前位置更新为：

(20) rb(t) = p′(t) - db·n(t)

其中db为有符号法向偏移量，n(t)为瞬时单元法向量。对于本研究采用的双层虚拟粒子排布方式，第一层的|db| = 0.5df，第二层的|db| = 1.5df，符号由结构壁面的内法向方向决定。对于由局部偏移向量描述的虚拟粒子，位置通过下式更新：

(21) rb(t) = p′(t) + E(t)dlocal

其中E(t)为由当前单元几何信息重建的局部坐标变换矩阵。

对于采用局部坐标格式更新的虚拟粒子，虚拟粒子的法向量并不直接赋值为所关联FEM单元的法向量。相反，在投影点p'(t)和虚拟粒子位置rb(t)更新之后，虚拟粒子的法向量通过下式重建：

(22) nb(t) = (rb(t) - p′(t)) / |rb(t) - p′(t)|

这种处理方式使得位于边缘和角落附近的虚拟粒子在结构变形过程中能够保持其相对于FEM表面的初始相对偏移方向。

为满足SPH框架内的无滑移边界条件，虚拟粒子的瞬时物理速度ub(t)通过分析直接由节点的运动学信息获得。考虑到虚拟粒子沿单元法向的局部偏移距离极小，可以合理地忽略单元瞬时旋转引起的线速度增量。因此，虚拟粒子的瞬时速度可近似为其投影点p′(t)的平动速度，该速度由单元的三顶点瞬时平动速度插值得到：

(23) ub(t) = Σ(i=1 to 3) λi·vnode,i(t)

其中vnode,i(t)表示节点i的瞬时平动速度。

在SPH求解器计算完作用在边界虚拟粒子上的力之后，这些力被传递回FEM结构。对于每个关联的三角形单元，作用于第k个虚拟粒子的力通过其重心权重分配到三个FEM节点上：

(24) Fnode,i = Σ(k=1 to Nb) λi,k·Fb,k, (i = 1, 2, 3)

其中Fnode,i为节点i接收的总等效力，Nb为绑定到该三角形单元上的虚拟粒子总数，λi,k为第k个虚拟粒子相对于节点i的重心坐标权重，Fb,k为作用在第k个虚拟粒子上的流体相互作用力与边界排斥力之和。通过上述双向数据交换过程，实现了SPH与FEM之间的耦合。

3.2. SPH-DEM耦合算法

本研究采用Du等人（2026b）开发的完全解析SPH-DEM耦合框架来处理流体-颗粒相互作用。在该框架中，固体颗粒被视为SPH流体域内的移动固体边界。不再使用经验或半经验阻力关联式，而是直接从流体粒子与布置在颗粒表面的边界虚拟粒子之间的局部SPH相互作用中评估作用在每个DEM颗粒上的水动力。因此，颗粒边界及其周围流体运动在SPH离散化中被显式表示。这种完全解析的处理方式要求流体粒子间距df充分小于特征颗粒尺寸D，其中D定义为DEM颗粒的等效直径。典型的分辨率准则为df < D/10。因此，本方法更适合包含适中数量、相对较大且粒径分布相对较窄颗粒的系统。对于涉及大量颗粒或宽粒径分布的系统，所需的SPH分辨率可能导致相当大的计算成本。

基于上述方法，前述虚拟粒子边界方法也应用于固体颗粒。作用在单个DEM颗粒上的总水动力Ff和总力矩Tf按下式计算：

(25) Ff = Σ(b) Σ(f) Ffb

(26) Tf = Σ(b) Σ(f) (rf + rb/2 - ro) × Ffb

其中Ffb表示流体粒子对虚拟粒子施加的力（下标b代表固体边界），包括压力梯度和粘性力，以及作用于边界附近流体粒子以防止非物理穿透的排斥机制反作用力。ro代表固体颗粒的质心位置。

3.3. DEM-FEM耦合算法

为描述离散固体颗粒与柔性结构之间的相互作用，本研究采用了双向耦合的DEM-FEM框架。柔性结构使用四面体FEM单元进行离散化，从四面体网格中提取的外部三角形表面单元在DEM求解器中被视为可变形固体边界。这样，颗粒-结构相互作用即表示为单个DEM颗粒与可变形结构的三角形表面单元之间的接触相互作用。

非球形颗粒使用多面体DEM模型来描述。多面体颗粒所采用的接触模型基于Feng等人（2012）提出的弹性接触能量守恒原理。接触检测和接触力计算算法的详细描述可参见我们之前的研究（Liu et al., 2021b, Xie et al., 2020）。双向DEM-FEM耦合通过力传递和边界几何更新来实现。在DEM到FEM一侧，固体颗粒与三角形表面单元之间的接触力首先在DEM求解器中计算。每个接触力然后根据有限元形函数在精确接触点处的评估值，分配到被接触三角形单元的三顶点上，这些顶点同时也是四面体FEM网格对应的表面节点。这些映射的节点力作为外部动力载荷施加在FEM求解器中，驱动节点位移、结构变形和柔性结构的瞬态运动学响应。

在FEM到DEM一侧，当柔性结构在颗粒冲击下发生变形后，更新后的表面网格节点坐标和速度被传递回DEM求解器。DEM求解器然后根据更新的节点信息重建可变形三角形表面边界，以便后续的颗粒-结构接触计算使用瞬时结构构型。该双向DEM-FEM耦合算法的更多细节可参见我们之前的研究（Ma et al., 2025a, Shen et al., 2024）。

3.4. SPH-DEM-FEM耦合框架

集成SPH-DEM-FEM框架的总体求解流程和数据传递在图4中示出。本框架采用显式分区交错式双向耦合策略。SPH、DEM和FEM求解器使用相同的时间步独立推进。几何信息、运动学信息和力信息在每个时间步通过共享内存进行交换，且在一个耦合周期内不执行子迭代。

图4. 耦合SPH-DEM-FEM框架示意图。

在该框架中，使用自研DEM软件DEMSLab求解固体颗粒的运动和接触相互作用，使用内部SPH求解器求解流体流动，柔性结构的FEM变形则使用RecurDyn求解。在每个耦合周期开始时，进行初始同步。SPH求解器从共享内存中读取最新的DEM颗粒信息，包括颗粒位置、朝向和速度，以及当前FEM表面节点坐标和速度。这些数据用于更新SPH域中DEM颗粒的移动边界表示和柔性结构的变形。同时，DEM求解器读取更新后的FEM表面网格信息以及SPH求解器在上一周期计算的水动力和力矩。FEM求解器接收来自SPH求解器的流体力以及来自DEM求解器的颗粒-结构接触力，这些力同样是在上一周期中获得的。

同步之后，三个求解器独立推进其计算。SPH求解器更新边界虚拟粒子，推进流体粒子，评估作用在DEM颗粒上的水动力和力矩，以及作用在结构边界虚拟粒子上的流体力。DEM求解器计算颗粒间接触力和颗粒-结构接触力，并将这些接触力与上一步同步中接收的水动力相结合，更新DEM颗粒的平动和转动运动。与此同时，FEM求解器将接收到的流体载荷和接触载荷作为外部动力载荷，计算柔性结构的瞬态变形和节点运动。在耦合周期结束时，进行二次同步以准备下一周期所需的数据。SPH求解器将新计算的水动力、力矩和结构流体载荷写入共享内存。同时，DEM求解器写入更新后的颗粒位置、朝向和速度，FEM求解器写入更新后的表面节点坐标和速度。这些数据将在下一个耦合周期中使用。

由于数据交换遵循显式交错顺序，交换的力和边界运动信息在耦合周期内被顺序使用。这引入了显式分区方案通常带来的一步数据延迟。在本研究的模拟中，通过使用足够小的耦合时间步和适当的同步间隔来控制这一效应。第4节中给出的验证结果表明，所采用的耦合策略为本研究所考虑的水坝溃决冲击问题提供了稳定且准确的结果。

4. 验证

为验证所开发的SPH-DEM-FEM框架，进行了水坝溃决试验以验证其在预测涉及柔性结构和非球形颗粒悬浮多相流的复杂行为方面的准确性和稳定性，并将模拟结果与实验数据进行对比。鉴于本框架中使用的完全解析SPH-DEM和DEM-FEM耦合子模型已在先前的工作中得到验证（Du et al., 2026b, Ma et al., 2025a），本节的验证工作将专门针对用于求解流固相互作用的SPH-FEM耦合模型和SPH-DEM-FEM耦合框架。

为系统地执行验证流程，针对两个涉及水坝溃决流冲击柔性板的标准场景进行了实验研究和数值模拟。这两个场景分别涉及纯流体条件下的水坝溃决流冲击柔性板以及颗粒悬浮多相条件下的水坝溃决流冲击柔性板。具体而言，纯流体水坝溃决冲击柔性板场景旨在验证双向SPH-FEM耦合框架的基准准确性（第4.3节）。随后，颗粒悬浮水坝溃决冲击柔性板场景被用来在复杂多相冲击条件下全面检验完全耦合SPH-DEM-FEM框架的可靠性和准确性（第4.4节）。在进行这些验证分析之前，两个案例的实验装置和相应的模拟条件分别在第4.1节和第4.2节中详细说明。

4.1. 物理实验

水坝溃决流冲击柔性板的实验在一个由高度透明亚克力玻璃构建的矩形水箱中进行。水箱的尺寸和实验装置的示意图如图5所示。水箱内部安装有一个可沿导轨槽滑动的可拆卸闸门，用于控制流体或流体-固体颗粒混合物的释放。为确保上游蓄水池的严格密封性并在闸门移除期间最小化滑动摩擦，在实验前将涂有凡士林的泡沫条附着在闸门的两个垂直边缘上。闸门的驱动装置设计类似于Lobovsky等人（2014）的设计，即通过定滑轮利用自由下落的重物为闸门提供驱动力。在本研究的所有实验中，闸门的平均提取速度约为2 m/s。同时，一块4 mm厚的聚氨酯（PU）柔性板垂直悬挂在下游区域。该板由两块间距4.0 mm的亚克力板在150 mm高度处紧密夹持固定，从而建立了固定支撑边界条件。此外，在整个实验过程中，水箱底部始终保持水平位置。

图5. 水坝溃决实验中带有柔性板的水箱尺寸和布置。

在每次实验运行前，上游蓄水池的液位始终保持在200 mm，而下游通道保持完全干燥。

针对载粒溃坝冲击柔性板的工况，进一步引入了非球形颗粒，在储水底部分布了总计100颗多面体颗粒，堆积高度约为50 mm（如图6所示）。多面体颗粒由树脂通过3D打印工艺制造，制造公差为0.1 mm。颗粒密度为1142 kg/m3，等效直径为16.95 mm。在实验过程中，待初始状态稳定后，释放自由落体配重以快速抽出闸门，从而触发溃坝流动。同时，为确保数据采集的时间同步性，在y轴正交方向布置一台工作帧率为500 fps的高速相机，与闸门运动同步触发，记录瞬时流动行为和柔性板的动态变形。此外，为减少随机误差并确保所观测现象的可重复性，每种实验工况（包括纯流体和载粒工况）均独立重复三次。

4.2 模型设置

对于流体相，初始流体粒子间距设为1.0 mm，光滑长度定义为1.2 df。所有固体壁面均施加无滑移边界条件。在纯流体溃坝工况中，流体粒子和边界虚粒子的总数分别约为200万和97万。类似地，在载粒多相溃坝工况中，这两个数值分别约为168万和115万。在溃坝释放前，允许流体粒子在重力作用下沉降，以耗散动能并达到稳定的初始状态。

如前所述，柔性板的材料选用聚氨酯（PU）。本质上，PU是一种粘弹性材料，表现出显著的应变率效应（Yi等，2006），其动态杨氏模量在更高载荷速率下会增大。尽管如此，已有研究表明，聚氨酯材料的主要刚度响应可在明确定义的变形范围内用等效线弹性参数表示。例如，Kanyanta等（2009）利用实验测定的弹性属性对PU在特定载荷条件下的行为进行建模。近期，Raf等（2026）报道，SLS制造的TPU的平衡应力-应变响应在应变约10%以内近似为线性，并相应地在该范围内使用线弹性模型确定有效弹性参数。在本实验中，柔性板的最大应变保持在5%以下，且在此范围内的应力-应变响应在各测试应变率下均近似为线性。因此，采用等效线弹性规律来表征当前短时冲击条件下柔性板在加载阶段的刚度。为考虑应变率依赖的刚化效应，杨氏模量根据实测的应变率响应和特征冲击应变率确定，而非直接规定为一般静态值。尽管等效线弹性处理对当前加载范围和主要响应量而言是一种合理的工程简化，但它不可避免地忽略了粘弹性阻尼和应变率相关的滞后效应。然而，这些效应对本研究所考虑的短时冲击响应的影响预计是有限的。

受Eyring应变率过程理论（Bauwens-Crowet等，1969，该理论描述了聚合物变形和屈服的应变率对数依赖性）的启发，此处采用半经验对数关系来拟合实测的应变率依赖杨氏模量。PU的宏观变形涉及克服与分子链运动相关的能垒，其力学响应通常表现出对应变率对数的依赖性（即材料在外动态载荷下变形的速率）（Mulliken和Boyce，2006）。据此，采用以下半经验关系确定等效动态杨氏模量：

E(ε?)=E?+K ln(ε?/ε??) ……(27)

其中E?为参考应变率ε??=10?3 s?1下的准静态杨氏模量，K为拟合材料常数，ε?为当前特征应变率。为确定该方程中的未知常数（E?和K），使用万能试验机（Somarathna等，2020）对PU材料进行标准单轴拉伸试验。具体而言，按照ASTM D412（Die C）标准制备哑铃形试样，在室温下以三种不同加载速率进行测试。图7给出了应变率分别为10?3、10?2和10?1 s?1时的应力-应变曲线。应注意，本实验溃坝条件下观测到的最大结构应变小于5%。因此，通过专门在5%以下的应变范围内进行线性拟合，计算出三种应变率下的杨氏模量分别为37.2 MPa、41.5 MPa和45.6 MPa。将这三组数据代入式(27)并进行线性拟合，确定常数E?和K的值分别为37.2 MPa和1.82 MPa。最后，为评估溃坝冲击期间的特征应变率，分析了无柔性板存在时的流动运动学。在柔性板的具体位置（x=0.154 m），提取初始冲击速度（即前锋波的峰值推进速度）为v=1.5 m/s。已知柔性板的宏观结构高度为H=0.15 m，估计冲击阶段具有代表性的宏观应变率为ε?=v/H=10 s?1。将此工作应变率代回标定后的经验公式(27)，得到等效动态杨氏模量约为54 MPa。该特定值随后用作后续仿真中弹性本构模型的输入参数。

与实验约束条件一致，柔性板的上端指定为固定支撑。PU材料密度设为实验测定的值1161.75 kg/m3。为反映PU材料的近不可压缩性，泊松比指定为0.49。此外，采用从应变率依赖拉伸试验和估计的特征冲击应变率中标定得到的等效动态杨氏模量54 MPa作为弹性本构模型。该标定的模量用于表征测试的短时冲击条件下PU板的主要刚度响应。

对于载粒多相溃坝工况，颗粒-颗粒、颗粒-水槽和颗粒-板之间的摩擦系数分别指定为0.35、0.2和0.25。类似地，这些相互作用对应的恢复系数分别定义为0.7、0.8和0.6。此外，对于刚性接触（颗粒-颗粒和颗粒-水槽界面），法向和切向弹簧常数分别设为7.5×10? N/m和5.5×10? N/m。为适当考虑PU材料的弹性，降低了颗粒-板相互作用的接触刚度，法向和切向弹簧常数分别指定为1.5×10? N/m和1.1×10? N/m。初始颗粒-流体构型通过初步的SPH-DEM沉降仿真获得。多面体颗粒与流体一起引入左侧储水池，在此初始化阶段保持闸门关闭。随后推进耦合的SPH-DEM计算，直至自由液面和颗粒堆积体均趋于稳定。所得到的颗粒堆积高度约为50 mm，相应的流体状态作为后续溃坝冲击仿真的初始条件。

所有仿真的积分时间步长统一设为8×10?? s。此外，在测试开始前，柔性板自由下端的最下方点（图5中以实心红色圆圈标出）与水槽左壁之间的水平距离为0.154 m。在后续分析中，该点的位移被选为两种工况的核心定量评价指标。除非另有说明，其余仿真参数与相应实验条件保持一致。

4.3 纯流体工况下的模型验证

本节主要关注下游推进前锋的演化以及柔性板的动态变形行为。图8首先给出了实验和仿真获得的流动行为快照。如图8中的流动快照所示，闸门抽出后，上游储水池中的流体在重力作用下向下游涌动。在流体对柔性板初始冲击时（t≈72 ms），动能不足以引起任何可视觉辨认的结构变形。因此，在板前发生轻微的流体堆积。在持续的水动力冲击驱动下，柔性板逐渐发生弯曲变形，其自由下端尖部逐步抬升。这种变形使流体能够迅速从板下方通过。与此同时，向下游推进的流体前锋在垂直于流动方向上呈现出相对光滑的剖面。在224 ms时，该高速流体前锋冲击远端下游垂直壁面，产生强烈的垂直涌浪。同时，另一部分流体逐渐沿柔性板上爬（160 ms）。当其接近上端夹持处时，上爬的流体受到重力影响，形成明显的反向卷曲坍塌破碎波，随后回落至底部流体层（224 ms和280 ms）。

为定量支持视觉观察到的流动形态，提取了前锋位移（误差线代表三次独立重复实验的标准偏差）以及柔性板下端尖部的水平和垂直位移的时间演化，分别如图9和图10所示。在图9中，实验数据以均值和标准偏差误差线表示，连续线代表数值结果。在初始阶段（t<72 ms），前锋在干渠中表现出稳定的非线性加速推进，如图9早期阶段所示。在冲击柔性板时（t≈72 ms），初始流体冲击不足以引起柔性板的宏观可见变形，导致结构前方出现短暂的流体堆积。这种物理流阻直接导致流体向下游传播出现瞬时减速，在图9前锋位移曲线上表现为轻微的平缓段。

在流体动量传递的直接驱动下，柔性板在t=72 ms后立即开始变形。随着流体不断堆积并持续向下游推挤，柔性结构发生快速弯曲。该加载过程反映在图10中上升的结构位移曲线上，柔性板最终达到最大挠度（水平方向约85 mm，垂直方向约31 mm）。随后，随着堆积流体逐渐从板下方通过，水动力载荷减小，柔性板进入卸载和回弹阶段（t>224 ms）。

为进一步检验SPH-FEM耦合处理的压力响应，图11在选定时刻展示了流体压力分布和柔性板的冯·米塞斯应力分布。在t=160 ms时，流体已冲击柔性板，并在结构前方形成高压区域。随着柔性板继续变形，压力场随局部流体运动和变化的结构构型而演化。在t=224 ms、280 ms和320 ms时，流体在变形板附近区域堆积并通过。

放大视图进一步展示了冲击和变形过程中流体-结构界面附近的压力演化。下载：下载高分辨率图像（341KB）下载：下载全尺寸图像图11.纯流体溃坝工况下流体的压力分布和柔性板的冯·米塞斯应力分布。

总之，所提出的SPH-FEM耦合框架在流体-结构相互作用整个过程中与实验测量结果表现出良好的一致性，如图9和图10中波前传播和柔性板位移的对比所示。定量地讲，在所有考虑的工况中，预测的波前时间演化和结构位移的平均误差均保持在3%以下。此外，图11中所示的压力分布在冲击过程中保持物理合理。在所选压力快照中未观察到明显的非物理压力集中现象，且压力场在所示阶段中平稳演化。这些结果表明，所提出的SPH-FEM耦合方法能够准确捕捉流体运动和结构变形，同时保持稳定的压力响应，展现了其在涉及柔性结构的流体-结构相互作用问题中的精度和数值稳定性。

4.4. 含颗粒工况下的模型验证

为了对多相流体-结构相互作用过程进行直观和定性的评估，图12展示了含颗粒溃坝流冲击柔性板的实验和模拟的连续快照。在流体-固体混合物释放后，流体相最初超过颗粒相向下游推进。与纯流体工况相比，含颗粒工况中前导流体剖面明显变薄，这是由于颗粒相的存在和影响所致。当这层衰减的前导流体最初接触柔性板时（t ≈ 76 ms），其动能不足以瞬间触发宏观可辨的结构变形。因此，在板前方形成局部流体积聚，且持续时间比纯流体工况更长。随后，当流体-固体混合物的主体冲击板时，水动力载荷逐渐增加并驱动柔性板弯曲。在变形早期阶段，下部自由尖端的抬升高度不足以让颗粒通过，导致固体颗粒在板前方短暂积聚（t = 160 ms）。相比之下，流体相能够更早地从部分抬升的板下方通过。随着柔性板底部持续抬升，固体颗粒可以从其下方通过，并在周围流体的水动力作用下加速。在t = 240 ms时，前导流体冲击水箱最右侧的垂直壁，形成一个明显的垂直射流，固体颗粒紧随其后在产生的冲击力作用下向上运动。在此阶段，由于颗粒之间、流体与刚性壁之间的剧烈碰撞，自由液面经历严重的变形、破碎和飞溅。与此同时，在板的夹持上端附近，一部分流体向上攀爬并最终形成重力驱动的坠落式碎波（t = 240 ms和280 ms），表现出与纯流体工况中观察到的类似动态特征。

下载：下载高分辨率图像（466KB）下载：下载全尺寸图像图12.含颗粒溃坝流冲击柔性板的实验和模拟连续快照。

为了定量追踪上述过程，图13展示了流体前锋位移和固体颗粒前锋位移的时间演化。图13(a)中流体前锋位移在最初与柔性板接触后呈现出明显的平台期。该平台期对应于图12中观察到的瞬态流体积聚现象。由于前导流体层比纯流体工况中的更薄，在板变形足以允许进一步传播之前，流体前锋被阻挡的时间更长。图14中纯流体工况与含颗粒工况的波前位移历史的对比进一步支持了这一解释。具体而言，该图揭示含颗粒工况中的波前平台期明显长于纯流体工况。此外，根据图9和图13，含颗粒工况中的实验误差棒比纯流体工况更宽，这可归因于初始颗粒堆积状态的不确定性以及颗粒运动和接触导致的随机波动增强。

下载：下载高分辨率图像（198KB）下载：下载全尺寸图像图13.流体和固体颗粒前锋位移的时间演化：(a) 流体前锋，(b) 固体颗粒前锋。下载：下载高分分辨率图像（121KB）下载：下载全尺寸图像图14.纯流体工况与含颗粒工况的波前位移历史对比。

图13(b)中固体颗粒前锋的变化展现了根本不同的物理响应。与流体相不同，颗粒前锋位移曲线没有明显的平台期或明显的减速现象。为了进一步探究这一现象的根本原因，图15中提取了柔性板下部自由尖端的位移时间历史。结合图13(b)和图15的观察，当颗粒前锋到达柔性板时（即前锋位移为54 mm、时间约为100 ms），下部自由尖端已经积累了相当程度的垂直位移。这是因为前导流体比颗粒更早到达，从而启动了柔性结构的加载过程。虽然这一初始位移不足以让颗粒通过，但这一预变形提前开始了开启过程，从而有效缩短了颗粒到达后柔性板达到足够间隙高度所需的时间。随后，在颗粒相和流体相的持续冲击下，柔性板底部尖端迅速抬升并快速超过颗粒尺寸（在约180 ms时达到约20 mm的抬升高度）。因此，颗粒在板前方仅经历非常短暂的积聚，随后从板下方通过并在流体驱动力作用下加速。此外，颗粒前锋的位移曲线表现出持续前进的趋势。

下载：下载高分辨率图像（176KB）下载：下载全尺寸图像图15.柔性板下部自由尖端的实验与模拟位移时间历史对比：(a) 水平位移，(b) 垂直位移。

此后，当颗粒顺利通过柔性板下方时，颗粒与板之间几乎没有直接接触。这表明柔性板的后续变形主要受流体流动影响。此外，本研究中的颗粒速度从整体来看相对较低。特别是当柔性板达到最大开度时，底部区域几乎完全被颗粒覆盖，此时颗粒速度应趋近于零。这意味着颗粒对流体流动行为的影响相对较弱。因此，纯流体工况与含颗粒工况中柔性板的最大位移应该相似。图16进一步展示了纯流体工况与含颗粒工况中板底部自由尖端的垂直位移对比。图16显示两种工况中下部自由尖端的最大位移没有显著差异，这进一步验证了上述推断。最后，随着柔性板左侧的水动力载荷逐渐减小，结构进入卸载和回弹阶段。

下载：下载高分辨率图像（95KB）下载：下载全尺寸图像图16.纯流体工况和含颗粒工况的柔性板位移完整时间历史。

为了进一步检验含颗粒溃坝工况下的耦合响应，图17展示了选定时刻流体的压力分布和柔性板的冯·米塞斯应力分布。在t = 160 ms时，含颗粒流体已达到柔性板，压力场开始在冲击区域附近集中。随着板的变形，局部压力分布随流动模式和结构构型的改变而演化。在t = 240 ms、280 ms和320 ms时，流体和颗粒继续与变形的板相互作用，而柔性结构的应力分布与其弯曲变形保持一致。放大视图进一步展示了耦合交互区域附近的局部压力和应力演化。

下载：下载高分辨率图像（343KB）下载：下载全尺寸图像图17.含颗粒溃坝冲击工况下流体的压力分布和柔性板的冯·米塞斯应力分布。

总体而言，含颗粒溃坝流冲击柔性板物理过程的实验与数值结果取得了良好的一致性，定性快照对比和定量分析均表明模拟与实验之间仅有小差异。压力分布和应力分布进一步表明，在所考虑的短持续时间冲击条件下，所提出的SPH-DEM-FEM框架能够提供稳定且物理合理的流体-颗粒-结构耦合响应描述。具体而言，所有工况中流体前锋、颗粒前锋以及结构位移的预测误差均保持在5.8%以下。这表明SPH-DEM-FEM模拟能够成功捕捉受柔性结构变形影响的含颗粒复杂多相流动行为。这些结果支持了所提出框架在本研究考虑的含颗粒溃坝冲击构型中的可靠性。

5. 结论

本研究开发了一个全耦合的SPH-DEM-FEM数值框架，用于模拟流体、非球形颗粒和柔性结构之间的复杂相互作用。在该框架中，SPH用于捕捉剧烈的自由表面流动，DEM用于求解非球形颗粒的运动学和接触力学，FEM用于计算结构变形。同时，引入了一种基于重心坐标系和局部坐标系的虚拟粒子更新算法，以实现可变形结构的实时边界更新。利用两种溃坝冲击工况——纯流体条件和含颗粒条件下涉及柔性板的溃坝冲击——对拟议框架的可靠性进行了评估。对于这两种测试构型，预测的波前和颗粒前锋演化、多相流动模式和柔性板变形与实验测量结果表现出令人满意的一致性。这些对比展示了该框架在考虑几何形状、边界条件和短持续时间冲击载荷下再现耦合响应能力。

尽管具备这些能力，当前方法仍面临各种需要进一步研究的局限性。例如，在结构建模方面，柔性结构的网格生成主要依赖四面体单元，缺乏单元类型的多样性。此外，流体-结构耦合算法要求流体粒子和结构网格具有可比的空间分辨率。因此，在模拟极薄的柔性结构时，流体域分辨率必须按比例细化，不可避免地导致总体计算成本显著增加。此外，当前的SPH-DEM耦合采用完全解析策略，在处理大量颗粒时会产生过高的计算成本。未来研究将专注于解决这些挑战，以进一步提高该框架在涉及柔性结构的大规模工程多相流模拟中的计算效率和适用性。

CRediT作者贡献声明

Yu Yihang：撰写-审阅与编辑，撰写-原稿，可视化，验证，软件，方法论，调查，形式分析，数据整理。Du Jiahang：撰写-审阅与编辑，验证，软件，调查，形式分析，数据整理。Zhang Jiabao：验证，数据整理。Ma Huaqing：撰写-审阅与编辑，项目管理，资金获取，概念化。Zhao Yongzhi：监督，软件，资源，项目管理，资金获取。

Yu Yihang | Du Jiahang | Zhang Jiabao | Ma Huaqing | Zhao Yongzhi

浙江大学能源工程学院先进装备研究所，杭州 310027，中国