Agnès Leroy
email: agnes.leroy@edf.fr

Damien Violeau
email: damien.violeau@edf.fr

Martin Ferrand
email: martin.ferrand@edf.fr

Christophe Kassiotis
email: christophe.kassiotis@enpc.fr

Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH

Keywords: SPH, projection method, incompressible, boundary conditions

come

using a stabilising procedure based on particle shifting. However, we consider an extension of this model to Reynolds-Averaged Navier-Stokes equations based on the k -ǫ turbulent closure model, as done in [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. The discrete SPH operators are modified by the new description of the wall boundary conditions. In particular, a boundary term appears in the Laplacian operator, which makes it possible to accurately impose a von Neumann pressure wall boundary condition that corresponds to impermeability. The shifting and free-surface detection algorithms have also been adapted to the new boundary conditions. Moreover, a new way to compute the wall renormalisation factor in the frame of the unified semi-analytical boundary conditions is proposed in order to decrease the computational time. We present several verifications to the present approach, including a lid-driven cavity, a water column collapsing on a wedge and a periodic schematic fish-pass. Our results are compared to Finite Volumes methods, using Volume of Fluids in the case of free-surface flows. We briefly investigate the convergence of the method and prove its ability to model complex free-surface and turbulent flows. The results are generally improved when compared to a weakly compressible SPH model with the same boundary conditions, especially

Introduction

Modelling incompressible flows with the Smoothed Particle Hydrodynamics (SPH) method has classically been done through weakly compressible SPH (WCSPH) models, as is thoroughly described in [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF]. In this case, the pressure is calculated through an artificial equation of state, which causes the pressure prediction to be noisy and, in many cases, inaccurate. To remedy this issue, truly incompressible approaches were developed in the framework of SPH. In particular, Cummins and Rudman [START_REF] Cummins | An SPH projection method[END_REF] adapted the projection method of Chorin [START_REF] Chorin | Numerical study of thermal convection in a fluid heated from below[END_REF][START_REF] Chorin | On convergence of discrete approximations to the Navier-Stokes equations[END_REF] to SPH by solving a discrete Poisson equation for pressure, leading to an incompressible SPH method (ISPH). Comparisons between ISPH and WCSPH were done by Lee et al. [START_REF] Lee | Truly incompressible approach for computing incompressible flows in SPH and comparisons with the traditional weakly compressible approach[END_REF], which showed that ISPH makes it possible to reduce the computational time while providing a better description of the pressure field than WCSPH. Several versions of the SPH projection method were proposed, the main three of them being: i ) the one proposed by Cummins and Rudman, which consists in maintaining zero divergence velocity, ii ) the one proposed by Shao and Lo [START_REF] Shao | Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free-surface[END_REF], which consists in keeping density invariance and iii ) the one proposed by Hu and Adams [START_REF] Hu | An incompressible multi-phase SPH method[END_REF], based on combining the two previously mentioned methods and thus solving two Poisson equations. In 2009, Xu et al. [START_REF] Xu | Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach[END_REF] made a comparative study between those methods and showed that each of them presented drawbacks. According to the latter authors, imposing the density invariance leads to noisy pressure fields, while imposing the zero velocity divergence gives very smooth pressure fields but leads to instabilities due to particle clustering. On the other hand, the method proposed by Hu and Adams [START_REF] Hu | An incompressible multi-phase SPH method[END_REF], though being stable and providing smooth pressure fields, suffers from very high computational times. Thus, Xu et al. [START_REF] Xu | Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach[END_REF] proposed a stabilising method for the ISPH model based on keeping divergence-free velocity field, which makes it possible to accurately estimate the pressure while keeping computational time smaller than WCSPH. This method consists in slightly shifting the position of the particles at each iteration so as to avoid highly anisotropic particle spacing. The hydrodynamic variables are then corrected by adding the advection term corresponding to the position shift. This method was improved by Lind et al. [START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF], who proposed an expression for the position shift based on Fick's law of diffusion. They also extended the shifting method to free-surface flows. The algorithm proposed by Lind et al. [START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF] seems stable and able to accurately model a great variety of complex flows. Yet, the problem of the pressure wall boundary conditions remains.

A classical way of imposing wall boundary conditions in SPH is the imposition of repulsive forces such as the Lennard-Jones potential [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] or Monaghan and Kajtar's method [START_REF] Monaghan | SPH particle boundary forces for arbitrary boundaries[END_REF]. These methods are easy to implement even for complex geometries and are computationally cheap, but lead to spurious behaviour, as pointed out by Ferrand et al. [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. In particular, the fluid does not remain still near the walls in a hydrostatic case. Besides, they make it difficult -if not impossible -to accurately prescribe a Neumann pressure wall boundary condition, which is a serious issue for ISPH. Most available ISPH models in the literature are thus based on ghost particles [START_REF] Shao | Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free-surface[END_REF][START_REF] Issa | Numerical assessment of the Smoothed Particle Hydrodynamics gridless method for incompressible flows and its extension to turbulent flows[END_REF][START_REF] Lee | Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method[END_REF] or mirror particles [START_REF] Oger | Specific pre/post treatments for 3-D SPH applications through massive HPC simulations[END_REF], for example in [START_REF] Lee | Truly incompressible approach for computing incompressible flows in SPH and comparisons with the traditional weakly compressible approach[END_REF][START_REF] Xu | Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach[END_REF][START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF]. These two techniques are widely used to impose wall boundary conditions in SPH. However, they present serious drawbacks. First, ghost particles are not easy to place in case of complex geometries, especially in three dimensions. Moreover, for nearly all the existing ISPH models combined to ghost or mirror particle methods, a homogeneous Neumann wall boundary condition is imposed on the pressure [START_REF] Lee | Truly incompressible approach for computing incompressible flows in SPH and comparisons with the traditional weakly compressible approach[END_REF][START_REF] Xu | Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach[END_REF][START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF]. This is done by manipulating the relevant entries in the linear system so that the value of the pressure is mirrored across the solid boundary. This is not an exact prescription of Neumann pressure wall boundary condition, and is a serious issue since the proper imposition of pressure boundary condition is crucial when solving the pressure Poisson equation. Yildiz et al. [START_REF] Yildiz | SPH with the multiple boundary tangent method[END_REF] proposed a new method for placing the ghost particles which seemed to improve the accuracy of the imposition of wall boundary condition, but still not exact, and their condition remained homogeneous.

However, in many cases the pressure gradient at a solid wall is non-zero, so that a homogeneous boundary condition is not appropriate. In [START_REF] Hosseini | Pressure boundary conditions for computing incompressible flows with SPH[END_REF], Hosseini Recently, other methods were proposed to model solid boundaries that account for the kernel truncation close to the wall, through the use of a wall renormalisation factor in the SPH discrete interpolation. Kulasegaram et al. [START_REF] Kulasegaram | A variational formulation based contact algorithm for rigid boundaries in two-dimensional SPH applications[END_REF] and De Leffe et al. [START_REF] De Leffe | Normal flux method at the boundary for SPH[END_REF] proposed approximate methods to calculate the renormalisation factor while Feldman and Bonet [START_REF] Feldman | Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems[END_REF] proposed an analytical method for simple wall shapes. In these works, the application of the renormalisation factor in the discrete SPH interpolation formula led to the application of a boundary force in the Navier-Stokes equations. In [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF] With this method the estimation of the fields is very accurate, even close to the walls. Associating the wall boundary conditions proposed by Ferrand et al. [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF] to an ISPH model would make it possible to impose exactly arbitrary Neumann (or Dirichlet) boundary conditions on the pressure, and thus to properly model flows involving complex boundary geometries while taking advantage of the ISPH method. From now on these boundary conditions will be referred to as USAW boundary conditions (unified semi-analytical wall boundary conditions).

