用于海底滑坡的交错网格显式-隐式两相两点方法（MPM）

《INTERNATIONAL JOURNAL OF MECHANICAL SCIENCES》：Staggered-grid explicit-implicit two-phase two-point MPM for submarine landslides

　　洪峰|尹镇宇
香港理工大学土木与环境工程系，香港九龙红磡

摘要
准确模拟海底滑坡对于地质灾害评估至关重要。两相两点物质点法（MPM）通过将耦合的土壤-水系统离散为两组相互作用的物质点，被广泛用于此类问题。然而，一个稳健的MPM框架仍在研究中。本研究提出了一种新颖的交错网格显式-隐式两相两点MPM模型。数学模型采用了两相流动公式，其中土壤用规范化的μ(I)非牛顿材料表示，水用牛顿材料表示。两组拉格朗日物质点用于离散土固相和液相，同时采用共享的欧拉交错网格，从而同时减轻了虚假压力模式，并实现了简单有效的土壤-水耦合方案。该解决方案在显式-隐式MPM框架内采用了算子分裂策略。在第一个预测阶段，显式计算土固相和液相的中间速度。第二个校正阶段隐式求解压力泊松方程并校正中间速度场。所提出的MPM模型进一步应用于一系列基准案例，包括盖驱动空腔流、浸没颗粒柱坍塌和海底滑坡。MPM结果与其他数值研究的结果以及实验测量结果具有良好的一致性，证明了所开发的MPM框架描述海底滑坡的能力。最后，将新方法与现有的两相两点MPM公式进行了比较，为选择适当的海底滑坡数值模型提供了实用指导。

1. 引言
海底滑坡是沿海底形成的斜坡沉积物的大规模重力驱动的质量运动。这些事件直接威胁到海底基础设施，并可能产生破坏性海啸，从而将风险扩展到沿海设施和人类安全[[1], [2], [3], [4], [5]]（见图1）。一个著名的例子是1963年10月9日瓦容特水库的海底滑坡。大约3亿立方米的物质位移产生了80米高的海啸，冲破了大坝，摧毁了隆加罗内镇，导致近两千人死亡[6]。还有许多其他由滑坡引起的海啸事件被记录在案，包括1936年挪威的Norddal峡湾滑坡[7]、1969年密西西比河三角洲事件[8]、1998年巴布亚新几内亚海底滑坡[9]、利图亚湾滑坡[10]和斯特龙博利岛滑坡[11]。因此，准确预测滑坡运动对于灾害评估和风险缓解至关重要。

下载：下载高分辨率图片（1004KB）
下载：下载全尺寸图片
图1. 海底滑坡和两相模型。(a) 具有复杂土壤-水相互作用的海底滑坡；(b) 使用体积平均技术将海底滑坡表示为两相模型。

为了实现上述目的，已经开发了几种数值方法。传统方法，包括极限平衡法（例如[12]）、极限分析法（例如[13]）和有限元法（例如[14]），被广泛用于评估滑坡稳定性和计算安全系数。然而，由于这些方法在小变形框架内运行，它们无法模拟滑坡失败后的行为。特别是，传统的FEM在变形较大时会出现严重的网格变形，导致数值不稳定和模拟终止。这一根本限制推动了先进的大变形数值方法的发展，这些方法能够模拟从开始到滑坡失败整个过程。

在大变形模拟框架内，采用了两种不同的海底滑坡概念框架。第一种是单相模型（例如[[15], [16], [17]]），将滑坡视为单一的均质化相。这种方法的控制方程从深度平均公式（如浅水或Boussinesq型方程）到更复杂的系统（如势流或完整的Navier-Stokes方程）不等。虽然这种方法可以以相对较低的计算成本有效模拟滑坡运动，但它不考虑土壤-水相互作用项。相反，两相土壤-水耦合框架将土壤骨架（在本研究中称为“土固相”）和水视为不同的但相互作用的相（例如[18]）。尽管比单相框架需要更多的计算成本，但两相框架可以更准确地捕捉土壤-水相互作用，从而提供更优越的滑坡演变和波浪传播表示。

基于用于土固相和液相的运动学描述（欧拉或拉格朗日），海底滑坡的两相数值模型可以分为三种主要类型：欧拉-欧拉、拉格朗日-拉格朗日和欧拉-拉格朗日（参见[19]的综述）。欧拉-欧拉方法在固定的欧拉参考框架内描述土固相和液相。它通常使用计算流体动力学（CFD）技术实现，主要是有限体积法（FVM）（例如[[20], [21], [22], [23]]）和有限差分法（FDM）（例如[24]）。虽然这种方法为模拟大变形流动提供了数值稳定性，并可以利用现有的商业或开源CFD软件（例如Fluent、OpenFOAM），但在准确捕捉演变的土壤-水界面和水自由表面方面仍存在重大挑战。精确跟踪这些界面对于可靠的海啸波预测至关重要。

拉格朗日-拉格朗日方法将土壤和水离散为移动的拉格朗日点或颗粒。该框架包括物质点法（MPM）（例如[[25], [26], [27], [28], [29], [30]]）、平滑粒子流体动力学（SPH）（例如[[31], [32], [33], [34]]）和粒子有限元法（PFEM）（例如[35,18]）等技术。其主要优势在于它能够处理大变形、自由表面和演变的土壤-水界面而不会遇到网格变形。作为一种混合框架，欧拉-拉格朗日方法结合了两种描述方法的优点。它通常使用拉格朗日方法（例如MPM、SPH）描述滑坡体，并结合欧拉求解器处理周围的水（例如[36,37]）。这种组合利用了高效的土壤建模和稳健的流体模拟，但在开发耦合算法时引入了挑战，通常导致相当大的编码复杂性[[38], [39], [40]]。

在MPM社区中，早期的研究使用两相一点公式模拟海底滑坡中的饱和多孔介质，其中一层物质点代表固液混合物[[41], [42], [43]]。随后，开发了两相两点MPM，其中两组独立的物质点分别代表土壤骨架和孔隙水，并在背景计算网格上解决它们的相互作用[[44], [45], [46]]。由于其捕捉大变形和土壤-水耦合的稳健性，这种两相两点MPM已被广泛用于模拟海底滑坡[[47], [48], [49]]。在两相两点框架内，已经开发了两种主要变体（见表1）：

(i) 显式MPM框架：早期的基础研究使用MPM求解Biot混合物模型（即耦合土壤变形与孔隙流体流动的饱和多孔系统），从而建立了第一个两相MPM公式[44,50,41,45,42,43]。在早期的两相两点模型[44,50,45,46]中，使用显式时间积分的MPM求解土壤骨架（土固相），其力学响应由速率无关的本构模型（例如Mohr-Coulomb或Drucker-Prager）描述（例如[52]）。同时，液体相也使用显式方案求解，使用状态方程描述水压。这种方案的关键特点是土壤和液体相在共享的计算背景网格上离散。因此，网格节点同时存储两种成分的速度场，允许在统一的数值框架内直接计算土壤-水耦合机制。

表1. 现有用于土壤-水耦合问题的两相两点MPM公式总结。

参考文献
固体相 | 流体相 | 空单元 | 网格 | 积分 | 本构模型 | 网格 | 积分 | 可压缩性
[44] | 共位 | 显式 | MCC | 共位 | 显式 | 可压缩
[50] | 共位 | 显式 | MCC | 共位 | 显式 | 可压缩
[51] | 共位 | 显式 | MCC | 共位 | 显式 | 可压缩
[53] | 共位 | 显式 | DPC | 共位 | 显式 | 可压缩
[54] | 共位 | 显式 | 应变软化DPC | 共位 | 显式 | 可压缩
[55] | 共位 | 显式 | DPC | 共位 | 显式 | 可压缩
[56] | 共位 | 显式 | MCC | 共位 | 显式 | 可压缩
[26] | 共位 | 显式 | 非结合 | MCC | 共位 | 显式 | 可压缩
[25] | 共位 | 显式 | 弹性塑性-μ(I) | 共位 | 显式 | 可压缩
[57] | 共位 | 显式 | DPS | 半交错 | 显式-隐式 | 不可压缩
[58] | 共位 | 显式 | DP | 半交错 | 显式-隐式 | 不可压缩
[27] | 共位 | 显式 | DP | 交错 | 显式-隐式 | 不可压缩

