Three-Dimensional Simulations of Scour around Pipelines of Finite Lengths

In the past few decades, there have been many numerical studies on the scour around offshore pipelines, most of which concern two-dimensional setups, with the pipeline infinitely long and the flow perpendicular to the pipeline. Based on the Ansys FLUENT flow solver, this study establishes a numerical tool to study the three-dimensional scour around pipelines of finite lengths. The user-defined functions are written to calculate the sediment transport rate, update the bed elevation, and adapt the computational mesh to the new boundary. The correctness of the model has been verified against the measurements of the conventional two-dimensional scour around a long pipe and the three-dimensional scour around a sphere. A series of computations are subsequently carried out to discover how the scour hole is dependent on the pipeline length. It is found that the equilibrium scour depth increases with the pipeline length until the pipeline length exceeds four times the pipe diameter.


Introduction
Pipelines are often regarded as the lifeline for the oil and gas industry. Local scour around offshore pipelines may leave pipelines to be suspended above the sea floor, making them vulnerable when exposed to currents, waves, shipment, and fishery activities. Against the backdrop of rapid advancements in offshore resource exploitation, the pipeline scour problem has attracted intensive research in the past few decades. Numerous empirical formulae have been proposed to estimate the threshold condition for the scour to occur, the equilibrium scour depth, and the time scale of its growth under different combinations of the pipeline size, initial embedment, flow condition, and bed sediment property. Underwater surveys on the wellbeing of on-bottom pipelines are often conducted after the passing of a tropical storm, and a risk assessment is then held based on these empirical formulae, pipeline strength, and importance, etc. to decide whether any scour remedial measures are required before the arrival of the next storm.
As early as the 1960s, Herbich et al. (1984) conducted experiments and proposed empirical methods to estimate the scour depth beneath offshore pipelines under the action of waves [1]. Mao et al. (1987) experimentally studied the two-dimensional scour processes around fixed or vibrating pipes subject to steady currents [2]. Truelsen et al. (2005) studied the local scour around a near-bed sphere in both wave and current conditions and analysed the self-bury phenomenon [3]. Among others, Gao et al. (2006) carried out sophisticated experiments to investigate the influence of vortex-induced vibrations on the scour beneath an elastically mounted pipe [4].
Early numerical studies are based on the potential flow theory, which can correctly predict the maximum scour depth and the shape of the upstream part of the scour hole. 2 of 16 However, these previous models cannot explain the gentle slope of the scour hole downstream of the pipe. It has been understood in Sumer et al. (1988) that the gentle slope of the scour hole is mainly due to the flow separation and vortex shedding in the wake [5], which is later confirmed by Li and Cheng (1999) in the numerical prediction of the equilibrium scour hole below pipelines [6]. Chiew (1990) established empirical formulae concerning the flow and sediment transport characteristics through the scour hole under unidirectional flow conditions [7]. Brørs (1999) reproduced the measured scour hole profiles in the numerical simulation of clear-water scour around an on-bottom pipe, although the vortex shedding phenomenon was absent from the simulation [8]. Shamloo et al. (2001) studied the flow and scour around a half-buried sphere, which could be regarded as a habitat structure, exposed to steady currents [9]. Their key findings were that the maximum scour occurred in front of the structure and the maximum scour depth was 0.67 times the structure diameter. An excellent review was conducted by Sumer and Fredsøe (2002) [10], which reported the key research and significant findings on pipeline scour. , and  developed a two-dimensional model that accurately captured the local scour phenomena around the fixed pipelines induced by currents and oscillatory flows, with the inclusion of the effects of both the suspended load and bedload [11,12]. The model emphasized the important role of streaming in unsteady sediment transport and was later applied to predict the scour around vibrating pipes, with the results summarized into nondimensional regression formulae to describe the relationship between the scour depth, vibration frequency, and amplitude [13]. Based on analytical/numerical data, Yang et al. (2005) studied the shear stress distribution on the seabed, the pressure difference between the front and rear of the pipeline, and the free-spanning mechanisms of pipelines [14]. Yang et al. (2006) examined the influence of the wake vortex oscillation on the scour [15]. According to these numerical results, they obtained the seepage field and hydraulic gradient to quantify the influence of the seepage flow. Hatton et al. (2007) used Flow-3D software to numerically simulate the three-dimensional flow around partially buried submarine pipelines subject to combined actions of waves and currents, but the bed was fixed and only the scour potential was inferred [16]. Zhang et al. (2012) relied on the commercial software, FLUENT, to obtain solutions to the Reynolds-averaged Navier-Stokes (RANS) equations and estimate the two-dimensional scour evolution around a horizontal cylinder [17]. The adaptive meshing method was used to fit new bed positions. The model predictions were found to agree well with measurements.
As can be seen, the majority of the existing scour models have focused on the vertical two-dimensional cases. Due to the limitation of the computational power, genuine threedimensional simulations of the scour phenomena are rare. This paper studies the scour development beneath finite-length pipelines under the influence of steady currents, which is inherently three-dimensional. The numerical simulations are based on the secondary development of the commercial software, ANSYS FLUENT, employing the adaptive meshing method and k-ε turbulence model. The User-Defined-Functions (UDFs) of this commercial flow solver allow the computations of the sediment transport and bed deformation. After verification, the model is used to perform parametric studies to reveal the dependence of the scour evolution on the pipe length. The present study offers a convenient way of establishing three-dimensional scour models, which are generic and can also be used to study the scour phenomena around other types of underwater structures.

