Towards a High Order Convergent ALE-SPH Scheme with Efficient WENO Spatial Reconstruction

Abstract

This paper studies the convergence properties of an arbitrary Lagrangian–Eulerian (ALE) Riemann-based SPH algorithm in conjunction with a Weighted Essentially Non-Oscillatory (WENO) high-order spatial reconstruction, in the framework of the DualSPHysics open-source code. A convergence analysis is carried out for Lagrangian and Eulerian simulations and the numerical results demonstrate that, in absence of particle disorder, the overall convergence of the scheme is close to the one guaranteed by the WENO spatial reconstruction. Moreover, an alternative method for the WENO spatial reconstruction is introduced which guarantees a speed-up of 3.5, in comparison with the classical Moving Least-Squares (MLS) approach.

1. Introduction

The Smoothed Particle Hydrodynamics (SPH) numerical method was originally conceived for the simulation of astrophysical phenomena [1,2]. Its meshless character and Lagrangian description make it very attractive for the simulation of complex flows with highly distorted moving interfaces, as it was first demonstrated by Monaghan [3]. Since then, SPH has been used in a large number of applications that range from environmental and coastal engineering to energy production [4]. Nowadays, SPH is widely used for Computational Fluid Dynamics [5,6] and applied to a number of problems in hydraulics, including wastewater works [7], turbine design [8], fish passage flows [9], interaction of free-surface flows with flexible structures [10], sloshing in partially filled tanks [11], analysis of Wave Energy Converters (WEC) [12], study of the impact of sea waves on structures [13], and Large Eddy Simulation (LES) modelling of turbulent flows for moderate [14] and, more recently, high [15] Reynolds numbers. For geometrically complex problems and/or multiphysics applications for which the creation of computational grids is a practical burden, SPH already represents a real alternative to more established mesh-based tools.

Nevertheless, the SPH Research and Engineering International Community (SPHERIC), https://spheric-sph.org, accessed on 5 April 2021, has identified a series of Grand Challenges [16] that need to be overcome to further extend the usage of SPH in industrial applications. In particular, the present work addresses the first of these Grand Challenges which is the low accuracy and convergence rate of the method [6]. Indeed, the theoretical second order spatial accuracy of the continuous SPH integral interpolation [17] cannot be achieved with the original SPH formulation [3] and usually the convergence rate practically attained is around first order. This is mainly due to the Lagrangian nature of the scheme, which inevitably leads to an irregular particle distribution and consequent poor accuracy of the SPH spatial operators [18]. Early attempts to increase the accuracy of SPH were based on re-normalisations and corrections [19,20,21] to enforce zeroth and first order consistency of the SPH operators. In a later work [22], this methodology was extended to higher orders of consistency; interestingly, the method kept the consistency also at boundaries where the SPH computational stencils have incomplete support. To alleviate the computational cost, Sibilla [23] derived an algorithm from the previous work [22] where the unknowns are only the derivatives up to the third order (hence fixing to second order the level of consistency attained on the first derivative). An alternative to achieve higher orders of consistency, known as Moving Least-Squares (MLS) hydrodynamics [24], substitutes the use of integral interpolations by the minimisation of a least-squares functional. Overall, the downside of all these methods is that the symmetry of the interactions between particles is broken and hence the attractive conservation properties of SPH are lost. To avoid this there is a second research path in the literature that aims at preventing the formation of irregularities in the particle distribution as the simulation evolves. This line of research encompasses techniques such as the periodic reinitialisation (or remeshing) of the particle locations [25] where the properties of the shifted particles are corrected (as they have been moved in an unchanged velocity field). Previously, the widely-known XSPH scheme [26] already updated the positions of the particles with slightly perturbed non-Lagrangian velocities, and later it was demonstrated that the gradient of the partition of unity of the kernel provides a shifting field more effective in regularising particle distributions [27]. Lind and Stansby [28] demonstrated that high-order convergence in SPH solutions can be achieved if particles are located on a Cartesian grid and kept fixed during the simulation (Eulerian approach). Despite this algorithm no longer featuring the Lagrangian description that had so far characterised the SPH method, this encouraging result has inspired other works that use particle regularisation strategies in conjunction with arbitrary Lagrangian–Eulerian (ALE) schemes to improve the convergence of quasi-Lagrangian SPH schemes (see [29] for a recent example). Finally, there is a third research path to improve convergence in (Lagrangian) SPH schemes, firstly proposed in [30], consisting in using high-order spatial reconstructions to increase the accuracy of particle-particle interactions modelled as Riemann problems (which have a great stabilization effect and do not require tuning of empirical coefficients). In line with this, Avesani et al. [31] introduced a polynomial Weighted Essentially Non-Oscillatory (WENO) reconstruction, improving the accuracy of this ALE-SPH scheme. In a recent work [32], a new Arbitrary Derivative in Space and Time (ADER) integrator has been proposed to avoid the expensive computation of the reconstruction at intermediate stages of a given time step. It is worth mentioning that other approaches to adopt WENO reconstructions in SPH schemes have been proposed in the literature. Contrary to the polynomial WENO above-mentioned, both Vergnaud et al. [33] and Zhang et al. [34] use 1D stencils along each interacting particle pair. Moreover, alternative methods to stabilize mesh-based schemes have been extended to SPH such as the a posteriori Multi-dimensional Optimal Order Detection (MOOD) technique [35].

