通过ISPH方法开发用于多相流模拟的介观尺度颗粒模型

《Computational Particle Mechanics》：Development of Interphase Particle for multiphase flow simulation by ISPH

　　鹤田直树|后藤隆文|渡边茂人|阿巴斯·海耶尔|后藤仁
日本横须贺港湾机场研究所海岸水力工程部

**摘要**
在传统的不压缩平滑粒子流体动力学（ISPH）中，通过保持每个粒子的粒子数密度恒定来确保体积守恒。然而，在多相流模拟中，不同密度的存在会破坏基于粒子数密度的连续性方程的兼容性。这通常会导致界面处的数值不稳定性和不准确性。为了解决这个问题，本研究提出了一种新的、稳定的多相方案，引入了一种称为“界面粒子”（IP）的数学构造。IP作为一种中介，负责不同相之间的粒子间作用力。与传统的基于插值的方法或网格外方法不同，IP直接集成到压力泊松方程矩阵中，使得所有相的压力可以在单一离散化空间内同时求解。通过几个基准测试验证了所提方案在能量守恒、界面再现性（包括壁边界）和数值稳定性方面的良好性能：两相静止流体柱、振荡液滴、上升气泡的合并以及溃坝现象。

1. **引言**
气液两相流在各种工程现象中起着关键作用。特别是在海岸和海洋领域，波浪破碎过程中，细小气泡的卷入显著影响波高波动和能量耗散[1]。Hossain等人（2024年）[2]量化了破碎波中气泡卷入对能量耗散的影响。Watanabe等人（2005年）[3]进行了LES模拟，研究了大尺度涡旋将气泡吸入海洋的过程，并证明了气相卷入决定了水中的涡旋形成和动量传输机制。此外，被包裹的气囊表现出“缓冲效应”，这通过Bagnold的压力-脉冲理论来描述沿海结构的峰值冲击压力[4]。Bredmose等人（2001年）[5]也研究了这种空气效应的建模。这种气液相互作用在开发高效内燃机时也是一个基本问题。在这种情况下，燃料雾化控制混合气的形成和随后的燃烧效率；因此，精确的喷雾动力学数值预测[6,7]对于优化性能至关重要。除了开放表面和喷雾流动外，管道中的气液流动在管道[8]和污水系统[9]等众多工程领域中也至关重要。

如前所述，气液界面在工程中具有相当大的重要性，需要精确跟踪以理解相关现象。然而，由于密度差异显著，每个相的流体在界面附近的行为有很大不同，形成了复杂的流场，这使得稳定模拟变得困难。图1展示了这一问题的一个典型例子。空气和水的密度比约为1000：1000，因此加速度比也约为1000，导致两者之间的运动有显著差异。处理界面的方法对于高精度和稳定地模拟这些复杂界面至关重要。已经开发了多种计算流体动力学模型来模拟这些复杂的气液界面。为了处理界面，广泛使用了水平集[10,11]和体积积分（VOF）方法[12,13]，因为它们能够处理重大的拓扑变化。VOF方法结合了自适应网格细化技术[14]或水平集方法[15]。浸没边界（IB）方法[16]也被用来有效地模拟流体与复杂固体几何形状之间的相互作用。然而，这些基于网格的方法往往难以保持质量守恒或捕捉高度碎裂的界面[17]。因此，基于粒子的方法（如SPH[18]、[19]、[20]）和MPS[21]方法作为模拟涉及极端变形的剧烈多相流的强大替代方案出现[22]。

**下载：**下载高分辨率图像（169KB）
**下载：**下载全尺寸图像

图1. 界面上方的加速度差异。

在SPH社区中，提出了多种基于粒子的多相方法来克服数值不稳定性。例如，Khayyer等人（2016年）[23]提出了一种一阶精确的密度平滑方案，用于MPS方法，以处理尖锐的密度不连续性并提高数值稳定性。该方案的性能优于SPH和基于MPS的多相流模拟中常用的零阶方法。处理目标相之间显著不同的密度界面至关重要，修改SPH的密度模型可以获得良好的性能[24]。为了提高稳定性，引入人工排斥力是有效的，从而使粒子均匀分布。例如，Sun等人（2019年）[24]和Ming等人（2017年）[25]对界面施加稳定力，并展示了它们在多相流再现性方面的良好性能。Rezavand等人（2020年）[27]将Riemann求解器扩展到两相流中，有助于稳定界面相互作用。他们还在轻相中使用了传输-速度公式[28]来消除非物理空洞或碎片化。在粒子分布方面，粒子移动[29,30]技术也可以成为多相流模拟的强大工具[31,32]。类似网格外方法的插值方法也对再现界面有效。Shimizu等人（2018年）[33]开发了一种使用自由表面边界模型的插值技术，即空间势粒子（PPP）[34]。PPP通过虚拟粒子建立了两相之间的连接，并给出了它们之间的相互作用。所提出的方案能够在没有压力噪声的情况下实现稳定模拟。

