IISPH-FLIP for incompressible fluids

IISPH-FLIP for incompressible fluids——2014 EG

Jens Cornelis    Markus Ihmsen    Andreas Peer    Matthias Teschner

University of Freiburg, Germany

本文是对原文的翻译，供有需求者参考，翻译不当之处请指正，欢迎互相交流。

Abstract

我们提出使用隐式不可压缩光滑粒子流体动力学（IISPH）进行流体 - 隐式 - 粒子（FLIP）求解器中的压强投影和边界处理，用于模拟不可压缩流体。这种新颖的结合解决了现有SPH和FLIP求解器的两个问题，即FLIP中的质量守恒以​​及SPH中的效率和内存消耗。首先，SPH部分可以模拟具有完美质量守恒的不可压缩流体。其次，FLIP部分有效地丰富了SPH部分的细节，其细节与具有相同数量粒子的标准SPH模拟相当，同时将性能提高了7倍，并显著降低了内存消耗。我们证明了所提出的IISPH-FLIP求解器可以模拟不可压缩流体，具有可量化的，难以察觉的低于0.1%的密度偏差。我们展示了使用仅15GB内存在单台台式PC上处理的多达1.6亿粒子的大规模场景。并展示了单向和双向耦合固体。

1. Introduction

流体 - 隐式 - 粒子（FLIP）概念广泛用于视觉效果，代表了高质量，多功能，强大和高效流体模拟的最新技术[ATW13]。 FLIP由Brackbill和Ruppel [BR86，BKR88]提出，并且由Zhu和Bridson [ZB05，Bri08]在图形界引入了改进版本。在FLIP中，流体用粒子表示，而投影步骤在辅助欧拉网格上执行。多功能效果已集成到FLIP中，例如，与固体[BBB07]和多相流[BB12]的双向耦合。就内存消耗和计算时间而言，与其他粒子方法相比，例如光滑粒子流体动力学（SPH），FLIP具有概念上的优势。首先，根据[ATT12]，FLIP处理的时间步长比SPH更大。其次，原始FLIP概念不需要查找，存储和处理相当大量的相邻粒子。第三，可以灵活地处理FLIP粒子。可以轻松生成或删除它们以保留流体薄片，拆分或合并以考虑自适应网格，或重新定位以解决不均匀的空间粒子分布[ATT12，ATW13]。所有这些操作都旨在尽可能减少对流体动力学的影响，同时它们会在SPH模拟中造成严重的稳定性问题。

现有的FLIP求解器使用粒子来处理对流。然而，压强投影，体力和边界处理[BBB07]通常在固定[ZB05，ATT12]或自适应欧拉网格[BR86，ATW13]上进行处理。作为欧拉网格的替代方案，我们建议在拉格朗日网格上进行压强计算，即SPH粒子使用IISPH [ICS* 13]。这种IISPH-FLIP组合的特点是精确的质量守恒，这是使用欧拉网格的FLIP求解器中的典型问题。在欧拉压力投影中，避免甚至量化由于质量损失导致的体积变化是相当困难的，例如， [CM12]，我们表明，在本文提出的求解器中，由于密度偏差导致的体积变化可以精确量化并有效地最小化。此外，[GB13]中讨论的由于FLIP粒子的重播引起的潜在质量增加不再是问题，因为质量守恒完全由IISPH压强投影处理。在所有呈现的IISPH-FLIP实验中，体积变化处于低于0.1%的难以察觉的水平。

本文提出的技术是FLIP概念的新颖变体。将FLIP粒子的速度插值到辅助拉格朗日网格上，即SPH粒子，其中通过用IISPH求解压强泊松方程来计算无散速度场。将来自连续步骤的速度差插值回去并添加到最终对流的FLIP粒子的速度。

由于提出的求解器将IISPH集成到FLIP求解器中，其实现遵循FLIP概念。然而，压强样本的位置，即SPH粒子位置不是固定的，而是以它们的速度对流。SPH和FLIP粒子之间的插值如FLIP中那样执行。压强投影用IISPH [ICS* 13]计算。

现有的FLIP方法在粗欧拉网格上使用大量FLIP粒子产生了的高度详细的模拟。对于IISPH-FLIP也有同样的效果，其中粗采样的SPH粒子和大量FLIP粒子之间的插值丰富了流体模拟。