Recently, Macià et al. [START_REF] Macià | A boundary integral sph formulation -consistency and applications to isph and wcsph[END_REF] applied the USAW boundary conditions to ISPH, but they focused on the prescription of Dirichlet boundary conditions on the pressure field, which is not appropriate in dynamic cases. Moreover, they did We will first describe the SPH interpolation in the frame of USAW wall boundary conditions. Then, the new ISPH model will be explained. We will see how the algorithm proposed by Lind et al. [START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF] can be adapted to the new boundary description, and see how to impose a non-homogeneous Neumann pressure boundary condition. We include laminar and turbulent (Reynoldsaveraged) flows in the same framework, our purpose being to unify all wall boundary treatment from [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF], including the Poisson equation and the k -ǫ model. In this way our method can deal with basic industrial turbulent flows with a reasonable quality of predictions, though the RANS (Reynolds-Averaged Navier-Stokes) approach is rather crude compared to LES (Large Eddy Simulation) models. Note that the k -ǫ model was applied for the first time to SPH by Violeau [START_REF] Violeau | One and two-equations turbulent closures for Smoothed Particle Hydrodynamics[END_REF] and by Shao [START_REF] Shao | Incompressible SPH simulation of wave breaking and overtopping with turbulence modelling[END_REF], but in these works the boundary conditions were not imposed properly on the turbulent fields.

We will also explain how the method proposed by Bonet and Feldman to compute analytically the wall renormalisation factor can be applied to our description of the solid boundaries, in order to reduce computational time. Finally, the results obtained on several 2-D validation cases will be described. We will investigate the convergence of the method as well as its ability to model complex free-surface and turbulent flows. Comparisons will be done with other numerical methods.

SPH interpolation in the frame of unified semi-analytical wall boundary conditions

SPH is now a well-known method, and we assume the reader is rather familiar with its basics. Thus, we will not describe the classical SPH interpolation and operators. For an extensive description and analysis of the SPH method see [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF][START_REF] Violeau | Fluid Mechanics and the SPH method[END_REF]. In this section, however, we will summarise the USAW boundary conditions used herein. In this work, fluid particles which do not belong to a boundary are called "free" particles. Solid boundaries in the method proposed by Ferrand et al. [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF] are modelled by vertex particles v ∈ V and segments s ∈ S (see Figure 1). The vertex particles are truncated fluid particles placed at the wall, with velocity imposed as equal to the wall velocity. They were introduced to compute more accurately the fields and their derivatives close to the walls.