如前所述，已经提出了多种方法来提高多相粒子模拟的数值精度和稳定性。通过解决具有大幅不同密度的相之间的压力来实现稳定模拟的挑战导致了额外程序的引入，包括人工模拟的力或移动技术。这些过程增强了粒子分布的规则性并达到了稳定状态。问题在于，大多数这些系统需要在紧凑的支持范围内具有平滑的界面。这对于稳定地施加强度或以明确区分所属相区的方式移动粒子是必要的。与上述旨在改善粒子分布的方法相反，另一种概念允许不同相之间存在不规则性，包括粒子的重叠。这种方法涉及为每个相应用不同的离散化空间，并使用网格外等方法进行插值，以建立每个离散化空间之间的连接。然而，鉴于准确再现界面力所需的可靠插值的必要性，其在剧烈界面流动情况下的应用受到限制。这种限制是由于在界面周围不可避免地会出现不稳定的截断核支持，特别是在涉及飞溅粒子的场景中。此外，基于投影的粒子方法面临着基于粒子数密度的连续性方程的根本挑战，因此，用于处理不同密度的底层算法也存在问题。这将在后续章节中详细讨论。这是成功再现多相流体物理状态所必须解决的关键方面。

为了进一步提高基于投影的粒子方法的适用性，本研究解决了传统ISPH（不可压缩SPH）[35]中多相流的一个关键问题。这个问题涉及使用基于粒子数密度的压力泊松方程进行体积守恒。为了解决这个问题，我们提出了一种新的虚拟粒子，即界面粒子（IP），它将不同相之间的粒子间作用力作为其目标粒子的代理。IP保持与目标粒子的质量-体积-密度一致性。通过在两相静止流体柱、振荡液滴、上升气泡的合并以及溃坝等基准测试中验证了所提出的IP方案的数值稳定性和准确性。

