Los movimientos de fondo representados por fronteras inmersas: experimentos numéricos usando las simulaciones de gran escala
José Eduardo Alamy Filho^{1} Harry Edmar Schulz^{2} André Luiz Andrade Simões^{3}
^{1}Faculty of Civil Engineering. Universidade Federal de Uberlândia. Uberlândia, MG, 38400902. Brazil. Email: zealamy@yahoo.com.br
^{2}Nucleus of Thermal Engg. and Fluids and Dept. of Hydraulics and Sanitary Engg. Universidade de São Paulo. São Carlos, SP, 13566590. Brazil. Email: heschulz@sc.usp.br
^{3}Dept. of Hydraulics and Sanitary Engg. Universidade de São Paulo. São Carlos, SP, 13566590. Brazil. Email: simoes@sc.usp.br
RESUMEN
Se utilizó el método de fronteras inmersas para evaluar el transporte de sedimentos en lechos deformables. Los procedimientos de la simulación de gran escala fueron utilizados para el tratamiento matemático de la turbulencia, y se utilizó la ecuación de adveccióndifusión para calcular la concentración de sedimentos. Se aplicó el método de diferencias finitas con malla desplazada para la solución numérica de las ecuaciones básicas (NavierStokes filtrado, continuidad y las ecuaciones de adveccióndifusión). Las derivadas espaciales fueron discretizadas mediante diferencias de segundo orden centrado. Se utilizó un esquema explícito de AdamsBashforth de segundo orden para la evolución del tiempo en la ecuación de adveccióndifusión, mientras que un esquema AdamsBashforth de cuarto orden fue utilizado para las ecuaciones de NavierStokes filtrado. La simulación numérica reproduce estructuras de flujo, como los remolinos grandes después de las crestas de las dunas y los vórtices contrarrotativos, que son importantes en el transporte de sedimentos. Los flujos de resuspensión y sedimentación (dependientes de la concentración de partículas) se calcularon utilizando las ecuaciones propuestas en este estudio. Las deformaciones del lecho por la erosión y la deposición pueden ser bien seguidas a través de los procedimientos aquí presentados, que muestran que esta metodología es adecuada para evaluar las modificaciones del lecho y el transporte de sedimentos en los flujos aluviales.
Palabras clave: Mecánica de fluidos computacional, método de fronteras inmersas, simulaciones de gran escala, transporte de sedimentos, lechos deformables.
ABSTRACT
The Immerse Boundary Method (IBM) was used to evaluate the sediment transport over deformable beds. Large Eddy Simulation (LES) procedures were used for the mathematical treatment of turbulence, and the advectiondiffusion equation was used to calculate sediment concentration. The Finite Differences Method with staggered grid was applied for the numerical solution of the governing equations (filtered NavierStokes, Continuity and advectiondiffusion equations). Spatial derivatives were discretized using second order centered differences. A second order explicit AdamsBashforth scheme was used for the time evolution in the advectiondiffusion equation, while a fourth order AdamsBashforth scheme was used for the filtered NavierStokes equations. The numerical simulation reproduced flow structures like large eddies after the dune crests and counterrotative vortices, which are important in sediment transport. Resuspension fluxes and sedimentation (dependent on particle concentration) were calculated using equations proposed in this study. The deformations of the bed caused by erosion and deposition may be well followed through the present procedures, showing that this methodology is adequate to evaluate bed modifications and sediment transport in alluvial flows.
Keywords: Computational fluid mechanics, immerse boundary method, large eddy simulation, sediment transport, deformable beds.
INTRODUCTION
Even if sediment transport is being studied over decades, its physical and mathematical description, in a more complete and useful form, is still a target focused by many researchers. Along these decades, different results were merged together to furnish the best description of sediment transport, involving conclusions for sedimentation velocity, general turbulent flows, formation and movement of eddies close to solid boundaries, among others ([79]; [18]).
Natural flows involving particle sedimentation and sediment transport are still too complex to be fully described by the available theoretical tools. Erosion and deposition imply in deformable boundaries, which quantification is still not completely understood. Numerical approximations thus seem to be the best way to predict the evolution of bed forms, but limitations also exist for such approximations, for example related to the computational effort. Ideally, Direct Numerical Simulation (DNS) using the flow and the sediment transport equations would solve the problem, but DNS is still impracticable for the usual scales and Reynolds numbers. Appropriate approximations may be obtained using adequate turbulence models for flow simulation (for example using Large Eddy Simulation, LES), together with adequate procedures to predict the bed deformation.
In the present context of sediment transport, two differential approaches are used in the literature: 1) the trajectories of the particles are followed using the so called "lagrangean models"; 2) the concentration of sediments in the flow is described along space and time. The lagrangean approach is still only interesting for low concentrations of sediments (volume fractions of the order of 10^{3}). However, natural flows involve higher volume fractions, which suggests the use of the advectiondiffusion equation for practical purposes. The continuous deformation of the boundaries of the flow imposes corrections in the domain of calculation after each time step of the simulation. This may also be treated in two ways: 1) the whole numerical grid is recalculated, adjusting it to the new frontiers (which may increase the computational costs), and 2) the numerical grid for the flow computations remains untouched, while a moving boundary, defined within the flow, is modified (Immerse Boundary Method, [12]).
In this study the bed deformations in sediment transport are calculated using the advectiondiffusion equations, the Immerse Boundary Method and Large Eddy Simulation.
MATHEMATICAL TREATMENT OF TURBULENCE
Water flows were simulated, allowing to assume incompressible flows. For LES, the filtered NavierStokes and continuity equations are:

(1) 

(2) 
and are the filtered velocity components and the pressure. i and j represent the directions in space. ? is the water viscosity and ?_{sg} is the subgrid viscosity. F_{i} is a body force that imposes the no slip condition at the immersed boundary. The Leonard and cross tensor were discarded due to the use of a second order advection scheme ([1516]). The subgrid viscosity was calculated here using the Second Order Velocity Structure Function (, which involves only velocity differences between adjacent cells), and is presented as:

(3) 
In equation (3), C_{k} is the Kolmogoroff constant, usually assuming the value 1.40 (for fully developed turbulence), ? is the grid dimension, and is calculated as:

(4) 
ADVECTIONDIFFUSION EQUATION
The advectiondiffusion equation applied to sediment concentration, and filtered accordingly to the LES procedures, assumes the form (without the cross and Leonard turbulent fluxes):

(5)

is the filtered concentration, D is the molecular diffusivity and D_{sg} is the subgrid diffusivity. w_{s} is the vertical sedimentation velocity. It is usually assumed that [1]:

(6) 
? is a factor which quantifies the influence of the particles in the equality. The turbulent Schmidt number, s, is sometimes assumed equal to 1.0 ([23], for example), but may also be a function of sediment and flow characteristics [5]. Equations (7a) and b are frequently used:

(7a) 

(7b) 
?_{s} is the sediment density, while u* is the shear velocity. The sedimentation velocity depends on the volume fraction of the sediment in the fluid, . For very dispersed suspensions, the Stokes law holds, or, expressing mathematically, , where d is the diameter of the particle. For concentrated suspensions, empirical studies as well as theoretical studies ([3]; [1314]; [2122], among others) present different equations for the sedimentation velocity. In this study the equation proposed by [14] was used, because it reproduces the form of other equations found in the literature through adequate small order truncations:

(8) 
w_{0} is the sedimentation velocity for ?=0 (Stokes sedimentation law). ?_{i} are constants, which were adjusted based on experimental data. A truncation on the third order term (n=3) allows to reproduce adequately experimental results for usual particle diameters (as those of reference [19], for example). For the set of data analysed by the authors, it was found that ?_{0}=35.123; [?_{1}=82.710; ?_{2}=52.646 and ?_{3}=2.071. Figure 1 shows the obtained adjustment.
Figure 1. Data of reference [19] adjusted to equation (8).
IMMERSED BOUNDARY METHOD
The immersed boundary method consists in using two different grids, a fixed and usually rectangular grid for the calculation domain (the so called eulerian grid), and a grid formed by a set of points following the form of the boundary surface of the flow (the so called lagrangean grid or immersed boundary). In this study the immersed boundary is located along the deformable bed of the flow. This grid deforms as consequence of erosion, sedimentation and resuspension. The no slip condition is imposed using the force field F_{i} shown in equation (1). Figure 2 shows both grids representing the flow domain and a wavy bed form.
Figure 2. The rectangle represents all the dominium of calculation (cartesian grid does not change along the time), with flows at both sides of the lagrangean grid which represents the bed form.
[10] and [11] present a description of the so called Virtual Physical Model used here to calculate the force on the lagrangean points. The force field F_{i} affects the fluid close to the immersed boundary (spreading o_{r}f the force) and is calculated using the function , given by equation (9).

(9) 
is obtained using the NavierStokes equations at the interfacial points, furnishing:

(10) 
represents the points of the eulerian grid. N is the number of dimensions of the problem. The function D of equation (9) governs the spreading of the force into the fluid. In this study, for a 3D flows, D was expressed by equations (11) and (12).

(11) 

(12) 
?x, ?y and ?z are the mesh sizes in the orthogonal directions and e is the normalized distance between a point in the boundary and a point in the flow.
EROSION AND DEPOSITION MODELS
The model described in [21] and [22] was used for the resuspension of sediments, and is shown in equation (13). Previous results using this equation are reported, for example, by [6], [20] and [23].

(13) 
In this equation 
and 


is the resuspension flux, u* is the shear velocity, and u*_{cr} is the critical value above which resuspension occurs. The critical value was obtained from the studies of [4], who suggested, based on experimental studies:

(14) 
Figure 3 shows the balance of mass for the cell just above the bed, in a 3D situation, used to calculate resuspension and deposition. represents advectives fluxes represents diffusive fluxes _{sed} represents the sedimentation flux and represents the resuspension flux, given by equation (13).
Figure 3. Sediment fluxes close to the bed which result in deposition or erosion.
The equation for the mass balance of Figure 3 is:

(15) 
Equation (15) allows to calculate the deformation of the bed, or in other words, the variation of the position of the lagrangean points, ?_{k}:

(16) 
?t is the time step used in the calculations, ? is the porosity of the bed. Positive and negative values of ?_{k} imply in deposition and erosion, respectively. The positions of the lagrangean points were corrected after each time step.
SIMULATION PARAMETERS AND PROCEDURES
A predefined initial wavy bed was used, and is shown in Figure 4. The dimensions of the bed undulations were based on observed values in a channel described by [17]. The eulerian grid was built with 121,500 points and the lagrangean grid with 11,041 points.
Figure 4. The "physical dominium" and the numerical dominium established for the calculations.
The upper and lower boundaries of the domain were defined as symmetry boundaries. The lateral walls were defined as rigid no slip boundaries. A constant velocity profile was imposed at the inlet; while normal derivatives of the velocity were set as zero at the outlet. The inlet horizontal velocity and a null transversal velocity were set as initial conditions for the velocity in all the domain. A null pressure was set as initial condition in all the domain. For the pressure corrections, zero normal derivatives were set at all boundary surfaces, with exception of the outlet, were a null pressure correction was set. The Reynolds number was 3,000 at the inlet section (H=0.0375m). For the concentration field, null fluxes were imposed at the upper and lateral surfaces. At the immersed boundary the conditions varied accordingly to the calculated resuspension fluxes. At the outlet, the normal derivatives were set as zero, and in the inlet a homogeneous sediment flux of 0.96 kg/m^{2}s was imposed. The characteristics of the sediment, needed to calculate the sedimentation velocity, were ?_{s}=2,650 kg/m^{3}and diameter of 0.12 mm (fine sand). The derivatives of velocities and concentrations were discretized using finite differences over a staggered grid (velocities stored at the faces and scalars stored at the center of the cells). Spatial derivatives were built using second order central differences. Considering the time evolution, a second order explicit AdamsBashforth scheme was used for the advectiondiffusion equation, while a fourth order AdamsBashforth scheme was used for the filtered NavierStokes equations. The time steps varied from the initial value of 10^{4}s to the final value of 10^{3}s, following a geometrical series with ratio of 1.05. Fractionated steps were used for the calculations. The computational code was built in Fortran Power Station 4.0. A simple PC computer, 3.2 GHz, was used for the calculations, without any parallel processing. The time needed for the calculations were approximately one day of calculation per second of flow.
RESULTS
Large flow structures and sediment transport
Figure 5 shows the vectors of the velocity field in the central longitudinal plane of the channel for the nondimensional times t^{+}=4.5, 9.0, 9.9 and 16.7, with t^{+}=Ut/H, U=8.0 cm/s and H=3,75 cm. It can be seen that the successive vortical structures indicated for t^{+}=4.5 suspend the sediment that, otherwise, would flow more close to the bed (as occurs for t^{+}=16.7, where the structures moved downward). The colored scale indicates the sediment concentration (see videos 1, 2 and 3).
Figure 5. Sediment concentration field [kg/m^{3}] and deposition on the channel bed. The vortical structures suspend the sediment. t^{+}=Ut/H, U=8.0 cm/s and H=3,75 cm.
The resuspension fluxes occur mainly right after the crests, where large vortical structures are formed. Figure 6 presents results of calculations made by [2] in which vorticity and concentration fields were compared, furnishing results comparable to the present study. The calculated twodimensional flow (Figure 6) presents recirculation patterns within the valleys of the wavy bed, and the detachment of large counter rotation structures from the crests. The concentration field is highly correlated with the vorticity field, and shows the relevance of the large flow structures for the sediment transport. A so called mushroom structure is also observed at the position x=1.2 m. This 2D experiment preceded the 3D experiment presented here.
Figure 6. Vorticity [1/s] and sediment concentration [kg/m^{3}] in a two dimensional flow, as calculated by [2]. t^{+}=Ut/H, where U=3.0 cm/s and H=28.0 cm.
Moving bed
The modifications in the bed form, produced by sedimentation and erosion, were reproduced by the evolution of the immersed boundary, without instabilizing the flow computations. Figures 7a trough k show the evolution of the bed form for the present calculations. Two sequential positions of the bed surface (immersed boundary) are superposed in each figure. The dark color is the "new" position of the immersed boundary, and shows the regions of deposition of sediments. The light color is the "old" position of the immersed boundary, and shows the regions of resuspension of sediments for the time interval considered. The sequence of figures allows to attest for the convenience of the methodology described here for the calculation of moving beds. Figure 7 shows the formation of ripples close to the channel inlet at the beginning of the experiment (Figures 7a and 7b), the predominant erosion of the crests (Figures 7b to e), the predominant deposition on the valleys (Figures 7b to e) and the general trend to elevate the bed position along all the channel, which is a consequence of the continuous input of sediment at the inlet section of the channel. For the conditions of this numerical experiment, the flow remained stable along all the calculations (see videos 4 and 5).
Figure 7. 3D evolution of the moving bed calculated with immersed boundaries, large eddy simulation and the advectiondiffusion equation.
CONCLUSIONS
The present study shows that 3D moving beds can be adequately calculated using immersed boundaries together with the advectiondiffusion equation for the evaluation of sediment transport and Large Eddy Simulation for the calculation of the turbulent flow. The joint use of these tools shows a high potential for practical applications. A predefined initial bed shape was modified along the calculations considering the results of the sediment fluxes for the first cell above the immersed boundary. The forces that guarantee the no slip condition at the immersed boundary were recalculated for each time step, accordingly the new position of the boundary determined by the mass fluxes. The turbulent flow was calculated using Large Eddy Simulation and a fixed eulerian grid for the flow domain during all the modifications suffered by the moving bed. As expected, erosion occurred mainly at the bed crests, while deposition occurred mainly at the bed valleys. As sediment was furnished continuously to the flow, the position of the bed was elevated along the calculations. The present numerical experiment was conducted in a 3.2 GHz personal computer, without any parallel processing. The numerical code was built in Fortran Power Station 4.0. For the numerical conditions of this experiment, one second of flow could be calculated in about one day.
ACKNOWLEDGEMENTS
The authors are indebted to CAPES, FAPESP, and CNPq, Brazilian research support foundations. The second author is indebted to CAPES, for the fellowship AEX N° 2201/062.
REFERENCES
[1] J.E. Alamy Filho. "Modelação numérica de processos de sedimentação em escoamentos turbulentos e análise de ressuspensão em canais". Doctoral thesis. School of Engineering at São Carlos. Universidade de São Paulo. Brasil. 2006.
[2] J.E. Alamy Filho and H.E. Schulz. "Modeling sediment transport for moving beds using the advectiondiffusion equation and the immersed boundary method". In: G.H. Merten, C. Poleto and A.L.O. Borges (Eds.). Selected papers of the VII ENESSediments: a multidisciplinary challenge. ABRH Associação Brasileira de Recursos Hídricos. Vol. 1, Chap. 4, pp. 127147. Porto Alegre, Brasil. 2007.
[3] T.E. Baldock, M.R. Tomkins, P. Nielsen and M.G. Hughes. "Settling velocity of sediments at high concentrations". Coastal Engineering. Vol. 51, pp. 91100. 2004. DOI: 10.1016/j,coastaleng.2003.12.004.
[4] A.R.J. Borges and J.A.G. Saraiva. "An erosion technique for assessing ground level winds". Wind Engineering. Oxford: Pergamon Press. Vol. 1, pp. 235242. 1980.
[5] N.L. Coleman. "Flume studies of the sediment transfer coefficient". Water Resources Research. Vol. 6, Issue 3, pp. 801809. 1970. DOI: 10.1029/WR006i003p00801.
[6] M. Garcia and G. Parker. "Entrainment of bed sediment into suspension". Journal of Hydraulic Engineering, ASCE. Vol. 117, Issue 4, pp. 414435. 1991.
[7] M.H. García, F. López and Y. Niño. "Characterization of nearbed coherent structures in turbulent open channel flow using synchronized highspeed video and hotfilm measurements". Experiments in Fluids. Vol. 19, pp. 1628. 1995. DOI: 10.1007/BF00192229.
[8] M. Garcia, Y. Niño and F. López. "Laboratory observations of particle entrainment into suspension by turbulent bursting". In: Ashworth, Bennett, Best y McLelland (Eds.). Coherent Flow Structures in Open Channels. Wiley & Sons. Chap. 3, pp. 6386. 1996.
[9] S.J. Kline, W.C. Reynolds, F.A. Schraub and P.W. Runstadler. "The structure of turbulent boundary layers". Journal of Fluid Mechanics. Vol. 12, pp. 741773. 1967.
[10] A.L.F. Lima e Silva. "Generating and implementing a new methodology for flow calculations over complex structures: immersed boundary method and virtual physical model" [in Portuguese]. Doctoral thesis. Universidade Federal de Uberlândia. Brasil. 2002.
[11] A.L.F. Lima e Silva, A. SilveiraNeto and J.J.R. Damasceno. "Numerical simulation of two dimensional flows over a circular cylinder using the immersed boundary method". Journal of Computational Physics. Vol. 189, pp. 351370. 2003. DOI: 10.1016/S00219991(03)002146.
[12] C.S. Peskin. "Numerical analysis of blood flow in the heart". Journal of Computational Physics. Vol. 25, pp. 220252. 1977. DOI: 10.1016/00219991(77)901000.
[13] J.F. Richardson and W.N. Zaki. "Sedimentation and fluidization: Part 1". Transactions of the Institution of Chemichal Engineers. Vol. 32, pp. 3553. 1954.
[14] H.E. Schulz and J.E. Alamy Filho. "Settling velocity for particle sets and its dependence on the sediment concentration  a study based on the kinetic energy quantification". Ciência & Engenharia (Science & Engineering Journal). Vol. 14, Issue 2, pp. 99106. 2005.
[15] S. Shaanan, J.H. Ferzinger and W.C. Reynolds. "Numerical simulation of turbulence in the presence of shear". Rep. TF6. Dept. of Mechanical Engineering. Stanford University. 1975.
[16] A. Silveira Neto, D. Grand, O. Métais and M. Lesieur. "A numerical investigation of the coherent structures of turbulence behind a backwardfacing step". Int. Journal of Fluid Mechanics. Vol. 256, pp. 125. 1993.
[17] L.B.S. Souza, H.E. Schulz, S.M. Villela, J.S. Gulliver and L.B.S. de Souza. "Experimental Study and numerical simulation of sediment transport in a shallow reservoir". Journal of Applied Fluid Mechanics. Vol. 3, Issue 2, pp. 921. 2010.
[18] A.J. Sutherland. "Proposed mechanism for sediment entrainment by turbulent flows". Journal of Geophysical Research. Vol. 72, pp. 61836194. 1967.
[19] M.R. Tomkins, T.E. Baldock and P. Nielsen. "Hindered settling velocity for sand". Abstract ICCEInternational Conference for Coastal Engineering. 2004.
[20] K. Trouw, J.J Williams and C.P. Rose. "Modelling sand resuspension by waves over a rippled bed". Estuarine, coastal and shelf science. Vol. 50, pp. 143151. 2000.
[21] L.C. van Rijn. "Sediment transport, Part I: bed load transport". Journal of Hydraulic Engineering, ASCE. Vol. 110, Issue 10, pp. 14311456. 1984.
[22] L.C. van Rijn. "Sediment transport, Part II: suspended load transport". Journal of Hydraulic Engineering, ASCE. Vol. 110, Issue 10, pp. 16131641. 1984.
[23] E.A. Zedler and R.L. Street. "Largeeddy simulation of sediment transport: currents over ripples". Journal of Hydraulic Engineering, ASCE. Vol. 127, Issue 6, pp. 444452. 2001.
Received: August 11, 2011 Accepted: September 3, 2012