They are specially important when dealing with turbulence, where the fields Throughout this work, we will use the 5th order Wendland kernel [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree[END_REF], defined by:

w h (q) = α W,2 h 2 1 - q 2 4
(1 + 2q) for q q max = 2

where α W,2 is its normalising constant, equal to 7/4π in 2-D, h is the smoothing length and q = |rr ′ | h with r and r ′ two position vectors. We impose h = 2δr for all the simulations.

With the present boundary conditions, the SPH discrete interpolation of a field A at particle a with position r a reads:

[A] γ a = 1 γ a b∈P V b A b w ab (2)
where V b is the volume of particle b and w ab = w h (r ar b ). γ a is the wall renormalisation factor mentioned in the introduction, defined as in [START_REF] Kulasegaram | A variational formulation based contact algorithm for rigid boundaries in two-dimensional SPH applications[END_REF] and [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]:

γ a = Ω∩Ωa w(r a -r ′ )d n r ′ ( 3
)
where Ω is the fluid domain, Ω a is the compact support of the kernel at particle a and n is the space dimension. Note that γ a is equal to 1 inside the fluid, due to the normalisation property of the kernel. On the other hand, γ a is inferior to 1 when the kernel support intersects the wall. In the method proposed by Ferrand et al. [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF], γ a is computed as:

dγ a dt = u R a • ∇γ a (4)
where u R a represents the particle's velocity with respect to the wall. However, we will present in section 3.4 another method to compute γ a , as accurate as the latter but which decreases computational time.

In this framework, the discrete SPH differential operators are different from the classical ones [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. The antisymmetric form of the discrete gradient reads:

(∇A) a ≈ G γ,+ a {A b } = 1 γ a b∈P V b (A a + A b ) ∇w ab - 1 γ a s∈S (A a + A s ) ∇γ as (5)
where m and ρ are the mass and density of particles. The latter being kept constant in ISPH, we will omit the particle subscript in its notation. ∇γ as is the contribution of segment s to the gradient of γ a , defined as:

∇γ as = ∂Ωs∪Ωa w(r a -r ′ )n s d n-1 r ′ (6)
∂Ω s is the boundary area spanned by segment s and n s is the inward normal to the wall on s (see Figure 1). The following property holds [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]:

∇γ a = s∈S ∇γ as (7)
It is also possible to define a discrete symmetric gradient:

(∇A) a ≈ G γ,- a {A b } = - 1 γ a b∈P V b A ab ∇w ab + 1 γ a s∈S A as ∇γ as (8)
where A ab = A a -A b and A as = A a -A s . In case the discrete gradient of a vector field is calculated, the formulae ( 5) and ( 8) remain unchanged except that A ab ∇w ab and A as ∇γ as are replaced by A ab ⊗∇w ab and A ab ⊗∇γ as respectively.

In the SPH literature, it is recommended to use the symmetric gradient for accurate estimation of the required quantities (e.g. velocity gradients) while the antisymmetric form is more relevant when estimating the pressure gradient in the momentum equation, for reasons of momentum conservation (see e.g. [START_REF] Price | Smoothed particle hydrodynamics and magnetohydrodynamics[END_REF]).

The symmetric form of the divergence operator reads:

(∇ • A) a ≈ D γ a {A b } = - 1 γ a b∈P V b A ab • ∇w ab + 1 γ a s∈S A as • ∇γ as (9)
Finally, the discrete Laplacian operator proposed by Ferrand et al. reads:

[∇ • (B∇A)] a ≈ L γ a {B b , A b } = 1 γ a b∈P V b (B a + B b ) A ab r 2 ab r ab • ∇w ab - 1 γ a s∈S [B s (∇A) s + B a (∇A) a ] • ∇γ as ( 10
)
where B is a diffusion coefficient for the field A, r ab = r ar b and r ab = |r ab |.

In case A is a vector, the Laplacian will be denoted by L γ a {B b , A b }. In case B = 1, it will be denoted by L γ a {A b }.

Formulation of the new ISPH model

Discrete Navier-Stokes equations for an incompressible flow

The Navier-Stokes equations in a continuous framework for an incompressible flow read:

∇ • u = 0 du dt = - 1 ρ ∇p + 1 ρ ∇ • (µ E ∇u) + g (11)
where the density, velocity, pressure, time, dynamic viscosity and gravity are noted respectively ρ, u, p, t, µ and g. We recall that µ = ρν, ν being the kinematic molecular viscosity. We defined µ E = µ+µ T where µ T is the turbulent dynamic viscosity, equal to zero for a laminar flow. We also defined:

p = p + 2 3 ρk (12)
which is used for turbulent flows in the context of Reynolds-averaged fields, k being the turbulent kinetic energy. The discrete SPH Navier-Stokes equations are then based on the discrete operators proposed in the previous section:

D γ a {u b } = 0 du a dt = - 1 ρ G γ,+ a {p b } + 1 ρ L γ a {µ E,b , u b } + g (13)
ν T,a = µ T,a ρ is the turbulent kinematic viscosity of particle a. It is a function of its turbulent kinetic energy k a and dissipation rate ǫ a , according to:

ν T,a = C µ k 2 a ǫ a (14)
where C µ is a constant defined in Table 1.

The ISPH model deals with the resolution of (13) through a procedure described in the next section. In the laminar case, k, ǫ and ν T are equal to 0 and it is only necessary to impose a Neumann boundary condition on the velocity in the viscous term of [START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and multigrid method[END_REF], which is done by writing [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]:

L γ a {µ, u b } = 2µ γ a b∈P V b u ab r 2 ab r ab • ∇w ab - 2µ γ a s∈S u R as • t as δr as t as |∇γ as | (15)
where u R as is the particle's velocity with respect to the segment s, and:

t as = u R as -u R as • n s |u R as -(u R as • n s ) | (16)
We also defined:

δr as = max(r as • n s , δr) (17)
where r as = r ar s .

The use of the k -ǫ model is the same as in WCSPH with the present boundary conditions (see [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]). More information about the k -ǫ model can be found in [START_REF] Launder | Mathematical models of turbulence[END_REF][START_REF] Pope | Turbulent flows[END_REF]. The quantities k a and ǫ a are calculated through a semi-implicit time-scheme involving the SPH form of the standard k -ǫ model (see [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]):

k n+1 a -k n a δt = P n a -ǫ n a k n+1 a k n a + 1 ρ L γ a µ + µ n T,b σ k , k n b ( 18
)
ǫ n+1 a -ǫ n a δt = ǫ n a k n a C ǫ1 P n a -C ǫ2 ǫ n+1 a + 1 ρ L γ a µ + µ n T,b σ ǫ , ǫ n b ( 19
)
where σ k , C ǫ1 , C ǫ2 and σ ǫ are constants described in Table 1, the superscripts n and n + 1 refer to the time iteration numbers and δt is the time step. P a is the production of turbulent kinetic energy of particle a and is calculated according to a semi-linear model [START_REF] Guimet | A linearised turbulent production in the k-ǫ model for engineering applications[END_REF]:

P a = min C µ k a S a , ν T,a S 2 a ( 20
)
where S a = √ 2S a : S a , with:

S a = 1 2 G γ,- a {u b } + G γ,- a {u b } T (21) Table 1: Values of the k -ǫ model constants [21] κ C µ C ǫ1 C ǫ2 σ k σ ǫ 0.41 0.09 1.44 1.92 1.0 1.3
The imposition of boundary conditions on u, k and ǫ is a crucial issue in the k -ǫ model. Our choices in terms of boundary conditions were based on the Code_Saturne Theory Guide [START_REF]Code_Saturne 3.0.0 theory guide[END_REF] (see Part IV-B), which describes the implementation of the k -ǫ model in a well-established Finite Volumes (FV) code. Thus, the following equations ( 22) to ( 29) can be considered as an SPH equivalent of the latter FV approach. Note that the Neumann wall boundary conditions are applied through the Laplacian operator like in mesh-based methods, whereas the Dirichlet boundary conditions are imposed at the vertex particles which are involved in the volumic terms. Thus, in the aforementioned equations the particles a belong to F , whereas the particles b belong to P = F ∪ V .

A non-homogeneous Neumann wall boundary condition on the velocity is applied in the viscous term of ( 13) by writing:

1 ρ L γ a {µ E,b , u b } = 2µ γ a ρ b∈P V b u ab r 2 ab r ab • ∇w ab - 2 γ a s∈S u 2 * ,as t as |∇γ as | (22)
where u * ,as is the friction velocity at the wall seen by particle a, computed through an iterative process solving the following implicit equation:

u R as • t as u * ,as = 1 κ ln δr as u * ,as ν + 5.2 (23)
where κ is the Von Karman constant (see Table 1).

On the other hand, the velocity at the vertex particles is left to evolve as a function of the viscous term:

u n+1 v = u n v + δt 1 ρ L γ v {µ E,b , u b } (24)
but its normal component is imposed to be equal to zero by projecting u n+1 v along the tangent to the wall.

A homogeneous Neumann wall boundary condition on the turbulent kinetic energy is applied in (18) by writing:

L γ a {µ + µ T,b σ k , k b } = 1 γ a b∈P V b 2µ + µ T,a + µ T,b σ k k ab r 2 ab r ab • ∇w ab (25)
A compatible Dirichlet boundary condition on k is imposed at all vertex particles v:

k v = 1 α v b∈F V b k b w vb (26)
As for the dissipation ǫ, it was necessary to improve the treatment of the diffusion boundary term in [START_REF] Kawamura | DNS of turbulence and heat transport in a channel flow with different Reynolds and Prandtl numbers and boundary conditions[END_REF] compared to what was proposed in [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. Indeed, the formulation they proposed did not give correct results close to the walls. Thus, the Neumann wall boundary condition on the dissipation rate is applied in [START_REF] Kawamura | DNS of turbulence and heat transport in a channel flow with different Reynolds and Prandtl numbers and boundary conditions[END_REF] by writing:

L γ a {µ + µ T,b σ ǫ , ǫ b } = 1 γ a b∈P V b 2µ + µ T,a + µ T,b σ ǫ ǫ ab r 2 ab r ab • ∇w ab + 4C µ σ ǫ γ a s∈S k 2 a δr as |∇γ as | (27)
A compatible Dirichlet boundary condition is imposed on ǫ at all vertex particles v:

ǫ v = ǫ s1 + ǫ s2 2 (28)
where:

ǫ s = 1 α s b∈F V b ǫ b + 4C 3/4 µ k 3/2 b κδr bs w sb (29)
A justification for this formulation of the boundary conditions on ǫ is given in the Appendix A. Note that the wall boundary conditions imposed on ǫ have a great impact on the flow representation.

New ISPH algorithm

The present model follows the structure of the one proposed by Lind et al. [START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF], which is based on the projection method proposed by Cummins and Rudman [START_REF] Cummins | An SPH projection method[END_REF], where the velocity field is maintained divergence-free, and a stabilising method consisting in a particle shift is included. First, the particles are moved to an intermediate position r * according to:

r * a = r n a + δt 2 u n a ( 30
)
An estimation of the velocity field is then done based on the momentum equation without the pressure gradient term, so that:

u * a -u n a δt = 1 ρ L γ a {µ E,b , u n b } + g (31)
u n a is the velocity at time n at particle a and u * a is the estimated velocity field. The second part of the momentum equation reads:

ũn+1 a -u * a δt = - 1 ρ G γ,+ a {p n+1 b } (32)
where ũn+1 a is the velocity calculated through the projection method at time n + 1. Applying the divergence operator to the continuous form of this equation while imposing that ∇ • ũn+1 a = 0 gives a pressure Poisson equation that has to be solved to calculate the pressure field at the next time-step. Considering that the density is constant in our work, this equation reads:

∇ 2 pn+1 a = ρ δt ∇ • u * a ( 33
)
Here we do not write it in a discretized form here since its discretization will be dealt with in section 3.3. After the pressure is calculated, the velocity field is corrected by applying equation [START_REF] Meyer | Assessment of implicit large-eddy simulation with a conservative immersed interface method for turbulent cylinder flow[END_REF], and the new position of the particles is calculated according to:

rn+1 a = r * a + δt 2 ũn+1 a (34)
To ensure the stability of the simulations, after the particles were moved according to equation [START_REF] Monaghan | SPH particle boundary forces for arbitrary boundaries[END_REF], they are slightly shifted according to a diffusion law:

r n+1 a = rn+1 a + δr a (35)
where:

δr a = -C shift h 2 ∇C a (36)
C shift is a diffusion coefficient, the value of which having been calibrated from various test cases and taken equal to 0.7 for the Wendland kernel. ∇C a is the gradient of the particle concentration. The following discrete gradient was used instead of the classical one described by Lind et al.:

∇C a ≈ G γ a {1} = 1 γ a b∈P V b ∇w ab - 1 γ a s∈S ∇γ as (37)
In this expression, the boundary term running over the segments s prevents the particles from leaving the domain when the diffusion is applied near a wall.

In their work, Lind et al. [START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF] observed that it was necessary to introduce an additional term in the concentration gradient in order to avoid particle clumping.

This was not the case in the present work due to the fact that we use a Wendland kernel, which is known to avoid particle clumping due to the positiveness of its Fourier transform [START_REF] Dehnen | Improving convergence in smoothed particle hydrodynamics simulations without pairing instability[END_REF]. Moreover, applying the particle shift close to the free-surface would lead to an unphysical behaviour of the particles due to the kernel truncation, which is not accounted for near the free-surface, even with our boundary conditions. To avoid this issue, the shift is not applied to the particles whose distance to the free-surface is lower than hq max /2 (see eqn.

(1) for the definition of q max ). This criterion was established by numerical experiments on the dam-break over a wedge case (section 4.3). It was expressed as a function of hq max in order to have the same number of particle layers not shifted near the free surface, regardless of the kernel choice. This method is equivalent to the one proposed by Lind et al. [START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF] applying a modified particle shift near the freesurface. After the particles' positions are shifted, their velocities are modified according to a first-order Taylor expansion:

u n+1 a = ũn+1 a + G γ,- a {ũ n+1 b } • δr a (38)
In case of turbulent flow, a similar process is applied to the turbulent kinetic energy and dissipation rate. This marks the end of a time-step and a new one begins with [START_REF] Mayrhofer | Study of differential operators in the context of the semi-analytical wall boundary conditions[END_REF].

To perform simulations including free-surfaces with ISPH, it is necessary to impose zero pressure at the free-surface (Dirichlet boundary condition). Thus, the particles belonging to the free-surface have to be tracked, which is done trough a criterion based on the value of the divergence of their position, D γ a {r b }. Indeed, ∇ • r should be equal to n in n dimensions, which is not exactly verified near the free-surface due to the kernel truncation. A particle is identified as belonging to the free-surface if D γ a {r b } ≤ 1.5 in 2-D [START_REF] Lee | Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method[END_REF]. With this tracking, however, some particles happen to cross the wall when they belong to thin jets impacting it with high velocity (typically 3-4 particles in the case of the jet impacting a wall in the triangular wedge case, see Section 4.3). This is because their pressure is set to zero while they reach the wall so that they end by crossing it. To solve this issue, for each free particle with divergence of the position lower than 1.5, a test is performed to check whether it will cross the wall at the next time-step, which is done through the following criterion:

r av + δt u a • r av r 2 av < hq max 8 (39)
If the latter relation holds, the free particle a and the vertex particle v are not identified as free-surface particles. This technique was tested on the triangular wedge case (Section 4.3). Other techniques exist to track the free-surface (see [START_REF] Marrone | Fast free-surface detection and level-set function definition in SPH solvers[END_REF] for example), but the present one proved sufficient to ensure the impermeability of the walls while solving properly the pressure Poisson equation.

In summary, the structure of the algorithm is almost unchanged compared to the one proposed by Lind et al. [START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF], but the differential operators used are different. Our model performs better near the walls without the use of ghost particles and includes turbulence treatment. We saw that the free-surface detection and the shifting algorithms were slightly modified, but the most important change concerns the Laplacian operator. We will now see how this change makes it possible to impose an accurate non-homogeneous Neumann pressure boundary condition.

Laplacian operator and imposition of wall pressure boundary conditions

In the pressure Poisson equation, the Laplacian operator ( 10) is involved with B = 1 and A = p. In [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF], the summation term involving the segments s is the boundary term. The treatment of the latter is crucial, since pressure wall boundary condition are applied through it. It involves the pressure gradient at the segments and at the fluid particles close to the wall. Here we assume that (∇p) a • n s ≈ (∇p) s • n s , which is justified by the fact that the pressure field does not vary much near the walls. Using the fact that ∇γ as is oriented along n s by the definition (6), we obtain:

L γ a {p b } = 2 γ a b∈P V b pab r 2 ab r ab • ∇w ab - 2 γ a s∈S (∇p) s • ∇γ as (40)
The notation L γ a from now on will refer to this expression instead of the one of equation [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF] when it is applied to the pressure. This formulation of the Laplacian led to instabilities on hydrostatic cases since it is not first order consistent.

To solve this issue, we use the equality:

∇ 2 p = ∇ 2 (p + ρgz) ( 41
)
where z is the vertical coordinate and g the gravity magnitude. We now solve a modified Poisson equation:

L γ a {p b + ρgz b } = ρ δt D γ a {u * b } (42)
which is an SPH form of [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF]. D γ a is given by [START_REF] Feldman | Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems[END_REF].

It is now necessary to define the pressure gradient term at the segments,

(
∇p n+1 v • n v = ρ δt (u * v -u wall v ) • n v (43)
The same applies to the segments since their velocity is calculated as the average of the velocities of the vertex particles at each of its vertices [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. Finally, the discrete pressure Poisson equation ( 42) can be written as:

2 γ a b∈P V b pab + ρgz ab r 2 ab r ab • ∇w ab = 2ρ γ a s∈S u * s -u wall s δt + g • ∇γ as + ρ δt D γ a {u * b } (44)
One can now check on a simple case that this pressure wall-boundary condition is physical. Let us consider the case of a fluid at rest with a free-surface in a rectangular tank. Following the steps of the projection method, we have:

u * s = δtg (45)
because the velocity at the initial time n is equal to zero. Then the condition imposed on the pressure gradient at the wall is:

∇(p n+1 s ) • n s = ρg • n s (46)
which is the expected non-homogeneous boundary condition. Thus we see that the condition [START_REF] Trask | Accuracy and performance of implicit projection methods for transient viscous flows using SPH[END_REF] provides the exact pressure condition in order to balance gravity forces on a horizontal bed. This condition is non-homogeneous in many cases since the right-hand side depends on viscous and external forces through u * . In Hosseini and Feng's paper [START_REF] Hosseini | Pressure boundary conditions for computing incompressible flows with SPH[END_REF], the same boundary condition was imposed on the pressure at inflow or outflow boundaries.

Equation ( 44) corresponds to a linear system:

A • p = B ( 47
)
where p is the vector of all particle pressures, B is the vector of right-hand side values at all particles and A is a non-symmetric sparse matrix corresponding to the discrete Laplacian operator. To solve this system, linear solvers like GMRES [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] or Bi-CGSTAB [START_REF] Van Der | Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems[END_REF] are used. In the case of confined flows, if no Dirichlet condition is imposed the system has an infinity of solutions, and the matrix A is not invertible. We make it invertible by adding a small perturbation through a slight reinforcement of the diagonal terms.

Reducing computational time: analytical computation of γ a with the 5th order Wendland kernel

To ensure stability, several conditions concerning the time-step value have to be imposed [START_REF] Morris | Modeling low Reynolds number incompressible flows using SPH[END_REF][START_REF] Violeau | Maximum time step for keeping numerical stability of viscous weakly compressible SPH[END_REF], namely the Courant-Friedrichs-Lewy (CFL) condition and others relative to viscous forces and acceleration, which reads:

δt = min (δt CFL , δt visq , δt force , δt γ ) (48)
where:

                     δt CFL = C CFL h u ref δt visq = C visq min a∈P h 2 ν E,a δt force = C force min a∈P h ua (49)
C CFL = 0.2, C visq = 0.125 and C force = 0.25 are constants which were found by numerical experiments (see Morris et al. [START_REF] Morris | Modeling low Reynolds number incompressible flows using SPH[END_REF] for the last two values). ua is the total acceleration of particle a and u ref is a reference velocity, which corresponds to the numerical speed of sound for a WCSPH model and to the maximum velocity in the fluid for an ISPH model. The speed of sound in WCSPH is usually fixed as 10 times the maximum flow velocity [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF]. In most simulations the CFL prevails, which leads to a time-step 10 times smaller with WCSPH than with ISPH. But when the calculation of γ a is done through (4), an additional condition on δt has to be imposed, which reads [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]:

δt γ = C t,γ 1 max a∈P,s∈S ∇γ as • u R a (50)
We recall that u R a was defined in eqn. (4). It was found by numerical experiments that C t,γ = 0.004 [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. This condition appeared to prevail in many cases, so that the time-step size should be the same for ISPH and WCSPH. This would be a major drawback for ISPH with these boundary conditions since the matrix inversion makes the method much slower than WCSPH for a given value of the time-step. This is why we propose a method to compute γ a analytically without solving (4), which avoids the condition [START_REF] Yildiz | SPH with the multiple boundary tangent method[END_REF]. It follows the idea proposed by Feldman and Bonet [START_REF] Feldman | Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems[END_REF], which consists in writing γ a as a boundary integral by applying Gauss's theorem to (3):

γ a = - ∂Ω W (|r a -r ′ |) • n (r ′ ) d n-1 Γ (r ′ ) ( 51
)
where W is defined as:

w (q) = ∇ • W (q) ( 52
)
Since w is a radial function, so is W:

W (q) = -ϕ (q) r (53)
where q = r h and r = r ar ′ . Then we have :

γ a = ∂Ω ϕ (q) r • n (r) d n-1 Γ (r) (54)
Here we only consider the case of a two-dimensional space. The calculations were done for the 5th order Wendland kernel (1), which yields:

ϕ (q) = 1 2πh 2 q 2 1 - q 2 5 1 + 5q 2 + 2q 2 for q 2 (55)
As pointed out in [START_REF] Feldman | Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems[END_REF], the function h 2 ϕ (q) presents a singularity in q = 0, so that the Gauss theorem invoked to obtain (51) is only valid if the integration Figure 2: Notations for the calculation of γas .

is done on ∂Ω ∪ ∂Ω ǫ , with Ω ǫ a small sphere of radius ǫh centred on q = 0. By decreasing ǫ to zero, it is possible to show that:

γ a = 1 - s γas (56)
with:

γas = -n s • s ϕ r h rd n-1 Γ (r) (57)
Recall n s is the inward unit normal on segment s. Let t s be the unitary vector tangential to s (see Figure 2). We have r = r 0 -r a +yt s so (57) can be simplified to give:

γas = n s • r a0 s ϕ r h d n-1 Γ (r) (58)
where r a0 = r ar 0 and r 0 is the orthogonal projection of r a on the segment s. Let y be the coordinate along t s , r a0 = |r a0 | the distance from the integration point to the segment and q a0 = r a0 /h. We define:

ζ i = min 1 2 q2 a0 + y 2 i h 2 , 1 (59)
for i = 1, 2, with qa0 = min (q a0 , 2). y spans the interval

[y 1 = r av1 • t s ; y 2 = r av2 • t s ].
It is then found that:

γas = sign (n s • r a0 ) 1 4π [sign (y 2 ) ψ (q a0 , ζ 2 ) -sign (y 1 ) ψ (q a0 , ζ 1 )] (60)
with:

ψ (q, ζ) = q ζ 2 -q 2 4        -4 3 ζ 5 + 7ζ 4 -5 12 q 2 + 14 ζ 3 + 7 3 q 2 + 5 ζ 2 -1 4 5 8 q 2 + 21 q 2 ζ + 7 6 q 4 + 35 6 q 2 -7        -5 8 q 2 + 21 q 5 16 arg cosh 2ζ q + 2 arctan 4ζ 2 q 2 -1 (61)
Let us now consider the particular case where a is located exactly on the straight line driven by the segment s. The limit of γas when q a0 tends to zero (the y i remaining different from zero) is:

lim qa0-→0 γas = 1 4 [sign (y 2 ) -sign (y 1 )] (62)
If the integration point is located inside the segment, y 1 and y 2 have opposite signs and y 2 is positive according to our notations, so we find γas = 1/2, as expected. On the contrary, if the point is located outside the segment, y 1 and y 2 have the same sign, and we obtain γas = 0. Thus, for a point located on a straight wall, (60) gives the expected result: γ a = 1/2.

The case where the integration point is located at the intersection of two segments corresponds to a singularity. Let us consider a point a belonging to the segment s 1 and getting closer of one of its extremities r v , in the direction of segment s 2 , which makes an internal angle θ v with s 1 at the point r v (see Figure 2). Let us assume that the lengths of the two segments are large enough so that only the segments s 1 and s 2 have a contribution. What we saw before

shows that γas1 = 1/2 for any value of the distance r av > 0. Making r av tend to zero we obtain:

γvs2 = 1 2 - θ v 2π (63)
Finally:

γ v = 1 -(γ vs1 + γvs2 ) (64) = θ v 2π
which is the expected result. If the shape of the wall changes quickly close to the vertex particle v, other positive or negative contributions can appear in (56), but there is no singularity problem. In Table 2 the techniques used to compute γas in all the situations are summarised.

It was checked that the computed results perfectly match the theoretical values of γ a in cases of 1) a straight infinite wall and 2) an arbitrary angle.

Validation of the model on laminar free-surface and confined flows

The ISPH algorithm itself, without USAW boundary conditions, is relatively well established [START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF], so that we don't present any validation on cases without Table 2: Summary of the various cases for the calculation of γas in 2-D.

Particle position γas

Near a wall eqn (60)

On a segment 1 2

On a vertex v 1 2 1 - θ v 2π
walls in this paper. We start the validation with laminar flows where reference results are widely available in the literature. In Section 5, validation on two confined turbulent flows will be presented, one of them being a 2-D turbulent Poiseuille channel flow, which is the standard case for validation of the k -ǫ model. In all simulations the reference density of the fluid is ρ = 1000 kg.m -3 .

According to what was said earlier, in the following the present model will be referred to as ISPH-USAW.

Lid-driven cavity

The lid-driven cavity test-case is classical in fluid dynamics and is much used to validate numerical models. It consists of a square closed cavity of size L whose lid slides laterally at a constant velocity, driving the fluid under the effect of the viscosity. For Reynolds numbers lower than about 7500 [START_REF] Peng | Transition in a 2-D lid-driven cavity flow[END_REF], it reaches a steady-state after some time. Then, it is possible to compare the results between different computational fluid dynamics (CFD) codes. In particular, the SPH results were compared to the ones obtained by Ghia et al. [START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and multigrid method[END_REF] by a multigrid simulation method, and to the ones obtained with Code_Saturne, a widely validated code based on FV [START_REF] Archambeau | Code_Saturne: a Finite Volume code for the computation of turbulent incompressible flows -industrial applications[END_REF]. The FV simulations were always done with 512 cells. Three Reynolds numbers were considered here: 100, 400 and 1000. The Reynolds number is based on the size of the cavity and the velocity of the lid:

Re = U L ν ( 65
)
The latter is progressively accelerated from 0 to U . We define the dimensionless

variables x + = x/L, z + = z/L, u + x = u x /U , u + z = u z /U , p + = 2p/ρU 2 .
A representation of the results obtained with the present ISPH-USAW model and FV after convergence for a Reynolds number of 1000 is presented Figure 3, which qualitatively shows that the two CFD codes give very similar results.

Simulations on this test-case showed that the impermeability of the walls is granted by the ISPH-USAW model.

For the Reynolds number 100, we compared ISPH-USAW results to Yildiz et al.'s results [START_REF] Yildiz | SPH with the multiple boundary tangent method[END_REF] based on an ISPH model with the multiple boundary tangent method (ISPH-MBT). A discretization of 120 × 120 particles was used in both methods. The velocity profiles in x + = 1/2 and z + = 1/2 are shown in Figure 4, where the same quality of results was obtained with both ISPH models compared to Ghia et al. and to FV results. We could not compare pressure results since there were none available in [START_REF] Yildiz | SPH with the multiple boundary tangent method[END_REF].

For the Reynolds number 400, we compared ISPH-USAW results to WC-SPH using USAW boundary conditions (WCSPH-USAW). A discretization of 200 × 200 particles was used in both methods. For WCSPH-USAW the numerical speed of sound was taken equal to 10U , and a background pressure was imposed, without which cavities appear in the flow which is in agreement with [START_REF] Lee | Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method[END_REF]. Besides, a Ferrari density correction [START_REF] Ferrari | A new 3D parallel SPH scheme for free surface flows[END_REF] was applied, which was adapted to WCSPH-USAW by Mayrhofer et al. [START_REF] Mayrhofer | Study of differential operators in the context of the semi-analytical wall boundary conditions[END_REF]. The velocity profiles are shown on the left side of Figure 5, where the same quality of results was obtained with ISPH-USAW and WCSPH-USAW compared to Ghia et al. to FV. The pressure profiles in z + = 1/2 and on the diagonal of the cavity, defined as that between the bottom-left and the top-right corners, are shown on the right side of Figure 5. It appears that WCSPH-USAW results are much inferior to ISPH-USAW results in terms of pressure prediction, even with a Ferrari density correction.

For the Reynolds number 1000, we compared our ISPH-USAW results to WCSPH-USAW and to the results obtained by Xu et al. [START_REF] Xu | Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach[END_REF] using ISPH with a classical ghost particles technique (ISPH-GP). A discretization of 240 × 240 particles was used in all methods. The velocity profiles are shown on the left side of Figure 6, where the same quality of results was obtained with both ISPH models compared to Ghia et al. and to FV. The velocity results obtained with WCSPH-USAW are slightly inferior to the two ISPH models. Both ISPH models are much better than WCSPH in terms of pressure prediction, as can be seen in Figure 7. Finally, the computational time with ISPH-USAW was smaller than with WCSPH-USAW as shown in Table 3, and FV performed faster.

For the three Reynolds numbers ISPH-USAW results are in good agreement with the ones obtained with FV and by Ghia et al. in terms of velocity and pressure, which shows that the boundary conditions are imposed satisfactorily for laminar flows. It is expected that ISPH-MBT and ISPH-GP perform well on this test-case where the geometry is simple. Though, no convergence study was presented in the two latter works, so that the order of convergence of those models is not known.

To quantify the error made with our ISPH model compared to the FV method, convergence studies were performed at a Reynolds number of 1000 where the results obtained with FV on a cavity discretized by 512 × 512 cells were taken as a reference. The particles positions, through:

L 2 = 1 V b∈P V b u sol b,x -u ref b,x u max 2 (66)
where V = b∈P V b is the total volume of the computational domain, u sol is the velocity obtained by the ISPH model, u ref is the velocity obtained with FV and u max = U is the maximum theoretical velocity of the flow. The results of the convergence study are shown on the right side of Figure 6, where it appears that the order of convergence of ISPH-USAW is close to 2, whereas WCSPH-USAW presents a convergence order less than one and an error about 10 times higher than with ISPH-USAW.

Infinite array of cylinders in a channel

The second confined laminar flow considered in this work consists of a very viscous flow around an infinite array of cylinders confined in a channel. This case was chosen in order to check that ISPH-USAW can accurately predict hydrodynamic forces on walls. The problem considered in this work is the same as in [START_REF] Liu | Viscoelastic flow of polymer solutions around a periodic, linear array of cylinders: comparisons of predictions for microstructure and flow fields[END_REF] and [START_REF] Ellero | SPH simulations of flow around a periodic array of cylinders confined in a channel[END_REF]. A cylinder of radius R c = 0.02m is placed at the halfheight of a channel, at z = z c = 0.04m. The latter is bounded by walls on its upper and lower sides and periodic boundary conditions are applied along body force is the one proposed in [START_REF] Mayrhofer | Study of differential operators in the context of the semi-analytical wall boundary conditions[END_REF]:

F n = F n-1 + v -2ṽ n-1 + ṽn-2 δt ( 67
)
where ṽn is the average longitudinal flow velocity in the unobstructed channel at time n, computed as:

ṽn = 1 N n c a∈F ∪Ωc u n x ( 68
)
where Ω c is a slice of the channel located at x = L c of width equal to the initial interparticular spacing δr, and N c is the number of fluid particles located in this slice at time n.

The total drag force per unit length acting on the cylinder, F D , was computed for several values of L. This force is oriented along the x-direction and was computed as:

F D = s∈Γ -p s n s + µ ∇u s + ∇u T s • e x S s (69)
where Γ is the boundary of the cylinder, S s is the length of the segment s and the gradient of velocity at the segments was computed as:

∇u s = 1 2 i=1,2 G γ,- vi {u b } (70)
where the v i are the vertices linked together by segment s. For the following comparisons, the dimensionless drag coefficient will be used which is defined as

C D = F D /µ v [8]
. Figure 8 shows the values of C D obtained with ISPH-USAW compared with the results of Liu et al. [START_REF] Liu | Viscoelastic flow of polymer solutions around a periodic, linear array of cylinders: comparisons of predictions for microstructure and flow fields[END_REF] for several lengths of the channel. These results were obtained with a Finite Elements Method (FEM).

The agreement is good for the three values of L considered.

Let us now consider only the case where L = 6. A comparison of velocity profiles was done with results obtained by Ellero et al. [START_REF] Ellero | SPH simulations of flow around a periodic array of cylinders confined in a channel[END_REF] with the Immersed Boundary Method (IBM) [START_REF] Meyer | Assessment of implicit large-eddy simulation with a conservative immersed interface method for turbulent cylinder flow[END_REF][START_REF] Meyer | A conservative immersed interface method for large-eddy simulation of incompressible flows[END_REF] and with WCSPH using mirror particles to model boundaries (WCSPH-MP). For the SPH simulations, a discretization of 120 particles along the height of the channel was used. We observe that the ISPH-USAW velocity profiles match quite well the ones obtained with IBM (see . They attributed this to the fact that the discretization error becomes predominant for lower resolutions but it does not seem to be a relevant explanation since we did not observe this phenomenon in our simulations. Nevertheless, our results show that the pressure prediction is more accurate with ISPH-USAW than with WCSPH-USAW.

Note that for this test-case the numerical stability is conditioned by the viscous force, so that the time-step is the same with WCSPH and ISPH. Thus, computational times are higher with the latter. They are presented in Table 3. To reduce computational times at low Reynolds numbers with ISPH a solution would be to treat the viscous term implicitly, as was presented in [START_REF] Trask | Accuracy and performance of implicit projection methods for transient viscous flows using SPH[END_REF] for example.

Dam-break over a wedge

This case was simulated in order to check that our new ISPH-USAW model can accurately represent violent free-surface flows. It consists of a schematic 2-D dam-break in a 2 meters long and 1 meter high pool, presenting a triangular wedge in the bottom. The geometry is the same as in [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. The initial interparticular spacing for the simulations with ISPH and WCSPH was taken equal to 10 -2 m and the kinematic viscosity to 10 -2 m 2 s -1 . In the case of the WCSPH method, a Ferrari density correction was used [START_REF] Ferrari | A new 3D parallel SPH scheme for free surface flows[END_REF] and the numerical speed of sound was taken equal to 20ms is useful to check the accuracy of our method. The results obtained with VoF were considered as a reference against which the ones obtained with SPH were compared. The comparison is presented Figure 11 in a qualitative way. The dimensionless time was defined as t + = t g/H where g is the magnitude of the gravity field and H is the initial fluid depth (H = 1m). The two methods give similar results. Differences appear between the models that can be due to the two-phase nature of VoF, while the SPH models are single-phase. Moreover, in the visualisation of VoF results, the free surface is considered as the locations where the volume fraction is 0.5, which can explain some of the differences appearing in Figure 11 to capture air inside the fluid in the two-phase VoF simulation, which does not happen with SPH. On Figure 11, one can observe that a consequent number of particles remains stuck to the walls during the SPH simulation. For example, this can be seen quite well at time t + = 3.13. This is due to the high viscosity of the fluid considered here. Furthermore, particle clumping is observed at the free-surface, which is well visible on the jet. This is due to the switch off for the diffusion shift close to the free-surface as mentioned in Section 3.2. In order to quantitatively compare the different methods, the evolution of the pressure force applied on the left side of the wedge during the simulation is plotted, as in [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. This normal force F was computed by integrating the pressure on the left side of the wedge, Γ, according to:

F = s∈S∪Γ p s S s (71)
where S s is the surface of the segment s. In this case all the surfaces of the segments are equal to δr. The results obtained with ISPH-USAW, WCSPH-USAW and VoF are compared in Figure 12. The peaks that appear on the VoF curve correspond to the collapse of trapped air bubbles, which hampers the convergence of the linear solver. The three methods give similar results.

However, the evolution of the value of the force is smoother with ISPH-USAW than with WCSPH-USAW. Besides, the prediction of the maximum value of the force is closer to the one obtained by VoF with ISPH-USAW than with WCSPH-USAW. When the pressure peek occurs, the effect of air is likely to be small, so that ISPH probably predicts that peek better than WCSPH.

On the other hand, simulations on this test-case showed that the impermeability of the walls is granted by the ISPH-USAW model even in the presence of strong impact of the water on a solid wall. For the latter, the computational time was smaller than for WCSPH-USAW, as shown in Table 3. VoF presented higher computational time than the two SPH models, which also happened on the next test-case (Section 4.4).

Water wheel

A water wheel case is now proposed in order to show that the new ISPH-USAW model is able to represent flows where complex free-surface shapes and complex wall boundaries are involved. The geometry of the problem is presented Figure 13. The wheel radius R was set to 1m. The wheel turns counterclockwise at π/2 rad.s -1 , driving the fluid. The viscosity was set to 10 -2 m 2 s -1 . Thus, the Reynolds number is about 300 and it is possible to assume that the flow is laminar. The latter is periodical along x, presents a free-surface and a horizontal Figure 13: Water wheel test-case: scheme of the geometry.

bottom along z = 0. The dimensionless time was defined as t + = t g/H

where H is the water height at the initial (here H = 0.9m). As for the dambreak case, the results obtained with ISPH-USAW are compared to the VoF where the volume fraction is 0.5, which gives the impression that there is water on the paddles. This is a drawback of the VoF method where the free-surface is fuzzy. A quantitative comparison was done by comparing the time evolution of the pressure force applied on the bucket P (in red in Figure 13) obtained with the three methods. The results are presented Figure 15, where we present smoothed results for the sake of readability, since they were very noisy with the three methods. With ISPH-USAW and VoF this is explained by the fact that it is hard for the pressure solver to converge. With VoF this is due to the rotating mesh, while with ISPH-USAW it is due to the few particles wetting the wheel arms when they are above the free-surface. Although some differences appear due to the fact that we are comparing a single-phase model with a two-phase one, ISPH-USAW and VoF results are in reasonable agreement. On the other hand, with WCSPH-USAW the pressure peaks present much greater amplitudes. The amplitude of the pressure force peaks is slightly higher with ISPH-USAW than with VoF because of the presence of air trapped between the wheel and the fluid. The air pockets provide an additional pressure on the wheel, but they also reduce the water level beneath it, which in the end reduces the force due to water on the paddle. In spite of this, the results obtained with ISPH-USAW are quite satisfactory and show that the new model is robust and accurate, even with complex walls. Besides, the computational time was lower with ISPH-USAW than with WCSPH-USAW and VoF performed slower than the two SPH models, as shown in Table 3 (all codes running on one CPU).

The very high computational time exhibited by VoF on this case is due to the difficulty the pressure solver had to converge due to the rotating mesh, which led to high numbers of solver iterations. where y is the distance to the lower wall. The friction Reynolds number is defined as:

Confined turbulent flows

Re * = u * e ν (73)
It is equal to the dimensionless vertical coordinate at the centre of the channel, e + , and was taken equal to 640, so that the molecular viscosity of the fluid was taken equal to 1.5625 × 10 -3 m 2 .s -1 . The results presented below were obtained with an initial interparticular spacing of 5 × 10 -2 m.

The results obtained with ISPH-USAW are presented in Figure 16 explains why the turbulent kinetic energy profile obtained with DNS is different from the ones obtained with FV and ISPH-USAW close to the wall.

Fish-pass

Let us now consider another turbulent case, more complex and closer to reality: a water flow through a periodical fish-pass system, which is the one considered in [START_REF] Violeau | Fluid Mechanics and the SPH method[END_REF][START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. It consists of a succession of pools communicating through vertical slots. When the number of pools is high enough, the flow can be considered as periodical and it is sufficient to study one of them. Experimental results [START_REF] Tarrade | Experimental approach to adapt the turbulent flow in the vertical slot fishways to the small fish species[END_REF] showed that the mean flow is approximately parallel to the bottom of the pool, the latter being inclined of an angle I ≈ 0.1 rad compared to the horizontal. Thus, the flow was modelled in two dimensions (top-viewed) and the variations along the vertical were neglected. The effect of gravity was not taken into account and the free-surface behaviour was not represented. Thus, this flow does not represent the real one, since turbulence is a three dimensional phenomenon and the free-surface cannot remain perfectly horizontal. For a complete description of the geometry of the fish-pass, see [START_REF] Violeau | Fluid Mechanics and the SPH method[END_REF]. In our simulations the flow was driven by a constant body force along the x axis of magnitude 1.885 m.s -2 . The Reynolds number is between 10 5 and 10 6 since the molecular viscosity of the fluid considered is ν = 10 -6 m 2 .s -1 , the characteristic length is the size of the slot, 0.3m and the characteristic velocity in the fluid is close to 1m.s -1 . The results obtained with the new ISPH-USAW model were compared to the ones obtains with FV and with WCSPH-USAW. In all cases the RANS equations were solved using a k -ǫ model, as presented in section 3.1. The SPH simulations were done with 58823 particles while the simulations with FV were done with 24632 cells. A qualitative comparison of the results obtained with ISPH-USAW and FV after 20s of physical time is presented in Figure 18.

A quantitative comparison of the three methods was done by comparing bottom-left part of the wall (profile P 4 in Figure 18). The results are shown in Figure 23, where we see that ISPH-USAW improves a lot the distribution of wall pressures. Note that the differences observed between the two SPH models and FV can be due to slight differences in the imposition of boundary conditions in the k -ǫ model. On this test-case, WCSPH performed faster than ISPH and FV performed faster than the SPH models (see Table 3). In summary, the new ISPH-USAW model makes it possible to accurately represent turbulent flows presenting complex wall boundaries, while such flows are very hard to model using ghost or mirror particles, due to the accuracy required regarding the imposition of a non-homogeneous Neumann boundary condition on p and ǫ.

Conclusions

In this work a new ISPH method was proposed, in which solid boundaries are modelled through the unified semi-analytical wall (USAW) boundary con- were presented to show that ISPH-USAW is able to accurately model complex laminar and Reynolds-averaged turbulent flows, even with complex geometries.

Convergence studies were done on a lid-driven cavity. The solution obtained with FV was taken as a reference and it was observed that the convergence order was close to 2 for the new ISPH-USAW model, which shows that the wall boundary conditions are imposed satisfactorily. In general, the results obtained with the proposed ISPH-USAW model were better than with WCSPH-USAW, especially regarding the pressure prediction, and were obtained in most cases with a similar computational time. To achieve this reduction of computational time in ISPH-USAW, the wall renormalisation factor γ a was computed analytically, extending the method proposed by Feldman and Bonet [START_REF] Feldman | Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems[END_REF]. The accuracy of the k -ǫ turbulence model combined to ISPH-USAW was checked on a turbulent channel flow where an excellent agreement between our results and DNS and FV results was observed. Besides, our results were in fairly good agreement with the ones obtained with FV in the case of the fish-pass. It should be noticed that with the USAW boundary conditions it is possible to apply the ISPH method to complex geometries, not easy to handle with the traditional SPH wall treatments like ghost particles. All the results presented in this paper concerned 2-D flows, but the extension of this work to 3-D, based on [START_REF] Mayrhofer | Unified semi-analytical wall boundary conditions in SPH: analytical extension to 3-D. Numerical Algorithms[END_REF], does not present any further theoretical issues. However, it requires parallel computing for efficiency reasons.

Acknowledgements

This work was partly funded by the French Research Agency (CIFRE agreement # 2011-0264).

  et al. tested a rotational projection scheme in SPH which makes it possible to impose a non-homogeneous Neumann pressure boundary condition by imposing a homogeneous boundary condition on the dynamic pressure. However, this technique does not make it possible to impose arbitrary non-homogeneous boundary condition.

  , Ferrand et al. proposed a different way of computing the renormalisation factor together with a new formulation of the differential operators. This formulation is similar to the one Kulasegaram et al. proposed for the pressure gradient, but the boundary terms are properly represented for all the differential operators. In this framework, the imposition of boundary conditions can be done in a natural way through the boundary term of the new Laplacian operator. This was applied in [10] to the k -ǫ turbulence model where Neumann boundary conditions could be prescribed exactly on k and ǫ for the first time in SPH, the condition on ǫ being non-homogeneous.

  not present any applications of their ISPH model to 2-D or 3-D. In the present work an ISPH model is developed, in which exact arbitrary Neumann boundary conditions can be prescribed when solving the pressure Poisson equation in 2-D.

Figure 1 :

 1 Figure 1: Sketch of the different sets of entities involved in the representation of the USAW boundary conditions.

Figure 3 :

 3 Figure 3: Lid-driven cavity case for Re = 1000: comparison of the results obtained after convergence with ISPH-USAW, on the left, and with FV, on the right.

Figure 4 :

 4 Figure4: Lid-driven cavity for Re = 100. Comparison of the velocity profiles in x + = 1/2 and z + = 1/2 between ISPH-USAW, ISPH-MBT[START_REF] Yildiz | SPH with the multiple boundary tangent method[END_REF], FV and the results of Ghia et al.[START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and multigrid method[END_REF].

Figure 5 :

 5 Figure 5: Lid-driven cavity for Re = 400. Velocity profiles (top), pressure profiles in z + = 1/2 (bottom-left) and pressure profiles on the diagonal (bottomright). Comparison between FV, WCSPH-USAW and ISPH-USAW. Velocity results are also compared to Ghia et al.'s results [13].

Figure 6 :Figure 7 :

 67 Figure6: Lid-driven cavity for Re = 1000. On the left: comparison of the velocity profiles in x + = 1/2 and z + = 1/2 between ISPH-USAW, ISPH-GP[START_REF] Xu | Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach[END_REF], WCSPH-USAW, FV and the results of Ghia et al.[START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and multigrid method[END_REF]. On the right: convergence studies with ISPH-USAW and WCSPH-USAW.

Figure 8 :

 8 Figure 8: Infinite array of cylinders in a channel: dimensionless drag force as a function of the inter-cylinder distance. Comparison between ISPH-USAW and the results obtained by Liu et al. [26].

Figure 9 )

 9 Figure 9). Ellero et al. obtained slightly better velocity profiles with WCSPH-MP, which can be explained by the fact that they used a ratio h/δr = 4.5,

- 1 . 5 Figure 9 : 6 .Figure 10 :

 159610 Figure 9: Infinite array of cylinders in a channel: velocity profiles for the case L = 6. Comparison between ISPH-USAW, WCSPH-MP and IBM [8].

Figure 11 :

 11 Figure 11: Dam-break over a wedge. Comparison of the free-surface shapes and pressure fields obtained with VoF (on the left) and ISPH-USAW (on the right) at different times.

Figure 12 :

 12 Figure 12: Dam-break over a wedge. Comparison of the evolution of the pressure force applied on the left-side of the wedge between VoF (6322 cells), ISPH-USAW (5881 particles) and WCSPH-USAW (5881 particles).

Figure 14 :

 14 Figure 14: Water wheel test-case. Comparison of the free-surface shapes and velocity fields between VoF on the left and ISPH-USAW on the right at t + = 66.

Figure 15 :

 15 Figure 15: Water wheel test-case. Evolution of the smoothed pressure force magnitude applied on the bucket P . Comparison between VoF, WCSPH-USAW and ISPH-USAW.

  Two validation cases were performed to assess the performance of the k -ǫ model in the SPH incompressible formalism. Let us recall that since we use a model based on the RANS formalism, only the mean quantities of the flows are modelled, which proves sufficient in many industrial studies. A more accurate model would need, e.g. LES, but this is not the purpose of the present work.

Figure 16 :Figure 17 :

 1617 Figure 16: Turbulent Poiseuille channel flow at Re * = 640. Comparison of the dimensionless velocity profiles obtained by ISPH-USAW, DNS and FV.

Figure 18 :

 18 Figure 18: Fish-pass after 20s. Comparison of the results obtained with ISPH-USAW (top) and FV (bottom).

Figure 19 :

 19 Figure 19: Fish-pass after 20s. Mean velocity profiles on P 1 (left), P 2 (middle) and P 3 (right) obtained with FV, ISPH-USAW and WCSPH-USAW.

Figure 20 :

 20 Figure 20: Fish-pass after 20s. Pressure profiles on P 1 (left), P 2 (middle) and P 3 (right) obtained with FV, ISPH-USAW and WCSPH-USAW.

Figure 21 :Figure 22 :

 2122 Figure 21: Fish-pass after 20s. Turbulent kinetic energy profiles on P 1 (left), P 2 (middle) and P 3 (right) obtained with FV, ISPH-USAW and WCSPH-USAW.

Figure 23 :

 23 Figure 23: Fish-pass after 20s. Velocity and pressure profiles on profile P 4 obtained with FV, ISPH-USAW and WCSPH-USAW.

  ditions. The major improvement compared to a classical ISPH model is the exact imposition of a non-homogeneous Neumann boundary condition on the pressure field to solve the pressure Poisson equation, which makes it possible to prescribe the impermeability condition on solid walls. Various test-cases

  ∇p) s •n s , through which we impose a von Neumann boundary condition. Let us consider a particle v belonging to the wall. It is not a Lagrangian particle since it does not move according to the Navier-Stokes equations. Instead, its Lagrangian velocity is imposed as equal to u wall

v

. This corresponds to both impermeability and no-slip conditions. Note that in case of turbulence, the Dirichlet imposed on the velocity at vertex particles in the viscous term is not used in the pressure Poisson equation. Projecting the second part of the momentum equation

[START_REF] Meyer | Assessment of implicit large-eddy simulation with a conservative immersed interface method for turbulent cylinder flow[END_REF]

onto the normal n v to the wall in v and substituting u n+1 v by its imposed value, one obtains:

Table 3 :

 3 Computational times of the various models on several test-cases. The calculations were performed on 1 CPU.

	Model	Number of cells/particles Time
	Lid-driven cavity (Re = 1000, 60s of physical time)
	FV ISPH-USAW WCSPH-USAW	512 × 512 200 × 200 200 × 200	38 h 31 h 32 h
	Infinite array of cylinders (80s of physical time)
	ISPH-USAW	12.659e 3	10h00
	WCSPH-USAW	12.659e 3	1h30
	Dam-break over a wedge (2s of physical time)
	VoF	6.322e 3	> 1h
	ISPH-USAW	5.881e 3	20 min
	WCSPH-USAW	5.881e 3	30 min
	Water wheel (30s of physical time)
	VoF ISPH-USAW WCSPH-USAW	≈ 8e 4 ≈ 3e 4 ≈ 3e 4	5 days h 18.5 h
	Fish-pass (20s of physical time)
	FV ISPH-USAW WCSPH-USAW	≈ 2.5e 4 ≈ 6e 4 ≈ 6e 4	26 h 76 h 55 h

Turbulent channel flow

In order to test the performance of the k -ǫ model associated to ISPH, a turbulent Poiseuille channel flow was modelled. The half-height of the channel, e, is equal to 1m and periodical conditions are applied along the horizontal. An external force of constant magnitude, f = 1.0 m.s -2 , is applied. The friction velocity, u * , can be calculated by writing a balance of the forces and is equal to √ f e = 1 m.s -1 . At the initial time, the particles are aligned along horizontal lines and they remain so during the simulation, even after 100s of physical time (about 60000 iterations), with either ISPH-USAW or WCSPH-USAW. The following dimensionless variables were defined:

Appendix A. Boundary conditions imposed on ǫ in the kǫ model

The Laplacian operator involved in the diffusion of ǫ reads:

The Neumann boundary condition is applied on ǫ by imposing the terms ∂ǫ ∂n s and ∂ǫ ∂n a . Since ǫ quickly varies close to the wall the same treatment as for the pressure or velocity fields, which consists in equalling these two terms, cannot be applied. Instead, we write:

where r a ′ = 1 2 (r a + r s ). We assume that the theory of turbulent boundary layer is valid and use the theoretical relations ǫ = Considering that k slowly varies close to the wall, (A.1) can be written:

On the other hand, the Dirichlet boundary condition is imposed at the vertex particles based on a FV formulation where the Dirichlet boundary condition on ǫ was 2nd order accurate in space on an orthogonal mesh.

Let us first consider a 1D situation with the same notations as before. We use the following Taylor series expansions:

The extension to 2D is done by interpolating ǫ s based on the value of the surrounding ǫ a through:

κδr bs w bs (A.8)

Finally, the Dirichlet boundary condition is imposed through the vertex particles by writing: