Robust Simulation of Small-Scale Thin Features in SPH-based Free Surface Flows

Robust Simulation of Small-Scale Thin Features in SPH-based Free Surface Flows——TOG 2015

Xiaowei He1, Huamin Wang2;Fengjun Zhang1, Hongan Wang1, Guoping Wang3, Kun Zhou4

1Institute of Software, Chinese Academy of Sciences

2The Ohio State University, USA

3Peking University, China

4Zhejiang University, China

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

Abstract

光滑粒子流体动力学（SPH）是一种高效，质量守恒，并且能灵活处理拓扑变化的方法。然而，由于许多鲁棒性和稳定性问题，在基于SPH的自由表面流中难以模拟小尺度薄特征。在本文中，我们从两个方面解决了这个问题：表面力的鲁棒性和薄特征的数值不稳定性。我们在扩散表面模型下提出了一种基于自由表面能量函数的新的表面张力方案。我们开发了一种不使用空气粒子来计算自由表面流动的空气压力的有效的方法。与以前的表面力公式相比，我们的公式在薄特征情况下对粒子稀疏性更具鲁棒性。为了避免薄特征上的数值不稳定，我们提出通过估计两个尺度的内部压强来调整内部压力，并使用几何感知的各向异性核函数来过滤力。我们的结果证明了我们的算法在处理各种小型薄液体特征方面的有效性，包括薄片，薄喷射和水溅。

Additional Key Words and Phrases: smoothed particle hydrodynamics(SPH), surface tension, diffuse interface model, liquid animation, thin feature

1. INTRODUCTION AND BACKGROUND

小型薄的特征，例如小水流和薄片，在基于物理的液体动画中提供了有趣的细节。但是如何防止它们被分辨率限制和数值不稳定性破坏在计算机图形学中具有挑战性。然而大多数研究工作都用于在基于网格的欧拉模拟器[Losasso et al. 2004; Irving et al. 2006; Kim et al. 2007; Sussman and Ohta 2009]和基于网格的拉格朗日模拟器[Thurey et al. 2010; Wojtan et al. 2010; Brochu et al. 2010; Zhang et al. 2011; Clausen et al. 2013]下解决这个问题，对光滑粒子流体动力学（SPH）及其模拟器的研究很少。实际上，SPH对自由表面流动中液体表面周围的粒子缺乏高度敏感，这使得小尺度的薄特征更难以模拟。由于SPH模拟器因其效率，质量守恒和处理拓扑变化的灵活性而受到许多应用的欢迎，我们认为有必要在基于SPH的自由表面流中稳健地模拟薄特征。

与最近关于基于粒子模拟的分辨率限制的工作不同[Adams et al. 2007; Solenthaler和Gross 2011; Ando et al. 2012; Ando et al. 2013]，我们的工作重点是小规模薄特征的数值方面。具体而言，我们感兴趣的是知道如何提高其鲁棒性，即使在没有足够的粒子的情况下。根据我们的经验，我们发现了与此问题相关的两个主要因素。

第一个因素是表面力，特别是表面张力。在现实世界中，表面张力在维持和破坏小型薄特征方面起着重要作用。在SPH框架下有两种典型的计算表面张力的方法：连续表面力（CSF）方法[Morris 2000;Muller et al. 2003; Hu and Adams 2006]和粒子间相互作用力（IIF）方法[Nugent and Posch 2000; Tartakovsky and Meakin 2005; Becker and Teschner 2007]。通过将表面张力定义为宏观水平的平均曲率流，CSF方法计算每个粒子的表面法线，然后使用光滑核函数来计算表面法线的散度。或者，IIF方法计算微观水平的表面张力作为两个粒子之间的分子间力。虽然这两种方法都足以满足大型水体的需要，但随着粒子数量的减少，它们的精度和鲁棒性都较差，无论表面张力系数如何，都会使薄特征难以保存。

