Visual Simulation of Detailed Turbulent Water by Preserving the Thin Sheets of Fluid

: When we perform particle-based water simulation, water particles are often increased dramatically because of particle splitting around breaking holes to maintain the thin ﬂuid sheets. Because most of the existing approaches do not consider the volume of the water particles, the water particles must have a very low mass to satisfy the law of the conservation of mass. This phenomenon smears the motion of the water, which would otherwise result in splashing, thereby resulting in artifacts such as numerical dissipation. Thus, we propose a new ﬂuid-implicit, particle-based framework for maintaining and representing the thin sheets and turbulent ﬂows of water. After splitting the water particles, the proposed method uses the ghost density and ghost mass to redistribute the difference in mass based on the volume of the water particles. Next, small-scale turbulent ﬂows are formed in local regions and transferred in a smooth manner to the global ﬂow ﬁeld. Our results show us the turbulence details as well as the thin sheets of water, thereby obtaining an aesthetically pleasing improvement compared with existing


Introduction
In the fields of smoke, fire, and water simulation, many techniques have been studied to improve the visualization detail of surfaces and flows. Expression of turbulence and vorticity in fluid flow is very important in Computer Graphics. Chen et al. [1] synthesized turbulence using a hybrid method. Zhao et al. [2] proposed random forcing method on an upsampled grid to synthesize fluid turbulence. The original vortex particles method [3] was extended by Yoon et al. [4]. The vortex particles are applied onto an upsampled grid with higher resolution. Kim et al. [5] proposed a method for buoyant sources by extending the vortex particle approach. They computed the source terms on the grid and the computed terms are mapped to the vortex particle system. Weissmann et al. [6] simulated vorticity using vortex filaments entirely, and source terms are used for obstacle-induced shedding. Barnat and Pollard [7] expressed vortex filaments based on Weissmann et al.'s method by using an improved re-sampling scheme.
The representation of turbulence details that arise from interactions between fluid and solid objects or large fluid bodies is a prevalent topic in computer graphics. In particular, the generation and preservation of small-scale details is an important aspect when animating turbulent fluids. Thus, physics-based [5,[8][9][10][11], noise-based, and texture-based [12][13][14] methods have been proposed for use in fluid simulations. There are adaptive-grid and multigrid methods for turbulent flows around solid objects [15][16][17] and self-turbulent flows [18,19]. It is possible to temporarily model the generation and evolution of small-scale vortices using the Navier-Stokes equations locally with a higher grid resolution. However, it becomes difficult to preserve the turbulence details if the vorticities leave the high-resolution region. Thus, various vortex particle methods [3,5] have been proposed based on the hybrid particle-grid framework to effectively express small-scale details of fluid.
In computer graphics, most existing particle-based fluid methods have tried to capture the details on the free surfaces of water, e.g., self-turbulent flows generated by a strong vortex [12,13,18,19] and turbulent flows around an object [15][16][17]20]; in contrast, this paper focuses on synthesizing the missing turbulence details while still preserving the fluid sheets. Recently, methods have been investigated for conserving thin fluid sheets, which represent small-scale details [21][22][23]. However, previous studies used methods that enforced mass conservation even when the number of water particles increased during the representation of thin fluid sheets [21,22]. Thus, they often assigned small mass values to regions with geometrically large volumes, so the stormy turbulent flows that are assumed to emerge in water simulation were oversmoothed because of the resulting low density (see Figure 1). To the best of our knowledge, modeling realistic turbulence formation and preserving thin fluid sheets in the global flow using a single method is still a challenging problem. In this study, we present a new method that achieves both of these goals in a particle-based fluid framework. Our method improves the simulation of turbulent water phenomena in a particle-based water simulator while using a fluid sheets solver. Small-scale turbulences are generated and maintained in local regions by adjusting the uneven mass distribution according to the volumes of the water particles. The creation of local turbulence and the preservation of thin sheets are combined in a novel manner. Thus, the ghost density of all the water particles, including newly inserted thin particles, is calculated to maintain the fluid sheets while preserving the turbulence details. The calculated ghost density of the increased water particles depends on the number of water particles that exist before checking the conservation of mass. To overcome any unsuccessful conservation of mass caused by the increased number of water particles, the ghost density is used to ensure mass conservation. Namely, the mass is geometrically redistributed according to the distribution of the water particles; thus, the turbulence details are preserved near the fluid sheets.
Problem Statement. In previous methods for preserving fluid sheets [21,22], every attribute except for the mass is linearly interpolated when the particles split. When splitting occurs to obtain a pair of particles, (i, j), the new mass is given by m = (m i + m j )/3. The new masses of particles i and j decrease to m new i = 2 3 m i and m new j = 2 3 m j after the split. The fluid velocity, u, and particle density, ρ, on the grid are affected by the changes in mass.
As shown in Figure 2, the water particles split near the fluid sheets to cover the breaking holes. From a visual perspective, the quality of the water surface increases but the mass of the water particles near the thin sheets decreases because of the aggressive mass conservation during the split. As the fluid sheets grow, the water particles split further and their masses decrease. As shown in the inset in Figure 2, the fluid motion appears to be oversmoothed around the region with turbulent motion, which is caused by the interaction between the fluid sheets and underlying water. To address this issue, we propose a fast and efficient method for synthesizing the missing water turbulence near the thin fluid sheets while still maintaining the conservation of mass. Hence, the proposed method obtains more aesthetically pleasing results compared with existing methods. Problem with a previously proposed thin sheet preservation algorithm [21,22] (scalar value: Mass of particles; resting mass: 1.0; gray rectangle: Interacting region between the underlying flow and thin fluid sheet areas; arrow: Velocity of particle).