在其最初开发之后，显式MPM框架已成功扩展到模拟各种岩土工程问题，包括滑坡引发的海啸[[54], [55], [56]]、溢坝[[59]]以及非饱和土壤中的斜坡失稳[[26]]。尽管这些研究表明了该方法在这些应用中的有效性，但现有的公式面临两个主要数值挑战。首先，使用近似不可压缩流体方程可能导致孔隙压力不稳定[60]。其次，依赖显式时间积分带来了严格的稳定性约束，需要较小的时间步长，从而导致较高的计算成本[61,62]。

(ii) 混合网格MPM框架：为了解决上述完全显式MPM框架中弱可压缩流体相的孔隙压力不稳定问题，开发了一种混合网格两相两点MPM框架（例如[[27], [57], [58]]。在这种混合框架中，土壤固相使用显式MPM方案和共位网格求解，而不可压缩流体相则使用交错或半交错网格的显式-隐式MPM公式处理。

与完全显式MPM框架类似，这种混合网格MPM框架中的土壤固相仍然使用Mohr-Coulomb和Drucker-Prager等本构模型显式建模。不同之处在于处理流体相不可压缩性的方法，可以考虑几种不同的解决方案（见图2）：
(i) 共位网格类型：所有速度分量和压力存储在相同的网格节点位置（例如[[63], [64], [65], [66], [67]]）。然而，当在这种共位网格上使用基于分数步长的显式-隐式方案时，需要稳定性技术（例如Rhie-Chow插值）来避免虚假的棋盘式压力模式。对于这种网格类型的整体隐式方案，通常需要额外的稳定技术（例如[63,64]）。
(ii) 交错网格类型：速度和压力存储在空间偏移的位置，以本质上增强数值稳定性。遵循经典的MAC（标记和单元）方法，水平速度通常存储在垂直单元面上，垂直速度存储在水平单元面上，压力存储在单元中心（例如[[68], [69], [70], [71], [72]]。
(iii) 半交错网格类型：这种方法在单元节点存储水平和垂直速度分量，在单元中心存储压力（例如[73]）。

基于上述不可压缩流体相解决方案，开发了两种主要类型的混合网格MPM框架：
(a) 半交错网格MPM框架（例如[57,58]）：例如，Gao等人[57]将半交错网格不可压缩MPM与共位网格显式MPM用于土壤耦合。不可压缩流体使用分数步长方法模拟，包括速度预测步骤和速度校正步骤。使用交错网格存储速度。土壤使用显式MPM建模，速度存储在共位网格上。由于土壤和水速度存储在相同的网格节点上，这种方法能够实现有效的土壤-水耦合。然而，半交错网格公式的一个显著限制是它容易受到沙漏模式的影响，表现为虚假的速度振荡，通常需要额外的数值阻尼来进行稳定。
(b) 交错网格MPM框架（例如[[27], [74], [75]]：在后续发展中，He等人[27]将交错网格不可压缩MPM与共位网格显式MPM用于固体耦合。虽然这种公式有效解决了压力不稳定问题，但它引入了显著的几何复杂性：每个相的离散变量存储在不同的网格位置。具体来说，固相速度定义在网格节点上，液相速度在单元面上，水压在单元中心。这种状态变量在计算网格上的分散使得直接计算土壤-水相互作用机制变得复杂。尽管使用插值技术在这些位置之间传递压力信息，但可能会导致物理上不一致的压力场。主要原因是固相压力是通过对基于网格的插值得到的，其对固相速度的影响也是通过插值施加的，而流体压力处理方式不同，不依赖于此过程。因此，一个稳健的两相两点MPM框架仍在开发中。

在这项工作中，开发了一种新的海底滑坡两相两点MPM框架。首先，我们采用两相流动公式，而不是Biot混合物理论，其中土固相和液相都建模为不可压缩流体：土固相遵循规范化的μ(I)流变学，液相遵循牛顿流变学。然后使用交错网格显式-隐式两相两点MPM求解该系统。土固相和液相分别用不同的拉格朗日物质点层离散，而两个相共享一个固定的欧拉交错网格。两个阶段都使用相同的Chorin型算子分裂方案来求解。该算法包括：(i) 显式预测步骤，通过不考虑压力项求解动量方程来获得中间速度；(ii) 隐式校正步骤，利用计算出的压力梯度更新速度以确保不可压缩性。通过这种方式，土壤-固体相和流体相共享相同的动量解框架、错位网格和时间积分算法。因此，我们将这种方法称为“错位网格显式-隐式两相两点MPM”。与我们之前的工作[76]相比，之前的工作主要开发和分析了一维干燥颗粒流动的新流变模型，而本研究提出并验证了用于海底滑坡的新两相两点MPM框架。与现有的两相两点MPM公式相比，所建立的框架具有四个主要优点。首先，它通过将流体视为不可压缩的来获得稳定的流体场（例如，避免压力不稳定性），而不是依赖于近似不可压缩的状态方程。其次，它采用了一个简单且经过良好验证的土壤本构模型，易于在MPM中实现，避免了现有两相两点框架中所需的复杂应力积分算法。第三，它导致了两相两点MPM的简单实现，因为土壤-固体相和流体相共享相同的动量解程序、相同的错位网格和相同的时间步长。最后，与替代公式相比，它降低了计算成本，因为将土壤和流体相都描述为类流体流变材料允许更大的允许时间步长。据作者所知，这是首次提出这种类型的两相两点MPM框架用于水-土相互作用问题。

本文的结构如下。第2节讨论了海底滑坡的两相流动数学模型。第3节详细介绍了所提出的错位网格显式-隐式两相两点MPM公式。第4节提供了数值实例来说明MPM框架的有效性。第5节提供了分析和讨论，第6节给出了结论。

2. 两相流动数学模型
描述海底滑坡及其产生的波浪的物理过程的数学模型对于建立合理的数值框架至关重要。本研究采用了已被广泛用于海底滑坡的两相流动数学模型[[77], [78], [79], [80], 16]。在本节中，将讨论该数学模型的基本假设、守恒方程和应力框架。

2.1 基本假设
在介绍具体公式之前，首先明确以下基本假设：
A1：土壤被视为由土壤-固体部分和流体部分（即本研究中的水）组成的饱和多孔介质；忽略孔隙中的空气存在。
A2：海底滑坡和产生的波浪被建模为使用体积平均技术的两相流动系统。图1说明了海底滑坡体积平均技术的基本原理。在这种方法中，土壤和水被平均为连续体，土壤-固体相和水的体积分别表示为：
(1)Vˉs = nsVs, ρˉs = nsρs
(2)Vˉw = nwVw, ρˉw = nwρw
其中下标s和w分别表示土壤-固体相和水相。对于每个相，V表示体积，ρ表示固有密度，n表示体积分数。平均后的对应量分别用Vˉ和ρˉ表示。
A3：假设土壤-固体相和水相都是不可压缩的。
A4：假设土壤-固体相和水相的体积分数是恒定的。这与传统的饱和多孔介质模型不同，在传统模型中孔隙率可以变化。然而，在临界状态的小变形范围内，这个假设是有效的。在临界状态的大变形下，松散砂和密实砂都不会进一步改变体积，因此体积分数假设是合理的。
A5：土壤-固体相和水相之间没有质量传递。
A6：通过固体-流体相互作用项考虑土壤-固体相和水相之间的动量交换。
A7：流体相被建模为不可压缩的牛顿流体。
A8：土壤-固体相被建模为不可压缩的非牛顿材料；其行为不由传统的速率无关的弹塑性材料描述（见第2.4小节）。

