https://orcid.org/0000-0002-9702-5612
Abstract
When a traditional Weakly-Compressible Smoothed Particle Hydrodynamics (WCSPH) model is used to simulate free surface flow with a large Reynolds number, an unstable numerical calculation due to high random pressure oscillations will result, while an accurate pressure field is of vital significance for simulating violent fluid-structure interactions. Riemann-based SPH and Delta-SPH are widely used to solve this problem. In this paper, to enhance computational efficiency, the SPH method is implemented on a General Processing Unit (GPU) platform using Compute Unified Device Architecture (CUDA). Parallelized SPH programs including the standard SPH method, Riemann-based SPH and Delta-SPH are verified by a dam break model with large Reynolds number and violent deformation of free surfaces. The results show that all SPH methods can vividly reflect the whole process of splashing, rolling and backward jet flow; both the Riemann-based SPH and the Delta-SPH methods are effective in alleviating the problem of inhomogeneous pressure distribution in the simulation process; the Riemann-based SPH method has better stability even with relatively large particle spacing, and it has higher accuracy in simulating impact pressure. When the number of particles reaches 100,000, compared with a single-thread Central Processing Unit (CPU) implementation, the speedups obtained with NVIDIA Titan V with high computing cores and Quadro K2200 with low computing cores are thousands and hundreds, respectively.
HIGHLIGHTS
For SPH methods parallelized on a GPU, the speedups obtained with high and low computing cores are thousands and hundreds compared with the CPU implementation.
All the traditional SPH, Delta-SPH and Riemann-based SPH can vividly reflect the whole process of dam break.
Riemann-based SPH has higher accuracy in simulating impact pressure and better stability even with relatively large particle spacing.
INTRODUCTION
Violent wave impact pressure is crucial to the design and operation safety of marine structures (Peregrine 2003; Greco et al. 2007) and hydraulic structures (Hou et al. 2014b; Tay & Thorpe 2014), and a dam break model is often used to describe its dynamic characteristics (Marrone et al. 2011). Due to the complexity of violent dynamics and large fluid deformations of free surfaces, physical model tests often struggle to reflect the interaction of multi-physical effects (Marrone et al. 2011), while numerical simulation can provide detailed flow characteristics without test limitation, which has extremely important today.
As a meshless Lagrangian method, Smoothed Particle Hydrodynamics (SPH) is essentially a method of interpolating function information and its derivatives through discrete point information (Liu & Liu 2003). The SPH method was first proposed to study three-dimensional celestial motion (Gingold & Monaghan 1977; Lucy 1977). Subsequently, it was extended to the field of fluid dynamics (Monaghan 1994; Monaghan & Kocharyan 1995; Liu & Liu 2010). In the simulation of large deformation hydrodynamic problems, it can avoid the defects due to grid distortion in traditional grid methods. Therefore, it is widely used in many problems such as free surface flow (Hu & Adams 2007; Adami et al. 2012; Hou et al. 2014a) and explosion phenomenon (Matsumoto et al. 2005).
In the Weakly-Compressible SPH (WCSPH) method, the artificial equation of state allowing weak compressibility is used to associate fluid pressure and density (Monaghan 1994; Morris et al. 1997). When simulating free surface flows with a large Reynolds number and violent fluid-structure interaction, the standard WCSPH method may produce large interfacial oscillation and false pressure, resulting in the instability of the numerical calculation (Lee et al. 2008; Antuono et al. 2012). To alleviate this problem, an artificial viscosity term was introduced in the momentum equation to dampen the pressure oscillations (Monaghan & Gingold 1983), since it would produce excessive dissipation which cannot reflect the real physical flow characteristics (Zhang et al. 2017). Subsequently, the Riemann-based SPH (Vila 1999; Moussa 2006) and Delta-SPH (Ferrari et al. 2009) methods were proposed to improve the large interfacial oscillation and false pressure, and many scholars have subsequently worked on improvements and applications for these methods (Roubtsova & Kahawita 2006; Molteni & Colagrossi 2009; Marrone et al. 2011; Zheng & Chen 2019; Meng et al. 2020).
In this paper, an efficient meshless SPH model based on Compute Unified Device Architecture (CUDA) parallel computing technology is developed. Riemann-based SPH and Delta-SPH are used to alleviate the pressure oscillation and to improve the calculation stability of violent impact flows. Taking a dam break with high Reynolds number and violent fluid-structure interaction as an example, the simulating stability and accuracy are compared and analyzed.
NUMERICAL METHOD
, 
, p represent the fluid density, velocity and pressure, respectively, 
is the viscosity term, and 
represents gravitational acceleration.SPH discretization of governing equations
indicates the density, velocity and pressure; mj is the mass of the particle j, xi is the position of the particle i, and xj is the position of the neighboring particle j; 
is the kernel function, and h is the smoothing length. When h tends to zero, the kernel becomes the Delta function. The integral of the kernel on its support is equal to 1; 
can simply be referred to as 
. In this paper, the used kernel function is a cubic spline function, which can better meet the accuracy requirements of the calculation (Yang & Liu 2012). Its expression is as follows:
is equal to 
, 
and 
in 1D, 2D and 3D, respectively.
can be expressed as: 
, where the discrete form of laminar viscosity is expressed as Equation (8), and the discrete form of turbulent viscosity is expressed as Equation (9) (Liu & Liu 2003):
represents water molecular viscosity, taken as 10-6m2/s, 
is the distance between the i- and j-th particles, coefficient η is 0.1h to avoid a zero denominator, and 
represents the Sub-Particle Scale (SPS) stress tensor of fluid particles. Dalrymple & Rogers (2006) introduced SPS into WCSPH, which can be obtained by the following formula:
, 
represents the strain tensor of fluid particles, k is the SPS turbulent kinetic energy, 
is a diagonal matrix with module equal to 1, 
is the Smagorinsky constant, usually taking as 0.12, CI = 0.0066; 
is the particle spacing.
Equation of state
is the reference pressure, 
is the reference density of water, γ = 7, and c is the numerical sound velocity. The compressibility is controlled by numerical sound velocity. In order to ensure that the density changes within 1% to meet the assumption of weak compressibility, a numerical sound velocity is usually selected at least 10 times larger than the estimated maximum fluid velocity, i.e., 
.In the WCSPH method, the fluid pressure is related to the density through the equation of state, allowing smaller compressibility. When the flow is inviscid or with a high Reynolds number, the standard WCSPH formula will produce false pressure oscillation, resulting in instability of the numerical calculation. Many researchers have adopted density filtering and a Riemann solver to alleviate this problem.
Delta-SPH method
is introduced into the continuous equation as a dissipative term to smooth the density oscillation in the calculation, so as to alleviate the pressure oscillation. The density dissipation term converges to zero as the particle size decreases. The essence of the Delta-SPH is the density filtering method, and researchers (Marrone et al. 2011; Zheng & Chen 2019) found that the Delta-SPH method can significantly reduce the numerical oscillation and effectively solve the numerical dissipation problem in the SPH calculation. After introducing the dissipation term into the continuity equation, it can be expressed as follows:
, 
represents the Delta-SPH coefficient, usually taken as 0.1.Riemann method
, subscripts L and R denote left and right states, respectively.
The construction and solution of a Riemann solver. (a) Construction of a Riemann problem along the interaction line of particles i and j. (b) The solution to the Riemann problem.
The construction and solution of a Riemann solver. (a) Construction of a Riemann problem along the interaction line of particles i and j. (b) The solution to the Riemann problem.
and pressure 
can be approximated as following:
and 
.
and 
, instead of average velocity in the continuity equation and average pressure in the momentum equation. The continuity equation and momentum equation of the standard SPH method are amended as follows (Zhang et al. 2017; Rezavand et al. 2019):
and 
. Replacing 
and 
in the above equations by 
and 
, respectively, gives
.Boundary conditions
Boundary treatment is a difficult problem in the SPH method. Effective boundary conditions can keep density and pressure stable when the fluid particles are close to the boundary, without producing particle leakage. When the SPH method calculates the particle characteristics near a wall, it needs to satisfy that there are enough and uniformly distributed particles in the support domain. Therefore, when dealing with the boundary, particles with physical quantities need to be arranged outside the wall.
represents boundary particle point pressure, 
represents fluid particle pressure, ΔV represents the particle volume. The density of virtual particles is obtained by the equation of state (Equation (12)) according to the pressure.Time marching scheme
GPU acceleration
When calculating the physical qualities of each particle, the SPH method needs to search all adjacent particles in its support domain to obtain the characteristic information of each particle by weighted average. When the number of particles is tens of thousands, the simulation speed of a traditional Central Processing Unit (CPU) using serial mode is very slow. When the number of particles reaches hundreds of thousands or even millions, then the use of a CPU parallel mode is also limited by the number of CPU cores, which makes it difficult to meet the calculation requirements. In the SPH calculation, different particles adopt the same steps, and this large-scale high synchronization logicless computing mode is very suitable for parallel computing. The neighbor search and particle summation are regarded as the key procedures of an SPH-based fluid simulator. Both the above processes are very time-consuming because the neighbor search is performed several times at each time step and the particle summation needs intensive force computations using Equations (5)–(11), which are over all the neighboring particles located in the support domain of the particle of interest. To focus on the above issues, we implemented all the main simulation steps on a modern GPU so as to fully parallelize the SPH method (Hou et al. 2017). The processor of the high-performance server used in the numerical calculation is an Intel Xeon [email protected] GHz, and the graphics cards are a NVIDIA Titan V and Quadro K2200. In this study, the acceleration effect of the Graphics Processing Unit (GPU) is analyzed by taking a dam break model as an example.
The performance results are presented in Figure 2, which exhibits the speedup achieved with the efficient GPU implementations in comparison to the single-threaded CPU. Compared with the single-threaded CPU implementation, the speedups calculated by the Titan V and K2200 increase with the increase of the particle numbers, and when the number of particles is over 10,000, the speedup of the K2200 increases slowly, while the speedup of the Titan V still increases with a high speed. For example, when the number of particles reaches 100,000, the speedups calculated with the Titan V and K2200 compared with the single-threaded CPU implementation are about 2500 and 550, respectively. The total achieved speedup on Titan V is about five times that of K2200 when the SPH method is implemented by the GPU code. In other words, the speedup calculated by the advanced GPUs with higher computing cores is obviously greater than that of the CUDA-enabled GPUs with low computing cores.
Speedup for different numbers of particles (N) achieved by the GPU implementation on NVIDIA Titan V and Quadro K2200 in.comparison to the serial CPU code.
Speedup for different numbers of particles (N) achieved by the GPU implementation on NVIDIA Titan V and Quadro K2200 in.comparison to the serial CPU code.
MODEL VERIFICATION AND COMPARISON
= 1.0 × 10−15 is added to avoid a zero denominator.
Arrangement of the dam break model.
Water flow pattern
Figures 4(a)–4(c) show the flow patterns and pressure distributions of the dam break evolution for five transient processes at t(g/H)0.5 = 2.11, 3.31, 4.93, 5.87, and 6.65 for the three SPH methods, respectively. The dimensionless time t* is taken as t(g/H)0.5, and Figure 4(d) shows the flow patterns of the dam break flow simulated by Rezavand et al. (2019). They vividly demonstrate the main flow characteristics of dam-break evolution. It can be seen from Figure 4 that the flow patterns of all three SPH methods are in general agreement with those simulated by Rezavand et al. (2019). At t(g/H)0.5 = 2.11, the water block collapsed rapidly under the action of gravity, and the water-front approached the right wall quickly along the horizontal solid wall. At t(g/H)0.5 = 3.31, the water-front encounters the vertical wall on the right side and collides violently, and then the water-front climbs upward along the vertical wall. At t(g/H)0.5 = 4.93, the climbing water-front began to turn left. When t(g/H)0.5 = 5.87, the rolling water-front impacts the downstream inflow water and forms a closed cavity. When t(g/H)0.5 = 6.65, the water-front collides with the water surface again, which forms a backward jet. Therefore, the SPH method can vividly reflect the splashing, tumbling and fusion process of the water-front during the dam break. Comparing the three SPH methods, it can be seen that the pressure distributions calculated by the Riemann-based SPH method and the Delta-SPH method are more stable and smoother, while the pressure distribution of the standard SPH method is disordered and a large number of negative pressures appear, which can lead to numerical instability and even terminate the simulation process. The Riemann-based SPH method and the Delta-SPH method can improve the problem of inhomogeneous pressure distribution significantly.
Flow pattern of dam break at different times. (a) Standard SPH method. (b) Delta-SPH method. (c) Riemann SPH method. (d) Simulated in Rezavand et al. (2019).
Flow pattern of dam break at different times. (a) Standard SPH method. (b) Delta-SPH method. (c) Riemann SPH method. (d) Simulated in Rezavand et al. (2019).
Impact pressure
Figure 5 shows the time histories of pressure at the measurement point P for the three SPH models when H/dx is 80. As can be seen at t(g/H)0.5 = 2.41, due to the nappe impact on the right side of the wall, the pressure at the measurement point suddenly increased and reached a first peak; the water then climbs upward along the vertical wall, at which the measured pressure basically remained unchanged; at t(g/H)0.5 = 6.41, the nappe turned over to the left and impacted on the downstream water to form a closed cavity, when the pressure reached its second peak value. It can also be seen that when using the standard SPH model, the high-frequency oscillation of the pressure at the measurement point is more serious, while the pressure oscillation problem is significantly reduced when the Riemann method and the Delta-SPH method are adopted. Therefore, both the Riemann SPH method and Delta-SPH method can effectively alleviate the pressure oscillation problem in the simulation process and improve the stability of the numerical calculation.
The time histories of pressure at the measurement point for H/dx= 80.
To study the convergence behavior, three different particle spacings (H/dx= 40, H/dx= 80 and H/dx= 100) were adopted. Figure 6 compares the calculated time history of the pressure with the experimental values. It can be seen that for the standard SPH model, the high-frequency pressure oscillation of the measuring point is serious, and with the decrease of particle spacing, the pressure oscillation is alleviated. Delta-SPH can effectively alleviate the pressure oscillation, but it is sensitive to the particle spacing. When particle spacing is large, although the pressure oscillation is significantly attenuated compared with the standard SPH model, the pressure oscillation cannot be ignored, which affects the stability of the numerical calculation. Besides, the peak pressure is significantly larger than the experimental value for the Delta-SPH method. When the Riemann method is adopted, the pressure time-history curves calculated under three particle spacings are relatively consistent, and the peak pressure at the measuring point is also consistent with the experimental value. Therefore, the Riemann solver based on switching function can more effectively alleviate the oscillation problem caused by violent fluid-structure impact flows, and relatively large particle spacing can generate more accurate and stable calculation results, which can effectively reduce the computational load.
Comparison of the time variations of dimensionless pressure for three SPH models under different particle spacing with experimental data (Buchner 2002). (a) Standard SPH method. (b) Riemann SPH method. (c) Delta-SPH method.
Comparison of the time variations of dimensionless pressure for three SPH models under different particle spacing with experimental data (Buchner 2002). (a) Standard SPH method. (b) Riemann SPH method. (c) Delta-SPH method.
Figure 7 gives the pressure time history obtained by the Riemann SPH method and Delta-SPH method together with that of Rezavand et al. (2019), with the particle spacing of H/dx = 80. It can be seen that the peak pressure using the Riemann-based method is closer to the test value, followed by the Delta-SPH method. Compared with the calculation results in this paper, the pressure oscillation of the SPH model used by Rezavand et al. (2019) is more significant due to the consideration of air influence, and the peak pressure is also significantly larger than the other results. Compared with the experimental measurement, the second pressure peak obtained by numerical calculation is later. The main reason is that the friction force on the wall is not considered in the numerical simulation. In addition, the simulation is two-dimensional, which is also different from the actual three-dimensional test.
Comparison of the time variations of dimensionless pressure for different methods with experimental data (Buchner 2002).
Comparison of the time variations of dimensionless pressure for different methods with experimental data (Buchner 2002).
Water-front motion
Figure 8 shows the time evolution of the water-front for the three SPH models under different particle spacing. It can be seen that the calculated results of the three SPH models converge as the particle spacing decreases. When t(g/H)0.5 = 1.7, the computational results agree well with shallow water theory (Ritter 1982). However, as observed by Colagrossi & Landrini (2003), the movement of the water-front in the calculation is relatively faster than the experimental result, which is mainly influenced by the turbulent wall boundary layer, wall roughness and experimental errors. Moreover, the shallow water assumption does not hold at the initial moments.
Time evolution of the water-front for the three SPH models under different particle spacings. (a) Standard SPH. (b) Riemann-based SPH. (c) Delta-SPH.
Time evolution of the water-front for the three SPH models under different particle spacings. (a) Standard SPH. (b) Riemann-based SPH. (c) Delta-SPH.
CONCLUSIONS
In this paper, based on CUDA parallel computing technology, an efficient meshless SPH model has been developed. A dam break problem with high Reynolds number was used to verify its effectiveness. The main conclusions are as following:
- (1)
Compared with a single-threaded CPU implementation, the speedups calculated by NVIDIA Titan V and Quadro K2200 increase with a increase of particle numbers, and when the number of particles is over ten thousand, the speedup of K2200 increases slowly, while the speedup of Titan V still increases with a high speed. When the number of particles reaches one hundred thousand, compared with the single-threaded CPU implementation, the speedups obtained with NVIDIA Titan V with high computing cores and Quadro K2200 with low computing cores are thousands and hundreds, respectively. The total achieved speedup on Titan V is about five times that of K2200.
- (2)
The numerical simulation results can vividly reflect the whole process of the splashing, tumbling and fusion process of the water-front. In order to reduce the false pressure oscillation in the standard WCSPH model and improve the instability of numerical calculation, this study adopts THE Riemann-based SPH method and Delta-SPH method, and the accuracy of the two SPH methods was verified through a dam break problem. The simulation results show that the pressure field distribution is stable and smooth in both the Riemann-based SPH method and Delta-SPH method.
- (3)
Through a convergence study, it was found that when using the standard WCSPH model, the pressure oscillation is alleviated with the decrease of particle spacing. The Delta-SPH method is sensitive to particle spacing. When the particle spacing is large, although the pressure oscillation attenuation is significant compared with the standard WCSPH model, the pressure oscillation is still intensive, and the peak pressure is significantly greater than the experimental value. When the Riemann-based SPH method is adopted, the pressure time-history curves obtained by three particle spacings are relatively consistent, and the peak pressure at the measuring point is consistent with the experimental test. Therefore, the Riemann SPH method based on switching function can more effectively alleviate the pressure oscillation caused by violent wave impact, and the relatively large particle spacing can obtain accurate and stable calculation results, which can effectively reduce the computational load.
ACKNOWLEDGEMENTS
This study has been partly funded by the National Key Research and Development Program of China (No. 2020YFC1807905), National Natural Science Foundation of China (No. 52079090) and Basic Research Program of Qinghai Province (No. 2022-ZJ-704).
CONFLICT OF INTEREST STATEMENT
The authors are not affiliated with or involved with any organisation or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this paper.
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.
%20cropped.png?versionId=7153)