SPH-ALE Scheme for Weakly Compressible Viscous Flow with a Posteriori Stabilization

A highly accurate SPH method with a new stabilization paradigm has been introduced by the authors in a recent paper aimed to solve Euler equations for ideal gases. We present here the extension of the method to viscous incompressible flow. Incompressibility is tackled assuming a weakly compressible approach. The method adopts the SPH-ALE framework and improves accuracy by taking high-order variable reconstruction of the Riemann states at the midpoints between interacting particles. The moving least squares technique is used to estimate the derivatives required for the Taylor approximations for convective fluxes, and also provides the derivatives needed to discretize the viscous flux terms. Stability is preserved by implementing the a posteriori Multi-dimensional Optimal Order Detection (MOOD) method procedure thus avoiding the utilization of any slope/flux limiter or artificial viscosity. The capabilities of the method are illustrated by solving one- and two-dimensional Riemann problems and benchmark cases. The proposed methodology shows improvements in accuracy in the Riemann problems and does not require any parameter calibration. In addition, the method is extended to the solution of viscous flow and results are validated with the analytical Taylor–Green, Couette and Poiseuille flows, and lid-driven cavity test cases.


Introduction
Smoothed Particle Hydrodynamics (SPH) is a widely used mesh-free method for Computational Fluid Dynamics. The method was originally developed to simulate problems in astrophysics, but it is currently applied to many engineering fields [1][2][3]. In recent years, SPH has gained many improvements in accuracy and efficiency increasing the attention of the scientific community. However, regardless of all the unquestionable advances, SPH has not attained a complete state of maturity yet, and it faces several challenges that require to be addressed by the scientific community [4].
Similar to the Finite Volume Method [5,6], in SPH there are two main approaches to model liquids. One is based on the incompressibility assumption of the Navier-Stokes equations. This assumption leads to the decoupling of the equations and the continuity equation can be considered as a constraint the velocity field has to satisfy. The methods are based on the solution of a Poisson equation for the pressure field, using the pressurecorrection idea from grid-based methods [7][8][9]. This approach is known as Incompressible SPH (ISPH). The second approach, introduced by Monaghan [10] is based on Weakly Compressible hypotheses (WCSPH). In this approach, the incompressibility is approximated by artificially allowing a slight flow compressibility. One advantage of this approach is that it avoids the need for solving a Poisson equation to compute the pressure field. The computation of pressure only requires the use of an equation of state (EOS). As the density of most liquids is nearly constant, a barotropic approximation is reasonable, and a linear EOS depending only on density is often used [11]. Both approaches has advantages and are comparable to the current standards in CFD. Thus, the proposed methodology is a competitive alternative to current state-of-art SPH methods for incompressible flows. In this work, we extend the method proposed in [22] to weakly compressible viscous flows. The usual approach to extend the SPH-ALE methodology to viscous flows is to consider the SPH-ALE discretization of the Euler equations supplemented with an additional term that accounts for the viscous effects [23][24][25]. Here, we propose a different approach based on the use of Moving Least Squares approximations (MLS). This approach does not increase the computational complexity of the scheme since the viscous terms are computed using the same reconstruction already calculated for the convective terms.

