Coupled CFD-DEM simulation of submarine landslide movement behavior

Submarine landslide is a common disaster geological phenomenon in the ocean, which can be destructive to underwater infrastructure. Therefore, it is necessary to simulate the evolutionary behavior of submarine viscous landslides. In this paper, computational fluid dynamics (CFD) and discrete element method (DEM) are used to establish the fluid-solid coupling model of water-particle interaction. Firstly, the inter-particle cohesion model is introduced to train the coupled CFD-DEM analysis of the kinematic evolution process, and at the same time, two typical cases are carried out to verify the analysis. In the experimental simulation, the kinematic and morphological characteristics of the submarine landslide were simulated considering the viscous effect and initial velocity of the landslide, and the influence mechanism of the landslide motion and evolution process was investigated in depth. The results show that the coupled method can better simulate the motion of submarine landslide, and the viscous effect of the landslide has a significant influence on its kinematic and morphological characteristics, and the initial velocity also significantly affects the evolution and distribution characteristics of the particle flow field of each part of the landslide in the process of motion. This study is important for the simulation and effective prediction of the kinematic evolution process of real submarine landslides.


Introduction
As the development of offshore oil and gas resources enters deep water, offshore oil and gas extraction wells, submarine pipeline and cable systems (including oil and gas pipelines, submarine fiber optic cables and cables) will face more and more complex marine geohazard problems [1].Among them, submarine landslides are one of the more common marine geological hazards, often occurring in offshore deltas and continental slopes, which are usually three to four orders of magnitude larger than onshore landslides and are very destructive, causing serious damage to underwater infrastructure and even triggering tsunami disasters once they occur [2].
At present, submarine landslides have attracted the attention of experts and scholars from various countries, especially the simulation and prediction of submarine landslide motion process has become a current research hotspot [1,3].In the analytical research, the mainstream numerical simulation methods include: finite difference method (FDM) [4,5], large deformation finite element (LDFE) [6], computational fluid dynamics (CFD) method [7], and the computational fluid dynamics (CDFE) method.computational fluid dynamics (CFD) [8], and material point method (MPM) [9].Using the above numerical simulation methods, scholars have carried out many valuable explorations on the motion process of submarine landslides, such as: Gauer et al [10] used CFD numerical method based on Bingham model to treat the landslide body as a plastic non-Newtonian fluid, and inverted the motion of the largest known Storegga submarine landslide in the world, and verified it by comparing with the measured landslide sediment profiles.Most of the above numerical methods assume that the landslide body is a single-phase continuous medium, ignoring the coupled water-soil exchange between seawater and the landslide body during the motion of the submarine landslide, which makes it difficult to restore the real evolution of the submarine landslide under the influence of the water environment.Therefore, it is necessary to combine the theoretical methods of continuous and discontinuous media [11] and develop CFD-DEM fluid-solid coupling analysis methods to simulate the motion evolution behavior of submarine viscous landslides.
In this paper, the open-source software Open FOAM, LIGGGHTS, and the coupling module CFDEM are used to describe the particle composition of landslides using the discrete element method (DEM) [12], and the CFD method is used to calculate the motion of ambient water and its interaction with particles [13,14].On this basis, the cohesive force model is introduced into the particle equation of motion and coupled with the open-source program Open FOAM to realize the simulation of the motion process of cohesive landslides and extend to the simulation of the kinematic and morphological characteristics of submarine cohesive landslides.This work is an inheritance and development of the existing coupled CFD-DEM analysis method, and provides a possibility for future in-depth investigation of the motion and evolution process of real submarine landslides and disaster prediction This work is a continuation and development of the existing coupled CFD-DEM analysis method, and also provides the possibility to investigate the motion and evolution of real submarine landslides and disaster prediction.

CFD-DEM equation
When The control equation of fluid motion is the modified Navier-Stokes equation.Its vector form of mass conservation, momentum conservation equations are [13].
Equation:  is the porosity. is the fluid density;  is the fluid velocity. is the pressure. is the fluid viscous stress tensor. is the particle being the interaction force between particles and fluid; g is the acceleration of gravity.
According to Newton's second law of motion, the advection and rotation of particle iControl equations are [15], respectively.
Equation:  is the particle mass;  , is the velocity of particle;  is the total external force acting on the particle [15]; is the rotational inertia; , is the angular velocity of particle;  and  are the tangential force and rolling friction respectively, moments generated by the rolling friction.
The interaction force between the particles and the fluid is denoted by  ,and its expression is as follows.Equation:  D, is the drag force;  LM, is the lift force;  is the number of particles in a calculation cell; ∆ is the volume of a fluid calculation cell.

𝑭
The CFD-DEM solution process is as follows.
(1) Establish a CFD-DEM water-soil coupled computational model for landslides under different operating conditions.
(2) Given the initial flow field, initial particle position, initial particle velocity and fixed step size.
(3) Introduce the cohesive force model in the particle equation of motion, and the DEM solver calculates the particle motion position and velocity.
(4) The particle position and velocity information is passed to the CFD solver.
(5) Discrete solution of the momentum equation based on the velocity and pressure fields, then solving the pressure correction equation to correct the pressure and velocity, and iterating until the flow field converges.
(6) Calculate the fluid volume fraction and the mean particle velocity in the fluid grid.( 7) Derive the fluid forces on each particle and solve for the particle-fluid momentum exchange source term.
(8) Considering the void fraction and the momentum exchange source term, the CFD solver solves for the flow field degree.
(9) After the convergence of the flow field, the solution of the particle phase is carried out, i.e., repeat steps (3) to ( 9) until the end of the simulation.
(10) coupled with the above process to calculate the force characteristics in the process of interaction between soil particles and seawater and obtain the dynamic change process of viscous sliding body movement in seawater.Particle control equation considering cohesion model.

Particle control equation considering cohesion model
When considering the cohesive force, the total external force Fi acting on the particle can be expressed as: C,  D,  LM,  CO, 6 Equation:  C, the inter-particle contact force. CO, is the cohesive force. CO, in Equation ( 6) is analyzed by JKR model.JKR model is one of the most widely used mathematical models in current research, where the cohesion increases with the contact area; Johnson et al. [16] assumes that the cohesion value between two particles  CO, is determined by the intrinsic cohesion energy density k and the effective contact surface area S, i. e.
CO,  7 Equation:  is defined as the cohesion per unit effective surface area, namely the cohesive energy density (J/ );  is the effective contact area.The expression of  is relatively complex, and its value can be based on experience or reference to [17], or the expression [18]  Equation: ∆ is the center distance of particles i and j;  and  are the radii of particle.and.Coupling (7), (8), and ( 9), the new cohesive force expression  CO, can be obtained, and then coupling ( 6) and (3), the controlling equation of particle motion containing the cohesive force model can be derived, thus realizing the CFD-DEM fluid-solid coupling simulation of cohesive slip body motion.It is worth noting that the cohesive force here is different from that of cohesive soils (clay, pulverized clay) in geomechanics, and the correlation or conversion relationship about it still needs to be studied.

Model setup and study conditions
In order to effectively simulate the evolution of the motion of submarine viscous slides.In this paper, the original model of simulation calculation is constructed.As shown in Figure 1, a seabed simulation area with 4 m length, 2 m width and 1 m depth is 5°, and a volume of 50 mm 50mm 50 mm sliding block is located in the fluid at one end of the depth of the simulation area, causing collapse and then underwater landslide.
The calculation conditions consider are detailed in Table 1.Other parameters were set as: particle diameter 1 mm, particle density

Numerical model validation
For the above CFD-DEM coupling equation and its numerical model, a progressive validation work was conducted using 2 cases.Firstly, the settling process of individual particles in water is simulated to verify the accuracy of the coupling force calculation to test the reliability of the CFD-DEM coupling procedure; subsequently, the collapse process of the particle accumulation is simulated, where the accumulation collapses to form a flow and eventually reaches stability.Finally, the numerical simulation of the evolution process of viscous granular flow transport is compared with the indoor model tests carried out by Trygve et al [20] in order to qualitatively verify the reliability of developing a fluid-solid coupling method for the simulation of viscous slip motion processes.

Verification of single particle sedimentation
Set up a fluid domain of 0.1 m x 0.1 m x 0.5 m, As shown in Figure 2. The main parameters used for the numerical simulations are listed in Table 2.The sedimentation process of individual particles in water was simulated to verify the accuracy of the coupling force calculation to demonstrate the reliability of the CFD-DEM coupling procedure.

Verification of collapse of non-cohesive particle accumulation
The collapse process of particle accumulation shows the collapse to form a flow and finally reach a stable.The original model for the numerical simulation is shown in  To verify the accuracy of particle flow-flow field bi-directional coupling calculation, the selected parameters are shown in Table 3.The distribution of particle velocity and flow field velocity was obtained at different moments, as shown in Figure 3.

Figure 3. Simulated particle and flow field velocity under different time scale
In conducting the simulation, the initial state is assumed to be a sliding body on a submerged slope capable of collapsing under its own weight, i.e., after the action of the trigger factor, it slides under the action of gravity, and the evolutionary process of its motion on the seafloor after sliding and the change of its morphological characteristics are examined, without considering the effects brought about by waves and bottom currents, etc. in this process.

Simulation of movement characteristics of viscous landslides
The viscous aggregation energy densities of different slides were set under different working conditions.The morphological characteristics and kinematic characteristics (motion speed and transport distance) of the slider with different viscosities were obtained by simulating the motion process of the slider under water.Different parts of the initial slider are represented by different colors, as shown in Fig. 4, where red represents the leading-edge part of the slider, green represents the middle part of the slider, and blue represents the trailing edge part of the slider.It can be seen from Figure .7 that the lateral width of the landslide body expands and the speed increases as the landslide body moves forward.
With the increase of viscosity, the velocity and transport distance of the landslide body gradually decrease.In addition, it can be seen that the average sliding distance of the landslide body increases linearly with time and decreases with the increase of viscous energy density at a certain time.By comparing and analysing the numerical results for shorter time periods, it is found that the sliding distance of non-viscous slides is quite farther than that of viscous slides, up to about 8%.In particular, it is noted that the kinematic characteristics (e.g., velocity, distance, etc.) of the landslide body at different locations during the movement of the landslide body are significantly different.In order to analyze the kinematic characteristics of the local location of the landslide body more visually and clearly, the movement velocity (which is the average velocity rather than the velocity of a characteristic particle) at each location of the trailing, middle and leading edges of the landslide bodyis real-time monitoring.
The motion velocity of all parts of the sliding body showed the basic characteristics of increasing rapidly with time and decreasing slowly after reaching the peak velocity.The central and trailing edge parts of the slip are the latest to reach the peak velocity, which is due to the viscous traction effect of the front part of the slip and the inertial force of the rear part of the slip (the acceleration path is longer and its gravitational potential energy is released more fully), so that a larger average velocity can be reached.With the friction of the seabed and the interaction between the slip and the ambient water, the velocity of each part of the slip gradually decreases slowly after reaching the peak, and its decay rate shows a gradual increase from the front to the back, and there is a significant normalization trend of the average velocity of each part of the leading, middle and trailing edges, as shown in the overall trend diagram in Figure

Simulation of the morphological characteristics of a viscous slip with initial velocity
For the D-F conditions in Table 1, the simulation of the motion process was carried out by setting up viscous slides with different initial velocities, focusing on the evolution of morphological features, and the results are shown in Figure 7, where different parts of the slides are represented by different colors.Analysis of Fig. 7 shows that the initial velocity has a great influence on the morphological characteristics of the viscous sliding body in motion.As the initial velocity increases, the length of the sliding body becomes smaller, the lateral expansion increases, the morphology of the sliding body changes significantly, and the center of gravity of the shape of the sliding mass becomes more aggregated.In contrast to the slip body without initial velocity, it can be seen that after a period of slip body motion with initial velocity, the particles are more fully mixed, and the particles in the middle and rear of the initial state also appear in the leading edge part of the slip body.The larger the initial velocity is, the more aggregated the shape of the sliding body is.By comparing the results of numerical analysis, it is found that the slip distance with initial velocity is more than 8% farther than the slip distance without initial velocity.8, in the initial stage of movement, with a larger initial velocity of the slide, due to the leading edge of the greater water resistance and friction at the bottom, making the slide in the state of resistance than the sliding force, its speed will gradually decrease; and no initial velocity of the slide in the state of the sliding force is greater than the resistance, it is gradually accelerated.After a period of movement, the velocity changes due to the collision and cohesion of the particles inside the sliding body.In short, the average sliding distance of the sliding body to meet, the greater the initial velocity, the farther the sliding distance of the law.Figure 9 shows the relationship between the dimensionless length λ A and time t for different initial velocities of the sliding body.The analysis shows that for a certain initial velocity, the length λ A increases linearly with the increase of time t, and then gradually decreases and stabilizes after reaching a certain value.For a certain moment, the dimensionless length λ A decreases significantly with the increase of the initial velocity of the slider.In addition, the larger the initial velocity of the slide, the lower and earlier the peak point of λ A appears.And Fig. 10 gives the variation of dimensionless width λ B with time.The analysis shows that the dimensionless width λ B increases with time; in particular, as the initial velocity increases, the value of λ B increases accordingly.As the initial velocity of the slider increases, the length of the slider decreases while the width increases.In the case of the change in the initial velocity of the slide, the greater this initial velocity, the greater the water resistance to the front of the slide, the slower its evolution, and thus the phenomenon of the rear pushing out the front; therefore, when comparing the effect of the size of the initial velocity, there is a special phenomenon that the length of the slide decreases and the lateral expansion increases with the increase of the initial velocity.

Conclusion
For submarine viscous landslides, a coupled CFD-DEM model and analysis method describing the interaction between environmental water and soil particles are established, and the characteristics of submarine landslide kinematic behavior and evolution process are thoroughly discussed.
The viscous cohesion model is introduced into the coupled CFD-DEM algorithm, which better reproduces the motion behavior and evolution process of submarine viscous landslides.
By simulating the kinematic behaviors of different viscous slides, the significant effects of slip viscosity on the kinematic characteristics (distance and velocity) and morphological characteristics (length, width, shape, etc.) of submarine landslides are revealed.

Figure 2 .
Figure 2. Model of particle accumulation collapse

Figure 4 .Figure 5 .
Figure 4.The movement characteristics of sliding mass under different working conditions

Figure 6 .
Figure 6.Velocity and time relationship of each part of sliding body

Figure 7 .Figure 8 .
Figure 7.The morphological characteristics of sliding mass movement at different initial velocities

Figure 9 . 9 Figure 10 .
Figure 9.The relationship between dimensionless length A and the time at different initial velocities [19]in the simplified Johnson-Kendall-Roberts (SJKR) model: Equation:  Is radius;  is Poisson's ratio;  is Young's modulus;  is energy per unit contact area, it can be obtained by combining numerical simulation and test.The expression of the concerned  in the SJKR model in LIGGGHTS is[19]

Table 2 .
Parameters for single particle sedimentation simulation