Formation of sediment patterns in channel flow: minimal unstable systems and their temporal evolution

TL;DR

This study numerically investigates sediment pattern formation in channel flow, identifying minimal unstable system sizes and analyzing the evolution and characteristics of ripples over different domain lengths.

Contribution

It determines the minimal domain size for unstable sediment patterns and explores the nonlinear evolution of ripples in large-scale simulations.

Findings

The cutoff length for pattern formation is 75-100 particle diameters.

Ripple amplitude growth has exponential and nonlinear regimes.

Ripple shape exhibits a power-law spectrum decay.

Abstract

The phenomenon of sediment pattern formation in a channel flow is numerically investigated by performing simulations which resolve all the relevant scales of the problem. The numerical approach employed and the flow configuration considered is identical to our previous study (Kidanemariam and Uhlmann J. Fluid Mech., vol. 750, 2014, R2), the only difference being the length of the computational domain. By successively reducing the streamwise length, the minimum box dimension which accommodates an unstable sediment bed is revealed, thus determining the lower threshold of the unstable modes. For the considered parameter point, the cutoff length for pattern formation lies in the range 75-100 particle diameters (3-4 times the clear fluid height). We also simulate the flow in a very long streamwise box with a size of 48 times the clear fluid height (featuring over one million particles), accommodating approximately 11 initial ripple units with a wavelength of 100-110 particle diameters. The amplitude of the patterns exhibits two regimes of growth: an initial exponential regime, with a growth rate independent of the chosen domain size, and a subsequent non-linear regime which is strongly constrained by the domain length. In the small domain cases, after the initial exponential regime, the ripples evolve steadily maintaining their shape and migration velocity, at a mean wavelength equal to the length of the domain. The asymmetric ripple shape is characterized by a spectrum which exhibits a power-law decay over the first few dominant non-dispersive modes propagating at the mean dune migration velocity. The particle flowrate and the mean interface shear stress exhibited an increase with increasing ripple dimensions. Nevertheless, the relationship between the two was observed to be approximately described by the empirical power law formula for sediment transport by Wong & Parker (2006).

Tables

Table 1. Table 1: Physical parameters of the simulations. In all cases, the values of the bulk Reynolds number R ​ e b 𝑅 subscript 𝑒 𝑏 Re_{b} , the particle-to-fluid density ratio ρ p / ρ f subscript 𝜌 𝑝 subscript 𝜌 𝑓 \rho_{p}/\rho_{f}  and the Galileo number G ​ a 𝐺 𝑎 Ga  are imposed. Derived physical parameters include the friction Reynolds number R ​ e τ 𝑅 subscript 𝑒 𝜏 Re_{\tau} , the Shields number θ 𝜃 \theta  as well as the length scale ratios H f / D subscript 𝐻 𝑓 𝐷 H_{f}/D , H b / D subscript 𝐻 𝑏 𝐷 H_{b}/D and D + superscript 𝐷 D^{+} . Note that three independent simulations are performed for cases H4 and H12 each for the purpose of ensemble averaging. The values of the relevant parameters are determined over a time interval at the end of the respective simulations (cf. table  2 ).

Case	$R{e}_b$	$R{e}_{\tau }$	${\rho }_p/{\rho }_f$	$Ga$	$D^{+}$	$H_f/D$	$H_b/D$	$\theta$
H3	3011	244.8	2.5	28.37	9.62	25.44	12.96	0.12
H41	3011	273.6	2.5	28.37	10.89	25.12	13.28	0.15
H42	3011	265.4	2.5	28.37	10.55	25.15	13.25	0.14
H43	3011	263.0	2.5	28.37	10.45	25.17	13.23	0.14
H6	3011	303.0	2.5	28.37	11.90	25.47	12.93	0.18
H7	3011	309.1	2.5	28.37	12.22	25.30	13.10	0.19
H121	3011	301.3	2.5	28.37	12.01	25.08	13.32	0.18
H122	3011	301.6	2.5	28.37	12.04	25.05	13.35	0.18
H123	3011	298.6	2.5	28.37	11.91	25.07	13.33	0.18
H48	3011	293.1	2.5	28.37	11.69	25.07	13.33	0.17

