Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges

Physics-based fluid simulation has played an increasingly important role in the computer graphics community. Recent methods in this area have greatly improved in generating complex visual effects and also in computational efficiency. Novel techniques have emerged to deal with complex boundaries, multiphase fluids, gas-liquid interfaces, and fine details. In parallel, the combined use of machine learning, image processing, and fluid control technologies has brought many interesting and novel research perspectives. In this survey, we provide an introduction to theoretical concepts underpinning physics-based fluid simulation and their practical implementation with the aim to serve as a guide for both newcomer and seasoned researchers for exploring the physics-based fluid simulation field, with a focus on recent developments in the last decade. Driven by the distribution of recent publications in the field, we structure our survey to cover physical background, discretization approaches, computational methods that address scalability, fluid interactions with other materials and interfaces, and methods for expressive aspects of surface detail and control. From a practical perspective, we overview existing implementations available for the above methods.


Introduction
Fluids are a crucial element in visual simulations given their ubiquitous existence in natural environments.Their versatile motion and complex behavior make fluids an attractive, but also tricky to describe and compute, target for graphics simulations.As such, the simulation of fluids has long been one of the most important subjects in computer graphics.The development of computer technology has made it possible to simulate complex fluid phenomena directly using the governing equations of fluid dynamics.Scores of physicsbased methods and techniques have been proposed to this end, ranging from simple but inaccurate models to progressively refined, complex, techniques that capture increasingly challenging dynamics of the interacting media.The diversity of such methods, with widely varying assumptions, modeling power, and underlying implementation techniques, has made understanding the current state of the art increasingly challenging, especially for practitioners and newcomers.
This survey aims to address the above understanding challenge by covering the major research topics of physics-based fluid simulation in computer graphics and new trends in those topics over the last decade.It discusses the different goals of fluid simulation, techniques proposed to address such goals, challenges of such techniques, and key findings in the field.We structure our survey in a top-down manner as follows.Section 2 presents our methodology used in collecting relevant work from the literature (with a focus on the last decade) and proposes a classification of fluid simulation into seven main topics.Section 3 introduces the physical background and discretization approaches used by fluid simulation required to understand the remainder of the survey.Sections 4-10 discuss the papers found in each of the above-mentioned seven topics, as follows.Section 4 presents the various classes of advanced adaptive, parallel, and data-driven computational approaches used to accelerate the implementation of the simulation models outlined in Sec. 3. The next three sections detail more specific and challenging simulation contextsfluid interaction with different materials (Sec.5), multiphase simulations (Sec.6), and gas-liquid interfaces (Sec.7).Sections 8 and 9 discuss artistic measures for improving the quality, respectively controlling the appearance and motion, of the simulated fluids.Section 10 looks into special fluids.Finally, Section 11 concludes our survey by identifying key directions for future research.Appendix A presents and discusses software implementations of fluid simulation covering numerical simulation, modeling, and rendering aspects.

Methodology for constructing the survey
Physics-based fluid simulation is a research area that has been active for many decades with input from fields as diverse as engineering, physics, mathematics, and computing science.Providing a survey of fluid simulation covering all these fields is too large a challenge for a single paper.Furthermore, we believe that the interests of typical researchers and practitioners in computer graphics focus on a subset of the above aspects, and structure our survey accordingly, as follows.
We selected as main information sources articles published in ACM Transactions on Graphics (TOG), IEEE Transactions on Visualization and Computer Graphics (TVCG), and Computer Graphics Forum (CGF), which are arguably the three most influential and representative computer graphics journals.As our survey aims to cover recent tendencies, we included all relevant papers from these journals published in the last decade (2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021)(2022).We further included some papers presented at key graphics conferences like ACM SIGGRAPH / Eurographics Symposium on Computer Animation (SCA), ACM SIGGRAPH (SIGGRAPH), and ACM SIGGRAPH Asia Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 3 (SIGGRAPH Asia).Finally, we included earlier papers that still significantly impact recent research.In total, we collected and further analyzed 327 such papers.The origins of physics-based fluid simulation in computer graphics can be traced back to the 1970s with the development of 'particle systems' [1].However, a fully-developed, reasonably stable physics-based system for fluid animation was not achieved until the end of the last century [2].In the first decade afterwards, various simulation strategies continued to evolve, focusing on improving stability, accuracy, and efficiency in fluid simulations.We discuss early development progress in this respect in Sec.3.3.Concurrently, research on specific fluid phenomena and the tailoring of simulated effects began to emerge and gain momentum.
As we entered the second decade of the 21st century, which is the primary focus of this survey paper, main research interests in fluid simulation shifted towards addressing specific effects that are challenging to achieve using conventional fluid simulation methods.Alongside this shift, advances in machine learning technologies have opened up new ways for integrating neural networks with simulation algorithms, pushing the boundaries of what can be accomplished in fluid simulations.As this survey aims to discuss the current advancements in this area comprehensively, we classify the collected papers into seven relevant topics that span the past decade based on our detailed analysis of these papers.Next, we select a subset of representative papers within each topic and discuss these in greater detail.
Figure 1 shows the seven identified topics at the first level of a hierarchy depicting our survey's structure.Further levels refine these into sub-topics.The seven main topics are as follows: • Advanced computational approaches: (Sec.4) Methods that aim to make full use of powerful parallel computing resources for fluid simulation; • Fluid coupling with multi-materials: (Sec.5) Methods that model the interaction between fluid and solid objects of various sorts of shapes and textures; • Multiphase liquids: (Sec.6) Methods for the simulation of liquid-liquid interaction effects of various phases; • Gas-liquid interfaces: (Sec.7) Methods dealing with scenarios where forces on the gas-liquid interfaces dominate the fluid motion; • Fine details enhancement: (Sec.8) Methods that concentrate on preserving/enhancing fluid motion on a detail level; • Fluid control: (Sec.9) Methods that allow visual designers to control the appearance and style of fluid simulations; • Special fluids: (Sec.10) Methods that simulate nonconventional fluids, e.g., highly viscous/thin, sensitive to magnetic fields, or targeting materials that are strictly speaking not fluids but behave like fluids.Our seven-topic classification aims first and foremost to identify salient trends in the past decade.As such, the emergence of these topics is based on a significant portion of the found papers that can be grouped within each topic -in other words, the topics reflect a 'data-driven organization' of how research in fluids in the past decade proceeded.This is in contrast with other surveys which group works based on a predefined, model-driven taxonomy proposed by the respective survey authors.
It is insightful to analyze how the found papers distribute over topics and how these topics evolve over time.Figure 2 shows the number of papers found for each of the seven identified topics.We see that, while variations exist (some topics being more popular than others), each topic has a significant number of papers, with 4-8 top-tier papers per year on average, thereby supporting the claim that our identified topics are a good way to organize research.Figure 3 refines the above insight, showing the number of articles per topic and per year.We see that, while some topics show an increase in publications (e.g.advanced computational approaches, the red curve in Fig. 3), all identified topics have been 'alive' over the past decade -another indication that they are suitable for organizing our survey.F l u i d c o u p l i n g w i t h m u l t i -m a t e r i a l s M u l t i p h a s e l i q u i d s G a s -l i q u i d i n t e r f a c e s

Fluid simulation overview
The development of fluid simulation in computer graphics is deeply rooted in the history of physics.In the 19th century, scientists such as Sir Isaac Newton and Claude-Louis Navier c o u n t y e a r A d v a n c e d c o m p u t a t i o n a l a p p r o a c h e s F l u i d c o u p l i n g w i t h m u l t i -m a t e r i a l s M u l t i p h a s e l i q u i d s G a s -l i q u i d i n t e r f a c e s F i n e d e t a i l s e n h a n c e m e n t F l u i d c o n t r o l S p e c i a l f l u i d s Fig. 3 Trends in number of papers published per topic over the past decade.contributed significantly to the understanding of fluid mechanics, paving the way for the Navier-Stokes equations.These equations, which govern fluid motion, form the foundation of modern fluid simulation algorithms.
This section offers a basic introduction to fluid simulation and provides background knowledge for the rest of the survey.For a more comprehensive understanding of fluid simulation, we refer to Bridson's book [3].For more specific knowledge about Lagrangian-based smoothed particle hydrodynamics and material point methods, we refer to the surveys of Koschier et al. [4] and Jiang et al. [5] respectively.Readers less familiar with this area are highly encouraged to read this section before diving into the seven topics described next.We first introduce relevant physical principles behind fluid simulation such as the continuum hypothesis (Sec.3.1) and Navier-Stokes equations (Sec.3.2).We next present the early development of this area (Sec.3.3), including a brief overview of the ideas behind different discretization strategies.
Table 1 summarizes the main notations used in our survey.We use these notations accompanied by subscripts, superscripts, and brackets with parameters to denote these quantities under various conditions, as explained further in context.

Fluid mechanics
Matter in nature is built up of atoms and molecules that are discrete and separated by space.Simulating fluid at the microscopic level to describe macroscopic phenomena is only doable on supercomputers with weeks, if not months, of computing time.Computer graphics pursues a balance between efficiency and fidelity.For this, fluid mechanics based on a continuum hypothesis is the level at which physical properties are modeled.
Fluid mechanics models an object with matter continuously distributed over its body, an approximation called the Stress tensor P a Table 1 Common notations (and their descriptions) used in this paper.For vectors and tensors, the 'Unit' column shows the unit of their norm values.
continuum hypothesis.This means that any infinitely small volume element in the fluid is seen as a continuous medium, also called a fluid parcel.As Landau and Lifshitz stated [6], a fluid parcel is 'very small compared with the volume of the body under consideration, but large compared with the distances between molecules'.
In fluid mechanics, the continuity equation describes the transportation of physical properties in space and time as ∂A (x, t) ∂t + ∇ • (A (x, t) u (x, t)) = s (x, t) , where A can be an arbitrary scalar, vector or tensor physical property, u is the velocity, and s is the source term for A, all described at time t and location x.Equation (1) states that the change rate of any physical property at a fixed position ∂A/∂t depends on the variation brought by the flux of Au and the source term s.
Lagrangian and Eulerian viewpoints.Considering the physical attribute A in Eqn.(1), a flow field can be analyzed from a Lagrangian or an Eulerian viewpoint, as follows.
The Eulerian viewpoint studies the physical field using fixed positions.The change rate of the physical value A at a given position x is the ∂A (x, t)/∂t term in Eqn.(1), which is caused by both the flux and source terms.While intuitive, this does not explicitly express the motion of the fluid parcel in the continuum hypothesis as parcels constantly travel through fixed locations at all times.
In contrast, the Lagrangian viewpoint studies the change rate of physical attributes with respect to the fluid parcel by Physics-based where D (•) /Dt, the so-called material derivative, is the change rate of A within a fluid parcel.In Eqn. ( 2), u and s are respectively the velocity, and the source term, of a specific fluid parcel.Hence, all positions x can be substituted with the parcel identifier ι.For brevity, we next omit the explicit mention of (x, t), (t), and ι unless required by the context.

Navier-Stokes equations
Numerous methods for calculating fluid motion have been developed, spanning from Lagrangian to Eulerian perspectives.However, the underlying physical principles for nearly all these approaches are rooted in the Navier-Stokes equations, which govern the dynamics of fluid flow and serve as a fundamental foundation for fluid simulations.We describe these briefly next.
Mass conservation.In a closed system, fluid mass is conserved over time.This principle is represented by the continuity equation (Eqn.( 1)).By letting A be the fluid density and setting s ≡ 0, Eqn.(1) can be written as For the case of incompressible flow, the density variation within the flow is conserved, i.e., Dρ/Dt = 0.This condition further implies a divergence-free velocity field, as expressed by where s m is the momentum source altering the speed of each fluid parcel and ⊗ represents the outer product.Following this, a basic form of the Navier-Stokes momentum equation for viscid compressible flow further specifies s m into three separate terms as where p is the pressure, g is the gravitational acceleration, and µ is the dynamic viscosity coefficient describing how viscous a fluid is.Equation (6) says that the velocity change rate for a fluid parcel is affected by three force terms given by pressure (−∇p), viscosity (µ∇ 2 u), and gravity (ρg).

Early developments
As computer technology advanced in the 20th century, numerical methods became popular for solving partial differential equations including the Navier-Stokes equations.With the advent of powerful computer hardware and software, computer graphics began to incorporate these physics-based algorithms, enabling increasingly realistic fluid simulations.
Dating back to the 1970s, William T. Reeves, a member of Lucasfilm's Computer Division, Computer Graphics Group, pioneered the development of 'particle systems' [1,7].These systems enabled the realistic depiction of elements such as smoke and fire in films, as seen in "Star Trek II: The Wrath of Khan."This breakthrough laid the foundation for early fluid simulation techniques in computer graphics.In the 1990s, physics-based fluid simulation began to gain traction.Wejchert and Haumann [8] used a simplified version of the Navier-Stokes equations to animate irrotational, incompressible linearized fluid flow, providing the physics foundation for their fluid animations.Subsequently, Stam and Fiume [9] incorporated the complete Navier-Stokes equations to create turbulent wind effects.
On the Lagrangian side, Desbrun and Cani [10] introduced Smoothed Particle Hydrodynamics (SPH) to the computer graphics field for simulating highly deformable bodies.On the Eulerian side, Foster and Metaxas [11] used the Navier-Stokes equations on fixed grids to simulate fluid motion.The study of fluid simulation reached a significant milestone at the end of the 20th century with Stam's Stable Fluids method [2].This finally made stable, three-dimensional, physics-based fluid simulation an attainable goal, capable of producing realistic fluid effects.It was the first unconditionally stable method for fluid simulation and introduced the concept of semi-Lagrangian advection, and was one of the earliest works to apply the idea of hybrid simulation into the field.
Hybrid methods in fluid simulation merge the strengths of both Lagrangian and Eulerian approaches, yielding more versatile and robust systems.Two foundational principles that underpin hybrid fluid simulation are Harlow's [12] particlein-cell (PIC) method and the refined fluid implicit particle (FLIP) method of Brackbill and Ruppel [13].These techniques have significantly contributed to the widespread success and adoption of hybrid fluid simulation in the 21st century.The field also saw a major advancement when Zhu and Bridson [14] applied the FLIP method to simulate incompressible flow.This moved hybrid fluid simulation to new heights as it enabled exploring complex fluid dynamics with enhanced precision and stability.The continuous evolution of hybrid fluid simulation techniques has had a profound impact on computer graphics, facilitating the creation of realistic and visually stunning effects.