In the present work, the same WENO spatial reconstruction adopted in [31] has been used for the flux computation in the ALE-SPH scheme. The complex interaction between the WENO reconstruction error and the SPH interpolation error is not well understood yet. For this reason, one of the main aims is to demonstrate whether the convergence properties of the global scheme are inherited from those of the WENO reconstruction as this would imply that increasing the order of the WENO polynomials is an effective way to obtain a high order version of the SPH method. For a synthetic 2D vortex case, it will be demonstrated that the discretization error in the SPH interpolators for divergencies in the governing equations prevents to exploit the high-order capabilities of the polynomial WENO reconstruction, which become apparent in a Eulerian simulation. On the other hand, an alternative way to compute the WENO reconstruction polynomials is proposed: the MLS fits used in [31] have been replaced with a corrected SPH interpolation [22] with the aim of remarkably improving the efficiency of the numerical scheme. Though boundary conditions are not included in the declared scope of this work, a high-order boundary treatment is also required to achieve global high order convergence of the scheme, as demonstrated for the case of Eulerian SPH [36]. For the implementation, the DualSPHysics open-source project [37] has been chosen as baseline for its efficiency and capability of exploiting the computational power of modern graphics processing units (GPUs). Recently, a new boundary treatment has been developed [38] with higher order capabilities (at least, for pressures). It should be noted as well that there is additional research on other high-order mesh-free methods such as the Local Anisotropic Basis Function Method (LABFM) [39] that interestingly collapses to SPH in the limiting low order case.

The rest of this paper is organised as follows: Section 2 provides first a concise but complete description of the ALE Riemann-based SPH scheme with polynomial WENO reconstruction; then summarizes the general procedure to restore consistency in SPH interpolations, and details the context where this correction is applied in the present work. Section 3 gives first the numerical evidence to support the newly proposed method for the generation of the WENO polynomials; then follows the convergence study of the WENO-SPH scheme for the 2D vortex case. Finally, Section 4 draws the conclusions and proposes topics for further research.

2. Materials and Methods

2.1. Arbitrary Lagrangian–Eulerian Riemann-Based SPH

The starting point to derive the expressions of the SPH scheme used in this work is the conservative form of the Euler equations for the conservation of mass and momentum [40]:

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}& =-\nabla ·\left(\rho \mathit{v}\right),\hfill \\ \hfill \frac{\partial }{\partial t}\left(\rho \mathit{v}\right)& =-\nabla ·\left(\rho \mathit{v}\otimes \mathit{v}+p\stackrel{=}{I}\right).\hfill \end{array}$$

(1)

In standard (Lagrangian) SPH schemes the fluid domain is discretized into small computational nodes called particles, with position ${\mathit{r}}_i$ and associated volume $V_i$, that move with the fluid velocity. On the contrary, in an ALE SPH scheme the velocities of these particles, known as the transport velocities (${\mathit{v}}_{0i}$), can be arbitrarily prescribed. The equations for the evolution of the mass and the momentum of particle i can be obtained from Equation (1) after adjusting the advective terms with the transport velocities:

$$\begin{array}{cc}\hfill \frac{1}{V_i}\frac{d\left(V_i{\rho }_i\right)}{dt}& =-{\left[\nabla ·\left(\rho \left(\mathit{v}-{\mathit{v}}_0\right)\right)\right]}_i,\hfill \\ \hfill \frac{1}{V_i}\frac{d\left(V_i{\rho }_i{\mathit{v}}_i\right)}{dt}& =-{\left[\nabla ·\left(\rho \mathit{v}\otimes \left(\mathit{v}-{\mathit{v}}_0\right)+p\stackrel{=}{I}\right)\right]}_i.\hfill \end{array}$$