Table 2. Table 2: Numerical parameters of the simulations. L i subscript 𝐿 𝑖 L_{i} is the domain length adopted in the i 𝑖 i -th direction. In all the cases, the spherical particles have a diameter D = 10 / 256 𝐷 10 256 D=10/256 and a uniform grid spacing Δ ​ x = Δ ​ y = Δ ​ z = 1 / 256 Δ 𝑥 Δ 𝑦 Δ 𝑧 1 256 \Delta x=\Delta y=\Delta z=1/256 is adopted yielding the above listed resolutions. N p subscript 𝑁 𝑝 N_{p}  is the number of spherical particles adopted. T o ​ b ​ s subscript 𝑇 𝑜 𝑏 𝑠 T_{obs} is the total simulation time interval staring from the instant at which the mobile particles are released, while T o ​ b ​ s s superscript subscript 𝑇 𝑜 𝑏 𝑠 𝑠 T_{obs}^{s} corresponds to the interval (at the end of the respective simulations) during which statistical averaging is performed to determine the relevant parameters listed in table  3 ( R ​ e τ 𝑅 subscript 𝑒 𝜏 Re_{\tau} , D + superscript 𝐷 D^{+} , H b / D subscript 𝐻 𝑏 𝐷 H_{b}/D , and θ 𝜃 \theta ) and to compute steady-state ripple quantities which are discussed in sections 5.2 , 5.3 , 5.4 and 5.5 . T b ≡ H f / u b subscript 𝑇 𝑏 subscript 𝐻 𝑓 subscript 𝑢 𝑏 T_{b}\equiv H_{f}/u_{b} is the bulk time unit.

Case	$\left[L_x\times L_y\times L_z\right]/D$	$D/\Delta x$	$\Delta x^{+}$	$N_p$	$T_{obs}/T_b$	$T_{obs}^s/T_b$
H3	$76.8\times 38.4\times 76.8$	10	0.96	65359	638	483
H41,2,3	$102.4\times 38.4\times 76.8$	10	1.06	86645	401/400/853	88/84/542
H6	$153.6\times 38.4\times 76.8$	10	1.19	127070	918	513
H7	$179.2\times 38.4\times 76.8$	10	1.22	150521	977	566
H121,2,3	$307.2\times 38.4\times 76.8$	10	1.20	263412	911/807/890	283/177/262
H48	$1228.8\times 38.4\times 76.8$	10	1.17	1053648	462	28

Equations

$$R{e}_b=\frac{u_b{H}_f}{\nu },$$

$$Ga=\frac{(({\rho }_p/{\rho }_f-1)|\mathbf{g}|D^3{)}^{\frac{1}{2}}}{\nu }=\frac{U_{g}D}{\nu },$$

$$\theta =\frac{u_{\tau }^2}{\left({\rho }_p/{\rho }_f-1\right)|\mathbf{g}|D}=\frac{u_{\tau }^2}{U_g^2}={\left(\frac{D^{+}}{Ga}\right)}^2,$$

$$h_b(x,t)=y|⟨{\varphi }_p{⟩}_z(x,y,t)=⟨{\varphi }_p{⟩}_z^{thresh}.$$

$$h_f(x,t)=L_y-h_b(x,t).$$

$$h_b^{\prime }(x,t)=h_b(x,t)-⟨h_b{⟩}_x(t).$$

$${\sigma }_h^2(t)=⟨h_b^{\prime }(x,t)\cdot h_b^{\prime }(x,t){⟩}_x.$$

$$R_h(\delta x,t)=\frac{⟨h_b^{\prime }(x,t)\cdot h_b^{\prime }(x+\delta x,t){⟩}_x}{{\sigma }_h^2(t)}.$$

$${\lambda }_h(t)=2x_{min}|\forall \delta x\in [0,L_x/2]:R_h(\delta x,t)\ge R_h(x_{min},t).$$

$${\sigma }_h/D=A\exp (Bt/T_b)$$

$$A_j(t)=A_{j0}\exp ({\sigma }_{j}t)$$

$${\tilde{H}}_b(\tilde{x})=⟨h_b(\tilde{x},t){⟩}_t,$$

