Smoothed Particle Hydrodynamics Techniques for the Physics Based Simulation of Fluids and Solids
Dan Koschier, Jan Bender, Barbara Solenthaler, Matthias Teschner
EEUROGRAPHICS 2019/ W. Jakob and E. Puppo
Tutorial
Smoothed Particle HydrodynamicsTechniques for the Physics Based Simulation of Fluids and Solids
Dan Koschier , Jan Bender , Barbara Solenthaler , and Matthias Teschner University College London, UK RWTH Aachen University, Germany ETH Zurich, Switzerland University of Freiburg, Germany
Disclaimer : This document is work-in-progress and will be maintained and updated over a longer period. Please find thecurrent version of the document under
InteractiveComputerGraphics.github.io/SPH-Tutorial . Abstract
Graphics research on Smoothed Particle Hydrodynamics (SPH) has produced fantastic visual results that are unique across theboard of research communities concerned with SPH simulations. Generally, the SPH formalism serves as a spatial discretizationtechnique, commonly used for the numerical simulation of continuum mechanical problems such as the simulation of fluids,highly viscous materials, and deformable solids. Recent advances in the field have made it possible to efficiently simulatemassive scenes with highly complex boundary geometries on a single PC [Com16b, Com16a]. Moreover, novel techniquesallow to robustly handle interactions among various materials [Com18,Com17]. As of today, graphics-inspired pressure solvers,neighborhood search algorithms, boundary formulations, and other contributions often serve as core components in commercialsoftware for animation purposes [Nex17] as well as in computer-aided engineering software [FIF16].This tutorial covers various aspects of SPH simulations. Governing equations for mechanical phenomena and their SPH dis-cretizations are discussed. Concepts and implementations of core components such as neighborhood search algorithms, pres-sure solvers, and boundary handling techniques are presented. Implementation hints for the realization of SPH solvers for fluids,elastic solids, and rigid bodies are given. The tutorial combines the introduction of theoretical concepts with the presentationof actual implementations.
Keywords:
Physically-based animation, Smoothed Particle Hydrodynamics, fluids, elastic solids, rigid bodiesCategories and Subject Descriptors (according to ACM CCS) : Computer Graphics [I.3.7]: Three-Dimensional Graphics andRealism—Animation
1. Introduction
The SPH concept is increasingly popular in a large variety of ap-plication areas that range from entertainment technologies to en-gineering. On the one hand, this popularity is based on the factthat Lagrangian approaches in general – and SPH in particular –can naturally handle scenarios that would be rather involved forEulerian approaches. A favorable example is a free-surface fluidwith geometrically complex and dynamic solid boundaries. Suchsettings are especially relevant for special effects productions inindustry. The scenario, however, has the same relevance in engi-neering, e.g. , for the analysis of vehicles in water passages, for theprediction of rain water evacuation on a vehicle with moving wipersor for the design of gear boxes with optimized lubrication.A second important aspect for the impressive advances in SPHbased techniques is the fact that various research communities con- tribute to different aspects of SPH simulations.
E.g. , the simulationcommunity has a strong focus on the accuracy of SPH discretiza-tions or on specific properties of the discretizations. Kernel func-tions and the effect of the size of kernel support domain are inves-tigated. Effects of the sampling quality onto SPH approximationsare analyzed, leading to concepts such as kernel gradient correc-tion, particle shift, ambient pressure or density diffusion, just toname a few. The computer science community – the graphics com-munity in particular – focuses on efficient algorithms for neigh-borhood searches, efficient pressure solvers, and flexible boundaryhandling. Also, pre- and post-processing is a typical graphics topic, e.g. , boundary sampling and visualization. The graphics commu-nity also experiments with combinations of different discretizationconcepts.
E.g. , some projects have started to investigate combi-nations of SPH and MLS discretizations which is less typical in c (cid:13) (cid:13)(cid:13)
E.g. , some projects have started to investigate combi-nations of SPH and MLS discretizations which is less typical in c (cid:13) (cid:13)(cid:13) a r X i v : . [ c s . G R ] S e p . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids the simulation community, where we currently see a strong focuson SPH within Lagrangian approaches with exclusive SPH confer-ences and SPH initiatives.Still, simulation and computer science are different communi-ties, but there is a growing acceptance of advances across commu-nities. Graphics papers use state-of-the-art kernel functions, rang-ing from cubic spline to Wendland kernel types. The kernel gradientcorrection is employed in a growing number of approaches. Viceversa, the simulation community adopts efficient data structures forneighborhood searches, concepts for non-uniformly boundary sam-plings, and efficient pressure solvers.This tutorial aims at a practical introduction of the SPH conceptand its application in the simulation of fluids, elastic solids, andrigid solids. It starts with a description of the SPH concept and itsusage for the interpolation of field quantities and for the computa-tion of spatial derivatives. Then, the governing equations for fluidsand solids are stated and the SPH concepts for the simulation of flu-ids and solids are outlined. In the following, various aspects of SPHsimulations are explained in more detail. One of these aspects is theneighborhood search that is required for all SPH computations, asthe interpolation of a quantity or a spatial derivative is always com-puted as a sum over adjacent particles. Another important aspectis incompressibility which is not only relevant for fluids, but also, e.g. , in the case of elastic solids. Next, boundary handling conceptsare explained, e.g. , the interaction for fluid particles at solid walls,at free surfaces, i.e. , at the interface between fluid and air, or theinteraction of particles from different fluids, i.e. , multiphase fluids.Other topics are viscosity, surface tension, and vorticity. Further,the SPH simulation of elastic solids and SPH-based contact han-dling between rigid bodies is described. Moreover, the techniquesfor the usage of SPH discretizations in data driven fluid simulationsare presented. Finally, SPlisHSPlasH, an open-source library forthe physically-based SPH simulation of fluids and solids, is intro-duced. The most important quantities that will be used throughoutthis tutorial are summarized in Tab. 1.
2. Foundations
In this section, we introduce the fundamental concept of SPH forthe phenomenological simulation of fluids and solids. The sectionis primarily based on the excellent work of Price [Pri12] and Mon-aghan [Mon05] but, moreover, includes important theoretical andpractical insights that we have gained over the years working onSPH based techniques.We first show how the SPH formalism discretizes spatial quan-tities using a set of particles equipped with a kernel function. Sec-ondly, we discuss the approximation quality that can be expectedand provide practical examples to illustrate the consequences forphysics-based simulations targeting computer graphics applica-tions. Thirdly, we show how 1 st - and 2 nd -order differential oper-ators are discretized and present specialized variants of the discreteoperators tailored to specific circumstances. Finally, we give a briefintroduction of the conservation law of linear momentum and theconcept of stress in order to derive the governing equations forfluids and elastic solids and present a simple approach to simu-late weakly compressible fluids using the knowledge that we havegained up to this point. Variable Description Unit d Spatial dimension – A Auxiliary function – t Time s ρ Volumetric mass density kgm − p Pressure Pa m Mass kg Ψ Pseudo-mass kg µ Dynamic viscosity Pas ν Kinematic viscosity m s − h Smoothing length m (cid:126)
Kernel support radius m˜ h Particle size m σ Kernel normalization factor m − x Position vector of material particle m r Distance vector between two material particles m u Displacement of a material particle m v Velocity vector of material particle ms − a Acceleration vector of material particle ms − ω Angular velocity vector of material particle rads − f Body force Nm − F Force N τ Body torque Nmm − τ Torque Nm Θ Microinertia m s − T Cauchy stress tensor Nm − P st Piola-Kirchhoff stress tensor Nm − J Deformation gradient – ε Strain tensor – E Strain rate tensor s − Table 1: Table of notation.
The concept of SPH can be generally understood as a method forthe discretization of spatial field quantities and spatial differentialoperators , e.g. , gradient, divergence, curl, etc . In order to under-stand the basic idea, we first have to introduce the Dirac- δ distri-bution and the corresponding Dirac- δ identity. δ is a generalizedfunction defined as δ ( r ) = ∞ if r = (cid:82) δ ( r ) dv = i.e. , m = : (cid:82) ρ ( x ) dv . However, if an idealized point mass isconsidered, the concept of a density function loses its meaning asthe point mass has no spatial extents. In this case the density cannot be described as a function, anymore, but collapses to the Dirac- c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids − . − . . . . x . . . . . N ( , σ ) σ = σ = σ = σ = σ = σ = Figure 1: Gaussian bell function (normal distribution) N ( , σ ) with varying variance σ . For σ → δ distribution. The arrow indicates a function value of ∞ .The function family has non-compact support. δ distribution scaled using the point mass. Another intuition of in-terpreting the Dirac- δ distribution is to understand it as the limit ofthe Gaussian normal distribution as the variance approaches zero(see Fig. 1).Now that we have understood the Dirac- δ distribution, we canapply the Dirac- δ identity as the basis for the discretization. Theidentity states that the convolution of a continuous compactly sup-ported function A ( x ) with the Dirac- δ distribution is identical to A itself, i.e. , A ( x ) = ( A ∗ δ ) ( x ) = (cid:90) A (cid:0) x (cid:48) (cid:1) δ (cid:0) x − x (cid:48) (cid:1) dv (cid:48) , (2)where dv (cid:48) denotes the (volume) integration variable correspondingto x (cid:48) . We will later approximate the integral of Eq. (2) using a sum for nu-merical quadrature. Since δ ( r ) is, however, neither a function norcan be discretized, we first make a continuous approximation to theDirac- δ distribution as a preparation to the discrete approximationof the integral. A natural choice to approximate δ is to use a nor-malized Gaussian since δ is equal to the normal distribution withzero variance. Consequently, convolving a field quantity A with aGaussian effectively smoothes A . We will later see that the Gaus-sian is, however, not an optimal choice due to its non-compact sup-port domain and will therefore consider more general smoothingfunctions W : R d × R + → R which we will refer to as kernel func-tions or smoothing kernels . Formally the continuous approximationto A ( x ) with W ( r , h ) is A ( x ) ≈ ( A ∗ W ) ( x )= (cid:90) A (cid:0) x (cid:48) (cid:1) W (cid:0) x − x (cid:48) , h (cid:1) dv (cid:48) , (3) where h denotes the kernel’s smoothing length. The smoothinglength controls the amount of smoothing and consequently howstrongly the value of A at position x is influenced by the valuesin its close proximity. This means the smoothing effect increaseswith growing smoothing lengths. The following properties are fur-thermore desired: (cid:90) R d W (cid:0) r (cid:48) , h (cid:1) dv (cid:48) = normalization condition) lim h (cid:48) → W (cid:0) r , h (cid:48) (cid:1) = δ ( r ) ( Dirac- δ condition) W ( r , h ) ≥ positivity condition) W ( r , h ) = W ( − r , h ) ( symmetry condition) W ( r , h ) = (cid:107) r (cid:107) ≥ (cid:126) , ( compact support condition) ∀ r ∈ R d , h ∈ R + , where (cid:126) denotes the support radius of the kernelfunction. Moreover, the kernel should be at least twice continuouslydifferentiable to enable a consistent discretization of 2 nd -order par-tial differential equations (PDEs). It is essential to use a kernel thatsatisfies the first two conditions (normalization and Dirac- δ ), in or-der to ensure that the approximation in Eq. (3) remains valid. Thepositivity condition is not strongly required (there are also kernelsthat do not have this property). However, in the context of physicalsimulations kernels that take negative values may lead to physi-cally inconsistent estimates of field quantities, e.g. , negative massdensity estimates, and should therefore be avoided. We will latersee that the symmetry condition ensures 1 st -order consistency ofthe continuous approximation. Finally, ensuring that the kernel iscompactly supported is a purely practical consideration that willcome into play after discretizing the continuous integral and willbe discussed later. To keep this tutorial practical, we refrain fromdiscussing how to construct SPH kernels and would like to refer thereader to the review of Liu and Liu [LL10] for a discussion on ker-nel construction and an overview over a range of smoothing kernelssuitable for SPH.A typical choice for the smoothing kernel is the cubic spline ker-nel W ( r , h ) = σ d ( q − q ) + ≤ q ≤ ( − q ) for < q ≤
10 otherwise , (4)with q = h (cid:107) r (cid:107) . The kernel normalization factors for the respectivedimensions d = , , σ = h , σ = π h , and σ = π h . Pleasenote that there exist different formulations of the cubic spline ker-nel throughout SPH literature that are differently parametrized withrespect to h . This kernel fulfills all of the discussed kernel proper-ties and has the particular advantage that its smoothing length isidentical to the kernel support radius, i.e. , h = (cid:126) , which helps toavoid confusions in the implementation. For a graphical illustra-tion please see Fig. 2. The plots demonstrate that the kernel is C -continuous. Therefore, derivatives of order > nd -order derivatives solely based on the kernel gradient. Otherwise, if c (cid:13) (cid:13)(cid:13)
10 otherwise , (4)with q = h (cid:107) r (cid:107) . The kernel normalization factors for the respectivedimensions d = , , σ = h , σ = π h , and σ = π h . Pleasenote that there exist different formulations of the cubic spline ker-nel throughout SPH literature that are differently parametrized withrespect to h . This kernel fulfills all of the discussed kernel proper-ties and has the particular advantage that its smoothing length isidentical to the kernel support radius, i.e. , h = (cid:126) , which helps toavoid confusions in the implementation. For a graphical illustra-tion please see Fig. 2. The plots demonstrate that the kernel is C -continuous. Therefore, derivatives of order > nd -order derivatives solely based on the kernel gradient. Otherwise, if c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids W ( r , h ) h = √ h = √ h = √ h = √ − − ∂ W ( r , h ) ∂ r · − h = √ h = √ − . − .
25 0 .
00 0 .
25 0 . r − − − ∂ W ( r , h ) ∂ r · − Figure 2: Graph of the cubic spline kernel (see Eq. (4)) and itsderivatives.desired, kernels of higher regularity can be found in the literature, e.g. , in the work of [LL10].Let us consider the field A : R d → R . In order to investigate theaccuracy of the continuous approximation, a Taylor series expan- sion of A in x (cid:48) about x can be applied, i.e. , ( A ∗ W )( x ) = (cid:90) (cid:20) A ( x ) + ∇ A | x · ( x (cid:48) − x ) + ( x (cid:48) − x ) · ∇∇ A | x ( x (cid:48) − x ) + O (( x (cid:48) − x ) ) (cid:21) W (cid:0) x − x (cid:48) , h (cid:1) dv (cid:48) (5) = A ( x ) (cid:90) W (cid:0) x − x (cid:48) (cid:1) dv (cid:48) + ∇ A | x · (cid:90) ( x − x (cid:48) ) W (cid:0) x − x (cid:48) (cid:1) dv (cid:48) + O (( x − x (cid:48) ) ) . (6)It is trivial to see that the approximation of ( A ∗ W ) to A is 1 st -orderaccurate if the integral in the first term of Eq. (6) becomes 1, and ifthe integral in the second term vanishes. The first condition is au-tomatically fulfilled if the kernel is normalized ( cf. , normalizationcondition). The second condition is met if the kernel is symmetric( cf. , symmetry condition). Consequently, given a normalized, sym-metric kernel we can expect that the approximation is (at least) ableto exactly reproduce functions up to 1 st -order. The remaining step to realize the SPH discretization is to replacethe analytic integral in Eq. (3) by a sum over discrete samplingpoints as follows: ( A ∗ W )( x i ) = (cid:90) A (cid:0) x (cid:48) (cid:1) ρ ( x (cid:48) ) W (cid:0) x − x (cid:48) , h (cid:1) ρ (cid:0) x (cid:48) (cid:1) dv (cid:48) (cid:124) (cid:123)(cid:122) (cid:125) dm (cid:48) (7) ≈ ∑ j ∈F A j m j ρ j W (cid:0) x i − x j , h (cid:1) = : (cid:104) A ( x i ) (cid:105) , (8)where F is the set containing all point samples and where all fieldquantities indexed using a subscript denote the field evaluated at therespective position, i.e. , A j = A (cid:0) x j (cid:1) . For improved readability, wewill drop the second argument of the kernel function and use theabbreviation W i j = W (cid:0) x i − x j , h (cid:1) in the remainder of this tutorial.The physical interpretation of this is that we keep track of a numberof points that "carry" field quantities. In this particular case, eachpoint j has a certain location x j and carries a mass sample m j anda field sample A j . It is not mandatory that the particle keeps trackof its density ρ j as this field can be reconstructed from its locationand mass as explained later. Due to the analogy to physical particlesthe term smoothed particle has been coined in the pioneering workof Gingold and Monaghan [GM77]. Nevertheless, we would liketo stress the fact that a set of SPH particles must not be misunder-stood as discrete physical particles but simply as a spatial functiondiscretization.Analogously to the brief error analysis for the continuous ap-proximation, a Taylor series expansion of (cid:104) A (cid:105) in x j about x i revealsthe accuracy of the discretization (cid:104) A ( x i ) (cid:105) = A i ∑ j m j ρ j W i j + ∇ A | x i · ∑ j m j ρ j ( x j − x i ) W i j + O (( x j − x i ) ) . (9) c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Due to the discretization the resulting approximation is only 1 st -order accurate if ∑ j m j ρ j W i j = ∑ j m j ρ j ( x j − x i ) W i j = . (10)Even presuming that a normalized symmetric kernel is used, theconditions are highly dependent on the sampling pattern leading tothe fact that not even a 0 th -order consistent discretization can beguaranteed. In practice, however, the approximation is sufficientlyaccurate to approximate physical field functions to obtain realisticsimulations. If desired, 0 th -order consistency can be easily restoredby normalizing the SPH approximation with ∑ j m j ρ j W i j or even 1 st -order consistency can be restored by the cost of a small matrix in-version (see [Pri12]).To give the reader a notion of the quality of the discrete approx-imation of functions, we have discretized a linear and a quadraticpolynomial as well as a trigonometric function using a fairly coarseSPH discretization equipped with the cubic spline kernel. The sam-pling pattern is illustrated in Fig. 3 while the the function and ap-proximation graphs are depicted in Fig. 4. In this example we have s Figure 3: Point sampling of rectangular domain. Test functions arediscretized using SPH. Function values are sampled along the redpath parametrized by s .used the cubic spline kernel with a smoothing length of h = . m i = i.e. , h = h . We, also recommend this heuristic to estimate a goodsmoothing length in practice. In three-dimensional discretizationsthis leads to a number of approx. 30 −
40 particles in a fully popu-lated neighborhood.Although no consistency can be strongly guaranteed in the ab-sence of certain particle configurations that strongly fulfill the con-ditions in Eq. (10), the graphs demonstrate that even a coarse sam-pling results in a discretization with good accuracy away from theboundary of the particle set. The phenomenon of decreasing ap-proximation quality in the close proximity of the domain bound-ary can be simply explained by the lack of sampling points outsidethe domain and is usually referred to as boundary deficiency . Inthe course of this tutorial practical solutions to this particular prob-lem will be discussed. We would also like to assure the reader thateven without further considerations to recover the consistency or-der, SPH based approaches are able to produce robust and highly- − A pp r o x i m a t i o n A ( x , y ) = ( x + y ) A ( x , y ) = ( x + y ) − A ( x , y ) = sin 5 x cos 3 y h A ih A ih A i s − . . . E rr o r A − h A i A − h A i A − h A i Figure 4: Comparison between analytic test functions and accord-ing SPH discretization using the sampling pattern illustrated inFig. 3.realistic results as demonstrated in countless publications that havebeen published within recent decades.
As previously mentioned, it is not required that the particles "carry"the mass density field as it can be reconstructed. Evaluating thedensity field at position x i using the SPH discretization in Eq. (8)results in ρ i = ∑ j m j W i j (11)and is therefore solely dependent on the sample position and themass field. Alternatively, the density can be tracked by discretizingthe mass density field using the SPH sampling and by numericalintegration of the continuity equation which describes the densityevolution, i.e. , ˙ ρ = − ρ ( ∇ · v ) . However, as also discussed by Ran-dles and Libersky [RL96], this approach is less robust and leadsto accumulating errors in the density field due to the errors of theunderlying numerical integration of the continuity equation.Note that the density can be reconstructed at any position byEq. (11) but the reconstructed density is typically underestimatedat the free surface due to particle deficiency ( cf. , Fig. 5). This mustbe considered when implementing a pressure solver as discussedin-depth in Section 4. c (cid:13) (cid:13)(cid:13)
As previously mentioned, it is not required that the particles "carry"the mass density field as it can be reconstructed. Evaluating thedensity field at position x i using the SPH discretization in Eq. (8)results in ρ i = ∑ j m j W i j (11)and is therefore solely dependent on the sample position and themass field. Alternatively, the density can be tracked by discretizingthe mass density field using the SPH sampling and by numericalintegration of the continuity equation which describes the densityevolution, i.e. , ˙ ρ = − ρ ( ∇ · v ) . However, as also discussed by Ran-dles and Libersky [RL96], this approach is less robust and leadsto accumulating errors in the density field due to the errors of theunderlying numerical integration of the continuity equation.Note that the density can be reconstructed at any position byEq. (11) but the reconstructed density is typically underestimatedat the free surface due to particle deficiency ( cf. , Fig. 5). This mustbe considered when implementing a pressure solver as discussedin-depth in Section 4. c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Figure 5: Particle deficiency problem. The green particle has a fullneighborhood while the red particle at the free surface has no neigh-bors in the air. Therefore, its density is underestimated by Eq. (11).
Besides the discretization of field quantities, it is usually necessaryto discretize spatial differential operators in order to numericallysolve physical conservation laws. In the remainder of this tutorial,we will assume that the smoothing length h is constant in space(and time). Based on the discrete SPH approximation in Eq. (8) thegradient of the underlying field can be approximated straightfor-wardly using ∇ A i ≈ ∑ j A j m j ρ j ∇ W i j . (12)Given discrete representations of higher-dimensional functions, e.g. , A : R d → R n , even more complex first-order spatial differ-ential operators can be directly discretized, e.g. , ∇ A i ≈ ∑ j m j ρ j A j ⊗ ∇ W i j (13) ∇ · A i ≈ ∑ j m j ρ j A j · ∇ W i j (14) ∇ × A i ≈ − ∑ j m j ρ j A j × ∇ W i j , (15)where a ⊗ b = ab T denotes the dyadic product. Unfortunately, these"direct" derivatives lead to a poor approximation quality and unsta-ble simulations. For this reason many discrete differential operatorshave emerged over time.In this tutorial, we will cover the two most widely used formula-tions for first order derivatives, i.e. , the difference formula and the symmetric formula . Difference Formula
Analyzing the error in the gradient based on Taylor series expan-sion (similar to the one carried out in Eq. (9)) reveals that the gradi-ent estimate is only 0 th -order (1 st -order) accurate if the first (both)of the following constraints are fulfilled: (cid:104)∇ (cid:105) = ∑ j m j ρ j ∇ i W i j = and ∑ j m j ρ j ( x j − x i ) ⊗ ∇ i W i j = (cid:49) . (16)In order to recover 0 th -order accuracy we can simply subtract thefirst error term of the Taylor series resulting in the improved ap- proximation ∇ A i ≈ (cid:104)∇ A i (cid:105) − A i (cid:104)∇ (cid:105) = ∑ j m j ρ j ( A j − A i ) ∇ i W i j . (17)In the rest of this thesis we will refer to this gradient estimateas difference formula . The same formula can be straightforwardlyapplied to the higher-dimensional first-order differential operatorspresented in Eqs. (13) to (15). This gradient estimate finally resultsin a more accurate discretization but keep in mind that we still ex-pect a linear error. However, linear accuracy is sometimes requiredand can be restored at the cost of solving a small linear equationsystem per evaluation, i.e. , ∇ A i ≈ L i (cid:32) ∑ j m j ρ j ( A j − A i ) ∇ i W i j (cid:33) L i = (cid:32) ∑ j m j ρ j ∇ i W i j ⊗ ( x j − x i ) (cid:33) − . (18) Symmetric Formula
Motivated from classical mechanics for hydrodynamical systems,a discrete formula for the pressure force/gradient, starting from thediscrete Lagrangian and the density estimate, can be derived. Thisresults in the following approximation ∇ A i ≈ ρ i (cid:32) A i ρ i (cid:104)∇ ρ (cid:105) + (cid:104)∇ (cid:18) A i ρ i (cid:19) (cid:105) (cid:33) = ρ i ∑ j m j (cid:32) A i ρ i + A j ρ j (cid:33) ∇ i W i j . (19)Please note, that we did not include the lengthy derivation as thisis out of the scope of this tutorial but kindly refer the reader to thereport of Price [Pri12].Since this formula does not satisfy the constraints in Eq. (16),it is clear that it is not able to exactly reproduce constant or lineargradient functions. However, the massive advantage of this is thatdiscrete physical forces using this particular gradient estimate ex-actly conserve linear and angular momentum which is an essentialcriterion for robust simulations.Deriving the criterion using Taylor series expansion of Eq. (19)reveals that the constant error of the symmetric gradient is governedby how much ∑ j m j (cid:32) ρ i + ρ j (cid:33) ∇ i W i j ≈ (20)deviates from 0. As noted by Price [Pri12], the symmetric formu-lation "cares" about the particle ordering and the discrete physicalforces will try to reorder the particle configuration until Eq. (20) isfulfilled. This is in contrast to forces formulated using the differ-ence formula.To summarize, the difference formula does indeed lead to a moreaccurate gradient estimate than the symmetric formula. In the con-text of physical forces the higher accuracy comes at the cost of aloss in momentum conservation and can therefore lead to unstablesimulations. For the stated reasons, we recommend to use gradient c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids estimates of the symmetric type when quantities are discretized thatdirectly affect particle trajectories, e.g. , physical forces, impulses,and to use the difference formula when 1 st -order differentials areestimated for indirect use, e.g. , the velocity divergence during pres-sure solves. Similar to the direct 1 st -order derivatives (Eqs. (13)-(15)) theLaplace operator can be directly discretized, i.e. , ∇ A i ≈ ∑ j m j ρ j A j ∇ i W i j . (21)This, however, leads again to a very poor estimate of the 2 nd -orderdifferential. A improved discrete operator for the Laplacian waspresented by Brookshaw [Bro85]: ∇ A i ≈ − ∑ j m j ρ j A i j (cid:107)∇ i W i j (cid:107)(cid:107) r i j (cid:107) . (22)The main idea leading to this particular formulation is to effectivelyuse solely a 1 st -order derivative of the kernel function and to realizethe second derivative using a finite-difference-like operation, i.e. ,dividing by the particle distance.2 nd -order derivatives of vectorial field quantities are realizedanalogously resulting in ∇ A i = − ∑ j m j ρ j A i j (cid:107)∇ i W i j (cid:107)(cid:107) r i j (cid:107) (23) ∇ ( ∇ · A i ) = ∑ j m j ρ j (cid:2) ( d + )( A i j · ˜ r i j ) ˜ r i j − A i j (cid:3) (cid:107)∇ i W i j (cid:107)(cid:107) r i j (cid:107) , (24)where d and ˜ r = r ij (cid:107) r ij (cid:107) denote the spatial dimension and the normal-ized distance vector between particles i and j , respectively. A prob-lem of the discrete Laplace operator defined in Eq. (23) in the con-text of physics simulations is that forces derived using this operator, e.g. , viscosity forces, are not momentum conserving. Fortunately,we get the following expression by adding together Eqs. (23) and(24): ∑ j m j ρ j ( A i j · ˜ r i j ) ˜ r i j (cid:107)∇ i W i j (cid:107)(cid:107) r i j (cid:107) ≈ ∇ ( ∇ · A i ) d + − ∇ A i ( d + ) . (25)This identity has the important consequence that in the case of adivergence-free vector field, i.e. , ∇ · A =
0, the Laplace operatorcan be discretized using ∇ A i ≈ ( d + ) ∑ j m j ρ j A i j · r i j (cid:107) r i j (cid:107) ∇ i W i j , (26)resulting in forces composed of terms that solely act along the "lineof sight" between two interacting particles i and j . This particu-lar choice has the advantage that derived physical forces recovermomentum conservation [Pri12]. Therefore, we recommend to useEq. (26) as discrete Laplace operator for divergence-free vectorfields. In order to improve readability, we will drop the differentia-tion index for differential operators in the remainder of this tutorial.We will use the convention that the spatial operators always differ-entiate with respect to the variable according to the first index suchthat e.g. , ∇ W i j ≡ ∇ i W i j . In order to simulate the dynamic behavior of fluids and solids, amathematical model that describes physical phenomena and mo-tion of the matter is required. In computer graphics related research,we are generally interested in the appearance of objects and flu-ids in motion on humanly perceivable scales which is dominantlygoverned by the matter’s macroscopic behavior. An important classof mathematical models that describe the macroscopic mechanicalbehavior of fluids and solids is based on continuum theory. Un-fortunately, we can not cover an introduction to continuum me-chanics as this is out of the scope of this tutorial. For a thoroughintroduction we would like to refer the reader to the works of Abe-yaratne [Abe12] and Lai et al. [LKR09]. Nevertheless, we wouldlike to informally describe the basic idea of continuum theoreticalmodels in the following.Physics teaches us that all matter is formed out of discrete par-ticles such as atoms, molecules, etc . Therefore, we know that thedistribution of mass within matter is not continuous but can ratherbe interpreted as a system of discrete mass points. Nonetheless, thevast majority of macroscopic mechanical phenomena can be ac-curately described when the corresponding matter is idealized as acontinuum, i.e. , a region of continuously distributed mass. This ide-alization then implies that a portion of matter can always be dividedinto smaller portions independent of the size of the regions. This inturn confirms the theoretical existence of a material particle , i.e. ,a portion of matter contained in an infinitesimal volume. Contin-uum theory then aims to model macroscopic physical phenomenaand neglects effects that can be observed on microscales. In thefollowing, we will summarize the most important local conserva-tion laws required for the numerical simulation of (in)compressiblefluids and solids. Continuity Equation
The continuity equation describes the evolution of an object’s massdensity ρ over time, i.e. , D ρ Dt = − ρ ( ∇ · v ) , (27)where D ( · ) Dt denotes the material derivative . This relation is espe-cially important when incompressible materials are modelled. Inthis particular case the constraint D ρ Dt = ⇔ ∇ · v = A note on the material derivative:
The material derivative de-scribes the time rate of change of a field quantity at a materialpoint. It is important to understand that the explicit form of thematerial derivative is dependent on the type of coordinates that areused to the describe the system.
Eulerian coordinates describe afield quantity at spatially fixed points in space, observing the mo-tion of the continuum as time passes. This type of coordinates isusually employed for mesh-based simulation techniques for fluids.In contrast,
Lagrangian coordinates "track" the individual mate-rial particles as they move through space and time. Lagrangian co- c (cid:13) (cid:13)(cid:13)
Lagrangian coordinates "track" the individual mate-rial particles as they move through space and time. Lagrangian co- c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids ordinates are commonly employed for the particle based simula-tion of fluids, such as SPH, or the mesh-based simulation of elasticsolids. Given the same field quantity once described in Euleriancoordinates A E ( t , x ) and Lagrangian coordinates A L ( t ) the materialderivative has the following explicit forms DA E Dt = ∂ A E ∂ t + v · ∇ x A E and DA L Dt = ∂ A L ∂ t . (29)The second term of the material derivative for Eulerian coordinatesis referred to as convection term or self-advection term. As opposedto some people’s beliefs, the convection term is non-existent whena quantity is described in Lagrangian coordinates. In the remain-der of this tutorial, we will exclusively describe quantities usingLagrangian coordinates. The conservation law of linear momentum can be interpreted as ageneralization of Newton’s second law of motion for continua andis also often called the equation of motion . It states that the rate ofchange of momentum of a material particle is equal to the sum ofall internal and external volume forces acting on the particle, i.e. , ρ D x Dt = ∇ · T + f ext , (30)where T denotes the stress tensor and f ext body forces – we un-derstand a body force as a force per unit volume. This equation isindependent of the material of the underlying matter as the mate-rial’s behavior is "encoded" in the stress tensor and described usingso-called constitutive laws. Navier-Stokes Equation
A typical constitutive relation for incom-pressible flow is T = − p (cid:49) + µ ( ∇ v + ∇ v T ) , (31)where p and µ denote the pressure and dynamic viscosity of thefluid. If the incompressibility is intended to be strongly enforced,the pressure p can be interpreted as a Lagrange multiplier that hasto be chosen such that the Eq. (28) is fulfilled. If strong enforce-ment of incompressibility is not required, the constitutive relationcan be instead completed by a so-called state equation that relatesgeometric compression (changes in mass density) with the pres-sure, i.e. , p = p ( ρ ) (see Section 4.4). A simple example for a stateequation is a variation of the ideal gas equation that linearly penal-izes deviations from a rest density ρ scaled by a positive stiffnessfactor k resulting in p ( ρ ) = k (cid:16) ρρ − (cid:17) .By plugging Eq. (31) into Eq. (30) we arrive at the incompress-ible Navier-Stokes equation ρ D v Dt = −∇ p + µ ∇ v + f ext . (32) Elasticity
The stress tensor of elastic solids is solely dependent onthe geometric deformation of an object, e.g. , T = T ( J ) , where J denotes the deformation gradient which will be later introduced inSection 10. Obviously, the constitutive model can be augmentedaccordingly, if viscoelastic, plastic, thermoelastic, or other defor-mation inducing phenomena have to be modeled. The previously introduced linear momentum conservation law(Eq. (30)) in combination with a constitutive relation, e.g. , Eq. (31),is a PDE in time and space that describes the motion of any objectcomposed of the material modeled by the constitutive law. In or-der to model a specific problem and to ensure a unique solution,initial conditions, i.e. , the initial shape and velocity of the objectat every point, and boundary conditions constraining the positionand/or velocity field have to be specified. As there is, in general, noknown analytic solution to the mixed initial-boundary value prob-lem in arbitrary scenarios, numerical solving is inevitable and re-quires discretization of the associated differential operators. In theprevious sections, we have seen several discrete differential opera-tors based on the SPH formalism that can be employed to discretizethe spatial differential operators. After spatial discretization, we areleft with a system of ordinary differential equations (this method-ology is often called method of lines ) that is typically discretizedusing standard time integration schemes such as the implicit or ex-plicit Euler method, Runge-Kutta schemes, etc . In the remainder ofthis tutorial, we will see several variations of these discretizationstailored to specific problems in physics based simulation.
Before we will discuss a simple example of a complete simulationloop, the concept of operator splitting is introduced. Its importanceis emphasized by the fact that the vast majority of today’s SPHbased simulators follow the concept. The basic idea is to decom-pose the underlying PDE, e.g. , the Navier-Stokes equation in thecase of fluids, into several sequential subproblems and to employindividual techniques for solving each subproblem. This simplifiesthe complexity of the overall problem and sometimes also decou-ples field variables such as velocity and pressure in the numeri-cal solver. It moreover allows us to use stable implicit updates forstiff subproblems while cheap explicit updates for the remainingterms can be used. A potential operator split for the incompress-ible Navier-Stokes equation (Eq. (32)) for low-viscous fluids withstrong enforcement of the incompressibility constraint (Eq. (28))might look as described in the following. Given the current geom-etry of the continuum x ( t ) and its velocity field v ( t ) at time t , wesplit the problem into a sequence of subproblems in order to obtain x ( t + ∆ t ) and v ( t + ∆ t ) :1. Update v by solving D v Dt = ν ∇ v + ρ f ext ,2. determine ∇ p by enforcing D ρ Dt = v by solving D v Dt = − ρ ∇ p and4. update x by solving D x Dt = v ,where ν = µ ρ denotes the kinematic viscosity. In this way, the"weaker" forces could be handled using explicit time integrationwhile we can solve for the pressure gradients using a more sophisti-cated implicit solver in order to keep the simulation robust for largetime steps. It should further be noticed that the individual steps arenot performed in parallel but the updated fields (in this case v and ∇ p ) are fed forward into the next substep resulting in a somewhatimplicit handling which has demonstrated to improve stability inpractice as also discussed by Bridson [Bri15] for grid-based fluidsimulation. c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids As previously described, an SPH discretization of the underlyingPDE leaves us with a system of ordinary differential equations(ODEs) in time following the method of lines. This, of course,requires us to discretize the ODE in time. Due the operator split-ting approach, as introduced in the previous section, each indi-vidual subproblem has to be numerically integrated in time. The-oretically, a different time integration scheme can be employedfor each individual step. In, practice most methods mainly rely onsimple and efficient explicit time integration schemes. The, by far,most frequently used scheme is the semi-implicit Euler scheme, as e.g. , employted in [BK17, IAAT12, SB12, ICS ∗ ∆ t † . However,we also understand that choosing an overly large time step widthresults in decreased accuracy of the numerical approximation andmay lead to a less stable simulation which might ultimately result ina breakdown of the simulation. In the context of computer graph-ics research, we care most about carrying out a robust and stablesimulation in a resource-efficient manner while the numerical ac-curacy is often of subordinate importance. This does not mean thatwe do not care about accuracy at all, as the realism of the resultinganimations often improves with better accuracy of the numericalapproximation. We are simply putting a higher priority on main-taining a robust and stable simulation under extreme conditions andin highly complex scenarios than on achieving the highest possibleaccuracy.In order to find a "good" time step width ∆ t that is as large as pos-sible to achieve high performance but sufficiently small to maintainstability, the vast majority of approaches adaptively estimate thetime step using a heuristic based on the Courant-Friedrichs-Lewy(CFL) condition. The CFL condition is a necessary condition forthe convergence of numerical solvers for differential equations and,as a result, provides an upper bound for the time step width, i.e. , ∆ t ≤ λ ˜ h (cid:107) v max (cid:107) , (33)where ˜ h , v max , and λ denote the particle size, the velocity at whichthe fastest particle travels and a user-defined scaling parameter, re-spectively. The intuition behind this condition is that all particlesare only allowed to move less than the particle diameter per timestep for λ =
1. As this is only a necessary but no generally sufficientcondition, the scaling parameter is heuristically chosen to keep thesimulation stable, i.e. , λ ≈ . † We will later see that larger time step widths not always result in bet-ter performance. This is especially true when iterative pressure solvers areemployed (see Sec. 4). dition typically leads to stable simulations [SP09, ICS ∗
14, BK17].Although obvious from Eq. (33), we would like to stress the factthat the maximally allowed time step decreases with higher veloci-ties and spatial resolution. We would further like to point out that itis in practice useful to specify global bounds, i.e. , a lower and upperbound, for the time step as we want to produce a certain number offrames per second and want to avoid that the simulation comes tohalt if a single particle moves with very high velocity.
Based on the knowledge that we have acquired up to this point, weare now able to implement a simple state-equation based simulatorfor weakly compressible fluids with operator splitting using SPHand symplectic Euler integration. for all particle i do Reconstruct density ρ i at x i with Eq. (11) for all particle i do Compute F viscosity i = m i ν ∇ v i , e.g. , using Eq. (23) v ∗ i = v i + ∆ tm i ( F viscosity i + F ext i ) for all particle i do Compute F pressure i = − ρ ∇ p using state eq. and Eq. (19) for all particle i dov i ( t + ∆ t ) = v ∗ i + ∆ tm i F pressure i x i ( t + ∆ t ) = x i + ∆ t v i ( t + ∆ t ) Algorithm 1: Simulation loop for SPH simulation of weakly com-pressible fluids.The few lines in Algorithm 1 are already enough to implementa simple fluid solver. However, the algorithm does, unfortunately,not handle boundary conditions. A practical workaround to modelboundaries in the discrete model is to sample the boundary geome-try with static (non-moving) fluid particles. The pressure forces willthen "push away" particles that attempt to penetrate the boundary. Amore consistent handling of boundary conditions will be discussedin Section 5.
3. Neighborhood Search
A major insight that we can gain from Algorithm 1 is that eval-uating the individual force terms is rather inefficient. It requiresto compute the previously defined discrete differential operatorswhich in turn require to compute a sum over all particles result-ing in a runtime complexity of O ( n ) , where n is the number ofparticles. If we, however, use a smoothing kernel that fulfills the compact support condition , most terms of the sums vanish sincethe kernel function and its derivatives for particles that are furtheraway from i than the kernel support radius (cid:126) vanish. Assuming thatwe have a list of neighbors for each particle i that lie within a ra-dius of (cid:126) around i , the algorithmic complexity reduces to O ( mn ) ,where m is the maximum number of neighboring particles. In prac-tice, m is usually bounded by a constant such that we can expectlinear runtime complexity, i.e. , O ( n ) .The problem of finding the neighbor list is commonly referred to c (cid:13) (cid:13)(cid:13)
A major insight that we can gain from Algorithm 1 is that eval-uating the individual force terms is rather inefficient. It requiresto compute the previously defined discrete differential operatorswhich in turn require to compute a sum over all particles result-ing in a runtime complexity of O ( n ) , where n is the number ofparticles. If we, however, use a smoothing kernel that fulfills the compact support condition , most terms of the sums vanish sincethe kernel function and its derivatives for particles that are furtheraway from i than the kernel support radius (cid:126) vanish. Assuming thatwe have a list of neighbors for each particle i that lie within a ra-dius of (cid:126) around i , the algorithmic complexity reduces to O ( mn ) ,where m is the maximum number of neighboring particles. In prac-tice, m is usually bounded by a constant such that we can expectlinear runtime complexity, i.e. , O ( n ) .The problem of finding the neighbor list is commonly referred to c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids (cid:126) (cid:124) (cid:123)(cid:122) (cid:125) (cid:126) Figure 6: Uniform grid approach to efficiently find neighboring par-ticles. Particles that are closer to the orange particle than the kernelsupport radius (cid:126) must lie within the one ring of the cell containingthe orange particle if the grid size is ≥ (cid:126) .as the fixed-radius near neighbor problem and is widely addressedin the computational geometry literature. The naïve approach, i.e. ,brute-force, has a computational complexity of O ( n ) and is there-fore not optimal. In this section, we will present an algorithm toapproach the problem in a computationally more efficient way, i.e. ,compact hashing [IABT11]. The basic idea of the approach is toplace a uniform grid over the domain spanned by the particles witha grid cell size equal to the kernel support radius (cid:126) . Assumingthat a particle is located in the grid cell represented by the tuple c = ( i , j , k ) , where i , j , and k denote row, column, and depth columnof the cell in the grid. Then, it is obvious that we only have to queryfor potentially neighboring particles in the cell c itself and its one-ring, i.e. , ( i − , j − , k − ) , ( i − , j − , k ) , ( i − , j − , k + ) , . . . , ( i + , j + , k + ) . The strategy then results in an algorithm witha computational complexity of O ( n ) for construction and O ( ) tofind the neighbors of a single particle implying O ( n ) to find the setof all the neighbors of all of the particles. Obviously, the grid-basedapproach can easily be generalized to higher dimensions. Please seeFig. 6 for a graphical illustration. As discussed before, the uniform grid based approach results ina good computational runtime complexity. However, there is stillpotential for optimizations in terms of memory consumption, cacheefficiency, and parallel processing. In this regard, the concept of compact hashing was proposed by Ihmsen et al. [IABT11] and willbe explained in the following. An open source
C++ implementationof a variant of this approach can be found online [Kos19].A particular disadvantage of spatial grids is that memory for allcells in the grid has to be allocated although only a small numberof cells might be occupied by particles. Due to the curse of dimen- H a s h t a b l e Secondarystructure1 2 ... n-1 n z-sorted ......1 1 1 1 ... ... k k k kCell-wiseparticle indices
Figure 7: Hash table of size m pointing to secondary data struc-ture (of n non-empty cells) to realize compact hashing. Each arrayto store the particle index lists in the secondary structure preallo-cates k entries to minimize the number of memory allocations inthe course of the simulation.sionality the memory requirements increase quickly with increas-ing domain size. It would be more memory efficient to only storethe populated cells and, hence, employ a sparse representation ofthe grid. Therefore, Ihmsen et al. suggest to store the grid cells ina hash map by hashing the index tuple c = ( i , j , k ) to a scalar indexfollowing [THM ∗ ( c ) = [( p i ) XOR ( p j ) XOR ( p k )] mod m , (34)where p = p = p = m is the hash table size. Pleasenote, that it generally cannot be avoided that several spatial cellsare mapped to the same hash value (hash collision). The effect ofoverpopulated entries in the hash table might lead to a slow-downof the neighborhood query. However, as suggested by Teschner etal. [THM ∗
03] the number of hash collisions can be reduced byincreasing the hash table size, i.e. , trading memory for speed. Asnoted by Ihmsen et al. the hash table is usually sparsely filled whenused in conjuction with SPH discretizations. Therefore, we wouldlike to avoid to unnecessarily preallocate a large amount of mem-ory. Moreover, the cache-hit rate of this approach can not expectedto be optimal as the cells that are spatially close are not necessarilyclose in memory.In order to reduce the frequency of allocations, Ihmsen et al.suggest to only store a handle per hash table entry that points toa secondary data structure – a contiguous array of the populatedcells (see Fig. 7). Each item of the secondary structure stores a listof the particle indices contained in the respective cell. In this waymemory for a used cell is only allocated if it contains particles andthe memory can be (optionally) deallocated if the cell gets empty.Each storage for the index arrays in the secondary data structurecan be further preallocated with the maximally expected number c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids of particles in a cell. To summarize, the memory consumption nowscales linear with the number of particles and not with the volumeof the simulation domain.As it lies in the nature of spatial hash tables to scatter data ac-cording to spatially close cells, the indirection to the secondary datastructure allows us to optimize for spatial locality in memory. Torealize this, Ihmsen et al. suggest to sort the non-empty cells ac-cording to a space-filling Z-curve. The cache hit-rate can further beoptimized by analogously sorting the per-particle data in the sameway. However, performing the actual sort ( O ( n log n ) ) causes com-putational overhead and since the particles are constantly movingthroughout space during the simulation, it is advised to update theZ-sort in fixed intervals, e.g. , after every 1000 th time step. This isjustified as the order is expected to be roughly maintained over asmall number of time steps due to temporal coherence.Finally, several operations such as the hash table construction,updates and neighborhood queries can be (partially) parallelized tofurther optimize performance. For further details on the approach,we would like to refer the reader to the according original pa-per [IABT11].
4. Pressure Solvers
Incompressibility is an essential aspect in realistic fluid simula-tions. The fluid volume should not noticeably oscillate or generallygrow or shrink over time. Fluid solvers preserve the fluid volume bycomputing a pressure acceleration − ρ ∇ p where the pressure p isproportional to the volume deviation. Then, the term − ρ ∇ p accel-erates particles from high pressure, i.e. , regions with large volumedeviations, to low pressure, i.e. , regions with small volume devi-ations. If there would be no volume deviation everywhere in thefluid, the pressure would be zero and the pressure gradient and thepressure acceleration would also be zero.Solver implementations typically distinguish pressure acceler-ation a p = − ρ ∇ p and all other non-pressure accelerations a nonp which improves the intuition of the incompressibility concept.First, a predicted velocity is computed with, e.g. , v ∗ = v ( t ) + ∆ t a nonp ( t ) . Then, pressure is computed from the volume deviationafter advecting the fluid with v ∗ . Finally, the respective pressureacceleration would be applied as, e.g. , v ( t + ∆ t ) = v ∗ + ∆ t a p ( t ) tominimize the volume deviation. This final velocity update is oftenreferred to as pressure projection which is related to the fact thatthe velocity change ∆ t a p ( t ) should be minimal. I.e. , the pressureacceleration should change the velocity field as little as possible.Conceptually, pressure is proportional to the volume deviation.However, there exist various alternatives to actually compute thepressure. First, the volume deviation can be explicitly computedfrom the density or the velocity divergence can be used to computea differential update of the volume deviation. Second, pressure canbe computed locally with a state equation or it can be computedglobally by solving a Pressure Poisson Equation (PPE). The firstaspect determines whether the fluid volume oscillates or continu-ously changes, while the second aspect influences the solver per-formance.
The volume deviation is typically deduced from the density devi-ation. Although SPH solvers can easily handle both formulations,it is probably due to historical reasons that the density formula-tion is preferred over the volume formulation. The SPH density ata particle i is computed with ρ i = ∑ j m j W i j and the deviation tothe rest density ρ is considered for the pressure computation. Notethat the density deviation is often clamped, e.g. , max ( ρ i − ρ , ) or max (cid:16) ρ i ρ − , (cid:17) , as a simple solution to the particle deficiencyproblem at the free surface (see Fig. 5). The continuity equation relates the time derivative of the densityto the velocity divergence: D ρ D t = − ρ ∇ · v . This fact can be usedto predict a particle density from its previous density and, e.g. , thepredicted velocity: ρ ∗ i = ρ i ( t ) − ∆ t ρ i ( t ) ∇ · v ∗ i . Here, ρ ∗ i is a pre-diction of the particle density after advecting the particles with v ∗ i for time ∆ t . If it is assumed that the current density equals the restdensity, the predicted density is computed as ρ ∗ i = ρ − ∆ t ρ ∇ · v ∗ i which means that ρ ∗ i − ρ = − ∆ t ρ ∇ · v ∗ i can be used as a measurefor the density deviation. It can be seen that minimizing the den-sity deviation is related to minimizing the velocity divergence. Theterm − ∆ t ρ ∇ · v ∗ i is a density change at a particle if the particlesare advected with v ∗ i for time ∆ t . Both forms imply challenges. If pressure accelerations are derivedfrom the explicit form of the volume deviation, the fluid volumeoscillates due to an over-correction of the pressure acceleration.These oscillations have to minimized. At least, they should not beperceivable. Using the differential form to compute the volume de-viation results in a drift of the fluid volume, typically a volumeloss. The differential form assumes that the current density is cor-rect. It minimizes density changes between simulation steps, butpotentially existing density deviations are not detected or corrected.Here, the challenge is to minimize the volume drift. Although vol-ume drift often occurs in Eulerian pressure solvers and volume os-cillations often occur in Lagrangian solvers, both issues are not re-lated to the Eulerian or Lagrangian perspective. If an SPH solverwas using the differential form to compute the density deviation, itwould suffer from volume drift. If a Eulerian solver, e.g. , FLIP, wasusing the explicit form for the density computation, it would sufferfrom oscillations.
State equations are used to compute pressure from density devi-ations. The density deviation can be computed explicitly or froma differential form. The deviation can be represented as a quo-tient or a difference of actual and rest density. One or more stiff-ness constants are involved. Some examples are: p i = k (cid:16) ρ i ρ − (cid:17) , p i = k ( ρ i − ρ ) or p i = k (cid:18)(cid:16) ρ i ρ (cid:17) k − (cid:19) . As ρ i < ρ is not consid-ered to solve the particle deficiency problem at the free surface, the c (cid:13) (cid:13)(cid:13)
State equations are used to compute pressure from density devi-ations. The density deviation can be computed explicitly or froma differential form. The deviation can be represented as a quo-tient or a difference of actual and rest density. One or more stiff-ness constants are involved. Some examples are: p i = k (cid:16) ρ i ρ − (cid:17) , p i = k ( ρ i − ρ ) or p i = k (cid:18)(cid:16) ρ i ρ (cid:17) k − (cid:19) . As ρ i < ρ is not consid-ered to solve the particle deficiency problem at the free surface, the c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids computed pressure is always non-negative. SPH fluid simulationsthat use a state equation to compute pressure are often referred to ascompressible or weakly compressible. In contrast, fluid implemen-tations that solve a PPE to compute pressure are known as incom-pressible. These terms basically indicate that it is more challengingto minimize compression with a state equation than with a PPE.It might look confusing that arbitrary pressure values can becomputed for a given density ρ i > ρ dependent on the state equa-tion and the stiffness constant(s). Here, it is interesting to note thatthe parameters do not govern the pressure, but the compressibilityof the SPH fluid. This can be seen in a simple example with a fluidat rest under gravity. In this case, the pressure acceleration at allparticles cancels gravity, i.e. , g − ρ i ∇ p i = . Discretizing the pres-sure gradient with SPH yields g = ∑ j m j (cid:18) p i ρ i + p j ρ j (cid:19) ∇ W i j . Using, e.g. , p i = k ( ρ i − ρ ) , yields g = k ∑ j m j (cid:18) ρ i − ρ ρ i + ρ i − ρ ρ j (cid:19) ∇ W i j . Itcan be seen that a variation of the stiffness constant k is related toa variation in the density deviation ρ i − ρ . This relation, however,is not simply k ( ρ i − ρ ) = const since the erroneous particle sam-pling for ρ i (cid:54) = ρ influences the SPH discretization of the pressuregradient. But generally, the stiffness constant in the state equationgoverns the density deviation. Larger values result in smaller devi-ations and require smaller time steps. Smaller values lead to largerdensity deviations, i.e. , less realistic simulations. Also, the bound-ary handling fails if the tolerated density deviation is too large. The general idea of the pressure computation is to end up with pres-sure accelerations that cause velocity changes that in turn cause dis-placements such that all particles are uncompressed, i.e. , have theirrest density or rest volume. PPE solvers or projection schemes solvea linear system to compute the respective pressure field.
We consider the predicted velocity after all non-pressure acceler-ations: v ∗ = v ( t ) + ∆ t a nonp ( t ) . If the particles would be advectedwith this velocity, we can use the continuity equation to estimatea predicted density ρ ∗ = ρ ( t ) − ∆ t ρ ( t ) ∇ · v ∗ . Now, the goal ofthe pressure computation is a pressure acceleration − ρ ( t ) ∇ p ( t ) that corresponds to a velocity change − ∆ t ρ ( t ) ∇ p ( t ) whose diver-gence −∇ · (cid:16) ∆ t ρ ( t ) ∇ p ( t ) (cid:17) corresponds to a density change pertime − ρ ( t ) ∇· (cid:16) ∆ t ρ ( t ) ∇ p ( t ) (cid:17) that cancels the predicted density de-viation per time ρ − ρ ∗ ∆ t , i.e. , ρ − ρ ∗ ∆ t − ρ ( t ) ∇ · (cid:16) ∆ t ρ ( t ) ∇ p ( t ) (cid:17) = ∆ t ∇ p ( t ) = ρ − ρ ∗ ∆ t . (35)Note that ∇ · ∇ p ( t ) = ∇ p ( t ) . In this equation, the pressure p isunknown. We have one equation per particle, resulting in a sys-tem with n equations and n unknown pressure values for n parti-cles. Various similar PPE forms can be derived more formally, e.g. ,starting with the continuity equation at time t + ∆ t : − ρ ( t + ∆ t ) ∇ · v ( t + ∆ t ) = D ρ ( t + ∆ t ) D t . The time derivative of the density is approx-imated with D ρ ( t + ∆ t ) D t = ρ ( t + ∆ t ) − ρ ( t ) ∆ t . The velocity is written as v ( t + ∆ t ) = v ∗ − ∆ t ρ ( t + ∆ t ) ∇ p ( t + ∆ t ) using an implicit update withthe pressure acceleration at time t + ∆ t . Imposing the constraint ρ ( t + ∆ t ) = ρ , we get − ρ ∇ · v ∗ + ∇ · ( ∆ t ∇ p ( t + ∆ t )) = ρ − ρ ( t ) ∆ t and finally ∆ t ∇ p ( t + ∆ t ) = ρ − ( ρ ( t ) − ∆ t ρ ∇ · v ∗ ) ∆ t . (36)If compared carefully, there are very minor differences betweenEqs. (35) and (36) in the computation of the predicted density andthe pressure acceleration. As ρ ( t ) ≈ ρ , however, these differencesare negligible. The biggest difference would actually be the SPHdiscretizations of the Laplacian. Eq. (35) works with the particleneighborhood at time t , while Eq. (36) requires the neighborhoodat time t + ∆ t . As this would require an additional neighbor searchper simulation step, it is generally ignored.Eqs. (35) and (36) make use of the density invariance as sourceterm in the PPE. Motivated by the continuity equation, the di-vergence of the predicted velocity could be used alternatively.To derive the respective form, we start with ∆ t ρ ( t ) ∇ p ( t ) = v ∗ − v ( t + ∆ t ) . Taking the divergence and imposing the constraint thatthe velocity field at t + ∆ t should be divergence free, i.e. , ∇ · v ( t + ∆ t ) =
0, we get ∇ · (cid:16) ∆ t ρ ( t ) ∇ p ( t ) (cid:17) = ∇ · v ∗ and ∆ t ∇ p ( t ) = ρ ( t ) ∇ · v ∗ . (37)We search for pressure p such that the pressure acceleration − ρ ( t ) ∇ p ( t ) corresponds to a velocity change − ∆ t ρ ( t ) ∇ p ( t ) whosedivergence −∇ · (cid:16) ∆ t ρ ( t ) ∇ p ( t ) (cid:17) cancels the divergence of the pre-dicted velocity, i.e. , −∇ · (cid:16) ∆ t ρ ( t ) ∇ p ( t ) (cid:17) + ∇ · v ∗ = e.g. , aFLIP issue, but only due to the typically used velocity divergenceas source term in the PPE. There exist various alternative discretizations with benefits anddrawbacks. Here, we discuss one option referred to as IISPH whichhas been proposed in [ICS ∗ i : ∆ t ∇ p i = ρ − ρ ∗ i . (38)If a quantity is considered at time t , e.g. , pressure, the time index isomitted. Using SPH, the source term is computed as ρ − ρ ∗ i = ρ − ρ i − ∆ t ∑ j m j (cid:0) v ∗ i − v ∗ j (cid:1) · ∇ W i j (39)with v ∗ i = v i + ∆ t a nonp i . The computation of the Laplacian is real-ized as the computation of the divergence of the velocity change c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids due to the pressure acceleration, i.e. , ∆ t ∇ p i = − ∆ t ρ i ∇ · (cid:0) ∆ t a p i (cid:1) = ∆ t ∑ j m j (cid:16) a p i − a p j (cid:17) · ∇ W i j (40)with pressure acceleration a p i = − ρ i ∇ p i = − ∑ j m j (cid:32) p i ρ i + p j ρ j (cid:33) ∇ W i j . (41)Using Eqs. (39) and (40), we can compute the left-hand and right-hand side of Eq. (38) at each particle using SPH sums over adjacentparticles. The IISPH discretization of Eq. (38) is ∆ t ∑ j m j (cid:16) a p i − a p j (cid:17) · ∇ W i j = ρ − ρ ∗ i . (42)Introducing the velocity change due to pressure acceleration v p i = ∆ t a p i , Eq. (42) can be written as ∆ t ∑ j m j (cid:16) v p i − v p j (cid:17) · ∇ W i j = ∆ t ρ i ∇ · v p i = ρ − ρ ∗ i . (43)It can be seen that we search for pressure values such that the pres-sure acceleration causes a velocity change v p i such that the diver-gence of v p i causes a density change that corrects from the predicteddensity ρ ∗ i to the rest density ρ , i.e. , ρ ∗ i + ∆ t ρ i ∇ · v p i = ρ . Eq. (42) is considered at all particles. We have n equa-tions with n unknown pressure values p i . Each equation does notonly contain an unknown pressure value p i , but also unknown pres-sure values at neighbors of i and - due to Eq. (41) - also unknownpressures at neighbors of neighbors of i . The set of Eq. (42) at allparticles forms a system Ap = s and Eq. (42) at particle i can bewritten as ( Ap ) i = s i : ( Ap ) i = ∆ t ∑ j m j (cid:16) a p i − a p j (cid:17) · ∇ W i j (44) s i = ρ − ρ ∗ i . (45) A is not a diagonal matrix, but at least sparse. In row i , matrix ele-ments a ii and also elements a i j are non-zero with j being a neighboror a neighbor of a neighbor of i .Extracting the matrix elements a i j would not be impossible, butrather tedious and error-prone. Fortunately, element extraction isnot required to solve the system as solver implementations typi-cally just require the computation of ( Ap ) i . The term ( Ap ) i canbe computed by first evaluating a p i with Eq. (41), followed by thecomputation of ∆ t ∇ p i using Eq. (40). An explicit notion of theelements of A is typically not required as, e.g. , for Conjugate Gra-dients. Some solvers, e.g. , Jacobi, require the diagonal elements a ii . Relaxed Jacobi scheme
In the following, we discuss the imple-mentation of a Jacobi variant to solve Ap = s . The method startswith an initialization of the pressure vector, e.g. , p ( ) = . Then,the weighted / damped / relaxed Jacobi scheme iteratively updates all pressure values using p ( l + ) i = ( − ω ) p ( l ) i + ω a ii (cid:32) s i − ∑ j (cid:54) = i a i j p ( l ) j (cid:33) (46)with l indicating the iteration and ω being a relaxation coefficient.The relaxation coefficient is typically set to ω = . a ii , but also elements a i j for neighboring parti-cles are required in an implementation. Interestingly, the update inEq. (46) can be rewritten as p ( l + ) i = ( − ω ) p ( l ) i + ω a ii (cid:16) s i − ( Ap ( l ) ) i + a ii p ( l ) i (cid:17) = p ( l ) i + ω a ii (cid:16) s i − ( Ap ( l ) ) i (cid:17) . (47)This update requires the computation of the term ( Ap ( l ) ) i whichcan be computed from Eqs. (40) and (41) without notion of theelements a i j . As IISPH only considers non-negative pressure, theactual update is p ( l + ) i = max (cid:18) p ( l ) i + ω a ii (cid:16) s i − ( Ap ( l ) ) i (cid:17) , (cid:19) . (48) Pressure clamping
The clamping has been proposed in [ICS ∗ [ , p max ] to, e.g. , [ − , p max − ] or [ , p max + ] by adding constant offsets to all values doesnot change the pressure gradient. SPH discretizations, however, be-have differently for negative and positive pressure values. Reasonshave not been investigated yet, but one could speculate that negativeand positive pressures result in different pressure gradients in caseof incomplete neighborhoods where missing contributions are im-plicitly assumed to have zero pressure. In IISPH simulations, min-imum pressure is generally zero being consistent with the implicitassumption of zero pressure for missing samples at free surfaces. Diagonal element
The update in Eq. (48) requires the diagonalelement a ii . The element can be calculated by accumulating all co-efficients of p i after inserting Eq. (41) into Eq. (40). Here, it isimportant to keep in mind that i is one of the neighbors of eachneighbor of i . Finally, the diagonal element is a ii = − ∆ t ∑ j m j (cid:32) ∑ j m j ρ j ∇ W i j (cid:33) · ∇ W i j − ∆ t ∑ j m j (cid:32) m i ρ i ∇ W i j (cid:33) · ∇ W i j . (49) Stop criterion
There is no agreement on when to stop the Jacobiiterations in Eq. (48). The iterations could be stopped after a fixednumber. Alternatively, a predicted density deviation is often con-sidered. This is motivated by the fact that ( Ap ( l ) ) i is a predicteddensity change at particle i due to the pressure field at iteration l . I.e. , ρ err , ∗ i = (( Ap ( l ) ) i − s i ) / ρ = (( Ap ( l ) ) i + ρ ∗ i − ρ ) / ρ is a c (cid:13) (cid:13)(cid:13)
There is no agreement on when to stop the Jacobiiterations in Eq. (48). The iterations could be stopped after a fixednumber. Alternatively, a predicted density deviation is often con-sidered. This is motivated by the fact that ( Ap ( l ) ) i is a predicteddensity change at particle i due to the pressure field at iteration l . I.e. , ρ err , ∗ i = (( Ap ( l ) ) i − s i ) / ρ = (( Ap ( l ) ) i + ρ ∗ i − ρ ) / ρ is a c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids predicted relative density error at particle i , if pressure accelera-tions according to the pressure field p ( l ) would be applied. Typ-ically, the average of all ρ err , ∗ i is taken as a stop criterion, e.g. , ρ avg_err , ∗ i = n ∑ i | ρ err , ∗ i | . In [ICS ∗ e.g. , ρ avg_err , ∗ i < . i.e. , the oscillation of the overall fluid volume isbelow 0 . Implementation
Alg. 2 shows the implementation of IISPH.Source term s i and diagonal element a ii are computed once. In iter-ation l , the pressure Laplacian ( Ap ( l ) ) i = ∆ t ∇ p i is computed intwo steps. First, the pressure acceleration ( a p i ) ( l ) is computed andstored at the particles. Then, the divergence of the velocity changedue to the pressure acceleration, i.e. , ( Ap ( l ) ) i , is computed. Neigh-borhood search, the computation of the predicted velocity v ∗ i andthe advection of particles are omitted in Alg. 2. for all particle i do compute diagonal element a ii with Eq. (49)compute source term s i with Eq. (39)initialize pressure p ( ) i = l = repeatfor all particle i do compute pressure acceleration ( a p i ) ( l ) with Eq. (41) for all particle i do compute Laplacian ( Ap ( l ) ) i with Eq. (40)update pressure p ( l + ) i with Eq. (48) l = l + until ρ avg_err , ∗ i < . In the following a pressure solver based on a predictor-correctorapproach is introduced [SP09].
Motivation
The density ρ i ( t + ∆ t ) can be estimated with ρ i ( t + ∆ t ) = ∑ j m j W i j + ∆ t ∑ j m j ( v i − v j ) · ∇ W i j + ∆ t ∑ j m j ( ∆ t a nonp i − ∆ t a nonp j ) · ∇ W i j + ∆ t ∑ j m j ( ∆ t a p i − ∆ t a p j ) · ∇ W i j . (50)Again, the time index is omitted for quantities at time t . Using thepredicted density ρ ∗ i = ∑ j m j W i j + ∆ t ∑ j m j ( v i − v j ) · ∇ W i j + ∆ t ∑ j m j ( ∆ t a nonp i − ∆ t a nonp j ) · ∇ W i j (51) and the constraint ρ i ( t + ∆ t ) = ρ , we can write ρ = ρ ∗ i + ∆ t ∑ j m j ( ∆ t a p i − ∆ t a p j ) · ∇ W i j . (52)When using the symmetric SPH formulation for the pressure accel-eration, the terms a p i and a p j are computed from unknown pressurevalues at particle i , at neighbors of i and at neighbors of neigh-bors of i . Now, PCISPH introduces approximations and simpli-fications to end up with only one unknown pressure value p i inEq. (52) [SP09]. Then, each pressure value p i can be computedfrom one equation. Solving a linear system is avoided. Simplifications
The pressure acceleration is discretized with thesymmetric formulation. It is assumed that the pressure p j at neigh-boring particles equals the pressure p i at particle i . Further, m j = m i for all neighbors and ρ i = ρ j ≈ ρ . Then, the pressure accelerationcan be written as a p i = − ∑ j m j (cid:32) p i ρ i + p j ρ j (cid:33) ∇ W i j (53) ≈ − m i p i ( ρ ) ∑ j ∇ W i j . (54)Using this approximation, Eq. (52) can be written as ρ = ρ ∗ i + ∆ t m i ( ρ ) ∑ j (cid:32) − p i ∑ j ∇ W i j + p j ∑ k ∇ W jk (cid:33) · ∇ W i j . (55)Using p i ≈ p j and approximating ∑ k ∇ W jk ≈ ∇ W ji , the equationcan further be simplified to ρ = ρ ∗ i + p i ∆ t m i ( ρ ) ∑ j (cid:32) − ∑ j ∇ W i j + ∇ W ji (cid:33) · ∇ W i j = ρ ∗ i − p i ∆ t m i ( ρ ) (cid:32) ∑ j ∇ W i j · ∑ j ∇ W i j + ∑ j ( ∇ W i j · ∇ W i j ) (cid:33)(cid:124) (cid:123)(cid:122) (cid:125) − k PCI . (56)In PCISPH, the coefficient k PCI is considered for a template par-ticle with perfect sampling, i.e. , ∑ j ∇ W i j · ∑ j ∇ W i j + ∑ j ( ∇ W i j ·∇ W i j ) = const. Finally, we have the following equation per parti-cle: p i = k PCI ( ρ − ρ ∗ i ) (57)with k PCI = − . ( ρ ) ∆ t m i · ∑ j ∇ W i j · ∑ j ∇ W i j + ∑ j ( ∇ W i j · ∇ W i j ) (58)for a template particle i with perfect sampling. This is a state equa-tion, where the stiffness constant is not user-defined, but motivatedby the fact that Eq. (50) should be satisfied, i.e. , the pressure shouldinduce pressure accelerations such that the particles have their restdensity at time t + ∆ t . c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Iterative refinement
The locally optimized state equation is oneimportant property of the PCISPH concept. A second significantcharacteristics is the iterative refinement of the pressure field.While this sounds expensive, it is motivated by large time stepscompared to simple state-equation solvers. PCISPH computes afirst estimate of the pressure p ( ) i with Eq. (57). This predictedpressure field is used to compute ( a p i ) ( ) with Eq. (53). Then, thepressure field is iteratively refined by p ( l + ) i = p ( l ) i + k PCI ( ρ − ρ ∗ i − ∆ t ∑ j m j ( ∆ t ( a p i ) ( l ) − ∆ t ( a p j ) ( l ) ) · ∇ W i j ) (59)in iteration l . The term ( ρ p i ) ( l ) = ∆ t ∑ j m j ( ∆ t ( a p i ) ( l ) − ∆ t ( a p j ) ( l ) ) · ∇ W i j (60)is one option to predict the density change due to the pressure accel-erations. If this density change does not cancel the predicted den-sity deviation ρ − ρ ∗ i , the predicted pressure is corrected. Similarlyto IISPH, the process is stopped if ( ρ − ρ ∗ i + ( ρ p i ) ( l ) ) / ρ is suffi-ciently small, e.g. , smaller than 0 . ρ p i from particle displacements due to thepressure accelerations, i.e. , ( ρ p i ) ( l ) = m i (cid:32) ( ∆ x i ) ( l ) · ∑ j ∇ W i j − ∑ j ∇ W i j · ( ∆ x j ) ( l ) (cid:33) . (61)If ( ∆ x i ) ( l ) = ∆ t ( a p i ) ( l ) and m i = m j , Eqs. (60) and (61) are identi-cal. We prefer Eq. (60) as it will also be used in the discussion ofthe relation between PCISPH and IISPH. Implementation
Alg. 3 shows the implementation of PCISPH.The stiffness constant k PCI is computed once at the beginning ofthe simulation. The same coefficient is used for all particles. Thepressure field is predicted with the state equation in Eq. (57). Theeffect of the respective pressure acceleration onto the density is es-timated. Remaining deviations from the rest density are used tocompute pressure corrections with Eq. (59). Neighbor search andthe advection of the particles are omitted in Alg. 3.
The PCISPH pressure solver updates the pressure with p ( l + ) i = p ( l ) i + k PCI ( ρ − ρ ∗ i − ∆ t ∑ j m j ( ∆ t ( a p i ) ( l ) − ∆ t ( a p j ) ( l ) ) · ∇ W i j ) (62)which simplifies to the state equation p ( ) i = k PCI ( ρ − ρ ∗ i ) (63)for the first update if the pressure is initialized with p ( ) i =
0. Thesame applies to IISPH. The solver updates pressure with p ( l + ) i = p ( l ) i + ω a ii (cid:16) s i − ( Ap ( l ) ) i (cid:17) (64) compute stiffness constant k PCI with Eq. (58) for all particle i do compute predicted density ρ ∗ i with Eq. (51)initialize pressure p ( ) i with Eq. (57) l = repeatfor all particle i do compute pressure acceleration ( a p i ) ( l ) with Eq. (53) for all particle i do compute density change ( ρ p i ) ( l ) with Eq. (60)update pressure p ( l + ) i with Eq. (59) l = l + until ρ avg_err , ∗ i < . p ( ) i = ω a ii s i = ω a ii ( ρ − ρ ∗ i ) (65)for the first update if the pressure is initialized with p ( ) i = I.e. ,if IISPH or PCISPH stop after one pressure update, they are state-equation solvers.Another remarkable aspect are the stiffness constants in PCISPH(Eqs. (62) and (63)) and IISPH (Eqs. (64) and (65) ). The PCISPHstiffness constant is k PCI = − . ( ρ ) ∆ t m i · ∑ j ∇ W i j · ∑ j ∇ W i j + ∑ j ( ∇ W i j · ∇ W i j ) . (66)The IISPH constant is ω a ii with ω = . a ii = − ∆ t ∑ j m j (cid:32) ∑ j m j ρ j ∇ W i j (cid:33) · ∇ W i j − ∆ t ∑ j m j (cid:32) m i ρ i ∇ W i j (cid:33) · ∇ W i j . (67)Applying the same assumptions as for PCISPH, i.e. , ρ i = ρ and m i = m j , the diagonal element simplifies to a ii = − ∆ t m i ( ρ ) (cid:32) ∑ j ∇ W i j · ∑ j ∇ W i j + ∑ j ( ∇ W i j · ∇ W i j ) (cid:33) . (68)Thus, k PCI = ω a ii . (69)Further, ∆ t ∑ j m j ( ∆ t ( a p i ) ( l ) − ∆ t ( a p j ) ( l ) ) · ∇ W i j ) = ( Ap ( l ) ) i (70)which means that the pressure update with Eq. (62) in the PCISPHsolver is equal to the pressure update with Eq. (64) in the IISPHsolver. There might be insignificant differences in the computations c (cid:13) (cid:13)(cid:13)
0. Thesame applies to IISPH. The solver updates pressure with p ( l + ) i = p ( l ) i + ω a ii (cid:16) s i − ( Ap ( l ) ) i (cid:17) (64) compute stiffness constant k PCI with Eq. (58) for all particle i do compute predicted density ρ ∗ i with Eq. (51)initialize pressure p ( ) i with Eq. (57) l = repeatfor all particle i do compute pressure acceleration ( a p i ) ( l ) with Eq. (53) for all particle i do compute density change ( ρ p i ) ( l ) with Eq. (60)update pressure p ( l + ) i with Eq. (59) l = l + until ρ avg_err , ∗ i < . p ( ) i = ω a ii s i = ω a ii ( ρ − ρ ∗ i ) (65)for the first update if the pressure is initialized with p ( ) i = I.e. ,if IISPH or PCISPH stop after one pressure update, they are state-equation solvers.Another remarkable aspect are the stiffness constants in PCISPH(Eqs. (62) and (63)) and IISPH (Eqs. (64) and (65) ). The PCISPHstiffness constant is k PCI = − . ( ρ ) ∆ t m i · ∑ j ∇ W i j · ∑ j ∇ W i j + ∑ j ( ∇ W i j · ∇ W i j ) . (66)The IISPH constant is ω a ii with ω = . a ii = − ∆ t ∑ j m j (cid:32) ∑ j m j ρ j ∇ W i j (cid:33) · ∇ W i j − ∆ t ∑ j m j (cid:32) m i ρ i ∇ W i j (cid:33) · ∇ W i j . (67)Applying the same assumptions as for PCISPH, i.e. , ρ i = ρ and m i = m j , the diagonal element simplifies to a ii = − ∆ t m i ( ρ ) (cid:32) ∑ j ∇ W i j · ∑ j ∇ W i j + ∑ j ( ∇ W i j · ∇ W i j ) (cid:33) . (68)Thus, k PCI = ω a ii . (69)Further, ∆ t ∑ j m j ( ∆ t ( a p i ) ( l ) − ∆ t ( a p j ) ( l ) ) · ∇ W i j ) = ( Ap ( l ) ) i (70)which means that the pressure update with Eq. (62) in the PCISPHsolver is equal to the pressure update with Eq. (64) in the IISPHsolver. There might be insignificant differences in the computations c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids of ρ ∗ i and ( Ap ) i between PCISPH and IISPH, but both solvers areessentially equal.Ihmsen et al. [ICS ∗
14] report significant performance differ-ences between PCISPH and IISPH. These differences are possiblydue to the fact that PCISPH computes one global stiffness constant k PCI for a template particle, while IISPH computes a ii for each par-ticle. In particular, ∑ j ∇ W i j · ∑ j ∇ W i j + ∑ j ( ∇ W i j · ∇ W i j ) is as-sumed to be constant in PCISPH, while it is computed per parti-cle in IISPH. The respective difference affects the performance. If k PCI < ω a ii , the convergence of PCISPH is worse than with IISPH.If k PCI > ω a ii , PCISPH could be unstable. Such potential insta-bilities are difficult to deal with and might result in the require-ment of significantly smaller time steps compared to IISPH. Thereare also other smaller differences between PCISPH and IISPH. E.g. , PCISPH computes ρ ∗ i from advected samples without updatedneighborhood, while IISPH uses the velocity divergence to estimate ρ ∗ i . Such differences, however, are probably less relevant for stabil-ity or convergence differences. In addition to PCISPH, there exist various other pressure solversthat are closely related to a PPE solver, e.g. , Local Poisson SPH[HLL ∗ C i = ρ i ρ −
1. The mag-nitudes of these constraints are equal to pressure values computedwith a state equation p i = k (cid:16) ρ i ρ − (cid:17) with stiffness constant k = i.e. , to compute pres-sure such that the resulting pressure accelerations minimize densitydeviations or the divergence of the velocity field. Nevertheless, thecomputed velocity field depends on the employed discretizations. E.g. , Fürstenau et al. [FAW17] compare three discretizations of thepressure Laplacian. They consider the following form of the PPEwith velocity divergence as source term: ∇ · (cid:18) ρ i ∇ p i (cid:19) = ∇ · v ∗ i ∆ t . (71)Three variants to compute the pressure Laplacian are analyzed. Thefirst variant is ∇ · (cid:18) ρ i ∇ p i (cid:19) = m i ρ i ∑ j ρ i + ρ j ρ i ρ j ( p j − p i )( x i − x j )( x i − x j ) + ε · ∇ W i j , (72)referred to as finite difference scheme, where ε is a small constantwhich is added to avoid singularities. The second variant is ∇ · (cid:18) ρ i ∇ p i (cid:19) = m i ρ i ∑ j (cid:18) ∇ p j ρ j − ∇ p i ρ i (cid:19) · ∇ W i j , (73) referred to as double summation scheme. This approximation isused, e.g. , in IISPH [ICS ∗ ∇ · (cid:18) ρ i ∇ p i (cid:19) = m i ρ i ∑ j ρ i + ρ j ρ i ρ j ( p j − p i ) ∇ W i j , (74)referred to as second derivative scheme. All these options can beused to compute the left-hand side of the PPE that contains thepressure Laplacian. If the PPE is solved with a Jacobi method, therequired diagonal element a ii varies accordingly.The three discretizations result in different solutions for the ve-locity field. It is difficult to come up with general conclusions, asall discretizations have benefits and drawbacks. E.g. , the double-summation approach used in IISPH seems to have an improvedsolver convergence compared to the finite-difference scheme andthe second-derivative scheme. On the other hand, the computed ve-locity field suffers from high-frequency noise. This noise, however,is rather low and it depends on the application whether this is anissue or not. In typical free-surface scenario, this noise is not an is-sue and the double-summation discretization is just faster than theother two options as less solver iterations are required for a speci-fied tolerated density deviation.Another degree-of-freedom is the form of the source term. Asalready shown and discussed in Section 4.5.1, the source term caneither represent the divergence of the predicted velocity or the de-viation of the predicted density to the rest density. The predictedvelocity is the velocity after applying all non-pressure accelera-tions. The predicted density is the estimated density after advectingall samples with the predicted velocity. The equivalence of bothformulations follows from the continuity equation. The PPE withdensity invariance as source term is ∆ t ∇ p i = ρ − ρ ∗ i (75)with source term s i = ρ − ρ ∗ i . (76)The PPE with velocity divergence as source term is ∆ t ∇ p i = ∆ t ρ i ∇ · v ∗ i (77)with source term s i = ∆ t ρ i ∇ · v ∗ i . (78)The properties of both variants have been analyzed in [CBG ∗ ∗ c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids advect the particles with the velocity from the PPE with densityinvariance as source term. This avoids the density drift. The ve-locity from the PPE with velocity divergence, however, could beused for the final velocities at the particles as these velocities haveless artificial viscosity and less noise. This has been actually doneby Bender and Koschier [BK15] who proposed to solve two PPEswith different source terms and to use one solution to compute theadvected particle positions and one solution as the final velocityfield. Solving two PPEs is obviously more expensive than a simpleIISPH solver. However, each PPE solve typically requires very fewiterations, i.e. , less than ten, and both PPEs share the same matrixwhich can be exploited to get an efficient solver as shown in thenext section. DFSPH [BK17] is a variant of the idea to solve two PPEs with dif-ferent source terms. The original DFSPH solver does not computepressure, but some stiffness parameter κ i per particle i . Accordingto the work of Band et al. [BGPT18], however, this stiffness param-eter is closely related to pressure with p i = κ i ρ i . In the followingwe introduce DFSPH using a pressure formulation and replace thestiffness parameter in order to get a formulation which is closer tothe ones of PCISPH and IISPH. In this way it is easier to see thesimilarities and the differences.DFSPH conceptually solves two PPEs, one with density invari-ance as source term, the other one with velocity divergence assource term. Instead of using one solution for the position updateand one solution as the final velocity field, DFSPH combines bothsolutions to compute the final velocity field which is used to ad-vect the particles. This combination is motivated by the fact that afirst PPE solve with density invariance computes particle positionsof an incompressible fluid state, but not necessarily a divergence-free velocity field. That’s why, a second PPE solve with velocitydivergence computes a divergence-free velocity field. Divergence-Free Solver
If we solve the PPE with velocity diver-gence as source term (see Eq. (77)), we can derive the followingequation for the corresponding pressure value p vi of a particle i : p vi = ∆ t D ρ i Dt · ρ i (cid:107) ∑ j m j ∇ W i j (cid:107) + ∑ j (cid:107) m j ∇ W i j (cid:107) (cid:124) (cid:123)(cid:122) (cid:125) k DFSPH i , (79)where the time derivative of the density is determined as D ρ i Dt = ∑ j m j ( v i − v j ) · ∇ W i j . (80)Note that this formulation is similar to the pressure computation ofPCISPH (see Eqs. (57) and (58)). But in contrast to PCISPH thesource term is the velocity divergence and not the density devia-tion. Moreover, DFSPH does not use a global factor k PCI that isdetermined for a template particle but computes the actual factor k DFSPH i for each particle i in each time step.Algorithm 4 shows the divergence-free solver. In each iterationfirst the divergence is updated for all particles. Then the pressurevalues are determined and the predicted velocity is updated accord-ingly. while (cid:16)(cid:16) D ρ Dt (cid:17) avg > η div (cid:17) ∨ ( iter < ) do for all particles i do D ρ i Dt = − ρ i ∇ · v ∗ i for all particles i do p vi = ∆ t D ρ i Dt k DFSPH i , p vj = ∆ t D ρ j Dt k DFSPH j v ∗ i : = v ∗ i − ∆ t ∑ j m j (cid:18) p vi ρ i + p vj ρ j (cid:19) ∇ W i j Algorithm 4: Divergence-free solver
Constant Density Solver
The constant density solver uses thePPE with density deviation as source term (see Eq. (75)). For thisPPE we get a pressure of p i = ∆ t ( ρ ∗ i − ρ ) k DFSPH i , (81)where the predicted density is determined by ρ ∗ i = ρ i + ∆ t D ρ i Dt = ρ i + ∆ t ∑ j m j ( v ∗ i − v ∗ j ) · ∇ W i j . (82)Note that the factor k DFSPH i is used for both, the divergence-freeand the constant density solver. Therefore, it has to be computedonly once per simulation step. Algorithm 5 demonstrates an imple-mentation of the constant density solver. while ( ρ avg − ρ > η ) ∨ ( iter < ) do for all particles i do compute ρ ∗ i for all particles i do p i = ρ ∗ i − ρ ∆ t k DFSPH i , p j = ρ ∗ j − ρ ∆ t k DFSPH j v ∗ i : = v ∗ i − ∆ t ∑ j m j (cid:18) p i ρ i + p j ρ j (cid:19) ∇ W i j Algorithm 5: Constant density solver
DFSPH Simulation Step
Algorithm 6 shows a simulation stepwith DFSPH and how both solvers are integrated in the time step.Note that the neighborhoods, the particle densities and the factor k DFSPH i are computed once at the beginning of the simulation forthe initial state and then updated once per time step. The algorithmfirst computes predicted velocities by integrating all non-pressureaccelerations. Then the density deviation is corrected using the con-stant density solver which yields new particle positions. Hence,the neighborhoods, the density values and the factors must be up-dated. After correcting the density deviation the velocity field istypically not divergence-free. This is corrected in the last step bythe divergence-free solver which gives us the final velocities. Notethat the order of the steps is a bit different than the order of othersolvers but in this way it is guaranteed that the density deviationsand the divergence error are both corrected at the end of a timestep. Moreover, in this way we have to update the factor k DFSPH i only once per time step but are able to use it twice: for the constantdensity solver and for the divergence-free solver.DFSPH solves two PPEs which is more expensive than solving c (cid:13) (cid:13)(cid:13)
Algorithm 6 shows a simulation stepwith DFSPH and how both solvers are integrated in the time step.Note that the neighborhoods, the particle densities and the factor k DFSPH i are computed once at the beginning of the simulation forthe initial state and then updated once per time step. The algorithmfirst computes predicted velocities by integrating all non-pressureaccelerations. Then the density deviation is corrected using the con-stant density solver which yields new particle positions. Hence,the neighborhoods, the density values and the factors must be up-dated. After correcting the density deviation the velocity field istypically not divergence-free. This is corrected in the last step bythe divergence-free solver which gives us the final velocities. Notethat the order of the steps is a bit different than the order of othersolvers but in this way it is guaranteed that the density deviationsand the divergence error are both corrected at the end of a timestep. Moreover, in this way we have to update the factor k DFSPH i only once per time step but are able to use it twice: for the constantdensity solver and for the divergence-free solver.DFSPH solves two PPEs which is more expensive than solving c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids for all particles i do compute non-pressure accelerations a nonp i adapt time step size ∆ t according to CFL condition for all particles i do predict velocity v ∗ i = v i + ∆ t a nonp i correct density error using algorithm 5 for all particles i do update position x t + ∆ ti = x i + ∆ t v ∗ i update neighborhoods for all particles i do update density ρ i update factor k DFSPH i using Eq. (79) correct divergence error using algorithm 4 for all particles i do update velocity v t + ∆ ti = v ∗ i Algorithm 6: Simulation step with DFSPHjust one like PCISPH or IISPH. However, the second solve is notthat expensive since the costly computation of the factor k DFSPH i hasto be performed only once per step. Experiments have shown thatsolving both PPEs leads to a better stability which enables largertime steps and therefore a faster simulation [BK17]. The perfor-mance can be further improved by using a warm start. More detailsabout this can be found in [BK17]. Iterative PPE solvers are more expensive to compute than EOSsolvers. Their utility, however, is motivated by the fact that PPEsolvers work with significantly larger time steps compared to EOSsolvers. Solenthaler and Pajarola [SP09] show an improved overallperformance of PCISPH compared to the EOS solver in [BT07].Ihmsen et al. [ICS ∗
14] show an improved performance of IISPHcompared to PCISPH and Bender and Koschier [BK17] show a per-formance gain of DFSPH compared to IISPH. So, DFSPH has thebest overall performance of all discussed PPE variants.Although PPE solvers work with big time steps, they do notreach their best overall performance for the largest possible timestep as discussed in [ICS ∗
14] and [IOS ∗ e.g. , one layer of fluid particleson a planar boundary, EOS solvers are faster. The overall numberof particles does not necessarily influence the solver performance.An EOS solver is more efficient than a PPE solver for one billionfluid particles in one layer on a plane, while a PPE solver is moreefficient than an EOS solver for one hundred particles on top ofeach other in one column under gravity.Independent from whether there is a performance gain of a PPE solver, they are more simple to handle than EOS solvers. In an EOSsolver, the stiffness constant has to be found to realize a desireddensity deviation and the time step has to be found to get a stablesimulation. In a PPE solver, the desired density deviation is ex-plicitly specified. The time step is also easier to estimate as it istypically rather larger, corresponding to CFL numbers close to one.
5. Boundary Handling
In order to complete the discretization of a mixed initial-boundaryvalue problem (see Section 2) the boundary of the simulation do-main has to be discretized and the corresponding boundary con-ditions must be enforced. In recent years, a wide variety of ap-proaches to represent boundary geometries and to enforce bound-ary conditions has been presented. The approaches can be roughlycategorized into particle-based approaches, e.g. , [AIA ∗ ∗
18, GPB ∗ e.g. ,[KB17, HKK07a, HKK07b, BLS12].The particle based strategy is probably the most popular repre-sentation type. The main idea is to sample the boundary geometryusing an additional set of so-called boundary particles equippedwith a (sometimes specialized) kernel function. The advantageshere are that the representation is consistent with the discretiza-tion of the fluid or solid. Modeling, either explicit/implicit bound-ary forces for weak satisfaction of boundary conditions or algo-rithms to strongly satisfy the constraints is probably more straight-forward than using implicit or mesh-based techniques. Most meth-ods, however, have the constraint that the particle size used to sam-ple the boundary has to be the same as the particle size of the con-tinuum discretization. The disadvantages are that even the repre-sentation of simple geometries, such as a plane, requires a largenumber of boundary particles that have to be accounted for dur-ing the neighborhood search and in the evaluation of field quan-tities, e.g. , Eq. (11). Moreover, determining "good" samplings isgenerally non-trivial. Too sparse samplings might not sufficientlycover the surface of the boundary leading to the issue of SPH par-ticles penetrating the boundary. Too dense samplings lead to an in-creased computational effort and to higher memory requirements.Also a somewhat "bumpy" sampling might lead to a bias in theparticle trajectories as the discretized surfaces’ smoothness suffersfrom sampling noise leading to unwanted perturbations in the sim-ulation ( cf. , [KB17]) if no additional considerations are made, suchas proposed by Band et al. [BGPT18].Implicit boundaries use an implicit function – typically a signeddistance field (SDF) – to represent the boundary geometry. Advan-tages of this type of approaches are that the boundary representa-tion is decoupled from the particle size. As a result, more flexibledata structures, e.g. , adaptive octrees with higher-order approxima-tions [KDBB17], can be used to memory-efficiently and accuratelyrepresent the boundary geometry. This circumstance also avoidsthe problem of noisy boundary samplings, the resulting bias, andunwanted perturbations in the particle trajectories. Typical disad-vantages are, that implicit representations do not directly integratewith particle based continuum discretizations. In order to coupleboth discretization types, special considerations have to be madeand the corresponding implementation is rather involved.In the remainder of this section, we will discuss approaches to c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids FluidparticleSolidparticles
Kernel support
Figure 8: Idea of the particle-based boundary handling. The bound-ary is represented with particles x i b . These particles are consideredin the computation of density ρ i , pressure p i , and pressure force F pi of nearby fluid particles x i . If a fluid particle moves closer tothe boundary, its density increases. If a fluid particle is too close,its density is larger than the rest density, i.e. , ρ i > ρ , which causespressure p i > F pi (cid:54) = thataccelerates the fluid particle away from boundary particles. Pleasenote that i b denote boundary neighbors of a fluid particle i . Accord-ingly, i f refer to indices of fluid neighbors of a fluid particle i .handle non-penetration of rigid boundaries using particle samplingapproaches. We gradually develop a formulation starting with asimple dense, uniform multilayer sampling of the boundary andshow how the method can be simplified to a uniform single layersampling and consequently even to a robust and consistent for-mulation using non-uniformly sampled boundaries. We moreoverdiscuss how the fluid-boundary coupling can be improved usingpressure mirroring or pressure extrapolation and how these tech-niques can be incorporated into the previously discussed pressuresolvers. For implicit or mesh-based boundary handling techniqueswe would like to refer the reader to the corresponding literature, e.g. , [KB17,HKK07a,HKK07b,BLS12,MFK ∗ ∗ In this concept, the boundary is represented with particles and theseboundary particles are incorporated into the computation of density ρ i , pressure p i , and F pi at nearby fluid particles x i . This is illustratedin Fig. 8. There are two main aspects to discuss. First: The bound-ary can be sampled in different ways. E.g. , , the boundary can besampled with particles of uniform size corresponding to the size ofa fluid particle. Or the boundary can be sampled with particles ofnon-uniform size. Also, the boundary can be sampled with severallayers of boundary particles as indicated in Fig. 8 or the boundarycan be sampled with just one layer of particles. The second aspectis the actual computation of density ρ i , pressure p i , and pressureforce F pi of nearby fluid particles x i . Such computations require in-formation from boundary samples, e.g. , pressure. Such pressure atboundary samples can be estimated in different ways. Here, typicalexamples are pressure mirroring, i.e. , it is assumed that the pres-sure at the boundary particles equals the pressure at adjacent fluidparticles. Alternatively, pressure can be extrapolated from the fluidinto the boundary. While the pressure extrapolation is theoreticallythe correct choice, this concept is challenging to realize due to thefact that the computation of the pressure gradient at a fluid particleclose to the boundary is error-prone. In the following, we discussdifferent combinations of the aforementioned variants. FluidSolid
Figure 9: Sample-based boundary representation with several lay-ers. The usage of several boundary layers avoids incomplete neigh-borhoods for fluid particles x i close to the boundary. If a fluidparticle x i is close to a boundary, it generally has fluid neighbors x i f and boundary neighbors x i b as illustrated in Fig. 9. All theseneighbors contribute to the density computation, i.e. , ρ i = ∑ i f m i f W ii f + ∑ i b m i b W ii b . (83)All samples x i , x i f , x i b have the same size and we take the boundarysamples as static fluid samples resulting in the same rest density ρ for all fluid and boundary particles. This also means that all massesare equal: m i = m i f = m i b . Thus, the density computation could alsobe written as ρ i = m i ( ∑ i f W ii f + ∑ i b W ii b ) . If the boundary samplesbelong to a rigid body, it is not perfectly intuitive to represent itsboundary with particles that have mass and rest density of a fluidparticle. This issue results from the fact that SPH formulations of-ten prefer to weight the contribution of a particle i with m i ρ i insteadof using its volume V i = m i ρ i . Now, in the boundary handling, one ac-tually works with the volume of boundary particles. Nevertheless,this volume is often represented with some rest density - typicallythe fluid rest density - and the respective mass.If the density ρ i in Eq. (83) is larger than the rest density ρ , apressure p i is computed at fluid particles. We have seen that p i canbe computed from a state equation, e.g. , p i = k (cid:16) ρ i ρ − (cid:17) or fromsolving a PPE.Now, the pressure force at the fluid particle can be computed as F pi = − m i m i ∑ i f (cid:32) p i ρ i + p i f ρ i f (cid:33) ∇ W ii f − m i m i ∑ i b (cid:32) p i ρ i + p i b ρ i b (cid:33) ∇ W ii b . (84)It can be seen in Eq. (84) that we basically compute a pressureforce component with respect to fluid neighbors and a pressureforce component with respect to boundary neighbors. Using sev-eral layers of boundary particles guarantees that the neighborhoodof a fluid particle is completely sampled, even if this particle is veryclose to the boundary. Fully filled neighborhoods keep the errorsdue to missing samples in Eqs. (83) and (84) small.The computation in Eq. (84) requires positions, densities andpressures of adjacent fluid and boundary particles. While thesequantities are known for fluid particles, density and pressure at c (cid:13) (cid:13)(cid:13)
Figure 9: Sample-based boundary representation with several lay-ers. The usage of several boundary layers avoids incomplete neigh-borhoods for fluid particles x i close to the boundary. If a fluidparticle x i is close to a boundary, it generally has fluid neighbors x i f and boundary neighbors x i b as illustrated in Fig. 9. All theseneighbors contribute to the density computation, i.e. , ρ i = ∑ i f m i f W ii f + ∑ i b m i b W ii b . (83)All samples x i , x i f , x i b have the same size and we take the boundarysamples as static fluid samples resulting in the same rest density ρ for all fluid and boundary particles. This also means that all massesare equal: m i = m i f = m i b . Thus, the density computation could alsobe written as ρ i = m i ( ∑ i f W ii f + ∑ i b W ii b ) . If the boundary samplesbelong to a rigid body, it is not perfectly intuitive to represent itsboundary with particles that have mass and rest density of a fluidparticle. This issue results from the fact that SPH formulations of-ten prefer to weight the contribution of a particle i with m i ρ i insteadof using its volume V i = m i ρ i . Now, in the boundary handling, one ac-tually works with the volume of boundary particles. Nevertheless,this volume is often represented with some rest density - typicallythe fluid rest density - and the respective mass.If the density ρ i in Eq. (83) is larger than the rest density ρ , apressure p i is computed at fluid particles. We have seen that p i canbe computed from a state equation, e.g. , p i = k (cid:16) ρ i ρ − (cid:17) or fromsolving a PPE.Now, the pressure force at the fluid particle can be computed as F pi = − m i m i ∑ i f (cid:32) p i ρ i + p i f ρ i f (cid:33) ∇ W ii f − m i m i ∑ i b (cid:32) p i ρ i + p i b ρ i b (cid:33) ∇ W ii b . (84)It can be seen in Eq. (84) that we basically compute a pressureforce component with respect to fluid neighbors and a pressureforce component with respect to boundary neighbors. Using sev-eral layers of boundary particles guarantees that the neighborhoodof a fluid particle is completely sampled, even if this particle is veryclose to the boundary. Fully filled neighborhoods keep the errorsdue to missing samples in Eqs. (83) and (84) small.The computation in Eq. (84) requires positions, densities andpressures of adjacent fluid and boundary particles. While thesequantities are known for fluid particles, density and pressure at c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Figure 10: Working with equal pressure p i at a fluid particle i andat its adjacent boundary samples results in a repulsion force in caseof p i > i.e. , ρ i b = ρ . Regarding the pressure p i b , we have already briefly mentioned pressure mirroring and pres-sure extrapolation. Here, the simplest idea is to mirror the pressurefrom a fluid particle to an adjacent boundary particle, i.e. , p i b = p i .These assumptions result in an adapted form of Eq. (84) for thepressure force where all required quantities at fluid and boundaryneighbors are known: F pi = − m i m i ∑ i f (cid:32) p i ρ i + p i f ρ i f (cid:33) ∇ W ii f − m i m i ∑ i b (cid:32) p i ρ i + p i ( ρ ) (cid:33) ∇ W ii b . (85)Working with equal pressure at a fluid sample and at its neighbor-ing boundary samples theoretically corresponds to a pressure gra-dient of zero which in turn results in a pressure acceleration of zero.In practice, however, the SPH derivative approximation always re-sults in the desired gradient with the respective repulsion force. Ifthe neighborhood of a fluid particle is completely filled, as shownin Fig. 9, the true gradient is slightly underestimated by SPH. Thisleads to a small, practically not relevant amount of penetration ofthe fluid into the boundary. In the other case, where a single fluidparticle without fluid neighbors is close to the boundary as depictedin Fig. 10, the contributions from missing fluid neighbors are im-plicitly assumed to be zero. Having a fluid neighbor with zero pres-sure or not having this fluid neighbor has the same effect on theSPH approximation.For a single fluid particle at a boundary, Eq. (85) simplifies to F pi = − p i m i ρ i ∑ i b ∇ W ii b for ρ i = ρ . The term − ∑ i b ∇ W ii b can beinterpreted as the surface normal of the boundary close to position x i . The pressure force at particle x i corresponds to this vector scaledwith the fluid particle pressure p i . So, if the density of the fluidparticle is larger than its rest density, we get a positive pressure anda repulsion force from the boundary into normal direction.In summary, computing the fluid density ρ i with Eq. (83), thepressure p i with a state equation or a PPE, and the pressure forcewith Eq. (85) realizes a boundary handling with pressure mirror-ing in case of a uniformly sampled boundary with several layers. FluidSolidMissing samples
Figure 11: One layer of uniform boundary samples is easier to gen-erate than several layers. The boundary layer is used to naturallydetect the proximity of a fluid sample to a boundary. Due to thesize of the kernel support, however, more than one layer of samplesmight be required in SPH computations. The contributions of thesemissing samples can be analytically estimated.The boundary handling works for fluid particles with a complete orincomplete neighborhood.
One layer of uniform boundary samples:
It can be difficult togenerate multiple layers of uniform boundary samples for arbitrar-ily shaped geometries. Also, if the relative position of a fluid sam-ple to the boundary can be determined, it is not necessarily requiredto explicitly represent the boundary particles. Instead, their contri-butions can analytically be estimated. Fig. 11 shows a setting withone layer of uniform boundary samples. These samples are used toestimate that a fluid particle is close to a boundary. For the com-putation of SPH approximations, however, more than one layer ofsamples might be required as indicated in Fig. 11.For a given position x i of a fluid particle close to a one-layerboundary, the density can be written as ρ i = m i ∑ i f W ii f + m i ∑ i b W ii b + m i ∑ i m W ii m . (86)The index i m refers to missing samples in the neighborhood of par-ticle i . While it is a natural way to encode the contributions of miss-ing samples as an offset, these contributions are typically encodedas a correcting factor of the contributions from boundary samplesinstead, i.e. , ρ i = m i ∑ i f W ii f + γ m i ∑ i b W ii b . (87)Correcting coefficients for contributions from missing boundarysamples have been proposed in [AIA ∗ γ in Eq. (87) depends on various as-pects, e.g. , kernel function, kernel support, dimensionality. In par-ticular, however, it depends on the position of a fluid particle rela-tive to the boundary. In practice, the correcting coefficient is deter-mined for a template particle in a perfect sampling pattern as shownin Fig. 11. If the neighborhood of a fluid particle at the bound-ary is fully filled, the kernel sum over all neighbors gives one overthe particle volume according to the general kernel properties, i.e. , ∑ i f W ii f + ∑ i b W ii b = V i . If there are samples missing, this equationdoes not hold and we introduce the correcting coefficient γ to ob- c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids tain the desired result: ∑ i f W ii f + γ ∑ i b W ii b = V i . (88)Solving this equation gives the desired value for the correcting co-efficient: γ = V i − ∑ i f W ii f ∑ i b W ii b . (89)The corrected density computation is the basis for the pressurecomputation which is typically not affected by missing boundarysamples. A state equation just works with the density of the fluidparticle itself. A PPE typically uses the densities of adjacent fluidparticles. So, boundary samples are not required in the pressurecomputation.In the subsequent computation of the pressure forces, however,the missing samples have to be accounted for. Similar to the densitycomputation, a correcting coefficient γ is introduced: F pi = − m i m i ∑ i f (cid:32) p i ρ i + p i f ρ i f (cid:33) ∇ W ii f − γ p i m i ρ i ∑ i b ∇ W ii b . (90)This coefficient is derived from a property of the kernel gradi-ent. If the neighborhood of a particle is perfectly sampled, thesum of the kernel gradient over the neighbors is zero: ∑ i f ∇ W ii f + ∑ i b ∇ W ii b = . In case of missing boundary samples, the sum is notzero and a correcting coefficient for the boundary contributions isintroduced to meet the constraint: ∑ i f ∇ W ii f + γ ∑ i b ∇ W ii b = . (91)Solving this equation results in γ = ∑ i f ∇ W ii f · ∑ i b ∇ W ii b ∑ i b ∇ W ii b · ∑ i b ∇ W ii b . (92)Similar to the coefficient γ , γ also depends on dimenionality, ker-nel function and support. It also depends on the position x i of fluidparticle i relative to the boundary. In practice, however, the coeffi-cient is determined for a template setting with perfect sampling asdepicted in Fig. 11.In summary: One-layer boundary representations are more sim-ple to generate than multi-layer boundaries. The contributions ofmissing samples in the computation of the density and the pressureforce at a fluid particle close to the boundary can be approximatedwith correcting coefficients. One layer of non-uniform boundary samples:
The next step toan even more flexible boundary representation is to work with sam-ples of arbitrary size as proposed in [AIA ∗
12] and illustrated inFig. 12. This sampling is motivated by the fact that its generationis really simple. Particles can be arbitrarily placed on a boundarygeometry, as long as each boundary particle is equal or smaller thana fluid particle. Even more than one boundary particle at the sameposition can be handled. This extreme case is certainly not optimalin terms of performance, but the boundary handling works.The basic idea of the boundary handling with non-uniformboundary samples is the consideration of the actual contribution,
FluidSolidMissing samples
Figure 12: The boundary geometry is sampled with one layer ofparticles of non-uniform size. The sampling should be sufficientlydense, i.e. , each boundary particle should be equal or smaller thana fluid particle. Otherwise, leakage could occur.
Perfect sampling Perfect sampling withmissing samples Arbitrary sampling withmissing samples
Figure 13: Derivation of the mass m i b of a boundary sample of ar-bitrary size. The artificial mass is computed from the volume of thesample and the rest density of the fluid. The left image shows thevolume for a filled neighborhood and uniform sample size. The vol-ume is ˜ h . It can also be computed as one over the sum of the kernelvalues from neighboring particles. In the middle image, we have re-moved all samples except the samples in a plane which representsthe boundary. We still know that the size of x i b should be ˜ h , butthe kernel sum would not provide the correct result. That’s why thecorrecting coefficient γ is introduced to get the expected result ˜ h . I.e. , γ accounts for missing neighbors in the volume computation.The volume computation V i b = γ ∑ ibb W ibibb also works for arbitrarysample sizes within the plane as shown in the image on the right.The artificial mass of a boundary sample is finally computed fromits volume and the rest density of a fluid particle. i.e. , the actual volume of each boundary sample. The density of afluid particle near the boundary is computed with ρ i = m i ∑ i f W ii f + ∑ i b m i b W ii b , (93)where m i b represents the contribution of boundary sample i b . If aboundary sample is bigger, its contribution is bigger and this con-tribution is encoded in the mass m i b . The contribution, i.e. , the ar-tificial mass of a boundary sample is deduced from its volume andthe rest density of the fluid as explained in Fig. 13. So, the mass m i b of a boundary sample is computed as m i b = ρ γ ∑ i bb W i b i bb . (94)The correcting coefficient γ actually accounts for missing contri-butions in two cases. It cancels missing contributions in the compu-tation of the mass m i b of a boundary sample. In parallel, it accountsfor missing contributions in the density computation of a nearby c (cid:13) (cid:13)(cid:13)
Figure 13: Derivation of the mass m i b of a boundary sample of ar-bitrary size. The artificial mass is computed from the volume of thesample and the rest density of the fluid. The left image shows thevolume for a filled neighborhood and uniform sample size. The vol-ume is ˜ h . It can also be computed as one over the sum of the kernelvalues from neighboring particles. In the middle image, we have re-moved all samples except the samples in a plane which representsthe boundary. We still know that the size of x i b should be ˜ h , butthe kernel sum would not provide the correct result. That’s why thecorrecting coefficient γ is introduced to get the expected result ˜ h . I.e. , γ accounts for missing neighbors in the volume computation.The volume computation V i b = γ ∑ ibb W ibibb also works for arbitrarysample sizes within the plane as shown in the image on the right.The artificial mass of a boundary sample is finally computed fromits volume and the rest density of a fluid particle. i.e. , the actual volume of each boundary sample. The density of afluid particle near the boundary is computed with ρ i = m i ∑ i f W ii f + ∑ i b m i b W ii b , (93)where m i b represents the contribution of boundary sample i b . If aboundary sample is bigger, its contribution is bigger and this con-tribution is encoded in the mass m i b . The contribution, i.e. , the ar-tificial mass of a boundary sample is deduced from its volume andthe rest density of the fluid as explained in Fig. 13. So, the mass m i b of a boundary sample is computed as m i b = ρ γ ∑ i bb W i b i bb . (94)The correcting coefficient γ actually accounts for missing contri-butions in two cases. It cancels missing contributions in the compu-tation of the mass m i b of a boundary sample. In parallel, it accountsfor missing contributions in the density computation of a nearby c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids fluid particle using Eq. 93. The pressure force is now computedwith F pi = − m i m i ∑ i f (cid:32) p i ρ i + p i f ρ i f (cid:33) ∇ W ii f − γ p i m i ρ i ∑ i b m i b ∇ W ii b , (95)which is very similar to the pressure force for uniform samples. Theonly difference is the consideration of the individual masses m i b ofthe boundary particles instead of the standard mass m i for samplesof uniform size. The correcting factor γ is γ = ∑ i f ∇ W ii f · ∑ i b ∇ W ii b ∑ i b ∇ W ii b · ∑ i b ∇ W ii b . (96)The same factor has already been used, motivated, and derived inthe case of a one-layer boundary with uniform samples.
6. Viscosity
Modeling and simulating viscosity is often vital for physics simula-tions as the phenomenon is responsible for various visually appeal-ing effects, such as buckling and coiling but also general energydissipation. In recent years, various methods for the realistic sim-ulation of low viscous fluids like water as well as highly viscousfluids like honey, mud, or dough were proposed. Fig. 14 shows dif-ferent examples of highly viscous materials.In this section we first introduce the term for the viscous force inthe Navier-Stokes equations. Then we discuss important methodsfor the simulation of low viscous flow and highly viscous materials.We present explicit approaches which are typically used for lowviscous fluids and we discuss implicit methods for highly viscousmaterials.
In the Navier-Stokes equations for incompressible fluids, the vis-cosity term is determined by a material parameter µ and the Lapla-cian of the velocity field ∇ v (see Eq. (32)). Before we discussdifferent approaches to compute this viscosity term, we first wantto show how this term is derived.The stress tensor of a Newtonian fluid is defined as T = − p (cid:49) + µ E , (97)where E is the strain rate tensor which is determined as E = ( ∇ v + ( ∇ v ) T ) . (98)When substituting the stress tensor in the conservation law of lin-ear momentum (see Eq. (30)) and considering the incompressibilityconstraint ∇ · v = ρ D v D t = −∇ p + µ ∇ · ∇ v (cid:124) (cid:123)(cid:122) (cid:125) ∇ v + µ ∇ · ( ∇ v ) T (cid:124) (cid:123)(cid:122) (cid:125) ∇ ( ∇· v )= + f ext , (99)where the right term of the strain rate tensor vanishes and the finalviscosity force is determined as f visco = µ ∇ v . (100) Recent viscosity solvers either use a strain rate based formulationto compute viscous forces or they directly determine the Laplacianof the velocity field. While the Laplacian formulation ensures thatthe right term of the strain rate tensor vanishes by definition, ap-proaches based on the strain rate have to enforce a divergence-freevelocity field, otherwise ∇ · ( ∇ v ) T (cid:54) = which leads to undesiredbulk viscosity [Lau11,PICT15]. In the following we introduce bothapproaches and discuss the advantages and disadvantages. In the Navier-Stokes equations for incompressible fluids the vis-cous force is defined by the Laplacian of the velocity field (seeEq. (100)). The standard SPH discretization of this Laplacian is de-termined as ∇ v i = ∑ j m j ρ j v j ∇ W i j . (101)However, this formulation has two major disadvantages. First, itis sensitive to particle disorder [Mon05, Pri12]. Second, typicallyGaussian-like kernel functions are used in SPH and their secondderivatives changes the sign inside the support radius (see Fig. 2).Different approaches were proposed to avoid this problem. First,instead of computing the second derivative directly, it can also bedetermined taking two first SPH derivatives [FMH ∗
94, WBF ∗ ∗ ∗
14] and for vector quan-tities [ER03, JSD04, Mon05, Pri12, WKBB18]. In the following wewill discuss this approach in more detail.In order to approximate the Laplacian of the velocity field wecombine an SPH derivative with a finite difference derivative whichyields the following equation [Mon05]: ∇ v i = ( d + ) ∑ j m j ρ j v i j · x i j (cid:107) x i j (cid:107) + . h ∇ W i j , (102)where x i j = x i − x j , v i j = v i − v j and d is the number of spatialdimensions. Note that a term 0 . h is introduced in the denom-inator to avoid singularities. The approximation of the Laplacianin Eq. (102) has some nice features. First, it is Galilean invariant.Moreover, it vanishes for rigid body rotation. This is an importantproperty since there is no friction if all particles rotate uniformly.Finally, the introduced formulation conserves linear and angularmomentum [Mon92].Instead of computing the Laplacian of the velocity field in or-der to simulate viscosity, in some works XSPH is used as artificialviscosity model ( e.g. , [SB12]). XSPH determines the smoothed ve-locity ˆ v i of a particle i asˆ v i = v i + α ∑ j m j ρ j ( v j − v i ) W i j , (103)where 0 ≤ α < c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Figure 14: Different examples for the simulation of viscous behavior. Left: The simulation of highly viscous fluids enables realistic bucklingeffects [WKBB18]. Center: A viscous dough interacts with a fast moving solid [PICT15]. Right: A melting Eiffel tower is simulated (themodel is courtesy of Pranav Panchal) [PT16].neighborhood. The advantage of this formulation is that no kernelderivative is required. The disadvantage is that α is not physicallymeaningful. Simulating the behavior of highly viscous materials implies that theviscosity coefficient is large. However, in this case explicit viscos-ity solvers tend to get unstable. Therefore, it is recommended to usean implicit method in order to simulate highly viscous fluids. In thissubsection we introduce some of the most important implicit vis-cosity solvers [TDF ∗ Takahashi et al. [TDF ∗ As discussed in the previous subsec-tion, one way to compute the second derivative of the velocity fieldis to take two first SPH derivatives. This approach is used by Taka-hashi et al. to formulate an implicit integration scheme for the vis-cosity term in the Navier-Stokes equations. The implicit integrationenables a stable simulation of highly viscous fluids. In each simula-tion step Takahashi et al. first determine the strain rate E i for eachparticle i using Eq. (98). Following the Navier-Stokes equations(see Eq. (99)) the authors then compute the divergence of the strainrate as ∇ · (cid:16) ∇ v i + ( ∇ v i ) T (cid:17) = ∑ j m j (cid:32) E i ρ i + E j ρ j (cid:33) ∇ W i j . (104)Using this formulation the implicit integration scheme can be de-rived as v ( t + ∆ t ) = v ∗ + ∆ t ρ µ ∇ · (cid:16) ∇ v ( t + ∆ t ) + ( ∇ v ( t + ∆ t )) T (cid:17) , (105)where v ∗ is the predicted velocity which is determined by integrat-ing all non-pressure forces except viscosity. Takahashi et al. substi-tute Eq. (104) in Eq. (105) and solve the resulting formula to getthe new velocities of the particles. The advantage of this implicitscheme is that highly viscous fluids can be simulated in a stableway while the viscosity is independent of the temporal and spatialresolution. However, in this formulation all second-ring neighborsof a particle have to be considered in order to compute one first-order SPH derivative after the other. This leads to many non-zeroelements in the system matrix which decreases the performancesignificantly. Peer et al. [PICT15, PT16]
Instead of using a classical implicittime integration scheme, Peer et al. propose to decompose andmodify the velocity gradient ∇ v . The goal of the authors is tomodify only the shear rate in order to simulate a viscous behav-ior. Hence, they exploit the fact that the velocity gradient can bedecomposed as ∇ v = R + V + S , (106)where R = ( ∇ v − ( ∇ v ) T ) is the spin rate tensor, V = ( ∇ · v ) (cid:49) the expansion rate tensor and S = E − V the traceless shear ratetensor. This decomposition enables to modify the traceless shearrate tensor without influencing the other components of the velocitygradient. Therefore, the authors define a target velocity gradient ∇ v target = R + V + ξ S (107)which reduces the shear rate by a user-defined factor 0 ≤ ξ ≤ v i ( t + ∆ t ) = ρ i ∑ j m j (cid:32) v j ( t + ∆ t ) + ∇ v target i + ∇ v target j x i j (cid:33) W i j . (108)This yields a linear system Av ( t + ∆ t ) = b , where the matrix entriesand the right hand side vector are defined as A i j = − m j W i j , (109) A ii = ρ i − m i W ii , (110) b i = ∑ j m j ∇ v target i + ∇ v target j x i j W i j . (111)This system can be decomposed to get three smaller linear sys-tems for the x-, y- and z-component of the velocity. Finally, theauthors propose to solve the three systems using a conjugate gradi-ent method.Later, Peer and Teschner [PT16] extended this method by sim-ulating vorticity diffusion in order to improve the rotational mo-tion. The diffusion process in a viscous fluid is described by D ω Dt = ν ∇ ω . The authors determine ω = ( ω x , ω y , ω z ) T from the spin ratetensor as R = − ω z ω y ω z − ω x − ω y ω x . (112) c (cid:13) (cid:13)(cid:13)
Instead of using a classical implicittime integration scheme, Peer et al. propose to decompose andmodify the velocity gradient ∇ v . The goal of the authors is tomodify only the shear rate in order to simulate a viscous behav-ior. Hence, they exploit the fact that the velocity gradient can bedecomposed as ∇ v = R + V + S , (106)where R = ( ∇ v − ( ∇ v ) T ) is the spin rate tensor, V = ( ∇ · v ) (cid:49) the expansion rate tensor and S = E − V the traceless shear ratetensor. This decomposition enables to modify the traceless shearrate tensor without influencing the other components of the velocitygradient. Therefore, the authors define a target velocity gradient ∇ v target = R + V + ξ S (107)which reduces the shear rate by a user-defined factor 0 ≤ ξ ≤ v i ( t + ∆ t ) = ρ i ∑ j m j (cid:32) v j ( t + ∆ t ) + ∇ v target i + ∇ v target j x i j (cid:33) W i j . (108)This yields a linear system Av ( t + ∆ t ) = b , where the matrix entriesand the right hand side vector are defined as A i j = − m j W i j , (109) A ii = ρ i − m i W ii , (110) b i = ∑ j m j ∇ v target i + ∇ v target j x i j W i j . (111)This system can be decomposed to get three smaller linear sys-tems for the x-, y- and z-component of the velocity. Finally, theauthors propose to solve the three systems using a conjugate gradi-ent method.Later, Peer and Teschner [PT16] extended this method by sim-ulating vorticity diffusion in order to improve the rotational mo-tion. The diffusion process in a viscous fluid is described by D ω Dt = ν ∇ ω . The authors determine ω = ( ω x , ω y , ω z ) T from the spin ratetensor as R = − ω z ω y ω z − ω x − ω y ω x . (112) c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Analogous to Eq. (107) the vorticity is reduced by solving the sys-tem ∇ ω target i = ξ ∇ ω . (113)The resulting vector ω target i is used to determine a target spin ratetensor R target i which is substituted in Eq. (107) before reconstruct-ing the velocity field using Eq. (108).The proposed methods are very efficient and enable a stable sim-ulation of highly viscous materials. However, these methods havealso some disadvantages. The reconstruction of the velocity fieldusing SPH is problematic as discussed in [BGFAO17] and intro-duces a significant damping. When simulating highly viscous flu-ids, this damping effect is not that crucial but this approach is notrecommended for the simulation of low viscous flow. Another dis-advantage of the methods is that the viscosity parameter ξ is notphysically meaningful and depends on the temporal and spatial res-olution. Bender and Koschier [BK17]
The authors of this work alsoreduce the strain rate by introducing a user-defined coefficientwhich is similar to the core idea of Peer et al. However, insteadof modifying the velocity gradient and reconstructing the veloc-ity field, Bender and Koschier define a velocity constraint function C i ( v ) = E i − γ E i with the user-defined coefficient 0 ≤ γ ≤
1. Theconstraint is defined as six-dimensional vector function where thevector contains the elements of the upper triangular part of the sym-metric strain rate tensor. Finally, the constraint is enforced by firstsolving the linear system (cid:32) ρ i ∂ E i ∂ v i (cid:18) ∂ E i ∂ v i (cid:19) T + ∑ j ρ i ∂ E i ∂ v j (cid:18) ∂ E i ∂ v j (cid:19) T (cid:33) µ i = E i − γ E i (114)for the Lagrange multiplier µ by Jacobi iterations. The final veloci-ties are then determined as v i ( t + ∆ t ) = v ∗ i + m i (cid:32) m i ρ i (cid:18) ∂ E i ∂ v i (cid:19) T µ i + ∑ j m j ρ j (cid:18) ∂ E j ∂ v i (cid:19) T µ j (cid:33) . (115)Details about the computation of ∂ E / ∂ v can be found in [BK17].The advantage of solving a constraint function instead of us-ing the velocity field reconstruction approach of Peer et al. is thatalso low viscous fluids can be simulated. The disadvantages of themethod are that solving six-dimensional constraints using Jacobi it-erations is computationally expensive and the introduced viscositycoefficient depends on the temporal and spatial resolution. Weiler et al. [WKBB18]
All implicit viscosity methods intro-duced so far, use a formulation based on the strain rate tensor E .The strain rate is determined by Eq. (98), where the velocity gra-dient ∇ v is computed using the following SPH discretization (seeSection 2.5): ∇ v i = ρ i ∑ j m j ( v j − v i ) ∇ W Ti j . (116)Weiler et al. found out that this SPH discretization is negativelyaffected by the particle deficiency problem at the free surface of afluid. For a rotational velocity field (see Fig. 15, left) the strain rate ms 1s Figure 15: When computing the strain rate tensor usingEq. (116), errors occur at the free surface due to particle defi-ciency [WKBB18]. Left: Velocity field of a rotational motion.Right: The corresponding strain rate should be zero for all parti-cles. However, the plot of the Frobenius norm of the tensors showsan error at the free surface.should be E = since a rotation is a rigid body motion which doesnot deform the body. However, when using the SPH discretiza-tion in Eq. (116), the strain rate is not zero at the free surface (seeFig. 15, right). In this experiment we can observe a significant errorat the boundary. The main problem is that the viscosity solver triesto counteract this erroneous strain rate which leads to ghost forces.These forces causes severe visual artifacts and a loss of angularmomentum which is discussed later in more detail.To solve this problem, Weiler et al. developed an implicit viscos-ity solver which directly determines the Laplacian of the velocityfield instead of using the strain rate. Their approach is based on theimplicit integration scheme v ( t + ∆ t ) = v ∗ + ∆ t ρ µ ∇ v ( t + ∆ t ) . (117)This is similar to the one of Takahashi et al. [TDF ∗
15] but uses theLaplacian of the velocity field instead of the divergence of the strainrate. To compute the Laplacian, the approximation in Eq. (102) isused. Since this approximation vanishes for rigid body rotationsand conserves linear and angular momentum [Mon92], the pro-posed approach solves the problems of the methods above.Eq. (117) is a linear system which has to be solved to get theunknown new velocities v ( t + ∆ t ) . Using the SPH discretization ofthe Laplacian in Eq. (102), we can rewrite this linear system as ( (cid:49) − ∆ t A ) v ( t + ∆ t ) = v ∗ , (118)where the matrix A contains a 3 × A i j for each pair ofneighboring particles i and j : A i j = − ( d + ) µ m i j ρ i ρ j ∇ W i j x Ti j (cid:13)(cid:13) x i j (cid:13)(cid:13) + . h , A ii = − ∑ j A i j . (119)Note that the average mass m i j = . ( m i + m j ) is used in orderto obtain a symmetric system. The resulting system can be solvedefficiently by a matrix-free conjugate gradient method. The con-vergence can be improved by a block Jacobi preconditioner wherethe preconditioner matrix is block diagonal with the 3 × (cid:49) − ∆ t A ii . Moreover, starting the conjugate gradient solver with an c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids initial guess of v ∗ + ∆ v using the velocity difference of the last step ∆ v = v ( t ) − v ∗ ( t − ∆ t ) further improves the performance.Weiler et al. propose to extend this viscosity formulation alsofor the boundary in order to simulate materials that stick to solidobjects. In SPH often a particle-based surface representation of theboundary is used. For such a surface representation the diagonalmatrix blocks in Eq. (119) and the right hand side of the system inEq. (118) have to be adapted as A ii = − ∑ j A i j + ( d + ) ∑ k µ b m k ρ i ∇ W ik x Tik (cid:107) x ik (cid:107) + . h , (120) b i = v ∗ i − ( d + ) ∆ t ∑ k µ b m k ρ i v k · x ik (cid:107) x ik (cid:107) + . h ∇ W ik , (121)where m k is the mass of the boundary particle k (see Section 5).Note that a reaction force has to be applied to the boundary particlesto get a consistent two-way coupling [AIA ∗ Comparison
In this part we want to compare all implicit viscositysolvers introduced above. All approaches except the one of Weileret al. [WKBB18] are based on a strain rate formulation. As dis-cussed above and shown in Fig. 15, this leads to errors at the freesurface due to particle deficiency. In practice this can lead to ar-tifacts at the surface (see Figs. 16a, 16b and 16c). Moreover, thestrain rate error at the free surface causes a significant loss of an-gular momentum which leads to a damped rotational motion (seeFigs. 17a, 17b and 17c).Weiler et al. analyzed these problems and proposed a newmethod which computes the Laplacian of the velocity field insteadof using the divergence of the strain rate. Moreover, they use anSPH approximation of the Laplacian that vanishes for rigid bodyrotation and conserves linear and angular momentum. In this waythe problems at the free surface can be solved (see Figs. 16d and17d).While the methods of Peer et al. [PICT15, PT16] and Benderand Koschier [BK17] use a viscosity parameter that depends onthe temporal and spatial resolution, Takahashi et al. [TDF ∗
15] andWeiler et al. [WKBB18] solve this problem by using a consistentimplicit time integration.Finally, when comparing the performance, the approach of Peeret al. [PICT15] is the fastest method. This is due to the fact thatthey can decompose their linear system in three smaller ones whilethis cannot be done for the approach of Weiler et al. Bender andKoschier use a Jacobi solver which converges slower and Takahashiet al. have to consider the second-ring neighbors which results inlarge computational overhead.
For the simulation of low viscous fluids an explicit viscosity for-mulation should be used since explicit methods are computation-ally less expensive. We recommend to compute the viscous forcein Eq. (100) by approximating the Laplacian of the velocity fieldusing Eq. (102). An alternative, which is computationally less ex-pensive but also less accurate, is to use XSPH as artificial viscosity. An implicit viscosity solver is recommended for the simula-tion of highly viscous fluids due to stability reasons. As discussedabove, the implicit strain rate based formulations suffer from an er-ror in the SPH approximation that causes visual artifacts and leadsto a loss of angular momentum. The method of Weiler et al. avoidsthe problem and therefore generates more realistic results. Finally,we think that it would be an interesting open problem for future re-search to find a better SPH approximation of the strain rate withoutproblems at the free surface. Such an approximation would solvethe problems of the strain rate based formulations.Note that all viscosity methods that were discussed in thissection are implemented in our open-source framework SPlisH-SPlasH [Ben19b].
7. Surface Tension
Surface tension is an important physical phenomenon which is aubiquitous effect in daily life. For example, surface tension forceskeep liquid molecules together when pouring water into a glass.The surface forces are the result of intermolecular attractive forcesat microscopic scales. The molecules attract each other inside ofa fluid while the molecules at the surface are pulled inwards.Therefore, surface tension minimizes the surface area which causesdroplets of water to form a sphere when external forces are ex-cluded. We typically speak of cohesion if molecules of the sametype attract each other while adhesion describes the attractive forcesbetween molecules of different types. Cohesion and adhesion areimportant effects when simulating surface tension.In recent years, various methods were proposed to simulate sur-face tension effects in SPH-based fluid simulations ( e.g. , [BT07,AAT13, HWZ ∗ Surface tension is the result of attracting forces between molecules.Methods that are based on a microscopic point of view aim to simu-late the intermolecular cohesive forces. However, since the smallestelement in an SPH simulation is a particle, these forces are deter-mined at the particle-level.Becker and Teschner [BT07] propose to compute the particle ac-celeration due to cohesion byD v i D t = − α m i ∑ j m j (cid:0) x i − x j (cid:1) W i j , (122)where α is a coefficient to control the surface tension of the fluid.This equation interpolates the position differences in the neighbor-hood of a particle i and computes an acceleration to attract theneighboring particles. Note that adhesion can also be simulated byEq. (122) when using a sum over all neighboring boundary parti-cles. c (cid:13) (cid:13)(cid:13)
Surface tension is an important physical phenomenon which is aubiquitous effect in daily life. For example, surface tension forceskeep liquid molecules together when pouring water into a glass.The surface forces are the result of intermolecular attractive forcesat microscopic scales. The molecules attract each other inside ofa fluid while the molecules at the surface are pulled inwards.Therefore, surface tension minimizes the surface area which causesdroplets of water to form a sphere when external forces are ex-cluded. We typically speak of cohesion if molecules of the sametype attract each other while adhesion describes the attractive forcesbetween molecules of different types. Cohesion and adhesion areimportant effects when simulating surface tension.In recent years, various methods were proposed to simulate sur-face tension effects in SPH-based fluid simulations ( e.g. , [BT07,AAT13, HWZ ∗ Surface tension is the result of attracting forces between molecules.Methods that are based on a microscopic point of view aim to simu-late the intermolecular cohesive forces. However, since the smallestelement in an SPH simulation is a particle, these forces are deter-mined at the particle-level.Becker and Teschner [BT07] propose to compute the particle ac-celeration due to cohesion byD v i D t = − α m i ∑ j m j (cid:0) x i − x j (cid:1) W i j , (122)where α is a coefficient to control the surface tension of the fluid.This equation interpolates the position differences in the neighbor-hood of a particle i and computes an acceleration to attract theneighboring particles. Note that adhesion can also be simulated byEq. (122) when using a sum over all neighboring boundary parti-cles. c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids (a) Takahashi et al. [TDF ∗
15] (b) Peer & Teschner [PT16] (c) Bender & Koschier [BK17] (d) Weiler et al. [WKBB18]
Figure 16: The SPH approximation of the strain rate tensor is erroneous at the free surface due to particle deficiency which leads to artifacts(a,b,c). The formulation of Weiler et al. avoids this strain rate computation by using the Laplacian of the velocity field which solves thisproblem (d) [WKBB18]. (a) Takahashi et al. [TDF ∗
15] (b) Peer & Teschner [PT16] (c) Bender & Koschier [BK17] (d) Weiler et al. [WKBB18]
Figure 17: When simulating the rope coiling effect with the introduced implicit viscosity solvers, we can observe a loss of angular momentumcaused by errors in the SPH strain rate approximation (a,b,c). Using a formulation based on the Laplacian of the velocity field instead, solvesthis problem and enables a uniform coiling (d) [WKBB18].
Akinci et al. [AAT13] present a method based on amacroscopic point of view. Instead of only considering cohesiveforces, also forces to minimize the surface area are determined.Their method computes the cohesive force of a particle as F cohesion i ← j = − α m i m j x i − x j (cid:107) x i − x j (cid:107) W cohesion (cid:0) (cid:107) x i − x j (cid:107) (cid:1) , (123)where i and j are neighboring particles and W cohesion is a specialcohesion kernel which is defined by W cohesion ( r ) = π h ( h − r ) r r > h ∧ r ≤ h ( h − r ) r − h r > ∧ r ≤ h n i = h ∑ j m j ρ j ∇ W i j (125)yields a surface normal pointing into the fluid. Note that the factor h is used to make the normal scale independent. The magnitude ofthe resulting vector is close to zero in the interior of the fluid andproportional to the curvature at the free surface. Hence, a symmet-ric force that counteracts the curvature can be defined as F curvature i ← j = − α m i ( n i − n j ) . (126)This force is used to minimize the surface area.Finally, both forces are combined as F surface tension i ← j = K i j ( F cohesion i + F curvature i ) , (127)where K i j = ρ ρ i + ρ j is a symmetric factor that amplifies the surfacetension forces at the free surface. At the surface ρ i and ρ j are under-estimated due to particle deficiency and K i j > K i j ≈ Adhesion
To simulate the attractive forces between fluid particlesand the boundary, an adhesion force is introduced as F adhesion i ← k = − β m i m k W adhesion ik x i − x k (cid:107) x i − x k (cid:107) , (128) c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Figure 18: Different surface tension effects are realized by approximating intermolecular attractive forces in the fluid [AAT13].where β is the adhesion coefficient and k denotes a neighboringboundary particle. The computation of the mass m k of a boundaryparticle is described in detail in Section 5. Akinci et al. propose touse another specialized kernel function for the computation of theadhesion forces which is defined as W adhesion ( r ) = . h . (cid:113) − r h + r − h r > h ∧ r ≤ h h / h areattracted by the boundary.
8. Vorticity
One of the most visually appealing phenomena in dynamic fluidsis the evolution of chaotic structures due to turbulences. Turbulentmotions are largely caused by the interaction of many unsteady vor-tices on various scales. A vortex is, moreover, defined as a localspinning motion in the fluid – mathematically spoken the vorticityis a vector field ω = ∇ × v . (130)In SPH fluid simulations turbulent details quickly get lost dueto numerical diffusion [dGWH ∗
15] or due to coarse sampling ofthe velocity field [IOS ∗
14, CIPT14] which negatively influencesthe visual liveliness of the flow. In order to facilitate the forma-tion of vortices in the simulation and to counteract numerical dif-fusion, a range of approaches has been proposed in the past. Mostof these approaches originate from research concerning Eulerian,grid-based discretizations and can roughly be categorized into vor-ticity confinement techniques, Lagrangian vortex methods, fluidup-sampling, and more recently micropolar models. In this sectionwe will discuss a simple vorticity confinement approach followingMacklin and Müller [MM13] and an SPH discretization of a mi-cropolar model that facilitates the formation of vortices as proposedby Bender et al. [BKKW18].
As already discussed, SPH discretizations tend to overly dissipateenergy in turbulent flow. Therefore, Macklin and Müller [MM13]employ a method based on vorticity confinement in order to coun-teract the dissipation by amplifying existing vortices. The tech-nique consists of three steps: 1. The vorticity ω i for each particle i is computed using a discretecurl operator, e.g. , ω i = ∇ × v i = − ∑ j m j ρ j v i j × ∇ i W i j . (131)2. A corrective force is computed and applied that amplifies thealready existing vortical motion, i.e. , F vorticity i = ε vorticity ( η (cid:107) η (cid:107) × ω i ) (132) η = ∑ j m j ρ j (cid:107) ω j (cid:107)∇ i W i j , (133)where ε vorticity denotes a small constant used to steer the amountof amplification.3. The velocity field is smoothed using XSPH (see Eq. (103)) inorder to enforce a coherent particle motion.Algorithmically, this correction is applied just before the advectionof the SPH particle positions.While this approach is very simple and effectively amplifies ex-isting vortices it has some drawbacks. It is hard to choose thecontrol parameter ε vorticity such that overamplification is avoided.Moreover, the method can, in the best case, only conserve existingvortices but does not facilitate the formation of new ones. The most prominent mathematical model describing the dynam-ics of Newtonian fluids is the Navier-Stokes model. As describedin Section 2 the model can be derived from the conservation lawof linear momentum (see Eq. (30)) and by presuming that the me-chanical stress is composed of isotropic pressure and a diffusingviscous term. An important assumption of the model is that the in-finitesimally small particles which compose a fluid continuum arenot subject to rotational motion. This also implies that the law ofangular momentum conservation is identically fulfilled if and onlyif the stress tensor is symmetric.In this section we introduce the concept of micropolar fluids andpresent a material model that generalizes the Navier-Stokes equa-tions for the simulation of incompressible, inviscid turbulent flowas proposed by Bender at al. [BKKW18]. Following the defini-tion of Łukaszewicz [Łuk99], a micropolar fluid follows constitu-tive laws modeled using a generally non-symmetric stress tensor. c (cid:13) (cid:13)(cid:13)
As already discussed, SPH discretizations tend to overly dissipateenergy in turbulent flow. Therefore, Macklin and Müller [MM13]employ a method based on vorticity confinement in order to coun-teract the dissipation by amplifying existing vortices. The tech-nique consists of three steps: 1. The vorticity ω i for each particle i is computed using a discretecurl operator, e.g. , ω i = ∇ × v i = − ∑ j m j ρ j v i j × ∇ i W i j . (131)2. A corrective force is computed and applied that amplifies thealready existing vortical motion, i.e. , F vorticity i = ε vorticity ( η (cid:107) η (cid:107) × ω i ) (132) η = ∑ j m j ρ j (cid:107) ω j (cid:107)∇ i W i j , (133)where ε vorticity denotes a small constant used to steer the amountof amplification.3. The velocity field is smoothed using XSPH (see Eq. (103)) inorder to enforce a coherent particle motion.Algorithmically, this correction is applied just before the advectionof the SPH particle positions.While this approach is very simple and effectively amplifies ex-isting vortices it has some drawbacks. It is hard to choose thecontrol parameter ε vorticity such that overamplification is avoided.Moreover, the method can, in the best case, only conserve existingvortices but does not facilitate the formation of new ones. The most prominent mathematical model describing the dynam-ics of Newtonian fluids is the Navier-Stokes model. As describedin Section 2 the model can be derived from the conservation lawof linear momentum (see Eq. (30)) and by presuming that the me-chanical stress is composed of isotropic pressure and a diffusingviscous term. An important assumption of the model is that the in-finitesimally small particles which compose a fluid continuum arenot subject to rotational motion. This also implies that the law ofangular momentum conservation is identically fulfilled if and onlyif the stress tensor is symmetric.In this section we introduce the concept of micropolar fluids andpresent a material model that generalizes the Navier-Stokes equa-tions for the simulation of incompressible, inviscid turbulent flowas proposed by Bender at al. [BKKW18]. Following the defini-tion of Łukaszewicz [Łuk99], a micropolar fluid follows constitu-tive laws modeled using a generally non-symmetric stress tensor. c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Moreover, the definition includes that the fluid consists of rigid,spherical (and therefore rotationally invariant) particles. Based onthe non-symmetric stress measures, the micropolar model addition-ally models rotating motions of the infinitesimal spherical parti-cles using an angular velocity field. Due to the additional rotationaldegrees of freedom, the generation of vortices is facilitated and awider range of potential dynamic effects are captured by the model.Please note, that for this section we will neglect any dissipationterms, such as viscosity, as the main goal is to generate undamped,highly turbulent flows. For the complete model, we would like tokindly refer the reader to the original paper [BKKW18].
Balance Law for Angular Momentum Conservation
Similar to the conservation law of linear momentum (see Eq. (30))a balance law for angular momentum can be derived, i.e. , ρΘ D ω Dt = T × + τ , (134)with [ T × ] i = ∑ j ∑ k ε i jk T jk , where ε i jk and τ denote the Levi-Civitatensor and external body torque. Further, Θ represents a scalar,isotropic microinertia coefficient. A physical interpretation for thisquantity is that each infinitesimal fluid particle has a certain inertialresistance against rotational accelerations. Bender et al. suggest tochoose Θ = s − based on experimentation. We would furtherlike to stress the fact that Θ is not at all related to the spatial extentsof an SPH particle as it is defined in the continuous setting. Constitutive Model
It is essential to understand that the classical model actually also re-spects angular momentum conservation (see Eq. (134)). It is in thiscontext, however, rarely explicitly mentioned as the balance law isusually identically fulfilled based on two assumptions. Firstly, theclassical approach does not model external torques, i.e. , τ ≡ . Sec-ondly, stress tensor T is usually chosen as a symmetric tensor andsuch that T × ≡ D ω Dt ≡
0. For this reason, the balancelaw of angular momentum is not particularly useful for symmetricstress measures.In order to account for the microstructured particles and to utilizethe balance law of angular momentum, Bender et al. propose to usethe following constitutive relation: T = − p (cid:49) − µ t ∇ v + µ t ω × , (135)with (cid:2) ω × (cid:3) = ∑ i ε jik ω i , where µ t denotes the "transfer coefficient".We will later discuss a physical interpretation of the term in thefinal PDE that is controlled using µ t . In order to ensure consistencywith the second law of thermodynamics µ t ≥ Equations of Motion
As the conservation laws and constitutive equation are now estab-lished, we can finally derive the augmented equations of motionthat build the basis for the numerical simulation. By plugging theconstitutive relation (135) into the conservation laws (30) and (134) and by applying the incompressibility condition (28), we arrive atthe following representation: D v Dt = − ρ ∇ p + ν t ∇ × ω + f ext ρ (136) Θ D ω Dt = ν t ( ∇ × v − ω ) + τρ , (137)where ν t denotes the kinematic transfer coefficient. Looking atEq. (136), we quickly notice that it is identical to the inviscidNavier-Stokes equation (Euler equation) but augmented by the term ν t ∇ × ω . A complementary term also governed by the transfer co-efficient resides in Eq. (137), i.e. , ν t ( ∇ × v − ω ) . The terms ef-fectively convert angular accelerations into linear accelerations andvice versa. Physically, they can be interpreted as dissipation-freefriction or as a dissipation-free viscosity coupling linear and rota-tional motion.In order to realize the model in the implementation, it is requiredto discretize Eqs. (136) and (137). This means that additional to adiscrete representation of x and v it is necessary to discretize thevorticity ω . Please note, that due to the assumption in the micropo-lar model, that the microstructure of the material particles is spher-ical we do not have to discretize and track the rotational field (onlythe angular velocity field) which makes the implementation lesscomplicated and delivers better performance. The transfer forcesand torques can then, following the splitting approach, simply beapplied in line with the non-pressure forces and integrated explic-itly as they are considerably less stiff than the pressure forces. It isfurther advised to "filter" the resulting velocity field using XSPH(see Eq. (103)) in order to ensure coherent particle motion. Algo-rithm 7 shows an exemplary pseudocode of the resulting method. for all particle i do Find neighbors for all particle i do Compute density ρ i for all particle i do Compute non-pressure forces F non-pressure i Compute transfer forces F transfer i = m i ν t ∇ × ω i Compute transfer torque τ transfer i = m i ν t ( ∇ × v i − ω i ) Compute time step size ∆ t according to CFL for all particle i dov ∗ i = v i + ∆ tm i ( F non-pressure i + F transfer i + F ext i ) for all particle i do Enforce incompressibility using pressure solverUpdate v ∗ i for all particle i dov i ( t + ∆ t ) = v ∗ i x i ( t + ∆ t ) = x i + ∆ t v i ( t + ∆ t ) ω i ( t + ∆ t ) = ω i ( t ) + ∆ tm i Θ i ( τ transfer i + τ ext i ) Algorithm 7: Simulation loop for SPH simulation turbulent microp-olar fluids.Bender et al. [BKKW18] demonstrated the effect of the transfercoefficient using an intuitive example (see Fig. 19). In the experi-ment a fluid flowing in a narrow channel was simulated. Moreover, c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Figure 19: Simulation of 4 .
7M turbulent fluid particles with threeobstacles and increasing transfer coefficient ν t . Top-down: ν t = . s − , ν t = . s − , ν t = . s − .Figure 20: 1M fluid particles interact with a fast rotating propellerresulting in highly turbulent flow.three obstacles were placed in the channel to provoke turbulenceswhile the transfer coefficient was continuously increased. In thetop image we can see that that the flow is only moderately turbu-lent for a transfer coefficient ν t = . s − . For larger values thevorticity significantly increases (middle) and even tends to get un-realistic for values greater than 0 . s − (bottom). Furthermore,they showcase the visual realism that can be achieved in turbulentscenarios (see Fig. 20). Discussion
Two methods to improve the behavior of the simulation in the pres-ence of turbulences have been explained. In this paragraph, wewould like to discuss the similarities and differences between vor-ticity confinement and the micropolar model.Both methods build on the concept of obtaining/maintaining avorticity field (angular velocity field) ω following Eq. (130). How- ever, the main idea of vorticity confinement is to merely identifyand amplify existing vortices. Moreover, the vorticity will alwaysbe derived from the linear field. In contrast, the micropolar ap-proach builds on the concept of angular momentum conservationand on modeling a constitutive model for turbulences. In this moresophisticated setting, the velocity field and the vorticity are dis-cretized independently and strongly coupled via the transfer termsin Eqs. (136) and (137). In this way there is a complex interactionbetween both physical quantities that not only conserves existingvortices better but also facilitates the formation of new vortices.This effect can be exemplified using the lid-driven cavity exper-iment carried out by Bender et al. [BKKW18] (see Fig. 21). Inthis experiment the "lid" (top-side) of a two-dimensional domainfilled with water is accelerated with constant velocity. Given suit-able model parameters, the velocity field is expected to stabilize ina big central vortex and three minor vortices rotating in the oppositedirection. This result shows that vorticity confinement successfullyamplifies the vortical motion but is not able to form the additionalcorner vortices. In contrast, the micropolar approach yields the ex-pected result.
9. Multiphase Fluids
The simulation of multiple immiscible and miscible fluids greatlyenhance visual effects in graphics. In contrast to Eulerian ap-proaches, the particle representation of SPH offers the advantagethat fluid interfaces are sharply defined. In this section, we firstpresent how the standard SPH equations can be adapted to modeldensity discontinuity across fluid interfaces, and introduce the re-sulting adapted force equations. We then discuss models for cap-turing complex mixing phenomena.
A simple approach to simulate multiple fluids with SPH is to as-sign different labels to particles of different phases, and assign-ing them with corresponding physical attributes such as massesand rest densities [MSKG05]. Typically, each particle’s rest vol-ume remains constant to ensure a uniform particle sampling, thus m a ρ a = m b ρ b for two fluid types a and b . The momentum equation canbe solved with the single flow SPH formulation presented in theprevious sections, while simply using the physical attributes storedon the particles. However, for high density ratios between phases,this can lead to instability problems that are not time step related.The desired density discontinuity across the interface is smootheddue to the nature of SPH of summing up contributions from particleneighbors. As a consequence, pressure and force fields are affected,which manifests as spurious interface tension [Hoo98, AMS ∗ δ i = ∑ j W i j was introduced and thestandard SPH equations were adapted accordingly [TM05, HA06,SP08]. The density of a particle i is then computed as˜ ρ i = m i δ i . (138) c (cid:13) (cid:13)(cid:13)
A simple approach to simulate multiple fluids with SPH is to as-sign different labels to particles of different phases, and assign-ing them with corresponding physical attributes such as massesand rest densities [MSKG05]. Typically, each particle’s rest vol-ume remains constant to ensure a uniform particle sampling, thus m a ρ a = m b ρ b for two fluid types a and b . The momentum equation canbe solved with the single flow SPH formulation presented in theprevious sections, while simply using the physical attributes storedon the particles. However, for high density ratios between phases,this can lead to instability problems that are not time step related.The desired density discontinuity across the interface is smootheddue to the nature of SPH of summing up contributions from particleneighbors. As a consequence, pressure and force fields are affected,which manifests as spurious interface tension [Hoo98, AMS ∗ δ i = ∑ j W i j was introduced and thestandard SPH equations were adapted accordingly [TM05, HA06,SP08]. The density of a particle i is then computed as˜ ρ i = m i δ i . (138) c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Figure 21: Velocity fields of the lid-driven cavity benchmark. Standard SPH (left) and vorticity confinement (middle) are only able toproduce one large central vortex. In contrast, the micropolar approach yields (right) the expected result, i.e. , one central vortex and threesmaller vortices in the corners which are rotating in the opposite direction.Like this, the density of particle i is not influenced by the mass ofits neighbors j , while still receiving the geometric contribution W i j from j . The state equation of Sections 4.4 can then be changed suchthat the pressure is computed with the adapted density as˜ p i = k (cid:32)(cid:18) ˜ ρ i ρ (cid:19) k − (cid:33) . (139)Solenthaler et al. [SP08] derived adapted forces by substituting˜ ρ i and ˜ p i into the Navier-Stokes equations and applying the SPHformalism. The resulting pressure force term is then given as F p = − ∇ ˜ p δ . By employing the quotient rule we then get ∇ ˜ p δ = ∇ ( ˜ p δ ) + ˜ p δ ∇ δ . After applying the SPH rule and replacing V by ρ ,the pressure force equation can be written as F pi = − ∑ j (cid:32) ˜ p j δ j + ˜ p i δ i (cid:33) ∇ W i j . (140)Similar derivations can be found in [TM05, HA06]. The viscosityforce (and other force terms) can be derived analogously and isgiven as F vi = δ i ∑ j µ i + µ j δ j ( v j − v i ) ∇ W i j . (141)Note that the above equations are identical to the standard SPHequations when applied to a single phase flow. For multiple flu-ids, however, the adapted method eliminates any spurious tensioneffects and notably increases stability. The method has been ex-tended with an incompressibility condition and solid-fluid cou-pling [AIA ∗
12, GPB ∗ ∗ ∗ Fluid mixing can be simulated by solving the diffusion equation ∂ C ∂ t = α ∇ C , which evolves the concentration C over time. WithSPH, this equation can be written as [MSKG05] ∂ C i ∂ t = α ∑ j m j C j − C i ρ j ∇ W i j , (142)where α defines the diffusion strength. Another SPH formulationfor computing the diffusion has been presented in [LLP11].More complex mixing effects can be simulated by taking the c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids flow motion and force distributions into account as demonstratedin the SPH-based mixture model of Ren et al. [RLY ∗ ρ m D t = ∂ρ m ∂ t + ∇ · ( ρ m v m ) = , (143)where ρ m is the rest density of the mixture and v m is the mixturevelocity, averaged over all phases. ρ m and v m are computed usingthe volume fraction α k of a phase k with rest density ρ k , i.e. , ρ m = ∑ k α k ρ k and v m = ρ m ∑ k α k ρ k v k . The momentum equation for themixture is given asD ( ρ m , v m ) D t = −∇ p + ∇ · ( τ m + τ Dm ) + ρ m g , (144)where τ m and τ Dm are the mixture’s viscous stress and diffusiontensors, respectively.In each simulation step, the drift velocity v mk = v k − v m is com-puted, which represents the relative velocity of phase k to the mix-ture. The equation can be rewritten using individual terms for slipvelocity due to body forces, pressure effects that cause fluid phasesto move from high to low pressure regions, and a Brownian diffu-sion term that represents phase drifting from high to low concentra-tion regions. The drift velocity is then used to calculate the diffu-sion tensor τ Dm and change in volume fraction D α k / D t . The SPHequations for the mixture model described above can be found inthe work of Ren et al. [RLY ∗ ∗
16] extended the mixture model to handle the in-teraction between fluid and solid phases, and demonstrated variouseffects including dissolution of solids, flows in porous media, andinteraction with elastic materials. Another extension has been pre-sented by Yang et al. [YCR ∗
15] where an energy-based model wasused. The approach integrates the Cahn-Hilliard equation that de-scribes phase separation, expanding the capability of a multi-fluidsolver and enabling incompressible flows.
10. Deformable Solids
The simulation of deformable solids is an active research topic incomputer graphics. The most popular simulation approaches in thisarea, like the finite element method (FEM) [KBT17, KKB18] andPosition-Based Dynamics (PBD) [BKCW14, BMM14, BMM17],are mesh-based. However, also meshless methods were investigatedlike the moving least squares (MLS) method [AW09]. In this sec-tion we show that SPH is also an interesting meshless method tosimulate deformable solids. An advantage of an SPH-based simu-lation of deformables is that this enables a simple coupling betweenfluids and solids in a unified framework.
In this subsection we first introduce a continuum mechanical for-mulation for linear elasticity. In the next subsection we then showhow to discretize the resulting equations using SPH. The deformation of a solid is defined by the function Φ ( X ) = X + u = x (145)which maps a point X in the reference configuration to its currentposition x in the deformed configuration, where u = x − X is thedisplacement vector. Differentiating this function with respect tothe reference position X gives us the deformation gradient J = ∂Φ ( X ) ∂ X = ∂ x ∂ X = (cid:49) + ∂ u ∂ X . (146)This quantity can be used to measure the strain of a deformed body.In computer graphics often a linear strain measure is used to avoidthe solution of a non-linear system of equations. For the same rea-son we introduce the linear infinitesimal strain tensor ε ( J ) = ( J + J T ) − (cid:49) . (147)The next step is to define a constitutive model for linear elastic-ity. We follow the work of Sifakis [Sif12] and define it in terms ofthe strain energy density: Ψ ( J ) = µ ε : ε + λ ( ε ) , (148)where µ and λ are the Lamé coefficients [Sif12]. The first Pi-ola–Kirchhoff stress tensor is determined by differentiating thestrain energy density with respect to the deformation gradient P ( J ) = ∂Ψ∂ J . (149)For our linear elasticity model this yields P ( J ) = µ ε + λ tr ( ε ) (cid:49) . (150)Note that the Lamé coefficients can be computed from Young’smodulus k and Poisson’s ratio ν by µ = k ( + ν ) , λ = k ν ( + ν )( − ν ) . (151)This is often more intuitive since Young’s modulus is a measureof stretch resistance while Poisson ratio is a measure of incom-pressibility. Finally, the elastic body forces are determined as thedivergence of the stress tensor f = ∇ · P . (152)In the following subsection we will discuss how this continuouselastic material model can be discretized using the SPH formula-tion. The deformation of an elastic body is determined with respect to itsreference configuration (see Eq. (145)). Typically the initial shapeof a body is used as reference configuration in an SPH simulation.Since the topology of an elastic body does not change during thesimulation, we store the neighborhood of each particle i in the ref-erence configuration. In the following this reference neighborhoodis denoted by N i . c (cid:13) (cid:13)(cid:13)
In this subsection we first introduce a continuum mechanical for-mulation for linear elasticity. In the next subsection we then showhow to discretize the resulting equations using SPH. The deformation of a solid is defined by the function Φ ( X ) = X + u = x (145)which maps a point X in the reference configuration to its currentposition x in the deformed configuration, where u = x − X is thedisplacement vector. Differentiating this function with respect tothe reference position X gives us the deformation gradient J = ∂Φ ( X ) ∂ X = ∂ x ∂ X = (cid:49) + ∂ u ∂ X . (146)This quantity can be used to measure the strain of a deformed body.In computer graphics often a linear strain measure is used to avoidthe solution of a non-linear system of equations. For the same rea-son we introduce the linear infinitesimal strain tensor ε ( J ) = ( J + J T ) − (cid:49) . (147)The next step is to define a constitutive model for linear elastic-ity. We follow the work of Sifakis [Sif12] and define it in terms ofthe strain energy density: Ψ ( J ) = µ ε : ε + λ ( ε ) , (148)where µ and λ are the Lamé coefficients [Sif12]. The first Pi-ola–Kirchhoff stress tensor is determined by differentiating thestrain energy density with respect to the deformation gradient P ( J ) = ∂Ψ∂ J . (149)For our linear elasticity model this yields P ( J ) = µ ε + λ tr ( ε ) (cid:49) . (150)Note that the Lamé coefficients can be computed from Young’smodulus k and Poisson’s ratio ν by µ = k ( + ν ) , λ = k ν ( + ν )( − ν ) . (151)This is often more intuitive since Young’s modulus is a measureof stretch resistance while Poisson ratio is a measure of incom-pressibility. Finally, the elastic body forces are determined as thedivergence of the stress tensor f = ∇ · P . (152)In the following subsection we will discuss how this continuouselastic material model can be discretized using the SPH formula-tion. The deformation of an elastic body is determined with respect to itsreference configuration (see Eq. (145)). Typically the initial shapeof a body is used as reference configuration in an SPH simulation.Since the topology of an elastic body does not change during thesimulation, we store the neighborhood of each particle i in the ref-erence configuration. In the following this reference neighborhoodis denoted by N i . c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Figure 23: Deformable solids simulated using an SPH formulation [PGBT17]. Left: An elastic bunny model is hit by a rigid capsule. Right:A deformable hairball interacts with water.
Deformation Gradient
A straightforward SPH discretization ofthe deformation gradient is given by J i ( x , X ) = ∑ j ∈N i V j x ji ∇ W (cid:0) X i j (cid:1) T , (153)where x ji = x j − x i and X i j = X i − X j . Since this SPH approxima-tion is determined in the reference configuration, the rest volume V j has to be used in the sum. However, this formulation fails tocapture rotational motion since it is not first-order consistent (seeSection 2.3). Kernel Gradient Correction
Bonet and Lok [BL99] have shownthat the gradient of the kernel has to fulfill the following conditionto ensure that the computation is first-order consistent and thereforecorrectly captures rotational motion: ∑ j ∈N i V j X ji ∇ W (cid:0) X i j (cid:1) T = (cid:49) . (154)Now we can formulate a corrected kernel gradient which satisfiesthe condition by construction:˜ ∇ W i (cid:0) X i j (cid:1) = L i ∇ W (cid:0) X i j (cid:1) , (155)where L i is a correction matrix that is defined as L i = ∑ j ∈N i V j ∇ W (cid:0) X i j (cid:1) X Ti j − . (156)Note that this correction matrix only depends on the rest volumeand the particle positions in the reference configuration. Therefore,this matrix is precomputed at the beginning of the simulation. If thematrix in Eq. (156) is singular and cannot be inverted, e.g. , due toa collinear or coplanar particle configuration, the Moore–Penroseinverse is used instead. Corotated Approach
Using the corrected kernel gradient we geta first-order consistent SPH formulation for the deformation gradi-ent: J i ( x , X ) = ∑ j ∈N i V j x ji ˜ ∇ W (cid:0) X i j (cid:1) T . (157) The deformation gradient is used to compute the linear infinites-imal strain tensor by Eq. (147). Note that we use a linear strainmeasure since in an implicit formulation it is more efficient tosolve a linear system than a non-linear one. However, the linearstrain tensor is not invariant under rotations. In computer graphicsa common solution for this problem is to use a corotational ap-proach [BIT09, KKB18]. The core idea of this approach is to ex-tract the rotation and to compute the strain measure in an unrotatedframe. In the following we will show how the rotation can be ex-tracted and introduce the computation of a corotated deformationgradient.The deformation gradient computed with the corrected kernelgradient (see Eq. (157)) is able to capture the rotation correctly.Hence, the per-particle rotation R i can be directly extracted from J i , e.g. , by using the efficient and stable method of Müller etal. [MBCM16].The extracted rotation matrix R is used to rotate the referenceconfiguration so that the resulting displacement vector u = x − RX contains no rotation. Since we rotate the reference configuration,we also have to rotate the corrected kernel gradient as it dependson the reference positions. This yields the rotated corrected kernelgradient ∗ ∇ W i (cid:0) X i j (cid:1) = R i L i ∇ W (cid:0) X i j (cid:1) . (158)Putting all together gives us the corotated deformation gradient J ∗ i ( x , X ) = (cid:49) + ∑ j ∈N i V j (cid:0) x ji − R i X ji (cid:1) ∗ ∇ W i (cid:0) X i j (cid:1) T . (159)Now we compute the strain tensor ε ( J ∗ ) using Eq. (147) and thestress tensor P ( J ∗ ) using Eq. (150). Finally, the force is determinedas the divergence of the stress tensor. In our SPH formulation thisyields [Gan15]: F i (cid:0) J ∗ (cid:1) = ∑ j ∈N i V i V j (cid:18) P i (cid:0) J ∗ i (cid:1) ∗ ∇ W i (cid:0) X i j (cid:1) − P j (cid:0) J ∗ j (cid:1) ∗ ∇ W j (cid:0) X i j (cid:1)(cid:19) . (160) c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids (a) reference configuration (b) deformed configuration Figure 24: Due to numerical errors the transformed reference vector J i X ji and the current deformed vector x ji are typically not equal inan SPH simulation. The error is defined by the difference e iji = J i X ji − x ji .If we simply add these particle forces to our system, we get anexplicit approach for the simulation of deformable bodies. How-ever, this approach is only conditionally stable and requires smalltime steps when simulating stiff solids. To improve the stability wewill introduce an implicit approach in the next subsection. The implicit method described in the following is based on the workof Peer et al. [PGBT17]. Since the elastic forces depend linearly onthe particle positions, an implicit formulation of the time step isstraightforward v t + ∆ t = v t + ∆ tm F (cid:16) J t + ∆ t (cid:17) , (161)where J t + ∆ t = J ∗ (cid:16) x t + ∆ t , X (cid:17) is the deformation gradient at the endof the time step. This means that we use the new particle positions x t + ∆ t to determine the elastic forces at the end of the time step.In this formulation we have unknown positions x t + ∆ t and un-known velocities v t + ∆ t . In the next step we substitute the positionsby x t + ∆ t = x t + ∆ t v t + ∆ t to get a linear system for the new velocities.Moreover, we split the computation of the force into F (cid:16) J t + ∆ t (cid:17) = F (cid:0) J t (cid:1) + F (cid:16) J ∆ t (cid:17) , (162)where J t = J ∗ (cid:0) x t , X (cid:1) is the deformation gradient at the beginningof the time step and J ∆ t = J ∗ (cid:16) ∆ t v t + ∆ t , (cid:17) is the deformation gra-dient that corresponds to the position change in one time step dueto the velocities v t + ∆ t . Now we can bring all terms that depend onthe unknown new velocities to the left hand side and the rest to theright hand side: v t + ∆ t − ∆ tm F (cid:16) J ∗ (cid:16) ∆ t v t + ∆ t , (cid:17)(cid:17) = v t + ∆ tm F (cid:0) J ∗ (cid:0) x t , X (cid:1)(cid:1) . (163)This yields a linear system for the velocities v t + ∆ t which can ef-ficiently be solved using a matrix-free conjugate gradient method.More details about the matrix-free solver can be found in the workof Peer et al. [PGBT17]. In SPH simulations zero-energy modes can occur which are sim-ilar to hour glass modes in finite element methods [Gan15]. Thedisplacement field in the neighborhood of a particle i has to be de-fined exactly by its corresponding deformation gradient J i . Thismeans that if we transform the vector X ji from particle i to oneof its neighbors j in reference space using the deformation gradi-ent J i X ji , this should give us the actual vector x ji in the deformedconfiguration. Hence, the vector e iji = J i X ji − x ji (164)should be zero. However, this is typically not the case due to nu-merical errors (see Fig. 24).Ganzenmüller [Gan15] proposes to compute a penalty force tominimize the error vector e ji : F HG i = − α k ∑ j ∈N i V i V j W (cid:0) X i j (cid:1) (cid:107) X i j (cid:107) (cid:32) e iji · x ji (cid:107) x ji (cid:107) + e ji j · x i j (cid:107) x i j (cid:107) (cid:33) x ji (cid:107) x ji (cid:107) , (165)where the coefficient α controls the amplitude of the zero-energymode suppression and k is the Young’s modulus. In this way thesystem gets more stable and hourglass modes are suppressed.
11. Rigid Solids
Recently, it has been demonstrated that SPH can also be used torealize a rigid body simulation with contact handling [GPB ∗ In the following we will discuss how to compute rigid-rigid contactforces F rr . We want to compute these forces similar to the fluid-rigid interface forces which were discussed in Section 5. Therefore,we first introduce an artificial rest density ρ r = r . Note that the magnitude of the rest density can be chosenarbitrarily since we are only interested in a density deviation. Ifthere is a contact, we get a density deviation of ρ r − ρ r > c (cid:13) (cid:13)(cid:13)
Recently, it has been demonstrated that SPH can also be used torealize a rigid body simulation with contact handling [GPB ∗ In the following we will discuss how to compute rigid-rigid contactforces F rr . We want to compute these forces similar to the fluid-rigid interface forces which were discussed in Section 5. Therefore,we first introduce an artificial rest density ρ r = r . Note that the magnitude of the rest density can be chosenarbitrarily since we are only interested in a density deviation. Ifthere is a contact, we get a density deviation of ρ r − ρ r > c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Figure 25: The SPH-based rigid body solver enables a strong two-way coupling between fluids and rigid bodies. In this scenario 43.8M fluidparticles interact with 50M static and 2.3M dynamic rigid body particles. Up to 90k simultaneous rigid-rigid contacts were resolved.Figure 26: Two creatures run in highly viscous mud, break through a wall of rigid bodies and collide with a deformable tree. This complexsimulation shows the power of a unified SPH solver.the corresponding particle r . In this case our goal is to determinecontact forces F rr such that ρ r = ρ r .Now we will derive an implicit method to compute the unknowncontact forces. We start with the continuity equation D ρ r Dt = − ρ r ∇ · v r , (166)where ρ r and v r are the density and the velocity of a rigid bodyparticle r , respectively. Then we use a backward difference timediscretization and introduce a constant density constraint ρ t + ∆ tr = ρ r to get ρ r − ρ r ∆ t = − ρ r ∇ · v t + ∆ tr , (167)where v t + ∆ tr is the velocity of the rigid body particle r at time t + ∆ t . This velocity vector can be written as v t + ∆ tr = v t + ∆ tR + ω t + ∆ tR × r t + ∆ tr , (168)where v t + ∆ tR and ω t + ∆ tR are the linear and angular velocity of therigid body R at time t + ∆ t , respectively, and r t + ∆ tr is the vector fromthe center of mass of the rigid body to the position of the particle r .The new velocities are determined by an Euler integration step v t + ∆ tR = v R + ∆ t m R (cid:32) F R + ∑ k ∈R F rr k (cid:33) (169) ω t + ∆ tR = ω R + ∆ t I − R (cid:32) τ R + ( I R ω R ) × ω R + ∑ k ∈R r k × F rr k (cid:33) , (170)where I R is the inertia tensor of the rigid body. The vectors F R c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids and τ R contain all forces and torques acting on the body exceptthe unknown rigid-rigid contact forces F rr k . R denotes the set of allparticles of rigid body R . Note that all quantities on the right handside are at time t . For improved readability we omitted the timeparameter for all quantities at the current time t .In the next step we substitute Eqs. (168)-(170) in Eq. (167) toget a linear system for the unknown rigid-rigid contact forces ρ r − ρ r ∆ t = − ρ r ∇ · (cid:18) v R + ∆ tm R F R (cid:19) − ρ r ∇ · (cid:32) ∆ tm R ∑ k ∈R F rr k (cid:33) − ρ r ∇ · (cid:16)(cid:16) ω R + ∆ t I − R ( τ R + ( I R ω R ) × ω R ) (cid:17) × r t + ∆ tr (cid:17) − ρ r ∇ · (cid:32) ∆ t (cid:32) I − R ∑ k ∈R r k × F rr k (cid:33) × r t + ∆ tr (cid:33) . (171)To simplify this system we use the approximation r t + ∆ tr = r r . More-over, we introduce the velocity vector v sr = v R + ∆ tm R F R + (cid:16) ω R + ∆ t I − R ( τ R + ( I R ω R ) × ω R ) (cid:17) × r r . (172)This vector determines the new velocity of a particle r after a timestep which considers all forces and torques except the unknowncontact forces. In this way we can write the right-hand side of ourlinear system in a compact form: s r = ρ r − ρ r ∆ t + ρ r ∇ · v sr . (173)The left-hand side contains all terms of Eq. (171) that depend onthe rigid-rigid contact forces F rr k . The resulting linear system hasthe form − ρ r ∇ · (cid:32) ∆ tm R ∑ k ∈R F rr k + (cid:32) ∆ t I − R ∑ k ∈R r k × F rr k (cid:33) × r k (cid:33) = s r . (174)The left-hand side can further be simplified by introducing the ma-trix K rk = m R (cid:49) − ˜ r r I − R ˜ r k , (175)where ˜ r r is the cross product matrix of r r to get − ρ r ∇ · (cid:32) ∆ t ∑ k ∈R K rk F rr k (cid:33) = s r . (176)Note that the matrix K rk is well-known in the area of rigid bodysolvers [Mir96, BET14].Our goal is to resolve the contacts by a pressure force. Therefore,we define F rr k = − V k ∇ p k , (177)where V k is an artificial volume of particle k and p k is an unknownpressure which is used to resolve the collision. This yields the finallinear system ρ r ∇ · (cid:32) ∆ t ∑ k ∈R V k K rk ∇ p k (cid:33) = ρ r − ρ r ∆ t + ρ r ∇ · v sr . (178)Solving the linear system gives us the unknown pressure values for all rigid particles. Note that the linear system contains one equa-tion for each rigid particle. However, if a particle r has no contactto a particle of another rigid body, we can remove the correspond-ing equation from the system and set p r = In the following we describe how the quantities in the derived linearsystem are computed.The artificial rest volume of a rigid particle r is determined as V r = . ∑ k ∈R W rk . (179)A detailed discussion about this computation can be found in thework of Gissler et al. [GPB ∗ ρ r = ρ r = ∑ k V r ρ r W rk , where the sum considers the particles k of allrigid bodies within the support radius of the kernel. Due to the sumover the particles of neighboring rigid bodies, we get a density de-viation of ρ r − ρ r > r as V r = ρ r V r ρ r .The divergence on the right-hand side of the linear system isdetermined as ∇ · v sr = ρ r ∑ k ∈R V k ρ k (cid:0) v sk − v sr (cid:1) · ∇ W rk . (180)Finally, we solve the linear system using a relaxed Jacobi solverand update the pressure in iteration l + p l + r = p lr + β RJ r b r (cid:32) s r − ρ r ∇ · (cid:32) ∆ t ∑ k ∈R V k K rk ∇ p lk (cid:33)(cid:33) , (181)where b r is the diagonal element of the linear system and β RJ r is therelaxation coefficient which is set to β RJ r = . . The introduced SPH-based rigid body solver enables a strong cou-pling between fluids and rigid bodies. As shown in Fig. 27 thesolver is able to accurately handle complex scenarios with thou-sands of simultaneous contacts. It can be easily extended to sim-ulate friction effects [GPB ∗
12. Data Driven Fluid Simulation
Using machine learning for fluid simulations is a largely unex-plored research area, but first results are promising and indicate thepotential of such data-driven approaches. In the Lagrangian con-text, the seminal work of Ladický et al. [LJS ∗
15] employed Re-gression forests to infer particle accelerations (and velocities) us-ing handcrafted, SPH-based features. We discuss this work in moredetail below. Um et al. [UHT18] presented a method to augment c (cid:13) (cid:13)(cid:13)
15] employed Re-gression forests to infer particle accelerations (and velocities) us-ing handcrafted, SPH-based features. We discuss this work in moredetail below. Um et al. [UHT18] presented a method to augment c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Figure 27: Closeup of the water gate scene in Fig. 25. The singlechain elements are connected only by rigid-rigid contacts.simulations with learned splashes from a high-resolution FLIP sim-ulation, but included also an example where SPH training data wasused. Somewhat related to SPH, Schenk et al. [SF18] proposed adifferentiable PBF solver [MM13] for deep neural networks. Theyhave presented convolution layers for summing up contributionsfrom neighbors and for fluid-object interaction, which potentiallycan be adapted to SPH fluids as well. Other work mainly focusedon Eulerian simulations, for example to substitute the pressureprojection step with a CNN [TSSP16], to synthesize flow simu-lations from a set of reduced parameters [KAT ∗ In the following, we give an overview of the regression forest ap-proach for SPH presented by Ladický et al. [LJS ∗ t + ∆ t is efficiently estimatedgiven the state at time t .As input to the regressor, a feature vector Φ x i is evaluated foreach particle. The features are designed such that they representthe individual forces and constraints of the Navier-Stokes equa-tions: the used features model pressure, incompressibility, viscos-ity and surface tension. In order to evaluate features without using Figure 28: The regression forest approach for SPH [LJS ∗
15] en-ables real-time simulations of over a million particles. At runtime,additional external forces can be added to mimic different materialproperties without retraining the model.an explicit neighbor search step, context-based integral features arecomputed that are defined as flat-kernel sums of rectangular regionssurrounding a particle. The different box sizes allow to capture thebehavior of both close and distant particles, and the features canbe evaluated in constant time and are robust to small deviations.More details on the computation of context-based integral featurescan be found in [LJS ∗
1. Learning naïve prediction:
The first approach directly learnsparticle accelerations a , given the evaluated features at state t . Theregression problem is formulated as a i ( t ) : = Reg ( Φ x i ) , (182)where Reg ( . ) is the learned regression function. Velocities and po-sitions are then integrated with v i ( t + ∆ t ) = v i ( t ) + a i ( t ) ∆ t (183) x i ( t + ∆ t ) = x i ( t ) + v i ( t ) + v i ( t + ∆ t ) ∆ t . (184)This formulation mimics standard SPH and hence does not considerincompressibility.
2. Learning prediction with hindsight:
The second strategy firstcomputes external forces, advects particles, and applies collisionhandling. Then, this intermediate state with particle positions x ∗ i isused to compute integral features. The regression learns a corrective c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids acceleration and is defined as a i ( t ) : = Reg ( Φ x ∗ i ) , (185)followed by advection as in Eqs. (183) and (184). Unlike the naïveapproach, the regressor is able to predict compressions and hence tocounteract those with a corrective acceleration. Conceptually, thisapproach mimics PCISPH [SP09].
3. Learning correction:
The third approach starts similarly as thesecond one, but instead of learning accelerations, corrective veloc-ities are computed. The regression problem is defined as ∆ v corri : = Reg ( Φ x ∗ i ) , (186)and positions and velocities are updated with v i ( t + ∆ t ) = v ∗ i ( t ) + ∆ v corri (187) x i ( t + ∆ t ) = x ∗ i ( t ) + ∆ v corri ∆ t . (188)This approach counteracts compressions as well, and conceptuallymimics PBF [MM13]. Unlike PBF, the regressor takes into accountinformation from a larger neighborhood, and hence does not requireseveral iterations to converge.For training the regression forest, 165 scenes consisting of 1-6 million particles and moving obstacles (sphere, box, cylinder)were randomly generated and computed for 6 seconds. The trainingtime was 4 days on 12 CPUs, and the size of the resulting modelwas about 40 MB. With the regression fluid approach it is possibleto simulate 1 to 1.5 million particles in real-time, and hence thisapproach represents an attractive alternative to traditional solversfor games and virtual reality applications (Fig. 28).With the naïve prediction, strong compression artifacts are vis-ible. The system cannot self-correct in the next frames since themodel has never seen such distorted states during the training. Bothprediction with hindsight and learning corrections can handle in-compressibility well, however the third approach seems to lead tosmaller errors compared to the ground truth data. Additionally, withthe second and third approaches, external forces can be added with-out retraining the model. This allows adding surface tension, fric-tion, or drag effects at runtime to mimic different material prop-erties as illustrated in Fig. 28. The disadvantage of the regressionfluids approach - and in fact of all machine learning based strate-gies - is that learning methods are not capable to extrapolate themodel far outside the training data ( e.g. , when domain size or fluidresolution change).
13. SPlisHSPlasH
In this section we want to introduce SPlisHSPlasH [Ben19b] whichis an SPH-based open-source library for the physically-based sim-ulation of fluids (see Fig. 29). The SPlisHSPlasH framework con-tains a reference implementation of many of the methods intro-duced in this tutorial and several simulations shown in the figureswere performed using this library. Therefore, we think that ouropen-source framework perfectly supplements these course notes.In the current version SPlisHSPlasH implements six of the mostpopular explicit and implicit pressure solvers [BT07, SP09, MM13, ICS ∗
14, BK15, WKB16] which enable the simulation of incom-pressible fluids with several million particles. Moreover, the ex-plicit [SB12, Mon92] and implicit viscosity methods [TDF ∗ ∗ ∗ ∗
12] to simulate the coupling between rigid bodies andfluids. The required surface sampling of the bodies is performedautomatically using a Poisson disk sampling. For the simulationof dynamic rigid bodies SPlisHSPlasH uses the open-source Po-sitionBasedDynamics library [Ben19a]. This library simulates therigid bodies using a position-based approach [DCB14]. Collisionsbetween the bodies are efficiently detected using signed distancefields [KDBB17] while the contacts are resolved using a projectedGauss-Seidel method [BET14]. Finally, SPlisHSPlasH implementsdifferent methods for the simulation of deformable solids [BIT09,PGBT17] using an SPH formulation (see Section 10). Since an SPHformulation is used, the two-way coupling between solids and flu-ids is simply handled by the implemented multiphase method.SPlisHSPlasH uses a neighborhood search based on the com-pact hashing approach of Ihmsen et al. [IABT11]. This approach isdiscussed in more detail in Section 3. The neighborhood search isimplemented in our open-source library CompactNSearch [Kos19].The SPlisHSPlasH framework has many more features like emit-ters, adaptive time-stepping or the support of different kernel func-tions. Moreover, the library has some useful tools like volume sam-pling of closed geometries or the export of particle data for Mayaor Houdini. New simulation scenarios can be created easily using aJSON-based scene file format. Finally, due to a modular concept itis simple to extend the library and to integrate own SPH methods.Therefore, we think that SPlisHSPlasH is a good starting point forall beginners in the area of SPH-based simulations.
14. Conclusion
This tutorial introduced state-of-the-art SPH techniques for thephysics based simulation of fluids and solids in graphics and pre-sented practical guidelines for implementations. Various conceptsthat are particularly relevant for graphics applications were dis-cussed. We showed that with the recent improvements SPH mod-els have matured and ultimately emerged as a competitive alter-native to Eulerian fluid simulations or hybrid approaches. Particu-lar challenges of SPH concerning neighborhood search algorithms,pressure solvers, or versatile fluid-solid interaction techniques havebeen overcome. With the improved robustness and efficiency, mil-lions of particles can today be simulated on a single desktop com- c (cid:13) (cid:13)(cid:13)
This tutorial introduced state-of-the-art SPH techniques for thephysics based simulation of fluids and solids in graphics and pre-sented practical guidelines for implementations. Various conceptsthat are particularly relevant for graphics applications were dis-cussed. We showed that with the recent improvements SPH mod-els have matured and ultimately emerged as a competitive alter-native to Eulerian fluid simulations or hybrid approaches. Particu-lar challenges of SPH concerning neighborhood search algorithms,pressure solvers, or versatile fluid-solid interaction techniques havebeen overcome. With the improved robustness and efficiency, mil-lions of particles can today be simulated on a single desktop com- c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids Figure 29: Screenshots of the SPlisHSPlasH fluid simulation framework. Left: Two-way coupling of dynamic rigid bodies and a color-codedfluid. Center: Two-phase double dam break simulation. Right: Highly viscous bunny collides with a static dragon model.puter. Accordingly, the Lagrangian SPH method reaches an un-precedented level of visual quality, where fine-scale surface effectsand flow details are reliably captured.The SPH community – in graphics as well as in other researchdisciplines – is very active and the field advances quickly. Eachcommunity contributes to different aspects of SPH simulations, andthe research often finds applications across disciplines. For graph-ics applications, it was especially important to efficiently enforceincompressibility on unstructured particles and hence to eliminatethe severe time step restrictions of standard SPH techniques. Wehave presented a practical introduction to various SPH conceptsthat enforce volume conservation and/or divergence-free velocityfields. A current difficulty is that these approaches render time step-ping more challenging, since the largest possible time step does notnecessarily result in the best overall performance. Future work iscertainly necessary to establish a CFL condition for these methods,as well as to overcome the current main performance bottleneckwhich is still the time step restriction especially when using mil-lions of particles.Using large particle numbers is one of the key componentsfor high visual quality and production level results. Such high-resolution simulations pose new challenges and existing conceptsmight need to be revisited. Especially speed, flexibility and con-trollability are core aspects, for which solutions are still largelymissing. This problem, however, affects not only the SPH field butthe entire fluid community in graphics likewise, and has triggeredresearch on pre- and post-processing methods or data-driven ap-proaches.Our tutorial introduced the SPH-based open-source librarySPlisHSPlasH that contains reference implementations of manyconcepts that we discussed. This implementation is an excellentstarting point for students, researchers and practitioners, and mayserve as a valuable tool for future research.
References [AAT13] A
KINCI
N., A
KINCI
G., T
ESCHNER
M.: Versatile surface ten-sion and adhesion for SPH fluids.
ACM Transactions on Graphics 32 , 6(2013), 1–8. 25, 26, 27, 37[Abe12] A
BEYARATNE
R.:
Continuum Mechanics: Volume II of LectureNotes on The Mechanics of Elastic Solids . techreport, MIT Departmentof Mechanical Engineering, 2012. 7 [AIA ∗
12] A
KINCI
N., I
HMSEN
M., A
KINCI
G., S
OLENTHALER
B.,T
ESCHNER
M.: Versatile rigid-fluid coupling for incompressible SPH.
ACM Transactions on Graphics 31 , 4 (July 2012), 1–8. 18, 20, 21, 25,30, 37[AMS ∗
07] A
GERTZ
O., M
OORE
B., S
TADEL
J., P
OTTER
D., M
INIATI
F., R
EAD
J., M
AYER
L., G
AWRYSZCZAK
A., K
RAVTSOV
A., M ON - AGHAN
J., N
ORDLUND
A., P
EARCE
F., Q
UILIS
V., R
UDD
D.,S
PRINGEL
V., S
TONE
J., T
ASKER
E., T
EYSSIER
R., W
ADSLEY
J.,W
ALDER
R.: Fundamental differences between SPH and grid methods.
Mon. Not. R. Astron. Soc. 380 , 3 (2007), 963–978. 29[AW09] A
DAMS
B., W
ICKE
M.: Meshless Approximation Methods andApplications in Physics Based Modeling and Animation. In
Proceedingsof the Eurographics conference (2009), EG ’09, Eurographics Associa-tion, pp. 213–239. 31[Ben19a] B
ENDER
J.: PositionBasedDynamics Library. https://github.com/InteractiveComputerGraphics/PositionBasedDynamics , 2019. 37[Ben19b] B
ENDER
J.: SPlisHSPlasH Library. https://github.com/InteractiveComputerGraphics/SPlisHSPlasH , 2019. 25, 37[BET14] B
ENDER
J., E
RLEBEN
K., T
RINKLE
J.: Interactive Simulationof Rigid Body Dynamics in Computer Graphics.
Computer GraphicsForum 33 , 1 (2014), 246–270. 35, 37[BGFAO17] B
ARREIRO
H., G
ARCÍA -F ERNÁNDEZ
I., A
LDUÁN
I.,O
TADUY
M. A.: Conformation constraints for efficient viscoelastic fluidsimulation.
ACM Transactions on Graphics 36 , 6 (2017), 221.1–221.11.24[BGI ∗
18] B
AND
S., G
ISSLER
C., I
HMSEN
M., C
ORNELIS
J., P
EER
A.,T
ESCHNER
M.: Pressure boundaries for implicit incompressible sph.
ACM Transactions on Graphics 37 , 2 (Feb. 2018), 14:1–14:11. 18[BGPT18] B
AND
S., G
ISSLER
C., P
EER
A., T
ESCHNER
M.: MLS pres-sure boundaries for divergence-free and viscous SPH fluids.
Computers& Graphics 76 (nov 2018), 37–46. 17, 18[BIT09] B
ECKER
M., I
HMSEN
M., T
ESCHNER
M.: Corotated SPH fordeformable solids. In
Proceedings of Eurographics Conference on Nat-ural Phenomena (2009), pp. 27–34. 32, 37[BK15] B
ENDER
J., K
OSCHIER
D.: Divergence-Free Smoothed Parti-cle Hydrodynamics. In
ACM SIGGRAPH/Eurographics Symposium onComputer Animation (2015), pp. 1–9. 17, 37[BK17] B
ENDER
J., K
OSCHIER
D.: Divergence-Free SPH for Incom-pressible and Viscous Fluids.
IEEE Transactions on Visualization andComputer Graphics 23 , 3 (2017), 1193–1206. 9, 17, 18, 24, 25, 26, 37[BKCW14] B
ENDER
J., K
OSCHIER
D., C
HARRIER
P., W
EBER
D.:Position-Based Simulation of Continuous Materials.
Computers &Graphics 44 , 1 (2014), 1–10. 31[BKKW17] B
ENDER
J., K
OSCHIER
D., K
UGELSTADT
T., W
EILER
M.: c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids A micropolar material model for turbulent sph fluids. In
ACM SIG-GRAPH/Eurographics Symposium on Computer Animation (July 2017),pp. 1–8. 37[BKKW18] B
ENDER
J., K
OSCHIER
D., K
UGELSTADT
T., W
EILER
M.:Turbulent micropolar sph fluids with foam.
IEEE Transactions on Visu-alization and Computer Graphics (2018). 27, 28, 29[BL99] B
ONET
J., L OK T.-S.: Variational and momentum preservationaspects of smooth particle hydrodynamic formulations.
Computer Meth-ods in Applied Mechanics and Engineering 180 , 1 (1999), 97 – 115. 32[BLS12] B
ODIN
K., L
ACOURSIÈRE
C., S
ERVIN
M.: Constraint fluids.
IEEE Transactions on Visualization and Computer Graphics 18 (2012),516–526. 16, 18, 19[BMM14] B
ENDER
J., M
ÜLLER
M., M
ACKLIN
M.: A Survey onPosition-Based Simulation Methods in Computer Graphics.
ComputerGraphics Forum 33 , 6 (2014), 228–251. 31[BMM17] B
ENDER
J., M
ÜLLER
M., M
ACKLIN
M.: A survey on posi-tion based dynamics, 2017. In
EUROGRAPHICS 2017 Tutorials (2017),Eurographics Association. 31[Bri15] B
RIDSON
R.:
Fluid Simulation for Computer Graphics, SecondEdition . Taylor & Francis, 2015. 8[Bro85] B
ROOKSHAW
L.: A method of calculating radiative heat diffu-sion in particle simulations.
Publications of the Astronomical Society ofAustralia 6 , 2 (1985), 207–210. 7, 22[BT07] B
ECKER
M., T
ESCHNER
M.: Weakly compressible SPH for freesurface flows. In
ACM SIGGRAPH/Eurographics Symposium on Com-puter Animation (2007), pp. 1–8. 18, 25, 26, 37[CBG ∗
18] C
ORNELIS
J., B
ENDER
J., G
ISSLER
C., I
HMSEN
M.,T
ESCHNER
M.: An optimized source term formulation for incompress-ible SPH.
The Visual Computer (Feb. 2018). 16[CIPT14] C
ORNELIS
J., I
HMSEN
M., P
EER
A., T
ESCHNER
M.: IISPH-FLIP for incompressible fluids.
Computer Graphics Forum 33 , 2 (may2014), 255–262. 27[CM99] C
LEARY
P. W., M
ONAGHAN
J. J.: Conduction modelling usingsmoothed particle hydrodynamics.
Journal of Computational Physics148 , 1 (1999), 227 – 264. 22[Com16a] C
OMPUTER A NIMATION , RWTH A
ACHEN U NIVERSITY :Divergence-free sph for incompressible and viscous fluids. , 2016. 1[Com16b] C
OMPUTER G RAPHICS , U
NIVERSITY OF F REIBURG : Ter-rain 2 - up to 500 million particles with PreonLab (FIFTY2). , 2016. 1[Com17] C
OMPUTER G RAPHICS , U
NIVERSITY OF F REIBURG : Shipunder attack. ,2017. 1[Com18] C
OMPUTER G RAPHICS , U
NIVERSITY OF F REIBURG : An im-plicit SPH formulation for incompressible linearly elastic solids. , 2018. 1[DCB14] D
EUL
C., C
HARRIER
P., B
ENDER
J.: Position-based rigid-body dynamics.
Computer Animation and Virtual Worlds 27 , 2 (Sept.2014), 103–112. 37[dGWH ∗ DE G OES
F., W
ALLEZ
C., H
UANG
J., P
AVLOV
D., D ES - BRUN
M.: Power Particles: An incompressible fluid solver based onpower diagrams.
ACM Transactions on Graphics 34 , 4 (2015), 50:1–50:11. 27[ER03] E
SPANOL
P., R
EVENGA
M.: Smoothed dissipative particle dy-namics.
Physical Review E 67 , 2 (2003), 026705. 22[FAW17] F
ÜRSTENAU
J.-P., A
VCI
B., W
RIGGERS
P.: A comparativenumerical study of pressure-Poisson-equation discretization strategiesfor SPH. In (2017). 16[FIF16] FIFTY2 T
ECHNOLOGY : PreonLab Promotion. , 2016. 1 [FLR ∗
13] F
ERRAND
M., L
AURENCE
D. R., R
OGERS
B. D., V
IOLEAU
D., K
ASSIOTIS
C.: Unified semi-analytical wall boundary conditionsfor inviscid, laminar or turbulent flows in the meshless SPH method.
International Journal for Numerical Methods in Fluids 71 , 4 (Feb. 2013),446–472. 19[FM15] F
UJISAWA
M., M
IURA
K. T.: An Efficient Boundary Handlingwith a Modified Density Calculation for SPH.
Computer Graphics Fo-rum 34 , 7 (2015), 155–162. 19[FMH ∗
94] F
LEBBE
O., M
UENZEL
S., H
EROLD
H., R
IFFERT
H.,R
UDER
H.: Smoothed Particle Hydrodynamics: Physical viscosity andthe simulation of accretion disks.
The Astrophysical Journal 431 (Aug.1994), 754–760. 22[GAC ∗
09] G
RENIER
N., A
NTUONO
M., C
OLAGROSSI
A., T
OUZÉ
D. L., A
LESSANDRINI
B.: An Hamiltonian interface SPH formulationfor multi-fluid and free surface flows.
Journal of Computational Physics228 , 22 (2009), 8380 – 8393. 30[Gan15] G
ANZENMÜLLER
G. C.: An hourglass control algorithm for la-grangian smooth particle hydrodynamics.
Computer Methods in AppliedMechanics and Engineering 286 (apr 2015), 87–106. 32, 33[GBP ∗
17] G
ISSLER
C., B
AND
S., P
EER
A., I
HMSEN
M., T
ESCHNER
M.: Generalized drag force for particle-based simulations.
Computers& Graphics 69 (dec 2017), 1–11. 37[GM77] G
INGOLD R. A ., M ONAGHAN
J.: Smoothed Particle Hydrody-namics: Theory and Application to Non-Spherical Stars.
Monthly No-tices of the Royal Astronomical Society , 181 (1977), 375–389. 4[GPB ∗
19] G
ISSLER
C., P
EER
A., B
AND
S., B
ENDER
J., T
ESCHNER
M.: Interlinked sph pressure solvers for strong fluid-rigid coupling.
ACMTransactions on Graphics 38 , 1 (Jan. 2019), 5:1–5:13. 18, 30, 33, 35[HA06] H U X., A
DAMS
N.: A multi-phase SPH method for macroscopicand mesoscopic flows.
Journal of Computational Physics 213 , 2 (2006),844–861. 29, 30[HKK07a] H
ARADA
T., K
OSHIZUKA
S., K
AWAGUCHI
Y.: Smoothedparticle hydrodynamics in complex shapes. In
Spring Conference onComputer Graphics (2007), pp. 191–197. 18, 19[HKK07b] H
ARADA
T., K
OSHIZUKA
S., K
AWAGUCHI
Y.: SmoothedParticle Hydrodynamics on GPUs. In
Computer Graphics International (2007), pp. 63–70. 18, 19[HLL ∗
12] H E X., L IU N., L I S., W
ANG
H., W
ANG
G.: Local PoissonSPH for Viscous Incompressible Fluids.
Computer Graphics Forum 31 (2012), 1948–1958. 16[Hoo98] H
OOVER
W.: Isomorphism linking smooth particles and embed-ded atoms.
Physica A: Statistical Mechanics and its Applications 260 , 3(1998), 244–254. 29[HS13] H
ORVATH
C. J., S
OLENTHALER
B.: Mass preserving multi-scale SPH. Pixar Technical Memo 13-04, Pixar Animation Studios,2013. 30[HWZ ∗
14] H E X., W
ANG
H., Z
HANG
F., W
ANG
H., W
ANG
G., Z
HOU
K.: Robust Simulation of Sparsely Sampled Thin Features in SPH-BasedFree Surface Flows.
ACM Transactions on Graphics 34 , 1 (2014), 7:1–7:9. 25, 37[IAAT12] I
HMSEN
M., A
KINCI
N., A
KINCI
G., T
ESCHNER
M.: Uni-fied spray, foam and air bubbles for particle-based fluids.
The VisualComputer 28 , 6-8 (2012), 669–677. 9[IABT11] I
HMSEN
M., A
KINCI
N., B
ECKER
M., T
ESCHNER
M.: AParallel SPH Implementation on Multi-Core CPUs.
Computer GraphicsForum 30 , 1 (Mar. 2011), 99–112. 10, 11, 37[IAGT10] I
HMSEN
M., A
KINCI
N., G
ISSLER
M., T
ESCHNER
M.:Boundary handling and adaptive time-stepping for PCISPH. In
VirtualReality Interactions and Physical Simulations (2010), pp. 79–88. 18[ICS ∗
14] I
HMSEN
M., C
ORNELIS
J., S
OLENTHALER
B., H
ORVATH
C.,T
ESCHNER
M.: Implicit incompressible SPH.
IEEE Transactions onVisualization and Computer Graphics 20 , 3 (2014), 426–435. 9, 12, 13,14, 16, 18, 37 c (cid:13) (cid:13)(cid:13)
IEEE Transactions onVisualization and Computer Graphics 20 , 3 (2014), 426–435. 9, 12, 13,14, 16, 18, 37 c (cid:13) (cid:13)(cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids [IOS ∗
14] I
HMSEN
M., O
RTHMANN
J., S
OLENTHALER
B., K
OLB
A.,T
ESCHNER
M.: SPH Fluids in Computer Graphics.
Eurographics (Stateof the Art Reports) (2014), 21–42. 18, 22, 27[JSD04] J
UBELGAS
M., S
PRINGEL
V., D
OLAG
K.: Thermal conductionin cosmological SPH simulations.
Monthly Notices of the Royal Astro-nomical Society 351 , 2 (2004), 423–435. 22[KAT ∗
19] K IM B., A
ZEVEDO
V., T
HUEREY
N., G
ROSS
M., S
OLEN - THALER
B.: Deep fluids: A generative network for parameterized fluidsimulations.
Computer Graphics Forum (2019). 36[KB17] K
OSCHIER
D., B
ENDER
J.: Density maps for improved sphboundary handling. In
ACM SIGGRAPH/Eurographics Symposium onComputer Animation (July 2017), pp. 1–10. 18, 19[KBT17] K
OSCHIER
D., B
ENDER
J., T
HUEREY
N.: Robust eXtendedFinite Elements for Complex Cutting of Deformables.
ACM Transac-tions on Graphics 36 , 4 (2017), 55:1–55:13. 31[KDBB17] K
OSCHIER
D., D
EUL
C., B
RAND
M., B
ENDER
J.: Anhp-adaptive discretization algorithm for signed distance field genera-tion.
IEEE Transactions on Visualization and Computer Graphics 23 ,10 (2017), 2208–2221. 18, 37[KKB18] K
UGELSTADT
T., K
OSCHIER
D., B
ENDER
J.: Fast corotatedFEM using operator splitting.
Computer Graphics Forum 37 , 8 (2018).31, 32[Kos19] K
OSCHIER
D.: CompactNSearch Library. https://github.com/InteractiveComputerGraphics/CompactNSearch , 2019. 10, 37[Lau11] L
AUTRUP
B.:
Physics of Continuous Matter . Taylor & Francis,2011. 22[LJS ∗
15] L
ADICKÝ
L., J
EONG
S., S
OLENTHALER
B., P
OLLEFEYS
M.,G
ROSS
M.: Data-driven fluid simulations using regression forests.
ACMTransactions on Graphics 34 , 6 (Oct. 2015), 199:1–199:9. 35, 36[LKR09] L AI M., K
REMPL
E., R
UBEN
D.:
Introduction to ContinuumMechanics . Butterworth-Heinemann, 2009. 7[LL10] L IU M., L IU G.: Smoothed Particle Hydrodynamics (SPH): anOverview and Recent Developments.
Archives of Computational Meth-ods in Engineering 17 , 1 (2010), 25–76. 3, 4[LLP11] L IU S., L IU Q., P
ENG
Q.: Realistic simulation of mixing fluids.
The Visual Computer 27 , 3 (2011), 241–248. 30[Łuk99] Ł
UKASZEWICZ
G.:
Micropolar Fluids . Modeling and Simula-tion in Science, Engineering and Technology. Birkhäuser Boston, 1999.27[MBCM16] M
ÜLLER
M., B
ENDER
J., C
HENTANEZ
N., M
ACKLIN
M.:A robust method to extract the rotational part of deformations. In
Pro-ceedings of ACM SIGGRAPH Conference on Motion in Games (2016),MIG ’16, ACM. 32[MFK ∗
15] M
AYRHOFER
A., F
ERRAND
M., K
ASSIOTIS
C., V
IOLEAU
D., M
OREL
F.-X.: Unified semi-analytical wall boundary conditionsin SPH: analytical extension to 3-D.
Numerical Algorithms 68 , 1 (Jan.2015), 15–34. 19[Mir96] M
IRTICH
B. V.:
Impulse-based dynamic simulation of rigidbody systems . PhD thesis, University of California at Berkeley, 1996.35[MM13] M
ACKLIN
M., M
ÜLLER
M.: Position Based Fluids.
ACMTransactions on Graphics 32 , 4 (2013), 1–5. 16, 27, 36, 37[MMCK14] M
ACKLIN
M., M
ÜLLER
M., C
HENTANEZ
N., K IM T.-Y.:Unified Particle Physics for Real-Time Applications.
ACM Transactionson Graphics 33 , 4 (2014), 1–12. 37[Mon92] M
ONAGHAN
J.: Smoothed Particle Hydrodynamics.
AnnualReview of Astronomy and Astrophysics 30 , 1 (1992), 543–574. 9, 22, 24,37[Mon05] M
ONAGHAN
J. J.: Smoothed Particle Hydrodynamics.
Reportson Progress in Physics 68 , 8 (2005), 1703–1759. 2, 22 [MSKG05] M
ÜLLER
M., S
OLENTHALER
B., K
EISER
R., G
ROSS
M.:Particle-based fluid-fluid interaction. In
ACM SIGGRAPH/EurographicsSymposium on Computer Animation (2005), p. 237. 29, 30[Nex17] N
EXT L IMIT : RealFlow Showreel 2017. , 2017. 1[PGBT17] P
EER
A., G
ISSLER
C., B
AND
S., T
ESCHNER
M.: An im-plicit sph formulation for incompressible linearly elastic solids.
Com-puter Graphics Forum (2017), n/a–n/a. 32, 33, 37[PICT15] P
EER
A., I
HMSEN
M., C
ORNELIS
J., T
ESCHNER
M.: AnImplicit Viscosity Formulation for SPH Fluids.
ACM Transactions onGraphics 34 , 4 (2015), 1–10. 22, 23, 25, 37[Pri12] P
RICE
D. J.: Smoothed particle hydrodynamics and magneto-hydrodynamics.
Journal of Computational Physics 231 , 3 (Feb. 2012),759–794. 2, 5, 6, 7, 22[PT16] P
EER
A., T
ESCHNER
M.: Prescribed velocity gradients forhighly viscous SPH fluids with vorticity diffusion.
IEEE Transactionson Visualization and Computer Graphics (2016), 1–9. 23, 25, 26, 37[RL96] R
ANDLES
P., L
IBERSKY
L.: Smoothed particle hydrodynam-ics: Some recent improvements and applications.
Computer Methods inApplied Mechanics and Engineering 139 , 1 (1996), 375 – 408. 5[RLY ∗
14] R EN B., L I C., Y AN X., L IN M. C., B
ONET
J., H U S.-M.:Multiple-Fluid SPH Simulation Using a Mixture Model.
ACM Transac-tions on Graphics 33 , 5 (2014), 1–11. 31[SB12] S
CHECHTER
H., B
RIDSON
R.: Ghost SPH for animating water.
ACM Transactions on Graphics 31 , 4 (2012), 61:1–61:8. 9, 22, 30, 37[SF18] S
CHENCK
C., F OX D.: Spnets: Differentiable fluid dynamics fordeep neural networks. In
CoRL (2018), vol. 87 of
Proceedings of Ma-chine Learning Research , PMLR, pp. 317–335. 36[SG11] S
OLENTHALER
B., G
ROSS
M.: Two-scale particle simulation.
TOG 30 , 4 (2011), 72:1–72:8. 30[Sif12] S
IFAKIS
E.:
SIGGRAPH 2012 Course Notes FEM Simulationof 3D Deformable Solids Part 1 . Tech. rep., University of Wisconsin-Madison, 2012. 31[SP08] S
OLENTHALER
B., P
AJAROLA
R.: Density Contrast SPH In-terfaces. In
ACM SIGGRAPH/Eurographics Symposium on ComputerAnimation (2008), pp. 211–218. 29, 30, 37[SP09] S
OLENTHALER
B., P
AJAROLA
R.: Predictive-corrective incom-pressible SPH.
ACM Transactions on Graphics 28 , 3 (2009), 40:1–40:6.9, 14, 18, 37[TDF ∗
15] T
AKAHASHI
T., D
OBASHI
Y., F
UJISHIRO
I., N
ISHITA
T.,L IN M.: Implicit Formulation for SPH-based Viscous Fluids.
ComputerGraphics Forum 34 , 2 (2015), 493–502. 22, 23, 24, 25, 26, 37[THM ∗
03] T
ESCHNER
M., H
EIDELBERGER
B., M
ÜLLER
M., P
OMER - ANTES
D., G
ROSS
M. H.: Optimized spatial hashing for collision de-tection of deformable objects. In
Vmv (2003), vol. 3, pp. 47–54. 10[TM05] T
ARTAKOVSKY
A., M
EAKIN
P.: Modeling of surface tensionand contact angles with smoothed particle hydrodynamics.
Physical Re-view E 72 , 2 (2005), 026301. 29, 30[TSSP16] T
OMPSON
J., S
CHLACHTER
K., S
PRECHMANN
P., P
ERLIN
K.: Accelerating Eulerian Fluid Simulation With Convolutional Net-works, jul 2016. 36[UHT18] U M K., H U X., T
HUEREY
N.: Liquid splash modeling withneural networks.
CGF 37 , 8 (2018), 171–182. 35[WBF ∗
96] W
ATKINS
S. J., B
HATTAL
A. S., F
RANCIS
N., T
URNER
J. A., W
HITWORTH
A. P.: A new prescription for viscosity in smoothedparticle hydrodynamics.
Astron. Astrophys. Suppl. Ser. 119 , 1 (1996),177–187. 22[WBT18] W
IEWEL
S., B
ECHER
M., T
HUEREY
N.: Latent-spacePhysics: Towards Learning the Temporal Evolution of Fluid Flow, feb2018. 36[WKB16] W
EILER
M., K
OSCHIER
D., B
ENDER
J.: Projective fluids. In
ACM Motion in Games (2016), pp. 1–6. 37 c (cid:13) (cid:13) . Koschier, J. Bender, B. Solenthaler & M. Teschner / SPH Techniques for the Physics Based Simulation of Fluids and Solids [WKBB18] W EILER
M., K
OSCHIER
D., B
RAND
M., B
ENDER
J.: Aphysically consistent implicit viscosity solver for sph fluids.
ComputerGraphics Forum 37 , 2 (2018). 22, 23, 24, 25, 26, 37[XFCT18] X IE Y., F
RANZ
E., C HU M., T
HUEREY
N.: tempoGAN:A temporally coherent, volumetric GAN for super-resolution fluid flow.
TOG 37 , 4 (2018), 95. 36[YCR ∗
15] Y
ANG
T., C
HANG
J., R EN B., L IN M. C., Z
HANG
J. J., H U S.-M.: Fast multiple-fluid simulation using helmholtz free energy.
ACMTransactions on Graphics 34 , 6 (Oct. 2015), 201:1–201:11. 31[YJL ∗
16] Y AN X., J
IANG
Y.-T., L I C.-F., M
ARTIN
R. R., H U S.-M.:Multiphase sph simulation for interactive fluids and solids.
ACM Trans-actions on Graphics 35 , 4 (July 2016), 79:1–79:11. 31 c (cid:13) (cid:13)(cid:13)