Discretization Strategies
As fluid simulations in computer graphics have evolved since the early development in the 20th century, the field has branched out into three distinct categories: Eulerian schemes, Lagrangian schemes, and hybrid schemes.Each of these approaches offers unique advantages and challenges, contributing to the comprehensive understanding of fluid dynamics in CG.
Eulerian schemes.These simulation methods use the Eulerian viewpoint introduced in Sec.3.1, i.e., compute property values at fixed points in the simulation domain.For this, the domain is typically divided into evenly-distributed cells.In a traditional collocated grid structure (Fig 4a), all physical values are evaluated at the center of each cell.To derive a continuous flow field with values at arbitrary positions, e.g., the grey dot in Fig. 4c, one can use a weighted interpolation of neighboring cell values.The staggered grid (Fig. 4b) stores physical values at cell edges and cell centers separately.Compared with collocated grids, staggered grids are currently more popular for simulating incompressible fluids given their higher stability.It is noteworthy that staggered grids are related to the marker-and-cell (MAC) method [15] which was used in the early days of computational fluid dynamics to solve incompressible flow problems.
Lagrangian schemes.In the Lagrangian framework, domain discretization is done by a set of particles moving with the fluid flow, each approximating the physical values of a fluid parcel.Hence, Lagrangian schemes conserve mass by construction.Since particle locations can be much more freely distributed over the computational domain than covering the same domain with a grid, Lagrangian schemes are also good at modeling complex free-surface details.
Currently, SPH is one of the most popular Lagrangian methods for fluid simulation, with origins in the works by Lucy [16] and Gingold and Monaghan [17].SPH has evolved significantly over time, with various advancements and improvements.
Figure 5 shows how SPH performs interpolation: Physical value A at the location x i of particle i is computed by using a smoothing kernel W as where h is called the smoothing length, V is the volume of (the parcel of) each particle, and j indicates all particles closer to i than the distance h.To compute higher-order quantities, e.g., pressure gradients, one can simply replace the kernel function W in Eqn. ( 7) by its higher-order counterpart (gradient in our example).Initially, the weakly-compressible SPH (WCSPH) [18] approach was introduced, where pressure computation was performed explicitly.Later, the Predictive-Corrective Incompressible SPH (PCISPH) [19] method was proposed, which introduced a prediction-correction scheme for implicit pressure computation.This technique improved the stability and accuracy of fluid simulations by enforcing incompressibility more effectively.Further developments led to the introduction of Implicit Incompressible SPH (IISPH) [20], which provided a more strictly incompressible simulation with increased computational efficiency.Most recently, the Divergence-Free SPH (DFSPH) [21] method has been developed, which further enforces the divergence-free condition within a simulation.Position Based Dynamics (PBD) is a versatile and efficient simulation method for handling various physical phenomena, including fluids, deformable solids, and cloth.PBD was first introduced by Müller et al. [22] as an alternative to traditional force-based dynamics, focusing on the direct manipulation of object positions instead of computing forces and accelerations.In the context of fluid simulation, the Position Based Fluids (PBF) method was proposed by Macklin and Müller [23], which builds upon the principles of PBF and enforces incompressibility by iteratively adjusting particle positions.
Hybrid schemes These schemes combine the advantages of Lagrangian and Eulerian schemes by representing the motion of the fluid flow with Lagrangian particles while computing dynamics (forces) on an Eulerian grid.
As Figure 6, shows, to combine particles and grids, physical values need to be mapped from particles to grids (P2G) and from grids to particles (G2P), before and after the dynamic simulation separately.A so-called shape function, similar to the kernel function W for the SPH method, performs these mapping procedures.
In the original PIC [24], only the momentum term is transferred between P and G.The later FLIP [13] transfers the differential of momentum to obtain better dynamic effects at the cost of stability.The Material Point Method (MPM) introduced by Sulsky et al. [25] is another extension of the original PIC.It adds a new dimension to fluid simulation by considering the deformation gradient information along with the momentum term, making it suitable for simulating a wide range of materials, including fluids, granular materials, and deformable solids.
Throughout the development of the PIC, FLIP, and MPM methods, these techniques have evolved and merged to form more advanced approaches.The Affine Particle-In-Cell (APIC) method [26] extends the MPM framework by incorporating affine velocity fields, which reduces numerical dissipation and offers improved stability compared to both PIC and FLIP.The Polynomial Particle-In-Cell (PolyPIC) [27] method takes the MPM framework a step further by incorporating higher-order polynomial velocity fields, building upon the advancements made by the APIC method.Finally, the Moving Least Squares Material Point (MLS-MPM) [28] method utilises moving least squares for grid interpolation and differentiation in MPM simulations, further enhancing the accuracy and robustness of the approach.
p Fig. 6 Particle-In-Cell method, hybrid scheme.This example uses a cell-centered grid so information is stored at the yellow points.During simulation, momentum and weight (which can be used to get the velocity u t I on the grid) are transferred from particles to the cell centers.In the next step, forces are applied to grid nodes to compute the new velocity u t+1 I .Finally, the velocity is re-transferred from the grid to particles.When particles get their new velocities, the new positions can be easily found by forward Euler integration.

Advanced computational approaches
Fluid simulation requires a high discretization resolution to reach high visual quality.Yet, more discrete particles or denser grids ask for more computing resources.This section surveys approaches from recent years for improving computational efficiency.We organise these into approaches that use adaptive time and/or space sampling (Sec.4.1), GPU or CPU parallelisation (Sec.4.2), and the more recent datadriven approaches (Sec.4.3).For a more extensive survey about this area, we refer to the work of Manteaux et al. [29].

Adaptive solutions
A stable, sufficiently accurate, and detailed simulation requires adequate temporal and spatial resolution.Time steps must be short enough to ensure stability, and high-resolution grids or dense particles are needed to capture fine details.However, computational cost increases with both spatial and temporal resolution, and an overall high resolution is not always needed.For example, time steps must be small for violent motion but can be longer when the overall movement is slow; high spatial resolution is needed to capture delicate splashes and sprays, but is less important deep inside the fluid where such detail is not visible.
As such, adaptivity uses high resolution only in the needed time and space instances and uses low resolution elsewhere to preserve accuracy and detail while reducing computing costs.Figure 7 illustrates this strategy with particles as an example.Adaptive methods can be categorized into temporal and spatial adaptivity.Temporal adaptivity dynamically changes the time step, either globally or locally for different parts of fluid.Spatial adaptivity adjusts the resolution for different fluid regions, or changes the method of discretization for a similar effect.These two approaches are described next.

Temporal adaptivity
Temporal adaptivity adjusts the time step length dynamically.A straightforward strategy is to adapt the time step globally, i.e., use the same time step for the entire simulation domain.The time step size is determined dynamically at each time step under a restriction.For further performance gains, different time steps can be used for different spatial domain zones, thereby reducing the total number of integration steps needed.

Global time step
The Courant-Friedrichs-Lewy condition (CFL) [31] is a well-known method for determining the time step size.Most current simulation methods compute a global time step according to the CFL condition at each time step.Generally, the CFL condition takes the form where ∥u c ∥ is the speed of information Determining an optimal value for C max often involves an extensive trial-and-error process tailored to a specific scenario.Sun et al. [32] addressed this issue in their work, where they considered metrics related to stability of MPM simulations, such as the deformation gradient.By using these metrics, they were able to more effectively identify the performance limits and improve the overall stability of simulations.
Asynchronous time integration.When dealing with scenarios involving both intense waves and calm regions, implementing a global time step restriction can be inefficient and wasteful.To address this, the concept of regional time stepping was initially introduced in the SPH method by Goswami et al. [33].This approach subdivides the simulation space into smaller regions, allowing each region to have its own independent time step.Recognizing the grid-based nature of the subdivided regions, Fang et al. [34] extended this idea to the MPM method.In their technique, a scheduler determines the order in which blocks are updated, while a buffer block is employed to handle boundaries between blocks with different time steps.This resulted in significant performance improvements, achieving speed-ups of 9.8 compared to traditional synchronous MPM implementations.Inspired by Fang et al. [34], Koike et al. [30] proposed an asynchronous time integrator for Eulerian liquid simulation, with an interpolation strategy to deal with boundaries between different-time-step zones, and an advection scheme to prevent seams at the boundaries.separate regions, they still necessitate synchronisation for all regions at simulation time barriers.Reinhardt et al. [35] presented a fully asynchronous time integration model for SPH fluid animation where each particlehas an individual time step and is processed using a priority queue.

Spatial adaptivity
These methods change the spatial resolution or the discretization method in different spatial regions so as to keep fine detail in some regions but use coarser (thus faster to compute) detail in less important regions.Spatial adaptivity methods are heavily dependent on the underlying discretization.We next detail different spatial adaptive approaches for Eulerian, Lagrangian, and hybrid methods.Eulerian grids.Grid-based Eulerian methods use adaptive grid structures to achieve dynamic spatial resolution.However, compared to uniformly distributed Eulerian grids, designing a stencil on an adaptive grid for pressure solving to attain highorder accuracy and forming a symmetric positive-definite linear system that can be efficiently solved on non-symmetric adaptive grids poses challenges.
The octree data structure is one of the approaches for grid adaptivity that allows for changing the resolution of axis-aligned structured grids.As shown in Fig. 9, each cell is divided into eight equal children by cutting it in half along each axis .Octrees have the advantage of regularity, support fast discretization, and are simple to implement.However, on the transition between different grid levels, octrees have T-junctions which cause challenging numerical issues.
Losasso et al. [36] proposed the first octree-based liquid solver using a set of symmetric differential operators, which enables solving the Poisson equation on unrestricted octree grids.In the octree, velocity is stored on cell faces while pressure is stored at cell centers.The velocity divergence ∇ • u at cell centers is computed considering all cell faces f as where V c is the cell volume, and n f , u f , and S f are the outward-pointing normal, velocity, and area of face f , respectively.The pressure gradient on each face is computed from the pressure of the two adjacent cells using where ∆x denotes the cell size; and subscripts 1 and 2 denote the adjacent cells to that face.Dynamically adjusted octree grids also present challenges in terms of modifying and accessing data.Setaluri et al. [37] proposed a sparse paged grid (SPGrid) data structure that constructs the octree as a hierarchy of sparsely populated regular grids instead of a standard pointer-based tree.Goldade et al. [38] recognized the limitations of the first-order accuracy of the velocity field for octrees and contributed a variational finite difference discretization method to it, enabling a more efficient viscous simulation.Ando and Batty [39] focused on using octree grids to enhance surface detail exclusively.This approach further reduces implementation complexity while retaining the benefits of octrees.While the particular attention to maintaining data order for efficient computation is advantageous, it also presents a challenge in system design.Shao et al. [40] noted this issue and identified an underutilized potential within the regular Cartesian grid structure.
They ingeniously integrated the single instruction, multiple data (SIMD) approach with a multigrid structure, aiming to streamline and minimize the needed multiplications.Their method showed significant speed-ups of 2.0 to 14.6 times compared to contemporary adaptive octree solvers found in commercial software for large-scale simulations.
Several approaches have been inspired by and extended from the octree grid concept.These works aim to improve efficiency and accuracy in various ways.Ferstl et al. [41] proposed a hexahedral finite element discretization multigrid solver on adaptive octree grids.By specially treating boundary conditions on the free surface, they achieved second-order accuracy on the surface.Aanjaneya et al. [42] focused on enhancing pressure projection on octrees.They used a finite volume power diagram to accurately recover irregular embedded boundaries that cross grids, reaching both second-order accurate and symmetric positive definite (SPD) conditions.Xiao et al. [43] introduced an adaptive staggered-tilted (AST) grid for conducting adaptive fluid simulations on a regular discretization.By adding a tilted grid to an octree structure, they avoided T-junctions and further improved the adaptivity of the simulation.Some methods in fluid simulation employ multiple grids with different resolutions or structures to simulate various parts of the fluid, later compositing these elements together.This approach contrasts with using a single adaptive grid for the entire fluid domain.Gao et al. [44] devised a technique that divides the domain into nested partitions with different resolutions, effectively handling multi-resolution fluid behavior.English et al. [45] used overlapping Cartesian grids with varying scales and rotations to represent the fluid domain, constructing a local Voronoi diagram for managing pressure projection near grid interfaces.Li et al. [46] introduced an adaptive relaxation method for kinetic approaches, enabling fluid sampling at arbitrary overlapping resolutions and providing efficient representation of fluid behavior across a wide range of scales.
While many adaptive methods that use single or multiple grid structures can disrupt the uniform data structure of the original Cartesian grid, some works have managed to strike a balance between maintaining uniformity and introducing adaptivity.Zhu et al. [47] used a uniform grid within a cubical region of interest, extending the grid into the far-field by stretching cells along an axis.This approach retains the benefits of a uniform grid while providing adaptivity in specific areas.Ibayashi et al. [48] proposed a technique for dynamically warping uniform grids, combining the advantages of both unstructured and structured grids.
Lagrangian methods.Particle-based Lagrangian approaches, such as SPH, achieve spatial adaptivity through defining a desired resolution for each particle with a sizing function.By adjusting particle sampling through local merging or splitting of particles (as shown in Fig. 7), these methods are able to dynamically change the resolution, offering more efficient and accurate simulations while focusing on areas of interest.
The early study of adaptive SPH can be traced back to the work of Adams et al. [49].They introduced a sizing function based on geometric local feature size that allows focusing computational resources on geometrically complex regions.Yet, adaptive particles yield density errors due to different resolution scales, which can lead to instabilities.To address this issue, Orthmann and Kolb [50] proposed a temporal blending technique to limit the temporal rate of resolution change, thereby significantly reducing the error.With the advent of more strictly incompressible implicit SPH approaches, the size difference between neighbouring particles needs to be kept even lower to avoid instability.Winchenbach et al. [51] achieved this by forming a continuous transition of particle resolution by tuning the splitting and merging pattern, and introducing mass redistribution between particles.A simplified version of temporal blending is also incorporated.In their method, splitting supports arbitrary 1 : n patterns.Merging uses an (n + 1) : n pattern where one particle is merged into the others.Mass redistribution divides the excessive mass m ex of one particle i equally among n particles.The physical attributes A of the mass receiving particles j are updated to A * j by Zhai et al. [52] took inspiration from this method to propose an adaptive scheme for Power Particles.structure for efficient neighbour search and ray tracing for adaptive SPH based on hash-maps.Among vortex methods, ways to solve the Poisson problem efficiently with adaptive data structures were also studied.The Poisson problem is an N -body problem, where the interaction between each object and the rest of the objects is considered.Naively solving this problem requires O(N 2 ) computations, so adaptive methods are used to reduce this complexity.The Fast Multipole Method (FMM) [57] uses an octree to solve the N -body problem approximately in O(N log η N ) where η ∈ {0, 1} by approximating interactions between far-away bodies by the body centers instead of computing all pairwise interactions.Zhang and Bridson [58] proposed a novel Particle-Particle Particle-Mesh method, which is easier to implement and parallelize on GPUs, and applied it on a vortex segment solver.Angelidis [59] used FMM with added support for non-uniform particle sampling to simulate incompressible smoke with vortices.combination of Lagrangian and Eulerian representations for fluid simulation.Ando et al. [63] applied the particle splittingcollapsing scheme, similar to traditional Lagrangian adaptive mechanisms, to adjust the granularity of fluid representation in relation to the distance to the fluid surface for FLIP.They used the finest particles to represent splashes and sheets.Ando et al. [64] later introduced an adaptive liquid solver on tetrahedral meshes, which combined a variant of FEM with FLIP advection, adapting both the size of the particles and that of tetrahedral meshes coherently for more efficient simulation.In highly stable situations, Yue et al. [65] explored the possibility of simulating the interior area as soft continuum materials to reduce computational costs at the solver level.
To further save on computation costs, more recent hybrid approaches aimed to 'hollow' the inner area of the simulated fluid by using Eulerian simulation only, with particles applied near the surface.Chentanez et al. [66] proposed coupling pure Lagrangian and Eulerian methods to simulate single fluid bulks, addressing the issue of fluid representation transition and coupling stability between two different fluid solvers.However, coupling two different solvers can still be prone to instability, so Ferstl et al. [67] later returned to adaptive FLIP simulation using FLIP particles within a narrow band of the fluid surface.Nakanishi et al. [60] focused on constructing a proper data structure to adaptively merge and split octree grids for FLIP, adapting the size of the background grid only near the fluid surface (Fig. 10 ).[67] to replace some of the fluid surface using level set surfaces via an introduced transition function.This idea represents a breakthrough in the field, enabling a seamless merging of physical simulation details and largescale representations.To address the setback of inefficient, highly-dissipative wave propagation during the transition of surface representation, Huang et al. [69] employed hybridization of volumetric and surface-based advection-projection discretizations, with the Boundary Element Method (BEM) applied for long-lasting waves.A single processing unit can be a multicore-CPU or a GPU.A CPU has few, but high-processing-power, cores.A GPU has weaker, but many more cores which can result in better performance for repetitive operations.

Parallelization
Parallelization of simulation algorithms offers a promising pathway to augment fluid simulation speed by capitalizing on the existing parallel computation capacities of modern GPUs and CPUs.We divide such methods into three classes: parallelization on a single processing unit (multi-core CPU or a GPU); parallel techniques with multiple processing units, especially multiple GPUs; and using distributed systems (see also Fig. 13).

Single processing unit
Parallelization can be most easily achieved on a single processing unit such as a multicore CPU or a GPU.Pure Eulerian and Lagrangian simulation techniques can be readily parallelized under these conditions.However, hybrid methods necessitate meticulous measures due to their composite nature of particles and grids.The Voxel Data Block (VDB) method pioneered by Museth et al. [70] offers a robust framework for manipulating sparse volumetric data.Wu et al. [61] successfully integrated this approach into the FLIP method by resolving the parallel particle-to-grid rasterization challenge inherent in hybrid methods (Fig. 11).Gao et al. [71]accomplished a similar integration for MPM based on the SPGrid structure [37], wherein the particle-to-grid transfer mechanism was a crucial aspect that needed resolution.Furthermore, Chu et al. [72] applied and optimized the Schur complement theory for FLIP simulations.They divided the simulation domain into multiple subdomains with face edges and cross-points to establish a parallel-friendly data structure.Aiming to streamline the development of parallel programs, Hu et al. [73] proposed a new data-oriented programming language -Taichi.This language facilitates efficient authorship, access, and maintenance of sparse data structures, and is accompanied by a compiler designed to automatically optimize and parallelize code on CPU or GPU platforms.Subsequently, Hu et al. [74] enhanced Taichi to include bit-level memory control over numerical data types.

Multiple processing units
Multi-GPU techniques are the approach of choice for scaling up the simulation size.In this case multiple GPUs cooperate with the CPU(s), together forming a heterogeneous computing structure.The most challenging point for this structure is to reduce both the frequency and the amount of the data exchange between CPU-GPU and GPU-GPU, which is very expensive.Liu et al. [75] initially introduced the Schur complement method to the graphics community as a strategy to tackle this issue.As previously outlined, the Schur complement method has the significant benefit of segmenting the entire simulation domain into several regions.Each region can be efficiently computed using a single GPU, with the interaction boundary between these regions elegantly managed by the CPU(s).This setup eliminates the need for data transmission from a GPU to the CPU and subsequently to another GPU.Meanwhile, Wang et al. [76] took a different approach by directly adapting the MPM algorithm structure to fit a multi-GPU framework.They accomplished this by developing a particle data structure to encourage coalesced memory access and circumvent atomic operations during particle-to-grid data writing.Additionally, they proposed a kernel fusion strategy to diminish the number of GPU kernel launches and reduce global memory requirements.Chen et al. [62] further optimized kinetic methods for scenarios containing complex solids (Fig. 12).They introduced a multi-kernel launch methodology for parallelism enhancement and a parametric cost model to improve performance optimization.This exploration of fluid simulation parallelization demonstrates the richness and versatility of strategies within this domain.

Distributed systems
In the pursuit of scalability in fluid simulation, the potential of distributed platforms is harnessed to delegate tasks across multiple computing nodes.These distributed simulations frequently utilize automatic task allocation to ensure efficient processing.Biddiscombe et al. [77] laid a crucial Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 13 groundwork by providing a GUI-based interface and analysis system for high-performance computing (HPC).Their innovative approach involved substituting the I/O layer in the Hierarchical Data Format Version 5 (HDF5) with a parallel data transfer driver, thereby enabling parallel simulation, analysis, and GUI operation concurrently on one or multiple devices.Mashayekhi et al. [78] developed a system that automatically distributed tasks across many multi-core cloud computing nodes to dynamically manage fluid partitions.Shah et al. [79] proposed a load balancing scheme for sparse fluid simulations.Qu et al. [80] outlined a simple yet effective solution to accelerate distributed fluid simulations which uses micro-partitioning to greatly improve load balance and communication performance.

Data-driven approaches
Although great progress has been made in space-time adaptive and parallel computing in recent years, the simulation of fluid by traditional physical methods still requires high computational resources.Strictly limited time steps are needed to ensure simulation stability when solving governing equations in a discretized space.Implementing efficient and accurate end-to-end fluid simulation pipelines is also technically highly challenging.Data-driven approaches provide an alternative solution for real-time interactive fluid simulation, which we outline next.

Model reduction
Model reduction is achieved by precomputing a set of simulation sequences to obtain a low-dimensional representation of fluid motion, allowing for efficient and fast re-simulation.
Treuille et al. [81] first introduced model reduction to fluid simulation via Galerkin projection.They constructed vector field basis functions for fluid dynamics based on principal component analysis (PCA) to generate real-time fluids.For this, they computed the Galerkin projection of the differential equation F onto the reduced dimensional space as The high-dimensional space vector v ∈ R n and the lowdimensional space vector r ∈ R m are transformed by the projection operator P (v) = r and its inverse Later, improvements such as the extension of the Galerkin projection to non-polynomial systems [82] and a multidimensional cubature method supporting semi-Lagrangian advection [83] were proposed.De Witt et al. [84] used Laplace eigenfunctions instead of PCA eigenvectors.The method performs a Galerkin projection of the vorticity form of the Navier-Stokes equations and is therefore not data-driven but physically driven.Liu et al. [85] further stabilized the method using a variational integrator, providing structure coefficients without artifacts.Zhai et al. [86] proposed a model reduction method based on empirical modal decomposition (EMD), which can decompose the flow field into various frequency components as the basis vectors for model reduction.The method can extract the characteristic parameters of the original fluid to achieve inverse modeling.To reduce the memory needed to store basis functions, Cui et al. [87] generalized the dynamics to Neumann boundary conditions using analytic eigenfunctions and the Fast Fourier Transform (FFT).This approach allows the use of thousands of basis functions to produce more convincing and fine-grained fluid dynamics.Using this, they proposed an analytical extension of the Laplace eigenfunction method [88].This spiral-spectral fluid simulation method is capable of producing realistic turbulent effects over a variety of radial domains, both surface and bulk.Mercier and Nowrouzezahrai [89] constructed an anisotropic vector field basis function that can accommodate curved boundaries and coupling with dynamic obstacles.The method sacrifices a physically accurate solution for a visually plausible simulation.A reduced model for fluids based on incompressible polynomial vector fields was proposed by Panuelos et al. [90] to reduce the computational cost of highly viscous fluids.

Machine learning
Machine learning methods have brought about a revolution in physics-based fluid simulation, with deep learning techniques particularly proving the possibilities of data-driven approaches.Ladicky et al. [91] expressed physics-based fluid simulation as a regression problem and used a regression forest algorithm to approximate the dynamic behavior of fluid particles.This method strongly generalizes to simulating large-scale scenes in real time.Raveendran et al. [92] used interpolation on existing fluid simulations to rapidly generate a large number of new simulation results.Thuerey [93] improved this method using a signed distance function to fully automate the matching process.Recently, Oh and Lee [94] proposed a temporal interpolation network based on optical flow and forward advection that can derive high frame rate smoke simulations from low frame rate simulations.
With the development of deep learning, Convolutional Neural Networks (CNNs) [96] and Artificial Neural Networks (ANNs) [95] were introduced to solve pressure calculations in fluid simulations and accelerate the pressure projection step (Fig. 14).Wiewel et al. [97] proposed an LSTM architecture to predict the evolution of fluids over time.They used CNNs to map the 3D fluid simulation to a low-dimensional latent Fig. 15 The position information of the fluid is extracted from example videos as a reference (left).The 3D simulation is projected onto the screen space (middle).The difference between the two is next iteratively reduced (right) [99].
space, thus greatly speeding up the simulation.They further improved the method using latent space subdivision [98], allowing for more stable predictions of complex long-term series.Takahashi and Lin [99] proposed a framework capable of extracting physical parameters from real fluid videos and applying them to new scenarios to generate the user's ideal fluid behavior (Fig. 15).Eckert et al. [100] created ScalarFlow, the first large-scale volumetric dataset for real smoke reconstruction for computer graphics and machine learning.ScalarFlow also makes an important contribution by providing reliable benchmark data and evaluation criteria.

Fluid coupling with multi-materials
A key topic frequently mentioned in fluid simulation is the coupling of fluids with their environment composed of different materials.Indeed, in computer graphics, fluids are mostly attractive due to the way they interact with their surroundings, as the attention of the spectator is arguably attracted by the interfaces, or boundaries, between fluids and the rest of the vir-tual world.Moreover, the behavior of fluid is strongly affected by how these surrounding factors themselves evolve.In this section, we explore this topic with the focus on recent works that aim to accurately and efficiently model the coupling with multiple complex materials.We split the discussion into three separate subtopics: meshless methods (Sec.5.1), mesh-based methods for handling fluid boundary conditions in the case of solid boundaries (Sec.5.2), and solutions designed to model more complex couplings emerging for multiple boundaries (Sec.5.3).

Meshless Methods
Particle based boundaries For most Lagrangian simulation approaches, the solid boundary sampled by so-called boundary particles is the main enabler of inter-particle interactions between fluid and other objects (Fig. 16).To couple fluid with solid boundaries, early methods used various approaches such as collision detection and ghost particles.Becker et al. [101] computed the contact point between the fluid and solid particles and controlled the normal and tangential velocities to impose boundary conditions.Yang et al. [102] facilitated the interaction between SPH fluid and nonlinear FEM deformable solid by sampling proxy particles across the boundary, and handled fluid-solid coupling with momentumconserving collisions.Schechter and Bridson [103] used the ghost particle method to generate a thin layer of ghost particles in the nearby solid and air, reducing numerical errors caused by non-uniform particle distributions near boundaries.He et al. [104] generated staggered particles between neighboring particles to resolve the problem of shape functions losing the Kronecker delta property, thereby supporting various slip boundary conditions.
∂Dx 16 Schematic diagram of fluid-solid coupling using meshless methods.Left: Discretization of the SPH approach in the coupling with materials.Fluid and solid particles both contribute to the boundary handling.Ω denotes the domain, f denotes fluid, s denotes solid, ∂Dx denotes the surrounding spherical neighborhood.Right: Solid domain in the neighbourhood of a fluid particle represents the contribution of the boundary density value at the particle position.
To further reduce computing time and numerical errors, Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 15 Akinci et al. [105] proposed a versatile and efficient method for SPH boundary handling without the need for collision detection or generating extra particles.The method resolved to directly handle the problem of uneven boundary sampling by evaluating a relative contribution of each boundary particle using Shepard interpolation (Fig. 17) given by to compute the correct density without depending on boundary sampling.In the subscripts, f and s respectively represent fluid and solid particles, ρ 0 denotes the fluid rest density, and Θ s computes the contributions of a boundary particle according to its volume.Compared to the volume of fluid particles with the same size V f , the volume of solid particles V s varies based on the local solid sampling distribution.Denser sampled regions have smaller solid volumes to maintain interaction stability.At the mere cost of a one-time evaluation of the above procedure, thin boundary geometries with only one layer of particles and non-manifold geometries can be supported.
Fig. 17 Large-scale fluid-solid interaction [105].A boat with ragdolls sails under a bridge (left).As the flow rate increases and the bridge is released, a second boat impacts the bridge.
The method of Akinci et al. [105] has gained widespread popularity due to its simplicity and efficiency.Macklin et al. [106] applied this method to unify various physical behaviours using a single particle system for real-time applications built on the PBD.Cornelis et al. [107] used the concept of [105] to couple high-resolution FLIP with a low-resolution implicit SPH method.Peer et al. [108] built an implicit formulation for the simulation for the incompressible linearly elastic solids embedded in the ISPH pressure solver which further enables a pressure-based boundary treatment using the method of [105].Takahashi et al. [109,110] integrated the approach of [105] into their multilevel particle-based solver, in which they adaptively assigned various roles to the particles, to guarantee the solvability of the linear system in a unified manner regardless of the arrangements of the particles.
Although Akinci et al. [105] eliminated a variety of artifacts of particle-based fluid-rigid coupling, numerical issues, such as penetration across the boundary, and the lack of higherorder accuracy of pressure due to the mirroring scheme of physical values from fluid to boundary particles, still limit the time step size and stability under drastic scenarios.Shao et al. [111] treated surface and inner boundary particles differently to prevent particle deficiency and penetration issues.Instead of using the state of density compression denoting the magnitude of the pressure force, Band et al. [112] introduced the notion 'volume compression' into coupling problems, where the rest volumes rather than the mass of boundary particles are applied to derive a continuous pressure force: which considers fluid samples f i with all fluid and boundary neighbours j of the fluid sample f i in the same way.The results of their experiments show a significant improvement of stability and a wider range of possible time step sizes.Gissler et al. [113] resolved the stability issue in [105] by interlinking an artificial pressure solver for boundary particles with the fluid pressure, achieving a fully dynamic two-way coupling (Fig. 18).Truong et al. [114] prevented penetration by treating particle collisions with particle merging and splitting.
Fig. 18 SPH fluid (43.8M particles) in a terrain (50M static rigid particles) is two-way coupled with a water wheel that is connected to a gate via gears and a chain [113].The gears, chain, and water gate are modeled with 2.3M dynamic rigid particles.Up to 90k simultaneous rigid-rigid contacts are handled.
Unsampled boundaries.The influence of a solid model upon the surrounding fluid particles can also be handled carefully using a mesh-based (Eulerian) solid representation for a stable and effective coupling.Vines et al. [115] coupled Lagrangian vortex particles to mesh-based solids by generating vorticity at the solid boundary.Fujisawa and Miura [116] considered the influence of triangle mesh boundaries on the integration of a kernel function for SPH without the need for boundary particles.However, the method cannot handle solid boundaries with complex geometry and is computationally expensive compared to [105].Chang et al. [117] extended [116] to support arbitrarily shaped solid boundaries by converting the volume integral inside the solid boundary to a surface integral.Koschier and Bender [118] presented the 'density maps' method, precomputing a continuous boundary density field to efficiently handle arbitrary boundary geometries.Bender et al. [119] also targeted the expensive renormalization process in [116] by storing the volume contribution from the boundary on a spatial grid that can be efficiently queried at runtime.

Mesh-based Methods
Another mainstream approach for fluid simulation is to use Eulerian and Lagrangian meshes as shown in Fig. 19.Yet, loss of mass and difficulties in handling drastic deformations often limit the practicability of these approaches.We organize mesh-based methods in several classes, as follows.
Lagrangian/Eulerian meshes Early on, Clausen et al. [120] used fully Lagrangian tetrahedral meshes to significantly reduce numerical viscosity in simulations with relatively low resolutions and long time steps.Azevedo and Oliveira [121] proposed a semi-Lagrangian method that introduced curvilinear grids and achieved more accurate boundary conditions in simulations with moderate resolutions.Besides these Lagrangian-meshed methods, Teng et al. [122] later incorporated the previous work into an Eulerian fluid solver and resolved complex contact scenarios between multiple solids and fluids.Recently, Takahashi and Lin [123] formulated the implicit viscosity integration as a minimization problem in which the volume fractions are consistently evaluated to handle sub-grid details.Focusing on the poor boundary conditions of irregular boundaries defined under coarse grids, the cut-cell approach became a major trend to improve convergence for Neumann boundary conditions.Fluid grid cells are clipped against the solid boundary represented by a triangle mesh, forming several distinct polyhedral sub-cells at each time step, allowing small details to be handled without refining or rotating the grid.Based on the multigrid scheme proposed by Chentanez and Mueller-Fischer [124] with a variational discretization compatible on all levels, Weber et al. [125] presented a cutcell-based multigrid scheme on staggered grids that is secondorder accurate for Neumann boundaries.To better capture the flow across thin solids and gaps, Azevedo et al. [126] further proposed a topology-preserving pressure projection scheme on cut-cell meshes.Following this, Zarifi and Batty [127] used cut-cell discretization for the two-way fluid-deformable interaction, enforcing the free-slip boundary condition at the actual interface.Their method computed the pressure based on the MAC grid of three dimensions x,y,z as where Ω contributes to the domain region and S(Ω) is the area of the domain region.The cut-cell is further extended by Chen et al. [128] to represent sub-grid structures on the free surface of the liquid, with a new iso-surface Poisson solver with desirable properties such as second-order accuracy and symmetric positive definiteness.Tao et al. [129] introduced an algorithm based on VEM for simulating fluid flow on cut-cell meshes, which effectively handles complex geometries and accurately captures intricate features, including thin tubes and extremely thin walls.