第二个因素是吸引力和排斥力所固有的数值不稳定性。与线性弹簧力不同，基于SPH的吸引力（包括表面张力和空气压力）在粒子靠近时更强，而在粒子分离时更弱。如果它们是唯一的力，它们会将粒子分成若干簇。这个问题通常称为拉伸不稳定性。以前关于拉伸不稳定性的研究主要集中在弹性固体的大变形上[Chen et al. 1999; Monaghan 2000]，所提出的技术不能直接应用于液体模拟。同时，当基于SPH的排斥力（例如​​内部压力）由于数值误差而略微偏离同一直线或平面时，往往将粒子推出。吸引力和排斥力都会导致薄液体特征破裂，例如图6所示的薄片示例。我们注意到数值不稳定性与现实世界表面张力不稳定性不同。它存在于自由表面流动中很大程度上是由于粒子仅定义在自由表面的液体侧。因此，如Schechter和Bridson [2012]所建议的那样，在空气侧添加幽灵粒子有助于减少数值不稳定性，但这需要更多的实现工作和计算成本。

基于这两个观察结果，我们做出了以下贡献，以鲁棒地模拟小尺度的薄特征：

• 在扩散界面模型下从表面能量函数导出的表面张力方案。即使没有足够数量的粒子，该方案也可以稳健地反映液体表面的局部几何形状。
• 仅基于液相中的液体粒子的气压公式。它可以产生各种气压效应，而且计算开销很小。
• 基于双尺度压强估计和几何感知各向异性过滤的内压力算法。它有效地减少了数值不稳定性，而不会影响水体的可压缩性。

我们实现了新方法并将它们集成到标准SPH / WCSPH求解器中。整个系统高效且兼容图形硬件加速。我们的实验表明，它可以逼真而稳健地模拟各种小尺寸的薄特征，如薄喷射（图9），薄膜（图10）和水溅（图8）。

2. PREVIOUS WORK

光滑粒子流体动力学（SPH）已被广泛用于计算物理学和计算机图形学中以模拟动态液体行为。以前的研究主要集中在一些问题上，包括人工粘度[Monaghan 1989; 1994]，不可压缩性[Becker和Teschner 2007; Solenthaler和Pajarola 2009]，边界条件[Muller et al。 2003; Schechter和Bridson 2012]，与其它流体和固体耦合[Monaghan 1994; Muller等人 2005; Solenthaler和Pajarola 2008; Ihmsen等 2010; Akinci等 2012]，以及可变尺寸的粒子[Adams et al. 2007; Solenthaler和Gross 2011;Ando等人 2012; Ando 等人 2013]。

在这些问题中，表面张力及其对小尺寸薄特征的影响是一个研究较少的问题。最初为多相流[Morris 2000]，Muller及其同事[2003]扩展了连续体表面力（CSF）方法，以处理自由表面情况。 Hu和Adams [2006]通过将表面张力表示为应力张量的散度而不是表面法线来改善CSF方法的鲁棒性。 Sirotkin和Yoh [2011]提出了一种新的光滑核函数和梯度校正项，以避免CSF方法中的压缩不稳定性。如Nugent和Posch [2000]所示，基于粒子的表面张力流也可以通过粒子间相互作用力（IIF）方法计算。利用排斥力和吸引力的结合，Tartakovsky和Meakin [2005]使用IIF方法模拟表面张力和流固耦合效应。 Becker和Teschner [2007]应用IIF方法计算自由表面流动中的表面张力。

不幸的是，CSF和IIF的准确度取决于足够数量的粒子。当它们处理具有较少粒子的薄特征时，它们的结果变得不太可靠且噪声更大。另外，Zhang [2010]和Andersson以及合作者[2010]提出重建液体表面以进行表面张力计算。 Yu及其同事[2012]将液体表面随着时间的推移明确地保持为三角形网格。两种方法都比基于粒子的表面张力方法更稳健，但它们需要额外的计算成本。由于自由表面流动中的许多问题不会出现在多相流中，因此直接的想法是在自由表面的空气侧创建幽灵粒子，正如Schechter和Bridson [2012]所示。当场景包含许多薄液体特征时，处理这些新粒子的计算开销可能很大。

