Turbulent Details Simulation for SPH Fluids via Vorticity Refinement
Sinuo Liu, Xiaokun Wang, Xiaojuan Ban, Yanrui Xu, Jing Zhou, Jiří Kosinka, Alexandru C.Telea
VVolume xx ( ), Number z, pp. 1–13
Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement
Sinuo Liu , Xiaokun Wang , Xiaojuan Ban , Yanrui Xu , Jing Zhou , Jiˇrí Kosinka , Alexandru C. Telea University of Science & Technology Beijing, Institute of Artificial Intelligence, China University of Science & Technology Beijing, School of Computer & Communication Engineering, China University of Groningen, Bernoulli Institute, Netherlands Utrecht University, Department of Information and Computing Sciences, NetherlandsWith vorticity refinementWithout vorticity refinement
Figure 1:
Our Vorticity Refinement (VR) solver applied to an DFSPH [BK15] simulation (1.18M particles). In this scene, a breaking damcollides with a board, creating turbulence. Zoom-ins compare the surface under DFSPH without (top) and with (bottom) our VR solver. Theresult shows that our method better captures turbulence details.
Abstract
A major issue in Smoothed Particle Hydrodynamics (SPH) approaches is the numerical dissipation during the projection pro-cess, especially under coarse discretizations. High-frequency details, such as turbulence and vortices, are smoothed out, leadingto unrealistic results. To address this issue, we introduce a Vorticity Refinement (VR) solver for SPH fluids with negligible com-putational overhead. In this method, the numerical dissipation of the vorticity field is recovered by the difference between thetheoretical and the actual vorticity, so as to enhance turbulence details. Instead of solving the Biot-Savart integrals, a streamfunction, which is easier and more efficient to solve, is used to relate the vorticity field to the velocity field. We obtain turbu-lence effects of different intensity levels by changing an adjustable parameter. Since the vorticity field is enhanced accordingto the curl field, our method can not only amplify existing vortices, but also capture additional turbulence. Our VR solver isstraightforward to implement and can be easily integrated into existing SPH methods.
CCS Concepts • Computing methodologies → Physical simulation; † Co-first author, contributed equally ‡ Corresponding author: [email protected]; [email protected] submitted to COMPUTER GRAPHICS
Forum (10/2020). a r X i v : . [ c s . G R ] S e p S. Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement
1. Introduction
Fluid simulation is a hot topic in computer graphics, with huge re-search and application demands. Within this context, the SmoothedParticle Hydrodynamics (SPH) method simulates fluids with largedeformations accurately and efficiently, showing abundant detailsand vivid motion. In the past decades, several solutions havebeen proposed to enforce incompressibility [SP09, ICS ∗
14, BK15].However, numerical dissipation problems still remain and causea significant loss of turbulence details [JST17, FGG ∗ ∗ ∗
19] looked into cor-recting the velocity field through the vorticity recovered from thekinematic viscosity by increasing the vorticity field proportionallyby the energy dissipated, which is based on the rotational kineticenergy. In contrast, we focus on getting the ideal vorticity field di-
Pressure projection ζ n = ∇ × 𝒗 n Vorticity
Velocity 𝒗 n 𝒗 𝑛𝑜𝑟𝑚𝑎𝑙 + 𝒗 𝑡𝑎𝑛𝑔𝑒𝑛𝑡𝑖𝑎𝑙 𝒗 𝑡𝑎𝑛𝑔𝑒𝑛𝑡𝑖𝑎𝑙 Δ ζ = ζ n+1 − ζ refine 𝒗 n+1 ζ = ∇ × 𝒗 𝑡𝑎𝑛𝑔𝑒𝑛𝑡𝑖𝑎𝑙 ∇ ⋅ 𝒗 𝑛𝑜𝑟𝑚𝑎𝑙 > 0 Advection Projection
Figure 2:
This diagram shows how the vorticity dissipates in theadvection-projection process. Linear velocity at time step n splitsinto normal (red) and tangential (green) components during theadvection procedure. Next, projection eliminates the normal com-ponent using the pressure force to keep the flow divergence-free. rectly from the curl of the Navier-Stokes equations, which is a morephysically reasonable model.Summarizing, the main contributions of this paper are: • A vivid turbulence-details generation method that recovers nu-merical dissipation through vorticity field correction; • A novel vorticity-based constraint and stream function solutionfor simulating turbulence; • An orthogonal solver for the SPH fluid framework with turbu-lence simulation that can be easily integrated into other particle-based methods and fluid solvers.The remainder of the paper is organized as follows. Section 2gives an overview of related work. Section 3 discusses the accuracyof numerical calculations in SPH. Our vorticity refinement schemeis presented in Section 4. Section 5 presents and discusses our ex-perimental results. Finally, Section 6 concludes the paper.
2. Related Work
Fluid simulation is a well researched topic in computer graphics.Early works on this topic include [Mon94, FM96, Sta99, MCG03].For recent overviews, we refer to Bridson’s book [Bri15] and thestate-of-the-art report of Koschier et al. [KBST19]. We further dis-cuss more specific work related to our context, namely SPH-basedfluid simulation (Sec. 2.1) and turbulence simulation (Sec. 2.2).
Monaghan simulated free surface flows with SPH [Mon94],which laid the foundation for fluid simulation. Later, Muller etal. [MCG03] proposed to simulate fluids using the ideal gas stateequation with surface tension and viscosity forces, but without fullincompressibility. An improved weakly-compressible SPH (WC-SPH) method was proposed by Becker and Teschner [BT07]. Theuse of the stiff equation of state (EOS) significantly increased re-alistic effects, but the efficiency of such methods is limited by thesize of the used time step. To further enforce incompressibility andimprove numerical accuracy, much effort has been invested intoimplicit pressure solvers. Previous approaches can be categorized submitted to COMPUTER GRAPHICS
Forum (10/2020). . Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement
Figure 3:
A breaking dam collides with a board (1.18M particles). (a) Only few vortex effects are formed using DFSPH. (b) The MP solverenhances the simulation result to some extent. (c) Our method greatly improves the turbulence details, and the surface details are richer thanwith the MP solver. The MP solver and our solver achieve different styles. as methods that project particle positions onto an incompressiblestate using iterative EOS solvers, and pressure projection meth-ods [IOS ∗ ∗
12] and Solenthaler and Pajarola [SP09] pro-posed predictive-corrective approaches that iteratively project par-ticle positions onto an incompressible state. This is also done inposition-based fluids (PBF) [MM13]. However, PBF avoids ac-cumulating pressure or pressure forces that eventually update thevelocity and the position. Ihmsen et al. [ICS ∗
14] proposed im-plicit incompressible SPH (IISPH) following the strategy of pres-sure projections. Separately, Bender and Koschier [BK15] pro-posed a method that enforces a low compression (0 .
01% ) and adivergence-free velocity constraint (DFSPH). Among all the vari-ants of the SPH method, the typical advection-projection modelsare PCISPH [SP09], IISPH [ICS ∗
14] and DFSPH [BK15]. In thispaper, we use the DFSPH approach as a baseline for comparisonsof computational efficiency and stability.It is well known that SPH approaches suffer from nu-merical dissipation problems, especially for coarse discretiza-tions [Mon94, dGWH ∗
15, BKKW18]. Ihmsen et al. [IOS ∗ Restoring high-frequency details has been an important challengein fluid simulation since its very beginning [KTJG08, JSMF ∗ ∗
15] successfully restoremost of the rotational motion using a hybrid method.Although the general simulation methods mentioned above canhandle numerical dissipation on a macroscopic level, both Eulerianand Lagrangian approaches face challenges when simulating high-frequency details such as turbulence. Therefore, methods specifi-cally designed for refining turbulent details have emerged. These submitted to COMPUTER GRAPHICS
Forum (10/2020).
S. Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement can be classified into three categories: up-res methods, vorticityconfinement methods, and Lagrangian vortex methods [BKKW18],as follows.
Up-res methods add high-frequency details over a coarse dis-cretization. Mercier et al. [MBT ∗
15] proposed a post-processingmethod to apply fine turbulence over particle-based fluid sur-faces. High-resolution surface points are seeded after curvatureevaluation, and the detailed surface waves are then evolved overcoarse particles. Edwards and Bridson [EB14] proposed an adap-tive volumetric-mesh method for grid-based fluids by using theadaptive discontinuous Galerkin method. Machine learning meth-ods such as Convolutional Neural Networks (CNNs) [CT17] havebeen applied in fluid simulation to synthesize high-resolution tur-bulence on rough simulation results based on a high-resolutionsource. However, training CNNs is time-consuming and often re-quires delicate hyperparameter tuning. Overall, up-res methods cantypically improve only surface effects.
Vorticity confinement methods aim to find existing vortices andrecover their dissipation. A new forcing term is added to increasethe velocity of target positions, and to enforce the rotation, of thevortex. Lentine et al. [LAF11] improved vorticity confinement tobe both energy conserving and momentum conserving. Jang etal. [JKB ∗
10] used multi-level vorticity confinement to acquire bet-ter results. Macklin and Muller [MM13] presented a simple methodto amplify the existing vorticity through accelerating particles us-ing SPH. Overall, vorticity confinement methods provide a simpleway for preserving vortices, but are in general unable to create ad-ditional turbulence details. Moreover, they are prone to adding ex-cessive energy to the system so that energy conservation is likely tobe violated, leading to unstable results.
Lagrangian vortex methods build on the vorticity representa-tion of the Navier-Stokes equations [PK05], which have less nu-merical dissipation and more divergence retention than vortic-ity confinement methods. These methods can be applied to par-ticles [WLB ∗ ∗
12] also treated boundaries of an Euleriangrid to solve this issue. A disadvantage of these methods is thatthe velocity field has to be recovered by solving the Biot-Savartintegrals or a vector-valued Poisson equation. Recently, Bender etal. [BKKW18] introduced the MicroPolar fluid solver (MP solver)for inviscid fluids in order to capture the micro-rotation of fluidparticles, achieving impressive visual turbulent features. Wang etal. [WLB ∗
20] proposed a turbulence refinement method based onthe Rankine vortex model for particle-based simulation. Zhang etal. [ZBG15] proposed an Integrated Vorticity of Convective Kine-matics (IVOCK) method to restore dissipated energy by measuringvorticity loss in advection. This method can cheaply capture muchof the lost details for smoke and fire, but does not work well forliquid simulations. In [ZBG15], only the vorticity dissipated dur-ing the advection step is considered. The refined linear velocity intheir paper is the velocity after the advection step. This velocityis then further affected by viscosity and the projection step. Vis-cosity may become another source of vorticity dissipation and the pressure force may introduce vorticity errors into the velocity fieldafter the projection step. Although it maintains an incompressibledensity field, it is not necessarily divergence-free.
Our method is inspired by the idea of stream functions [ZBG15],extended to Lagrangian fluid simulations. This allows us to effi-ciently derive velocity refinement from the vorticity field. Recov-ering turbulence from the curl form of the Navier-Stokes equationshas a long history. In 2005, Park and Kim [PK05] gave the gov-erning equations of the vortex method and introduced the conceptof the stream function. [ZBG15] and our work, among many othervortex methods, utilize this concept to reduce numerical dissipationduring simulation. In our method, we derive the dissipated vorticityduring the whole advection-projection step in the SPH approach.This can be easily done with little extra computation overhead.With respect to the concept of the stream function, we carry outthe Biot-Savart summation process within smoothing length, whichmakes it less accurate but more efficient than [ZBG15]; we showthis to be sufficient to maintain stability. This is because, theoret-ically, the refined velocity is the curl of the stream function, andany curl of a vector field is divergence-free. Moreover, we imple-mented our method using DFSPH (Divergence-free SPH), whichincludes an extra divergence-correction solver, thereby eliminatingpossible errors caused by the summation process. Moreover, we donot need to solve the Biot-Savart integrals or a vector-valued Pois-son equation. In contrast to the MP solver [BKKW18], in which themotion equation is obtained from the MicroPolar model and dis-cretized with SPH, we derive the vorticity equation from the curl ofthe Navier-Stokes equations, and recover velocity from the vortic-ity field using stream functions. Our results show that our methodcan not only enhance existing vortices, but also generate turbulenceat potential locations of new vortices.
3. SPH Discretization for Fluid Simulation
Traditional Lagrangian-based fluid simulations use the fluid gov-erning equations, the Navier-Stokes equations, to solve for the po-sition and velocity of each fluid particle. The acceleration of thefluid particles is obtained by the combination of pressure a pres , vis-cous force a vis , and gravity a g as DvvvDt = a pres + a vis + a g = − ρ ∇ p + ν v ∇ vvv + g , (1)where D denotes the material derivative, ρ is the density of the fluid, p represents pressure, vvv is velocity, ν v is the kinematic viscositycoefficient, a value that characterizes various fluid types (set to ν v = .
05 in our experiments), g is the gravitational acceleration, and ∇ denotes the Laplace operator.The SPH approach can be used to discretize the Navier-Stokesequations to numerically solve them. The continuous physical val-ues in space can be discretized using a smooth kernel W as in A ( xxx i ) = ∑ (cid:107) xxx i − xxx j (cid:107)≤ h m ( xxx j ) A ( xxx j ) ρ ( xxx j ) W (cid:0) xxx i − xxx j , h (cid:1) (2)with A ( xxx i ) being a certain quantity associated with particle i atlocation xxx i . This quantity can be interpolated from the values ofneighbour particles, indexed by j , within a support radius h . The submitted to COMPUTER GRAPHICS Forum (10/2020). . Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement DFSPH MP solver ( ν t = .
05) MP solver ( ν t =0.2) Our method ( α = . Figure 4:
Comparison of vorticity in a 2D scene extracted from the 3D scene in Fig. 3. This visualization corresponds to the second columnin Fig. 3. Color shows the vorticity magnitude of the particles, thereby allowing one to compare the vorticity of DFSPH, the MP solver, andour VR method. quantities m and ρ stand for mass and density, respectively. To sim-plify notation, we next use the shorthand A i to denote the quantity A evaluated at position xxx i .The density of a fluid can be derived by simply replacing A by ρ . In our work, we use the cubic spline kernel [Mon85]: W i j = π h − x + x ≤ x ≤ ( − x ) ≤ x ≤ x ≥ , where x = (cid:107) xxx i − xxx j (cid:107) / h and W i j is a short form of W (cid:0) xxx i − xxx j , h (cid:1) . Toobtain a better accuracy of the approximation of the divergence ofvelocity, the gradient and the curl of velocity, we apply the differ-ence form of the SPH discretization as: ∇ ⊗ A = ∑ j m j ρ j (cid:0) A j − A i (cid:1) ⊗ ∇ W i j , (3)which expresses the gradient ( ∇ A ), divergence ( ∇ · A ), and curl( ∇ × A , in which case the right hand side is negative) of A . Sincethe second derivative is often sensitive to particle disorder and signchanges inside the support radius h , we use artificial viscosity toapproximate the Laplacian as follows [KBST19]: ∇ A ( xxx i ) = ( d + ) ∑ j m j ρ j A i j · xxx i j xxx i j + . h ∇ W i j , (4)where d is the space dimension (in our case equal to 3), xxx i j = xxx i − xxx j , and A i j = A i − A j .Simulating incompressible fluids in DFSPH follows severalsteps, including advection and projection, and an extra diver-gence correction step which is applied to keep the velocity fielddivergence-free. The whole process is summarized in Algorithm 1,where ∆ t denotes the size of one time step, a adv = a vis + a g , and a pro j and a correct are the change rate of velocity derived formthe implicit pressure field to satisfy the incompressibility and divergence-free conditions accordingly. Further, ρ is the rest den-sity of the fluid, and ρ err , div err , n , and n (cid:48) are user-specified scalarvalues as thresholds. ALGORITHM 1:
Advection-projection with divergence correction
Advectiom process:compute a adv ˜ vvv : = vvv n + ∆ t a adv ˜ xxx : = xxx n + ∆ t ˜ vvv ˜ ρ : = positionBasedDensity ( ˜ xxx ) Projection process: while ( ˜ ρ − ρ ) > ρ err || numberO f Iterations < np : = positionBasedPressure ( ˜ xxx ) a proj : = pressureBasedForce ( p ) ˜ vvv : = ˜ vvv + ∆ t a proj ˜ ρ : = positionBasedDensity ( ˜ xxx + ∆ t ˜ vvv ) xxx n + = ˜ xxx Divergence correction process: while ( ∇ · ˜ vvv ) > div err || numberO f Iterations < n (cid:48) p : = velocityBasedPressure ( ˜ vvv ) a correct : = pressureBasedForce ( p ) ˜ vvv : = ˜ vvv + ∆ t a correct vvv n + : = ˜ vvv
4. Vorticity Refinement Model for Turbulence Simulation
Our method is closely related to Lagrangian vortex methods,namely it restores the velocity field through vorticity. In ourmethod, besides velocity vvv , each particle has a vector vorticity at-tribute ζζζ defined as ζζζ = ∇ × vvv . (5)In a particle system, vorticity is a quantity used to describe the ro-tation of a particle. For the vorticity at the position of particle i , thevalue can be derived using Eqn. 3 as: ζζζ i = ζζζ ( vvv i ) = ∇ × vvv i = ∑ j m j ρ j ( vvv i − vvv j ) × ∇ W i j . (6) submitted to COMPUTER GRAPHICS Forum (10/2020).
S. Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement (a) DFSPH(b) MP solver ( ν t = . α = . α = . Figure 5:
Comparison of DFSPH, the MP solver ( ν t = . ), andour method ( α = . and α = . ) for a simulation of a stick (rod)mixing water. Besides turbulence enhancement, our method standsout in keeping the flow trail visible and maintaining the stabilityof the surface. As visible, the MP solver and our solver can obtaindifferent enhancement results. Similarly to the divergence error issue [BK15], vorticity dissipa-tion can also hinder the performance of a simulation. Recent SPHapproaches [ICS ∗
14, BK15] for fluid animation can only correctnegative divergence of the velocity field. As a result, the kinetic en-ergy from the vorticity field is still allowed to be transformed intopositive divergence during simulation, causing the loss of surfacedetails and overall dynamic motions, effectively violating (the dis-crete version of) Eqn. 5. Given that the numerical dissipation of vorticity occurs betweentime steps, an ideal non-dissipative rate of change of vorticity isrequired to know the exact vorticity loss in each projection step. Weachieve this through the curl of the Navier-Stokes equation (Eqn.1)as: ∇ × (cid:18) DvvvDt (cid:19) = D ζζζ Dt = ζζζ · ∇ vvv + ν v ∇ ζζζ , (7)where ζζζ · ∇ vvv is the stretching term, which is vital for physicallymeaningful turbulence motion evolution. We use Eqn. 7 to obtainthe exact non-dissipative vorticity change of fluid particles betweentime steps, including boundary treatment [AIA ∗ ζζζ · ∇ vvv in Eqn. 7 is a vector, which we compute, percoordinate, using the difference form of the SPH approximation(Eqn. 3) via ∇ v { x , y , z } i = ∑ j m j ρ j (cid:16) v { x , y , z } j − v { x , y , z } i (cid:17) ∇ W i j , (8)where v xi is the x component of the velocity of particle with index i , and similarly for y and z . For the particle with index i , the vector ζζζ i · ∇ vvv i can be thus derived as ζζζ i · ∇ vvv i = ζζζ i · ∇ v xi ζζζ i · ∇ v yi ζζζ i · ∇ v zi . (9)The Laplacian of ζζζ in Eqn. 7 is derived using the artificial approx-imation analogous to Eqn. 4. Hence, for the particle with index i , ν v ∇ ζζζ i can be derived as ν v ∇ ζζζ i = ( d + ) ν v ∑ j m j ρ j ζζζ i j · xxx i j xxx i j + . h ∇ W i j . (10)According to Eqn. 7, the ideal change of the vorticity field withrespect to time, i.e., from time t n to t n + , is: ζζζ n + = ζζζ n + ∆ t D ζζζ n Dt , (11)and the dissipative vorticity update is given by: ∆ζζζ = ζζζ n + − ∇ × ˜ vvv , (12)where ˜ vvv is the (intermediate) velocity, as in the last line of Algo-rithm 1.We next explain how we apply the update of Eqn. 12. Assumethat we know the velocity and position of all fluid particles at time t n , and that the velocity at this time step is non-dissipative. We thenget the velocity and position at time t n + using the DFSPH ap-proach. Next, we compute the vorticity at the current time t n andthe next time t n + , denoted ζζζ n and ˜ ζζζ , respectively, from the ve-locity field using Eqn. 6. By our assumption, ζζζ n is ideal, but ˜ ζζζ is dissipative due to numerical integration. Thus the ideal vortic-ity value for a fluid particle at t n + , denoted ζζζ n + , is computedbased on ζζζ n and the vorticity equation (Eqn. 7). Hence, the dissipa-tive vorticity value for this particle in Eqn. 12 can be converted to ∆ζζζ = ζζζ n + − ˜ ζζζ . The dissipated vorticity is used to refine the velocityusing the stream function, as explained next. submitted to COMPUTER GRAPHICS Forum (10/2020). . Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement
Figure 6:
A breaking dam collides with a ship and 2 static pillars (1.70M particles). The key area is marked in white and zoomed in on. Ourmethod and the MP solver have more details than the DFSPH solution. In the second column, the MP solver becomes unstable and someparticles explode, while our method enhances the turbulence effect in a more stable and realistic way.
Inspired by [ZB14], we express the relationship between the veloc-ity vvv and the vorticity ζζζ using the stream function ψψψ as: vvv = ∇ × ψψψ , ∇ ψψψ = − ζζζ . (13)Greenâ ˘A ´Zs function provides a semi-analytical solution for thestream function. The derivation from the stream function to linearvelocity can be solved using Eqn. 3. Generalized by the Helmholtzdecomposition, the stream function is the vector potential ψψψ of thevelocity field vvv , which can be defined as ψψψ ( xxx ) = (cid:90) R ∇ × vvv ( yyy ) π (cid:107) xxx − yyy (cid:107) , (14)that is, the stream function ψψψ at position xxx is computed by integrat-ing the curl of velocity vvv at position yyy over the three-dimensionalspace R . Using Eqn. 5, we next discretize Eqn. 14 to get the streamfunction at the local position of particle with index i as: ψψψ i = π ∑ (cid:107) x i − x j (cid:107)≤ h ∆ζζζ j V j (cid:13)(cid:13) xxx i − xxx j (cid:13)(cid:13) , (15)where V j stands for the volume represented by the particle withindex j . Ideally, V j should be infinitely small and all distances (cid:107) xxx i − xxx j (cid:107) should be considered in the summation in Eqn. 15. However,to limit computational overhead and its adaptability to SPH, we only include neighbouring particles within a smoothing radius h inEqn. 15. This is justified by the fact that the influence of neighbourparticles shrinks with distance. Although the approximation couldpotentially induce instability and dissipation, our results show thatthis improves performance without sacrificing turbulent details, asalready observed e.g. in [MCG03].With the stream function obtained for each particle, the refinedvelocity for the particle with index i is derived as ∆ vvv i = ∑ j m j ρ j (cid:16) ψψψ i − ψψψ j (cid:17) × ∇ W i j . (16)To extend the flexibility of our method, we introduce an adjust-ment parameter α ∈ R , with the default value of 1 representingthe ideal vorticity refinement. It controls the amount of turbulenceadded to every simulation time step. Therefore the refined linearvelocity at t n + is expressed as vvv n + = ˜ vvv + α∆ vvv . (17)Since the divergence of the curl of any field is zero, the correctionof linear velocity due to vorticity does not cause any further diver-gence deviations. Hence, our method does not contradict any SPHprinciples, making it easier to implement into standard Lagrangianapproaches. Algorithm 2 summarizes our method, integrated withthe DFSPH technique for SPH simulation; see also Fig. 2. submitted to COMPUTER GRAPHICS Forum (10/2020).
S. Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement
ALGORITHM 2:
Our vorticity refinement (VR) solver
Compute current vorticity field: ζζζ n = ∇ × vvv n Advection-projection: vvv adv = advectPro ject ( vvv n ) Correct divergence field: ˜ vvv = correctDivergence ( vvv adv ) Vorticity through linear field: ˜ ζζζ = ∇ × ˜ vvv Compute vorticity equation: ζζζ n + = ζζζ n + ∆ t D ζζζ n Dt (Eqn. 11) Dissipation of vorticity: ∇ × vvv ( yyy ) = ζζζ n + − ˜ ζζζ (Eqn. 12) Compute stream function: ψψψ = (cid:82) R ∇ × vvv ( yyy ) π (cid:107) xxx − yyy (cid:107) (Eqn. 14) Refinement of linear velocity: ∆ vvv = ∇ × ψψψ (Eqn. 13) Refine linear velocity: vvv n + = ˜ vvv + α ∆ vvv (Eqn. 17)
5. Results and Discussion
We next test our novel Vortex Refinement (VR) method on sev-eral scenes, comparing it with the state-of-the-art micropolar (MP)model and classical SPH approaches.Both the VR and the MP method are integrated with DFSPH inthe following experiments to show the applicability of our method.We used the boundary handling method proposed by Akinci et al.[AIA ∗ α in our method, there isa scalar ν t in the MP method to control it. Based on the mechanismof the MP method [Eri66, BKKW18], ν t greater than ν v can poten-tially violate the second law of thermodynamics. In all experimentsbelow we set ν v = .
05. We therefore choose ν t = .
05 as a nat-ural refinement for the MP solver, which corresponds to α = α greater than 1 and ν t greater than 0 .
05; seeFigs. 5 and 9. As stated in [BKKW18], fluids are reasonably stablewhen ν t ≤ . To show the effectiveness of our approach numerically, we exe-cuted two breaking-dam experiments, and we executed two otherexperiments for parameter discussion and energy comparison withother methods, as follows.
Breaking dam with a board.
In Figs. 1 and 3, a board collides witha breaking dam which only allows fluid to go through the so-createdgap. Figure 3 shows the results with 1.18M particles. Only few vor-tex effects can be seen using DFSPH. Water flushes through the gapand dissipates quickly without clear turbulence effects. Comparedto DFSPH, our solver generates several realistic vortices aroundthe board and corners. The MP solver also improves the visual re-sult, but not as obviously as our method. Since our method refinesparticle velocity based on the vorticity field, vortices are naturallypreserved and turbulence is generated from the dissipated energyin a realistic way. In Fig. 4, the vorticity magnitude of all parti-cles is visualized. The comparison shows that both our method andthe MP method yield higher energy values than DFSPH. The MPsolver adds energy in a natural way, while our method recovers en-ergy from numerical dissipation more effectively and is thus ableto simulate more details.
Breaking dam with three obstacles.
As shown in Fig. 11, a break-ing dam scenario with static obstacles was tested using 457K fluidparticles. The fluid flows in from the left and hits the wall on theright. Several waves are generated in the process, which then comeback and interact with three rigid bodies. Desirable turbulence canbe observed over the surface. We compared our method with theDFSPH and MP solvers. In DFSPH, the fluid seems to go aroundthe pillars and forms splashes, but scarcely any complex turbulenceeffects. In contrast, our solver creates small-scale vortices insteadof just the fluid smoothly flowing around the pillars. Since thesesmall vortices cannot sustain a self-spinning state, they quicklybreak down into turbulence. Compared to the MP solver in thisscene, our method seems to generate more turbulent details butsmaller vortices. The MP solver and our method can achieve dif-ferent visual effects.
Energy Comparison.
An energy comparison of a breaking damexperiment (see Fig. 8) is shown in Fig. 7. The left plot shows theenergy comparison, while the right plot shows the energy increaseratio relative to DFSPH. When t ∈ [ , ] , the fluid keeps flowingand forming turbulence. If the energy is larger than that of DFSPH,then energy is recovered (or added) successfully. After the watersurface calms down (after about 10s), the scene should contain onlypotential energy (no kinetic energy). The energy of the traditionalDFPSH method can be used as a benchmark: If a method gener-ates, at this time point, more energy than DFSPH, then this methodis considered to create additional energy. In this comparison ex-periment, the energy values after 10s for both the MP solver with ν t = . , . α = . , . , . α = . ν t = . , . , . α = . Breaking dam with a hemisphere: parameter influence.
In thisexperiment (see Fig. 9) we flush a hemisphere obstacle with a fixedvolume of fluid. This means only limited kinetic energy is involvedin this scenario (from gravitational potential energy). We simulatedthe flow using DFSPH, our method with α = α = .
2, andthe MP solver with ν t = . ν t = . ν t = .
4. When com-paring the DFSPH approach with our method with α = ν t = .
05, both methods are able to increase the turbu-lence performance, but our result is more pronounced than the MPone. To obtain more obvious turbulence effects, we increase theturbulence control parameters in the two methods, which meansthat more energy is added to the simulation. The renderings showthat our method with α = . v t = .
2. To keep ourmethod in line with the underlying physics, as explained for the ear-lier example, we do not use higher parameter values. The MP solveradds more turbulence in this scene. The obtained results are visu-ally more salient for large parameter values, e.g. ν t = .
2. However, ν t cannot be increased indefinitely. For example, if we set ν t = . submitted to COMPUTER GRAPHICS Forum (10/2020). . Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement Energy
T i m e ( s )
D F S P H V R 1 . 0 0 V R 1 . 1 0 V R 1 . 2 0 V R 1 . 3 0 M P 0 . 0 5 M P 0 . 1 0 M P 0 . 2 0 M P 0 . 3 0 M P 0 . 4 0
Percentage
V R 1 . 0 0 V R 1 . 1 0 V R 1 . 2 0 V R 1 . 3 0 M P 0 . 0 5 M P 0 . 1 0 M P 0 . 2 0 M P 0 . 3 0 M P 0 . 4 0
Figure 7:
Comparison of energy changes for different methods using different parameter values in the breaking dam scene shown in Fig. 8.Left: direct energy comparison. Right: energy increase ratio relative to DFSPH. When the fluid is flowing (t ∈ [ , ] ), our method and the MPsolver are able to add energy to the scene and enhance the visual effect. However, the MP solver with ν t = . , . , . and our method with α = . do not converge after t = due to the excessive energy added. Figure 8:
A breaking dam scenario for the evaluation of energychanges (461K fluid particles); see Fig. 7. can be applied to scenes that are more sensitive to physics laws,such as adding more details to a relatively stably-flowing scene. Incontrast, the MP method can be used in scenes where one wants tocreate a stronger visual impact, such as collapses or violent shocks.Overall, this experiment shows that the MP solver and our solvercan achieve different turbulence effects. Our method achieves bet-ter turbulence results without adding energy sources. In contrast,the MP solver can add small vortices, but when increasing its pa-rameter values, energy sources will pop up and prevent the fluidfrom calming down.
To further demonstrate the turbulence quality of our method, wesimulated several complex scenarios with dynamic boundary con-ditions and compared them with the MP solver.
Spinning Propeller.
A propeller is slowly submerged into water,after which it starts spinning at 3 radians per second. Fig. 12 showsthe results of this simulation using 1.29M fluid particles for DF-SPH, MP, and VR (our method). Observe that neither the complexflow nor strong turbulence effects are produced and preserved using DFSPH. Both our method and the MP method enhance the visualeffect. In contrast to the MP method, our method adds energy in aphysically reasonable way (no turbulence in front of the propeller)and creates vivid turbulent details over the free fluid surface. Thekey areas are zoomed in on. Also, a vortex is observed with ourmethod after the propeller has stopped spinning (see also the sup-plementary video).
Boat-sinking.
In this scenario, a boat and two columns interactwith a breaking dam. Figure 6 shows the results using 1.7M fluidparticles. The potential energy of the fluid transforms into the ki-netic energy of the fluid particles and the boat. The water is firstviolently displaced when it hits the column and the boat, and nextgradually calms down as time goes by, finally reaching a stablestate. We see that the DFSPH method produces relatively weaklyturbulent details, which get lost quickly due to numerical dissipa-tion. In contrast, our method and MP server shows more naturaldynamics with realistic turbulent effects on the fluid surface. Thefluid gradually calms down as time goes on. Our method and theMP method achieve different styles.
Stirring water.
In Fig. 5, a cylindrical stick was inserted into atank of water, and stirred at a uniform speed for several seconds.The water splashed around due to the quick movement of the stick.Observe that the trace left on the surface lasts longer in our methodthan with the MP method, which is a critical point for boat-sailinganimation scenarios. After the stirring process, the stick is pulledout of the fluid, and the water starts to calm down. The DFSPH ap-proach calms the fluid down quickly due to numerical dissipation.The surface details are clearer and sharper in our method. Also, wenotice a disturbance wave in the MP method, caused by the fact that ν t exceeds the kinematic viscosity.The above three scenarios show that our method can keep sta-bility when dealing with extreme conditions like strong collisions,while physically preserving energy. Moreover, in the accompany- submitted to COMPUTER GRAPHICS Forum (10/2020). S. Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement S. Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement (a) DFSPH(b) MP solver ( n t = . n t = . n t = . a = . a = . Figure 9:
Comparison of DFSPH, the MP solver with n t = . , . , . , . , and our method with a = . , . in the breaking dam scenariowith a static spherical obstacle placed to the right. The water splashed around due to the quick movement of the stick.Observe that the trace left on the surface lasts longer in our methodthan with the MP method, which is a critical point for boat-sailinganimation scenarios. After the stirring process, the stick is pulled out of the fluid, and the water starts to calm down. The DFSPH ap-proach calms the fluid down quickly due to numerical dissipation.The surface details are clearer and sharper in our method. Also, we submitted to COMPUTER GRAPHICS
Forum (9/2020).
Figure 9:
Comparison of DFSPH, the MP solver with ν t = . , . , . , . , and our method with α = . , . in the breaking dam scenariowith a static spherical obstacle placed to the right. ing video it can be seen that our method not only amplifies existingvortices but also generates new ones. Computational overhead.
The computational overhead of ourmethod is negligible compared to the whole SPH simulation pro-cedure. Table 1 shows the computing times for DFSPH, the MPsolver, and our method for different simulation scenes. The dif-ferent computation times are explained as follows. Compared toDFSPH, both turbulence methods (MP and ours) need to compute the vorticity field, i.e., solve for the Laplacian ∇ ζζζ . Further, ∇ × ζζζ (in the MP solver) and ∇ vvv (in our method) also need to be solvedfor. The difference is that our method needs to compute ψψψ ( ζζζ ) and ∇ × ψψψ to get the refined velocity, but as Table 1 shows, the extracomputational effort is negligible. submitted to COMPUTER GRAPHICS Forum (10/2020). . Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement Percentage
T i m e ( s )
V R 1 . 0 0 V R 1 . 2 0 M P 0 . 0 5 M P 0 . 2 0 M P 0 . 4 0
Figure 10:
Similar to Fig. 8, we compared the energy increase ratioof different methods to the DFSPH method using the scene shown inFig. 9. We see that all methods produce more energy than DFSPHwhen the fluid is flowing (t ∈ [ , ] ). Note that the MP solver with ν t = . , . does not converge after t = . Experiment Fig. Particles ∆ t (ms) Steps DFSPH (m) MP (m) VR (m)Board 3 1.18M 2.4 9542 2401.2 2565.9 2565.1Stirring 5 1.39M 2.4 9542 2399.9 2864.4 2693.4Sphere 9 899.8K 2.4 8375 1657.4 1715.2 1827.2Pillars 11 457K 3 6667 156.4 188.3 218.5Propeller 12 1.29M 3 7334 1782.1 2309.1 2338.9 Table 1:
Total time comparisons of three methods: DFSPH, MP,and our method (VR) over five simulations. ∆ t, in milliseconds, isthe time step used in the experiments, and the total computationtimes, in minutes, include the costs of the density solver and thedivergence-free solver in DFSPH.
6. Conclusion and Discussion
We have presented a particle-based turbulence refinement methodthat recovers lost velocity from the difference between the theoreti-cal and the actual vorticity value. Our method can not only increaseexisting vortices significantly by recovering numerical dissipation,but also generates new turbulence at potentially different locations.The turbulence-enhancement parameter of our method has a theo-retically optimal value α = (a) DFSPH(b) MP solver(c) Our method Figure 11:
A comparison of DFSPH, the MP solver and ourmethod in a breaking dam scene with 3 static pillars as obstacles.The water hits the cylinders and the right wall and bounces back,forming turbulence in the process. Compared to DFSPH, our solvergives rise to tiny vortices instead of the water simply going aroundthe pillars. Compared to the MP solver, our method seems to pro-vide more turbulent details but smaller vortices in this scene. ods. Our method can simulate typical turbulent scenes efficientlyand is relatively stable even for scenarios with highly turbulentflow. At the same time, we should note that some vorticity is lostin such cases. While this small amount of loss does not affect thegeneral visual quality, decreasing it is an open topic for future re-search, which can be expected to lead to even more realistic fluidsimulations.In the future, we aim to investigate merging our method with mi-crostructural models, since these models show great potential forrough simulation conditions and also have a close relationship withviscosity. Improving computation accuracy is another potential fu-ture research direction. Finally, increasing the computational scal-ability of our method by e.g. efficient and effective parallelizationis attractive for making our method directly applicable to complexreal-world and/or interactive simulations.
References [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 (TOG) 31 , 4 (2012), 62.[AN05] A
NGELIDIS
A., N
EYRET
F.: Simulation of smoke based submitted to COMPUTER GRAPHICS
Forum (10/2020). S. Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement (a) DFSPH(b) MP solver(c) Our method
Figure 12:
A propeller interacts with 1.29M fluid particles using DFSPH, the MP solver, and our method. DFSPH (top row) is not able toproduce a complex flow. Under stable conditions, the MP solver (middle row) can only add limited turbulence effects to DFSPH, while ourVR method is able to generate realistic vortices (bottom row). The improved turbulence performance can be seen clearly in the insets. on vortex filament primitives. In
Proceedings of the 2005 ACMSIGGRAPH/Eurographics symposium on Computer animation (2005),ACM, pp. 87–96.[BK15] B
ENDER
J., K
OSCHIER
D.: Divergence-free smoothed par-ticle hydrodynamics. In
Proceedings of the 14th ACM SIG-GRAPH/Eurographics Symposium on Computer Animation (2015),ACM, pp. 147–155.[BKKW18] B
ENDER
J., K
OSCHIER
D., K
UGELSTADT
T., W
EILER
M.:Turbulent micropolar sph fluids with foam.
IEEE Transact. on Visual-ization and Computer Graphics (2018).[Bri15] B
RIDSON
R.:
Fluid simulation for computer graphics . AK Pe-ters/CRC Press, 2015.[BT07] B
ECKER
M., T
ESCHNER
M.: Weakly compressible sphfor free surface flows. In
Proceedings of the 2007 ACM SIG-GRAPH/Eurographics symposium on Computer animation (2007), Eu-rographics Association, pp. 209–217.[CT17] C HU M., T
HUEREY
N.: Data-driven synthesis of smoke flowswith cnn-based feature descriptors.
ACM Transactions on Graphics(TOG) 36 , 4 (2017), 69.[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 Trans. Graph. 34 , 4 (2015), 50–1. [EB14] E
DWARDS
E., B
RIDSON
R.: Detailed water with coarse grids:combining surface meshes and adaptive discontinuous galerkin.
ACMTransact. on Graphics (TOG) 33 , 4 (2014), 136.[Eri66] E
RINGEN
A. C.: Theory of micropolar fluids.
Journal of Mathe-matics and Mechanics (1966), 1–18.[EWPT17] E
BERHARDT
S., W
EISSMANN
S., P
INKALL
U., T
HUEREY
N.: Hierarchical vorticity skeletons. In
Proceedings of the ACMSIGGRAPH/Eurographics Symposium on Computer Animation (2017),ACM, p. 6.[FGG ∗
17] F U C., G UO Q., G
AST
T., J
IANG
C., T
ERAN
J.: A poly-nomial particle-in-cell method.
ACM Trans. Graph. 36 , 6 (Nov. 2017),222:1–222:12.[FM96] F
OSTER
N., M
ETAXAS
D.: Realistic animation of liquids.
Graphical models and image processing 58 , 5 (1996), 471–483.[FSJ01] F
EDKIW
R., S
TAM
J., J
ENSEN
H. W.: Visual simulation ofsmoke. In
Proceedings of the 28th conference on Computer graphicsand interactive techniques (2001), ACM, pp. 15–22.[GNS ∗
12] G
OLAS
A., N
ARAIN
R., S
EWALL
J., K
RAJCEVSKI
P.,D
UBEY
P., L IN M.: Large-scale fluid simulation using velocity-vorticitydomain decomposition.
ACM Transactions on Graphics (TOG) 31 , 6(2012), 148. submitted to COMPUTER GRAPHICS
Forum (10/2020). . Liu & X. Wang & X. Ban et al. / Turbulent Details Simulationfor SPH Fluids via Vorticity Refinement ∗
12] H E X., L IU N., L I S., W
ANG
H., W
ANG
G.: Local poissonsph for viscous incompressible fluids. In
Computer Graphics Forum (2012), vol. 31, Wiley Online Library, pp. 1948–1958.[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.[IOS ∗
14] I
HMSEN
M., O
RTHMANN
J., S
OLENTHALER
B., K
OLB
A.,T
ESCHNER
M.: Sph fluids in computer graphics.[JKB ∗
10] J
ANG
T., K IM H., B AE J., S EO J., N OH J.: Multilevel vorticityconfinement for water turbulence simulation.
The Visual Computer 26 ,6-8 (2010), 873–881.[JSMF ∗
18] J
ESCHKE
S., S
K ˇRIVAN
T., M
ÜLLER -F ISCHER
M., C
HEN - TANEZ
N., M
ACKLIN
M., W
OJTAN
C.: Water surface wavelets.
ACMTrans. Graph. 37 , 4 (July 2018), 94:1–94:13.[JSS ∗
15] J
IANG
C., S
CHROEDER
C., S
ELLE
A., T
ERAN
J., S
TOM - AKHIN
A.: The affine particle-in-cell method.
ACM Transactions onGraphics (TOG) 34 , 4 (2015), 51.[JST17] J
IANG
C., S
CHROEDER
C., T
ERAN
J.: An angular momen-tum conserving affine-particle-in-cell method.
Journal of ComputationalPhysics 338 (2017), 137–164.[KBST19] K
OSCHIER
D., B
ENDER
J., S
OLENTHALER
B., T
ESCHNER
M.: Smoothed particle hydrodynamics techniques for the physics basedsimulation of fluids and solids.[KLLR05] K IM B., L IU Y., L
LAMAS
I., R
OSSIGNAC
J. R.:
Flowfixer:Using bfecc for fluid simulation . Tech. rep., Georgia Institute of Tech-nology, 2005.[KTJG08] K IM T., T
HÜREY
N., J
AMES
D., G
ROSS
M.: Wavelet tur-bulence for fluid simulation. In
ACM Transactions on Graphics (TOG) (2008), vol. 27, ACM, p. 50.[LAF11] L
ENTINE
M., A
ANJANEYA
M., F
EDKIW
R.: Mass and mo-mentum conservation for fluid simulation. In
Proceedings of the 2011ACM SIGGRAPH/Eurographics Symposium on Computer Animation (2011), ACM, pp. 91–100.[MBT ∗
15] M
ERCIER
O., B
EAUCHEMIN
C., T
HUEREY
N., K IM T.,N
OWROUZEZAHRAI
D.: Surface turbulence for particle-based liquidsimulations.
ACM Transactions on Graphics (TOG) 34 , 6 (2015), 202.[MCG03] M
ÜLLER
M., C
HARYPAR
D., G
ROSS
M.: Particle-based fluidsimulation for interactive applications. In
Proceedings of the 2003 ACMSIGGRAPH/Eurographics symposium on Computer animation (2003),Eurographics Association, pp. 154–159.[MM13] M
ACKLIN
M., M
ÜLLER
M.: Position based fluids.
ACM Trans-actions on Graphics (TOG) 32 , 4 (2013), 104.[Mon85] M
ONAGHAN
J.: Particle methods for hydrodynamics.
Com-puter Physics Reports 3 , 2 (1985), 71–124.[Mon94] M
ONAGHAN
J. J.: Simulating free surface flows with sph.
Jour-nal of computational physics 110 , 2 (1994), 399–406.[NZT19] N
ARAIN
R., Z
EHNDER
J., T
HOMASZEWSKI
B.: A second-order advection-reflection solver.
Proceedings of the ACM on ComputerGraphics and Interactive Techniques 2 , 2 (2019), 1–14.[PK05] P
ARK
S. I., K IM M. J.: Vortex fluid for gaseous phenomena. In
Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium onComputer animation (2005), ACM, pp. 261–270.[SP09] S
OLENTHALER
B., P
AJAROLA
R.: Predictive-corrective incom-pressible sph. In
ACM transactions on graphics (TOG) (2009), vol. 28,ACM, p. 40.[Sta99] S
TAM
J.: Stable fluids. In
Proceedings of the 26th annual con-ference on Computer graphics and interactive techniques (1999), ACMPress/Addison-Wesley Publishing Co., pp. 121–128.[WLB ∗
20] W
ANG
X., L IU S., B AN X., X U Y., Z
HOU
J., K
OSINKA
J.:Robust turbulence simulation for particle-based fluids using the rankinevortex model. In (2020), IEEE, pp. 657–658. [WP10] W
EISSMANN
S., P
INKALL
U.: Filament-based smoke with vor-tex shedding and variational reconnection. In
ACM Transactions onGraphics (TOG) (2010), vol. 29, ACM, p. 115.[XBP ∗
19] X U Y., B AN X., P
ENG
Y., W
ANG
X., L IU S., Z
HOU
J.: Tur-bulence enhancement for sph fluids visualization. In
International Con-ference on Cooperative Design, Visualization and Engineering (2019),Springer, pp. 254–260.[ZB05] Z HU Y., B
RIDSON
R.: Animating sand as a fluid.
ACM Trans-actions on Graphics (TOG) 24 , 3 (2005), 965–972.[ZB14] Z
HANG
X., B
RIDSON
R.: A pppm fast summation method forfluids and beyond.
ACM Transactions on Graphics (TOG) 33 , 6 (2014),206.[ZBG15] Z
HANG
X., B
RIDSON
R., G
REIF
C.: Restoring the missingvorticity in advection-projection fluid solvers.
ACM Transactions onGraphics (TOG) 34 , 4 (2015), 52.[ZNT18] Z
EHNDER
J., N
ARAIN
R., T
HOMASZEWSKI
B.: An advection-reflection solver for detail-preserving fluid simulation.
ACM Transac-tions on Graphics (TOG) 37 , 4 (2018), 1–8.[ZYF10] Z HU B., Y
ANG
X., F AN Y.: Creating and preserving vorticaldetails in sph fluid. In
Computer Graphics Forum (2010), vol. 29, WileyOnline Library, pp. 2207–2214. submitted to COMPUTER GRAPHICS