2.2 质量守恒方程
连续性方程强制质量守恒。这个原则在不同学科中有不同的含义。在固体力学中，变形体的质量是恒定的。在流体力学中，控制体积内的质量流入量减去质量流出量等于体积内的质量变化率。这个关系表示为：
(3)DmDt = ∫∫∫V(DρDt + ρ(?·u))dV = 0
其中DmDt是物质导数；ρ是密度；u是速度；V是控制体积。
进一步假设密度梯度为零，即?ρ=0，上述方程简化为：
(4)DρDt + ρ(?·u) = ?ρ?t + (?ρ)·u + ρ?·u = ?ρ?t + ρ?·u = 0
在两相流动模型中，土壤-固体相和流体组分都满足上述连续性方程。因此，可以使用体积平均密度写出这两个相的质量守恒方程：
(5)?(nsρs)?t + ?·(nsρsus) = 0
(6)?(nwρw)?t + ?·(nwρwuw) = 0
其中us和uw分别是土壤-固体相和流体相的速度。
假设土壤-固体相和流体相都是不可压缩的且体积分数恒定，我们得到以下质量守恒方程：
(7)ns?·us = 0
(8)nw?·uw = 0

2.3 动量守恒方程
流体的动量方程源自柯西动量方程，进而得到众所周知的纳维-斯托克斯动量方程，表示为：
(9)DuDt = 1ρ?·σ + b
其中σ是柯西应力张量；b是作用在连续体上的加速度，在这里我们只考虑重力。
上述方程用欧拉形式表示。忽略物质导数中的对流项后，可以用拉格朗日形式重写，这在固体力学中广泛使用：
(10)ρ?u?t = ?·σ + ρb
在两相流动模型中，土壤-固体相和流体组分都满足上述动量方程。采用拉格朗日形式并使用体积平均密度和体积平均应力，这两个相的动力量守恒方程可以写为：
(11)nsρs?us?t = ?·σˉs + nsρsb
(12)nwρw?uw?t = ?·σˉw + nwρwb
其中σˉs和σˉw分别是土壤-固体相和流体相的平均应力张量。
当进一步考虑土壤-固体相和流体相之间的动量交换时，上述动量方程变为：
(13)nsρs?us?t = ?·σˉs + nsρsb + fd
(14)nwρw?uw?t = ?·σˉw + nwρwb ? fd
其中fd是固体-流体相互作用项，即阻力。
对于阻力模型，可以使用线性达西关系或非线性达西-福希海默关系。达西-福希海默关系在高流速下提供了更准确的固体-流体相互作用力的描述。然而，本研究采用达西定律来简化当前的MPM公式。未来的工作可能将达西-福希海默关系纳入所提出的MPM模型中。达西定律表示为：
(15)?p = ηwκq
其中ηw表示流体的粘度；κ和q表示绝对渗透率和流量；表示为：
(16)κ = d21/50nw3(1?nw)2
(17)q = nw(uw ? us)
其中d是土壤-固体颗粒直径。
压力降?p代表固体-流体相互作用力fd。结合方程(15)-(17)，我们可以得到以下具体方程：
(18)fd = ηwnw2/150d2(1?nw)2(uw ? us)

