Unresolved CFD and DEM Coupled Simulations on Scour around a Subsea Pipeline

: In this paper, numerical studies were carried out on scour around a subsea pipeline. A coupled solver between computational ﬂuid dynamics (CFD) and discrete element method (DEM) was selected to simulate ﬂuid ﬂow and particle interactions. To select and validate the numerical model parameters in the solver, angles of repose and incipient motion were simulated. From the validation studies, the selected coefﬁcient of rolling friction with spherical particles could predict the behavior of non-spherical particles. The ﬂuid ﬂow around the subsea pipeline was simulated, and the motion of individual soil particles was tracked. Particle motions were generated by the drag force, due to a high velocity. Three scour development process, such as onset of scour, tunnel erosion, and lee-wake erosion, were studied and discussed. The scour depth evolution showed good agreement with the experimental data. It was conﬁrmed that the selected solver, with numerical model parameters, predicted the scour process around a subsea pipeline well. lower than the experimental data in the onset of scour and tunnel erosion. While the scour depth was in good agreement with the experimental data in the lee-wake erosion, the selection of the coefﬁcient of rolling friction by the validation test in this study more accurately inﬂuenced the prediction in the onset of scour and tunnel erosion.


Introduction
When subsea structures are placed on a sea floor, vortices form due to currents and waves occur around foundations and subsea pipelines. The formation of a vortex structure is determined by the geometry and flow [1]. The vortices are related scour, which results in the loss of the soil and erodes the seabed. The scour process around a subsea pipeline is expressed by the interaction between the fluid flow and the seabed [2]. The fluid flow around the subsea pipeline may cause a free span. Excessive bending moment and fatigue load, due to the free span, can reduce the stability of the pipeline, leading to a structural failure [3]. Pipelines installed on the sea bottom can be classified into two main purposes: submarine cables for power supply and communication systems and subsea pipelines that transport oil or gas to offshore or onshore processing facilities. When the structure is damaged, it causes economic damage, due to the interruption of continuous resource production work. In addition, it takes enormous cost and time to recover the polluted marine environment caused by the oil spill. For this reason, it is important to accurately predict scour around the subsea pipeline, which requires an understanding of the interaction between the fluid flow and seabed.
Many researchers have done a series of experiments to study the behavior of scour around the subsea pipeline [3][4][5][6]. Mao [3] concluded that the difference in pressure between the upstream and downstream of the pipeline caused the erosion of sand particles beneath the pipeline. Chiew [7] used the term piping to indicate a dominant factor in the early stage of scour. For factors influencing the scour process, stream velocity, grain size, pipe diameter, initial gap between a pipeline and a seabed, and flow depth were experimentally defined, with various environmental conditions [3,5,7]. Empirical equations for the depth of scour beneath the subsea pipeline was presented [5,8,9].
Various numerical models have been proposed, such as single-and two-phase models, to simulate a scour process. Single-phase models could be subdivided into a rigid seabed approaches and deformable seabed approaches. In the rigid seabed approaches, scour was predicted by the shear stress magnitude at the seabed boundary [10,11]. On the other hand, in the deformable seabed approaches, the elevation of the seabed was changed by a morphological model, which is considered a sediment transport model [12][13][14]. To consider the deformable seabed, the periodic updates of the deformed mesh were required [15,16]. Zhao et al. [17] investigated a local scour around various pipelines with combined wave and current conditions by a fully coupled numerical model. Liang et al. [18] established a three-dimensional model to investigate scour around underwater structures and found that the equilibrium scour depth increased with the pipeline length until the pipeline length exceeded four times the pipe diameter. Above all, the single-phase models were unable to consider interactions between particles.
The two-phase models could be subdivided into two approaches depending on the fluid and seabed modeling method: Euler-Euler approach and Euler-Lagrangian approach. The Euler-Euler approach was based on an interpenetrating continuum for the fluid and particles. The movement of the particles that made up the seabed was explained by the sediment concentration. Yeganeh-Bakhtiary et al. [19] employed the Euler-Euler two-phase model to simulate the tunnel erosion beneath a pipeline and successfully predicted the bed profile and flow behavior. Fraga et al. [20] investigated scour beneath piggyback pipelines with current flow by a two-phase Euler-Euler model. The Euler-Euler two-phase model had limitations, in that it was difficult to provide information about particles, such as particle and particle interactions, as well as particle position and velocity [21]. On the other hand, the Euler-Lagrangian approach treated the soil as a discrete phase that allowed tracking of individual particles. The Euler-Lagrangian approach has been applied in various engineering problems to simulate fluid and particle flow.
Computational fluid dynamics (CFD) and discrete element method (DEM) coupled models were commonly used in the Euler-Lagrangian approach. DEM had the advantage of being able to predict the exact behavior of particles and handle a large amount of particles [22,23]. CFD and DEM coupled models could be divided into resolved or unresolved approaches, depending on how fluid and particle interactions were considered. The resolved CFD and DEM models were fully resolving the flow around each particle and directly calculating the drag force. A grid that was at least 8-10 times smaller than a particle diameter was required, indicating huge computational requirements [24]. The resolved methods were adequate for a dilute flow system or small system. In the resolved CFD and DEM model, a dense grid system was needed to compute the fluid force acting on a particle; whereas, in the unresolved CFD and DEM coupled model, a dense grid was not required to obtain an accurate drag [23,25,26]. The unresolved CFD and DEM coupled model presented high computational efficiency for flows with bulk particles [23,27]. However, there was the grid dependency that yields valid computational results for fluid and particle interactions [28]. Various numerical methods have been proposed to increase the reliability of the results, regardless of the grid size ratio [29][30][31][32][33]. Song and Park [33] used a kernel function-based averaging method to reduce the grid dependency.
The CFD and DEM coupled models have been applied to geotechnical engineering problems, e.g., the seepage flow in soils, sediment transport, and sand pile formation [34][35][36]. Simulations for the scour process around a subsea pipeline were recently carried out, and the applicability of the CFD and DEM coupled model was discussed [37,38]. Zhang et al. [37] used the CFD and DEM coupled model to investigate the onset of scour, and the motion of each particle was considered in detail. Yang et al. [38] applied the CFD and DEM coupled model to simulate the scour process and provided detailed information to better understand the scour mechanism below a pipeline. Hu et al. [39] adopted the CFD and DEM coupled model to study the effect of gap ratio and incipient velocity on a local scour around two pipelines in tandem and revealed that a scour development was closely related to the individual particle behavior. Liu and Tau [21] employed the CFD-DEM coupled model to simulate a local scour behavior around a bridge pier with clear-water conditions. Yazdanfar et al. [40] proposed a novel upscaling methodology, by a microscale CFD and DEM coupled model, to predict a live-bed scour. However, studies on local scour around the pipeline are still insufficient. More research on improved CFD and DEM coupled models is needed to overcome existing constraints, such as high computational cost, inaccuracies, due to the use of spherical particles, and grid dependency [21,28,37,38]. Recently, Song and Park [33] presented an improved unresolved CFD and DEM coupled solver for particulate flow. The solver showed that it can reduce grid dependency and predict particle sedimentation for the single-particle settlement.
The purpose of this paper is to investigate the applicability of the CFD and DEM coupled solver, developed by Song and Park [33], to simulate fluid flow and seabed interactions around a subsea pipeline. For the simulations, the unresolved CFD and DEM coupled solver using the kernel-based averaging method was used [33]. The novelty and purpose were to suggest a proper selection procedure for the numerical model parameter and apply on scour around the pipeline. From the selected rolling friction coefficient, scour with non-spherical particles could be predicted by spherical particles. The selected solver was developed based on the open source CFD (OpenFOAM) and DEM (LIGGGHTS) libraries. To verify and validate the numerical methods, angles of repose and an incipient motion were simulated.
The present paper is organized as follows. Section 2 describes computational methods including governing equations for the unresolved CFD and DEM coupled solver and numerical methods. Section 3 presents the problem description. Section 4 shows numerical model parameters. Section 5 shows the results and discussion for the angle of repose, incipient motion of particles, and scour process around a pipeline. Finally, in Section 6, some concluding remarks are provided.