液体动画中薄特征的存在也依赖于液体表面重建过程。 Blinn [1982]提出的blobby球方法使用径向基函数之和从标量场提取等值面。 Zhu和Bridson [2005]后来通过增加对局部粒子密度变化的补偿来改进这种方法以减少人工碰撞和压痕。Adams和合作者[2007]提出随时间跟踪粒子表面距离，以便对非均匀粒子的液体表面可以平滑地重建。 Yu和Turk [2010]使用基于局部粒子分布的各向异性光滑核函数，而不是各向同性光滑核函数，以减少表面碰撞而不破坏薄特征。 Bhatacharya及其同事[2011]将液体表面重建表示为约束优化问题，并使用水平集方法来最小化液体表面的薄板能量。 Akinci和合作者[2012]应用网格操作来减少碰撞并有效地提高重建质量。虽然我们的工作专注于数值模拟，但我们的系统也可以从这些液体表面重建技术的使用中受益，以获得更鲁棒的薄特征效果。

3. SURFACE FORCES

在本节中，我们提出了处理基于SPH的自由表面流动的表面张力和气压的新技术。这两种技术都基于扩散界面模型，其历史可以追溯到范德瓦尔斯的早期工作[1893]。扩散界面模型背后的基本思想是假设液体表面具有小而有限的厚度，物理量可以快速且平稳地从一个相位变化到另一个相位。扩散界面中的表面能量可以定义为亥姆霍兹(Helmholtz)自由能函数[Cahn和Hilliard 1958]：
$$\epsilon ={\int }_V[f(c)+\frac{\kappa }{2}|\nabla c{|}^2]dV\text{\hspace{1em}(1)}$$

其中$V$是液体体积，$\kappa$是平方梯度能量系数，$f(c)$是体积自由能密度，$c$是凝结场。通常，如果点在体积内，则凝结值$c$为1，如果一个点在体积之外，则其为0；当一个点穿过表面时，其从1到0光滑变化。直观地，$|\nabla c|$指示接口的位置以及$c$的变化速度。方程1中的平方梯度能量项是表面张力能量
$${\epsilon }^S={\int }_V\frac{\kappa }{2}|\nabla c{|}^{2}dV\text{\hspace{1em}(2)}$$

这与表面积成正比。然后可以将该能量的梯度表示为使表面积最小化的表面张力。扩散界面模型自然地与基于粒子的表示兼容，因为它不需要显式表示液体表面。

3.1 Surface Tension Force

为了使用SPH框架下的扩散界面模型计算表面张力，我们简单地在每个粒子上定义$c=1$，称为色场[Morris 2000; Muller等 2003]。然后我们将${\nabla }_{i}c$计算为：
$${\nabla }_{i}c=\frac{{\sum }_j{V}_j{c}_j{\nabla }_i{W}_{ij}^h}{{\sum }_j{V}_j{W}_{ij}^h}\text{\hspace{1em}(3)}$$

其中$V_j$是粒子$j$的体积，$W_{ij}^h=W(r_{ij},h)$是具有半径$h$的平滑核函数，并且$r_{ij}$是粒子$i$和$j$之间的距离。这里使用分母${\sum }_j{V}_j{W}_{ij}^h$来补偿自由表面流中缺少的空气粒子。使用等式3，我们可以在每个粒子$i$处获得$\frac{\kappa }{2}|{\nabla }_{i}c{|}^2$。通过将其视为每个粒子的光滑能量密度并忽略其他粒子对其的影响，我们假设可以通过分别最小化每个能量密度来最小化${\epsilon }^S$。这个假设允许我们使用能量密度的梯度来定义表面张力：
$$F_i^S=V_i{\nabla }_i(\frac{\kappa }{2}|{\nabla }_{i}c{|}^2)=\frac{\kappa }{2}\sum _j{V}_i{V}_j|\nabla c_j{|}^2{\nabla }_i{W}_{ij}^h\text{\hspace{1em}(4)}$$

