SPH Fluids in Computer Graphics

(1)

SPH Fluids in Computer Graphics

Markus Ihmsen¹, Jens Orthmann², Barbara Solenthaler³, Andreas Kolb², Matthias Teschner¹

1University of Freiburg

2University of Siegen

3ETH Zurich

Abstract

Smoothed Particle Hydrodynamics (SPH) has been established as one of the major concepts for fluid animation in computer graphics. While SPH initially gained popularity for interactive free-surface scenarios, it has emerged to be a fully fledged technique for state-of-the-art fluid animation with versatile effects. Nowadays, complex scenes with millions of sampling points, one- and two-way coupled rigid and elastic solids, multiple phases and additional features such as foam or air bubbles can be computed at reasonable expense. This state-of-the-art report summarizes SPH research within the graphics community.

Keywords:Physically-based animation, fluid animation, Smoothed Particle Hydrodynamics

Categories and Subject Descriptors (according to ACM CCS): Computer Graphics [I.3.7]: Three-Dimensional Graphics and Realism—Animation

1. Introduction

This section starts with a compact overview of SPH fluids.

Sec.1.1introduces the underlying equations, in particular the momentum equation. Sec.1.2explains how to use SPH for the interpolation of fluid quantities and for the approximation of spatial derivatives in the momentum equation. Fi- nally, Sec.1.3presents a first simple SPH fluid solver and introduces its components.

The remainder of this state-of-the-art report discusses various solver components in detail. Sec.2explains approaches to estimate the neighborhood of a particle. Although recent research in that direction focuses on uniform grids, various realizations based on different concepts and considering different architectures such as GPUs are presented. Sec.3dis- cusses various ways to compute pressure forces, either using a state equation (EOS), EOS-based iterative refinement, or pressure projection. Benefits, drawbacks, and issues with performance comparisons are discussed. Sec 4 discusses methods for boundary handling with a focus on boundaries that are represented with particles. Sec.5discusses variants that employ adaptive time or space discretization to opti- mize performance, while Sec.6 discusses multiphase fluids. Sec.7discusses various techniques to reconstruct interface representations from particles and specific rendering approaches. Sec.8presents approaches to efficiently add

small-scale detail, e.g., foam, spray, and tiny air bubbles, to pre-computed SPH simulations in order to enhance the visual quality. Finally, Sec.9presents a few aspects that could be addressed in future research.

1.1. Governing Equations

We consider a fluid that consists of a set of small moving fluid elements, i.e., particles. Each particleihas a massmi

and carries attributes such as densityρi, pressurepior vol- umeVi. Over timet, particle positionsxiand the respective attributes are advected with the local fluid velocityv_i:

dx_i

dt =vi. (1)

As the particles move with the fluid flow, the time rate of change of the velocity v_iis governed by the Lagrange form of the Navier-Stokes equation (see AppendixAfor a comparison of Lagrange and Euler formulations):

dv_i dt =−1

ρi

∇p_i+ν∇²v_i+F^other_i mi

. (2)

The term−¹

ρi∇p_idescribes the particle acceleration due to pressure differences in the fluid. The respective pressure

c

The Eurographics Association 2014.

(2)

force generally dominates all forces. It is responsible to preserve the volume of the fluid. As the particle mass is constant, preserving the fluid volume corresponds to preserving the density. Small and preferably constant density deviations are important for high-quality simulations. Otherwise, a per- ceivable and disturbing bouncing of the free surface occurs.

The termν∇²v_irepresents the acceleration due to friction forces between particles with different velocities. Although the kinematic viscosityνis known, e.g.,ν≈10⁻⁶m²·s⁻¹ for water, larger user-defined values are typically preferred to improve the stability of SPH simulations. Viscosity is also realized with alternative methods such as artificial viscosity [Mon92] or XSPH [Mon89]. The term ^F

other i

mi describes other accelerations such as gravity.

1.2. SPH

The SPH concept is used to interpolate fluid quantities at arbitrary positions and to approximate the spatial derivatives in Eq. (2) with a finite number of sample positions, i.e., adjacent particles.

Interpolation: A quantityAiat an arbitrary position x_iis approximately computed with a set of known quantitiesAj

at neighboring particle positionsx_j: Ai=

∑

j

mj

ρj

AjWi j (3)

withWi jbeing a kernel function of the form Wi j=W

kx_i−x_jk h

=W(q) = 1

h^d f(q) (4) wheredindicates the number of dimensions andhis the so- called smoothing length. Kernel functions should be close to a Gaussian [Mon92], but with a compact support that typically ranges fromhfor the bell-shaped function [Luc77]

to 3hfor the quintic spline function [Mor96]. The number of adjacent particles that are considered in the summation of Eq. (3) depends on the dimensionalityd, the support of the kernel function and the particle spacing which is typically close toh[Mor96,Mon05]. The choice of the kernel function, the number and the disorder of considered particles influence the accuracy of the summation in Eq. (3). A typical function for the kernel shown in Eq. (4) in 3D would be the cubic spline [Mon92]:

f(q) = 3 2π







2

3−q²+¹₂q³ 0≤q<1

1

6(2−q)³ 1≤q<2

0 q≥2

. (5)

There is, however, no consensus on the optimal kernel with respect to the trade-off between accuracy and computa- tional cost. E.g., [APKG07,OK12,OHB^∗13] use a set of ker- nels proposed in [DC96,MCG03,MSKG05]. [ZYF10] uses a quintic spline [Mor96]. [BT07,AIS^∗12,SB12] use the cubic spline (Eq. (5)).

Spatial derivatives: Spatial derivatives can be computed in various ways. In order to address issues of the original formulations ∇Ai=∑j

mj

ρjAj∇Wi j and

∇²Ai=∑j mj

ρjAj∇²Wi j, various alternatives have been in- vestigated and currently, the following approximations are preferred [Mon92,MFZ97]:

∇A_i = ρi

∑

j

m_j Ai

ρ²_i +Aj

ρ²_j

!

∇W_{i j}, (6)

∇ ·A_i = −1 ρi

∑

j

mjA_{i j}· ∇Wi j, (7)

∇²Ai = 2

∑

j

mj

ρj

Ai j

xi j· ∇Wi j

xi j·xi j+0.01h², (8) with A_{i j} = A_i−A_j, A_{i j} = A_i−A_j, x_{i j} = x_i−x_j and

∇Wi j= _∂W

i j

∂xi,x,^∂W_∂x^{i j}

i,y,^∂W_∂x^{i j}

i,z

T

. Eq. (6) and Eq. (8) can be used in the computation of particle accelerations in Eq. (2).

Eq. (7) can, e.g., be used to predict density changes from the divergence of the velocity field based on the continuity equation [ICS^∗13]. While pressure forces are preferably computed with the formulation in Eq. (6), e.g. [BT07,SP09b, RWT11], there exist various alternative forms, e.g. [MCG03, APKG07,LD09].

