Numerical prediction of bed-load and surface deformation on a granular bed sheared by a turbulent boundary-layer

TECHNICAL PAPERS

FLUID MECHANICS

Numerical prediction of bed-load and surface deformation on a granular bed sheared by a turbulent boundary-layer

Erick de Moraes Franklin

Itajubá Federal University - Unifei. Institute of Mechanical Engineering Itajubá, MG, Brazil. [e-mail]

ABSTRACT

The entrainment of solid particles by a fluid flow is frequently found in nature and in industry. A better knowledge of this is of importance to understand the physical nature of the phenomenon and to improve industrial processes. Bed-load occurs if the shear stresses exerted by the fluid on the granular bed are bounded to some limits: a mobile granular layer takes place over the fixed part of the bed. When the fluid is a liquid, the thickness of this mobile layer is a few grain-diameters. Under these conditions, an initially flat granular bed may be unstable, giving rise to ripples and dunes. Some examples are the dunes seen in deserts, but also the ripples appearing in petroleum pipelines conveying sand. This article presents a mathematical model for the prediction of bed-load under a turbulent boundary-layer of a liquid, and its numerical implementation. The model, kept as simple as possible, focuses on the mobile layer of the granular bed, reducing then the computation domain. It is able to capture the pertinent physics involved, such as the growing of instabilities on the granular bed, so that, if employed together with the stability analysis of Franklin (2010), it selects and predicts the evolution of most unstable modes.

Keywords: mathematical modeling, numerical simulations, turbulent boundary-layer, bedload, instability

Introduction

The granular media is of great importance in our quotidian, and their transport by a fluid flow is frequently found in nature and in industry. It is present in the erosion of coastal regions, in the displacement of desert dunes, in sewerage systems and in hydrocarbon pipelines conveying sand, for instance. A better knowledge of this kind of transport is then of great importance to both understand nature and improve related industrial processes.

If the force due to the shear stresses exerted by the fluid flow on the granular bed is small compared to the weight of individual grains, but is able to move some grains, the flow cannot transport grains as a suspension. Instead, a mobile granular layer, known as bed-load, takes place over the fixed part of the bed. If the fluid is a liquid, the thickness of this mobile layer is a few grain-diameters (Bagnold, 1941; Raudkivi, 1976).

Motivated by environmental and industrial issues, many studies concerning the bed-load were conducted in the last decades. Perhaps the most comprehensive is Bagnold (1941), which discusses the mechanisms of granular entrainment, proposes semi-empirical expressions for the estimation of the flow rate of grains and discusses the effects of the granular mobility on the fluid flow (feedback effect). One of the main difficulties concerning bed-load is that the mobile granular layer is not a continuous media, so that the great number of contact points determine the behavior of the system. Dealing with this great number of contact points is not feasible and, on the other hand, a rheology for the granular matter is not known at present. As a consequence, there are many different formulations for the bed-load flow rate and there isn't a real consensus about it.

The growth of bedforms is another difficulty related to the modeling of the bed-load. These forms, initially two-dimensional, may grow and generate patterns such as ripples or dunes. In nature, one example is the aquatic ripples observed on river beds, which create a supplementary friction between the bed and the water, affecting then the water depth and being related to flood problems. In industry, examples are mostly related to closed-conduit flows conveying grains, such as sewerage systems or the hydrocarbon pipelines conveying sand. In these cases, the bedforms generate supplementary pressure loss and pressure and flow rate transients (Kuru et al., 1995; Franklin, 2008). In order to understand the growing of bedforms, many works on the stability of granular beds sheared by a fluid were made in the last decades (Kennedy, 1963; Reynolds, 1964; Engelund, 1970; Richards, 1980; Elbelrhiti et al., 2005; Claudin and Andreotti, 2006; for instance). Engelund and Fredsoe (1982) present a remarkable overview on this subject.

Recently, Franklin (2010) discussed some physical aspects of the instability, explained its mechanisms and presented a linear stability analysis. The stability analysis was made in the specific case of granular beds sheared by turbulent boundary-layers of liquids, without free surface effects. Different from previous linear stability analyses, Franklin (2010) showed that the length-scale of the initial bedforms varies with the fluid flow conditions. In order to understand the wavelength of ripples and dunes observed in nature, which usually have their wavelength predicted by linear analyses though they are clearly in a nonlinear phase, Franklin (2011) presented a nonlinear stability analysis in the same scope of Franklin (2010). The employed approach was the weakly nonlinear analysis (Landau and Lifchitz, 1994; Schmid and Henningson, 2001; Drazin and Reid, 2004; Charru, 2007), useful whenever a dominant mode can be proved to exist. Franklin (2011) showed that the bed instabilities saturate after the initial exponential growth (linear phase), i.e., they attenuate their growth rate and maintain the same wavelength.

In both Franklin (2010) and Franklin (2011) the results of the analyses were compared to experimental data concerning ripples in closed-conduit flows (which is a case where free surface effects are absent). In particular, the dependence of the bedform wavelength on the fluid flow conditions and the saturation of the bedform amplitude were confirmed by experimental results.

Bed-load numerical simulations are another approach employed in an effort to understand the instabilities appearing on granular beds. The main problem is then the modeling of the granular media: it is a discrete media for which a lagrangian description is not feasible given the large number of discrete elements. To solve this problem, it is usual to employ some semi-empirical expressions for the bed-load flow rate, such as those proposed by Bagnold (1941) and Meyer-Peter and Mueller (1948). Another difficulty is the computation of the fluid flow over a deformable bed, which perturbs in its turn the fluid flow.

Amongst the works on mathematical modeling and numerical simulations of bed-load under turbulent boundary-layers, there are Kroy et al. (2002a), Kroy et al. (2002b) and Hersen (2004). These models were all conceived for turbulent gas flows (aeolian case), and are based on four equations: (i) the shear stress caused by the fluid flow on the bed surface; (ii) a semi-empirical equation for the bed-load flow rate; (iii) an equation accounting for the relaxation between the fluid flow and the granular flow; and (iv) the mass conservation of the granular bed.

Because the semi-empirical expressions for the bed-load flow rate are based on the shear stress on the bed, Kroy et al. (2002a) and Kroy et al. (2002b) proposed the direct employment of an equation for it, instead of computing the entire fluid flow above the bed. If the deformation keeps a small aspect ratio, the shear stress caused by the fluid flow on a deformed bed can be obtained analytically by Perturbation Methods. For example, this was done in the case of turbulent boundary-layers by Jackson et al. (1975), Hunt et al. (1988) and Weng et al. (1991). Kroy et al. (2002a) and Kroy et al. (2002b) simplified the results of Jackson et al. (1975), Hunt et al. (1988) and Weng et al. (1991) and obtained an expression containing only the dominant physical effects of the shear stress perturbation. Concerning the bed-load flow rate, they employed an expression derived from the aeolian bed-load model of Sauermann et al. (2001), whose constants were adjusted from experimental data. The employed equation accounting for the relaxation between the flow rate of grains and the fluid flow was also obtained from Sauermann et al. (2001). The model was then implemented in a numerical code and some simulations of the aeolian case were performed. From different initial conditions, the model was able to show the evolution from an initial bump to dunes. However, the authors were not interested in the initial instabilities, and for that they did not try to find a most unstable wavelength, which is the one expected to give origin to dunes.

The model and the simulations presented in Hersen (2004) are based on the same equations of Kroy et al. (2002a) and Kroy et al. (2002b). The main difference in Hersen (2004) is that the work was focused on the formation of the crescentic shape barchan dunes.

This article presents a mathematical model for the transport of grains as bed-load by a turbulent boundary-layer, when the fluid is a liquid. The model is kept as simple as possible, but it is able to capture the pertinent physics involved, such as the growing of instabilities on the bed. Concerning the fluid flow, it employs an analytical equation for the shear stress on the bed given by Perturbation Methods, which limits the model to low aspect ratios of the bedforms. Another restriction is that the model is applied only to liquids. This restriction was not imposed in previous models of this kind (Kroy et al., 2002a; Kroy et al., 2002b; Hersen, 2004), although the bed-load characteristics in liquids are different from the aeolian case (Raudkivi, 1976). To compute the bed-load flow rate, a semiempirical formula is employed. Concerning the boundary condition, it is assumed the existence of a flat granular bed at the upstream end of the domain. As initial condition, the model needs an initial bedform, which is taken as an undulated granular bed of low aspect ratio. The implementation of the model in a numerical code is presented and the results of some simulations are confronted to the linear stability analysis of Franklin (2010) and to experimental data previously published.

The main objective of this paper is to propose a simple model that, if employed with the stability analysis of Franklin (2010) as initial condition, is able to predict the bed-load flow rate and the initial evolution of granular beds submitted to turbulent boundarylayers of liquids. The initial evolution of the bedforms, based on the most unstable wavelengths, is completely absent in current numerical codes on this subject.

The next section discusses the involved physics and presents the equations composing the model. The following sections describe the numerical implementation of the model in a computational code, and the main results from the numerical simulations. They are followed by the conclusions section.

Nomenclature

A = constant B = constant C = constant c = phase velocity, m.s^-1 d = mean grain diameter, m F = Fourier transform F ^-1 = inverse Fourier transform g = acceleration of gravity, m.s^-2 H = amplitude of the initial bedform, m h = local height of the granular bed, m i = imaginary number k = wave-number, m^-1 L = length-scale, m N = number of iterations Q = volumetric flow rate of grains by unit of width, in the basic
state, m².s^-1 q = local volumetric flow rate of grains by unit of width, m².s^-1 u_* = shear velocity, m.s^-1 US = grain settling velocity, m.s^-1 Re_s = settling Reynolds number S = ratio between the specific masses S_d = Standard deviation, m t = time, s X = longitudinal length of the computation domain, m x = horizontal (longitudinal) coordinate, m y = vertical coordinate, m y₀ = rugosity height, m

