Simulation of Water-Soil-Structure Interactions Using Incompressible Smoothed Particle Hydrodynamics

: In the present work, an incompressible smoothed particle hydrodynamic (SPH) method is introduced to simulate water-soil-structure interactions. In the current calculation, the water is modelled as a Newtonian fluid. The soil is modelled in two different cases. In the first case, the granular material is considered as a fluid where a Bingham type constitutive model is proposed based on Mohr-Coulomb yield-stress criterion, and the viscosity is derived from the cohesion and friction angle. In addition, the fictitious suspension layers between water and soil depending on the concentration of soil are introduced. In the second case, Hooke’s law introduces elastic soil. In ISPH, the pressure is evaluated by solving the pressure Poisson equation using a semi-implicit algorithm based on the projection method and an eddy viscosity for water is modelled by a large eddy simulation with the Smagorinsky model. In the proposed ISPH method, the pressure is stabilized to simulate the multiphase flow between soil and water. Numerical experiments for water-soil suspension flow of Louvain erosional dam break with flat soil foundation, is simulated and validated using 3D-ISPH method. Coupling between water-soil interactions with different solid structures are simulated. The results revealed that, the suspension layers with the Bingham model of soil gives more accurate results in the experiment as compared to the case of the Bingham model without suspension layers. In addition, the elastic soil model by the Hooke’s law can simulate soil hump accurately as compared to the Bingham model. From the simulations, avoiding erosion behind the structure for preventing the structure break during flood are investigated by using an extended structure or a wedge structure.

the effective stress, that is Terzaghi's concept of effective stress. They used this approach to model the saturated soil problem with SPH method. On the other hand, Ulrich et al. [Ulrich, Koliha and Rung (2011)] modelled the water/soil interaction by SPH method as a multiphase flow and the layers between water and soil are introduced as fictitious suspension layers, which deserve more attention in several cases. In the current work, we follow Ulrich et al.'s [Ulrich, Koliha and Rung (2011)] technique for treating the interaction between water, soil and suspension layers. The fluids are assumed to be Newtonian and an eddy viscosity is modelled by means of large eddy simulation using the Smagorinsky model. The soil model considers the granular material as a fluid where a Bingham type constitutive model is proposed based on Mohr-Coulomb yield-stress criterion and the viscosity is derived from the cohesion and friction angle. A concentration based approach to mimic the stresses inside the fictitious suspension layer is introduced which is derived from a Chezy-relation between the shear stresses and the local flow velocity as proposed by Fraccarollo et al. [Fraccarollo and Chapart (2002)]. Recently, Ren et al. ] derived dual-support SPH (DS-SPH) in the field of solid mechanics. The main advantage of DS-SPH method is the easy formulation of tangent stiffness matrix. Dai et al. [Dai, Ren, Zhuang et al. (2017)] introduced different support domain and dual-support of SPH method for treating elastic mechanics. Ren et al. ] developed dual-support smoothed particle hydrodynamics (DS-SPH) method with variable smoothing length. The current DS-SPH is applied into weakly compressible flow including water droplet flow and 2D dam break over dry bed. It is reported that DS-SPH can reduce the computational cost compare to conventional SPH. Ren et al. ] presented the dual-horizon peridynamics (DH-PD) formulation for simulations of crack paths with variable horizon and particle sizes. Ren et al. [Ren, Zhuang, Cai et al. (2016)] developed dual-horizon peridynamics which includes variations on the horizon sizes. Rabczuk et al. ] proposed the novel nonlocal operator theory based on the vibrational principle for solving the partial differential equations. The novel nonlocal operator can stabilize meshes methods by the dual-concept and allows implementation of complete implicit methods. In the present study, the pressure is stabilized by introducing the source term which contains both contributions from velocity-divergence free and density invariance conditions Sonoda (2011a, 2011b); Asai, Aly, Sonoda et al. (2012)]. In addition, the eddy viscosity based on the Smagorinsky model is introduced. To simulate soil hump for seawall, Bingham flow approach with suspension region and solid approach for soil, are integrated, then the technique is applied to a seawall collapse simulation during a tsunami. The modification in the original Bingham flow model [Ulrich, Koliha and Rung (2011)], which is based on weakly compressible approach for water, is the extension to the incompressible formulation by ISPH. First, the water-soil suspension flow of Louvain erosional dam break [Fraccarollo and Chapart (2002)] with flat soil foundation is simulated using 3D-ISPH method. This simulation is validated by comparing it to the experimental results. Second, several numerical tests for fluidstructure soil foundation interactions are discussed.

Mathematical analysis
The governing equations of transient compressible fluid flow include the equations for conservation of mass and momentum are: where subscripts α and β refer to the spatial coordinates, t is the time, g is the gravitational acceleration, v is the velocity vector.

