Divergence-Free Smoothed Particle Hydrodynamics

Divergence-Free Smoothed Particle Hydrodynamics——SCA 2015

Jan Bender∗           Dan Koschiery

Graduate School CE

TU Darmstadt

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

Abstract

在本文中，我们介绍了一种有效且稳定的隐式SPH方法，用于不可压缩流体的物理模拟。在计算机图形领域，最有效的SPH方法仅关注密度误差的校正以防止体积压缩。然而，不可压缩流动的连续性方程也需要无散的速度场，大多数方法都忽略了这一点。虽然有几种方法考虑了速度散度，但它们要么很慢，要么具有可感知的密度波动。

我们的新方法使用两个压强求解器的有效结合，实现低体积压缩（低于0.01%）和无散速度场。这可以看作是在位置级别和速度级别上强制执行不可压缩性。第一部分对于真实的物理行为至关重要，而无散状态显著增加稳定性并减少求解器迭代次数。此外，它允许更大的时间步长，从而获得相当大的性能增益，因为粒子邻域更新频率减小。因此，我们的无散SPH（DFSPH）方法比不可压缩流体的现有最先进的SPH方法明显更快且更稳定。我们在数百万个快速移动的粒子的模拟中证明了这一点。

Keywords: fluid simulation, Smoothed Particle Hydrodynamics, divergence-free fluids, incompressibility, implicit integration

1 Introduction

在过去几年中，光滑粒子流体动力学（SPH）成为计算机图形学中模拟复杂水效果的常用方法。 SPH是一种无网格拉格朗日方法，它通过仅考虑一组有限的相邻粒子来计算给定点处的流体量。该领域中最具挑战性的问题之一是强制执行不可压缩性，这对于模拟真实的物理行为至关重要。

一般来说，不可压缩流体的SPH模拟是使用拉格朗日坐标系中的不可压缩，等温Navier-Stokes方程进行的。
$$\frac{D\rho }{Dt}=0\iff \nabla \cdot \mathbf{v}=0\text{\hspace{1em}(1)}$$
$$\frac{D\mathbf{v}}{Dt}=-\frac{1}{\rho }\nabla p+\nu {\nabla }^2\mathbf{v}+\frac{\mathbf{f}}{\rho }\text{\hspace{1em}(2)}$$

其中$D(\cdot )/Dt$表示物质导数，$\rho ,p,\nu ,\mathbf{v}$和$\mathbf{f}$分别表示密度，压强，运动粘度，速度和体力。注意，由于不可压缩性，密度相对于时间的偏导数为零，即$\partial \rho /\partial t=0$，这意味着(1)中两个方程的等价。根据这些方程，不可压缩流体必须满足无散条件$\nabla \cdot \mathbf{v}=0$，这意味着流体必须具有无散速度场。根据连续性方程$\frac{D\rho }{Dt}=-\rho \nabla \cdot \mathbf{v}$和无散条件，密度必须随时间保持恒定$\frac{D\rho }{Dt}=0$。因此，理论上，具有无散速度场的流体具有恒定的密度，因此是不可压缩的。然而，在实践中，在模拟中强制实现无散条件不足以保证不可压缩性。数值时间积分的不可避免的误差导致密度偏差，这总结了模拟。由于无散条件没有考虑所产生的密度误差，因此无法避免由于数值误差引起的体积压缩。为了校正数值密度误差，必须满足第二个条件$\rho -\rho 0=0$，我们称之为恒定密度条件。总之，不可压缩流体的Navier-Stokes方程需要无散速度场，另外需要稳定形式的恒定密度条件来抵消数值误差。

最近，在计算机图形领域，隐式压强求解器在模拟不可压缩流体的SPH模拟中变得流行。这些求解器确定粒子模型的压力，以防止体积压缩。目前，用于不可压缩流体的最有效的SPH压强求解器仅考虑了仅依赖于粒子位置的恒定密度条件。然而，得到的速度场通常不是无散的（见图2，左），这是不可压缩流动的连续性方程所要求的。到目前为止，只有少数SPH方法考虑速度偏差，但它们要么不能强制执行低密度误差，要么它们相对较慢。

