Implicit Incompressible Porous Flow using SPH

一句话总结

本文提出首个基于 SPH 的隐式不可压多孔流（porous flow）求解器：通过一套考虑孔隙率的新密度估计让流体与固体粒子直接”重叠”，把孔隙率纳入压力求解器以严格保持不可压性，并用强耦合的隐式非压力力（拖曳、毛细、浮力）取代传统的显式 Darcy 近似，从而稳定地模拟海绵吸水、挤压排水、毛细上升等复杂多孔交互。

研究背景

• 领域现状：自然界大多数固体（木头、土壤、羊毛、海绵）都含有可让流体进入并流动的孔隙，这种多孔流会显著影响两相的宏观行为。SPH 作为纯粒子的拉格朗日方法，天然适合统一离散化各种材料与交互力，是多物理场仿真的常用工具。
• 核心痛点：图形学中已有的 SPH 多孔流方法（如 Lenaerts、Ren 等）通过在固-流界面处”缩放流体粒子体积”来表示吸收——粒子跨越界面时不断缩小/增大。这带来两个严重问题：其一，界面处相邻粒子质量差异巨大，导致隐式对偶压力求解器条件数很差；其二，内部流动只用 Darcy 定律显式求解，不保证两相之间的力平衡，多孔效应靠压力场的启发式手段拼凑，缺乏动量守恒。快速多孔流场景下粒子体积剧变还会引发”爆炸”式失稳。
• 本文 idea：既然缩放粒子是万恶之源，就干脆不缩放。作者观察到，只要允许流体与固体粒子占据同一空间（重叠域方法），并重新推导一个把孔隙空间考虑进去的密度估计，就能用标准 SPH 压力求解器直接对重叠区域强制不可压。在此基础上，把拖曳、毛细、浮力等交互力都写成对下一时刻速度线性、可组装进一个强耦合隐式线性系统的形式，配合固体弹性与流体黏性一起求解，从而兼顾稳定性、动量守恒与大时间步。

方法

孔隙密度估计

方法的地基是重新定义”密度”。采用连续介质视角，用孔隙率 \(\phi \in [0,1)\) 度量固体内部自由空间：\(\phi=0\) 为无孔材料，\(\phi \to 1\) 为纯空隙。

标准 SPH 中粒子密度为核加权求和 \(\rho_i = \sum_j m_j W_{ij}\)，其中 \(W_{ij}=W(x_i-x_j,h)\) 是光滑核。作者的关键改动在于区分两类粒子的采样体积含义（见正文对体积定义的说明）：固体粒子 \(s\) 的采样体积 \(V_s^0\) 包含空隙，但只承载不含孔的固体质量 \(m_s=(1-\phi)\rho_s^0 V_s^0\)。据此得到多孔固体密度：

\[\rho_s = (1-\phi)\rho_s^0 \sum_{t\in\mathcal{N}_s^S} V_t^0 W_{st}\]

它逼近的是多孔物体整体密度 \((1-\phi)\rho_s^0\)，而非固体相本身密度。对不可压多孔物体 \(\rho_s\) 应保持常量；对可压材料（如海绵），孔隙可被挤出，允许 \(\rho_s\) 增大，但绝不超过固体相密度 \(\rho_s^0\)（此时孔隙全部塌陷，再压就是压固体本身，方法禁止之）。

流体密度则不同：流体采样体积 \(V_i^0\) 完全被流体填满。当流体与固体重叠时，需扣除固体相自身占据的空间。借鉴 number density 思路，但只把固体的”无孔部分”（乘因子 \(1-\phi\)）计入流体密度估计：

\[\hat{\rho}_i = \rho_i^0 \sum_{j\in\mathcal{N}_i^F} V_j^0 W_{ij} + \rho_i^0 \sum_{t\in\mathcal{N}_i^S} (1-\phi)V_t^0 W_{it}\]

若附近无固体，第二项消失并退化回标准公式，因此 \(\hat{\rho}_i\) 在多孔域内外都保持一致。孔隙率 \(\phi\) 正好限制了流体能与固体重叠多少，从而决定了流体粒子在多孔体内的间距（吸收上限）。

压力求解器

通过求解压力泊松方程强制不可压：

\[\Delta t \nabla^2 p = \frac{\rho^0 - \rho^*}{\Delta t}\]

其中 \(\rho^*\) 是仅由非压力力预测的密度。求解器需要两样东西：由速度得到的密度变化率，以及压力梯度产生的加速度。前者对 \(\rho_s\) 与 \(\hat{\rho}_i\) 取时间导数得到（本质是速度散度的 SPH 近似，与连续性方程一致）：

\[\frac{D\hat{\rho}_i}{Dt} = -\rho_i^0 \sum_{j\in\mathcal{N}_i^F} V_i^0 v_{ji}\cdot\nabla W_{ij} - \rho_i^0 \sum_{t\in\mathcal{N}_i^S}(1-\phi)V_t^0 v_{ti}\cdot\nabla W_{it}\]

压力加速度方面，流体把固体当作边界粒子处理，采用对称梯度形式并用修正后的采样体积：

\[a_i^{press} = -\rho_i^0 \sum_{j\in\mathcal{N}_i^F} V_j^0 \left(\frac{p_i}{\rho_i^2}+\frac{p_j}{\rho_j^2}\right)\nabla W_{ij} - \rho_i^0 \sum_{t\in\mathcal{N}_i^S}(1-\phi)V_t^0 \frac{p_i}{\rho_i^2}\nabla W_{it}\]

固体粒子此阶段不考虑来自流体的压力力（浮力单独作为耦合力处理，以便与弹性求解器直接耦合）。利用 SPH 常见的压力钳制（只把压缩视为违反连续性、右端为负时压力才激活），作者巧妙地用”替换右端静密度项”的方式在可压/不可压多孔物体间切换：可压材料允许 \(\rho_s < \rho_s^0\)，不可压材料则只允许 \(\rho_s < (1-\phi)\rho_s^0\)。整体用改进的 IISPH 求解，松弛 Jacobi 迭代、每步钳制压力为正。

多孔交互力

固体粒子用带孔隙耦合力 \(f_s^{pore}\) 的 Cauchy 动量方程，流体粒子则把孔隙力加到 Navier-Stokes 动量方程：

\[\frac{Dv_i}{Dt} = -\frac{1}{\hat{\rho}_i}\nabla p + \frac{\mu_i}{\hat{\rho}_i}\nabla^2 v + \frac{1}{m_i}f_i^{pore} + \frac{1}{m_i}f_i^{ext}\]

孔隙力被分解为最常见的三种效应之和：毛细 \(f^{cap}\)、拖曳 \(f^{drag}\)、浮力 \(f^{buo}\)。由于多孔体内拖曳往往很强、显式求解极不稳定，所有非压力力都写成对 \(v^{n+1}\) 线性的隐式形式。

毛细作用：简化自 Jeske 等人的黏附模型，用毛细系数 \(C^{cap}\) 表达流固间吸引：

\[f_i^{cap} = -\sum_{t\in\mathcal{N}_i^S} C^{cap}(S_t^n)\frac{\bar{m}_{it}}{\bar{\rho}_{it}} x_{it}^{n+1} W_{it}^n\]

该向量指向固体邻居最密集的方向，把流体拉入多孔体。饱和度 \(S_s\) 度量局部孔隙被流体填充的程度：

\[S_s = \frac{1}{\phi}\frac{\sum_{j\in\mathcal{N}_s^F} V_j^0 W_{sj}}{\sum_{t\in\mathcal{N}_s^S}\frac{m_t}{\rho_t}W_{st}}\]

并钳制到 \([0,1]\)。毛细势随饱和度按落差因子 \(\eta^{cap}\in[0,1]\) 衰减：

\[C^{cap}(S_s) = C^{cap}_0 (1 - S_s \eta^{cap})\]

其中 \(C^{cap}_0\) 决定毛细可达的最大高度，\(\eta^{cap}\) 则控制下部饱和区被拉入多少流体。为写成速度线性，用 \(x^{n+1}=x^n+\Delta t\,v^{n+1}\) 展开，把力拆成依赖 \(v^{n+1}\) 的项与只依赖已知量的右端项 \(b^{cap,n}\)。固体侧的力由动量守恒 \(f_{i\leftarrow s}^{cap}=-f_{s\leftarrow i}^{cap}\) 得到。

拖曳：借鉴 Bui 和 Nguyen 把拖曳当作黏性摩擦，写成相对速度拉普拉斯的近似：

\[f_i^{drag}(v^{n+1}) = \tilde{d}\sum_{t\in\mathcal{N}_i^S}\frac{\mu^{por}}{1-\phi}\frac{m_i m_t}{\rho_i^n \rho_t^0} v_{it}^{n+1}\cdot\frac{x_{it}^n}{\|x_{it}^n\|^2+0.01h^2}\nabla W_{it}^n\]

其中 \(\mu^{por}\) 为多孔黏性系数，\(\tilde{d}=2(d+2)\)，\(d\) 为空间维数；形式类似 Navier-Stokes 黏性项，并改造 Weiler 等人的黏性法以在不同质量粒子间保持力对称。

