Anne-Éléonore Paquier

Thibault Oudart

Caroline Le Bouteiller

Samuel Meulé

Philippe Larroudé
email: philippe.larroude@legi.grenoble-inp.fr

Robert A Dalrymple

Anne-Eléonore Paquier

3D numerical simulation of seagrass movement under waves and currents with GPUSPH

Keywords: Numerical seagrass movement, GPUSPH, Fluid structure interaction, Waves, Current

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Traditionally, to study fluid-structure interactions, structures usually are described with a Lagrangian formulation while fluids often are represented through an Eulerian formulation. The formulations generally are coupled using an Arbitrary Lagrangian-Eulerian formulation for the fluid. Numerous fluid-structure interactions have been studied according to such formulations, including valve springs [START_REF] Rugonyi | On finite element analysis of fluid flows fully coupled with structural interactions[END_REF], the interaction of compressible and incompressible fluids with structures [START_REF] Bathe | Finite element developments for general fluid flows with structural interactions[END_REF], and the absorption of hydro-elastic shocks (Le [START_REF] Tallec | Fluid structure interaction with large structural displacements[END_REF].

The Lagrangian Smoothed Particle Hydrodynamics (SPH) method is a different approach to modeling in which fluids and/or solids are represented by a set of particles.

Each solid particle has defined individual material properties and moves under the effect of the fluid, itself represented by another set of particles. The use of a Lagrangian formulation for fluids can be more convenient for certain types of applications, such as fluid-free surface studies or large movements along a fluid-solid interface. Indeed, the SPH method does not require specifically processing the fluid-free surface (no surface tension) and the two groups of particles (fluid and solid) can be tracked simultaneously.

The SPH method has already been successfully used to study interactions between fluids [START_REF] Monaghan | Gravity currents descending a ramp in a stratified tank[END_REF] and solids [START_REF] Gotoh | On the state-of-the-art of particle methods for coastal and ocean engineering[END_REF][START_REF] Gray | SPH elastic dynamics[END_REF][START_REF] Liang | Incompressible SPH Simulation of Solitary Wave Interaction with Movable Seawalls[END_REF][START_REF] Van Liedekerke | Solving microscopic flow problems using Stokes equations in SPH[END_REF][START_REF] Zhang | Smoothed particle hydrodynamics and its applications in fluid-structure interactions[END_REF].

The SPH method was developed in 1977 [START_REF] Gingold | Smoothed particle hydrodynamics: theory and application to non-spherical stars[END_REF][START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF] for applications in astrophysics. J. Monaghan (an Australian researcher in applied mathematics) contributed greatly to the development of the method. Since 1985, the method has been used for other applications. In the early 1990s, the SPH method was used for fluid simulations including free surface flows, waves [START_REF] Monaghan | Simulating free surface flows with SPH[END_REF][START_REF] Randles | Smoothed particle hydrodynamics: some recent improvements and applications[END_REF][START_REF] Dalrymple | Numerical Modeling of Water Waves with the SPH Method[END_REF], multiphase fluids [START_REF] Monaghan | SPH simulation of multi-phase flow[END_REF] and weakly compressible flows [START_REF] Morris | Modeling low Reynolds number incompressible flows using SPH[END_REF]. In the field of solid mechanics, simulations have mainly focused on the study of hypo-elastic solids [START_REF] Gray | SPH elastic dynamics[END_REF]. Recent developments in science and engineering have increased the need for modeling Fluid-Structure Interactions (FSI) to analyze multi-physics processes like retroactive loops, in which the pressure and viscous stresses of the fluid create deformation in the solid, which in turn affects stress, pressure, and velocity in the fluid.

Modeling FSI can, for instance, be useful for aerodynamics studies, biomechanics simulations, or airbag design [START_REF] Agamloh | Application of Fluid Structure Interaction Simulation of an Ocean Wave Energy Extraction Device[END_REF][START_REF] Farhat | Provable second-order time-accurate loosely-coupled solution algorithms for transient nonlinear computational aeroelasticity[END_REF]. The earliest models of stimulant-elastic solid and fluid interactions using the SPH method appeared a little more than ten years ago and focused on the laws of interaction [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF][START_REF] Amini | A new model to solve fluid-hypo-elastic solid interaction using the smoothed particle hydrodynamics (SPH) method[END_REF] and the definition of the solid interface [START_REF] Ha | Shape design optimization of SPH fluid-structure interactions considering geometrically exact interfaces[END_REF]. To date, the SPH method has not often been applied to fluid-structure interactions in environmental studies [START_REF] Zhang | Smoothed particle hydrodynamics and its applications in fluid-structure interactions[END_REF]. Here, it is proposed to use this method to evaluate the influence of submerged aquatic vegetation on waves and currents in natural hydrosystems.

Seagrasses form underwater meadows in the coastal zone. These meadows are known to attenuate waves [START_REF] Bradley | Relative velocity of seagrass blades: Implications for wave attenuation in low-energy environments[END_REF][START_REF] Fonseca | A preliminary evaluation of wave attenuation by four species of Seagrass[END_REF][START_REF] Koftis | Wave damping over artificial Posidonia oceanica meadow: A large-scale experimental study[END_REF][START_REF] Paul | Spatial and seasonal variation in wave attenuation over Zostera noltii[END_REF][START_REF] Kobayashi | Wave Attenuation by Vegetation[END_REF], and currents [START_REF] Fonseca | Influence of the seagrass, Zostera marina L., on current flow[END_REF][START_REF] Fonseca | A comparison of canopy friction and sediment movement between four species of seagrass with reference to their ecology and restoration[END_REF][START_REF] Peterson | IGF-I and IGF-II mRNA expression in slow and fast growing families of USDA103 channel catfish (Ictalurus punctatus)[END_REF][START_REF] Widdows | Effects of seagrass beds (Zostera noltii and Z. marina) on near-bed hydrodynamics and sediment resuspension[END_REF], and to significantly influence sediment dynamics [START_REF] Boscutti | Distribution modeling of seagrasses in brackish waters of Grado-Marano lagoon (Northern Adriatic Sea)[END_REF][START_REF] Madsen | The interaction between water movement, sediment dynamics and submersed macrophytes[END_REF].

Depending on the biometric characteristics of the species (shoot density, leaf area index, plant stiffness, submergence ratio) and the morphological characteristics of the meadows (extension, fragmentation), the effect on hydrodynamics will vary [START_REF] Fonseca | Flow in seagrass canopies: The influence of patch width[END_REF][START_REF] Koftis | Wave damping over artificial Posidonia oceanica meadow: A large-scale experimental study[END_REF][START_REF] Paul | Spatial and seasonal variation in wave attenuation over Zostera noltii[END_REF][START_REF] Stratigaki | Large-scale experiments on wave propagation over Posidonia oceanica[END_REF]. The nonlinear interactions of these parameters with hydrodynamics variables (wave height, wave period, water level, presence of a current) can also influence the attenuation of waves, the modification of currents and sediment dynamics [START_REF] Bradley | Relative velocity of seagrass blades: Implications for wave attenuation in low-energy environments[END_REF][START_REF] Paquier | Wind-Induced Hydrodynamic Interactions With Aquatic Vegetation in a Fetch-Limited Setting: Implications for Coastal Sedimentation and Protection[END_REF][START_REF] Paul | Spatial and seasonal variation in wave attenuation over Zostera noltii[END_REF].

Large-scale hydrosystem models can integrate the role of vegetation through a change in the roughness coefficient. This can be done empirically through the modification of Manning's n value (for example, in the ADCIRC model [START_REF] Luettich | ADCIRC: an advanced three-dimensional circulation model for shelves coasts and estuaries, report 1: theory and methodology of ADCIRC-2DDI and ADCIRC-3DL[END_REF]). In other models, such as Xbeach [START_REF] Roelvink | Modelling storm impacts on beaches, dunes and barrier islands[END_REF] or Telemac [START_REF] Hervouet | Hydrodynamics of Free Surface Flows, Modelling with the Finite element Method[END_REF], a specific vegetation roughness, as expressed by [START_REF] Mendez | An empirical model to estimate the propagation of random breaking and nonbreaking waves over vegetation fields[END_REF], can be used. The calculation of this parameter is based on the geometric and physical characteristics of the vegetation meadow. This approach was very innovative in 2004, and is now widely used in modeling. However, it does not account for vegetation flexibility, which can be very important in subaquatic vegetation, especially when canopies are sparse [START_REF] Paquier | Wind-Induced Hydrodynamic Interactions With Aquatic Vegetation in a Fetch-Limited Setting: Implications for Coastal Sedimentation and Protection[END_REF]. As seagrass meadows are declining around the world [START_REF] Waycott | Accelerating loss of seagrasses across the globe threatens coastal ecosystems[END_REF], it seems essential to better take into account vegetation flexibility in hydrodynamics modeling and to better understand the role of declining meadows or, more generally, sparse canopies. This could help to investigate seagrass meadows' limitations as coastal buffers protecting against erosion and storm surges, and could potentially enhance the protection of these declining ecosystems.

Sub-Aquatic Vegetation (SAV) is known to play a role in wave and storm surge mitigation, and, thus, in coastal protection. The impacts of different species under different conditions have already been studied (see introduction of [START_REF] Paquier | Wind-Induced Hydrodynamic Interactions With Aquatic Vegetation in a Fetch-Limited Setting: Implications for Coastal Sedimentation and Protection[END_REF] for a review of the different processes observed). SAV has the ability to attenuate waves and currents and to modify turbulence. However, no study to date has compared the impacts of the different existing species. Such studies are nearly impossible to do in the field considering the variety of environments in which the different species live.

Furthermore, experiments in flumes with real marine plants often are forbidden to protect the marine plant species. A numerical study could, therefore, help remedy the lack of in-situ studies. However, the modification in hydrodynamics generated by the plants themselves is still difficult to estimate since plant movements under waves and currents are still not well simulated and no exact numerical model exists [START_REF] Mendez | An empirical model to estimate the propagation of random breaking and nonbreaking waves over vegetation fields[END_REF][START_REF] Sánchez-González | Wave attenuation due to Posidonia oceanica meadows[END_REF]). Generally, existing numerical studies use rigid bodies to simulate plants with different roughness coefficients [START_REF] Dalrymple | Wave Diffraction Due to Areas of Energy Dissipation[END_REF][START_REF] Wallace | Effects of Seagrass on Nearshore Current and Wave Dynamics[END_REF]. However, among the numerical studies dealing with plants, [START_REF] Beudin | Development of a coupled wave-flow-vegetation interaction model[END_REF] describe a three dimensional (3D) modeling framework that simulates the interactions between waves, flow, and flexible vegetation. [START_REF] Beudin | Development of a coupled wave-flow-vegetation interaction model[END_REF] evaluated a plant drag model against field measurements in eelgrass canopies and investigated the hydrodynamics in an idealized shallow-water setting.