Governing Equations
Navier-Stokes equations are derived by imposing the physical principles of conservation of mass, momentum, and energy to a Newtonian fluid system. Adopting an ALE approach the system of equations can be expressed in a differential conservative form by where w w w f rame stands for the velocity of the Lagrangian frame and U U U is the vector of conservative variables. The operator L w w w f rame is called the transport operator linked to w w w f rame . Its application over U U U is designated by L w w w f rame (U U U) and results in L w w w f rame (U U U) = ∂ t U U U + ∇ · (w w w f rame ⊗ U U U). We denote with F F F the non-viscous flux tensor (convective and pressure terms), F F F v represents the viscous tensor and vector S S S T contains the source terms. For two-dimensional cases, the vectors and tensors previously introduced are given by where the fluxes F F F x and F F F y are the rows of the flux tensor F F F , namely, The fluid velocity vector and its components in x and y directions are denoted by u u u = (u, v) T . Density and pressure are designed by ρ and p. We use E for the specific total energy defined as the sum of the internal energy (e) and the kinetic energy according to E = e + 1 2 (u 2 + v 2 ). Moreover, the total enthalpy is defined as H = E + p/ρ. For the diffusive terms, τ ij denotes the viscous tensor component and q i the thermal conduction flux component. For an incompressible Newtonian fluid, τ ij can be expressed ) being µ the dynamic fluid viscosity. Similarly, the thermal flux could be expressed in terms of temperature gradients and thermal conductivity according to ). Finally, in the source term vector (S S S T ), f x and f y represent external force components per unit mass, whereasq h is a volumetric heat source.
In order to model a weakly-compressible flow, we consider two different equations of state (EOS): Tait and Tammann EOS.
Tait EOS models a barotropic fluid, and the pressure only depends on the density, that is, p = p(ρ). Tammann EOS is more general and relates pressure with both the density and the internal energy, that is, p = p(ρ, e). Tait EOS keeps the energy equation decoupled from the momentum equation and can lead to computational cost savings when energy effects on the flow are negligible. However, when shock waves are present in the flow, the Tammann EOS is a more convenient choice. Table 1 shows the expressions to evaluate pressure and acoustic sound speed for any of the EOS adopted in this work. Caloric equations are also provided although its inclusion is optional for a barotropic fluid. The last row in the table presents the whole set of constants values required to properly set each EOS. These values are fixed case by case in the validation section. Zero subindex in Tait equation means the constant is associated to the reference state of the fluid. Tait EOS [26] Tammann EOS [27] Pressure

SPH-ALE Formulation
The computational domain Ω is discretized by a set of N particles at positions r r r i = (x i , y i ) T . Index i is used to label particles and ranges from 1 to N. Each particle has an effective volume V i . Besides the volume, particles have other properties that we refer with a generic variable φ i . Moreover, each particle i has n i neighboring particles inside its compact support domain D i as schematically represented in Figure 1. The SPH-ALE formulation was introduced by Vila and Ben Moussa [14,15] to increase the accuracy and stability of SPH methods in nonlinear systems of conservation laws. Vila and Ben Moussa applied this formulation to the Euler equations and presented a system of equations in semidiscrete form that has many similarities with the finite volume formalism. Since then, several developments were made by incorporating consolidated techniques used in mesh based methods like approximate and partial Riemann solvers, MUSCL, MOOD, and WENO [19,22,28].
In the SPH-ALE formulation, the interaction of each neighboring particle j with the particle i admits a representation as a flux at the midpoint ij located at r r r ij = 1 2 (r r r i + r r r j ). Fluxes are computed from solutions to one-dimensional moving Riemann problems. Thus, we can associate the particle i as the left state, particle j as the right state and the moving interface with the midpoint ij. Figure 2 shows the definition of one of the moving Riemann problems. Unit vector n n n ij points from particle i to particle j. We use index ij − and ij + to denote reconstruction values at the interface from the left and from the right. The kernel gradient can be expressed in terms of the unit vector n n n ij . Kernel functions which depend only on the distance between particles can be expressed as W ij = W(r r r i − r r r j , h ij ) = W(q ij ), The gradient of the kernel function is given by ∇W ij = | ∂W ∂q ij  Several authors have proposed the SPH-ALE formulation for the Navier-Stokes equations [23][24][25]. In these works, the viscous term was discretized using an approximation of the Laplacian based on a hybrid SPH gradient by means of a first-order finite difference scheme [29]. In this work, we propose a different discretization for the viscous term of the Navier-Stokes equations. Observation of the Navier-Stokes in ALE form Equations (1)- (4) suggests that the viscous terms can be computed in the form of a diffusive flux, following a similar approach as the one used for the convective terms.
Thus, the proposed resulting semi-discretized form of the Navier-Stokes equations is given by where Equation (5) expresses the evolution of the conservative variables and Equation (6) describes the evolution of the effective volume associated to particle i. In the above equations, W ij = W(r r r i − r r r j , h ij ) is a kernel function and h ij = 1 2 (h i + h j ) is the smoothing length. The length h i is linked to the particle volume via equation where β is a constant and D is the dimension of the computational domain. Throughout this work, we have set β = 2. Note that, as the interparticle distance dx can be estimated as dx = V 1 D i , we adopt the practical consideration of linking the smoothing kernel length by a constant factor β to the interparticle distance dx.
Moreover, we choose the cubic spline kernel proposed by Monaghan and Lattanzio [30] as kernel function. The support radius for cubic spline kernel is given by i . This implies that, for an initial uniform distribution of particles, R = 4dx, resulting in 9 neighbors for 1D tests and 49 neighbors for 2D problems.
Tensors G G G ij and F F F v ij in Equation (5) account for the Euler and viscous fluxes in the interface ij, respectively. The terms appearing with minus sign inside the parenthesis (H H H i and F F F v i ) are tensors evaluated at the position of particle i that assure at least zero order consistency at discrete level [19].

Discretization of the Convective Terms
The numerical flux tensor G G G ij is computed using the Rusanov flux in the co-moving frame according to The term S + ij is the maximum eigenvalue of the Jacobian matrix which in the Arbitrarian Lagrangian-Eulerian (ALE) framework reads ij is the jump of the reconstructed conservative variables. Moreover, the term w w w ij is the velocity of the reference frame at the interface ij. On an Eulerian frame, w w w ij = 0, whereas on a Lagrangian frame w w w ij = u u u ij = (u u u + ij + u u u − ij )/2. We note that, despite the known diffusive behavior [31] of the Rusanov flux, it can be easily used with different EOS, so it is a convenient choice for the problems addressed here. Tensor In the SPH-ALE formulation, each particle i is associated with a velocity frame w w w i and a material velocity u u u i . The velocity frame w w w i can be freely chosen and determines the evolution of particle positions. For the Eulerian approach of the method we set w w w i = 0 and particles are fixed in space. For the Lagrangian version, we set w w w i = u u u i and perform a weighted average interpolation of the velocity [32] to update particle positions. Therefore, the evolution of the particle position must satisfy for Eulerian/Lagrangian frame Equation (7).

High-Order Reconstruction and Moving Least Squares Approximations
One way to increase the accuracy of the resulting scheme is to compute the reconstruction of the variables at each integration point ij using a high-order approximation. For a given variable φ, which is known on each particle, we can compute the reconstructed variable at integration point, φ ij , by means of Taylor series as where the first and successive derivatives are computed, following our previous work [22], using MLS approximations.

Viscous Terms Discretization
In this work, we propose to extend the same discretization that is typically performed in the Finite Volume method to SPH-ALE discretizations.
is the viscous tensor computed as a function of the state of the i-th particle and the gradients obtained at integration point ij as where the derivatives required to evaluate the viscous tensor are computed with MLS reconstruction. Note that since the MLS reconstruction is performed for the convective terms, it does not require any additional reconstruction procedure to obtain a highly accurate discretization of the fluxes.

The Multi-Dimensional Optimal Order Detection (MOOD) Method
The a posteriori MOOD paradigm, introduced by Clain et al. [21], is used in this work to determine the optimal order of the polynomial reconstruction iteratively by building a candidate solution U * for time t n+1 based on the t n solution. The candidate solution is then run by a series of detectors that check if the solution has a certain set of desirable properties. If any of the particles is flagged as invalid, the candidate solution at that particle is discarded and recomputed from the original solution at t n but using a more dissipative scheme by lowering the polynomial reconstruction degree.

Mood Loop
The MOOD approach is composed of a Particle Polynomial Degree (PPD) and a chain of detectors. The PPD refers to the actual polynomial degree used to compute the candidate solution U * . We evaluate the flux at the midpoint ij between particles i and j, taking the minimum of the respective PPDs for the polynomial reconstruction as PPD ij = min(PPD(i), PPD(j)). The chain of detectors controls the validity of the resulting solution and the particle's PPD is decremented where any of the detectors flag the solution as invalid.
The MOOD loop iterates through the PPD map, initialized with maximal order d max = 3, decreasing the order of the particles that present a non-physical or invalid solution. In this work only orders 3 and 1 are used, being the order 1 called the parachute scheme that, by definition, fulfills all the detectors requirements.
If the PPD map is modified, the candidate solution is declared not valid and therefore must be recomputed. Note that only the particles where the PPD map is modified are recomputed.

Chain Detectors
In order to obtain a stable solution within the SPH formulation, a chain of detectors is used to assess whether the solution is admissible or not. In this work, we employ two detectors: Physical Admissibility Detector (PAD): it checks that all the particles in the solution have positive density and pressure at all times. It also accounts for NaN (Not a Number) values that arise in the candidate solution.
Numerical Admissible Detector (NAD) [33]: relaxed version of the Discrete Maximum Principle (DMP) [21]. It checks that the solution is monotonic and thus, no new extrema are created. It compares the candidate solution with the solution obtained in the previous Runge-Kutta step.
where V i is the set of closest particles and the tolerance δ is defined following [33] as

Numerical Tests
We present the numerical tests selected to assess the ability of the SPH-ALE-MOOD scheme to produce accurate and robust approximations. All the numerical examples have been computed using a third-order Runge-Kutta scheme for time integration.

1d Riemann Problems
The first test cases are devoted to assess the stability and diffusive properties of the SPH-ALE-MOOD scheme. Here, we consider several one-dimensional tests.
The first test case is the 1D Riemann problem (R1) which is one of the four proposed by Marongiu in [34]. In the context of SPH, the works presented in [11,18] also simulate this 1D Riemann problems with Tait EOS. The fluid is water modeled with Tait EOS (ρ 0 = 1000 kg/m 3 , c 0 = 1466.0 m/s, γ = 7 and p 0 = 0 Pa). The domain is [0, 0.1] m and the initial condition is defined as A discretization of 100 particles is used and the solution is advanced up to t f inal = 10 −5 s. The exact solution consists of a rarefaction wave traveling to the left and a shock wave traveling to the right. As a reference solution, we use the exact solution obtained with the algorithm given in [35] applied to the Tait EOS as indicated in [36]. Figure 3 plots the pressure and velocity profiles obtained with the base scheme (firstorder SPH-ALE scheme) [15] and with the SPH-ALE-MOOD model. The SPH-ALE base scheme smears the solution in the shock and rarefaction wave. As expected, for the same number of particles the SPH-ALE-MOOD provides a better representation of the shock front. The front and tail of the rarefaction wave provided for the SPH-ALE-MOOD are also accurately captured and are free of overshoots near discontinuities. Both Eulerian and Lagrangian versions of the scheme produce very similar results, so we only plot here the results obtained with the Lagrangian description.
We consider a second one-dimensional Riemann problem (R2). In this case, the liquid is assumed to follow the Tamman EOS (γ = 7.15 and p c = 3 × 10 8 Pa). This test was proposed by Ivings and Toro in [36] and has been also presented in [37] with a SPH-ALE code with MUSCL reconstruction using a minmod limiter and a finer particle resolution. The initial conditions for this problem are  In Figure 4 we plot the results for pressure, velocity, density, and internal energy at the final time t f inal = 7 × 10 −5 s obtained with the SPH base scheme and the SPH-ALE-MOOD using a Lagrangian description. The SPH-ALE-MOOD improves the results of the SPH base scheme in all the salient features present in the flow. In the density and internal energy plots, it is observed that the resolution of the contact discontinuity is not as sharp as the one obtained for the shock front. We note that the smearing in the contact discontinuity is inherent to the approximations made in the derivation of the Rusanov flux as reported in previous works [19,22].