(2)

Using the identity:

$$\nabla ·\phi =\nabla ·\phi +\phi ·\nabla 1,$$

(3)

where $\phi$ is a generic vector or tensor field, and applying standard SPH operators for the divergence discretization [41], the following expressions are obtained:

$$\begin{array}{cc}\hfill \frac{d\left(V_i{\rho }_i\right)}{dt}& =-V_i\sum _j{V}_j\left({\rho }_i\left({\mathit{v}}_i-{\mathit{v}}_{0i}\right)+{\rho }_j\left({\mathit{v}}_j-{\mathit{v}}_{0j}\right)\right)·\nabla W_{ij},\hfill \\ \hfill \frac{d\left(V_i{\rho }_i{\mathit{v}}_i\right)}{dt}& =-V_i\sum _j{V}_j\left({\rho }_i{\mathit{v}}_i\otimes \left({\mathit{v}}_i-{\mathit{v}}_{0i}\right)+p_i\stackrel{=}{I}+{\rho }_j{\mathit{v}}_j\otimes \left({\mathit{v}}_j-{\mathit{v}}_{0j}\right)+p_j\stackrel{=}{I}\right)·\nabla W_{ij}.\hfill \end{array}$$

(4)

Here $W_{ij}$ is the so-called smoothing kernel function. Following the approach of mesh-based Godunov methods, Vila [30] introduced the idea of computing the fluxes exchanged within a pair of interacting particles via the solution of a Riemann problem located at the midpoint of the pair and that moves at speed ${\overline{\mathit{v}}}_{0ij}=\left({\mathit{v}}_{0i}+{\mathit{v}}_{0j}\right)/2$. These fluxes appear naturally when interpreting Equations (4) as a space centred approximation of:

$$\begin{array}{c}\hfill \frac{d\left(V_i{\rho }_i\right)}{dt}=-V_i\sum _j{V}_j\left(2{\left[\rho \left(\mathit{v}-{\mathit{v}}_0\right)\right]}_{ij}\right)·\nabla W_{ij},\\ \hfill \frac{d\left(V_i{\rho }_i{\mathit{v}}_i\right)}{dt}=-V_i\sum _j{V}_j\left(2{\left[\rho \mathit{v}\otimes \left(\mathit{v}-{\mathit{v}}_0\right)+p\stackrel{=}{I}\right]}_{ij}\right)·\nabla W_{ij}.\end{array}$$

(5)

Two algorithmic elements, that will be discussed in the next subsections, are used to supply the input required for the equations in (5):

• A Riemann solver, that can be exact or approximate;
• A reconstruction procedure to infer the values of mass and momentum fluxes at the midpoint of the $ij$ pair, for both the left (L) and the right (R) states of the Riemann problem.

Furthermore, an extra equation is needed to link the volume variations with the arbitrarily prescribed transport field. Considering (from Continuum Mechanics) that the divergence of a transport field provides the corresponding local rate of variation of a unitary volume, and applying standard SPH operators for the divergence, the following expression is obtained:

$$\frac{d{V}_i}{dt}=V_i\sum _j{V}_j\left({\mathit{v}}_{0j}-{\mathit{v}}_{0i}\right)·\nabla W_{ij}.$$

(6)

Finally, a weakly compressible model is adopted, coupling the two equations in (5) via the equation of state

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

(7)

where $c_0$ is the speed of sound, and ${\rho }_0$ is a reference density. To obtain a larger computational time step, $c_0$ is chosen as 10 times the maximum expected fluid velocity (a common compromise in weakly compressible SPH schemes [3]). The equations in (5) and (6) are explicitly integrated in time with a symplectic scheme, under a CFL condition (see [37] for details).

2.2. Rusanov Flux

The Rusanov flux is used here to (approximately) solve the Riemann problem that appears in (5), as it is simple and computationally efficient. For notation convenience, the conserved variables are grouped together in a column vector

$$Q=\left(\begin{array}{c}\rho \\ \rho \mathit{v}\end{array}\right),$$

(8)

and the mass and momentum fluxes are arranged into a matrix

$$F\left(\rho ,\mathit{v},p,{\mathit{v}}_0\right)=\left(\begin{array}{c}\rho \left(\mathit{v}-{\mathit{v}}_0\right)\\ \rho \mathit{v}\otimes \left(\mathit{v}-{\mathit{v}}_0\right)+p\stackrel{=}{I}\end{array}\right).$$