Related Work
In this section, we briefly review previous researches related to the central issues addressed in this study, i.e., turbulence details and thin sheet details of fluids.
Turbulence Details. Stam [24] proposed an unconditionally stable semi-Lagrangian method for turbulence details that is used widely in the fluid simulation despite numerical dissipation. To express the small-scale features of smoke and overcome numerical dissipation, Fedkiw et al. [25] proposed the vorticity confinement method. However, it cannot calculate the amount of confinement force to be injected without causing instability. To ensure the stability, He et al. [26] proposed the adaptive vorticity confinement method by adjusting the coefficients for vorticity according to the vorticity magnitude.
To express the details of turbulence in smoke simulations, the vortex particle method is often used. To compute the velocity of turbulence, Park and Kim [27] solved the vorticity-related Poisson equation for a stream function. They evaluated the velocity field by tracking the gradient of the stream function. To express the details of turbulent flow, Selle et al. [3,5] integrated vortex particles using background grids. Angelidis and Neyret [28,29] used a vortex structure based on filament primitives to express the detailed features. Pfaff et al. [11] proposed a vortex sheet method being applied on a high-resolution triangular mesh to synthesize the detailed turbulent flow between interfaces. Some previous studies have used high-level turbulence models to express subgrid details. For example, Schechter and Bridson [30] captured energy cascading with a subgrid turbulence model. Kim et al. [13] expressed the turbulent flow using combination of the curl-noise function and a Kolmogorov wavelet form. Pfaff et al. [14] expressed the turbulence features near the boundaries of obstacles using a new energy model. They simulated vortex shedding using the distributed vortex particles near the boundary layers. They also extended this method to simulate anisotropic details of turbulence [12]. Mauricio et al. [31] accurately simulated vortex shedding by creating the vorticity directly on solid surfaces without modeling the boundary. Kim et al. [9] also expressed the stretching and wiggling motions of liquid with vortex stretching. Hong et al. [8] expressed the cellular patterns and wrinkling on the surface of a flame using the combination of detonation shock dynamics and the Navier-Stokes equation. Recently, Mercier et al. [10] substantially enhanced the details of the turbulent surface by adding a wave simulation to water surfaces. However, this method has been suffered from numerical dissipation causing the missing turbulence because it only improves the details of the original water surface just like a bump mapping method. The method adds a wave simulation near the surface particle to express the motion of the water including the valley with strong flow or turbulence. Advecting particles in this way makes it difficult to represent thin sheets or numerically lost turbulence due to motion dependent on the wave simulation, so it is not appropriate to apply directly to our method. Thin Sheet Details. The level-set method [32] is frequently used in the field of fluid simulation. In this method, the distance to the nearest surface, which is defined implicitly as a zero contour, is saved for each node. This method has been refined to obtain many other novel approaches; for example, Enright et al. [33], Wang et al. [34], and Mihalef et al. [35] attempted to reduce the numerical loss using Lagrangian particle approaches.
Bargteil et al. [36] increased the accuracy of fluid simulations by improving the signed distance field calculated from the reconstructed surface. Heo and Ko [37] preserved the details of thin surfaces using the spectrally refined level set and the high-order reinitialization method. However, these methods are insufficient for fully preserving the ruptures on a thin surface. As an alternative approach for solving this issue, mesh-based surface tracking has been explored in the fields of medical image analysis and fluid dynamics [23,38,39]. In these methods, explicit surface elements are typically advected via underlying flows to represent the thin features, but it is still difficult to address self-intersections and the complex topology changes that occur near interfaces. These problems can be temporarily addressed by resolving only the colliding parts into resampling meshes, but the algorithm becomes highly complex because of the diversity of possible situations.
There are various methods for reconstructing surfaces from the point cloud [40][41][42][43] used in particle-based fluid simulation methods, the fluid surfaces are defined implicitly based on the distances among particles and their density. However, these methods obtain poor results and blobby ruptures when the particles are sparse. Recently, Ando and Tsuruno [22] proposed an algorithm for preserving thin sheets in a fluid-implicit particle (FLIP)-based water simulation. In this method, particles are inserted into breaking holes to fully restore the thin sheets. They also extended this method to adaptively sample anisotropic particles and represent thin sheets more efficiently [21]. These methods well represent thin sheets of fluid. However, in the process of splitting water particles, inappropriate mass redistribution criterion causes water particles to have too small a mass, resulting in the problem of turbulence being lost.