$$⟨q_p⟩-⟨q_{p,min}⟩\approx \beta {\Phi }_b{u}_D{H}_D,$$

$$R_{ht}(\delta x,\delta t)=⟨h_b^{\prime }(x,t_1)\cdot h_b^{\prime }(x+\delta x,t_1+\delta t){⟩}_x\forall \delta t\in \left[0,t_2-t_1\right]$$

$$r_x^m(\delta t)=x_{max}(\delta t)|\forall \delta x\in [-L_x/2,L_x/2]:R_{ht}(\delta x,\delta t)\le R_{ht}(x_{max},\delta t),$$

$$\beta =\frac{1}{H_D{\lambda }_h}{\int }_0^{{\lambda }_h}({\tilde{H}}_b(x)-\min ({\tilde{H}}_b))\mathrm{d}x.$$

$${\phi }_j(t)=atan2(\Re ({\hat{h}}_{bj}(t)),\Im ({\hat{h}}_{bj}(t)))\forall j=0,...,N_x/2$$

$$c_j(t)=-\frac{1}{{\kappa }_j}\frac{\mathrm{d}{\phi }_j}{\mathrm{d}t}(t).$$

$$-\Big{\langle}\frac{{\rm d}p}{{\rm d}x}\Big{\rangle}(L_{y}-y)=\underbrace{\rho_{f}\nu\partial_{y}\langle u\rangle-\rho_{f}\langle u^{\prime}v^{\prime}\rangle+\int_{y}^{L_{y}}\langle f_{x}\rangle{\rm d}y}_{\langle\tau_{tot}\rangle}\;$$

$$q_p(t)=\frac{\pi D^3}{6L_x{L}_z}\sum _{l=1}^{N_p}{u}_p^{(l)}(t),$$

$$⟨q_p⟩/q_{ref}=A{\tilde{\theta }}^{\alpha },$$

Full text

In the small domain cases, after the initial exponential regime, the ripples evolve steadily maintaining their shape and migration velocity, at a mean wavelength equal to the length of the domain. The asymmetric ripple shape is characterized by a spectrum which exhibits a power-law decay over the first few dominant non-dispersive modes propagating at the mean dune migration velocity. The rate of particle transport and the mean interface shear stress exhibited an increase with increasing ripple dimensions. Nevertheless, the relationship between the two was observed to be approximately

described by the   empirical power law formula for sediment transport by Wong & Parker (2006).

1 Introduction

Sediment patterns, also commonly termed as dunes, ripples or simply bedforms, can be abundantly observed in deserts, coastal areas, natural streams or man-made canals. These sediment features, besides being simply fascinating, have important implications in many fields of science and engineering. For instance, bedforms greatly influence the rate of sediment transport as well as the stability of hydraulic structures in a given river. Thus, fundamental understanding of the mechanisms which are behind their formation as well as predicting their characteristics is crucial. However, this task has been very challenging due to the complex interaction between the sediment particles and the driving turbulent flow.

Sediment patterns in general can be broadly classified as ‘aeolian’ or ‘subaqueous’, the former being driven by wind forces while the latter, which the present study is focused on, form under water or due to the action of a similar fluid.

The sediment-to-fluid density ratio of subaqueous patterns is much smaller than that of the aeolian ones. As a result, although the basic formation process is the same, there are fundamental differences regarding the relevant controlling mechanisms of their generation in both categories (cf. Bagnold, 1941; Sauermann et al., 2001; Andreotti & Claudin, 2013; Duran et al., 2014).

Concerning subaqueous sediment patterns, one can further distinguish between those that form under the action of a steady flow, such as in rivers and canals, and those that are worked by an unsteady oscillatory flow, for instance by the action of a sea wave in coastal regions (cf. Sleath, 1976; Blondeaux, 1990; Vittori & Blondeaux, 1990; Blondeaux et al., 2015). Henceforth we will restrict our attention to the study of patterns occurring in steady flows. Although subaqueous bedforms are three-dimensional objects, under certain circumstances they tend to be statistically invariant in the cross-stream direction (Yalin, 1977).

Indeed, a large part of the previous theoretical modeling and experimental measurement has focused on two-dimensional subaqueous bedforms (Best, 2005), and this is the approach we take in the present work.