Material Point Method
Since its introduction to computer graphics, the Material Point Method (MPM) has garnered significant attention.By integrating features of Lagrangian particle representation and Eulerian grid representation, MPM offers a powerful technique for coupling fluid and solid simulations.However, conventional MPM solvers have drawbacks such as computational inefficiency and limited capability to handle self-contact collisions, despite their physical realism and geometric convenience.To improve upon these, Gao et al. [130] presented an adaptive Generalized Interpolation Material Point (GIMP) method with extensively optimized particle-grid transfer memory efficiency and parallelism.Hu et al. [28] proposed Compatible Particle-In-Cell (CPIC) which enables the handling of discontinuous material points and infinitely thin boundaries by leveraging the relative positions Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 17 of grid nodes and particles.They also embedded the Moving Least Squares (MLS) method into MPM to double the computation speed.However, such MPM approaches do not address the inconsistent tangential velocities at the interface between multiple materials, leading to visually unpleasant artificial stickiness.To alleviate this, Fang et al. [131] presented a ghost matrix operator-splitting scheme for monolithic coupling between incompressible fluids and elastic solids and designed a novel interface quadrature cut-cell MPM formulation for free-slip boundary condition.Subsequently, Cao et al. [132] extended some ideas of [131] from incompressible to compressible flow.
Monolithic schemes Monolithic solvers simulate various materials and their interactions within a unified system that includes boundary conditions.These schemes naturally ensure a more robust interface of large density ratios and enable large time steps.They not only occur in SPH methods [113] but also in the mesh-based (e.g., MPM) methods [131,132].Aanjaneya [133] proposed a monolithic solver for efficiently simulating the interaction between rigid bodies and incompressible fluids.The solver remains robust even in poorly conditioned scenarios with large density ratios between the solid objects and the fluid.Lai et al. [134] introduced a V-cycle of the Full Approximation Scheme (FAS) multigrid method to solve the linear complementary problems to achieve better scalability and efficiency compared to previous methods.Takahashi and Batty [135] proposed a monolithic pressureviscosity-contact solver to simulate the complex interactions between rigid bodies and liquids, efficiently managing incompressibility and offering the option for implicit viscosity integration in liquids.It also addresses contact resolution for rigid bodies and handles mutual coupling.
Partitioned schemes Compared to the monolithic scheme, partitioned schemes can deal more flexibly with multiple co-existing solvers by alternating between the solid and fluid while applying suitable boundary conditions.Akbay et al. [137] employed fluid and solid solvers as independent components with restricted interfaces, promoting modularity and facilitating code reusability.Lee et al. [138] used a partitioned methodology to connect a coarse background Eulerian grid with a fine ALE mesh in their simulation framework for character-water and hair-water interactions.
Recently, several coupling approaches for the special treatment of fluid have emerged.Brandt et al. [139] modeled fluids and deformable objects as incompressible media avoiding expensive operations such as interface tracking and boundary condition handling.Ruan et al. [136] used a three-way coupling method, employing a thin liquid membrane to model contact between solid objects and fluid driven by strong surface tension.(Fig. 20).

Coupling with complex boundaries
Besides interacting with the dynamic or static boundaries of surroundings composed of rigid or elastoplastic materials, fluids often interact with other solid matter with complex and diverse physical attributes.Such cases require intricate boundary models and corresponding solvers.Next, we discuss recent advances in fluid coupling with thin film surfaces and/or porous materials including hair, cloth, sponges, and sand.

Strands and cloths
Coupling between hairs and fluid is complex due to the wetting of hairs.Such phenomena are shown in Fig. 21.Rungjiratananon et al. [140] used an Eulerian approach to capture hair porosity and wetting effects, and a Lagrangian approach to simulate individual hair strands and their interactions, resulting in a detailed and dynamic hair simulation.Chen et al. [141] proposed a real-time painting system that aimed to generate realistic paintings by simulating the interactions among the brush bristles, paint, and canvas.Fei et al. [142] proposed a multi-component framework to model wet hair.PIC and Kirchhoff Rods were applied to model fluid and the hair separately.A height-field was introduced to represent the liquid volume around each hair strand, considering the wet condition.Fei et al. [143] next extended the fluid attributes to compressible, shear-dependent liquids.A modified second-order Coulomb cone model was also designed to capture cohesion and friction during strand collisions.Lee et al. [144] used a tetrahedral volume mesh to embed hair, enabling the hairs to adhere to their embedded positions, and facilitating simulations with millions of hairs during water-hair interactions.
For the interaction between cloths and fluid, Huber et al. [145] proposed an efficient method for two-way interaction between particle-based fluid and thin triangular meshes, enabling cloth-fluid coupling even at large time steps.Jiang et al. [146] created an anisotropic hyperelastic model that distinguishes the response to manifold strain, shearing, and compression in orthogonal directions.This model facilitates the coupling of various materials such as elastic surfaces, curves, fluids, and granular materials.Fei et al. [147] introduced a method for simulating the intricate dynamics of woven or knit fabrics, both partially and fully saturated, interacting with liquids, using the method of Jiang et al. [146] for contact and collisions.To simulate stain formation and evolution on cloths, Wang et al. [148] developed a pigmented solution by utilizing a homogenization process that combines inhomogeneous and/or anisotropic properties into bulk anisotropic diffusion tensors.Zheng et al. [149] formalized the spreading of stain in woven fabric into in-yarn diffusion and cross-yarn diffusion, and introduced a triple-layer model to manage wetting and wicking calculations.[150] simulated liquid flow within the porous object with a deforming unstructured mesh and modeled liquid diffusion based on saturation, as well as allowing the liquid to be absorbed by, or leave, the solid.

Sponge-like porous materials Patkar and Chaudhuri
Thin film surfaces.Real-life situations often involve surface flow phenomena, such as rainwater cascading down a tree trunk or the gradual progression of a water front in a shower room.For such flows, Vantzos et al. [151] proposed a triangle mesh model for simulating the motion of a thin viscous fluid film on a curved surface.The model includes discretization for curvature and advection operators to ensure accurate simulation results.Ren et al. [152] expanded the standard shallow-water flow model to accommodate general triangle meshes.They introduced a feature-based bottom friction model, allowing for the capture of non-viscous flow motion along edges and creases on detailed 3D meshes.Granular materials Such materials, e.g., sand, can themselves exhibit flow-like behavior.In addition to their porous behavior, simulating dissolution of granular materials also attracted much attention.Yan et al. [154] combined a hypoplastic model with SPH to simulate granular materials diffusing into fluid.Yang et al. [153] integrated the phasefield method to simulate liquids and multiple types of solids with dissolution achieved by evolving the granular particles' concentration and phase (Fig. 22).Tampubolon et al. [155] used continuum mixture theory to simulate water-sand coupling in MPM where different phases are coupled by momentum exchange.Gao et al. [156] modeled the motion of solid sediment particles inside fluids with MPM.Sediment was modeled by Drucker-Prager elastoplasticity with two-way coupling between sediment and fluid being achieved.He et al. [157] proposed position-based constraints for granular flows, using cohesion and friction models that vary across space, with cohesion affected by water saturation.Takahashi and Batty [158] simulated two-way coupling between rigid bodies and continuum granular materials or liquids with a monolithic solver that combines pressure, friction, and contact interactions.Gao et al. [159] used a hybrid scheme to accurately simulate the behavior of discontinuous fluid-like substance.This approach integrated an affine particle-in-cell solver with density fields, enabling transformations across granular particles, dust clouds, powders, and their mixtures within a unified framework.

Multiphase liquids
The real world is replete with complex fluid phenomena such as dissolution, dispersion, and Rayleigh-Taylor instabilities,

Non-mixing fluids
Recent methods for simulating non-mixing fluids include particle-based methods, e.g., SPH [160], PBF [161], MPM [162], and mesh-based methods [163][164][165][166]. Phases in non-mixing fluids typically do not merge together.Nonmixing fluid phases inherently resist fusion, creating challenges in (a) tracking the interfaces between different phases and (b) accurately calculating force interactions between phases at liquid-liquid or liquid-solid interfaces.The unique properties of gas-liquid interfaces are discussed separately in Sec. 7. SPH simulations of multiple fluids, particularly those with high density ratios, can produce erroneous interface tension and a non-physical separation between the fluids.In response to this issue, Solenthaler and Pajarola [160] modified the standard SPH equations to account for the density discontinuity across the interface.They used an SPH density interpolation formulation ρ i = m i j W (x i − x j , h) instead of the standard one ρ i = j m j W (x i − x j , h).
Alduán et al. [161] proposed a versatile simulation framework based on PBF methods designed to address VFX production demands.On the phase interface, they used an SPH density interpolation formulation which is similar to the method in [160].Subsequent PBF calculations were also modified using this uniform-density formulation.High viscosity was handled by breaking the XSPH calculation into multiple lower-viscosity stable iterations.The surface tension effect was bounded for stable artistic controls.
Fig. 24 Rayleigh-Taylor instability appearing at the interface between a high-density medium atop a low-density medium due to gravity [165].
Yan et al. [162] extended the MPM method to cope with interacting solid-fluid simulations.Upon the detection of a solid-fluid collision along the phase boundary, an antipenetration force is applied to both the fluid and solid particles in opposite directions along the surface normal.The same strategy is used to simulate non-mixing fluids where the antipenetration force is calculated between each distinct phase and all the other phases collectively (Fig. 23).
Mesh-based methodologies offer the advantage of explicit interface tracking using high-resolution geometric structures.Da et al. [164] considered a special emphasis on topologicalchange handling in 2D surface tracking for non-mixing fluid simulations (Fig. 23).They constructed highly distorted interfaces featuring thin sheets and tiny regions.Misztal et al. [163] , on the other hand, aimed to prevent mismatches between phase occupancy regions and the simulation quantity storage grid.They achieved this by discretizing each phase region using unstructured 3D tetrahedral grids, tracking the deformation of the tetrahedra.They managed topology changes and mesh quality enhancements using a modified 3D deformable simplicial complex method.In contrast, Li et al. [165] avoided complex remeshing operations by combining mesh-based tracking and surface reconstruction from distance fields (Fig. 24).They reconstructed the surface meshes between each phase using an unsigned distance function and an indicator function.The mesh was stored as later interpolation reference in a semi-Lagrangian updating of the two functions in the next time step.This approach was extended in [166] for surface tracking of more than three phases.In this method, special care is taken to handle mesh penetrations and ensure consistency between meshes and the regional level-set functions on the multi-fluid interfaces.An extended triangulation template strategy is also proposed to handle triple junctions which standard marching cube algorithms fail to reproduce.

Mixing fluids
Another category of multiple-fluid flows involves miscible or dispersed fluid mixtures where interfaces can be challenging to track continuously or may not exist at all.In contrast to the non-mixing case (Sec.6.1), different phases now always co-exist at the same spatial position (Fig. 25).A key problem is to calculate how the local volume fractions c V k of each phase k change during the simulation, i.e., solving the multiphase continuity equations.Other challenges include simulating diffusive behavior, incompressibility enhancements, and stability improvements.We summarize several typical methods which can be integrated into SPH [167,168], PBF [153,169], and MPM [170] solvers.Ren et al. [167] used a multiphase mixture model for complex multiphase mixing and unmixing effects.Their SPH-based approach (called WCSPH) solves on each particle the mixture continuity and momentum equations given by

Mixture model
where is the mixture velocity, p m is the mixture's pressure, and T m , T Dm are the mixture viscous stress and diffusion tensors, respectively.Phase velocities are assumed to be different from each other.For each phase k, a drift velocity u mk = u k − u m is analytically computed at the start of each time step.These phase-wise drift velocities are used to calculate the phase volume fraction change Dc V k /Dt as well as the mixture diffusion tensor T Dm .Following this, the aggregate particle motion, individual phase velocities, and the phase volume fraction changes on the particles during the simulation are solved.Yan et al. [154] extended this mixture model to cope with solid phases.Ren et al. [171] further introduced a virtual phase concept for multiphase simulations containing porous solids considering the absorbed and nonabsorbed parts of a single phase to be two virtual phases that can be universally handled by the mixture model.The result is a unified algorithm framework for multiphase flows inside and outside porous solids.To alleviate the incompressibility issue of the WCSPH framework [167], Jiang et al. [168] used volume-weighted mixture velocities u m = k c V k u k , thereby ensuring a divergence-free mixture velocity field solvable by an iterative incompressible SPH solver for the single-fluid case.To capture multiphase fluids with highly dynamic relative motions, Jiang and Lan [172] presented a dynamic mixture model which abandoned the local equilibrium condition.This method also allowed for fluid control in the multiphase environment by solving the Navier-Stokes equations for each phase flexibly.In contrast, Ren et al. [173] used the deformation gradient to construct a set of linear equations that matched the local volume change resulting from the momentum-equation-solved velocities that resulted from the continuity-equation-solved fraction changes, and solved these equations for enhanced incompressibility.

Non-trivial-diffusion Traditionally, phase-mixing effects in fluid simulation have been modeled using the diffusion equation Dc
Dt = α∇ 2 c, where α is the diffusion coefficient and c is the concentration, which assumes uniform phase velocities and movement in accordance with the aggregate motion.This approach has been employed in various works such as Im's [174] diffusive dissolved air transfer model for calculating bubble distribution in freezing ice blocks and He et al.'s [175] two-phase diffusive model for simulating diffusive appearances with varying sharpness in materials like ink and bubbles.Other researchers have advanced this field with more sophisticated models.Yang et al. [169] integrated the Cahn-Hilliard equation into multiphase simulation using an energybased model to capture complex multiphase effects such as unmixing and extraction.They computed the change of the mass fraction c m k of each phase k as where M is a degenerate mobility constant and ϕ k is the k-th phase's chemical potential relying on the derivative of a case-specified Helmholtz free energy function at the current concentration composition.This model was able to capture complex multiphase effects such as unmixing and extraction (Fig. 26).Yang et al. [153] extended this model by using a unified Helmholtz free energy form to handle both solid and liquid phases, thereby expanding the capacity of the PBF multiphase solver.Chen et al. [170] proposed a moving least square reproducing kernel particle method for better precision and stability of particle-based simulations.Using an advanced interpolation scheme, they integrated the Cahn-Hilliard equations into MPM solvers and achieved good mass conservation, stability, and sub-grid details in multiphase fluids.
Xue et al. [176] modeled anisotropic diffusive effects using non-Fourier diffusion which was integrated into a phase field formulation using an MPM solver.The resulting constitutive model is given by where q is the associated diffusion flux, q C and q F represent Cattaneo and Fourier diffusion, respectively, τ is the relaxation time with respect to the flux, α is the diffusion coefficient, and X is the quantity being diffused.G T is an adimensional parameter that represents the weight between Cattaneo-type diffusion and Fourier-type diffusion.Their method reproduced complex folding effects of poroelastic materials during wetting and also directional diffusive transportation effects.Su et al. [177] adopted the anisotropic diffusive model in [176] for temperature transport in an extended MPM phase change solver allowing the simulation of richer phenomena.They also introduced an integration scheme that provides second-order accuracy with only first-order algorithmic overhead.
Additional works studied other mixing-related phenomena in recent years.Stomakhin et al. [178] proposed an MPM approach to solving heat-induced phase change of various materials.A carefully designed projection solver allowed them to simulate nearly incompressible phase-changing materials in MPM.Hochstetter and Kolb [179] presented an SPH method to simulate evaporation and condensation of liquids.Their technique utilized particles to signify the liquid phase, while the grid primarily served as a medium for simulating the air phase and facilitating water vapor transport.This method used Fourier's law as a basis for heat transfer between grid cells and particles, thereby advancing the understanding of multiphase heat and mass transfer phenomena.Liquid Bubbles Foam Fig. 27 Schematic diagram of bubbles and foam.Gas accumulates under water and forms bubbles.As bubbles rise to the water surface, their volume increases due to a pressure decrease.When reaching the surface, bubbles may form foam.