(9)

With this, the Rusanov flux has the following expression:

$$F_{ij}^{Rusanov}=\frac{1}{2}\left(F\left({\rho }_{ij}^L,{\mathit{v}}_{ij}^L,p_{ij}^L,{\overline{\mathit{v}}}_{0ij}\right)+F\left({\rho }_{ij}^R,{\mathit{v}}_{ij}^R,p_{ij}^R,{\overline{\mathit{v}}}_{0ij}\right)\right)-\frac{c_{ij}}{2}\left(Q_{ij}^R-Q_{ij}^L\right)\otimes {\mathit{n}}_{ij},$$

(10)

where ${\mathit{n}}_{ij}$ is the unit vector pointing from particle i to particle j, and $c_{ij}$ is the maximum value of speed of sound in the pair. Here, since Equation (7) has been adopted as the equation of state, the speed of sound is uniform for all the particles (and constant in time). It is worth mentioning that the second term on the right-hand side of Equation (10) provides the necessary upwinding to numerically stabilise the scheme avoiding the use of artificial viscosity (with the benefit of not having to tune unphysical coefficients for each application).

2.3. High-Order Polynomial Weighted Essentially Non-Oscillatory Reconstruction

To compute the fluxes of Equation (10), it is first necessary to get the values of the conserved variables, $\rho$ and $\rho \mathit{v}$, for the left and right states of each individual Riemann problem. The $ij$ interface is located at the midpoint of the corresponding particle pair; however, at a given time step, the fluid variables are known only at the specific locations of the particles (i.e., a piece-wise constant distribution of data). Hence, it is needed to locally reconstruct the fluid fields from the particle data. It is well known that the piece-wise constant reconstruction (where the variables of the left and right states of the Riemann problem are taken equal to the variables of the i and j particles, respectively) originates too much numerical diffusion [42,43].

In this work it is proposed to use a polynomial WENO (non-linear) reconstruction. Originally extended from mesh-based methods to SPH in [31], the method is suitable for approximating functions that contain discontinuities or exhibit large gradients, preserving the stability of the numerical scheme. It is assumed that the reconstructed functions are polynomials of a specified degree. For a given (central) particle i, the corresponding local polynomial reconstruction has the following form:

$$Q_i\left(\mathit{r}\right)=Q_i+\sum _{k=1}^N{C}_i^k{B}_i^k\left(\mathit{r}-{\mathit{r}}_i\right),$$

(11)

where N is the size of the polynomial basis (that depends on the polynomial degree and on the space dimension), $B_i^k\left(\mathit{r}-{\mathit{r}}_i\right)$ are the associated basis functions, and $C_i^k$ are the (unknown) polynomial coefficients.

For each i-th particle, several candidate stencils for the reconstruction are initially considered, each of them composed by a subset of the particles in the neighbourhood. Figure 1 shows a sketch of the candidate stencils for a sample randomized 2D particle distribution. The so-called central stencil is composed by all the particles at distance lower or equal to the length of the kernel support. On the other hand, there are several lateral one-sided stencils, all of them composed by particles at distance lower or equal to the double of the length of the kernel support. In 2D, there are 8 lateral one-sided stencils: a neighbouring j-th particle belongs to the s stencil ($1\le s\le 8$) when the polar $\theta$ coordinate of the vector ${\mathit{r}}_{ij}={\mathit{r}}_j-{\mathit{r}}_i$ fulfils the expression:

$$\theta \in \left[\left(s-1\right)\frac{\pi }{4},s\frac{\pi }{4}\right].$$

(12)

It should be noted that the particles belonging to the central stencil are also part of one of the lateral one-sided stencils. In general, SPH particles move in time; hence, at each time step, the particles that form a given candidate stencil change.

For each candidate stencil, a (local) candidate polynomial reconstruction is built. Differently from what was done in [31], a corrected SPH interpolation has been used in the present work to compute the polynomial coefficients, as will be explained in Section 2.4.

In this stencil-adaptive WENO method, the final reconstructed polynomial is a weighted average of all the candidate stencil polynomials,

$$Q_i\left(\mathit{r}\right)=\sum _s{\omega }_s{Q}_i^s\left(\mathit{r}\right).$$

(13)

To determine the weights ${\omega }_s$ that appear in (13), the following unnormalized weights are computed first:

$${\tilde{\omega }}_s=\frac{{\lambda }_s}{{\left(\epsilon +{\sigma }_s\right)}^4},$$

(14)

where ${\lambda }_s=1$ for the one-sided lateral stencils ($1\le s\le 8$), and ${\lambda }_0={10}^5$ for the central stencil; $\epsilon$ is introduced to avoid dividing by zero (valid values are ${10}^{-7}$ and ${10}^{-14}$ for single and double precision, respectively); and ${\sigma }_s$ is an indicator of the smoothness of the candidate polynomial based on the values of its coefficients:

$${\sigma }_s=\sum _{k=1}^N{\left(C_{i,s}^k\right)}^2.$$

(15)

It should be noted that Equation (14) favours the contribution of the central polynomial candidate, and effectively filters out the contributions of polynomial candidates that are oscillatory. Finally, a normalization of the weights is done via the expression:

$${\omega }_{s_t}=\frac{{\tilde{\omega }}_{s_t}}{{\sum }_s{\tilde{\omega }}_s}.$$

(16)

Once a local polynomial reconstruction is defined for every particle, the left state of the Riemann problem for the $ij$ pair is obtained from the reconstruction local to particle i, and the right state is obtained from the reconstruction local to particle j.

2.4. Corrected SPH Estimation of Derivatives

In this work, the consistency correction proposed in [22] is applied to compute the coefficients of the polynomials involved in the WENO reconstruction. For completeness, the following paragraphs explain the process to restore consistency in SPH.

The objective is to use SPH summations to compute at the location of a given particle i the values of a function, $f_i$, and its derivatives of successive order, $f_{i,\alpha }$, $f_{i,\alpha \beta }$, …(the indexes $\alpha$, $\beta$, …traverse the different dimension variables of the space). For convenience, a Taylor expansion of f around particle i is used to express the value of f at the neighbouring particle j,

$$f_j=f_i+f_{i,\alpha }\left({\mathit{r}}_j^{\alpha }-{\mathit{r}}_i^{\alpha }\right)+\frac{\left({\mathit{r}}_j^{\alpha }-{\mathit{r}}_i^{\alpha }\right)\left({\mathit{r}}_j^{\beta }-{\mathit{r}}_i^{\beta }\right)}{2!}{f}_{i,\alpha \beta }+\cdots ,$$

(17)

where ${\mathit{r}}_i^{\alpha }$ represents the $\alpha$ component of the position vector of particle i. It should be noted that the Einstein summation convention is adopted.

Without loss of generality, it may be assumed that the objective is to achieve first order consistency (or second order accuracy) in the SPH estimation. Then, all the terms in Equation (17) containing second or higher order derivatives can be neglected, and the remaining unknowns are $f_i$ and the first order derivatives $f_{i,\alpha }$. Multiplying the resulting expression by the product of the volume of particle j times the value of the smoothing kernel, $V_j{W}_{ij}$, and summing for all the neighbouring particles j, a first equation is obtained:

$$\sum _j{f}_j{V}_j{W}_{ij}=f_i\sum _j{V}_j{W}_{ij}+f_{i,\alpha }\sum _j\left({\mathit{r}}_j^{\alpha }-{\mathit{r}}_i^{\alpha }\right)V_j{W}_{ij}.$$

(18)

Taking similar actions but using instead the first order derivatives of the smoothing kernel, $W_{ij,\beta }$, leads to analogous equations (as many as the number of dimensions of the space):

$$\sum _j{f}_j{V}_j{W}_{ij,\beta }=f_i\sum _j{V}_j{W}_{ij,\beta }+f_{i,\alpha }\sum _j\left({\mathit{r}}_j^{\alpha }-{\mathit{r}}_i^{\alpha }\right)V_j{W}_{ij,\beta }.$$

(19)

Equations (18) and (19) form a linear system that can be solved, for instance, with an LU decomposition [44].

Please note that, if achieving a higher order of consistency is required, then higher order derivatives of the kernel function can be considered in the Taylor expansion of Equation (17), and this process can be readily generalised.

Generation of polynomials for the Weighted Essentially Non-Oscillatory Reconstruction

In the WENO reconstruction it is necessary to estimate the polynomial coefficients $C_i^k$ that appear in Equation (11). This is a computationally expensive process that needs to be done for each of the candidate stencils (e.g., nine times in a 2D space). Avesani et al. [31] use a least-squares fit to carry out this task. In the present work we propose to use a corrected SPH interpolation described in Section 2.4 for the estimation of the polynomial coefficients (equivalent to the estimation of the reconstructed function derivatives by virtue of a Taylor expansion). The benefits in terms of computational efficiency and accuracy will be supported with data subsequently.

