Review of smoothed particle hydrodynamics: towards converged Lagrangian flow modelling

1. Introduction

Computational fluid dynamics (CFD) is an established capability in most areas of engineering: aerospace, automotive, nuclear, chemical, offshore, marine, hydraulic, etc. The methods are mesh based with a range of turbulence models generally available with various codes, e.g. ANSYS and OpenFOAM. For many applications multi-processing is standard, and for complex geometries mesh generation is a time-consuming preliminary activity. Methods for handling interface problems have been developed, through volume-of-fluid and level-set methods, but such methods are not naturally suited to mass conservation. Meshless methods have also been developed with either vortex or mass particle tracking as Lagrangian formulations. Vortex methods have the advantage of inherent vorticity conservation in two dimensions, but in three dimensions, this is not the case and there are greater requirements for a numerical solution. They are also not well suited to multi-phase problems. Methods with particles based on mass have however developed rapidly, showing versatility for diverse problems. The most common form is known as smoothed particle hydrodynamics (SPH) and is the subject of this review. This evolved from astrophysics in the 1990s for engineering applications, showing the versatility and robustness of the method, until recently without emphasis on accuracy and convergence. The early method was of weakly compressible form.

Since 2016, six review papers focusing on different elements of SPH have been published, indicating its increasing maturity and application [1–6]. These papers show the diversification of SPH in terms of both its applications and the challenges in securing the underlying theoretical basis. Based on these reviews and the experience of the authors and many others mentioned in this article, the advantages were that:

• (i)

  A free surface or interface is automatically followed and does not require special treatment; for example, for breaking waves in the ocean.

• (ii)

  The method is inherently mass conservative.

• (iii)

  Dynamic body interaction is handled automatically since the formulation is Lagrangian.

• (iv)

  Fluid properties associated with different phases may be easily incorporated.

• (v)

  The method has well-characterized numerical stability.

There are a number of perceived weaknesses, as listed below. These areas remain the subject of intense research, and in the relevant sections of this review the latest work in addressing these weaknesses is highlighted and discussed:

• (i)

  low convergence rate (§§5 and 6),

• (ii)

  approximation of boundary conditions (§7),

• (iii)

  lack of turbulence modelling (§8),

• (iv)

  large pressure fluctuations due to spurious pressure waves in the untreated weakly compressible formulation (§9),

• (v)

  formation of voids (§10),

• (vi)

  high computational cost (§§6, 11 and 12).

In addition to investigating natural phenomena benefitting from the Lagrangian description of motion (e.g. Monaghan [7] on astrophysics, Vacondio et al. [8] on flooding and Xenakis et al. [9] on landslides and tsunamis), SPH is increasingly used for applications where safety is an important concern, such as nuclear and renewable energy, civil and environmental engineering, aerospace and automotive industries, and chemical engineering and food production. Existing simulation tools in these fields of application are based on the finite-element method (FEM), finite-volume method (FVM), finite-difference method and spectral methods which are mature algorithmically and are well understood with a theoretically sound foundation and with well-established high-order convergence. However, these industries are increasingly considering the use of particle methods, such as SPH, to overcome certain limitations of these methods to address previously intractable problems. Naturally, this has led to questions about the accuracy of SPH and most importantly its convergence properties.

In many applications, demonstration of higher-order convergence is an essential pre-requisite to widespread use. The historical precedent took place in the 1980s and 1990s with numerical methods finding widespread use in the aerospace industry. Uniform reporting of the convergence behaviour of numerical methods became a critical issue [10], leading to the emergence of standardized methods such as the grid convergence indicator to present the convergence behaviour of simulation techniques in a consistent and comparable manner. Academic journals went even further, implementing editorial policies whereby manuscripts submitted would not be sent to reviewers if there was no convergence study [11]. The same has yet to happen for SPH [2].

As will be explained in §5, research into assessing and improving the convergence and accuracy of SPH is in full swing. Since the range of applications of SPH is so broad, the convergence of SPH must take into account all the factors contributing to the behaviour of the numerical solutions. Figure 1 shows the factors considered by the authors to be the most important to make progress towards higher-order convergence of SPH. These include physical aspects from boundary conditions, viscosity and turbulence, the basic choice of SPH formulation for fluids (either incompressible or weakly compressible) and novel numerics with higher-order kernels and manipulation of the truncation error terms to pragmatic developments essential for application to real problems, including adaptivity, coupling and hardware acceleration.

This review focuses on progress towards converged solutions. The general formalism of SPH is first described, then the weakly compressible form of a solution is defined, followed by the incompressible form. Formulations to high order, increasing accuracy, are then summarized. Particle adaptivity for variable resolution is described, although this has received limited attention to date. The next section describes the range of boundary condition formulations. The vital aspects of viscosity and turbulence follow. The effective method of stabilization of the weakly compressible form due to diffusion in the continuity equation is then described, followed by the particle regularization method (known as shifting), which is required for stabilization of the incompressible form and generally improving accuracy of both forms. Coupling of complex inner flow domains by SPH with efficient outer domain solvers for simpler flow physics is then described. The important topic of run-time acceleration by parallel processing on CPUs and use of GPUs (graphics processing units developed for the computer games industry) follows. Finally, some simulations are presented showing converged solutions, and conclusions are drawn.

2. Basic smoothed particle hydrodynamics concepts

In order to understand the challenges involved in improving the convergence properties of the SPH method presented in §§5–13, the basics of the SPH methodologies are presented herein. For a full description, the reader is referred to texts such as Violeau [12].

The basic principle of the SPH formulation is the integral representation of a function f, which may represent a numerical or physical variable defined over a domain of interest Ω at a point x. The integral approximation or kernel approximation reads [13]

where h is defined as the smoothing length that characterizes the size of the support domain of the kernel and W is the weighting or kernel function. The kernel function is chosen to be a smooth, isotropic and even function with compact support (i.e. the finite radius of influence around x). Note that this is an important advantage over vortex methods.

In some basic analysis using a Taylor expansion, Monaghan [13] showed that, with positive kernels, the error of the integral approximation in equation (2.1) is O(h²)