直观地，表面张力能量密度$\frac{\kappa }{2}|{\nabla }_{i}c{|}^2$可以被认为是粒子$i$的局部表面积的近似值。通过总计一组吸引力，表面张力试图使其最小化。在内部粒子表面积为零且表面粒子显式指定的理想情况下，我们可以简单地将表面张力视为相邻表面粒子引起的吸引力之和，如图1所示。为了确保实际中的动量守恒，我们计算了两个表面张力能量密度的平均值，并将其用于以下表面张力公式：
$$F_i^S=\frac{\kappa }{4}\sum _j{V}_i{V}_j(|\nabla c_i{|}^2+|\nabla c_j{|}^2){\nabla }_i{W}_{ij}^h\text{\hspace{1em}(5)}$$

我们的方法的主要优点是它对粒子稀疏性的鲁棒性，这在薄特征上是常见的。与依赖于$\nabla c$来确定法线方向和平均曲率的CSF方法不同，我们的方法使用$|\nabla c{|}^2$来仅估计局部表面积。因此，当正常估计成为问题时，例如由单个粒子层组成的液体薄片，我们的方法仍然可以精确地计算表面张力。图2在光滑核函数的支持域包含少于20个粒子的情况下，比较了我们的方法与CSF方法（[Muller等人 2003]）和IIF方法（[Becker和Teschner 2007]）。它表明我们的方法在2D和3D方面都更具鲁棒性。

图1. 三种表面情况下的表面张力。在这些情况下，我们可以将表面张力模型化为表面粒子之间的吸引力之和。它试图将凸面和凹面变形成平面，使表面张力能量最小化。

图2. 三种表面张力方法的2D和3D比较例子。我们通过将液体粒子限制在密闭容器中并施加表面张力和空气压力来模拟凹面实例（将在3.2小节中讨论）。

3.2 Air Pressure Force

由于空气压力，水不能自由地离开固体表面，也不能占据现实世界中的气泡体积。使用液体和空气粒子在多相流中模拟空气压力是很简单的。对于单相自由表面流动，Schechter和Bridson [2012]提出通过在液体表面周围添加幽灵空气粒子以类似的方式计算空气压力。由于使用空气粒子需要更多的内存和计算成本，一个有趣的问题是：我们能否不使用空气粒子模拟空气压力？

要回答这个问题，我们首先假设空气粒子仍然存在。液体粒子$i$应该被空气粒子和液体粒子包围，如图3所示。设$p_{atm}$为每个空气粒子$k$的气压。粒子$i$的空气压力可以计算为：
$$F_i^a=-V_i{p}_{atm}\sum _k{V}_k{\nabla }_i{W}_{ik}^h\text{\hspace{1em}(6)}$$

假设空气粒子和液体粒子光滑均匀分布，我们有：
$$\sum _j{V}_j{\nabla }_i{W}_{ij}^h+\sum _k{V}_k{\nabla }_i{W}_{ik}^h=\nabla 1=0\text{\hspace{1em}(7)}$$

其中第一项对液体粒子求和，第二项对空气粒子求和。使用方程7，我们可以替换方程6中的求和得到：
$$F_i^a=V_i{p}_{atm}\sum _j{V}_j{\nabla }_i{W}_{ij}^h\text{\hspace{1em}(8)}$$

直观地，方程8通过分配相邻的具有负气压的液体粒子，将气压力表示为吸引力的总和。

为了比较该方法与Schechter和Bridson [2012]的ghost SPH方法的性能，我们创建了一个实心球体示例，如图4所示。这个例子表明两种方法都可以产生流动效果（与我们在3.1小节中的表面张力公式一起），其中水在固体表面上流动并在球体底部合并。但由于我们的方法不需要空气粒子，因此它比ghost SPH方法更快，并且其计算成本与薄特征无关。此外，我们的方法在计算气压力时不需要任何额外的存储成本，而ghost方法在该示例中增加了大约25%的存储开销。

图3. 一个表面粒子附近有空气粒子和液体粒子。通过在每个相邻液体粒子处定义负气压，我们可以在不显式定义空气粒子的情况下计算气压力。

图4. 流动的水。虽然ghost SPH方法和我们的方法都可用于模拟实心球体上的流动水，但我们的方法不需要空气粒子并且运行速度更快。

4. NUMERICAL INSTABILITY

