- Email: [email protected]

S0045-7930(15)00023-7 http://dx.doi.org/10.1016/j.compfluid.2015.01.014 CAF 2793

To appear in:

Computers & Fluids

Received Date: Revised Date: Accepted Date:

12 July 2013 28 December 2014 27 January 2015

Please cite this article as: Nguyen, Y.Q., Wells, J.C., Modeling bedform development under turbulent flows using Large-Eddy-Simulation and Immersed-Boundary-Method, Computers & Fluids (2015), doi: http://dx.doi.org/ 10.1016/j.compfluid.2015.01.014

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Modeling bedform development under turbulent ows using Large-Eddy-Simulation and Immersed-Boundary-Method

Y Q. Nguyen , J.C. Wells 1,

Graduate School of Science and Engineering, Ritsumeikan University 1-1-1 Nojihigashi, Kusatsu, Shiga, Japan

Abstract

A numerical model was built to study the mechanism of sedimentary bedform development in hydraulically smooth turbulent ows. The model consisted of a module for ows, a module for sediment transport, and a module for bed surface evolution. The ow was unsteady, three-dimensional and modeled by a Large-Eddy-Simulation (LES) method coupled with an ImmersedBoundary-Method (IBM). Governing equations were discretized on a xed Cartesian grid by Finite Dierence Method. Sediment (bedload) transport was estimated by Van Rijn formula corresponding to bed shear stress distribution obtained from the ow solution. The bed surface evolution was adapted to the sediment ux and described by the Exner-Polya equation. Updated bed surface was then used as the boundary for solving the ow eld in the next time step. Time-advancement was discretised by the AdamsBashforth method. The model was rst validated by two test cases of bed

Corresponding author

(Y Q. Nguyen), 1 Present address: Dept. of Civ. Eng., HCMC University of Technology, Vietnam.

Email addresses: [email protected] [email protected] (J.C. Wells)

Preprint submitted to Computers and Fluids

February 4, 2015

shear stress and turbulence statistics over xed sinusoidal bed surfaces. It was then employed to study initiation and development of bedforms from an initially at bed to fully developed forms. Formation of newly and successively formed bedforms, growing and downstream propagating process of existing bedforms were very close to experimental observations. Instantaneous bed shear stress and corresponding sediment ux around evolving bedforms, which were dicult to observe in experiments, were also well produced by this model. Keywords: sediment transport, bedforms, turbulent ows, Large-Eddy-Simulation, Immersed-Boundary-Method

2

1. Introduction

In hydraulics, sediment transport and sedimentary bedform are of fundamental problems. They are known to strongly aect dynamics of rivers, such as erosion and deposition of river bed and banks, as well as transport capacity of the ow. Particularly, the bedform increases ow resistance, hence, in uences living environment of aquatic creatures and enhances hydraulic disasters, such as draughts or oods. The bedform developing in hydraulically smooth turbulent ow is named ripple whose development and mature dimensions are independent of the