Preserving the Thin Sheet and Turbulence Details of Water
In the proposed method, we generate thin fluid sheets and turbulence details using three-dimensional water particles with an existing particle-based water solver, i.e., FLIP [42]. Furthermore, we express the detailed turbulent flow around fluid sheets. This algorithm operates as follows.

1.
Water particles are advected and their density is calculated.

2.
Singular value decomposition (SVD) is used to extract thin particles from water particles.

3.
The distances and relative velocities between thin particle pairs are used to find candidate positions. 4.
The candidate positions are used to insert and break new water particles. 5.
To restore the missing turbulence, the ghost density is calculated for all of the water particles, including the newly inserted water particles. The ghost mass is then calculated based on the density to ensure the conservation of mass. 6.
The ghost mass is used to rasterize the velocities of particles on an Eulerian grid and advect the water particles using a FLIP solver. 7.
The fluid surfaces are reconstructed.

Extracting Thin Particles
The proposed method is based on Ando's methods [21,22] for turbulent flow modeling, which we briefly review here. First, the particle distribution is analyzed using anisotropic kernels. The particle directions are then determined by calculating the weighted average covariance per particle, C i , as follows: where d 0 denotes the initial space between particles, α 1 denotes a scaling constant for distance d 0 , and The covariance matrices, C i , calculated using the formulas above are used to obtain the eigenvectors and eigenvalues by SVD, as follows, from which the stretch and direction among adjacent particles are extracted: where e n denotes the principal axes ordered by the variance and σ n denotes the degree of stretch. The thin particles are extracted by examining σ 3 ≤ α 2 σ 1 , where α 2 denotes a threshold that controls the degree of thinness. Figure 3 shows the selected thin particles, and the difference between the eigenvalues increases as the water particles collide with the ground and spread out. As a result, the more strongly the particles are spread, the wider the selection range of the thin particles.

Finding Candidate Positions
The extracted thin particles are used to add particles to the breaking holes. We look for pairs of particles, (i, j), that connect the breaking holes. The midpoint of a pair, (p i + p j )/2, is considered as a candidate position for splitting. We select a pair that satisfies all of the following conditions: where α 3 and α 4 denote constants that determine the minimum and maximum distances between candidate positions, respectively, and u is the velocity of a particle. Equation (5) checks whether the two particles are located at a appropriate distance from each other.Equation (6) checks whether the candidate midpoints are sparsely distributed to add a new particle at radius α 3 . Equation (7) checks whether the pair of particles are moving away from each other over time. After these pairs have been identified, the midpoints are placed into list S without duplication.