Bedforms can be characterized with respect to their spatial dimensions and their dynamics (Julien, 1998). Classically, a distinction is made between ripples, dunes and antidunes (Engelund & Fredsoe, 1982; García, 2008). Ripples and dunes exhibit some similarities, such as an asymmetric triangular-like shape with a gentle slope on the upstream side of their crest and a steeper slope on the downstream side. The principal distinction between ripples and dunes is the separation of length scales; the wavelength of observed ripples is commonly believed to scale with the grain size, and that of larger scale dunes is supposedly controlled by the flow geometry, typically the flow depth (Yalin, 1977). Experiments show that at low flow rates small-scale undulations appear out of an initially flat sediment bed, at first having no clear-cut shape (Coleman & Melville, 1994). When increasing the flow rate, these ripples give way to the much larger dunes, which are characterized by flow separation and recirculation in their downstream faces, and which propagate downstream. At even higher flow rates antidunes appear, which tend to propagate upstream. In the case of a flow with a free surface, both dunes and antidunes dynamically interact with surface waves (see e.g. Best, 2005). Nevertheless, bedforms are equally observed in closed conduits which possess no free surface, and thus the existence of a free surface is not necessary for their formation (Yalin, 1977; Bridge & Best, 1988; Langlois & Valance, 2007; Ouriemi et al., 2009; Cardona Florez & Franklin, 2016). The problem of bedform and free surface wave interaction is not within the scope of the present study.

In most previous theoretical work on pattern formation, the background turbulent flow is typically represented by a Reynolds averaged Navier-Stokes equations (RANS) model, while the sediment bed evolution is described by the sediment continuity equation, i.e. the Exner equation. The hydrodynamic and morphodynamic problems are then coupled by an algebraic expression for the particle flux as a function of the local bed shear stress. Hydro-morphodynamic linear (or weakly nonlinear) stability analysis is then performed in order to determine the stability of the sediment bed, the controlling parameter of the instability as well as the initial unstable pattern wavelength (e.g. see reviews by Seminara, 2010; Charru et al., 2013). However, there is no clear consensus among the different predictions from these approaches, and, when compared to experimental observations, the outcome of these models is often unsatisfactory. For instance, the prediction for the initial pattern wavelength can be off by an order of magnitude (Raudkivi, 1997; Langlois & Valance, 2007; Coleman & Nikora, 2009a; Ouriemi et al., 2009). This inadequacy can be linked to, among others, the insufficient predictive capability of the adopted algebraic expressions for the particle flux.

The subsequent bedform evolution is an even more complicated issue which includes aspects like coarsening (bedform wavelength and amplitude growth), asymmetry, coalescence (fusion of bedforms), 3D patterning (loss of bedform translational invariance). Even in simplified and controlled experiments where bedforms are essentially two-dimensional (see e.g. Betat et al., 2002), the evolution process is a result of nonlinear interaction mechanisms which are far from our understanding. Certainly, linear stability theories are not adequate in describing the bedform transient processes such as the evolution of the amplitude or the evolution of bedform morphology towards the asymmetric shape (which is even observed in the absence of flow separation). Weakly nonlinear stability theories have been proposed (see e.g. Colombini & Stocchino, 2008), but this approach is still unable to reliably predict the observed bedforms in reference experiments (e.g. Ouriemi et al., 2009).

There exists a very limited amount of experimental studies available reporting on the initial formation and development stages of patterns (see e.g. Coleman & Melville, 1994; Betat et al., 2002; Coleman et al., 2003; Langlois & Valance, 2007; Ouriemi et al., 2009; Cardona Florez & Franklin, 2016). Even these studies, limited by measurement difficulties, often fall short of providing details with respect to the very first instants of the bed instability, such as the dispersion relation of the unstable modes. Such quantities are crucial when it comes to assessing the validity of the various proposed theoretical models. On the other hand, although there are a number of numerical studies available on the problem of subaqueous sedimentary patterns, most of these studies are based on the continuum description of the sediment bed, which again rely on the use of (semi-)empirical relationships to couple the hydro-morphodynamic problem (see e.g. Chou & Fringer, 2010; Khosronejad & Sotiropoulos, 2014). In the granular flow community, a number of studies have been performed with the aid of a Lagrangian description of the sediment bed (based on discrete-element models), but without coupling to a resolved turbulent background flow (Durań et al., 2012; Schmeeckle, 2014; Maurin et al., 2015).

