Abstract
To achieve small-scale ocean scene simulation in a marine simulator, we present an incompressible SPH algorithm by relaxing the density-invariant condition for incompressible fluids. As there are larger density errors of the fluid particles near or on the boundary, more iteration numbers are required. Taking boundary handling to modify the density of the fluid particles, the relaxation method is used to optimize the iterative method to solve the density-invariant condition, which can reduce the iteration numbers and average density deviation and improve the accuracy. Our proposed approach can achieve incompressible SPH and solve the problem of the particle deficiency. While ensuring the stability, the approach can allow large time steps and control the density deviation below 0.01% and improve the efficiency by reducing the iteration numbers and optimizing the calculating procedure and the initial value selection.
1. Introduction
Ocean wave simulation is an important part of fluid simulation. The realistic and real-time characteristics and interactivity of the simulation directly affect the degree of realism of the sea scene in a marine simulator. At present, physically-based simulation methods are used to enhance environmental realism by means of the methods of computational fluid dynamics (CFD). The methods based on the Navier-Stokes equation can be divided into two categories: the Euler method, which is a grid-based method, and the Lagrange method, which is a particle-based method. Compared with the Euler method, the Lagrange method can effectively avoid the problem of numerical dissipation because there is no direct convection term, and it can suitably address large deformation and free surface problems. In the Lagrange method, the smoothed particle hydrodynamics (SPH) method is more realistic, stable, and efficient because it is easy to deal with the simulation of surface flow and the fluid-structure coupling of large deformation. At present, this method has become one of the most promising simulation methods. Desbrun et al. [1] first introduced the standard SPH (SSPH) to computer graphics (CG) for simulating deformable bodies. Since then, Müller et al. [2] have applied SSPH to simulate water with a free surface.
Currently, there are two main methods to simulate incompressible fluid by SPH: one is to regard incompressible fluid as a slightly compressible fluid and introduce a state equation to simulate weak compressible fluid, which is called weakly compressible SPH (WCSPH). The other way is to use a semi-implicit method to enforce incompressibility, called incompressible SPH (ISPH). The WCSPH method was proposed by Monaphan et al. [3] and Becker al. [4], which simulated the free surface flow. The equation of state was employed to obtain the fluid pressure, which is an explicit method, but smaller density fluctuations can lead to a sudden increase in pressure. It is often necessary to introduce a larger sound speed to limit the density fluctuation rate to less than 1%, which requires a smaller time step and affects the efficiency of the algorithm. A slightly larger time step causes instability.
The ISPH method is a better way to simulate incompressible fluids using the projection method [5–7], called the semi-implicit method. The fluid pressure is obtained by solving the pressure Poisson equation (PPE). This method can allow a larger time step and strictly enforce incompressibility, but it needs to solve large linear equations. The coefficient matrix is very large, resulting in high computational cost. The methods for solving the PPE are usually grid-based methods and meshless particle methods. Grid-based methods [7–10] are traditional, usually taking a conjugate gradient, a successive over-relaxation (SOR), a multiple grid, or a Voronoi diagram method to solve the PPE on the background grid.
The meshless particle method is mainly used to solve the PPE by iteration methods such as the predictor-corrector method and Jacobi-style method to realize a density-invariant condition to simulate incompressible fluid. Solenthaler and Pajarola [11] proposed the PCISPH method, which uses the pressure field of the current time step to obtain the particle position of the next time step and then uses the density estimation to correct the current pressure field. The density error of the next time step is less than the given threshold to achieve the density-invariant condition by iterative method. The implicit incompressible SPH (IISPH) method proposed by Ihmsen et al. [12] is also strictly a semi-implicit incompressible SPH method. The relaxed Jacobi method is used to solve the PPE. The divergence-free SPH (DFSPH) method proposed by Bender et al. [13] is also a predictor-corrector method for solving the PPE. However, in addition to the density-invariant condition, the velocity divergence-free condition is used as the source term of the PPE, and two pressure solvers are combined to enforce incompressibility. Subsequently, Bender et al. [14] extend the method and applies it to the simulation of viscous flow.
There are two ISPH versions: CFD version and CG version. The two ISPH versions have obviously different features and scope of applications. In CFD version, the approach is developed for solving complex hydrodynamic problems such as modeling complex free surface flows, multiphase flows, and fluid-structure interactions and have widely applied in a wide range of hydrodynamic problems (e.g., water-body interactions, water impacts, water entries and exits, and flow fast bodies) and explosion problems (e.g., underwater explosions, contact explosions, etc.) [15]. In CG version, the approach is widely applied in computer graphics (e.g., the simulation of navigation, salvage, flying, virtual naval battlefield, game, and film), which is faster and easier for modeling the incompressible fluid. Different from CFD version, they simplified the physical model to save the computational cost. Their goal is to achieve a balance between realism and real-time. Raveendran et al. [16] proposed that solving the PPE on a coarse grid enforces a divergence-free velocity field and generates pressure forces that are then passed along to the WCSPH or predictive-corrective incompressible SPH (PCISPH) solver to achieve incompressible flow. This approach can be seen as the combination of Euler and SPH approaches. Losasso et al. [17] combined the particle level set (PLS) and SPH method, which uses the PLS method to model dense liquid volumes and the SPH method to simulate small-scale phenomena, such as spay and foam. The SPH solver also needs to solve the PPE on a background grid, not only to maintain the divergence-free flow field but also to maintain a predefined target density.
In ocean engineering, the breaking wave with spay and foam happened in the interaction of the high speed ship and wave. Marrone and Colagrossi [18] employed a 2D+t SPH approach to study the breaking wave pattern of a fast ship for different forward ship speeds. The vortical structures due to plunging of the bow wave are tracked using passive markers which are initially positioned on the undisturbed free surface. Sun and Colagrossi [19] presented algorithms for the detection of Lagrangian Coherent Structures in viscous flows through the Finite-Time Lyapunov Exponents (FTLEs). FTLEs introduced to Hydrodynamics are used for capturing the bow breaking wave feature. In order to make it easier, they adopt a simplified ship hull. In computer graphics, and Losasso et. al. [17] extended the SPH solver to simulate diffuse phenomena such as mixtures of spray and air, but foam is not based on physical model and a large number of spay particles can only be rendered offline. Ihmsen et al. [20] present a postprocessing model to simulate spay, foam, and air bubbles. This method enables efficient processing and versatile handling, but they ignored the effect of diffuse material on the fluid simulation and rendered the scene with mental ray 3.9, based on the ray casting and cannot meet the requirements of real-time rendering.
Macklin and Müller [21] presented a position based fluids (PBF) method which integrates an iterative density solver into the Position Based Dynamics framework (PBD) and applies to fluid simulation. Similar to the PCISPH method, the PBF method employs the predictor-corrector method to constrain the fluid position so as to enforce incompressibility. The position change is iteratively computed for each particle rather than the pressure so that PBF method allows larger time step, achieves better stability, and is suitable for real-time rendering. Because the particle position updating is based on geometry, PBF method is not physical simulation in its strict sense. In particular, collision detection and response of PBD are adopted for boundary handling and there is a larger density error at the interface between fluid particles and the ground or wall, which requires more iteration numbers, thus affecting the efficiency of the solver. In this work, the boundary particle method is introduced into the PBF method to correct the density error of fluid particles near or on the boundary. For density-invariant conditions, the relaxation method, optimal initial value, and iterative termination condition are adopted to accelerate convergence, reduce the number of iterations, improve the efficiency of predictor-corrector iteration, and increase the time step for improving the overall performance.
2. Position-Based Fluids Method
The incompressibility of fluid requires that the fluid particle densities remain constant. The idea of the PBF method is to modify the particle position when the density error of fluid particles exists, and the particle position is iteratively corrected so that the predicted density error is lower than a predefined value. The density constraint of a fluid particle is defined aswhere , , and denote the position, density, and rest density of a fluid particle , respectively.
The fluid particle density , can be solved by density summation formula in SPH method as follows:where and denote the mass and position of a neighbor particle , respectively, is a kernel function with support radius , and .
In order to enforce incompressibility, the fluid density must be kept constant which is realized by correcting the fluid position. The position changes should satisfy the following constraint:From (3) using the first-order Taylor expansion, we can get the following:Considering that the constraints do not affect the motion state, the direction of partial derivative of with respect to is perpendicular to the original direction of motion. The direction of is the same as the direction of . Importing a Lagrange multiplier , is expressed as follows:Insert (5) into (4), and we obtainSo (5) changes intoThe position change and factor of a fluid particle are computed asDerivative is given byThe denominator in (9) causes instability when the neighbors of the support domain are insufficient so we use a parameter to avoid singularity and set is . The factor is now
To speed up the updating of the particle position changes, the position change is updated by the factor of the neighbors in the support domain of particle asThe whole algorithm is shown in Algorithm 1.
| ||||||||||||||||||||||||||||||||||||
3. Incompressible SPH Method
3.1. Relaxation Method
With the increase of the accuracy requirement of fluid incompressibility and the number of fluid particles, the computational cost time for solving the PPE is dramatically increased. To solve the PPE quickly and accurately and save the computational cost time, we adopt relaxation method to accelerate convergence. So the predicted position is adjusted bywhere is the relaxation factor. If is chosen properly, the convergence can be greatly accelerated. We use different relaxation factors to test the maximum number of iterations when the maximum volume compression is set to in 0.1% or 0.01%, respectively. The maximum number of iterations in different relaxation factors is shown in Table 1.
When is set to 0, our method is equivalent to the PBF method and the maximum iteration number is 47. When is greater than 0.5, although the maximum number of iterations is still decreasing, the simulation tends to be unstable. When , the maximum number of iterations decreases with the increase of . Therefore the relaxation factor is set to 0.5.
3.2. Boundary Handling
Recently, the boundary handling methods of SPH fluid simulation mainly include penalty method, ghost particle method, and boundary particle method. In penalty method, distance-based penalty forces, such as Lennard-Jones forces [3] and normal repulsion forces [22], are exerted for preventing the fluid particles at the fluid-solid interface from penetrating the solid boundary. But a large penalty force is needed, which results in large pressure fluctuation of the fluid particles on the boundary and momentum conservation is not guaranteed. Becker et al. [23] employed direct force to correct the position and velocity of the fluid particles for two-way fluid-rigid coupling and non-penetration. Different slip conditions are realized by including a nonsymmetric friction force. But penalty force and direct force method have the problem of the fluid particles adhering to the boundary, which is caused by the insufficient number of the fluid particles at the fluid-solid interface.
In the ghost particle method [24, 25], ghost particles are dynamically generated, which is symmetrical to the fluid particle outside the boundary. A ghost particle gets the same density, pressure, and velocity at opposite direction of its corresponding fluid particle, which participates in the density estimation. This method was successfully employed to simulate the different slip conditions and ensure conservation of momentum, but it is difficult to deal with the complex boundary problem. Aiming at the defect of standard ghost particle method in modeling complex geometries, Marrone et al. [18] develop the fixed ghost particle method. Differently from the standard ghost particle method where at each time step any particle nearby the solid boundary is mirrored into a ghost with respect to that boundary; the method fixed the mirror particles to the outside of the fluid at the boundary. The interpolation of these particles is carried out by using the method of higher order difference to compute the physical quantities of these particles. Finally, the corresponding fixed ghost particles are assigned by the mirror rule. This method solves the problem of particle defects and overlap caused by standard ghost particle method and is suitable for dealing with complex geometries and gives quite good and accurate results. Subsequently, Sun et al. [26] introduced it into the numerical simulation of the self-propulsive motion of a fishlike swimming foil. However, this method is complicated in computing the physical quantities attributed to each ghost particle; an interpolation point is associated with it. This interpolation point is obtained by mirroring the position of the fixed ghost particle into fluid domain, which results in high computational cost and requires high accuracy for interpolation methods.
In boundary particle method [12, 27, 28], layer, two-layer, or multi-layer boundary particles are sampled at the solid boundary. The number of boundary particle layer is determined by the support domain. The boundary particles have the same rest density of the fluid particles. The velocity of the boundary particles is zero for the static boundary, while the velocity is the same as the velocity of the dynamic boundary. The density and pressure of the boundary particles also participate in the interpolation of the density, pressure, and internal force of the fluid particles.
Therefore, we employ the boundary particle method to solve the problem of the particle deficiency, sample the solid triangles by Poisson disk sampling method [25], and generate single layer boundary particles which participate in the density estimation of the fluid particles. In order to reduce the error of density estimation caused by unreasonable mass or nonuniform distribution of boundary particles in single layer boundary. The mass of boundary particles is adaptively changed according to the distribution of boundary particles. Assuming that the density of a boundary particle b is equal to the density of a fluid particle , the mass of a boundary particle is computed aswhere denotes the volume of a boundary particle. Its value is , and the number density of a boundary particle . After the correction of the continuity equation, the predicted density of a fluid particle is computed as
In this way, the density fluctuation caused by under sampling or over sampling can be corrected to improve stability. Similarly, after the correction of the momentum equation, the position change of a fluid particle is computed as
The whole algorithm of incompressible fluid simulation based on relaxing density-invariant condition is illustrated in Algorithm 2. We modify the PBF algorithm as follows: first of all, the convergence speed is accelerated by relaxation method and the relaxation factor is set to 0.5 (line (16) of Algorithm 2). Then, collision detection and response are replaced with boundary handling (line (14) of Algorithm 2). In (2), the density equation may lead to density underestimation and cause particle clustering near a free surface where kernel truncation occurs. We solved this issue simply by clamping density to rest density of a fluid particle. Considering the effect of relative velocity on density computing, (15) and (16) are used to replace (2) and (12) for computing the density and position change. Finally, if solverIterations is set too small (line (7) of Algorithm 1), the average density error in simulation may not reach the required target but the iteration terminates, which makes the density estimation inaccurate. The limit of the maximum number of iterations is replaced with the average density error less than the predefined threshold as the termination condition of iteration (line (7) of Algorithm 2), so that the average density error can satisfy the requirements of accuracy.
| ||||||||||||||||||||||||||||||||||||||||
4. Foam
The effect of broken waves, such as foam and spray, directly affects the visual realism of the whole scenario. The key to simulate foam and spray is the determination of foam generation area, the advection, and dissolution of foam and spray. The potential value of each fluid particle calculated by geometric characteristics, kinetic energy, and relative velocity is compared with a given threshold to determine the area of foam generation. The detailed small-scale effects are achieved by a postprocessing technology and modeled with unified particles which are called diffuse particles. We classify the particle into foam and spay. In our experiments, the diffuse particles with less than 8 fluid neighbors are considered as spray particle and, in all other cases, the particles are considered to be foam. Physically motivated rules are employed to generate and advect and dissipate spay and foam, but interparticle forces and the influence of spay and foam particles onto the fluid particles are neglected. The rendering of the foam and spay is employed by the screen space foam rendering [29].
4.1. Generation
The probability of a fluid particle generating diffuse particles is expressed by the control function which is defined aswhere , is the minimum and maximum threshold and is the potential for the generation of diffuse material.
Geometric Characteristics. The wave crest is unstable and easy to produce diffused particles. At the wave crest, the surface curvature is relatively high and the local surface is convex. The surface curvature can be approximated withwhere is a normalized surface normal. In order to distinguish convex from concave regions, considering the vector product between and , wave crests are identified usingwithThe condition of the identification of wave crests is adjusted as
Finally, the likelihood of a particle to be at the crest of a wave is computed aswith user-defined and .
Kinetic Energy. The kinetic energy is used as a measurement for the generation of diffuse particles. The higher the velocity of a fluid particle, the greater the energy and the easier the diffusive particles will be generated at the position. The potential to generate diffuse particles due to kinetic energy is computed aswith user-defined and .
Relative Velocity. The relative velocity is used to determine regions where air is potentially trapped since these are large for impacts and vortices. The larger the relative velocity of fluid particles is, the easier it is to generate diffuse particles. The potential to generate diffuse particles due to relative velocity is computed aswith user-defined and , where is the with the normalized relative velocity and normalized distance vector , is a radially symmetric weighting function defined as
Therefore, the number of diffuse particles generated by a fluid particle is computed as
Similar to [20], we sample the diffuse particles in a cylinder which becomes a rectangle in two dimension (see Figure 1). The position and velocity of a diffuse particle d are computed asandwhere three uniformly distributed random variables , the distance to the cylinder axis , the azimuth , the distance , and normalized vector is perpendicular to the velocity direction of the fluid particle.
4.2. Advection
Spay particles are affected by gravity and external force , and using the Euler method, the velocity and position are updated as
But foam particles are mainly influenced by the fluid particles of the supporting domain. The velocity of the foam particle d is computed by the average velocity of the fluid particles of the supporting domain which is computed aswhere . Then the position of the foam particles is updated as
4.3. Dissipation
Similar to the particle system, the diffuse particles’ lifetime is initialized with a predetermined value. In each simulation step, the time step is subtracted from the lifetime of foam particles, whereas the lifetime of spray particles is not reduced. Foam particles are finally deleted when the lifetime is smaller or equal to zero.
5. Results
To verify the validity of the algorithm proposed in this paper, we performed simulations of several small-scale ocean scenes, such as breaking dam, wave heaving, rolling, and breaking using 262k fluid particles on two Intel (R) Core (TM) i5-6400 processors with 2.2 GHz CPU, 8 g memory, and NVIDIA GeForce GTX960 graphics. The visual effect is shown in Figures 2 and 3. The boundary handling was performed using our proposed method, the neighborhood search algorithm using the space hash function [30], smoothing kernels function with the cubic B-spline function [31], viscous force with the method proposed by Schechter and Bridson [25], and vorticity confinement with the method of Macklin and Müller [21]. To render the fluid surface, we employed the screen space fluid rendering method [32].
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
6. Discussion
6.1. Convergence
To compare the convergence and overall performance of the proposed approach and PBF method, we simulated a breaking dam scenario with 4500 fluid particles and 18630 boundary particles. The particle radius r is 0.025 m and the average density change rate set to 0.01% or 0.1%, and the fixed time step is 0.005 s. The average iteration numbers of each physical time step can be found in Figure 4.
(a)
(b)
When the maximum density error is set to 0.1%, the average iteration of PBF method is 10.6 times and our approach is 6.2 times, and when the maximum density error is set to 0.01%, the average iteration number of PBF method is 29.0 times and our approach is 11.7 times, which greatly improved convergence. Solenthaler and Pajarola [11] indicated that the peak of average iteration number appears when fluid particles collide with the ground and wall. The reason for this is that there exist boundary deficiencies between the fluid particles and the interface of the ground and wall, and thus a large density error is predicted. In our proposed approach, a layer of boundary particles is placed in the boundary position and participates in the density calculation to correct the boundary deficiency, and the peak of the average iteration number appears in the position of larger surface curvature changes.
A convergence comparison of one physical time step at the time of 0.2 s is shown in Figure 5.
(a)
(b)
When the maximum density error is set to 0.01%, the maximum number of iterations for the PBF method is 47 and the maximum average density change rate is 3.91%, while the maximum number of iterations for our approach is 16, and the maximum average density change rate is 2.90%; when the maximum density error is set to 0.1%, the maximum number of iterations for the PBF method is 17 and the maximum average density change rate is 5.83%, while the maximum number of iterations for our approach is 10, and the maximum average density change rate is 5.49%. Solenthaler and Pajarola [11] showed that the maximum density change rate can be controlled within 1% and, at the time of 0.2 s, the maximum average density change rate is 6.2%. On one hand, in our approach, a boundary particle is adopted to correct the fluid particle density near or on the boundary, which greatly reduces the predicted maximum density error; on the other hand, it uses the relaxation method and optimizes the selection of the initial value to significantly decrease the maximum iteration number and improve convergence of iteration.
6.2. Performance
For predictive-corrective iteration, as long as the iteration number is sufficient, the precision of the result can be guaranteed; however, the convergence speed is sometimes slow, which causes a heavy calculation burden. Its convergence also depends on the selection of the initial value. We adopt the relaxation method and optimize the selection of the initial value with 0.5 times the value from the last time step to reduce iteration numbers and accelerate the convergence speed. We used different fixed time step to control the density change rate within the range of 0.1% or 0.01%, and the tested pressure solver and the overall time were used to illustrate the time step and convergence effect on the overall performance. The measurements are shown in Tables 2 and 3.
As shown in Tables 2 and 3, under the condition of a larger time step, our proposed approach exhibits better performance than the PBF method. As shown in Table 2, by adopting our approach, the maximum density error is controlled within 0.1%, the pressure solver obtained a speedup of 1.9, and the total time obtained a speedup of 1.6 times that of the PBF method. Table 3 shows that the maximum density error is controlled by 0.01%, the pressure solver obtained a speedup of 2.6 times, and the total time obtained a speedup of 2.3 times that of the PBF method. The improvement in the overall performance is related not only to the reduction of iterations but also to the decrease of the consuming time of each physical time step which saves data storage and operation times.
7. Conclusions and Future Work
In this paper, we have proposed an incompressible fluid simulation approach based-on relaxing density-invariant condition, realized the simulation of small-scale ocean scenes, such as breaking dam and wave heaving, rolling, and breaking, and solved the problem of the boundary particle deficiency problem. The execution efficiency is improved by reducing the number of iterations and optimizing the calculation process. Compared with the PBF method, the pressure solver and the total calculation of our approach obtained a speedup of 2.6 times and 2.3 times, respectively, when the maximum density error is controlled within 0.01%.
The density filtering using the well-known MLS interpolation [33, 34] or the Shepard kernel interpolation [35] is an more efficient way to solve the problem of the particle deficiency near the free surface. It is our next research. Based on the ISPH by relaxing the density-invariant condition method, we plan to simulate the interaction between a vessel underway and water, including spray and foam, on the bow of the vessel in a marine simulator.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was funded by the Ministry of Science and Technology of People’s Republic of China via National 863 Project No. SS2015AA010504. Funding was also provided by the China Department of Education of Fujian Province under Award no. JA15269.