3. Results and Discussion

3.1. Moving Least-Squares vs. Corrected SPH Interpolation for the Generation of Reconstruction Polynomials

This section focuses on the comparison of the computational performance and the accuracy of an MLS fit used in [31] and the corrected SPH interpolation described in Section 2.4 to compute the coefficients $C_i^k$ of the polynomial reconstruction (11). With this aim, both methods will be used to reconstruct the following 2D scalar function:

$$f\left(x,y\right)=sin\left(\pi x\right)·sin\left(\pi y\right),$$

(20)

in the squared domain $\left[-1,1\right]\times \left[-1,1\right]$ with periodic boundary conditions. The (exact) value of f is known at the particle positions, and two different particle arrangements have been adopted: (a) particles distributed in the vertices of a Cartesian grid with horizontal and vertical spacing equal to a certain length, $\Delta x$, and (b) particles distributed randomly by applying a uniformly distributed random perturbation (of maximum length $\Delta ={\epsilon }_{rand}·\Delta x$, with ${\epsilon }_{rand}=10\%$) to the x and y coordinates of an initial Cartesian grid configuration. For each particle i, a second order polynomial reconstruction:

$$\begin{array}{cc}\hfill P_i\left(x,y\right)=f\left(x_i,y_i\right)+C_i^1\left(x-x_i\right)+& C_i^2\left(y-y_i\right)+\hfill \\ & C_i^3{\left(x-x_i\right)}^2+C_i^4{\left(y-y_i\right)}^2+C_i^5\left(x-x_i\right)\left(y-y_i\right)\hfill \end{array}$$

(21)

is produced from the data of the particles in a neighbourhood of radius 4 times $\Delta x$. To evaluate the reconstructions, a Cartesian reconstruction grid is created initially superimposed to the above-mentioned particle arrangement (a), then the grid is shifted by the distance $0.1·\Delta x$ both in positive x and y directions. For the computation of the error, each reconstructing polynomial will be evaluated at the reconstruction grid point placed closest to the centre of the reconstruction, and compared to the exact value of f. Finally, for the corrected SPH interpolations, the popular cubic spline kernel is used [37], with a smoothing length $h=2\Delta x$.

Figure 2 shows the convergence of the L2 norm of the error in the reconstruction of Equation (20) for the MLS and for the corrected SPH interpolation methods (annotations of the order of convergence are also included). The MLS reconstructions have been computed via a Singular Value Decomposition (SVD) algorithm (see, e.g., chapter 15 in [45] for details). Both above-mentioned uniform and randomized particle distributions have been used. A third order convergence rate is achieved as one would expect from the second order approximation of (21). However, in both particle arrangements, the accuracy of the reconstruction when using corrected SPH interpolations is higher for a given resolution.

More importantly, the average measured computational cost for the method using the corrected SPH interpolation is 3.5 smaller than for MLS. Considering that the WENO method requires one of these polynomial reconstructions per candidate stencil (e.g., 9 polynomials in a 2D case), it is expected that the associated speed-up in the complete SPH scheme will be remarkable.

3.2. Weakly Compressible 2D Vortex

A weakly compressible 2D vortex, as described in [46], has been chosen for the evaluation of the SPH scheme used in this work. This is an inviscid, stationary problem that has circular streamlines. Interestingly, an analytical solution given by the following expressions is available:

$$\begin{array}{cc}\hfill \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]\hfill \\ \hfill u_{\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}}}\hfill \end{array},$$

(22)

where $c_0$ and ${\rho }_0$ have already been defined in (7), ${\rho }_M$ is the asymptotic value of density, and $r_0$ marks the position where the (tangential) velocity is maximal. Please note that the analytical solution depends only on the radial distance r. The actual values of these parameters used in this work are: $c_0$ = 1.0 m/s, ${\rho }_0$ = 1.0 kg/m², ${\rho }_M$ = 1.01 kg/m², and $r_0$ = 0.1 m. The size of the circular simulation domain is R = 0.7 m. There are several reasons to select this test case for the preliminary studies of the present WENO-SPH scheme. Firstly, boundary conditions can be easily enforced with the appropriate number of rows of dummy particles around the simulation domain (where the analytical solution is imposed). This allows to avoid the situation where some of the candidate reconstruction stencils lack of sufficient neighbouring particles to build the corresponding polynomials. On the other hand, the test case is challenging as in general vortices are difficult to capture for Lagrangian particle methods.