where f_I is the function interpolated by equation (2.1) and f_E is the exact value. This is the important feature of the SPH interpolation process that conventionally has imposed a theoretical upper bound on the convergence order of 2. From observation in the literature, reviewed values of 1.1–1.8 are quite common (e.g. [14–16]). This is primarily due to the errors introduced by the discretization process described below and the irregular particle distribution [17]. In Lagrangian simulations, the particle distribution is inevitably irregular to some degree and state-of-the-art treatments will be reviewed in §§9 and 10.

In a discrete domain, equation (2.1) can be approximated by using an SPH summation in the form of

where V is the volume of the particle expressed as the ratio of the mass m to density ρ and N is the number of particles within the support. Throughout this paper, the subscript i denotes the interpolating particle and j refers to the neighbouring particles. The 〈…〉 symbol denotes an SPH interpolation and will be dropped for simplicity. The final form of the particle approximation in a discrete form is

where W_ij = W (x_i − x_j, h) and f_j = f(x_j). The discretization process of the above equation introduces another source of error (see Quinlan et al. [17] for an in-depth investigation).

For gradients, the following form of the first-order gradient operator is widely used:

where ϕ is a general variable and ∇_iW_ij is the gradient of the smoothing kernel. A similar expression is obtained for the divergence operator by decomposing the operator into a series of gradient operators.

Higher-order derivatives are formulated in terms of combinations of first-order SPH kernel derivatives and finite differences. The Laplacian operator employed for a quantity b is that suggested by Morris et al. [18] and will be referred to as the Morris operator,

where η = 0.001 h to avoid singularity as r_ij → 0 and r_ij is distance. Here, a is a scalar variable.

The Morris operator is generally employed for internal flows. However, it is well known that the Morris operator is inaccurate near free surfaces and requires improved formulations necessary to improve the convergence properties. Near a free surface, the Schwaiger [19] operator has been shown to provide more accurate interpolation by correcting for the truncated kernel support. Specifically

where Γ is defined as

where superscripts α and β denote the spatial coordinates. More details on the performance of the operator can be found in Lind et al. [20] and Schwaiger [19]. It should be noted that both Laplacian operators are non-conservative.

There are two forms of SPH solver for the Navier–Stokes equations, weakly compressible SPH (WCSPH) and incompressible SPH (ISPH), the same as for mesh-based CFD. Espanol & Revenga [21] provide a theoretical connection between the incompressible and compressible contribution to the viscous term and readers are recommended to read this seminal paper.

3. Weakly compressible smoothed particle hydrodynamics

The Navier–Stokes equations can be written in Lagrangian form as

and

where u is the fluid velocity vector, ρ is fluid density, t is time, P is fluid pressure, ν is kinematic fluid viscosity and g is gravity.

In the WCSPH formulation, a stiff equation of state (EoS) is often used for closure of the Navier–Stokes equations. The EoS relates fluid pressure to changes in density via

where γ is the fluid polytropic index so that small variations in density lead to large variations in pressure, ρ_o is the reference density and c_o is the reference numerical speed of sound. c_o is typically set to c_o = 10U_max, where U_max is the maximum expected flow velocity such that the overall flow Mach number Ma ≤ 0.1, giving an expected density variation of less than 1%. If γ is set to a value of 1, then the EoS simplifies to p = c_o²(ρ − ρ_o). First used in SPH by Morris et al. [18], this has been found to be more amenable to developing theoretical concepts and understanding within the WCSPH approach (e.g. [22]).

The discrete Lagrangian WCSPH Navier–Stokes equations can be written as [13]

and

where Π_ij denotes the viscous terms in the discrete SPH formulation. The final term in equation (3.5), δ_i, is an artificial density diffusion term used to limit spurious density fluctuations. Such density fluctuations can cause spurious pressure fluctuations and various damping schemes have been applied but these so-called δ-SPH equations are becoming the preferred option (see §9).

4. Incompressible smoothed particle hydrodynamics

The Navier–Stokes equations in Lagrangian form for an incompressible fluid become

and

To enforce incompressibility and ensure a divergence-free velocity field, the projection method of Cummins & Rudman [23] is generally used. Note that incompressibility implies no volumetric variations; hence, density and mass are constant.

Firstly, particles are advected to an intermediate position x* with velocity uⁿ using

An intermediate velocity u* is calculated at position x* based on the viscous forces of the momentum equation

The pressure at time n + 1 is obtained by the implicit solution of the Poisson equation, which has the following form:

The discretization of the left-hand side of the Poisson equation (4.5) is performed with either the Morris or Schwaiger operator. Next the velocity at time n + 1 is obtained by the projection of the intermediate velocity onto divergence-free space such that

resulting in a divergence-free velocity field. Finally, the particles are advected to the n + 1 position by a central difference scheme

In ISPH, simulations can quickly become unstable as accurate computation produces clustering, notably near stagnation points (e.g. [24]). Particle regularization or shifting (described later in §10) is essential, producing accurate converged solutions. It should be noted that this is an entirely numerical device separate from the underlying physics and that this is a distinction from δ-SPH in WCSPH.

5. High-order and highly accurate smoothed particle hydrodynamics

While theoretical general proofs of convergence are being established (e.g. [25,26]), the practical convergence properties of particle methods and SPH started to receive attention in the mid-1990s with seminal works by Belytschko et al. [27], Liu et al. [28,29], Dilts [30,31], Bonet & Lok [32], Welton [33], Vila [34] and others looking at the role of the truncation error and the use of correction techniques to ensure polynomial consistency up to nth order and convergence. While these approaches certainly provided an option to correct for the SPH error and reproduce functions consistently, in practice they have been rarely used because SPH was viewed wholly as a Lagrangian method, with construction and solution of the correction matrices a costly operation at each time step. Furthermore, simulations were far more difficult to keep stable. The reasons for this are still not entirely understood (although the papers of Violeau & Leroy [35,36] have provided some insight into possible time integration issues). During the 2000s, focus shifted towards expanding the range of applications for SPH, so the convergence properties received limited attention (with a few exceptions, e.g. [37]). It was not until the mid-2010s, with the SPH community in possession of increasing experience in its application, that the limitations of the SPH convergence behaviour started to receive renewed attention [38,39], in particular the effect of particle disorder.