1.3. Concept of an SPH-based Fluid Solver

The basic building blocks of SPH-based fluid solvers are:

neighborhood search, pressure computation and time integration. The neighborhood search is typically accelerated by a spatial access structure, e.g. a uniform grid [THM^∗03, LD08,GSSP10,IABT11], with a cell size that is equal to the kernel support, e.g., 2hfor the kernel in Eq. (5). Details are discussed in Sec.2.

In order to compute the pressure gradient in Eq. (2), pressure piis computed from the densityρi. While Sec.3dis- cusses a variety of alternatives, state equations are a simple and popular choice and, in this context, the formulation

pi=k ρi

ρ0

7

−1

!

(9) seems to be preferred [Mon94,BT07,LD09,YT10,RWT11, Mon12,YT13]. The valueρ0 is the desired rest density of the fluid, kis a stiffness constant that scales the pressure and, thus, the pressure gradient and the respective pressure forces. In practice, a larger stiffness constant reduces the compressibility of the fluid, but demands smaller integration time steps.

Alg.1 shows an example of a simple SPH simulation step. As the algorithm employs a state equation, it can be referred to as state equation SPH (SESPH) [ICS^∗13]. In an implementation, the kernel functionWi j has to be specified, e.g., cubic with a support of 2h(see Eq. (4) and Eq. (5)).

The particle massmihas to be specified, e.g.,mi=h³ρ0as

c

(3)

Algorithm 1SPH with state equation.

for allparticle ido find neighbors j for allparticle ido ρi=∑jmjWi j

computepiusingρi(e.g. Eq. (9)) for allparticle ido

F_i^pressure=−^mⁱ

ρi∇pi(e.g. Eq. (6)) F^viscosity_i =miν∇²vi(e.g. Eq. (8)) F^other_i =mig

F_i(t) =F^pressure_i +F^viscosity_i +F^other_i for allparticle ido

vi(t+∆t) =vi(t) +∆tFi(t)/mi

x_i(t+∆t) =x_i(t) +∆tv_i(t+∆t)

in [ICS^∗13] ormi= (²₃h)³ρ₀as in [SB12] resulting in different numbers of neighbors that contribute to the SPH sums of a particle.

In terms of numerical integration schemes, the semi- implicit Euler, also referred to as symplectic Euler or Euler- Cromer, is regularly used, e.g. [IAAT12,SB12,ICS^∗13].

The time step∆tis governed by the Courant-Friedrich-Levy (CFL) condition, e.g.,∆t≤λ_kvmax^h kwithλ≈0.4 [Mon92],h being the particle diameter andv^maxbeing the maximum velocity of all particles. Forλ=1, this constraint states that all particles move less than the particle diameter per time step.

A detailed discussion of further aspects that can be incorpo- rated into the time step analysis can be found in [IAGT10].

Boundary handling forces as discussed in Sec.4are not considered in Alg.1.

2. Neighborhood Search

SPH requires the computation of sums over dynamically changing sets of neighboring particles. The search of these neighborhood sets is generally accelerated by spatial data structures that should be efficiently generated and queried, preferably in a parallelizable way. While the neighborhood search has similarities with, e.g., collision detection or inter- section tests in raytracing, it is additionally characterized by the fact that more than one space cell has to be queried to find the neighbors of a particle, i.e., cells adjacent to the cell of a particle have to be accessed. The data structure should also be efficient for arbitrary, sparsely filled simulation domains with a non-uniform particle distribution.

Although hierarchical data structures such as kd-trees are used in multi-resolution scenarios with a variable kernel support [KAD^∗06,Kei06,APKG07,SBH09], there seems to be consensus to employ uniform grids for standard SPH with a fixed kernel support. Grids are built inO(n)and particles are accessed inO(1), while hierarchical data structures are typi-

cally built inO(nlogn)and accessed inO(logn)[HKK07b].

Therefore, this section focuses on uniform grids.

In addition to the implementation of the grid, there are further choices that influence the efficiency of an SPH solver. E.g, it is not obvious whether to store the neighborhood set for reuse as, e.g., in [IABT11], or not as, e.g., in [HKK07b,ZSP08,GSSP10]. This mainly depends on the number of neighborhood queries per simulation step. While this number is low for non-iterative SPH solvers (see Sec.3.1 and Sec.3.2), it can be large for iterative solvers (see Sec.3.3 and Sec.3.4). Another performance issue is the frequent re- computation of neighborhood sets. This issue can be allevi- ated by Verlet lists [Ver67,Hie07,PH10], where a set of potential neighbors is computed within a distance that is larger than the actual kernel support. Actual neighbors are computed from the set of potential neighbors which is updated only every n-th simulation step depending on the ratio between kernel support and query distance. However, Verlet lists are prohibitively memory-intensive for complex scenarios and slower than the strategies presented in the following.

The following discussion is limited to various implementations of uniform grids in the context of SPH applications with uniform kernel support.

2.1. Uniform Grid

Space is subdivided into cubic cells and each particle is associated to one cell in the construction stage. To find all relevant neighbors, the cell of a particle and all adjacent cells are queried. If the cell size is equal to the kernel support, 3³cells have to be queried in 3D, which is optimal according to [IABT11]. The number of particles associated with a cell and the number of neighboring particles depend on the initial particle distance, e.g., approximately 8 and 40, re- spectively, for a kernel support twice as large as the initial particle distance. While the query step is easy to parallelize, grid construction is not due to potential write conflicts (race conditions).

2.1.1. Index Sort

In order to avoid race conditions in the parallel construction of uniform grids, particles can be sorted with respect to a key that is unambiguously assigned to each cell [PDC^∗03, Gre08,KS09]. Instead of storing references to all particles, a cell only stores one reference to it’s first particle in the sorted array. Parallel reduction is commonly employed to compute these references [KS09]. By sorting particles according to their spatial cell, particles in the same cell are close in memory which improves memory coherence. However, particles in neighboring cells are not necessarily close in memory.

2.1.2. Z-index Sort

Efficient algorithms have to enforce low memory transfers such that threads can perform the operations with almost

(4)

no latency. The transfer rate decreases for higher cache-hit rates, i.e., the percentage of accesses where the requested data is already present in the cache. The cache-hit rate of any SPH implementation can be optimized by mapping the spatial locality of particles onto memory. This can be achieved by employing a space-filling Z-curve for computing cell indices [GSSP10,IABT11].

Rather than sorting an n-dimensional space one di- mension after another, a Z-curve orders the space by n-dimensional blocks of 2ⁿ cells. This ordering pre- serves spatial locality due to the self-containing (recur- sive) block structure. Consequently, it leads to a high cache-hit rate while indices can be computed fast by bit- interleaving [PF01]. As shown in [IABT11], the Z-curve increases the cache-hit rate and, thus, improves the performance for the query and processing of particle neighbors.