Greek Symbols

Δ = interval (step) θ = Shields number

κ = von Kármán constant λ = wavelength of the initial instabilities, m µ = dynamic viscosity, Pa.s ρ = specific mass, kg.m^-3 σ = growth rate, s^-1 τ = shear stress, Pa χ = mean of the Gaussian function, m ω = frequency, s^-1 φ = concentration (volume fraction) of grains

Subscripts

d

= relative to deposition

drag = inertial scale k = Fourier space max = relative to the most unstable (amplified) mode 0 = relative to the basic state p = relative to grains s = relative to settling sat = relative to the saturated regime x = real space

Superscripts

^

= perturbation

Mathematical Model

The present model follows the lines of Kroy et al. (2002a), Kroy et al. (2002b) and Hersen (2004), and is also based on the following four equations: (i) the shear stress caused by the fluid flow on the bed surface; (ii) a semi-empirical equation for the bed-load flow rate; (iii) an equation accounting for the relaxation between the fluid flow and the granular flow; and (iv) the mass conservation of the granular bed. However, these equations are worked at differently here, mainly with respect to relaxation effects, and with respect to the inclusion of gravity effects. Also, differently from Kroy et al. (2002a), Kroy et al. (2002b) and Hersen (2004), the present model is applicable to turbulent liquid flows, and is focused on the early stages of the granular bed instabilities.

As the interest here is in the first stages of the bed instabilities, the model is two-dimensional. This is justified by the Squire's Theorem, which states that the most unstable modes in parallel flows are two-dimensional (Schmid and Henningson, 2001; Drazin and Reid, 2004; Charru, 2007). The physical concepts and the equations employed in the model are presented next.

-Shear stress on the bed

The perturbation of a turbulent boundary-layer by a hill with small aspect ratio was analytically found by Jackson and Hunt (1975) and by Hunt et al. (1988). Their results were later applied to forms with higher aspect ratio by Weng et al. (1991). Jackson and Hunt (1975), Hunt et al. (1988) and Weng et al. (1991) found that the perturbed shear stress is shifted upstream the dune crest. Kroy et al. (2002a) and Kroy et al. (2002b) simplified the results of Weng et al. (1991) and obtained an expression containing only the dominant physical effects of this perturbation, making clearer the reasons for this upstream shift. For a two-dimensional hill with a height h(x),a surface rugosity y₀ and a length 2L between the half-heights (total length ≈ 4L), they showed that the perturbation of the longitudinal shear stress, in dimensionless form and in the Fourier space is

where k= λ^-1 is the longitudinal wave-number (λ is the wavelength), i is the imaginary number, the subscript k is related to the Fourier space, F is the Fourier operator, x is the longitudinal direction and A and B are considered as constants, as they vary with the logarithm of L/y0. Equation (1) was obtained for H/L < 0.05, but Carruthers and Hunt (1990) showed that reasonable results are obtained when Eq. (1) is applied to slopes up to H/L = 0.3.

The fluid flow over the bed can be written as a basic state flow, unperturbed, plus a flow perturbation. The basic state is defined as the fluid flow, and the corresponding bed-load flow rate, over a flat granular bed. The shear stress on the bed surface can then be written, in the real space, as

where τ0 is the shear stress caused by the basic state flow on the bed and F^-1 is the inverse Fourier operator.

The shear stress caused by the fluid flow on the bed surface can be obtained directly from Eqs. (1) and (2). The great advantage of this method is that there is no need to compute the fluid flow in regions far from the surface, as it would be necessary with, for example, a RANS (Reynolds Average Navier-Stokes) method.

Flow rate of grains in the basic state

In the basic state, the fluid flow and the transport of grains are in equilibrium: their flow rates are in a steady state regime without spatial variations, so that the term "saturated" is often employed. However, as pointed out in the introduction, there are many different formulations for the flow rate of grains, and they are all semi-empirical. The one employed here is from Meyer-Peter and Mueller (1948), based on exhaustive experiments with water streams

where q_sat is the volumetric flow rate of grains by unit of width (saturated), S is the ratio between the grains specific mass ρ_p and the fluid specific mass ρ, g is the acceleration of gravity, d is the mean grain diameter, θ is the Shields parameter (ratio of the entraining force to the resisting force):

and θt is the threshold Shields parameter (Buffington and Montgomery, 1997).

-Relaxation between the flow rate of grains and the fluid flow

For fluid flows over undulated beds, the shear stress on the bed is a function of the position. It is then expected some inertial (or relaxation) effect between the grains and the fluid, so that the flow rate of grains will lag some distance with respect to the fluid flow. This distance is a characteristic length usually called "saturation length", L_sat. A simplified expression taking into account the relaxation effect can be obtained from the erosion-deposition model of Charru et al. (2004):

Kroy et al. (2002a), Kroy et al. (2002b) and Hersen (2004) considered that the saturation length has an inertial origin, and is proportional to the traveling distance of individual grains, given by L_drag=d ρ_p/ρ. According to this expression, L_drag is an inertial lengthscale obtained when the density of the granular material is many times larger than the density of the fluid, ρ_p >> ρ. It is then pertinent when the fluid is a gas. When the fluid is a liquid, however, ρ_p ≈ ρ and it has been argued by Charru (2006) and Franklin (2010) that this lengthscale can no longer be applied. Instead, a relaxation length based on the deposition of an individual grain must be used, L_d:

where u*₀ is the shear velocity of the basic state flow, defined by and Us is the grain settling velocity. The order of magnitude of L_sat, O(L_sat), has been reported as O(L_d) < O(L_sat) < 10.O(L_d) (Charru et al., 2004; Franklin, 2008), depending on the particular conditions of the liquid flow: O(L_sat) = O(L_d) in hydraulic smooth regimes and O(L_sat) = 10.O(L_d) in hydraulic rough regimes. C_sat is then an adjustable constant that, following the inertia of the boundary-layer at the grain scale, gives the proportion between L_sat and L_d. For a large range of the settling Reynolds number (R_eS = ρUSd/µ), the settling velocity Us may be evaluated as

where CD is the drag coefficient that may be evaluated by the Schiller-Neuman correlation when ReS < 800 (Clift et al., 1978).

-Mass conservation of grains

The two-dimensional mass conservation applied to the grains can be written as

where φ is the concentration (volume fraction) of grains in the granular bed, considered as constant in the model.

-Gravity effects

As discussed in Franklin (2010), gravity weakens the transport of grains over positive slopes (upstream the crests) and facilitates it over negative slopes (downstream the crests), being inversely proportional to the slope of the bed. It can then be incorporated in the constant B of the shear stress perturbation (Eq. (1)).

-Summary

Computations employing the model are relatively fast and easy. First, there is no need to compute the fluid flow in regions above the bed because the model needs only the shear stress on the bed surface, whose analytical solutions are known from Perturbation Methods. Second, the set of equations to be solved are uncoupled, and two of the four basic equations have analytical solutions (the others are first order differential equations). If an explicit scheme is employed in Eq. (8), numerical solutions can be searched by evaluating in sequence the presented equations, given an upstream condition (the boundary condition), and an initial bedform (initial condition). The structure employed to perform the numerical computations is given next.

Numerical Implementation

The model was implemented in a numerical code whose structure is given below.

1. Entry:

1. Entry: a.

Fluid and grain properties;

b. Initial condition (initial h(x)); c. Boundary condition (u*₀);

d.

Numerical parameters; e. Computation of L_sat by Eq. (6).

2.

Loop (until the desired number of iterations, or the total time is achieved):

a. Fourier transform of h(x) and computation of by Eq. (1); b. Inverse Fourier transform of and computation of τ by Eq. (2); c. Computation of q_sat(x) by Eq. (3); d. Computation of q(x) by the numerical solution of Eq. (5); e. Computation of the values of h(x) at the new time step by solving Eq. (8); f.

Storage of the data at the desired instants.

3. Post-processing: a. Computation of the bedform growth rate;

b.

Computation of the bedform celerity. A brief explanation of some of the numerical steps is given next.

-Fluid and grain properties (1.a): these concern the specific masses of the liquid and of the grain, the dynamic viscosity of the liquid, the mean diameter of the grains and the threshold Shields number (Buffington and Montgomery, 1997).

-Initial condition (1.b): it is the initial form of the bed h(x). For all the simulations presented here the initial bedform was a Gaussian Function

where H is the initial amplitude (crest), ξ is the Gaussian mean and s_d is the standard deviation. This form has the advantage of tending to zero when the domain is large enough, meaning that the basic state is expected to exist at the boundaries of the domain. Another advantage is that it can be easily adjusted to different lengths, positions and aspect ratios by varying sd, ξ and H. In order to be coherent with Jackson and Hunt (1975), Hunt et al. (1988) and Weng et al. (1991), it is considered here that L = sd so that the total length is approximately 4s_d.

-Boundary condition (1.c): is the condition far upstream the initial bedform, over a flat bed, corresponding then to the basic state.

-Numerical parameters (1.d): they correspond to the numerical constants necessary to the code functioning. They are: the spatial resolution Δx, the time step Δt, the total number of iterations N and the numerical scheme to be employed in the numerical solution of Eq. (8).

-Computation of τ(x) (2.a and 2.b): it is divided in four steps. First, a Fast Fourier Transform (FFT) operation is made on h. Then, is evaluated by Eq. (1) and an Inverse Fast Fourier Transform (IFFT) operation is made on . Finally, τ(x) is evaluated by Eq. (2).

-Computation of q(x) (2.d): q(x) is evaluated by the numerical integration of Eq. (5). An Upwind scheme is employed, with the boundary condition q(x = 0) = q_sat(x = 0).