Discrete Element Method
In the DEM solver, each particle was individually considered. The collisions between particles and walls, as well as particles and particles, were simplified to a linear springdashpot model [41]. An individual particle was described by Newton's second law. The governing equations for translational and rotational motions of a particle were, respectively, written as follows: where m p , u p , f p,p , and f p,w are the particle mass, particle translational velocity, contact force between particle and particle, and contact force between particle and wall, respectively. m p g and f p, f are the gravitational force and particle and fluid interaction force, including the drag force ( f d ), and buoyance force ( f b ), respectively. I p , ω p , and T t are the inertia moment of a particle, angular velocity of the particle, and torque due to the rolling resistance, as well as the collision with the particle or the wall, respectively. The Hertzian law of contact was used in conjunction with Coulomb's law of friction to describe the particle and particle contacts [23]. The particle and wall interaction were modeled in the same way as the particle and particle interaction, and the wall is treated as a particle with an infinitely large radius. The contact force ( f p,p ) between particle and particle can be divided to the normal ( f c,n ) and tangential ( f c,t ) forces, as follows: 4 of 20 where subscripts n and t indicate the normal and tangential directions, respectively. k and γ are the stiffness and damping coefficients, respectively. µ s is the sliding friction coefficient. δ and v are the displacement and relative velocity between two contacting particles, respectively. Since spherical particles were used in the simulation of non-spherical particles, such as sand and rock, and the effect of the shape of the actual particles should be considered by using a rolling friction model in the DEM solver [42,43]. In the present study, the directional constant torque rolling friction model was used to consider particle non-sphericity [23,42,44].
where M r is the total rotational torque. f c,n is the contact normal force. R r is rolling radius defined by R r = r i r j / r i + r j . r i and r j are the radius of the two spherical particles in contact, and µ r is the coefficient of rolling friction.