2. **方法**
ISPH[35]是一种基于以下基本方程的基于投影的粒子方法：
(1) \(D\rho\dot{t}+\rho\nabla\cdot u = 0\)
(2) \(\dot{u}\dot{t}=-\nabla p\rho+\nu\nabla^{2}u+g\)
其中\rho是密度，\(t\)是时间，\(u\)是速度向量，\(p\)是压力，\(\nu\)是动态粘度系数，\(g\)是重力加速度向量。利用连续性方程，得到压力泊松方程：
(3) \(\frac{1}{\rho}D\rho\dot{t}=-\nabla\cdot u^{*}=-\nabla^{2}p\rho\Delta t\)
其中下标\(*\)表示在计算压力泊松方程之前的投影方法的第一步后的临时状态。这些方程在ISPH模拟中通过对每个计算粒子进行局部加权平均来计算。目标粒子\(i\)的局部平均密度\(\langle\rho\rangle\)计算如下：
(4) \(\langle\rho\rangle_{i}=\sum_{j} w_{ij}v_{ij}=\sum_{j}\rho_{j}V_{j}w_{ij};\quad w_{ij}=w(q_{ij}),\quad q_{ij}=|\r_{j}-r_{i}|/h\)
其中\(m\)是粒子质量，\(V\)是粒子体积，\(w\)是加权平均操作的核函数，\(r\)是坐标向量，\(h\)是平滑长度，下标\(i\)和\(j\)分别表示目标粒子\(i\)及其邻近粒子\(j\)。在本研究中，均匀应用了Wendland核[36]。
(5) \(w_{q}=\alpha Ds\{1-\frac{q^{2}}{1+\frac{q}{2}}\quad(0\leq q\leq2),\quad 2其中\(\alpha ds\)是一个归一化参数。下面定义的粒子数密度\(n\)对于构建方程(4)至关重要：
(6) \(n_{i}=\sum_{j}v_{j}w_{ij}\)

在单相流的情况下，如果每个粒子的设定密度\(\rho\)相等，则需要始终保持每个粒子的局部平均密度恒定以确保体积守恒。在方程(4)中，通过控制方程(6)中的粒子数密度\(n\)保持恒定可以实现这一点。在传统的isph单相流算法中，通过以下方式使用粒子数密度\(n\)来实现体积守恒：
(7) \(\frac{1}{\rho}d\rho\dot{t}=\frac{1}{n}d_{n}\dot{t}\)
使用这种替换后，方程(3)重写为：
(8) \(\frac{1}{\rho}d\rho\dot{t}=-\nabla^{2}p\rho\delta t\to\frac{1}{n}d_{n}\dot{t}=-\nabla^{2}p\rho\delta t\)

多相流模拟的问题在于粒子密度\(\rho\)是变化的。因此，不能简单地将密度\(\rho\)替换为粒子数密度\(n\)，也就不能用方程(7)和(8)。图2展示了这样一个典型例子。尽管粒子数密度保持恒定，从而实现了体积守恒，但局部平均密度仍然不同。

**下载：**下载高分辨率图像（258kb）
**下载：**下载全尺寸图像

图2. isph中密度与粒子数密度之间的关系。

为了适应多相流的粒子分布，方程(1)需要在右侧有一个未知值来考虑变化的密度。这意味着传统的基于投影的isph算法不能应用于多相流。这在实施多相流模拟时是一个关键问题，因为如图1所示，由于密度差异显著，会导致数值不稳定。

为了解决这个根本问题，引入了伪粒子来表示不同相之间的过渡界面。此后，这些伪粒子被称为“界面粒子”（ip）。需要注意的是，ip是一个表示粒子边界或连接的数学构造，它不对应任何物理实体，也不存在物理上。如图3所示，ip被设置为与目标粒子\(i\)具有相同的位置、速度、质量和密度。考虑到粒子间作用力，ip不仅与目标粒子\(i\)连接，还与具有不同密度的其他ip粒子\(k\)连接。ip将所有邻近粒子\(k\)视为目标粒子\(i\)的单一粒子。而目标粒子\(i\)则通过其自身的ip考虑来自不同相的粒子\(k\)的粒子间作用力。

**下载：**下载高分辨率图像（267kb）
**下载：**下载全尺寸图像

图3. 界面粒子的图形表示。

目标粒子\(i\)的ip核权重\(w\)由不同相中邻近粒子\(k\)的总权重给出：
(9) \(w_{iip}=sum_{k \in k}v_{k}ip_{wiip_{k}ip_{k}\)
其中\(k,k = \{i\neq k\}\)
由于每个ip的位置和体积与其目标粒子相同，ip之间的核权重如下：
(10) \(v_{kip}_{wiip_{k}}=\frac{v_{k}w_{ik}}{v_{k}}\)
利用这个关系和方程(4)，可以得到局部平均密度\(\rho\)：
(11) \(\langle\rho\rangle_{i}=\frac{\rho_{i}ip_{wiip_{i}}+\sum_{l \in l}v_{l}w_{il}}_{\text{相同相}}+\sum_{k \in k}v_{k}w_{ik}ip_{k}}_{\text{其他相}}+\frac{\rho_{i}sum_{j \in v_{j}w_{ij}}{l=\{r_{i}=\rho_{i}\text{且}l\neq ip},j=k\cup l}\)
图4说明了方程(11)中的关系。该方程使得仅使用相同的密度\(\rho_{i}\就可以计算目标粒子\(i\)的局部平均密度。物理上存在的粒子\(i\)和\(k\)可以通过类似于单相流模拟中使用方程(7)的方式保持粒子数密度恒定，从而确保体积守恒，而不考虑邻近粒子的不同密度。 其中\(\alpha ds\)是一个归一化参数。下面定义的粒子数密度\(n\)对于构建方程(4)至关重要： (6) \(n_{i}=\sum_{j}V_{j}w_{ij}\) 在单相流的情况下，如果每个粒子的设定密度\(\rho\)相等，则需要始终保持每个粒子的局部平均密度恒定以确保体积守恒。在方程(4)中，通过控制方程(6)中的粒子数密度\(n\)保持恒定可以实现这一点。在传统的isph单相流算法中，通过以下方式使用粒子数密度\(n\)来实现体积守恒： (7) \(\frac{1}{\rho}d\rho\dot{t}=\frac{1}{n}D_{n}\dot{t}\) 使用这种替换后，方程(3)重写为： (8) \(\frac{1}{\rho}d\rho\dot{t}=-\nabla^{2}p\rho\Delta t\to\frac{1}{n}d_{n}\dot{t}=-\nabla^{2}p\rho\Delta t\) 多相流模拟的问题在于粒子密度\(\rho\)是变化的。因此，不能简单地将密度\(\rho\)替换为粒子数密度\(n\)，也就不能用方程(7)和(8)。图2展示了这样一个典型例子。尽管粒子数密度保持恒定，从而实现了体积守恒，但局部平均密度仍然不同。 **下载：**下载高分辨率图像（258kb） **下载：**下载全尺寸图像 图2. isph中密度与粒子数密度之间的关系。 为了适应多相流的粒子分布，方程(1)需要在右侧有一个未知值来考虑变化的密度。这意味着传统的基于投影的isph算法不能应用于多相流。这在实施多相流模拟时是一个关键问题，因为如图1所示，由于密度差异显著，会导致数值不稳定。 为了解决这个根本问题，引入了伪粒子来表示不同相之间的过渡界面。此后，这些伪粒子被称为“界面粒子”（ip）。需要注意的是，ip是一个表示粒子边界或连接的数学构造，它不对应任何物理实体，也不存在物理上。如图3所示，ip被设置为与目标粒子\(i\)具有相同的位置、速度、质量和密度。考虑到粒子间作用力，ip不仅与目标粒子\(i\)连接，还与具有不同密度的其他ip粒子\(k\)连接。ip将所有邻近粒子\(k\)视为目标粒子\(i\)的单一粒子。而目标粒子\(i\)则通过其自身的ip考虑来自不同相的粒子\(k\)的粒子间作用力。 **下载：**下载高分辨率图像（267kb） **下载：**下载全尺寸图像 图3. 界面粒子的图形表示。 目标粒子\(i\)的ip核权重\(w\)由不同相中邻近粒子\(k\)的总权重给出： (9) \(w_{iip}=sum_{k \in k}v_{k}ip_{wiip_{k}ip_{k}\) 其中\(k,k=\{i\neq k\}\) 由于每个ip的位置和体积与其目标粒子相同，ip之间的核权重如下： (10) \(v_{kip}_{wiip_{k}}=\frac{V_{k}w_{ik}}{V_{k}}\) 利用这个关系和方程(4)，可以得到局部平均密度\(\rho\)： (11) \(\langle\rho\rangle_{i}=\frac{\rho_{i}IP_{wiIP_{i}}+\sum_{l \in l}v_{l}w_{il}}_{\text{相同相}}+\sum_{k \in k}v_{k}w_{ik}ip_{k}}_{\text{其他相}}+\frac{\rho_{i}sum_{j \in v_{j}w_{ij}}{l=\{r_{i}=\rho_{i}\text{且}l\neq ip},j=K\cup l}\)>其中\(\alpha ds\)是一个归一化参数。下面定义的粒子数密度\(n\)对于构建方程(4)至关重要：
(6) \(n_{i}=\sum_{j}v_{j}w_{ij}\)

在单相流的情况下，如果每个粒子的设定密度\(\rho\)相等，则需要始终保持每个粒子的局部平均密度恒定以确保体积守恒。在方程(4)中，通过控制方程(6)中的粒子数密度\(n\)保持恒定可以实现这一点。在传统的isph单相流算法中，通过以下方式使用粒子数密度\(n\)来实现体积守恒：
(7) \(\frac{1}{\rho}d\rho\dot{t}=\frac{1}{n}d_{n}\dot{t}\)
使用这种替换后，方程(3)重写为：
(8) \(\frac{1}{\rho}d\rho\dot{t}=-\nabla^{2}p\rho\delta t\to\frac{1}{n}d_{n}\dot{t}=-\nabla^{2}p\rho\delta t\)

多相流模拟的问题在于粒子密度\(\rho\)是变化的。因此，不能简单地将密度\(\rho\)替换为粒子数密度\(n\)，也就不能用方程(7)和(8)。图2展示了这样一个典型例子。尽管粒子数密度保持恒定，从而实现了体积守恒，但局部平均密度仍然不同。

**下载：**下载高分辨率图像（258kb）
**下载：**下载全尺寸图像

图2. isph中密度与粒子数密度之间的关系。

为了适应多相流的粒子分布，方程(1)需要在右侧有一个未知值来考虑变化的密度。这意味着传统的基于投影的isph算法不能应用于多相流。这在实施多相流模拟时是一个关键问题，因为如图1所示，由于密度差异显著，会导致数值不稳定。

为了解决这个根本问题，引入了伪粒子来表示不同相之间的过渡界面。此后，这些伪粒子被称为“界面粒子”（ip）。需要注意的是，ip是一个表示粒子边界或连接的数学构造，它不对应任何物理实体，也不存在物理上。如图3所示，ip被设置为与目标粒子\(i\)具有相同的位置、速度、质量和密度。考虑到粒子间作用力，ip不仅与目标粒子\(i\)连接，还与具有不同密度的其他ip粒子\(k\)连接。ip将所有邻近粒子\(k\)视为目标粒子\(i\)的单一粒子。而目标粒子\(i\)则通过其自身的ip考虑来自不同相的粒子\(k\)的粒子间作用力。

**下载：**下载高分辨率图像（267kb）
**下载：**下载全尺寸图像

图3. 界面粒子的图形表示。

目标粒子\(i\)的ip核权重\(w\)由不同相中邻近粒子\(k\)的总权重给出：
(9) \(w_{iip}=sum_{k \in k}v_{k}ip_{wiip_{k}ip_{k}\)
其中\(k,k = \{i\neq k\}\)
由于每个ip的位置和体积与其目标粒子相同，ip之间的核权重如下：
(10) \(v_{kip}_{wiip_{k}}=\frac{v_{k}w_{ik}}{v_{k}}\)
利用这个关系和方程(4)，可以得到局部平均密度\(\rho\)：
(11) \(\langle\rho\rangle_{i}=\frac{\rho_{i}ip_{wiip_{i}}+\sum_{l \in l}v_{l}w_{il}}_{\text{相同相}}+\sum_{k \in k}v_{k}w_{ik}ip_{k}}_{\text{其他相}}+\frac{\rho_{i}sum_{j \in v_{j}w_{ij}}{l=\{r_{i}=\rho_{i}\text{且}l\neq ip},j=k\cup l}\)
图4说明了方程(11)中的关系。该方程使得仅使用相同的密度\(\rho_{i}\就可以计算目标粒子\(i\)的局部平均密度。物理上存在的粒子\(i\)和\(k\)可以通过类似于单相流模拟中使用方程(7)的方式保持粒子数密度恒定，从而确保体积守恒，而不考虑邻近粒子的不同密度。>另一方面，IP并不是一种物理上存在的粒子，并且在每个时间步长都会被重新分配；因此，它们的误差会得到更新，而不会影响物理上存在的粒子。下载：下载高分辨率图片（311KB）下载：下载全尺寸图片图4. 粒子与IP之间的密度关系。压力泊松方程（PPE）中的源项和压力的拉普拉斯算子是通过IP计算的。速度散度计算如下：
(12)??·u?i=uiiIP?=0ViIP·?wiiIP?ownIP+∑l∈LuilVl·?wil?samephase=∑l∈LuilVl·?wil?samephase
对于粒子i的IP，有如下方程：
(13)??·u?iIP=uiIPi?=0Vi·?wiIPi?particlei+∑k∈KuiIPkIPVkIP·?wiIPkIP?otherphase
IPs=∑k∈KuikVk·?wik?otherphase
将方程(12)和(13)相加，对应的就是不使用IP技术时粒子i处的速度散度。对于压力的拉普拉斯算子，粒子i的设置如下：
(14)??2pρ?i=?·piiIPρ￣iiIPViIP?wiiIP?ownIP+∑l∈L?·pilρ￣ilVl?wil?samephase=∑k∈K?·piiIPρiVk?wik?ownIP+∑l∈L?·pilρiVl?wil?samephase
其中ρil上的上划线表示粒子i和l的平均密度。注意，方程右侧的第一项是由单个IP表示的；因此，该值被分配给拉普拉斯压力系数矩阵中的单个分量。而IP的设置如下：
(15)??2pρ?iIP=?·(piIPiρ￣iIPiVi?wiIPi)?particlei+∑k∈K?·(piIPkIPρ￣iIPkIPVkIP?wiIPkIP)?otherphase
IPs=∑k∈K?·(piIPiρiVk?wik)?particlei+∑k∈K?·(piIPkIPρ￣ikVk?wik)?otherphase
与(14)类似，方程(15)右侧的第一项也被分配给矩阵的单一分量。使用方程(12)、(13)、(14)、(15)在PPE中，粒子i及其自身的IP之间的相互作用是分别求解的，然后相互连接。
从PPE中获得粒子压力后，通过将IP的相互作用反映回目标粒子i来计算压力梯度项：
(16)??pρ?i=piiIPρ￣iiIPViIP?wiiIP?ownIP+∑l∈Lpilρ￣ilVl?wil?samephase+piIPiρ￣iIPiVi?wiIPi?particlei+∑k∈KpiIPkIPρ￣iIPkIPVk?wik?otherphases
viaIP?viaownIP
这种算法意味着IP作为目标粒子与其他相中的物理存在粒子相互作用的代理。这些粒子通过它们在PPE中的自身IP在数学上与其他相的粒子相连。这个过程中的关键是IP是一个数学构造，代表了不同相粒子的互连，并且它也被作为一种额外的粒子纳入PPE中。因此，包括IP在内的所有粒子的压力都在相同的离散化空间内同时求解。这与基于插值技术的网格外方法完全不同，如图5所示。IP是一个伪粒子，代表了不同相之间的过渡界面。下载：下载高分辨率图片（238KB）下载：下载全尺寸图片图5. 网格外方法与IP的比较。
以下是IP的分布规则：
- 物理上存在的粒子通过它们在PPE中的自身IP在数学上相互连接，从而在不同相之间产生粒子间力。
- IP不仅连接到目标粒子i，还连接到具有不同密度的相邻物理存在粒子k的IP。
- IP也被作为额外的粒子纳入矩阵中，并在PPE中获得压力。
- IP设置在目标粒子i的位置，具有相同的速度、质量和密度。
- IP在每个时间步长都会被重新分配。
- 每个存在于界面附近的物理目标粒子都会生成一个IP。然后，通过IP获得的粒子间力被加到现有的目标粒子上，以计算压力梯度力，即它们被积分。
图6显示了IP应用的流程图。下载：下载高分辨率图片（457KB）下载：下载全尺寸图片图6. IP的流程图。
所提出的方法涉及生成替代粒子（IPs）以分解PPE中的相互作用，并重新整合它们以获得压力梯度项。这种基于IP的分解满足了ISPH连续性方程，该方程是通过用粒子数密度替换密度从均匀密度场推导出来的。此外，由于大的密度差异在PPE中产生的误差在每个时间步长通过重置IP而被锁定和消除；因此，这些误差不会影响物理上存在的粒子。
在本研究中，液体和空气相都被视为不可压缩流体。然而，通过考虑可压缩形式的PPE [37]，IP方案可以很容易地扩展到不可压缩-可压缩多相计算中。此外，通过考虑方程(10)、(11)、(12)、(13)、(14)、(15)、(16)、(9)，IP方案可以方便地扩展到多分辨率多相流动，因为在推导出的公式中考虑了每个单独粒子的体积。
3. 仿真结果和讨论
为了检验所提出的IP方案的性能，实施了针对两相静水流体柱、在中心力场下具有两种流体的振荡液滴 [38]、两个上升气泡的合并 [39] 以及气-液溃坝的仿真。为了确保实际适用性，在仿真中将Courant数αdt设置为0.1或0.05。CISPH-HS-HL-LBP-DS-ADS [[40]、[41]、[42]、[43]、[44]、[45]] 被用作基本设置以提高数值精度和稳定性。DS仅应用于壁粒子，而ADS应用于流体粒子之间的相互作用。
3.1. 两相静水流体柱
进行了两相静水流体柱的仿真，以检验所提出的IP方案的基本性能。两相流体柱各高1.0米，深2.0米，放置在一个宽2.0米的容器中，如图7所示。密度设置为下部柱ρd = 1000.0 kg/m3，上部柱ρl = 1.0 kg/m3。粘度系数分别设置为下部柱μd = 0.001 kg/(m·s)，上部柱μl = 0.000018 kg/(m·s)。计算粒子直径设置为d = 0.025米。在标准模型中，壁粒子具有与下部柱相同的属性。所有基准测试都使用了传统的ISPH壁边界模型。固定粒子作为壁粒子分布，并以类似流体粒子的方式纳入PPE矩阵中。IP方案应用于壁粒子和流体粒子。数值设置详情列在表1中。下载：下载高分辨率图片（176KB）下载：下载全尺寸图片图7. 静水两相流体柱仿真的初始条件。表1. 静水两相流体柱仿真的数值设置。
表1. 静水两相流体柱仿真的数值设置。
粒子间距 d [m] 0.025
平滑长度比 h/d 1.5
表面阈值 β [ni/n0] 0.95
密度比 ρd/ρl 1000
核心 Wendland核
时间步长规则 CFL条件（Courant数 αdt = 0.1）
PPE求解器类型 ICCG
PPE容差 1 × 10?8
PPE的平均迭代次数
无IP时：83.8，有IP时：111.7
增强方案 CISPH-HS-HL-LBP-DS-ADS
图8显示了仿真的粒子快照，以及速度矢量、密度和压力场（上图：无IP；下图：有IP）。由于两柱之间的压力差异显著，因此在每个仿真的结果中显示了不同颜色范围的压力场快照。从快照中可以看出，无IP时仿真中产生了非物理的循环流动。壁（ρd = 1000.0 kg/m3）和较高柱（ρl = 1.0 kg/m3）之间的密度差异导致了非物理的力，将较轻的柱从侧壁和较低柱推开。因此，上部柱的自由表面线上的粒子运动以及两柱之间的界面处的压力分布出现了轻微的振荡。相比之下，IP有效地抑制了非物理的循环流动和振荡。这种相位之间的稳定框架也对壁有效。
图8. 两相静水流体柱仿真中t = 20.0s时的速度矢量（左图：无IP；右图：有IP）和不同颜色范围的密度（左图：0 – 10 N/m2，右图：0 – 1000 N/m2）和压力场的快照。
图9显示了仿真结果中容器中心线（x = 1.0m）上的压力时间序列（上图：无IP；下图：有IP）。图10显示了容器中心位置（x = 1.0m）处每个位置的压力标准差。在没有IP的仿真结果中，较轻柱的每个测量位置都有明显的噪声。此外，由于两者之间的密度差异，两柱界面周围的压力的连续扰动在整个仿真中都很明显。另一方面，IP有效地抑制了所有线上的噪声，包括y = 1.0m处的界面周围的噪声。
图9. 两相静水流体柱仿真中的压力时间序列（上图：无IP；下图：有IP）。
图10. 两相静水流体柱仿真中的压力标准差。
3.2. 在中心力场下的振荡液滴
为了检验所提出的IP方案对于更动态的自由表面变形的有效性，实施了广泛模拟的基准测试——在中心力场下包含两种流体的2D振荡液滴 [38]。这个基准测试被广泛尝试 [46]、[47]、[48]、[49]]。外力fex设置为fex = –Ω2r；Ω = 1.2。为了检验收敛性，粒子大小d分别设置为d = 0.005、0.01和0.02米。Courant数设置为αdt = 0.1。初始条件如图11所示。密度设置为：ρl = 1.0 kg/m3，ρd = 1000.0 kg/m3。初始速度由以下公式给出：
(17)u=σx
(18)v=?σy
其中σ设置为σ = 0.4。其他数值设置详情列在表2中。下载：下载高分辨率图片（617KB）下载：下载全尺寸图片图11. 在中心力场下振荡液滴仿真的初始条件。表2. 在中心力场下振荡液滴仿真的数值设置。
表2. 在中心力场下振荡液滴仿真的数值设置。
粒子间距 d [m] 0.020、0.010、0.005（3种情况）
平滑长度比 h/d 1.5
表面阈值 β [ni/n0] 0.95
密度比 ρd/ρl 1000
核心 Wendland核
时间步长 CFL条件（Courant数 αdt = 0.1）
PPE求解器类型 ICCG
PPE容差 1 × 10?8
PPE的平均迭代次数
无IP时：43.4（d = 0.020m），70.1（d = 0.020m），88.2（d = 0.010m），131.1（d = 0.010m），155.1（d = 0.005m）
262.0（d = 0.005m）
增强方案 CISPH-HS-HL-LBP-ADS
图12显示了仿真结果中粒子的快照及其密度和压力分布。为了清晰显示密度较大和较轻粒子的压力，在每个快照中应用了两种颜色范围。从快照中可以看出，尽管由于密度差异导致压力突然变化，但密度和压力场仍然保持平滑。下载：下载高分辨率图片（1MB）下载：下载全尺寸图片图12. 在中心力场下振荡液滴仿真中t = 20.0s时的粒子快照（左图：密度；右图：不同颜色范围的压力场（左图：0 – 170 N/m2，右图：0 – 0.55 N/m2）。
图13显示了仿真结果中动能和势能及其总和的时间序列，表3显示了总能量的最大相对漂移。PPE中的相间连接方法两者之间存在显著差异。然而，它们的性能几乎相同。
图13. 在中心力场下振荡液滴仿真中的动能和势能时间序列（上图：动能；中图：势能；下图：它们的总和）。表3. 在中心力场下振荡液滴仿真中总能量的最大相对漂移。
表3. 在中心力场下振荡液滴仿真中的总能量的最大相对漂移。
3.3. 两个上升气泡的合并
为了检验IP在再现两相流动中的大型和复杂界面变形方面的性能，模拟了两个上升气泡的合并 [39]。图14显示了初始条件和带有密度场的粒子快照。当它们通过较密集的流体上升时，两个气泡合并。较密集流体的密度为ρd = 1000.0 kg/m3，较轻流体的密度为ρl = 100.0 kg/m3。粘度系数分别为μd = 0.156 kg/(m·s) 和 μl = 0.078 kg/(m·s)。这里没有考虑表面张力，因此Bond数Bo设置为Bo = ∞。粒子直径d设置为d = 0.005米。为了提高稳定性，应用了粒子位移技术 [29]。数值设置详情列在表4中。由于空间分辨率的原因，在初始配置中圆形气泡的表面出现了不规则性。通过急剧跟踪界面，这些不规则性可能会加剧。为了防止这种情况，将平滑长度设置为一个足够大的值。
图14. 两个上升气泡合并仿真中的粒子快照 [39]（黑线：Level-set方法 [39]，蓝色粒子：当前仿真）。表4. 两个上升气泡合并仿真的数值设置。每个快照都显示了模拟结果：左侧是含IP的情况，右侧是不含IP的传统方法。不含IP时，一些颗粒会在两相界面附近显著分散，这意味着上升气泡的轮廓无法被准确捕捉。相比之下，IP能够提供精确的轮廓。

3.4 堤坝溃塌
在本节中，通过针对坝坝溃塌的二维模拟来检验IP对剧烈多相流的数值稳定性，这些模拟包含了强烈且大幅变形的自由表面流动以及溅射现象。如图15所示，较轻和较重的流体分别由密度为ρl = 1.0 kg/m3和ρd = 1000.0 kg/m3的颗粒组成。粘度系数分别设置为μl = 0.000018 kg/(m·s)和μd = 0.001 kg/(m·s)。较轻的颗粒被视为不可压缩的，并且忽略了表面张力。颗粒大小d设置为d = 0.25m。为了提高稳定性，将颗粒移动技术[29]应用于CISPH-HS-HL-LBP-DS-ADS方法[[40], [41], [42], [43], [44], [45]]。其他数值设置细节列在表5中。
下载：下载高分辨率图像（160KB）
下载：下载全尺寸图像
图15. 堤坝溃塌模拟的初始条件。
表5. 堤坝溃塌模拟的数值设置。
颗粒间距 d [m] 0.25
平滑长度比 h/d 1.5
表面阈值 β [ni/n0] 0.95
密度比 ρd/ρl 1000
核函数 Wendland核函数
时间步长 CFL（Courant数 αdt = 0.1）
PPE求解器类型 ICCG
PPE容忍度 1 × 10??
PPE平均迭代次数 单相：55.1，含IP的多相：297.1
增强方案 CISPH-HS-HL-LBP-DS-ADS-PS

图16显示了含IP与否的单相和多相流模拟中颗粒及其速度向量和密度的快照。图17显示了模拟中的压力分布。不含IP的多相流模拟结果显示，在t = 8.0秒时，当较重流体喷射到静止表面产生冲击力时，压力分布中存在显著噪声。不稳定状态在x = 40米附近出现非物理性的空洞，速度扰动在整个计算域中扩散。即使在t = 11.0秒压力稳定后，这样的非物理性空洞仍然存在于较重流体表面附近，以及远离表面的水箱角落。另一方面，即使在包含散射的轻质和重质颗粒的破碎波的情况下，IP也能实现稳定的模拟。在t = 7.0秒和t = 8.0秒时，溅射的较重颗粒产生了附近轻质颗粒的剧烈局部流动。尽管界面附近的速度场发生了突然变化，压力分布仍然平滑，周围区域保持宁静。当喷射流在x = 40米左右撞击较重相的表面时，轻质流体团块被喷射流捕获并压缩，逐渐增加压力。这延缓了喷射流的崩塌。与单相流模拟相比，较重表面的颗粒之间更连续且紧密地连接在一起。这可以在t = 8.0秒和t = 9.0秒的溅射尖端看到。在多相流模拟中，较轻相和较重相之间的界面不被视为自由表面，而是将其附近的颗粒纳入PPE中，它们的体积被视为连续体的一部分。值得注意的是，模拟使用基于Courant数αdt = 0.1的实际时间步长稳定实现，且界面周围没有明显的数值扩散。

下载：下载高分辨率图像（5MB）
下载：下载全尺寸图像
图16. 堤坝溃塌模拟中颗粒的速度向量和密度快照（左：单相流；中：不含IP的两相流；右：含IP的两相流）。
下载：下载高分辨率图像（5MB）
下载：下载全尺寸图像
图17. 堤坝溃塌模拟中颗粒的压力快照（左：单相流；中：含IP的多相流；右：不含IP的多相流）。

4. 结论
传统的ISPH通过在整个计算过程中保持恒定的颗粒数密度来确保体积守恒。这一框架基于这样一个事实：在压力泊松方程（PPE）中，局部平均密度可以用颗粒数密度来替代。然而，多相流中的不同密度破坏了兼容性，导致与连续性方程不一致。
为了解决这个问题，本研究提出了一种稳定的ISPH多相方案，以实现稳定和准确的模拟。考虑到密度与颗粒数密度之间的关系，引入了一种名为“Interphase Particle”（IP）的伪颗粒来考虑不同相之间的颗粒间力。IP在目标颗粒与其他相中实际存在的邻近颗粒的相互作用中起中介作用。其他相颗粒的相互作用通过IP来考虑。同时，目标颗粒仅考虑与其密度相同或具有相同密度的邻近颗粒之间的颗粒间力。这样可以类似于单相流模拟，通过颗粒数密度来确保体积守恒。IP是一个数学构造，代表了它们之间的相互作用，并且它也被作为实际存在的颗粒纳入PPE矩阵中，以便在同一离散空间内同时求解每个相中流体颗粒的所有压力。这与基于插值技术的网格外方法完全不同。
通过针对两相静水流体柱、在中心力场下振荡的液滴[38]、两个上升气泡的合并[39]以及坝坝溃塌的模拟，展示了IP的增强性能。不同相之间的流体运动在墙壁附近以及流体颗粒附近是不稳定的。由于IP也可以应用于墙壁颗粒，因此这种不稳定性可以得到有效抑制，如在两相静水流体柱的模拟中可以看出。IP在能量守恒方面的性能在振荡液滴的基准测试中得到了验证。尽管相之间的相互作用过程与传统框架完全不同，IP在能量守恒和压力分布方面几乎得到了相同的结果。IP在合并两个上升气泡的基准测试中显示出界面的高再现性。虽然传统方法涉及在不同相的界面周围分散颗粒，但IP保持了连续且平滑的界面轮廓，没有这种分散。IP的数值稳定性在坝坝溃塌模拟中得到了验证。模拟包括了剧烈流动，例如颗粒被破碎波溅射。IP即使在实际时间步长下也能确保模拟的稳定性和可靠性。
对于未来的工作，需要一个表面张力模型[24,50]和一个额外的框架来考虑气相的可压缩性，以实现更真实的气液可压缩-不可压缩多相流[37,51]的物理再现。

CRediT作者贡献声明：
鹤田直树：撰写 - 审稿与编辑，撰写 - 原始草稿，可视化，验证，软件，项目管理，方法论，概念化。
后藤孝文：撰写 - 审稿与编辑，验证，方法论。
渡边茂男：可视化，验证。
阿巴斯·哈耶尔：撰写 - 审稿与编辑，验证，监督。
後藤仁：撰写 - 审稿与编辑，监督，方法论，概念化。