Lind & Stansby [40] showed the very high accuracies possible for fixed (Eulerian) particles in distributions that minimize the error of particle disorder. High-order convergence (fourth and sixth order) was achieved using new non-monotonic kernels, easily implementable and enabling errors that hit machine precision readily for idealized viscous fluid flows (such as Taylor–Green vortices). The attractive Lagrangian features of SPH can still be called upon, and Lind & Stansby [40] offered a mechanism (demonstrated with a propagating wave test case) where a fixed (or Eulerian) region of particles may transition through an arbitrary Lagrangian–Eulerian (ALE)-like approach to become fully Lagrangian at the free surface. It has since been demonstrated that the robustness of SPH offers further flexibility in how one may transition between Eulerian and Lagrangian regions: even a sharply discontinuous Eulerian–Lagrangian interface (including treatments such as shifting in the Lagrangian region) provides a stable and readily implementable Eulerian–Lagrangian SPH scheme of superior accuracy to the fully Lagrangian equivalent [41].

Mixing Lagrangian and Eulerian discretization is not new in CFD, e.g. the marker-and-cell method, but the advantage here is that calculations are done over a node (particle) of one data type that exists in a single layer (without overlapping or buffer zones), with programming and book-keeping greatly simplified. Nevertheless, in the search for improved accuracy, the use of background meshes for the most error-prone parts of particle-based algorithms is being revisited, most recently in projection-based particle methods [42]. Furthermore, rather than a Cartesian grid, Voronoi cells are increasingly being considered as a means to improve accuracy. In Fernandez-Gutierrez et al. [43], for example, a Voronoi discretization replaces that of SPH near boundaries in order to improve the consistency and implementation of boundary conditions, although this approach does require an overlap or buffer zone between SPH and Voronoi regions, with interpolation between the two regions that is not necessary in wholly SPH particle approaches like Fourtakas et al. [41]. The iterative shifting approach of Vacondio & Rogers [44] is emerging as the most likely means of achieving high-order accuracy in domains where uniform or Cartesian distributions are difficult to apply, as well as enabling high-order quasi-Lagrangian simulations. Shifting iteratively per time step to the point where particle distributions yield high (fourth)-order convergence has been demonstrated recently for select internal flow problems using ALE-WCSPH [16]. The extension to high-order free-surface flows is currently being explored.

As an alternative to controlling particle distributions, some studies have recently opted for construction of local stencils between particles upon which weighted essentially non-oscillatory (WENO) schemes are applied for increased accuracy in calculating particle interactions [45]. While accuracy increases are available through high-order WENO schemes, the associated SPH interpolants have so far remained low order, and so global high-order error reduction is restricted. Understanding the complex interaction between WENO reconstruction error and SPH interpolation error will be necessary for further accuracy gains in schemes such as this. Indeed, as high-order approaches begin to increase in popularity in SPH (including those available via implicit consistency correction schemes, such as Sibilla [46]), greater insight into SPH error characteristics will be needed to maximize benefit. Compared with algorithm development, there is a dearth of work in this area, but some progress is being made. Violeau et al. [47] undertook an analysis of the spectral properties of the SPH Laplacian (Morris) operator, and obtained results useful for the theoretical analysis of convergence of the SPH Poisson equation, as applied in ISPH schemes. Underpinning any numerical analysis in SPH is consideration of the smoothing error, and, recently, this has been determined exactly as a function of kernel standard deviation dependent on the kernel Laplace transform [48]. The technique has been applied to key SPH operators, e.g. gradients and Laplacians, and will no doubt be a key stepping stone to formal SPH error analysis in more practical contexts.

Accuracy increases through the use of high-order kernels or WENO reconstructions may be viewed as the SPH equivalent of ‘p-refinement’, to borrow terminology from high-order numerical schemes, such as hp-FEM. However, when seeking high-order solutions, it is important to bear in mind that there may be regions of the flow where high-order convergence is not demonstrable or even possible given singularities in higher derivatives of the primitive variables or geometry. In this event, a return to particle (or ‘h’) refinement is the obvious way to increase spatial accuracy, and doing this in an efficient and cost-effective way, e.g. only when and where needed, has resulted in a number of important studies on particle refinement adaptivity. The next section briefly considers some recent and important work in this area.

6. Adaptivity

Implementations of mesh/grid adaptivity are mature in conventional CFD, but the equivalent based on particles is relatively new in SPH. In this context, we consider only spatial adaptivity where particle number and/or volume are required to differ in different regions of the flow according to simulation requirements (e.g. error minimization or computational savings). Basic implementations of this have been around for some time, where, in known regions of interest, particles of smaller size are introduced prior to the simulation to improve resolution (see, for example, Omidvar et al. [49]). The real goal here, however, is a dynamic adaptivity to rival conventional CFD, including minimization of key error measures in two and three dimensions. The most appropriate error measures and criteria for fully dynamic adaptivity are still being explored, but Vacondio et al. [50] have made important headway in this area. Extending the two-dimensional (2D) dynamic particle splitting and coalescing study of Vacondio et al. [14], optimal three-dimensional (3D) particle distributions of daughter particles—those particles that have formed from the splitting of a larger ‘mother’ particle—have been determined. In this case, density error is minimized while mass and momentum are conserved during the particle-splitting process. Of course, lessons can be learned from mesh-based methods, and Chiron et al. [51] employ concepts from adaptive mesh refinement in the development of a robust adaptive particle refinement (APR) scheme that improves computational time while retaining a data structure suitable for parallelization. However, interpolation between particles of different refinements uses a Shepard filter—a low-order interpolant. The effect on accuracy in this context has not been fully explored, and, indeed, the use of APR with Shepard filtering for ALE-SPH schemes is not recommended.

In the majority of cases, adaptivity is used to increase particle resolution near boundaries (both solid and free), as this is often where the most interesting and challenging physics resides. Therefore, in the pursuit of converged and high-order solutions, even with careful adaptivity and locally high particle resolutions, accurate boundary conditions (at second order and above) are essential, and this remains a subject of intense research.

7. Boundary conditions

Boundary treatments affect the convergence behaviour of all SPH codes where boundaries should be imposed. The original SPH formulations were derived for astrophysics simulations where boundaries do not exist. For engineering applications, the truncation of the SPH kernel support in the vicinity of boundaries has presented a major challenge.