这就引出了对所提出的求解器的第二个观点。它也可以被视为SPH或IISPH扩展，其中SPH采样中增加了额外的FLIP粒子和FLIP插值策略。实验结果表明，在粒子数相同的情况下，IISPH-FLIP的整体流动和细节与IISPH相当，但与IISPH相比，性能显著提高，内存需求降低。在具有15GB内存的单台PC上计算了具有高达1.6亿个FLIP粒子的场景（图1）。

图1：具有1.6亿个FLIP粒子的场景。粒子对速度进行颜色编码。

2. Related work

在图形应用程序中，将不同组件组合成一个新解算器相当流行，并且通过一个解决器的优点来证明其在某些情况下解决了其他解算器的缺点，反之亦然。一种流行的方法是SPH粒子与网格求解器的组合，例如[LTKF08，LHK09，GLHB09，ZYF10]。这些方法在适当的情况下用SPH表示流体域的部分，而域的其他部分在网格上求解。两个域通常是双向耦合的。与这些方法相反，我们提出用SPH粒子表示整个流体域，并用FLIP粒子丰富整个域。

其他方法结合了多种分辨率。例如，[LZF10]在粗网格上计算压强，如果需要，在精细网格上细化。在SPH的背景下，[RWT11]也在粗网格上计算压强。对于高分辨率SPH粒子，所产生的压力与基于WCSPH [BT07]或PCISPH [SP09]的压力相结合。与这些方法相反，我们提出用IISPH计算SPH粒子上的低分辨率压强场，该压强场与高分辨率FLIP粒子结合。

我们的方法适用于两种空间分辨率。然而，它与空间自适应方法的关系不大，例如[APKG07，SG11，ZLC* 13]，其中域的不同部分用不同的细节水平表示。相比之下，我们的方法使用低分辨率SPH粒子和高分辨率FLIP粒子对整个流体进行采样。

我们的方法与FLIP密切相关，例如[BR86，ZB05，ATW13，GB13]，我们提出了一种新的变体。我们表明，IISPH可用于压强投影，有效解决标准FLIP中的质量和体积守恒问题。我们的方法也与SPH密切相关，例如[Mon92，DC96，MCG03，Mon05，BT07，SP09，HLWW12，ICS* 13，IOS* 14]。在这里，我们展示了我们的方法有效地丰富了标准的IISPH的细节。

众所周知，当需要强制执行小密度偏差时，基于状态方程的SPH求解器的性能较低，如[SP09]中所述。然而，[SP09，HLWW12]中提出的预测矫正公式缓解了这一问题，[ICS* 13]进一步表明IISPH优于预测矫正公式[SP09，BLS12]。3.3节讨论了进一步激发其在FLIP压强投影中使用的IISPH属性。

3. Basics

我们提出在FLIP框架内使用IISPH概念进行压强计算，这也可以看作是带有FLIP粒子的SPH扩展。为了完整起见，简要介绍了这两个部分。

3.1. FLIP

标准FLIP [BR86]使用FLIP粒子和辅助粗欧拉网格。使用相邻FLIP粒子速度的加权平均值加上由于体力引起的速度变化来初始化网格节点处的速度。然后，压强投影计算速度变化，从而在网格上产生无散速度场。相应的压强泊松方程通常用例如改进的不完全Cholesky预处理器的共轭梯度求解[Bri08，GB13]（MICCG），而边界处理可以与[BBB07]结合。由于体力和压强投影引起的速度变化从网格插值回FLIP粒子并且加到它们之前的速度上。这样得到的粒子速度通常与[Bri08]中提出的插值网格总速度混合。然后，对流FLIP粒子。粒子运动在固体边界处被剪除。

3.2. IISPH

IISPH [ICS* 13]是一个不可压缩的SPH变体，它在SPH中实现了标准投影方法[Cho68]。在SPH粒子速度上加上由于体力和粘度引起的速度变化。通过求解压强泊松方程并应用由压强梯度产生的相应速度变化来校正该中间速度。最后，SPH粒子以校正的速度对流。 IISPH在泊松方程的源项中使用密度不变性条件以避免体积漂移。由于体积漂移在IISPH中不是问题，因此少量的求解器迭代就足够了。边界处理基于[AIA* 12，IAGT10，ACAT13]并入IISPH。