尽管第3节中提出的公式可以稳健地计算表面力，但我们仍然可以看到薄的液体特征受到数值不稳定性的影响。这个问题与SPH框架下的吸引力和排斥力有关。

与基于SPH的吸引力相关的不稳定性问题，包括我们的表面张力和我们的空气压力，被称为拉伸不稳定性。为了理解这个问题，让我们考虑一个包含一个可移动粒子$a$和两个固定粒子$b$和$c$的一维情况，如图5a所示。假设粒子具有相同的大小并且它们仅受吸引力的影响。如果$a$恰好位于$b$和$c$的中间，它受到的合力为零并且它可以保持静止。然而，如果$a$由于数值误差而稍微靠近$b$，那么$b$施加在$a$上的吸引力将更大并且$c$施加在$a$上的吸引力将更小。因此总力不平衡并且将$a$推向更接近$b$的位置。在模拟中，这将导致粒子形成一组簇。根据Swegle和合作者[1995]的观点，拉伸不稳定性的存在可以在数学上被确定为一个充分条件：$\sigma W^{"}>0$，其中$\sigma$是应力状态，$W^{"}$是光滑核函数到粒子距离的二阶导数。

虽然基于SPH的排斥力不具有拉伸不稳定性问题，但它们具有其自身的不稳定性问题，如图5b所示。在该示例中，当仅存在排斥力时，粒子$a$可以在$b$和$c$之间保持静止。但是，如果$a$稍微偏离直线，则排斥力将其进一步推出。结果，排斥力不能保持自由表面流动中的薄液体特征，包括薄喷射流和薄板。我们注意到这个问题不会出现在多相流中，因为周围的空气粒子会阻止液体粒子轻易地逃逸薄特征。

由吸引力（或称为拉伸不稳定性）引起的不稳定性可以在任何地方发生，而由排斥力引起的不稳定性仅在薄的特征上发生。在实践中，我们没有注意到水体中的拉伸不稳定性，因为排斥力可以避免粒子任意接近。然而，拉伸不稳定性在薄特征上成为问题，其中粒子较稀疏并且大部分现有模拟器低估了内部压力。

基于前面的分析，我们得出了如下数值不稳定性的解决方案。为了减少拉伸不稳定性，我们首先使用两个光滑核函数增加薄特征上的内部压力。一旦排斥的内部压力变大，它们的不稳定性问题就会变得明显，并且也需要解决。然后，我们对内部压力应用各向异性过滤器，使其影响在薄特征内受到限制。通过以这种方式计算内部压力，我们的系统可以稳健地保持薄的特征，例如图6中的薄板。

图5. 基于SPH的吸引力和排斥力的数值不稳定性问题的一维示例。

图6. 薄片。如果没有通过我们的方法计算内部压力，薄片在小扰动后会破裂成水滴，如(a)所示。使用我们的方法，薄片保持静止，如(b)所示。

4.1 Two-Scale Pressure Estimation

当使用大光滑核函数时，大多数算法无法区分薄特征上的粒子稀疏性和低内压区域中的粒子稀疏性。结果，在薄特征上经常低估内部压力并且不能防止发生拉伸不稳定性。使用小核函数可以降低拉伸不稳定性，但由于计算中涉及的粒子较少，估计的压强噪声更多且不太可靠。

为了可靠地估算水体和薄特征的内部压强，我们的想法是将小型和大型光滑核函数一起使用。设$h=R$是大核函数的半径，$h=r$是小核函数的半径。我们首先使用Solenthaler和Pajarola [2008]提出的多相方法计算密度：
$${\rho }_i^r=\alpha m\sum _j{W}_{ij}^{r}and{\rho }_i^R=m\sum _j{W}_{ij}^R\text{\hspace{1em}(9)}$$