-Computation of h(x) for the new time step (2.e): h(x) is evaluated by the numerical solution of Eq. (8). As q(x) and q_sat(x) are already known, an explicit numerical scheme can be easily employed to determine h(x) in the new time step. In the code, one can choose from four schemes that were implemented: LAX + FTCS, Upwind + LAX, Upwind + First Order Euler, FTCS + First Order Euler (Press et al., 1992). Due to the very small phase shifts between q(x) and h(x) in some cases, the First Order Euler schemes work better than the LAX schemes: the schemes based on LAX become more dependent on the spatial resolution because they employ spatial averages on the preceding time steps. This phase shift is the responsible for the stability or the instability of the bedform (Franklin, 2010). The use of an explicit scheme for the numerical solution of Eq. (8) means that the Courant condition must be verified. This was assured for all the simulations.

-Computation of the bedform growth rate and of the bedform celerity (3.a and 3.b): the growth rate σ and the bedform celerity c are evaluated from the transversal and the longitudinal displacement velocities of the bedform crest, respectively. Those velocities are evaluated by fitting the displacements as linear functions of the time.

Main Results

Simulations were performed employing the numerical code described in the previous section, and the main results are presented here. For the present simulations, the aspect ratio of the initial bedform, the mean diameter of the grains and the fluid properties (water) were fixed. In order to compare the numerical model with the linear stability analysis of Franklin (2010), both the length of the initial bedform and the shear velocity were varied. The choice of these parameters is based on the main conclusion of Franklin (2010) that the length-scale of the most unstable mode varies with the fluid flow, which is different from previous stability analyses. Effects due to gravity (slope) and to the grains diameter are not discussed in this paper. The bedform evolution was then computed and its growth rate and celerity were evaluated.

Figure 1 shows three examples of bedform development predicted by the numerical code. The initial form, shown in dashed line, was the same for the three cases: a Gaussian function with s_d = 0.1 m and aspect ratio H/(4sd) = 0.1. For the simulation presented in Fig. 1, the total number of iterations was N= 10⁵ with a time step Δt = 10^-5s, corresponding then to a total time of 1s. The total domain in the x direction was X= 2 m and it was discretized in intervals of Δx = 10^-3m. The saturation constant is C_sat = 1, corresponding to a hydraulic smooth regime. Figure 1.a corresponds to u*₀ = 0.04 m/s, Fig 1.b corresponds to u*₀ = 0.08 m/s and Fig. 1.c corresponds to u*₀ = 0.16 m/s.

From Figs. 1.b and 1.c it is clear that the simulated cases are unstable: the bedform grows as the time increases, while the form is displaced to the right, showing a positive celerity. The same occurs in Fig. 1.a, however the time-scale is too large (slow growth rate) to be clear in the figure.

The slope of the upwind face (upstream the crest) decreases, while the slope of the lee face (downstream the crest) increases. This indicates that the initial Gaussian bedform tends to a transverse ripple, which has an upstream face of small slope, and a high slope lee face (in fact an avalanche slope). This is exactly what is verified experimentally (Kuru et al., 1995; Franklin, 2008), so that the model is able to predict the tendency of the bedform evolution, in the early stages of the instability (or stability). Also from Fig. 1, it is clear that the growth rate and the celerity of the bedforms increase with the shear velocity of the basic flow.

Franklin (2010) showed from the mass conservation of grains that there is erosion in regions where the gradient of the bed-load flow rate is positive, and deposition where it is negative, so that the phase lag between the flow rate of grains and the bedform is a stability criterion. In other words, if the maximum of the flow rate of grains occurs upstream the bedform crests, deposition must occur at the crests: the bed is unstable and bedforms will grow. On the contrary, the bed is stable and bedforms tend to disappear. Franklin (2010) then discussed three reasons for the phase lag. One is the fluid flow perturbation caused by the shape of the bed, which is shifted upstream the bed and is an unstable mechanism (Jackson et al., 1975; Hunt et al., 1988; Weng et al., 1991). Another one concerns the relaxation effects related to the transport of grains: the grains lag some distance (or time) with respect to the fluid flow, so that the bed-load flow rate is shifted downstream the fluid flow and is a stable mechanism (Valance and Langlois, 2005, and Charru, 2006, in the case of viscous flows; Franklin, 2010, in the case of turbulent flows). The last one is related to gravity (slope effect): the bed-load flow rate is larger over negative slopes, so that it is shifted downstream the fluid flow and is another stable mechanism.

In the present simulations we do not analyze slope effects, so that the relaxation is the only stable one. In order to verify if the model is able to predict also stable solutions for the same initial bedform, the relaxation effects can be increased by increasing the value of C_sat. From experiments, C_sat can be estimated to be of order 10 in the transition to the hydraulic rough regime (Kuru et al., 1995; Franklin, 2008). Figure 2 shows the bedform development predicted by the numerical code when Csat = 50. The initial form, shown in dashed line, was the same as that of Fig. 1: a Gaussian function with s_d = 0.1 m and aspect ratio H/(4s_d) = 0.1. Figure 2a corresponds to u*₀ = 0.04 m/s, Fig 2.b corresponds to u*₀ = 0.08 m/s, Fig. 2.c corresponds to u*₀ = 0.16 m/s and all the other parameters are the same as in Fig. 1.