Inserting and Removing Particles
When the candidate positions are very dense, they are not suitable for inserting new particles. In this section, we describe how to filter unnecessary particles and insert a minimal number of particles. I is the final list of the indices of inserted candidates. In this list, the first particle, I 1 , is the most remote candidate within S: where ρ j is the density of particle j. After determining I 1 , from S, we omit every nearby candidate around p I 1 that exists within a radius of α 3 d 0 . In the next step, we find nearby particles around p I 1 in S that exist beyond a radius of α 3 and add them to the list N I 1 . If we find a particle, the closest candidate is inserted within N I 1 : If the search is formulated as I n+1 = search(I n ), which means that I 2 can be found in I 1 , then this task can be repeated to find I 3 , I 4 , ..., I n+1 in a serial manner. If the next search fails, Equation (8) is used instead.
As mentioned in the problem statement, the attributes of particles are subject to interpolation after splitting. Thus, the following formula is used to calculate the mapped velocity, u, and particle density, ρ, in the Eulerian grid: where α 5 denotes a scaling constant for velocity u, and: In contrast, inserted particle i is removed if it satisfies either of the following conditions: where α 6 denotes the maximum density coefficient and α 7 denotes the minimum distance between particles. A particle is removed if it satisfies Equations (13) and (14). When a particle is removed, its mass is added back to the source particle before splitting.

Synthesizing Turbulence Details
The algorithm described in previous section fully represents the fluid sheets at breaking holes. However, the original mass of the water particles is distributed by aggressive conservation of mass whenever the particles split. A problem with this method is that the frequent splitting of water particles near the fluid sheets significantly decreases their mass. Thus, if the fluid sheets are larger, the water particles split more frequently and the mass values decrease. To illustrate this issue, the fluid sheets generated by stretching and their masses are colored in Figure 4, which shows that the mass value is smaller in the center because of the continuous insertion of thin particles caused by stretching. As shown in Figure 4a, the total mass value is increased due to the newly inserted orange particles, and the mass value around the orange particles is kept small in order to satisfy the mass conservation condition. In Figure 4b, blue particles have mass values close to zero, but they are expressed as a voluminous fluid surface, often resulting in awkward results.  Therefore, the stormy turbulence of water, which occurs when the fluid sheets collide with the underlying water (see Figure 5a), cannot be represented. In existing methods, the region where fluid sheets are generated has a smaller mass (see blue particles in Figure 5a), so any interactions with the underlying water barely change the mass, thereby hindering the generation of turbulence details and leading to oversmoothing. In contrast, the proposed method accurately captures the variations in mass at the point of interaction (see Figure 5b). Though the water particles composed of thin sheet particles are clustered, the intensity of the motion is weak because of their small mass values. However, due to the mass redistribution, our results well represent the motion of the water due to the interaction. We introduce ghost density and ghost mass to address this challenge. First, to calculate the mass while considering the newly inserted water particles, the ghost density ρ G is obtained as follows: where the mass term is omitted to determine the density according to the distribution of water particles. Thus, Equation (15) yields a weight that is used to provide a larger mass to the region where many particles are distributed geometrically. This value is then used to calculate the ghost mass m G : where m 0 denotes the initial mass of particles, which is set at 1.0, and σ target is a value used for the conservation of mass. Therefore, the conservation of mass is maintained through σ target during the simulation even though the number of particles changes. Equations (16)- (19) are used to distribute the mass value according to the distribution of water particles, thereby assigning a relatively larger mass to regions where water particles are congested; thus, the mass near the fluid sheets will not differ significantly. This distribution method aims to avoid ignoring the motions of particles with relatively low masses when the difference in the mass of the water particles is significant. The ghost density and ghost mass facilitate the robust restoration of turbulent water flows, which are blurred by previously proposed methods. Ando's method simply splits water particles to represent only thin sheets shapes [21,22]. The mass redistribution criteria used in the method cause turbulence to disappear in areas where thin sheets are interacted with the underlying water or coupled to the object. To reduce this problem, using the averages that you have suggested can be an alternative. However, it is likely that the turbulence flow will not be restored properly because only the mass of newly added thin sheets particles is updated in the thin sheets area. Equations (16)- (19) are formulas that redistribute mass while satisfying mass-conservation using the total number of particles including particles used to represent thin sheets. Figure 6 shows the distributed masses of the water particles obtained using the existing approach and the proposed method. The existing method simply distributes the original mass during the splitting process, so the thin particles have relatively low masses (see the aqua line in the figure). In contrast, the proposed method obtains refined masses by redistributing the mass according to the density of the water particles, thereby synthesizing the mass of the reduced thin particles (see the grey line in the figure). In the chart above, the 4th, 5th, 6th, and 7th particles in aqua color have relatively small mass values because they are particles split around thin sheets. Because this problem dissipates turbulence, the ghost density p G using Equation (15) and the mass-conservation Equation (15) are used to redistribute the mass. As a result, we avoided the problem of allocating a small mass value by distributing the mass at the density based on the number of particles, rather than simply splitting it. In the above chart, the mass value of the 4th, 5th, 6th, and 7th particles marked with grey color is increased because the number of particles near the thin sheets is relatively larger than that of the other regions.
In order to represent the turbulence of the lost water, the mass value is redistributed adaptively, but the total mass of the fluid is not changed by Equation (16)(17)(18)(19). In summary, while retaining the total mass of the fluid, it restores turbulence that has been lost only by local changes in the mass ratio. Figure 6 shows an example of changing the mass ratio of each particle while maintaining the total mass of the mass. Figure 6 is an intuitive illustration of the change in mass, and the results of the actual algorithm can vary. In Figure 5, the turbulence is lost in the aqua line in Figure 6, and the restored turbulence is the grey line in Figure 6.
In the proposed method, m G is used when rasterizing the velocity of water particles in the Eulerian grid and calculating the density of the water particles. According to this formula, stormy turbulence is well-represented when the fluid sheets collide with water, as shown in Figure 5b. (a) (b) Figure 5. Thin sheets and water-particle mass represented using the existing approach and the proposed method. Particles are colored according to their masss (red: High; blue: Low). (a) Existing method [21,22]; (b) Proposed method.

