A Symmetric Particle-Based Simulation Scheme towards Large Scale Diffuse Fluids

: We present a symmetric particle simulation scheme for diffuse ﬂuids based on the Lagrangian Smoothed Particle Hydrodynamics (SPH) model. In our method, the generation of diffuse particles is determined by the entropy of ﬂuid particles, and it is calculated by the velocity difference and kinetic energy. Diffuse particles are generated near the qualiﬁed diffuse particle emitters whose diffuse material generation rate is greater than zero. Our method ﬁts the laws of physics better, as it abandons the common practice of adding diffuse materials at the crest empirically. The coupling between diffuse materials and ﬂuid is a post-processing step achieved by the velocity ﬁeld, which enables the avoiding of the time-consuming process of cross ﬁnding neighbors. The inﬂuence weights of the ﬂuid particles are assigned based on the degree of coupling. Therefore, it improved the accuracy of the diffuse particle position and made the simulation results more realistic. The approach is appropriate for large scale diffuse ﬂuid, as it can be easily integrated in existing SPH simulation methods and the computational overhead is negligible.


Introduction
Diffuse fluids include spray, foam and bubble [1,2].They are widely presented in metallurgical engineering, industry analysis, virtual reality game, and many other fields.Diffuse materials usually attach to the fluid's boundary or splash around.These phenomena indicate the shape of the boundary and the tendency of the flow.The addition of these details can make the fluid simulation more realistic (Figure 1).However, diffuse materials have many complicated physical natures.It is extremely hard to simulate diffuse fluids exactly.On one hand, diffuse materials are usually generated by the rushing fluid flow which is hard-to-describe.In addition, they need to be deleted dynamically because of the uncertain survival time.On the other hand, the coupling between diffuse materials and fluid is difficult to be modeled.
In fluid simulation field, the Euler method based on the point of flow field space is difficult to model fluids with great deformation due to the limitation of fixed grid [3].However, the Lagrange method based on fluid particle is suitable for simulating splash scenes with diffuse materials [4].In this paper, a symmetric particle-based method named Smoothed Particle Hydrodynamics (SPH) [5] is used to handle this procedure.SPH method is one of the Lagrange methods.There are two main difficulties to simulate diffuse fluids.One is how to model the generation and termination of diffuse particles [6].Due to the unstable state of diffuse materials, they need to be defined as dynamically surviving particles.The complex generating conditions and deleting conditions should be established for diffuse fluids.
Another is how to handle the coupling between diffuse materials and fluid efficiently [7].Small-scale scenes have higher requirements on the surface characteristics [8,9], while large-scale scenes pay more attention to the accuracy of the description of the movement tendency [10].At the same time, the number of particles in a large scene increases exponentially, and the process of simulation will consume much time [11].A dam burst hit a ball from the right side.The 60th frame of the experiment.Render result without diffuse materials (left).Our method (right).When the dam impacts the ball, curved water will be formed around the ball.Compared with the non-diffuse-material results, the rendering result with diffuse materials can capture the direction, shape and other details of fluid more clearly, and the simulation result is more realistic.
A new efficient diffuse fluid simulation scheme based on SPH is proposed in this paper.The generation of diffuse particles is modeled by the entropy of the flow field, which is calculated from the velocity difference and kinetic energy.When the disorder, or randomness of the flow field is large, the diffuse material generation rate of certain fluid particles will be greater than zero, and these particles will be designated as diffusible particle emitters.Diffuse particles are generated near these emitters randomly.To enhance efficiency, the coupling between diffuse materials and fluid is achieved by the velocity field.The weight of the coupling depends on the degree of the interaction between fluid and diffuse materials.This step is a post-processing step, ignoring the extremely small forces exerted by the diffuse material on the fluid, and reduces the time spent on the process of finding neighbors.