From Fig. 2b we observe that, due to the increased relaxation effect (increased phase shift between q and q_sat, as shown by Eq. 5), the growth rate is smaller when C_sat = 50 in comparison with C_sat = 1 (Fig. 1.b). Figure 2.c shows that for u*₀ = 0.16 m/s the relaxation effects are even larger, being stronger than the unstable mechanism and being responsible for the stabilization of the bed. The effect of a greater relaxation (smaller growth rate) also occurs in Fig. 2.a (in comparison to Fig. 1.a); however, the time-scale is too large (slow growth rate) to be clear in the figure.

The initial growth rate and the celerity of the bedform were computed from the evolution predicted by the numerical code. For the present simulations, the growth rate σ was evaluated as the transversal displacement velocity of the crest divided by the total length of the initial bedform, and the bedform celerity c was evaluated as the longitudinal displacement velocity of the crest. Table 1 shows the celerity c and the growth rate σ obtained in the simulations of Figs. 1 and 2.

From Table 1, when C_sat = 1, variations of the shear velocity of the basic flow u*₀ by a factor j imply variations in both the bedform celerity c and in the growth rate σ by a factor j³. The variation of c and of σ, as u*₀ 3, indicates that c and σ vary directly with the saturated flow rate of grains in the basic state (Franklin, 2010): the effect of the saturation length Lsat is too weak to affect the celerity or the growth rate. In fact, as discussed in Franklin (2010), the modulation of length-scale, celerity and growth rate depends on the phase shift between the fluid flow and the bed morphology (unstable effect) and on the phase shift between the bed-load flow rate and the fluid flow (stable effect). When Csat = 1, the latter, although exists, seems too small (many times smaller than the former) to affect c and σ. Then, c and σ vary only with the saturated flow rate of grains. When the value of Csat is increased to Csat = 50, the phase shift between the bed-load flow rate and the fluid flow is increased (Eqs. (5) and (6)) and both c and σ show variations that are lower than u*₀³, indicating that relaxation effects are present.

The present model can be directly compared to stability analyses if the evolutions of bedforms with different initial lengths are computed. The reason for this is that stability analyses are made on a continuum spectrum of modes. For example, this was done by Franklin (2010), who found a most unstable mode (most amplified growth rate) of a given wavelength, and predicted the variations of the growth rate and of the celerity for this specific wavelength. Franklin (2010) also showed the existence of a cut-off wavelength, below which the bed is always stable.

Figure 3 compares the celerity predictions from the linear stability analysis of Franklin (2010) with the celerity computed with the present model. The lines correspond to Eq. (19) of Franklin (2010) and the symbols to the celerity computed with the model. The dashed line and the circles correspond to u*₀ = 0.04 m/s and the continuous line and the asterisks correspond to u*₀ = 0.08 m/s. Figure 3.a presents the celerity c as a function of the wave-number k and Fig. 3.b presents the celerity normalized by the settling velocity c/Us as a function of the wave-number normalized by the saturation length k.L_sat. Because of the Courant condition, the total time (or the number of iterations) had to be larger in order to obtain the bedform celerity, and this was done for only some wave-numbers (but within a large range in the spectrum). The celerity computed numerically from the bedform evolution agrees with the stability analysis predictions: the simulated points c(k) follow the analytical curve c(k) given by Franklin (2010). Also, for all the modes (including the most unstable one) the numerical model gives the scaling c ~ u ², *0 as predicted by Franklin (2010).

Figure 4 compares the growth rate predictions from the linear stability analysis of Franklin (2010) with the growth rates computed with the present model. Figure 4.a presents the growth rate σ as a function of the wave-number k in dimensional form. The lines correspond to the Eq. (17) of Franklin (2010) and the symbols to the growth rates computed with the model. The dashed line and the circles correspond to u*₀ = 0.04 m/s and the continuous line and the asterisks correspond to u*₀ = 0.08 m/s.

The bedform evolution predicted by the model agrees with the stability analysis, showing the existence of a most unstable mode (the growth rate has a maximum) and that it is long-wave instability: bedforms increase (instability) if their initial wave-numbers are smaller than the cut-off wave-number, and bedforms decrease (stability) otherwise. The value of the cut-off wave-number predicted by the stability analysis is slightly lower, with the curve σ(k) decreasing (after the maximum) faster than the points computed with the model. This difference comes from the simplifications in the flow rate of grains done in Franklin (2010): (i) the multiplier factor was made equal to 1 and the threshold stress was neglected in the equation for the flow rate of grains (Eq. (5) of Franklin, 2010); and (ii) the latter was linearized (Eq. (7) of Franklin, 2010). Concerning the most unstable mode the agreement is good: the differences between the wave-number for the most unstable mode k_max predicted from the stability analysis and the computed numerically are less than 6% when u*₀ = 0.04 m/s and less than 5% when u*₀ = 0.08 m/s. Also, the value of the growth rate of the most unstable mode σ_max is multiplied by a factor 2 when the shear velocity u*₀ is doubled. This means that the numerical model captures the scaling σ_max ~ u*₀ predicted by Franklin (2010).