其中$m$是粒子质量，$\alpha$是比例因子，三维空间中为$(r/R{)}^3$，二维空间中为$(r/R{)}^2$。我们通常设置$R=2.5d$，$r=d$，其中$d$是两个粒子之间的预期参考距离。然后我们使用局部泊松方法[He et al. 2012]将两个密度分别转换为两个压强$p_i^r$和$p_i^R$。为了提高效率，我们只应用积分域的半径减少为零的特殊公式。对于水体中的粒子，压强$p_i^R$更准确，但对于薄特征上的粒子则低估了压强。同时，压强$p_i^r$不太可靠，但它对薄特征没有低估的问题。为了提供从一个到另一个的平滑过渡，我们计算粒子$i$的最终内部压强为：
$$p_i=max(p_i^R,\beta p_i^r)+p_{atm}\text{\hspace{1em}(10)}$$

系数$\beta$在方程10中有两个目的。首先，它确保水体中粒子的压强不受可能含有噪声的$p_i^r$的影响。其次，它控制了薄特征上的排斥内压力的大小，因此它们不会抑制有吸引力的表面力，尤其是表面张力。由于表面张力与表面张力能量密度有关，我们使用经验方程来定义$\beta$：
$$\beta =\gamma \underset{i}{\max }(\kappa |{\nabla }_{i}c{|}^2)/\underset{i}{\max }(p_i^r)\text{\hspace{1em}(11)}$$

其中$\gamma$需要位于$[0,0.5]$的范围内，以获得合理的模拟结果。 $\beta$的使用可以被认为是表面张力效应和薄特征之间的平衡。较小的$\beta$使表面张力效应更显著，但破坏了薄特征，而更大的$\beta$保留了薄的特征，但削弱了表面张力效应。我们注意到$\beta$完全是从算法的角度引入的，它没有物理意义，但它为我们提供了一种调整视觉效果的便捷方法。我们将通过在结果部分中将$\gamma$设置为不同的值来演示$\beta$如何影响流体行为。

4.2 Anisotropic Pressure Filtering

在我们使用4.1小节中的双尺度压强估算方法修正拉伸不稳定性之后，我们现在必须面对由内部压力增加引起的不稳定性。这种不稳定性问题经常被$p_i$中的噪声夸大，即使在方程10中使用小$\beta$之后也是如此。

受Yu和Turk [2010]提出的各向异性表面重建方法的启发，我们通过在内部压力上应用各向异性过滤器来解决这种不稳定性。令$C_i={\sum }_j(x_j-x_i)(x_j-x_i{)}^T{W}_{ij}^R$是在粒子$i$处定义的各向异性协方差矩阵。我们提出张量矩阵$T_i$为：
$$T_i=\frac{p_i^R}{p_i}I+(1-\frac{p_i^R}{p_i})\frac{C_i}{||C_i|{|}_2}\text{\hspace{1em}(12)}$$

对于水体中的粒子，$p_i=p_i^R$，$T_i$是单位矩阵。对于薄特征上的粒子，$p_i$大于$p_i^R$并且$T_i$变得更加各向异性。使用这个张量矩阵，我们将内部压力表示为：
$$F_i^p=-\frac{1}{2}\sum _j{V}_i{V}_j(p_i{T}_i+p_j{T}_j)\cdot {\nabla }_i{W}_{ij}^h\text{\hspace{1em}(13)}$$

直观地，方程13减小了垂直于薄特征的方向上的内部压强效应。以这种方式，内部压力消除了拉伸不稳定性，而不会由于其自身的不稳定性问题而破坏薄特征。对于水体中的粒子，方程13将简化为没有各向异性过滤的标准公式。

5. RESULTS

（请参考动画结果的补充视频。）我们实现了我们的方法，并将它们集成到标准的SPH / WCSPH求解器中，所有的光滑核函数都是根据Muller及其同事[2003]的工作选择的。我们使用Ihmsen及其同事[2011]提出的并行索引排序算法来构建统一网格并加速邻域搜索。我们使用Bhatacharya及其合作者[2011]提出的水平集方法进行液体表面重建。为了建模固体物体，我们均匀地创建靠近固体表面的固体粒子，并将它们视为幽灵液体粒子。在这种情况下，我们忽略固 - 液表面张力并将固 - 空表面张力视为液 - 空表面张力。我们在具有6GB内存的四核Intel Xeon W3550 3.07GHz工作站上测试我们的系统。我们在模拟中将时间步长设置为$\Delta t=0.001s$。表I显示了我们的例子的系数和时间（每个时间步长），包括粒子数N，表面张力有效$\kappa$，表面力的平均计算时间$t^s$，内部压强的平均计算时间$t^P$和总计算时间$t^{tot}$。