The research presented herein aims to reproduce plant movement under waves and/or current using the SPH method. The method was tested by building numerical two dimensional (2D) and 3D flumes with the GPUSPH model (GPUSPH official website. www.gpusph.org. Accessed: 2017-03-31.;[START_REF] Gómez-Gesteira | User Guide for the SPHysics Code v2.0[END_REF][START_REF] Hérault | SPH on GPU with CUDA[END_REF][START_REF] Hérault | Atteindre la meilleure précision dans une mise en oeuvre SPH[END_REF], then comparing the simulations with the available data.

After a brief description of the SPH method (definition of the two numerical approaches used for plants and water), several tests are presented of the FSI simulations using the proposed model and these are compared with other simulations and with experimental data (the FSI was calibrated and evaluated in the open source GPUSPH code [START_REF] Hérault | GPUSPH (Version 2.0) (Software[END_REF].

The numerical simulation results are compared with experimental results based on real or plastic seagrass meadows. The Young modulus is evaluated for real or fake plants for both the experimental and numerical approaches. The wave attenuation due to the seagrass meadow in the simulation also compared with experimental data. Finally, a first simulation also is run with field data, collected in a bay occupied by a Zostera Noltei meadow (another species of marine plant) in a large lagoon in Southern France.

The SPH method and the FSI approach

The SPH method is based on the theory of interpolation [START_REF] Benz | Smooth Particle Hydrodynamics: A Review[END_REF][START_REF] Liu | Mesh Free methods: Moving beyond the finite element method[END_REF][START_REF] Monaghan | Why particle methods work[END_REF][START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF][START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF]. Its formulation often is divided into two parts: (i) the integral representation, and (ii) the particle approximation. The equation for the integral representation generally used in SPH is given in Appendix A.

Equations for modeling the fluid

The fluid was assumed to be only slightly compressible, isothermal, and governed by the continuity equation:

where 𝜌, 𝑣 𝑖 , 𝜎 ij , 𝑎𝑛𝑑 𝑏 𝑖 are the density, velocity in i direction, stress tensor in ij directions and force vector, respectively; t is time, x i is the position in the i direction.

Stress tensor 𝜎 ij is defined as:

𝜎 ij = -𝑝𝛿 ij + 𝑆 ij ( 4
)
where 𝑝 is pressure, 𝛿 ij is the Kronecker tensor and 𝑆 ij is the deviatoric shear stress both in the ij directions.

Equations for modeling the solid

The solids are governed by the same equations as the fluid. The stress tensor for a solid can similarly be decomposed into its isotropic pressure p and deviatoric shear stress S, as follows:

𝜎 ij = -𝑝𝛿 ij + 𝑆 ij (5)
where 𝑝 = -𝜎 kk 3 ⁄ is the isostatic pressure, 𝛿 ij is the Kronecker tensor, and 𝑆 ij is the deviatoric shear stress (as for the fluid).

Equation 2 can, thus, be written as:

Dv 𝑖 Dt = 1 𝜌 𝜕𝑆 ij 𝜕𝑥 𝑖 - 1 𝜌 𝜕𝑝 𝜕𝑥 𝑖 + 𝑏 𝑖 (6)
The linear relation that connects the stress and strain tensors can be derived from the time equation, 𝑆 𝑖𝑗 . This is done by applying the Jaumann rate and gives the following corrected Hooke's law:

DS ij Dt = 2𝜇 (𝐷 ij - 1 3 𝐷 kk 𝛿 ij ) + 𝑆 ik 𝛺 jk + 𝛺 ik 𝑆 kj ( 7
)
where 𝜇 is the shear modulus, defined as:

𝜇 = 𝐸 2(1+𝜈) (8)
and where 𝐸 represents Young's modulus and 𝜈 represents Poisson's ratio. 𝐷 𝑖𝑗 is the strain rate tensor:

𝐷 ij = 1 2 ( 𝜕𝑣 𝑖 𝜕𝑥 𝑗 + 𝜕𝑣 𝑗 𝜕𝑥 𝑖 ) (9)
and 𝛺 ij is the rotation tensor:

𝛺 ij = 1 2 ( 𝜕𝑣 𝑖 𝜕𝑥 𝑗 - 𝜕𝑣 𝑗 𝜕𝑥 𝑖 ) (10)
In the case of large deformations, energy is not conserved. However, this method has the advantage of describing both the fluid and the solid dynamics in terms of pressure and velocity.

Pressure can be defined in the same way for both solids and fluid by applying the linearized equation of state:

𝑝 = 𝑐 0 2 (𝜌 -𝜌 0 ) (11)
with the speed of sound being 𝑐 0 = √ 𝜀 𝜌 0

for the fluid and 𝑐 0 = √ 𝐾 𝜌 0 for the solids.

𝜀 is the modulus of the fluid studied, 𝜌 0 a reference density, and 𝐾 is the modulus of elasticity of the solid, which can be expressed as a function of Young's modulus, as follows:

𝐾 = 1 3 𝐸 (1-2𝜈) (12)

Approximation of the motion equations through the SPH method

The fluid particles

The following equation is obtained for the velocity of the fluid particles:

Dv 𝑖 𝑞 Dt = -∑ 𝑚 𝑏 𝑏 ( 𝑝 𝑎 𝜌 𝑎 2 + 𝑝 𝑏 𝜌 𝑏 2 ) 𝛿 ij 𝛻 𝑎 𝑊 ab + 𝑔 𝑖 (13)
with 𝑔 𝑖 = gravity, m b the masse of the b particle, and the classical indices i, j, q, a, and b of the SPH method .

The solid particles

The following equation is obtained for the velocity of the solid particles:

Dv 𝑖 𝑞 Dt = ∑ 𝑚 𝑏 𝑏 ( 𝑆 ij 𝑎 𝜌 𝑎 2 + 𝑆 ij 𝑏 𝜌 𝑏 2 -( 𝑝 𝑎 𝜌 𝑎 2 + 𝑝 𝑏 𝜌 𝑏 2 ) 𝛿 ij ) .𝛻 𝑎 𝑊 ab + 𝑔 𝑖 (14)
However, Eq. 14 suffers from numerical instability and requires two terms to be introduced: one term for artificial viscosity and another term to correct tensile instability in the solid traction.

Artificial viscosity, 𝛱 𝑎𝑏 , was proposed by Monaghan [START_REF] Monaghan | Shock simulation by the particle method SPH[END_REF][START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] and is used to smooth the oscillations caused by speed and to prevent the particles from being non-uniformly distributed in space when they are too close to each other.

Dv 𝑖 𝑎 Dt = ∑ 𝑚 𝑏 𝑏 (-( 𝑝 𝑎 𝜌 𝑎 2 + 𝑝 𝑏 𝜌 𝑏 2 ) 𝛿 ij + 𝑆 ij 𝑎 𝜌 𝑎 2 + 𝑆 ij 𝑏 𝜌 𝑏 2 + 𝛱 ab 𝛿 ij ) . 𝜕𝑊 ab 𝜕 𝑎 𝜕𝑥 𝑗 + 𝑔 𝑖 (15) 𝛱 𝑎𝑏 { -𝛼𝑐 𝑎𝑏 ´𝜇𝑎𝑏 𝜌 𝑎𝑏 ´𝑣𝑎𝑏 𝑟 𝑎𝑏 < 0 0 𝑣 𝑎𝑏 𝑟 𝑎𝑏 > 0 (16)
with

𝜇 𝑎𝑏 = ℎ(𝑣 𝑎𝑏 𝑟 𝑎𝑏 ) |𝑟 𝑎𝑏 | 2 +0.1ℎ , 𝑟 𝑎𝑏 = 𝑟 𝑎 -𝑟 𝑏 , 𝑣 𝑎𝑏 = 𝑣 𝑎 -𝑣 𝑏 , 𝑐 𝑎𝑏 ´= 1 2 (𝑐 𝑎 + 𝑐 𝑏 ) (see Appendix A).
This last term introduces a numerical dissipation, which must be minimized so as not to affect the solution. [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF] proposed a compromise value between oscillation and dissipation.

Correcting for speed

To move the particles, the GPUSPH code uses the XSPH technique [START_REF] Monaghan | On the problem of penetration in particle methods[END_REF]. For each particle, a term containing the average speed of the surrounding particles is added to the original equation. Particle velocities, 𝑣 ⃗, were corrected according to the XSPH scheme to smooth out unwanted numerical peaks.

𝑣 𝑎 ⃗⃗⃗⃗⃗ ~= 𝑣 𝑎 ⃗⃗⃗⃗⃗ + 𝜃 ∑ 𝑚 𝑏 𝑝 𝑎𝑏 𝑏 (𝑣 𝑏 ⃗⃗⃗⃗⃗-𝑣 𝑎 ⃗⃗⃗⃗⃗)𝑊 𝑎𝑏 ∑ 𝑚 𝑏 𝑝 𝑎𝑏 𝑏 𝑊 𝑎𝑏 (17)
The value of 𝜃 depends on the medium: 𝜃 = 0.08 for the fluid which is increased to 1 to support high deformation for the solid.

Validation of solid elasticity

As proposed by [START_REF] Gray | SPH elastic dynamics[END_REF] and [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF], the sensitivity of the artificial viscosity and the XSPH factor can be tested by simulating the free oscillations of an elastic plate having one end free and the other held in place. The plate is initially horizontal (Fig. 1a) and an initial velocity distribution is assigned according to the analytical expression of the free oscillations of a thin plate:

𝑣 𝑧 = 𝑣 𝑙0 𝑐 0 𝑓(𝑥) 𝑓(𝐿) (18)
where L is the length of the elastic plate and 2). The small difference between the results can be attributed to several factors, but it's difficult to compare 2D simulations and the current 3D simulations, since the mesh is different. The difference between the current results and those presented by [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF] and [START_REF] Gray | SPH elastic dynamics[END_REF] probably come from the current use of 3D simulation.

As proposed by [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF], this experiment allowed the artificial viscosity coefficient, 𝛼, to be tested. In agreement with [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF], the current numerical experiment showed that if 𝛼 = 0, the simulation is unstable; when 𝛼 > 1 the free oscillations damping reduces. Therefore, it was decided to use 𝛼 = 1, which caused some slight energy dissipation but guaranteed good stability. [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] proposed a tensile instability corrective formula specifically for 2D simulations, but which does not seem to be needed for 3D simulations. It can be replaced by simply increasing the XSPH coefficient. For 0.5 < 𝜃 < 2.5, similar results always were obtained. For simulations with higher initial velocities of the free end of the plate (𝑣 𝑙0 = 0.05 and 𝑣 𝑙0 = 0.1), good high oscillations without breakage were obtained for 𝜃 = 1 but the plate broke or become unstable for lower 𝜃 values (see Fig. 1b), and permanent deformations occurred for higher 𝜃 values.

Fluid-Solid Interactions

To simulate seagrass motion under waves and wave induced-currents, solid particles must be defined to simulate the plants (a few millimetres thick) and the flume (a few meters long). It was decided to create different particle sizes depending on the nature of the element (water, interface, solid (see Fig. 1c)) and to adjust some interaction laws.

In order to use smaller particle sizes for the solid, interface particles were constructed from the solid particles. They were then initialized with the sum of the masses of the neighboring solid particles. The interface particles interact with the fluid particles according to the different laws previously presented. The stresses on fluid particles are directly transmitted to the interface particles, and for solid particles, the stresses are transmitted proportionally to the mass ratio between the interface particles and the solid particles. This treatment prevents a solid particle from being subjected to significant stress due to its low mass. It also makes it possible to smooth and evenly distribute the force on the solid particles.

Three different forces were tested to determine which one had the best fit for the current simulations. The computational methods developed for FSI simulations with deformable structures in the context of particle methods can be classified into several categories, including FSI solvers. The latter were developed by coupling particle methods with other computational methods or by integrating Lagrangian particles into solvers for both the fluid and the solid. Several examples can be found in the reviews on SPH methods done by [START_REF] Gotoh | On the state-of-the-art of particle methods for coastal and ocean engineering[END_REF] and [START_REF] Zhang | Smoothed particle hydrodynamics and its applications in fluid-structure interactions[END_REF].

Following [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF], two types of boundary conditions were tested. First, the Lennard-Jones boundary condition as used by [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF] was tested. The well-known Lennard-Jones formulation is a purely repulsive force principally used for molecule interactions. The force, 𝑓(𝑟), per unit mass between a solid and a fluid particle separated by a distance 𝑟 has the form:

𝑓(𝑟) = 𝐷 (( 𝑟 0 𝑟 ) 𝑝 1 -( 𝑟 0 𝑟 ) 𝑝 2 ) . 𝑟 𝑟 0 𝑖𝑓 𝑟 > 𝑟 0 𝑓(𝑟) = 0 𝑖𝑓 𝑟 ≤ 𝑟 0 (20)
The constants 𝑝 1 and 𝑝 2 (Eq. 20) must satisfy the condition 𝑝 1 > 𝑝 2 to be repulsive; we set 𝑝 1 = 12 and 𝑝 2 = 6, were selected as proposed by [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF].

The mass of an interface particle is calculated as the sum of the masses of the nearest solid particles within its radius of action (a solid particle being bound to a single interface particle). The interaction is then calculated between the interface particles and the fluid particles, and the force calculated on the interface particles then is transmitted to the linked solid particles, taking into account the mass ratios. The distance used in the interaction laws between the interface particles and the fluid particles depends on the size ratio between the interface particles and the solid particles at initialization so that the interaction volume equals the volume of the solid particles.

Coefficient D depends on the water depth, 𝐻, and it is usually taken as 𝐷 = 5𝑔𝐻. D is applied to all the particles, and, therefore, does not generate any difference between particle types.

A formulation proposed by [START_REF] Kajtar | SPH simulations of swimming linked bodies [J][END_REF] also was tested, i.e., the same correction, as explained for (Eq. 20), was applied to the length scale. Both boundary forces yielded similar results in the calculation of solid and fluid interactions.

Finally, the second boundary condition formulation that was tested concerned extending the fluid and solid momentum and continuity equations to all particles regardless of their nature. This formulation gave good results for the experiments of [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF]. However, in cases with large differences in density between the fluid and the solid, when impacts occurred, this formulation created local high pressure. Studies have questioned or perfected these different approaches to boundary conditions [START_REF] Crespo | Boundary Conditions Generated by Dynamic Particles in SPH Methods[END_REF][START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. In the current case, the standard options in the GPUSPH code worked well enough for the current needs.

Calibration and comparison with experimental data

Antoci et al. ( 2007) did a laboratory experiment and compared the results with SPH simulations. The experiment consisted in deforming an elastic bundle acting as a dam in a water column. The configuration of the experiment is described in [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF]. Table 1 lists the size experiment parameters. The elastic properties of the beam can be determined from [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF]. For the current simulations, a Young's modulus that provided the results closest to the experimental results was chosen. Liquid and elasticity properties are listed in Table 2.

To calibrate and validate the fluid-structure interaction (FSI) model used for plant movement, the current results were compared for the vertical (Fig. 3a) and horizontal movements (Fig. 3b) at the end of the plate with [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF]. The numerical results are consistent with the experimental results found by [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF].

A fluid compressibility below the theoretical value was used to reduce the time step, as was done by [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF]. Despite this difference in compressibility, the displacements obtained were within the expected range.

In the current study, the Lennard-Jones [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF], and [START_REF] Kajtar | SPH simulations of swimming linked bodies [J][END_REF] boundary conditions were used, and better results were found with the latter formulation. The proposed multi-scale approach also was validated with [START_REF] Kajtar | SPH simulations of swimming linked bodies [J][END_REF] boundary conditions for solid particles until they were ten times closer together than the fluid particles and with an interface particle similar to the one described in the previous section. A multi-scale approach such as this is needed to simulate small elastic solids in a fluid domain.

Results and discussion

The plant movements predicted with the proposed SPH simulations were compared with both experimental results (Le [START_REF] Bouteiller | Sediment transport and shear stress partitioning in a vegetated flow[END_REF][START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] and field observations [START_REF] Paquier | Wind-Induced Hydrodynamic Interactions With Aquatic Vegetation in a Fetch-Limited Setting: Implications for Coastal Sedimentation and Protection[END_REF].

Different species of seagrass colonize the shallow coastal waters around the world, often forming large underwater meadows. Seagrass meadows are known to reduce wave height and current velocity near the seabed, but the simulation meadow (see Fig. (2019) worked on real Zostera noltei.

Comparison using the initial 2D approach

The comparison was done on the basis of [START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] study. Plant modeling relies on solid knowledge of the physical properties of the plant (height and thickness but also density and elastic modulus). In addition, it is complicated to numerically simulate a very thin leaf. For the species Posidonia Oceanica, [START_REF] Folkard | Hydrodynamics of model Posidonia oceanica patches in shallow water[END_REF] gives clear leaf properties (Table 3). Poisson's ratio is not given, but generally ranges from 0.3 to 0.4 for this type of solid.

Several limitations were encountered during the simulation process of the Posidonia Oceanica leaf thickness. First, the thickness of the plant leaf is naturally 0.1 mm and the particle size used for the simulations was 5 mm (a ratio of 1/150). To deal with this issue two choices were available: either divide the particle size by 150 or modify the plant properties while conserving a comparable flexibility. The 2D code as used originally [START_REF] Oudart | SPH model to simulate movement of grass meadow of Posidonia under waves[END_REF] does not accept multi-sized particles, which would induce a decrease in the water particle size. Furthermore, reducing the particle size would have dramatically increased the number of particles, and, thus, the computation time. The memory needed to run such a calculation is much too large to allow the use of the SPH method. Therefore it was decided to modify the plant properties. Plant thickness, Young's modulus and fluid compressibility were modified.

The impact of density was not investigated; and another study should be done on this point. The final modified parameters of the plants are listed in Table 3.

Some seagrass species' leaves (including Posidonia Oceanica) grow from shoots in small groups (the number of leaves per shoot depends on the species and can vary from one shoot to another). These shoots cannot be simulated in a 2D model. However, several leaves close to each other were simulated in order to represent a natural meadow as well as possible. Figure 5 shows the movement of three simulated leaves (as a first approach to meadow simulation) under a breaking wave. The numerical representation of the deformation of the plants in a meadow is therefore possible. [START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] did laboratory experiments to study the flow-induced reconfiguration of flexible subaquatic vegetation. The plant mimics used were inspired by Posidonia Oceanica plants; the drag and posture of the mimics were within the natural range for stiffness and buoyancy. [START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] then compared the experimental results with a theoretical model. Two current velocities were tested: one lower and one higher than the restoring forces, the latter flow bending the leaves. The current simulation was compared to [START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] experiments from a qualitative point of view. [START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] built two kinds of leaves:

-one made of silicone foam (Young's modulus, E = 500 kPa;  = 330 kg/m 3 with a seawater density assumed to be 1025 kg/m 3 ; thickness = 1.9 mm);

-one made of high-density polyethylene (HDPE) (Young's modulus, E = 0.93 GPa;  = 50 kg/m 3 ; thickness = 0.4 mm).

Figure 6a shows [START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] experimental set up and the comparison with the current simulations (Figs. 6b,c and Figs. 6d,e). For both current velocities, the SPH simulation gives the correct inclination of the plant. For the low current velocity (U = 0.16 m/s), an inclination angle of 45° was obtain versus 38° angle obtained in the experiment with silicon foam leaf blades (see Figs. 6b,c).

For the higher current velocity (U = 0.32 m/s), [START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] measured an angle of 21° while an angle of 26° was simulated. The shape of the numerically simulated blade is in reasonable agreement with the experimental results shown in black on the picture (see Fig. 6e).

In the experiment, [START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] also built artificial leaf blades from HDPE.

Again, the SPH simulation showed results similar to the experimental results. The difference [START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] observed in the angle taken for the HDPE leaf compared to the silicon leaf is well represented by the SPH modeling.

Comparing simulated 2D and 3D plant effects on waves with experimental results

The simulated 2D results were compared with the experimental results of [START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] and [START_REF] Stratigaki | Large-scale experiments on wave propagation over Posidonia oceanica[END_REF]; both researchteams worked on Posidonia Oceanica mimics. The simulated 3D results were compared with those of [START_REF] Paquier | Wind-Induced Hydrodynamic Interactions With Aquatic Vegetation in a Fetch-Limited Setting: Implications for Coastal Sedimentation and Protection[END_REF], who worked in the field on real Zostera noltei plants (a smaller species than

Posidonia Oceanica).

The 3D modeling approach used in the current study allowed simulation of plant movements under waves. The results for different h 0 /h plant ratios (where h 0 is the reference water level in the channel and h plant is the height of the plant) were compared for one plant located at a distance of x = 1.2 m (Fig. 7) from the piston paddle. The height of the wave generated by the piston paddle was influenced at x = 0.6 m from the paddle. This boundary condition is too close to the seagrass meadow and the computational domain had to be increased [START_REF] Oudart | SPH model to simulate movement of grass meadow of Posidonia under waves[END_REF][START_REF] Oudart | Two numerical approaches: 2D and 3D SPH model to simulate extreme waves over a barrier island[END_REF].

The wave heights at x = 1 m and x = 1.3 m were nearly identical for all the simulations: wave height increased at the approach to the plant and decreased significantly afterwards. However, when approaching the area of the channel slope (at x = 1.7 m), there was another increase in wave height, which varied according to the ratio h 0 /h plant . The presence of the plant generated a decrease in the wave height at this location, due to a loss of energy resulting from friction. Reflecting this situation numerically would require complementary simulations with the foot of the beach slope far from the end of the sea grass meadow to avoid any slope influence on wave attenuation.

Different scenarios were simulated: a channel without plants (Case 0), with one elastic plant (Case 1), with one rigid plant (Case 2), and with three elastic plants located at x = 1 m, x = 1.2 m, x = 1.4 m from the wave maker (Case 3). Plant height was set at 0.1 meter (h 0 /h plant= 0.5 m, h 0 = 0.2 m); the changes in wave height relative to h 0 along the cross section is represented in Fig. 7. In the presence of a rigid plant (Case 2), wave height decreased then increased again after the plant due to the recirculation generated by the rigid obstacle. With a three-plant meadow (Case 3), the decrease in wave height was higher than with a single plant (Case 2) and significantly higher than without any plants (decrease 24% greater for Case 3 than for Case 0). [START_REF] Stratigaki | Large-scale experiments on wave propagation over Posidonia oceanica[END_REF] observed a decrease of 30 to 35% in wave height in their experimental study. These higher rates can be explained by a longer friction zone in the experiment than in the numerical simulations shown in Fig. 7. [START_REF] Antoci | Numerical simulation of fluid-structure interaction by SPH[END_REF][START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF][START_REF] Stratigaki | Large-scale experiments on wave propagation over Posidonia oceanica[END_REF]. These experiments are consistent with a 2D approach because very little flow can pass on each side of the seagrass. Indeed, the plant mimics in the foregoing studies were as wide as the channel or the meadow took up the whole width of the experimental set up. Obviously, to simulate experiments with sparser seagrass meadows or marine plants from coastal zones other than the ones studied by [START_REF] Stratigaki | Large-scale experiments on wave propagation over Posidonia oceanica[END_REF] and [START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF], a 3D simulation is necessary. The proposed multi-scale 3D FSI simulations with the SPH method are able to reproduce the shape and phasing of leaf movements resulting from the wave passing over the plant or plant meadow. This is the first step toward developing a more quantitative approach. It also was possible to simulate this movement for an entire 3D meadow of 200 plants (Fig. 8a) and for a 3D meadow with two separate zones of seagrass (Fig. 8b).

With a 3D method, different particle sizes can be used for water and plants. This makes it possible to (i) create plants with a more realistic thickness, (ii) to simulate several plants (200 in the current case), and (iii) to simulate several densities of seagrass meadow with realistic bathymetry (see Fig. 9). Fig. 9d shows the numerical results with and without a plant meadow on a real bathymetry (the "Point de Berre" beach in Berre Lagoon, France, [START_REF] Paquier | Wind-Induced Hydrodynamic Interactions With Aquatic Vegetation in a Fetch-Limited Setting: Implications for Coastal Sedimentation and Protection[END_REF]). This first qualitative approach with the 3D SPH method clearly shows the attenuation of the waves by the plants. The simulated plant movements (Figs. 9a, b andc) are also qualitatively consistent with those observed in the study area (video available in the supplementary material in [START_REF] Paquier | Wind-Induced Hydrodynamic Interactions With Aquatic Vegetation in a Fetch-Limited Setting: Implications for Coastal Sedimentation and Protection[END_REF].

The next research step will be to calibrate the model with in-situ wave measurements.

Comparison of simulated flow properties (Manning's coefficient) with experimental results

The aim of the TEC21 post-doc program at the University of Grenoble was to numerically reproduce the experiment of Le [START_REF] Bouteiller | Sediment transport and shear stress partitioning in a vegetated flow[END_REF], which had been done on mimics of Zostera marina in an inclined flume (15 m long and 1 m wide) located in the River Dynamics Laboratory at Simon Fraser University in Canada (Figs.

10a,b). The central section of the flume between 5 and 11 m was covered with a vegetated bed (see Fig. 10a) consisting of staggered artificial blades. The blades were made of a low density polyethylene film (920 kg/m 3 ) with a Young's modulus of 2x10 8

Pa, a length of 25.5 cm, a width of 0.75 cm and a thickness of 0.2 mm (see Fig. 10b).

The material and dimensions of the blades were chosen for their similarity to natural Zostera marina plants to maintain the ratio of flexural rigidity to buoyancy forces [START_REF] Ghisalberti | Mixing layers and coherent structures in vegetated aquatic flows[END_REF]. Two densities of plants were used in the experiments: 130 blades (low density) and 800 blades (high density) per square meter.

To simulate the interactions between the low-density meadow and hydrodynamics, several lines of plants were built. For the low-density meadow experiment, the lines were 6.25 cm from each other and the plant mimics were spaced 12.5 cm apart on each line. The alignment of the plants on one line was shifted by 6.25 cm from the previous line to form an alternating network. The nearest surrounding plants were, therefore, around 8.8 cm away. For the high-density meadow experiment, the lines were separated by 2.5 cm while the plant mimics on each line were affixed every 5 cm. Again, an alternating network was built. The nearest surrounding plants were, therefore, around 3.56 cm away. In the simulation presented here, the low-density meadow case was considered (see Fig. 10c). The simulated movements of the plants correspond to what was observed during the actual flume experiment (see Figs. 10b,c,and d), even if a better quantitative evaluation of the simulated plants' movements is still needed. A Manning's n coefficient similar to the one found experimentally by Le [START_REF] Bouteiller | Sediment transport and shear stress partitioning in a vegetated flow[END_REF] based on a uniform flow hypothesis was evaluated. For the low-density meadow, Le [START_REF] Bouteiller | Sediment transport and shear stress partitioning in a vegetated flow[END_REF] found a Manning's n coefficient value of 0.023, which remained constant under increasing flow velocity. In the current study, a

Manning's coefficient equal to 0.025 was obtained. A better parametrization of the Manning's n evaluation is still needed but it does seem appropriate to use the GPUSPH simulation to calculate the Manning's n coefficient. The evaluation of the Manning's n following the SPH simulation gave good results for the mesh (number of particles) used in the test cases. In future studies, the mesh should be increased to refine the results.

What is pointed out in the current study is that it is feasible to use the same Manning's n evaluation methodology for 3D SPH simulations and real experiments.

Conclusions

The results obtained for the various FSI simulations with the SPH method were reasonable, even though improvements can be made in the future. First, in a deformation test case of an elastic plate under a column of water, it was shown that it is possible to model the movement of a solid submitted to water effects with the SPH method and FSI. Plant movements also were successfully simulaed in a wave flume.

The different curves taken by the plant under the passing waves show the role of the plant in wave height modification. In these initial simulations, the decrease in wave height was significant. The influence of plant flexibility (in opposition to a rigid body) also was investigated. Further numerical experiments are needed to include larger and/or higher meadows that have been studied in physical experiments and in-situ in the literature.

Current modeling limitations include the sizing problems of the particles related to the cost in computation time but these limitations will undoubtedly be overcome in the future. Even so, the SPH method is still time-consuming in terms of computation for simulations over large areas. It is, however, a good way to evaluate an equivalent viscosity that can be applied in other methods and models. Once wave damping in-situ can be quantitatively simulated, researchers will be able to simulate the role of any plant and to determine an equivalent viscosity depending on each plant's properties.

Improvements will likely be made to the methodology presented here for comparing SPH models with existing experiments and the computational time required for 3D calculations will tend to decrease with better computer technology. For the time being, however, it is difficult to do quantitative comparisons between simulations and existing experimental results due to the differences in the initial parameters themselves and to the limitations of the SPH method (especially regarding computational time). The next step in this research will be to build simpler theoretical cases (for which the dimensioning of the elements and the quantification of the fluid and solid movements are known) through SPH simulations and to evaluate the need to add a coefficient to the law of interaction.

It was also possible to simulate long plastic blades with high flexibility based on the experimental flume data from Le [START_REF] Bouteiller | Sediment transport and shear stress partitioning in a vegetated flow[END_REF]. This allowed a

Manning's n coefficient to be calculated consistent with the one experimentally evaluated by Le [START_REF] Bouteiller | Sediment transport and shear stress partitioning in a vegetated flow[END_REF]. The next step in this study will be to ensure the quantification of this numerical reproduction and to improve the accuracy of the estimated Manning's n roughness coefficient for a variety of seagrass species. There is a real need for these values for the coastal management community's use in integrated models.

Interest in the role of nature-based defences is growing (Bridges et al., often protected, marine plants; therefore, the authors hope that the virtual FSI channel will contribute to solving these environmental objectives. Finally, a layer of sediment [START_REF] Ghaïtanellis | A sph elastic-viscoplastic model for granular flows and bed-load transport[END_REF] or gravel [START_REF] Kazemi | SPH modelling of turbulent open channel flow over and within natural gravel beds with rough interfacial boundaries[END_REF] could also be taken into account, which, when added to the plant/wave interaction we simulated in our study, would be able to numerically reproduce a complete coastal system.

  𝑓(𝑥) = (𝑐𝑜𝑠(𝑘𝐿) + 𝑐𝑜𝑠ℎ(𝑘𝐿))(𝑐𝑜𝑠ℎ(𝑘𝑥) -𝑐𝑜𝑠(𝑘𝑥)) + (𝑠𝑖𝑛(𝑘𝐿) -𝑠𝑖𝑛ℎ(𝑘𝐿))(𝑠𝑖𝑛ℎ(𝑘𝑥) -𝑠𝑖𝑛(𝑘𝑥)) (19) while 𝑣 𝑙0 = 0.01 determines the initial velocity of the free end and 𝑘𝐿 = 1.875 is the fundamental mode. For a plate with the following features: 𝐿 = 0.2 m, 𝐻 = 0.02, m ℎ = 0.2 m, 𝜌 = 1000 kg/m 3 , 𝐾 = 3.25 × 10 6 N/m 2 , 𝜇 = 715000 N/m 2 , the analytical non-dimensional amplitude (A) and period (T) of the first oscillation are A/L = 0.115 and 𝑇 𝑐 0 𝐿 ⁄ = 72.39. Gray et al. (2001) obtained A/L=0.125 and 𝑇 𝑐 0 𝐿 ⁄ = 82 and Antoci et al. (2007) obtained A/L=0.124 and 𝑇𝑐 0 𝐿 ⁄ = 81.5. For the proposed model the following was obtained A/L=0.116 and 𝑇 𝑐 0 𝐿 ⁄ = 73.34 (Fig.

  4), could be also a problem for the simulations. The different datasets used for the current comparisons were collected on real or mimicked plants of different species.[START_REF] Luhar | Flow-induced reconfiguration of buoyant and flexible aquatic vegetation[END_REF] and[START_REF] Stratigaki | Large-scale experiments on wave propagation over Posidonia oceanica[END_REF] worked on mimics of Posidonia oceanica, Le[START_REF] Bouteiller | Sediment transport and shear stress partitioning in a vegetated flow[END_REF] workedon mimics of Zostera marina and Paquier et al.

Fig. 1 .

 1 Fig. 1. a) Set up of the plate, b) example of high deformation in a longer plate, and c)

Fig. 2 .

 2 Fig. 2. Comparison between the 3D numerical simulation and the theoretical vertical

Fig. 3 .

 3 Fig. 3. Comparison between numerically simulated (sim) and experimental (exp)

Fig. 4 .

 4 Fig. 4. Seagrass meadow of Posidonia Oceanica (left) and (right) schematic SPH

Fig. 5 .

 5 Fig. 5. Movement of three plants under a breaking wave using FSI simulated with SPH.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Comparison between the experiment with artificial plants made of foam in Luhar

Fig. 8 .

 8 Fig. 8. Movement of the sea grass meadow under waves at different depths: a) with one

Fig. 9 .Fig. 10 .

 910 Fig. 9. GPUSPH 3D simulation of the real bathymetry in the Berre Lagoon: a), b) and c)

Table 1 .

 1 Dimensions in the experiment.

	Dimensions in the experiment
	A(m)	0.1
	H(m)	0.14
	B(m)	0.1
	B*(m)	0.098
	L(m)	0.079
	s(m)	0.005

Table 2 .

 2 Fluid and solid properties and elastic modulus used

	Fluid properties
	Density, ρ(kg/m 3 )	1,000
	Compressibility modulus, ε (N/m 2 )	2x10 6
	Solid properties
	Density, ρ s (kg/m 3 )	1,100
	Young's modulus, E (N/m 2 )	1.2x10 7
	Shear's modulus, µ (N/m 2 )	4.27x10 6
	Bulk modulus, K (N/m 2 )	2x10 7

Table 3 .

 3 Properties of real Posidonia Oceanica and modified properties used for numerical simulations, from[START_REF] Folkard | Hydrodynamics of model Posidonia oceanica patches in shallow water[END_REF]

	Property	Real Plant	Simulated Plant
	Density, ρ s (kg/m 3 )	910 ± 110 kg/m 3	910 kg/m 3
	Young's modulus, E	4.7 ± 0.6 × 10 8 N/m 2	2.1 × 10 4 N/m 2
	(N/m 2 )
	Thickness (mm)	0.1 mm	15 mm
	Bulk modulus, K (N/m 2 )	3.92 ± 0.6 × 10 8 N/m 2	1.7 × 10 4 N/m 2
	Shear's modulus, µ (N/m 2 )	1.81 ± 0.5 × 10 8 N/m 2	8.7 × 10 3 N/m 2

Acknowledgements

The authors gratefully acknowledge the support of the NVIDIA Corporation who donated the GTX 780 GPU used for this research. This research was partially supported by the LabEx Tec 21 program (Investissement d'Avenir -grant agreement ANR-11-LABX-0030) and by

financial support from the "Agence de l'Eau RM&C" through the CANOPé research program.

AEP was partially supported by the French National Research Agency (ANR) through ANR @RAction chair medLOC (ANR-14-ACHN-0007-01project leader Thomas Stieglitz). The authors also thank M. Luhar and H. Nepf (from MIT) for the exchange of information and data.

APPENDIX A:

The integral representation of an arbitrary function, 𝐹(𝑟), is defined as follows:

where 𝛿(𝑟 -𝑟′) represents the Dirac delta function, 𝑉 is the area of integration (in two dimensions) and 𝑟 is the position vector. For the SPH method an approximation to the delta function is:

where ℎ is the smoothing length defined by the area of influence of the smoothing function (or weight function), 𝑊.

In the SPH method, the system is represented by a number of particles, each carrying a mass and other fluid parameters and defined in a specific area of space. The particle approximation is used to discretize the domain. If 𝐹(𝑟′) is known in 𝑁discrete points 𝑟 1 , 𝑟 2 ,...𝑟 𝑁 , then it can be approximated by:

where 𝑏 corresponds to a particle neighbor. The kernel function must satisfy [START_REF] Benz | Smooth Particle Hydrodynamics: A Review[END_REF][START_REF] Liu | Mesh Free methods: Moving beyond the finite element method[END_REF][START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF]:

The first condition is the normalization condition, the second is the Delta function property (observed when the smoothing length approaches zero), and the third is the compact condition where 𝜅 is a constant related to the length of smoothing and which defines the effective area of the smoothing function.

The stability of the SPH algorithm depends heavily on the second derivative of the kernel [START_REF] Morris | Modeling low Reynolds number incompressible flows using SPH[END_REF]. There are different kernel functions, the most well-known being a Gaussian function. In the current study, the Wendland function primarily was used:

, and 𝛹 is the normalization constant, equal to 10 7 ⁄ 𝜋 for two dimensions.

The accuracy of the SPH method can be improved by increasing the number of neighbors considered, which also significantly increases computational time. In the case of 2D modeling, a number of neighboring particles between 20 and 30 is a good compromise between accuracy and computational time; this number needs to be increased to 80 for 3D modeling.

There are two forms used to calculate a gradient of a scalar F in SPH. One is mainly used when the particle density varies only slightly: