Abstract
The smoothed particle hydrodynamics (SPH) method has been popularly applied in various fields, including astrodynamics, thermodynamics, aerodynamics, and hydrodynamics. Generally, a high-precision interpolation is required to calculate the particle physical attributes and their derivatives for the boundary treatment and postproceeding in the SPH simulation. However, as a result of the truncation of kernel function support domain and irregular particle distribution, the interpolation using conventional SPH interpolation experiences low accuracy for the particles near the boundary and free surface. To overcome this drawback, stable regularized moving least-squares (SRMLS) method was introduced for interpolation in SPH. The surface fitting studies were performed with a variety of polyline bases, spatial resolutions, particle distributions, kernel functions, and support domain sizes. Numerical solutions were compared with the results using moving least-squares (MLS) and three SPH methods, including CSPH, K2SPH, and KGFSPH, and it was found that SRMLS not only has nonsingular moment matrix, but also obtains high-accuracy result. Finally, the capability of the algorithm coupled with SRMLS and SPH was illustrated and assessed through several numerical tests.
1. Introduction
The smoothed particle hydrodynamics (SPH) method is a Lagrangian meshless method that was developed for the astrophysical flow simulation in 1970s [1, 2]. It has been a promising numerical technique for the complex fluid flow simulations as well as for various fields, including thermodynamics and aerodynamics in recent years [3–7]. Given the widespread using of SPH presently, a stable and high-order interpolation method is indispensable to overcome the lack of accuracy of the standard SPH approximation for the boundary treatment, field refreshing, and postproceeding, especially in the zone with the truncation of support domain.
Some researchers try to restore SPH consistency and investigate the truncation error near the free surface and boundary and significant progress has been done over the years to improve the accuracy with the SPH interpolation [8–12]. Most of their approximation results highly depend on the characteristics of algebraic matrices. For example, the corrected SPH (CPSH) [13] and its extensions K2SPH [14] and KGFSPH [15], both of them have the second-order accuracy theoretically, and the approximation results are solved with linear systems which are built up with the particles in support domain, namely, the results rely on the vicinity particle and special treatment should be applied for irreversible algebraic matrix, and the order of interpolation decreases consequently. To avoid the tensile instability problem inherently existed in SPH, a more accurate alternative to the SPH gradient formulation, Lagrangian gradient smoothing method (L-GSM), was introduced by Mao [16–18], and it could handle the incompressible flows with complicated free surfaces effectively and easily.
Moving least-squares (MLS) method is a high-order mesh-free interpolation algorithm, and the combination of MLS and SPH method has been motived by the fact that mesh-free and mesh-adaptive discretizations are always better to cope with large geometric changes of interest domain, such as free surface and large deformation of elastic structure, than the classical grid-based methods. Currently, MLS is widely used for boundary treatment [19] and field refreshing in SPH simulation. Although the moment matrix in MLS may be singular for ill-quality particle sets, several techniques have been proposed to deal with ill-posed least-squares problem such as coordinate transformation, perturbation of nodal position, and matrix triangularization [20–22].
In this paper, a stable regularized moving least-squares (SRMLS) interpolation whose moment matrix can be always nonsingular was introduced to interpolate the particle properties for the phases of boundary treatment, density re-initialization, and postproceeding in SPH simulation. Different from previous studies [23–25], to highlight the advantages of SRMLS, the properties of this method with the nonuniform particles that were generated by a SPH pretreatment method [26] are presented and discussed. The surface fitting solutions were compared with the results using moving least-square (MLS) method and three SPH methods, including CSPH, K2SPH, and KGFSPH, and it was found that SRMLS not only has nonsingular moment matrices, but also obtains high accuracy result. Then the capability of the SRMLS was also illustrated and assessed through some numerical tests with SPH.
This paper is organized as follows. The SPH and stable regularized moving least-squares (SRMLS) methods are reviewed in the second section. With the discussion of surface fitting results, the properties of SRMLS are presented in the third section. Afterwards, the coupled algorithm with SRMLS and SPH method is checked with a series of 2D simulations in the fourth section. In the first case (Taylor–Green vortex flow), the algorithm is tested with periodic boundary, and the method is checked with external and internal moving boundary in the second simulation (Lid-Dirven cavity flow), and the third case (moving square inside a rectangular box), respectively. The algorithm performs for a real application with free surface in the last case (solitary wave slamming). For all simulations, verifications are performed by comparing the computed solutions with available reference results. Finally, the conclusions are summarized in fifth section.
2. Mathematical Models
In this section, the mathematical models, SPH as well as stable regularized moving least-squares (SRMLS) methods, used in this paper will be briefly recalled. For SPH approach, the weakly compressible flow model was adopted and the pressure of fluid depends only on the density.
2.1. SPH Methodology
With the above assumption, the Lagrangian form of the Navier–Stokes equation can be written aswhere is fluid density, u is particle velocity, t is time, p is pressure, is the body force, and υ is kinematic viscosity. With the hypothesis of weak compressibility, the pressure can be computed by an explicit expression:where c0 is the artificial sound speed, it can be taken as 10 times of maximum velocity, and γ = 7 for simulation with water.
To get rid of the pressure noise, this paper adopted a variant SPH, δ-SPH method [27], which can effectively improve the stability of flow pressure with density correction. The discretization form of the δ-SPH equations are written aswhere the subscript i and j are the particle indexes, h is the smooth length, , is the distance difference vector, is called as kernel function or weight function value, is the body force associated to the generic particle I, and Vj is the volume of particle j.
The renormalized density gradient, , should be calculated firstly to get Ψij in equation (6), it is defined as
In δ-SPH equations, δ and α are the coefficients of density and velocity diffusion terms, respectively. α = 0.01 and δ = 0.1 in this paper. More details about the δ-SPH and its extensions and applications have been documented in the work of Marrone et al. [19] and Bouscasse et al. [27]. As only δ-SPH method is discussed in present work, it is referred to as SPH in latter discussion for convenience.
2.2. Stable Regularized Moving Least-Squares Method
The domain Ω is discretized with n particles, . is the target interpolation node. To approximate the distribution of a given function f in Ω, over these particles, the MLS obtains the shape function by minimizing an error function between the nodal values of the given function and the local approximation function as
The polynomial can be computed by a direct solution of a least-squares problem and then used to interpolate or estimate its derivatives:where is the unknown coefficient and is the basis function. In SRMLS, the stable basis function in 2D can be chosen as follows:(i)Linear Basis(ii)Quadratic Basis In 3D space, the basis can be chosen as follows:(iii)Linear Basis(iv)Quadratic Basis
in which .
Equation (8) can be written in matrix form aswhere
By taking the derivate of E respect to a as
Then one can obtainwhere moment matrix, A, is equal to and . In mesh-free method, it may be time consuming and not even obtain desired accuracy when A is ill-conditioned. To overcome this defect, regularization is introduced for SRMLS, namely, adding slightly disturbances in P and W:where β is the disturbance constant.
Finally, the local approximation and corresponding derivatives can be computed aswhere is the shape function and the index following comma represents spatial derivative.
The shape function is written aswhere , , . The derivatives of shape function can be found in Appendix A.
2.3. Coupling Algorithm
By coupling the two approaches, it is possible to obtain an algorithm that retains the best for each single approach when simulating flow with domains characterized by different kernel functions, support domain sizes, and spatial resolutions. Generally, we want to exploit the ability of SPH method to simulate the free surface flows with violent deformation, complicated boundary geometrics, and get rid of high-frequency pressure vibration through SRMLS interpolation.
The contribution of SRMLS in SPH method is providing high-quality interpolation results for particle’s unknown physical characteristics, instead of using conventional SPH kernel approximation, in the scheme of boundary treatment, and post-proceeding. As an example, a sketch of interpolation for the value of fictitious particle i locates in support domain (the pink zone) with the 12 vicinity particles is shown in Figure 1, SRMLS can be used in this procedure. The boundary treatment method used in this paper can be found in the work of Chen et al. [28].
SRMLS was also applied for another important scheme, the density refreshing, in the present study. To restrain the vibrations of pressure, the density field is refreshed every 200 steps for all the numerical tests. In terms of postproceeding, it is necessary to project the results to a background grid to draw the streamline and capture the free surface [29], or to sample result of arbitrary position in the interest fluid domain, SRMLS method has the ability to deal with the interpolation or sampling in such postproceeding schemes easily. The flow chart of the coupled algorithm is shown in Figure 2.
3. Surface Fitting
In this section, comprehensive studies of surface fitting with SRMLS were carried out. Different from previous work [23, 24], this paper focuses on the error analyzation of SRMLS results with the nonuniform nodes for surface fitting. Specially, the irregular nodes were generated with the SPH pretreatment method [26] in order to imitate particles in SPH.
Franke’s function [30], which is a standard test function for two-dimensional scattered data fitting since the seminal survey of Franke, is considered as test function:defined over the domain and . The shape and surface curvature of Franke’s function are shown in Figure 3.
The 2D domain was discretized with sets of three type particles (random, irregular, and regular particles). Figure 4 presents the computational domain, local snapshots, and irregular particles-distribution with dx = 4.00E − 2, in which, the value of ghost particle is defined as the analytical result of corresponding domain particle. Note that the irregular particles were generated with random particles by using SPH pretreatment method. It is obviously that the irregular particles-distribution was quasiuniform in Figure 4(c).
To analyze the interpolation results, the Root Mean Square Error (RMSE) and Relative Error (RE) between numerical and analytical results were calculated with NI = 2,500 (50 × 50) sampling points distributed uniformly in the square. The error estimators are written as
In consideration of the co-operation with SPH method, this paper is interested in the interpolation results of function value and first-order derivatives with various kernel functions, support domain sizes, spatial resolutions which are the sensitive parameters in SPH simulations. Surface fitting solutions were compared with the results of CSPH, K2SPH, KGFSPH, and MLS, and the properties of SRMLS will be discussed in the following subsections. Simultaneously, the second-order derivatives of the numerical exercises with quadratic basis were also presented to highlight the advantages of SRMLS. Cases in the surface fitting studies are summarized in Table 1.
3.1. Dependence of the Error on Disturbance Constant
It is necessary to investigate the influence of disturbance constant, β, on SRMLS results for guiding the choosing of disturbance constant value. 120 (4 × 15 × 2) cases with four disturbance constants, β = 1.0E – 1, 1.0E – 3, 1.0E – 5, and 1.0E – 7, and fifteen sets of irregular particles, dx = 4.00E – 2, 2.00E – 2, 1.33E – 2, 1.00E – 2, 8.00E – 3, 6.67E – 3, 5.71E – 3, 5.00E – 3, 4.00E – 3, 3.33E – 3, 2.50E – 3, 2.00E – 3, 1.67E – 3, 1.25E – 3, and 1.00E – 3, linear and quadratic basis were carried out for surface fitting. All the numerical tests were performed with the renormalized Gauss kernel function [31] and h = 1.0dx in this subsection.
Figure 5 shows the comparison of the RMSE for function value and first-order derivatives with different disturbance constants and bases. For the results with linear basis, it is obvious that SRMLS and MLS obtain same order of accuracy for both function value and first-order derivatives. No significant influence of disturbance constant is observed in the SRMLS result with linear basis. However, the disturbance constants affect the SRMLS results with quadratic basis clearly, especially with β equal to 1.0E – 1. The results are almost same when disturbance constant less than 1.0E – 3.
(a)
(b)
(c)
The comparison of the RMSE for second-order derivatives with different disturbance constants and the relative error contours of SRMLS with dx = 2.50E – 3 and β = 1.0E – 5 are shown in Figure 6. Furthermore, the convergence rates of function value and derivatives with different disturbance constants and bases are presented in Table 2. For all methods, the convergence rate corresponds to the slope resulting from least-squares fitting of the RMSE data points. It should be noted that the divergent results were not included. For linear basis, both MLS and SRMLS have second order of accuracy for function value interpolation. Not only the function value but also all derivatives approximations of SRMLS with quadratic basis have about one order higher than the solutions with linear basis.
(a)
(b)
(c)
Generally, the RMSE of SRMLS decreases with the particle spacing. MLS has lower stability than SRMLS as excepted. The results with disturbance constant less than 1.0E – 3 are close to the others. Finally, disturbance constant, β, was chosen as 1.0E − 5 in next numerical exercises.
3.2. Dependence of the Error on Kernel Function
Kernel function plays an important role in SPH simulation, and it is better to choose the same weight function for both SPH and SRMLS to simplify the coupled algorithm. In order to see how the type of the kernel function can influence the SRMLS interpolation accuracy in this subsection, the interpolations were repeated with five weight functions which are widely applied in SPH, they are the quintic spline function [5]:where for 2D case and for 3D case, .
The cubic spline function [3]where for 2D case and for 3D case.
The quadratic smooth function [5]where for 2D case and for 3D case.
The renormalized Gauss kernel functionwhere for 2D case and for 3D case.
The Wendland C4 function [32] is given bywhere for 2D case and for 3D case. They are denoted as Quintic, Cubic, Quadratic, Gauss, and Wendland, respectively, in this paper.
In this subsection, the 2D square was discretized with 40,000 irregular particles which were generated by using SPH pretreatment. The particle interval dx = 5.00E − 3 and the smooth length h = 1.0dx. 35 (5 × 7) cases were carried out.
Figure 7 gives the comparison of the RMSE of function and derivate values with various kernel functions and the SRMLS relative error contours with renormalized Gauss kernel function and linear basis. As expected, the interpolation of function value has higher order of accuracy than the approximation of derivative. For the function value interpolation, the worst results are obtained with CSPH, and the orders of accuracy with K2SPH, KGFSPH, and linear basis MLS and SRMLS are nearly the same with different kernel functions. In terms of the first-order derivate results, in addition to RMSE of linear basis MLS with Wendland C4 kernel function, the approximation derivates in each direction with CSPH, K2SPH, KGFSPH, and linear basis MLS and SRMLS are almost the same. For numerical exercises with quadratic basis, only the accuracy of SRMLS result improves. The SRMLS relative error contours are also presented in Figure 7. It can be observed that more accurate results could be obtained in the zones with small surface curvature, and high relative errors always locate in the fields with high surface curvature or close to boundary.
(a)
(b)
(c)
Benefit from larger support domain, namely, more particles contribute to the interpolation, a higher order of accuracy obtained with quintic spline and renormalized Gauss kernel function than the others in the calculation with SRMLS and quadratic basis. Although, Wendland kernel function is highly acceptable to SPH researchers, it was not chosen as weight function in this paper because of its poor performance in both MLS and SRMLS. As a result of good performance for interpolation and SPH simulation, the renormalized Gauss function was applied in this study.
3.3. Dependence of the Error on Support Domain Size
With the conclusion in last subsection, better results were obtained with the quantic spline and renormalized Gauss kernel function. This may benefit from larger support domain. So it is essential to study and assess the influence of the support domain size on interpolation stability and precision.
Same with the setups in the last subsection, the domain was filled with same irregular particles, namely, total particle number is 40,000; the renormalized Gauss kernel function was employed. The effect of support domain size was investigated with six smooth lengths, h = 0.9dx, 1.0dx, 1.1dx, 1.2dx, 1.3dx, 1.4dx. 42 (6 × 7) cases were carried out. The comparison of the RMSE with different smooth lengths and the SRMLS relative error contours with h = 1.2dx and quadratic basis are presented in Figure 8.
(a)
(b)
(c)
Same conclusions that the interpolation precision of function value is higher than that of derivative and the function value approximations with K2SPH, KGFSPH, and linear basis MLS and SRMLS are nearly the same with various smooth lengths can be drawn. Besides the results of MLS and SRMLS with quadratic basis, the interpolation precision slightly decreases with the smooth length. While the interpolation precision with MLS and quadratic basis increases with the smooth length, none of the derivate results are acceptable. As expected, the best results are calculated with SRMLS and quadratic basis, and the interpolation precision slightly decreases with the smooth length for function value estimation.
To further show the advantages of the SRMLS method with quadratic basis, relative error contours with h = 1.2dx and renormalized Gauss function are also plotted in Figure 8. It is clear that relative error is highly influenced by surface curvature.
The corresponding maximum condition numbers of the moment matrices with different smooth lengths are summarized in Figure 9 and Table 3, where is the average number of particle number in support domain. Method with quadratic basis can obtain larger condition number than that with linear basis, and stable regularized method could decrease the condition number significantly. All the maximum condition numbers decrease with respect to the smooth length. Finally, h = 1.2dx will be used in the coming calculations to ensure enough particles in support domain, acceptable interpolation accuracy, and time consumption.
3.4. Dependence of the Error on Spatial Resolution
Generally, verification should be taken out for all numerical simulations, and multi cases should be tested with various domain discretization resolutions, namely, the grid or particle sizes for grid-based or meshless methods, respectively, to estimate the discretization error in numerical uncertainty study. Consequently, it is important to check the dependence of error on spatial resolution.
Fifteen sets of particle spacings, dx = 4.00E − 2, 2.00E − 2, 1.33E − 2, 1.00E − 2, 8.00E − 3, 6.67E − 3, 5.71E − 3, 5.00E − 3, 4.00E − 3, 3.33E − 3, 2.50E − 3, 2.00E − 3, 1.67E − 3, 1.25E − 3, and 1.00E − 3, with the renormalized Gauss kernel function and h = 1.2dx were used in this subsection. The spatial resolution parameters are summarized in Table 4.
The maximum condition numbers of moment matrices are presented in Figure 10. Simultaneously, the comparison of the RMSE with different particle intervals and the SRMLS relative error contours with renormalized Gauss kernel function, dx = 2.50E − 3, and quadratic basis are shown in Figure 11. The right-hand side of Figure 11 is the relative error contour with irregular particle, renormalized Gauss kernel function, dx = 2.50E − 3, and quadratic basis. It can be found that the error distribution is similar with contours in Figure 8.
(a)
(b)
(c)
Furthermore, to highlight the advantages of SRMLS method with quadratic basis and check its performance with more violent condition, the surface fittings were repeated with random particles. Figure 12 shows the maximum condition numbers of moment matrices with various particle spacings and three particle distributions. The maximum condition numbers with regular nodes are close to those with irregular particles. Maximum condition numbers keep the same with respect to the spatial resolution for SRMLS with both regular and irregular nodes; however, they increase with the spatial resolution for MLS. Significant oscillations can be observed with random distributions in MLS and SRMLS.
Apart from the results with MLS, all the interpolation precisions increase with the spatial resolution. Same conclusions that K2SPH, KGFSPH, and linear basis MLS and SRMLS obtain almost same order of accuracy for both function value and first-order derivatives interpolation can be drawn. In MLS, the maximum condition numbers larger than 1010 for these cases with quadratic basis when particle interval less than 1.33E − 2, which means the moment matrices may be singular and could not be inversed. That is the reason for the RMSE increase with spatial resolution when particle spacing less than 1.33E − 2. The maximum condition numbers of SRMLS are acceptable and almost unchanged respect to spatial resolutions; consequently, the corresponding results have the highest accuracy with expectation.
Table 5 lists the calculated convergence rates of function value and first-order derivatives obtained from the different interpolation methods with irregular particles. For all methods, the convergence rate corresponds to the slope resulting from a least-squares fitting of the RMSE data points. However, the convergence rate of RMSE with MLS and quadratic basis is not calculated as the results are not converged. It can be observed that SMLRS with quadratic basis obtains fastest convergence rates for both function value and first-order derivatives estimations. Conversely, CSPH exhibits very poor convergence rate. K2SPH, KGFSPH, MLS, and SRMLS with linear basis achieve same order of accuracy.
Figure 13 shows the comparison of the RMSE between SRMLS and MLS methods with quadratic basis for different particle spacings and patterns. Function value, first-order and second-order derivatives are illustrated. It implies that the numerical results with MLS and quadratic basis diverge from the actual values when particle spacing is less than 2.00E – 2. Refer to the maximum condition numbers in Figure 12; SRMLS and MLS obtain same accuracy for regular and irregular particles. The random distributions increase the venture of irreversible and decrease the interpolation accuracy; however, the SRMLS can overcome this drawback.
(a)
(b)
(c)
(d)
(e)
(f)
The convergence rates of the second-order derivatives with SRMLS and quadratic basis for different particle patterns are summarized in Table 6. The results show that SRMLS with quadratic basis exhibit approximately first-order convergence rates on second-order derivatives for different particle patterns.
To highlight the advantages of SRMLS, the function value, first-order, and second-order derivatives approximations of another test function are presented in Appendix B. Same performance can be seen in those results. Some explanation should be stated for that the performance of SRMLS is the same for regular and irregular particle distributions. The reason is the special production method of the irregular particles in this paper. The irregular particles are generated from random particles by using SPH pretreating method. Generally, this SPH pretreating method is used to create quasiuniform initial particles in domain with complicated boundaries. Actually, the irregular particles are quasiuniform. Their patterns are different from regular particles, so they are referred as to irregular particles.
4. Numerical Results and Discussion
Four 2D test cases are reported to access the adaptability and convergence properties and to validate the coupled algorithm in comparison with analytical solutions and reference results obtained with other solvers in this section.
In the first test (Taylor–Green vortex flow) the algorithm was computed with periodic boundary; then the method was checked with external and internal moving boundary in the second simulation (Lid-Dirven cavity flow) and the third case (moving square inside a rectangular box), respectively. The last case (solitary wave slamming) was carried out to verify the robustness and reliability for the actual problem with violet free surface deformation.
For each test case, a smooth length h = 1.2dx has been chosen, except where otherwise specified. The renormalized Gauss kernel function was applied as the high performance in both SRMLS and SPH applications. To maintain the stability of SPH simulation, a particle shifting technique [33] was adopted in all simulations. Then interpolation in the schemes of boundary treatment, density re-initialization, and postproceeding were dealt with SRMLS with quadratic basis.
4.1. Taylor–Green Vortex Flow
The Taylor–Green vortex test is an unsteady flow of decaying vortices, which has an exact closed form solution of the incompressible Navier–Stokes equations. The domain is a square [0, L] × [0, L], and the Reynolds number is defined as Re = UL/υ. Periodic boundary conditions are enforced on the domain boundaries and the 2D solution is given bywhere u and are the components of the velocity u and U is initial the maximum velocity. Integrating in space the equation (29), one obtains the theoretical kinetic energy decay:where E0 is the initial kinetic energy.
In this test, L = 1 m, U = 1 m/s, and the corresponding Reynolds number is Re = 100. The reference length and time are, respectively, L, L/U; the pressure is made nondimensional through . Four particle sizes with dx = 2.00E – 2, 1.00E – 2, 5.00E – 3, and 2.50E – 3 m were carried out. The corresponding computing settings are summarized in Table 7, and the initial condition of Taylor–Green flow is shown in Figure 14.
(a)
(b)
In Figure 15, the convergence for the kinetic energy and max velocity magnitude decay predicted by the SPH is compared against the analytical solution for Re = 100. The discrepancy between the numerical simulation and the analytical solution is very small already for dx = 5.00E – 3 m.
(a)
(b)
Figure 16 shows the comparison of the velocity magnitude and pressure predicted by the present SPH for different time instants. A particle spacing of dx = 1.00E – 2 m has been adopted in order to highlight the particle configuration during the simulation. In those plots, the pressure fields are free of noises and smooth, and the symmetry results of both velocity magnitude and pressure exhibit the high performance of this coupled algorithm for the simulation with periodic boundary. Simultaneously, the particle shifting technique was also tested in present simulation, and no particle collection and cavitation were observed during the case running.
(a)
(b)
4.2. Lid-Dirven Cavity Flow
Lid-Dirven cavity flow is a widely used validation test case used for other numerical codes, e.g., finite volume, finite difference, and spectral element. This makes the problem ideal for assessing the ability of various SPH techniques in comparison with conventional methods where the results are well known. However, this case is particularly challenging for Lagrangian particle methods since the solution is widely ruled by the presence of vortices which are particularly difficult to capture. Indeed, particles moving exactly along the flow trajectories tend to create Lagrangian coherent structures, and this effect worsens as the Reynolds number gets higher.
The problem consists of a 2D rectangular cavity filled with fluid that is covered by a horizontal moving lid as shown in the schematic in Figure 17. Here, L = 1 m, U = 1 m/s, υ = 0.001 m2/s, then the corresponding Reynolds number is 1,000. The sensitivity of solution to the spatial resolution was investigated using dx = 1.00E – 2, 5.00E – 3, and 2.50E – 3 m. The case computing parameters are summarized in Table 8. For comparison, the case was also simulated with two grid-based finite difference methods, Marker and Cell method (MAC) [34] and Constrained Interpolation Profile method (CIP) [35], in which 100 × 100 Cartesian grids and center difference were applied.
Figure 18 shows the solutions of the horizontal velocity along the central vertical line and the vertical velocity along the central horizontal line with three spatial resolution levels being presented and compared to the solutions from Ghia et al. [36] and mentioned mesh-based methods. As expected, SPH solutions are in good agreement with the reference solutions from Ghia and CIP methods, better than MAC with center difference.
(a)
(b)
Note that the contours projected from Lagrangian solution with a spatial resolution 400 × 400 particles are plotted in Figure 19. They provide a view of the steady flow obtained with SPH. The left figure presents the velocity magnitude contour where streamlines are plotted and significant vortex can be found in the bottom corners. While the right plot displays the pressure contour, the pressure field is smooth.
(a)
(b)
4.3. Moving Square inside a Rectangular Box
The flow across a square cylinder is an important test case with interior moving boundary for the validation of separated flow especially in the turbulent regime. In the past, it has been regularly used to validate Large Eddy Simulation (LES) models. The test proposed is a tough benchmark for Lagrangian solvers like SPH as well. The main goal is to prove that the coupled particle solver adopted is able to generate diffuse and convect correctly the unsteady vorticity. Due to the presence of sharp edges, the vorticity generated is quite intense. Even though the geometry is very simple, the flow can be very complex; because of increasing the Reynolds number, the vorticity decays at lower scales. The Reynolds number of 100 is here analyzed.
The simulation domain adopted is represented in Figure 20, the red zone is the solid square while the black boundary of the square is the initial position of the moving box at time t = 0. The square moves rightward for t > 0 with a prescribed motion which is presented in Figure 21, starting from rest and accelerating to a final steady maximum velocity. This test is characterized only by no-slip penetrable solid boundary conditions on all walls of the tank and the square.
The convergence study was carried out with three particle spacings, dx = 2.50E − 2, 1.67E − 2, and 1.25E – 2 m. The computing parameters are summarized in Table 9. The viscous and pressure components (Cvd and Cpd) of drag coefficient are plotted against time in Figure 22. The drag coefficient , where Fx is the total horizontal force. The black circle is the result of the FDM Navier–Stokes solver from Colicchio et al. [37]. It can be observed that the forces do not present oscillations due to penetrable solid boundary, and results are converged with the spatial resolution, while the peak of viscous drag coefficient is slightly larger than reference result.
Furthermore, the snapshots of predicted velocity magnitude and vorticity fields at t = 5.0 s (left) and 8.0 s (right) are shown in Figure 23. The vorticity contour is more noisy than velocity magnitude because of the error in the evaluation of the curl and of the irregular distribution of the particles. Nonetheless, the features of the field are well captured.
(a)
(b)
4.4. Solitary Wave Slamming
Solitary wave is a kind of extreme wave phenomenon in coastal and ocean engineering problems which is used to represent certain characteristics of the tsunami wave. It is a good choice to study the phenomenon of wave propagation and transformation over a slope [38]. In consideration of the real application and free surface problem that SPH does well in, generally, solitary wave slamming was carried out in this section. The solitary wave propagation over a constant water depth and impact on a solid wall with inclination angle of 120° was considered to verify the effectiveness and stability of the coupled SPH method.
A piston wave maker is in the left side of the tank, and the motion of this plate can refer to Ma and Zhou [39]:where , , and . The analytical solution for the wave profile can be derived from the equation:in which is the solitary wave velocity, η is the water surface elevation, λ is the wave amplitude, and d is the water depth.
The wave simulation physical model is shown in Figure 24, and a model with inclination angle α = 120° on the right wall was simulated. The common basic information of the model includes that the length of the tank model is 10 m, the depth of the water is 0.25 m, a wave gage is located at x = 2.0 m, and two measurement points are located on the right wall at a distance of y = 0.05 m (P0) and y = 0.15 m (P1) to monitor the pressure time history for comparison with experimental results that was finished in Harbin Engineering University.
In this part, the solitary wave amplitude is λ = 0.15 m, the particle spacing dx = 1.00E – 2 m, so the particle size is 1.00E – 4 m2, total fluid particle number is 24,820, and the time step is equal to 1.00E – 4 s. The simulation was repeated with four smooth lengths, h = 1.00dx, 1.10dx, 1.20dx, and 1.30dx. The computing settings are summarized in Table 10.
The comparison of solitary wave profiles at different time is shown in Figure 25, and the time history for numerical and experimental wave elevation at wave gage is presented in Figure 26. According to the results, the numerical wave elevations meet well with both analytical and experimental results.
In order to illustrate the adaptability of coupled SPH algorithm for free surface problem clearly, Figure 27 gives the comparison of experimental and numerical pressure at P0 and P1. The results are converged with smooth length. Furthermore, Figure 28 gives the comparison of present pressure contour near the end with result of Chen et al. [28] at t = 5.4 s. Comparing with the numerical results of Chen et al. [28], which were without SRMLS, the coupled SPH method performs better in the solitary wave slamming. Figure 29 presents some snapshots of the pressure distribution with streamlines at t = 4.0 s and 5.8 s. From the comparison, it is obvious that the pressure is smooth.
5. Conclusions
In this paper, a stable regularized moving least-squares (SRMLS) interpolation whose moment matrix can be always nonsingular is introduced and coupled with the SPH method. The properties of SRMLS with the irregular particles which were generated by a SPH pretreating method for surface fitting aare presented and discussed. Numerical solutions were compared with the results using moving least-squares (MLS) and three SPH methods, including CSPH, K2SPH, and KGFSPH. It was found that SRMLS not only has nonsingular moment matrix, but also obtains high accuracy result.
Afterwards, the capability of the coupled method with SRMLS and SPH was also illustrated and assessed through some numerical tests, including Taylor-Green vortex flow, Lid-Dirven cavity flow, moving square inside a rectangular box, and solitary wave slamming, wherein, four type boundary conditions, periodic boundary, free surface and external and internal moving walls, were checked with the presented algorithm. All the results show convergence with particle size or smooth length and obtain good agreement with available reference solutions.
The main conclusions of this paper lie in the following aspects. Firstly, SRMLS method always obtains higher accuracy results with quadratic basis than that with linear basis as expected. Secondly, a large support domain can improve the stability of SRMLS method, especially for the interpolation with quadratic basis, while the precision would decrease slightly. Thirdly, SRMLS have better results with the renormalized Gauss kernel and quintic spline functions than the others, as a result of the lager support domain. Fourthly, stabilization and regularization could decrease condition number of moment matrix significantly, namely, the interpolation stability would be improved and the singular matrix may be prevented as well. Finally, the quadratic basis, the renormalized Gauss kernel function and h = 1.2dx are recommended for SRMLS method.
Instead of conventional SPH kernel approximation, SRMLS deals with the interpolations for particle’s unknown physical characteristics in the schemes of SPH boundary treatment, density refreshing, and postproceeding. With the comparison of reference results, the coupled algorithm shows high adaptability with various boundary conditions.
Furthermore, it should be noted that there are some limitations in the present works as follows. (1) The inlet and outlet boundary which is widely used in ocean engineering is not tested with the coupled method. (2) The interpolation plays an important role in overlap particle technique which is meaningful for future large-scale simulation; therefore, it is worth to investigate the application of SRMLS for the overlap zone. (3) The SRMLS method can not only deal with mentioned interpolations, but also provide high-order approximations to improve SPH algorithm. (4) All the tested cases are limited in 2D, a 3D study will be carried out in the future. (5) SPH method is widely applied for free surface flow with large deformation in ocean engineering and coastal engineering such as green water problem, slamming, and wave run-up; more cases should be carried out for free surface flow applications with the coupled method.
Appendix
A
A.1. Property of SRMLS Method
Mirzaei [24] has proved that the MLS approximation can be implemented in a more stable fashion, if a shifted and scaled polynomial basis function is used as a basis for P. Here, this appendix proves the moment matrix, , in SRMLS is nonsingular.where , . Since E is a diagonal matrix with nonzero diagonal elements, one can haveas
It leads to
From equation (18), the rank of matrix is always m. Finally, it proves that matrix will be nonsingular.
A.2. Derivate of Shape Function
The shape function in SRMLS approximation is given as
The first and second-order derivatives can be calculated byin which, , , , , and .
Some operations are required to get the derivatives of . It should be noted that . Taking the derivative of both sides, one can have
Then it yields
At last, the second-order derivate of is written aswhere the index following comma is a spatial derivative.
B
To highlight the advantages of SRMLS, numerical exercises were carried out for the second test function:defined over the domain and . The shape of this function is shown in Figure 30.
Figure 31 shows the comparison of the RMSE between SRMLS and MLS methods with quadratic basis for different particle spacings and patterns, function value, first-order, and second-order derivatives being illustrated. The convergence rates of the function value, first-order and second-order derivatives approximations with SRMLS and quadratic basis for different particle patterns are summarized in Table 11. All methods perform same with the tests of Franke’s function.
(a)
(b)
(c)
(d)
(e)
(f)
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 there are no conflicts of interest regarding the publication of this paper.
Authors’ Contributions
H. Jiang and Y. Chen made the computations and data analysis. S. Jin finished the preliminary code in his master study. X. Zheng and Q. Ma guided the research work and proof-read the paper. All authors contributed to the work.
Acknowledgments
The study was supported by National Natural Science Foundation of China (Grant no. 51909039) and Shandong Provincial Key Laboratory of Ocean Engineering (Grant no. 201807).