在本文中，我们介绍了一种新的高效稳定的SPH方法，与大多数以前的方法相比，它实现了无散条件和恒定密度条件。这是通过两个压强求解器的结合来实现的，这两个压强求解器分别考虑了散度误差和密度误差。我们的第一个求解器强制实现无散速度场（见图2，右），这在SPH模拟中具有几个优点。首先，模拟变得更加稳定，特别是在快速移动粒子的模拟中（见图1）。增加的稳定性允许我们执行更大的时间步长，这增加了整体性能，因为不用频繁地计算昂贵的邻域搜索。此外，可以通过取决于最大粒子速度的CFL条件更精确地确定最大时间步长。最后，强制无散速度场显著减少了纠正密度误差的第二个求解器所需的迭代次数。已经存在不同的恒定密度求解器，它们可以与我们的无散解算器结合，例如，预测 - 纠正不可压缩SPH（PCISPH）[Solenthaler和Pajarola 2009]和隐式不可压缩SPH（IISPH）[Ihmsen等 2014a]。然而，在这项工作中，我们介绍了一种新方法，它完全适合我们的无散解算器，因为它使用模拟技术，因此可以重用预先计算的系数。与当前最先进的SPH方法相比，使用我们的新方法进行的模拟速度提高超过20倍。

图1：我们新的SPH方法可以稳定模拟高速不可压缩流体，同时保持无散速度场。左侧的模拟中显示了240万个流体粒子和600万个边界粒子。此外，我们的方法明显快于目前最先进的SPH方法，能够在5秒的时间步长模拟由500万个流体粒子和4000万个边界粒子组成的复杂场景，最大体积压缩率为0.01%（右）。

图2：在125k粒子的溃坝模拟中，IISPH（左）和DFSPH（右）的速度散度之间的比较。散度误差是颜色编码的，其中白色是最小值，红色是最大值。

2 Related Work

在计算机图形学领域中，开发了用于模拟不可压缩流体的各种方法。在本节中，我们将简要介绍相关方法。对于一般调查，我们参考了Bridson [2008]的工作，而Ihmsen等人[2014b]最近的最新报告专门讨论了用于流体模拟的SPH方法。

Monaghan [1994]使用状态方程通过使用刚度加权压强项惩罚密度误差来弱强制不可压缩性。他们根据声速选择刚度系数，以保持可观的可压缩性。 Desbrun和Gascuel [1996]利用基于状态方程（EOS）的求解器来模拟高度可变形的物体。后来，Muller等人 [2003]在计算机图形学界引入了基于EOS的流体模拟方法，而Adams等人 [2007]通过空间自适应性扩展了这种方法。 Becker和Teschner [2007]提出了一种EOS方法，通过预定的，依赖于场景的刚度系数来限制最大压缩。然而，小密度容差导致高刚度系数导致刚性微分方程，因此导致严格的时间步长限制。

为了进一步加强密度守恒，可以使用在没有压力的情况下对流粒子后确定的中间密度来计算压强，概念上称为分裂[Chorin 1968; Bridson 2008]。 Solenthaler和Pajarola [2009]采用了这种分裂概念，并通过迭代压强求解器对其进行了扩展，以便将最大密度误差保持在用户定义的容差范围内。后来，该方法通过刚性流体耦合[Akinci et al. 2012]和一种新型表面张力模型[Akinci et al. 2013]扩展。 Macklin和Muller [2013]提出了基于位置的流体（PBF）—— 一种类似的概念，它迭代地细化粒子位置以强制执行不可压缩性。 Bodin等人 [2012]在速度水平上解决了由刚体力学引起的密度误差的完整约束。然而，由于数值近似，正则化参数和时间积分误差，他们报告压缩误差高达17%。所有这些基于EOS的分裂迭代方法都只能在速度或密度水平上强制执行不可压缩性，而我们的方法同时考虑了这两种。最近，Kang和Sagong [2014]提出了一种方法，通过额外满足无散条件来扩展Macklin和Muller [2013]的工作。然而，随着粒子位置随着速度投影而更新，它们不能保证在每个时间步之后的速度场无散。此外，他们报告的运行时间与Macklin和Muller的方法相似或更差，而我们的方法产生的加速倍数高达一个数量级，特别是对于大时间步长（见表2）。