The most common boundary conditions are termed ‘dynamic’, where dummy particles are fixed in the surface, typically with two rows to provide kernel support for the fluid particles, and treated as fluid particles except that they are given the body velocity (zero if stationary) [52]. This is robust and convenient but does not accurately predict the surface shear stress or zero normal pressure gradient. For many water wave problems with inertia dominant, this may be acceptable. Adami et al. [53] presented a useful modification of this approach that takes advantage of the ease of implementation of ‘dynamic boundary particles' but with values imposed on the boundary particles that take into account the flow properties, such as the pressure field, leading to improved behaviour.

Mirror particles, analogous to images in vortex methods, are quite accurate as body particles are created that mirror those in the flow close enough to the surface to be within their support domain. For no-slip boundary conditions, the velocities are equal and opposite and the normal pressure gradient is zero; for slip boundary conditions, the tangential velocities are simply equal. This is, however, numerically cumbersome as an uncertain number of new particles are created at each time step. Furthermore, mirror particles are difficult to implement in complex domains, and have traditionally been restricted to comparatively simple geometries. Some local mirror techniques have been developed [54–56], but these also present challenges for broader generalization for 3D complex geometries. A hybrid approach is to keep the body particles fixed, as dummy particles, but provide flow properties from the interpolation points in the flow [57]. Fourtakas et al. [16] deployed a Local Uniform STencil (LUST) per fluid particle to complete support near boundaries, which employs symmetry arguments in much the same way as mirror particles, but avoids common mirror particle problems such as boundary particle overlap. The locality of the stencil enables complex geometries to be readily modelled, and suits parallelization on GPUs as well as extension to three dimensions [16].

Well-established methods in conventional CFD are increasingly being considered for SPH. For example, immersed boundary methods, quite popular in conventional CFD, are relatively new in SPH, but have recently been developed for the Lagrangian and Eulerian forms with promising results for thin slender bodies [58]. Principles from finite-difference approaches may also be used, and Zheng et al. [59] have employed a simplified finite-difference interpolation scheme for discretizing and imposing boundary conditions in ISPH without the use of mirror particles. Semi-analytical approaches that calculate factors associated with the truncated kernel at boundaries offer perhaps the most mathematically rigorous way of imposing boundary conditions [15,60], but at added complexity and computational expense. Implementations exist in two dimensions and Mayrhofer et al. [61] presented a 3D approach. While very accurate, these approaches are quite restrictive in that mathematical derivations are kernel dependent, restricting generality of application. Nevertheless, the accuracy offered is attractive and recent works have applied the semi-analytical approach to open as well as solid boundaries in weakly compressible [62] and incompressible [63] flows. Open boundaries in coastal and offshore applications will invariably contain a free surface, and it is here that one has to return to fully numerical (rather than semi-analytical) approaches. The recent approach of Tafuni et al. [64] is gaining popularity owing to its robustness, simplicity and flexibility in treating open boundaries. Inflows, outflows and mixed open boundaries, with or without a free surface, can all be modelled. Based on a buffer region approach with first-order accurate SPH interpolation, the method has been straightforwardly extended to three dimensions as well as parallelized on both CPUs and GPUs. This is built upon in work by Verbrugghe et al. [65], who, for the generation of nonlinear waves, include velocity corrections at the inlet and outlet to help reduce reflections in the flow domain from a ‘sponge’ layer that cannot absorb waves totally.

For the water surface, no special boundary condition is required for WCSPH as pressure → 0 as density → 0. Colagrossi et al. [66], however, showed that the choice of the SPH operators employed at the surface has important consequences for the convergence properties and theoretical considerations for the free-surface boundary condition. For ISPH, a zero or constant pressure condition is imposed for the Poisson solver for particles identified as being on the free surface by using the divergence $\mathrm{\nabla }.\mathbf{\text{r}}$ of position r for each particle with threshold values of less than 1.5 in two dimensions or 2.5 in three dimensions, 2 and 3 being values of internal flow where each particle has a complete kernel support. Lind et al. [67] coupled WCSPH with ISPH for air–water flows (termed incompressible–compressible SPH (ICSPH)), and, in this approach, water/air interface velocities are calculated across both domains and pressure from the compressible air is imposed on the water interface for the Poisson solver. In this way, velocity and pressure are continuous across the interface, as required without surface tension effects.

For offshore and coastal engineering applications, wave maker motion may be imposed directly, or velocities on a vertical plane may be imposed, with absorption in a damping zone or sponge layer over at least one wavelength [68]. Inflow and outflow conditions for currents are similar to those for mesh-based approaches [64]. Flows of this type are a popular application for SPH and are of characteristically high Reynolds number and often turbulent. The above discussion around boundary conditions is of course critical for fundamental studies of boundary layers and turbulent flow in SPH, but due consideration must also be given to the numerical representation of viscosity terms, as well as turbulence models more generally. With this in mind, some recent and important approaches to viscosity and turbulence modelling in SPH are now considered.

8. Viscosity and turbulence

Turbulence modelling and direct numerical simulation (DNS) remain key challenges in CFD research generally. Many leading DNS approaches are based on high-order grid-based methods (e.g. [69]), and it is in this area that high-order and converged SPH offers great promise in the longer term by providing accurate turbulent flow solutions for more diverse and practical applications involving complex geometries and free surfaces. For moderate to low Reynolds numbers, the laminar viscosity model as detailed above (equation (3.4)) is an appropriate choice to model viscous effects directly; the operator is a direct SPH discretization of the Laplacian, and quite straightforward to implement. Similar discretizations exist that include corrections to reduce the error: Fatehi & Manzari [70] provide a comparative study for SPH second derivatives, as well as a new scheme showing order of magnitude decreases in error. In recent years, studies around viscosity in SPH have focused on its use as eddy viscosity in turbulence modelling, but not exclusively. The particulate nature of SPH has allowed insights into the viscosity of fluids not offered by conventional schemes. Colagrossi et al. [71] consider a Hamiltonian particle system and use SPH to formulate a continuum approximation for the full Newtonian viscous stress, including shear and bulk viscosities. Based on momentum arguments, the system derived from the Hamiltonian implies that the bulk viscosity may be positive and in contradiction to the Stokes hypothesis. It is also found that predicted values of bulk viscosity appear close to experimental values for common gases and liquids, including water. The relationship between SPH and the underlying Hamiltonian or Lagrangian is often overlooked during practical engineering simulation, but its exploration is frequently worthwhile, offering deep, sometimes surprising, physical insights. There are indications that this holds true for turbulence modelling with SPH—generally, one of the great challenges of computational fluid mechanics.

