1 Introduction

Fluid effects play an important role in the computer games industry: these are the complex interplay of various phenomena, such as convection, diffusion, turbulence and surface tension. Considerable research has taken place to improve the behaviour and performance and with the ever-increasing performance of hardware, simulations of the underlying physics can better approximate natural dynamics for representing water, waves, fire and gas. Hence, currently researches propose to simulate the interaction of multiple materials (fluid and solid objects). Physically, the Navier-Stokes Equation describes the state of the fluid, and there are many methods for solving incompressible flow by this equation.

1.1 Physical Equation

It is well known in Computational Fluid Dynamics (CFD) that the Navier-Stokes Equation that precisely describes the fluid’s acceleration by governs the behaviour of the fluid

$$\begin{aligned} \nabla \cdot \varvec{v} = 0 \end{aligned}$$
(1)
$$\begin{aligned} \frac{\partial \varvec{v}}{\partial t}=-\left( \varvec{v} \cdot \nabla \right) \varvec{v} -\frac{1}{\rho }\nabla p + \mu \nabla ^{2} \varvec{v} + \varvec{f} \end{aligned}$$
(2)

where \(\varvec{v}\) is the fluid’s velocity at a grid’s centre, p is pressure, \(\rho \) is density, \(\mu \) is the coefficient of viscosity and \(\varvec{f}\) includes any other external forces such as gravity or boundary confinement. These equations give a precise description of the evolution of a velocity field over time, and tell exactly how the velocity will change over an infinitesimal time-step as a function of the four terms (advection, pressure, viscosity and external-forces) in Eq. 2, but these equations do not consider the mass evolving in the computational physics transferring.

1.2 Related Work

The Eulerian grid-based fluid implementations have been popular for real time solutions as they provide a better description of the fluid properties. However, they do have a major disadvantage in that the grid must be fixed in the space, so making it difficult to track and depict the detailed behaviours (such as mist, foam and spray). The Navier-Stokes Equation was first used in [7] to animate gases and it produced good results on relatively coarse grids. In order to increase the details with lower memory consumption for grid-based simulation, [17] used a dynamic grid that had a low memory footprint when representing a high-resolution Level Set method. In [12], all quantities at cell centres were stored for an iterative solver, although a staggered MAC grid is robust and can make it easier to define boundary conditions [11, 20]. In order to reduce the effect of numerical dissipation caused by the use of an implicit semi-Lagrangian integration scheme [6], physical based approaches have considered to use momentum conservation in incompressible flow and the work in [14] can conserve the energy with the semi-Lagrangian method used in the simulation.

For the coupling between fluid and solid objects, there exist alternative concepts to incorporate boundary conditions for Eulerian fluids. [2] improved the FLIP method of [21] for two-way solid-fluid interaction. In that approach, the intra-pressure is formulated as a kinetic energy for the coupling problem. [10] presented a two-way coupling for deformable and solid thin objects, the algorithm using ray-casting to avoid fluid leaking through thin solids represented by mesh. [15] presented a GPU approach for the semi-Lagrangian scheme of [19] with arbitrary boundary conditions for the fluid simulation being generated and traced directly in Real-time.

For the GPU acceleration, NVIDIA’s Compute Unified Device Architecture (CUDA) has been applied to a large number of GPU-based fluid dynamics implementations, primarily in the engineering and scientific computing fields [9]. [13] proposed a CPU-GPU multigrid Poisson solver that exploits both the CPU and GPU to improve the performance and accuracy of the advection step. Later, the work in [3] presented a hybrid grid of two kinds of cell composition, and [4] described a novel gas simulation system that dynamically translates the fluid simulation domain to track the object and fluid surface.

1.3 Contribution

In this paper, we propose to maximize the performance of fluid simulation, focusing on a grid-based solution that can be efficient with parallel acceleration. For the Navier-Stokes Equation, we introduce a stable handling for the viscosity-friction term while updating the mass and momentum. The pressure-gradients term can be obtained by the Ideal Gas Law as suggested by [16]. Different from the traditional backwards advection method, we propose to use the forward-tracing-based motion and the evolving of fluid component (A) can be obtained by

$$\begin{aligned} A' = \frac{\partial A}{\partial t} + \left( A \cdot \nabla \right) A = -\frac{1}{\rho }\nabla p + \mu \nabla ^{2} A \end{aligned}$$
(3)

Consider the velocity evolving in Navier-Stokes Equation, the transmissions of velocities would be replaced by the transporting of momentum and mass (explain in Sect. 2). Further, we propose a fast component collection method for the fluid transport handling. This idea is designed for parallel implementation and the processing time taken is \(O\left( n \log n \right) \) which is suitable for real-time simulations. In addition, the kinetic energy transferring and conservation of momentum can calculate the coupling status, when solid object interaction is considered during fluid simulation.