与非迭代和迭代EOS求解器相比，可以求解压强泊松方程（PPE）将中间速度投影到无散状态。特别是，在基于网格的欧拉方法中，这是常见的做法[Bridson 2008]。因此，Cummins和Rudman [1999]使用多重网格方法求解PPE，其中粒子被视为最精细的网格。后来，Foster和Fedkiw [2001]采用半拉格朗日方法，其中在模拟网格上强制执行不可压缩性，同时使用水平集和自由移动的无惯性粒子抵消质量耗散。在接下来的几年中，开发了使用粒子和背景网格的进一步混合方法[Zhu and Bridson 2005; Raveendran等 2011; Ando等 2013]。 Losasso等人[2008]提出了一种先进的方法，其中PPE在背景网格上求解，不仅保持无散场，而且保持预定义的目标密度。 Sin等人 [2009]引入了另一种有趣的方法，他们在每个步骤中的点样本构建的Voronoi图上求解PPE。但是，所有这些方法都必须维护并有时重建额外的网格数据结构，从而导致高内存消耗，特别是对于大型域。我们的方法完全避免了背景网格，因此消耗较少的内存，特别是对于如图1所示的大型场景。

一些方法直接求解无网格离散化的PPE。 Hu等人[2007]提出了一种用于多相流的不可压缩SPH方法，其中获得速度的时间积分步骤被细分为两个半步。在第一个半步期间，使用位置变化迭代梯度下降求解器消除了密度波动。在第二个半步期间，速度散度的误差以类似的方式解决。然而，他们将该方法仅应用于多达14,400个粒子的非复杂二维场景。He等人 [2012]通过局部泊松求解器扩展了Solenthaler和Parajola [2009]的工作，实现了无散速度场和恒定密度。不幸的是，他们必须在每个求解器迭代中确定其泊松求解的粒子邻居，而积分域不一定相等或者是局部核函数支持的子集。由于搜索半径增加，这导致相当大的计算量。此外，与Solenthaler和Parajola的基础方法相比，他们仅报告了大约1.5的小加速倍数，而我们使用DFSPH实现了超过20倍的加速。在Ihmsen等人[2014a]最近的一项工作中，使用PPE以便将密度偏差迭代地降低至0.1%或甚至更小。但是，与我们的方法相比，他们没有考虑消除速度散度导致大量的求解器迭代，甚至可能在特殊情况下导致伪影（见第4节）。

3 Fluid Simulation

在我们的模拟中，我们使用不可压缩的等温NavierStokes方程（参见第1节）。但是，我们跳过方程(2)右侧的第二个粘性项，并使用Schechter和Bridson [2012]提出的XSPH变量来模拟粘度。我们使用如下所述的SPH方法在空间上离散化方程(2)，同时通过分别实现3.2和3.3节中描述的无散条件和恒定密度条件来满足方程(1)所表示的不可压缩性条件。

使用SPH概念，位置${\mathbf{x}}_{\mathbf{i}}$处的量通过一组相邻粒子${\mathbf{x}}_{\mathbf{j}}$的值近似[Monaghan 1992]。SPH公式中最重要的量之一是密度，可以使用以下概念近似：
$${\rho }_i=\sum _j\frac{m_j}{{\rho }_j}{\rho }_j{W}_{ij}=\sum _j{m}_j{W}_{ij}$$

其中$m_j$是粒子$j$的质量，$W_{ij}=W(x_i-x_j,h)$是具有支撑半径$h$的类高斯核函数。在计算密度之后，流体的压强场通常由以下状态方程（EOS）确定：
$$p_i=\frac{\kappa {\rho }_0}{\gamma }((\frac{{\rho }_i}{{\rho }_0}{)}^{\gamma }-1)$$

其中$\kappa$和$\gamma$是刚度参数，${\rho }_0$是静止密度。在我们的工作中，我们考虑特殊情况
$$p_i=\kappa ({\rho }_i-{\rho }_0)\text{\hspace{1em}(3)}$$

其中$\gamma =1$是计算机图形学中的常见选择[Desbrun and Gascuel 1996; Muller等 2003; Ihmsen等 2014b]。

在标准SPH方法[Monaghan 1992]中，由EOS计算的压强场直接用于确定压力来抵消体积压缩。然而，减小压缩流体需要很高的刚度系数$\kappa$，这限制了时间步长并因此限制了整体性能。因此，开始研究隐式压强求解器来模拟不可压缩的流体，例如[Solenthaler和Pajarola 2009; Ihmsen等 2014a]。这些求解器通常通过求解线性系统来计算压强场，该线性系统允许更大的时间步长并因此获得显著的性能增益。