3.3. IISPH properties

IISPH在泊松方程的源项中使用密度不变性条件，该条件能够校正之前引入的密度或体积偏差。这与无散度条件形成对比，其不考虑之前时间步长的偏差。对于无散条件，产生的体积漂移只能通过相当大量的求解器迭代来抵消。相比之下，体积漂移在IISPH中不是问题，并且少量的求解器迭代就足够了。 [ICS* 13]讨论了具有多达4000万个压强样本（SPH粒子）和0.025秒的时间步长的场景，其中IISPH求解器平均仅需要4.1次迭代来求解具有难以察觉的密度偏差的泊松方程。 IISPH可以无矩阵方式实现，每个粒子的内存占用量为7个标量值。这些性质激励IISPH成为FLIP中有效的体积守恒压强投影的一个有趣候选者。第5节讨论了与标准FLIP [Bri08]中压强预测的性能比较。

4. IISPH-FLIP

在欧拉流体模拟中，FLIP作为一个处理对流项的概念被引入。因此，标准FLIP中的投影步骤是在欧拉网格上计算的。然而，也可以采用一组拉格朗日样本位置，例如SPH粒子，用于FLIP中的压强计算，如[BKR88]中简要提到的那样。这种组合显然不是很有用，因为在没有FLIP粒子的拉格朗日方法中可以平凡地处理对流。尽管如此，拉格朗日SPH和FLIP粒子的组合具有上述关于FLIP中的质量守恒和SPH中细节的有效计算的益处。在下文中，我们描述了所提出的IISPH-FLIP求解器的概念，然后是实现细节。

4.1. Concept

我们的方法适用于位置为$x_F$，速度为$v_F$的高分辨率FLIP粒子，位置为$x_S$，速度为$v_S$的低分辨率SPH粒子。根据[ZB05]，每个压强样本，即SPH粒子，使用8个FLIP粒子。如果不相关，则省略时间索引。此外，使用具有紧凑支持的核函数$W_{ij}=W(||x_i-x_j||,2h)$，其依赖于FLIP或SPH粒子i和j的距离。我们采用支持$2h$的三次样条函数 [Mon92]。 SPH粒子的间距是$h$。

对于每个时间步，通过使用核函数$W_{i_S{j}_F}$插值相邻FLIP粒子$j$的速度$v_{j_F}$来重新初始化SPH粒子$i$的中间速度$v_{i_S}^{\ast }$：
$$v_{i_S}^{\ast }=\frac{1}{{\sum }_{j_F}{W}_{i_S{j}_F}}\sum _{j_F}{v}_{j_F}{W}_{i_S{j}_F}\text{\hspace{1em}(1)}$$

与场量的标准SPH插值相比，求和必须归一化。这是因为FLIP粒子的实际体积或密度未知。由于SPH粒子携带速度，因此不一定需要这种初始化。但是，该过程减少了IISPH模拟的数值耗散。此外，产生的内聚效果类似于标准FLIP，因此我们更喜欢这种额外的插值。

现在，IISPH处理SPH粒子。根据投影或分裂概念，IISPH通过求解泊松方程和速度$v_{i_S}(t+\Delta t)$来计算SPH粒子的压强。边界处理已合并到IISPH中。在我们的场景中，容许密度偏差始终设置为0.1%。与标准IISPH相比，不需要人工粘度。

在下一步中，将SPH粒子的速度插值回FLIP粒子上。这里，我们遵循标准技术，例如[ZB05，Bri08，ATT12，GB13]，其中混合了FLIP和PIC插值。当FLIP插值速度差时，PIC会将总速度插值到FLIP粒子上。考虑到核函数支持中的所有SPH粒子，FLIP粒子处的PIC速度插值计算为
$$v_{i_F}^{PIC}=\frac{1}{{\sum }_{j_S}{W}_{i_F{j}_S}}\sum _{j_S}{v}_{j_S}(t+\Delta t)W_{i_F{j}_S}\text{\hspace{1em}(2)}$$

于是，FLIP速度插值计算为
$$v_{i_F}^{FLIP}=v_{i_F}(t)+\frac{1}{{\sum }_{j_S}{W}_{i_F{j}_S}}\sum _{j_S}(v_{j_S}(t+\Delta t)-v_{j_S}^{\ast })W_{i_F{j}_S}\text{\hspace{1em}(3)}$$