There has been a promising gain of momentum in recent years on the number of direct numerical simulation (DNS) studies on the problem of sediment transport in which the detail of the flow, even at the boundaries of each individual particles, is faithfully resolved (see e.g. Kidanemariam & Uhlmann, 2014b, a; Vowinckel et al., 2014; Derksen, 2015). We have recently performed novel DNS of the flow over an erodible bed of spherical sediment particles in a statistically unidirectional channel flow configuration in both the laminar and turbulent regimes (Kidanemariam & Uhlmann, 2014a). These simulations are, to the best of our knowledge, the first to successfully simulate the evolution of a bed of mobile sediment particles (leading to pattern formation) by means of DNS (Colombini, 2014). The results from these simulations have shown that the chosen streamwise length of the computational domain ($12H_{f}$, where $H_{f}$ is the channel fluid height) was able to accommodate a few integer multiples of the initial pattern wavelength, which allowed us to address some of the outstanding questions on sedimentary patterns such as the initial dominant wavelength and pattern amplitude growth (Kidanemariam & Uhlmann, 2014a). Ideally, in order to accurately capture the natural selection mechanism of the unstable wavelength and its subsequent evolution, it is necessary to consider a computational domain size which is much superior than the anticipated pattern wavelength. However, the domain sizes considered by Kidanemariam & Uhlmann (2014a) might be marginal when compared to the observed dominant wavelengths. It was noted that the discreteness of the numerical harmonics could influence the initial wavelength as well as the time evolution processes. In light of this aspect, assessing the influence of the domain size, on the selection and evolution process is crucial.

Moreover, various theoretical and experimental studies have shown that there is a lower threshold of the unstable wavelengths (see for instance Charru, 2006; Franklin & Charru, 2011). Determination of this value is important, since the lower bound of unstable modes is believed to be linked to some relevant length scale which controls erodible bed instability (Hersen et al., 2002; Claudin & Andreotti, 2006; Andreotti et al., 2011).

In the present work, we have carried out a series of simulations of the same flow configuration as in our previous study (Kidanemariam & Uhlmann, 2014a). We have kept all the physical and numerical parameters of the simulations identical, except for the streamwise length of the computational box. The latter was varied between $3H_{f}$ and $48H_{f}$, in order to determine whether at the cosidered parameter point, there exists a cutoff length for pattern formation. To this end, we have successively reduced the domain size until it is smaller than an unknown threshold and cannot accomodate an unstable bed. Moreover, in order to assess the influence of the domain size on the selection of the initial ripple wavelength and its subsequent evolution, we have chosen a computational domain size which is four times larger than in our previous study. It should be noted that a very large number of spherical particles (approximately 1.1 million in total) are considered to represent the mobile bed. This simulation is the first of its kind to break the $\mathcal{O}(10^{6})$ fully resolved particle milestone. Furthermore, based on the analysis of the DNS data, we address the relationship between the evolving patterns, their migration velocity and the particle flow rate at the present parameter point.

2 Numerical method

In the present work we have used the same numerical procedure as Kidanemariam & Uhlmann (2014a, b). The numerical treatment of the fluid-solid system is based upon the immersed boundary technique of Uhlmann (2005), wherein the incompressible Navier-Stokes equations are solved with a second-order finite-difference method throughout the entire computational domain $\Omega$, adding a localized force term which serves to impose the no-slip condition at the fluid-solid interface. The particle motion is obtained via integration of the Newton-Euler equations for rigid body motion, driven by the hydrodynamic force (and torque) as well as gravity and the force (torque) resulting from solid-solid contact. The collision process between the immersed particles is described through a discrete element model (DEM) based on the soft-sphere approach. A pair of particles is defined as ‘being in contact’ when the smallest distance between their surfaces, $\Delta$, becomes smaller than a force range $\Delta_{c}$. The resulting contact force is then the sum of an elastic normal component, a normal damping component and a tangential frictional component. The elastic part of the normal force component is a linear function of the penetration length $\delta_{c}\equiv\Delta_{c}-\Delta$, with a stiffness constant $k_{n}$. The normal damping force is a linear function of the normal component of the relative velocity between the particles at the contact point with a constant coefficient $c_{n}$. The tangential frictional force (the magnitude of which is limited by the Coulomb friction limit with a friction coefficient $\mu_{c}$) is a linear function of the tangential relative velocity at the contact point, again formulated with a constant coefficient denoted as $c_{t}$. A detailed description of the collision model and extensive validation can be found in Kidanemariam & Uhlmann (2014b).