2.4 应力框架
在两相流动数学模型[[81], [82], [83]]中，土壤骨架和流体相的应力分别表示为：
(19)σˉs = nsσs
(20)σˉw = nwσw
其中σs和σw分别是土壤-固体相和水相的柯西应力张量；每个张量包含一个各向同性和一个偏应力分量：
(21)σs = ?psI + τs
(22)σw = ?pwI + τw
其中ps和τs分别是土壤-固体相的各向同性和偏应力张量；pw和τw分别是水相的各向同性和偏应力张量。这四个术语将在下文中讨论。
流体相的本构关系包括各向同性和偏应力分量。各向同性应力对应于水压力，通过将水视为不可压缩并满足无散度速度场来隐式求解。因此，只有粘性偏应力需要本构建模，由牛顿流变模型提供：
(23)τw = ηw(?uw + (?uw)T
土壤-固体相的本构关系可以归类为三类主要模型：(i) 类固体型（例如Mohr-Coulomb或Drucker-Prager）[84,85]，其中体积分数是可变的。这些模型将各向同性应力与体积分数相关联，以捕捉临界状态前的收缩和膨胀。它们准确描述了速率无关的小应变行为。(ii) 类流体型（例如μ(I)模型），其中体积分数是恒定的。它们假设颗粒介质在达到临界状态时保持恒定体积，并使用有效的粘性流变来描述偏应力。(iii) 相变型（Feng等人的Elastoplasticity-μ(I)模型[25]），它将颗粒介质视为体积可变的非线性材料。这些模型在小应变下表现出类固体行为，在大变形下转变为类流体行为。由于本研究假设体积恒定，因此不能使用类固体模型或相变模型。只采用类流体模型。
接下来是广泛使用的μ(I)模型[86,87]，以及我们团队[76]开发的正则化μ(I)流变模型。在这个模型中，考虑无散度速度场来计算各向同性应力，粘性应力表示为：
(24)τs = ηsγ˙ = ηs(?us + (?us)T
(25)ηs = μ(I)|γ˙|ps = (μs + μd ? μsI0 + II)ps|γ˙|
其中ηs是土壤-固体相的表观粘度；剪切率张量表示为γ˙，其二阶不变量表示为|γ˙|。μd和μs分别是动态和静态摩擦系数。
由于模型在低剪切率下粘度变得无限大，因此直接应用该模型会导致数值不稳定。因此，我们通过将表观粘度分解为静态和增量部分来正则化模型：
(26)ηs = ηsta + ηΔ
(27)ηsta = psμs|γ˙|
(28)ηΔ = ps(μd ? μs)I0ps/ρg + dg|γ˙|
其中ηsta和ηΔ分别是土壤-固体相的静态和增量表观粘度。
使用L1正则化技术[76]对静态表观粘度项进行正则化，可以表示为：
(29)ηsta1 = psμs|γ˙| + λ
其中λ是正则化参数，在本研究中的值为1 × 10??，基于之前的敏感性分析[76]。
上述两相流动模型的应力框架与Biot的混合理论不同，后者的讨论见附录A。另外，请注意，上述正则化的μ(I)模型适用于颗粒材料。对于粘土土壤，应使用适当的流变模型（例如Bingham或Herschel-Bulkley）。

3. 错位网格显式-隐式两相两点MPM
本节主要讨论错位网格显式-隐式两相两点MPM，包括算子分裂方法、显式-隐式MPM算法、显式预测步骤和隐式校正步骤。

3.1 算子分裂方法
采用Chorin的算子分裂方法来求解上述两相流动数学模型[88,76]。该算法包括一个预测步骤，通过不考虑压力求解动量方程来计算中间速度；以及一个校正步骤，使用计算出的压力梯度更新速度以确保不可压缩性。
第2节中提出的数学模型可以使用算子分裂方法重新表述（见图3），其步骤如下：
(i) 预测步骤：不考虑压力场ps和pw，方程(13)-(14)表示为：
(30)nsρsus* ? ustΔt = ?·τs + nsρsb + fd
(31)nwρwuw* ? uwtΔt = ?·τw + nwρwb ? fd
其中us*和uw*分别是土壤-固体相和水相的中间速度。
(ii) 校正步骤：然后使用未知的压力场（即方程(32)-(33)）校正上述中间速度。方程(32)-(33)和方程(30)-(31)的结合形式恢复了第2节中提出的动量方程。
(32)nsρsust + 1 ? us*Δt = ??pst + 1
(33)nwρwuwt + 1 ? uw*Δt = ??pwt + 1
其中ustin*和uwt*分别是土壤-固体相和水相的最终速度。

同时，使用第2.2小节中的质量守恒方程，并将方程(7)-(8)的梯度代入方程(32)-(33)的梯度中，我们得到以下泊松压力方程：
(34)Δtnsρs?2pst + 1 = ?·us*
(35)Δtnwρw?2pwt + 1 = ?·uw*

3.2 显式-隐式MPM算法
上述数学模型的数值解需要考虑空间离散框架、时间积分方案和网格类型。采用两相两点MPM作为空间框架。在这种方法中，土壤-固体相和水相由两组独立的拉格朗日物质点离散化，这些点在静止网格内移动。与MPM概念一致，所有物理状态变量都由物质点携带，而计算网格仅作为求解控制方程的临时支架。信息通过形状函数在点与网格之间传递，类似于SPH中的核近似。
关于时间积分，预测-校正算法内同时使用了显式和隐式方案。预测步骤采用显式MPM公式来求解两个相的动量方程(30)-(31)。随后，校正步骤使用隐式MPM公式来解决与压力场相关的控制方程（公式（32）-（35））。在求解压力泊松方程时，对于网格类型，采用相同的交错（MAC）网格排列，以消除虚假压力模式。在这种统一设置中，土壤-固体和土壤相的压力定义在单元中心，而它们的相应速度定义在单元边界。这种交错排列在速度预测和压力求解阶段始终得到应用，从而提高了计算 robustness 和效率。基于上述考虑，这种方法构成了交错网格显式-隐式双相两点MPM。3.3. 显式预测步骤预测步骤使用显式时间积分。基本弱形式公式和MPM离散化在下面讨论。采用Galerkin方法将公式（30）和（31）的强形式转换为弱形式，表示如下：（36）∫Ωsnsρsas*·δusdΩ=?∫Ωsτs:?δusdΩ+∫Ωs(nsρsbs+fd)·δusdΩ+∫ΓNstn·δusdΓ（37）∫Ωwnwρwaw*·δuwdΩ=?∫Ωwτw:?δuwdΩ+∫Ωw(nwρwbw?fd)·δuwdΩ+∫ΓNwtn·δuwdΓ其中Ωs和Ωw分别是土壤-固体和流体相的体积；δus和δuw分别是土壤-固体和流体相的虚拟测试函数；这些函数在Dirichlet边界处为零。上述弱形式方程的离散化在标准MPM框架内进行。在此框架中，连续体被表示为一组粒子，表示如下：（38）Ωs=∑ps=1NpsVps（39）Ωw=∑pw=1NpwVpwin其中Ωs和Vps分别是土壤-固体相的连续体和粒子的体积，Ωw和Vpw对水相具有相同的意义。网格的物理信息从材料点获得，表示如下：（40）□Ia,st=∑spNIa,spt□spt（41）□Ia,wt=∑wpNIa,wpt□wpt其中□spt和□wpt分别代表土壤-固体和流体粒子的物理信息；□Ia,st和□Ia,wt是网格的物理信息。NIa,spt是土壤-固体粒子与网格之间的形状函数；NIa,wpt是流体粒子与网格之间的形状函数。在这项工作中，选择了二次B样条形状函数，表示如下：（42）N(x)={34?x2,0≤|x|<1212(32?|x|)2,12≤|x|<320,32≤|x|其中x表示粒子-网格距离；N(x)表示形状函数的值。根据上述MPM离散化方案，显式预测步骤包括以下阶段（见图4）：(a) 更新粒子位置图4. 交错网格显式-隐式双相两点MPM。首先在每个步骤更新粒子位置，以确保速度场无散度（见[76]）。(43) xspt+1=xspt+uspt+1Δt(44) xwpt+1=xwpt+uwpt+1Δt其中xspt+1和xspt分别是时间t + 1和t时的土壤-固体粒子的位置；uspt+1是时间t + 1时的土壤-固体粒子的速度；xwpt+1和xwpt分别是时间t + 1和t时的水粒子的位置；uwpt+1是时间t + 1时的水粒子的速度。(b) 粒子到网格粒子信息通过形状函数插值传输到网格。对于映射过程（即P2G和G2P），交错网格和共位网格之间的区别在于：（i）在交错网格中，信息在单元边界和粒子之间传输，而在共位网格中，信息在网格节点和粒子之间传输；（ii）在交错网格中，速度分量根据各自的坐标方向进行映射，而在共位网格中没有考虑这种方向性分离。因此，土壤-固体相的质量和速度在a方向上映射到网格单元中，表示如下：(45)mIa,st=∑spNIa,sptmsp(46)uIa,st=∑spNIa,sptmspusptmIa,st其中msp和uspt分别是时间t时固体相的质量和速度；mIa,st和uIa,st表示时间t时沿a方向的单元边界的质量和速度，a表示边界方向，I表示单元；NIa,spt表示形状函数，通过NIa,spt=N(xIa?xspt)计算得到。类似地，流体单元的质量和速度表示如下：(47)mIa,wt=∑wpNIa,wptmwp(48)uIa,wt=∑wpNIa,wptmwpuwptmIa,wt(c) 计算粒子内部应力土壤-固体和流体粒子的速度梯度通过公式（49）和（50）计算。注意，在计算内部应力之前必须强制网格上的边界条件。(49)?uspt=∑a(?usp,at)=∑a(∑I?NIa,sptuIa,st)(50)?uwpt=∑a(?uwp,at)=∑a(∑I?NIa,wptuIa,wt)其中?uspt和?uwpt分别表示土壤-固体粒子和水粒子的速度梯度；?usp,at和?uwp,at分别表示在某个方向上土壤-固体粒子和水粒子的速度梯度。土壤-固体相的粘性剪应力通过规则化μ(I)模型计算（见第2节），表示如下：(51)ηspt=f(uspt)(52)τspt=ηspt(?uspt+(?uspt)T)其中f(uspt)表示本构关系。流体相的粘性剪应力通过第2节讨论的牛顿流变模型计算，表示如下：(53)τwpt=ηw(?uwpt+(?uwpt)T)(d) 计算中间速度中间速度基于内力、外力和阻力计算。土壤-固体和流体相的内力表示如下：(54)fIa,sint,t=?∑spτspt·?NIa,sptVspt(55)fIa,wint,t=?∑wpτwpt·?NIa,wptVwpt其中fIa,sint,t和fIa,wint,t分别是土壤-固体和水网格的内力。土壤-固体和水网格的外力只考虑重力，表示如下：(56)fIa,sext,t=∑spNIa,sptmsptb(57)fIa,wext,t=∑wpNIa,wptmwptb其中fIa,sext,t和fIa,wext,t分别是土壤-固体和流体网格的外力。阻力更新如下：(58)fIa,swd,t=ηw(nIa,wt)2150d2(1?nIa,wt)2(uIa,wt?uIa,st)其中fIa,swd,t是阻力；nIa,wt是网格节点的水体积分数。使用以下标准确定是否应在网格上计算固-流相互作用力：(59)mIa,st>0,mIa,wt>0由于使用相同的交错网格，土壤-固体和水的网格速度使用以下公式计算：(60)aIa,s*=fIa,sint,t+fIa,sext,t+fIa,swd,tmIa,st(61)aIa,w*=fIa,wint,t+fIa,wext,t?fIa,swd,tmIa,wt(62)uIa,s*=uIa,st+aIa,s*Δt(63)uIa,w*=uIa,wt+aIa,w*Δt其中aIa,s*和aIa,w*分别是土壤-固体和水节点的中间加速度。3.4. 隐式校正步骤本小节介绍了隐式步骤，包括求解压力泊松方程和计算最终速度。注意，压力泊松方程使用有限差分方法离散化和求解，可以视为嵌入在整个MPM框架内的有限差分框架。首先，土壤-固体相的压力泊松方程表示如下：(64)(?2pst+1)(i,j)=ps,(i+1,j)t+1+ps,(i?1,j)t+1?2ps,(i,j)t+1(Δx)2+ps,(i,j+1)t+1+ps,(i,j?1)t+1?2ps,(i,j)t+1(Δy)2(?·us*)(i,j)=us,(i+0.5,j)*?us,(i?0.5,j)*Δx+vs,(i,j+0.5)*?vs,(i,j?0.5)*Δy其中ps,(i,j)t+1是单元(i, j)的固体相压力；ps,(i+1,j)t+1, ps,(i?1,j)t+1, ps,(i,j+1)t+1和ps,(i,j?1)t+1是包围该单元的土壤-固体相压力；us,(i+0.5,j)*, us,(i?0.5,j)*, vs,(i,j+0.5)*, 和 vs,(i,j?0.5)*代表包围该单元的中间速度。类似地，流体相的压力泊松方程表示如下：(65)(?2pwt+1)(i,j)=pw,(i+1,j)t+1+pw,(i?1,j)t+1?2pw,(i,j)t+1(Δx)2+pw,(i,j+1)t+1+pw,(i,j?1)t+1?2pw,(i,j)t+1(Δy)2(?·uw*)(i,j)=uw,(i+0.5,j)*?uw,(i?0.5,j)*Δx+vw,(i,j+0.5)*?vw,(i,j?0.5)*Δy将上述离散方程应用于整个域，得到椭圆型偏微分方程，表示如下：(66)AI,spI,st+1=UI,s(67)AI,wpI,wt+1=UI,ww其中AI,s、pI,s和UI,s分别是从方程（64）综合得到的固体相的系数矩阵、未知压力向量和中间速度向量；AI,w、pI,w和UI,w分别是从方程（65）综合得到的水相的系数矩阵、未知压力向量和中间速度向量。采用共轭梯度方法求解上述椭圆型偏微分方程。随后，使用两个步骤过程校正中间速度：(f) 校正中间速度，和 (g) 网格到粒子。(f) 校正中间速度我们可以使用以下方程得到最终速度：(68)uIa,st+1=uIa,s*?Δtnsρs?pIa,st+1(69)uIa,wt+1=uIa,w*?Δtnwρw?pIa,wt+1(g) 网格到粒子在将信息从网格传输到粒子之前，必须强制网格上的边界条件。然后，使用以下插值方程获得粒子速度。(70)usp,at+1=usp,at+∑I(uIa,st+1?uIa,st)NIa,spt(71)uwp,at+1=uwp,at+∑I(uIa,wt+1?uIa,wt)NIa,wpt其中usp,at+1和usp,at分别是时间t + 1和t时的固体粒子速度；uwp,at+1和uwp,at分别是时间t + 1和t时的水粒子速度。总结来说，阶段(a)-(g)定义了显式-隐式MPM循环。每个时间步骤都从更新粒子位移和重置背景网格开始（即，所有节点速度信息被清除，而网格保持不变）。临界时间步长的计算在附录B中给出。此外，采用GPU加速技术来提高计算效率，详见附录C。在这项研究中，考虑了Neumann边界条件，即?pIa,wt+1·n=0,?pIa,st+1·n=0（其中n是单位向外法向量）。同时，使用基于网格的技术施加无滑移和无穿透条件，即所有固体边界上的水平和垂直速度都设为零[76,27]。4. 数值示例在本节中，进行了许多2D数值模拟，以验证所提出的MPM框架，重点是其模拟纯流体流动和耦合土壤-水相互作用的能力。土壤相求解器的验证已在之前的研究中单独进行（见[76]）。因此，本研究评估了流体求解器和耦合算法。选定的基准案例及其具体的验证目标总结在表2中。所有示例的计算设置是一致的：使用方形背景网格，共轭梯度迭代求解器的收敛容忍度设置为1 × 10?12，并使用我们自家的MPM代码在配备Intel Core i9–13,900 K处理器和NVIDIA GeForce RTX 4080 GPU的工作站上执行模拟。表2. 本研究或我们之前研究中数值示例的总结。方面性能验证参考文献土壤模拟密集颗粒流[76]流体模拟不可压缩流4.1 肖节 土壤-水耦合问题模拟海底滑坡4.2 肖节 土壤-水耦合问题模拟海底滑坡和产生的波浪4.3 肖节 土壤-水耦合问题模拟海底滑坡和产生的波浪4.44.1. 情况1：在盖驱动腔体中的牛顿流4.1.1. 模型描述首先使用经典的盖驱动腔体基准来验证牛顿流的流体相MPM求解器。在这种模拟中，流体最初在方形腔体内由于顶部边界施加的水平速度而形成涡流。图5显示了MPM模型配置，对应于Ghia等人[89]的经典盖驱动腔体基准和Wang等人[90]的最新设置。腔体为1 m × 1 m的正方形。左侧、右侧和底部墙壁施加无滑移边界条件。对于顶部墙壁，施加水平速度ulid，并将法向自由度设置为零。域使用矩形元素的结构化网格离散化。网格尺寸的选择遵循Wang等人[90]之前的研究。尽管Wang等人[90]的网格敏感性分析表明dx=0.024 m可以得到满意的结果，但我们采用更细的分辨率dx=0.005 m以确保更高的精度。每个元素初始化四个材料点，共计160,400个材料点。时间步长设置为10–3 s，满足所需的临界时间步长（见附录B），以保证计算稳定性。为了达到雷诺数Re=400 [89,90]，流体参数设置如下：密度ρw=1000 kg/m3，盖速度ulid= 0.4 m/s，动态粘度ηw= 1 Pa·s；模拟时间设置为80 s。对于另一个雷诺数1000，流体参数设置如下：密度ρw=1000 kg/m3，盖速度ulid=1.0 m/s，动态粘度ηw= 1 Pa·s。模拟时间设置为50 s。4.1.2. 模型结果图6首先展示了Re=400时的速度场演变。从t = 5 s到20 s，形成一个主涡流，其核心从右上区域向腔体中心迁移。随后，从t = 20 s到80 s，涡流位置在中心稳定，进一步移动很小。这种涡流的初始形成、迁移和最终稳定与Ghia等人[89]和Wang等人[90]的基准结果非常吻合。这种一致性成功验证了当前MPM求解器模拟不可压缩牛顿流的能力。图7进一步提供了定量验证，比较了MPM模拟和建立的基准数据的标准化水平和标准化垂直速度剖面。沿着水平中心线（y = 0.5L）和垂直中心线（x = 0.5L）的标准化垂直速度与Ghia等人[89]和Wang等人[90]的SPH/FVM模拟的结果进行了对比。如图7(a)和7(b)所示，所有数据集之间的密切一致性验证了当前不可压缩MPM公式对于牛顿流的准确性和有效性。在Re=400时，对盖子驱动的方形腔体流动的MPM结果、SPH结果和FVM解进行比较。(a) y = 0.5L处的归一化垂直速度；(b) x = 0.5L处的归一化水平速度。图8显示了Re=1000时MPM情况从t = 2.5 s到t = 40 s的流线图。主涡旋从右上区域向腔体中心移动。大约在t ≈ 20 s后，涡旋结构逐渐稳定，形成了一个主导的中心涡旋和两个较小的角涡旋。图9将MPM结果与参考解[89,91]进行了比较。这种高度的一致性证明了所提出的MPM模型在捕捉不可压缩腔体流动方面的有效性。下载：下载高分辨率图片（2MB）下载：下载全尺寸图片

图8. Re=1000时盖子驱动的方形腔体流动的MPM结果流线的演变。(a) t = 2.5 s；(b) t = 5 s；(c) t = 10 s；(d) t = 20 s；(e) t = 30 s；(f) t = 40 s。下载：下载高分辨率图片（381KB）下载：下载全尺寸图片

图9. Re=1000时盖子驱动的方形腔体流动的MPM模型结果与其他研究结果的比较。(a) y = 0.5L处的归一化垂直速度；(b) x = 0.5L处的归一化水平速度。

4.2. 情况2：浸没颗粒柱堆塌 [92]
虽然逻辑上的进展应该包括验证固相MPM求解器，但由于在之前的研究[76]中已经详细介绍了该单独的土相求解器的验证过程，因此这里不再重复。因此，本小节重点验证所提出的MPM模型在模拟浸没颗粒柱堆塌方面的有效性。

4.2.1. 模型描述
浸没颗粒柱堆塌是验证土-流体相互作用问题的数值模型的经典基准，并且已经通过实验和数值模拟得到了广泛的研究（例如[93,94,92]）。在这个测试中，最初用可移动的挡板将一个饱和沙子柱限制住。突然移除这个挡板会触发一次海底滑坡，从而可以观察到堆塌过程。我们采用了Wang等人[92]提出的平板上浸没颗粒柱堆塌实验作为基准，该实验已经使用SPH[94]、MPM[26]和通用粒子动力学[95]等方法进行了广泛的模拟。几何形状如图10所示。该设置包括一个6厘米×8厘米的柱子，其中填充了玻璃珠，浸没在一个长50厘米、宽10厘米的水箱中。

图10. 平板上浸没颗粒柱堆塌的示意图[92]。

图11展示了用于模拟浸没颗粒柱堆塌实验[92]的MPM模型。其几何形状复制了原始实验的尺寸。网格尺寸的选择遵循了Feng等人[26]之前的MPM研究。该域使用均匀的背景网格离散化，矩形元素的网格尺寸为dx=0.005米，固体和流体相都在这个共同网格上进行了求解。每个元素初始化了36个物质点，从而形成了6887个固体点和72,600个流体点的系统。模拟时间为1秒，足够颗粒物质沉淀。所有墙壁（即底部、左侧和右侧）都施加了无滑移边界条件。时间步长设置为5 × 10??秒，这是低于临界时间步长的（见附录B），以确保数值稳定性。

图11. Wang等人[92]提出的浸没颗粒堆塌的MPM设置。

材料参数定义如下：对于流体相（水），密度为ρw= 1000 kg/m3，动态粘度为ηw=1 × 10?3 Pa·s。对于固相（玻璃珠），我们采用Wang等人[92]提供的属性：颗粒密度ρs= 2500 kg/m3，固体体积分数ns= 0.55（孔隙率e ≈ 0.82），颗粒直径d = 2.5 × 10?3 m。土壤的参数包括：参考惯性数I0=0.279，静止摩擦系数μs=0.975，动态摩擦系数μd=1.63，以及正则化参数λ=1 × 10??。

4.2.2. 模型结果
图12显示了海底滑坡过程中土相速度大小的发展。在t = 0.1秒时，不稳定性开始出现，尤其是在柱子的外缘。到t = 0.2秒时，靠近自由表面的颗粒开始向下塌陷，而柱子的核心部分基本保持静止。在t = 0.2秒到0.8秒之间，颗粒物质向前移动。最终，在t = 1.0秒时，由于底部摩擦作用，所有运动停止。这种初始失败、塌陷、传播和最终停止的现象与其他使用[26,94]进行的数值研究中的定性行为非常吻合。这种一致性验证了所提出的交错网格显式-隐式MPM模型在模拟海底滑坡方面的有效性。

图12. 浸没颗粒堆塌过程中土-固相速度大小的发展。(a) t = 0.1 s；(b) t = 0.2 s；(c) t = 0.8 s；(d) t = 1.0 s。

图13展示了颗粒堆塌过程中流体相自由表面和压力场的演变。在t = 0.1秒时，由于固-流体相互作用，靠近土壤柱的自由表面开始波动。这种扰动在t = 0.2秒时被放大。随后，从t = 0.2秒到1.0秒，产生的波浪向正x方向传播，与滑坡移动的方向一致。此外，压力场保持静水压力稳定，从表面的0逐渐过渡到底部的ρwgh。波浪动力学和静水压力的物理真实表示证实了所开发的MPM模型正确捕捉了流体行为，进一步证明了所提出的MPM模型的有效性。

为了提供定量验证，我们将我们MPM模拟的最终沉积剖面与实验数据[92]、SPH结果[92]和显式MPM结果[26]进行了比较。如图14所示，我们模型预测的剖面与实验基准和数值参考结果均吻合良好。这种紧密匹配从定量上证明了建立的MPM框架能够准确模拟水-土耦合问题。

然而，所有数值模型（SPH、显式MPM和交错网格显式-隐式MPM）都显示与实验结果存在不匹配。这种差异主要是由于所使用的本构模型不同。具体来说，本研究采用了类流体的速率依赖模型，而其他方法使用的是类固体的速率无关模型。这两种方法都无法完全捕捉颗粒流中的复杂固-流体相变[25]，从而导致了一些与实验结果的偏差。尽管存在这种不匹配，但开发的显式-隐式MPM仍然提供了相对准确的预测。为了解决这种差异，未来的工作将致力于开发和实现一个非牛顿流变学模型，以考虑在所提出的两相两点MPM框架内的残余临界状态应力（见小节5.2）。

4.3. 情况3：倾斜平面上的海底滑坡 [96]
本小节验证了所提出的交错网格显式-隐式MPM框架，包括其在模拟海底滑坡和生成波浪方面的性能。

4.3.1. 模型描述
本小节采用了Rzadkiewicz等人[96]的实验进行验证，该实验研究了倾斜平面上的海底滑坡。如图15所示，实验设置包括一个长4.0米、宽1.5米的水箱。一个倾斜角度为45°、高度为1.5米的基座位于距离左侧水箱壁1.5米的位置；其上下长度分别为1.0米和2.5米。在一个倾斜面上放置了一个0.65米×0.65米的三角形沙质块，距离右侧水箱壁1.0米。滑坡和随后的波浪生成是由突然移除挡板引起的。

图15. 倾斜平面上海底滑坡的模型设置示意图[96]。

MPM模型复制了实验的几何形状。网格尺寸的选择遵循了我们之前的MPM研究[25]。离散化使用均匀网格（dx= 0.008米），每个单元有四个物质点，共生成了13,117个固体点和236,407个流体物质点。为了稳定性，采用了5 × 10??秒的时间步长（与第4.2节一致）。所有限制壁（倾斜壁、底部、左侧和右侧）都施加了无滑移边界条件。模拟时间为1秒。

材料参数定义如下：流体相（水）的密度为ρw=1000 kg/m3，动态粘度为ηw=1 × 10?3 Pa·s。对于固相（玻璃珠），我们采用了Zhu等人[34]和Feng等人[25]提供的属性：颗粒密度ρs=2650 kg/m3，初始固体体积分数ns=0.62（孔隙率nw=0.38），平均颗粒直径d = 3 × 10?3 m。颗粒材料的参数包括：μs=0.775（对应于Zhu等人[34]中的内摩擦角37°），μd= 0.997，I0= 0.279，以及λ= 1 × 10??。

图16展示了土相速度大小的变化。在t = 0.1秒时，未观察到显著的运动，但颗粒物质显示出沿倾斜表面移动的趋势。随后，在重力作用下开始塌陷，最高的速度出现在流动前沿。在t = 0.4秒到0.8秒之间，整个质量开始移动。最终，在t = 1.0秒时，由于底部摩擦作用，所有运动停止。这种观察到的初始失败、塌陷、传播和最终停止的现象与使用[26,94]进行的其他数值研究中的定性行为非常吻合。这种一致性验证了所提出的交错网格显式-隐式MPM模型在模拟海底滑坡方面的有效性。

图17 进一步展示了海底滑坡过程中土-固相速度大小的变化。(a) t = 0.1 s；(b) t = 0.3 s；(c) t = 0.4 s；(d) t = 0.8 s。

图17 进一步展示了海底滑坡过程中的流体压力。总体而言，这种情况呈现了一个稳定且线性的压力场，与小节4.2中的压力场相似。波浪的生成和传播直接对应于图16中显示的颗粒质量的运动。最初，在t = 0.1秒时，没有明显的波浪。当滑坡在t = 0.3秒和0.4秒时加速向下移动时，水面上形成了一个凹陷。到t = 0.8秒时，表面形成了一个明显的三段式轮廓，特征是一个中央凸起平台两侧是两个凹区域。这种流体自由表面的发展与Feng等人[25]使用显式MPM获得的数值结果一致，进一步验证了当前方法捕捉海底滑坡动态的能力。

图18 比较了当前MPM模型预测的颗粒质量自由表面轮廓与实验结果[96]以及其他数值方法（PFEM [35]、SPH [34]和显式MPM [25]）的解。在t = 0.4秒时，当前结果与基准数据总体吻合。值得注意的是，所提出的MPM模型预测的滑出距离比其他模拟更接近实验测量值。到t = 0.8秒时，当前模型仍提供了最准确的滑出距离预测。然而，预测的沉积高度低于实验和其他数值模型中观察到的高度，这表明虽然滑出距离被很好地捕捉到了，但最终的沉积高度被低估了。下载：下载高分辨率图像（503KB）下载：下载全尺寸图像图18. 从当前MPM结果、实验数据以及已发表的SPH、PFEM和显式MPM解决方案得到的滑坡自由表面。(a) t = 0.4秒；(b) t = 0.8秒。这些观察结果可以这样解释：首先，本研究建立的MPM采用了一个不可压缩的正则化μ(I)模型，该模型忽略了由于膨胀和收缩导致的体积变化。因此，所提出的MPM模型很好地预测了滑出距离，但低估了沉积厚度。其次，作为一个类流体本构模型，正则化μ(I)流变学本质上忽略了传统弹塑性模型[25,35,34]所捕捉到的类固体速率独立行为。因此，所提出的MPM模型产生的自由表面轮廓与使用传统弹塑性公式得到的轮廓有显著差异。为了解决这一限制，应该考虑一个能够考虑残余临界状态应力的非牛顿流变学模型（见第5节的讨论）。图19进一步比较了当前显式-隐式MPM在不同时间间隔的水体自由表面与实验数据和其他数值方法的结果。在t = 0.4秒时，当前模型预测的自由表面轮廓与实验数据相似，比PFEM [35]和SPH [34]的结果更为接近。到了t = 0.8秒时，MPM的结果与测量数据以及Feng等人的显式MPM模拟结果相符，并且比SPH和PFEM的预测更准确。这些结果进一步证明了所提出的显式-隐式MPM在捕捉海底滑坡波浪生成方面的能力。下载：下载高分辨率图像（348KB）下载：下载全尺寸图像图19. 滑坡产生的波浪比较：当前MPM结果与实验和其他数值模拟结果。(a) t = 0.4秒；(b) t = 0.8秒。4.4. 情况4：斜面上的海底滑坡[97]本小节进一步验证了所提出的MPM框架在模拟海底滑坡及其产生的波浪方面的性能。与小节4.3不同，这里省略了对特定状态变量（例如颗粒速度和水压）的详细分析，以避免重复。4.4.1. 模型描述如图20所示，第四个验证案例基于Grilli等人的实验[97]。原始实验[97]的设置包括一个6270毫米×330毫米的水池和一个位于35°斜面上的三角形颗粒柱（高42毫米，长84毫米）。波浪测量仪用于记录波高的时间变化。与Feng等人[25]一致，本研究采用了缩小的水池几何尺寸，长度为5340毫米，宽度为43毫米，以降低计算成本。实验测量了位于滑坡趾部下游600毫米和1600毫米处的波浪高度，这与Feng等人[25]采用的位置相匹配。下载：下载高分辨率图像（153KB）下载：下载全尺寸图像图20. 斜面上的海底滑坡示意图[97]。图20显示了MPM模型中的网格尺寸。网格尺寸的选择遵循了我们之前的MPM研究[25]。域使用均匀的MPM网格进行离散化，单元格大小为dx=0.005米。初始配置每个流体单元格分配16个材料点，每个固体单元格分配16个材料点，总共269,704个水材料点和3225个固体材料点。考虑了2秒的总模拟时间。使用5 × 10??秒的时间步长以确保数值稳定性，这与前面的小节一致。所有墙壁都应用了无滑移边界条件。材料参数定义如下：流体相（水）的特性是密度ρw= 1000 kg/m3和动态粘度ηw= 1 × 10?3 Pa·s。对于固体相，玻璃珠的参数如下：颗粒密度ρs= 2970 kg/m3，初始固体体积分数ns=0.634（孔隙率nw=0.366）；平均颗粒直径d = 3 × 10?3 m；μs=0.70；μd= 0.80，I0= 0.279，λ=0.0001。4.4.2. 模型结果图21定性地展示了由开发的交错网格显式-隐式MPM模拟的土壤剖面的演变。颗粒质量最初从t = 0.0秒到0.02秒在重力作用下向下移动。随后，从t = 0.02秒到0.47秒，它沿着斜面流动，形成了一条明显的带状轮廓。到t = 0.62秒时，颗粒物质到达水平面。该图显示了当前模拟、实验数据[97]以及显式MPM模型[25]在所有观察时间间隔内的良好一致性，验证了该模型在模拟海底滑坡动力学方面的有效性。下载：下载高分辨率图像（584KB）下载：下载全尺寸图像图21. 从当前MPM结果、实验数据和显式MPM解决方案得到的滑坡自由表面。(a) t = 0.0秒；(b) t = 0.02秒；(c) t = 0.17秒；(d) t = 0.32秒；(e) t = 0.47秒；(f) t = 0.62秒。图22进一步提供了当前MPM结果、实验数据[97]、FVM结果[80]和显式MPM结果[25]之间的自由表面高度的定量比较。这些数据集之间的一致性再次验证了所提出模型在捕捉波浪生成方面的有效性。下载：下载高分辨率图像（429KB）下载：下载全尺寸图像图22. 本研究中MPM结果与实验结果[97]、FVM结果[80]、显式MPM结果[25]之间的自由表面高度比较。(a) 传感器1；(b) 传感器2。5.1. 比较分析在验证了所提出的交错网格显式-隐式MPM模型之后，我们现在进一步分析和讨论这个MPM模型。首先，我们将它的模拟结果与使用两种替代的两相两点MPM公式（显式MPM模型和混合网格MPM模型）得到的结果进行比较。然后我们讨论了该模型的局限性和未来的工作。5.1. 比较分析我们首先使用小节4.2中的海底滑坡案例，将建立的显式-隐式MPM框架与显式MPM进行比较。图23展示了所提出的MPM和Feng等人[26]的显式MPM的模拟结果，揭示了两个关键差异：(i)滑坡剖面的差异：两种MPM模型都成功捕捉到了海底滑坡的整个过程。然而，所提出的交错网格显式-隐式MPM预测的土壤自由表面轮廓与显式MPM模型的略有不同。这种差异源于它们根本不同的本构框架：所提出的模型采用了一个依赖于速率的非牛顿流变学模型，而共位网格显式模型使用了一个与速率无关的Drucker-Prager模型。使用不同的材料参数进一步导致了预测滑坡剖面的差异。Feng等人[25]提供了关于不同本构模型差异的详细讨论。(ii)水压的差异：两种MPM模型都描述了滑坡产生的波浪的演变。图23清楚地表明，交错网格减少了虚假的压力模式，并产生了一个稳定的压力场，而共位网格导致了一个不稳定的压力场。这种改进主要归因于压力解算法：显式模型采用了弱可压缩公式，而当前方法利用了专为不可压缩流动设计的显式-隐式方案。下载：下载高分辨率图像（249KB）下载：下载全尺寸图像图23. 交错网格显式-隐式MPM与显式MPM之间的水压比较（黑线代表滑坡剖面）。(a) 交错网格显式-隐式MPM；(b) 显式MPM [26]。然后，将所提出的模型与小节4.3中的海底滑坡案例相结合，使用混合网格MPM进行比较。图24展示了所提出的MPM和He等人[27]的混合网格MPM的模拟结果。与图23类似，一个关键差异在于预测的海底滑坡滑出路径。不同的本构模型导致了不同的运动模式。关于水压，虽然He等人[27]没有报告，但两种算法固有的不可压缩公式表明它们都能够产生稳定的压力场。下载：下载高分辨率图像（238KB）下载：下载全尺寸图像图24. 交错网格显式-隐式MPM与混合网格MPM之间的土壤-固体相剖面比较。(a) 交错网格显式-隐式MPM；(b) 混合网格MPM [27]。表3总结了所提出的交错网格显式-隐式MPM、显式MPM和混合网格MPM模型之间的主要区别。除了它们在土壤和流体相方面的固有算法差异外，还有五个主要对比点：(i)计算效率：交错网格显式-隐式MPM的关键时间步长主要受渗透性控制，而混合网格MPM还必须满足土壤相的Courant-Friedrichs-Lewy (CFL)条件。显式MPM进一步受到流体相CFL条件的限制。因此，交错网格显式-隐式MPM可以采用最大的允许时间步长，因此在三种公式中具有最低的计算成本。这些约束的方程在附录B中有详细说明。(ii)滑坡运动：由于用于土壤-固体相的不同算法，交错网格显式-隐式MPM通常采用依赖于速率的本构模型，而显式MPM和混合网格MPM可以采用与速率无关的本构模型。(iii)水压稳定性：交错网格显式-隐式和混合网格MPM公式在模拟稳定的不可压缩流体流动方面表现出优势。然而，显式MPM倾向于产生不稳定和振荡的水压场。(iv)本构模型复杂性：交错网格显式-隐式MPM需要最简单的土壤本构公式，因为它避免了其他两种框架中所需的复杂应力积分算法。(v)网格系统复杂性：交错网格显式-隐式和显式MPM使用简化的网格结构，对固体和流体相都使用相同的网格。混合网格MPM使用独立的网格，具有最高的网格复杂性。(vi)数学模型：交错网格显式-隐式MPM采用不可压缩的两相流公式，这与Biot的混合理论不同。这种区别在附录A中有详细说明，可能导致对海底滑坡的不同预测。表3. 不同类型的两相两点MPM框架之间的比较。方面项目交错网格显式-隐式MPM显式MPM混合网格MPM固体相时间积分显式-隐式显式显式本构模型非牛顿模型类固体模型（例如，MC，DP）类固体模型（例如，MC，DP）网格类型交错网格共位网格共位网格流体相时间积分显式-隐式显式显式-隐式本构模型牛顿模型状态方程牛顿模型网格类型交错网格共位网格交错网格关键时间步长Δtcri = ΔtkΔtcri = min(Δtk,ΔtCFL,s, ΔtCFL,w)Δtcri = min(Δtk,ΔtCFL,s)关键时间步长交错网格MPM>混合网格MPM>显式MPM滑坡运动交错网格MPM≠（显式MPM，混合网格MPM）水压稳定性（交错网格MPM≈ 混合网格MPM）> 显式MPM土壤本构复杂性交错网格MPM<（显式MPM，混合网格MPM）MPM网格复杂性（交错网格MPM≈ 显式MPM）< 混合网格MPM备注1此比较描述了每种MPM方法的基本差异、独特优势和固有弱点。这些见解可以用于选择这些模型来解决海底滑坡问题。5.2. 局限性和未来工作尽管第5.1节中概述了优势，但下面还讨论了几项限制和未来的工作。(i)恒定体积假设。当前工作对土壤采用了恒定体积假设。虽然这个假设在临界状态下的大变形时是有效的，但它未能捕捉到海底滑坡开始时的体积行为（即小变形范围），此时松散的沙子会收缩，而密集的沙子会膨胀。因此，这个假设可能会引入错误，例如低估沉积高度。为了解决这一限制，应该开发并在所提出的两相两点MPM框架内实现一个能够考虑残余临界状态应力的非牛顿流变学模型（参见[80]中引用的模型）。(ii)三维大规模真实世界海底滑坡模拟。本工作采用了二维框架。未来的工作应该将模型扩展到三维，并针对真实世界中的海底滑坡案例进行验证。(iii)表面侵蚀模拟。尽管本工作使用了所提出的两相两点MPM模型来模拟海底滑坡，但该模型在模拟表面侵蚀（即土壤表面的颗粒运动）方面具有巨大潜力，其中体积变形可以忽略不计，剪切流变学占主导地位（参见[98]）。为了进一步验证这一能力，未来的工作应该将模型扩展到其他条件，如溢出、局部冲刷和海岸侵蚀。6. 结论在这项研究中，我们提出了一种新颖的交错网格显式-隐式两相两点MPM框架，用于模拟海底滑坡。该MPM框架基于一个两相流数学模型，其应力框架不同于Biot的混合理论。使用两组材料点来离散化水和土壤-固体相，并使用共享的交错网格来解决它们的控制方程。它们之间的相互作用是根据来自相同网格节点的物理信息计算得出的，从而实现了土壤与水之间的平滑且有效的耦合。该方法采用了显式-隐式MPM框架内的算子分割策略，并分为两个阶段进行。在第一阶段，预测步骤明确计算土壤-固体相和流体相的中间速度；在第二阶段，校正步骤隐式求解压力泊松方程，并利用其解来更新这些中间速度。通过模拟一系列实验进一步评估了该模型的性能。该模型对海底滑坡及其产生的波浪的预测结果与现有实验数据以及其他研究的结果高度吻合，验证了所开发的交错网格显式-隐式MPM模型的有效性。为了突出其独特优势，将这一经过充分验证的MPM模型与两种现有的双相两点MPM公式（完全显式模型和混合网格MPM模型）进行了比较：首先，通过将流体视为不可压缩介质来生成稳定的流体场；其次，它使用了一个简单的土壤本构模型，该模型在MPM中易于实现，避免了弹性塑性模型所需的复杂应力积分算法；第三，它简化了双相两点MPM的实施过程，因为土壤相和流体相共享相同的动量求解方法、相同的交错网格和相同的时间步长；最后，与其它公式相比，它降低了计算成本，因为它允许更大的时间步长。尽管本研究在MPM框架内采用了不可压缩的双相流模型，但该框架也可以应用于其他数值方法，如MPS和SPH。

数据可用性声明：本研究期间生成和/或分析的数据集可应请求向相应作者获取。

作者贡献声明：
冯航：撰写-审稿与编辑、撰写-初稿、可视化、验证、方法论、概念化。
尹振宇：撰写-审稿与编辑、监督、资源提供、数据调查。