Related Work
Diffuse materials play an important role in fluid simulations, metallurgical industry and other fields.Although the regular movements of fluid can be successfully achieved with both Eulerian and Lagrangian techniques [4].The Lagrange method is more suitable for simulating the fluid with great deformation because its particle-based system does not have the limitation of fixed mesh [12].The schematic diagrams of Eulerian and Lagrangian method shown in Figure 2. In this paper, one of the Lagrange methods named Smoothed Particle Hydrodynamics (SPH) is extended to deal with diffuse fluid.SPH method was first used in the field of astrophysics [13,14].It was later introduced by Desbrun et al. [15] into the field of computer graphics to simulate the deformable body.In 2003, Müller et al. [12] used SPH method for real-time interactive simulation of fluid.However, since the ideal gas equation was employed to solve the pressure, the incompressibility of the fluid cannot be guaranteed.In the field of Computational Fluid Dynamics (CFD), Monaghan [16] proposed a SPH method for free surface flow, replacing the ideal gas state equation with Tate equation.This method achieved a weak compressibility of the fluid by assuming a large sound speed and stiffness coefficients, which is very robust and easy to program.In 2007, Becker et al. [17] introduced this idea into the field of fluid animation.Weakly Compressible SPH (WCSPH) significantly increases the realism of the simulation, but it is at the expense of limiting the time step, which affects the efficiency of the algorithm.Physics-based fluid simulation and Computational Fluid Dynamics seem to be similar in mathematical expression, but the two are essentially different.The most significant difference is that fluid simulation as a branch of computer animation pursues realistic picture effects and computational efficiency [18], while CFD as a branch of numerical calculation pursues computational accuracy [19].In 2009, Solenthaler et al. [20] proposed the Predictive-Corrective Incompressible SPH (PCISPH) approach, using a predictive-corrective method to adjust the pressure by density which avoiding solving complex equations.PCISPH method improved the fidelity and computational efficiency, and it made a significant progress in the field of fluid simulation based on SPH.When updating the position of the particle, the explicit method was originally used for the calculation [12,17,20], which is, using the speed of this time-step to solve the position of following time-step.In 2014, Ihmsen et al. [21] proposed an implicit method to calculate the particle's position of the next time-step with the particle's velocity of the next time-step which making the calculation more accurate.At the same time, the paper gives a new iterative method, and the simulation efficiency has been further improved.In 2015, Bender et al. [22] drew on the idea divergence-free velocity method, enforcing incompressibility of fluid both on position level and velocity level.This method was first proposed by Cummins et al. [23] in 1999, projecting the velocity field onto a divergence-free space by solving a pressure Poisson equation.To improve the ISPH modeling capacity, Gui et al. [24] combined the standard density-invariant and the velocity divergence-free formulations in a weighted average form and proposed a mixed source term model.This model can obtain more accurate impact pressure and force as compared with the results obtained by using either the density-invariant or the velocity divergence-free ISPH model.In 2017, Winchenbach et al. [25] proposed an adaptive SPH method, which changed the previous fixed-mass setting of particles and set the best mass value according to the distance of the particles from the free surface, so that the simulation result effectiveness increased.
The SPH method has been continuously revised to improve its symmetry and accuracy.In the 1980s, Monaghan [26] presented a symmetrization formula to achieve better results.The method is applied to the derivation of the SPH equations which conserve linear and angular momentum exactly.In 1996, Johnson et al. [27] set forth a normalized smoothing function (NSF) algorithm.This approach can obtain exact normal strain rates for conditions of constant linear velocity distributions.In 1999, Bonet et al. [28] presented a discrete variational SPH framework which ensures the balance of linear and angular momentum.In 2008, Khayyer et al. [29] proposed a corrected incompressible SPH (CISPH) method to track water-surface in breaking waves accurately.Corrective terms are derived for the incompressible SPH formulations based on a variational approach, which ensures the preservation of angular momentum.At present, the SPH method can well simulate general fluid dynamics problems.
As a particle-based method, the SPH method is often used to simulate turbulence and sea waves.Back in 2004, Gotho et al. [30] combined the incompressible Smoothed Particle Hydrodynamics (SPH) method with a Large Eddy Simulation (LES) model to study the wave transmission and reflection by a half-immersed curtain breakwater.For the plunging waves, Shao et al. [31] presented a 2-D large eddy simulation (LES) modelling approach in 2006.The computations are in good agreement with the documented data, especially the computed turbulence quantities under the breaking waves.In 2015, Shadloo et al. [32] simulated tsunami runup on idealized geometries for the validation and exploration of three-dimensional flow structures in tsunamis.This study proves that SPH is able to reproduce the runup of long waves for different initial and geometric conditions.In June of 2017, Bender et al. [33] introduced a novel micropolar material model for the simulation of turbulent inviscid fluids.This model is linear and angular momentum conserving, and its computational overhead is negligible.
For diffuse materials, there are many simulation models focused on the behavior of microscopic bubbles [8,34].In 2011, Ihmsen present a particle-based model for the complex bubble flow with respect to deformation and merging of bubbles and volume-dependent buoyancy [35].Chentanez used a hybrid grid representation composed of regular cubic cells on top of a layer of tall cells, and focused effort on the area near the surface where it most matters [36].In these methods, air is considered as a separate phase and the behavior of diffuse fluids is understood as interaction between multiphase flows.However, the simulation of large-scale macro-diffusible fluids remains to be studied.In 2003, Takahashi et al. [37] proposed a method for foam simulation, which combining grid and particle methods to add foam at the crest of wave formation.Although this method had certain visual effects, it ignored the physical reality and do not modeling from the perspective of bubble formation.Due to the lack of theoretical support, the resulting effect has some limitations.In 2009, Mihalef et al. [38] used kinetic energy and surface energy to determine the generation of small-scale splash, which could simulate droplets and bubbles realistically.In 2016, Meric proposed to analyze the acceleration and curvature values to identify the fast and complicated water flows, and foam particle is advected by its classified type based on its velocity [39].This paper is inspired by those methods mentioned above and a symmetric SPH-based simulation model for diffuse fluid is proposed.Experiments demonstrate that the addition of diffuse materials makes the results more vivid.The simulation results of our method can show the shape and movement trends of fluid more clearly with a negligible computational overhead.