For convenience, all the simulation results are expressed in non-dimensional variables, marked by a hat symbol (⌃). These variables are defined as follows:

$$\begin{array}{ccccc}\hfill \widehat{x}& =\frac{x}{r_0},\hfill & \hfill \widehat{z}& =\frac{z}{r_0},\hfill & \hfill \widehat{t}=\frac{t·v_{\theta }\left(r=r_0\right)}{r_0},\\ \hfill \widehat{p}& =\frac{p}{{\rho }_0·v_{\theta }^2\left(r=r_0\right)},\hfill & \hfill \widehat{v}& =\frac{v}{v_{\theta }\left(r=r_0\right)}.\hfill \end{array}$$

(23)

Please note that the non-dimensional time above can be interpreted as the angle (in radians) that the radius vector of the fastest particles (i.e., the ones at $r=r_0$) traverses in that time.

Regarding the SPH scheme, second order polynomials are used for the WENO reconstruction. The kernel is a cubic spline (see definition in [37]), with a smoothing length $h=2\Delta x$ for all the simulations. The CFL number has been set equal to 0.2. Unless explicitly stated otherwise, the simulations are performed adopting a Lagrangian transport velocity (${\mathit{v}}_0=\mathit{v}$).

Figure 3 shows the pressure distribution along the line $\widehat{z}=0$ and at time $\widehat{t}=0.37$, computed with a standard SPH scheme with artificial viscosity (see details in [37]), and two Riemann-based SPH (RSPH) schemes: one with a piece-wise constant reconstruction and the other one with the polynomial WENO reconstruction. The resolution of the simulations is $R/\Delta x=39$, and the exact (stationary) solution is provided for comparison. The value chosen for the artificial viscosity coefficient ($\alpha =0.05$) is typical of applications in the hydraulics field. Please note that the selected physical time corresponds to a sufficiently early instant in the simulation which guarantees that the fluid has not experienced large deformations yet. This is needed to obtain meaningful results in Lagrangian simulations without a dedicated treatment for irregular particle distributions. It may be seen that the artificial viscosity used in the standard SPH scheme introduces a spurious spike around $\widehat{x}=0$. On the other hand, the RSPH scheme is very good at stabilizing the results without spurious effects. However, the unavoidable numerical diffusion that it introduces is too large when using the piece-wise constant reconstruction, affecting the accuracy of the solution. Finally, the high-order reconstruction in the RSPH-WENO scheme largely improves its accuracy.

Figure 4 shows the distribution of particles of the 2D vortex simulation (with the RSPH-WENO scheme) at a later physical time, $\widehat{t}=3.7$, with the magnitude of the gradient of partition of unity superimposed:

$$∥{\nabla }_{\widehat{x}}{C}_i∥=r_0∥\nabla C_i∥=r_0∥\sum _j{V}_j\nabla W_{ij}∥,$$

(24)

In this plot, the resolution is $R/\Delta x=92$. The high-order spatial reconstruction makes the particles move accurately along their Lagrangian trajectories which in its turn leads to anisotropic particle distributions. An immediate consequence is the (local) increase of the gradient of partition of unity which translates into low accuracy of the standard SPH approximations of divergence operators used in Equations (5) and (6).

Aiming at further comparing the different schemes herein adopted, Figure 5 shows the convergence of the L2 norm of the pressure error at $\widehat{t}=0.37$ for the following five different RSPH schemes:

• RSPH with a piece-wise constant reconstruction. This is the less accurate scheme due to the very large diffusion introduced by the Riemann problem with a trivial reconstruction.
• RSPH with a MUSCL reconstruction. Uses linear reconstructions based on gradients estimated by standard SPH operators, with slope limiters (to prevent oscillatory fields) as first depicted by [30]. The MUSCL reconstruction improves the accuracy of the scheme to the point that anisotropic distributions (like the ones displayed in Figure 4 for the RSPH-WENO scheme) appear, especially for high resolutions. The result is an increase in the discretization error that saturates the decrease in the global error.
• RSPH with WENO reconstruction: The high-order spatial reconstruction further improves the accuracy of the trajectories described by the particles which, in conjunction with the usage of standard divergence SPH operators in Equations (5) and (6), ruins the beneficial effect of the high-order reconstruction.
• RSPH with WENO reconstruction and Eulerian transport. In these simulations, particles remain fixed on the vertices of a Cartesian grid (${\mathit{v}}_0=0$). Hence, the gradient of partition of unity is exactly zero in the whole field, unleashing the full potential of the WENO reconstruction to reach a maximum convergence rate of ∼2.7. Considering that the theoretical convergence for the second order reconstructing polynomials used is third order, this result suggests that the overall convergence of the scheme is guided by the order of the WENO polynomials.

