Numerical Investigation of the Influence of Water Jumping on the Local Scour beneath a Pipeline under Steady Flow

Rigid-lid approximation is usually used to replace the free surface in scour simulation. The influence of the rigid lid assumption on the prediction precision of scour hole in steady flow is studied in this paper. Firstly, a local scour model was constructed based on the open sources Computational Fluid Dynamics (CFD) model OpenFOAM, where both the bed load and suspended load were considered. In the present model, the bed shear stress was calculated by the Newton shear stress formula, instead of the traditional calculation method with the assumption that the flow velocity in vertical direction complies with a logarithmic distribution. The Volume of Fluid (VOF) method was used to capture the free surface and a moving-mesh method was used to track the change of bed surface. Then, several experiments were chosen to validate the model, and the modeling results fitted well with the measured data. Lastly, the effect of the rigid lid assumption on surface elevation, bed shear stress and the profile of the scour hole in steady flow are studied. The result shows that the surface elevation suffers a drop above the pipeline, and the difference of surface elevation between the upstream and downstream increases with decreasing dimensionless depth. Compared with the free surface condition, the bed shear stress and scour hole depth computed with the rigid lid approximation were underestimated.


Introduction
When a pipeline is placed on an erodible sand bed, a scour hole will appear beneath it.This scour hole does great harm to the safety of the pipeline.Many numerical models and physical experiments have been carried out to study the hydrodynamic characteristics and the transport of sediment near the pipeline during the past decades.With respect to studies of the local scour caused by pipelines under steady flow, Mao [1] conducted laboratory experiments and pointed out that the whole pipeline scouring process can be divided into three stages.The first stage is the onset of erosion, which is caused by the action of vortices and the underflow.The second stage is jet scour.As a substantial amount of water goes through the gap under the pipeline, the hole becomes deeper and extends.With the increase in the scour hole, the third stage occurs, namely, the lee-wake scour, which is important to the development of the downstream scour hole.In addition, Jensen et al. [2] studied the erosion of pipeline under waves and found that the Keulegan-Carpenter (KC) number plays a significant role in affecting the slope of scour holes.Matteo et al. [3,4] studied the sediment particle dynamics and turbulence structures around a submarine pipeline under wave with Particle Tracking Velocimetry (PTV) and the relationships between dimensionless scour depth (S/D) and the Keulegan-Carpenter (KC) value.Liang et al. [5,6] studied the erosion of three-dimensional pipelines under steady current and wave-current conditions with physical models, where the incoming flow was not perpendicular to the pipeline.Zhang et al. [7] adopted physical models to study the scour hole of a pipeline under unsteady flows.Burak et al. [8] studied the scour beneath pipelines under irregular wave attacks using physical experiments.
On the other hand, the local scour numerical models also underwent considerable development.Liang et al. [9] studied the scour under a pipeline in a steady flow based on a finite difference method model.Fuhrman et al. [10] simulated the backfilling process caused by waves with a numerical model.Zhao et al. [11] studied the effect of the distance between two pipelines on the scour hole under them.Liu [12] studied nonlinear wave-induced scour under the pipeline on a slope, and Bjarke et al. [13] simulated the erosion under a pipeline under wave-current coupled conditions.
Water jumping is a common free surface phenomenon when flow meets the pipeline, while in most scour simulations [9,11,13,14], the rigid lid assumption was adopted to simplify the simulation.As for using the rigid lid assumption to replace free surface, Baykal et al. [14] and Roulund et al. [15] pointed out that rigid lid approximation can provide reasonable results when the Froude number is less than 0.2.Zhao et al. [16] studied the scour under a pipeline under the vibration condition, and the free surface model was considered in his study.He found that the surface elevation drops above the pipeline, and the difference of water levels between the upstream and downstream increases with an increase in the Froude number.The maximum difference reaches 0.3D when Froude number Fr = 0.29 and ratio of depth to the radius equals 2.5.Liu [17] studied the free surface effects on the scour under waves and found that it is unreasonable to use the rigid-lid approximation as the wave-induced scour around large piles are simulated.A Froude number Fr = 0.2 is usually regarded as a critical parameter to distinguish whether the rigid lid assumption is suitable for scour simulations.However, the errors of results induced by rigid lid approximation when the Froude number is around 0.2 have never been studied.
To study this, a numerical scour model is constructed based on the open source CFD model (OpenFOAM) firstly, where the bed load and suspended load sediment transport model are fully coupled with the flow model.The finite volume method is used to discrete the equations, the moving mesh method is used to simulate the change of bed surface, the free surface is captured by VOF method and the Finite Area Method (FAM) is used to map the parameters from a three-dimensional hydrodynamic model to a two-dimensional morphology model for solving the Exner equation.The sand slide model is used to smooth the bed surface, and more details can be seen in Fan et al. (2017) [18].Then, the bed shear stress, suspended load model and scour model are validated by several experiments.The result of validation shows that the constructed numerical model simulated the scour accurately.The effect of dimensionless depth on bed shear stress and the errors caused by using rigid lid assumption are studied.The result shows that the dimensionless depth plays a significant role in impacting the computed results.For the channel without pipeline, the errors caused by rigid lid approximation under steady flow increase with the improvement of the Froude number.When the Froude number reaches 0.3, the errors between rigid lid assumption and free surface are 6%.As for the channel with a pipeline, the surface elevation suffers a drop above the pipeline, and the difference of surface elevation between the upstream and downstream increases with the decreasing of dimensionless depth (H/D), namely, the ratio of depth (H) to the diameter of pipeline D. Compared with the free surface condition, the bed shear stress and scour depth with rigid lid assumption were underestimated.

Flow Model
The basic governing equations of flow simulation are incompressible continuity equation and momentum equation, namely, the Reynolds average Navier-Stokes equations [17]: where u is the velocity field in Cartesian coordinates, ρ is the multiphase flow density, P is the total pressure, g is the gravitational acceleration, µ is the dynamic molecular viscosity, τ is the Reynolds stress tensor.K γ is the curvature of the interface, σ T is the surface tension constant, and γ is the volume fraction of fluid.

Turbulence Closure Model
For the standard k-ε turbulence closure model, where a two-phase density item was considered, can be described as below: where ρ is density of two-phase flow, α w is the volume rate of water, ρ w is the density of water, ρ a is the density of air, k is the turbulence kinetic energy, ε is the dissipation rate of turbulence, C 1ε = 1.44 , C 2ε = 1.99, σ ε = 1.3 and µ t is turbulence viscosity.

Sediment Transport Equations
The sediment transport modes include bed load and suspended load.The bed load is a thin layer moving in continuous contact with the bed, while the suspended load is in the water transported by the flow.

Bed Load Transport Equations
During the past decades, many bed load transport formulas have been proposed, such as those by Van Rijn [19], Simith et al. [20], Mayer et al. [21] and Soulsby et al. [22].In the present work, the one proposed by Huai et al. [23] is adopted, where the buoyancy caused by flow is considered. ) θ cr = τ cr g(s − 1)d 50 (10) Water 2017, 9, 642 where q b is the bed load transport rate per unit width; θ b is the shields number; g is gravitational acceleration; d is the median grain diameter; s is the relative density of sediment; and θ cr is the critical shields number, which can also be calculated by Equation (11): where θ cro is the critical shields number for a flat bed; φ is the angle between the bed's steepest slope and the bed's shear stress direction; β is the angle between the bed and horizontal plane; µ s is the static friction coefficient; and D * is the dimensionless grain size.

Suspended Load Transport Equations
For the suspended load model, its governing equation is a classical convection-diffusion equation [24].
where c is the sediment concentration; U is the flow velocity field; w s is sediment fall velocity in still water; and ν t is diffusivity of the sediment, which equals turbulence eddy viscosity in the present work.

Boundary Condition
As for the bottom boundary condition of sediment concentration, the suspended concentration at a reference height, which has been recognized as the boundary of bed load and suspended load, is specified [25].The sediment concentration at this reference height can be computed as below: where C b is the sediment concentrate at a reference level, ρ s is the density of sand; ∆ b is the distance from the reference level to the bed, T is the non-dimensional excess shear stress, which can be calculated by Equation (16).
where µ sc is the efficiency factor current, τ sc is the bed shear stress caused by current, and τ cr is the critical bed shear stress.Other details about these parameters can be found in the user manual of Delft3D-Flow [26].

Scour Model
The change of bed elevation can be calculated by solving the Exner equation [27]: where η is the bed elevation; n is the porosity of sand; q b is the bed load transport rate; and D is the deposition rate; E is the erosion rate.
Water 2017, 9, 642 5 of 23 There are several empirical models for calculating the erosion rate; all of them are perform well.In this paper, the model by van Rijin [20] was used, which can be described as below: where C b is the sediment concentration near the bed at a specific location, η is the vector direction that is perpendicular to the bed surface.

Dynamic Mesh Model
As mentioned before, the dynamic mesh method is used to capture the change of bed elevation.After the change of bed elevation is solved by Equation ( 17), the Laplacian equation will be solved to compute the displacement of inner nodes of the mesh: where ν 1 is the displacement field of mesh points; γ is the diffusion coefficient, which is used to control the mesh motion; the inverse distance diffusion coefficient is chosen in this paper, which is more suitable for scour simulation.

Numerical Schemes of the Flow Model
The governing equations of flow simulation are the Navier-Stokes (N-S) Equations, namely, continuity equation and moment equations.For solving the N-S Equations, a turbulence closure model like k-epsilon, k-omega, RNG k-epsilon, LES model and so on is necessary.To discrete the N-S Equations, the finite volume method was used, which is good for mass conservation.The upwind scheme was used to discrete the divergence items, and the Gauss linear scheme was used to discrete the Laplace item.
Nowadays, there are three approaches for solving the pressure-velocity coupling formula due to the improvement in numerical methods.They are PISO (pressure implicit splitting of operators), SIMPlE (semi-implicit method for pressure linked equation) and PIMPLE (Coupling method of PISO and SIMPLE).In the present studies, the PIMPLE scheme was used.
To solve the bed boundary layer, the mesh grid near the bed needs to be fine enough in common CFD codes since the boundary layer is quite thin usually.In the present scour model, the wall function is used to resolve the boundary layer problem.The computational zone is divided into two parts: one is the boundary layer zone, which is resolved by empirical equations.The other layer is the turbulent zone, which is calculated through the turbulence closure model and N-S Equations.With this method, the first mesh layer does not need to be seated in the boundary layer.

Numerical Schemes for Sediment Transport and Scour Model
The bed load and suspended load are independent variables, and they are solved successively in the sediment transport model.The bed shear stress, shields number and the bed load transport rate are computed by flow velocity.When flow is simulated, the bed load can be calculated in a straightforward manner.The suspended load equation is a convection-diffusion equation that is transported by flow field.The finite volume method is used to solve the equation.The erosion rate is calculated through an empirical formula when the flow model computation is finished.The deposition rate of a suspended load can be computed after the convention-diffusion equation is solved.
The Exner Equation will be solved to obtain the bed surface based on the parameters of sediment transport calculated by former calculations.The Exner Equation is a two-dimensional problem, while the flow and sediment transport parameters are calculated in three dimensional.Therefore, the information needs to be mapped from 3D to 2D.The Finite Area Method (FAM) is used to solve this problem.
The time step of flow simulation is small.If the morphology is solved with the same time step, the bed elevation will change very slightly.Therefore, morphology evolution needs much computational time.In order to enhance the modeling speed, the time marching scheme is used.The time step of the morphology model is 10 times that of the flow simulation.
Figure 1 shows the flow chart of scour model constructed in present paper, and the details can be summarized as follows: (A), N-S Equations are solved, which are used to obtain the flow field information, and bed shear stress, suspended load concentration and bed load transport rate are calculated.
(B), the erosion rate and deposition rate of suspended load are computed, and the Exner equation is solved with FAM to obtain the change of the bed surface.
(C), the Laplace equation is solved to obtain the movement of mesh vertices, and then mesh points are moved with assigned values.
(D), the sand slide model is applied to adjust the bed surface when the bed angle is larger than the reposing angle.
(E), the mesh points are moved to new location, which is corrected by step (D), until all angles are smaller than the reposing angle.
(F), go to the next time step and repeat steps (A)-(E) until the specified time has been reached.
Water 2017, 9, 642 6 of 23 The time step of flow simulation is small.If the morphology is solved with the same time step, the bed elevation will change very slightly.Therefore, morphology evolution needs much computational time.In order to enhance the modeling speed, the time marching scheme is used.The time step of the morphology model is 10 times that of the flow simulation.
Figure 1 shows the flow chart of scour model constructed in present paper, and the details can be summarized as follows: (A), N-S Equations are solved, which are used to obtain the flow field information, and bed shear stress, suspended load concentration and bed load transport rate are calculated.
(B), the erosion rate and deposition rate of suspended load are computed, and the Exner equation is solved with FAM to obtain the change of the bed surface.
(C), the Laplace equation is solved to obtain the movement of mesh vertices, and then mesh points are moved with assigned values.
(D), the sand slide model is applied to adjust the bed surface when the bed angle is larger than the reposing angle.
(E), the mesh points are moved to new location, which is corrected by step (D), until all angles are smaller than the reposing angle.
(F), go to the next time step and repeat steps (A)-(E) until the specified time has been reached.

Bed Shear Stress
Bed shear stress is one basic essential factor that affects the precision of scour models directly.It is difficult to obtain reasonable results when the bed shear stress is uncorrected.In order to validate the bed shear stress model, an oscillation flow experiment with a high Reynolds number [28] was simulated.Figure 2 shows the schematic profile of the experiment.It can be seen that the workplace shape is a rectangle, 10 m in length, 0.28 m in height and 0.39 m in width, and the oscillation flow is controlled by a pneumatic system.More details can be found in Jensen et al. (1989) [28].

Bed Shear Stress
Bed shear stress is one basic essential factor that affects the precision of scour models directly.It is difficult to obtain reasonable results when the bed shear stress is uncorrected.In order to validate the bed shear stress model, an oscillation flow experiment with a high Reynolds number [28] was simulated.Figure 2 shows the schematic profile of the experiment.It can be seen that the workplace shape is a rectangle, 10 m in length, 0.28 m in height and 0.39 m in width, and the oscillation flow is controlled by a pneumatic system.More details can be found in Jensen et al. (1989) [28].In the scour model, the bed friction velocity was no longer achieved based on the parabolic distribution assumption.By contrast, the bed shear stress τb should be computed before the friction velocity u* is calculated.
For the Newtonian fluid under laminar conditions, the bed shear stress can be calculated according to the definition of the stress in Newtonian fluid.
While for the Reynolds-averaged Navier-Stokes equations (RANS) turbulence condition, the Reynolds stress should be considered in the shear stress model, as illustrated by Equation (23).The viscosity coefficient μ is 1 × 10 −6 , and the dynamic viscosity coefficient ν is derived from the computed result of the turbulence model: To simply the simulation, only the working area was simulated, and a two dimensional flow model was constructed, which is 10 m in a horizontal direction and 0.28 m in a vertical direction.To ensure the accuracy of the result, a fine mesh 5 mm in the horizontal and 1 mm in the vertical direction was generated.The inflow velocity at any specific moment can be computed by U = U0 × sin (w × t + Q); U0 is the maximum inlet flow velocity (0.32 m/s), ω is the frequency oscillatory flow that can be calculated by ω = 2π/T, T is the period of oscillation flow, i.e., 2.7 s, and Q is the initial phase of model, which is 0 in this model.The Reynolds number is 37,000, which is in the turbulence regime; therefore, the standard k-e model is used.The efficiency roughness height in the model is 2.5 times the sand diameter, which is around 0.0005 m.The bed friction velocity at x = 5 m is collected for comparing with the measured data, and the comparison result is shown in Figure 3.In the scour model, the bed friction velocity was no longer achieved based on the parabolic distribution assumption.By contrast, the bed shear stress τ b should be computed before the friction velocity u * is calculated.
For the Newtonian fluid under laminar conditions, the bed shear stress can be calculated according to the definition of the stress in Newtonian fluid.
While for the Reynolds-averaged Navier-Stokes equations (RANS) turbulence condition, the Reynolds stress should be considered in the shear stress model, as illustrated by Equation (23).The viscosity coefficient µ is 1 × 10 −6 , and the dynamic viscosity coefficient ν is derived from the computed result of the turbulence model: To simply the simulation, only the working area was simulated, and a two dimensional flow model was constructed, which is 10 m in a horizontal direction and 0.28 m in a vertical direction.To ensure the accuracy of the result, a fine mesh 5 mm in the horizontal and 1 mm in the vertical direction was generated.The inflow velocity at any specific moment can be computed by U = U 0 × sin (w × t + Q); U 0 is the maximum inlet flow velocity (0.32 m/s), ω is the frequency oscillatory flow that can be calculated by ω = 2π/T, T is the period of oscillation flow, i.e., 2.7 s, and Q is the initial phase of model, which is 0 in this model.The Reynolds number is 37,000, which is in the turbulence regime; therefore, the standard k-e model is used.The efficiency roughness height in the model is 2.5 times the sand diameter, which is around 0.0005 m.The bed friction velocity at x = 5 m is collected for comparing with the measured data, and the comparison result is shown in Figure 3.

Zero Entrainment Experiment
The experiment conducted by Wang and Ribeerink et al. [29] demonstrated a case where sands fall from a sand feed to water and is then transported with flow.The main purpose of simulating this experiment is to test the convection-diffusion equation.In the experiment, the bed was made up of a perforated plate and there were some compartments under the plate.Therefore, once the sands fall to the bed, re-suspension is difficult.The schematic profile of the experiment can be seen in Figure 4.The water depth of the channel is 0.215 m and the inlet mean velocity is 0.56 m/s, and the median partical diameter is 0.25 m.To save computational time, the inlet boundary location was set at the beginning of the perforated plate(x = 0).The measured suspended concentration profile at x = 0 was used as the inlet boundary condition, while other hydrodynamic parameters like flow velocity profile, pressure, eddy viscosity and so on were obtained from the last simulation.This means that the model should be simulated twice.First, a hydrodynamic model was computed solely, and the average velocity and pressure were used as the inlet boundary condition.Then, the hydrodynamic parameters at a specific place, where the turbulence was full developed, were abstracted as the boundary condition of the following simulation.The bed roughness height in models is 0.005 m, and the sediment falling velocity has a uniform value of 0.022 m/s.

Zero Entrainment Experiment
The experiment conducted by Wang and Ribeerink et al. [29] demonstrated a case where sands fall from a sand feed to water and is then transported with flow.The main purpose of simulating this experiment is to test the convection-diffusion equation.In the experiment, the bed was made up of a perforated plate and there were some compartments under the plate.Therefore, once the sands fall to the bed, re-suspension is difficult.The schematic profile of the experiment can be seen in Figure 4.The water depth of the channel is 0.215 m and the inlet mean velocity is 0.56 m/s, and the median partical diameter is 0.25 m.

Zero Entrainment Experiment
The experiment conducted by Wang and Ribeerink et al. [29] demonstrated a case where sands fall from a sand feed to water and is then transported with flow.The main purpose of simulating this experiment is to test the convection-diffusion equation.In the experiment, the bed was made up of a perforated plate and there were some compartments under the plate.Therefore, once the sands fall to the bed, re-suspension is difficult.The schematic profile of the experiment can be seen in Figure 4.The water depth of the channel is 0.215 m and the inlet mean velocity is 0.56 m/s, and the median partical diameter is 0.25 m.To save computational time, the inlet boundary location was set at the beginning of the perforated plate(x = 0).The measured suspended concentration profile at x = 0 was used as the inlet boundary condition, while other hydrodynamic parameters like flow velocity profile, pressure, eddy viscosity and so on were obtained from the last simulation.This means that the model should be simulated twice.First, a hydrodynamic model was computed solely, and the average velocity and pressure were used as the inlet boundary condition.Then, the hydrodynamic parameters at a specific place, where the turbulence was full developed, were abstracted as the boundary condition of the following simulation.The bed roughness height in models is 0.005 m, and the sediment falling velocity has a uniform value of 0.022 m/s.To save computational time, the inlet boundary location was set at the beginning of the perforated plate(x = 0).The measured suspended concentration profile at x = 0 was used as the inlet boundary condition, while other hydrodynamic parameters like flow velocity profile, pressure, eddy viscosity and so on were obtained from the last simulation.This means that the model should be simulated twice.First, a hydrodynamic model was computed solely, and the average velocity and pressure were used as the inlet boundary condition.Then, the hydrodynamic parameters at a specific place, where the turbulence was full developed, were abstracted as the boundary condition of the following simulation.The bed roughness height in models is 0.005 m, and the sediment falling velocity has a uniform value of 0.022 m/s.
The comparison between modeled and measured data is shown in Figure 5.It can be seen that the sediment concentration changes with the location of the reference point.It is obvious that the maximum sediment concentration value in the reference vertical section declines with the increase of the x coordinate, and a slight change occurs in the parabolic distribution of the concentration profile.Overall, the modeled result and measured data match very well; though there is still some error between them, the maximum concentration near bed is overestimated slightly.
Water 2017, 9, 642 9 of 23 The comparison between modeled and measured data is shown in Figure 5.It can be seen that the sediment concentration changes with the location of the reference point.It is obvious that the maximum sediment concentration value in the reference vertical section declines with the increase of the x coordinate, and a slight change occurs in the parabolic distribution of the concentration profile.Overall, the modeled result and measured data match very well; though there is still some error between them, the maximum concentration near bed is overestimated slightly.

Net Entrainment Experiment
The suspended load model is validated by the van Rijn's [30] net erosion experiment.The schematic view of the experiment is shown in Figure 6.It can be seen that the total experiment area is divided into two parts: a rigid bed and a loose sand bed.For the rigid bed area, the water is clean and there is no sediment from the bed.As for the loose bed area, the suspended load will be entrained from the bed.The flow inlet is seated in the rigid bed area.The water depth of the experiment is 0.25 m, the inlet velocity is 0.67 m/s, the sand particle diameter in the loose bed is 0.2 mm, the corresponding setting velocity is 0.022 m/s, and the critical shields number is 0.042.The sediment concentration of the bottom boundary plays a crucial role in the suspended load model, which affects the result of convection and diffusion equation directly.The concentration at the reference height above the bed can be calculated by Equation (15).In this model, the reference height in this simulation is 0.01 m, and the bed roughness height is 2.5d50, which is approximately 0.005 m.Therefore, the concentration at reference height Cb is 4.6 kg/m 3 , which is larger than Liang's [19] result.In OpenFOAM, the basic cell is the mesh volume cell and variations like velocity, pressure, concentration and so on, are stored at the center of the volume.The mesh size in the vertical direction near the bed is 0.01 m, which means the concentration at the mesh center is Cb.In this simulation, the scour model is not included.As for the flow model, the k-epsilon closure model is used to simulate the turbulence and the VOF method is used to capture the water surface.For the turbulence to be

Net Entrainment Experiment
The suspended load model is validated by the van Rijn's [30] net erosion experiment.The schematic view of the experiment is shown in Figure 6.It can be seen that the total experiment area is divided into two parts: a rigid bed and a loose sand bed.For the rigid bed area, the water is clean and there is no sediment from the bed.As for the loose bed area, the suspended load will be entrained from the bed.The flow inlet is seated in the rigid bed area.The water depth of the experiment is 0.25 m, the inlet velocity is 0.67 m/s, the sand particle diameter in the loose bed is 0.2 mm, the corresponding setting velocity is 0.022 m/s, and the critical shields number is 0.042.

Net Entrainment Experiment
The suspended load model is validated by the van Rijn's [30] net erosion experiment.The schematic view of the experiment is shown in Figure 6.It can be seen that the total experiment area is divided into two parts: a rigid bed and a loose sand bed.For the rigid bed area, the water is clean and there is no sediment from the bed.As for the loose bed area, the suspended load will be entrained from the bed.The flow inlet is seated in the rigid bed area.The water depth of the experiment is 0.25 m, the inlet velocity is 0.67 m/s, the sand particle diameter in the loose bed is 0.2 mm, the corresponding setting velocity is 0.022 m/s, and the critical shields number is 0.042.The sediment concentration of the bottom boundary plays a crucial role in the suspended load model, which affects the result of convection and diffusion equation directly.The concentration at the reference height above the bed can be calculated by Equation (15).In this model, the reference height in this simulation is 0.01 m, and the bed roughness height is 2.5d50, which is approximately 0.005 m.Therefore, the concentration at reference height Cb is 4.6 kg/m 3 , which is larger than Liang's [19] result.In OpenFOAM, the basic cell is the mesh volume cell and variations like velocity, pressure, concentration and so on, are stored at the center of the volume.The mesh size in the vertical direction near the bed is 0.01 m, which means the concentration at the mesh center is Cb.In this simulation, the scour model is not included.As for the flow model, the k-epsilon closure model is used to simulate the turbulence and the VOF method is used to capture the water surface.For the turbulence to be The sediment concentration of the bottom boundary plays a crucial role in the suspended load model, which affects the result of convection and diffusion equation directly.The concentration at the reference height above the bed can be calculated by Equation (15).In this model, the reference height in this simulation is 0.01 m, and the bed roughness height is 2.5d 50 , which is approximately 0.005 m.Therefore, the concentration at reference height C b is 4.6 kg/m 3 , which is larger than Liang's [19] result.In OpenFOAM, the basic cell is the mesh volume cell and variations like velocity, pressure, concentration and so on, are stored at the center of the volume.The mesh size in the vertical direction Water 2017, 9, 642 11 of 23 near the bed is 0.01 m, which means the concentration at the mesh center is C b .In this simulation, the scour model is not included.As for the flow model, the k-epsilon closure model is used to simulate the turbulence and the VOF method is used to capture the water surface.For the turbulence to be fully developed, the length of the rigid bed and loose bed is 20H and 60H, respectively.The suspended concentration of inlet is zero for the clear water condition.The concentration profile at x = 4H, 10H, 20H and 40H is abstracted to validate the simulated result, and x = 0 is the location of the beginning of the loose bed. Figure 7 shows that the sediment concentration profile is parabolic.It can be seen that the sediment concentration decreases rapidly with the height distance from the bed surface, and it becomes a parabolic profiles.In a word, the modeled result matches well with the experimental data.
Water 2017, 9, 642 11 of 23 fully developed, the length of the rigid bed and loose bed is 20H and 60H, respectively.The suspended concentration of inlet is zero for the clear water condition.The concentration profile at x = 4H, 10H, 20H and 40H is abstracted to validate the simulated result, and x = 0 is the location of the beginning of the loose bed. Figure 7 shows that the sediment concentration profile is parabolic.It can be seen that the sediment concentration decreases rapidly with the height distance from the bed surface, and it becomes a parabolic profiles.In a word, the modeled result matches well with the experimental data.

Local Scour under a Pipeline
Mao [1] conducted a laboratory experiment to study the local scour under a pipeline under steady flow.Here, the scour model will be calibrated by the clear water case, where the average income velocity is 0.35 m/s, the depth of the experiment is 0.35 m, the pipeline diameter is D = 0.1 mm, and the median particle diameter of sediment d = 0.36 mm.
The meshing process usually affects the result of numerical simulation.In present work, Bjarke's [13] meshing strategy is adopted for scour simulation, and he pointed out that increasing the mesh resolution near the bed and pipeline can solve the problem of mesh distortion as the scour hole broadens and deepens.Therefore, a gradually changed mesh is used, the minimum mesh size near the bottom is 0.5 mm and that near the pipeline is 0.3 mm.To prevent the change of mesh topology,

Local Scour under a Pipeline
Mao [1] conducted a laboratory experiment to study the local scour under a pipeline under steady flow.Here, the scour model will be calibrated by the clear water case, where the average income velocity is 0.35 m/s, the depth of the experiment is 0.35 m, the pipeline diameter is D = 0.1 mm, and the median particle diameter of sediment d = 0.36 mm.
The meshing process usually affects the result of numerical simulation.In present work, Bjarke's [13] meshing strategy is adopted for scour simulation, and he pointed out that increasing the mesh resolution near the bed and pipeline can solve the problem of mesh distortion as the scour hole broadens and deepens.Therefore, a gradually changed mesh is used, the minimum mesh size near the bottom is 0.5 mm and that near the pipeline is 0.3 mm.To prevent the change of mesh topology, a sinusoidal hole with 0.1D is set under the pipeline.To ensure the turbulence and wake vortex are developed fully, the horizontal span of the numerical model ranges from −25D to 25D.The mesh condition near the pipeline can be seen in Figure 8. a sinusoidal hole with 0.1D is set under the pipeline.To ensure the turbulence and wake vortex are developed fully, the horizontal span of the numerical model ranges from −25D to 25D.The mesh condition near the pipeline can be seen in Figure 8.Before the simulation, we performed a mesh convergence study in which we simulated five different mesh size cases of 1 cm, 5 mm, 1 mm, 0.5 mm and 0.1 mm.We then abstracted the friction velocity of the bed at x = 0D and x = 2D.D is the pipeline diameter, and the results are shown in Table 1.We found that the modeled results resembled more closely the measurements when the mesh size was finer, but there was no significant improvement with mesh sizes finer than 0.5 mm.Thus, in this work, we adopted a mesh size of 0.5 mm in the next simulation.To obtain the income boundary condition of the flow model, a separated long flume with the same mesh size as the scour model was simulated first.The bed efficiency roughness height ks = 2.5D and the RNG k-ε turbulence model were used to simulate the turbulence characteristic.As for the sediment transport model, the critical shields number for a flat bed is 0.048; the reposing angle is 34 degrees, and the sediment concentration at the inlet boundary is 0 due to the clear water condition.
The scour profile at 10 min, 30 min, 200 min and 370 min was selected to verify the scour model and the moving mesh method.The comparison between measured and modeled data is shown in Figure 9.As illustrated in Figure 9, the scour hole deepens and extends dramatically during the first 30 min, and the hole depth almost keeps steady during the period from 200 min to 360 min.The maximum equilibrium scour depth is about 0.54D, and the maximum deposition height is around Before the simulation, we performed a mesh convergence study in which we simulated five different mesh size cases of 1 cm, 5 mm, 1 mm, 0.5 mm and 0.1 mm.We then abstracted the friction velocity of the bed at x = 0D and x = 2D.D is the pipeline diameter, and the results are shown in Table 1.We found that the modeled results resembled more closely the measurements when the mesh size was finer, but there was no significant improvement with mesh sizes finer than 0.5 mm.Thus, in this work, we adopted a mesh size of 0.5 mm in the next simulation.To obtain the income boundary condition of the flow model, a separated long flume with the same mesh size as the scour model was simulated first.The bed efficiency roughness height ks = 2.5D and the RNG k-ε turbulence model were used to simulate the turbulence characteristic.As for the sediment transport model, the critical shields number for a flat bed is 0.048; the reposing angle is 34 degrees, and the sediment concentration at the inlet boundary is 0 due to the clear water condition.
Water 2017, 9, 642 13 of 23 The scour profile at 10 min, 30 min, 200 min and 370 min was selected to verify the scour model and the moving mesh method.The comparison between measured and modeled data is shown in Figure 9.As illustrated in Figure 9, the scour hole deepens and extends dramatically during the first 30 min, and the hole depth almost keeps steady during the period from 200 min to 360 min.The maximum equilibrium scour depth is about 0.54D, and the maximum deposition height is around 0.48D.Though there is a slight difference in the dune's profile between modeled and measured results, the precision of the modelling result can satisfy the following study's need.The velocity vector and scour profiles at 0.1 min, 0.5 min, 1 min and 5 min are abstracted to study the developing process of the scour hole, as shown in Figure 10.It is clearly seen that the hole under the pipeline deepens dramatically and a small dune appears in the downstream at the beginning of the simulation.Then, the dune moves toward the downstream slowly; its height increases gradually, and at the same time the scour hole depth increases gradually.The velocity vector and scour profiles at 0.1 min, 0.5 min, 1 min and 5 min are abstracted to study the developing process of the scour hole, as shown in Figure 10.It is clearly seen that the hole under the pipeline deepens dramatically and a small dune appears in the downstream at the beginning of the simulation.Then, the dune moves toward the downstream slowly; its height increases gradually, and at the same time the scour hole depth increases gradually.The velocity vector and scour profiles at 0.1 min, 0.5 min, 1 min and 5 min are abstracted to study the developing process of the scour hole, as shown in Figure 10.It is clearly seen that the hole under the pipeline deepens dramatically and a small dune appears in the downstream at the beginning of the simulation.Then, the dune moves toward the downstream slowly; its height increases gradually, and at the same time the scour hole depth increases gradually.

The Channel without a Pipeline
A separated long channel without pipelines is simulated to study the difference of bed shear stress results under free surface and rigid lid conditions.The model parameters are shown in Table 2, where Froude numbers 0.1, 0.2 and 0.3 are considered.Every Froude number case includes seven different Reynolds number conditions to prevent the effect of individuals.Figure 11 shows the bed shear stress computed by rigid lid approximation and the free surface model.The x axis represents the bed shear stress value calculated by the rigid lid model, and the y axis is the value computed under the free surface condition.The linear coefficients of the result are calculated through the minimum square method.It can be observed that the linear coefficients of the fitted line decline with increasing Froude number.When the Froude number reaches 0.3, the slope coefficient of the fitting line is 0.94, which means the difference of the result of the free surface and rigid lid model is approximately 6%.Therefore, in the simulation of the channel without obstructions, even if the Froude number is larger than 0.2, the results computed by the rigid lid assumption are also acceptable-smaller than 10%.

The Channel without a Pipeline
A separated long channel without pipelines is simulated to study the difference of bed shear stress results under free surface and rigid lid conditions.The model parameters are shown in Table 2, where Froude numbers 0.1, 0.2 and 0.3 are considered.Every Froude number case includes seven different Reynolds number conditions to prevent the effect of individuals.Figure 11 shows the bed shear stress computed by rigid lid approximation and the free surface model.The x axis represents the bed shear stress value calculated by the rigid lid model, and the y axis is the value computed under the free surface condition.The linear coefficients of the result are calculated through the minimum square method.It can be observed that the linear coefficients of the fitted line decline with increasing Froude number.When the Froude number reaches 0.3, the slope coefficient of the fitting line is 0.94, which means the difference of the result of the free surface and rigid lid model is approximately 6%.Therefore, in the simulation of the channel without obstructions, even if the Froude number is larger than 0.2, the results computed by the rigid lid assumption are also acceptable-smaller than 10%.

Validation of the Free Surface
An experimental case on the free surface flow around a half-circular cylinder is simulated first, which was conducted by Frobes and Scheawartz [31].The computational parameters of the model are the same as those used in Zhao's model [16].The ratio of the water depth to the cylinder diameter is 1.25, the Reynolds number defined by Re = U × D/v is 1.92 × 10 5 , and the Froude number Fr = 0.373.Figure 12 shows the comparison between the modeled and measured free surface profiles.The change in free surface elevation above the half cylinder is well predicted by the numerical model.

Validation of the Free Surface
An experimental case on the free surface flow around a half-circular cylinder is simulated first, which was conducted by Frobes and Scheawartz [31].The computational parameters of the model are the same as those used in Zhao's model [16].The ratio of the water depth to the cylinder diameter is 1.25, the Reynolds number defined by Re = U × D/v is 1.92 × 10 5 , and the Froude number Fr = 0.373.Figure 12 shows the comparison between the modeled and measured free surface profiles.The change in free surface elevation above the half cylinder is well predicted by the numerical model.

The Surface Elevation and Bed Shear Stress
Frobes (1982) studied the free-surface profile as the fluid flow over a semicircular obstacle on the bottom [31], and two parameters, Froude number Fr and dimensionless circle radius, were considered in his study.It was found that the ratio of the diameter of the circle to the depth also has a significant effect on the water level, and a train of nonlinear stokes wave is observed in the downstream when the dimensionless circle radius and Froude number reach certain values.
To study the impact of the ratio of the depth to the pipeline diameter on the surface elevation, four cases with the same Froude number (Fr = 0.2 and Fr = 0.3) and different dimensionless depths (H/D) are simulated, and the pipeline diameter in those simulations are kept constant, i.e., approximately 0.1 m.The ratios of the depth to the pipeline diameter are 2, 3, 4 and 5, respectively.The incoming flow velocity of those cases can be calculated based on the Froude number formula.
Figure 13 illustrates the distribution of the free surface of the Froude number Fr = 0.2 and Fr = 0.3.The x axis is the dimensionless distance, and the y axis is the ratio of surface elevation to the diameter of the pipeline.The free surface experiences a slight drop in the area above the pipeline, and the difference between the surface elevation upstream and downstream decreases with an increase in the value of dimensionless depth.For the condition Fr = 0.2, the maximum drop of the surface elevation is 0.07D as the dimensionless depth H/D = 2.As for the condition of Fr = 0.3, this value reaches 0.22D.Besides, when Fr = 0.3 and H/D are smaller than or equal to 3, a regular wave will appear in the downstream region.

The Surface Elevation and Bed Shear Stress
Frobes (1982) studied the free-surface profile as the fluid flow over a semicircular obstacle on the bottom [31], and two parameters, Froude number Fr and dimensionless circle radius, were considered in his study.It was found that the ratio of the diameter of the circle to the depth also has a significant effect on the water level, and a train of nonlinear stokes wave is observed in the downstream when the dimensionless circle radius and Froude number reach certain values.
To study the impact of the ratio of the depth to the pipeline diameter on the surface elevation, four cases with the same Froude number (Fr = 0.2 and Fr = 0.3) and different dimensionless depths (H/D) are simulated, and the pipeline diameter in those simulations are kept constant, i.e., approximately 0.1 m.The ratios of the depth to the pipeline diameter are 2, 3, 4 and 5, respectively.The incoming flow velocity of those cases can be calculated based on the Froude number formula.
Figure 13 illustrates the distribution of the free surface of the Froude number Fr = 0.2 and Fr = 0.3.The x axis is the dimensionless distance, and the y axis is the ratio of surface elevation to the diameter of the pipeline.The free surface experiences a slight drop in the area above the pipeline, and the difference between the surface elevation upstream and downstream decreases with an increase in the value of dimensionless depth.For the condition Fr = 0.2, the maximum drop of the surface elevation is 0.07D as the dimensionless depth H/D = 2.As for the condition of Fr = 0.3, this value reaches 0.22D.Besides, when Fr = 0.3 and H/D are smaller than or equal to 3, a regular wave will appear in the downstream region.
Compared with other conditions, the drop of surface elevation between upstream and downstream regions is the largest as the dimensionless depth equals 2. Therefore, the bed shear stress of Fr = 0.2, H/D = 2 and Fr = 0.3, H/D = 2 under free surface and rigid lid assumption are calculated.Figure 14 shows the distribution of the dimensionless friction velocity under the pipeline before the scour model works.It can be seen that the bottom friction velocity calculated by free surface model is larger than that computed under the rigid lid assumption.This is because the surface elevation above the pipeline will drop if the free surface is considered in the simulation, which means that the area of that section will decrease, leading to the average velocity of that section increasing.The differences of modeled results mainly appear under the pipeline and downstream.When Fr = 0.3, the maximum difference of dimensionless friction velocity between the rigid lid assumption and free surface occurs at x/D = 0, i.e., approximately 12%.
The x axis is the dimensionless distance, and the y axis is the ratio of surface elevation to the diameter of the pipeline.The free surface experiences a slight drop in the area above the pipeline, and the difference between the surface elevation upstream and downstream decreases with an increase in the value of dimensionless depth.For the condition Fr = 0.2, the maximum drop of the surface elevation is 0.07D as the dimensionless depth H/D = 2.As for the condition of Fr = 0.3, this value reaches 0.22D.Besides, when Fr = 0.3 and H/D are smaller than or equal to 3, a regular wave will appear in the downstream region.Compared with other conditions, the drop of surface elevation between upstream and downstream regions is the largest as the dimensionless depth equals 2. Therefore, the bed shear stress of Fr = 0.2, H/D = 2 and Fr = 0.3, H/D = 2 under free surface and rigid lid assumption are calculated.Figure 14 shows the distribution of the dimensionless friction velocity under the pipeline before the scour model works.It can be seen that the bottom friction velocity calculated by free surface model is larger than that computed under the rigid lid assumption.This is because the surface elevation above the pipeline will drop if the free surface is considered in the simulation, which means that the area of that section will decrease, leading to the average velocity of that section increasing.The differences of modeled results mainly appear under the pipeline and downstream.When Fr = 0.3, the maximum difference of dimensionless friction velocity between the rigid lid assumption and free surface occurs at x/D = 0, i.e., approximately 12%.

The Scour Hole
According to the study above, the surface elevation near the pipeline will experience a slight drop, which also affects the precision of the scour simulation.Here, the impact of the rigid lid model on the local scour will be studied based on Mao's experiment, where the dimensionless depth is the same (H/D = 3.5) and the Fr = 0.187, Fr = 0.214 and Fr = 0.267.The mesh adopted in this section is the same as that used in the scour validation case.
Figures 15-17 show the scour hole profile modeled at different times of the three cases.Through comparison with the result computed by the free surface, the hole depth with the rigid lid approximation was underestimated, and there was also a slight difference in the scour hole profile.The main difference occurs at the location under the pipeline and in the downstream area.A higher Froude number will enlarge the errors of local scour simulation under the rigid lid assumption.When Fr = 0.267 and the dimensionless depth H/D = 3.5, the maximum difference between the rigid lid model and measured data is approximately 14%.
On the other hand, both the results computed with the rigid lid assumption and a free surface match very well with the measured data.Although the scour profile is more accurate, the Central

The Scour Hole
According to the study above, the surface elevation near the pipeline will experience a slight drop, which also affects the precision of the scour simulation.Here, the impact of the rigid lid model on the local scour will be studied based on Mao's experiment, where the dimensionless depth is the same (H/D = 3.5) and the Fr = 0.187, Fr = 0.214 and Fr = 0.267.The mesh adopted in this section is the same as that used in the scour validation case.
Figures 15-17 show the scour hole profile modeled at different times of the three cases.Through comparison with the result computed by the free surface, the hole depth with the rigid lid approximation was underestimated, and there was also a slight difference in the scour hole profile.The main difference occurs at the location under the pipeline and in the downstream area.A higher Froude number will enlarge the errors of local scour simulation under the rigid lid assumption.When Fr = 0.267 and the dimensionless depth H/D = 3.5, the maximum difference between the rigid lid model and measured data is approximately 14%.On the other hand, both the results computed with the rigid lid assumption and a free surface match very well with the measured data.Although the scour profile is more accurate, the Central Processing Unit (CPU) time of the free surface model is almost 1.5 times that for the model with the rigid lid approximation. (c)