Fluid Simulation Model
The addition of diffusion material can enrich the abundance of fluid simulation scene and enhance the fidelity of simulation results.In this paper, the calculation of diffusion material is a post-processing step.The physics-based fluid simulation process needs to be modeled first.This section expounds the physical constraints of the fluid and the SPH simulation method.

Physical Constraint
Smoothed Particle Hydrodynamics is a particle-based method for fluid simulation.In this scheme, fluid is modeled by a number of particles with attribute such as mass m, position x, and velocity v.Each moving particle is constrained by three conservation laws, including conservation of mass, conservation of momentum and conservation of energy.These constraints are expressed as governing equations of fluid in the field of fluid mechanics: where ρ denotes the density of fluid, t denotes the time, v is the velocity, p is the pressure, σ is the viscous stress tensor, and a ext is the acceleration generated by the external force, including gravity, surface tension, adhesion and other forces.Equation ( 1) is derived from the law of conservation of mass to ensure the continuity of fluid.Since mass is an inherent property of the particle and the number of particles in the scene remains constant, the law of conservation of mass is automatically satisfied.
Equation ( 2) is the Navier-Stockes equation (N-S Equation) derived from the law of conservation of momentum, and it is used to describe the fluid's motion.It can be seen that the acceleration of a particle is mainly contribution by the pressure acceleration a pressure = − 1 ρ ∇p , the viscous acceleration a viscous = 1 ρ ∇σ, and external acceleration a ext [40].The external force is applied to the particles directly, the viscous force is often calculated by a user-defined artificial viscosity, and the pressure is calculated by the equation of state derived from the law of conservation of energy: where γ = 7, and c denotes the velocity of light.Equation ( 3) is the Tate equation for ensuring the incompressibility of fluid, which can be calculated efficiently.Compared with the ideal gas state equation p = k(ρ − ρ 0 ), this equation solves the problem of high compressibility.Although this equation requires smaller time steps compared to the Poisson equation ∇ 2 p = ρ ∇v ∆t , the computation per time step is significantly faster and the implementation easily fits into the SPH framework [17].
The stress state of any point in the continuous flow field can be obtained by solving the above three equations.In addition, to implement fluid simulations, the discretization method is also needed to solve the stress information at each discrete point.

Smoothed Particle Hydrodynamics
The SPH simulation scheme is essentially an interpolation method.Each particle will carry attributes including mass, position, velocity, etc (Figure 3).The continuous flow field is dispersed as a large number of moving fluid particles.The interpolation result A i of particle i at position x i is contributed by all the particles j in the neighborhood, with the weight W related to the position: where m is the mass which is constant, W ij = W(x i − x j , h) is the smooth kernel function with homogeneity, and h is the support radius.In this paper, we take W as a B-cubic spline function to ensure the stability, accuracy and efficiency of SPH method.Fluid is considered to be composed of movable particles that carry physical attributions.The physical attributions of a particle are weighted in summation by its neighbors.
By substituting density into Equation ( 4), we can get the density at location x i , which is the last unsolved attribute: The SPH method has the inherent problem of force asymmetry.Muller [12] solve the problem by using the idea of arithmetic averaging.The new interpolation function is as follows: Equation ( 6) takes the property values of particle i into account, ensuring the force symmetry.In addition, the gradient of A i is written as Discretize the N-S Equation ( 2) by Equation ( 6), and the force exerted by neighbor particles on one particle can be obtained.
In general, the conservativeness of governing equations of fluid is guaranteed by the symmetry of space translation, time translation and rotation transformation.The symmetry of force in the SPH method is guaranteed by the combination of interpolation method and arithmetical mean.Dispersing the flow field into a finite number of moving particles and solving the governing equations of fluid by SPH method, the state of flow field will be obtained.The calculation of diffuse materials is a post-processing step, which will be elaborated in the next section.