Index sort variants are considered to be one of the fastest spatial acceleration methods. However, in these schemes, the memory consumption scales with the simulation domain [Gre08]. In order to represent infinite domains with low memory consumptions, spatial hashing can be employed.

2.1.3. Hashing

In spatial hashing [THM^∗03], the effectively infinite domain is mapped to a finite list. The hash function that maps a po- sitionx= (x,y,z)to a hash table of sizemhas the following form:

c= hjx

d k

·p1

xor

jy d k

·p2

xor

jz d k

·p3

i

%m, (10) withp₁,p₂,p₃ being large prime numbers that are chosen as 73856093, 19349663 and 83492791 [THM^∗03], respec- tively. Different spatial cells can be mapped to the same hash cell (hash collision), slowing down the neighborhood query. The number of hash collisions can be reduced by increasing the size of the hash table, trading memory for speed. [THM^∗03] suggests to reserve memory for a certain numberkof entries in all hash cells on initialization in order to avoid frequent memory allocations. However, for SPH fluids, the hash table is generally sparsely filled. Thus, a significant amount of memory is unnecessarily pre-allocated.

Furthermore, cells that are close in memory are necessarily not close in space which reduces the cache-hit rate for the neighborhood query. These issues are addressed by the compact hashing method.

2.1.4. Compact Hashing

Compact hashing [IABT11] uses a secondary data structure which stores a compact list of non-empty (used) cells. Hash cells just store a handle to their used cell. Memory for a used cell is allocated if it contains particles and deallocated if the cell gets empty. Thus, constant memory is consumed for the hash table and additional memory for the list of used cells.

Thereby, the memory consumption scales with the number of particles and not with the simulation domain.

Generally, the hash function is designed to abolish spatial locality in order to minimize the number of hash collisions. As this would result in an increased memory transfer and longer query times compared to index sort methods, [IABT11] proposes to reorder the particles and the compact list of used cells according to a Z-curve. As particles show a high temporal locality, sorting might not be required in each simulation step. Depending on the performance of the sorting algorithm and the number of particles that move between cells, reordering can be performed in each [DRF12]

to every 100th [IABT11] simulation step in order to yield the optimal performance. Interestingly, when reordering is performed in each step, the compact list just needs to store a reference to the first particle in the sorted array and the number of entries, instead of references to all particles in the cell. This further reduces the memory consumption of the compact list.

2.1.5. GPUs

Since there are almost no data dependencies, SPH methods generally map well to the streaming architecture of today’s Graphic Processing Units (GPUs) [HKK07b,Gre08,ZSP08, YWH^∗09,GSSP10,OK12,MM13]. On GPUs, neighborhood computations can be carried out via two dual mapping operations: either by scattering particle contributions onto points of evaluation, or by gathering of particle data at sampling positions.

Scattering operations are efficiently realized in combination with the rasterization pipeline of programmable graphics hardware [KLRS04,KSW04,AIY^∗04,KC05,HCM06, HKK07a,HKK07b,ZSP08]. Field functions are mapped to fragment shaders which blend particle contributions into several texture slices using render-to-texture mechanisms [KC05]. The absence of acceleration structures and neighborhood queries makes scattering very attractive for interactive rendering of SPH data [vdLGS09,FGE10,FAW10].

In contrast, gathering operations better exploit paral- lelism in combination with generic programming APIs such as CUDA or OpenCL, which is advantageous for simulation [Gre08,GSSP10,OK12,MM13] and surface reconstruc- tion [AIAT12,AAIT12]. However, a vital criterion for the efficiency of gathering approaches is thread coherence including coherent memory access between concurrent threads.

Thus, neighborhood search is in general realized with an index sort [Gre08] or Z-index sort [GSSP10,ZGHG10], employing a data-parallel radix sort [SHG09] and stream com- paction [SHZO07] as main building blocks. Furthermore, it is advantageous to avoid data transfer between CPU and GPU, favoring fully GPU-based systems which due to memory limits trade particle numbers for computation speed.

However, in combination with domain decomposition techniques, multi-GPU solvers [VBDRC12,ZSL^∗13,RBH^∗13]

can be utilized in order to increase the particle resolution.

c

(5)

3. Incompressibility

Enforcing incompressibility is essential for realistic SPH fluid simulations. While oscillations of the free surface due to compressibility are less prominent in small scenarios, oscillations become more significant for larger scenes with, e.g., millions of particles [ICS^∗13].

The computation of particle states with low compressibility, usually between 0.1% and 1% [MM97,SP09b,ICS^∗13], is one of the most expensive steps in SPH. Either the computation per time step is efficient, albeit at the expense of a small time step, or the computation is expensive for a larger time step.

This survey classifies all discussed incompressibility approaches into four classes. First, non-iterative approaches are discussed that employ an equation of state (EOS). Sec- ond, EOS solvers with splitting are outlined. Third, iterative EOS solvers are explained. Fourth, solvers based on a pressure Poisson equation are outlined.

3.1. Non-iterative EOS Solvers

State equations are used in the context of Alg.1to compute pressure from density. For example,pi=c²sρihas been proposed in [MM97]. Here,csis often referred to as the speed of sound. This is motivated by the fact that fluids with a Mach number greater than 0.3 are of considerable compressibility.

As a given constantcsresults in a certain compressibility, the respective simulation does not represent a real fluid, but an artificial one with a reduced sound speed [Mon94]. In practice,cs is a stiffness constant that scales the pressure and, thus, the pressure gradient and the respective pressure forces.

A larger value reduces the compressibility of the fluid, but also limits the time step. Therefore, the stiffness constants for different variants of the state equation are all denoted withkin the following.

Alternatively to [MM97], pressure is often related to ρi−ρ0, the deviation of the actual density to the rest density, i.e., the density error. Here, slightly varying forms are employed, e.g., pi =k(ρi−ρ0) [DC96,MCG03] or pi=k(ρ_i/ρ₀−1)[APKG07]. However, as already discussed in Sec.1, the formulation pi=k((ρi/ρ0)⁷−1)is widely used [Mon94,BT07,LD09,YT10,RWT11,YT13]. [ZYF10]

proposespi=k((ρi/ρ0)²−1).

Differences in terms of stability and performance between the various EOS formulations are rarely analyzed in the literature. Generally, the incorporation of any of the discussed state equations seems to be less efficient compared to iterative EOS solvers (Sec.3.3) and pressure projection (Sec.3.4). Nevertheless, it would be interesting to see the performance of a state equation in combination with a non- iterative splitting approach as given in Alg.2.

3.2. Non-iterative EOS Solvers with Splitting

