Direct numerical simulation of pattern formation in subaqueous sediment
DDirect numerical simulation of pattern formation insubaqueous sediment
Aman G. Kidanemariam ∗ and Markus Uhlmann †∗ Institute for Hydromechanics, Karlsruhe Institute of Technology76131 Karlsruhe, Germany
Abstract
We present results of direct numerical simulation of incompressible fluid flow over a thickbed of mobile, spherically-shaped particles. The algorithm is based upon the immersed bound-ary technique for fluid-solid coupling and uses a soft-sphere model for the solid-solid contact.Two parameter points in the laminar flow regime are chosen, leading to the emergence of sedi-ment patterns classified as ‘small dunes’, while one case under turbulent flow conditions leadsto ‘vortex dunes’ with significant flow separation on the lee side. Wavelength, amplitude andpropagation speed of the patterns extracted from the spanwise-averaged fluid-bed interfaceare found to be consistent with available experimental data. The particle transport rates arewell represented by available empirical models for flow over a plane sediment bed in both thelaminar and the turbulent regimes.
The process of erosion of particles from an initially flat subaqueous sediment layer and theirdeposition at certain preferential locations leads, under certain circumstances, to the amplificationof small perturbations and gives rise to wave-like bed shapes which are commonly described asripples or dunes. These sedimentary patterns are commonly observed in river and marine flows, aswell as in various technical applications involving shear flow over a bed of mobile sediment particles.From an engineering point of view it is highly desirable to be able to predict the occurrence andthe nature of this phenomenon, since the bed-form significantly influences flow characteristics suchas resistance, mixing properties and sediment transport.Most of the previous theoretical work on the formation of sediment patterns is based uponthe notion that a flat bed is unstable with respect to perturbations of sinusoidal shape. It wasKennedy (1963) who first studied this instability problem by considering a potential flow solution,and over the years the concept has been applied by a number of researchers for a variety of flowconditions, in the laminar (Charru & Mouilleron-Arnould, 2002; Charru & Hinch, 2006) and tur-bulent regime (Richards, 1980; Sumer & Bakioglu, 1984; Colombini, 2004; Colombini & Stocchino,2011). Invoking a disparity in time scales between the flow and the bed shape modification, mostof the approaches have considered the bed shape as fixed for the purpose of the analysis. Thehydrodynamic stability problem is then complemented by an expression for the particle flux as afunction of the local bed shear stress at a given transversal section of the flow.It is now generally accepted that the mechanism which destabilizes a flat sediment bed is thephase-lag between the perturbation in bed height and the bottom shear stress as a consequence offluid inertia. A balance between this destabilizing mechanism and other stabilizing effects such asgravity (Engelund & Fredsoe, 1982) or phase-lag between bottom shear stress and the particle flowrate (Charru, 2006) is believed to result in instability of the bed at a certain preferred wavelength. ∗ [email protected] † [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] M a y a ) i j T c ( i,j ) T c ( i,j ) δ ( i,j ) c F c ( i,j ) n F c ( i,j ) t − F c ( i,j ) t − F c ( i,j ) n ∆ c ( b ) xyz g L y L z h f = 2 Hh b I ( c ) xyz g L y L z h f = Hh b I Figure 1: ( a ) Schematic showing the collision force F c ( i,j ) and torque T c ( i,j ) acting on particlewith index i during contact with particle j . The subscripts n and t indicate the normal andtangential components, respectively; ∆ c denotes the force range, and δ ( i,j ) c is the overlap length.( b ) Computational domain and coordinate definition in cases LC1 and LC2. ( c ) The same for caseTO1. The computational domains are periodic along the x - and z -directions with periods L x , L z ,respectively.Linear stability analysis is often applied to the problem in order to predict the most unstable wave-length; compared to experimental observations, however, predictions resulting from this approachcan be broadly described as unsatisfactory, sometimes predicting pattern wavelengths which areoff by an order of magnitude (Raudkivi, 1997; Langlois & Valance, 2007; Coleman & Nikora, 2009;Ouriemi et al. , 2009).Most available experimental studies report wavelengths of the developed bed-forms after theyhave undergone a coarsening process (the temporal evolution of the initial patterns to their ‘steady-state’ form), or possibly after they have coalesced with other bed-forms. There are several exper-imental studies which report on the initial wavelength and its development (Coleman & Melville,1994; Betat et al. , 2002; Coleman et al. , 2003; Langlois & Valance, 2007; Ouriemi et al. , 2009).However, the reported data is widely dispersed. Today it is still a challenge to capture the three-dimensional nature of the individual particle and fluid motion within the bed layer in a laboratoryexperiment, which in turn has hindered the assessment of the various theoretical approaches.In the present work we numerically investigate the development of subaqueous patterns in astatistically uni-directional channel flow configuration both in the laminar and turbulent regimes.A sufficiently large number of freely-moving spherical particles are represented such that they forma realistic sediment bed in the simulation. To our knowledge, no attempt to numerically simulatethe evolution of a bed of mobile sediment particles (leading to pattern formation) by means ofdirect numerical simulation (DNS), which resolves all the relevant length and time scales of theturbulent flow as well as the individual sediment particles, has been reported to the present date.The present study focuses on aspects related to the initial bed instabilities and their subsequentshort-time development. The numerical treatment of the fluid-solid system is based upon the immersed boundary tech-nique. The incompressible Navier-Stokes equations are solved with a second-order finite-differencemethod throughout the entire computational domain Ω (comprising the fluid domain Ω f and thedomain occupied by the suspended particles Ω s ), adding a localized force term which serves toimpose the no-slip condition at the fluid-solid interface. The particle motion is obtained via in-tegration of the Newton equation for rigid body motion, driven by the hydrodynamic force (and2ase Re b Re τ ρ p /ρ f Ga D + H/D
Θ [ L x × L y × L z ] /D ∆ x + D/ ∆ x N p LC1 700 32.31 2.5 2.42 1.20 26.92 0.25 307 . × . × . . × . × . . × . × . Re b and Re τ are the bulk andfriction Reynolds numbers, respectively, ρ p /ρ f is the particle-to-fluid density ratio, Ga is theGalileo number, D is the particle diameter, H is the equivalent boundary layer thickness, Θ theShields number, L i is the domain size in the i th coordinate direction, ∆ x is the uniform meshwidth and N p is the number of particles.torque) as well as gravity and the force (torque) resulting from inter-particle contact. Furtherinformation on the extensive validation of the direct numerical simulation (DNS) code on a wholerange of benchmark problems can be found in Uhlmann (2005), Uhlmann & Duˇsek (2014) andfurther references therein. The code has been previously employed for the simulation of variousparticulate flow configurations (Uhlmann, 2008; Chan-Braun et al. , 2011; Garc´ıa-Villalba et al. ,2012; Kidanemariam et al. , 2013).In the present case, direct particle-particle contact contributes significantly to the dynamics ofthe system. In order to realistically simulate the collision process between the immersed particles,a discrete element model (DEM) based on the soft-sphere approach is coupled to the two-phaseflow solver. The DEM used in the present work employs a linear mass-spring-damper systemto model the collision forces, which are computed independently for each colliding particle pair.Any two particles are defined as ‘being in contact’ when the smallest distance between theirsurfaces, ∆, becomes smaller than a force range ∆ c as illustrated in figure 1( a ). The collisionforce is computed from the sum of an elastic normal component, a normal damping componentand a tangential frictional component. The elastic part of the normal force component is a linearfunction of the penetration length δ c ≡ ∆ c − ∆, with a stiffness constant k n . The normal dampingforce is a linear function of the normal component of the relative velocity between the particles atthe contact point with a constant coefficient c n . The tangential frictional force (the magnitude ofwhich is limited by the Coulomb friction limit with a friction coefficient µ c ) is a linear function ofthe tangential relative velocity at the contact point, again formulated with a constant coefficientdenoted as c t . A detailed description of the collision model and extensive validation tests withrespect to available experimental data of a single particle colliding with a wall in a viscous fluidand in the case of bedload transport under laminar shear flow has been recently provided byKidanemariam & Uhlmann (2014).The four parameters which describe the collision process in the framework of this model ( k n , c n , c t , µ c ) as well as the force range ∆ c need to be prescribed for each simulation. From an analyticalsolution of the linear mass-spring-damper system in an idealised configuration (considering abinary collision of uniformly translating spheres in vacuum and in the absence of external forces),a relation between the normal stiffness coefficient k n and the normal damping coefficient c n canbe formed by introducing the dry restitution coefficient ε d . This latter quantity is a materialproperty, defined as the absolute value of the ratio between the normal components of the relativevelocity post-collision and pre-collision. In the present simulations, ∆ c is set equal to one gridspacing ∆ x . The stiffness parameter k n has a value equivalent to approximately 17000 times thesubmerged weight of the particles, divided by the particle diameter. The chosen value ensuresthat the maximum overlap δ c over all contacting particle pairs is within a few percent of ∆ c .The dry coefficient of restitution is set to ε d = 0 . 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 µ c = 0 . e b , R e p i p e Ga (2 H/D ) Figure 2: Different regimes of sediment bed patterns obtained in the pipe flow experiment ofOuriemi et al. (2009), shown in the parameter plane ( Re b , Ga (2 H/D ) ): ‘flat bed in motion’ ( (cid:4) );‘small dunes’ ( ◦ ); ‘vortex dunes’ ( (cid:78) ). For the pipe flow data the Reynolds number Re pipe basedupon the pipe diameter d pipe and the bulk velocity q f /d pipe is used. The data points in theturbulent channel flow experiment of Langlois & Valance (2007) are indicated by: D = 100 µm ( (cid:77) ); D = 250 µm ( (cid:67) ); D = 500 µm ( (cid:66) ). The following symbols refer to the present simulations: • , LC1; • , LC2; • , TO1.process and those of the smallest flow scales, the Newton equation for particle motion is solvedwith a significantly smaller time step than the one used for solving the Navier-Stokes equations(by a factor of approximately one hundred), while keeping the hydrodynamic contribution to theforce and torque acting upon the particles constant during the intermediate interval. The flow considered in this work is horizontal plane channel flow in a doubly-periodic domain asshown in figure 1( b, c ). Three simulations are performed, henceforth denoted as case LC1, case LC2(both in the laminar flow regime) and case TO1 (in the turbulent regime). In cases LC1 and LC2the domain is bounded in the vertical direction by two solid wall planes, whereas for reasonsof computational cost in case TO1 an open channel is simulated, i.e. only the lower boundaryplane corresponds to a no-slip wall, while a free-slip condition is imposed at the upper boundaryplane. As shown in figure 1 the Cartesian coordinates x , y , and z are aligned with the streamwise,wall-normal and spanwise directions, respectively, while gravity acts in the negative y -direction.The flow is driven by a horizontal pressure gradient at constant flow rate q f (per unit spanwiselength) which results in a shear flow of height h f over a mobile bed of height h b ; spatial averages (cid:104) h f (cid:105) x and (cid:104) h b (cid:105) x of both quantities are defined in § (cid:104)·(cid:105) t . The bulk Reynolds number of theflow is defined as Re b = 2 Hu b /ν , where u b ≡ q f / (cid:104) h f (cid:105) xt is the bulk velocity, H is the equivalentboundary layer thickness (i.e. H = (cid:104) h f (cid:105) xt / H = (cid:104) h f (cid:105) xt in case TO1, cf.figure 1 b, c ), and ν is the kinematic viscosity. Similarly, the friction Reynolds number is defined as Re τ = u τ (cid:104) h f (cid:105) xt /ν , where the friction velocity u τ is computed by extrapolation of the total shearstress to the fully-developed value of the wall-normal location of the average fluid-bed interface (cid:104) h b (cid:105) x . Further physical parameters are the ratio of particle to fluid density, ρ p /ρ f , the Galileonumber Ga = u g D/ν (where u g = (( ρ p /ρ f − | g | D ) / and D the particle diameter), the Shieldsnumber Θ = u τ /u g and the length scale ratio H/D ; these together with the chosen numericalparameters are shown in table 1. The present simulations consumed a total of approximately5 million core hours on the computing system SuperMUC at LRZ M¨unchen. Typical runs ofcase TO1 were carried out on 576 cores.Figure 2 shows the three parameter points of the present simulations in the plane spanned by4 a ) yD ( b ) yD x/D ( c ) yD ( d ) yD t u b /H Figure 3: ( a ) Close-up showing an instantaneous distribution of particles (colored according totheir streamwise velocity) in case TO1. ( b ) The spanwise-averaged two-dimensional solid volumefraction (cid:104) φ p (cid:105) z ( x, y, t ) corresponding to the snapshot shown in ( a ) displayed in greyscale. Theextracted fluid-bed interface location h b ( x, t ) is shown with a magenta colored line. ( c ) Timeevolution of the bed height: (cid:104) h b (cid:105) x ( t ) is shown as solid lines, min x ( h b ) and max x ( h b ) are shown asdashed lines. Black color is used for case LC1, red color for case LC2. ( d ) The same quantities asshown in ( c ), but for case TO1(see also the movies available as supplementary material). Re b and Ga (2 H/D ) ) in comparison to the laboratory experiments of Ouriemi et al. (2009) andLanglois & Valance (2007). Note that the former experiment was performed in pipe flow, whereasthe latter was in plane channel flow. It can be seen that the cases LC1 and LC2 fall into theregime where the formation of ‘small dunes’ is expected while ‘vortex dunes’ can be anticipatedin case TO1. The simulations were initiated as follows. In a first stage the initial sediment bed was generatedby means of a simulation of particles settling (from random initial positions) under gravity andunder solid-solid collisions but disregarding hydrodynamic forces. The result is a pseudo-randomlypacked bed of initial bed thickness h b ( t = 0) above the bottom wall. Then the actual fully-coupledfluid-solid simulation is started with all particles being initially at rest. In cases LC1 and LC2, theinitial fluid velocity field is set equal to a laminar Poiseuille flow profile with the desired flowratein the interval h b ( t = 0) ≤ y ≤ L y and zero elsewhere. After starting the simulation, individualparticles are set into motion due to the action of hydrodynamic force/torque, and erosion takesplace. In case TO1 the fluid-solid simulation was first run with all particles held fixed in order todevelop a fully-turbulent field over the given sediment bed. After approximately 100 bulk timeunits the particles were released, and the bed started to evolve away from its initial macroscopicallyflat shape, as discussed in the following. The location of the interface between the fluid and the sediment bed has been determined in thefollowing way. First, a solid phase indicator function φ p ( x , t ) is defined which measures unity if x islocated inside any particle and zero elsewhere. Spanwise averaging then yields (cid:104) φ p (cid:105) z ( x, y, t ) whichis a direct measure of the instantaneous, two-dimensional solid volume fraction. The spanwise-averaged fluid-bed interface location h b ( x, t ) is finally extracted by means of a threshold value,chosen as (cid:104) φ p (cid:105) threshz = 0 . h b ( x, t ) = y | (cid:104) φ p (cid:105) z ( x, y, t ) = (cid:104) φ p (cid:105) threshz , (1)5 a ) t u b / H x/D ( b ) t u b / H x/D ( c ) t u b / H x/D Figure 4: Space-time evolution of the fluctuation of the fluid-bed interface location, h (cid:48) b ( x, t ) = h b ( x, t ) − (cid:104) h b (cid:105) x ( t ), normalized with the particle diameter D : ( a ) case LC1; ( b ) case LC2; ( c ) caseTO1.as illustrated in figure 3 (b) . The corresponding spanwise-averaged fluid height is then simply givenby h f ( x, t ) = L y − h b ( x, t ). Figure 3 (c,d) shows the time evolution of the streamwise average ofthe bed height defined in (1), (cid:104) h b (cid:105) x ( t ), as well as its minimum and maximum values. It can beobserved that after a small initial dilation (cid:104) h b (cid:105) x quickly reaches an approximately constant valuein all three cases. Contrarily, the maximum and minimum values of the bed elevation continueto diverge over the simulated interval of approximately 1000 bulk time units, not reaching anequilibrium state. Space-time plots of the sediment bed height fluctuation with respect to the instantaneous stream-wise average, defined as h (cid:48) b ( x, t ) = h b ( x, t ) − (cid:104) h b (cid:105) x ( t ), are shown in figure 4. The emergence ofdune-like patterns can be clearly observed, with a streamwise succession of alternating ridges andtroughs forming shortly after start-up in all three cases. In the two simulations in the laminarregime we obtain similar fluctuation amplitudes. However, the propagation velocity is signifi-cantly larger in case LC2 (i.e. at larger Shields number) than in case LC1. The turbulent caseTO1, on the other hand, is found to exhibit a higher growth rate, rapidly leading to enhancedfluctuation amplitudes as compared to both laminar cases. Furthermore, these space-time plotsshow the occasional occurrence of dune mergers with a subsequent increase of wavelength and an y / D ( a ) x/D ( b ) x/D Figure 5: Close-up of the instantaneous, spanwise-averaged solid-volume fraction (cid:104) φ p (cid:105) z (plottedin greyscale) as well as the streamlines computed from a spanwise-averaged instantaneous flowfield in the final phase of the simulated interval: ( a ) case LC1; ( b ) case TO1.6 a ) λ a v / D t u b /H ( b ) (cid:104) h (cid:48) b h (cid:48) b (cid:105) / x / D t u b /H Figure 6: ( a ) Time evolution of the mean wavelength of the sediment bed height normalized withthe particle diameter. The dashed lines indicate the wavelengths of the second to sixth streamwiseharmonics in the current domain. ( b ) Time evolution of the r.m.s. sediment bed height. Thedashed line shows the fit obtained by Langlois & Valance (2007, figure 6 a ) for their case with Re H = 15130, Θ = 0 . H/D = 35. Note that this fit was obtained for t ≥ u b /H (their firstdata point is indicated by the symbol ‘ (cid:66) ’). In both graphs solid lines with the following colorscorrespond to the present cases: LC1, LC2, TO1.apparent decrease of the propagation speed. For times t (cid:38) H/u b the sediment bed patterns inthe turbulent case TO1 (cf. figure 4 c ) remain roughly invariant with two distinct dunes featuringsomewhat different elevation amplitudes.A visualization of the fluid-bed interface and the streamlines of the spanwise averaged flow fieldtowards the end of the simulated intervals is shown in figure 5. It is found that the patterns inthe laminar cases indeed correspond to ‘small dunes’ in the terminology of Ouriemi et al. (2009),and to ‘vortex dunes’ with significant separation on the lee-side in the turbulent case (the graphfor case LC2 is similar to case LC1 and has been omitted). These results are, therefore, consistentwith the regime classification based upon the channel (or pipe) Reynolds number proposed bythese authors (cf. figure 2).The instantaneous two-point correlation of the bed height fluctuation as a function of stream-wise separations r x , defined as R h ( r x , t ) = (cid:104) h (cid:48) b ( x, t ) h (cid:48) b ( x + r x , t ) (cid:105) x , exhibits a clear negative min-imum in all of the present cases (figure omitted). Therefore, we can define an average patternwavelength λ av as twice the streamwise separation at which the global minimum of R h occurs. Thetime evolution of the mean wavelength λ av , normalized by the particle diameter, is shown in fig-ure 6( a ). Also indicated with horizontal dashed lines in the graph are the wavelengths of the secondto sixth harmonics in the present computational domain (recall that L x /D = 307 . λ av in the turbulent case TO1exhibits several fast oscillations between L x / L x / L x / λ av /D ≈ tu b /H ≈
250 theaverage wavelength then grows at an increasing rate, settling again at λ av ≈ L x / . D untilthe end of the simulated interval. Contrarily, the two laminar cases have a less oscillatory initialevolution. Case LC2 first settles into a plateau-like state (with λ av ≈ L x /
3) after approximately250 elapsed bulk time units. Subsequently the wavelength corresponding to the second harmonic( λ av ≈ L x /
2) grows in amplitude and becomes dominant after approximately 750 bulk units.Case LC1 does not appear to settle into any of the harmonic wavelengths of the current domain,exhibiting an average wavelength in the range of 100 −
150 particle diameters for the most partof the simulated interval. In turbulent channel flow experiments, Langlois & Valance (2007) havedetermined values of the initial pattern wavelength of λ/D ≈ − − . H , except for a few initial bulk7 a ) (cid:104) q p (cid:105) x / q v ( b ) (cid:104) q p (cid:105) x / q i t u b /H ( c ) (cid:104) q p (cid:105) x t / q i (cid:104) q p (cid:105) x t / q v , −1 −3 −2 −1 Θ Figure 7: Analysis of the volumetric particle flow rate (per unit spanwise length), q p ( x, t ). ( a ) Thetime evolution of the streamwise average value, (cid:104) q p (cid:105) x ( t ) in cases LC1 (black color) and LC2(red color) is shown with solid lines. The viscous scale q v = Ga ν is used for the purpose ofnormalization. Additionally, the two dashed lines in each case indicate the respective minimumand maximum values (min x q p , max x q p ). ( b ) The same quantity in case TO1, using the inertialscale q i = u g D for normalization. ( c ) The value of the time average (cid:104) q p (cid:105) xt over the final part ofthe simulations plotted versus the Shields number Θ: • , LC1; • , LC2; • , TO1. Note the differentscalings ( q v in the laminar cases, q i in the turbulent case). The open circles are for featurelessbedload transport in laminar flow (Kidanemariam & Uhlmann, 2014); the dashed line is the fit (cid:104) q p (cid:105) xt /q v = 1 .
66 Θ . from that reference. The solid line is the Wong & Parker (2006) version (cid:104) q p (cid:105) xt /q i = 4 .
93 (Θ − . . of the Meyer-Peter & M¨uller (1948) formula for turbulent flow.time units. This range is comparable to the range found for the initial wavelength of the ‘smalldunes’ in pipe flow reported as 2 . − . H by Ouriemi et al. (2009, figure 7 b ) and for the ‘vortexdune’ data shown in their figure 3( c ), where an initial wavelength of 4 H is observed.Let us turn to the amplitude of the sediment patterns. The evolution of the r.m.s. value of thefluid-bed interface location is plotted in figure 6( b ). The fact that no saturation is observed by theend of the simulated intervals is consistent with experimental observations, where it was found thateven after an order of magnitude longer times the amplitude of the patterns continues to grow (theintervals simulated in the current work correspond to roughly one minute in the experiments ofLanglois & Valance, 2007, conducted over more than one hour). Both of the present laminar casesexhibit growth at an approximately constant rate (with a slope of 6 · − in the units of figure 6 b ).Contrarily, the time evolution in the turbulent case TO1 is quasi-linear with different slopes indifferent time intervals. The initial growth of the turbulent case (for times up to tu b /H ≈ (cid:46) tu b /H (cid:46)
500 are approximately ten times higher than thegrowth in the laminar cases, while in the remaining two intervals the growth rate is comparable tothe laminar value. As can be seen in figure 6( b ), the time evolution in the final period of case TO1is roughly equivalent to the one determined by Langlois & Valance (2007) in turbulent flow atcomparable parameter values ( Re H = 15130, Θ = 0 . H/D = 35).The propagation speed of the patterns can be determined from the shift of the maximum ofthe two-point/two-time correlation of the fluid-bed interface fluctuation h (cid:48) b ( x, t ). It turns out thatthe patterns in case LC1 propagate at a relatively constant speed of approximately 0 . u b , whilethe propagation velocity decreases with time during the coarsening process in cases LC2 and TO1,reaching values of 0 . u b and 0 . u b , respectively, in the final period of the current simulations.The latter number is consistent with the range of values reported for ‘vortex dunes’ by Ouriemi et al. (2009, figure 3 b ).The volumetric particle flow rate (per unit spanwise length), q p ( x, t ), is analyzed in figure 7.The solid lines in figure 7( a, b ) show the temporal evolution of the streamwise average (cid:104) q p (cid:105) x ( t ),8hich is observed to reach approximately constant values after a few hundred bulk time units inall cases. The continuous growth of the pattern amplitudes (cf. figure 6 b ) seems to have only amild influence upon the total particle flow rate, irrespective of the flow regime. These graphs alsoshow for each instant the maximum and minimum values (in space) of the particle flow rate, drawnas dashed lines. Although these extrema curves are noisier, it can be observed that the maximacontinue to grow until the end of the simulations, consistent with the increase in the amplitude ofthe propagating patterns. Of particular interest in view of applications is the scaling of the particletransport rates, typically expressed as a function of the Shields number Θ. Figure 7( c ) shows thespace-averaged values, additionally averaged in time over the final part of the simulations, denotedas (cid:104) q p (cid:105) xt . It is found that the present values of (cid:104) q p (cid:105) xt /q v (where q v = Ga ν ) in the two laminarcases are only slightly below the (approximately) cubic power law fitted by Kidanemariam &Uhlmann (2014) to their simulation data for featureless bedload transport. It is obviously notpossible to infer scaling from two data points. However, if a power law of the particle flow rate asa function of the Shields number is assumed, the present laminar data suggest a cubic variation.Turning to the turbulent case TO1, figure 7( c ) shows that the value for (cid:104) q p (cid:105) xt /q i (with the inertialscaling q i = u g D ) is very close to the value given by the empirical law of Wong & Parker (2006),which in turn is a modified version of the Meyer-Peter & M¨uller (1948) formula for turbulent flow.Wong & Parker’s formula is valid for plane sediment beds. The fact that the present data agreeswell with that prediction together with the observed mild variation in time (cf. figure 7 b ) showsthat the presence of ‘vortex dunes’ does not strongly affect the net particle transport rate. We have performed direct numerical simulation of the flow over an erodible bed of sphericalsediment particles above both the threshold for particle mobility and for pattern formation. Twocases in laminar flow (with different Galileo and Shields numbers) lead to the formation of ‘smalldunes’, while one case under turbulent flow conditions exhibits ‘vortex dunes’, consistently withthe regime classification of Ouriemi et al. (2009). The reconstruction of the fluid-bed interfacefrom a spanwise-averaged solid volume fraction (involving a threshold value) has allowed us toanalyze the length scales, amplitude and propagation velocity of the sediment patterns in detail.In all three respects, the results of the present simulations are found to be consistent with availableexperimental data.We have observed that the continuing growth of the dune patterns, which have not reached astatistically stationary state after approximately 1000 bulk time units, does not strongly affect thenet volumetric particle transport rates. In the two laminar cases the particle flow rate (per unitspan) is consistent with a cubic power law as a function of the Shields number; these values arefound to be not far from those obtained in featureless bedload transport. The value pertaining tothe turbulent case is very well predicted by the transport law of Wong & Parker (2006) which isderived for turbulent flow in the presence of a plane mobile bed. The present results therefore seemto suggest that the presence of ‘small dunes’ as well as that of ‘vortex dunes’ up to the amplitudesencountered in the present simulations does not lead to a modification of the net particle transportrate which would require a correction of the respective transport laws. This conclusion should bereassessed in the future when much longer temporal intervals can be covered.The present work demonstrates that the DNS-DEM approach to sediment pattern formation isfeasible today. Although still costly in terms of computational resources, it is already possible toaddress some of the outstanding questions with this method. Some aspects which are of importancein geophysical applications (such as reaching the fully rough turbulent regime, guaranteeing anasymptotically large computational domain and integrating over asymptotically long temporalintervals) still present a considerable computational challenge.As a next step, the streamwise length of the computational domain should be extended in order Note that in (Kidanemariam & Uhlmann, 2014) the Shields number (termed Θ
Pois therein) was defined basedupon the assumption of a parabolic fluid velocity profile for consistency with the reference experiment. As a result,the fit represents even the data points at larger values of the Shields number Θ
Pois with good accuracy.
9o reduce the influence of the discreteness of the numerical harmonics upon the pattern wavelength.Conversely, shrinking the box length will allow to reveal the smallest amplified wavelength of thesedimentary patterns. Finally, an in-depth investigation of the flow field which develops overthe time-dependent sediment bed can be carried out based upon the simulation data. Preliminaryvisualization suggests that in the turbulent case the coherent structures leave their footprint in thebed shape, visible as longitudinal ridges and troughs superposed on the roughly two-dimensionaldune patterns. Such an analysis is left for future work.This work was supported by the German Research Foundation (DFG) through grant UH 242/2-1. The computer resources, technical expertise and assistance provided by the staff at LRZM¨unchen (grant pr58do) are thankfully acknowledged.Movies of the particle motion are available as ancillary files along with the arXiv submissionand from the following URL:
References
Betat, A., Kruelle, C. A., Frette, V. & Rehberg, I.
Eur. Phys. J. E. Soft Matter (5), 465–76. Chan-Braun, C., Garc´ıa-Villalba, M. & Uhlmann, M.
J. Fluid Mech. , 441–474.
Charru, F.
Phys.Fluids (12), 121508. Charru, F. & Hinch, E. J.
J. Fluid Mech. , 111–121.
Charru, F. & Mouilleron-Arnould, H.
J. Fluid Mech. , 303–323.
Coleman, S. E., Fedele, J. J. & Garcia, M. H.
J. Hydraul. Eng. (December), 956–965.
Coleman, S. E. & Melville, B. W.
J. Hydraul. Eng. (4),544–560.
Coleman, S. E. & Nikora, V. I.
Water Resour. Res. (4), n/a–n/a. Colombini, M.
J. Fluid Mech. ,1–16.
Colombini, M. & Stocchino, a.
J. Fluid Mech. ,121–131.
Engelund, F. & Fredsoe, J.
Annu. Rev. Fluid Mech. (1),13–37. Garc´ıa-Villalba, M., Kidanemariam, A. G & Uhlmann, M.
Int. J. Multiph. Flow , 54–74. Kennedy, J. F.
J. FluidMech. (4), 521–544. 10 idanemariam, A.G. & Uhlmann, M. Int. J. Multiphase Flow (submitted) . Kidanemariam, A. G., Chan-Braun, C., Doychev, T. & Uhlmann, M.
New J. Phys. (2), 025031. Langlois, V. & Valance, A.
Eur. Phys. J. E. Soft Matter (3), 201–8. Meyer-Peter, E. & M¨uller, R.
Proc. 2nd Meeting ,pp. 39–64. IAHR, Stockholm, Sweden.
Ouriemi, M., Aussillous, P. & Guazzelli, ´E.
J. Fluid Mech. , 295–319.
Raudkivi, A. J.
J. Hydraul. Eng. (1), 58–64.
Richards, K. J.
J. Fluid Mech. (3), 597–618. Sumer, B. M. & Bakioglu, M.
J. FluidMech. , 177–190.
Uhlmann, M.
J. Comput. Phys. (2), 448–476.
Uhlmann, M.
Phys. Fluids (5), 053305. Uhlmann, M. & Duˇsek, J.
Int. J. Multiphase Flow , 221–243. Wong, M. & Parker, G.
J. Hydr. Eng.132