两个速度都用正则化参数$\alpha =min(\frac{6\Delta tv}{h^2},1)$混合，该参数与时间步长$\Delta t$，SPH粒子间距$h$和流体粘度$\nu$相关，在我们的场景中设置为0.01：
$$v_{i_F}(t+\Delta t)=\alpha v_{i_F}^{PIC}+(1-\alpha )v_{i_F}^{FLIP}\text{\hspace{1em}(4)}$$

最后，使用更新速度的欧拉步更新FLIP和SPH粒子的位置：
$$x_{i_F}(t+\Delta t)=x_{i_F}(t)+\Delta t{v}_{i_F}(t+\Delta t)\text{\hspace{1em}(5)}$$
$$x_{i_S}(t+\Delta t)=x_{i_S}(t)+\Delta t{v}_{i_S}(t+\Delta t)\text{\hspace{1em}(6)}$$

该算法如算法1所示。

4.2. Implementation

邻域搜索：我们使用z-索引排序和紧凑哈希[IABT11]的组合来加速邻域搜索。所有FLIP和SPH粒子都存储在相同的均匀网格中，边长为$2h$，因为在IISPH中使用相同的$2h$支持的相同三次样条核以及FLIP和SPH粒子之间的速度插值。在压强计算之前，针对每个SPH粒子查询FLIP和SPH邻居。不存储~240个FLIP邻居，因为它们仅在等式1中使用一次。在迭代的IISPH求解器中需要存储约~40个SPH邻居。在压强投影之后，第二个查询找到每个FLIP粒子的~40个SPH邻居。同样，这些邻居不会被存储，因为它们仅用于等式2，3中速度插值。

边界处理：我们用粒子表示固体边界，并使用[AIA* 12]的概念来处理SPH粒子与边界粒子的单向和双向耦合。 FLIP粒子的边界处理通过方程4中的空间自适应$\alpha$来实现。如果FLIP粒子的相邻SPH粒子处于边界粒子的影响半径内，我们设置$\alpha =1$。该技术缓解了FLIP粒子在边界处的聚集。这与标准FLIP形成对比，其中FLIP粒子的路径通常在边界处被剪除。

粒子属性：SPH粒子具有IISPH所需的标准属性。然而，FLIP粒子需要较少的属性。特别是，不存储对FLIP或SPH邻居的引用。这与无矩阵的IISPH实现一起，实现了高效存储，例如， 15GB处理1.5亿粒子（图1）。

孤立的粒子：如果SPH粒子没有FLIP邻居，则方程1无法计算。相应地，$v_{i_S}^{\ast }$设定为孤立的SPH粒子的速度。如果FLIP粒子没有邻居，则将重力作用于其速度。

5. Result

我们实现了IISPH，标准FLIP和IISHP-FLIP，并在5.1节中比较模拟结果，在5.2节中比较性能。IISPH和IISPH-FLIP使用三次样条核函数[Mon05]。在IISPH中，我们另外模拟[Mon05]中的粘滞力和[BT07]中的表面张力。注意，我们不需要在IISPH-FLIP中应用额外的粘度。表面张力是可选的，没有应用于所提出的IISPH-FLIP模拟。我们对标准FLIP求解器的实现遵循[ZB05]，每个压强样本使用8个FLIP粒子，并采用变分方法进行边界处理，如[BBB07]。使用MICCG [Bri08]求解压强系统。这些图像是用Houdini的Mantra渲染器[Sid13]渲染的。随附的视频以每秒50帧的速度编码。

为了比较视觉特性和实现的求解器的性能，我们建立了[KFV* 05]中规定的溃坝场景。

5.1. Visual comparison

我们首先将标准FLIP与IISPH进行比较，以证明两个求解器的关系和特性。

我们相信这有助于评估所提出的IISPH-FLIP方法的好处。用所有考虑的方法模拟的溃坝场景的相应帧如图2所示。

图2：障碍物允许质量损失的视觉比较。对于所考虑的方法，可以感知到速度场的明显差异。粒子对速度进行颜色编码。

IISPH vs FLIP