Diffuse Material Simulation System
Diffuse materials are usually generated by the rushing fluid flow.When the disorder, or randomness of the system become larger, the air will be dragged out and resulting in diffuse fluids.Although current SPH-based simulation methods could describe regular movements of fluid, details such as diffuse materials were often ignored or simplified during modeling process.There have been many studies focusing on the process of bubble in small scenes.However, there are still some problems to be solved for the macroscopic movement of diffuse materials in large scenes.
In this paper, an efficient diffuse fluid simulation scheme based on SPH is presented.Diffuse materials are considered when the diffuse material generation rate calculated by the entropy of the flow field is greater than zero.They are generated by the qualified diffuse particle emitters in fluid particles.The coupling between diffuse materials and fluid is a post-processing step achieved by the velocity field.Three coupling parameters are used to indicate the degree of coupling between the diffuse material and the fluid.Therefore, when the calculation of the flow field by the model above is completed, the generation and termination of diffuse materials and the coupling with other materials will be processed.This approach avoids the time-consuming process of finding neighbors and enhance the efficiency.

Generation of Diffuse Material
Diffusion Entropy.In some early approaches, diffuse materials are added at the wave crest empirically based on visual effects.However, these methods do not obey the physical laws.In nature, air will be dragged into turbulent fluid flow, resulting in spray, foam and bubble.Slowly flowing streams are always crystal clear, while a huge splash will be sparked when the wave crash into the rock in high speed.Therefore, we consider fluids with more kinetic energy and more chaotic internalities to have a greater ability to produce diffuse materials.The degree of disorder can be expressed by the velocity difference between particles.The simulation of diffuse fluid is processed from the perspective of physics.
The concept of entropy is introduced to describe the potential of flow field to generate diffuse materials.Entropy is a measure of disorder within a macroscopic system or of the availability of the energy in a system to do work [41].In statistical mechanics, entropy is a measure of the number of ways in which a system may be arranged, often taken to be a measure of "disorder".The higher the entropy, the higher the disorder [42].For fluid, entropy is the degree of disorder within the flow field and the energy it contains.Therefore, the diffusion entropy S D of the fluid particle i is related to the kinetic energy E and velocity difference v di f f : The kinetic energy of particle i is as follows: The velocity difference of particle i can be calculated by interpolation: where v ij = v i − v j is the velocity difference between particle i and particle j, W D is the kernel function with normality: where x ij = x i − x j is the position difference vector, h is the support radius.W is an even function whose function image is symmetric about the longitudinal axis.Equation ( 10) can calculate the velocity difference between particles i and its neighbors, but it cannot accurately represent the generation probability of the diffuse material.Under the same velocity difference, more air will be mixed in as the two particles move toward each other.Conversely, if two particles are moving away from each other, the amount of diffuse material produced will be less.Therefore, Equation (10) is not sufficient enough to express the amount of bubbles generated by the velocity difference.
In the SPH scheme, vij • xij can be used to indicate whether the two particles are moving away from each other or close to each other, where vij = is the direction vector of the velocity difference vector, and xij = is the direction vector of the position difference vector (Figure 4).The value of vij • xij is negative when particles are in opposite directions, and positive when particles are backwards.Therefore, we use the parameter λ to determine the additional weight of the velocity difference: The corrected velocity difference equation is as follows: The diffusion entropy of each fluid particle is calculated by the product of kinetic energy and velocity difference.Obviously, the higher the diffusion entropy is, the higher the probability of generating diffuse particles in the vicinity of the fluid particles is.Therefore, the number of diffuse particles generated around each fluid particle is related to the entropy.
Number of Diffuse Particles.In order to solve the number of generated diffuse particles, the mapping function φ is used to map diffusion entropy S D to [0, 1] interval to obtain the diffuse material generation rate R D .The mapping function for a scalar quantity M is defined as: where ϕ min is the upper limit, ϕ max is the lower limit.Velocity difference and kinetic energy are mapped respectively to obtain accurate calculation results.The diffuse material generation rate R D is the product of these two value: where are the user-defined parameters.Throughout the flow field, the diffuse material generation rate R D of each moving particle is calculated.When the R D is greater than zero, the particle is considered as a diffuse particle emitter.Diffuse particles are generated in the vicinity of these emitters randomly.The number of diffuse particles produced in each time step is: where ∆t is the time step, k max is the maximum sampling rate, that is, the maximum number of diffuse particles that can be generated in each time step.Diffuse Particle Emitter.Diffuse particles are randomly generated near the fluid particle selected as diffuse particle emitter.In one time-step, the fluid particles move from position x f (t) to position x f (t + ∆t).In addition, the diffuse particles are generated in the cylinder centered on the line connecting x f (t) and x f (t + ∆t).Its default height is ∆x f = v f • ∆t , and the default bottom radius r is related to the support radius of fluid particles.
In order to make the generation of diffuse material more realistic, the initialization parameters need to be set randomly.Firstly, a vector n 1 perpendicular to vector ∆x f is chosen randomly.Then, three random variables k h , k r , k θ are set to produce another vector n 2 that is not linearly related with n 1 and ∆x f .The exact location of diffuse particle i is determined by the following parameter: where h i is the height, θ i is the azimuth, r i is the distance to the cylinder (Figure 5).When the values of k h , k r and k θ are between 0 and 1, the diffuse particles are generated inside the cylinder.In experiments, when the support radius and time-step is small, the size of the cylinder becomes very small.Newly generated diffuse particles will accumulate together, causing undesirable visual effects.Therefore, the value of the above three parameters can be adjusted according to experiments in order to achieve the best results.The vector n 2 is calculated by Rodrigues' rotation formula [43]: where v f is the direction vector of v f .The initial position and velocity of a diffuse particle are as follows: where, n 1 and n 2 act as position in Equation ( 21), and velocity in Equation (22).In summary, the number of newly generated diffuse particles per unit time-step is determined by the entropy of the fluid particles, which is calculated from the velocity difference and kinetic energy of the fluid particles.The diffuse material emitter will be designated, and its diffuse material generation rate will be determined by mapping the entropy to [0, 1] interval.The maximum sampling rate is set to control the number of diffuse particles generated.Random parameters are used to initialize diffuse particles, making the distribution of particles more random and experimental results more realistic.

Interaction with Fluid
The traditional force analysis of microscopic diffuse fluids simulation is complex.Many surface properties are considered including surface tension, boundary handling, etc.In a small-scale scenario, a two-way coupling model needs to be established for the interaction between the fluid and the air, including the description of the surface characteristics at the two-phase interface.The time-consuming cross-neighbor lookup process between fluid particles and foam particles is unavoidable.However, large-scale scenes emphasize the tendency of macroscopic movements, while the requirements for surface characteristics are not high.Therefore, we proposed a simplified one-way coupling model designed to quickly predict macroscopic motion of diffusible materials.
Diffuse particles are essentially air particles.The mass carried by these particles is much smaller than the fluid particles nearby.Therefore, it is difficult for diffuse materials to possess enough momentum to have a significant effect on the movement of the fluid.The force that the diffuse material exerts on the fluid is ignorable.In our method, the calculation process ignores the negligible effect of the diffuse material on the fluid and only solves the effect of the fluid on the diffuse particles.Besides, the number of particles in large scenes increases exponentially, and the process of two-way coupling will consume much time.Our model addresses macroscopic simulations of large-scale scenes and deals with the interaction of diffusible materials with fluids through a one-way coupling model.The coupling between the fluid with the diffuse materials is achieved by the velocity field.In our model, only the fluid neighbors of the fluid and the fluid neighbors of the diffuse materials need to be searched, which greatly saves the simulation time.This process does not require density calculation, avoiding the problem of distortion of the simulation result due to the density instability at the two-phase interface in two-way coupling models.In addition, it is a post-processing step in which the movement of the diffuse materials is not processed until after the state of the fluid has been calculated.Therefore, our method can be easily integrated in existing SPH simulation methods.
The velocity of a diffuse particle is affected by external forces, buoyancy and fluid coupling forces.In the broken wave, part of diffuse particles detaches from the particle cluster as independent flying particles, known as spray.Sometimes, air is trapped in fluid and forms bubbles that rise gradually under buoyancy.Bubbles are also transported by the fluid due to the frictional force.The foam is formed when rising bubbles reach the surface, and it is advected by the water.Therefore, we set three parameters k ext , k b , k cp to represent the weight of these forces on the velocity of the diffuse particle.The weighted average velocity of the fluid particles around a diffuse particle d at x d is: where x f is the position of fluid particles, v f (t + ∆t) = x f (t+∆t)−x f (t) ∆t . The velocity of a diffuse particle is: where k ext = 1 consistently, k b is the rising parameter to determine the speed of the rise of diffuse particles, k cp is the following parameter to indicate how much the diffuse particles are affected by the neighboring fluid particles.When k cp = 1, diffuse particles move completely with the fluid particles.The position of the foam particle is updated as follows: Due to the different states of spray, foam, and bubble, the three parameters in the velocity formula have different values.As shown in Figure 6, the diffuse particles are divided into three categories according to the number of neighbor fluid particles.When the number of fluid neighbors is between 6 and 20, the diffuse particle is considered as foam.Diffuse particles with too few fluid neighbors are classified as spray, while others are considered as bubbles.
Spray is the free-flying particle that breaks away from the fluid during exercise.Spray particles move along a parabola by gravity, and their trajectories are unaffected by fluid particles (Figure 7).The velocity of the spray particle is only affected by the acceleration of the external force  Unlike spray, foam is on the surface of the water.The buoyancy and the gravity of foam particles are balanced, and there is no need to account for buoyancy and gravity changes.Therefore, k ext = k b = 0.However, the foam contacts the free surface of the fluid, and its velocity is affected by the fluid neighbors (Figure 8).The velocity of a foam particles is update as: The movement of bubble is more complicated.Bubble particles are not only affected by the neighbor fluid particles, but also lifted by the buoyancy (Figure 9).Therefore, the velocity update formula of bubble particles is expressed as: In this paper, the velocity field is used to deal with the coupling between the fluid and the diffuse materials.It avoids the searching for the neighbor of diffuse particles and greatly reduces the consumption of computing resources.Diffuse particles are classified into three categories when dealing with the interaction of diffuse materials with fluid.Spray, foam and bubble have different properties, and their interactions with the fluid is significantly different.The coupling parameters are set to different values according to the degree of coupling with the fluid, and the physical properties of different diffuse materials have been fully described.