Fluid model
For the simple Newtonian fluids, the stress tensor σ is given by: where p is the pressure, αβ δ is the unit tensor and αβ τ is the stress tensor and represents the viscous stresses which depend on an isotropic viscosity * µ and gradient of velocity as: where, ϑ is viscosity coefficient and strain rate-tensor is defined as: And the effective dynamic viscosity * µ is composed from a viscosity µ and an eddy viscosity T µ as: In this paper, it is assumed that the eddy viscosity is modeled by the static Smagorinsky For incompressible fluids . 0; ∇ = v the viscous stress tensor is: The total stress tensor is therefore: Note that, in the most general incompressible flow approach, the density is assumed by a constant value with its initial value 0 ρ . Then, the governing equations for an incompressible Newtonian fluid are summarized:

Projection method
In the projection method [Cummins and Rudman (1999)], the velocity-pressure coupling problem has been solved separately for velocity and pressure. Here, all the state variables may update from a previous time step to current time step. In this below, superscripts (n) and (n+1) indicate previous and current time step, respectively. In the first predictor step, the intermediate distate without pressure gradient is assumed and its velocity field is indicated by * v . The intermediate velocity field can be evaluated by solving the following equation: (Predictor): Then, the following corrector step introduces an effect of remaining 'current' pressure gradient term as follow: (Corrector): where * α ∆v indicates the incremental velocity from the predicted velocity * α v .
The key point here is the evaluation of 'current' pressure value. By taking the divergence of correction step (Eq. (14)) as: From the incompressible condition (Eq. (9)), which leads: By substitute Eqs. (17) into (16), this leads to the following pressure Poisson equation (PPE): where, which is calculated from SPH approximation.
The above corrector step can be implemented by substituting the pressure gradient with the solution of PPE. The flow chart for solving steps of fluid flow using ISPH method has been introduced in Fig. 1. n=n+1

Input initial model with conditions
Step 4: update position Check END

Start End
Projection scheme: Step 1: predict the velocity

Soil model
In this section, the two different approaches of soil modelling are discussed. In the first approach, the soil model considers the granular material as a fluid with a variable viscosity, where a Bingham type constitutive model is proposed based on Mohr-Coulomb yield-stress criterion [Ulrich, Koliha and Rung (2011)]. The other approach [Bui, Sako and Fukagawa (2007)] is the use of nonlinear material constitutive model in the framework of solid mechanics.

Bingham flow
Here, the soil particles are treated as a viscous material with a variable viscosity. Where, the viscosity is based on Mohr-Coulomb yield-stress criterion for granular material and is derived as follows: For non-Newtonian fluids with the yield strength, the relation between shear stress τ and shear strain rate γ  is given by: where, τ y is the yield shear strength. For shear stresses below the yield stress, a Bingham fluid has behavior such as a rigid body and doesn't deform but when the shear stress surpasses the yield stress, the flow failure occurs resulting in large deformations.
where, so µ is the viscosity after yield, max s µ is the maximum viscosity for a given soil.

Water-soil suspension flow
The suspension layer is nested between soil and water regions depending on the concentration of soil. Particles which reside inside a fictitious suspension layer are assigned to a viscosity that is derived from a Chezy-relation between the stress and the local flow velocity along a route outlined by Fraccarollo et al. [Fraccarollo and Chapart (2002)].

Elastic soil
In the current study, the elastic model will be introduced to describe the soil behavior. Hooke's law is used as the constitutive model: where, is the deviatoric shear strain rate tensor; K is the elastic bulk modulus, which relates to the shear modulus G and Poisson's ratio λ through the following equations:

SPH formulation 3.1 SPH concepts
In the SPH method, a physical scalar function ( ) where, the subscripts i and j indicate positions of labeled particle, and j m means representative mass related to particle j. Note that the triangle bracket i φ 〈 〉 means SPH approximation of a function φ . W is a weight function called by smoothing kernel function in the SPH literature. In the smoothing kernel function, are the distance between neighbor particles and smoothing length, respectively. The divergence of the scalar function can be assumed by using the above defined SPH approximation as follows: Also, the other expression for the gradient can be represented by: In this paper, quintic spline function is utilized as a kernel function as follows: ( ) where, d β is 2 3 7 478 and 3 358 h h π π in two and three dimension space.

Discretization of projection method
Here, the projection method for incompressible fluid and Bingham flow problems is discretized into particle quantities based on the SPH methodology. For this purpose, the gradient of pressure and the divergence of velocity are approximated as follow: The second derivative of velocity for the viscous force and the Laplacian pressure have been proposed by Morris et al. [Morris, Fox and Zhu (1997)] by an approximation expression as follows: where η is a parameter to avoid a zero dominator, and its value is usually given by  (36) and the strain rate tensor of a particle has discretized into the SPH formulations as:

Artificial viscosity
In the numerical solutions, unphysical oscillations are appearing if the dissipative term is not introduced into the governing equations. To improve the numerical stability and to damp out such undesirable oscillations, artificial viscosity ij Π is introduced into the momentum equation as follows: The artificial viscosity ij Π is derived by [Monaghan (1992)] as: In the above equation, α Π and β Π are constants and are chosen according to particular applications; C is the sound speed in soil, which is calculated from elastic bulk modulus and density as

Tensile instability and artificial stress method
In the case of applied SPH in solids, the SPH particles mimic the behavior of the atoms. The instability, which is strictly related to the interpolation technique of the standard SPH method [Rabczuk, Belytschko and Xiao (2004)], is especially noticeable when simulating tension states in solids. The SPH particles forming clumps and causing non-physical fractures in the material. To avoid particles clumping, the artificial stress was proposed by Gray et al. [Gray, Monaghan and Swift (2001)] to eliminate the effects of tensile instability.
The key idea of artificial stress is to introduce a small repulsive force between neighboring particles to avoiding particles clumping. Rabczuk et al. [Rabczuk and Belytschko (2007)] avoided instabilities in simulating the cracks by 3D-meshfree particle method. In the current work, the artificial stress is introduced in three dimensions as follows: where, n is the exponent dependent on the smoothing kernel and ij f is the repulsive force term and is specified, according to Monaghan [Monaghan (2000)], in terms of the kernels as: where d ∆ is the initial particle spacing. h is assumed to be constant in the current work.

Results and discussions
In this section, the numerical examples have been introduced to validate the current scheme. In addition, several numerical tests for water-soil-solid interactions were performed.

Water/soil-suspension flow
The Louvain erosional dam break experiment presented by Fraccarollo et al. [Fraccarollo and Chapart (2002)] is used to validate the suspension model. Fig. 3 introduces the initial diagram of the Louvain erosional dam break experiment. The model has the following dimensions; width of the water is Lw=1 m and height Hw=0.1 m. The width of the soil bed is Ls=2 m and height Hs=0.6 m. The ratio of the density between soil and water is 1.54. The collapse of the water column induces a surge leading to erosions of the soil. A suspension layer form between the pure soil and fluid whose evolution has been tracked in the experiments. Fig. 4 shows the evolution of the three phases in the experiment and the simulations at times 0.25, 0.5, 0.75 and 1.0 sec, respectively. Two simulations are considered in the comparison. The first simulation refers to suspension layers nested between water and soil. The second simulation is performed without a special treatment of the suspension layer. From the current investigation, the simulation with suspension layer gives more accurate validations compare to the experiment. While, the simulation without suspension layer has delay on the wave front of the water over soil compare to the experimental results and suspension case.  [Fraccarollo and Chapart (2002)] and simulations using current ISPH method at times 0.25 and 1.0 sec

Water-soil-solid interactions
Here, water, soil and solid interactions are introduced. To simulate soil hump, Bingham flow for seawall is not enough since the soil is modeled as a fluid with a variable viscosity, which leads to very soft soil even if the physical parameters are chosen with high values. Solid approach for soil is applied to simulate soil hump. The initial models for coupling between water-soil interactions with fully structure body and structure body over soil hump were introduced in Fig. 5. Here, the density of the structure body is taken as 3 2.8 ρ = gm/cm 3 . In Fig. 6, the suspension layers are formed between the water and soil. After impact to the structure body, the splash waves of the water are formed over the structure. Later, the fluid flow makes an erosion behind the structure. Similar tendencies occur between the fluid-soil and structure over soil hump in Fig. 7. The only difference here appears after erosion in the soil hump. This simulation can predict the effects of the soil erosion in a dam break analysis. Additional numerical tests were presented in Figs. 8-10. The effects of the water flow velocity on both of the shape of the waves over structure and formed erosion around structure are presented in Fig. 8. It is observed that the shape of the wave over structure is strongly affected by the initial velocity of the fluid flow. Moreover, an extra erosion occurs near to the structure at low initial velocity of the fluid flow. One simulation trying to prevent erosion near to the structure is showed in Fig. 9. In this simulation, an extended structure behind the main structure was added to prevent erosion. Another simulation which tries to prevent the impact of the erosion by adding a wedge for the structure is shown in Fig. 10. In this case, the wedge can prevent the effects of the erosion around the structure during the fluid impacts.

Conclusion
In this study, an ISPH method has been used to simulate water, soil and solid structure interactions. The fluid is modeled as a Newtonian fluid and the soil is modeled in two different cases depending on the nature of the simulation. Firstly, the soil is simulated by Bingham model. In this model, the granular material is taken as a fluid with derived viscosity from the cohesion and friction angle. Moreover, the Bingham type constitutive model is proposed based on Mohr-Coulomb yield-stress criterion. From the validation with the experimental results, the simulation with suspension layer gives more accurate validations compare to the experiment.
Secondly, the soil is modeled by the solid mechanics, and the soil constitutive model is based on the Hooke's law of linear elasticity. The nested suspension layers between water and soil are formed in both of the two cases. From the current simulations, the following points are reported: • The shape of the wave over structure is strongly depend on the initial velocity of the fluid flow. • An extra erosion occurs near to the structure at a low initial velocity of the fluid flow.
• An extended structure behind the main structure try to prevent erosion.
• Adding a wedge can prevent the structure from the impact of the erosion. Finally, the solid structure is taken as a rigid body in this study and as a future work, the solid structure will be taken as a concrete material with deformations.