Figure 4.b presents the growth rate σ as a function of the wavenumber k in dimensionless form: σ is normalized by the deposition time td = d/Us and k by the saturation length Lsat. The symbols and the lines are the same as in Fig. 4.a. As the grains properties are the same for all the analyzed cases, the deposition time is the same for both cases (td = 8.1.10^-3s), so that the ratio of the magnitudes of σ is maintained when comparing u*₀ = 0.04 m/s to u*₀ = 0.08 m/s. The effect of normalizing k by L_sat is to bring the two cases to the same positions in the normalized abscissa, including the normalized wave-numbers for the most unstable mode and for the cut-off. This indicates that the most unstable mode scales with the saturation length L_sat, in agreement with Franklin (2010), but also with the experimental results of Kuru et al. (1995) and Franklin (2008).

Experiments in the scope of this model were presented in Kuru et al. (1995) and in Franklin (2008). Kuru et al. (1995) performed experiments on a 7 m long, 31.1 mm diameter horizontal pipe, and employed mixtures of water and glycerin as the fluid media and glass beads as the granular media. Franklin (2008) performed experiments on a 6 m long, horizontal closed-conduit of rectangular cross-section (120 mm wide by 60 mm high), made of transparent material, and employed water as the fluid and glass and zirconium beads as the granular media. In both works, the authors measured the wavelengths of the initial bedforms (ripples) appearing on the granular bed. Given the small time scales of the problem and the presence of high uncertainties, the celerity and the growth rate were not reported.

In order to directly compare the experimental data with the present computations, Fig. 5 presents the dimensionless mean wavelength λ/d of initial ripples as a function of the dimensionless shear velocity u*/U_s. From this figure, we can observe that the dimensionless wavelength increases with the fluid flow conditions, and that the experimental data seems dispersed around a line. If we take into account the relatively high uncertainties, often present in measurements of bedforms instabilities, the alignment of the experimental data in Fig. 5 supports the results of the proposed model (but also of Franklin, 2010): that λ ~ L_sat ~ d u*/U_s for the most unstable mode.

Finally, the numerical simulations with the present model agree with the stability predictions concerning the growth rate, the celerity, the most unstable wavelength and the cut-off wavelength. However, differently from the stability analysis, the model does not select the most unstable wavelength, so that it needs an initial bedform as an entry. If an appropriate routine based on the stability analysis is included in the model, it can select the most unstable mode as an entry. As the growth is exponential in the initial phase of the instability, the most unstable mode prevails, so that the model needs to compute the evolution of this sole mode.

Conclusions

This paper presented a simple and fast computing model for the transport of grains as bed-load by a turbulent boundary-layer of a liquid. Although simple, the model is able to capture the pertinent physics involved, such as the growing of instabilities. The main differences from previous models are that it concerns specifically liquid flows and that it is focused on the early stages of the bedform evolution. If employed together with the stability analysis, it can predict the evolution of the most unstable mode as long as the aspect ratio is small (smaller than 0.3): as the growth is exponential in the initial phase of the instability, the most unstable mode prevails, so that the model needs to compute the evolution of this sole mode. The proposed model agrees with the stability analysis of Franklin (2010) and compares well with the available experimental data of Kuru et al. (1995) and Franklin (2008). The initial evolution of the bedforms, based on the most unstable wavelengths, is completely absent in current numerical codes on this subject. The proposed model can then be included in current numerical codes in order to correctly predict the very first forms appearing on the bed, giving then a complete picture of the bedforms evolution.

Acknowledgements

The author is grateful to Petrobras for the financial support to write this article (contract number 0050.0045763.08.4).

Paper received 28 April 2011.

Paper accepted 03 November 2011.

Technical Editor: Francisco Cunha