浮力：某些多孔材料因周围流体压力而漂浮。作者不把这个力放进流体压力求解器（否则流体会无视弹性力直接推开固体），而是把固体视为固定、用最近一次压力解得的流体压力梯度镜像出作用在固体上的力 \(b_s^{buo,n}\)，因此流体本身不再额外受浮力（\(f_i^{buo}=0\)）。

非耦合力与孔隙固体效应

流体黏性沿用 Weiler 等人的隐式求解器（用流体相密度 \(\hat{\rho}\)），固体弹性用 Peer 等人的共旋线性弹性求解器（用 \(V_s^0\) 和 \(\rho_s\)）。为体现饱和引起的力学变化，作者对应力张量做扩展：

\[\sigma_s = 2\mu(S_s)\epsilon_s + \lambda(S_s)\mathrm{tr}(\epsilon_s)\mathbf{1} + \sigma_s^{bloat}\]

膨胀（bloating）项按饱和度加入，模拟吸水后孔隙膨胀或吸湿溶胀：

\[\sigma_s^{bloat} = -\eta^{bloat} S_s \mathbf{1}\]

Lamé 系数也随饱和度调整，以刻画木头吸水变软等现象：

\[\mu(S_s) = (1+\eta_m S_s)\mu_0,\quad \lambda(S_s) = (1+\eta_l S_s)\lambda_0\]

变化因子为正增大、为负减小抗变形能力（后者需钳制保正）。

求解器整体流程

每个时间步先算密度（三套密度公式）与饱和度；再把所有非压力力组装成对 \(v^{n+1}\) 线性的系统：

\[M v^{n+1} - \Delta t\, f(v^{n+1}) \approx M v^n + \Delta t\, b^n\]

其中 \(M\) 为质量矩阵，\(f\) 与右端 \(b\) 是各交互力与非耦合力对应项之和。该系统用无矩阵共轭梯度法求解，并用上一步的速度增量作为初值。解得速度后预测密度、经压力加速度修正，再用辛欧拉法积分得到新位置。

实验与结论

• 稳定性：在 1ms 时间步、\(\mu^{por}=10\,\mathrm{Pa\,s}\) 下，显式拖曳无法充分减速流体、粒子直接穿透多孔块；增大系数又导致失稳。隐式力即便在较小系数下也能明显减速并保持稳定。与 Ren 等人方法对比，在快速多孔流场景中对方因界面处密集粒子团的压力尖峰频繁”爆炸”，本文方法因完全避免了粒子剧烈缩放而在快慢两种情形都稳定。代价是本文更慢（慢速例每步约 194ms vs 对方 126ms），因为对方省去了被吸收流体的压力计算且只做单向耦合。
• 验证：梯形坝内流体输运实验中，本文无需毛细即可复现与 Casagrande 解析解一致的自由面，而 Ren 等人方法难以独立控制流速与毛细高度、且底部流体被过度压缩；多孔球下落实验中，方法依固体相密度正确再现漂浮/下沉行为。
• 参数研究：吸收流体量与孔隙率 \(\phi\) 近似线性；多孔黏性 \(\mu^{por}\) 控制均匀渗入速度；毛细 \(C^{cap}_0\) 控制吸收高度与速度，\(\eta^{cap}\) 控制下部湿区被拉入的流体量。含黏附的例子最耗时（每步约 373.7ms，其中压力求解占 69%、强耦合非压力力占 20%）。
• 视觉展示：海绵吸水后被挤压排水、fusilli 形多孔物吸水变软、多孔兔子在螺旋桨扰动湍流中互动（配合微极模型）、以及毛细驱动流体沿字母上行等复杂艺术场景。方法已实现于 SPlisHSPlasH 库（C++，AVX+OpenMP），作者计划开源。
• 局限：目前仅支持共旋线性弹性固体，未来可换非线性弹性或颗粒材料求解器；尚未研究非线性拖曳模型；压力求解器是速度瓶颈，改用 divergence-free SPH 或可提速但需支持多孔吸收/排出时的预期速度散度。

启示

• 把”孔隙率”直接编进密度估计与压力求解器，是让两相合法重叠而不牺牲不可压性的关键，避开了粒子缩放这一长期不稳定源头——一个”改地基”胜过”打补丁”的范例。
• 将拖曳、毛细、浮力、黏性、弹性统统写成对下一时刻速度线性的形式并组装进单一强耦合隐式系统，是隐式方法能承受大拖曳力、支持大时间步的核心工程手段。
• 用动量守恒的物理耦合力替代 Darcy 定律的启发式压力场，在物理自洽性上带来实质提升，代价是压力求解成为性能瓶颈，指明了后续优化方向。