Termination of Diffuse Material
In nature, diffuse material only has a short uncertain lifetime.Therefore, diffuse particles are defined as dynamically surviving particles in this paper.In order to describe the random disappearance of diffuse materials, each diffuse particle is assigned a random remaining lifetime t remain during initialization.Within each time step, t remain will decrease ∆t.The diffuse particles are deleted dynamically when the remaining lifetime runs out.
It is not hard to find out that foams disappear more slowly where they accumulate.In addition, bubble particles live longer than foam particles.Thus, the parameter k t is used to denote the weight of the reduction in each time step.The value of k t depends on the type of diffuse materials.For the foam particles, k t = 1 consistently.Experimental results will be better when the parameters of the spray and bubble are less than 1.The remaining lifetime for a diffuse particle is:

Implementation
In the SPH scheme, we employ the cubic spline kernel [5] as the smooth kernel in interpolations.We use compact hashing proposed in [44] to find neighbor particles.For the adaptive time-stepping schemes, we use the method explained in [45].We employed the fluid-solid coupling method proposed in [46] and Bullet [47] for simulating rigid bodies.In following experiments, surface is reconstructed by anisotropic method proposed in [48].The fluid is reconstructed into grid file according to the position of particles, and it is rendered as transparent liquid using the software named Blender [49].Foam particles are assigned to a number of tiny spheres and rendered as tiny particles according to their position.The application of our symmetric particle-based simulation scheme towards large scale diffuse fluids is presented in Algorithm 1.

Experimental Section and Results
The following experiments were performed in order to verify the effectiveness of the algorithm.The platform for these experiments is Intel (R) Xeon (R) CPU E5-2637W v2 @ 3.50 GHz 80.0 GB memory (Intel, Santa Clara, CA, USA), 64-bit Windows operating system, rendering platform for Intel (R) Xeon (R) CPU E5-2687W v4 @ 3.00 GHz 72.0 GB memory (Intel, Santa Clara, CA, USA), 64-bit Windows operating system.