在本文中，我们提出了一种新的隐式SPH方法用于不可压缩流体。我们的方法分别使用两个求解器来满足无散条件和恒定密度条件。无散解算器的目标（见3.2节）是获得无散速度场。但是，如第1节所述，这不足以保证实际中的不可压缩性，因为无法校正由数值误差引起的密度偏差。因此，我们采用恒定密度求解器（参见第3.3节）作为消除密度误差的稳定器。求解器的关键思想是确定每个邻域的刚度系数${\kappa }_i$（见式(3)），用于局部求解相应的条件。为了满足全局条件，求解器以并行的Jacobi方式处理邻域，这是SPH求解器中的常见选择[Solenthaler和Pajarola 2009; Macklin和Muller2013; Ihmsen等 2014a]。单独调整刚度系数相当于压强场的隐式计算。然而，这使我们能够得到一个简单而一致的公式来表示恒定密度和无散力。

3.1 Simulation Step

在模拟步骤结束时，必须满足恒定密度条件和无散条件。无散求解器计算压力，该压力被积分一次以获得满足相应条件的速度变化。由恒定密度求解器确定的压力必须积分两次，以获得所需的位置变化。因此，作为副作用，它也会改变速度。因此，首先，我们执行恒定密度求解器并修改速度和位置，然后执行无散解算器，校正得到的速度以获得无散状态。由于两个步骤都是在一个循环中执行的，因此在计算无散速度场之前执行密度稳定不会产生任何限制。

算法1详细描述了我们的新方法的模拟。首先，我们使用紧凑哈希来确定粒子邻域$N_i$ [Ihmsen et al. 2011]。然后，对于每个粒子，我们计算密度${\rho }_i$和因子${\alpha }_i$（参见第3.2节）。${\alpha }_i$是两个求解器的公因子，分别是校正密度误差和散度误差所必需的。由于该因子仅取决于当前位置，因此在执行求解器之前对其进行预计算，这显著降低了两个求解器的计算量。

模拟循环的第一步是确定所有非压力$F^{adv}$，如重力，表面张力和粘度，并根据Courant-Friedrich-Levy（CFL）条件$\Delta t\le 0.4\frac{d}{||v_{max}||}$调整时间步长 [Monaghan 1992]，其中$d$是粒子直径，$v_{max}$是最大粒子速度。在12行中，非压力用于计算预测速度$v_i^{\ast }$。恒定密度求解器使用该预测速度和预先计算的因子${\alpha }_i$来确定每个邻域的压力，以便校正密度误差${\rho }_i^{\ast }-{\rho }_0$（参见第3.3节）。然后，这些位置将及时更新。因此，必须更新邻域$N_i$，密度${\rho }_i$和因子${\alpha }_i$。在第21行中，无散解算器计算压力以满足$\frac{D{\rho }_i}{Dt}=0$，以使速度场无散（见第3.2节）。最后，所得到的速度变化用于更新粒子速度。

注意，每个模拟步骤仅确定一次$N_i$，${\rho }_i$和${\alpha }_i$。但是，我们不会像在之前的工作中那样在模拟循环的开始计算这些值，因为我们的两个求解器在不同的时间点执行。相反，必须在第一个模拟步骤之前第2-6行中初始化这些值。然后在第16-20行的每个时间步骤中更新它们。

3.2 Divergence-Free Solver

我们的无散求解器满足条件$\frac{D\rho }{Dt}=0$，这意味着密度不随时间变化。这相当于无散条件（见方程(1)），因此得到无散速度场。

如果对于粒子$i$不满足条件$\frac{D{\rho }_i}{Dt}=0$，则求解器计算粒子及其邻域的一组压力，用于校正散度误差。粒子$i$的压力由下式确定
$$F_i^p=-\frac{m_i}{{\rho }_i}\nabla p_i\text{\hspace{1em}(4)}$$

压强梯度通过微分方程(3)w.r.t来计算。$x_i$使用SPH公式[Ihmsen et al. 2014b]：
$$\nabla p_i={\kappa }_i^v\nabla {\rho }_i={\kappa }_i^v\sum _j{m}_j\nabla W_{ij}$$

其中${\kappa }_i^v$是我们想要确定的刚度参数。此外，我们考虑来自粒子$i$作用于相邻粒子$j$上的压力$F_{j\leftarrow i}^p$，以获得满足条件$F_i^p+{\sum }_j{F}_{j\leftarrow i}^p=0$的一组对称压力。这意味着所有内力总和为零，这是动量守恒所需的。力$F_{j\leftarrow i}^p$类似于等式(4)计算，除了我们将压强相对于相邻位置$x_j$进行求导：
$$F_{j\leftarrow i}^p=-\frac{m_i}{{\rho }_i}\frac{\partial p_i}{\partial x_j}=\frac{m_i}{{\rho }_i}{\kappa }_i^v{m}_j\nabla W_{ij}\text{\hspace{1em}(5)}$$