Gas-liquid interfaces
In fluid simulations, the influence of gas is often ignored.However, numerous real-world fluid phenomena, including the formation of water droplets and bubbles, cannot be accurately represented without considering the role of gases.The phenomena formed by gas-fluid interactions are complex and diverse.In this section, we discuss gas-liquid interface phenomena grouped into three categories: free surface fluids (Sec.7.1), bubbles, foam and glugging (Sec.7.2), and spray and splashing (Sec.7.3).
In free surface fluid simulation, the emphasis is typically placed on calculating the fluid surface and accounting for surface tension, rather than explicitly modeling the presence of air or gases.We introduce some typical methods such as contact angle, surface tracking, and continuous surface force.We then discuss bubbles, foam and glugging together because of their similar characteristics, which are the effects of a small amount of gas being wrapped by the closed fluid.The bubbles we are discussing here are those formed by air gathering in water, while foam is formed by bubbles rising to the surface of the fluid (Fig. 27).Glugging occurs e.g. when a liquid is rapidly poured from a bottle with a narrow opening (Fig. 28) and is a multiphase phenomenon where bubbles are generated automatically.Finally, spray and splashing are formed by free liquids in the gas.These are usually produced by violent collisions of fluids and require more accurate simulation methods.Fig. 28 The glugging effect [180].

Free surface fluids
Physically correct or even at least plausible gas-liquid interface modeling is challenging.This is largely due to the fact that, while scalar fields such as pressure can be approximated well using particles or grids at macroscopic scales, surface tension (and similar) effects are the result of microscopic intermolecular forces (Fig. 29).This makes introducing surface tension effects into standard Lagrangian and Eulerian solvers (see Sec. 3.3) non-trivial.
Wang et al. [181] introduced the contact angle to calculate surface tension.This angle exists at the junction of solid, liquid, and gas, which indicates the hydrophilicity and hydrophobicity of solid materials (Fig. 30).They used signed distance fields to represent such surfaces and constructed a Liquid Water molecule Molecular force virtual surface below the solid one to replace the real solidfluid interface.The distance field can be modified by the virtual surface.Following this, the stable contact angle θ s can be obtained to estimate the surface tension from where γ sa , γ la , and γ ls are the interfacial tension coefficients for the solid-air, liquid-air, and liquid-solid surfaces, respectively.However, this method uses a grid to represent the internal volume of the fluid, which requires significant memory and computing time.
Air Liquid Solid γ la γ sa γ ls θ s Fig. 30 Schematic diagram of stable contact angle θs.γsa, γ la , and γ ls are respectively the interfacial tension coefficients for solid-air, liquid-air, and liquid-solid surfaces.When a drop of liquid rests on a solid surface in equilibrium, the angle between the solid-liquid interface and the gas-liquid interface is called the stable contact angle.
Water drop animation is the main focus in Zhang et al. [182].The crucial part of their Lagrangian system that allows efficient simulations of water drop motions is the reduction of volumetric fluid dynamics over the whole liquid volume to a deformable surface model.While also using the contact angle method like [181], their model focuses only on the surface and is as such more computationally efficient.
Da et al. [183] proposed a surface-only model that avoids having to deal with degrees of freedom inside liquids and (often) far away from their surface.This is the first such model for 3D liquids with the first advection-projection scheme for surface-based liquids, albeit partly limited to bodies dominated by surface tension and inertia, although still capable of modeling effects such as crown splashing.
Akinci et al. [184] employed SPH as a method to model surface tension and adhesion forces.Theirs is the first method Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 23 Fig. 31 Simulation of high surface energy liquids including shape change and two-way coupling with solids [188].
that correctly handles large surface tension (and adhesion) without the need for ghost particles or artificial pressure forces.The cohesion force is described as where γ represents the surface tension coefficient and B denotes a spline function.The method is simple to integrate with existing SPH solvers and can simulate effects such as water crown formation and rolling water droplets.By assigning each particle a value corresponding to an estimate of its surface area, leading to an implicit definition of the free surface of the fluid, Orthmann et al. [185] achieved conservative transport within and between surfaces, including correct handling of thin sheets and other singularities.This allows for effective simulations of detergents, cleansing, and coating.
Yang et al. [186] used a pairwise-force model called PF-SPH which relies on larger support radii than traditional SPH.The method improves the accuracy of the surface tension calculation by using anisotropic filtering to scale neighbouring particle interaction forces.
Energy-based methods have also been used to simulate free surface fluids.He et al. [187] modified earlier surface tension and air pressure formulations for SPH-based free surface flows, building on the diffuse interface model.They introduced a modified surface tension energy formulation E s as where V represents the volume of the liquid, κ is a coefficient associated with squared gradient energy, and l denotes the condensation field.The surface tension energy E s is directly related to the surface area of a fluid interface.Its gradient can be computed to determine the surface tension force acting on the interface.This improves the robustness of the model vs particle sparsity and in turn leads to increased stability.The model can simulate delicate surface tension effects such as water/milk crowning.Classical methods often struggle when it comes to relatively high coefficient/parameter values such as those controlling surface energy.To address this, Hyde et al. [188] developed an implicit Lagrangian formulation.This formulation specifically targets liquids with significant surface energy, such as liquid metals (Fig. 31).By treating discrete forces as gradients of the potential energy that are proportional to the surface area of the liquid, this approach enables more accurate and stable simulations.Chen et al. [189] proposed an MPM approach which generalizes [188] by improving resampling via new types of temporary 'balance' particles which achieve the perfect conservation of grid linear and angular momenta.

Bubbles, foam, and glugging
Indeed, the influence of gas on fluid simulation extends beyond just the free surface of the fluid.It encompasses the behavior of the fluid interior and involves more intricate interaction processes.Single-phase liquid simulations typically struggle to capture phenomena like bubbles, foam, and glugging effects, which necessitates the modeling of gas and liquid as two-phase flows.
Patkar et al. [190] presented a hybrid Lagrangian-Eulerian scheme for converting between small (i.e., sub-grid and under-resolved) Lagrangian bubbles and larger well-resolved bubbles modeled with an Eulerian approach based on level sets.Their framework includes a bubble seeding mechanism to realistically simulate fluid structure interaction with complex (moving) objects.Cho and Ko [191] combined the volume of fluid (VOF) with the sub-grid refinement of the level set method to simulate moving interfaces in two-phase flows.
Goldade et al. [192] developed a model for immersed bubble simulation which avoids advection and projection inside bubbles.It is based on constraint-based incompressible bubbles (with zero density) and affine fluid regions (to account for non-zero density coefficients).The simulation region is divided into a fluid region Ω f , solid region Ω s , and air region Ω a .Any enclosed and continuous region filled with air is treated as a bubble.Linear velocity constraints are imposed on each bubble via where Ω ai is the continuous region of bubble i. n f and n s are fluid's normal and solid boundary's normal, respectively.Additional special bubble simulations exist.Paddilla et al. [193] modeled bubble rings via vortex filaments of variable thickness assuming that advective inertial forces are small with respect to viscous forces.Filaments are expressed as a configuration manifold on which the equations of motion are geodesic.Langlois et al. [194] introduced a set of techniques aimed at generating sound representations for intricate twophase liquid animations.They extended the open-source Gerris solver [195] (a finite-volume-based multigrid solver) to achieve audio-visual fluid (and bubble) simulations.
Although some of the above methods can handle foams to some extent, specialized methods exist for this.Busaryev et al. [196] animated bubble interactions in liquid foams by treating (small) bubble particles as sites of Voronoi cells in a weighted diagram.Their framework handles bubble-bubble, bubble-liquid and also bubble-solid interactions, giving rise to foam simulations with bursting and coalescing.Kim et al. [197] modeled foam waves using the FLIP solver.Foam particles projected to 2D give rise to depth and acceleration maps, making the method efficient.The method provides the option to art-direct the foam effects using sketches and level-of-detail controls.Recently, Wretborn et al. [198] presented a realistic model for white-water simulation (Fig. 32).Their method enhances simulations with (tiny) bubble and foam detail by a stable coupling scheme between bubbles and water, a novel bubble emission scheme, and manifold advection for accurate foam tracking.
For the simulation of glugging, Boyd and Bridson [199] proposed the MultiFLIP method that extended the FLIP method to two-phase flows.They treated not only liquid, but also air, as incompressible phases, both modeled via particles.This (re)produces, among others, the glugging effect.
Ando et al. [200] introduced a stream function-based solver as a FLIP variant.In this approach, the stream function ψ is used to determine the divergence-free velocity field u, given by u = ∇ × ψ.Interestingly, their work shows that solvers based on stream functions are just as viable as regular pressure solves.The method is able to simulate glugging without modeling the second phase (air) explicitly.

Spray and splashing
Spray and splashing are very common phenomena in fluid scenes (Fig. 33).For scenes with intense collisions like turbulence, the final visual effect largely depends on the fidelity of the spray and splash simulation.Nielsen and Østerby [201] modeled spray as a two-way coupled two-continua with different volume fractions to achieve realistic spray motion.However, a grid-based density field cannot capture the motion of a single droplet.In contrast, Jones and Southern [202] focused on efficient physics-based droplet interaction.They introduced coalescence, separating, and fragmenting collision outcomes into a novel particle interaction model to simulate droplets.This provides a ballistic particle system for liquid droplets and spray.
Yang et al. [203] focused on spray simulation such as arising from high-speed/violent liquid streams.Similar to [190], they also used a hybrid Lagrangian-Eulerian model (with FLIP components in their case) to model mixture phenomena with high fidelity.Their efficient CUDA implementation allows modeling droplet and spray effects such as arising in waterfall and fountain simulations.
Guo et al. [204] addressed the stability challenges encountered in the two-phase lattice Boltzmann model (TP-LBM) by introducing a novel density-aware sub-grid-scale model.Their approach can uniformly simulate different gas-liquid phenomena, allowing for realistic and visually compelling representations of gas-liquid flow dynamics.
Li et al. [180] proposed a multiphase flow method to simulate complex effects such as bubbling, glugging, wetting, and splashing.A single model captures all these effects by building on the kinetic-based Lattice Boltzmann Phase-Field (LBM-PF) method.The interface motion is governed by the conservative phase-field equation where the phase field c Φ represents the percentage of the flow phase; the mobility M controls the degree of interface splitting, ξ denotes the interface width, and n inter corresponds to the interface normal.This method is highly general and versatile and can produce results comparable to industrystandard CFD simulations.
Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 25 In a similar vein, Li et al. [205] introduced a kinetic approach to multiphase fluids.Their scheme incorporates an accurate collision model and is able to robustly capture intricate and visually-appealing behaviors such as the injection of gas into the liquid (Fig. 34).
Fig. 34 The injection of gas into a glass of water [205].

Fine details enhancement
Capturing high-frequency details of fluid surfaces, such as vortices, waves, and turbulence, is key to enhancing the realism or artistry of fluid simulations.The use of higherorder advection numerical methods or finer discretization can alleviate the numerical dissipation problems inherent in fluid simulation.Yet, this creates very high memory and computation time needs.To cope with this, several methods were designed to specifically add fine details to a coarse fluid simulation in the last decade.We group these into three classes: reduced-dimensional simulation on the fluid surface only (Sec.8.1); dynamic methods to combat numerical dissipation (Sec.8.2); and data-driven methods (Sec.8.3).

Reduced-dimensional simulation on the fluid surface only
These techniques decouple the surface simulation from the volume simulation, allowing a secondary model with highresolution surface features to be added to a coarse (thus fast to compute) volume fluid model.We further split methods in this class into embedding techniques, 2D water wave simulation, and surface tracking and reconstruction, as follows.

Embedding techniques
The Closest Point Method (CPM) [206] is a numerical method for solving PDEs on surfaces.Unlike 2D surface parametrization, CPM typically uses a 3D Cartesian grid to discretize narrow spatial bands around the surface, which allows it to scale according to the complexity of the surface rather than the volume.Auer et al. [207] used CPM to simulate fluid effects on static surfaces in real-time.Auer and Westermann [208] followed up with a semi-Lagrangian CPM that alleviated some technical limitations of previous applications of CPM to deformed surfaces.Their method is unconditionally stable for surface deformation.Kim et al. [209] used CPM to explicitly perform high-resolution wave simulations on the liquid surface.They used the iWave algorithm [210] to produce more realistic water waves than the traditional wave equation, which can be expressed as where H is the fluid height, g is the gravity constant, and √ −∇ 2 is a fractional Laplacian operator.Mercier et al. [211] added a sub-grid wave model to particle-based liquid simulations to enhance such simulations with additional turbulence.Goldade et al. [212] worked to eliminate sub-grid errors in underlying surfaces and reduced artifacts in narrow bands around surfaces.Morgenroth et al. [213] used CPM to efficiently compute high-resolution 2D simulations on rough surfaces.Their method is similar to [208] but adds mass and momentum conservation and can produce interesting effects such as oil films on water surfaces and thermal convection on a hemisphere.

2D water wave simulation
To reduce computational complexity while keeping surface detail, some researchers have investigated the simulation of water waves on fluid surfaces.Water wave simulations are independent of degrees of freedom, so high-frequency visual detail can be created without increasing the simulation resolution.
In response to the inability of the shallow water equations (SWE) to capture motion details such as wakes, Pan et al. [214] used a 2D discrete vortex method to capture the wake behind a moving rigid body.Their method requires only a small number of wake particles and is fast enough for real-time applications.Yet, this method fails to handle the complex wake patterns caused by vortex stretching and tiling in 3D flow.Azencot et al. [215] used a scalar vorticity function on 2D domains to describe the vortex behaviour of fluid surfaces, greatly simplifying the analysis and simulation of fluids.Later, Azencot et al. [216] simulated the complex behaviour between multiple waves, including annihilation, recreation, splitting, and merging, by solving the EPDiff [217] on arbitrary triangle meshes using an explicit structurepreserving numerical scheme.
A key challenge in wave simulation is handling the coupling between waves and obstacles.Canabal et al. [218] generated rich water waves using a dispersion kernel as the spatially variant filter and further simulated interactions between waves and static or moving obstacles by modulating this dispersion kernel.The dispersion relation under the Airy wave theory [219] defines the propagation speed of each wave u c as where w is the wave number, g is the gravity constant, ρ and γ are the density and surface tension of the fluid, respectively, and H is the fluid height.Jeschke and Wojtan [220] simulated the movement and interaction of a large amount of waves by a wavefront tracking algorithm with multivalued function interpolation.Their method can model the dispersion, refraction, reflection and diffraction of waves well, but only handles scenes with static obstacles.Later, they introduced the concept of wave packets [221] which can handle the interaction of water waves and moving objects.They used an improved Lagrangian particle method to simulate the diffusion of water waves to add more visual detail.However, this method cannot be extended to moving 3D fluid simulations of surfaces.The same problem was addressed by Skrivan et al. [222] who decoupled the wave resolution from the simulated resolution using Lagrangian wave packets.This method significantly increases the visible detail on the fluid surface as a post-processing step.
Creating large open water animations and adding wave detail is a common requirement for a variety of interactive or offline applications.Implementing this requires carefully balancing visual quality vs computational resources.Nielsen et al. [223] proposed a wave synthesis technique based on the Fourier transform to enhance the details of wave animation.However, wave-obstacle interactions are difficult to incorporate into the spectral solver.Keeler and Bridson [224] proposed an efficient surface-only simulation of deep ocean waves and used a new indirect boundary integral equation to deal with wave-solid boundary interactions.The Method of Fundamental Solutions (MFS) was also used to generate realistic waves behind moving obstacles.Schreck et al. [225] proposed a novel discretization for MFS using wavelets and achieved naturally-looking wave interactions with complex boundaries (Fig. 35).Their method obtained impressive results on a large-scale ocean scene.Jeschke et al. have successively developed two interactive systems for the simulation of large ocean scenes that can handle detailed wave details [226] and coupled interactions with complex terrain [227], respectively.Recently, Schreck and Wojtan [228] proposed a coupled method of 3D liquid simulation and 2D wave propagation to simulate infinitely large bodies of water and fine surface wave detail.An empirically driven error compensation method is Fig. 35 Water waves can accurately interact with complex boundaries [225].also used to remove coupling errors from the simulation to achieve a seamless transition between 2D and 3D.