Sea Water Hit Static Rigid Body
The experimental scenes were based on a dam burst collides static rigid body.Rigid bodies of different shapes were selected to validate different aspects of our algorithm.In these experiments, the fluid is rendered as dark blue transparent liquid.Diffuse particles are rendered as tiny white spheres.Rigid bodies are rendered as white solids.In this subsection, we compare the visual effect of no-diffuse fluid, our method and adding foam at the wave crest.We also verified the influence of three coupling parameters on the experimental results.
The result of Figure 1 (page 2) shows the effects of diffuse materials on the simulation results.Seawater flows from right to left and hits the ball.In non-diffuse-material results, the boundary of the fluid is not clear, and the flow tendency is hard to determine.With the addition of diffuse materials, the curved water ring appearing around the sphere can be observed.As the experimental results show, the rendering result with diffuse materials can clearly capture details such as direction and shape of the fluid flow, and the simulation results are more realistic compared with the other method.
Figure 10 shows the comparison between the diffuse materials generated at location with larger entropy and the diffuse materials generated at wave crests.The water was set on the right side of the scene, colliding with the rigid body and scene boundaries in turn and bouncing to form waves. Obviously, the amount of diffuse materials generated by this method is more than that in previous methods, and the initial positions of diffuse particles are more accurate.The parameters in the experiment are shown in Table 1.
Figure 11 verifies the effects of three coupling parameters.In the experimental results without considering the coupling parameters, the distribution of diffuse particles is scattered.The calculation of the diffuse particle velocity is inaccurate, and the resulting position has obvious deviations.In the experimental results obtained, diffuse materials adhere to the fluid boundary.The trend of the flow and the boundary shape of fluid can be observed clearly.The phenomenon of fluid interacting with movable rigid bodies and creating diffuse materials is often observed in nature.Two movable rigid cylinders were employed to impact with a dam break to demonstrate that our model can simulate this effect.Water flows from below to above, hitting two cylinders and interacting with them to form splash and droplets.The method described in [46] is used to deal with the fluid-solid coupling, so the rigid body needs to be sampled.The sampling result is shown in Figure 12.In this experiment, the fluid is rendered as light blue transparent liquid.Diffuse particles are rendered as tiny translucent spheres.Figure 12 shows the experimental results of this scene.Fluid and rigid bodies collide to produce clusters of spray, foam and bubbles.As time move forward, diffuse particles splash around and eventually disappear.In this scenario, we found the following thresholds which produce expected results:

Diffuse Generated by Moving Body
Our diffuse module is a post-process step applied to particle-based fluid.Diffuse materials are divided into spray, foam and bubbles.In this experiment, nine rigid bodies of different shapes fall into the water under the action of gravity and produce a large number of diffuse particles.Different diffuse materials are marked with different colors in the rendering result (Figure 13).We can clearly observe the difference between them.The visual effect of the experimental results can be controlled by limiting the maximum number of diffuse particles.The number of diffuse particles as a function of time is shown in Figure 14, with a max diffuse particle number of 5.5 × 10 5 .In life, bubbles will be buoyant and reach the surface of the water, and the droplets will fall back to the surface under the influence of gravity.Both of the two diffuse materials eventually become floating foam.Our experimental results show that all diffuse materials will become floating foam and eventually disappear over time.Therefore, our model is in line with objective laws.In this scenario, we found the following thresholds which produce expected results: k  3.   We establish a one-way coupling model to handle the coupling process between fluids and diffusible materials.The interaction between the two is achieved through the velocity field, where the movement of the fluid is not affected by the diffuse materials.In two-way coupling, a cross neighbor finding is performed between the fluid particles and the diffuse particles.This process is very time consuming.In our model, only the fluid neighbors of the fluid and the fluid neighbors of the diffuse materials need to be searched, which greatly saves the simulation time (Figure 15).In order to verify the efficiency of our model, we compared the computational time of fluid simulation without diffuse materials, our method and the two-way coupling algorithm (Figure 16).In the experiment, the time of the two-way coupling algorithm is obtained by calculating the time of cross-neighbor finding.When the diffusible particles have not yet been generated, the calculation time of the three methods is generally the same.As the number of diffusible particles decreases, the time required for biphase coupling gradually decreases.Experimental results show that our method can handle fluid coupling with diffusible materials more efficiently.The computational overhead of our method is negligible.Comparison of the processing time of fluid simulation (yellow), our method (blue) and two-phase coupling method (green).When dealing with the coupling of fluids with diffusible materials, the computational overhead of our method is negligible.The calculation time required for the two-phase coupling algorithm is positively related to the number of diffusible particles.