求解器必须确定改变粒子速度的压力，以满足条件$\frac{D{\rho }_i}{Dt}=0$。粒子$i$的当前密度变化率是通过采用SPH散度公式计算的[Ihmsen等 2014b]：
$$\frac{D{\rho }_i}{Dt}=\sum _j{m}_j(v_i-v_j)\nabla W_{ij}\text{\hspace{1em}(6)}$$

在应用为粒子$i$确定的对称压力之后，该值应为零。压力导致速度变化$\Delta v_i=\Delta t{F}_i^p/m_i$，$\Delta v_j=\Delta t{F}_{j\leftarrow i}^p/m_i$。在等式(6)中插入这些项得到：
$$\frac{D{\rho }_i}{Dt}=\Delta t\sum _j{m}_j(\frac{F_i^p}{m_i}-\frac{F_{j\leftarrow i}^p}{m_i})\nabla W_{ij}\text{\hspace{1em}(7)}$$

在等式(7)中使用等式(4)和(5)得出刚度参数${\kappa }_i^v$的等式：
$$\frac{D{\rho }_i}{Dt}=\Delta t\sum _j{m}_j(\frac{F_i^p}{m_i}-\frac{F_{j\leftarrow i}^p}{m_i})\nabla W_{ij}$$
$$\frac{D{\rho }_i}{Dt}=-\frac{\Delta t}{{\rho }_i}\sum _j{m}_j({\kappa }_i^v\sum _j{m}_j\nabla W_{ij}+{\kappa }_i^v{m}_j\nabla W_{ij})\nabla W_{ij}$$
$$\frac{D{\rho }_i}{Dt}=-{\kappa }_i^v\frac{\Delta t}{{\rho }_i}(|\sum _j{m}_j\nabla W_{ij}{|}^2+\sum _j|m_j\nabla W_{ij}{|}^2)$$

解出${\kappa }_i^v$:
$${\kappa }_i^v=\frac{1}{\Delta t}\frac{D{\rho }_i}{Dt}\cdot \frac{-{\rho }_i}{\underset{{\alpha }_i}{\underbrace{|\sum _j{m}_j\nabla W_{ij}{|}^2+\sum _j|m_j\nabla W_{ij}{|}^2}}}\text{\hspace{1em}(8)}$$

其中${\alpha }_i$是仅取决于当前位置的因子。用刚度参数${\kappa }_i^v$计算的压力恰好满足条件$\frac{D{\rho }_i}{Dt}=0$，这意味着粒子$i$附近的速度场是无散的。然而，由于相邻粒子的刚度参数彼此依赖，因此它们是迭代确定的。注意，分母${\alpha }_i$可能导致粒子$i$具有较少邻居的特殊情况下的不稳定性。为了解决这个问题，如果它太小，我们约束分母。在我们的模拟中，我们使用的阈值为$10^{-6}$，这不会导致任何视觉伪影。

最后，我们确定粒子$i$的合力$F_{i,total}^p$，包括来自相邻粒子$j$的力
$$F_{i,total}^p=F_i^p+\sum _j{F}_{i\leftarrow j}^p=-m_i\sum _j{m}_j(\frac{{\kappa }_i^v}{{\rho }_i}+\frac{{\kappa }_j^v}{{\rho }_j})\nabla W_{ij}$$

其中$F_{i\leftarrow j}^p$类似于等式(5)计算。注意，该压力等价于Monaghan [1992]引入的对称压力。

我们的求解器使用Jacobi迭代并行确定压力，以使完整的速度场无散。由于因子${\alpha }_i$仅取决于当前位置，因此它们可以在迭代过程之前预先计算，并且不必在每个迭代步骤中更新。这使得迭代步骤的计算更简单，因为${\alpha }_i$是我们的求解器中使用的最复杂的项。注意，如果粒子$j$不是动态的，则$F_{j\leftarrow i}^p=0$，因此必须相应地针对静态边界粒子调整${\kappa }_i^v$的等式。

算法2概述了我们的无散求解器。求解器执行至少一次迭代，当平均密度变化率小于用户定义的阈值${\eta }^v$时结束。通过执行求解器的“热启动”可以显著提高收敛性。这意味着我们为每个粒子计算刚度值${\kappa }_i^v$。然后在下一个时间步骤中，我们首先使用结果值为每个粒子计算第7行，然后启动无散解算器。