Surface tracking and reconstruction
Tracking and reconstructing fluid surfaces is important for generating realistic animation effects.This is difficult to achieve due to the complex shapes and frequent topological changes of fluids.
For Eulerian fluid simulations, robust handling of surface triangle mesh splitting and merging can remove visual artifacts to preserve important surface features.Bojsen-Hansen et al. [229] proposed a method for tracking the topological evolution of surfaces that can solve the wave equations on lower resolution fluid surfaces to synthesize high-frequency details.Later, Bojsen-Hansen and Wojtan [230] presented a novel physics-based surface fairing method that solved the physical and topological artifacts arising from coupling high-resolution surface trackers with low-resolution fluid simulations by introducing an error metric and surface correction force.Edwards and Bridson [231] presented a new approach to adaptive fluid simulation.They track explicit triangulated mesh surfaces and use the p-adaptive Discontinuous Galerkin (DG) method [232] within detailed cut cells near the surface.The treatment of dynamics is guaranteed to be physically consistent while reducing computational costs by using coarse-grid fluid simulations.Chentanez et al. [233] devised a grid-free surface tracking method that deals with topological changes by removing overlapping triangles and performing effective triangulation of the generated holes.This method can be used in both mesh-based and particle-based simulations.
Inaccurate detection of free surface particles in particlebased fluid simulations can lead to unrealistic artifacts.Also, irregularly distributed particles can make the reconstructed surface bumpy.To address the fact that particles do not keep connectivity information, Yu et al. [234] proposed to periodically project surface meshes to match implicit surfaces Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 27 defined by fluid particles.This method allows the simulation of high resolution surface waves without the limitation of particle resolution.Later, Yu and Turk [235] proposed a method for reconstructing surfaces in particle-based fluid simulations.They utilized a stretched anisotropic smooth kernel to represent each simulated particle, resulting in a greatly improved surface quality (Fig. 36).Sandim et al. [236] proposed a fast free surface detection method that only requires the positions of particles to identify surface particles without using kernel functions or normal vectors.This method is applicable to cases with non-uniform particle distributions and complex free surface deformations.Dagenais et al. [237] used an explicit mesh projection method based on signed distance fields to preserve surface detail and introduced a new topology matching operation to maintain consistency between explicit surface and particle behaviour.

Dynamical methods for reducing numerical dissipation
The advective-projection method leads to numerical dissipation which results in kinetic energy decay and suppression of motion such as vortices and turbulence.Bulk enhancement methods aim to improve the whole fluid volume rather than only its free surface.We further group such methods that aim to improve system energy conservation and detail preservation by reducing numerical dissipation in vorticity confinement, vortex-based methods, and various variants of dynamics solvers, as described next.
Vorticity confinement methods.Such methods are based on the principle of vorticity conservation, which adds a vorticity control term to restrain the diffusion of vortices, thus simulating fluid dynamics problems such as turbulence and vortex streets without dissipation.Fedkiw et al. [238] first applied the vorticity confinement method to smoke simulations by adding an additional force field to maintain the airflow vorticity.Second-order vorticity confinement (VC2) further improved this method by ensuring momentum conservation.The confinement term of VC2, F conf , is given by where ∆x is the grid size, α and ζ are the positive and negative diffusion coefficients, respectively, ω is the angular velocity, and m is the harmonic mean of the local vorticity stencil.Lentine et al. [239] improved the vorticity confinement model by allocating global momentum to ensure momentum conservation.Jang et al. [240] applied the multi-level vorticity confinement method to simulate water turbulence in order to capture large/small-scale vortices and complex flow details.He and Lau [241] proposed adaptive adjustment of the positive diffusion term to balance constraints, which broadened the stability conditions of the VC2 method and made it capable of generating highly turbulent flows.However, the vorticity constraint method can only enhance existing vortices or turbulence and may not be effective for other types of flows, such as laminar flow.Moreover, computational costs can increase significantly due to additional constraints, especially for large-scale complex fluid scenarios.

Vortex-based methods
The potential vorticity field can be effectively modeled using vortex-based methods.These methods simulate the vorticity of the velocity field rather than the velocity field itself, so this automatically guarantees a divergence-free velocity field and removes numerical dissipation.Most such methods are Lagrangian and model the vorticity form ω v of the Navier-Stokes equation as This formulation represents the vorticity distribution as a superposition of singularities.Zhang et al. [242] proposed a new scheme named IVOCK, which aims to solve the errors and energy loss caused by the X.K. Wang et al.
self-advection step through compensating the vorticity error.However, this method is only applicable to fluid simulation on uniform grids.Liu et al. [243] extended this idea to particlebased turbulent detail simulation (Fig. 37).Recently, Xiong et al. [244] proposed a vortex segment method to simulate flows with strong anisotropic vortical features.Compared with existing Lagrangian vortex particle methods, this method can more vividly model complex phenomena such as the splitting and reconnection of two vortex tubes or vortex shedding near a solid boundary.
The main challenge of vortex methods is the handling of fluid-solid coupling and creating vortices at this coupling boundary.For this, Golas et al. [245] proposed a combination of Eulerian simulation and vortex singularity bases.By using Lagrangian vortex elements inside the fluid and enforcing boundary conditions in the Eulerian mesh, robust interaction of free surfaces and non-rigid obstacles can be achieved.Zhang et al. [246] used a FLIP approach to solve the Navier-Stokes equation using viscous particle strength exchange, handling the momentum transformation at the solid boundary effectively.Liao et al. [247] proposed a new wall-bounded turbulent smoke simulation method, which introduced particle-particle interactions to traditional vortex filament mesh calculations to accurately capture the vortices and thin turbulence generated by smoke-obstacle interactions.

Various variants of dynamics solvers
Additional methods improve on current advection-projection solvers or extend classical dynamical methods to achieve detail enhancement.In the advection-projection step, detail and energy preservation are greatly improved by the introduction of detail capture and shape correction techniques [248], the use of energypreserving reflection operators [249], and feature mapping with convectors [250].
Yang et al. [251] first introduced the Clebsch wave function as a system scaling variable to evolve Eulerian flow fields, which significantly improves the ability to generate and sustain vorticity in simulations of various gases and liquids.Later, in order to solve the numerical instability problem of Clebsch's method near dynamic interfaces, Xiong et al. [252] proposed a new wave function correction scheme and extrapolation algorithm, which achieved detailed simulation of various vortex structures on free surfaces.Recently, Feng et al. [253] proposed a numerical method for solving the Navier-Stokes equations based on the pulse gauge transformation which can generate rich vortical details by treating the fluid pulse as an auxiliary variable.Liu et al. [254] and Li et al. [255] created two turbulence simulation methods using an adaptive multi-relaxation scheme and statistical mechanics, respectively.Later, Lyu et al. [256] further improved the boundary treatment of dynamic solids, enabling the simulation of fluid-solid coupling between thin structures and turbulent fluids (Fig. 38).
The micropolar fluid model is an extension of the classical Navier-Stokes equation which takes into account not only the linear but also the angular velocities of the fluid particles, thus enhancing the eddy and turbulence details of the fluid.Bender et al. [257] first used the micropolar fluid model to simulate the turbulence phenomenon of non-viscous fluids.Subsequently, they post-processed the foam phenomenon on this basis to significantly improve the realism of the visual effect.

Data-driven methods for detail enhancement
Texture synthesis Geared toward the production of artistic effects based on fluid elements, texture synthesis is the technique of choice for adding detail to surfaces.As a postprocessing method, texture synthesis achieves detailed surface features through patch or style transfer based on deep learning.
Patch-based texture synthesis maps image textures or simulated features to the target flow field, improving the appearance of the source simulation by matching the target dataset features.Jamriška et al. [258] used per-pixel best-fit search to achieve rich visual effects through 2D input image appearance transfer.However, this method is limited to image space synthesis and is difficult to extend to 3D fluid surfaces.Gagnon et al. [259] proposed a temporally coherent patch-based texture synthesis method to handle scenes with significant deformation and topological changes.This approach aims to maintain a Poisson disk distribution of patches on a free surface to find the optimal parameter values and locations of time-varying patches.For the ghosting problem of overlapping patches in [259], Gagnon et al. [260] followed up with a solution Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 29 scheme based on patch erosion.They used feature-aware erosion to remove patch distortion textures to ensure realism of the mapping.
In contrast to surface texture synthesis, deep neural networks perform various stylization tasks on volumetric data by extracting 2D image features.Sato et al. [261] proposed a style transfer method that migrates high resolution turbulent details to low resolution flow fields, which speeds up the fine surface detail simulation by nearly 30 times.In addition, they used an optimized texture synthesis method to solve the problem of discontinuity at the patch boundary.Kim et al. [262] first proposed a transport-based neural style transfer algorithm that enables automatic conversion of the semantic structure of 2D images into 3D smoke simulations.The method achieves complex artistic effects by optimizing the transport of smoke to the desired stylized velocity field at each time step.However, this method cannot handle the transfer of color information.Therefore, they further redefined the method in a Lagrangian setting to ensure better temporal consistency and support for color stylization [263].Unlike stylization methods which focus on fluid simulation shapes, Guo et al. [264] proposed a Stylizing Prediction Network (SKPN) aimed at stylizing physical color appearances.The method can easily generate the user's desired color appearance without complex parameter tuning.Ten years ago, Zhang and Ma [266] proposed a spatiotemporal extrapolation technique that enables high-resolution flow features on coarse grids.Chu and Thuerey [267] enhanced the turbulence detail of the smoke simulation on the coarse grid by using local patch descriptors.Um et al. [268] proposed a deep neural network to capture smallscale splashed droplet details from low-resolution liquid simulations.Xie et al. [265] used a conditional generative adversarial network with a temporal discriminator to directly generate advected quantities with highly detailed and temporally coherent features for smoke simulation (Fig. 39).CNNs have been used to create matching models to correct the shape of low-resolution smoke simulations [269] and estimate physical parameters to guide the reconstruction of high-resolution velocity fields [270].Bai et al. [271] used a dictionary-based neural network for fluid upsampling.Yet, the choice of the training set was somewhat limited and the spatial and temporal consistency of the results could not be guaranteed.They next significantly improved the prediction quality of the network by adding filtering to the training process [272].With the development of Deep Neural Networks (DNNs) in recent years, Roy et al. [273] proposed a DNN-based method for improving the resolution of coarse particle liquid simulations.

Fluid control
The wealth of methods for fluid simulation surveyed so far leads next to a key question: How to control such simulations?Even though a method can technically produce highly accurate results, its ultimate target is to enable its users to steer the method toward the desired results.We next group and discuss fluid control methods in three classes using different control perspectives: scenario editing (Sec.9.1), artificial effects (Sec.9.2), and media-directed formation (Sec.9.3).

Scenario Editing
Scenario editing steers the fluid by generating new simulation scenarios based on existing simulation results without losing the characteristics of the original simulations (Fig. 40).On the positive side, such control is technically the closest to how a simulation works, so it can steer the fluid most directly.On the negative side, this control requires the end users' advanced skills and an understanding of the underlying simulation and overall fluid flow technicalities.We further divide scenario editing into three sub-types based on the implementation approach: target guided editing, adjustable editing, and camera-based editing.
Target guided editing.Such methods modify an existing fluid simulation to match a given target, e.g., a higher resolution.Gregson et al. [274] connected low-resolution smoke capture with its velocity field.They treated the pressure projection as a proximal operator and tracked the fluid by estimating its velocity.Through advection, their method obtained a high-resolution re-simulated smoke.Forootaninia and Narain [275] successfully guided high-resolution smoke flow by replacing its low-frequency component with a given guiding field.Generally, the guiding task is seen as an optimization problem that minimizes errors.This optimization problem was solved efficiently by using a fast primal-dual method [276].To achieve a more desirable artistic effect, it is essential to guide smoke animation in such a way that it aligns with one or multiple target density keyframes provided by the artist.In order to address this control problem, Pan and Manocha [277] formulated it as a space-time optimization.They employed an Alternating Direction Method of Multipliers (ADMM) optimizer [278] to derive a dense sequence that forces the smoke to attain the desired target shape.Tang et al. [279] proposed an advanced method that effectively addresses both the issues of unconstrained optimization and high-dimensionality of the parameter space simultaneously.
Adjustable editing.These methods aim to control, edit, or resize fluid simulations during implementation.Raveendran et al. [280] focused on control with an emphasis on the liquid surface and proposed a method for creating high-quality fluid animations that provides the animator with multiple control levels.Further on, Raveendran et al. [92] proposed a smooth blending method to interpolate between two or more existing pre-computed liquid simulations.With their method, one can generate hundreds of different plausible results at interactive rates with potential applications in games and virtual reality.Sato et al. [281] proposed a smoke blending method to help animators synthesize the desired fluid animation.Velocities at the boundaries are interpolated by minimizing an energy function.This approach significantly reduces computational costs by reusing existing flow data instead of creating realistic fluid animations by numerical simulation.Fluid carving is another way for editing fluid simulations.By utilizing seam carving, efficient and effective resizing of 4D fluid simulation data can be achieved [282].Flynn et al. [283] proposed a lattice-guided seam computation method which broke through the limitation of rectangular boundary and reduced calculation time.
Camera-based editing.In the context of large-scale scene simulation, due to computational cost considerations, it is often necessary to perform coarse-grained simulations of the entire world, and subsequently integrate finer-grained details into the scene.To integrate two simulated fluid scenes seamlessly, Bojsen-Hansen and Wojtan [284] presented a fluid modification approach with 'non-reflecting' boundary handling.They extended the simple Perfectly-Matched Layers (PMLs) method to handle coupling inflow/outflow boundaries with varying spatial and temporal conditions.The boundary can be modified easily during the simulation and it handles the multi-resolution combination.Stomanuykhin and Selle [285] introduced the Flow-Animated Boundary (FAB) method, in which the boundary can have a custom shape and vary over time, with materials outside the boundary dynamically removed using volume flux.

Artificial effects
One often needs to artificially edit and control fluids in order to achieve specific artistic effects.Different from scenario editing, artificial effects add new characteristics to a simulation by artificially guiding the formation of fluid shapes or movements during the process (Fig. 41).Due to the complex motion of fluids, keyframe animation, which involves controlling the flow of fluids to match keyframes, is a commonly used method for fluid control to reduce unrealistic effects in simulations.However, manually designed frames often lack volume preservation and exhibit excessive smoothness, resulting in the loss of simulation details.
Pan et al. [286] proposed a local control method instead of globally manipulating the entire fluid, allowing users to edit and control fluid shapes in specific regions using a brush-like tool.However, controlling the simulation process between keyframes is challenging.To address this issue, Lu et al. [287] drew inspiration from skeletal animation techniques.They introduced a method that controls fluid motion by manipulating Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 31 a point cloud with rigid body motion and incompressible deformations, and subsequently performing skinning operations on the point cloud.Similar to [287], Yan et al. [288] applied conditional generative adversarial networks to generate fluid splashes based on simple user-defined sketch input (Fig. 42).Control particles with attractive forces provide an efficient way to reproduce complex motion by using pre-computed templates [289].In this work, a set of shape-constraint particles are seeded and a repulsion force field is computed to control the shape of the final result.