Figure 6 displays the distribution of particles of the 2D vortex simulation at $\widehat{t}=3.7$ with the field of unitary volume variation superimposed

$$\widehat{V}-{\widehat{V}}_0=\frac{V-\left(\Delta x·\Delta x\right)}{\Delta x·\Delta x},$$

(25)

computed by the RSPH-WENO scheme. In this plot, the resolution is $R/\Delta x=92$. It may be seen that the algorithm keeps volume variations small, below a 3% (in absolute value).

Finally,

Table 1. Computational time of Riemann-based schemes relative to a standard SPH scheme (for a given resolution).

Riemann-Based Scheme	CPU Time / CPU Time Standard SPH
MUSCL	2.5
WENO	10.9

details the approximate computational time for two different Riemann-based SPH schemes relative to that of a standard SPH scheme with artificial viscosity (see details in [37]), for a given resolution. The large number of candidate polynomials reconstructions in WENO makes this method notably more expensive. However, with the high-order method the same level of accuracy can be reached with a lower resolution, as it may be seen in

Table 2. Performance comparison of standard SPH and Riemann-based with WENO reconstruction SPH, for the 2D vortex case. Values of the L2 norm of the pressure error at $\widehat{t}=0.37$ are provided.

SPH Scheme	$R/\Delta x$	$\parallel err\left(\widehat{p}\right){\parallel }_2$	CPU time / CPU Time Standard SPH
Standard SPH ($\alpha =0.6$)	160	0.020532	1
WENO-RSPH	52	0.020878	0.75

that compares the performance of standard SPH and the Riemann-based SPH scheme with WENO reconstruction. With the much lower resolution required by the WENO SPH scheme, a 25 percent speed-up is actually achieved. Note that the value of the artificial viscosity coefficient $\alpha =0.6$ has been chosen aiming at obtaining a similar amount of numerical diffusion (i.e., a similar depth of the central area in Figure 3) compared to that of the solution computed with the WENO SPH scheme. However, it is not possible to match exactly both the numerical diffusion and the L2 norm of the pressure error so a compromise value has been adopted.

4. Conclusions

In this paper, the convergence properties of a polynomial WENO ALE-SPH scheme have been investigated. For the case of a weakly compressible 2D vortex, it has been demonstrated first that the accuracy of the Lagrangian trajectories of the SPH particles improves with the WENO reconstruction, to the point that anisotropic particle configurations easily form. As a result, the truncation error in the discrete SPH operators becomes dominant and produces the saturation of the global error with increasing resolutions. It has also been shown that, with the aid of a Eulerian simulation with particles fixed on the vertices of a Cartesian grid, this truncation error can be eliminated, obtaining high order convergence (∼3). More importantly, the fact that the order of convergence achieved is very close to the order of the adopted WENO polynomials suggests that this type of schemes is a viable path to reach high order convergence in SPH (which is one of the main goals of the present work).

Furthermore, a novel efficient procedure based on a corrected SPH interpolation for the estimation of the coefficients of the underlying WENO (candidate) polynomial reconstructions has been proposed. Compared to the MLS approach used in [31], this method not only increases the accuracy of the results but brings a speed-up in computational time by a factor of 3.5.

In connection with the convergence analysis of the WENO SPH scheme, and thanks to the ALE framework of the SPH formulation used in this work, it should be noted that it is also possible to tackle the issue of the anisotropic particle distributions by means of particle regularisation techniques (as firstly proposed by [47]) that perturb the Lagrangian transport velocities with a small shifting velocity field that keeps the particle distribution regular. Contrary to the Eulerian simulation performed in the present work, this approach would allow to (approximately) keep the Lagrangian nature that has characterized so far the SPH method (very useful, e.g., in problems with free surfaces). Another interest of this technique in an ALE framework (which is left for future research) is that both conservation and consistency of the scheme would be kept.

The formulation of general boundary conditions (e.g., walls and free surfaces) in the framework of the polynomial WENO-SPH scheme has not been tackled in the present work. The challenge lies in the potential incompleteness of some of the candidate WENO reconstruction stencils when boundaries are present. The appropriate treatment that would keep high order convergence of the numerical solution also close to boundaries remains to be investigated.