Implementation
The proposed method was implemented using an Intel i7-2600k 3.40 GHz CPU with 16 GB RAM and an NVIDIA GeForce GTX 480 graphics card. A FLIP-based [22] water simulation was executed for the underlying water, where every step of the water solver was performed in parallel.
All the physical quantities were stored for the FLIP grid using the staggered marker-and-cell method [44]. An additional grid was employed for surface reconstruction. We compute the level-set values from the simulated water particles [45]. Next, the triangular meshes were reconstructed using the marching cubes algorithm [46]. A variational framework was used for the water-solid boundary conditions [47]. The water-air interface was approximated considering the minimum distance from the air cell center, x air to each particle position x in the water cell, min ( x − x air − r), where r denotes the radius of a particle. We used a preconditioned conjugate gradient algorithm to solve this linear system.

Results
We tested various scenarios-dam breaking, a sphere dropping onto underling water, a rotating box in a pool of water, and a water box falling over a box cliff-to verify the preservation of details for thin sheets and turbulent flow while using the proposed method. Figures 7 and 8 show the visual improvement obtained by using the proposed method in the dam breaking simulation. The proposed method successfully generated smooth thin sheets and a dynamic water surface. The existing method [21,22] smoothed the water features around the thin sheets (see Figure 7a), whereas the proposed method achieved a better representation of the thin sheets and captured the water surfaces in more detail (see Figure 7b). As you can see, thin sheet surfaces, represented by the addition of thin sheet particles, can represent a clear fluid surface with no holes. However, turbulence around the thin sheet surface is lost because the mass value of the particle converges to a value close to zero. As a result, the smoothing surface is represented in Figure 7a, and our results show that we have successfully recovered the lost turbulence. Figure 8 shows a comparison of scenes where the thin sheets interacted with the underlying water. The existing method provided a good representation of the fluid sheets, but the turbulent details emerging around the thin sheets were not represented. As shown in Figure 8a, no secondary motion could be seen despite the interaction, whereas the proposed method represented both the fluid sheets and turbulent details simultaneously (see Figure 8b).  Figure 9 shows the simulation results obtained via the proposed method. In this simulation, a water ball was dropped into a shallow water pool. Observing the particle views shown in the inset image of Figure 9, the proposed method nearly resolved the numerical dissipation problem caused by mass redistribution, as the particles were advected more turbulently. The advected particles directly affect the turbulence details on the water surfaces, so the visual output was improved. Figure 10 shows a stirring box rotating continuously in a pool of water. Observing Figure 10a, the existing method failed to geometrically redistribute the mass according to the amount of particles; thus, the turbulent flows generated by the thin fluid sheets were missing, whereas the proposed method yielded improved turbulent flow. Figure 11 shows a box of water falling over the edge of a cliff. It should be noted that the jumping water sheets moving from the cliff edge hit the pool without any rupture. The proposed method greatly improved the generation of water sheets and turbulent splashing.
These applications to different scenes demonstrate that the proposed method performs well in terms of thin water sheets and turbulence details by redistributing the mass value so it visually decreases, as well as rectifying the awkward motions of thin sheets vanishing into water. The proposed method incurs a slightly greater computational cost compared with existing methods, but it produces a more detailed simulation of water sheets and turbulence. The parameters and configurations in this paper are summarized in Tables 1 and 2.   Maximum insertion distance 3.5 α 5 Velocity radius scale 1.0 α 6 Maximum thin particle density 0.2 α 7 Maximum thin particle distance 0.2 ∆t Time-step 0.006