Conclusions
In this paper, a symmetric particle-based simulation scheme is proposed for diffuse materials.The entropy of flow field is used to estimate the number of new generated diffuse particles.The velocity differences and kinetic energy are considered to determine the entropy of fluid particles.Comparing to original methods based on experience, our model is more in line with the laws of physics.The coupling between diffuse materials and fluid is achieved by the velocity field.This post-processing step avoid the time-consuming process of finding neighbors and improve the efficiency.Equations with different weights are used to describe the degree of coupling between diffuse materials and fluid, which improved the accuracy of the position of diffuse particles and improved the quality of simulation significantly.It is found that the diffuse material can be well displayed in all cases, and the proposed SPH algorithm is stable and effective.
We further believe when more details like the temperature are added to model, one could begin to simulate phenomena of air-liquid conversion, e.g.water boils at high temperatures.

Figure 1 .
Figure 1.A dam burst hit a ball from the right side.The 60th frame of the experiment.Render result without diffuse materials (left).Our method (right).When the dam impacts the ball, curved water will be formed around the ball.Compared with the non-diffuse-material results, the rendering result with diffuse materials can capture the direction, shape and other details of fluid more clearly, and the simulation result is more realistic.

Figure 2 .
Figure 2. Schematics of two methods for discretizing the flow field.Grid-based Euler method (left).Particle-based Lagrange method (right).

Figure 3 .
Figure 3. Schematic for the symmetric particle-based method, Smoothed Particle Hydrodynamics (SPH).Fluid is considered to be composed of movable particles that carry physical attributions.The physical attributions of a particle are weighted in summation by its neighbors.

Figure 4 .
Figure 4.The effect of vij • xij .When two fluid particles move towards each other, vij • xij is less than zero (left).When the particle moves backwards, vij • xij is greater than zero (right).

Figure 5 .
Figure 5. Diffuse particle emitter.Diffuse particles are generated by the qualified fluid particle in the cylinder with radius r and height ∆x f .The initial result of diffuse particles is determined by the velocity and position of the fluid particles.

Figure 6 .
Figure 6.The different between spray, foam and bubble.

Figure 7 .
Figure 7.The trajectory of a spray particle.

Figure 8 .
Figure 8. Interaction between a foam particle and fluid particles.

Figure 9 .
Figure 9. Interaction between a bubble particle and fluid particles.

Figure 10 .
Figure10.Comparisons of the experimental results of our method (upper row) and foam at the wave crest (lower row).A dam burst collides a static cylinder.From left side to the right side in order of the 30th, 42th, 60th frames.

Figure 11 .
Figure 11.Comparisons of the experimental results without considering coupling parameters (left) and our method (right).The 80th frame of a dam burst collides a static cylinder.Considering the coupling parameters, the diffuse material is more in line with the fluid boundary, the experimental effect is more realistic.
k b = 5.6, k cp = 0.5, 33.The emission range of diffuse material emitters is magnified twenty times, k h = k r = 20.The parameters in the experiment are shown in Table2.

Figure 12 .
Figure 12.A dam burst collide two moving cylinders.The 31th frame of this experiment (rifht).The grid file (left) and the sample file (middle) of the solid cylinder.
33.The emission range of diffuse material emitters is magnified twenty times, k h = k r = 8.5.The parameters in the experiment are shown in Table

Figure 13 .
Figure 13.Add diffuse materials to the fluid simulation by post-processing steps.The left one shows the fluid simulation and its rendering results.The right one shows the experimental results with the addition of diffuse materials.Different diffuse materials are marked with different colors.

Figure 14 .
Figure 14.The number of diffuse particles over time.The max diffuse particle number is set to 5.5 × 10 5 .As time passes, all diffuse materials will become floating foam and eventually disappear.

Figure 15 .
Figure 15.Comparison of the neighbor finding process in two-phase coupling and our method.

Figure 16 .
Figure 16.Comparison of the processing time of fluid simulation (yellow), our method (blue) and two-phase coupling method (green).When dealing with the coupling of fluids with diffusible materials, the computational overhead of our method is negligible.The calculation time required for the two-phase coupling algorithm is positively related to the number of diffusible particles.

Table 1 .
The setting and statistics of a dam burst colliding static cylinder simulation.

Table 2 .
The setting and statistics of a dam burst colliding movable cylinder simulation.

Table 3 .
The setting and statistics of nine rigid bodies fall into the water.