Large eddy simulation of water flow over series of dunes

Water Science and Engineering, 2011, 4(4): 421-430doi: 0.3882/.issn.1674-2370.2011.04.006http://www.waterjoumal.cne-mail: wse2008@vip.163.comLarge eddy simulation of water flow over series of dunesJun LU'2, Ling-ling WANG*, Hai zHU', Hui-chao DAI'1. College of Water Conservancy and Hydropower Engineering, Hohai Universit, Nanjing 210098, P. R. China2. Zhangjiagang Water Conservancy Bureau, Zhangjiagang 215600, P: R. China .Abstract: Large eddy simulation was used to investigate the spatial development of open channelflow over a series of dunes. The three-dimensional filtered Navier-Stokes (N-S) equations werenumerically solved with the fractional-step method in sigma coordinates. The subgrid-scaleturbulent stress was modeled with a dynamic coherent eddy viscosity model proposed by the authors.The computed velocity profiles are in good agreement with the available experimental results. Themean velocity and the turbulent Reynolds stress afected by a series of dune-shaped structures werecompared and analyzed. The variation of turbulence statistics along the flow direction affected bythe wavy bottom roughness has been studied. The turbulent boundary layer in a complex geographicenvironment can be simulated well with the proposed large edy simulation (LES) model.Key words: large eddy simulation (LES); dunes; turbulent boundary layer;flow separation1 IntroductionIn sandy rivers, dunes and ripples are the most common river bed structures. Flowseparation and recirculation are enbanced because of the presence of dunes in river beds,which alter the overall flow resistance and, consequently, affect the transport of sediment andcontaminants within the river. Flow separation begins at the dune crest. On the lee side ofdunes, there exists a series of recirculation eddies, and on the opposite stoss side, that is, onthe upstream of the crest, the flow reattaches (Fourniostis et al. 2009).The turbulent open channel flow over a river bottom with dunes has been experimentallyand numerically studied under the assumption of a fixed bed without sediment movement. Lyn(1993) reported experimental research on the mean flow and turbulent characteristics overartificial space-periodic one-dimensional river bed structures using a laser Dopplervelocimetry (LDV). Wiberg and Nelson (1992) conducted experiments under unidirectionalflow over much smaller asymmetric and symmetric riverbed structures, including high-angleand low-angle ripples. Similar experiments have been carried out by many other researcherssuch as Nelson and Smith (1993) and Ojha and Mazumder (2008). Recently, numerical studiesThis work was supported by the Nationar Natural Science Foundation of China (Grant No. 51179058), theNational Science Fund for Distinguished Young Scholars (Grants No. 51125034 and 50925932), the SpecialFund for Public Welfare of the Water Resources Ministry of China (Grant No. 201201017), and the 111Project (Grant No. B12032).*Corresponding author (e-mail: wanglingling@hhu. edu.cn). Received Jan. 7, 2011; accepted Aug. 20, 2011中国煤化工MHCNMHGon the developing turbulent flow using a solution of the Reynolds-averaged Navier-Stokes.equations (RANS) have been carried out by researchers (Huai et al. 2010; Peric et al. 1988).For example, Huai et al. (2010) used the realizable k-ε model to simulate the buoyant walljet, and Peric et al. (1988) numerically simulated the flow over a typical dune using the k-&turbulence model. The influence of sand grain roughness was taken into account with the wallfunction approach. Comparisons with experiments show that the computed separated andreattached flow over a dune is in good agreement with the experimental data. Similar studieswere reported by Mendoza and Shen (1990), who presented an algebraic stress model withwall functions in place of the k-ε model. They were able to obtain quite realistic predictionsof detailed pressure, velocity, and turbulence profiles. Johns et al. (1993) employed theone-equation turbulence model, with the turbulence model length scale prescribed from anempirical correlation. The comparisons show that the near-wall velocity and turbulence data,particularly the wall shear stress, are not in good agreement with experimental data. Yoon et al.(1995) simulated the flow over a fixed dune using the k- -0 turbulence model of Wilcox(1993). Comparisons of model prediction results with measured velocity and turbulence fields,as well as the pressure and friction distributions along the dune, show good agreement. Lu andWang (2009) compared three different k-ε models, the standard k-ε model, therenormalization group (RNG) model, and the V2F model, for simulating the separated flowpassing sills. The results show that, of the three models, the performance of the V2F model isthe most encouraging. The commercial computational fluid dynamics (CFD) code FLUENTwas employed by Fourmiostis et al. (2009) to simulate sub critical, turbulent, and open-channelflows over a bottom with five dunes using the free surface treatment method based on therigid-lid approximation. They reported that the numerical prediction results of the meanvelocity and turbulence are in good agreement with available experimental data.Although the RANS model can provide satisfying prediction results of flow features overdunes, it cannot calculate the power spectrum of hydrodynamic turbulence, which is veryimportant for turbulence research and a rather interesting parameter for some specificengineering applications and refined turbulence structures. In recent years, large eddysimulation (LES) has been used to study hydrodynamic turbulence. For example, the separatedflow passing sills have been studied by Lu and Wang (2008). The computed results show thatLES is very powerful and encouraging. Therefore, it is possible to use the LES technique toobtain a detailed numerical simulation of open channel flow over a series of dunes.2 Governing equations and numerical methods2.1 Governing equations in σ coordinates and turbulence modelThe LES approach in σ coordinates was used by Lin and Li (2002). The equations forthe large-scale motion can be obtained by integrating a spatial filter with the Navier-Stokesequations (indicated by an over-bar). With a top-hat filter based on the Boussinesq assumptionJun LU et al. Water Science and Engineering, Dec. 2011, Vol. 4, No. 4, 421-430中国煤化工MHCNMHGand the principle of chain differentiation, the governing equations for the incompressible fuidin general coordinates can be written as品(可)=0(1)a(云)("1)=5s-132Otpag(r5o)+vg(r5g )+ 2(5)(2)where t is time, 5* and ζ' are the coordinate directions in the transformed space,5 =a5*/ax,，x is the Cartesian coordinate, J-' is the Jacobin factor of the .transformation, u, is the Cartesian component of the velocity field, f; is the gravitationalforce, P is the water density, p is the static pressure, v is the kinematic viscosity, and动is the subgrid stress.In natural open channel flow, the free surface elevation varies with time and the bottom isuneven. These cause certain difficulties when conventional Cartesian coordinates are used inthe disretization of the domain in the vertical direction. To solve the uneven physical domain,the vertical σ -coordinate transformation is used as follows:ζ=x=x, 52=x=y, ζ°=x=σ=z+h(3)h+nwhere η is the surface elevation, z is the vertical direction in Cartesian coordinates, andh isthe static water depth. The subgrid stess云can be decomposed into the sum of atrace-free factor 司and a diagonal tensor动:动-δ。些=-2vS,(4)1{u ou,(5)2(ax, ax,where动/3 is absorbed in the pressure term and δ; is the Kronecker sign. In this study,we modeled the eddy viscosity v, with the coherent eddy model (Lu and Wang 2009):v, =(C,4)|sa|(6)|so|=a25,5; +(1-a)2Qwhere 2Q=2y -Sy , and 2y=(u/ax,-qxu/ax)/2; a is the weight factor; s is alength scale defined hereas (4 +4 +4})" (Bardina et al.1980), where 4, 4,and sare the control volume dimensions in the x, y, and z directions, respectively;and C, isa coefficient.The coffcient C. is computed by a dynamic procedure. The initial dynamic constant,C, is calculated as follows (Germano et al. 1991; Lilly 1992):C=_LyMj(8)MyMyJun LU et al. Water Science and Engineering, Dec. 2011, Vol. 4, No. 4, 421-430423中国煤化工MYHCNMHGLy=渝一公司,(9)M, =-(可)|同|5 +893,(10)where the test filter width 0= k'd , and the coefficient k' can be optimized. The coefficientk' is taken tobe √5 , and a is taken to be 0.5. C? can exhibit long- time negative valueswhich will generate numerical instability. To solve this problem, we first calculated the valueof C2 with the coherent eddy model proposed by Lu'et al. (2010), and then filtered it inspace with a box filter as follows:c: =C(x,1)= [G(x-5)(x,)d5(11)where G(x-5) is a smooth function and C (x,t) is the dynamic constant calculated by .Eq. (8). Although the smooth function can be in many forms, a box filter function was used forconvenience in this study. Therefore, Eq. (11) can be written ascq=C=( LyM,(12)MM,JIn fact, we can repeat the averaging procedure above to make the model coefficientsmoother. In addition, the condition v, +v≥0 was imposed. This condition ensured that thetotal resolved dissipation remained positive or at zero.2.2 Numerical methods and boundary conditionsThe splitting operator approach was used to numerically solve the govening equations.The momentum equations were split into three steps in each time interval: advection, diffusion,and pressure propagation. The advection step was solved using a combination of the quadraticbackward characteristic method and the Lax-Wendroff method. The central difference methodwas used to solve the diffusion step. The pressure propagation step was used to solve pressureind gravitational forces. In order to satisfy the divergence-free condition as imposed by thecontinuity equation, the projection method was employed to calculate the pressure andvelocity fields to obtain the updated velocity field. The conjugate gradient stabilized(CGSTAB) method was used to solve the above equations.The goverming equations may be solved only when adequate boundary conditions areprovided. Several types of boundary conditions are usually imposed in open channel flowproblems. The no-slip boundary condition is imposed on the bottom wall and a zero gradientboundary condition is imposed on the two side walls. At the inflow boundary, the inflow ratewith a predetermined velocity distribution with added Gaussian distribution random signals isspecified and the gradient of the water surface elevation is assumed to be zero. At the outflowboundary, a convective boundary condition is imposed. A Lagrange Euler method is used tolocate the free surface elevation. Details of these boundary conditions can be found in Lin andLi (2002), Li and Ma (2003), and Lu and Wang (2009).424Jun LU et al. Water Science and Engineering, Dec. 2011, Vol. 4, No. 4, 421-430中国煤化工MYHCNMHG.3 Study caseFig.1 is a typical domain of open channel flow over a series of twelve identical dunes.The numbers and shape of the dunes shown here correspond to those used in the experimentsof Ojba and Mazumder (2008). The experimental channel has a length of 10 m, a width of0.5m,andadepthof0.5m.TheduneshaveameanlengthofL=32cmandameanheightofH = 3 cm at the crest The angles of the stoss side and lee side slopes of the dunes were 6° and50°, respectively (Fig. 1). Velocity profiles were measured by an acoustic Doppler velocimeter(ADV). The mean flow depth h is 30 cm, the mean horizontal velocity U is 0.5 m/s, and thedischarge Q is 0.04 m/s. The corresponding Reynolds number is Re =Uh/v=1.5x10' andthe Froude numberis Fr =u/Vgh =0.29 .I Flow@甲⑨⑧回②④回⑧@°②4➊◎③❺◎❺③⑥囱回❹卤③30mm- 320 mmFig. 1 Sketch of computed domain and dune profiles of Ojha and Mazumder (2008)The simulation domain was carefully chosen in order to properly set up the inflow andoutflow bouindaries. To ensure that the inlet flow fully developed, the length of the mainchannel upstream of the first dune was extended to 0.5 m. The length of the main channeldownstream, of the last dune was also extended to 0.5 m to avoid the influence caused by thedownstream outflow. A non-uniform grid of451 x 11 x 45 nodes in the x, y, and z directions,respectively, was used to discretize the computational domain (5 m long, 0.1 m wide, and 0.3m deep). The expansion ratio of the grid did not exceed 1.01. The time step was 0.0002 s. Thegrid and time step were small enough to obtain grid-convergent results. The total computingtime was 40s.4 Results and discussionFig. 2 shows the computed horizontal mean velocity profiles at the trough and crestpoints, respectively, of each dune. Corresponding profiles of the vertical mean velocity areshown in Fig. 3. In the figures, odd and even numbers inside circles represent the positions ofvelocity profiles at the trough and crest points, respectively (Fig. 1). From Fig. 2 and Fig.3, itcan be seen that the computed velocity results agree well with the experimental data. The mostevident feature is that the horizontal mean velocity near the trough points is negative except atthe trough of the first dune, which means that reversal flow exists at the trough points. TheJun LU et al. Water Science and Engineering, Dec. 2011, Vol. 4, No. 4, 421-430425中国煤化工MHCNMH G.computed results clearly show that the turbulent boundary layer is developing from the firstdune to the seventh dune and it reaches a quasi-steady state after the seventh dune.-- Computed。 Experimental60cm/s②①⑥⑨回⑧四皿⑧⑧②④.1.0p0.6 ti 0.4+).2 t010 12 1401214x/(a)At trough points(b) At crest pointsFig. 2 Horizontal velocity U over dunes (dashed lines are where U= 0)一Computed 。 Experimentalo①⑤①⑨四田四D⑨四四四0[0..8↑: 0..6.” 0..4 t)22 14(田) At trough pointsFig, 3 Vertical velocity W over dunes (dashed lines are where w= 0)To keep the figure readable, the velocity field of two adjacent dunes is shown as in Fig. 4.As can be seen, the spatial development of the turbulent open channel flow over the dunes hasthe characteristic that the flow separates at the dune crest and reattaches on the stoss side ofthe next dune. A similar phenomenon was shown by Foumiostis et al. (2009) in theircomputed results of five dunes. In Fig. 5, a closer view of the velocity vector between theninth and tenth dunes in the fully developed region is shown to provide a better view of theflow field over a complete trough region. The figure contains both the experimental andcorresponding computational velocity vector fields. Agreement is seen to be good within therecirculation zone. According to the flow pattemn, three distinct layers are observed: (1) theinternal layer at z/h S 0.05 , where the mean horizontal velocity always stays negative; (2) theadvecting and diffusing layer at z/h S 0.15, where the mean streamwise velocity becomespositive but is still much smaller than that in turbulent channel flow without dunes at the sameposition; and (3) the outer flow layer at z/h> 0.15, where the mean streamwise velocity isalmost the same as in turbulent channel flow with a flat bottom (Ojha and Mazumder 2008).426JunLU et al. Water Science and Engineering. Dec. 2011, Vol. 4, No.4, 421-430中国煤化工MYHCNMHG.. 40 cm/s0.30.0。业0.5 0.6 0.7 0.8 09 1.0 1.1 1.2 1.3x(m)Fig. 4 Computed velocity field一40cm/s.3 12.2t89 9919.2939.4 9.5 9.6 9.7 98899919.2939.49.59.69.798(a) Experimental results(b) Computational resultsFig. 5 Velocity field between ninth and tenth dunesFig. 6 shows the spatial development of the Reynolds stress -(u'w) along the flowover dunes. From Fig. 6, it can be seen that the flow characteristics vary up to the seventhdune, beyond which the entrance effect disappears. That is to say, the turbulence fullydevelops, because the increment of Reynolds stress intensity (the thickness in Fig. 6) growsmuch smaller after the seventh dune. Again, this is in qualitative agreement with theexperimental observations.The reattachment point is determined as the location closest to the bottom at which themean velocity changes sign. The predicted values at each reattachment point are shown inFig. 7. From Fig. 7, it can be seen that the first two computed reattachment lengths are muchlarger than the remaining values. Flow reaches the fully developed state after the seventh dune,which is in agreement with the result mentioned above. The predicted reattachment length isX,/Hs=6.0 at the quasi-steady state (as shown in Fig. 5), which is a lttle larger than theexperimental value of 5.8 (Fourniostis et al. 2009), where X, is the reattachment length.Kasagi and Matsunaga (1995) conducted an experiment of the flow over a backward-facingstep with a Reynolds number of 5540, giving a value of 6.5. With the increasing of theReynolds number in the range of 1 200 < Re < 6600 , the reattachment position tends to moveupstream.Dejoan and Leschziner (2004) reported a computed value of about 7.0. It seemsthat the experimental result underestimates the location of the reattachment point for purelytwo- dimensional conditions.Jun LU et al. Water Science and Engineering, Dec. 2011, Vol. 4, No. 4, 421-430427中国煤化工MHCNMH G.Reynolds shear stress (m2/s2)-1.3.21012246一1012x/Fig. 6 Distribution of Reynolds shear stressFig. 7 Rcattachment length at each duneThe free surface plays a significant role in the open channel over dunes. The normalizedfree surface level (n/h ) over dunes is shown in Fig. 8. From the figure, it can be seen that theaverage water level rises before the first dune because the existence of the dunes increases theroughness of the channel bottom. In the middle of the channel, free surface elevationfluctuates in accordance with the distribution of dunes on the bottom. The average water levelis low due to the reduction of the crosS-section area. Near the outflow region where dunesdisappear, the free surface tends to increase just like the outlet of a sluice.0.04 [.03 t.02 t0.01 t-001-0.02 512 14 16Fig. 8 Free surface elevation over dunesIn order to analyze the head loss caused by dunes, the energy equation of thecross-section between the inflow and outflow boundaries can be expressed asz. + fi+ qm.... + Pou + CaouCou+h.(13)pg 2ggwhere Zm is the inflow water elevation; Zom is the outlow water elevation; am and aouare cofficients, where ci.=1 and 0%ou=1; Pm and Pout are the dynamic water pressuresat the inlet and outlet positions, respectively; U;n and Uow are the mean stream-wisevelocities at the inlet and outlet positions, respectively; g is the gravitational acceleration;and hw is the water head loss.From Eq. (13), we find that the head loss between the inflow and outflow sections isabout 0.014 m.5 ConclusionsLarge eddy simulation of open channel flow over a series of dunes in the sigmacoordinates were carried out. The subgrid stress was modeled with the dynamic coherent eddy428Jun LU et al. Water Science and Engineering, Dec. 2011, Vol. 4, No. 4, 421-430中国煤化工MHCNMHGmodel proposed by the authors. The computed velocity profiles are in good agreement with theavailable experimental data. It was found that turbulence does not reach a fully developedstate until at the seventh dune. The computed results show that the length of the separationzone at the fully developed region is about 6Hg, which is a lttle larger than those in theexperiment. The spatially mean free surface level decreases in the flow direction. The headloss caused by dunes is about 0.014 m. With the employment of the sigma coordinates andspltting operator method, the dynamic coherent eddy model proposed by the authors has greatadvantages in simulating turbulent flows in complex geographic environments, which is ofgreat importance and engineering value to river dynamics simulation.ReferencesBardina, J, Ferziger, J. H, and Reynolds, W. C. 1980. Improved subgrid-scale model for large eddysimulation. Proceedings of the 13th Fluid and Plasma Dynamics Conference. Snowmass: AmericanInstitute of Aeronautics and Astronautics.Dejoan, A, and Leschziner, M. A. 2004. Large eddy simulation of periodically perturbed separated flow overa backward-facing step. Fluid and Heat Flow, 25(4), 581-592. [doi:10. 1016/j.jheatfluidflow.2004.Fourniostis, N. T, Toleris, N. E., and Demetracopoulos, A. C.2009. Numerical computation of turbulencedevelopment in flow over sand dunes. Advances in Water Resources and Hydraulic Engineering,Proceedings of 16th IAHR-APD and 3rd IAHR-ISHS, 943-848. Beijing: Tsinghua University Press.[doi:10.1007/978-3-540-89465-0_ 148]Gemano, M, Piomelli, U, Moin, P, and Cabot, W. H.1991. A dynamic subgrid-scale eddy viscosity model.Physics of Fluids, 3(7), 1760-1765. [doi:10.1063/1.857955]Huai, W. X, Sheng, Y. P, and Komatsu, T. 2003. Hybrid finite analytic solutions of shallow water circulation.Applied Mathematics and Mechanics (English Edition), 24(9), 1081-1088. [doi: 1000-0887(2003)09-0956-07]Huai, W. X., Li, Z. W, Qian, Z. D., Zeng, Y. H, Han, J., and Peng, W. Q.2010. Numerical simulation ofhorizontal buoyant wall jet. Journal of Hydrodynamics, 22(1), 58-65. [doi:10.1016/S1001-6058(09)60028-7]Johns, B., Soulsby, R. L.. and Xing, J. 1993. A comparison of numerical model experiments of free surfaceflow over topography with flume and ficld obscrvations. Journal of Hydraulic Research, 31(2), 215-228.[doi:10.1080/00221689309498846]Kasagi, N, and Matsunaga, A. 1995. Three- dimensional particle tracking velocimetry measurement ofturbulence statistics and energy budget in a backward-facing step flow. Fluid and Heat Flow, 16(6),477-485. [doi:10.1016/0142-727X(95)00041-N]Li, C. w., and Ma, F. x.2003. Large eddy simulation of difusion of a buoyancy source in ambient water.Applied Mathematical Modeling, 27(8), 649-663. [doi:10.1016/S0307-904X(03)00073-8]Lilly, D. K.1992. A proposed modification of the Germano subgrid-scale closure method. Physics of Fluids,4(3), 633-635. [doi:10.1063/1.858280]Lin, P. Z, and Li, C. W. 2002. A σ-coordinate three -dimensional numerical model for surface wavepropagation. International Journal for Numerical Methods in Fluids, 38(11), 1045-1068. [doi: 10.1002/1d.258]Lu, J, and Wang, L. L. 2008. Numerical study of large eddy structures-separated flows passing sls. Advancesin Water Resources and Hydraulic Engineering, Proceedings of 16th IAHR-APD and 3rd IAHR-ISHS,1795-1799. Bejjing: Tsinghua University Press. [doi:10.1007/978-3-540-89465-0_ 309]Lu, J, and Wang, L. L.2009. Comparison of several turbulent models for calculating separated flows passingon sill. Advances in Water Science, 20(2), 255 260. (in Chinese)Jun LU et al. Water Science and Engineering. Dec. 2011, Vol.4, No. 4, 421-430429中国煤化工MHCNMH G.Lu, J,, Tang, H. W, and Wang, L. L. 2010. A novel dynamic eddy model and its application to LES ofturbulent jet with free surface. Science in China, Ser. G 53(9), 1671-1680. [doi: 10.1007/s11433-010-4077-z]Lyn, D. A. 1993. Turbulence measurement in open channel flows over artificial bedforms. Journal ofHydraulic Engineering, 119(3), 306-326. [doi:10.1061/(ASCE)0733-94291993)119:3(306)]Mendoza, C, and Shen, H. W. 1990. Investigation of turbulent flow over dunes. Journal of HydraulicEngineering, 116(4), 459-477. [doi:10. 1061/(ASCE)0733-9429(1990)1 16:4(459)]Nelson, J. M., and Smith, J. D.1993. Mean flow and turbulence over two-dimensional bed forms. WaterResources Research, 29(12), 3925-3953. [oi:10.1029/93WR01932]Ojha, s. P., and Mazumder, B. S. 2008. Turbulence characteristics of flow region over a series of 2-D duneshaped structures. Advance in Water Resources, 31(3), 561-576. [doi:10.1016/jadvwatres.2007.12.001]Peric, M, Ruger, M., and Scheuerer, G 1988. Calculation of the Two-dimensional Turbulent Flow over a SandDune Model. Erlangen: University of Erlangen.Wiberg, P. L, and Nelson, J. M. 1992. Unidirectional flow over asymmetric and symmetric ripples.. Geophysica! Research, 97(8), 12745- 12761. [doi:10.1029/92JC01228]Wilcox, D. C.1993. Turbulence Modeling. for CFD. La Canada: DCW Industries.Yoon, J. Y, Patel, V. C, and Ettema, R. 1995. Numerical model of flow in ice-covered channels. Journal ofHydraulic Engineering, 122(1), 19-26. [doi;10.1061/(ASCE)0733-9429(1996)122:1(19)]430Jun LU et al. Water Science and Engineering, Dec. 2011, Vol. 4, No. 4, 421430中国煤化工MHCNMH G.