As an interesting alternative to Alg.1, the pressure could be computed with the density that is obtained after advecting the particles without pressure forces. This concept is known as splitting, e.g. [Cho68,Bri08], and the basis of various iterative solvers, e.g. [PTB^∗03,SP09b,SBH09,HLWW12, ICS^∗13,MM13]. Therefore, accelerations in Eq. (2) are split, i.e., considered in two different steps. First, an intermediate velocity v^∗ is computed from all non-pressure accelerations: (v^∗−v(t))/∆t=ν∇²vi(t) +F^other_i (t)/mi. Then, the acceleration from the pressure gradient is computed for the resulting velocity:(v(t+∆t)−v^∗)/∆t=−(1/ρ^∗_i)∇pi. Pressurepiis computed usingρ^∗_i, whileρ^∗_i is computed after advecting the particles withv^∗_i. The respective pressure forces are intended to project the intermediate velocityv^∗ onto a divergence-free velocity field, i.e., to minimize density deviations. Alg.2shows a simple implementation of this concept. While non-pressure forces compete with pressure forces in Alg.1, Alg.2considers the effect of non-pressure forces in the computation of pressure forces. Alg. 1 and Alg.2have not been compared in the literature yet. How- ever, since iterative SPH solvers based on splitting generally outperform Alg.1, e.g. [SP09b,HLWW12], simple non- iterative splitting as in Alg.2is certainly a promising concept.

Algorithm 2SPH with state equation and splitting.

for allparticle ido find neighbors j for allparticle ido

F^viscosiy_i =miν∇²v_i(e.g. Eq. (8)) F^other_i =m_ig

v^∗_i =v_i(t) +∆t^F

viscosity

i +F^other_i

mi

for allparticle ido

ρ^∗_i =∑jmjWi j+∆t∑j(v^∗_i −v^∗_j)· ∇W_{i j} computepiusingρ^∗_i (e.g. Eq. (9)) for allparticle ido

F_i^pressure=−_ρ^m∗ⁱ

i ∇pi(e.g. Eq. (6)) for allparticle ido

v_i(t+∆t) =v^∗_i +∆tF_i^pressure/m_i xi(t+∆t) =xi(t) +∆tvi(t+∆t)

3.3. Iterative EOS Solvers with Splitting

This concept is also based on the splitting concept illustrated in Alg.2. Intermediate velocities and positions are predicted using all non-pressure forces. To minimize density errors at the intermediate state, pressure forces are calculated. As an extension to Alg.2, these pressure forces are iteratively refined. In each iteration, the pressure forces lead to updated intermediate positions and velocities with a new density error that is considered in the force computation of the following iteration. If the density errorρerr- either the average

(6)

or the maximum - is below a thresholdη, the algorithm ter- minates. In contrast to non-iterative EOS solvers in Sec.3.1 and Sec.3.2, the stiffness constant in the state equation is generally computed and not user-defined. Instead, iterative EOS solvers are parameterized by a more appropriate and intuitive density errorη. Smaller values forηresult in more iterations, larger values are more efficient as less iterations are required.

Alg. 3 illustrates the concept. The notation is closely aligned with the iterative local Poisson approach (LPSPH) in [HLWW12], but can be mapped to other iterative EOS solvers, e.g. [SP09b] as well. While Alg.3is not an optimal implementation, it is a good illustration of the concept.

Algorithm 3Iterative EOS solver with splitting.

F^viscosity_i =m_iν∇²v_i(e.g. Eq. (8)) F^other_i =mig

v^∗_i =v_i(t) +∆t^F

viscosity

i +F^other_i

mi

x^∗_i =x_i(t) +∆tv^∗_i repeat

for allparticle ido computeρ^∗_i usingx^∗_i

computepiusingρ^∗_i, e.g.pi=k(ρ^∗_i −ρ0) computeρerr

for allparticle ido F^pressure_i =−^mⁱ

ρ^∗i ∇p_i(e.g. Eq. (6)) v^∗_i =v^∗_i +∆t^F

pressure i

mi

x^∗_i =x^∗_i +∆t²^F

pressure i

mi

untilρerr<η for allparticle ido

v_i(t+∆t) =v^∗_i xi(t+∆t) =x^∗_i

In addition to [HLWW12], other approaches implement this concept to enforce incompressibility with interesting variations. In predictive-corrective SPH (PCISPH) [SP09b], pressure is accumulated instead of pressure forces.

[HLWW12] recomputes particle neighbors in each iteration, [SP09b] does not. Also, different stiffness constants are used. PCISPH uses

k= 2ρ²₀

m²_i·∆t²∑j∇W_{i j}⁰·∑j∇W_{i j}⁰+∑j(∇W_{i j}⁰· ∇W_{i j}⁰) (11) withW_{i j}⁰being kernel values for a prototype particle. So,k can be mainly precomputed and only changes with the time step. In contrast, LPSPH uses

k= ρ^∗_ir²_i

2ρ0∆t² (12)

with initial particle distance 2ri which is derived from the pressure Poisson equation [HLWW12].

In each iteration, PCISPH and LPSPH compute intermediate positions. This is also done in position-based fluids (PBF) [MM13] and to some extent in [CBP05]. Following the splitting concept, PBF uses non-pressure forces to predict intermediate positions and velocities. Then, positions are iteratively refined to enforce incompressibility. There- fore, the state equationp_i=k(ρ_i/ρ₀−1)withk=1 is employed. Instead of accumulating pressure or pressure forces, PBF accumulates distances in each iteration with

∆x_i=−1 ρ0

∑

j

pi

βi

+pj

βj

∇W_{i j}, (13) where ∆x_i is the position change per iteration.βi and βj

are precomputed constants. The position update in Eq. (13) is closely related to an SPH pressure force. However, PBF [MM13] avoids accumulating pressure or pressure forces that eventually update the velocity and the position.

Another closely related approach is presented in [BLS12].

Splitting is employed and then, velocities and positions are iteratively updated. The same state equation as in [MM13]

is employed. In contrast to PBF [MM13] and similar to PCISPH [SP09b], pressure is accumulated. From an SPH perspective, [BLS12] is difficult to read. Density is computed with SPH using the regular kernel functionWi j. The kernel gradient, however, is denoted asfi jrˆi j. The approach solves for Lagrange multipliersλithat are related to the neg- ative pressure as can be seen, e.g., from the formulation of the constraint force:f_i=∑j(m_iλi+m_jλj)f_{i j}ˆr_{i j}.

Performance: The performance of iterative EOS solvers is commonly characterized by the maximum possible time step and the required number of iterations for a specified density error. [SP09b] shows that PCISPH allows for time steps that are up to two orders of magnitude larger than in a non-iterative EOS solver [BT07]. The average number of iterations is between three and five for density errors of 0.1%