ow depth [1, 2, 3, 4]. Ripple is observed to be initiated from small local disturbances on an initially at bed [4, 5, 6, 7]. The disturbances are about a few grain diameters [6, 7] and formed by random pile-ups of the sediment on the bed [6] or by random actions of high turbulent velocities [7, 8]. They continue to grow as more sediment is trapped, quickly become two-dimensional and disturb the ow eld [6, 7]. As a result, downstream bed shear stress distribution is caused to deviate from standard values in a region about 1100 the disturbance height [4, 7]. Bed surface is then deformed accordingly and the rst sand wave appears on the bed surface. The wave continues to grow as more sediment is deposited at its downstream side and, at a critical height, a ow separation becomes visible. Successive waves appear one by one downstream and form a ripple train [7]. While growing, the sandwaves are propagating downstream at speed inversely proportional to their heights [2]. This yields the `coalescence' process where the stronger and faster sand waves overtake the weaker ones [2, 7, 9]. After the coalescence process, existing system may be two- or three3

dimensional depending on the applied shear stress compared to the critical one [9]. At higher bed shear stress, the sandwaves are more two-dimensional. A coalesced wave is observed to keep growing until it reaches a mature dimension at which the sandwave simply propagates downstream without growing further [2, 3, 9]. In hydraulically smooth turbulent ows with the particle Reynolds number, de ned based on the grain diameter and friction velocity, less than 2.5, the ripple length is about from 700 wall units to 2500 wall units [2, 3] while the model ripple height is about 38 wall units [10]. There have been many attempts to study the process of bedform development numerically. However, just few of them focused on the hydraulically smooth region. Kennedy [11] is considered as the rst one working on this subject with a two-dimensional potential ow over an erodible bed. Later, Richards [12] added viscous eects to the ow model with an one-dimensional turbulence model for ows in hydraulically rough regions. This work was then extended to the hydraulically smooth region by Sumer & Bakioglu [13, 14] for analyzing ripple formations. In this work, ripple development from an initially at bed in a hydraulically smooth turbulent ow is modeled with a highly accurate solution for unsteady three-dimensional ows by Large-Eddy-Simulation (LES) method. Methods with Reynolds-Averaged Navier-Stokes equations (RANS) may be preferred in most engineering applications but in this case their performance is not satisfactory in solving ow separations behind the sandwaves which are critical to formation of the bedform [15, 16, 17], hence are not employed here. Direct-Numerical-Simulation (DNS) method can also be a good candidate, but requires more computational cost than LES [17]. 4

In addition, as the sandwaves grow, the bed surface changes in time. A curvilinear body- tted grid requires regridding, associating with interpolations of the ow eld, over time. To facilitate this problem, ImmersedBoundary-Method (IBM) [18, 19] is coupled into LES, hence bedforms are allowed to freely develop in a xed Cartesian grid system. In our previous short papers [20, 21], we reported preliminary results of bedform development simulations obtained from the model. We used a preliminary version of the model to examine initiation of a bedform on an initially at bed [20]. The model was then revised to be able to simulate bedforms up to mature state [21]. However, in those two papers, details of the model were not presented. In this paper, we add details of the model and discuss about its assumptions and limitations. More details of the equations and user-de ned parameters, such as the maximum bed slope and the adaption length, as well as eects of their selected values on the results are presented. In addition, the validation tests have been revised for better matching of boundary conditions and for more grid resolutions. We also discuss more about eects of the ow separations behind the bedforms on their development and mature dimensions.

5

2. Numerical models

To build the computational model, the following assumptions were used, following suggestions in the literature.

In hydraulically smooth ows, bed roughness and dynamics of individual grains have negligible eects on the ow eld [3, 4, 13, 14]. Therefore, in this study, the bed surface was treated as a smooth, continuum one. Eects of the free surface were ignored because development and mature dimensions of the bedform in the hydraulically smooth ows are independent of the ow depth and ow surface [1, 2, 3]. Accordingly, the free surface was treated as a xed free-slip one. Only bedload was considered for sediment transport as it was observed to be dominant in the process of ripple formation [5, 6, 9, 10]. The local bedload ux does not adapt instantaneously to the local bed shear stress, hence there is a lag time, or equivalently, a lag distance between the two. This lag distance is shown to play a central role on formations of the sand waves [11, 22]. Time scale of ow development is much smaller than that of the bedform development [12]. Accordingly, bed surface was treated as a xed one while hydrodynamic equations were being solved.

Based on the assumptions, a numerical model was built with the following components whose computational procedure is shown in Figure 1: 6

For the ow eld: the unsteady, three-dimensional (3D) Navier-Stokes equation was solved with a Large-Eddy-Simulation (LES) method coupled with an Immersed-Boundary-Method (IBM). For sediment transport: the bedload formula proposed by Van Rijn [23], which is valid for low particle Reynolds numbers, was employed. For bed elevation evolution: the continuity equation, Exner-Polya formula [4] which relates the gradient of the bedload ux and time variation of the bed surface, was selected.

The sediment transport and bed elevation were evaluated in two-dimensional (2D) space. They were assumed to change in time and only in the longitudinal direction while homogeneous on the spanwise direction. A 2D model for the bed evolution facilitates the computational considerably, but should be sucient to describe important characteristics of the bedform development process. It is reported that for very ne sands with very low particle Reynolds numbers and high bed shear stress, ripples are almost two-dimensional [2, 9]. For the assumption of 2D bedforms to be validated in this model, the particle Reynolds number was assumed to be less than 1.0. 2.1. Module of hydrodynamic equations

The governing equations for LES coupled with IBM for viscous, incompressible uids become [18, 24]: Buj 0 Bxj

7

(1)

Bui u Bui 1 Bp Bij B ui f Bt j Bxj Bxi Bxj Bxj Bxj i

(2) where: i,j=1,2,3 indicating streamwise, vertical, and spanwise directions, respectively; indicates ltered quantities of the velocity ui and the pressure p; and are the uid density and viscosity, respectively; fi is the arti cial body force in the i direction accounting for the solid portion; ij is the sub-grid scale (SGS) stress, which was computed by the Shear-Improved Smagorinsky eddy viscosity model [25]. This SGS stress model was selected for its simplicity, gentle restriction on time step, moderate computational cost, and performance comparable to the dynamic Smagorinsky model for this type of ow [25]. 2

ij u iuj i uj u

1 3

ij ij kk 2T Sij T

where

pCs q p|S| |xSy|q 2

(3) (4) (5)

1 p Bui Buj q 2 Bxj Bxi is strain-rate tensor. |S| 2Sij Sij is its magnitude. In Eq. (5), Cs equals 0:16, as recommended for this type of ow [25]. x y is the ensemble average which was evaluated along the homogeneous direction, i.e. the spanwise, in this model. For IBM, the direct forcing method [18, 24] was employed. In this method, the arti cial body force fi was implemented by applying no-slip condition on the solid surface. Velocities at cells inside but nearest to the solid surface were Sij

8

interpolated from the solved velocities at cells outside but nearest to the solid surface to satisfy zero-velocity on the solid boundary. Linear interpolation was used in this code. 2.2. Module of sediment transport

Once the 3D ow eld was solved, bed shear stress was calculated from velocity gradient along the bed surface and then averaged in the spanwise, z, direction. Equilibrium bedload ux was then estimated from the Van Rijn formula [23]: qb rps 1qgs : d : 0 5

where

1 5

0:053 Td :

:

2 1

0 3

p s 1qg { d d

1 3

2

is the particle parameter;

u u u

(6) (7)

(8) is the transport stage parameter; s is the speci c density of the particle (s{); g is the gravitational acceleration; d is the particle diameter; u and u are the bed-shear velocity and the critical bed-shear velocity for incipient bedload transport, respectively. The bedload ux was then modi ed due to eects of non-zero bedslope as [22]: T

2

2

c

2

c

c

B h qeq qb 1:0 Bx 1

in which h is the local bed elevation.

9

(9)

Fredsoe & Deigaard [26] showed that the bedslope also aects the critical condition for bedload motion as:

h{Bx 1 Btan

(10) with u : critical friction velocity for zero-slope bed which, in this study, was determined from White's experiments [27]; s: friction angle whose mean value is about 32 degrees with tan s 0:63 [12]. uc uc0

1

s

c0

2.3. Model for bed surface evolution

[4]:

Evolution of local bed surface was described by the Exner-Polya equation p1 nq BBt BBq

(11) where and are the local tangential and normal directions of the bed surface, respectively; q is the local bedload transport rate, and n is the porosity of the bed material. Since in this model, the bed surface was assumed as a continuum, n 0. The bedload ux q in (11) was related to the qeq calculated by Eq. (9) as [12, 22]: where

Bq q qeq B leq leq 3 d d: T

(12)

(13) is the adaption length, or the lag distance, taken as the average saltation step length proposed by Van Rijn [23]. 0 6

10

:

0 9

Giri & Shimizu [16], simulating sand-dunes in hydraulically rough ows and high particle Reynolds numbers, used an adaption length of order of 100 particle diameters. In this study, with low particle Reynolds numbers, leq was in order of 10 particle diameters. Tests with various adaption lengths, from 5 to 20 particle diameters, were conducted for the same other computational conditions. The dierences in the development paths and developed dimensions of the sandwaves among them were negligible. When sand-waves have reached certain heights, if the value of Bh{Bx is equal to or greater than tan s, Eq. (10) will be invalid. To prevent such irregularity and, also physic-contradicted, behavior, where the absolute value of the slope of the wave is greater than the friction angle s, the slope of the sandwave is xed to a certain value f ixed, which is slightly smaller than s. Accordingly, the following two constraints have to be satis ed: 1

$ &

Bh{Bx f ixed % ³ h:dx const: 1

1

(14)

The second constraint, together with employment of the Exner-Polya continuity equation of bedload ux, ensures conservation of the sediment mass when the slope is adjusted to f ixed. Tests with various values of f ixed which were around 20 degrees to 30 degrees show that they had negligible eects on bedform formations (see Figure 9). f ixed 30 was selected in our simulations. The above equations were discretized by a fourth order central nite difference in space and second order Adams-Bashforth method in time [28] on a xed staggered Cartesian grid. In the streamwise and spanwise directions, 11

periodic conditions were used. In the vertical direction, uniform grid spacing of about 0.9 wall unit was applied up to x {H (or y{H ) =0.2, where H is the ow depth, and then a hyper-tangential grid distribution was used. In the streamwise and spanwise directions, uniform grid spacings of 8.8 wall units were used. Compared to the considered range of the particle Reynolds number of less than 1.0, these grid resolutions oered ne and converged ow and sediment solutions. For stable simulations, the nondimensional time step t~ ptu {H q was 1:10 . According average Courant number [28] was less than 0.04. As shown in the computational procedure, Figure 1, from initial conditions of ow eld and sediment properties, Eq. (1) and (2) were solved. The bed shear stress was then obtained from the velocity gradient. The bedload

ux was estimated by Eq. (6) to (10). The bed surface was updated by Eq. (11) to (14). Over the updated bed surface, the hydrodynamic module was solved in an interval T N t, which was also the timestep for the bed evolution. Employment of T was to ensure that the ow had enough time to adapt to the new bed pro le. Values of N r10; 100s were con rmed to have minor eects on the results. 2

2

4

0

For simplicity, afterward, , , and are used to indicate streamwise, vertical, and spanwise directions, respectively. 2

x

y

z

12

3. Validation tests 3.1. Shear stress distribution over a sinusoidal bed surface

The hydrodynamic module was validated with a test case of bed shear stress distribution along a sinusoidal bed surface. The results were compared to those conducted with Direct-Numerical-Simulation (DNS) method by De Angelis et al [29] on an body- tted grid. The computational conditions are shown in Table 1. The computed results with three grid resolutions are presented in Figure 2, together with the DNS results. The LES+IBM did not perform satisfactorily with the coarse grid resolution of 64 72 64. The results of LES+IBM with the two ner grids are seen to converge and to agree well with the DNS ones. Compared to the body- tted DNS, the LES+IBM uses higher grid resolutions. This is mainly because the IBM requires ne grid spacings near the solid surface for implementation of arti cial force. Balaras [19] discussed about this issue and argued that the IBM still has more ecient CPU-time per node. 3.2. Mean ow eld over a sinusoidal bed surface

In this test case, the present code was tested for the ability of properly solving the mean outer ow eld. The test case was a ow over a sinusoidal bed surface, which was conducted by Ohta et al [30] with a DNS body tted method. The computational conditions are presented in Table 2. The Reynolds number is 300 based on the mean friction velocity corresponding to the total drag, u , Reu , and is 4300 based on the bulk velocity Ub, ReU . Figure 3 shows the statistics of the mean streamwise velocity, vertical velocity, turbulent kinetic energy, and the Reynolds stress at four dier0

0

b

13

ent positions along a wavelength, together with the DNS results. Excellent agreements between them can be seen here.

14

4. Bedforms evolving from an initially at bed

Initial conditions (inputs) for the numerical model were: (a) domain size, computational grid, and at bed level. (b) ow conditions: mean ow Reynolds number, Re H u H { , (c) sediment conditions: mean initial particle Reynolds number, Rep u d{ and mean initial Shields number u {ps 1qgd. Figure 4 shows an example of bedform development from an initially

at bed surface with H 300, Rep 0:1, and 1:9. Bed pro le is plotted at dierent nondimensional time t~, up to t~ 360. The initial at bed surface was introduced with a small disturbance of a sinusoidal wave with an amplitude less than one wall unit. From this disturbance, a small bedform appeared [20] and triggered successive visible bedforms as observed in this gure. The successive bedforms were not initiated simultaneously but appeared one by one. An existing bedform triggered initiation of a new one downstream. The newly formed bedform grew in time. When it gained a critical height, it then again triggered another new one downstream. The process was repeated. After a certain interval, a chain of bedforms, which was identical to the ripple train observed in experiments [7], was seen on the bed surface at t~ r30 80s in Figure 4. Once a bedform appeared, it already had a certain height and length which continued increasing in time to fully developed dimensions. At this time, the bedform simply migrated downstream without growing further. Examples of developed bedforms are the rst and the second ones from the right at t~ 360 in Figure 4. 0

2

0

0

0

0

15

0

Figure 5 shows variation of the bulk velocity normalized by the initial mean friction velocity, hU i {u , and the Reynolds number Re hU i H { , as the bedform developed, for the test case presented in Figure 4. Associating with the bedforms, form drag increased in time, while the driving pressure gradient was constant. As a result, the bulk velocity decreased in time. Up to t~ 100, the bulk velocity decreased rapidly as during this time, the bedforms developed rapidly. Later, the dropping rate was steady as the bedforms already lled up the domain. At t~ 320, the bulk velocity dropped to almost 60% of the initial one. Such drops also observed in the experiments by Coleman et al [5]. The process of initiation and migration of the ripple train is now examined in two categories: (a) How a new bedform is initiated, and (b) How an existing bedform migrates downstream. Figure 6 shows an example of instantaneous ow eld and bed shear stress, which were averaged in the spanwise direction, around a newly formed bedform. A new bedform is identi ed when a new front is formed downstream of an existing bedform. First, the existing upstream bedform caused behind it a separation region associated with a strong positive gradient of bed shear stress and corresponding bedload ux. Accordingly, the bed was eroded. Further downstream, the bed shear stress reached a peak and then experienced a negative gradient before returning to the base value. With the negative gradient of the bed shear stress, the sediment was deposited and formed a bump. The more the sediment was eroded, the deeper that area was; hence 0

16

the stronger the gradient of the bed shear stress was. Therefore, the bump was continuously fed and grew gradually. When the bump gained a certain threshold height, its downstream part turned into a clear front whose slope was still smaller than that of the xed value, f ixed. Examination of the bed shear stress showed that this happened after the bed shear stress, or the Shields number, , behind the bump was smaller than the critical one, c, for the bedload transport. Once c, the bedload ux was zero; hence there was discontinuity of the sediment transport on that area. Because the scoured sediment was deposited right before that discontinuity point, it helped the bump to become a visible front. At this instant, the bed shear stress downstream of the front also dropped to zero. A new separation was going to be created associated with strong positive gradient of the bed shear stress behind this point. The bed was then scoured more. The front became more visualized. Tests with dierent Rep1s and 1s showed that the height of such a newly formed bedform was u { r12:0 13:0s, thus was approximately equal to the thickness of the viscous sublayer, 11:6 [31]. Therefore, a newly formed bedform was just extruded above this layer and disrupted its continuity. The bedform now started directly aecting the turbulent core

ow. It then created a separation region [3]. Figure 7 helps to clarify the question of how an existing bedform grew and propagated downstream. Similar to the observations for the newly formed bedforms, separations before and after the bedform produced strong variation of the bed shear stress. Adapting to the bed shear stress, the bedload ux at the two ends of the bedform, where was less than c, was zero; hence the 0

0

17

sediment transport was limited within the surface of the bedform. Adapting to the gradient of the bed shear stress, the sediment was scoured over the upstream lee, particularly just behind the reattachment point, where the gradient was positive, and deposited over the downstream lee, where the gradient was negative. This cyclic process allowed the bedform to grow and to propagate downstream. Phase lag between the bed shear stress and the bed surface also helped the bedform gain height. The peak of the bed shear stress distribution was at a position upstream of the peak of the bed pro le. Downstream of the peak of the bed shear stress, the bed shear stress gradient turned negative and allowed the sediment to deposit from top of the bedform afterward. This deposition made the bedform height increase. Figure 8 shows dimensions of a mature bedform. The height of a fully developed bedform was about r45 50s wall units. The length was about r950 1000s wall units. These dimensions agreed well with the experimental ones [2, 3, 10]. The dimensions were almost invariant with various values of Rep 's and 's. For a Rep , higher only shortened the development time of the bedform to mature dimensions. With the mature height around r45 50s wall units, the sandwave was completely immersed in the viscous wall region [31]. While the sediment was transported only in the bedload mode in which the sediment was moved in a thin layer on the bed surface by the bed shear stress, the dominance of the turbulence above this region might prohibit the wave from growing further and forced it to translate downstream at this equilibrium height. Once the wave-height was xed, the wavelength could not increase further. First, the 0

0

18

amount of supplied sediment was already con ned within a single sandwave. Second, once the upstream sandwave reached the equilibrium dimensions, its downstream ow separation no longer grew; hence the adapting bed shear stress and bedload ux only uctuated around almost xed values. The sandwave then gradually reached the equilibrium or mature dimensions from this point. To check eects of selected f ixed in Eq. (14) on bedform formations, tests with two values of f ixed of 20 degrees and 25 degrees had been conducted. The bedform developments are plotted in Figure 9. Because the constraints in Eq. (14) were applied only when the slope of the downstream side of a bedform became higher than the friction angle s, bedforms obtained with two values of f ixed were almost the same up to t~ 40:0. Later, the bedforms diered on the downstream slopes, which were equal to the f ixed . Their developed lengths, heights, and downstream migrating speeds were identical to each other. Other tests with f ixed up to 30 degrees also yielded similar observations. Therefore, we had selected f ixed 30 degrees in our simulations.

19

5. Summaries

Development of sedimentary bedforms from an initially at bed in an hydraulically smooth turbulent ow was studied numerically. The results showed that the rst bedform was triggered by gradient of bed shear stress and sediment ux around an disturbance on the bed. Successive bedforms were then formed downstream of the rst one. Once a bedform reached a certain height, sediment transport was con ned within its surface and helped the bedform gain height and length. The dimensions of a mature sandwave was about [45-50] wall units in height and [950-1000] wall units in length, which agreed well with experimental observations in the literature.

20

Acknowledgements

The authors thank for support from Japanese Government Scholarship (MEXT). Valuable comments from the two anonymous reviewers are greatly appreciated. References

[1] J.R.L. Allen. Current ripples. North-Holland Publishing Co., 1968. [2] Raudkivi A.J. Ripples on stream bed. Journal of Hydraulic Engineering, 123(1):58{63, 1997. [3] M.S. Yalin. On the determination of ripple length. Journal of Hydraulics Division, 103(HY4):439{442, 1977. [4] M.S. Yalin. tion, 1977.

Mechanics of Sediment Transport.

Pergamon Press, 2 edi-

[5] S.E. Coleman, J.J. Fedele, and M.H. Garcia. Closed-conduit bedform initiation and development. Journal of Hydraulic Engineering, 129(12):956{965, 2003. [6] S.E. Coleman and B.W. Melville. Initiation of bed forms on a at sand bed. Journal of Hydraulic Engineering, 122(6):301{310, 1996. [7] P.B. Williams and P.H. Kemp. Initiation of ripples on at sediment beds. Journal of the Hydraulics Division, Proceedings of the ASCE 97(HY4):505{522, 1971. 21

[8] A. Gyr and A. Schmid. The dierent ripple formation mechanism. Journal of Hydraulic Research, 27(1):61{74, 1989. [9] P.A. Mantz. Bedforms produced by ne, cohesionless, granular and

akey sediments under subcritical water ows. Sedimentology, 25:83{ 103, 1978. [10] P.A. Mantz. Cohesionless, ne-sediment bed forms in shallow ows. Journal of Hydraulic Engineering, 118(5):743{764, 1992. [11] J.F. Kennedy. The mechanics of dunes and antidunes in erodible-bed channels. Journal of Fluid Mechanics, 16:521{544, 1963. [12] K.J. Richards. The formation of ripples and dunes on an erodible bed. Journal of Fluid Mechanics, 99(3):597{618, 1980. [13] B.M. Sumer, M. Bakioglu, and A. Bulutoglu. Ripple formation on a bed of ne, cohesionless, granular sediment. In Proceedings of Euromech 156: Mechanics of Sediment Transport, pages 99{110, 1982. [14] B.M. Sumer. On the formation of ripples on an erodible bed. of Fluid Mechanics, 144:177{190, 1984.

Journal

[15] Y.S. Chang and A. Scotti. Modelling unsteady ows over ripples: Reynolds-averaged Navier-Stokes equations (RANS) versus large-eddy simulation (LES). Journal of Geophysical Research, 109:C09012, 2004. [16] S. Giri and Y. Shimizu. Numerical computation of sand dune migration with free surface ow. Water Resources Research, 42:W10422, 2006. 22

[17] C.J. Keylock, R.J. Hardy, D.R. Parsons, R.I. Ferguson, S.N. Lane, and K.S. Richards. The theoretical foundations and potential for large-eddy simulation (LES) in uvial geomorphic and sedimentological research. Earth-Science Reviews, 71:271{304, 2005. [18] R. Verzicco, J. Mohd-Yusof, P. Orlandi, and D. Haworth. Large eddy simulation in complex geometric con gurations using boundary body forces. AIAA Journal, 38(3):427{433, 2000. [19] E. Balaras. Modeling complex boundaries using an external force eld on xed Cartesian grids in large-eddy simulations. Computers and Fluids, 33:375{404, 2004. [20] Q.Y Nguyen and J.C. Wells. Numerical modelling of bedform development; formation of sand-wavelets in hydrodynamically smooth ow over an erodible bed. Annual Journal of Hydraulic Engineering, JSCE, 52:163{168, 2008. [21] Q.Y Nguyen and J.C. Wells. A numerical model to study bedform development in hydraulically smooth turbulent ows. Annual Journal of Hydraulic Engineering, JSCE, 53:157{162, 2009. [22] M.D. Bui, T. Wenda, and W. Rodi. Numerical modeling of bed deformation in laboratory channels. Journal of Hydraulic Engineering, 130(9):894{904, 2004. [23] L.C. Van Rijn. Sediment transport, part i: Bed load transport. Journal of Hydraulic Engineering, 110(10):1431{1456, 1984. 23

[24] E.A. Fadlun, R. Verzicco, P. Orlandi, and J. Mohd-Yusof. Combined immersed-boundary nite-dierence methods for three-dimensional complex ow simulation. Journal of Computational Physics, 161:35{ 60, 2000. [25] E. Leveque, F. Toschi, L. Shao, and J.P. Bertoglio. Shear-improved Smagorinsky model for large-eddy simulation of wall-bounded turbulent

ows. Journal of Fluid Mechanics, 570:491{502, 2007. [26] J. Fredsoe and R. Deigaard. Mechanics of coastal sediment transport. Advanced series on ocean engineering 3. World Scienti c, 1994. [27] A.J. Raudkivi. Loose Boundary Hydraulics. Pergamon Press, 3 edition, 1990. [28] B.J. Geurts. Elements of Direct and Large-Eddy Simulation. Edwards, 2004. [29] V. De Angelis, P. Lombardi, and S. Banerjee. Direct numerical simulation of turbulent ow over a wavy wall. Physics of Fluid, 9(8):2429{2442, 1997. [30] T. Ohta, Y. Miyake, and T. Kajishima. Direct numerical simulation of turbulent ow in a wavy channel. JSME International Journal Ser. B, 41(2):447{453, 1998. [31] S.P. Pope. Turbulent ows. Cambridge University Press, 2000.

24

Table 1: Computational conditions of De Angelis et al's [29] test case. : ow depth; : bed wave-length; : bed wave-height; Re 0 u 0H { where u 0 is the mean friction velocity; Lx : computational domain in streamwise direction; Lz : computational domain in spanwise direction, Nx ; Ny ; Nz : grid sizes in the streamwise, vertical, and spanwise directions, respectively. H

DNS

LES+ IBM 1.04H 0.05H Re 170 Lx 6 4 Lz Nx { 21 64 or 128 or 256 Ny Nz 64 65 72 or 96 64 0

25

Table 2: Computational conditions of Ohta et al's [30] test case. Symbols are the same as in Table 1. ReU UbH { where Ub is the bulk velocity. b

DNS LES+ IBM 3.84(H ) 0.05(H ) Re 300 ReU 4300 Lx 4 1 Lz 0:5 Nx { 64 128 Ny Nz 128 64 72 64 0

b

26

Initial conditions

Large-Eddy-Simulation Immersed-Boundary-Method (three-dimensional)

Hydrodynamic equations

10 to 100 time steps

Calculate bed shear stress fixed bed

van Rijn (1984) Bedload formula (two-dimensional)

Calculate bedload flux

Exner-Polya equation (two-dimensional)

Update bed surface

END

Figure 1: Computational procedure.

27

LES+IBM (64 x 72 x 64) LES+IBM (128 x 72 x 64) LES+IBM (256 x 96 x 64) DNS De Angelis et al (1997)

4

τ₀/<τ₀>

3 2 1 0 -1

0

0.2

0.4

x/H

0.6

0.8

1

Figure 2: Shear stress distribution over a xed sinusoidal bed: LES+ IBM with three dierent grid resolutions vs. body- tted DNS results by De Angelis et al [29]. Black area is the bed surface.

28

y/ (Hα) ↵)

1 0.8 0.6

H H

0.4 0.2 22α↵

0 0

1

0

1

0

1

0

1

↵) α)

0.6

y/ (H

0.4 0.2 0

-0.04

0

-0.04

0

0

0.04

0

0.04

↵) α)

0.6

y/(H

0.8

0.4 0.2 0 0

0.01

0

0.01

0

0.01

0

0.01

2 k/U k/Ubb

1

↵) α)

0.6

y/ (H

0.8

0.4 0.2 0

0 0.004

0 0.004

0 0.004

0 0 hu v i /Ub 2

0 0.004

Figure 3: Statistical quantities of a ow over a sinusoidal bed. : LES+IBM (1287264), : DNS by Ohta et al (1998). From top to bottom: streamwise velocity

^1 , vertical velocity hvi u^2, turbulence kinetic energy k 1{2 uy hui u i u^iu^i, iu and Reynolds stress hu1v1i uz 2 u^1 u^2, where ^ indicates a time-averaged quan1u 29 tity. solid line

circle

x/H y0 H 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0 0.2 0.1 0

1.5

3

4.5

6

7.5

9

10.5

12

13.5

t˜ 360 340 320 300 280 260 240 220 200 180 160 140 120 100

80 60 50 40 30 20 10

0

Figure 4: Development of bedforms from an initially at bed, from bottom to top. The

ow is from left to right.

30

20 5000 16

12 3000

8 2000

4

0 0

Re = !U " H/ν

!U "/ u τ 0

4000

1000

0 50

100

150

200

250

300

350

t˜

Figure 5: Variation of the bulk velocity due to development of bedforms, for the test case in Figure 4.

31

0.5

1

1.5

2

2.5

3

3.5

4

0.1625

2

0.125

1

0.0875

0

0.1625

2

0.125

1

0.0875

0

0.1625

2

0.125

1

0.0875

0

y/H

0.05 0.2

−1 3

0.1625

2

0.125

1

0.0875

0

0.05 0.2

y/H

τ0 / ⟨τ0 ⟩

−1 3

τ0 / ⟨τ0 ⟩

y/H

0.05 0.2

τ0 / ⟨τ0 ⟩

−1 3

−1 3

0.1625

2

0.125

1

0.0875

0

0.05 0.2

τ0 / ⟨τ0 ⟩

y/H

0.05 0.2

y/H

4.5 3

τ0 / ⟨τ0 ⟩

0

−1 3

0.1625

2

0.125

1

0.0875

0

0.05 0

0.5

1

1.5

2

2.5

3

3.5

4

τ0 / ⟨τ0 ⟩

y/H

0.2

−1 4.5

x/H

Figure 6: Successive ow elds (arrows), bed shear stress distributions (circles), and bed pro les (solid line), with a time interval of T~ 1:5, from bottom to top.

32

3

τ0 / !τ0 "

2 1 0

−1

5

q˜/ !˜ q"

4 3 2 1 0

y/H

0.2 0.15 0.1 0.05 20.5

21

21.5

22

22.5

23

23.5

x/H

Figure 7: Instantaneous bedload ux, bed shear stress, and bed pro les around an evolving bedform at three successive time steps with a time interval of T~ 10. : bed shear stress; middle: nondimensional bedload ux, q~ q{rps 1qgs0:5 d1:5; bottom: bed pro le. Three lines indicate the distributions at each time step: at time t~, at time t~ T~, at time t~ 2T~. Top

33

y/H

1 0.75 0.5 0.25 0 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6

x/H

y+

60 40 20 0 0

200

400

600

800

1000

1200

1400

1600

1800

x+

Figure 8: Mature dimensions of a bedform, nondimensionalized by: the total ow depth, , in non-distored scales; viscous lengthscale, , with the vertical direction stretched for legibility. above

bellow

34

˜tt ~

0.2

60.0

0.1 0 0.2

25deg.

0.1

20deg.

50.0

0 0.2

40.0

y/H

0.1 0 0.2

30.0

0.1 0 0.2

20.0

0.1 0 0.2

10.0

0.1 0 0.2

0.0

0.1 0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6

6.5

7

x/H x/H

Figure 9: Development of bedforms with two values of f ixed: 20 degrees and 25 degrees.

35

Highlights

o We modeled bedform development in turbulent flows. o We used Large-Eddy-Simulation coupled with Immersed-BoundaryMethods. o Initiations, growing, and mature dimensions of bedforms were well reproduced.