Moving beds represented by immersed boundaries : numerical experiments using large eddy simulation

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 advection-diffusion 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 Navier-Stokes, Continuity and advection-diffusion equations). Spatial derivatives were discretized using second order centered differences. A second order explicit Adams-Bashforth scheme was used for the time evolution in the advection-diffusion equation, while a fourth order Adams-Bashforth scheme was used for the filtered Navier-Stokes equations. The numerical simulation reproduced flow structures like large eddies after the dune crests and counter-rotative 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 1 Faculty of Civil Engineering. Universidade Federal de Uberlândia. Uberlândia, MG, 38400-902. Brazil. E-mail: 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, 13566-590. Brazil. E-mail: heschulz@sc.usp.br 3 Dept. of Hydraulics and Sanitary Engg. Universidade de São Paulo. São Carlos, SP, 13566-590. Brazil. E-mail: simoes@sc.usp.br Alamy Filho, Edmar Schulz y Andrade Simões: Moving beds represented by immersed boundaries: numerical experiments...


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 ( [7][8][9]; [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 advection-diffusion 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 advection-diffusion 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 Navier-Stokes and continuity equations are: u i and p are the filtered velocity components and the pressure.i and j represent the directions in space.m is the water viscosity and m 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 ( [15][16]).The subgrid viscosity was calculated here using the Second Order Velocity Structure Function ( F 2 , which involves only velocity differences between adjacent cells), and is presented as: In equation ( 3), C k is the Kolmogoroff constant, usually assuming the value 1.40 (for fully developed turbulence), D is the grid dimension, and F 2 is calculated as:

ADVECTION-DIFFUSION EQUATION
The advection-diffusion equation applied to sediment concentration, and filtered accordingly to the LES procedures, assumes the form (without the cross and Leonard turbulent fluxes): c 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]: F 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: r 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, c / ϕ ρ = . 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]; [13][14]; [21][22], 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: w 0 is the sedimentation velocity for f=0 (Stokes sedimentation law).b 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 b 0 =35.123;b 1 =82.710;b 2 =52.646 and b 3 =2.071.Figure 1 shows the obtained adjustment.

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 of the force) and is calculated using the function f x t ( , ) k , given by equation (9).
k is obtained using the Navier-Stokes equations at the interfacial points, furnishing: x k 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 3-D flows, D was expressed by equations ( 11) and (12).
Dx, Dy and Dz 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].

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.
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 r 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 Adams-Bashforth scheme was used for the advection-diffusion equation, while a fourth order Adams-Bashforth scheme was used for the filtered Navier-Stokes 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 Figure 4.The "physical dominium" and the numerical dominium established for the calculations.
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.

Large flow structures and sediment transport
Figure 6 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).The resuspension fluxes occur mainly right after the crests, where large vortical structures are formed.Figure 7 presents results of calculations made by [2] in which vorticity and concentration fields were compared, furnishing results comparable to the present study.The calculated two-dimensional flow (Figure 7a) 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 2-D experiment preceded the 3-D experiment presented here.

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.Figure 8 shows the formation of ripples close to the channel inlet at the beginning of the experiment (Figures 8a and 8b), the predominant erosion of the crests (Figures 8b to e), the predominant deposition on the valleys (Figures 8b 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).

CONCLUSIONS
The present study shows that 3D moving beds can be adequately calculated using immersed boundaries together with the advection-diffusion 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.

Figure 3
Figure 3 shows the balance of mass for the cell just above the bed, in a 3-D situation, used to calculate resuspension and deposition.a represents advectives fluxes a u c i ( ) = ⋅ , d represents diffusive fluxes

Figure 3 .'
Figure 3. Sediment fluxes close to the bed which result in deposition or erosion.
Figures 8 a 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 8 .
Figure 8. 3D evolution of the moving bed calculated with immersed boundaries, large eddy simulation and the advection-diffusion equation.
followed through the present procedures, showing that this methodology is adequate to evaluate bed modifications and sediment transport in alluvial flows.