Extended Applications of the $\mathbf{\delta }$-SPH Model for the Numerical Study of Fluid–Soil Interactions

1. Introduction

The fluid–soil interaction problem occurs frequently in many natural phenomena and has wide engineering applications. Particularly, under the ocean and coastal engineering circumstances, it is inevitable that architectures and structures like seawalls, wind turbines, and pile lines are exposed frequently to complex environments like earthquakes, storm waves, and submarine landslides. Not only the significant disturbance between the fluid and soil aroused by extreme weather or natural disasters, but also the regular fluid–soil interactions like scour are serious problems during the layout and maintenance procedure of the ocean and coastal engineering. Therefore, it’s necessary to monitor the environmental load condition to ensure workplace safety and resource development efficiency. Noteworthily, the relevant loading analysis can be abstracted to a coupling question including fluid and soil. Thus, aiming to improve the loading forecast ability, the investigation of the interaction between fluid and soil phases is attached to great importance.

To better understand the mechanism of the fluid–soil interactions, scholars have made efforts to extend the studies by launching a series of experiments under different conditions (Spinewine and Zech [2002; 2007]). Although the movement process can be imaged with fast digital cameras and the velocity fields are available by tracking the movements, the loading situations including the stress of both phases are often difficult to capture. Furthermore, besides lacking the comprehensiveness of information, the expensive experimental cost is also worth mentioning.

Fortunately, owing to the rapid development of computer science and technology, numerical simulations have been significantly improved, which serves as a good alternative to investigate fluid–soil interactions. Compared to the experimental method, the interface characteristics and the stress acting of both phases can be extracted instantly through numerical approaches. Furthermore, numerical modeling offers an effective option which enables to permit different scenarios to be investigated quickly at a lower cost.

It is generally recognized that, under the continuum hypothesis, numerical methods are divided into mesh-based and meshless methods. Typical mesh-based numerical methods, like Finite Volume Method (FVM) and Finite Element Method (FEM), have been widely employed to solve hydrodynamic problems or geotechnical problems. However, the mesh-based methods have some limitations in large deformations or large relative motions on the interface, during the formation of scour and erosion induced by the fluid–soil interactions. By contrast, the problems above are not the case for meshless methods that have more advantages in avoiding mesh distortion under the conditions of large deformations or boundary movements. As a meshless particle method, the Smooth Particle Hydrodynamics (SPH) method possesses a harmonious combination of particle properties, Lagrangian properties, and adaptive properties that makes it suitable for solving problems involving high nonlinearity, free surface, interface variation and large deformation. In this case, the SPH method is competent to handle complex motions caused by fluid–soil interactions.

To solve astrophysical problems, Lucy [1977] and Gingold and Monaghan [1977] primarily proposed the SPH method, which was generalized to solve a wide range of engineering problems covering ocean, hydraulics, geo-environmental and geo-technical. Hence, in fluid mechanics and geotechnical fields, an extensive study of the theory and applications of SPH has been done in the last decades. In the fluid field, Monaghan [1994] first introduced the SPH method to the simulation of hydrodynamic problems, introduced the assumption of weak compressibility of fluid when dealing with the motion of free liquid surfaces, and formed the representative Weakly Compressible SPH (WCSPH) method. In terms of multi-phase flow, Morris [2000] proposed a surface tension model, Colagrossi and Landrini [2003] developed a multiphase SPH model for numerical simulation of air entrapment in violent fluid–structure interactions. Hu and Adams [2006] proposed a multiphase flow formula based on volume approximation. The $\delta$-SPH algorithm was then proposed to decrease the pressure noise, being applied to the free surface flow with high accuracy (Marrone et al. [2011]). Some research was dedicated to the improvement of the algorithm later, the Particle Shifting Technique (PST) was introduced to the $\delta$-SPH context to improve the particle distribution by Sun et al. [2017], and the density diffusion term was improved by Fourtakas et al. [2020], which reduced the computational costs. As for geotechnical engineering, Maeda et al. [2006] first applied SPH to the simulation of the infiltration failure of soil. Later, Bui et al. [2007] tried to use SPH to simulate the problem of jet breaking and modified the boundary defects of SPH particles located on the boundary between soil and water. The ideal elastoplastic model of soil was constructed under the SPH framework by Bui et al. [2008], which successfully realized the simulation of soil collapse. Extending from this model, the simulations of some seepage flow issues including embankment failures have been published in the next few years (Bui and Fukagawa [2013]; Bui and Nguyen [2017]). When it turns to the fluid–soil interaction under the coastal and ocean engineering circumstances, the investigations on fluid–soil interactions that have been carried out can be summarized into three categories (Luo et al. [2021]), that is subaqueous debris flow, sediment laden flow and the scouring or erosion of sediment. For the subaqueous debris flow, Shi et al. [2016] adopted the elastoplastic model to simulate the landslide-generated waves under slow and fast conditions, and Ghaïtanellis et al. [2021] considered the soil as granular rheological material and simulated the lake landslide and the subsequent tsunami. In terms of the sediment laden flow, Shi et al. [2017] simulated the sand cloud settling in free surface flows, where the fluid and soil were treated as two miscible fluids in a mixture model. When considering the sediment scouring or erosion problems, the soil is always treated as a slightly compressible pseudo-Newtonian fluid. The CFD-based models (Bui and Nguyen [2021]) like the elastic-viscoplastic model (Ulrich et al. [2013]; Ghaïtanellis et al. [2018]), Herschel–Bulkley–Papanastasiou (HBP) model (Fourtakas and Rogers [2016]; Zubeldia et al. [2018]) and $\mu (I)$-rheological constitutive model (Shi et al. [2019]) are used to describe the behavior of soil. The general elastoplastic constitutive model is also applied in the submerged soil excavation simulation by Yuan et al. [2019].

In general, based on the theory and application improvement, the numerical model of fluid–soil interaction generated a considerable scale of development. However, the current models lack sufficient loading information. Since an accurate prediction of stress is very important for some practical applications. Hence, in order to capture the fluid–soil coupling motion interaction better and offer adequate loading information, in this study, a fluid–soil interactions model will be constructed. In this model, the fluid and soil are treated as the immiscible continua, an elastic–perfectly plastic model is adopted to describe the behavior of soil, and the $\delta$-SPH method will be extended from the fluid phase to the soil phase to improve the accuracy of the stress field. Then the model will be validated on three benchmarks from simple to complex circumstances. After validation, the present model will be applied to simulate the fluid–soil coupling problems under high-velocity conditions.

The remainder of this paper is arranged as follows. In Sec. 2, a brief introduction to the fluid–soil interaction model is given. Particularly, the $\delta$-SPH method, stress diffusive term, two-phase interaction term and the boundary conditions are described. Afterward, in Sec. 3, numerical results reproducing the benchmark cases are presented, which ensures the validity and accuracy of the present model. In Sec. 4, the fluid–soil interaction model will be applied to simulate the case of water jet scouring, followed by the conclusions in Sec. 5.

2. Numerical Model

2.1. Fluid phase

2.1.1. Governing equations

The governing equations of fluid are described by the Navier–Stokes equations. In this paper, the fluid control equations are based on the weak compressibility hypothesis (Colagrossi and Landrini [2003]). Greek superscripts $\alpha$ and $\beta$ denote coordinate directions by employing Einstein’s summation. Therefore, in Lagrangian form, the continuity and momentum equations can be written as

where $\rho$ is the density, v denotes the velocity vector, $\sigma$ means total stress and f refers to the gravitational force. By substituting the SPH equation derived from the particle approximation theory into the governing equation, the governing equation in the SPH formalism can be obtained : where the subscript i denotes the interpolating particle and j refers to the neighboring particles and ${\delta }^{\alpha \beta }$ is the Kronecker’s delta. Aiming to improve numerical stability and suppress undesirable oscillation, the artificial viscosity term ${\mathrm{\Pi }}_{ij}$ is introduced to the momentum equation, which is defined by Monaghan [1994] with the expression of where ${\alpha }_{\mathrm{\Pi }}$ and ${\beta }_{\mathrm{\Pi }}$ are constants with the values of 0.1 and 1.0 in this study, respectively. The definition of the other parameters is given as follows:

$${\mathrm{\Phi }}_{ij}=\frac{h_{ij}{v}_{ij}^{\alpha }\cdot x_{ij}^{\alpha }}{{(x_{ij}^{\alpha })}^2+0.01h_{ij}^2},$$(4)

$${\rho }_{ij}=\frac{{\rho }_i+{\rho }_j}{2},c_{ij}=\frac{c_i+c_j}{2},h_{ij}=\frac{h_i+h_j}{2},$$(5)

$$x_{ij}^{\alpha }=\frac{x_i^{\alpha }-x_j^{\alpha }}{2},v_{ij}^{\alpha }=\frac{v_i^{\alpha }-v_j^{\alpha }}{2}.$$(6)

In the above equations, x is the position of the particle, c is the artificial sound speed derived according to the specific phase and h is the smoothing length set as 1.2 times the particle spacing under the two-dimensional (2D) conditions in this work.

As for the momentum equation of Eq. (2), the total stress tensor $\sigma$ can be divided into the isotropic pressure p and the shear stress $\tau$, which can be written in the general form as

In this paper, based on weakly compressible assumption, the isotropic pressure p is linked to density by means of the simplified Tait Equation of State (Antuono et al. [2010; 2015]) :

where $c_{f0}$ is the artificial sound speed of fluid and its value should satisfy the following restriction : where $p_{max}$ and $u_{max}$ are the maximum predicted pressure and velocity within the problem domain (Antuono et al. [2015]), respectively.

The shear stress in the momentum equation, in terms of Newtonian fluids, is considered to be proportional to the strain rate, which can be calculated as follows :

where $\mu$ is the viscosity, ${\dot{e}}^{\alpha \beta }$ is the deviatoric strain rate tensor written as

2.1.2. Density diffusive term

Due to the fact that the continuity equation in WC doesn’t guarantee strict mass conservation (Liu and Liu [2003]), there can be certain pressure noise phenomena. To address such an issue, the $\delta$-SPH method was proposed by Molteni and Colagrossi [2009]. Hence, the continuity equation in Eq. (2) can be written in the form of

where the coefficient $\delta$ is set as 0.1 for the density diffusive term (Marrone et al. [2011]), which converges to zero with the increase of the particle resolution. Besides, ${\psi }_{ij}$ is defined as where the term 0.01$h^2$ is used to prevent the numerical divergence when $x_{ij}$ approaches zero.

In addition, the $\delta$-SPH method affects the two phases separately. On this account, only the specified phase domain is taken into consideration for the computation of the density diffusive term.

2.2. Soil phase

2.2.1. Constitutive model

The governing equations of soil phase are also used in the form of Eq. (2) with a similar expression of stress tensor as Eq. (7). Instead of the CFD-based approaches, this study chooses the elastic–perfectly plastic constitutive model, which has been applied in simulating large deformation problems (Bui et al. [2008]; Feng et al. [2021]). In the absence of the equation of state, the elastoplastic model combines the yield surface and the plastic potential function to build the connection between the strain rate and stress rate that is used to evolve stress tensors over time. In this work, the Drucker–Prager yield criterion is adopted to determine the state of soil, which is expressed as

where $I_1$ and $J_2$ are the first and second invariants of the stress tensor, respectively, ${\alpha }_{\varphi }$ and $k_c$ are the Drucker–Prager constants associated with the Coulomb material constants c (cohesion) and $\varphi$ (internal friction). Based on the plane strain condition, the Drucker–Prager constants can be derived as

$${\alpha }_{\varphi }=\frac{tan\varphi }{\sqrt{9+12{tan}^2\varphi }},$$(15)

$$k_c=\frac{3c}{\sqrt{9+12{tan}^2\varphi }}.$$(16)

As for the plastic potential function g, the nonassociated flow rule is taken into consideration and the expression is obtained as

where $\psi$ is the dilatancy angle, which is assumed to be zero in this study.

When considering large deformations of materials, the stress rate concerning the rotation of the rigid body should be used. In this case, introducing the Jaumann stress rate, the stress rate can be rewritten as

Considering the above conditions, the constitutive relation in the SPH formalism can be finally obtained :

where G and K are the shear modulus and bulk modulus, respectively, ${\dot{e}}_i^{\alpha \beta }$ is defined in Eq. (11), ${\dot{\epsilon }}_i^{\gamma \gamma }$ is the sum of the normal strain rate components, $s_i^{\alpha \beta }$ is the deviatoric stress and ${\dot{\lambda }}_i$ is the rate of plastic multiplier given by

The detailed derivation of Eqs. (19) and (20) can be referred to Bui et al. [2008]. The definition of strain rate tensor ${\dot{\epsilon }}^{\alpha \beta }$ and the spin rate tensor ${\dot{\omega }}^{\alpha \beta }$ in Eqs. (19) and (20) is given as

In addition, the stress return mapping algorithm involving tension cracking treatment and stress scaling procedure is adopted to correct the stress back on the yield surface, which avoids the tension instability problem. Since this study uses the Leap-Frog scheme for integration, the above algorithm is applied at half a time step and more details of the algorithm can be found in Bui et al. [2008].

2.2.2. Stress diffusive term

When analyzing the stress distribution in geotechnical problems, numerical noise with high frequency can be observed in the stress field under large deformations (Feng [2021]) which is similar to the pressure oscillations observed in the WCSPH simulation of fluid (Molteni and Colagrossi [2009]). As observed in the experiments carried out herein, such a phenomenon affects the accuracy of stress distribution and the model’s ability to predict motion evolution.

To address this issue, a stress diffusive term is formulated by Feng et al. [2021], which significantly reduces the numerical noises and supplies a smooth stress field. After introducing the diffusive term, the stress rate equation should be written as

where stress diffusive term $D_i^{\alpha \beta }$ is similar to the density diffusive term developed by Fourtakas et al. [2020] with the expression as where $\zeta$ is the diffusion coefficient with a value of 0.1. And $c_{s0}$ denotes the artificial sound speed of soil, which is calculated by Young’s modulus and density as follows :

The diffusive operator ${\mathrm{\Gamma }}_{ij}^{\alpha \beta }$ in Eq. (24) is related to the stress components (Fourtakas et al. [2019]) and its definition is given as

where $K_0$ is the coefficient used to calculate the transverse earth pressure at rest and is defined as follows : with the friction angle $\varphi$ that is given by the material properties of the soil.

2.3. Two-phase interaction

To simulate the erosion process during the fluid–soil interactions, the force connected to the relative velocity between two phases was considered in the previous studies (Bui et al. [2007]; Fourtakas and Rogers [2016]). This force is derived from Darcy’s law, which is called seepage force with an expression as

where ${\gamma }_f$ is the unit weight of fluid, n is the porosity, $v_f$ and $v_s$ refer to the velocity of fluid and soil, respectively, and k is the coefficient of soil permeability. After imposing the seepage force, in the SPH formalism, the momentum equation can be rewritten as

In this study, it should be emphasized that we follow the assumption applied in the scouring simulation conducted by Fourtakas and Rogers [2016], that is, the seepage force is only applied to the particle pairs of different phases at the interface.

2.4. Boundary conditions

In this work, the boundary condition is modeled by the “fixed ghost particles”, which consist of several layers of ghost particles. In this case, the variables of the ghost particles can be obtained from the inner particles. In the previous research, the moving-least-square (MLS) interpolation developed by Marrone et al. [2011] has a wide range of applications, but it presents some difficulties for corners. In this work, we have adopted the Shepard interpolation to assign values to boundary particles proposed by Bui and Nguyen [2021]. The properties of the boundary particles can be extrapolated as

where ${\widehat{W}}_{ij}$ is the corrective kernel operator to address the truncate issue of the support domain when the particle approaches the solid boundary, which can be expressed in the form of Chen et al. [1999] :

The above equations create the nonslip boundary conditions, which are adopted under the soil–boundary interaction in this work. When it turns to the water–boundary interaction, we impose the free-slip condition by considering the viscous stress between fluid and ghost particles as zero (Adami et al. [2012]).

2.5. Time integration

As the boundary conditions and initial state are given, the standard time integration method is used to integrate the SPH form equations to obtain the field values of the discrete particles. Due to the low memory requirements and high efficiency, the Leap-Frog (LF) algorithm is employed to solve the equations above. Field variables such as velocity, density, stress tensor and particle position are updated by the following equations :

The stability of the LF integration scheme is restricted by the Courant–Friedrichs–Levy (CFL) condition, which requires the time step to be proportional to the minimum smooth length (Liu and Liu [2003]). In this study, the following criteria were used to determine the size of the time step :

where $C_{\mathrm{\text{cour}}}$ is the Courant coefficient, which is set as 0.25 in this study, $c_{s0}$ and $c_{f0}$ are the sound speed of the soil and fluid phases, respectively, which have been defined above.

3. Benchmark Validations

3.1. 2D dam break

First, we choose a traditional case to simulate the dam break flow against the vertical wall. As shown in Fig. 1, the water area is in the shape of a rectangle with a length of $1.2$m and a height of 0.6m. The probe is placed at the height of 584mm on the vertical wall to record the pressure, and the diameter of the probe is set as 90mm. The particle distance is $\mathrm{\Delta }x=0.015$m with a total of 3200 real particles to form the water area. The initial density of the water particles is ${\rho }_f=1000$kg/m³ with the initial pressure condition derived from the hydrostatic state.

The pressure comparison with the experiment data (Marrone et al. [2011]) is shown in Fig. 2, where we can see that the numerical simulation result is consistent with the experimental result. Despite the peak value of the numerical result being slightly higher than the experiment data, the other regions keep in line with the experimental curve. Hence, it is indicated that the present model is competent to capture the brief information of the fluid movement.

3.2. 2D soil collapse

To further validate the present model, a simulation of soil collapse under noncohesive conditions is carried out. The arrangement of the simulation follows the experimental set-up, which is illustrated in Fig. 3 and the soil particles are labeled to track their position development. In this case, 5000 real particles are used to form the initial soil area in the shape of a rectangle with a length of 0.2m and a height of 0.1m. The particle distance is $\mathrm{\Delta }x=0.0025$m and the material properties follow the experiment of Bui et al. [2008]. The initial density of the soil phase particles is ${\rho }_s=2683$kg/m³, with Young’s modulus is $E=0.84$MPa, and the Poisson ratio is $\nu =0.3$. The friction angles are chosen as $19.8^{\circ }$ and the initial pressure condition is calculated by the static state of the soil.

Figure 4(b) shows the position of soil particles after approaching the stable state, the collapse comes to the end near 0.6s, after that the slope toe remains stable while the upstream has a slight disturbance. Compared with the experiment result in Fig. 4(a), the numerical result achieves a similar distribution, while the toe of the slope expends not so far as the experiment. The distribution of accumulated plastic strain is given in Fig. 4(c), where we can see a clear division exists between the plastic deformation zone and the nonplastic deformation zone. As for the failure distribution result, it provides a failure distribution similar to the previous simulation.

Surface and failure line can be obtained from the above distribution results and their comparisons between the experiments are given in Fig. 5. As can be seen, the two surface curves develop in a similar trend and match well with each other especially from 0. 2m to 0.35m. In terms of the failure line, the two curves also achieve a good agreement though the end of the numerical failure line is located slightly further than the experimental one.

Figure 6 shows the vertical stress contour at the representative times, where we can see the stress distribution is even overall with a stable evaluation. The tensile instability problem was not observed on the free surface of the slope. When the flow develops to 0.1s, a slight degree of noise occurs at the toe of the flow, where several particles present a high-stress state. Since the toe of flow is located outside the failure line, the particles within that region have large relative motion with the boundary under nonslip condition, which leads to the local high stress on the boundary. This phenomenon can be improved by enhancing the resolution, and in this study, it gradually disappears at the following observation time instants with the relevant area becoming stable. The above simulation results show a reasonable behavior of soil, which proves that the current model can offer accurate information under the large deformation of soil.

3.3. 2D Erodible dam break

The simulation above presented a good reproduction of the fluid and soil motion, respectively. In order to evaluate the feasibility of the model in the fluid–soil interactions, we set a case of dam break wave propagating over a movable bed. This study is based on the small-scale experiments carried out by Spinewine and Zech [2007].

In this case, we adopt the same initial configurations as the flat sediment bed case in the experiment (Spinewine and Zech [2007]). The arrangement of the case is shown in Fig. 7, where the initial soil area has a length of 6m and a height of 0.1m and the initial fluid area has a length of 3m and a height of 0.35m. As for the physical parameters, we refer to the settings of Ghaïtanellis et al. [2018]. The initial density of the fluid particles is ${\rho }_f=1000$kg/m³ and the viscosity coefficient is $\mu =10^{-6}$Pa⋅s. The initial density of the soil particles is ${\rho }_s=2683$kg/m³, Young’s modulus is $E=0.8$MPa, the Poisson’s ratio is $\nu =0.3$, the friction angle is set to $\varphi =30^{\circ }$, the porosity is set to $n=0.47$ and the permeability coefficient is set to $k=0.05$cm/s.

The particle resolution is defined by the height of water H and the spacing of the particles $\mathrm{\Delta }x$, here we set the resolution test first including the cases of $H/\mathrm{\Delta }x=35$, $H/\mathrm{\Delta }x=70$ and $H/\mathrm{\Delta }x=140$, where the interface distribution at $t=500$ms is selected to compare with the experimental data. As shown in Fig. 8, the flow toe gets a more sufficient development when $H/\mathrm{\Delta }x=140$ with a closer agreement to the experimental data than the other groups of interfaces. Thus, we choose $H/\mathrm{\Delta }x=140$ for further comparison. In the following simulation results, the particle position distribution will be presented to show the movement of water and soil and be compared with the results of the small-scale dam break experiment at representative times.

According to Fig. 9, we can see that the flow development in the experimental can be reproduced well, the dynamics of the fluid and soil achieve a qualitative agreement between the numerical model and experimental observations. One thing to note here is that the state transition of the erodible and immobile sediment is judged by the displacement magnitude, which is represented the eroded state when the value is greater than the initial particle distance. In this case, all the interfaces can be captured directly from the particle distribution.

In order to estimate the two-phase interaction more directly, Fig. 10 shows a detailed comparison between experimental data (dotted line) and numerical results (solid line), which shows the representative interfaces including erodible–immobile sediment interface (red), the fluid–sediment interface (green) and the fluid free surface (blue). The Root Mean Square Error (RMSE) for each pair of curves is also calculated to supply the comparison, and the definition of RMSE is given as

where $Y_i$ is the experimental interface, $f(x_i)$ refers to the SPH interface location, and N means the total number of observation points.

As shown in Fig. 10, the simulation results match well with the experiment data. The free surface and the fluid–sediment interface present a consistent evaluation of the experiment at the early stages while seeming to lag behind the experimental result as the flow spreads further. Compared to the other two interfaces, the erodible–immobile sediment interfaces achieve the best agreement with the experiment result with the lowest RMSE at the representative times. Considering the interfaces matching results above, despite the acceptable error at the toe of flow, the current model is satisfactory in capturing the movement of fluid–soil interaction.

In this study, we also recorded the stress state in this case, the distribution is given in Fig. 11. From the contour, we can see that the majority of stress distribution is generally consistent with the actual situation. However, at the interface of fluid and soil, a different situation appeared at $t=250$ms: in the region without erosion, the stress distribution is relatively smooth; in the erosion development region, the stress fluctuation occurs at the beginning of the erosion, spreading along the improvement of the flow. When the later observation times are taken into consideration, the majority of stress noise is still located downstream, where the fluid–soil interaction behaves more violently as the flow develops. According to the above results, the current model is also capable of obtaining a reasonable stress distribution, especially under the low relative velocity condition.

4. Applications with Fluid–Soil Interactions

4.1. Vertical water jet scouring

As an important application of coastal and ocean engineering, water jet scouring plays a significant role in several circumstances including the removal of estuarial deposits and the laying and maintenance of the undersea pipeline. Because of the high speed of the water jet, the fluid–soil interactions behave violently, which induces the extremely large deformation of soil. To evaluate the present model for the fluid–soil interactions, a jet erosion case will be conducted.

The initial condition setting of the vertical water jet scouring case is shown in Fig. 12, where the soil area is in the shape of a rectangle, the length of which is 0.5m, and the height is 0.2m. The water particles move vertically downward with a velocity of $v_{\mathrm{\text{jet}}}=10$m/s and the width of the jet is 2cm. The particle spacing is $\mathrm{\Delta }x=0.005$m and the material parameters setting of both phases follow the erodible dam break case above.

First, the jet erosion process is given and the representative times from $t=0.005$s to $t=0.02$s are selected. From the results of Fig. 13, we can see the scour occurs soon after jet flow crushes the soil, and then the soil particles are splashed out of the pit by the overflow of jet flow with the tensile appearance along the outer edge of the pit. The development of erosion matches reasonably with the results obtained by the simulation of Bui et al. [2007], though the morphology of the scour pit is a bit different, which may be caused by the difference in the material properties.

Besides the particle position, the vertical stress distribution is given in Fig. 14. Through observing the stress development, the shock wave is generated from the initial jet impact point and expands deep into the soil phase in a short time. The stress shock wave reaches the bottom at a moment slightly earlier than 0.1s, the following interaction with the solid boundary makes it rebound and go upward, which enhances the soil stress dramatically. Under the circumstances of sufficiently deep soil layers, it can be inferred that the shock wave will continue to spread further without blocking. When focusing on the source of the shock wave, a minor disturbance exists on the impact position and extends over the interface between fluid and soil. Considering the location of the two phases, we can infer that the stress noise is primarily distributed in the fluid phase. The existence of negative pressure can be attributed to the phenomenon of uneven particle distribution at the interface, some voids within the water flow affect the stress field. Since this study focuses on the $\delta$-SPH scheme, it is expected that such a phenomenon will be improved with the implementation of the ${\delta }^{+}$-SPH scheme (Sun et al. [2017]) in the future. This is also an issue we are going to explore in our upcoming research. On the contrary, the stress field of soil phase behaves in a more uniform state.

4.2. Oblique water jet scouring

The initial condition of the oblique water jet scour case is shown in Fig. 15, where all the parameters follow the vertical water jet scouring case, but the particle distance is $\mathrm{\Delta }x=0.0025$m and the incident angle of the water jet has a deflection angle of $45^{\circ }$.

Figure 16 displays the erosion evolution of soil induced by the impact of the water jet at the representative times, which presents some phenomena, unlike the vertical case. Along the incoming direction of the water jet, the soil particles tend to be dragged upward more violently on the impactive side of the erosion pit. In this region, the jet flow runs over the pit and drags up the soil particles like the vertical case, which also amplifies the tensile appearance of the soil. On the contrary, thanks to avoiding the brunt of the water jet, the soil particles located on the weak side of the water jet perform a smaller deformation. Since the majority of fluid flows to the impactive side, no overflow occurs on the weak side, only the initial jet impact on the edge of the pit causes erosion with soil particles being splashed into the jet flow.

Compared with the vertical water jet scouring case, the present results seem reasonable for both soil and water behavior as the scour pit develops in an asymmetrical shape according to the inflow of the water jet. Similar to the vertical water jet scouring case, the distribution of particles at the interface or the free surface is not perfectly uniform. Some small voids still exist though increasing the resolution can improve this phenomenon. Thus, in subsequent studies, we will also make efforts to improve the particle distribution by adopting the approach mentioned in the vertical water jet scouring case.

5. Conclusions

Aiming to improve the capacity of the numerical simulation of fluid–soil interactions, a fluid–soil interaction model based on the SPH framework is applied. The water phase is treated as a Newtonian fluid with weak compressibility while the elastic–perfectly plastic model with Drucker–Prager yield criterion and nonassociated flow rule is used to describe the behavior of soil. It is worth noting that the application of $\delta$-SPH algorithm has extended from water to the soil phase, which apparently makes a contribution to improving the pressure distribution with the conjugation of stress diffusive term. In addition, to reproduce the erosion during the fluid–soil interactions, the seepage force is introduced to establish the connection between the two phases.

As for the validation, the simulation tests from the single-phase flow to the two-phase coupling have been executed, the simulation results achieve a reasonable agreement with the experiment data including the profiles of the representative interfaces or free surfaces.

After passing the basic validation, the present model is applied to simulate the vertical and oblique water jet scour case under a violent water impact and successfully perform the extremely large deformation of soil involving the splash phenomenon. The particle method shows strong robustness without any extra numerical techniques, which demonstrates the advancement of SPH in simulating fluid–soil interactions. Based on the performance, we can infer that SPH has a promising future with more extensive application, which covers the local scour of pile foundations for wind power or jacket platforms, erosion of coasts and a wide range of coastal and ocean engineering applications.

Acknowledgments

This research is supported by the National Key Research and Development Program of China (Grant No. 2022YFC2806300), the Fundamental Research Funds for the Central Universities, Sun Yat-sen University (Grant No. 231gbi023), and the GHfund A (Grant No. 202302014084).