and 1%, resulting in an overall speedup factor of fifty compared to [BT07]. [HLWW12] presents a slightly improved performance with a speedup of 1.5 compared to PCISPH.

PBF tolerates significantly larger time steps than PCISPH, but requires more iterations, resulting in a similar overall performance. In [BLS12], speedup factors of 25 to 250 are reported compared to non-iterative EOS solvers. Typically, 5 to 15 iterations are employed, while a density error of 1% is obtained with 100 iterations. However, as discussed in [ICS^∗13,MM13] a thorough performance analysis of iterative EOS solvers is rather involved. This is due to the fact that these solvers do not necessarily reach their optimal performance for the largest time step. On one hand, the number of iterations grows with the time step. On the other hand, the neighborhood search has to be performed per simulation step.

Discussion: Non-iterative EOS solvers use different state

c

(7)

equations, in particular they use different stiffness constants.

They also differ in terms of the quantity that is accumulated. Pressure, pressure forces or distances are accumulated to compute final particle positions. It would be interesting to analyze these aspects: Which stiffness constant is optimal?

Which quantity should be updated? Is it an option to analyze individual stiffness constants per particle?

3.4. Pressure Projection

As an alternative to state equations, pressure can also be computed by solving a pressure Poisson equation (PPE).

Following the splitting concept, intermediate velocitiesv^∗_i are predicted by applying all non-pressure forces. Then, pressure pi is computed from a discretized PPE in order to correct (project) the intermediate velocities to a divergence-free state vi(t+∆t) =v^∗_i −(1/ρ^∗_i)∇pi. This technique is commonly applied in grid-based approaches, e.g. [FF01,Bri08,CM11]. In SPH approaches, the pressure Poisson equation is used with different source terms, either using the divergence of the intermediate velocity fieldv^∗_i, e.g. [CR99,PTB^∗03], as

∇²pi= ρ0

∆t∇ ·v^∗_i (14) or the compressionρ0−ρ^∗_i after advecting the particles with v^∗_i, e.g. [SL03,KGS09], as

∇²pi= ρ0−ρ^∗_i

∆t² . (15)

Some authors propose combinations of both source terms [HA07,LTKF08]. The density can also be replaced by the number densityδi=∑jWi j [PTB^∗03,HA07,LTKF08].

The concept of pressure projection, referred to as incompressible SPH (ISPH), is illustrated in Alg.4. It is very similar to Alg.2, but instead of computing pressure per particle with the state equation, pressure is computed by solving the linear system that is described with the PPE, Eq. (14) or Eq. (15).

Recently, an alternative ISPH solver has been presented by [ICS^∗13]. The approach is referred to as implicit incompressible SPH (IISPH). Although solving a linear system seems to be expensive at first sight, IISPH outper- forms PCISPH [SP09b] and a standard ISPH variant [SL03].

IISPH employs the density invariance condition as source term in Eq. (15). Further, IISPH combines a discretized form of the continuity equation and an SPH form of the pressure force to a discretized form of the PPE. The employed form of the pressure force is equal to the pressure force that is used in the velocity update. This concept is different to approaches that either directly discretize the Laplace operator, e.g. [CR99,SL03,HA07,KGS09], or to approaches that use a background grid for the discretization [LTKF08,YLHB09,RWT11].

PPEs can be solved in various ways, e.g., using successive

Algorithm 4SPH with pressure projection (Incompressible SPH).

F^viscosiy_i =miν∇²v_i(e.g. Eq. (8)) F^other_i =mig

v^∗_i =v_i(t) +∆t^F

viscosity

i +F^other_i

mi

for allparticle ido

ρ^∗_i =∑jmjWi j+∆t∑j(v^∗_i −v^∗_j)· ∇Wi j

solve the PPE in Eq. (14) or Eq. (15) for allparticle ido

F_i^pressure=−^mⁱ

ρ^∗i ∇p_i(e.g. Eq. (6)) for allparticle ido

vi(t+∆t) =v^∗_i +∆tF_i^pressure/mi

x_i(t+∆t) =x_i(t) +∆tv_i(t+∆t)

Figure 1:A breaking dam simulated with 20 million SPH particles. The volume is preserved up to an error of 0.1%.

Velocities are color-coded.

over-relaxation (SOR) [FM96], conjugate gradient [CR99, FF01] or multigrid techniques [CM11]. In IISPH, relaxed Jacobi is employed to iteratively compute the pressure field.

[ICS^∗13] shows that the solver can be implemented in a very efficient, matrix-free way. Only seven scalar values are stored per particle and only two particle loops are required per iteration. Further, comparatively few iterations are sufficient to obtain small density deviations of down to 0.01%.

Compared to PCISPH, a speedup of up to six is presented for the computation of the pressure field.

Discussion: Although IISPH needs to solve a linear system, the actual implementation corresponds to iterative EOS solvers that accumulate pressure as in PCISPH. Basically, pressure is iteratively updated. In each iteration, PCISPH

(8)

scene particles neighborhood pressure IISPH iterations ∆t spacing

breaking dam (Fig.1) 20M 8.5s 21.5s 10 0.0034s 0.1m

fountain (Fig.2) 17M 7.0s 10.0s 4 0.0008s 0.0125m

ships (Fig.3) 19M 8.0s 24.0s 12 0.017s 0.5m

Table 1:Measurements given for a single simulation step performed on a standard six-core desktop computer.

Figure 2:A fountain simulated with 17 million SPH particles. Velocities are color-coded.

loops three times over all particles, IISPH requires only two loops. The memory footprint of seven scalar values per particle in IISPH is negligible. In terms of the time step, [ICS^∗13]

shows that the optimal time step does not necessarily correspond to the maximum possible time step in IISPH and PCISPH. [ICS^∗13] also illustrates the relation between time step and particle size. Time steps of up to 0.05s are presented for particles with a radius of 1m.

Figures 1, 2, 3 show example scenarios that have been computed with compact hashing [IABT11] and IISPH [ICS^∗13] on a standard six-core computer. Perfor- mance measurements are summarized in Table1.

3.5. Performance Comparison

Existing publications indicate that iterative EOS solvers, e.g., PCISPH and PBF, are more efficient than non-iterative EOS solvers and that pressure projection, e.g., IISPH, is more efficient compared to iterative EOS solvers. As can be deduced from the detailed performance analysis in [ICS^∗13], a thorough comparison is rather complex. This is due to various reasons.

Criteria: There seems to be an agreement to measure the overall computation time for a scenario to account for the different characteristics of existing solvers. Non-iterative

Figure 3:Three ships sailing at a speed of 60 kilometers per hour. The 19 million SPH particles are color-coded according to velocity.

EOS solvers are fast per simulation step with rather small time steps, while iterative EOS solvers and pressure projection schemes are expensive per simulation step, but allow for large time steps. As the overall computation time of all solvers largely depends on the obtained incompressibility, average or maximum density errors are considered to spec- ify the simulation quality.

Parameters: In non-iterative EOS solvers, the performance depends on one or more stiffness constants in the EOS. These stiffness constants govern the obtained density error, i.e. the simulation quality. Larger constants require smaller time steps, while the computation time per simulation step is constant. The desired density error can not be specified, but only tested. This is in contrast to iterative EOS solvers or pressure projection which are parameterized by a desired density error. Here, the computation time per simulation step grows for smaller specified density errors. On the other hand, larger time steps can be used. A performance comparison of non-iterative and iterative EOS solvers is only useful, if the density error in both approaches is compara- ble. Thus, the stiffness constant in one approach has to be mapped to the specified density error in the other approach.

Optimal Performance: As discussed in [ICS^∗13], SPH simulations with iterative EOS solvers or pressure projection schemes do not reach their maximum overall performance for the largest possible time step. I.e., there exists an optimal time step, where the combination of neighbor search and pressure computation reaches its best performance. While

c

(9)

this optimal time step is mentioned in [ICS^∗13], it has not been analyzed and it is not explained how to find it without testing.

Test scenarios:Depending on the test scenario, arbitrary performance differences can be obtained. E.g., for one SPH particle, a non-iterative EOS solver is the fastest solution, as iterative EOS solvers and pressure projection commonly perform a minimum number of iterations. Non-iterative EOS solvers can even be optimal for arbitrarily large scenes with millions of particles, if these particles form rather shallow water. Iterative EOS solvers and pressure projection, on the other hand, are more efficient for complex scenes. The term

"complex" is rarely defined, but usually refers to large scenarios with a certain fluid depth, e.g. a breaking dam. It is more expensive to obtain a desired density error for growing fluid depths. Additionally, with growing fluid depth the tol- erable density error decreases to avoid visible oscillations of the free surface, making the pressure computation even more expensive.

4. Boundary Handling

The interaction of the fluid with rigid boundaries requires special consideration in order to prevent penetration of objects. In most SPH implementations, solids are sampled with particles which exert forces on fluid particles. [Mon94, Mon05,MK09] compute distance-based penalty forces, e.g., Lennard-Jones forces which scale polynomially with the distance to the fluid particle. Although this approach has been adapted to realize the interaction of compressible SPH fluids with particle-sampled deformable meshes [MST^∗04, LAD08], it causes large pressure variations in the fluid which further restricts the time step for weakly compressible fluids. The main difficulty in penalty-based approaches is controlling the stiffness parameter which has to be bal- anced such that penetrations of rigid boundaries are avoided, while not causing pressures that are too high. Generally, these methods require small integration time steps to pro- duce smooth pressure distributions.

In order to overcome the issues of penalty-based methods and to have more control on the boundary condition, direct forcing has been proposed in [BTT09]. In this method, one- and two-way coupling of rigid bodies and fluids are realized by computing control forces and velocities using a predictor- corrector scheme. Thereby, different slip-conditions can be modeled, while non-penetration is guaranteed. Compared to penalty-based methods larger time steps can be used.

A phenomenon that occurs in distance-based penalty schemes and in direct forcing is sticking of fluid particles to the solid boundary. This results fromparticle deficiency at the interface with the solid boundary. Here, the support domain is not sufficiently sampled and, thus, field variables cannot be well approximated with the SPH interpolation concept. Harada et al. addressed this problem for the

distance-based approach [HKK07b] by employing a wall- weight function which based on the distance adds a precomputed contribution of the boundary to the fluid density.

This significantly reduces the stacking of particles, but introduces irregular density distributions at the boundary as shown in [IAGT10].

In order to obtain smoother transitions at the solid interface, boundary particles should contribute to the recon- struction of field variables. Therefore, solid objects are either pre-sampled with particles, e.g., [KAD^∗06,SSP07, FG07,IAGT10,SB12,AIS^∗12], or sampled on the fly, e.g., [HA06]. Furthermore, different strategies have been proposed how boundary particles are taken into account. For example in [MM97,SSP07,IAGT10], boundary particles are treated like fluid particles which are advected by the velocity of the solid, i.e., for each boundary particle a unique density and pressure value is computed. Another strategy is to mirror quantities of neighboring fluid particles onto boundary particles, e.g., a boundary particle gets the same pressure as its neighboring fluid particle. Mirroring quantities has been suc- cessfully employed to simulate different slip conditions for straight [HA06] and curved [MM97,SB12,AIS^∗12] boundaries.

The sampling distance of boundaries influences the numerical stability and the quality of the simulation significantly. Simple objects like a cube can be equidistantly sampled with particles, but for complex shapes with convex and concave regions, irregular samplings cannot be pre- vented. Furthermore, sampling only the surface of solid objects with a single layer is desirable for performance reasons and for thin objects. These challenges have been addressed in [AIS^∗12], where inhomogeneous samplings are handled by dynamically computing the relative contribution of a boundary particle. This yields smooth reconstructions for densely sampled, complex boundaries, including dynamically changing contacts with multiple objects as illustrated in Fig.4. Non-penetration is guaranteed even when sampling objects with a single particle layer, enabling plausible interactions with lower dimensional objects. Two-way coupling has been demonstrated even for large density ratios of up to 1000. As the approach purely builds on fluid quantities such as pressure or viscosity to handle collisions and friction, it can be easily integrated into any SPH implementation, e.g., IISPH as presented in [ICS^∗13]. An extension of [AIS^∗12]

for realizing one- and two-way coupling with deformable solids and cloth has been proposed in [ACAT13].

As shown in [OHB^∗13], the two-way coupling concept can be employed to model transport, i.e., the exchange of quantities between rigid objects and fluids. Thereby, convection is extended by a robust surface transport.

5. Adaptivity

The resolution of a fluid has a large impact on the resulting visual quality. Surface complexity is increased and dissi-

(10)

Figure 4:Simulation of 20 million SPH particles. Two-way fluid-solid coupling is realized with [AIS^∗12].

pation is reduced with increasing number of particles. Sim- ulations in the order of tens of million particles are, however, very challenging to compute within a given time frame on a desktop computer. Therefore, adaptive sampling techniques have been explored that follow the idea to allocate resources to visually interesting regions only. Adaptive sampling mechanisms either adapt the spatial resolution or adapt the sampling distance in time. In the following, we will sum- marize both, space and time adaptive methods.

5.1. Adaptive Spatial Discretization

Homogeneous fluid regions are often discretized with unnecessarily large particle numbers. Space adaptive methods thus try to resample particles in order to more accurately resolve regions with larger flow activity employing non- uniform particle radii. In order not to introduce high pressure variations, sampling techniques must maintain particle reg- ularity while staying computationally efficient. This survey classifies these level-of-detail approaches into two classes:

dynamic particle refinement methods and multi-scale techniques.

5.1.1. Dynamic Particle Refinement

The idea of a particle refinement is to dynamically exchange particle sets, either globally [CPK02] or locally [FB07], in order to increase the resolution in regions of complex flow. To introduce as little sampling errors as possible, a local error functional must be minimized [FB07], which in general requires Lagrange multipliers [LB95]. However, in computer graphics, simpler but faster sampling operators are used [DC99,KW02,LQB05,KAD^∗06,APKG07,ZSP08, OK12]. In high resolution regionsH, particles split to a fixed number of equally sized child particles, and vice versa merge to a single parent particle in low resolution regionsL. Reso- lution levels directly contribute to each other yielding a con- sistent fluid representation.

However, non-uniform smoothing radii lead to difficulties in reproducing quantity fields at level boundaries [BOT01].

Sampling errors are reduced by indirect interaction between particle levels [KAD^∗06] or by slowly changing smoothing radii over several layers [APKG07]. Still, sampling operators do not consider contributions from neighboring particles. SPH differentials are quite sensitive to the resulting irregular particle structure [BT07]. Thus, in order to maintain integration time steps equal to non-adaptive fluid solvers, an error-bound temporal blending of resolution levels is advantageous [OK12]. An example for an adaptively sampled SPH fluid using smooth blending between particle levels is outlined in Alg.5.

Algorithm 5Dynamic refinement SPH.

whileanimatingdo

compute physics forLandH

blend quantities between transitioning particles update blend weights using an error estimation identify split/merge regions

split particles inL/merge particles inH

To avoid temporal field discontinuities, new particles are smoothly blended in, while old particles are faded out using blend weights for all particles in transition. Blend weights are updated with respect to the introduced sampling error, leading to small transition regions as shown in Fig.5.

5.1.2. Multi-scale Methods

As an alternative to dynamically merging and splitting particles, level-of-detail can be achieved by using multiple resolution levels which are simulated in separate but coupled simulations [SG11,HS13]. The particle radius is fixed in each level, but arbitrarily large differences of particle sizes can be used between levels. The basic idea is presented in [SG11] for two levels and is summarized in Alg.6. The complete fluid is discretized and simulated in the coarse level L. A subset region of the fluid is simulated with higher resolution in the fine levelH. The boundary condition for the subset inHis defined by boundary particles. They contribute

c

(11)

Algorithm 6Two-scale SPH.

whileanimatingdo compute physicsL determine regions inL

transfer region information fromLontoH add / delete boundary particles inH

interpolate quantities fromLonto boundary inH fornSubstepsdo

compute physicsH

advect boundary particles inH update parent particle

interpolate feedback fromHontoL

to the physics of nearbyHparticles but are advected by the velocity of the closest particle inL, referred to as parent particle. Boundary particles can transition into the active region ofH, thus a relaxation scheme is necessary to restore the correct sampling density and to stabilize the system. Bound- ary particles are dynamically emitted and removed from the simulation. A feedback force transfers the velocities back fromHtoL. An example is illustrated in Fig.6.

This approach is extended in [HS13] to multiple resolution levels. A different boundary emission scheme is used to ensure conservation of the total mass, hence retaining the advantage of single-scale SPH solvers. The relaxation scheme is extended to support stability with complex collision boundaries.

Figure 5:In [OK12], resolution is smoothly changed via a temporal blending of particle levels (orange). An error estimation stabilizes integration time steps during blending.

Figure 6: Two-scale simulation [SG11]: dark blue corresponds to L, light blue to H (view frustum), and red to the boundary layer of H. Bottom row shows an overlay of the particle rendering and the surface computed with [SSP07].

Compared to the single-scale solution, computation time is decreased by a factor of 3 to 12. The overall performance gain depends on the size of the visually interesting subset and hence the particle count in the high-resolution region.

5.2. Adaptive Time Discretization

In contrast to an adaptive space discretization, the time domain can be adaptively sampled as well. In the following, time-adaptive methods are classified into globally adaptive schemes using a single but dynamic integration time step for all particles and a locally adaptive time stepping employing individual time steps for each particle.

5.2.1. Globally Adaptive Time Steps

The convergence of SPH simulations is bound to the CFL condition (see Sec.1). The idea of a globally adaptive time discretization [DC96,IAGT10,GP11] is to adjust the integration time step in each simulation step with respect to the

(12)

CFL condition. Thereby, the time step automatically adapts to the flow dynamics, hence avoiding complicated time step tuning.

However, in iterative EOS methods, the maximum force F^maxor maximum velocityv^maxis not known a priori. This might lead either to too restrictive time steps or to stability problems in case the dynamics are underestimated. To cope with this, [IAGT10] proposes to change the time step max- imally by 0.2% per simulation step. In some cases, the rela- tively small time step change might still not be sufficient to resolve large density fluctuations, requiring a separate shock handling mechanism.

5.2.2. Locally Adaptive Time Integration

Applying the same time step globally to all particles may restrict computation speed if precise integration is only locally required. Thus, the idea of locally adaptive methods [DC96,DC99] is to use individual time steps for each particle at the cost of non-symmetric particle interactions, ignoring the action equals reaction principle. In [DC96], each particle evaluates forces at the largest possible time step which satisfies stability restrictions. Individual particle states are then synchronized at∆tintervals.

In [GP11], barely active fluid regions are completely de- activated. Particles are divided into active particles and inac- tive particles which are static until the velocity magnitude or distance to the fluid surface triggers reactivation. Compared to non-iterative EOS solvers with constant time stepping, a speedup of 6.7 has been reported in [GP11] for scenarios with up to 1 million particles.

6. Multiphase Fluids

The SPH algorithm can be adapted for multi-fluid flows (e.g., [MSKG05,TM05b,HA06,SP08,SB12]). In contrast to mesh-based methods where diffusion is introduced at the interface, the particle representation offers the advantage that the interface between two fluids is sharply defined. In this survey, we first discuss the literature on liquid-liquid interaction, and then present the works on modeling the liquid-air interface. Further, we outline the literature on surface tension forces.

6.1. Liquid-Liquid Interface

Multiple fluids with small density ratios, i.e., up to a factor of 10, can be simulated with a standard SPH solver (Sec.1.2) by only changing the mass and rest density of the particles according to the fluid type [MSKG05]. The particle’s rest volume remains constant, thusm^a/ρ^a₀=m^b/ρ^b₀for two flu- idsaandb. The rest volume is invariant to changing the mass, provided that the rest density is adapted accordingly.

With the one-fluid formulation for the density and force computation, e.g., [Mon92,MCG03], discontinuities across

Figure 7:Standard SPH results in spurious behavior introduced by discontinuities across the interface of two fluids with density ratios (left). Adapted equations for multiple fluids according to [SP08] eliminate the artifacts (right).