2d Blast Explosion
The first two-dimensional test considered here is an extension of the one-dimensional shock tube R1 assuming cylindrical symmetry. The fluid is water modeled with Tait EOS (ρ 0 = 1000 kg/m 3 , c 0 = 1466.0 m/s, γ = 7 and p 0 = 0 Pa). The computational domain is a circle of radius R = 0.1 m centered at the origin and the initial conditions are given by The configuration mimics an explosion with a shockwave traveling outwards and a rarefaction moving towards the origin. The reference solution is obtained by using a one dimensional finite volume code with a very fine mesh as explained in [35].
The evolution of the flow is simulated until t f inal = 10 −5 s with the SPH-ALE-MOOD scheme. We consider three different particle initializations to evaluate the effects of the initial positions of the particles on the quality of the numerical results. A radial distribution disposing particles in rings, a Delaunay distribution, which places particles in barycenters of triangles and finally, the third initial layout of particles is the result of applying a random displacement to the Delaunay distribution. The number of particles of the radial distribution (∼90,000) is slightly higher than the one of the Delaunay and Random distribution (∼75,000). Figure 5 shows the initial distributions considered.  Figure 6 shows the density at final time for the three initial particle layouts. The reference solution is represented with a black solid line. It is observed that the results preserve the radial symmetry of the physical problem.

Taylor-Green Flow
The Taylor-Green flow is a classical test for numerical methods for the simulation of viscous flows. It provides an exact solution of the incompressible Navier-Stokes equations in a periodic domain [38,39]. The flow involves the decay of four counter-rotating vortices within the periodic region of size L × L as shown in Figure 8  The exact solution is given in [40] and reads where Re is the Reynolds number of the flow, defined as Re = ρ 0 U 0 L µ 0 . U 0 is a reference velocity magnitude, and ρ 0 and µ are constant values for the density and viscosity of the fluid, respectively.
A global decay kinetic energy factor, denoted by r(t), is defined as the ratio of the overall kinetic energy at time t (E k ) and the corresponding one to initial time E k0 Evaluation of (15) with the velocity field given by Equations (12) and (13) results in an exponential decay according to and by integration in the domain results E k0 = 1 4 U 2 0 L 2 . For the simulations presented in this test, we consider a Taylor-Green flow with L = 1 m, U 0 = 1 m/s and ρ 0 = 1 kg/m 3 . According to the weakly compressible approach, we assume that the fluid obeys the Tait equation with parameters ρ 0 = 1 kg/m 3 , c 0 = 10 m/s, γ = 7 and p 0 = 0 Pa. The case is simulated for Re = 10, Re = 100 and Re = 1000 with the Eulerian version of the SPH-ALE-MOOD scheme. Fluid particles are disposed inside the square domain on a Cartesian arrangement with dx = dy = 1/100 as shown in Figure 9.
The initial conditions are computed using Equations (12)- (14) with the value of the density obtained from the Tait EOS for the analytical pressure.     Figure 11 shows the time evolution of global decay of the kinetic energy, r(t), and the maximum velocity modulus obtained using the SPH-ALE-MOOD method and the corresponding reference incompressible solution for three different Reynolds numbers: Re = 10, Re = 100, and Re = 1000. The numerical results are in close agreement with the analytical solution, showing that the high-order reconstruction of the proposed scheme allows achieving a low dissipation scheme which is accurate for a wide range of Re numbers. Following the work in [41], we have measured the time evolution of the pressure at the center of the domain for Re = 100 and Re = 1000 cases. The results are compared in Figure 12 with the theoretical solution and the solutions obtained with the δ-SPH and δ + -SPH presented by Sun et al. [41,42]. We note that the proposed scheme shows better agreement with the reference solution for all particle resolutions. Moreover, it is remarkable the reduced amount of pressure oscillations, even for coarse discretizations. In Figure 13, the pressure field along y = 0.5 L is shown at t * = 6 for Re = 1000. The results closely follows the theoretical solution even for the coarser discretization. Concerning the computational cost of the proposed scheme, Figure 14 plots the CPU time consumed for different particle discretizations for a simulation until a final time of t * = 2 for Re = 100. As expected, the Eulerian scheme is faster than the Lagrangian method. Then, a possible way for improving the efficiency of the proposed method is to combine Eulerian and Lagrangian particles. This idea has been explored previously in the context of ISPH [43] and fits very naturally in the proposed formulation.

Couette and Poiseuille Flows
Couette and Poiseuille flows are special configurations of the incompressible Navier-Stokes equations that have analytical solution [44]. In both cases, a Newtonian fluid moves between two infinite parallel plates. The Couette flow is driven by the movement of one of the plates whereas the Poiseuille flow is driven by a pressure gradient. As velocity does not vary along the flow direction, a finite length of the plates is considered with periodic boundary condition in left and right sides. Figure 15 shows the geometry model and boundary conditions considered for the simulations. In both configurations the fluid is initially at rest. In the Couette flow, the time-dependent exact solution for the fluid velocity in the x-direction can be expressed as [44,45] u(y, t) = u p L y + where u p is the horizontal velocity of the upper plate and ν is the kinematic viscosity of the fluid.
Similarly, for the Poiseuille flow the transient exact solution is given by [44,45] u(y, t) = f x 2ν where f x denotes a force source term in the x-momentum equation and, as such, it must be included in the source term S S S T of the system of equations defined in (4). The force source term f x and the steady peak velocity in the midplane of the channel u peak are related by expression u peak = 1 8ν f x L 2 . For both problems, the Reynolds number is defined as Re = u max L ν considering the distance between plates L and the maximum velocity u max as the reference length and velocity scales.
In this work, we conduct Couette and Poiseuille simulations for Re = 10. The same value was adopted in the works of Chiron [23], Ferrand [46], and Fourtakas [47]. The fluid is modeled using the Tait equation with ρ 0 = 1 kg/m 3 , γ = 7, c 0 = 10 m/s, and p 0 = 0 Pa. The kinematic viscosity considered for the fluid is ν = 0.1 m 2 /s. The distance between plates is set to L = 1 m and half of this distance is considered for the periodic length in the flow direction.
In the Couette flow, u p is set to 1 m/s, which leads to Re = 10. For the Poiseuille flow, we impose the force source term as f x = 0.8 m/s 2 to produce the same Re number. Figure 16 shows the arrangement of the particles employed for both the Couette and Poiseuille tests. The number of fluid particles between walls and periodic zones is 40 and 20, respectively, resulting in a squared arrangement with distance between particles dx = L/40.
In addition to the fluid particles, we need to incorporate ghost particles for implementing the periodic and wall boundary conditions. For the wall ghost particles we follow the technique used in [13]. A schematic representation of this technique is shown on the right of Figure 16. Dirichlet boundary conditions for velocity on the wall require that ghost particles update their velocity u u u G = (u G , v G ) T following the vector equation.
where u u u F = (u F , v F ) T is the velocity of the mirroring fluid particle and u u u W = (u W , v W ) T is the velocity vector of the wall. In case of fixed walls u u u W = (0, 0) T and for the top moving wall in Couette flow u u u W = (u p , 0) T . As we already have commented, one of the main advantages of the SPH-ALE method is the ability to use either Eulerian or Lagrangian description, and both configurations are able to obtain accurate solutions for this test case. Figure 17 shows the velocity profiles obtained with the SPH-ALE-MOOD scheme for Poiseuille flow at Re = 10 for the Eulerian and Lagrangian description. The exact solution is computed using Equation (18). For t = 20 s the flow is practically in the steady state condition and the obtained numerical solutions agree almost perfectly with the exact solution.   Figure 18 shows the velocity profiles obtained with the SPH-ALE-MOOD scheme for the Couette flow at Re = 1, 10, 100, 1000. The exact solution is computed using Equation (17). At t = 0, the velocity of the moving plate changes abruptly from rest to an horizontal velocity u p forming a sharp discontinuity in the velocity field. A short time after that event, the obtained numerical results slightly deviates over the exact solution. This effect increases with the Reynolds number. For the last time instant, displayed for each Reynolds in Figure 18, the flow has practically reached the steady state, and the velocity profile is linear. The obtained results in the steady state are in close agreement with the exact solution for all the Reynolds numbers computed in this test case. In Figure 18, the deviation from the reference solution observed in the first time instants of the simulations for Re = 100 and Re = 1000, are due to a lack of particles. For these Reynolds numbers, the spatial discretization is not able to capture the abrupt change in the velocity. To verify this, we plot in Figure 19 the results obtained for Re = 100 at t = 0.2 (left) and Re = 1000 at t = 2 (right) for different particle resolutions. It is seen that as the particle resolution increases the deviation is reduced, as expected.