Computational Fluid Dynamics
The fluid satisfies the mass and momentum conservation equations, and the motion of the particles was described by the volume fraction transport equation. The momentum equation, related to the coupling with the DEM solver, was expressed as: where α f , u f , p, τ f , and g are the fluid volume fraction, fluid velocity, pressure, and stress tensor of the fluid, and gravity, respectively. R f ,p represents the volumetric momentum exchange between the fluid and particle. α f and R f ,p could be expressed as: where N p , V p,k , V c , and w p,k are the number of the particles, volume of the kth particle, volume of a single mesh, and weighting factor of the kth particle to a mesh, respectively. The weighting factor was expressed by the gaussian kernel function. < u p > and f p, f are the averaged particle velocity, based on each mesh, and forces acting on the particles by the fluid, respectively. In f p, f , the drag force ( f d ) and buoyance force were considered the main interaction forces between particle and fluid. More details for the unresolved CFD and DEM coupled solver are described in Song and Park [33].

Numerical Methods
The time derivative terms were discretized using a second-order accurate implicit scheme. The gradient at the cell centers was calculated using the Euler method. The convection terms were discretized using a second-order accurate total variation diminishing (TVD) scheme with the Sweby limiter [45]. The velocity and pressure coupling and overall solution procedure were based on the pressure implicit with the splitting of the operator (PISO) algorithm. For the turbulence closure, the k-ε turbulence model was used. The discretized algebraic equations for the pressure and velocity were solved using a pointwise Gauss-Seidel iterative algorithm, while an algebraic multi-grid method was employed to accelerate solution convergence.
The time step in the CFD solver was set to be bigger than that in the DEM solver. In every simulation, the time step in the DEM solver was set to be 20% smaller than the Rayleigh critical time step [33]. The simulation time was advanced when the normalized residuals for the solutions had dropped by six orders of magnitude.

Problem Description
The computational domain for the CFD and DEM simulation is shown in Figure 1.
Here, D is the pipeline diameter. The lengths from the center of the pipeline to the inlet and outlet boundaries were 10 D and 12 D, respectively. The water depth (H) and seabed depth (h) were 5 D and 1 D, respectively. The soil particles for the DEM simulation were randomly distributed and then settled under gravity to form the seabed. The pipeline was initially laid just above the seabed without an initial gap and fixed during the scour process [3]. The length of the seabed was set to 20 D, and the lengths of the inlet and outlet directions from the pipeline center were 9 D and 11 D, respectively. The width of the seabed was set to 2.16 mm in the out-of-plane direction. Due to the high computing cost, the length of the seabed was shorter than that of Mao's experiment; however, Yang et al. [38] confirmed that the development of scour could be sufficiently considered.
(TVD) scheme with the Sweby limiter [45]. The velocity and pressure coupling and overall solution procedure were based on the pressure implicit with the splitting of the operator (PISO) algorithm. For the turbulence closure, the k-ε turbulence model was used. The discretized algebraic equations for the pressure and velocity were solved using a pointwise Gauss-Seidel iterative algorithm, while an algebraic multi-grid method was employed to accelerate solution convergence.
The time step in the CFD solver was set to be bigger than that in the DEM solver. In every simulation, the time step in the DEM solver was set to be 20% smaller than the Rayleigh critical time step [33]. The simulation time was advanced when the normalized residuals for the solutions had dropped by six orders of magnitude.

Problem Description
The computational domain for the CFD and DEM simulation is shown in Figure 1.
Here, D is the pipeline diameter. The lengths from the center of the pipeline to the inlet and outlet boundaries were 10 D and 12 D, respectively. The water depth (H) and seabed depth (h) were 5 D and 1 D, respectively. The soil particles for the DEM simulation were randomly distributed and then settled under gravity to form the seabed. The pipeline was initially laid just above the seabed without an initial gap and fixed during the scour process [3]. The length of the seabed was set to 20 D, and the lengths of the inlet and outlet directions from the pipeline center were 9 D and 11 D, respectively. The width of the seabed was set to 2.16 mm in the out-of-plane direction. Due to the high computing cost, the length of the seabed was shorter than that of Mao's experiment; however, Yang et al. [38] confirmed that the development of scour could be sufficiently considered. The simulation conditions are presented in Table 1. The pipeline diameter (D) was 0.05 m. The shape of soil particles in the simulation was a uniform spherical particle with a diameter ( ) of 1 mm. The Shields parameter ( ) was 0.33 [3], which meant that livebed scour occurred when the mean bed shear stress in the upstream was larger than the threshold value required to move the soil particles. The large Shields parameter ( ) was chosen to accelerate scour development. The particle diameter was larger than that of Mao's experiment. To keep the Shields parameter ( = 0.33) to be the same as that in Mao's experiment, the mean flow velocity ( ) was increased to be 1.21 m/s. Based on the diameter of the pipeline, the Reynolds ( ) and Froude ( ) numbers were 6.05 × 10 and 1.71, respectively.   The simulation conditions are presented in Table 1. The pipeline diameter (D) was 0.05 m. The shape of soil particles in the simulation was a uniform spherical particle with a diameter (d p ) of 1 mm. The Shields parameter (θ) was 0.33 [3], which meant that live-bed scour occurred when the mean bed shear stress in the upstream was larger than the threshold value required to move the soil particles. The large Shields parameter (θ) was chosen to accelerate scour development. The particle diameter was larger than that of Mao's experiment. To keep the Shields parameter (θ = 0.33) to be the same as that in Mao's experiment, the mean flow velocity (U) was increased to be 1.21 m/s. Based on the diameter of the pipeline, the Reynolds (Re D ) and Froude (Fr D ) numbers were 6.05 × 10 4 and 1.71, respectively. On the inlet boundary, the Dirichlet condition was applied for the velocity, turbulence properties, and volume fraction, while the Neumann condition was applied for the pressure. On the outlet boundary, the Neumann condition was applied for the velocity, turbulence properties, and volume fraction, while the Dirichlet condition was applied for the pressure. A zero normal gradient was applied for the Neumann conditions of all variables. The no-slip condition was applied to the pipeline and bottom surfaces. At the free surface, the symmetry boundary conditions were applied. The fluid inlet boundary was specified with a logarithmic velocity profile, based on the friction velocity and bed roughness [19]. The mean velocity was calculated from the integral of the flow velocity over the depth at the inlet boundary. The mean flow velocity specified at the inlet was increased linearly from zero to 0.87 m/s over the first 6.4 s, in order to obtain the numerical stability [37,46].
The flow velocity profile was reproduced when there were the particles in the seabed as shown in Figure 1. The flow velocity profile by the selected unresolved CFD and DEM coupled solver was compared with that of a CFD solver. In the coupled solver simulation, the seabed particles were considered, while, in the simulation of the CFD solver, the seabed particles were not considered. Figure 2 shows the flow velocity profiles simulated by both solvers at a central location in a domain without the pipeline. The height of 0 m means the top of the sand surface. It was confirmed that the desired flow velocity profile could be obtained, even with the particles. On the inlet boundary, the Dirichlet condition was applied for the velocity, turbulence properties, and volume fraction, while the Neumann condition was applied for the pressure. On the outlet boundary, the Neumann condition was applied for the velocity, turbulence properties, and volume fraction, while the Dirichlet condition was applied for the pressure. A zero normal gradient was applied for the Neumann conditions of all variables. The no-slip condition was applied to the pipeline and bottom surfaces. At the free surface, the symmetry boundary conditions were applied. The fluid inlet boundary was specified with a logarithmic velocity profile, based on the friction velocity and bed roughness [19]. The mean velocity was calculated from the integral of the flow velocity over the depth at the inlet boundary. The mean flow velocity specified at the inlet was increased linearly from zero to 0.87 m/s over the first 6.4 s, in order to obtain the numerical stability [37,46].
The flow velocity profile was reproduced when there were the particles in the seabed as shown in Figure 1. The flow velocity profile by the selected unresolved CFD and DEM coupled solver was compared with that of a CFD solver. In the coupled solver simulation, the seabed particles were considered, while, in the simulation of the CFD solver, the seabed particles were not considered. Figure 2 shows the flow velocity profiles simulated by both solvers at a central location in a domain without the pipeline. The height of 0 m means the top of the sand surface. It was confirmed that the desired flow velocity profile could be obtained, even with the particles.  Table 2 shows the simulation parameters. A total number of 134,285 particles was used to form the seabed. The density of the particles was 2600 kg/m3. Young's modulus and Poisson ratio of the particles were 5 × 10 Pa and 0.45, respectively. The friction and rolling friction coefficients were 0.6 and 0.125, respectively. The time step for the CFD and DEM solvers were 2 × 10 s and 5 × 10 s, respectively. Every forty steps, both CFD and DEM results were coupled. The computational time for a single test took 10 days on a workstation equipped with Intel Xeon(R) Gold 6226R CPU @2.9 GHz processors.  Table 2 shows the simulation parameters. A total number of 134,285 particles was used to form the seabed. The density of the particles was 2600 kg/m3. Young's modulus and Poisson ratio of the particles were 5 × 10 6 Pa and 0.45, respectively. The friction and rolling friction coefficients were 0.6 and 0.125, respectively. The time step for the CFD and DEM solvers were 2 × 10 −4 s and 5 × 10 −6 s, respectively. Every forty steps, both CFD and DEM results were coupled. The computational time for a single test took 10 days on a workstation equipped with Intel Xeon(R) Gold 6226R CPU @2.9 GHz processors.

Angle of Repose for Glass Beads
To verify the numerical methods, the angle of repose for glass beads was simulated. Figure 3 shows the geometrical specifications of the baseplate and hopper. The diameters of top and bottom of the hopper were 250 mm and 50 mm, respectively. The height of the hopper was 400 mm, and the diameter of the baseplate (D b ) was 500 mm. The hopper moved slowly and vertically up from the initial position with a certain speed. The particles discharged from the hopper, accumulated on the baseplate, and gradually expanded to the edge [47].

Angle of Repose for Glass Beads
To verify the numerical methods, the angle of repose for glass beads was simulated. Figure 3 shows the geometrical specifications of the baseplate and hopper. The diameters of top and bottom of the hopper were 250 mm and 50 mm, respectively. The height of the hopper was 400 mm, and the diameter of the baseplate (Db) was 500 mm. The hopper moved slowly and vertically up from the initial position with a certain speed. The particles discharged from the hopper, accumulated on the baseplate, and gradually expanded to the edge [47]. The properties of the particles [47] are summarized in Table 3. The diameter ( ), density ( ), Poisson's ratio ( ), and Young's modulus ( ) were 10 mm, 2500 kg/m , 0.24, and 5.5 × 10 Pa, respectively. The total number of 13,000 particles were used. The measured sliding friction coefficient of the glass beads was 0.142, and the coefficient of rolling friction was assumed to be very small [47].  The properties of the particles [47] are summarized in Table 3. The diameter (d p ), density (ρ p ), Poisson's ratio (ν p ), and Young's modulus (E) were 10 mm, 2500 kg/m 3 , 0.24, and 5.5 × 10 6 Pa, respectively. The total number of 13,000 particles were used. The measured sliding friction coefficient of the glass beads was 0.142, and the coefficient of rolling friction was assumed to be very small [47].  Figure 4 shows the particle motions in the simulation of the angle of repose. The hopper filled with the particles elevated at 5 mm/s, and the gravity particles was discharged. The hopper was stopped at 110 mm height from the baseplate. The simulation was continued until the particles in the hopper were completely released and the pile was in a steady state. Then the number of particles of the stabilized pile (N s ) were obtained. The slope of the stabilized pile was not symmetrical because the released particles from the hopper were not symmetrical, causing the relative motion between the particles at the bottom of the pile and baseplate [48]. steady state. Then the number of particles of the stabilized pile ( ) were obtained. The slope of the stabilized pile was not symmetrical because the released particles from the hopper were not symmetrical, causing the relative motion between the particles at the bottom of the pile and baseplate [48]. The angle of repose ( ) could be calculated as: where H indicates the height of the pile, and Db is the baseplate diameter. Table 4 lists the predicted angle of repose and number of the particles of the stabilized pile. The simulation result showed that the angle of repose and number of particles of the stabilized pile agreed well with the experimental data [47].

Angle of Repose for Sand Particles
In the experiment, non-spherical particles were used. While spherical particles were used in the present study. To simulate the behavior of non-spherical particles by the spherical particles, the coefficient of rolling friction should be adjusted. To decide the rolling friction coefficient ( ) used in the scour simulation, a parametric study for the angle of repose for sand particles was carried out. The sand particle diameter was 1 mm. The density of the particles was 2600 kg/m 3 . Figure 5 shows the sand particle motions in the simulation of the angle of repose. The particles were discharged from the hopper [34]. Table 5 lists the rolling friction coefficients for a sliding friction coefficient. For the sliding friction coefficient, 0.6 was selected. This is because that sliding friction coefficients higher than 0.4 did not show a dominant effect on the angle of repose [48]. The angle of repose increased with increasing rolling friction coefficient. When the sliding friction coefficient was 0.6 and rolling friction coefficient was 0.125, the angle of repose ( ) was 32.02°, which The angle of repose (α) could be calculated as: where H indicates the height of the pile, and D b is the baseplate diameter. Table 4 lists the predicted angle of repose and number of the particles of the stabilized pile. The simulation result showed that the angle of repose and number of particles of the stabilized pile agreed well with the experimental data [47].

Angle of Repose for Sand Particles
In the experiment, non-spherical particles were used. While spherical particles were used in the present study. To simulate the behavior of non-spherical particles by the spherical particles, the coefficient of rolling friction should be adjusted. To decide the rolling friction coefficient (R µ ) used in the scour simulation, a parametric study for the angle of repose for sand particles was carried out. The sand particle diameter was 1 mm. The density of the particles was 2600 kg/m 3 . Figure 5 shows the sand particle motions in the simulation of the angle of repose. The particles were discharged from the hopper [34]. Table 5 lists the rolling friction coefficients for a sliding friction coefficient. For the sliding friction coefficient, 0.6 was selected. This is because that sliding friction coefficients higher than 0.4 did not show a dominant effect on the angle of repose [48]. The angle of repose increased with increasing rolling friction coefficient. When the sliding friction coefficient was 0.6 and rolling friction coefficient was 0.125, the angle of repose (α) was 32.02 • , which corresponded to the angle of repose of the sand [3,49]. From the simulation results, the sliding friction coefficient of 0.6 and rolling friction coefficient of 0.125 were selected.

Incipient Motion of Sand Particles
The incipient motion of sand particles was simulated to validate interactions between the fluid flow and seabed particles. The computational domain for the unresolved CFD and DEM coupled solver is described in Figure 6. The sand particles with a diameter of 1 mm were randomly distributed and densely packed in the seabed region of 4 D length, 1 D width, and 0.25 D depth. The inlet boundary with 4 D height was specified with a uniform flow. The inlet velocity was raised up linearly to target velocity over first 1 s for a numerical stability. The results of the particle motion were presented, at the same time, after the first 1 s. All simulations for the incipient motion were carried out with the same seabed particle arrangement.

Incipient Motion of Sand Particles
The incipient motion of sand particles was simulated to validate interactions between the fluid flow and seabed particles. The computational domain for the unresolved CFD and DEM coupled solver is described in Figure 6. The sand particles with a diameter of 1 mm were randomly distributed and densely packed in the seabed region of 4 D length, 1 D width, and 0.25 D depth. The inlet boundary with 4 D height was specified with a uniform flow. The inlet velocity was raised up linearly to target velocity over first 1 s for a numerical stability. The results of the particle motion were presented, at the same time, after the first 1 s. All simulations for the incipient motion were carried out with the same seabed particle arrangement. The incipient motion of the sand particles was determined using the criteria of Kramer [50]. Figure 7 shows the particle motion for various inlet velocities ( ) in the region between x/D = 2.125 and 4. Figure 7a shows no particle motion at = 0.3 m/s. The motion of very few particles was initiated at isolated spots at = 0.4 m/s. At = 0.5 m/s, several particle motions occurred. The particle motions were increased with increasing inlet velocity, and the particles were moved in groups in the flow direction. According to The incipient motion of the sand particles was determined using the criteria of Kramer [50]. Figure 7 shows the particle motion for various inlet velocities (u ∞ ) in the region between x/D = 2.125 and 4. Figure 7a shows no particle motion at u ∞ = 0.3 m/s. The motion of very few particles was initiated at isolated spots at u ∞ = 0.4 m/s. At u ∞ = 0.5 m/s, several particle motions occurred. The particle motions were increased with increasing inlet velocity, and the particles were moved in groups in the flow direction. According to Kramer [50], a level of weak movement was recommended as the incipient motion of the particles [51]. The weak movement represented the motion of small particles at isolated spots, which corresponded to the result in Figure 7b. In detail, the particle motion within the weak movement level is shown in Figure 8. As the inlet velocity (u ∞ ) increased from 0.4 m/s to 0.52 m/s, and the number of moving single particles increased, forming groups with surrounding particles. When single particles sporadically formed groups, it could be considered a transition from the weak level to the medium level.
The Hjulstrom diagram [52] was employed to define the critical states of the particle motion, which could be used to distinguish between the transport, erosion, and deposition of sediments using the particle size and mean flow velocity. The inlet velocity (u ∞ ) was not appropriate to characterize the flow situation for the sediment transport because the velocity profile changed over the length scale [53]. The mean flow velocity was defined as the depth-  Figure 9 shows the computed depth-averaged velocities with the particle diameter of 1 mm. The simulation results showed that it was predicted close to the critical velocity for the incipient motion of the particle. Kramer [50], a level of weak movement was recommended as the incipient motion of the particles [51]. The weak movement represented the motion of small particles at isolated spots, which corresponded to the result in Figure 7b. In detail, the particle motion within the weak movement level is shown in Figure 8. As the inlet velocity ( ) increased from 0.4 m/s to 0.52 m/s, and the number of moving single particles increased, forming groups with surrounding particles. When single particles sporadically formed groups, it could be considered a transition from the weak level to the medium level.  The Hjulstrom diagram [52] was employed to define the critical states of the particle motion, which could be used to distinguish between the transport, erosion, and deposition of sediments using the particle size and mean flow velocity. The inlet velocity ( ) was not appropriate to characterize the flow situation for the sediment transport because the Kramer [50], a level of weak movement was recommended as the incipient motion of the particles [51]. The weak movement represented the motion of small particles at isolated spots, which corresponded to the result in Figure 7b. In detail, the particle motion within the weak movement level is shown in Figure 8. As the inlet velocity ( ) increased from 0.4 m/s to 0.52 m/s, and the number of moving single particles increased, forming groups with surrounding particles. When single particles sporadically formed groups, it could be considered a transition from the weak level to the medium level.  The Hjulstrom diagram [52] was employed to define the critical states of the particle motion, which could be used to distinguish between the transport, erosion, and deposition of sediments using the particle size and mean flow velocity. The inlet velocity ( ) was not appropriate to characterize the flow situation for the sediment transport because the Here, is the top of the seabed surface. The mean flow velocity (V), which corresponded to the weak movement, ranged from = 0.4 m/s to 0.48 m/s and was calculated and indicated in the Hjulstrom diagram. Figure 9 shows the computed depth-averaged velocities with the particle diameter of 1 mm. The simulation results showed that it was predicted close to the critical velocity for the incipient motion of the particle.

Results and Discussion
Simulations for scour around a subsea pipeline were carried out using the numerical methods used in the simulations of the angles of repose and incipient motion of the par-

Results and Discussion
Simulations for scour around a subsea pipeline were carried out using the numerical methods used in the simulations of the angles of repose and incipient motion of the particles. A grid sensitivity study was performed using three different cases in Table 6. Here, 1, 2, and 3 represent coarse, medium, and fine grids, respectively. The total number of the grid was increased by 2 times. A representative grid size (∆x) was used [54]. The refinement factor for the representative grid size was √ 2. The particle size (d p ) was set consistently to 1 mm for the three cases. As the grid system was fined, the scour depth and drag force acting on the pipeline were converged. The drag coefficient showed the difference less than 5% [55]. The scour depth was converged in the case 2. The evolution of scour depth with time is shown in Figure 10. The depth of the scour hole (S) was normalized by the pipeline diameter (D). The simulation results in the early stage of scour slightly underestimated the scour depth, in which the inlet velocity was gradually increased [38]. In other words, the onset of scour was postponed, due to a low velocity, which corresponded to about t = 4.5 s. The scour depth increased sharply after the onset of scour, which corresponded to the tunnel erosion (about t = 7 s). After the tunnel erosion, the scour depth increased at a relatively lower rate, which corresponded to the leewake erosion (about t = 25 s). Additionally, the scour depth reached the equilibrium state. This trend was consistent with the experimental results [3]. In Yang et al.' simulation [38], the scour depth predicted lower than the experimental data in the onset of scour and tunnel erosion. While the scour depth was in good agreement with the experimental data in the leewake erosion, the selection of the coefficient of rolling friction by the validation test in this study more accurately influenced the prediction in the onset of scour and tunnel erosion.
J. Mar. Sci. Eng. 2022, 10, x FOR PEER REVIEW experimental data in the lee-wake erosion, the selection of the coefficient of rolling by the validation test in this study more accurately influenced the prediction in t of scour and tunnel erosion.  Figure 11 shows the evolution of a bed profile with time. The onset of sc occurred at about = 4.5 s. The onset of scour was initiated by the seepage flow the pipeline and pressure difference between the upstream and the downstream the pipeline [7]. In the single-phase models, and some two-phase models, it was c ing to reproduce the onset of scour without a very small gap between the seabed and the pipeline [13,19]. After the onset of scour, the tunnel erosion proceeded t  Figure 11 shows the evolution of a bed profile with time. The onset of scour was occurred at about t = 4.5 s. The onset of scour was initiated by the seepage flow beneath the pipeline and pressure difference between the upstream and the downstream around the pipeline [7]. In the single-phase models, and some two-phase models, it was challenging to reproduce the onset of scour without a very small gap between the seabed bottom and the pipeline [13,19]. After the onset of scour, the tunnel erosion proceeded to S/D < 0.3 (about t = 7 s), due to an increase in fluid flow entering the opening beneath the pipeline. The eroded particles beneath the pipeline formed a sand dune about half the diameter of the pipeline downstream. During the scour process, a large amount of particle movement at the rear side of the pipeline occurred at the lee-wake erosion (about t = 25 s). In addition, the scour hole was deepened and expanded downstream, while the sand dune moved farther and flattened in the direction of flow. These phenomena were consistent with the experimental observations [3,7].  Figure 11 shows the evolution of a bed profile with time. The onset of sco occurred at about = 4.5 s. The onset of scour was initiated by the seepage flow b the pipeline and pressure difference between the upstream and the downstream the pipeline [7]. In the single-phase models, and some two-phase models, it was ch ing to reproduce the onset of scour without a very small gap between the seabed and the pipeline [13,19]. After the onset of scour, the tunnel erosion proceeded to 0.3 (about = 7 s), due to an increase in fluid flow entering the opening bene pipeline. The eroded particles beneath the pipeline formed a sand dune about diameter of the pipeline downstream. During the scour process, a large amount of movement at the rear side of the pipeline occurred at the lee-wake erosion (about s). In addition, the scour hole was deepened and expanded downstream, while t dune moved farther and flattened in the direction of flow. These phenomena we sistent with the experimental observations [3,7].     Figure 12 shows the contours of the fluid velocity normalized by the inlet velocity. In the onset of scour, the small velocity of the flow in the opening beneath the pipeline appeared. In the tunnel erosion, the small opening beneath the pipeline expanded into the tunnel, due to the increased flow. In the lee-wake erosion, the flow was affected by the seabed, and the deformed seabed region could be observed. Figure 13 shows the velocity of the individual particles. The particles beneath the pipeline started to move at the onset of scour (t = 4.5 s). In the tunnel erosion (t = 7 s), the movement of the particles and number of moving particles around the pipeline were increased, indicating that the seabed under the pipeline was exposed to a strong erosion. In the lee-wake erosion (t = 25 s), the particles in the scour hole showed small movement, while the particles in the downstream of the scour hole moved in the flow direction.
The drag force is the dominant force in the fluid and particle interaction. Figure 14 shows the drag force acting on the individual particles in the flow direction. Until the tunnel erosion, the drag force increased mainly beneath the pipeline, causing the high scour depth rate. In the lee-wake erosion, the drag force applied to the seabed led to the expansion of the scour hole and the movement of many particles. The particle transport capacity in the scour process increased with the increase of the drag force. Figure 15 shows the drag force acting on individual particles in the vertical direction. The drag force acting in the vertical direction on the front and rear sides of the scour hole increased at the early stages and drove the motion of particles on the inclined bed inside the scour hole. From the results, the drag force acting on the particles played an important role in the transport of the sand particles.  Figure 13 shows the velocity of the individual particles. The particles beneath the pipeline started to move at the onset of scour ( = 4.5 s). In the tunnel erosion ( = 7 s), the movement of the particles and number of moving particles around the pipeline were increased, indicating that the seabed under the pipeline was exposed to a strong erosion. In the lee-wake erosion ( = 25 s), the particles in the scour hole showed small movement, while the particles in the downstream of the scour hole moved in the flow direction. The drag force is the dominant force in the fluid and particle interaction. Figure 14 shows the drag force acting on the individual particles in the flow direction. Until the tunnel erosion, the drag force increased mainly beneath the pipeline, causing the high scour depth rate. In the lee-wake erosion, the drag force applied to the seabed led to the expansion of the scour hole and the movement of many particles. The particle transport capacity in the scour process increased with the increase of the drag force. Figure 15 shows the drag force acting on individual particles in the vertical direction. The drag force acting in the vertical direction on the front and rear sides of the scour hole increased at the early stages and drove the motion of particles on the inclined bed inside the scour hole. From the results, the drag force acting on the particles played an important role in the transport of the sand particles.  Figure 14. The x-directional particle drag force ( , ). Figure 14. The x-directional particle drag force ( f d,x ).  Figure 16 shows the total force of the seabed particles. As the particles moved during the scour process, the interaction between the particles increased. In the region where the particle transport was most active, the interaction force between the particles was also the most intense. In Figures 13-15, it was shown that the particle transport was most affected by drag force acting on the particles. However, it should be noted that the particle and particle interaction was also active in the rear side of the sand dune, where the drag force was not dominant. It concluded that the interaction between the particles should be considered in the development process of the scour hole and sand dune in the tunnel erosion. These results show why the CFD and DEM coupled model is more suitable than the single-phase model for scour simulation.  Figure 16 shows the total force of the seabed particles. As the particles moved during the scour process, the interaction between the particles increased. In the region where the particle transport was most active, the interaction force between the particles was also the most intense. In Figures 13-15, it was shown that the particle transport was most affected by drag force acting on the particles. However, it should be noted that the particle and particle interaction was also active in the rear side of the sand dune, where the drag force was not dominant. It concluded that the interaction between the particles should be considered in the development process of the scour hole and sand dune in the tunnel erosion. These results show why the CFD and DEM coupled model is more suitable than the single-phase model for scour simulation.

Conclusions
In this study, the simulations were carried out on scour around the seabed pipeline, in order to investigate the applicability of the CFD and DEM coupled solver [33]. The unsolved CFD and DEM coupled solver was used to track the particle motion caused by the fluid flow. For the verification of the numerical methods, the angles of repose were simulated and compared with the experimental data. The predicted angles of repose were in good agreement with the experimental data, and the rolling friction coefficient of the sand particles was decided for the scour simulation. The rolling friction coefficient should be used to compensate for non-spherical sand particles and determine the proper value for application to scour simulation. To validate the selected numerical model parameters, the incipient motion was simulated and discussed the critical state to validate the interaction between the fluid flow and seabed particles. The incipient motion of the sand particles was predicted close to the critical velocity in the Hjulstrom diagram.
The scour simulations around the subsea pipeline were carried out and compared with the experimental data [3]. From the grid dependency test, the CFD and DEM coupled solver could be applied without grid dependency constraints. The flow fields and interactions between the flow and particle were presented in the three scour stages. At the

Conclusions
In this study, the simulations were carried out on scour around the seabed pipeline, in order to investigate the applicability of the CFD and DEM coupled solver [33]. The unsolved CFD and DEM coupled solver was used to track the particle motion caused by the fluid flow. For the verification of the numerical methods, the angles of repose were simulated and compared with the experimental data. The predicted angles of repose were in good agreement with the experimental data, and the rolling friction coefficient of the sand particles was decided for the scour simulation. The rolling friction coefficient should be used to compensate for non-spherical sand particles and determine the proper value for application to scour simulation. To validate the selected numerical model parameters, the incipient motion was simulated and discussed the critical state to validate the interaction between the fluid flow and seabed particles. The incipient motion of the sand particles was predicted close to the critical velocity in the Hjulstrom diagram.
The scour simulations around the subsea pipeline were carried out and compared with the experimental data [3]. From the grid dependency test, the CFD and DEM coupled solver could be applied without grid dependency constraints. The flow fields and interactions between the flow and particle were presented in the three scour stages. At the onset of the scour stage, the sand particles beneath the pipeline started to move, due to the seepage flow and pressure difference between the upstream and downstream around the pipeline. During the tunnel erosion stage, the sediment transport was significantly increased with the increasing drag force, indicating that the seabed beneath the pipeline was exposed to the strong erosion. The strong erosion by the tunnel flow also expanded the scour hole and formed the sand dune. It was confirmed that the interaction between the particles was most active beneath the pipeline and in the sand dune during the tunnel erosion stage. With entering the lee-wake erosion stage, the particle movement was very small in the scour hole but very active when moving downstream. The main particle motions around the pipeline were generated by the drag force, due to the high flow velocity and the interaction between the particles.
The present CFD and DEM coupled solver required a lot of computational time, but all of the above findings can be explained from the information of each individual particle without using empirical formulas or solving constitutive equations, including sand concentration. The CFD and DEM coupled model showed great potential to simulate particle and particle interactions, as well as the fluid and particle interactions in previous research and the current study.