Flow Equations
The whole mathematical model includes the hydrodynamic submodel and the morphological submodel. It has been well acknowledged that scour around pipelines is initiated by the failure of the underneath soil owing to the pressure difference between the upstream and downstream sides of the pipe. The analysis of the scour initiation requires the study of the seepage flow in the soil. The present study focuses on how scour will proceed once it has initiated, so the onset of scour is bypassed in this study by creating a small opening between the pipeline and the bed. Hence, only the hydrodynamics above the bed need to be considered. The three-dimensional Reynolds-Averaged Navier-Stokes (RANS) and continuity equations, together with k-ε turbulence closure equations, form the basis of the hydrodynamic submodel.
When considering incompressible flows, the three-dimensional RANS equations can be written in the Cartesian coordinates as: The continuity equation can be written as: where x is the spatial coordinate; i,j = 1, 2, 3 are the coordinate indices; t is time; u is the Reynolds-averaged velocity; ρ is density of the fluid; p is pressure; µ is the dynamic viscosity coefficient; τ is the Reynolds stress. The Boussinesq assumption is used to determine the turbulent stress tensor τ ij by introducing the eddy viscosity coefficient µ t : where k is the turbulent kinetic energy, δ . ιj is Kronecker delta, and S ij = 1 is the mean strain rate. The eddy viscosity coefficient, µ t , is considered isotropic and thus, a scalar. The semiempirical standard k − ε model requires the following two additional transport equations to be solved for the turbulence kinetic energy, k, and the turbulence kinetic energy dissipation rate, ε: where G k is the production rate of the turbulence kinetic energy.

Sediment Transport and Bed Deformation Equations
For the convenience of analyses, sediment transport may be conceived to occur in two distinct modes. The bedload refers to the sediment particles transported by rolling, sliding, and saltating along the bed. These moving particles are confined in a thin layer above the bed, and their weight is mainly balanced by the underlying particles. The suspended load refers to the sediment transported by suspension in the moving fluid. The particles in suspension are dispersed over the whole water column, and they are entirely supported by the surrounding fluid. The Van Rijn formula is used to estimate the bedload transport rate [18]: where T is the nondimensional transport stage parameter; q b is the volumetric bedload transport rate per unit width; u * is the friction velocity; u * cr is the critical bed friction velocity for the incipience of sediment motion; s is the specific weight of the sediment material; g is gravitational acceleration 9.8 m/s 2 ; d is the particle diameter; D * is the nondimensional particle size: where v = µ/ρ is the kinematic viscosity coefficient of the water taken to be 10 −6 m 2 /s. The friction velocity, u * , is determined by the bed shear stress, τ w , and the fluid density. It is introduced in the near-wall boundary layer theory mainly for the convenience of expression nondimensional groups.
where τ w = µ(∂u/∂y) y = 0 can be obtained from the computed velocity distribution immediately above the bed. The critical frictional velocity, u * cr , corresponding to the sediment incipience can be estimated by means of the nondimensional Shields parameter, The critical Shields parameter can be estimated with the following empirical equation: where θ cr0 is the critical Shields parameter on a horizontal bed. On an inclined bed, the critical Shields parameter, θ cr0 , can be modified using Allen's approach [19]: where α is the angle between the bed and the horizontal plane; Φ is the repose angle of the submerged sediment. Equation (13) indicates that the critical Shields parameter is increased for sediment moving up a slope, while decreased for sediment moving down a slope. When the slope angle, α, is close to the repose angle of sediment, Φ, the modified critical Shields number, θ cr , will tend to be 0, which may result in numerical instability or mistakes when using Equation (8). To solve this problem, Soulsby replaced the critical Shields parameter with the effective critical Shields number, θ ce [20]: Correspondingly, the Shields parameter, θ, which is related to the friction velocity, u * , is also modified into an effective quantity θ be : cos ω x cos ξ x + cos ω y cos ξ y + cos ω z cos ξ z where ω χ , ω y , and ω Z are the three angles between the flow direction and the three axes, and ξ x , ξ y , and ξ z are the three angles between the maximum slope direction of the bed and the three axes. Then, Equation (8) can be modified as: When the effective critical Shields parameter, θ ce , is exceeded by, θ be , sediment transport takes place. With the above modifications, the calculation of the transport stage parameter, T, is numerically stabilized. From the conservation of mass of the sediment material in the bed, the bed elevation, h b , can be updated as: where po is the porosity of the bed; q b . ι is the component of the bedload transport rate, q b , in the x i direction. For vertical two-dimensional problems, the bed is a curve, and thus, the spatial index is i = 1. For three-dimensional problems, the bed is a surface, and thus, i = 1,2. In Equation (17), De is the flux of the suspended sediment deposition to the bed, and En is the rate of sediment entrainment from the bed. This study only considers the effect from the bedload and neglects the influence of the suspended load.

Adaptive Meshing and Sand-Slide Models
The present study adopts triangular mesh in two-dimensional simulations and tetrahedral mesh in three-dimensional simulations. When using the three-dimensional tetrahedral mesh, the bed surface is divided into triangular elements.
The bed elevation change with time leads to the dynamic change of the bottom boundary of the computational domain. Adaptive meshing enables the computational grid to be updated to fit the new boundary. The adaptive meshing strategy adopted in this paper includes two steps. In the first step, the elevation changes at the centers of the boundary elements are transferred to the boundary nodes. Stresses are stored at the interfaces of the computational cells in FLUENT, hence, the bedload sediment transport rates. Similarly, the bed elevation change, ∆h b , calculated by Equation (17) at each time step, ∆t, corresponds to the value at the center of the triangular element on the bed, which needs to be transferred to the boundary nodes. This study uses the area-weighted method to interpolate the bed elevation changes.
where m is the number of elements surrounding a node; S j is the area of the element j; λ j is the weight for element j; ∆h bnode and ∆h b j are the elevation changes at the node and at the center of the element, j, respectively. In the second step, the internal nodes and elements are updated to adapt to the new bed shape. In this process, we combine the spring smoothing method and the remeshing method. The spring smoothing method assumes that each two neighboring nodes are connected by a fictitious spring, which all together form a spring network. When a boundary grid point incurs an initial displacement as a result of the bed elevation change, this local movement will spread out to its adjacent grid points by means of the elasticity of the springs, which enables a smooth transition of the local mesh. However, using the spring-smoothing method alone may give rise to entanglement when the displacement at the boundary exceeds the element size. Moreover, the quality of the grid may significantly deteriorate by only using the spring-smoothing method. To solve this problem, FLUENT can mark out the grid cells that are distorted too much and then replace them by undertaking a local remeshing. The combination of the spring smoothing and local remeshing methods can ensure the high grid quality during the bed morphological change, which helps improve the accuracy and stability of the computation. Because of the complicated behavior of the sediment transport and its nonlinear interaction with the flow field, positive feedback often occurs in the bed elevation update, which may then lead to unrealistic bed profiles and numerical instabilities. An anomaly often occurs when the predicted slope of the scour hole exceeds the repose angle of sediment. Different kinds of rectification techniques, such as backfilling and sand-slide models, have been proposed to overcome this difficulty [21,22]. A sand-slide model is employed in this study when the local bed slope exceeds the repose angle of sediment. The coordinates of the corresponding seabed nodes are adjusted so that the local bed slope becomes equal to the angle of repose while the centroid of the grid cell and the direction of the steepest slope of the bed remain unchanged. For any node, i, Z old i is adjusted to Z new i : where C is the centroid of a bed triangular element, and → C i represents the vector from the centroid to the unadjusted node, i; A is the direction vector of the fastest bed elevation change before the adjustment. Figure 1 demonstrates a comparison of the local mesh before and after introducing the sand-slide model in a typical example of scour beneath a sphere. There exist some unreasonably sharp bed irregularities without using the sand-slide model, which may lead to instabilities in the ensuing computations. After applying the aforementioned sand-slide model, the bed shape becomes smoother and more realistic. smoothing method. To solve this problem, FLUENT can mark out the grid cells that are distorted too much and then replace them by undertaking a local remeshing. The combination of the spring smoothing and local remeshing methods can ensure the high grid quality during the bed morphological change, which helps improve the accuracy and stability of the computation. Because of the complicated behavior of the sediment transport and its nonlinear interaction with the flow field, positive feedback often occurs in the bed elevation update, which may then lead to unrealistic bed profiles and numerical instabilities. An anomaly often occurs when the predicted slope of the scour hole exceeds the repose angle of sediment. Different kinds of rectification techniques, such as backfilling and sand-slide models, have been proposed to overcome this difficulty [21,22]. A sand-slide model is employed in this study when the local bed slope exceeds the repose angle of sediment. The coordinates of the corresponding seabed nodes are adjusted so that the local bed slope becomes equal to the angle of repose while the centroid of the grid cell and the direction of the steepest slope of the bed remain unchanged. For any node, , is adjusted to : where is the centroid of a bed triangular element, and ⃗ represents the vector from the centroid to the unadjusted node, ; is the direction vector of the fastest bed elevation change before the adjustment. Figure 1 demonstrates a comparison of the local mesh before and after introducing the sand-slide model in a typical example of scour beneath a sphere. There exist some unreasonably sharp bed irregularities without using the sand-slide model, which may lead to instabilities in the ensuing computations. After applying the aforementioned sandslide model, the bed shape becomes smoother and more realistic.

Numerical Method and Computational Procedure
The software, FLUENT, uses the finite volume method to solve the governing equations, with the computational domain divided into non-overlapping control volumes. The second-order upwind scheme is chosen to approximate spatial derivatives, while the firstorder implicit scheme is chosen to approximate temporal derivatives. Within a grid cell, the unknown quantities are assumed to vary linearly, and the spatial derivative and gradient are calculated based on the least-square method. The coupling of pressure and velocity is achieved with the SIMPLEC (semi-implicit method for pressure-linked equations consistent) algorithm, which is one of the SIMPLE's improved algorithms proposed by Van Doormal and Raithby [23].
The UDFs are a set of functions provided to the user, which serve as the interference to the standard computational procedures in FLUENT. In this study, the UDFs are used to construct the sediment transport and scour models. In calculating the sediment transport load, the UDFs use the F_STORAGE_R_N3V (f, t, SV_WALL_SHEAR) function to extract the bed shear stress information acquired by the flow solver. After solving the