the interface are smoothed due to the SPH nature of sum- ming up values from neighboring particles inside a compact support. Hence, the density of a particle computed with Eq. (3) is spuriously influenced by the density of the fluid on the other side of the interface. Consequently, pressure fields and the acceleration of a particle are affected, resulting in spurious interface tension effects [Hoo98] and a large gap between the fluids [AMS^∗07]. The problem is intensified if larger density ratios (>10) are used, introduc- ing instability issues which are unrelated to the time step size [CL03,SP08]. This is illustrated in Fig.7.

To account for the density discontinuity across the interface, the standard SPH equations can be adapted in terms of the number densityδi=∑jWi j [TM05b,HA06,SP08].

The density of a particleiis then computed as ˜ρi=miδi. This formulation offers the advantage that the density of particleiis not influenced by the mass of its neighbors j, even thoughireceives the geometric contributionW_{i j} from j. Hence, ifiand jbelong to different fluids, the computed density is not unphysically affected and sharp density pro- files are achieved. The pressure ˜piis then computed from ˜ρi

with Eq. (9).

In [SP08], ˜ρi and ˜pi are substituted into the pressure term of the Navier Stokes equations, resulting in the pressure force F^pressure = −∇p/δ. Employing the quotient˜ rule [Mon92] yields∇p/δ˜ =∇(p/δ) +˜ p/δ˜ ²∇δ. After applying the SPH rule and replacingV by 1/˜ρ, the pressure force equation is written as

F_i^pressure=−

∑

j

(p˜j

δj2+ p˜i

δi2)∇Wi j. (16) Though derived differently, Eq. (16) is applied in a similar form in [TM05b,HA06]. The viscosity force for multi-fluid systems can be derived similarly [SP08]. Note that, applied to a single fluid, the equations correspond to the standard one-phase flow formulation. By applying them to multiple fluids, however, spurious tension effects are avoided and stability is improved. The system can be further extended with an incompressibility condition [HA07].

c

(13)

Alternatively, the continuity equation can be used to compute the density in multi-phase flows [OS03,MR13]. The latter work is based on a Lagrangian with the continuity equation as a constraint. The continuity equation evolves the density over time bydρ/dt=−ρ∇ ·v[Mon94]. This formulation avoids the inherent issues of the summation equation, and hence is also applicable to free surface problems.

They typically require higher-order time stepping schemes and careful considerations of time step sizes to avoid accumulation of integration errors and thus drift from true mass conservation [SP08,SB12]. In [GAC^∗09], the density summation equation is used in combination with a Shepard kernel to accurately preserve the discontinuity at the interface.

The presented kernel is based on the volume distribution and the rate of change of the volume estimated by the continuity equation.

In the computer graphics literature, the density summation equation is preferred over evolving the density with the continuity equation. Avoiding error accumulation and thus drift from the rest volume is especially relevant for applications that require large time steps and robust results even if a comparably large time frame is simulated.

6.2. Liquid-Air Interface

In contrast to the liquid-liquid interface (Sec.6.1), the problems at the liquid-air interface are intensified due to the lack of air particles. Thus, the modified equations for multiphase fluids have no effect. Instead of discretizing the sur- rounding air completely, methods have been proposed to seed air particles only in areas close to the free surface, e.g., [MSKG05,IBAT11,SB12].

One-way coupling of air and water particles is presented in [SB12]. A narrow layer of air particles is dynamically generated, and quantities, such as the mass and velocity, are extrapolated. The air density is set toρ0 considering the p=0 surface boundary condition. Air particles only contribute to the density field and otherwise do not interact with fluid particles. They are passively advected with the velocity field of the fluid. Thus, particle resampling is necessary to reduce drift away from the liquid. However, large smoothing radii are required in order to smooth sampling errors. Air particles in [SB12] are only employed to account for particle deficiency and, thus, to correct the density estimates at the free surface. However, air particles are not used to model the entrainment of air.

Although large density ratios can be modeled with [SP08], small buoyant volumes are damped as soon as the fluid approaches the equilibrium state [LSRS99,SP08]. It is ar- gued that in this situation pressure forces arrange the particles in a stable lattice configuration which is difficult to break up. In order to model fast rising air bubbles, a buoy- ancy force can be added to counteract this effect, artificially simulating the high density ratio between water and air,

e.g., [MSKG05,CPPK07,IBAT11]. In these methods, particles are dynamically generated and deleted at the liquid- air interface which can be computed based on the gradient of the color field∇c[MSKG05]. In contrast to [MSKG05], two-way coupling of air and water particles in [IBAT11] is not realized via the pressure field, but according to the velocity field. This allows to simulate realistic drag and lift effects at comparatively large time steps which would cause numerical instabilities for both, standard SPH [MSKG05] and the modified multi-phase method [SP08].

Besides convection, physically-based formulation of reaction-diffusion processes at free surfaces requires stable modeling of phase singularities. The key concept of the implicit definition as given by [OHB^∗13] is to assign to each particle a value estimating its surface area, yielding a narrow band of surface particles as shown in Fig.8. Like this, quantity fields defined on the surface can be consistently represented, enabling anisotropic diffusion effects and robust handling of thin fluid sheets as they arise in combination with multiple rigid objects [AIS^∗12].

6.3. Surface Tension

In this survey, we distinguish between macroscopic and mi- croscopic views of the surface tension model.

6.3.1. Macroscopic Models

The macroscopic models, also referred as the continuum surface force model (CSF) [BKZ92], are based on a color fieldc which has a jump across the interface of two phases [BKZ92, Mor00,MCG03,MSKG05]. The color value of a particle iis interpolated with (3). Then,∇cis computed, yielding the surface normalnpointing into the fluid. It can either be computed with the original SPH formulation for the gradient [MCG03,MSKG05], by applying Eq. (6) [GAC^∗09], or by considering the color differenceci jbetween neighboring particles to get more accurate estimates of the surface nor- mal [Mor00,SP08]. The curvatureκis given asκ=−∇²c.

The resulting tension force then acts perpendicular to the interface of two phases, trying to minimize the curvatureκ.

In [MCG03,MSKG05,SP08], forces are only considered if the length of the normals are large enough. In order to ensure that tension forces are acting only on the interface of multiple liquids and not at the free surface, the color field has to be normalized [SP08].

Alternatively, the surface tension can be expressed as the divergence of the stress tensor as in [HA06,GAC^∗09]. The stress tensorTis again given by the color field and can be computed asT=σ1/|∇c|(1/3I|∇c|²− ∇c∇c^T)[HA06].

In [GAC^∗09], the formulation provides surface tension effects only between the two liquids but not at the surface.

However, the approximation of the Laplacian with SPH for computing surface normals is error prone. As shown