2 Solution Method

Commonly, the Navier-Stokes Equation in Eq. 2 can describe the velocity field, but it is not enough to describe the density field only by the changing velocity. [19] provided an improved solution to this problem and implemented a gas simulation with the density being described by the Navier-Stokes Equation. We extend this idea into both the momentum and mass fields. We know that the relation between mass and velocity can be considered as momentum, and the law of conservation of momentum describes the energy transfer in nature. Thus, we use mass and momentum instead of the velocity then apply to Eq. 3 as follows,

$$\begin{aligned} \varvec{E'} = -\frac{1}{\rho }\nabla p + \mu \nabla ^{2} \varvec{E} \end{aligned}$$
(4)
$$\begin{aligned} m' = -\frac{1}{\rho }\nabla p + \mu \nabla ^{2} m \end{aligned}$$
(5)

where \(\varvec{E}=m\varvec{v}\) is momentum. In these equations, we focus on the internal (viscosity, pressure) status and external coupling forces as described in Sect. 2.3. After solving, the new velocity can be obtained by the new momentum and mass as

$$\begin{aligned} \varvec{v}^{+\varDelta t} = \frac{\varvec{E}^{+\varDelta t}}{m^{+\varDelta t}} \end{aligned}$$
(6)

2.1 Viscosity-Friction

The viscosity is an internal friction force that describes a fluid’s internal resistance to flow. This resistance results in diffusion of the momentum (also velocity, density, etc.). It causes the fluid’s components to move towards the neighbourhood balance. To solve the mass and momentum fields we use:

$$\begin{aligned} \varvec{E}^{+\varDelta t} = \varvec{E} + \int \varvec{E'} dt \end{aligned}$$
(7)
$$\begin{aligned} m^{+\varDelta t} = m + \int m' dt \end{aligned}$$
(8)

where \(\varvec{E'}=\partial \varvec{E}/\partial t\) and \(m'=\partial m/\partial t\). However, this method is unstable when the viscosity is large, so we refer to an efficient method in [19] for a discussion on the Gauss-Seidel Relaxation iterative technology, given by,

$$\begin{aligned} \varvec{E}^{+\varDelta t} = \frac{\varvec{E}_i + \mu \left( \varDelta x\right) ^2 \cdot \varDelta t \cdot \sum _J \varvec{E}_j}{1 + \mu \left( \varDelta x\right) ^2 \cdot \varDelta t \cdot \left| J\right| } \end{aligned}$$
(9)
$$\begin{aligned} m^{+\varDelta t} = \frac{m_i + \mu \left( \varDelta x\right) ^2 \cdot \varDelta t \cdot \sum _J m_j}{1 + \mu \left( \varDelta x\right) ^2 \cdot \varDelta t \cdot \left| J\right| } \end{aligned}$$
(10)

where J is the set of neighbouring grids of current i, j is the element index of J, \(\left| J\right| \) is the number of elements in set J. This method can avoid the density becoming a negative value, and the new velocity can be more stable and realistic than that obtained by only solving the velocity field directly in Eq. 6.

2.2 Pressure-Gradients

Commonly, the pressure is computed by solving an adequately built equations system using a Conjugate Gradient solver [20], which is heavy in computational load and memory consumption. Focusing on the fluid sample taken on a cell, pressure can be sampled on the centre of the cell. Such a scheme is chosen due to it having better stability properties than a scheme where the samples are taken from the same location. Thus, we invoke the Ideal Gas Law to calculate the recent pressure. It can be obtained from the current density as

$$\begin{aligned} p = k \rho \end{aligned}$$
(11)

where k is the gas stiffness that depends on the temperature. We know that there is a constant rest density in the material and this state is more obvious in liquid behaviour. Thus, we use a modified Ideal Gas Law suggested in [16], where the fluid internal pressure can be obtained from the current density as

$$\begin{aligned} p = {\left\{ \begin{array}{ll} k \left( \rho - \rho _{ rest } \right) , &{} \rho > \rho _{ rest }\\ 0, &{} otherwise \end{array}\right. } \end{aligned}$$
(12)

where \(\rho _{ rest }\) is the rest density in the material. The pressure will be minimized as the density approaches to this value. The condition \(\rho >\rho _{ rest }\) recognizes that the pressure is an internal repulsion force, and ignores the attraction force between the two nearest grids: this result/effect is more obvious in liquid behaviour. Consequently, the fluid should exhibit some internal cohesion formation, resulting in attraction-repulsion force as

$$\begin{aligned} \varvec{f}_i= -\frac{1}{\varDelta x \cdot \rho _i} \nabla p_i \end{aligned}$$
(13)

2.3 Energy Transfer

In addition to model the boundary conditions, we must handle the momentum along both the fluid and solid object. Since the conservation of the total momentum demands that, the total momentum before the collision is the same as the total momentum after the collision,

$$\begin{aligned} \varvec{E}^{+\varDelta t}_{ fluid } + \varvec{E}^{+\varDelta t}_{ rigid } = \varvec{E}_{ fluid } + \varvec{E}_{ rigid } \end{aligned}$$
(14)

In order to make sure the overall energy is conserved during coupling, we have to find the relation between momentum and energy items. Moreover, the conservation of the total kinetic energy can be expressed by

$$\begin{aligned} \frac{1}{2}\varvec{E}^{+\varDelta t}_{ fluid }\varvec{v}^{+\varDelta t}_{ fluid } + \frac{1}{2}\varvec{E}^{+\varDelta t}_{ rigid }\varvec{v}^{+\varDelta t}_{ rigid } = \frac{1}{2}\varvec{E}_{ fluid }\varvec{v}_{ fluid } + \frac{1}{2}\varvec{E}_{ rigid }\varvec{v}_{ rigid } \end{aligned}$$
(15)

Specially, note that gravity should directly add momentum to fluid as potential energy then be converted into kinetic energy. The bottom boundary forces can offset it symmetrically. Thus, we ignore the momentum changes based on gravity in these equations. By solving Eqs. 14 and 15, the evolved fluid velocity \(\varvec{v}^{+\varDelta t}_{ fluid }\) in the new frame of reference can be determined by

$$\begin{aligned} \varvec{v}^{+\varDelta t}_{ fluid } = \frac{\varvec{E}_{ fluid } + \varvec{E}_{ rigid } - C\left( m_{ fluid }\varvec{v}_{ rigid } - \varvec{E}_{ fluid }\right) }{m_{ fluid } + m_{ rigid }} \end{aligned}$$
(16)

where C is the coefficient of the restitution and slip belonging to the normal and tangential directions at the collided face. We can demonstrate the scale of coefficient (\(C \in \left[ 0.0,1.0 \right] \)) for the handling various boundary conditions, as well as complex external constraints such as perfect elastic and free slip effects.

3 Implementation

In order to do simulation in real-time, we prefer to implement the fluid model in a GPU to accelerate our computation by parallelism.

3.1 Fluid Diffusion and Distribution

In the Eulerian grid-based approach, the fluid behaviour is described in terms of a fixed grid. Fluid components cannot be transported by moving these grids directly. To consider the cell size and diffusing the fluid mass (m) would lose the details in vorticity, so we change to diffuse the density (\(\rho =m/\varDelta x^d\), d is the dimension) instead. Consider the advection term \(-\left( \varvec{v}\cdot \nabla \right) \varvec{v}\) in Eq. 2, the minus symbol means these fluid components should be found by backward tracing [19]. However, this method may cause some mass loss. If mass (also momentum) is lost and there are few areas of coupling to exchange energy to, changes in energy may cause undesirable noise. To avoid this issue we introduce a forward tracing method that can conserve the total \(\rho \) throughout the simulation.

Fig. 1.
figure 1

These figures illustrate the forward tracing steps: after these components have been transported, they should be assigned to their nearest 4 (in 2D, 8 in 3D) grids.

Full size image

As shown in Fig. 1 in a 2D case, fluid material moves (forward) to a new position. These source components (density, momentum) have been divided by Linear-Interpolation and added to the nearest four grids (e.g. \(\rho _a\) will distribute to grid 8, 9, 13 and 14). Assume \(\rho _a\) is the current density in grid 21, the source \(\rho _a\) will divide to \(\rho _a^8\), \(\rho _a^9\), \(\rho _a^{13}\) and \(\rho _a^{14}\) respectively, with \(\rho ^8_a + \rho ^9_a + \rho ^{13}_a + \rho ^{14}_a = \rho _a\). Additionally, the final density of grid 8 should be \(\rho ^8_a + \rho ^8_b + \rho ^8_c\). Every grid must repeat this operation to obtain the final state.

3.2 Fluid-Grid Relation Table

In the forward tracing method, all grids must wait for fluid movement to finish, then each grid should evolve to the sum of the fluid components that have arrived, but it is difficult to distribute momentums and density to the grids efficiently. In our proposal, we define a pair of index arrays to represent the Fluid-Grid relation between the grid and the arrived fluid after diffusion (see Fig. 2).

Fig. 2.
figure 2

Flow chart of sorting and grouping the fluid density in the same grid.

Full size image

As shown in Fig. 2, array (a) stores the index of the grid that receives the fluid density and array (b) stores the sorted index of the density relevant to the corresponding grids. Sorting is by the parallel radix algorithm in the GPU environment and the performance is \(O \left( n \log n \right) \). According to these two arrays, we must provide the (First, Number) pair for each cell coordinate to record the table index of the first related density and the number of related densities in the sorted table.

Finally, every cell and its received density are grouped into the same region of the table, thus there is only one loop of the sum of densities and we know exactly how many densities must be queried in our method. This makes it suitable for parallel programming design and implementation. Note that using our conservative advection with diffusing can conserve the total mass throughout the simulation.

4 Experiment Results and Discussion

These experiments has been implemented on a PC with 2 GB memory with video card NVIDIA Quadro 6000 GPU. These simulations are implemented by CUDA 7.5. The scene was rendered by OpenGL 4.5, GLSL 4.5 in 2D and OptiX 3.9 in 3D.

4.1 Results

As shown in Fig. 3, the smoke is running in 2D and the cycle time-step is less than 1 / 60 s with the Courant-Friedrichs-Lewy (CFL) condition. Note that there is no gravitational effect in smoke, thus external forces only include the boundary and coupling forces. We use Eq. 11 for computing the smoke pressure; it illustrates that while the smoke rises by buoyancy, the potential energy will decrease with increasing kinetic energy and that will apply to the energy transference. However, there is no mass loss in our method, but it does cause some areas to contain undesirable dispersion-density/noise that is obvious in the smoke simulation. To alleviate this problem we should adjust neighbouring smoke-grid (For example in Fig. 1, very low density of smoke \(\rho _c^7\) in grid 7 has prorated to its neighbouring smoke-grid 2 and 8).

Fig. 3.
figure 3

These figures show the density of smoke at different time-steps. This simulation uses our method with the smoke injected from below. It processes at a resolution of \(500\times 900\) grids scale and coupling with two rotating stars.

Full size image
Fig. 4.
figure 4

These figures show the interaction between water and solid balls in different time-steps. This simulation processes at resolution \(128\times 128\times 128\) grids scale in 3D, and the fluid surface is constructed by the marching cube algorithm [8]. (Color figure online)

Full size image

As shown in Fig. 4, water is running in a 3D scene. We use Eq. 12 for computing the water pressure, and we also add some vorticity confinement [18] to the system for momentum conservation. Note that a particle level set method [5] is used to advect and treat the boundary tracing. For forward advection we need to treat the level set as a solid object, since there is no guarantee that the particle level set method had been also used in both time t and \(t^{+\varDelta t}\). Furthermore, the water would be affected by gravity and the upper grids energy will increase so that the energy becomes very large at the bottom. This is non-conservative from a global perspective. To solve this we produce a reaction vertical momentum after the gravity is added to these bottom grids.

These GPU-based fluid simulations with momentum conservation can be simulated in real time. The results show that our method can be used in gas and liquid with conservation of energy. The time performance of the experimental results is shown in Table 1.

Table 1. Performance results of GPU-based of fluid simulation with momentum conservation.
Full size table

4.2 Comparison

For rendering purposes, we improve the visual quality of the simulation by generating particles that penetrate the level set surface, then applying a local momentum conservation force to slightly alter the level set velocities in fluid surface regions. The comparison between our method, the methods using the GPU-based stable fluid in [1], and the mass and momentum conservation fluid [14] is shown in Fig. 5.

Fig. 5.
figure 5

Comparison of the simulation speed in FPS.

Full size image

Compared with [14], a speedup of around 10% to 20% can be achieved. Although the speed is slower than in [1], a reason is that the fluid surface rendering is also time-consuming, but there is not much difference and our method can avoid the mass loss and converse the energy of interaction.

5 Conclusion

In this paper, we have presented a novel momentum conservation method for simulating a wide variety of fluid behaviours. Our algorithm can avoid mass loss during the energy transfer, and the solution, designed for parallel computing, can be implemented in a Graphics Processing Unit (GPU). Moreover, a fast forward-tracing-based motion has been presented to handle the distribution of the fluid components. With the momentum conservation, we also considered the conservation of energy for various boundary conditions, including perfect elastic and free slip effects. The experiment results showed that our method can be around 10% to 20% faster than the previous momentum conservation fluid simulation and this is significant.