IISPH中的压强计算在SPH粒子上执行，而FLIP在辅助网格上计算压强。因此，在IISPH中，压强样本之间的距离对应于SPH粒子间距$h$，而在FLIP中，它对应于的网格点的距离$\Delta x_{grid}$。FLIP模拟中的细节通过所采用的FLIP粒子的额外自由度明显增加。我们通过在(i)相同分辨率的压强场和(ii)相同数量的FLIP粒子比较IISPH和FLIP来强调这种效果。

通过将SPH粒子间距$h$设定为FLIP中的网格间距$\Delta x_{grid}=1.85cm$，获得相同分辨率的压强场。因此，粒子的数量不同：IISPH中96k SPH粒子，FLIP中772k FLIP粒子。虽然压强是以相同的分辨率计算的，但低分辨率的IISPH模拟具有明显更高的数值耗散，从而导致阻尼。这可以从较低的涡度和较慢的波对流中观察到。

通过将SPH粒子间距减半至$h=0.925cm$，将SPH粒子的数量增加至与FLIP模拟中相同数量的FLIP粒子。对于该分辨率，IISPH模拟的数值粘度水平与FLIP模拟相当，可以通过比较不同帧的速度场来评估。波浪似乎也在以相同的速度移动。

图2还说明了FLIP遭受压缩。相比之下，IISPH在两种分辨率下都保留了流体体积。对于IISPH，容忍的体积压缩设置为0.1%。

IISPH-FLIP vs FLIP

对于相同数量的FLIP粒子和相同分辨率的压强场，IISPH-FLIP的模拟结果看起来类似于FLIP。 IISPH-FLIP中的数值耗散水平与标准FLIP相当。然而标准FLIP遭受质量损失并因此遭受体积压缩，但我们的IISPH-FLIP变体保留了流体体积，具有0.1%的不可感知的误差。

IISPH-FLIP vs IISPH

尽管与高分辨率IISPH模拟相比，IISPH-FLIP中仅使用了压强样本数量的$1/8$，但可以观察到类似的视觉细节和涡度。相比之下，低分辨率IISPH遭受更高的数值耗散，更少的细节和速度场平滑。因此，与低分辨率的IISPH模拟相比，IISPH-FLIP可以获得更好的模拟结果，同时粒子数相同的情况下，性能比高分辨率的IISPH模拟高出7倍。

Interpolation radii

IISPH-FLIP使用三种不同的支持半径：$r$用于IISPH核函数，$r_{F\to S}$用于从FLIP粒子到SPH粒子的速度插值，$r_{S\to F}$用于从SPH粒子返回到FLIP粒子的速度插值。根据IISPH核支持$r$，我们遵循[Mon92]并使用$r=2h$结合三次样条函数和粒子间距$h$。

与核函数支持$r$相比，速度插值的支持半径$r_{F\to S}$和$r_{S\to F}$对稳定性不太重要，但会影响性能和视觉结果。我们使用图3所示的溃坝场景对$r_{F\to S}$和$r_{S\to F}$的各种组合进行了实验。

对于$r_{F\to S}=2h$，从FLIP粒子到SPH粒子的速度插值需要0.82s并且对于$r_{F\to S}=h$减小到0.60s。对于从IISPH粒子到FLIP粒子的所得速度的插值也是如此，其中$r_{S\to F}=2h$的计算需要0.97s且$r_{S\to F}=h$仅为0.71s。每帧的总计算时间范围为，对于$r_{F\to S}=r_{S\to F}=h$，2.62s，对于$r_{F\to S}=r_{S\to F}=2h$，3.1s。

如图3所示，对于较小的插值半径，视觉外观不太乐观，因为FLIP粒子往往粘附到它们最接近的SPH粒子上，导致FLIP粒子的爆米花状模式。因此，我们更偏向使用$r_{F\to S}=r_{S\to F}=2h$进行速度插值，以获得视觉上令人满意的结果。

图3：FLIP和SPH粒子之间速度插值的支持半径的不同组合。

5.2. Performance comparison

我们使用所有线程在6核Xeon X5680 3.33 GHz，24GB RAM上模拟了的所有求解器的溃坝场景。性能测量总结在表1中。

表1：两种不同分辨率的IISPH、FLIP和IISPH-FLIP的性能比较。粒子数指的是用于IISPH的SPH粒子和用于FLIP和IISPH-FLIP的FLIP粒子。