Numerical Method and Computational Procedure
The software, FLUENT, uses the finite volume method to solve the governing equations, with the computational domain divided into non-overlapping control volumes. The second-order upwind scheme is chosen to approximate spatial derivatives, while the first-order implicit scheme is chosen to approximate temporal derivatives. Within a grid cell, the unknown quantities are assumed to vary linearly, and the spatial derivative and gradient are calculated based on the least-square method. The coupling of pressure and velocity is achieved with the SIMPLEC (semi-implicit method for pressure-linked equations consistent) algorithm, which is one of the SIMPLE's improved algorithms proposed by Van Doormal and Raithby [23].
The UDFs are a set of functions provided to the user, which serve as the interference to the standard computational procedures in FLUENT. In this study, the UDFs are used to construct the sediment transport and scour models. In calculating the sediment transport load, the UDFs use the F_STORAGE_R_N3V (f, t, SV_WALL_SHEAR) function to extract the bed shear stress information acquired by the flow solver. After solving the bed deformation equation, the bed elevation changes are obtained, based on which the nodes on the bed are moved using the macro, DEFINE_GRID_MOTION. Subsequently, the aforementioned mesh adaptation and sand-slide models are applied to update the mesh.
At the inlet of the domain, the flow is assumed to be fully developed and hydraulically smooth: where the friction velocity is u * = κu 0 / ln (δ/z 0 ); δ is the boundary layer thickness; z 0 = 2.5d/30 is the bottom roughness length; κ = 0.42 is the von Karman constant; u ∞ is the free stream velocity. The turbulence kinetic energy, turbulence length scale, and turbulence dissipation rate at the inlet are: The no-slip wall condition is specified on the solid surfaces inside the domain, including the sand bed and the pipe surface. The outlet flow condition is set at the downstream boundary, and the rigid lid condition is assumed at the top of the computational domain. Due to the different time scales of the flow field and the bed-shape changes, the hydrodynamic model and the morphological model are coupled in a quasi-steady manner. Different time steps, ∆t flow and ∆t bed , are selected in the hydrodynamic and morphological models to reduce the computational cost. The flow field is assumed to be steady between every two consecutive updates of the bed profile. The computational procedure is shown in Figure 2.
bed deformation equation, the bed elevation changes are obtained, based on which the nodes on the bed are moved using the macro, DEFINE_GRID_MOTION. Subsequently, the aforementioned mesh adaptation and sand-slide models are applied to update the mesh.
At the inlet of the domain, the flow is assumed to be fully developed and hydraulically smooth: where the friction velocity is * = к / ln ( / ); is the boundary layer thickness; = 2.5 /30 is the bottom roughness length; к = 0.42 is the von Karman constant; is the free stream velocity. The turbulence kinetic energy, turbulence length scale, and turbulence dissipation rate at the inlet are: The no-slip wall condition is specified on the solid surfaces inside the domain, including the sand bed and the pipe surface. The outlet flow condition is set at the downstream boundary, and the rigid lid condition is assumed at the top of the computational domain. Due to the different time scales of the flow field and the bed-shape changes, the hydrodynamic model and the morphological model are coupled in a quasi-steady manner. Different time steps, Δtflow and Δtbed, are selected in the hydrodynamic and morphological models to reduce the computational cost. The flow field is assumed to be steady between every two consecutive updates of the bed profile. The computational procedure is shown in Figure 2.

Vertical Two-Dimensional Scour around an Infinitely Long Pipe
The computation is conducted in a vertical two-dimensional domain of 20D in width and 3.5D in height. The pipe is located 8D downstream from the inlet and 12D upstream of the outlet. The pipe sits on the bed, with its bottom at the same level as that of the initial bed. The pipe has a diameter of 100 mm. The diameter of the sediment grain, d, is 0.36 mm. The sediment density is 2650 kg/m 3 . The repose angle, Φ, is 30 • . The freestream current velocity, U ∞ , is 0.4 m/s. As a result, the Reynolds number, Re, based on the pipe diameter and the incoming flow velocity, is 40,000 according to Equation (25): where v is the kinematic viscosity coefficient of water. The undisturbed Shields parameter, θ ∞ , can be calculated to be 0.048 by assuming a logarithmic velocity distribution above bed: where u * ∞ is the friction velocity sufficiently upstream of the pipe. The sediment incipient condition can be estimated to be θ cr = 0.034, which is slightly smaller than the undisturbed Shields parameter, indicating a live-bed scour scenario and the dominance of bedload in sediment transport. The depth of the initial scour hole is set to be 0.1D, which is much smaller than the equilibrium scour depth. When using the local mesh reconstruction algorithm in FLUENT to fit the updated bed boundary, the height of the near-bed grid cells is roughly equal to 0.19 mm, giving the value of the nondimensional distance to the wall, y + , of around 5.0. The grid resolution near the bed and pipeline is finer than elsewhere, in order to capture the steep velocity gradient and flow separation. The total number of grid cells is 51,083. Before allowing the bed elevations to change, the calculation is first carried out for a while to obtain a rational flow field corresponding to the initially fixed boundaries. Figure 3 shows the snapshots of the flow field and bed positions near the pipe at different stages of the scour: from the early rapid-development stage to the final steadyequilibrium stage. For clear visualization, only a third of the velocity vectors are plotted. As shown in the figure, the water flow velocity through the small scour hole is around 0.55 m/s, which is larger than the freestream velocity at the inlet. When t = 10 min, a sand hill with a height of 0.3D appears behind the pipeline. Downstream of the pipe, the flow has a large upward component due to the presence of the sand hill. The steady flow circulation in the wake is restrained and also pushed upwards by the upslope of the scour hole. No vortex shedding is formed. Behind the sand hill, a relatively small clockwise vortex is observed due to the sudden expansion of the flow. As the scour develops, the depth of the scour hole increases and the velocity of the flow through the hole decreases. The length of the scour hole increases faster than the scour depth, so overall the bed slope decreases over time. The wake vortex behind the pipeline becomes gradually horizontal and maintains a length of about 0.5D. The height of the sand hill continually decreases as the sand hill moves downstream. As a result, the vortex behind the sand hill gradually weakens and eventually disappears.

Vertical Two-Dimensional Scour around an Infinitely Long Pipe
The computation is conducted in a vertical two-dimensional domain of 20 in width and 3.5 in height. The pipe is located 8D downstream from the inlet and 12D upstream of the outlet. The pipe sits on the bed, with its bottom at the same level as that of the initial bed. The pipe has a diameter of 100 mm. The diameter of the sediment grain, , is 0.36 mm. The sediment density is 2650 kg/m . The repose angle, , is 30°. The freestream current velocity, , is 0.4 m/s. As a result, the Reynolds number, Re, based on the pipe diameter and the incoming flow velocity, is 40,000 according to Equation (25): where is the kinematic viscosity coefficient of water. The undisturbed Shields parameter, , can be calculated to be 0.048 by assuming a logarithmic velocity distribution above bed: where * is the friction velocity sufficiently upstream of the pipe. The sediment incipient condition can be estimated to be = 0.034, which is slightly smaller than the undisturbed Shields parameter, indicating a live-bed scour scenario and the dominance of bedload in sediment transport. The depth of the initial scour hole is set to be 0.1D, which is much smaller than the equilibrium scour depth.
When using the local mesh reconstruction algorithm in FLUENT to fit the updated bed boundary, the height of the near-bed grid cells is roughly equal to 0.19 mm, giving the value of the nondimensional distance to the wall, , of around 5.0. The grid resolution near the bed and pipeline is finer than elsewhere, in order to capture the steep velocity gradient and flow separation. The total number of grid cells is 51,083. Before allowing the bed elevations to change, the calculation is first carried out for a while to obtain a rational flow field corresponding to the initially fixed boundaries. As shown in the figure, the water flow velocity through the small scour hole is around 0.55 m/s, which is larger than the freestream velocity at the inlet. When t = 10 min, a sand hill with a height of 0.3D appears behind the pipeline. Downstream of the pipe, the flow has a large upward component due to the presence of the sand hill. The steady flow circulation in the wake is restrained and also pushed upwards by the upslope of the scour hole. No vortex shedding is formed. Behind the sand hill, a relatively small clockwise vortex is observed due to the sudden expansion of the flow. As the scour develops, the depth of the scour hole increases and the velocity of the flow through the hole decreases. The length of the scour hole increases faster than the scour depth, so overall the bed slope decreases over time. The wake vortex behind the pipeline becomes gradually horizontal and maintains a length of about 0.5D. The height of the sand hill continually decreases as the sand hill moves downstream. As a result, the vortex behind the sand hill gradually weakens and eventually disappears.    Figure 4 compares the predicted, represented by solid lines, and measured, represented by square symbols, bed profiles at different times. It can be clearly observed that the bed deformation is quick, and the scour hole rapidly grows during = 10~30 min, corresponding to the large sediment transport rate around the pipeline. A sand hill is formed as a result of the sand accretion. During this period, the sand hill is seen to move along the flow direction, and the sand accretion height increases. With the increase of time beyond = 30 min, the bed deformation slows down. At = 200 min, the scour hole continues to deepen, and the depth tends to the maximum value and the sand hill nearly diminishes, although a small backward-facing step at the bed can still be noticed. At = 370 min, the scour hole approaches the steady state, and the local bed deformation is negligible. Considering the complexity of the sediment transport and scour phenomena, the simulation results can be regarded to agree well with the experimental results, with only small discrepancies in the scour depth beneath the pipe and the sand hill shape downstream of the pipe.   Figure 4 compares the predicted, represented by solid lines, and measured, represented by square symbols, bed profiles at different times. It can be clearly observed that the bed deformation is quick, and the scour hole rapidly grows during t = 10 ∼ 30 min, corresponding to the large sediment transport rate around the pipeline. A sand hill is formed as a result of the sand accretion. During this period, the sand hill is seen to move along the flow direction, and the sand accretion height increases. With the increase of time beyond t = 30 min, the bed deformation slows down. At t = 200 min, the scour hole continues to deepen, and the depth tends to the maximum value and the sand hill nearly diminishes, although a small backward-facing step at the bed can still be noticed. At t = 370 min, the scour hole approaches the steady state, and the local bed deformation is negligible. Considering the complexity of the sediment transport and scour phenomena, the simulation results can be regarded to agree well with the experimental results, with only small discrepancies in the scour depth beneath the pipe and the sand hill shape downstream of the pipe.   Figure 4 compares the predicted, represented by solid lines, and measured, represented by square symbols, bed profiles at different times. It can be clearly observed that the bed deformation is quick, and the scour hole rapidly grows during = 10~30 min, corresponding to the large sediment transport rate around the pipeline. A sand hill is formed as a result of the sand accretion. During this period, the sand hill is seen to move along the flow direction, and the sand accretion height increases. With the increase of time beyond = 30 min, the bed deformation slows down. At = 200 min, the scour hole continues to deepen, and the depth tends to the maximum value and the sand hill nearly diminishes, although a small backward-facing step at the bed can still be noticed. At = 370 min, the scour hole approaches the steady state, and the local bed deformation is negligible. Considering the complexity of the sediment transport and scour phenomena, the simulation results can be regarded to agree well with the experimental results, with only small discrepancies in the scour depth beneath the pipe and the sand hill shape downstream of the pipe.      A sinusoidally shaped initial scour hole, which is 1.0 long and 0.1 deep, is initially introduced beneath the pipeline when the undisturbed Shields parameter, , is more than 0.04. A relatively deep initial scour hole is introduced to avoid computational difficulties and speed up the scour process to reach the equilibrium stage. However, such a deep initial scour hole may lead to exaggerated scour depth in the early stage, especially in situations with small equilibrium scour depths. If the undisturbed Shields parameter, , is smaller than 0.04, then the initial scour hole length and depth are reduced to 0.8 and 0.03 , respectively, considering the relatively shallow equilibrium scour holes in      A sinusoidally shaped initial scour hole, which is 1.0 long and 0.1 deep, is initially introduced beneath the pipeline when the undisturbed Shields parameter, , is more than 0.04. A relatively deep initial scour hole is introduced to avoid computational difficulties and speed up the scour process to reach the equilibrium stage. However, such a deep initial scour hole may lead to exaggerated scour depth in the early stage, especially in situations with small equilibrium scour depths. If the undisturbed Shields parameter, , is smaller than 0.04, then the initial scour hole length and depth are reduced to 0.8 and 0.03 , respectively, considering the relatively shallow equilibrium scour holes in A sinusoidally shaped initial scour hole, which is 1.0D long and 0.1D deep, is initially introduced beneath the pipeline when the undisturbed Shields parameter, θ ∞ , is more than 0.04. A relatively deep initial scour hole is introduced to avoid computational difficulties and speed up the scour process to reach the equilibrium stage. However, such a deep initial scour hole may lead to exaggerated scour depth in the early stage, especially in situations with small equilibrium scour depths. If the undisturbed Shields parameter, θ ∞ , is smaller than 0.04, then the initial scour hole length and depth are reduced to 0.8D and 0.03D, respectively, considering the relatively shallow equilibrium scour holes in those conditions. The roughness height of the bed is equal to the sediment particle diameter.

Three-Dimensional Scour around a Sphere
As for the case of the pipeline scour at the early stage, sediment is eroded quickly from the small gap between the underneath of the pipe and the bed and a small sand hill is formed as the eroded sediment deposits immediately behind the sphere. An illustration of the scour hole shapes observed in the experiment and predicted in the simulation is given in Figure 6, where the flow is to the right as indicated. Since the information about the lapse of time is missing when taking the photograph of scour, it only allows a qualitive comparison with the simulation result. The computational model successfully captures the general shape of the scour hole beneath the sphere and the sediment deposition behind it. those conditions. The roughness height of the bed is equal to the sediment particle diameter.
As for the case of the pipeline scour at the early stage, sediment is eroded quickly from the small gap between the underneath of the pipe and the bed and a small sand hill is formed as the eroded sediment deposits immediately behind the sphere. An illustration of the scour hole shapes observed in the experiment and predicted in the simulation is given in Figure 6, where the flow is to the right as indicated. Since the information about the lapse of time is missing when taking the photograph of scour, it only allows a qualitive comparison with the simulation result. The computational model successfully captures the general shape of the scour hole beneath the sphere and the sediment deposition behind it.  Figure 7 shows a quantitative comparison between the predicted and measured bed profiles along the centerline of the computational domain when = 0.12. The calculated result agrees precisely with the experimental result in the regions beneath and in front of the sphere. However, the simulated height of the sand hill downstream of the sphere is larger than the measured value, which may be attributed to the neglection of the suspend load in the simulation. When the undisturbed Shields parameter is as large as 0.12, the suspended load has a relatively large influence on the sediment transport in the wake region. The model proposed in this paper ignores the effect of suspended sediment by assuming the bedload dominates sediment transport, which works fine when the Shields parameter is close to the threshold value; thus, it is not surprising that the simulated shape of the sand hill behind the sphere demonstrates some degree of inaccuracy in this case. The bed profiles upstream and downstream of a sphere are more symmetrical than those of a pipeline, which can be explained by the differences between flows around the threedimensional and two-dimensional structures.  Figure 8 shows the dependence of the equilibrium scour depth on the undisturbed Shields parameter. The simulation results generally agree with the experimental results. As flow intensifies, both the scour hole depth and width increase. In the clear water scour case, the equilibrium scour depth increases rapidly with the increase of the undisturbed Shields parameter. When is 0.04∼0.045, the growth rate of the equilibrium scour  Figure 7 shows a quantitative comparison between the predicted and measured bed profiles along the centerline of the computational domain when θ ∞ = 0.12. The calculated result agrees precisely with the experimental result in the regions beneath and in front of the sphere. However, the simulated height of the sand hill downstream of the sphere is larger than the measured value, which may be attributed to the neglection of the suspend load in the simulation. When the undisturbed Shields parameter is as large as 0.12, the suspended load has a relatively large influence on the sediment transport in the wake region. The model proposed in this paper ignores the effect of suspended sediment by assuming the bedload dominates sediment transport, which works fine when the Shields parameter is close to the threshold value; thus, it is not surprising that the simulated shape of the sand hill behind the sphere demonstrates some degree of inaccuracy in this case. The bed profiles upstream and downstream of a sphere are more symmetrical than those of a pipeline, which can be explained by the differences between flows around the three-dimensional and two-dimensional structures. those conditions. The roughness height of the bed is equal to the sediment particle diameter.
As for the case of the pipeline scour at the early stage, sediment is eroded quickly from the small gap between the underneath of the pipe and the bed and a small sand hill is formed as the eroded sediment deposits immediately behind the sphere. An illustration of the scour hole shapes observed in the experiment and predicted in the simulation is given in Figure 6, where the flow is to the right as indicated. Since the information about the lapse of time is missing when taking the photograph of scour, it only allows a qualitive comparison with the simulation result. The computational model successfully captures the general shape of the scour hole beneath the sphere and the sediment deposition behind it.  Figure 7 shows a quantitative comparison between the predicted and measured bed profiles along the centerline of the computational domain when = 0.12. The calculated result agrees precisely with the experimental result in the regions beneath and in front of the sphere. However, the simulated height of the sand hill downstream of the sphere is larger than the measured value, which may be attributed to the neglection of the suspend load in the simulation. When the undisturbed Shields parameter is as large as 0.12, the suspended load has a relatively large influence on the sediment transport in the wake region. The model proposed in this paper ignores the effect of suspended sediment by assuming the bedload dominates sediment transport, which works fine when the Shields parameter is close to the threshold value; thus, it is not surprising that the simulated shape of the sand hill behind the sphere demonstrates some degree of inaccuracy in this case. The bed profiles upstream and downstream of a sphere are more symmetrical than those of a pipeline, which can be explained by the differences between flows around the threedimensional and two-dimensional structures.    Figure 8 shows the dependence of the equilibrium scour depth on the undisturbed Shields parameter. The simulation results generally agree with the experimental results. As flow intensifies, both the scour hole depth and width increase. In the clear water scour case, the equilibrium scour depth increases rapidly with the increase of the undisturbed Shields parameter. When θ ∞ is 0.04∼0.045, the growth rate of the equilibrium scour depth is significant. The increase rate reaches maximum when θ ∞ is approximately equal to the critical Shields parameter for the sediment motion, θ cr , and at this point the value of the normalized scour depth, S/D, exceeds 0.2. A further increase in the flow speed does not lead to a significant increase in the equilibrium scour depth, as it is the change of the sediment transport rate, rather than the sediment transport rate itself, that decides the scour.
depth is significant. The increase rate reaches maximum when is approximately to the critical Shields parameter for the sediment motion, , and at this point the of the normalized scour depth, / , exceeds 0.2. A further increase in the flow does not lead to a significant increase in the equilibrium scour depth, as it is the cha the sediment transport rate, rather than the sediment transport rate itself, that decid scour. When the value of the undisturbed Shields parameter is 0.04, the calculated eq rium scour depth shows an apparent difference from the measured value in the e ment, which may be due to the relatively deep scour hole prescribed at t = 0 in the lation. However, the comparison is still acceptable considering the large uncertain the experiments and empirical formulae used in the computation. In the case of Shields numbers, = 0.045, 0.06, 0.08, 0.12, the calculated equilibrium scour dept in excellent agreement with the experimental results, with errors no more than ± The validity of the adaptive-meshing and time-marching methods stated in Sectio well as the correctness of the UDFs, can, thus, be confirmed.

Scour around Pipes of Finite Lengths
This section focuses on the application of the established model to investiga scour around pipes of finite lengths. Figure 9 shows the closeup of an example me jacent to the pipe at the beginning of a simulation. The localized refinement of the near the cylindrical and end surfaces of the pipe can be clearly seen. A cuboid-s computational domain is created. A fully developed boundary layer flow, with lo velocity distribution along the positive -axis, enters the computational domain alo positive -axis from the left. A pipeline with a finite length is located inside the d with its bottom at the same level as the original bed level. The length of the pipe from 2 to 6 . The flow around the pipe is symmetric about the middle − pl the pipe, so the computational domain can be halved using the symmetric boundar dition. The inlet boundary is placed 5 upstream of the pipe and the outflow bou is placed 12 downstream of the pipe, to enable the sufficient length for the flow scour to evolve around the pipe. The height of the domain along the -axis is 3.5 width of the domain along the -axis is 5 , which guarantees an over 2D-wide g tween the end of the pipe and the side boundary to limit the sidewall effects on th and scour simulations. The initial manmade scour hole is 0.1D deep and 1D long, extends beyond the end of the pipe by 0.5D, as shown in Figure 9. The total num grid cells ranges from 320,000 to 420,000, depending on the length of the pipe. When the value of the undisturbed Shields parameter is 0.04, the calculated equilibrium scour depth shows an apparent difference from the measured value in the experiment, which may be due to the relatively deep scour hole prescribed at t = 0 in the simulation. However, the comparison is still acceptable considering the large uncertainties in the experiments and empirical formulae used in the computation. In the case of other Shields numbers, θ ∞ = 0.045, 0.06, 0.08, 0.12, the calculated equilibrium scour depths are in excellent agreement with the experimental results, with errors no more than ±0.01D. The validity of the adaptive-meshing and time-marching methods stated in Section 2, as well as the correctness of the UDFs, can, thus, be confirmed.

Scour around Pipes of Finite Lengths
This section focuses on the application of the established model to investigate the scour around pipes of finite lengths. Figure 9 shows the closeup of an example mesh adjacent to the pipe at the beginning of a simulation. The localized refinement of the mesh near the cylindrical and end surfaces of the pipe can be clearly seen. A cuboid-shaped computational domain is created. A fully developed boundary layer flow, with log-law velocity distribution along the positive z-axis, enters the computational domain along the positive y-axis from the left. A pipeline with a finite length is located inside the domain with its bottom at the same level as the original bed level. The length of the pipe varies from 2D to 6D. The flow around the pipe is symmetric about the middle Y − Z plane of the pipe, so the computational domain can be halved using the symmetric boundary condition. The inlet boundary is placed 5D upstream of the pipe and the outflow boundary is placed 12D downstream of the pipe, to enable the sufficient length for the flow and scour to evolve around the pipe. The height of the domain along the z-axis is 3.5D. The width of the domain along the x-axis is 5D, which guarantees an over 2D-wide gap between the end of the pipe and the side boundary to limit the sidewall effects on the flow and scour simulations. The initial manmade scour hole is 0.1D deep and 1D long, and it extends beyond the end of the pipe by 0.5D, as shown in Figure 9. The total number of grid cells ranges from 320,000 to 420,000, depending on the length of the pipe.    From Figures 10 and 11 and simulations performed on the pipelines of other lengths, multiple aspects of the scour hole development can be seen. At the beginning of the scour, the scour depth beneath the pipe increases rapidly. The sand hill behind the pipe is formed as a result of the sediment erosion beneath the pipe and the subsequent deposition. The sand hill is generally higher and migrates downstream faster near the symmetric − plane, as compared to the sand hill near the end of the pipe. An exception is observed around a short pipe at the early stage, as evidenced in Figure 10a,b, where the most significant sand hill is seen to appear near the pipe end. Later, however, this portion of the sand hill is slowly washed away by the end vortex, and the bed slope in this region becomes mild. As seen in Figure 10c,d, its height later becomes lower than that near the − symmetric plane. When the pipe length is long, the sand hill near the symmetric plane has a similar height to that near the pipe end. The general height of the sand hill varies slightly with the length of the pipe and remains at about 0.25 for the two cases presented in these figures.
As the length of the pipe increases, the rate of the scour hole development increases, including the scour hole depth and width. In addition, the sand hill migrates more rapidly downstream of the pipe as the pipe length increases. When the length of the pipe is 2 , the scour and the flow are strongly three-dimensional around the pipe, leading to the variation of the scour hole shape along the span direction. For a longer pipe, the scour depth and width are relatively unchanged along the centerline of the pipe except near the pipe end, indicating largely vertical two-dimensional characteristics in the middle sections of the pipe. Figure 12 shows the trend of equilibrium scour depth versus the pipe length with = 0.048. With the increase of the pipe length, the equilibrium scour depth increases, and the effect of the pipe end on flow and scour near the symmetric plane gradually diminishes. When the length of the pipe exceeds 4D, the equilibrium scour depth is close to the two-dimensional results, as reported in Mao (1988) [2], and the scour depth then becomes independent of the pipe length. From Figures 10 and 11 and simulations performed on the pipelines of other lengths, multiple aspects of the scour hole development can be seen. At the beginning of the scour, the scour depth beneath the pipe increases rapidly. The sand hill behind the pipe is formed as a result of the sediment erosion beneath the pipe and the subsequent deposition. The sand hill is generally higher and migrates downstream faster near the symmetric y − z plane, as compared to the sand hill near the end of the pipe. An exception is observed around a short pipe at the early stage, as evidenced in Figure 10a,b, where the most significant sand hill is seen to appear near the pipe end. Later, however, this portion of the sand hill is slowly washed away by the end vortex, and the bed slope in this region becomes mild. As seen in Figure 10c,d, its height later becomes lower than that near the y − z symmetric plane. When the pipe length is long, the sand hill near the symmetric plane has a similar height to that near the pipe end. The general height of the sand hill varies slightly with the length of the pipe and remains at about 0.25D for the two cases presented in these figures.
As the length of the pipe increases, the rate of the scour hole development increases, including the scour hole depth and width. In addition, the sand hill migrates more rapidly downstream of the pipe as the pipe length increases. When the length of the pipe is 2D, the scour and the flow are strongly three-dimensional around the pipe, leading to the variation of the scour hole shape along the span direction. For a longer pipe, the scour depth and width are relatively unchanged along the centerline of the pipe except near the pipe end, indicating largely vertical two-dimensional characteristics in the middle sections of the pipe. Figure 12 shows the trend of equilibrium scour depth versus the pipe length with θ ∞ = 0.048. With the increase of the pipe length, the equilibrium scour depth increases, and the effect of the pipe end on flow and scour near the symmetric plane gradually diminishes. When the length of the pipe exceeds 4D, the equilibrium scour depth is close to the two-dimensional results, as reported in Mao (1988) [2], and the scour depth then becomes independent of the pipe length. From Figures 10 and 11 and simulations performed on the pipelines of other lengths, multiple aspects of the scour hole development can be seen. At the beginning of the scour, the scour depth beneath the pipe increases rapidly. The sand hill behind the pipe is formed as a result of the sediment erosion beneath the pipe and the subsequent deposition. The sand hill is generally higher and migrates downstream faster near the symmetric − plane, as compared to the sand hill near the end of the pipe. An exception is observed around a short pipe at the early stage, as evidenced in Figure 10a,b, where the most significant sand hill is seen to appear near the pipe end. Later, however, this portion of the sand hill is slowly washed away by the end vortex, and the bed slope in this region becomes mild. As seen in Figure 10c,d, its height later becomes lower than that near the − symmetric plane. When the pipe length is long, the sand hill near the symmetric plane has a similar height to that near the pipe end. The general height of the sand hill varies slightly with the length of the pipe and remains at about 0.25 for the two cases presented in these figures.
As the length of the pipe increases, the rate of the scour hole development increases, including the scour hole depth and width. In addition, the sand hill migrates more rapidly downstream of the pipe as the pipe length increases. When the length of the pipe is 2 , the scour and the flow are strongly three-dimensional around the pipe, leading to the variation of the scour hole shape along the span direction. For a longer pipe, the scour depth and width are relatively unchanged along the centerline of the pipe except near the pipe end, indicating largely vertical two-dimensional characteristics in the middle sections of the pipe. Figure 12 shows the trend of equilibrium scour depth versus the pipe length with = 0.048. With the increase of the pipe length, the equilibrium scour depth increases, and the effect of the pipe end on flow and scour near the symmetric plane gradually diminishes. When the length of the pipe exceeds 4D, the equilibrium scour depth is close to the two-dimensional results, as reported in Mao (1988) [2], and the scour depth then becomes independent of the pipe length.

Conclusions
This study establishes a generic three-dimensional computer model for investigating the scour around underwater structures. The predicted scour processes around an infinitely long pipe and around a sphere are found to be in agreement with those observed in laboratory experiments. Numerical experiments are then conducted to examine the influences of the length of the pipe on the scour depth and sandhill development. The entire model has shown encouraging ability for predicting the local seabed evolution around pipelines and spheres.
As the scour hole deepens, the sand hill downstream of the pipe first grows and then gradually decreases in height as it travels downstream. As the length of the pipe increases, the rate of the scour development increases, and the speed of the sand hill migration downstream also increases. The maximum height of the sand hill appears first near the pipe end for short pipes, while the height of the sand hill is almost constant from the middle symmetric plane to the pipe end when the pipe is long. This can be explained by the flow field characteristics, especially the structure of the wake and end vortices. It has been found that the maximum scour depth increases with the pipeline length until the pipeline length exceeds four times the pipeline diameter. The variation of the local scour width along the spanwise direction is more apparent when the pipe length is smaller. The bed slope is generally milder near the pipe end than that near the middle symmetric plane. It should be noted that these conclusions are obtained from limited parametric studies. More simulations are needed to confirm these findings, especially using the parameters in the prototype conditions. Using UDFs, this paper introduces a quick way of building a versatile scour model based on a commercial flow solver, FLUENT. A drawback of this strategy is the relatively high computational cost, as the computations of sediment transport and bed deformation and the updates of the computational boundary and mesh are achieved through external interventions to the existing solution procedure of the flow field, rather than seamlessly integrating different submodels at the original source code level.