The four parameters which describe the collision process in the framework of this model ($k_{n}$, $c_{n}$, $c_{t}$, $\mu_{c}$) as well as the force range $\Delta_{c}$ need to be prescribed for each simulation. Note that the normal stiffness coefficient $k_{n}$ and the normal damping coefficient $c_{n}$ can be related by introducing the dry restitution coefficient $\varepsilon_{d}$, defined as the absolute value of the ratio between the normal components of the relative velocity post-collision and pre-collision.

In the present simulations, $\Delta_{c}$ is set equal to one grid spacing $\Delta x$. The stiffness parameter $k_{n}$ has a value equivalent to approximately $17000$ times the submerged weight of the particles, divided by the particle diameter. The chosen value ensures that the maximum overlap $\delta_{c}$ over all contacting particle pairs is within a few percent of $\Delta_{c}$. The dry coefficient of restitution is set to $\varepsilon_{d}=0.3$ which together with $k_{n}$ fixes the value for $c_{n}$. Finally, the tangential damping coefficient $c_{t}$ was set equal to $c_{n}$, and a value of $\mu_{c}=0.4$ was imposed for the Coulomb friction coefficient. This set of parameter values for the contact model is the same as used by Kidanemariam & Uhlmann (2014b).

Since the characteristic collision time is typically orders of magnitude smaller than the time step of the flow solver, the numerical integration of the equations for the particle motion is carried out adopting a sub-stepping technique, freezing the hydrodynamic forces acting upon the particles between successive flow field updates (Kidanemariam, 2015).

3 Flow configuration and parameter values

We have performed a total of ten independent simulations of the development of bedforms over a subaqueous sediment in an open channel flow configuration. As shown in figure 1 a Cartesian coordinate system is adopted such that $x$, $y$, and $z$ are the streamwise, wall-normal and spanwise directions, respectively. Mean flow and gravity are directed in the positive $x$ and the negative $y$ directions respectively. The computational domain is periodic in the streamwise and spanwise directions. A free-slip condition is imposed at the top boundary while a no-slip condition is imposed at the bottom wall. The simulations are labeled H3, H4, H6, H7, H12 and H48, indicating the approximate streamwise box length in terms of the mean fluid height $H_{f}$. The mean fluid height $H_{f}$ and the corresponding mean sediment bed thickness $H_{b}$ are computed by performing streamwise and time averaging of the spanwise-averaged instantaneous fluid height $h_{f}$ and sediment bed height $h_{b}$ (note that the definition of $h_{f}$ and $h_{b}$ will be made more precise in section 4.1). In order to perform ensemble averaging, case H4 and H12 are performed three times each adopting different initial conditions. Each simulation is performed independently following the simulation start-up procedure as detailed in Kidanemariam & Uhlmann (2014a).

In all cases, the channel is driven by a horizontal mean pressure gradient which is adjusted at each time step in order to impose a constant flow rate $q_{f}$. This results in a shearing flow of fluid height $H_{f}$ over a mobile bed of height $H_{b}$ (cf. sketch in figure 1). As is shown in table 3, the bulk Reynolds number of the flow, which is defined based on the mean height of the fluid as

[TABLE]

where $u_{b}\equiv q_{f}/H_{f}$ is the bulk velocity and $\nu$ is the kinematic viscosity, is set at a value such that the flow is fully turbulent. The friction Reynolds number $Re_{\tau}$, which is similarly defined based on the friction velocity $u_{\tau}$, is a posteriori determined by evaluating the total shear stress at the wall-normal location of the mean fluid-bed interface $y=H_{b}$ (cf. section 5.4 for details of the determination of $u_{\tau}$).