Discussion
A numerical scour model has been developed to study the effect of rigid lid assumption on scour under a pipeline under steady flow when the Froude number is close to 0.2.In the scour model, a new method to calculate bed shear stress was adopted, and a three-dimensional sand slide was constructed to smooth the mesh.A time marching scheme was used to accelerate the simulation, FAM was used to solve the Exner Equation, and a moving mesh method was used to capture the change of the bed surface.The effect of a free surface was considered in the local scour model.All of this work was well validated with experimental data.Last but not least, the scour under the steady flow under the free surface condition was studied with our model, and the study results show that: (1) Only using the Froude number to distinguish whether the rigid lid assumption is suitable for scour simulation is unreasonable, and the dimensionless depth should be considered.For the condition of the channel without a pipeline, even if the Froude number reaches 0.3, the difference of bed shear stress computed by the rigid lid assumption and free surface is only 6%.
(2) For the channel with a pipeline, the dimensionless depth (H/D) has a significant impact on the surface elevation simulation.The difference of surface elevation between the upstream and downstream increases with decreasing dimensionless depth under a fixed Froude number.When Froude number Fr = 0.3 and the dimensionless depth is smaller than 3, a train wave will occur in the downstream

Discussion
A numerical scour model has been developed to study the effect of rigid lid assumption on scour under a pipeline under steady flow when the Froude number is close to 0.2.In the scour model, a new method to calculate bed shear stress was adopted, and a three-dimensional sand slide was constructed to smooth the mesh.A time marching scheme was used to accelerate the simulation, FAM was used to solve the Exner Equation, and a moving mesh method was used to capture the change of the bed surface.The effect of a free surface was considered in the local scour model.All of this work was well validated with experimental data.Last but not least, the scour under the steady flow under the free surface condition was studied with our model, and the study results show that: (1) Only using the Froude number to distinguish whether the rigid lid assumption is suitable for scour simulation is unreasonable, and the dimensionless depth should be considered.For the condition of the channel without a pipeline, even if the Froude number reaches 0.3, the difference of bed shear stress computed by the rigid lid assumption and free surface is only 6%.
(2) For the channel with a pipeline, the dimensionless depth (H/D) has a significant impact on the surface elevation simulation.The difference of surface elevation between the upstream and downstream increases with decreasing dimensionless depth under a fixed Froude number.When Froude number Fr = 0.3 and the dimensionless depth is smaller than 3, a train wave will occur in the downstream region of the pipeline.
(3) Comparing with the result computed with the rigid lid approximation, the result where the free surface is considered is closer to the experimental data, while the cost of CPU time is also expensive.

Figure 1 .
Figure 1.Flow chart of the local scour model.

Figure 1 .
Figure 1.Flow chart of the local scour model.

Figure 3 .
Figure 3.The distribution of friction velocity of modeled values in different phases.

Figure 4 .
Figure 4. Schematic description of Zero Erosion experiment.

Figure 3 .
Figure 3.The distribution of friction velocity of modeled values in different phases.

Figure 3 .
Figure 3.The distribution of friction velocity of modeled values in different phases.

Figure 4 .
Figure 4. Schematic description of Zero Erosion experiment.

Figure 4 .
Figure 4. Schematic description of Zero Erosion experiment.

Figure 6 .
Figure 6.Schematic view of the entrainment experiment.

Figure 6 .
Figure 6.Schematic view of the entrainment experiment.

Figure 6 .
Figure 6.Schematic view of the entrainment experiment.

Figure 8 .
Figure 8.The mesh distribution around the pipeline in the scour test.

Figure 8 .
Figure 8.The mesh distribution around the pipeline in the scour test.

Figure 11 .
Figure 11.The correlation of the free surface modeled data and the rigid lid model (the circle is the modeled data, the line is the fitted line); (a) Fr = 0.1; (b) Fr = 0.2; (c) Fr = 0.3.

Figure 12 .
Figure 12.The calculated free surface profile result of the half-circular cylinder case.

Figure 12 .
Figure 12.The calculated free surface profile result of the half-circular cylinder case.

Figure 14 .
Figure 14.The distribution of dimensionless friction velocity under the pipeline at the initial moment.

Figure 14 .
Figure 14.The distribution of dimensionless friction velocity under the pipeline at the initial moment.

Figure 15 .
Figure 15.The scour profile of the free surface model and rigid lid assumption at different times in the case of Fr = 0.187, (a) T = 600 s; (b) T = 1800 s; (c) T = 12,000 s.Figure 15.The scour profile of the free surface model and rigid lid assumption at different times in the case of Fr = 0.187, (a) T = 600 s; (b) T = 1800 s; (c) T = 12,000 s.

Figure 15 .Figure 16 .Figure 16 .
Figure 15.The scour profile of the free surface model and rigid lid assumption at different times in the case of Fr = 0.187, (a) T = 600 s; (b) T = 1800 s; (c) T = 12,000 s.Figure 15.The scour profile of the free surface model and rigid lid assumption at different times in the case of Fr = 0.187, (a) T = 600 s; (b) T = 1800 s; (c) T = 12,000 s.

Figure 16 .Figure 17 .
Figure 16.The scour profile of the free surface model and rigid lid assumption at different times in the case of Fr = 0.214, (a) T = 390 s; (b) T = 1890 s; (c) T = 4860 s.

Figure 17 .
Figure 17.The scour profile of the free surface model and rigid lid assumption at different times in the case of Fr = 0.267, (a) T = 90 s; (b) T = 300 s; (c) T = 1800 s.

Table 1 .
Simulated bottom velocity response to the mesh size.

Table 1 .
Simulated bottom velocity response to the mesh size.

Table 2 .
The parameters of the study for hydrodynamic tests.

Table 2 .
The parameters of the study for hydrodynamic tests.