• Bagnold, R.A., 1941, "The physics of blown sand and desert dunes", Ed. Chapman and Hall, London, United Kingdom, 320 p.
• Buffington, J.M. and Montgomery, D.R., 1997, "A systematic analysis of eight decades of incipient motion studies, with special reference to gravelbedded rivers", Water Resour. Res., Vol. 33, pp. 1993-2029.
• Carruthers, D.J. and Hunt, J.C.R., 1990, "Fluid mechanics of airflows over hills: turbulence, fluxes, and waves in the boundary layer", Atmospheric Processes over Complex Terrain, Vol. 23, pp. 83-108.
• Clift, R., Grace, J. and Weber, M., 1978, "Bubbles, drops and particles", Ed. Dover Publications, Inc., New York, United States of America, 400 p.
• Charru, F., Mouilleron-Arnould, H. and Eiff, O., 2004, "Erosion and deposition of particles on a bed sheared by a viscous flow", J. Fluid Mech., Vol. 519 pp. 55-80.
• Charru, F., 2006, "Selection of the ripple length on a granular bed sheared by a liquid flow", Physics of Fluids, Vol. 18 (121508).
• Charru, F., 2007, "Instabilités hydrodynamiques", Ed. EDP Sciences, Les Ulis, France, 386 p.
• Claudin, P. and Andreotti, B., 2006, "A scaling law for aeolian dunes on Mars, Venus, Earth, and for subaqueous ripples", Earth Plan. Sci. Lett., Vol. 252, pp. 20-44.
• Drazin, P.G. and Reid, W.R., 2004, "Hydrodynamic stability", Ed. Cambridge University Press, Cambridge, United Kingdom, 605 p.
• Elbelrhiti, H., Claudin, P. and Andreotti, B., 2005, "Field evidence for surface-wave induced instability of sand dunes", Nature, Vol. 437 (04058).
• Engelund, F., 1970, "Instability of erodible beds", J. Fluid Mech., Vol. 42, pp. 225-244.
• Engelund, F. and Fredsoe, J., 1982, "Sediment ripples and dunes", Ann. Rev. Fluid Mech., Vol. 14, pp. 13-37.
• Franklin, E.M., 2008, "Dynamique de dunes isolées dans un écoulement cisaillé", (in french), Ph.D. Thesis, Université de Toulouse, Toulouse, France, 166 p.
• Franklin, E.M., 2010, "Initial instabilities of a granular bed sheared by a turbulent liquid flow: length-scale determination", J. Braz. Soc. Mech. Sci. Eng., Vol. 32, pp. 460-467.
• Franklin, E.M., 2011, "Non-linear instabilities on a granular bed sheared by a turbulent liquid flow", J. Braz. Soc. Mech. Sci. Eng., accepted.
• Hersen, P., 2004, "On the crescentic shape of barchan dunes", Eur. Phys. J. B, Vol. 4, pp. 507-514.
• Hunt, J.C.R., Leibovich, S. and Richards, K., 1988, "Turbulent shear flows over low Hills", Quart. J. R. Met. Soc., Vol. 114, pp. 1435-1470.
• Jackson, P.S. and Hunt, J.C.R., 1975, "Turbulent wind flow over a low hill", Quart. J. R. Met. Soc., Vol. 101, pp. 929-955.
• Kennedy, J.F., 1963, "The mechanics of dunes and antidunes in erodible-bed channels", J. Fluid Mech., Vol. 16, pp. 521-544.
• Kroy, K., Sauermann, G. and Herrmann, H.J., 2002a, "Minimal model for sand dunes", Phys. Rev. Letters, Vol. 88 (054301).
• Kroy, K., Sauermann, G. and Herrmann, H.J., 2002b, "Minimal model for aeolian sand dunes", Phys. Rev. E, Vol. 66 (031302).
• Kuru, W.C., Leighton, D.T. and McCready M.J., 1995, "Formation of waves on a horizontal erodible bed of particles", Int. J. Multiphase Flow, Vol. 21, No. 6, pp. 1123-1140.
• Landau, L.D. and Lifchitz, E.M., 1994, "Physique théorique: mécanique des fluides", Ed. Ellipses (traduction française), Poitiers, France, 752 p.
• Meyer-Peter, E. and Mueller, R., 1948, "Formulas for bed-load transport", Proceedings of the 2nd Meeting of the International Association for Hydraulic Research, Stokholm, Sweden.
• Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P., 1992, "Numerical Recipes in fortran: the art of scientific computing", Ed. Cambridge University Press, New York, United States of America, 933 p.
• Raudkivi, A.J., 1976, "Loose boundary hydraulics", Ed. Pergamon, Oxford, United Kingdom, 397 p.
• Reynolds, A.J., 1964, "Waves on the erodible bed of an open channel", J. Fluid Mech., Vol. 22, pp. 113-133.
• Richards, K.J., 1980, "The formation of ripples and dunes on an erodible bed", J. Fluid Mech., Vol. 99, pp. 597-618.
• Sauermann, G., Kroy, K. and Herrmann, H.J., 2001, "A continuum saltation model for sand dunes", Phys. Rev. E, Vol. 64 (31305).
• Schmid, P.J. and Henningson, D.S., 2001, "Stability and transition in shear flows", Ed. Springer, New York, United States of America, 556 p.
• Valance, A. and Langlois, V., 2005, "Ripple formation over a sand bed submitted to a laminar shear flow", Eur. Phys. J. B, Vol. 43, pp. 283-294.
• Weng, W.S., Hunt, J.C.R., Carruthers, D.J. , Warren, A., Wiggs, G.F. S., Livingstone, I., and Castro, I., 1991, "Air flow and sand transport over sand-dunes", Acta Mechanica, pp. 1-21.

Publication Dates

• Publication in this collection
  10 Apr 2012
• Date of issue
  Mar 2012

History

• Received
  28 Apr 2011
• Accepted
  03 Oct 2011