3.3 Constant Density Solver

虽然上一节中描述的求解器使速度场无散，但我们的恒定密度求解器使密度误差最小化，密度误差由实际密度与静止密度的偏差$\rho -{\rho }_0$确定。目前存在各种使密度偏差最小化的压强求解器，它们可以与我们的无散解算器结合，例如， PCISPH [Solenthaler和Pajarola 2009]或IISPH [Ihmsen等 2014a]。然而，在这项工作中，我们引入了一个新的求解器，它与我们的无散解算器类似，因此具有可以重用因子$\alpha$的优点。这显著地减少了计算工作量，因为如果因子$\alpha$已知，则所提出的迭代方法特别快。我们的新恒定密度求解器采用预测校正器方案，以便在校正密度误差的时间积分之后获得粒子位置。这个方案的关键思想类似于PCISPH。但是，与PCISPH相比，我们不使用带有填充邻域的预计算的原型配置来求解系统。

我们使用等式(6)对密度执行欧拉积分步骤，以便计算密度误差${\rho }_i^{\ast }-{\rho }_0$的预测：
$${\rho }_i^{\ast }={\rho }_i+\Delta t\frac{D{\rho }_i}{Dt}={\rho }_i+\Delta t\sum _j{m}_j(v_i^{\ast }-v_j^{\ast })\nabla W_{ij}$$

类似于等式(7)，我们通过求解下式来确定校正该密度误差的压力：
$${\rho }_i^{\ast }-{\rho }_0=\Delta t^2\sum _j{m}_j(\frac{F_i^p}{m_i}-\frac{F_{j\leftarrow i}^p}{m_i})\nabla W_{ij}\text{\hspace{1em}(9)}$$

这将得到刚度参数：
$${\kappa }_i=\frac{1}{\Delta t^2}({\rho }_i^{\ast }-{\rho }_0){\alpha }_i$$

算法3中概述了隐式压强求解器。注意，与无散解算器类似，我们执行热启动来提高求解器的收敛性。

3.4 Kernel

SPH模拟中使用的核函数是一个高斯近似。在以前的工作中，引入了几个核函数，例如poly6核函数，尖峰核函数和三次样条核函数。在一些工作中，甚至使用不同的核函数来计算$W_{ij}$及其梯度$\nabla W_{ij}$，例如，[Muller et al. 2003]。但是，在我们的预测器 - 校正器方案中，对两者使用相同的核函数很重要，否则预测和校正步骤将无法适配。在我们的工作中，我们使用三次样条核函数[Monaghan 1992]。

在SPH模拟中，通常使用具有紧凑支持的核函数，其在有限距离处消失，也称为支持半径$h$。通常，核函数可以写为$W_h(q(x))$，其中$q=\frac{||x||}{h}$。这意味着核函数仅在$0\le q<1$时非零。因此，这种核函数的一阶空间导数$\nabla W_h(q(x))$具有相同的紧凑支持。

特别是核函数梯度的计算是模拟步骤中最耗时的任务之一，因为必须在两个求解器的每个迭代步骤中确定每个粒子的整个邻域的梯度。而且，它们需要不同的非压力并计算因子${\alpha }_i$。但是，存储所有邻域的梯度需要大量内存，不建议用于大型场景。为了加速模拟，我们提出通过使用预计算的查找表来更快地计算核函数及其梯度。使用查找表进行快速函数计算的想法并不新鲜。然而，就我们所知，它尚未在SPH模拟中使用。

由于核函数$W_h$是具有紧凑支持的标量函数，因此通过常规采样可以很容易生成查找表。但是，梯度的处理方式不同。我们不是在所有三维中对矢量函数$\nabla W_h$进行采样，而是引入标量函数$g(q)$来减少计算量和内存需求：
$$\nabla W_h(q(x))=x\cdot g(q)withg(q)=\frac{\partial W_h}{\partial q}\cdot \frac{1}{h||x||}$$

还可以定期对函数$g(q)$进行采样以获得相应的查找表。最后，通过单个查找与x相乘确定梯度。

查找表的使用是一种简单但有效的技巧，也可用于加速其他SPH方法。在第4节中，我们讨论了有关采样距离，近似误差和性能增益的详细信息。

4 Results

本节中的所有时间均在两台Intel Xeon E52697处理器（2.7 GHz，每处理器12个内核和64 GB RAM）上进行测量。我们使用OpenMP并行化流体模拟。在我们的模拟中，我们使用Ihmsen等人[2011]的并行方法进行邻域搜索，使用Akinci等人 [2012]的刚体 - 流体耦合进行边界处理。为了模拟流体的粘度，我们采用了Schechter和Bridson[2012]的XSPH变体。我们结合Becker and Teschner [2007]，Akinci等人 [2013]以及He等人 [2014]的表面张力模型成功地测试了我们的方法。然而，对于本文中的结果，我们仅使用Akinci等人 [2013]的表面张力模型。自由表面的粒子缺乏是SPH模拟中众所周知的问题，其引起粒子聚集。在模拟中，我们通过将负压强约束到零来解决这个问题，这是一种常见的解决方案，参见例如[Ihmsen等 2014a]。在所有模拟中，我们将平均密度误差控制在小于0.01%，由于密度变化率小于0.1%，密度误差也控制在0.01%以内。除非另有说明，否则我们根据CFL条件在模拟中使用自适应时间步长（参见第3.1节）。

性能 我们用125k粒子进行了溃坝模拟，以比较我们的新方法与IISPH，PBF和PCISPH的性能。粒子半径为0.02米，我们使用不同的固定时间步长。表1总结了一秒钟内模拟的性能测量结果。在表2中，显示了加速因子。注意，在[Ihmsen等人 2014a]中类似的场景用于PCISPH和IISPH之间的性能比较，并测量了可比较的结果。

通过同时满足无散条件和恒定密度条件，在模拟期间保持较小的密度误差，这显著减少了所需的求解器迭代次数。在迭代求解器中使用热启动是减少迭代次数的另一种方法。在我们的模拟中，我们使用上一个时间步长的刚度值之和初始化DFSPH求解器。在溃坝场景中，这减少了大约3倍的迭代次数。虽然DFSPH在完全热启动情况下表现最佳，但当将最后一步的解乘以0.5时，IISPH具有最佳性能[Ihmsen et al. 2014a]。由于我们的模拟中的无散速度场和求解器的热启动，与IISPH相比，我们测量的加速为6.9倍，与PCISPH相比，高达23.9倍，时间步长为4 ms。 DFSPH求解器分别仅需要4.5和2.8次迭代来校正密度误差和散度误差，而第二好的方法IISPH已经需要50.5次迭代。在我们的实验中，我们测量了较小时间步长的较小加速比，因为对于DFSPH，通常使用最小迭代次数。但是，对于具有更多粒子的场景，其中迭代次数明显高于最小值，对于小步长，加速也更大。在我们的实验中，我们注意到，当使用时间步长使迭代次数的范围在2到20次迭代之间时，我们的方法表现最佳。

如表1所示，我们的方法在比IISPH，PBF和PCISPH更大的时间步长方面具有最佳的整体性能。这具有另一个优点：当使用更大的时间步长时，不需频繁地执行计算上昂贵的邻域搜索。

在溃坝模拟中，我们使用3.4节中介绍的核函数优化测量了大约30%的性能增益。对于核函数及其梯度的采样，我们使用了1000个采样点并测量了小于$10^{-11}$的最大局部误差。由于实现非常简单且误差可以忽略不计，因此这是我们方法的一个很好的扩展。

为了测量DFSPH在大量流体粒子模拟中的性能，我们首先模拟了一个具有两个龙模型和一个移动墙的溃坝场景（见图3）。该模型由230万个流体粒子和70万个边界粒子组成。此外，我们模拟了一个真正的溃坝，其中有500万个粒子流过一个由4000万个边界粒子取样的峡谷（见图1，右图）。对于两种模拟，算法1中主要步骤的平均计算时间如表3所示。

表1：使用不同固定时间步长的DFSPH与IISPH，PBF和PCISPH的比较，用于具有125k粒子的溃坝场景（参见图2）。该表显示了所需迭代的平均数量，求解器的总计算时间以及包括模拟超过一秒的邻域搜索的总时间。最低的总计算时间标记为粗体。注意，对于DFSPH，我们使用了无散解算器和恒定密度求解器所需的时间之和，因为两个求解器的迭代步骤几乎需要相同的时间。对于DFSPH，迭代计数的列包含恒定密度求解器（cd）和无散解算器（df）的值。

表2：该表显示了基于表1中的测量值，与IISPH，PBF和PCISPH相比DFSPH的加速因子。

表3：邻域搜索，$\alpha$的计算，所有非压力$F^{adv}$的计算，恒定密度求解器和用于溃坝模拟和峡谷模拟（见图1，右）的无散解算器（参见图3）的每个模拟步骤的平均计算时间（以秒为单位）。

图3：溃坝模型（230万个流体粒子），有两个龙模型和一个移动墙。速度场采用颜色编码：蓝色为最小值，白色为最大值。

内存要求 我们的方法的一个优点是内存需求很低，特别是在模拟具有数百万个粒子的大规模场景的情况下。每个粒子我们只需要存储两个求解器使用的标量值${\alpha }_i$。执行热启动时，我们必须为每个求解器存储一个额外的标量。为了进行比较，IISPH需要七个标量值用于求解器，一个用于热启动。因此，我们的方法比IISPH需要更少的内存。

稳定性 在下文中，我们表明当前最先进的压强求解器没有显式地强制实现无散速度场，不能满足连续性方程所要求的条件$\frac{D\rho }{Dt}=0$。这可能导致模拟中的不稳定性，尤其是当粒子具有高速度时。

图2比较了在125k粒子的溃坝模拟中IISPH和DFSPH的速度散度误差。该图表明，与没有校正速度散度的方法相比，我们的方法保持了无散速度场。IISPH的最大局部散度误差为108.3，而DFSPH的最大局部散度误差仅为1.9。在另一个散度比较中，我们模拟了具有80k流体粒子的静止液柱（见图4）。由于CFL条件返回静止粒子的任意大值，因此我们将时间步长限制为最大5 ms。在IISPH模拟中高达72.9的大散度误差与大时间步长相结合导致视觉伪影：流体粒子有时会因误差而跳起（参见图4，左）。DFSPH的最大局部散度误差为2.5，这允许稳定的模拟而没有伪影。

我们用落向地面上的27k快速移动粒子的立方体进行了稳定性测试，以显示它如何影响稳定性（见图5）。根据CFL条件选择时间步长。在这个测试中，我们比较了PCISPH和DFSPH。用DFSPH进行的测试进行两次。在第一次模拟中，我们停用了无散解算器并在第二次模拟中激活它，以显示模拟在无散速度场下更稳定。在模拟中，PCISPH由于撞击而变得不稳定，没有无散解算器的DFSPH显示伪影，使用两个求解器的DFSPH保持稳定。PCISPH模拟中的不稳定性甚至导致几个流体粒子穿过边界。为了使用PCISPH稳定模拟这种情况，必须大大减小时间步长，这显著降低了整体性能。

为了证明我们的方法在动态边界模拟中的稳定性，我们在我们的模拟器中集成了Bullet物理库[Coumans 2014]。在另一个稳定性实验中，我们在具有330k粒子的溃坝中扔下了几个具有不同速度的刚体。结果如图6和随附视频所示。最后，图1（左）表明DFSPH甚至允许使用240万个快速移动的粒子进行稳定模拟。

图4：具有80k粒子的静止液柱的顶部。 IISPH模拟中的大散度误差导致跳跃伪影（左）。 DFSPH保持无散速度场，因此允许稳定的模拟而没有伪影（右）。散度误差使用颜色编码：白色是最小值，红色是最大值。

图5：模拟中的稳定性比较，其中具有27k粒子的立方体以高速度落在地面上。 PCISPH（左）由于撞击而变得不稳定，并且几个流体粒子穿过边界，没有无散度解算器的DFSPH（中间）显示伪影，有无散度解算器的DFSPH保持稳定。使用与图3中相同的颜色编码。

图6：330k粒子与动态刚体的双向耦合。

5 Conclusion and Future Work

在本文中，我们提出了一种新的隐式SPH模拟方法，用于不可压缩流体，防止了体积压缩并强制实现无散速度场。与现有技术相比，无散度配置使隐式压强求解器的收敛显著加快，从而获得显著的加速。此外，我们证明了我们的方法能够稳健地处理具有数百万个快速移动粒子的场景，并产生令人信服的物理行为。

我们的方法也有一些局限性。在SPH模拟中，低估了自由表面附近的密度，这导致不自然的粒子聚集伪影。在我们的实现中，通过将负压强限制到零来解决该问题。然而，更好的解决方案是引入幽灵粒子，如Schechter和Bridson [2012]所提出的，以防止粒子缺乏，改善流体的物理行为。此外，在没有压强限制的情况下，可以采用更复杂的求解算法，如共轭梯度法，并且可以进一步提高收敛速度。这是我们未来研究的目标。在此背景下，我们还计划研究DFSPH是否可以提高高密度对比的多相模拟的稳定性。