Conventional CFD typically based on finite-volume or spectral methods remains mainstream for turbulence modelling, but SPH is making headway where it is hoped the intrinsic characteristics of SPH will offer insights unavailable to grid-based methods. As discussed in Di Mascio et al. [72], large eddy simulation (LES) is a natural option for SPH, given the similarity in frameworks between filtering and interpolating. Indeed, Di Mascio et al. [72] show that SPH smoothing can be viewed as LES Lagrangian filtering, provided terms additional to the LES convolution term are included in the governing equations. Mayrhofer et al. [73] showed that a minimum of 20 particles are required across a coherent turbulent structure to predict its evolution with LES to avoid incorrect reproduction of velocity–pressure interactions with the collocated nature of SPH. The affinity between SPH interpolation and LES filtering is not the only reason that SPH may be particularly well suited to turbulence modelling. Once again, promise is offered through energy considerations of the underlying discrete particle system—in this case the Lagrangian. By considering an SPH smoothed velocity in the Lagrangian of the particle system, Monaghan [74] offers an SPH model for turbulent flows that bears similarity to the Lagrangian-averaged Navier–Stokes (LANS) turbulence model. This approach offers an attractive particle-based model of turbulence that is shown to recover energy spectra well. Furthermore, the Lagrangian nature of the particles means the approach is well suited to the study of turbulent mixing and free-surface turbulent flows characterized by very large free-surface deformation [75].

The presence of viscosity in flows invariably has benefit in increasing stability and long-term accuracy. Indeed, this motivated development of the original artificial viscosity formulation [76], which converges to physical viscosity as particle spacing is refined. However, in recent years, it has become apparent that diffusion operators within the momentum equation are not the only way in which one may effectively stabilize a simulation. Applying an equivalent operator to the conservation of mass or continuity equation has proven particularly fruitful for many applications in WCSPH—originally known in the community as delta-SPH or δ-SPH.

9. Density diffusion schemes in weakly compressible smoothed particle hydrodynamics

Ever since its initial implementation in the form of a diffusion term given in equation (3.4) [77], and then the δ-SPH [22,78,79] scheme, density diffusion has become increasingly popular while undergoing a number of improvements. The first key improvement was the inclusion of corrections to normalize gradients in the density diffusion term to improve longer term accuracy [57].

In the δ-SPH method, as suggested by Antuono et al. [78], the density diffusion term δ_i is on the right-hand side of the continuity equation (3.5), which is given by

and

where a_d is the diffusion parameter, which is typically set to 0.1.

A second significant improvement to the accuracy and applicability of the δ-SPH scheme has been its use in combination with particle shifting, first explored in the context of variable resolution by Vacondio et al. [14]. The δ + SPH model [80], as the methodology has become known, uses δ-SPH in conjunction with the shifting algorithm, closely based on Lind et al. [20], to help eliminate effects such as the tensile instability and noise in velocity gradients (which persisted in the original δ-SPH scheme owing to any irregularity in particle distributions).

The δ + SPH scheme has since been applied to a wide range of flows, internal and free surface, and for a wide range of Reynolds numbers. Indeed, the effectiveness of combining δ-SPH diffusion terms with shifting is such that it has also been used successfully within the finite-particle method [81] for a range of viscous flows. Sun et al. [82] were able to approach high Reynolds number flows with δ-SPH, with shifting and a background pressure switch proving effective in modelling high Reynolds number flows past stationary and flexible bodies. Through the inclusion of an APR algorithm based on Vacondio et al. [14], Sun et al. [83] showed that δ + SPH can also be applied to a range of test cases for vortex-induced vibration involving multi-body interactions and large body movement. Again, good agreement with reference solutions is demonstrated. The δ parameter controls the diffusion within the conservation of mass equation and this accompanies the α parameter for artificial viscosity within the momentum equation [7]. While the choice of δ and α has been remarkably consistent and fixed for a wide variety of flows, a number of recent studies have looked more closely at the role and relationship of these diffusive terms. Meringolo et al. [84] determine that the parameters α and δ are closely related and, when used correctly, even quite different sets of parameter values result in similar total energy dissipation. Building on their inter-dependence, by employing an LES closure model, for turbulent flows at least, the parameters become dynamically determined, and in this framework may no longer be taken as free parameters. This is promising and indicative that δ-SPH, while often viewed as direct numerical diffusion of density, may at some level be a physical model for turbulent energy dissipation.

Motivated by a similar desire to remove empirical parameters in diffusion, Green et al. [85] showed that δ-SPH is exactly equivalent to a Riemann solver applied only to the continuity equation. Studies of internal energy are used to determine the most appropriate limiter, and, as desired, empirical parameters are avoided. The broad similarity in mathematical structure of this formulation to the δ-SPH scheme is unmistakable. Indeed, arguments have already been provided by Cercos-Pita et al. [86] that show equivalence of diffusion terms in the continuity equation with the Riemann solver SPH in their comparative study of density diffusion terms. This is an important observation that perhaps helps explain a revisiting of Riemann solver schemes in SPH in recent years (such as [87]). Needless to say, the work of Cercos-Pita et al. [86] helps to place δ-SPH-type schemes in an important mathematical context that is familiar to many in CFD more broadly.

The need to filter spurious pressure oscillations using δ-SPH or Riemann solvers arises in part from SPH formulations that rely on a very stiff empirical EoS alongside unphysical sound speeds (usually some orders of magnitude lower than the physically correct value). Kiara et al. [88] provide a comprehensive analysis of these oscillations and determine that, in some cases, the dynamics of these spurious acoustic solutions can scale with the speed of sound squared. Smoothing pressures using the treatments described above of course introduce artificial physical diffusion that may be unwanted. Removal of pressure fluctuations in post-processing has received comparatively little attention, but may offer a route to analysis of accurate pressures without the simulation dynamics being overly diffusive. Temporal filtering is suggested as an effective approach by Kiara et al. [88], and Meringolo et al. [89] consider the filtering of unwanted acoustic noise using wavelet transforms with some success in the analysis of simple prototype problems as well as violent sloshing flows. Comparatively understudied in the community, there is still much scope here to combine effective post-processing with controlled diffusion to deliver results of real use in industry.

The considerable improvement in results when using particle shifting in addition to δ-SPH highlights the importance of particle distribution in the accuracy and stability of simulations. SPH is traditionally thought of as a Lagrangian method, but increasingly there is a realization that controlling distributions in some way (and thereby directly or indirectly deviating from a true Lagrangian description of particle movement) is often very beneficial for a wide range of situations. Some recent advances in this area are reviewed in the next section.

10. Shifting and arbitrary Lagrangian–Eulerian approaches

Particle shifting is an increasingly standard technique used to prevent particle clustering and/or disordered particle arrangements, which can cause numerical instability. The notion of ‘shifting’ particles off their streamlines in order to improve distribution was first considered in SPH by Xu et al. [24], who based their formulation on the work of Nestor et al. [90] for the finite-volume particle method. This was reformulated by Lind et al. [20] to include free surfaces with a basis closely resembling the physics of Fick's law of diffusion. As the flow field evolves, SPH particles are displaced slightly into a more uniformly distributed configuration. A number of techniques (e.g. [91,92]) have been developed since, based on Xu et al. [24] or the diffusion method of Lind et al. [20].

In the majority of the particle-shifting techniques, a concentration gradient is used to evaluate a directional measure of particle disorder through

The concentration gradient is then used to diffuse the particles into a more uniform distribution following Fick's law,

where D_i is a numerical diffusion coefficient and Δx_s is a shifting displacement vector. From a stability analysis highlighted in [20,93], the maximum theoretical limit for the diffusion coefficient required to avoid instability is

The numerical diffusion coefficient can be rewritten as

where λ is a dimensionless parameter that defines the maximum shifting magnitude for a given concentration gradient per time step. From comparing (10.4) and (10.5), λ < 0.5 is required to avoid instability.

To improve the accuracy of the SPH approximation of derivatives for cases with disordered particle distributions, a correction for the kernel gradient was proposed by Bonet & Lok [32]. This is evaluated as

where ∇^cW_ij is the corrected kernel gradient and

These works effectively treated the particle shift as a treatment akin to re-meshing in conventional CFD, with shifted particles corrected for velocity, given that their position had now moved in an unchanged velocity field. In practice, this velocity correction was often not necessary as the shifting distance was very small compared with particle spacing, and so it became apparent quite quickly that the shifting of particles can be considered as being due to a non-Lagrangian transport velocity. Updating particles with a non-Lagrangian velocity is not a new idea, with X-SPH [94] using a particle transport velocity that is perturbed slightly by a weighted average of the velocities of surrounding particles. However, the approach based on Fick's law employed first gradients of the kernel, a measure that directly quantifies anisotropy in particle distribution, which particle shifting then aims to minimize.

Following Lind et al. [20], the incorporation of a transport velocity based on the first kernel gradient began to appear in SPH schemes. An example is that of Zhang et al. [87], which extends the transport velocity formulation of Adami et al. [95] to include a variable background pressure and localized smoothing lengths. However, to update particles with a non-Lagrangian velocity, care must be exercised if particle velocities are not subsequently corrected. Careful consideration of Reynolds transport theorem is, strictly speaking, required to accommodate potential mass changes that occur theoretically with non-Lagrangian particle movement. ALE schemes for SPH were first presented by Vila [34]; these schemes use conservation of fluxes between particles and have been widely used since, in conjunction with Riemann solvers, to remove pressure field fluctuations [96,97]. Oger et al. [98] presented a consistent ALE formulation with a transport velocity based on Fickian shifting. With the potential for mass to vary along a non-Lagrangian particle path, this formulation was consistent with the Reynolds transport theorem, improved particle distributions and thereby gave markedly improved results over the equivalent fully Lagrangian WCSPH scheme. Nevertheless, even if one does not use a formal ALE scheme, as the shifting distance/velocities are intentionally small, very good results can still be achieved. For example, recent work by Sun et al. [83] has shown that, through a quasi-Lagrangian scheme, the results remain accurate for a range of test cases, including Taylor–Green vortices and shallow water sloshing problems. Here, ‘quasi-Lagrangian’ means that additional transport velocity terms are included in the continuity equation, but, unlike a full ALE scheme, the mass per particle is assumed to remain approximately constant.

The importance of particle distribution for accuracy and stability is undeniable. While all the above studies build on the shifting approach of Lind et al. [20] and employ a single ‘shift’ per particle per time step, in recent years several works have seen benefit in applying the shifting procedure more frequently to further improve distributions. Litvinov et al. [39] applied a shifting procedure iteratively in a periodic domain to demonstrate that the particle distributions could be achieved that allow high accuracy and consistency, with near-ideal convergence of the sort only previously seen on a uniform Cartesian grid. However, this work was limited in that it remained quite idealized and without application to real fluid flow. The work of Vacondio & Rogers [44] from the 12th SPHERIC Workshop proceedings provided the initial framework for iterative shifting that is useful for practical fluid flow problems, with multiple shifting steps undertaken per time step to improve flow accuracy. Khayyer et al. [99] have since applied this iterative shifting scheme to multi-phase flows in high-density ratios with good results. This work also employed the free-surface shifting treatment of Khayyer et al. [92] to replace the free-surface shifting parameters as defined in Lind et al. [20].

Uniform or Cartesian particle distributions provide the most direct means of achieving high accuracy and convergence. Shifting, and to a greater extent iterative shifting, tries to recover this accuracy within a traditional SPH framework where one expects particles to move. Given the ALE context, we are motivated to re-think the meaning of an SPH particle, and it becomes reasonable to ask whether one should move SPH particles at all. In internal single-phase flows, for example, there are strong arguments for particles to remain fixed in controlled or pre-defined arrangements for which accuracy is potentially very high. Indeed, it was on fixed particles that the first exploration of high-order (above second) SPH was conducted [40] with accuracies demonstrated that rivalled spectral methods for practical fluid flow problems. The potential benefits of high-order SPH, including computational savings as well as increased accuracy, are considerable, and, as discussed above, it is an active area of current research. As well as allowing high-order convergence, viewing SPH particles in an ALE context as transitional entities between a fixed interpolation point and a material element of constant mass provides a clearer route for the development of rigorous coupling schemes with traditional grid-based methods. As with high-order SPH implementations, there are significant benefits to be had here, not least with regard to improved computation times. In the next section, recent approaches in the coupling of SPH with other methods are discussed.