2d Lid-Driven Cavity Flow
The final test presented to assess the behavior of the proposed method is the 2D flow inside a square lid-driven cavity of length L. A schematic setup of the geometry is shown in Figure 20. The lateral and bottom walls are stationary, while the top wall moves horizontally to the right at speed u w . Different Reynolds numbers, namely, Re = 100, Re = 400, and Re = 1000, are studied and results are compared to Ghia [48]. As in the previous case, the velocity of the frame is set to zero adopting the Eulerian version of the SPH-ALE-MOOD scheme.  Figure 21 shows the layout and the type of particles. A Cartesian layout is adopted using a discretization of 100 particles on each side. Lid-driven cavity does not introduce any new type of boundary conditions, but the wall corners need to be taken into account to properly update ghost particle information, as schematically presented in Figure 21 on the top right corner. In order to set the velocity for the corner particle u u u GC , we use a similar technique to the one proposed by Szewc [49]. Focusing on the top-right wall corner and considering the nearest four particles. We have one fluid particle with velocity u u u F , a ghost particle attached to the moving wall with velocity u u u GM , a ghost particle attached to the fixed wall with velocity u u u GF , and a ghost particle in the corner with velocity u u u GC . u u u F evolves with the governing equations, and that u u u GM and u u u GF are updated according to Equation (19). In order to set the velocity for the corner ghost particle u u u GC we impose the velocity in the vertex of the corner as the average of the four particles. Figure 22 shows the horizontal and vertical velocity profiles along the vertical and horizontal center line for Re = 100, Re = 400, and Re = 1000. Simulations were run for a t f inal = 500 s, clearly a time much longer than the one needed to reach the steadystate condition. Results are in good agreement with the reference solution [48] for all the Reynolds number considered. Moreover, the obtained solutions are compared with the ones obtained by Lee et al. [50]. We note that the two schemes use the same number of particles for Re = 400. For Re = 1000, the ISPH scheme from [50] uses a finer discretization (160 2 particles).
The contours of the velocity magnitude superposed with the streamlines after 500 s are shown in Figure 23. It is seen that the scheme is able to reproduce the primary and secondary vortices of the flow.

Conclusions
In this work, we have presented a new high-accurate SPH-ALE method for weakly compressible flow, which can deal with discontinuities. A new approach to compute the viscous flows terms is presented, where Moving Least Squares approximations are used to increase the accuracy and compute the derivatives needed for viscous fluxes. The performance of the proposed scheme is validated with a series of 1D and 2D benchmark problems, and compared with other WCSPH and ISPH schemes from the literature. The proposed method alleviates some of the known drawbacks of weakly compressible schemes. Thus, pressure oscillations are reduced compared with the δ-SPH scheme, and there is no need for an special initialization of the particles. The proposed scheme obtains accurate solutions for any initial distribution of the particles. Moreover, the proposed formulation simplifies the possible combination of the proposed method with grid-based methods (such as the finite volume method, which is closely related to the proposed scheme) or the combination of Eulerian and Lagrangian particles. Even though this issue is not addressed here, it is a very interesting research that will be pursued in forthcoming works.