Media-directed formation
Pure physics simulations often suffer from significant computational time requirements, making them impractical for real-world production applications.Media-directed formation aims to set up simulation scenarios based on real-life videos or images of the fluid, reproducing real-world scenes (Fig. 43).These methods estimate fluid properties such as volume, density, motion, and style from visual data.Newtonian fluids Following Newton's viscosity law, the viscosity of a Newtonian fluid (incompressible and isotropic) can be expressed by a material parameter µ called dynamic viscosity.Using this, the inner viscous stress tensor T vis is computed as where E is a symmetric strain rate tensor that describes the shear strain rate.This equation indicates that the viscous stresses of Newtonian fluids are at every point linearly correlated to the local strain rate.Using the spatial derivatives of the velocity field, the strain rate tensor E can be defined as Substituting Eqn.(31) and Eqn. ( 32) into the viscosity force field formulation f vis = ∇ • T vis and adding the incompressible condition ∇ • u = 0, the viscosity force field f vis can be computed as This gives the viscosity term in the momentum equation, corresponding to the viscosity term in the form of kinematic viscosity in Eqn.(6).
Discretizing the viscosity term a challenging problem for SPH based methods.Takahashi et al. [294] proposed an implicit Euler integration to solve the viscosity term separately, using two SPH first derivatives to discretize the strain tensor and the divergence of the stress tensor respectively.This method supports a larger range of viscosity and time step values, but a second-ring neighbor computation is required which impacts efficiency.Peer and Teschner [295] decompose the velocity gradient into three tensors -spin rate, expansion rate, and shear rate.A user-defined viscosity parameter modifies the shear component which describes the dissipation because of viscosity.This leads to a target velocity gradient that is used to reconstruct the final velocity field with a first-order Taylor approximation.Since the velocity gradient field is decomposed, the linear system of the velocity field can be solved separately.However, shear viscosity does physically affect the rotation component in the velocity gradient because of the tangential component in rotation.Peer et al. [296] extended this method using a vorticity diffusion scheme.The spin rate tensor in the target velocity gradient is modified by solving a vorticity diffusion process, which uses the viscosity parameter in [295], so that vorticity damping is introduced to achieve more realistic effects.Instead of using the strain rate, Weiler et al. [297] introduced an implicit viscosity solver based on the Laplacian of the velocity field in Eqn.(33).With a symmetric form of the approximation discretization of the viscosity term [298], an implicit linear system for a new velocity field can be obtained.Fig. 45 Coiling and bulking effects of highly viscous fluids [299].
Fig. 46 Change the viscosity to achieve the effect of different viscous fluids such as cream, jam and chocolate sauce [300].
The above-mentioned solvers separate the solving of pressure and viscosity, which reduces accuracy and cannot generate free surface details.Larinov et al. [299] introduced a unified pressure-viscosity solver based on implicit variational unsteady Stokes flow problems for grid-based methods, where the inertial force is considered to improve accuracy and achieve a wider viscosity range (Fig. 45).Combining the semi-implicit equation of correlation pressure (SIMPLE) method with SPH, Liu et al. [300] used the result of the pressure Poisson equation to improve the pressure in the viscosity solver in an iterative process, which converges to a globally optimal solution (Fig. 46).
Even with a stable solver, mimicking the viscosity coefficient of the target fluid is essential in order to realistically simulate highly viscous fluids.Takahashi and Lin [99] proposed a framework to find the required viscosity parameter from real videos of highly viscous fluids by minimizing an objective function that evaluates the difference between the silhouettes extracted from video frames and those obtained from the simulation.The Carreau-Yasuda model is a well-known method to simulate non-Newtonian fluids by defining a shear-rate-related viscosity µ as

Non-Newtonian fluids
where µ 0 is the zero-shear viscosity, µ ∞ is the infinite viscosity, ε is the shear rate, R is the relaxation time that scales the shear rate, n is a power law index, and a models the transition smoothness between the Newtonian plateau and the power law regime.When n = 1, this model becomes a Newtonian fluid with dynamic viscosity µ 0 .For shear-thinning fluids (n < 1), as the strain rate increases, viscosity will vary from µ 0 to µ ∞ .For shear-thickening fluids (n > 1), the viscosity increases as the shear rate increases.Bingham plastic fluids are another typical non-Newtonian fluid that behaves as a rigid body at low stress but flows as a Newtonian viscous fluid once the yield stress is exceeded.For many viscoplastic fluids, the stress curve of the flowing part is nonlinear when the shear rate exceeds a critical value.To capture shear thickening/thinning, µ is modeled by the Herschel-Bulkley model given by where k is the consistency coefficient, σ 0 y the yield stress, and ε0 the critical shear rate.Similar to Eqn. (34), when n = 1, this model describes the ideal Bingham plastic: n > 1 models shear thickening Bingham plastic, and n < 1 models the shear thinning Bingham plastic fluid.
Recent work used the above two models to simulate non-Newtonian fluids.Zhu et al. [293] simulated various codimensional features of different non-Newtonian fluids, e.g, shear thinning and thickening for Bingham plastics, and elastoplastics.The Carreau-Yasuda model for non-Newtonian fluids was used on a multi level-set model; semi-implicit methods were used for elasticity and variable viscosity.On the rims of thin fluid sheets, viscosity got an improved treatment to yield twisting motion.Yue et al. [301] used the non-Newtonian Herschel-Bulkley model to simulate dense foams composed of microscopic bubbles using MPM.They also proposed a particle resampling method for MPM and a tearing model to simulate tearing/connectivity recovery by explicitly han-dling the weakening regions detected from space.Mixtures of non-Newtonian fluids were studied by Nagasawa et al. [302].Using the Herschel-Bulkley model, a nonlinear blending model that satisfies the five blending laws [303], along with mass conservation, was proposed to capture non-Newtonian fluids' mixture behavior.For viscoelastoplastic materials, a constraint-based method [304] extended position-based dynamics to simulate elastoplastic and highly viscous fluids by recasting a constitutive model of viscoelasticity which defined governing equations for a conforming tensor.Fig. 47 Schematic representation of ferrofluids in interaction with a magnetic field.Sub-figure (a) depicts the superposition of a uniform vertical magnetic field and the magnetic field originating from an ellipsoid magnetization.There is a pronounced discontinuity in the magnetic field on the surface, most significantly at the extremities where there is a sharp escalation in the field strength.Sub-figure (b) illustrates a surface perturbation in the magnetized ferrofluid, resulting in a localized concentration of magnetic induction lines.The ferrofluid is drawn towards this bump due to a heightened field strength, which consequently enhances the gradient, amplifying the perturbation and leading to the formation of a spike.

Ferrofluids
The dynamic interactions of ferrofluids -liquid media responsive to external magnetic fields -have emerged as an area of considerable interest within the computer graphics research community.These magnetically active fluids, originally conceived by NASA to facilitate fuel transfer in spacecraft under microgravity conditions, derive their magnetic properties from the incorporation of nanoscale magnetic particles.Upon exposure to an external magnetic field, these dispersed particles within the ferrofluid polarize, thereby generating a distinct internal magnetic field.This induced magnetic field, working synergistically with the external one, is pivotal to the magnetization process of the ferrofluid, as shown in Fig. 47.
The computer graphics simulation of these captivating fluids, however, had not received significant attention, until the very recent pioneering work of Michels and his team, which directed attention towards this specialized field [305].We reference several key insights from their contributions in the subsequent discussion.The interrelationships of these spatial magnetic fields, both internally produced and externally imposed, conform to the well-established principles outlined in Maxwell's equations where B is the magnetic flux density describing the spatial magnetic field; H is the total magnetic field intensity; J is the free current density; and D is the electric displacement field.The ferrofluids further discussed in this section comply with the following model where µ E 0 is the vacuum permeability, M E is the magnetization field describing the density of the magnetic moment induced by the total magnetic field, H ext is the external magnetic field, and H int is the internal magnetic field generated by the ferrofluid.The magnetization field is a function of the known external magnetic field and the internal magnetic field, and the current internal magnetic field can be obtained by solving Poisson's equation with the assumptions of ∇ × H ext = J and ∂D ∂t = 0 in Eqn.(36), since the free current density J is not influenced by H int and the electric displacement produced by the flow of ferrofluid is not strong enough to influence the system.
Under the effect of magnetic force, fluid particles gather on tiny bumps near the surface where the synthesized magnetic field is stronger and pull the fluid to form spikes with attractive visuals.The spike shapes are also influenced by gravity and surface tension.
The past few years have witnessed the proposal of various methods for simulating ferrofluids within the field of computer graphics.Huang et al. [305] presented the smooth magnets method which utilizes Lagrangian particles embedded with magnetic nanoparticles to discretize fluids, as shown in Fig. 48.This method was the first one in the area of computer graphics to address the first-principle-based macroscopic simulation of ferrofluids.Their proposed magnetization model, along with the magnetic field's Poisson Fig. 48 A ferrofluid climbs up a steel helix under a strong external upward magnetic field to create surface spikes [305].equation, can be discretized using a smooth kernel function akin to that used in SPH.
Approaching from a Lagrangian perspective, the volumetric Kelvin force model is utilized to characterize the magnetic force interactions between particles.Standard Smoothed Particle Hydrodynamics (SPH) computations employ particles, enabling the enforcement of incompressibility and surface tension within ferrofluids.However, the implementation of a Kelvin force model engenders an unanticipated outwarddirected force on the surface, prompting particles to exhibit levitation near this surface.To address this phenomenon, Shao et al. [306] introduced a modification, replacing the Kelvin force model with a current loop model, thereby creating an inward force.This alteration permitted the integration of the magnetic model into Implicit SPH models.Consequently, it enhanced system stability and facilitated the use of larger time-step increments.
From an Eulerian perspective, Ni et al. [307] presented a level-set method for simulating various magnetic bodies including ferrofluids.The interplay between the magnetic field and the mechanical system was addressed as an interfacial issue, and a weighted average of the internal and external magnetic fields was calculated to manage discontinuities.The resulting magnetic force, coupled with surface tension, was integrated into the Navier-Stokes equations to direct the dynamics of the ferrofluids.
Advancing the research further, Huang and Michels [308] introduced the concept of surface-only ferrofluids, a singular study to date, tested successfully against real ferrofluids.Unlike previous approaches that incorporated magnetic force as an additional term in the momentum equation, their method infuses the discontinuity of magnetic pressure into the Dirichlet boundary condition within the pressure-projection of the Galerkin Boundary Element Method (BEM) based surfaceonly liquid solver [183], specifically at the fluid-air interface.This pioneering surface scheme enhances the Helmholtz decomposition step in surface-only fluid solvers via a more precise analytic integration process.Additionally, it aug-Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 35 ments the accuracy of the pressure projection step within surface-only fluid solvers through a Galerkin BEM.
From the hybrid discretization perspective, Sun et al. [309] utilized the MPM structure to further simulate nonlinear magnetic substances in pursuit of a more general magnetic description.They used the physically-realistic Langevin's nonlinear magnetization model so that the magnetic force between magnetic micro elements is bounded without additional numerical approximation.Following the concept of MPM, this method uses particles to carry microscopic magnetic quantities, and solves Poisson's equation of the magnetic fields and the Kelvin force on Cartesian grids.Without integrating the surface tension, this method cannot form stable spikes when simulating ferrofluids.However, thanks to the versatility of MPM, it can achieve a unified simulation and coupling of different magnetic materials.

Thin-film
Thin films and bubbles are fascinating phenomena that gained special attention.A common example is a soap bubble floating in air.Bubbles produced from pure water are usually few, small, and disappear quickly due to gravity, pressure, and strong surface tension.To produce more, larger, and longerlasting bubbles, surfactants are added to the water, e.g., fatty acids common in soaps.With surfactants interspersed among water molecules, surface tension is reduced so that larger bubbles appear.The tensile deformation of the film will recover due to the difference in surface tension working like elasticity, the so-called Marangoni effect (Fig .49).While the Marangoni effect models the resilience given by inconsistent surface tension, the Young-Laplace equation describes the capillary pressure difference ∆p caused by surface tension between air and the fluid as where γ is the surface tension and n is the surface outwardpointing normal.Eqn.(38) relates the pressure difference to film shape.Batty et al. [310] developed discrete viscous sheets by building on (Lagrangian) elastic thin shells.This reduceddimensional technique describes the sheets using triangular meshes with local thickness and used the area-based surface tension derived from the mid-surface of the (thin) shell.Wang et al. [311] enhanced this to capture the surface tension flow using moving-least-squares particles.Mixed Lagrangian-Eulerian approach models not only volumetric phenomena (3D) but also those arising from thin shells (2D), filaments (1D), and even individual points (0D).Surface tension (and other forces) are handled in a unified way across all (co)dimensions, including codimensional transitions.This enables simulating complex scenarios requiring careful (surface) tension treatment, such as two water jets colliding and forming a thin sheet.Similarly to [311], the codimensional surface tension flow of Zhu et al. [312] also relies on simplicial complexes and transitions between elements of different (co)dimensions, covering thin fluid sheets, filaments, and surface tension effects.Wang et al. [313] extended [311] to yield a thin-film SPH fluid model (Fig. 50).Films are modeled as codimension one (surface) particles with local thickness estimates.This particle setup interacts with a Lagrangian flow simulation using a thin film shape description.This physically couples aggressive surface deformations and strong tangential flows.A process of transformation from codimension one to codimension zero is used to simulate rupture.
Focusing on viscous thin films, Vantzos et al. [314] proposed a numerical scheme to simulate the thin film equation (modeled as a height function) on a planar domain, including gravity and other forces.They add a novel quadratic term to the governing equation to stabilize flow while keeping visual fidelity.Their scheme is fully local so it allows an efficient GPU implementation that leads to real-time simulations, such as of honey flowing through honeycombs.
Da et al. [315] studied soap films and foams whose dynamics are captured by a Lagrangian vortex sheet model with an emphasis on circulation.Surfaces are represented using multi-material triangular meshes supporting topological changes, and their tension forces lead to a circulation update rule based on mean curvature.Fig. 51 Cluster of bubbles [316].After the top-left bubble bursts at frame 21, the geometry of the remaining bubbles gradually transitions to the next equilibrium state, following Plateau's laws.
By Plateau's laws, a steady-state film consists of constant mean curvature parts and minimal surfaces (vanishing mean curvature).Ishida et al. [316] used this to model evolving foams via hyperbolic geometric flow, a type of mean curvature flow (Fig. 51).
Most existing soap film models [315,316] assumed that a film is infinitesimally thin and has no influence on its evolution.Ishida et al. [317] extended such methods to use film thickness, modeled on non-manifold meshes, as a reduced degree of freedom in the Navier-Stokes equations, and derive the motion equations.This provides an incompressible fluid solver for 2.5D films.
In addition to dynamic effects, bubbles also produce rich color effects due to the light interference caused by uneven film thickness.Besides soap bubbles, other fluid dynamics simulations are also possible on the thin film for additional visual effects.
Two further methods focus on spherical bubbles and fluid around them.Hill and Henderson [318] efficiently simulated fluids on a sphere surface.They handle poles/singularities of spherical coordinates -which would otherwise render the motion equations complex if used naively on the sphere.Their method also enables vector and scalar controls for art-directed spherical flows.Huang et al. [319] focused on chemicalmechanical simulation of soap film flows on spherical bubbles using lubrication theory.Considering the Marangoni effect and the capillary pressure difference, the stress condition at the film surface is where T c is the Cauchy's stress tensor, n is the outward unit normal vector at the surface, p a is the air pressure, and ∇ s is the 2D gradient operator.The surface tension γ is defined by a linear model where γ 0 is the surface tension of pure water, Γ is the surfactant concentration, and γ r is a constant that describes the Marangoni elasticity of the film.The surfactant concentration Γ is advected by an advection-diffusion equation.Huang et al.'s advection scheme used spherical coordinates; using local frames removed some artifacts of [318].This last method also proposed a physically accurate shader for real-time rendering under environmental lighting.Recently, Deng et al. [320] introduced MELP (Moving Eulerian-Lagrange Particles), a novel mesh-free method for incompressible fluids on moving foams and thin films.Their approach, including multi-MELP for interfacial flow, is able to model both large-scale surface deformations as well as detailed flows.