11. Coupling with other methods

For many practical engineering problems, fluid flows often interact with solid objects and structures which themselves may deform according to a governing equation for the solid. Of course, fully coupled fluid–structure interaction (FSI) is the more general physical problem (rather than fluid interaction with perfectly rigid boundaries, say), and FSI is often crucial for problems in naval hydrodynamics, offshore and coastal engineering and aerospace, to name only a few areas. For FSI, the coupling of SPH with other numerical methods can provide significant benefit. Usually, the aim of the coupling is increased fidelity of the physical model, particularly for the solid, using a method or code with a well-tested modelling capability for the material and dynamics in question. A desire for computational savings, given that stiff elastic solids require small time steps, for example, is also an important factor. In the case of multi-scale fluid simulations, where SPH may be coupled with another fluid solver, computational time becomes the key motivation, with SPH typically being the more expensive method and so restricted to as small a region of the domain as possible, where its use is most advantageous, e.g. in regions of high free-surface deformation.

For FSI, the coupling of SPH with the FEM is a popular choice, probably because of the considerable success of the latter in modelling solid mechanics problems. Fourey et al. [100] considered the coupling of SPH (for the fluid) with an FEM (for the structure) applied to violent FSI problems. Code Aster is used for the finite-element (FE) scheme, coupled in a weak sense to an SPH code, and studied problems include dam breaks through an elastic gate. In a subsequent paper, greater attention is given to energy considerations in the SPH–FE coupling, with interfacial energy compared and different structural models considered [101]. This work has since been extended to a full 3D SPH–FE coupling applied to complex applications, including tyre hydroplaning on rough ground, with good agreement shown against available analytical and experimental solutions [101]. Other recent popular couplings have involved using physics engines (such as Project Chrono) for moving structures with multiple members and SPH for the fluid [102].

SPH–FE coupling for fluid–fluid multi-scale flows has also been investigated in recent years. Often employed within coastal/offshore or naval engineering, SPH is typically used in wave-breaking regions with a potential flow solution employed externally, in regions where deformation is smaller and the flow is likely to remain singly connected. For example, Fourtakas et al. [103] undertook a coupling of an ISPH algorithm with the potential flow FEM QALE-FEM, with accuracy maintained and significant computational savings achieved compared with ISPH-only formulations. For problems in naval hydrodynamics, Serván-Camas et al. [104] coupled SPH with a far-field FEM solver. Although the far-field solution provided diffraction–radiation solutions, rather than a full description of the potential flow, experimental agreement with the coupled model was good.

A key issue in coupling SPH solvers with potential flow solvers is reconciling the irrotational and inviscid nature of the potential solver as it transitions into SPH (and vice versa). Typically, this is undertaken using a numerical buffer zone (as in [103,105,106]), and, while effective, it is understandable that many would wish to remove this inconsistency and solve the full Navier–Stokes equations on a mesh for certain cases. Coupling with the FVM is a popular choice here; finite volume having long been a stalwart underpinning many conventional CFD software packages. Chiron et al. [51] couple SPH with an FVM, and take care to include several important features, including correct net mass transfer between SPH and finite-volume domains, as well as accurate transport of vorticity and free surfaces. This work directly built on that of Marrone et al. [107], where the coupling was comparatively simple and the finite-volume solver only applied in internal flow and solid boundary regions, away from free surfaces. Similarly, Neuhauser & Marongiu [108] used an overlapping zone of SPH and FVM to interpolate data between the two methods, applying it to Taylor–Green vortices and flow past a NACA0012 aerofoil.

While we have so far mentioned SPH coupled with mesh-based approaches, it is of course possible that SPH can be coupled with particle methods of various types—including itself. Total Lagrangian particle formulations for the structure are sensible as this approach alleviates known numerical instabilities in computational solid mechanics. Han & Hu [109] considered an SPH FSI model, employing SPH with a transport velocity formulation for the fluid, combined with a total Lagrangian formulation for the structure—the coupling here is, of course, made relatively straightforward as both governing equations are solved within the SPH framework. A fluid–structure dam break problem has been considered by Sun et al. [83] that similarly combines SPH with a total Lagrangian particle method for the structure. The fluid phase also includes multi-phase physics and adaptive particle resolution, making the scheme quite attractive and capable of high material deformations in both fluid and solid. Falahaty et al. [110] considered nonlinear elastic structure–fluid interaction, using an incompressible particle method for the fluid and stress point integration for the structure based on a moving least-squares discretization. This coupling approach outperforms an equivalent scheme that employed nodal integration (rather than stress point) over a wide range of test cases. Coupling SPH with SPH for fluid–solid interaction offers a number of benefits in terms of code development, but also physical and mathematical consistency regarding interpolation error and conservation of quantities such as mass and momentum. Khayyer et al. [111] successfully linked ISPH for the fluid with elastic SPH for the solid, considering a wide range of test cases including collapse of a water column on an elastic plate, dam breaks with an elastic gate, sloshing with baffles and violent hydro-elastic slam.

Of course, one may couple two different SPH schemes in the modelling of fluid–fluid, rather than fluid–structure, problems, especially in cases where the material characteristics of the fluids are distinct enough to warrant separate governing equations and/or discretization techniques. One such example is the work of Lind et al. [67], who couple a compressible SPH fluid solver to an ISPH fluid solver for water-compressible air flows. Advantages include physically consistent compressibility in both phases and retention of a sharp discontinuous material interface. Of course, such an approach can be considered a subset of those approaches available for multi-phase flows in SPH more generally.

Both coupled SPH and high-order SPH can offer significant gains in computation time, and will generally require far fewer floating point operations and/or less computer memory to reach a desired error than low-order Lagrangian SPH models of the same problem. However, as indicated above in the discussions around turbulence modelling, as problems increase in size and complexity, parallelization still becomes a necessity, especially in three dimensions. Accordingly, in the next section, recent advances in the run-time acceleration of the SPH method are discussed.

12. Run-time acceleration

Traditionally, SPH has been considered to be a computationally expensive method. This is primarily due to two factors. (i) Around each particle, the kernel support has a far larger number of neighbours than the equivalent stencil for mesh-based schemes. In two dimensions, there are typically 20–50 neighbours, and in three dimensions, there are 120–400 neighbours depending on the ratio of h to particle size dp. (ii) For hydrodynamics simulations, the time step is limited by the speed of sound in WCSPH and by the maximum velocity in ISPH. This leads to time steps of the order of 10⁻⁵ s or smaller. So, simulations typically require 1 million time steps for 1 s of physical time. Most applications in three dimensions require 10–100 million particles. The computational expense to run an SPH simulation must, therefore, not be underestimated and has required significant development to accelerate simulations to be practicable.

As a particle method with a compact support around each particle, SPH appears to lend itself naturally to parallelization where vector lists of particle neighbours and particle interactions can be constructed in the code and executed on different parallel architectures. However, accelerating simulations has presented a considerable challenge to the SPH community, in particular creating SPH codes that are portable across the range of heterogeneous architectures we now have at our disposal. The main efforts have focused on two areas: (i) developing massively parallel SPH codes to run on supercomputing clusters comprising thousands of compute nodes with each node containing multiple cores and (ii) exploiting emerging hardware such as GPUs, which were developed for a different purpose of high-speed graphics manipulation. A further challenge has been accelerating the solution of SPH formulations that involve solving large sparse matrices such as those encountered with the pressure Poisson equation (PPE).

13. Smoothed particle hydrodynamics convergence in practice

Following the developments discussed above, in this section some practical demonstrations of convergence in SPH are shown. This is mainly in the context of ISPH to avoid possible contamination from spurious pressure oscillations in WCSPH. As mentioned, Lagrangian SPH typically converges between first and second order for commonly used values of particle spacing, and this is quite consistent across various SPH methods and applications. This behaviour is clearly demonstrated, for example, in the ISPH simulations of Chow et al. [121] for the case of the impulsively moving plate (figure 2). This is quite a challenging test case, not only because of the impulsive nature of the dynamics but also because the flow is inviscid (but stabilized through particle shifting). Though challenging, the impulsive plate test is useful as it provides a free-surface analytical solution via a perturbation analysis allowing well-defined convergence studies to be undertaken. In the absence of analytical solutions, one may be able to work with experimental data or the finest resolution simulation available to get some indication of convergence. In these cases, however, convergence rates may not be constant or even monotonic. Of course, the experimental data against which comparisons are made are invariably subject to measurement errors of various types, and the uniform convergence often seen when comparing against analytical solutions should not be expected. This sort of behaviour is clearly shown in figure 3 for the case of air–water plate-wave slam [122], where convergence to experimental peak pressure values on the impacting plate is demonstrated, but with the error hovering above and below the first-order convergence line.

Ultimately, first-order discretization error terms will always arise as a result of particle non-uniformity when using a fully Lagrangian particle update. Trying to reduce this error motivated the exploration of Eulerian and ALE-SPH schemes, and, as previously mentioned, this has led to promising high-order SPH formulations. In Eulerian form, Lind & Stansby [40] showed near-perfect second-, fourth- and sixth-order convergence for velocity and pressure for solutions of the full unsteady Navier–Stokes equations when considering the Taylor–Green vortex test case using Cartesian particle distributions (see, for example, figure 4). This test case was somewhat idealized in being periodic and without solid walls, which, at the time, were only available to be at best second-order accurate and, therefore, disruptive to ideal high-order convergence. This restriction is being resolved, however, with promising high-order boundary conditions being developed that are able to model internal wall-bounded flows, also in Eulerian form [123]. This new approach provides high-order approximations to no-slip and no-flux boundary conditions allowing globally high-order solutions to the Navier–Stokes equations (see, for example, the Taylor–Couette flow and convergence plots shown in figure 5). High-order convergence using ALE-SPH schemes is the natural next step following Eulerian SPH, and, as mentioned, research in this area continues apace. Currently, this is the clearest pathway towards truly converged Lagrangian flow solutions.

14. Conclusion

The aim of SPH is to provide a general capability for complex incompressible or compressible flows that may be in complex domains, with free surfaces, interfaces and multiple phases. Importantly, flow solutions must be converged and of high accuracy. Recalling the list of perceived weaknesses in the Introduction, various observations on progress towards this goal may be drawn from this review.

• (i)

  The convergence rate of SPH is the subject of much research; rates of approximately unity are now common and the standard theoretical limit of 2 is increasingly possible with particle distribution regularization, known as shifting. However with recent application of high-order kernels rates of 4 and higher are becoming available for internal flows. The target is to achieve such convergence rates in the presence of free surfaces and other complex physics. Recent work has started to provide rigorous proofs on the nature of convergence.

• (ii)

  Solid boundary conditions have improved beyond convenient dummy (fixed) solid particles and more accurate mirror particles. This is an active research area and a range of options is opening up.

• (iii)

  Adaptivity with particles of constant mass is automatic but variable particle mass will reduce particle numbers further with larger particles for areas of small flow gradients. This is particularly important for turbulence modelling with small-scale structures to be resolved, for example, near boundaries.

• (iv)

  Spurious pressure fluctuations in the weakly compressible form are largely overcome either by using a Riemann solver or by inclusion of a numerical diffusion term in the continuity equation (δ-SPH) and particle shifting.

• (v)

  Particle regularization through shifting works well by preventing voids and unphysical instabilities, and is developing.

• (vi)

  SPH is well suited for parallel processing, although the ISPH form requires a pressure Poisson solver. The use of GPUs is increasing, rapidly providing fast and inexpensive accessible parallel processing. Coupling is becoming an increasingly viable option for reducing computational cost and adding new physics, even when SPH is coupled with itself. For example, coupling ISPH for water with compressible SPH for air is enabling realistic speeds of sound in both phases.

Ultimately, the aim is that with converged and high-order solutions errors in a simulation will only be due to the physical assumptions, for example the stress model in turbulence or non-Newtonian fluids, not due to numerical discretization.

Data accessibility

This article has no additional data.

Authors' contributions

S.J.L., B.D.R. and P.K.S. contributed equally and gave final approval for publication and agree to be held accountable for the work performed herein.

Competing interests

We declare we have no competing interests.

Funding

Support of EPSRC grant nos. EP/H018638/1, EP/L014890/1, EP/R005729/1 and the Royal Academy of Engineering and Leverhulme Trust Senior Research Fellowship (LTSRF1718∖14∖37) for D.R. is gratefully acknowledged.