水射流。图9比较了使用不同表面张力系数的模拟效果。在这个例子中，水射流分解成具有相同体积但表面积较小的小水滴。通常称为Plateau-Rayleigh不稳定性，当表面张力系数从左向右增加时，这种效应更可能发生。

在小兔子上的巧克力。图7展示了不同表面力的模拟效果。在不使用表面张力的情况下，液体粒子可以自由地离开小流并形成许多液滴，如图7a所示。在不使用空气压力的情况下，小流不会在固体表面上流动，如图7b所示。如图7c所示，通过使用两种表面力，我们获得了物理上更合理的效果。

水冠。图8比较了使用不同$\gamma$值的流体行为。在$\gamma =0$的极端情况下，由于表面张力引起的数值不稳定，一束小液滴将被夹断。在$\gamma =1$的另一个极端情况下，这意味着压强将被过度校正，流体不能显示合理的表面张力行为。为了平衡表面张力效应和薄特征，一个合适的选择是设置$\gamma =0.2$，我们能够得到合理的模拟结果。

球面上的水。在这个例子中，我们提出另一个例子来证明我们的内部压力算法对于保留薄特征的重要性，如图10所示。在这里我们直接设置$\gamma =0.2$并忽略气压力使视觉差异更明显。在不调整内部压力的情况下，水膜迅速破裂成碎片。使用我们的技术后，在模拟中水膜一直保持，直到它被分辨率限制破坏。

图7. 兔子上的巧克力。此示例演示了使用不同表面力的效果。在大多数例子中，我们使用表面张力和气压力。

图8. 水冠。此示例显示了通过将$\gamma$设置为不同的值，压强校正如何影响流体行为。

图9. 水射流。该示例显示了使用不同表面张力系数（以$N/m$为单位）的Plateau-Rayleigh不稳定性效应。

图10. 球体上的水。在这个例子中，我们比较了不使用（顶部）和使用（底部）内部压力算法之间的差异。该算法有效地避免了数值不稳定性问题，但它无法阻止薄特征被分辨率限制破坏。

6. LIMITATIONS

我们的系统没有解决分辨率限制问题，因此它无法处理尺寸小于单个粒子的薄特征。它不能用于为大型水体中的气泡设置动画，因为它确实考虑了空气的可压缩性。它既不考虑固体表面性质也不考虑润湿效应，因此它不能模拟复杂的固体对液体的影响，例如疏水效应。尽管我们的表面张力公式比CSF和IIF方法更具鲁棒性，但它仍然是近似值，并且它不像ghost SPH方法那样准确[Schechter和Bridson 2012]。我们算法中的一些系数（包括$\kappa$和$\beta$）不是基于物理学的，它们需要针对不同的例子进行调整。最后，如何在液体表面重建过程中保留薄特征仍然是一个难题，因为薄的特征可能被错误地识别为噪声。使用较小的光滑核函数可以保留薄特征，但它可能不足以消除实际噪声，正如我们在一些示例中所发现的那样。

7. CONCLUSION AND FUTURE WORK

在本文中，我们发现表面力和数值不稳定性是影响基于SPH的自由表面流动中小尺度薄特征的两个主要因素。我们证明了自由表面能函数在形成表面张力方面的作用，我们研究了在不使用空气粒子的情况下处理气压效应的可能性。我们提出了一种计算薄特征内部压力的新算法，有效地减少了数值不稳定性问题。

我们接下来的计划是测试我们的系统与图形硬件加速的兼容性。我们也有兴趣将它与基于网格的跟踪方法相结合[Yu et al. 2012]，模拟我们当前系统无法处理的效果。从长远来看，我们希望探索利用自由表面流动为大型水体中复杂的液 - 固相互作用和气泡制作动画的可能性，其中大部分只能通过过去的多相流模拟。