Discussion and Additional Explanation
In this paper, we deal with the numerical dissipation problem occurring in the process of discretizing and numerically modeling Navier-Stokes equation, which is a physically based equation. The purpose of this research is to restore the thin thin sheets and turbulence flow of fluid, which was difficult to express by visual simulation due to this problem. So we provide a framework to improve fluid visual simulation by numerically modeling fluid thin sheets and turbulent flow.
First, we used FLIP solver [42], which is particle-based fluids, to calculate underlying water simulation, and extracted thin particles anisotropically to find the region where water particles move to irregular (see Section 3.1.1). Based on the position of the water particle, we calculated covariance matrix and PCA (Principal Component Analysis), and found orientation and stretching of water particles. This process roughly finds the region where holes occur, and finally finds the water particles present in the hole region through the conditions proposed in Section 3.1.2. In order to add particles from a region with a relatively large hole, we break a hole by adding particles from a region with a low density of water particles (see Section 3.1.3). At this time, the total mass of fluid is increased by newly added water particles, the constraint of mass conservation is not maintained, which is physically incorrect. In order to solve this problem, we used a simple scheme to divide mass into halves each time water particles are split. Using this mass transfer method, as shown in the Figure 4a, the water particle disappears because the water particle can not correctly express the turbulence flow due to the mass value reduced by the splitting scheme, despite the presence of many water particles. In order to solve this problem, we use Equation (15) to calculate ghost density ρ G . This value is the density calculated independent of the mass of the water particle, and density is calculated based on the volume of the particle, so you can restore the density value reduced by mass. In this process, ρ G of particles may differ from position to position, and mass conservation cannot be guaranteed. Therefore, we use ρ G to calculate ghost mass m G for maintaining mass conservation via Equations (16)- (19)). This formula redistributes the mass based on the ratio of ρ G calculated based on the particle volume and the sum of the ghost mass is normalized to be equal to the sum of the initial mass of fluid and the mass conservation condition is satisfied. As a result, ghost mass influences the calculation process of fluid's velocity, pressure, etc (see Equations (10) and (11)), so we can restore the turbulence flow lost in the splitting process.

Conclusions
In this study, we produce realistic turbulent flows around thin fluid sheets by redistributing the masses of particles. The proposed method considers the newly inserted particles to calculate the ghost density, which is then used to ensure the conservation of mass and to redistribute the mass of the particles. The calculated ghost mass affects the velocity step where it is mapped according to the density and the grid (see Equations (10) and (11)), thereby generating the turbulent flows and thus the motion of stormy water.
The proposed mass redistribution method is not based on physics, but it mathematically satisfies the conservation of mass to obtain turbulent flows. The proposed method handles turbulent flows by transferring the masses of water particles to particles near the fluid sheets, which causes the turbulent flows to decrease in some other regions. However, this issue is hardly visible because the mass is taken from a region where the density of the particles is low. The proposed hybrid particle-grid solver could be applied to wispy smoke, viscoelastic materials, and liquid ligaments to improve subgrid vortical details; these applications will be considered in future research.