Conclusion
Physics-based fluid simulation has been successful in games, film, and animation.Current advances enable high levels of control in production and have shown increasing acceptance and potential in virtual and augmented reality and other real-time graphics-intensive applications.Recent advances in physics-based fluid simulation methods rely on a complex mix of computational efficiency, realism, controllability, and ability to simulate diverse scenarios.This survey presented an in-depth overview of these methods over the last decade.We discussed the different goals in this field, techniques proposed to address these goals, and challenges of these techniques.
Our survey found seven major themes present in around 300 fluid simulation papers in the computer graphics community in the last decade -advanced computational approaches, interaction with materials, multiphase simulations, gas-liquid interfaces, enhancing fine details, simulation control, and special fluids.These themes structure our survey to outline important developments in this period and community.We Advanced computational approaches.Adaptive solutions and GPU parallelization are ubiquitous in computer graphics, but there is much room for developing such new methods for physics-based fluids.For instance, adaptivity usually makes an implementation complex and hard to be applied on GPU hardware.Spatial and temporal resolution are often related aspects.In space discretization, interactions between grids or particles at different scales may cause instability and fidelity loss, which can also create energy diffusion.Hence, in addition to improving computational efficiency and reducing overhead, proposing methods to reduce such unwanted effects on a wide range of resolution scales is important for future research.Using neural networks to learn fluid dynamic behavior will be a hot research topic in the future, as it has strong prospects for real-time simulation and industrial control.However, only summarizing physical laws from a large number of training data lacks underlying logical support.As such, more attention is likely to be paid next to how to inject physical prior knowledge into deep learning models.
Fluid coupling with multi-materials.The key challenge for such simulations is to keep accuracy, stability, and efficiency when coupling multiple materials in one scenario at one time.Using different solvers vs materials limited the diversity of fluid animation in the past.Current research moves from merging multiple solvers to developing monolithic ones.Monolithic solvers can simulate different materials and their interaction in a single framework which can eliminate stability issues.Yet, such solvers currently demand high computational resources.An open challenge is thus mixing hardware and algorithm designs to better support such solvers.
Multiphase liquids.Future work can explore many interesting aspects.One challenge is that current methods for incompressible fluid simulation do not handle high-density ratios well.Linear systems become ill-conditioned under highdensity contrasts so that Jacobi-like solvers fail to converge.Currently, many parameters of mixing fluids are manually adjusted.How to control parameters more intuitively to get desired visual results is worth further study.Modeling temperature, chemical reactions, elasticity blending, and optical blending are equally important open aspects in this area.
Liquid-gas interaction.Challenges for liquid-gas interaction include simulating gas-liquid phase transition and modeling its effect on surface tension, supporting non-manifold thin film structure, handling the transition between different codimensionalities, producing realistic surface colors for bubbles, reducing the simulation complexity while preserving accuracy for liquid-air coupling, and adding fine splash details that are not limited by particle size.The topic of achieving richer gas-liquid interaction phenomena while ensuring stability and efficiency is a subject of ongoing and future long-term research.
Fine details enhancement.Many existing detail enhancement techniques achieve detailed fluid surfaces without using high-resolution discretization.Realistic and fine-grained appearance representations are the primary pursuit in this direction.Energy conservation and detail preservation in accordance with physics laws are also important goals.How to improve the efficiency and scalability of detail enhancement is an open challenge for the end aim of providing real-time, interactive, and generalizable tools for artists.Future work will likely use deep learning techniques with innovative explorations in style transfer, high-resolution reconstruction, and detail generation.
Fluid control.Issues still exist in fluid control -high memory and computation costs, sensitive parameters, and manual feature labeling.A key difficulty of control methods is to achieve precise control.Achieving precise control, along with high accuracy and efficient computation, is -as for the other topics studied in this survey -one of the grand challenges of fluid simulation research.Special fluids.The challenge of simulating fluids with special forms lies in unifying governing equations and other physical characteristics in a consolidated system, which means integrating different solvers together.The targeted formulation of a monolithic-style solver is required to perform highperformance simulation results.
We hope this survey helps to introduce the theoretical concepts underpinning physics-based fluid simulation and their practical implementation to serve as a guide for researchers and practitioners as well as facilitate future works to exploit on the basis of recent developments.

A Practical resources
For beginners, how to quickly build their own fluid models and achieve convincing rendering results may be of the most concern and interest.In this appendix, we recommend several popular libraries, frameworks, as well as other software tools to help beginners get their foot in the door and get excited about learning.In addition, these resources are also useful to benchmark third-party fluid simulations.

A.1 Focus on simulation:
OpenMaelstrom is an open-source library for the simulation and rendering of fluids based on the SPH method.It provides many pressure solvers (IISPH, DFSPH, and Interlinked SPH for strong fluid-rigid coupling) and boundary handling methods.Moreover, it features spatial adaptation and full GPU support.MantaFlow includes a wide range of Navier-Stokes solver variants, and its parallelized C++ solver core, Python scene definition interface, and plugin system allow for quickly prototyping and testing new algorithms.An advantage of MantaFlow is that it can be easily integrated into Blender for end-to-end fluid simulation and rendering.In addition, MantaFlow can also be coupled with TensorFlow which helps the development of traditional physical simulation combined with deep learning techniques.PhysIKA (Physics-based Interactive Kinematics Architecture) is a node-based architecture targeted at real-time simulation of versatile physical materials.Currently, it supports simulating physical phenomena ranging from fluids, elastic objects and fractures, etc. PhysIKA is highly modularized and can also help the research community develop novel algorithms.PositionBasedDynamics supports the physically-based simulation of mechanical effects for elastic rods, deformable solids, rigid bodies, and fluids.SPlisHSPlasH, by the same author as PositionBasedDynamics, simulates complex fluid effects based on SPH methods.It includes several SPH solvers (including WCSPH, PCISPH, IISPH, and DFSPH) and provides different methods to simulate viscosity, surface tension, vorticity, and multiphase fluid interaction.PhysBAM is a multiphysics simulation library, capable of simulating rigid and deformable bodies, compressible and incompressible fluids, coupled solids and fluids, fracture, fire, smoke, hair, cloth, muscles, as well as many other natural phenomena.The PhysBAM library has a component called OpenGL Viewer, which displays and analyzes 1D, 2D, and 3D simulation data generated from PhysBAM projects in order to facilitate the fast verification and debugging of simulation methods.
SOFA is an open-source framework targeting real-time simulation, with an emphasis on medical simulation.It is based on about 15 years of research in physics simulation.SOFA is being used in many different projects, such as solid mechanics for the simulation of brain, ear, bones, heart, and liver; and fluid dynamics for the simulation of fat filling and blood flow in aneurysms.FleX is a particle-based simulation technique for real-time visual effects.It uses a unified particle representation for all object types, it enables new effects where different simulated substances can interact with each other seamlessly.It supports the simulation of rigid bodies, deformable objects, phase transitions, fluids, gases, and other phenomena.The goal of FleX is to use the power of GPUs to bring the capabilities of offline applications to real-time computer graphics.Realflow is a professional commercial fluid dynamics simulation software.It can be well connected with other 3D software, such as Maya, 3ds MAX, Cinema 4D, etc.It is powerful enough to simulate the fluid effects typically seen in high-speed macro photography; and to simulate viscous liquids such as cream, chocolate, oil, honey, and many others.Moreover, it allows the simulation of sand, snow, ice, and many other materials.It also provides a highly-optimized CPU and GPU particle solver where different types of materials are simulated within the same framework and are able to interact with each other.

A.2 Focus on modeling and rendering:
Highly accurate fluid simulation results require highly accurate rendering to convey the obtained simulation details.End-to-end systems cover most aspects of a graphics pipeline, including modeling, animation, and actual rendering.Here are some commonly used software for modeling and rendering.Houdini is a flexible node-based workflow following a dataflow computing model and allowing users to reuse computing nodes or even entire (sub)networks.Houdini has its own programming language, VEX, to deal with geometry for free development.The downside of Houdini is slow rendering.Blender is an efficient 3D modeling, rendering, and animation software.Its key advantages are hundreds of open source addons and extensive Python API.Every tool is available for scripting and customization.Both MAYA and 3ds MAX are Autodesk products that enable modeling, mapping, binding, animation, rendering, and more.MAYA's strength lies in animation and its special effects, mostly used for film and television animation; 3ds MAX has the advantage of being easy to model with and easy to learn.Both Maya and 3ds MAX are highly professional, accelerate workflows and provide stunning visuals.
Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 39 With the rise of the metaverse concept, NVIDIA Omniverse simulation platform is poised to launch the next wave of digitalization.NVIDIA Omniverse is an extensible open platform built for virtual collaboration and physically accurate real-time simulation.The advantage of Omniverse is that it facilitates real-time collaboration between users and applications, simplifying workflows by updating, iterating, and changing in real time without having to prepare data.It could change the way designers around the world collaborate and become the foundation of the metaverse.
a n c e d c o m p u t a t i o n a l a p p r o a c h e s

F i n eFig. 2
Fig. 2 Number of studied papers per identified topic.For each topic, the top (pure color) bar segment represents papers published in TOG, CGF, and TVCG.The bottom (shadowed) bar segment represents other key papers considered in our survey.

)
Navier-Stokes momentum equation.To further describe the motion of incompressible fluid flow, one can analyze the momentum of each fluid parcel.By introducing a momentum term ρu in Eqn.(1) and next using Eqn.(3), we obtain that ∂ρu ∂t

Fig. 4
Fig. 4 Schematic diagram of Eulerian grids.(a) Collocated grid where physical quantities are stored in the cell centres (yellow points).(b) Staggered grid where different variables are stored at different locations; in this example, pressure is stored at the cell centres (black points); velocity is split into its two Cartesian components and stored at the centres of the vertical cell edges (red and blue points).Subscripts I and J denote spatial indices.(c) Using bilinear interpolation to get the value of a physical quantity at any position.

Fig. 7
Fig. 7 Schematic diagram of particle-based adaptivity.Particlebased adaptivity adjusts particle size dynamically to reduce cost and preserve detail at the same time.

Fig. 9
Fig. 9 Schematic diagram of an octree represented in 2D where each cell has four children.3D octrees have eight children per cell.

Fig. 11 A
Fig.11A massive and sparse FLIP scene simulated using an efficient GPU parallel implementation[61] achieves one frame/second on an NVIDIA® Quadro GP100 GPU with 2 million particles and a 3360 × 160 × 2272-cell grid.

Fig. 13
Fig.13 Schematic diagram of different parallelization types.A single processing unit can be a multicore-CPU or a GPU.A CPU has few, but high-processing-power, cores.A GPU has weaker, but many more cores which can result in better performance for repetitive operations.

Fig. 14
Fig.14 Using Artificial Neural Networks (ANNs) to solve the pressure projection[95].The network inputs feature vectors extracted from training data to output pressure as close as possible to the ground-truth one.

Fig. 19
Fig.19 Schematic diagram of fluid-solid coupling using meshbased methods.Left: A sample MAC grid used in the fluid-solid coupling.Grey dots denote nodes of the solid regions.Fluid pressure values are stored at cell centers.Fluid velocity components are stored on cell faces.Right: A coupled solid-fluid system with the MPM method which allows to discretize the system with particles on a Cartesian grid and update the pressure based on DOFs on the grid.

Fig. 21
Fig. 21 Schematic diagram of wet hairs.Left: Wet hairs with a thin liquid layer flowing over their surface, and their motion is influenced by the surrounding liquid flow.Middle: The proximity and collision between wet hairs cause the adhesive force and contact force between hairs.Right: Cohesion of wet hair.Surface tension creates liquid bridges between closely positioned wet hairs.

Fig. 22
Fig. 22 Soluble and insoluble materials show different phenomena when coupled with fluid.Soluble and wettable granular materials [153].

Fig. 23
Fig. 23 Schematic diagram of non-mixing fluids.The lower-left part shows the anti-penetration force (F anti ) which is applied to particles in opposite directions along the surface normal.The embedding 3D fluid mesh can be re-tessellated in order to accommodate vertex displacement and produce changes in the interface topology.

c 1 c 2 Fig. 25
Fig. 25 Schematic diagram of mixing fluids.The lower-left part shows the fractions ci (volume fraction, mass fraction or concentration) of different phases in each particle during the diffusion process.The upper-right part shows the 3D fluid grid; black arrows represent grid forces due to diffusion.

Fig. 26
Fig. 26 Fluid extraction.From left to right, we see the change of solution due to greater solubility of the green liquid within the blue liquid than the red liquid.(a) At the beginning, the blue liquid and yellow mixture are put in a separating funnel.(b) The funnel is toppled.(c) Shaking vigorously mixes the fluids.(d) Turning the funnel upright results in a clear interface between the red liquid and a cyan mixture [169].

Fig. 29
Fig.29 Schematic diagram of surface tension.Liquids have forces between the same-kind molecules (cohesion) and between differentkind molecules (adhesion).The molecular force on the liquid surface is unbalanced, resulting in surface tension effects.

Fig. 32
Fig.32 River.A river flowing around sharp creases of a winding canyon creates a combination of calm and dynamic regions, including waterfalls and backdrafts.[198]©Weta FX.

Fig. 33
Fig. 33 Schematic diagram of spray and splashing.Spray and splashing are different in simulation scale.Compared with splashing, spray is composed of finer droplets.

Fig. 36
Fig. 36 Schematic diagram of isotropic and anisotropic particles.In contrast to isotropic particle fluids (left), anisotropic particle fluids (right) have a smoother surface.

Fig. 37
Fig. 37 Schematic diagram of the causes of vorticity dissipation and methods of refinement.

Fig. 38
Fig.38 Airflow over the delta wing of an aircraft showing a realistic and complex vortex distribution[256].

Fig. 39
Fig.39 Using TempoGAN[265] to generate high resolution smoke volumes from low resolution inputs.

Fig. 40
Fig. 40 Schematic diagram of scenario editing based on changing an existing simulation to achieve desired effects.

Fig. 41
Fig. 41 Schematic diagram of artificial effects: Artist-directed control to achieve non-physical effects.

( a )Fig. 42
Fig.42 Example of artificial effects: Fluid splashed on a butterfly shape generated by user sketching[288].

Fig. 43 1 Fig. 44
Fig. 43 Schematic diagram of media-directed formation: Images or videos are used to reproduce real-world scenes.
Such fluids do not follow Newton's viscosity law but rather show a non-linear relation between shear stress and the strain rate.For example, the viscosity of a Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 33 shear-thickening or dilatant fluid (e.g.starch paste) increases when the shear rate increases; the opposite happens for a shearthinning or pseudoplastic fluid (e.g.ketchup).Some non-Newtonian fluids have properties of solids, such as Bingham plastic fluids like toothpaste.Hence, it is impossible to use a constant viscosity for non-Newtonian fluids.Appropriate constitutive models are required to simulate such fluids.

Fig. 49
Fig.49 Schematic diagram of thin-film.Surfactants gather on both sides of the film.The tensile deformation yields an inconsistent surfactant concentration which leads to the surface tension gradient as shown in the figure.The film recovers its shape due to this gradient.
Physics-based fluid simulation in Computer Graphics: Survey, research trends, and challenges 37 also surveyed existing implementations used to compare and assess the quality of fluid simulations of various types; see also Appendix A.Several open challenges and directions for future work have emerged while analyzing the above seven themes, outlined next.