所有的模拟计算均使用自适应时间步长保持CFL数为1。因此，实际时间步长受压强采样点的距离限制：对于IISPH和IISPH-FLIP为$h$，对于FLIP为$\Delta x_{grid}$。因此，对于$\Delta x_{grid}=h=1.85cm$的IISPH-FLIP，FLIP和低分辨率IISPH，使用相同的时间步长，而对于$h=0.925cm$的高分辨率IISPH模拟，它被减半。

IISPH-FLIP使用相同数量的压强样本（SPH粒子）为IISPH模拟增加了计算开销，因为必须将量从SPH粒子插值到FLIP粒子并返回。对于给定的场景，使用IISPH-FLIP计算帧的时间比在低分辨率的IISPH模拟中长2.5倍。 IISPH-FLIP每帧占用3.1s。该时间0.24s用于邻域搜索，0.82s用于从FLIP粒子到SPH粒子的速度插值，1.07s用于IISPH压强投影，0.97s用于将所得速度插值回FLIP粒子。相比之下，低分辨率的IISPH每帧需要1.31s。

但是，如5.1节所述，附加程序显著增加了视觉细节，使模拟甚至可与高分辨率的IISPH结果相媲美。与具有相同数量粒子的IISPH模拟相比，IISPH-FLIP显著提高了性能，因为压强场的复杂性降低了8倍，实际时间步长可以加倍。对于使用772k SPH粒子的IISPH模拟，每帧需要21.4s，这比IISPH-FLIP所需的3.1s几乎慢7倍。

IISPH-FLIP优于我们的标准FLIP实现。这是因为FLIP求解器的性能受到许多方面的影响，这些方面在我们的实现中尚未优化。但是，应该注意的是，IISPH平均仅在9次迭代中执行压强计算，而使用MICCG的标准FLIP实现平均需要110次迭代。

5.3. Large-scale scenario

IISPH的一个主要限制是其内存消耗。我们的方法克服了这个限制，因为FLIP粒子不需要存储与压强投影相关的任何属性。不存储FLIP粒子邻居。

图1中的场景由1.6亿个FLIP粒子组成。在我们的IISPH实现中，这个场景的内存消耗为48GB，超过了当前标准工作机器的限制。相比之下，我们的IISPH-FLIP实现在这个场景中仅消耗15GB，在16核3.46 GHz Intel i7上平均每帧花费10分钟来计算。在可接受的时间内不可能模拟具有相同数量的SPH粒子的场景。

5.4. Two-way coupling

如4.2节所述，无需在固体边界处约束FLIP粒子轨迹。隐式处理边界碰撞，并且不必查询FLIP粒子的邻域以获得边界粒子，可以有效地执行水动力的计算而无需任何额外的邻域查询。在图4所示的场景中用40M FLIP粒子演示了耦合固体的实例。虽然没有明确考虑FLIP粒子的固体相互作用，但没有FLIP粒子穿透固体边界，并且船逼真地漂浮在流体粒子上。

图4：双向耦合船，逼真地漂浮在40M FLIP粒子上。粒子没有穿透固体物体。

6. Conclusion and future work

我们已经展示了如何将压强投影和边界处理的IISPH合并到标准FLIP模拟器中。通常，质量守恒是FLIP模拟中的一个问题。由于SPH的质量守恒，我们提出的IISPH-FLIP可以解决这个问题并模拟不可压缩的流体。另一方面，SPH通常会受到数值耗散的影响，特别是低分辨率模拟。添加FLIP粒子有助于保持涡度并为IISPH添加细节。虽然IISPH-FLIP显示的行为和细节程度与具有一半压强采样距离的IISPH模拟相当，但在这样的模拟中它的性能高出7倍。

我们已经表明，通过在空间调整流体的粘度可以隐式地处理FLIP粒子的边界处理。我们不需要在边界处约束粒子运动，只需依靠SPH粒子即可处理固体物体的双向耦合。

由于FLIP粒子不需要存储压强投影所需的属性，并且可以在线计算邻域而不用为存储邻居和距离消耗内存，因此我们的IISPH-FLIP模拟具有较低的内存占用。我们展示了一个在一台台式PC上计算1.6亿个FLIP粒子的场景，只消耗了15GB的内存。

由于我们同时使用了FLIP和SPH，因此可以研究其他技术，例如用于SPH的涡度约束或用于FLIP粒子的粒子重播。