1 Introduction

Smoothed particle hydrodynamics (SPH) is a numerical method originally developed for astrophysical simulations [1, 2] which has been applied to several different fluid dynamics phenomena. Thanks to its meshless and Lagrangian character, the method is suited to simulate topologically complex moving interfaces and large (nonlinear) deformations. This includes many applications in hydraulics involving complex highly dynamic free surfaces [3,4,5]. However, despite its potentially attractive features, SPH schemes present some relevant limitations, in comparison with classical Eulerian mesh-based numerical methods. To this end, the SPH rEsearch and engineeRing International Community (SPHERIC) has identified the current weaknesses that prevent widespread use of it in industry, classifying them into five categories or Grand Challenges (GCs) [6]. In this frame of reference, this paper aims to address GC1, which concerns the low accuracy and poor convergence rate of the method. Indeed, the SPH approximation of differential operators, differently from classical Eulerian schemes such as finite difference or finite volume, presents two classes of approximation: one related to the so-called smoothing length and the other caused by the spatial discretisation. The former appears at continuous level, whereas the latter is due to the approximation of the convolution integral with a discrete summation. For a comprehensive analysis of the effects of those errors on the accuracy and convergence rate of the SPH approximation of differential operators, the reader might refer to [7, 8].

In practice, the original SPH scheme [3] exhibits convergence rates around first order (rather than the theoretical second order of the continuous SPH spatial operators). Historically, the problem was addressed by introducing corrected SPH operators which ensure polynomial consistency (initially zeroth and first order [9,10,11], later up to nth order [12,13,14]) regardless of the particle distribution. The popularity of these corrections was relative as they pose difficulties to the stability of simulations and, importantly, break the appealing global conservation properties that have characterised SPH since its inception [15, 16]. Alternative approaches have targeted particle distributions to achieve a reduction of the discretisation error of SPH spatial operators. For example, Lind and Stansby [17] demonstrated that high-order convergence can be achieved for SPH in Eulerian simulations with particles distributed in the vertices of a Cartesian grid. In this context, it is also worth to mention the recent development of the Local Anisotropic Basis Function Method (LABFM) [18] as an interesting alternative to SPH schemes thanks to its framework to generate local high-order difference operators for arbitrary node distributions.

For applications with free-surface flows, where the Lagrangian description is desired, the strategy of periodic reinitialisation (or remeshing) of particle positions [19] appeared first, though it resulted too diffusive. Later the shifting technique based on Fick’s law of diffusion was developed [20, 21], where particle positions are slightly shifted towards regions of less particle concentration at each time step obtaining a remarkable increase in accuracy. However, updating particle positions in an unchanged velocity field while keeping constant their masses is not consistent with the Reynolds transport theorem [22]. A recent approach to the problem lies in arbitrary Lagrangian–Eulerian (ALE-SPH) schemes [23] with transport velocities prescribed (basing on Fickian shifting) to regularise the particle distribution [22], and where particles interact exchanging fluxes (computed via the solution of Riemann problems). Contrary to the traditional weakly compressible (WC-SPH) formulation where artificial viscosity [3] and diffusive terms [24, 25] are required for stability, in ALE-SPH schemes the intrinsic numerical dissipation of the Riemann solver guarantees a smooth pressure field (avoiding the need to tune empirical coefficients). It is remarkable that these schemes are consistent with the Reynolds transport theorem and, interestingly, they conserve global mass and momentum as well. Additional gains in the regularity of the particle distribution can be obtained with implicit shifting approaches such as the iterative method of Rastelli et al. [26, 27]. ALE-SPH schemes are, however, affected by excessive diffusion [23]. Avesani et al. [28] addressed this issue by introducing a formulation to compute the fluxes with a Riemann solver based on a Weighted Essentially Non-Oscillatory (WENO) reconstruction. Subsequently, an Arbitrary Derivative in Space and Time (ADER) integrator has been added [29] to skip the expensive reconstruction process at intermediate stages of a time integration step. More recently, Antona et al. [30] improved the efficiency of the WENO reconstruction in SPH by adopting corrected SPH operators [12] rather than moving least squares (MLS) [31] to markedly speed up the generation of the reconstruction polynomials. Very recently Vergnaud et al. [32] proposed a different approach based on a 1-D WENO reconstruction of the conserved variables at the particle–particle interface with the aim of avoiding the use of 2-D stencil in the WENO spatial reconstruction. However, despite all efforts in [28,29,30] dedicated to reduce diffusion in the ALE-SPH formulation, the analysis of the global error showed that the convergence rate of the scheme is still lower than second order and experiences a saturation effect.

Motivated by the above, this paper presents a novel ALE-SPH scheme where fluxes are computed by an approximate Riemann solver starting from a WENO polynomial reconstruction and two independent approaches have been included to improve the consistency of the SPH approximation of the divergence operators. The first one is obtained by applying a kernel correction procedure, and the second one is added by using a formulation for the transport velocity based on the explicit Fickian shifting. It is important to remark that, in this work, it has been preferred to study the effect of the latter two features separately to better assess the effects obtained by each one. The objective of the present work is to achieve high-order (above second) convergence in Lagrangian SPH simulations, particularly for test cases where highly anisotropic particle distributions develop easily. As suggested in [30], it is expected that improving the consistency of the divergence operators will decrease the associated discretisation error and, consequently, unleash the high-order capabilities of the WENO reconstruction. Finally, the improvement in accuracy will also be discussed in the presence of shock waves, generated by strong initial discontinuities in the pressure field.

The rest of the paper is organised as follows. It begins with a discussion of the basic methodology for ALE-SPH based on Riemann solvers in Sect. 2.1. The algorithm to compute SPH interpolations that guarantee the desired level of consistency is detailed in Sect. 2.2. Polynomial WENO reconstructions on sets of disordered particles are explained in Sect. 2.3. Details on the particle regularisation algorithm in ALE-SPH basing on Fickian shifting are given in Sect. 2.4. Section 3.1 tests ALE-SPH with the WENO reconstruction for the propagation of a small pressure perturbation in a fluid initially at rest, analysing the benefits of the scheme in terms of stability and low numerical diffusion. Then, Sect. 3.2 discusses accuracy and convergence for a 2-D vortex case followed by Sect. 3.3 that shows the robustness of the algorithm in the presence of strong discontinuities of a circular blast wave, and finally a standing wave has been simulated in Sect. 3.4 comparing the numerical dissipation obtained with the WENO and the MUSCL spatial reconstructions. Section 4 discusses the main results obtained and draws some conclusions.

2 Methodology

In this section, the high-order weakly compressible Lagrangian SPH scheme used in this work is explained including details on its numerically stable and low-dissipative formulation based on Riemann solvers; on the mesh-free WENO reconstruction adopted; and on the consistency enhancements introduced by using either corrected kernel interpolations or a particle regularisation algorithm in an ALE framework.

2.1 Arbitrary Lagrangian–Eulerian weakly compressible SPH based on Riemann solvers

The basic elements of the SPH formalism adopted in this paper, firstly established by Vila [23], can be reduced to the use of an ALE framework to lay out the equations governing the fluid motion and to the introduction of Riemann solvers to model the interactions between particles.

Following then the ALE formalism, the Euler equations can be written as:

$$\begin{aligned} \frac{\partial \textbf{Q}}{\partial t} + \nabla \cdot \textbf{F} = \textbf{0}, \end{aligned}$$
(1)

where t is time, and the conserved quantity \(\textbf{Q}\) and the associated flux \(\textbf{F}\) are a vector and a matrix, respectively, with the following expressions:

$$\begin{aligned} \textbf{Q} = \begin{pmatrix} \rho \\ \rho \textbf{v} \end{pmatrix}, \qquad \textbf{F} = \begin{pmatrix} \rho \textbf{v} \\ \rho \textbf{v} \otimes \textbf{v} + p \textbf{I} \end{pmatrix}. \end{aligned}$$
(2)

Here, \(\textbf{v}\) is fluid velocity, \(\rho \) is density, p is pressure and \(\textbf{I}\) represents the unit tensor. Note that, for the purposes of this paper, no source term originated from body forces (e.g. gravity) has been included in Eq. (1). This equation emphasises the conservation of some fluid properties in a control volume fixed in space. Aiming at building an expression to track the changes to the fluid state in a moving (not necessarily with the fluid velocity) frame of reference, Eq. (1) can be integrated on a generic volume \(\Omega \left( t \right) \) that moves and deforms in time:

$$\begin{aligned} \int _{\Omega \left( t \right) } \left( \frac{\partial \textbf{Q}}{\partial t} + \nabla \cdot \textbf{F} \right) \,\textrm{d}V = \textbf{0}, \end{aligned}$$
(3)

and then transformed on the basis of the Reynolds theorem to obtain:

$$\begin{aligned} \frac{d}{dt} \int _{\Omega \left( t \right) } \textbf{Q} \,\textrm{d}V - \oint _{\partial \Omega \left( t \right) } \left( \textbf{v}_0 \cdot \textbf{n} \right) \textbf{Q} \,dA + \int _{\Omega \left( t \right) } \nabla \cdot \textbf{F} \,\textrm{d}V = \textbf{0},\nonumber \\ \end{aligned}$$
(4)

where V is the volume and A is the area.

The second term in the expression above corresponds to a surface integral over the closed boundary \(\partial \Omega \left( t \right) \) of the compact volume \(\Omega \left( t \right) \). Here, \(\textbf{n}\) is the outward-pointing unit normal at each point on the boundary \(\partial \Omega \left( t \right) \), and \(\textbf{v}_0\) is the velocity of the corresponding area element (not the flow velocity). It is possible to express this surface integral in terms of an integral over the volume enclosed \(\Omega \left( t \right) \) if one considers first the identity:

$$\begin{aligned} \left( \textbf{b} \cdot \textbf{n} \right) \textbf{a} = \left( \textbf{a} \otimes \textbf{b} \right) \cdot \textbf{n}, \end{aligned}$$
(5)

where a and b are two generic vectors, then applies the divergence theorem. After regrouping terms, the following integral expression for the changes in the moving frame of reference is finally obtained:

$$\begin{aligned} \frac{d}{dt} \int _{\Omega \left( t \right) } \textbf{Q} \,\textrm{d}V = - \int _{\Omega \left( t \right) } \nabla \cdot \left( \textbf{F} - \textbf{Q} \otimes \textbf{v}_0 \right) \,\textrm{d}V \end{aligned}$$
(6)

On the other hand, Eq. (2) shows that the continuity and momentum conservation have been considered. In this paper, the system of equations is closed thanks the following equation of state:

$$\begin{aligned} p = c_0^2 \left( \rho - \rho _0 \right) , \end{aligned}$$
(7)

where p is pressure, \(c_0\) is the speed of sound, and \(\rho _0\) is a reference density.

The discretisation of the scheme is attained by first partitioning the problem domain into a finite set of N computational nodes or particles, \({\mathcal {P}}_i \left( 1 \le i \le N \right) \), located at positions \(\textbf{r}_i\) and with volumes \(V_i\). Any discrete quantity \(\textbf{g}_i \left( t \right) \) associated with the particle \({\mathcal {P}}_i\) is then computed as a space average over the volume of the particle:

$$\begin{aligned} \textbf{g}_i \left( t \right) = \frac{1}{V_i \left( t \right) } \int _{V_i \left( t \right) } \textbf{g} \left( \textbf{r}, t \right) \,\textrm{d}V. \end{aligned}$$
(8)

The semi-discrete form of Eq. (6) results consequently as:

$$\begin{aligned} \frac{d}{dt} \left( V_i \textbf{Q}_i \right) = - V_i \left[ \nabla \cdot \left( \textbf{F} - \textbf{Q} \otimes \textbf{v}_0 \right) \right] _i \end{aligned}$$
(9)

Note that, in the ALE-SPH scheme, the transport velocity \(\textbf{v}_{0i}\) associated with the trajectory described by the particle \({\mathcal {P}}_i\) differs in general from the corresponding fluid velocity \(\textbf{v}_i\). To predict the resulting change in the volume of the particle, it should be recalled the (purely geometric) relation between the divergence of a velocity field and the associated volumetric dilation rate:

$$\begin{aligned} \frac{\textrm{d}V_i}{\textrm{d}t} = V_i \left[ \nabla \cdot \textbf{v}_0 \right] _i. \end{aligned}$$
(10)

To complete the spatial discretisation in Eqs. (9) and (10), SPH summations are needed to approximate the divergence operators. Using standard zeroth-order consistent SPH approximations [33], the following expressions are obtained:

$$\begin{aligned}{} & {} \frac{d}{dt} \left( V_i \textbf{Q}_i \right) = - V_i \sum _j \left( \textbf{F}_i - \textbf{Q}_i \otimes \textbf{v}_{0i} + \textbf{F}_j - \textbf{Q}_j \otimes \textbf{v}_{0j} \right) \nonumber \\{} & {} \quad \cdot \nabla W_{ij} V_j, \end{aligned}$$
(11)
$$\begin{aligned}{} & {} \frac{\textrm{d}V_i}{\textrm{d}t} = V_i \sum _j \left( \textbf{v}_{0j} - \textbf{v}_{0i} \right) \cdot \nabla W_{ij} V_j, \end{aligned}$$
(12)

where \(W_{ij} = W \left( \textbf{r}_i - \textbf{r}_j, h \right) \) is the kernel function and h is the smoothing length that determines the (finite) size of the support domain of particle \({\mathcal {P}}_i\). In this paper, a computationally efficient cubic spline [15] has been chosen as kernel function:

$$\begin{aligned} W \left( \textbf{r}, h \right) = \frac{\kappa }{h^d} \left\{ \begin{matrix} 1 - \frac{3}{2} q^2 + \frac{3}{4} q^3, &{} 0 \le q \le 1, \\ \frac{1}{4} \left( 2 - q \right) ^3, &{} 1 \le q \le 2, \\ 0, &{} q \ge 2. \end{matrix} \right. \end{aligned}$$
(13)

Here, q is defined as the ratio \(\frac{\Vert \textbf{r} \Vert }{h}\), d is the dimension of the space, and \(\kappa \) is a normalisation factor (\(\frac{10}{7 \pi }\) in 2-D). For radially symmetric kernel functions, it can be easily demonstrated that \(\nabla \textbf{W}_{ij} = \psi \left( q \right) \textbf{n}_{ij}\) with \(\textbf{n}_{ij} = \frac{\textbf{r}_i - \textbf{r}_j}{\Vert \textbf{r}_i - \textbf{r}_j \Vert }\). Note then the symmetry of the form used for the SPH average in Eq. (11) which ensures global mass and momentum conservation of the scheme.

The SPH method is in general explicitly integrated in time, and therefore, the time step \(\Delta t\) is computed accordingly to the following condition:

$$\begin{aligned} \Delta t = CFL \min \left( \Delta t_f, \Delta t_c \right) \end{aligned}$$
(14)

with

$$\begin{aligned} \begin{aligned} \Delta t_f&= \min _i { \sqrt{ \frac{m_i h}{\Vert \frac{d \left( m_i \textbf{v}_i \right) }{dt}\Vert } } } \\ \Delta t_c&= \min _i { \frac{h}{\max \left( c_0, \, 10 \Vert \textbf{v}_i \Vert \right) } } \end{aligned} \end{aligned}$$

where \(m_i\) is the mass of particle i (obtained as the product \(V_i \rho _i\)) and CFL is the Courant number.

Being as well approximately centred in space, SPH schemes require some algorithmic artefact to ensure numerical stability. In this paper, it has been chosen to introduce some upwinding via the use of 1-D finite difference fluxes [23]. Similarly to mesh-based Godunov methods, each interaction between the central particle \({\mathcal {P}}_i\) and a neighbouring particle \({\mathcal {P}}_j\) in Eq. (11), \(\textbf{F}_i - \textbf{Q}_i \otimes \textbf{v}_{0i} + \textbf{F}_j - \textbf{Q}_j \otimes \textbf{v}_{0j}\), is interpreted as a centred approximation of the numerical flux associated with a Riemann problem posed in the midpoint (\(\bar{\textbf{r}} = \frac{\textbf{r}_i + \textbf{r}_j}{2}\)), along the direction \(\textbf{n}_{ji}\), and moving with velocity \(\bar{\textbf{v}}_{0ij} = \frac{\left( \textbf{v}_{0i} + \textbf{v}_{0j} \right) }{2}\). Naming such flux as \(\textbf{G}_{ij}\), Eq. (11) can be written as:

$$\begin{aligned} \frac{d \left( V_i \textbf{Q}_i \right) }{dt} = - V_i \sum _j 2 \textbf{G}_{ij} \cdot \nabla W_{ij} V_j. \end{aligned}$$
(15)

In the present work, a simple and computationally efficient Rusanov flux [34] is used to approximate the resolution of the moving Riemann problem:

$$\begin{aligned} \textbf{G}_{ij}= & {} \frac{1}{2} \left( \textbf{F} \left( \textbf{Q}_{ij}^L \right) - \textbf{Q}_{ij}^L \otimes \bar{\textbf{v}}_{0ij} + \textbf{F} \left( \textbf{Q}_{ij}^R \right) - \textbf{Q}_{ij}^R \otimes \bar{\textbf{v}}_{0ij} \right) \nonumber \\{} & {} - \frac{c_{ij}}{2} \left( \textbf{Q}_{ij}^R - \textbf{Q}_{ij}^L \right) \otimes \textbf{n}_{ji}, \end{aligned}$$
(16)

where \(\textbf{Q}_{ij}^L\) and \(\textbf{Q}_{ij}^R\) are the vectors of conserved quantities for the left and the right lateral states, and \(c_{ij}\) is the maximum value of speed of sound between the two particles that make up the pair. For the equation of state used in this work, Eq. (7), the field of speed of sound is actually both uniform in space and constant in time. The values of the vectors \(\textbf{Q}_{ij}^L\) and \(\textbf{Q}_{ij}^R\) for the lateral states are inferred from the fluid fields in the neighbourhood of the particles \({\mathcal {P}}_i\) and \({\mathcal {P}}_j\), respectively. Specifically, a high-order polynomial WENO reconstruction is here adopted to do such interpolation of the information coming from the surrounding fluid fields. (Further details are given in Sect. 2.3.) Note the second term in the right-hand side of Eq. (16) which naturally introduces the abovementioned upwinding (and the associated large stabilisation effect).

The system formed by Eqs. (12) and (15) is completed with the discrete version of Eq. (7):

$$\begin{aligned} p_i = c_0^2 \left( \rho _i - \rho _0 \right) \end{aligned}$$
(17)

and the expression needed to advance the positions of the particles:

$$\begin{aligned} \frac{d \textbf{r}_i}{dt} = \textbf{v}_{0i} \end{aligned}$$
(18)

These equations are integrated in time with an explicit second-order accurate symplectic scheme [16], where the size of the time step is bounded by Eq. (14).

A numerical value for the speed of sound \(c_0\) in Eq. (17) is chosen as only ten times the maximum expected fluid velocity in the simulation. This is a common assumption in so-called weakly compressible SPH schemes [3] that allows for much larger time steps at the expense of a small (\(\sim 1 \%\)) compressibility of the fluid.

2.2 High-order divergence operator

It is not uncommon for SPH practitioners (see, for example, [28, 29]) to use zeroth-order consistent SPH approximations for differential operators. In this work, the procedure devised by Liu and Liu [12] to build a higher-order consistent SPH divergence operator is used in Eqs. (12) and (15), and their effect in the global accuracy of the present scheme is studied.

The algorithm to increase the consistency of an SPH approximation is described here. For simplicity, the explanation is restricted to the case of guaranteeing first-order consistency of a two-dimensional function, though the process can be readily generalised for higher orders of consistency (and higher dimensions of the functional space).

A Taylor expansion of a scalar field g centred in the position of particle \({\mathcal {P}}_i\), \(\textbf{r}_i = \left( x_i, y_i \right) \), can be used to compute the values of the field at a generic location \(\textbf{r} = \left( x, y \right) \) in the vicinity with the following expression:

$$\begin{aligned} g \left( x, y \right)= & {} g_i + \left( x - x_i \right) \partial _x g_i + \left( y - y_i \right) \partial _y g_i \nonumber \\{} & {} + O \left( \Vert \textbf{r} - \textbf{r}_i \Vert ^ 2 \right) , \end{aligned}$$
(19)

where \(g_i\), \(\partial _x g_i\), and \(\partial _y g_i\) are the values of the scalar field and of its first derivatives (\(\frac{\partial g}{\partial x}\) and \(\frac{\partial g}{\partial y}\)) at the location of particle \({\mathcal {P}}_i\). Evaluating the expansion at the location of a neighbouring particle \({\mathcal {P}}_j\), and neglecting second- and higher-order terms, the expression takes the form:

$$\begin{aligned} g_j = g_i + x_{ji} \, \partial _x g_i + y_{ji} \, \partial _y g_i. \end{aligned}$$
(20)

Here, \(g_j\) is the value of the scalar field at the location of particle \({\mathcal {P}}_j\), and the notation \(x_{ji}\) and \(y_{ji}\) stands for \(x_j - x_i\) and \(y_j - y_i\), respectively. Multiplying Eq. (20) by the product of the value of the kernel function for the pair times the volume of particle \({\mathcal {P}}_j\), \(W_{ij}V_j\), and performing a summation for all the neighbouring particles, the following equation is obtained:

$$\begin{aligned} \sum _j g_j W_{ij} V_j= & {} g_i \sum _j W_{ij} V_j + \partial _x g_i \sum _j x_{ji} W_{ij} V_j \nonumber \\{} & {} + \partial _y g_i \sum _j y_{ji} W_{ij} V_j, \end{aligned}$$
(21)

where the unknowns are \(g_i\), \(\partial _x g_i\), and \(\partial _y g_i\). Two additional equations can be obtained from Eq. (20) taking similar steps but using instead the first derivatives of the kernel function, \(\partial _x W_{ij}\) and \(\partial _y W_{ij}\):

$$\begin{aligned}{} & {} \sum _j g_j \, \partial _x W_{ij} \, V_j \nonumber \\{} & {} \quad = g_i \sum _j \partial _x W_{ij} \, V_j + \partial _x g_i \sum _j x_{ji} \, \partial _x W_{ij} \, V_j \nonumber \\{} & {} \quad + \partial _y g_i \sum _j y_{ji} \, \partial _x W_{ij} \, V_j, \end{aligned}$$
(22)
$$\begin{aligned}{} & {} \sum _j g_j \, \partial _y W_{ij} \, V_j \nonumber \\{} & {} \quad = g_i \sum _j \partial _y W_{ij} \, V_j + \partial _x g_i \sum _j x_{ji} \, \partial _y W_{ij} \, V_j \nonumber \\{} & {} \quad + \partial _y g_i \sum _j y_{ji} \, \partial _y W_{ij} \, V_j. \end{aligned}$$
(23)

Therefore, Eqs. (21)–(23) constitute a linear system:

$$\begin{aligned} \textbf{A}_i \textbf{x}_i = \textbf{b}_i, \end{aligned}$$
(24)

where

$$\begin{aligned} \textbf{A}_i= & {} \begin{pmatrix} \sum _j W_{ij} V_j &{} \sum _j x_{ji} W_{ij} V_j &{} \sum _j y_{ji} W_{ij} V_j \\ \sum _j \partial _x W_{ij} \, V_j &{} \sum _j x_{ji} \, \partial _x W_{ij} \, V_j &{} \sum _j y_{ji} \, \partial _x W_{ij} \, V_j \\ \sum _j \partial _y W_{ij} \, V_j &{} \sum _j x_{ji} \partial _y W_{ij} \, V_j &{} \sum _j y_{ji} \, \partial _y W_{ij} \, V_j \end{pmatrix}, \end{aligned}$$
(25)
$$\begin{aligned} \textbf{b}_i= & {} \begin{pmatrix} \sum _j g_j W_{ij} V_j \\ \sum _j g_j \, \partial _x W_{ij} \, V_j \\ \sum _j g_j \, \partial _y W_{ij} \, V_j \end{pmatrix}, \end{aligned}$$
(26)

and the vector of unknowns is defined as:

$$\begin{aligned} \textbf{x}_i = \begin{pmatrix} g_i \\ \partial _x g_i \\ \partial _y g_i \end{pmatrix}. \end{aligned}$$
(27)

The linear system in Eq. (24) can be solved, for instance, with an LU decomposition algorithm [35]. Note that the method relies on the fact that \(\textbf{A}_i\) is in general diagonally dominant (provided there is near complete support) and hence it can be inverted [12, 30].

The procedure described is applied to compute the corrected SPH approximations of the two divergence operators that appear in Eqs. (12) and (15). A generic vector field in a two-dimensional space, \(\textbf{g} = \left( g_x, g_y \right) \), has first-order derivatives \(\left[ \partial _x g_{x} \right] _i\), \(\left[ \partial _y g_x \right] _i\), \(\left[ \partial _x g_y \right] _i\), \(\left[ \partial _y g_y \right] _i\) at the location of particle \({\mathcal {P}}_i\). Please note the notation \(\left[ \dots \right] _i\) to specify here properties at particle \({\mathcal {P}}_i\), for the sake of clarity. The corresponding value of the divergence operator at that same location is obtained from the definition of the differential operator:

$$\begin{aligned} \left[ \nabla \cdot \textbf{g} \right] _i = \left[ \partial _x g_x \right] _i + \left[ \partial _y g_y \right] _i. \end{aligned}$$
(28)

This approach is used in this paper to obtain high-order consistent versions of the divergence SPH operators. For example, Eq. (12) can be rewritten as:

$$\begin{aligned} \frac{\textrm{d}V_i}{\textrm{d}t} = V_i \left( \langle \left[ \partial _x v_{0x} \right] _i \rangle _C + \langle \left[ \partial _y v_{0y} \right] _i \rangle _C \right) \end{aligned}$$
(29)

where \(\langle \dots \rangle _C\) is the SPH-corrected spatial operator.

Then, by solving the linear system of Eq. (24) in which the generic function \(g_i\) in Eq. (27) is equal to the x-component of the transport velocity \(\left[ v_{0x} \right] _i\) for particle i, it is possible to obtain \(\langle \left[ \partial _x v_{0x} \right] _i \rangle _C\). Similarly, \(\langle \left[ \partial _y v_{0y} \right] _i \rangle _C\) can be obtained by solving the linear system in Eq. (24) a second time where \(g_i = \left[ v_{0y} \right] _i\). Please note that the matrix \(\textbf{A}_i\) does not depend on the scalar function \(g_i\) and therefore is unique for each particle.

The same procedure can also be adopted for any other divergence operator that appears in Eq. (15) (one for the continuity equation and two for each component of the momentum equation) and therefore all different SPH divergence operators adopted in the model can be linearly consistent.

Note that the generalisation to achieve a higher order of consistency implies firstly to retain higher-order terms in Eq. (20), and secondly to pose the necessary number of additional equations by repeating the procedure to obtain Eqs. (21)–(23) with the aid of higher derivatives of the kernel function as well. In the present work, a second-order consistent formulation is adopted.

2.3 High-order polynomial Weighted Essentially Non-Oscillatory Reconstruction

As previously mentioned in Sect. 2.1, in the present work the interaction between a pair of neighbouring particles \({\mathcal {P}}_i\) and \({\mathcal {P}}_j\) is modelled via a moving Riemann problem posed in the midpoint of the pair. To compute the lateral states (\(\textbf{Q}_{ij}^L\), \(\textbf{Q}_{ij}^R\)) of such Riemann problem, it is necessary, at each time step, to reconstruct smooth functions that represent the spatial variation of the solution within the support domain of particles \({\mathcal {P}}_i\) and \({\mathcal {P}}_j\); these reconstructions can then be evaluated at \(\bar{\textbf{r}} = \frac{\textbf{r}_i + \textbf{r}_j}{2}\). It is well known that, unless an accurate reconstruction procedure is used, any algorithm to compute the solution of the Riemann problem originates too much numerical diffusion, rendering the overall scheme very inaccurate. Particularly, it has been reported [36, 37] that a piece-wise constant approximation (where the lateral states are taken as those of the particles \({\mathcal {P}}_i\) and \({\mathcal {P}}_j\)) provides poor results. In this paper, a high-order nonlinear polynomial WENO reconstruction is used. Firstly proposed by Avesani et al. [28], this method is robust to approximate functions that have large gradients or discontinuities. The algorithm consists in three steps: construction of the candidate stencils; computation of a polynomial function per stencil; and nonlinear combination of the candidate polynomial functions. These are detailed in the subsequent paragraphs.

Construction of the candidate stencils Figure 1 shows a sketch of all the candidate stencils considered for a sample randomised two-dimensional particle distribution. A characteristic length \(h_\textrm{weno}\) is first chosen for the size of the support domain used in the reconstruction method. The main reconstruction stencil, known as the central stencil \({\mathcal {S}}_{i, 0}\), considers all the particles \({\mathcal {P}}_j\) at a distance sufficiently close to the centre of the reconstruction \({\mathcal {P}}_i\):

$$\begin{aligned} {\mathcal {S}}_{i, 0} = \left\{ {\mathcal {P}}_j :\left\| \textbf{r}_{ij} \right\| \le h_\textrm{weno} \right\} , \end{aligned}$$
(30)

where \(\textbf{r}_{ij} = \textbf{r}_i - \textbf{r}_j\). To confer directional (or adaptive) behaviour on the reconstruction method, several one-sided lateral stencils are defined. For two-dimensional problems, this paper uses the eight lateral stencils proposed in [28]. Here, a neighbouring particle \({\mathcal {P}}_j\) is assigned to a lateral stencil \({\mathcal {S}}_{i, s}\) (\(1 \le s \le 8\)) depending on the value of the angular coordinate \(\theta \in \left[ 0, 2\pi \right] \) of the vector \(\textbf{r}_{ji}\) in a polar coordinate system:

$$\begin{aligned}{} & {} {\mathcal {S}}_{i, s} = \left\{ {\mathcal {P}}_j :\left\| \textbf{r}_{ij} \right\| \le 2h_\textrm{weno}, \, \theta \in \left[ \left( s - 1 \right) \frac{\pi }{4}, s \frac{\pi }{4} \right] \right\} \nonumber \\{} & {} \quad \left( 1 \le s \le 8 \right) . \end{aligned}$$
(31)

Note that the radius of the circular sectors defining the envelope of the lateral stencils doubles the radius of the central one. This is needed to have a sufficient number of particles per stencil to carry out the second step of the reconstruction.

Fig. 1
figure 1

Sketch of the central and the one-sided lateral stencils for a sample randomised 2-D particle distribution

Full size image

Computation of a polynomial function per stencil For each candidate stencil \({\mathcal {S}}_{i,s}\) of the particle \({\mathcal {P}}_i\), a polynomial function \(\textbf{P}_{i,s}\) of a specified degree M is produced. These polynomials, defined in local coordinates \(\hat{\textbf{r}} = \left( \hat{x}, \hat{y} \right) = \frac{\textbf{r} - \textbf{r}_i}{h_\textrm{weno}}\), have the following form:

$$\begin{aligned} \textbf{P}_{i,s} \left( \hat{\textbf{r}} \right) = \textbf{Q}_i + \sum _{1 \le \alpha _1 + \alpha _2 \le M} \textbf{C}_{i, s}^{\alpha _1, \alpha _2} \; \hat{x}^{\alpha _1} \; \hat{y}^{\alpha _2}, \end{aligned}$$
(32)

where \(\alpha _1\) and \(\alpha _2\) are nonnegative integer indices, and \(\textbf{C}_{i,s}^{\alpha _1, \alpha _2}\) are the unknown polynomial coefficients that must be determined by interpolating the information provided by the corresponding stencil. Note that \(\hat{\textbf{r}}\) remains \(O \left( 1 \right) \) in the neighbourhood of \({\mathcal {P}}_i\) which avoids issues with computer precision irrespective of the method used for the interpolation. In [28, 29], it is proposed first to construct an overdetermined linear system of equations by forcing the polynomial to pass through all the points in the stencil, and then to solve it in a least-squares sense [38]. In this paper, it has been preferred to identify Eq. (32) with a truncated Taylor expansion of the field \(\textbf{Q}\) around the position of particle \({\mathcal {P}}_i\), which provides the following alternative expression for the unknown polynomial coefficients:

$$\begin{aligned} \textbf{C}_{i, s}^{\alpha _1, \alpha _2} = \frac{1}{\alpha _1! \alpha _2!} \left[ \frac{\partial ^{\alpha _1 + \alpha _2} \textbf{Q}}{\partial \hat{x}^{\alpha _1} \partial \hat{y}^{\alpha _2}} \right] _i. \end{aligned}$$
(33)

The spatial derivatives \(\left[ \frac{\partial ^{\alpha _1 + \alpha _2} \textbf{Q}}{\partial \hat{x}^{\alpha _1} \partial \hat{y}^{\alpha _2}} \right] _i\) can then be computed by means of consistent kernel interpolations following the procedure previously detailed in Sect. 2.2. (Note that the order of consistency chosen when assembling the matrix \(\textbf{A}_i\) in Eq. (24) must be equal to the order M of the polynomial reconstruction.) As demonstrated in [30], this approach provides higher accuracy and computational efficiency than the MLS method proposed in [28]. Finally, if a given stencil does not provide near complete support and the interpolation procedure fails as a consequence, the corresponding polynomial is ignored in the third step of the reconstruction. Note that the method relies on the fact that at least one candidate polynomial can be computed.

Nonlinear combination of the candidate polynomial functions The stencil-adaptive behaviour of the method is achieved by assigning a different weight \(\widetilde{\omega }_s\) to each of the candidate polynomials, which has the form:

$$\begin{aligned} \widetilde{\omega }_s = \frac{\lambda _{s}}{\left( \varepsilon + \sigma _s \right) ^4}. \end{aligned}$$
(34)

Here, the numerator provides a larger contribution to the central stencil by choosing \(\lambda _0 = 10^5\) and \(\lambda _s = 1\) for \(s \ge 1\). On the other hand, the smoothness indicator \(\sigma _s\) that appears in the denominator gets larger with the oscillatory character of the corresponding polynomial, hence penalising its contribution. As suggested in [28], this indicator is computed from the polynomial coefficients with the following expression:

$$\begin{aligned} \sigma _s = \sum _{1 \le \alpha _1 + \alpha _2 \le M} \left( \textbf{C}_{i, s}^{\alpha _1, \alpha _2} \right) ^2. \end{aligned}$$
(35)

Finally, the parameter \(\varepsilon \) is introduced in Eq. (34) to avoid null denominators. (Hence, it must be assigned a small value such as \(10^{-7}\) or \(10^{-14}\) for single or double precision, respectively.) Once the following normalised weights are computed:

$$\begin{aligned} \omega _s = \frac{\widetilde{\omega }_s}{\sum _r \widetilde{\omega }_r}, \end{aligned}$$
(36)

the final polynomial function that reconstructs the field of \(\textbf{Q}\) in the neighbourhood of particle \({\mathcal {P}}_i\) is obtained as a convex combination of all the candidate polynomials:

$$\begin{aligned} \textbf{P}_i \left( \hat{\textbf{r}} \right) = \sum _s \omega _s \; \textbf{P}_i^s \left( \hat{\textbf{r}} \right) . \end{aligned}$$
(37)

2.4 Particle regularisation in arbitrary Lagrangian–Eulerian frameworks

SPH schemes usually incorporate some sort of technique to slightly perturb the Lagrangian trajectories of the particles preventing the formation of anisotropies in the particle distribution [20, 21]. The strategy followed in this work exploits the ALE framework introduced in Sect. 2.1 where the transport velocity \(\textbf{v}_{0i}\) of particle \({\mathcal {P}}_i\) can differ from the corresponding fluid velocity \(\textbf{v}_i\). Oger et al. [22] demonstrated that this approach allows to preserve both conservation and consistency in a stable SPH scheme.

The key idea to maintain a relatively equispaced particle distribution as the simulation evolves consists in updating the transport field at the end of each time step according to the following expression:

$$\begin{aligned} \textbf{v}_{0i} = \textbf{v}_i + \delta \textbf{v}_i, \end{aligned}$$
(38)

where \(\delta \textbf{v}_i\) represents only a small perturbation to the Lagrangian velocity \(\textbf{v}_i\), and such that it induces a motion from regions of high concentration to regions of low concentration of particles. For a time step of size \(\Delta t\) (determined by the CFL condition of Eq. (14)), the perturbation in velocity \(\delta \textbf{v}_i\) can be written as the quotient

$$\begin{aligned} \delta \textbf{v}_i = \frac{\delta \textbf{r}_i}{\Delta t}, \end{aligned}$$
(39)

where the unknown \(\delta \textbf{r}_i\) is the so-called shifting displacement vector. As proposed by Lind et al. [20], \(\delta \textbf{r}_i\) can be computed with an explicit formula similar to Fick’s law of diffusion that provides the magnitude and direction of the position shift:

$$\begin{aligned} \delta \textbf{r}_i = - {\mathcal {D}} \nabla C_i. \end{aligned}$$
(40)

Here, \({\mathcal {D}}\) is a (strictly) numerical diffusion coefficient, and \(\nabla C_i\) is a measure of the gradient of particle concentration at the position of particle \({\mathcal {P}}_i\). Following [20], \(\nabla C_i\) is estimated basing on a standard SPH approximation of the gradient of partition of unity:

$$\begin{aligned} \nabla C_i = \sum _j V_j \left( 1 + f_{ij} \right) \nabla W_{ij}. \end{aligned}$$
(41)

The term \(f_{ij}\) is added to prevent vanishing values of \(\nabla C_i\) when two particles get extremely close to each other. The expression to compute \(f_{ij}\), firstly proposed in SPH literature to alleviate tensile instability issues [39], has the form:

$$\begin{aligned} f_{ij} = R \left( \frac{W_{ij}}{W \left( \Delta x \right) } \right) ^n, \end{aligned}$$
(42)

with parameters \(R = 0.2\) and \(n = 4\), and where \(\Delta x\) is the initial particle spacing. With regard to the diffusion coefficient in Eq. (40), the following form is used in this work:

$$\begin{aligned} {\mathcal {D}} = A_sh^2, \end{aligned}$$
(43)

where h is the smoothing length already defined in Sect. 2.1, and the dimensionless parameter \(A_s\) has an upper bound of 0.5 to avoid numerical instability [20, 21]. For the results in this paper, \(A_s = 0.5\) has been taken as suggested in [20]. Finally, it is known [20, 22] that Eqs. (40)–(43) do not guarantee that the magnitude of the shifting displacement vector \(\delta \textbf{r}_i\) obtained remains much smaller than the smoothing length. As suggested in [20], in the implementation used in this work, it has been chosen to cap the magnitude of \(\delta \textbf{r}_i\) to a 20% of h to prevent too large transport velocity that might generate stability issues.

3 Numerical results and discussion

This section tests the accuracy and convergence of the WENO SPH formulation enhanced with either corrected SPH operators for divergence terms as illustrated in Sect. 2.2 (hereafter referred as corrected WENO SPH), or with a transport field regularised with the explicit shifting algorithm illustrated in Sect. 2.4 (hereafter referred as regularised WENO SPH). For comparisons, the baseline polynomial WENO SPH of Sects. 2.1 and 2.3 will be referred simply as WENO SPH. Additionally, the WENO SPH scheme with Eulerian transport (i.e. with particle positions regularly spaced in a Cartesian grid and fixed in time) will be referred as Eulerian WENO SPH. Finally, the ALE-SPH scheme based on a piece-wise constant approximation (hereafter referred as piece-wise constant ALE-SPH) and on a MUSCL reconstruction (hereafter referred as MUSCL SPH) will be considered for completeness.

With this goal in mind, four different test cases have been chosen:

  1. 1.

    the propagation of a small pressure perturbation in a fluid initially at rest,

  2. 2.

    a two-dimensional vortex where strong anisotropies in the particle distribution appear easily, consequence of large fluid deformations,

  3. 3.

    a circular blast wave problem that tests the algorithm behaviour in the presence of pressure discontinuities, and

  4. 4.

    a standing wave which involves a dedicated treatment for free surfaces and a solid walls.

For the SPH scheme, the ratio of the smoothing length and the initial particle spacing remains fixed at value \(h / \Delta x = 2\), and the CFL coefficient is set to 0.2. The characteristic length for the WENO reconstruction is chosen the double of the smoothing length \(h_\textrm{weno} = 2h\), and second-order polynomials are used. Additionally, when corrected SPH interpolations are used to approximate a divergence operator, second-order consistency is enforced.

The computer code used to produce all the results in this section has been implemented basing on the DualSPHysics open-source project [40]: an efficient state-of-the-art SPH code for free-surface flows with capabilities for modern Graphics Processing Units (GPUs).

3.1 Propagation of a small pressure perturbation

A problem to test the propagation of a quasi-one-dimensional small pressure perturbation is constructed in 2-D by specifying a squared fluid domain with periodic boundary conditions in both directions and a region inside \(x \in \left[ x_1, x_2 \right] \) (with \(x_1 < x_2\)) initially at a slightly higher pressure than the rest. By virtue of Eq. (7), this condition can be expressed in densities as follows:

$$\begin{aligned} \rho \left( \textbf{x}, t=0 \right) = \left\{ \begin{matrix} \rho _H, &{} x_1 \le x \le x_2, \\ \rho _0, &{} \text {otherwise}, \end{matrix} \right. \end{aligned}$$
(44)

where \(\rho _0\) is the reference density previously introduced in Eq. (7), and \(\rho _H\) is the density in the high pressure region (\(\rho _H > \rho _0\)). Both regions are initially at rest, \(\textbf{v} = \textbf{0}\). The initial discontinuity splits into two smooth pressure waves travelling horizontally in opposite directions. In this paper, the values of the problem parameters are \(\rho _0 = 1.0\) kg/m2, \(\rho _H = 1.001\) kg/m2, \(x_1 = 1.1\) m, \(x_2 = 1.2\) m, and the sound speed is set to \(c_0 = 1.0\) m/s. The domain is a square box \(\left[ 0, 2 \right] \times \left[ 0, 2 \right] \) with periodic boundary conditions in both Cartesian coordinates [41].

For the presentation of the results, the following dimensionless variables (marked by a hat symbol, ⌃) are specified:

$$\begin{aligned} \begin{aligned} \hat{x}&= \frac{x}{D},&\hat{z}&= \frac{z}{D},&\hat{t}&= \frac{t \, c_0}{D},&\hat{p}&= \frac{p}{c_0^2 \left( \rho _H - \rho _0 \right) }, \end{aligned} \end{aligned}$$
(45)

where the characteristic length of the problem is defined as \(D = x_2 - x_1\). To assess the accuracy of the computations, the numerical solution obtained for a very fine resolution (\(D / \Delta x = 64\)) with the baseline WENO SPH scheme is used as reference. Note that in this particular test case the particles are initially arranged on a Cartesian grid and no significant particle displacement takes place. This means that there is no inconsistency of the SPH spatial operators generated by the particle disorder. All the subsequent results presented in this section correspond to time \(\hat{t} = 4\).

Figure 2 shows the pressure distribution in the \(\hat{z} = 10\) section computed with standard SPH [40] using two different values of artificial viscosity: \(\alpha = 0.1\) and \(\alpha = 2\). The resolution is \(D / \Delta x = 32\), and the reference solution is displayed superimposed in black. For the lower value of \(\alpha \), close to those typically seen in hydraulic engineering applications, the numerical solution with standard SPH presents a highly oscillatory character, specifically in between the two pressure waves. Increasing \(\alpha \) to 2 enhances numerical stability which helps to mitigate the oscillations. However, the size of the pressure oscillations around the location of the initial discontinuities remains similar (or even worse). Moreover, the associated numerical dissipation introduced, which translates into a net reduction of the magnitude of the pressure distribution maximum, is significant.

Fig. 2
figure 2

Pressure distribution for the small pressure perturbation test computed with standard SPH for two different values of artificial viscosity: \(\alpha = 0.1\) and \(\alpha = 2\)

Full size image

To show the benefits of using Riemann solvers for numerical stabilisation, Fig. 3 displays the same pressure profile computed now with the piece-wise constant ALE-SPH scheme and for different increasing resolutions. This scheme completely erases any spurious oscillation in the pressure field, also around the location of the initial discontinuities. In addition, it can observed that the numerical results converge to the reference solution. However, in the absence of (non-trivial) spatial reconstruction to compute the lateral states of the Riemann problem, the numerical dissipation associated with the Riemann solver is prohibitively large.

Fig. 3
figure 3

Pressure distribution for the small pressure perturbation test computed with the piece-wise constant ALE-SPH scheme for four different resolutions: \(D / \Delta x = 4, 8, 16, 32\)

Full size image

Results are greatly improved in terms of low numerical diffusion when using the high-order WENO spatial reconstruction, as shown in Fig. 4. Thanks to the WENO reconstruction, the scheme is able to resolve both the minimum and the maximum of the pressure profile, with convergent results for increasing resolutions.

Fig. 4
figure 4

Pressure distribution for the small pressure perturbation test computed with the WENO SPH scheme for four different resolutions: \(D / \Delta x = 4, 8, 16, 32\)

Full size image

In terms of computational cost, the ALE-SPH scheme in its current implementation is two times as expensive as the standard SPH scheme, while the WENO high-order formulation is about 10 times more expensive than the ALE-SPH scheme. It is important to note that these computational times are specific to the CPU implementation, and further work is needed to optimise the WENO high-order formulation for both CPUs and GPUs.

3.2 Two-dimensional vortex

The weakly compressible 2-D vortex introduced in [42] is an inviscid and stationary problem that presents circular streamlines. Considering polar coordinates \(\left( r, \theta \right) \) with origin in the centre of the vortex, the analytical solution for density \(\rho \) and tangential velocity \(v_{\theta }\) (the radial velocity \(v_r\) is null) depends only on the radial distance r and is given by:

$$\begin{aligned} \begin{aligned} \rho \left( r \right)&= \rho _0 + \left( \rho _M - \rho _0 \right) \left[ 1 - \frac{2r}{r_0} e^{- \frac{2r}{r_0}} - e^{- \frac{2r}{r_0}} \right] , \\ v_{\theta } \left( r \right)&= \frac{2r}{r_0} c_0 \sqrt{\frac{\rho _M - \rho _0}{\rho \left( r \right) } e^{- \frac{2r}{r_0}}}, \end{aligned} \end{aligned}$$
(46)

where \(\rho _0\) and \(c_0\) were defined when introducing Eq. (7), \(\rho _M\) is the asymptotic value of density when \(r \rightarrow \infty \), and \(r_0\) marks the radial distance at which the fluid (tangential) velocity is maximal. In this paper, the values set for these parameters are \(\rho _0\) = 1.0 kg/m2, \(c_0 = 1.0\) m/s, \(\rho _M = 1.01\) kg/m2, and \(r_0 = 0.1\) m. The simulation domain is a circle of radius \(R = 0.7\) m. To impose boundary conditions, the necessary number of rows of dummy particles around the simulation domain is added, with density and velocity calculated from Eq. (46). For simulations with a Lagrangian transport field (particles move with the local fluid velocity, \(\textbf{v}_{0i} = \textbf{v}_i\)), boundary particles change their momentum according to the centripetal force of the circular motion at constant speed that they describe. For the case of a Eulerian transport field (particles remain fixed in space, \(\textbf{v}_{0i} = \textbf{0}\)), all temporal derivatives for the boundary particles are zero as the problem is stationary.

For the presentation of the results, the following dimensionless variables (marked by a hat symbol, ⌃) are defined:

$$\begin{aligned} \begin{aligned} \hat{x}&= \frac{x}{r_0},&\hat{z}&= \frac{z}{r_0},&\hat{v}&= \frac{v}{c_0 \sqrt{\frac{\rho _M - \rho _0}{\rho _0}}}, \\ \hat{t}&= \frac{t \, c_0 \sqrt{\frac{\rho _M - \rho _0}{\rho _0}}}{r_0},&\hat{p}&= \frac{p}{c_0^2 \left( \rho _M - \rho _0 \right) }. \end{aligned} \end{aligned}$$
(47)

Plots displaying the spatial distribution of the approximation error of these magnitudes are provided. This error can be defined for a generic scalar property g of particle \({\mathcal {P}}_i\) as

$$\begin{aligned} err ( g )_i = \vert g_i^\textrm{exact} - g_i^\textrm{estimate} \vert . \end{aligned}$$
(48)

Moreover, convergence studies of the \(L_2\) error with resolution refinement are carried out. The \(L_2\) error for a generic scalar field g is computed by the expression:

$$\begin{aligned} \left\| err ( g ) \right\| _2 = \sqrt{\frac{1}{N} \sum _i^N \vert g_i^\textrm{exact} - g_i^\textrm{estimate} \vert ^2}. \end{aligned}$$
(49)

3.2.1 Results using corrected SPH operators for divergence terms

Figure 5 shows a comparison of the approximation error in pressures and horizontal velocities computed with the baseline WENO SPH and with the corrected WENO SPH schemes following the procedure described in Sect. 2.2. The measurement is taken at time \(\hat{t} = 2\), for a resolution \(r_0 / \Delta x = 28\), and errors are displayed superimposed to the distribution of particles. Note that the same error scales are used for both schemes to ease the comparison, and that results for the vertical velocity \(\hat{v}_z\) are omitted due to the similarity with those of the horizontal velocity \(\hat{v}_x\). The superior accuracy of the corrected WENO SPH scheme is apparent, especially for the pressure field. Note that the improvement in accuracy is obtained irrespective of the actual particle distribution that presents strong anisotropies in the central area of the vortex.

Fig. 5
figure 5

Effect of consistency correction in accuracy of the WENO SPH scheme for the 2-D vortex simulation

Full size image

To better quantify the gains in accuracy, a convergence analysis of the \(L_2\) errors with resolution refinement is shown in Fig. 6. Curves for the following schemes are displayed:

  • Piece-wise constant ALE-SPH. The very large numerical diffusion introduced by the Riemann solver in the absence of any spatial reconstruction leads to inaccurate results, both in pressure and velocity fields.

  • WENO SPH. The high-order spatial reconstruction largely diminishes the numerical diffusion of the Riemann solver enhancing the accuracy of the scheme. However, as the Lagrangian trajectories followed by the particles are more precise, non-uniform particle arrangements form easily which increase the truncation error of the (standard) SPH interpolators. The latter is responsible for the saturation effect seen in the reduction of the \(L_2\) error with increasing resolutions.

  • Eulerian WENO SPH. The regularity of a distribution formed by particles fixed in the vertices of a Cartesian grid eliminates the truncation error in the SPH interpolations. As a consequence, the saturation effect in the associated convergence curve vanishes. The order of convergence exhibited for pressure \(L_2\) errors reaches a maximum value of 2.5 at the highest resolutions. This is well above second order, moving towards the theoretical third-order convergence of the WENO reconstruction used.

  • Corrected WENO SPH. The use of corrected SPH interpolators for divergencies brings in a very large reduction in the saturation previously seen for the WENO SPH curve. In fact, this Lagrangian scheme exhibits convergence above second order for pressures.

Fig. 6
figure 6

Convergence of \(L_2\) error with resolution refinement in the corrected WENO SPH scheme for the 2-D vortex simulation: a comparison study

Full size image
Fig. 7
figure 7

Effect of the regularised transport in the field of gradient of partition of unity (magnitude) for the WENO SPH scheme in the 2-D vortex simulation

Full size image
Fig. 8
figure 8

Approximation error of the regularised WENO SPH scheme for the 2-D vortex simulation

Full size image
Fig. 9
figure 9

Convergence of \(L_2\) error with resolution refinement in the regularised WENO SPH scheme for the 2-D vortex simulation: a comparison study

Full size image

3.2.2 Results using a regularised transport field via explicit shifting

The alternative to achieve better consistency of the SPH approximation of differential operators exploits the ALE framework of the present scheme by adopting Eq. (38) for the transport velocity. Figure 7 compares two snapshots of the particle distribution obtained with the baseline WENO SPH and the regularised WENO SPH schemes (see Sect. 2.4. As for the results of the previous section, the measurement is taken at time \(\hat{t} = 2\), and for a resolution \(r_0 / \Delta x = 28\). To ease the comparison, the magnitude of the gradient of partition of unity is displayed superimposed (using the same scale in the two snapshots):

$$\begin{aligned} \left\| \nabla _{\hat{x}} C_i \right\| = r_0 \left\| \nabla C_i \right\| = r_0 \left\| \sum _j V_j \nabla W_{ij} \right\| . \end{aligned}$$
(50)

The effectiveness of the regularised transport field in the prevention of particle clustering can be appreciated by simple visual inspection. This is confirmed by the lower values of the magnitude of the gradient of partition of unity.

Figure 8 displays the approximation error in pressures and horizontal velocities computed now with the regularised WENO SPH scheme, at time \(\hat{t} = 2\) and for resolution \(r_0 /\Delta x = 28\). Note that the scales have the same values as in Fig. 5 to ease the comparison. It is clear that the regularised transport field can also be used to improve the accuracy of the computations of the WENO SPH scheme. The reduction in the error displayed is not as large as for the corrected WENO SPH. However, for different applications where the particle distribution might get very distorted, the use of some kind of particle distribution regularisation might be mandatory to preserve the stability of the scheme. Moreover, in this paper the shifting displacement vector \(\delta \textbf{r}_i\) is computed with the formula in Eq. (40) which, due to its explicit nature, does not guarantee obtaining a negligible \(\nabla C\) for all particles in the domain.

Fig. 10
figure 10

Field of particle unitary volume variation in the regularised WENO SPH scheme for the 2-D vortex simulation

Full size image

A convergence analysis of \(L_2\) errors (in pressures and velocities) with resolution refinement for the regularised WENO SPH scheme is displayed in Fig. 9. To perform a comparison study, the same curves corresponding to the baseline WENO SPH and the Eulerian WENO SPH simulations previously displayed in Fig. 6 have been reproduced in the graph as well. It may be seen how the regularised transport algorithm moderates the saturation of the \(L_2\) error reduction with increasing resolutions. The same reasons explained in the previous paragraph justify that this effect is weaker here than what was shown before for the corrected WENO SPH scheme in Fig. 6.

Finally, Fig. 10 displays the field of particle unitary volume variation defined as:

$$\begin{aligned} \Delta V_\textrm{unit} (\%) = 100 \times \frac{V - \left( \Delta x \right) ^2}{\left( \Delta x \right) ^2}, \end{aligned}$$
(51)

computed using the regularised WENO SPH scheme, at \({\hat{t}} = 2\) and for the resolution \(r_0 /\Delta x = 28\). The shifting algorithm used in this paper, in conjunction with the WENO scheme, is able to keep volume variations at a low level (below a 3%), which favours the accuracy of SPH operators.

Fig. 11
figure 11

Pressure and velocity distributions in the \(\hat{z} = 0\) section for the circular blast wave, computed with three different schemes: standard SPH with artificial viscosity (\(\alpha = 0.1\)), baseline WENO SPH (without corrected operators or particle regularisation), and Eulerian WENO SPH

Full size image
Fig. 12
figure 12

Pressure and velocity distributions in the \(\hat{z} = 0\) section for the circular blast wave, computed by the corrected WENO SPH scheme

Full size image
Fig. 13
figure 13

Effect of the regularised transport in the field of gradient of partition of unity (magnitude) for the WENO SPH scheme in the circular blast wave simulation

Full size image
Fig. 14
figure 14

Pressure and velocity distributions in the \(\hat{z} = 0\) section for the circular blast wave, computed by the regularised WENO SPH scheme

Full size image

3.3 Circular blast wave

This is a two-dimensional problem that models a cylindrical explosion which is simulated to verify the capability of the proposed numerical scheme to reproduce shock waves. Periodic boundary conditions have been adopted in both directions. The initial condition consists of a region at high pressure inside of a circle, surrounded by a region at low pressure. By virtue of Eq. (7), this condition can be expressed in densities as follows:

$$\begin{aligned} \rho \left( \textbf{r}, t=0 \right) = \left\{ \begin{matrix} \rho _I, &{} 0 \le \left\| \textbf{r} \right\| \le r_0, \\ \rho _0, &{} \left\| \textbf{r} \right\| > r_0, \end{matrix} \right. \end{aligned}$$
(52)

where \(r_0\) is the circle radius, \(\rho _0\) is the reference density previously introduced in Eq. (7), and \(\rho _I\) is the density in the high pressure inner region (\(\rho _I > \rho _0\)). Both regions are initially at rest. The circular discontinuity in pressures that initially joins the two regions splits into a shock travelling towards the outer boundary of the fluid domain and a rarefaction travelling towards the centre of the circle. In this paper, the values of the problem parameters are \(\rho _0 = 1.0\) kg/m2, \(\rho _I = 1.1\) kg/m2, \(r_0 = 0.5\) m, and the sound speed is set to \(c_0 = 1.0\) m/s. The simulation domain is a circle of radius \(R = 1.5\) m. All simulations have been carried out with a resolution \(r_0 / \Delta x = 120\).

As previously done for the 2-D vortex, the effect of the Liu & Liu [12] correction applied to the divergence operators (see Sect. 2.2) and the particle regularisation (see Sect. 2.4) is analysed. The results are reported using dimensionless variables defined as follows:

$$\begin{aligned} \begin{aligned} \hat{x}&= \frac{x}{r_0},&\hat{z}&= \frac{z}{r_0},&\hat{v}&= \frac{v}{c_0 \sqrt{\frac{\rho _I - \rho _0}{\rho _0}}}, \\ \hat{t}&= \frac{t \, c_0}{r_0},&\hat{p}&= \frac{p}{c_0^2 \left( \rho _I - \rho _0 \right) }. \end{aligned} \end{aligned}$$
(53)

To assess the accuracy of the computations, the solution of equivalent 1-D equations (with a geometrical source term) obtained under the assumption of cylindrical symmetry is used as reference (see Toro [34, Chapter 17] for details). As a preliminary step, simulations have been carried out using three different schemes:

  1. 1.

    standard SPH with artificial viscosity (\(\alpha = 0.1\)) [40],

  2. 2.

    baseline WENO SPH (without corrected operators for divergence terms or a regularised transport field), and

  3. 3.

    Eulerian WENO SPH.

Figure 11 shows the pressure and velocity distributions in the \(\hat{z} = 0\) section computed with the three schemes. The measurement is taken at time \(\hat{t} = 0.66\), and the reference solution is displayed superimposed. Note that plots for the vertical velocity \(\hat{v}_z\) are omitted due to the similarity with those of the horizontal velocity \(\hat{v}_x\). The results computed with the three schemes converge to the reference solution. However, both pressures and velocities computed with standard SPH are very noisy in the whole area that connects the rarefaction and the shock. On the contrary, WENO SPH provides smooth pressure and velocity fields thanks to the great stabilising effect of the Riemann solver. Moreover, the second-order WENO reconstruction provides the accuracy needed to resolve both the rarefaction and the shock. Nevertheless, a spurious wrinkle in the pressure curve can be observed around the radial location where the discontinuity was initially placed. This wrinkle is no longer present in the computations of Eulerian WENO SPH (where particles are arranged in a regular Cartesian grid), which suggests irregularities in the particle distribution are the cause. Next sections aim at mitigating the inaccuracy present in the pressure field for the baseline WENO SPH scheme.

3.3.1 Results using corrected SPH operators for divergence terms

Figure 12 displays the distributions of pressure and horizontal velocity in the \(\hat{z} = 0\) section obtained from the corrected WENO SPH scheme. It may be appreciated the effectiveness of this consistency correction that completely erases the spurious wrinkle in the pressure distribution previously shown in Fig. 11c for the computation with the WENO SPH scheme. For velocities, there is no substantial change, which is coherent with the use of a weakly compressible model.

3.3.2 Results using a regularised transport field via explicit shifting

Figure 13 compares two snapshots of the particle distribution with the magnitude of the gradient of partition of unity superimposed, obtained with the baseline WENO SPH and the regularised WENO SPH schemes. Different from the 2-D vortex, in this acoustic wave propagation problem it is not easy to extract conclusions by simple visual inspection of the particle distribution. Nevertheless, the field of gradient of partition of unity (displayed with the same scale in the two snapshots) shows how the regularised transport field achieves a more regular particle distribution around the radial location of the spurious pressure wrinkle (\(\hat{r} \simeq 1\)).

Figure 14 displays the pressure and horizontal velocity distribution in the \(\hat{z} = 0\) section computed now with the regularised WENO SPH scheme. Despite that the presence of the spurious pressure wrinkle is still evident, the regularised transport field manages to decrease its sharpness. As mentioned before in the analysis of the 2-D vortex, the partial effectiveness of the regularised transport field in erasing the spurious pressure wrinkle is attributed to the explicit nature of the computation of the shifting displacement vector \(\delta \textbf{r}_i\).

To complete the analysis, Fig. 15 presents the field of particle unitary volume variation, defined by Eq. (51), computed using the regularised WENO SPH scheme. It may be seen that for \(\hat{r} \simeq 1\) the volume variations are up to a 16%, notably larger than the values seen for the 2-D vortex. This result further supports the claim that there is big room for improvement in the shifting algorithm.

Fig. 15
figure 15

Field of particle unitary volume variation in the regularised WENO SPH scheme for the circular blast wave simulation

Full size image
Fig. 16
figure 16

Pressure and velocity distributions at time \(\hat{t} = 3\) for the standing wave computed with the WENO SPH scheme

Full size image

3.4 Standing wave

In this section, the evolution of a standing wave in a 2-D rectangular tank is studied for an inviscid fluid. For this test, the wave should oscillate periodically without any change for an infinite number of periods. In this test case, we can therefore evaluate the numerical dissipation of the WENO SPH scheme herein proposed, in comparison with the standard ALE-SPH formulation and with the MUSCL SPH formulation. Moreover, this basic test involves a free surface and wall boundary conditions and thus represents a first test towards the applicability of the present SPH algorithm to real engineering applications.

Specifically, the simulation domain is delimited by a tank of width L and sufficient height. Initially, the free surface is horizontal with a constant water height H. By setting an appropriate initial distribution of velocities (as specified below), and neglecting viscosity in the fluid, the linear theory predicts that the standing wave evolves indefinitely with a constant amplitude [43]. If this wave amplitude, D, is small compared to the water height H, the time evolution of the velocity field \(\textbf{v}\) is determined by the gradient of a potential function \(\psi \):

$$\begin{aligned} \textbf{v} = \nabla \psi , \end{aligned}$$
(54)

where \(\psi \) is given by:

$$\begin{aligned} \psi \left( x, z, t \right) = \psi _0 \left( x, z \right) \cos \left( \omega t \right) . \end{aligned}$$
(55)

Here, \(\omega \) is the angular frequency of the wave oscillations which for a given wave number k can be calculated with the formula:

$$\begin{aligned} \omega ^2 = gk\tanh \left( kH \right) , \end{aligned}$$
(56)

with g being the gravity acceleration. The wave number k in its turn is linked to the length \(\lambda \) of the wave simulated by the formula:

$$\begin{aligned} k = \frac{2 \pi }{\lambda }. \end{aligned}$$
(57)

On the other hand, the velocity potential at \(t=0\), \(\psi _0\), is obtained from the following analytical expression:

$$\begin{aligned} \psi _0 \left( x, z \right) = - \epsilon \frac{Hg}{2\omega } \frac{\cosh \left( k z \right) }{\cosh \left( kH\right) } \cos \left( kx \right) , \end{aligned}$$
(58)

where the origin of the coordinate system chosen corresponds to the bottom left corner of the tank, and \(\epsilon \) is defined as the ratio:

$$\begin{aligned} \epsilon = \frac{2D}{H}. \end{aligned}$$
(59)

Considering that the time derivative of the velocity potential \(\psi \) given by Eq. (55) is null for the initial instant, it can be shown that a hydrostatic distribution approximates the initial pressure field with sufficient accuracy [43]. In the present analysis, the values set for these parameters are \(L = 1\) m, \(H = 1\) m, \(D = 0.025\) m, \(g = 9.81\) m/s2 and \(\lambda = 1\) m. The reference density of the fluid is taken as \(\rho _0 = 1000\) kg/m2, and the sound speed is calculated from the expression \(c_0 = 10 \sqrt{gH}\). Two different types of boundary conditions have been adopted. For the bottom limit, a wall boundary condition has been imposed by discretising the bottom edge with the Modified Dynamic Boundary Conditions (mDBC) [44], whereas periodic boundary conditions have been adopted for the lateral sides of the tank. The free surface does not require a specific treatment in the present WENO SPH scheme: as explained in Sect. 2.3, the lateral reconstruction stencils without almost complete support are simply ignored in the nonlinear combination of candidate polynomial functions.

The results of the analysis are reported using dimensionless variables (marked by a hat symbol, ⌃) defined as follows:

$$\begin{aligned} \begin{aligned} \hat{x}&= \frac{x}{H},&\hat{z}&= \frac{z}{H},&\hat{v}&= \frac{v \, 2 \, \omega }{\epsilon \, H \, g \, k}, \\ \hat{t}&= \frac{t \, \omega }{2 \pi },&\hat{p}&= \frac{p}{\rho _0 \, g \, H}. \end{aligned} \end{aligned}$$
(60)

All the simulations have been carried out for a resolution of \(D / \Delta x = 1.19\). This resolution has been chosen to demonstrate both the applicability of WENO SPH for practical cases as well as the benefits of the WENO reconstruction over other reconstruct methods available in the literature. Note that a finer resolution should be used for this case if more accurate results are needed [45].

Figure 16 displays the pressure and velocity distributions (both horizontal and vertical) after three oscillation periods of the standing wave, computed with the WENO SPH scheme. Note that in this case the particle distribution remains pretty uniform. The results with Corrected WENO SPH and regularised WENO SPH schemes have been thus omitted as they do not introduces any significant improvement. Once again it can be seen that WENO SPH produces a smooth pressure field, which in this case approximates a hydrostatic distribution. For velocities, although no physical viscosity has been modelled, it is expected that the numerical dissipation introduced by the scheme itself will reduce the amplitude of the wave oscillations to zero. However, for the low-dissipative WENO reconstruction used, it can be observed that after the three-period oscillation the magnitude of the velocities is still well in the order of unity.

In fact, the present WENO algorithm can be compared with other reconstructions available in the literature in terms of the dissipation they introduce in the simulation of the standing wave. For this purpose, For this purpose, the total kinetic energy, \(E_{kin}\), at a given time can be calculated for the different ALE-SPH schemes using the following expression:

$$\begin{aligned} E_{kin} = \sum _i^N \frac{1}{2} \, V_i \, \rho _i \, \textbf{v}_i^2. \end{aligned}$$
(61)

Figure 17 shows the time evolution of the total kinetic energy for the standing wave simulation computed with three different schemes: the piece-wise constant ALE-SPH scheme; an ALE-SPH scheme that uses the well-known MUSCL reconstruction [23] (hereafter referred as MUSCL SPH); and the present WENO SPH scheme. Note that the magnitude displayed in the vertical axis of the graph, \(\hat{E}_{kin}\), is obtained by scaling the total kinetic energy with its initial value. The first thing to notice is the very large dissipation introduced by the piece-wise constant ALE-SPH scheme, which completely dampens the wave in less than a quarter of a period. In contrast, the use of the MUSCL reconstruction significantly reduces the numerical dissipation of the standing wave oscillations, which are still present after ten periods. Finally, WENO SPH is clearly a superior scheme, preserving the amplitude of the wave oscillations to a magnitude close to three times larger than that obtained with MUSCL SPH.

Fig. 17
figure 17

Time evolution of the total kinetic energy for the standing wave simulation: a comparison study

Full size image

4 Conclusions

In this paper, a WENO spatial reconstruction for ALE-SPH scheme is completed with two additional alternative approach to address the saturation in the convergence affecting the classical SPH operators. As demonstrated in the literature [30], the high-order WENO spatial reconstruction improves the accuracy of the Lagrangian trajectories which leads to a significant anisotropy in the particle distribution. In the present work, two additions have been proposed to enhance the consistency of the scheme in Lagrangian simulations:

  1. 1.

    high-order divergence SPH operators [12] that approximate differential operators in the governing equations, and

  2. 2.

    a regularised transport field in the ALE framework [22].

Firstly, the propagation of a small pressure perturbation has been used to demonstrate that an SPH formulation based on Riemann problems in conjunction with a high-order WENO reconstruction not only prevents the formation of spurious oscillations in the results but introduces less numerical dissipation than the artificial viscosity. Then, basing on the analytical solution for the (smooth) fluid variable distributions of a 2-D vortex, it has been proved that the two independent methods proposed to enhance consistency in WENO SPH effectively improve the accuracy and mitigate the saturation in global error convergence. For the scheme with corrected SPH interpolation of divergence terms, accuracies similar to those of Eulerian simulations in regular Cartesian grids are obtained in good part of the resolution range, demonstrating convergence above second order in pressures (for second-order WENO polynomials). The scheme has been tested also against the 2-D circular blast wave simulation with the aim of assessing the capability of preserving stability in the presence of strong discontinuities (preventing the formation of spurious pressure noise). The corrected SPH interpolation of divergence terms has been decisive to accurately resolve all the features of the flow. As in the case of the 2-D vortex, the regularised transport field also provides improvements in accuracy though these are smaller. Finally, a inviscid standing wave has been reproduced adopting the WENO and the MUSCL spatial reconstruction, demonstrating that the proposed approach is able to reduce the numerical dissipation also for free-surface flows. It has been claimed as well that the accuracy of the scheme with regularised transport via explicit shifting used in this work could be further improved. Though the implementation of an alternative implicit formulation [26] has been left out of the scope of this paper, it is expected that such scheme would largely improve the results. Contrary to the use of corrected SPH operators for divergence terms, the appeal of the approach with a regularised transport field is the capability to improve the accuracy of the method while keeping strict conservation of global mass and momentum. Future works include a further investigation of the numerical scheme assessing the best type and size of the adopted kernel function, an analysis on larger test cases to assess the efficiency of the proposed scheme in comparison with other existing approaches and, additionally, an investigation adapting the algorithm for boundary conditions to simulate arbitrarily complex solid walls while preserving the same convergence rate.