Numerical simulation of aeolian saltation within the sediment transport layer using granular kinetic theory

Juan Pedro Marval Lugo¹, Luis Rojas-Solorzano² And Jennifer Curtis³

1 Universidad Nacional Experimental Francisco de Miranda. Department of Energetic. jpmarval@gmail.com

2 Universidad Simon Bolivar. Department of Conversion of Energy. rrojas@usb.ve 

3 Universidad de Florida. Department of Chemical Engineering. jcurtis@che.ufl.edu 

ABSTRACT

Most granular flows at environmental conditions are unsteady and exhibit a complex physical behavior. Dune formation and migration in the desert are controlled not only by the flow of saltating particles over the sand bed, but also by turbulent atmospheric airflow. In fact, sediments are transported by the atmospheric airflow within a thin layer only a few centimeters above the sandy surface. These jumping particles reach a maximum sediment mass flux level at a certain delay time (known as the "saturation time") after the initial movement by sliding and rolling begins. Unlike sediment transport in water where the particles are lifted by the turbulent suspension, the saltating particles are kept alive in the layer mainly due to particle-particle and particle-bed collisions. In order to model this Aeolian transport of sand, Jenkins and Pasini (2005) proposed a two-fluid model (one-dimensional and steady state) using Granular Kinetic Theory (GKT) to describe the solid-phase stress. The present work extends the original idea of Jenkins and Pasini by using a more robust model of GKT for the kinetic/collisional contributions to the solid-phase stress tensor, together with a friction model for sustained contacts between particles. In addition, a standard k-ε turbulence model for the air and a drag model for the interaction between the phases is employed. A rectangular 2D geometry was chosen with a typical logarithmic profile for the inlet air velocity, along with an initial amount of sand at rest in the lower part of the simulation domain, resembling the particle saltating flow commonly seen in the vertical middle plane within Saltation wind tunnels. This model is validated with experimental data from Liu and Dong (2004) and the results given by Pasini and Jenkins (2005) and Marval et al. (2007). A good estimation for the particle erosion and mass flux in the Saltation layer is predicted, even though the profiles of mass flux and concentration within the transport layer are very thin and close to the bed.

Keywords: Saltation layer, Granular kinetic theory, CFD

Simulacion numérica de la saltación aérea dentro de la capa de transporte de sedimentos usando teoria cinetoca granular

RESUMEN

La mayoría de los flujos granulares que ocurren en condiciones ambientales son inestables y exhiben un complejo comportamiento físico. La formación y migración de dunas en los desiertos son controladas no sólo por el flujo de partículas en saltación sobre el lecho arenoso, sino también por el flujo turbulento de aire atmosférico. Estas partículas en saltación alcanzan un flujo másico de sedimentación máximo luego de cierto tiempo luego de la movilización inicial por desplazamiento y rodamiento. A diferencia del transporte de sedimentos en agua, en el cual las partículas son elevadas por la suspensión turbulenta, las partículas en saltación son avivadas en la capa principalmente por los choques existentes de partícula-partícula y partícula-lecho arenoso. Para el modelado de la teoría de Transporte Eólico de Arena, Jenkins y Pasini (2005) propusieron un modelo bifásico (un estado estable y unidimensional) usando la Teoría Cinética Granular (TCG) para describir los esfuerzos en la fase sólida. El presente trabajo extiende la idea original de Jenkins y Pasini usando un modelo más robusto de TCG para las contribuciones cinéticas/colisionales a los esfuerzos tensores de la fase sólida, junto a un modelo de fricción para los contactos continuos entre partículas. Fue seleccionada una geometría 2D rectangular con un perfil logarítmico típico para la velocidad de entrada del aire. Este modelo es validado con datos experimentales de Liu y Dong (2004) y los resultados dados por Pansini y Jenkins (2005) y Marval (2007). Una buena estimación para la erosión de las partículas y el flujo másico de saltación es predicha.

Palabras Clave: Saltation layer, Granular kinetic theory, CFD

Recibido: marzo 2008 Recibido en forma final revisado: octubre 2011

INTRODUCTION

Computational Fluid Dynamics (CFD) has recently impacted many fields of research (e.g. atmospheric flow, microfluidics) where before its application was unthinkable mainly due to the problems associated to the wide range in length and time scales, as well as the great computational effort necessary to obtain reliable results. Saltating flows are a type of multiphase flow where CFD can now have a great impact. The flow of sediment in the saltating layer is a mixture of two different phases. Air is the continuum phase or primary phase and sand is the dispersed phase. When the sand particles are considered as an immiscible "fluid", it is possible to use theories developed for two-phase flow. In this case, the dispersed phase of particles and the fluid phase are interpenetrating continua. Recently, the saltation layer has been modeled with mixture models (Ji et al. 2004; Alhajraf, 2004) neglecting the collisional effects intrinsic to the motion of the particles and the coupling of the particle interactions with the air stream. However, a complete description of sediment transport should consider three different flow regions present in the phenomenon and their interactions. The regions are as follows: (1) a bottom region where the particles are at rest (sand bed); (2) a region above the sand bed where the particles are jumping (saltating layer) and (3) a top region where there is only air. Therefore, the mixture model is inadequate, for modeling this phenomenon because it does not permit a description of the prevailing stratification. Multi-Eulerian models enable the simulation of the three co-existing regions and consider the movement of the cloud of particles with equations that are separate from the gas equations. This approach involves to incorporate a description for the solid-phase stress tensor using the Granular Kinetic Theory (GKT).

The movement of sand dunes is often a problem in desert areas. These dunes move onto communication routes like highways, freeways or other places of interest (e.g. cities, lakes, and parks). The dune fields in southern Morocco (Sauermann, 2005) and Dunes of Coro in Venezuela are examples of these problems. In these cases, the sand must be shoveled daily from the upwind side to the downwind side of the road using large bulldozers, obstructing the flow of traffic. The situation associated with the Venezuela dunes is more complex because the highway runs onto a thin portion of land (called an isthmus) which makes it impossible to construct another freeway without passing through the dunes. Thus, only a better understanding of Aeolian transport on a micro- and meso-scale, and generating improved mathematic models, will aid in predicting with reasonable precision the flow behavior on macro-scale. These models enable the design and development of control devices for dunes motion or in the construction of highways, overpasses, tunnels.

The present work introduces a numerical model for predicting the flow of sand particles in a saltating state over a sand bed using an Eulerian model for both phases. The gas phase turbulence closure model is the standard k-ε model, and the solid-phase stress is described using GKT, with a frictional model which is applied for sustained particle contacts. The mathematical model is implemented in the commercial code ANSYS-Fluent® and computed results are compared with experimental data given by Liu & Dong (2004) and the model predictions given by Pasini &Jenkins (2005) and Marval et al. (2007).

OVERVIEW OF SALTATION

When sand particles are dragged by a low velocity wind, two different mechanisms of transport are related to the suspension. These two mechanisms are known as "Saltation" and "Reptation" (Andreotti et al. 2006), recognized by Bagnold (1941) more than 60 years ago in his classical book "The physics of blown sand and desert dunes". Saltation and Reptation together contribute to the particle mass flux over the sandy surface and control the dynamics of a dune given its particular shape characteristics in the desert. Based on previous research, it is known that the particles interact with each other and with the sand bed (Andreotti et al. 2002). From the early work of Bagnold’s until now, much research effort, both experimental and theoretical/numerical modeling has been devoted to trying to describe the complex phenomenology of saltation.

Many experimental works has given insight into details about the Saltation layer (Liu & Dong, 2004; Nalpanis et al. 1993; Spies et al. 1995; Butterfield, 1999; Stout, 1998; Iversen & Rasmussen, 1999; Rasmussen & Sorensen, 1999, 2005; Okoli, 2003; Nishimura & Hunt, 2000; Ni et al. 2002; Dong et al. 2003, 2004; Bauer et al. 2004). Research results are presented in terms of two length scales. The small length scale is on the order of the particle diameter and consists of measurements of the ejection and impact velocities, ejection and impact angles, trajectories of the different types of particles (snow, mustard seed, quartz) and particle size. The mean ejection velocity is proportional to the friction velocity u* (which characterizes the flow and has dimensions of speed, but is a measure of shear, Eqn. 1) and does very much with the grain size.

On average, the sand particles are lifted in the airflow at an angle between 34˚ and 41˚ and descend with an impact angle between 11˚ and 14˚. In general, these angles decrease slightly when the particles are larger and increase significantly when the friction velocity is slightly raised. On the large scale, the following measurements are reported: particle mass flux, particle velocity and concentration profiles, saltation layer height, length of saturation and others macroscopic flow features. The majority of researchers agree that the particle mass flux and concentration of particles in saturation decrease exponentially with the height above the sand bed, unlike the particle velocity profile which increases exponentially with the height (White, 1996). Few studies have focused on the saturation length; however the general argument is that the equilibrium state is reached in a short time for a determinate length. Only recently, Dong et al. (2004) concluded that a completely saturated state within the Saltation layer is achieved after a few meters, depending on the friction velocity.

Figure 1 shows a summary of the main characteristics of the Saltation layer based on experimental observations. All the sketched profiles are plotted on a logarithmic scale in the vertical direction. The wind blows from left to right with the well-known logarithmic profile (Schlichting, 1979) (and erodes the sand bed starting at point A. Between points A and B, the air velocity profile is modified as momentum is transferred from the air to jumping particles until equilibrium is reached after point B. Basically, the saturation length depends on the friction velocity and particle diameter. After the saturation length, the saltating particle flow reaches a fully developed state and all flow variables only depend on the y-direction. Recent experimental work (Dong et al. 2003) has confirmed that inside the Saltation layer the air profile is modified by the presence of the particles and that the air profile has a convex shape near to the sand bed (y<δmom). Towards the top of the Saltation layer (y>δmom), (δmom is the layer where the appreciable momentum interchange occurs between both phases) the air profile is weakly influenced by the grains. In fact, the particle concentration and mass flux near to the top of the Saltation layer tend to zero.

Figure 1. Saltation layer skew.

The main purpose of the researchers has been to obtain reliable expressions that allow to calculate the total mass flux under specific conditions. Bagnold (1941) was the first to report an expression for the saturated mass flux over the sand bed, which shows dependency on the third power of the friction velocity. This equation does not take into account that the mass flux on the Saltation layer is only present when the friction velocity is greater than a threshold velocity for the initiation of particle movement or when the friction velocity is greater than the impact velocity for grains already in a saltating state. Experimental research has improved our understanding of saltating flow, but these studies have mainly enabled the development of mathematical expressions (analytical and experimental correlations) for predicting the mass flux at saturation conditions. Sorensen (2004) has reported an interesting summary of equations for predicting the saturation mass flux, this compilation shows mathematical expressions with a dependency according to the following relationship:

 Where V is the ratio between the friction velocity and threshold velocity, and f(V) is determined experimentally. For Bagnold’s equation this functions is a constant, which gives an untrue mass flux when the friction velocity is smaller than the pickup velocity which is the velocity required to start the movement of the particles. Each of those expressions for calculating the mass flux shares in common that the saltating particles travel in the wind direction totally free without obstacles or walls that blocks it along of its trajectory. This feature was described previously to indicate that these kinds of macro-expressions are invalid for predicting a flux of saltating particles interacting with a kind of static solid surface. In any case, the total mass flux over a vertical area can be calculated by mean of the following expression:

If we want to apply the Eqn. [3], it is necessary to know a priori the sand velocity and sand concentration profiles previously modified by the presence of obstacles and the form of sand bed eroded, and therefore, only a multiphase numerical simulation could be able to give a result close to the reality.

AEOLIAN SALTATION CLOUD MODELS

The Aeolian process begins when the aerodynamic force (i.e. drag) dislodges a few grains, which roll and slide over the sand bed until they enter into Saltation. The entrained grains are accelerated by the wind along their trajectory mainly by the drag force before they once again impact the bed. The particles which impact the bed probably rebound or dislodge other particles, creating a grain chain reaction. The visualization of particle collisions in experiments (Nalpanis et al. 1993; McEwan & Willets, 1993) was derived statistically to get the "splash process" (the interaction between an impacting grain and the bed) and this is modeled by the splash function (Ungar & Haff, 1985). The traveling grains extract momentum from the air and, consequently, the air decelerates and the process of direct aerodynamic entrainment decreases. This negative feedback mechanism reaches equilibrium in the saturation time, (the time at which the air will hold the maximum amount of sand particles in saltation) and, then, the grains are kept alive in the Saltation state only by their collisions with the sand bed. The rate of particle transport becomes constant when equilibrium is reached. The statistical models are based on the four sub-processes previously described: aerodynamics entrainment, the grain trajectory, the grain/collision, and the wind velocity modification. A common feature of all the models is that they consider saltating flow in a free condition of movement without any restrictions or obstacles. This premise makes the statistical model useless for the prediction of the patterns of erosion and/or deposition around solid objects caused by the grain-wall interaction.

CFD MODELS AND GRANULAR KINETIC THEORY

Currently, a CFD methodology for the appropriate treatment of the sediment transport by Aeolian Saltation does not exist. Recently, Ji, et al. (2004) and Alharaf (2004) proposed interesting models to simulate the large-scale Aeolian sediment transport, without considering in detail the fluid/ particle interaction and without taking into account the collisional interaction between particles as a fundamental parameter in their considerations. In both works, they used a homogenous model with a single set of equations, and the disperse phase was modeled with a slip velocity based on a drag function. Although, the effect of particles on the wind profile is weak, the air velocity profile within the Saltation layer is modified by the momentum transfer between the two phases. Therefore, in order to consider the mixture model, it would involve knowing a priori the concentration profile, velocity profile and other characteristics of the mixture, that are commonly unavailable. Therefore, although these previous works present significant advantages in the use of CFD to simulate the Saltation layer, they do not take into account some fundamental aspects present in the flow. Ji, et al. (2004) utilized the logarithmic and total entrained particle flux profiles observed in experiments in order to include functional relationships into his model. Meanwhile, Alhajraf (2004) focused his work on an algorithm for treating a moving boundary (erosion-deposition) based on handling an artificial computational domain.

On the other hand, there is some uncertainty about the particle suspension mechanisms present in sheet flows (water) and saltating flows (air), while experimental results only offer some qualitative characteristics of the two types of flows. Jenkins & Hanes (1998) implemented, for the first time, a simple one-dimensional and steady two-fluid model using Granular Kinetic Theory for describing the interaction among particles in the sediment transport layer both in water and air. Their results showed a qualitative good estimation of the solids erosion, which suggests that the Granular Kinetic Theory could be used to simulate this phenomenon with success. Pasini & Jenkins (2005) extended the previous work by adding an extra term in the gas momentum equations, originating from a second averaging process (given in Hsu et al. 2004) which described an additional mechanism of suspension due to turbulence effects together with the granular pressure gradient. Furthermore, Pasini & Jenkins (2005) considered a flow regime between the Saltation regime and turbulent suspension regime called the collisional regime while, at same time, they recognized that this regime had not yet been observed experimentally. Their results gave an overestimation of the sand mass flux in comparison with experimental results from particles in the Saltation state.

However, there is enough evidence from quantitative experimental results and analyses based on the ratio of the particle terminal velocity and the vertical turbulence fluctuation velocity (Sauermann, 2001; Pye & Tsoar, 1990) demonstrating that in the Saltation layer the turbulent suspension is negligible. Instead, the particles are only kept in movement in the Saltation layer by bed-particle collisions and momentum transfer from the air stream to the particles. In any case, the Saltation layer represents a dilute regime where the maximum volume fraction is around 5.10-4 (Liu & Dong, 2004), indicating that the saltating particles contribute only to the kinetic contribution to the solids viscosity and not the collisional regime like Jenkins & Hanes (1998) and Pasini & Jenkins (2005) suggest.

MATHEMATICAL MODEL

According to Dong, 2004, three well-known regions are considered in our simulations, according to the solid volume fraction and velocity fields. One region corresponds to the sand bed, where the solids packing level is high, and the gas and solid velocities are equal to zero. This zone cannot move like a fluid; it is a completely solid region, which is superficially eroded by the wind action. The second region is few centimeters above the sand bed (Saltation layer), where particles travel quickly, driven by a turbulent shearing airflow. In the third region, outside of the Saltation layer, there is only air moving parallel to the bed. In this work, we proposed an Eulerian model (inhomogeneous model) for both phases based on the fact that there are three zones with different concentration and velocities (similar to a stratified flow). Furthermore, the Stokes number for the Saltation layer condition is estimated to be greater than 1 according to the user’s manual of FLUENT®-ANSYS

Researchers agree on the use of the governing equations originally proposed by Anderson & Jackson (1967) and later on, by Jackson (1997), for the description of the gas-solid flow using an Eulerian approach. However, van Wachen et al. (2000) affirm that the difference with the Ishii (1975) equation sets is negligible on an engineering scale. We have adopted the Ishii equations applied to gas-solid flow by Enwald et al. (1996). These equations assume a suspension of identical spherical particles in an incompressible flow, where particles are characterized by a radius "a", and "L" is the shortest significant macroscopic length-scale related to the particle motion. For L>>a, both scales separate and define a process by averaging length-scale, l, given by the condition: a << l << L. This means that the variables are averaged over a region that is larger than the particle size, but smaller than the characteristic system length. Table 1 show the equations used in this simulation. The governing equations are expressed in Eqns. (4-10).

The closure model for the solid phase is based on the Granular Kinetic Theory, which comes from an analogy with the gas kinetic theory developed by Chapman & Cowling (1970). The main idea of this model consists of introducing a granular temperature, which is a measure of the energy level of the particle velocity fluctuations (Ogawa, 1978). Through the solution of the granular temperature equation, Eqn. (12), it is possible to obtain the pressure for the solid phase, Eqn. (13), and the entire transport coefficient, Eqns. (15-21). Lun et al. (1984) derived all these expressions considering the inelastic nature of particle collisions, the particles as identical spheres of diameter ds, composed of a material of density ρs, the interaction between particles occurring by instantaneous binary collisions, and the granular solid stress tensor resulting from kinetic and collisional contributions.

A constant restitution coefficient is used to take into account the energy dissipated by inelastic particle collisions. Equation (22) shows the radial distribution function, which can be interpreted as a correction factor that modifies the probability of collision when the solid volume fraction is high (dense regime). Regarding the solid viscosity, the model has been improved since the first proposal done by Lun et al. (1984). Gidaspow (1994) and Syamlal et al. (1993) implemented some changes in the original equation of Lun et al. (1984). However, Gidaspow (1994) did not account for the inelasticity of the particles in the kinetic contribution, and secondly, Syamlal et al. (1993) neglected the kinetic contribution in the dilute regime. Both equations are very similar, except when the solid volume concentration is low (dilute regime), where Gidaspow’s equation gives an inadequate finite value for the solids viscosity, and Syamlal et al.’s equations yield the correct trend in the solids velocity, but underestimate it, as the solids fraction approaches zero. Hrenya& Sinclair (1997) solve this problem by starting from the equation of Lun et al. (1984) and introducing a ratio between the mean free path and the characteristic length of the system. Hence, Eq. (17) conserves two important features: the solid viscosity is zero when the solid volume fraction is zero, and the solids viscosity is adequately predicted in the dilute regime. The Figures 2, 3 show comparisons among the conductivity and solids viscosity in low concentration for the Syamlal model and the Hrenya & Sinclair model. The most important features is that Syamlal´s curves behavior give a tendency to zero when the concentration is zero, but the Hrenya and Sinclair curves tend to Syamlal et al. curve when the mean free path is very low. In this case, the mean free path has an important role in the simulation due to very low concentration in the Saltation layer and the effects of small changes in the conductivity and solid viscosity which could affect considerably the results of the simulation.

Figure 2. Comparison of solid granular viscosities from different kinetic theory models αmax=0.65, e=0.9

The rapid granular flow regime is properly described by the GKT (Van wachen et al. 2000). However, when the solid volume fraction is very high (quasi-static flow regime) the GKT underestimates the solid viscosity. In fact, the granular temperature is very low as the result of the high solids packing and weak velocity fluctuations; therefore, this situation can be considered like a fluid-solid phase change. This new regime is characterized by long-lasting contacts between particles. Johnson and Jackson (1987) and Johnson et al. (1990) modeled a granular flow based on the work of Lun et al. (1984), but added the Coulomb  friction stress to the solid-phase stress in the momentum equation for quantifying this phenomenon. However, we use the frictional stress Eqn. (23) described by Syamlal et al. (1993), originally proposed by Shaeffer (1987), instead of the Coulomb equation. Johnson et al. (1990) proposed the functional form of the semi-empirical equation for the frictional pressure Eqn. (24) as a function of the solid volume fraction and, maximum and minimum solids packing. The value of the constants (Fr, n and p) used in this equation originates from Ocone et al. (1993) due to the similarity existing between the particle properties (diameter and density) used in their work and the properties of desert sand. A classical k-ε model for turbulence is used as a closure relation for the gas phase, Eqns. (28-32).

Figure 3. Comparison of solid granular conductivity from different kinetic theory models amax=0.65, e=0.9

The simulation was carried out in ANSYS-FLUENT®, which uses a control-volume-based technique to convert a general scalar transport equation to an algebraic equation that can be solved numerically. This control volume technique consists of integrating the transport equation about each control volume. Fluent uses by default second-order accuracy for the viscous terms, meanwhile for the convection terms a second order scheme was chosen. The SIMPLE (Patankar, 1980) method is used as a coupling algorithm between pressure and velocity, and for the temporal discretization an implicit second-order accuracy scheme was chosen.

GEOMETRY, GRID, BOUNDARY CONDITIONS AND INITIAL CONDITIONS

The geometry chosen for the simulation matches the experimental setup and results obtained by Liu & Dong (2004), which were made in a saltation wind tunnel with a total length of 37.78m and a working cross-section of 0.8m x 1 m. A two-dimensional channel with dimensions Lx x Ly was used for representing the computational domain. The vertical dimension Ly was chosen following two criteria: (a) it should be greater than the experimentally reported boundary layer thickness inside the tunnel; and (b) higher than the height for which the sand mass flux over the bed is almost zero at the maximum test velocity. Thus, based on the previous considerations, Ly was taken equal to 0.6m, while Lx has a value of 1m (Figure 4).

Gambit® v2.2.30 was used as the geometry and grid generator, with a uniform mesh in the x-direction and non-uniform mesh in the y-direction. A finer mesh was placed in the vertical direction near to the sand bed, where larger gradients are expected. Meanwhile, a coarser mesh was chosen for the top wall, where a small gradient is expected. The mesh refinement analysis was performed by comparison of the horizontal sand velocity at height of 0.062m from the bottom wall. Meshes of 15360, 24000 and 34560 cells were compared. A difference of 6.38% in the horizontal sand velocity were encountered between the 15460 and 24000 cells meanwhile a difference of 3.44% were encountered between the 24000 and 34560 cell. According to those results, the intermediate mesh with 24000 elements was adopted for the rest of the study. Details of the selected mesh are displayed on Figure 5. After the mesh sensitivity analysis was carried out, the computed results were compared with experimental data (Dong, 2003) and other results of simulations (Marval et al. 2007). The mesh contains 24,000 cells, depicting an aspect ratio with respect to the particle diameter in an approximate range between 5-60.

Although the chosen benchmark experiments have several free-stream velocities, we have taken into consideration the highest reported value of U∞=18 ms-1 due to good behavior found in the experimental curves of sand concentration and sand mass flux. Equation (33) shows the logarithmical air profile used in the simulation inlet boundary velocity, which was fitted based on air-wind profile data, without sand transport, presented by Liu & Dong (2004). The logarithmical air profile begins just over the sand bed and finishes in the limit of the boundary layer. Below the boundary layer the air velocity is zero, while above the boundary layer it is equal to U∞. The origin of the coordinate system is located on the bottom left corner, For the outlet boundary condition, a gauge pressure equal to zero was specified. Even though the top wall in the Saltation tunnel is located 1.0 m from the bottom, the top boundary condition for the numerical simulation is set at 0.6m (Ly) in order to reduce the number of cells and the total simulation time. Hence, a free slip (irrotational flow) condition for both phases is set at the top wall (instead of a no-slip boundary condition). A non-slip condition on the bottom wall for both phases was used. A criterion based on the hydraulic diameter and turbulence intensity was applied for the turbulent boundary conditions at the inlet and backflow outlet, with values of 0.6 m and 0.05%, respectively.

Transient simulations are strongly dependent on the initial condition, especially in this case where a quasi-steady state in the Saltation layer must be achieved in a reasonable interval of time (t1, see Figure 6) in order to match experimental data. The experimental procedure used in the wind Saltation tunnel indicates that the air flow over the sand bed is fully developed before the initiation of the drag on the particles by the air. This condition was obtained in experiments by covering the sand bed with a geo-textile fabric and then quickly uncovering the bed after the airflow was already fully developed. An equivalent procedure was implemented in ANSYS-Fluent® in order to mimic these conditions, i.e., the sand bed remained at rest, while the air flow was fully developed. This condition was obtained by temporarily not solving the equations of granular temperature, volume fraction, and, furthermore, by deactivation of the drag and lift forces. Fully developed air flow was obtained in a real time of only 0.08 seconds with a time step of 1.10-3 seconds. An initial sand bed with a thickness of 10mm and a concentration of 0.62 was placed on the bottom part of the geometry, which is equivalent to a weight of 16.43 kg/m of sand. The Table 1 shows the parameters used in the simulation.

Figure 4. Geometry and boundary conditions for the simulation.

Figure 5. Mesh samples used in the simulation: a) Top left corner, b) Bottom left corner, c) Whole view

Figure 6. Schematic of the temporal evolution of dimensionless mass flux and dimensionless total mass in the simulation

Table 1. Parameters used in the simulation

RESULTS AND DISCUSSION

Unquestionably, the Saltation layer represents a complex problem where the time dependence, small-scale variables (e.g. concentration profiles, mass flux profiles, etc.), and large-scale variables (mass total and mass flux evolution time) are closely inter-related. The strategy followed in this study consisted in obtaining a temporal condition along with a large-scale condition that allowed to make the comparison between the small-scale variables with experimental results. Due to the simulation configuration, there is always a sand mass loss in the domain from the beginning of the simulation until the moment the sand is depleted; therefore, there exists only a short period of time (t2-t1) where a match with experimental data can be made. This condition is called a quasi-steady state, which is shown schematically in Fig. 6, based on the dimensionless total mass and the dimensionless mass flux. From t=0s to t=t1 the erosion begins over the sand bed; after t=t1 is reached, the saturation in the Saltation layer and the mass flux is kept constant until t=t2. The last part of the process occurs between t=t2 and t=t3; in this interval of time, the mass flux is reduced to zero. The previously described process also occurs in the experimental tests and the characteristic times are mainly a function of the strength of the air-stream. The simulation total time was around of 48 CPU-hours for 4 seconds of real time.

Temporal evolution of total mass and mass flux

Figure 7 shows the sensitivity of the sand bed to the erosion according to the chosen frictional pressure model. The GKT model uses the kinetic/collisional pressure Eqn. (13) in the calculation of the frictional viscosity Eq. (23). This model led to nearly 50% sand bed erosion in 2 seconds, unlike results based on the frictional pressure model of Ocone et al. (1993) for which the sand bed is hardened due to the high solid pressure added in the solid momentum and granular temperature equations; therefore, the loss of sand mass is smaller in the latter case. We can also see that the variation in the FPL for the model of Ocone et al. (1993) produces only changes during the initial period of erosion, but the slope for both curves is the same. The rest of the simulations presented as the results in this work use the Ocone’s model like frictional pressure model.

The Figure 8 shows a similar behavior of the dimensionless total mass for three values of the mean free path. Over the 4 seconds the erosion is around 30%, which compared with the erosion obtained by Liu & Dong (2004) is high for the same interval of time with Saltation layer in total saturation (18.12% of erosion in 4 seconds). The difference in total erosion between the curves with a mean free path of λ=0.001 and λ=0.1 at 4 second, is 5%, depicting the best behavior among all cases. Previous results shown in figure 8 indicate that it is possible to adjust the Hrenya & Sinclair (1997) model through an appropriate selection of the mean free path to control the erosion. Also, we can see that around 2 seconds, all the curves in the Figure 8 change the slope, with a strong acceleration of the erosion.

Figure 7. Comparison of temporal evolution of dimensionless total mass from different frictional pressure models. Solid Viscosity and conductivity by Syamlal et al. (1993)

Figure 8. Comparison of temporal evolution of dimensionless total mass from different mean free path of the granular viscosity and conductivity by Hrenya and Sinclair´s model. All curves use a frictional pressure by Ocone et al. (1993) αmin,s=0.5

The dimensionless mass flux at the outlet in figure 9 shows sensitivity to the mean free path using the Hrenya and Sinclair model. Here, we see that around 2 seconds the dimensionless flux mass suffers a sudden increment until it reaches a value of 2.5, which is closely related with the value observed in figure 8, where also there is a change in the slope of the curves. In base to previous results, we can estimate the time t1 to be around 2 seconds, while a time of 4 seconds represents an intermediate time between t1 and t2, which can be used for comparing the small-scale variables.

The curve with a mean free path of λ=0.01 shows a more stable behavior (without oscillations), but at the same time, it presents the greater rate of erosion among all curves. However, although there exists uncertainty about the presence of oscillations for a mean free path equal to 0.001, this behavior can be associated with the changes in the conductivity and solid viscosity in the Saltation layer given by λ. The saturation mass flux used to do the dimensionless curves in the figure 9 was of 0.744 kg s-1 per wide, which corresponds with a total mass flux in saturation, as reported by Liu & Dong (2004) for a free stream air of 18m/s.

Solid volume fraction profiles

The profile of sand concentration within the Saltation layer presented in Figure 10 shows a comparison between the simulation results and experimental data. Liu & Dong (2004) data show an exponential decrease of the solids concentration with height over the sand bed, unlike the numerical results, the general behavior until α = 0.0002 is exponentially decreasing, the three curves descend linearly with height, producing a much higher concentration compared with the experimental data. The increase of solid viscosity with the reduction of the mean free path in the Hrenya & Sinclair (1997) model improves the predictions of the concentration profile causing more particles to scatter above the sand bed. In figure 11 a sand bed thickness of 0.007 m approximately, is obtained after 4 seconds, using a criterion of 0.5 for the solid volume fraction because this is the limit value for the activation of the frictional. model according to Johnson et al. in equation 24.

Figure 9. Temporal evolution of dimensionless mass flux at the outlet from different mean free path of the granular viscosity and conductivity by Hrenya and Sinclair´s model. All curves use a frictional pressure model based on Ocone et al. (1993) αmin,s=0.5

Figure 10. Comparison of the solid volume fraction above the sand bed with experimental data

Figure 11. Solids volume fraction within the sand bed and transport layer

Dong et al. (2004) gave evidence of the accumulative effect of the mass flux and, consequently, the distribution of sand particles in the Saltation state within a wind tunnel. They demonstrated that the mass flux in the Saltation layer increases with distance until it reaches a total saturation, after which, the concentration of particles and mass flux remains constant. This behavior can be seen qualitatively in the figure 12, where the concentration of sand increases over the sand bed along the x-direction. The thickness for constant saturation is reached near to the 0.5m over the horizontal length (Lsaturation), and the Saltation layer thickness (δsal) is around 0.08m, considering a criterion for the solid fraction volume of 1.10-6 which is the minimal fraction when the concentration becomes constant.

Velocities

Marval et al. (2007) presents a sand velocity profile reaching a thickness scarcely of 0.05m with a mean free path for the Hrenya and Sinclair model of λ=5. Instead, in our new results, the thickness of particles with velocity different to zero over the sand bed reached just more than 0.12m with the same model of Hrenya and Sinclair (1997) and a mean free path of 0.01 (Figure 13). This result puts in evidence that the reduction of the mean free path contributes with the scattering of particles from the sand bed towards the sediment transport layer. Of course, was recognized that the particle should be more spread with a velocity nearer to the air velocity over the height needed for this air free stream velocity (it should be around 0.4m). However, the influence of the conductivity and solid viscosity can help remarkably to improve this behavior if they are chosen with adequate values or a combination of them. In Figure 14, it is possible to see the sand bed in complete rest (bottom part), and over it there is a cloud of particles dragged by the action of the wind from left towards the right. These clouds of particles increment its thickness along the horizontal direction, like it is also seen in Figure 12 for the concentration contours.

It is very important to mention the temporal evolution of this cloud of particles before they reach this quasi-steady condition. In the first milliseconds, it begins a small, slight and thin layer of particles to move over the sand bed, while in a parallel form, from upwind over the sand bed, it begins a process of particles heating (increment of their granular temperature) due to the drag of wind, constituting a kind of little cloud of particle lifting. From this time, the initial cloud of particles grows in the vertical- and much more in the horizontal direction, propagating downwind until they reach the outlet of the domain. This process is achieved in 2 seconds, justly when is observed that the slope of the

Figure 12. Contour of solids volume fraction at 4 seconds, Hrenya and Sinclair´s model λ=0.001

Figure 13. Velocity profiles at the outlet at 4 second, Hrenya and Sinclair´s model λ=0.01

Figure 14. Contour of sand x-velocity at 4 seconds, Hrenya and Sinclair´s model λ=0.001

dimensionless total mass curve changes suddenly (more mass loss per second, Figure 8) and when the total mass flux reaches 2 seconds (see Figure 9). Obviously, when this cloud travels across the outlet boundary, the large scale variables are affected, and then these variables reach the quasi-steady state of saturation in the Saltation layer. In Figure 10 the oscillations in all the curves, with the λ=0.01 exception, shows a better behavior according to the description presented permanently for de addition in the flux in the outlet

Mass flux profiles

Figure 15 shows the mass flux profiles at the outlet. The mass flux occurs in a thin layer close to the sand bed, over a thickness of approximately 0.08m for λ=0.01, meanwhile the results of Marval et al. (2007) only give a thickness of 0.013m. On the scale of Figure 13, the dimensionless mass flux of Liu & Dong (2004) appears like a vertical line with approximately zero dimensionless mass flux, once again indicating that the sand particles in the simulation need to spread away from the sand bed in order to reduce the mass flux near the bed. The influence of the mean free path is put again in evidence when the mass flux profile is extended over more height from the sand bed. Here, it is important to indicate that the Liu & Dong data (2004) was obtained with mechanical traps over the sand bed, which only can capture particles jumping from 1cm above the bed, and thus, there exists uncertainty about the quantity above the sand that can have passed below this height and that was not registered by the sand traps. In any case, the simulation shows a result compatible with a logical reasoning, due to that the sand velocity is almost zero in the proximity to bed sand, while its concentration is very high, and when we evaluate the mass flux in agreement with the argument of equations [2] and we obtain values for the mass flux profile relatively high immediately, over the sand bed, but it decreases suddenly with the height. However the authors recognize that some adjustments has to be realized in the computational model which allows the reproduction with more details of the real behavior of this complex flux according to the results presented in the experimentation.

The curve of Pasini & Jenkins (2005), showed in the Figure 16, has a similar behavior with our numerical results. However, Pasini & Jenkins’s predictions show a slight peak of mass flux near to the bed, although, it has almost the same distribution over the height. From this comparison we can assure that the turbulent suspension added by Pasini & Jenkins (2005) in the y-momentum equation did not have an influence in the simulation, and that the unique mechanism for suspension is due to the solid pressure gradient, which is closely related to the variation of the granular temperature. The Equation 34 shows the Shield Parameter which is represented in the figure 16.

Figure 15. Comparison of dimensionless mass flux in the sediment transport layer with experimental data over the outlet at 4 seconds

Figure 16. Comparison of dimensionless mass flux above sand bed with numerical results by Pasini and Jenkins (2005). Height of sand bed: 0.007m

Solid viscosity, conductivity and granular temperature

The figures 18 and 19 (solid viscosity and conductivity profiles) show similar results in agreement with the behavior of the granular temperature profile presented in the figure 17. The granular temperature is a measure of fluctuation of sand velocity, and this variable into the sand bed must be zero because the sand bed is totally in rest.

Over the sand bed, the particles begin to acquire granular temperature, having the larger value near to 0.12m (same value of thickness found for the sand horizontal velocity). Above this height, the granular temperature is zero due to the absence of particles.

Based in the fact that the conductivity and solid viscosity have, for definition, a proportional dependency with the granular temperature, we can explain some features observed in the figures 18 and 19. Firstly, the solid viscosity in figure 18 shows a high value into the sand bed (y<0.01m) because the frictional viscosity is present right there. Over the sand bed, the frictional viscosity does not contribute with the solid viscosity, and thus, only the kinetic/collisional contribution is present. Alike the granular temperature, the solid viscosity is very low immediately above the sand bed, and then it is easier for the wind to drag the particles. Following the granular temperature profile, the solid viscosity increases accordingly with the height until the granular temperature ceases. The last description comes to explain the over mass flux near to the sand bed due to the low viscosity present in this zone.

With respect to the conductivity and its effects over the simulation above the sand bed (y>0.01), it is completely clear that the viscosity helps the scattering of particles into the Saltation layer, due to the greater conduction of granular temperature and subsequently, it provides a better pressure gradient in the layer that lifts more particles resting onto the sand bed towards the cloud.

The behavior of the solid phase has a duality like-gas and like-liquid (or solid) behavior. Into the sand bed, the solid phase behaves like a liquid, because the granular temperature is low, but the viscosity is high, meanwhile, in the Saltation layer, where the concentration is low, the solid phase behaves like a gas (following the Kinetic Theory). In the transport layer, if the granular temperature is low, the viscosity is high and conversely.

Figure 17. Comparison of granular temperature profile over the outlet at 4 seconds

Figure 18. Comparison of solid viscosities profile over the outlet at 4 seconds

Figure 19. Comparison of solid conductivity profile over the outlet at 4 seconds

CONCLUSIONS

From this numerical study, several conclusions may be drawn:

1. The work of Pasini & Jenkins (2005) on the Aeolian transport using GKT is extended, considering slightly inelastic particle-particle collisions and incorporating an improved two-dimensional transient model with a frictional sub-model to describe the sustained contacts between particles.

2. The simulation well described the solid-like characteristic of the sand bed, a sand transport layer over the sand bed and an air free flow outside the Saltation zone. These results (qualitative) show the ability of our model to describe in future simulations the interaction between the Saltation layer and the sand bed with obstacles, and moreover, to obtain the sand and air flow patterns around obstacles, including details about deposition and erosion.

3. Large-scale results over predicted the erosion as indicated by mass fluxes over 120%.

4. Small-scale results showed a sediment transport layer that is thin in comparison with the experimental data, which is susceptible to the conductivity, solid viscosity, distribution of granular temperature and the mean free path of the Hrenya and Sinclair model.

5. Based on these results, we can affirm that the GKT provides a relatively good description of the Saltation layer given the modification of the solid viscosity (due to the kinetic contribution) and conductivity, which greatly influence the mass flux, velocities and concentration.

Future work

The Aeolian sediment transport is influenced by many variables and, therefore, it is necessary to perform studies where the influences of many parameters are considered. The model proposed in this work allowed us to analyze the influence of the following parameters on the phenomenon: free stream velocity, slope of the sand bed, particle diameter and sand density. Also, we plan to extend the length of the computational domain in the x-direction in order to analyze the effects of particle accumulation on the mass flux and concentration; this could lead to better comparisons of simulated results with experimental data.

In addition, a more in-depth study of the effects caused by the granular pressure and temperature on the particle suspension should be conducted. It is expected that modifications in the granular temperature equation could improve the distribution of sand particles over the sand bed, based on consideration of the dissipation of granular energy and granular conduction.

Finally, including obstacles like cylinders, flat surfaces and aerodynamic profiles may enrich this study by evaluating their interaction with the sediment transport layer.

Acknowledgments

The author Juan P. Marval acknowledges the financial support given by the Fondo Nacional de Ciencia y Tecnologia (FONACIT, S3-2006000474) and by the Universidad Nacional Experimental Francisco de Miranda to cover the expenses during the research visit to the University of Florida. Also, authors want to recognize the valuable help of Eng. Basily Vasilacus by arranging the information and developing some figures and graphics.

References

1.Alhajraf, S. (2004). Computational Fluid Dynamic Modeling of Drifting Particles at Porous Fences. Environmental Modelling and Software 19, 163-170.         [ Links ]

2.Anderson, T.B. & Jackson, R. (1967). Fluid Mechanical Description of Fluidized Beds. Ind. Eng. Chem. Fundam. 6, 527–539.         [ Links ]

3.Andreotti, B., Claudin, P., Douady, S. (2002) Selection of dune shapes and velocities. Part 1: Dynamics of sand, wind and barchans. The European Physical Journal B 28, 321-339.         [ Links ]

4.Bagnolds, R. A. (1941). The Physics of Blown Sand and Desert Dunes, Chapman and Hall, London.         [ Links ]

5.Bauer, B., Houser, C., Nickling, W. (2004). Analysis of velocity profile measurements from wind-tunnel experiments with saltation. Geomorphology 59, 81-98.         [ Links ]

6.Butterfield, G.R. (1999). Near-Bed mass flux profiles in Aeolian sand transport: High-resolution measurements in a wind tunnel. Earth Surf. Process. Landforms 24, 293-412.        [ Links ]

7.Chapman, S. & Cowling, T. (1970). The Mathematical Theory of Non-Uniform Gases. Cambridge Mathematical Library, Cambridge, 3rd edition.         [ Links ]

8.Dong, Z., Liu, X., Wang, H., Wang, X. (2003). Aeolian sand transport: a wind tunnel model. Sedimentary Geology 161, 71–83.         [ Links ]

9.Dong, Z., Wang, H., Liu, X., Wang, X. (2004). The blown sand flux over a sandy surface: a wind tunnel investigation on the fetch effect. Geomorphology 57, 117-127.         [ Links ]

10.Drew, D. A. & Lahey, R.T. (1993) .In Particulate Two- Phase Flow. Butterworth-Heinemann, Boston, 509-566.         [ Links ]

11.Enwland, H., Peirano, E., Almstedt, A.-E. (1996). Eulerian Two-Phase Flow Theory Applied to Fluidization. Int. J. Multiphase Flow 22, suppl. 21-66.         [ Links ]

12.Garside, J. & Al-Dibouni, M. R. (1977). Velocity-Voidage Relationship for Fluidization and Sedimentation. Ind. Eng. Chem. Proc. Des. Dev. 16, 206.         [ Links ]

13.Gidaspow, D. (1994). Multiphase Flow and Fluidization, Academic Press Inc., Boston, 1st edition.         [ Links ]

14.Hrenya, C. & Sinclair, J. (1997). Effects of Particle-Phase Turbulence in Gas-Solid Flows. AIChE Journal 43, 853- 869.         [ Links ]

15.Hsu, T. J., Jenkins, J. T., Liu, P. L.-F. (2004). On Two-Phase Sediment Transport: Sheet Flow of Massive Particles. Proc. R. Soc. A 460, 2223-2250.         [ Links ]

16.Ishii, M. (1975). Termo-Fluid Dynamic Theory of Two- Phase Flow, Direction des Etudes et Recherches d`electricé de France, Eyrolles, Paris.         [ Links ]

17.Iversen, D. & Rasmussen, K. (1999). The effect of wind speed and bed slope on sand transport. Sedimentology 46, 723-731.         [ Links ]

18.Jackson, R. (1997). Locally Averaged Equations of Motion for a Mixture of Identical Spherical Particles and a Newtonian Fluid. Chemical Engineering Science 52, No. 15, 2457-2469.         [ Links ]

19.Jenkins, J. T. & Hanes, M. (1998). Collisional Sheet Flows of Sediment Driven by a Turbulent Fluid. Fluid Mech. 370, 29-52.         [ Links ]

20.Ji, S.B., Gerber, A.G., Sousa, A. (2004). A Convection- Diffusion CFD Model for Aeolian Particle Transport. Int. J. Numer. Meth Fluids 45, 797-817.         [ Links ]

21.Johnson, P. & Jackson, R. (1987). Frictional-Collisional Constitutive Relations for Granular Material, with Application to Plane Shearing. Journal of Fluid Mechanics 176, 67-93.         [ Links ]

22.Johnson, P., Nott, P., Jackson, R. (1990). Frictional- Collisional Equations of Motion for Particle Flows and their Application to Chutes. Journal Fluid Mechanic 210, 501-535.         [ Links ]

23.Liu, X. & Dong, Z. (2004). Experimental Investigation of the Concentration Profile of a Blowing Sand Cloud. Geomorphology 60, 371-381.         [ Links ]

24.Lun C. K. K., Savage, S.B., Jeffrey, D.J., Cherpuniy, N. (1984). Kinetic Theories for Granular Flow: Inelastic Particle in Couette Flow and Slightly Inelastic Particles in a General Flowfield. Journal of Fluid Mechanics 140, 223-256.         [ Links ]

25.Marval, J.P., Rojas, L., Curtis, J. (2007). Two-Dimensional Numerical Simulation of Saltating Particles Using Granular Kinetic Theory. FEDSM2007-37654. ASME/ JSME Fluids Engineering Division Summer Meeting. San Diego, California, USA.         [ Links ]

26.Mcewan, I. & Willets, B. (1993). Adaptation of Near- Surface Wind to the Development of Sand Transport. Journal of Fluid Mechanics 252, 99-115.         [ Links ]

27.Nalpanis, P., Hunt, J., Barret, C. (1993). Saltating Particles Over Flat Bed. Journal of Fluid Mechanics, 251, 661-685.         [ Links ]

28.Ni, J.R., Li, Z.S., Mendoza, C. (2002). Vertical profiles of Aeolian sand mass flux. Geomorphology 49, 205-218.         [ Links ]

29.Nishimura, K. & Hunt, J.C. (2000). Saltation and incipient suspension above a flat particle bed below a turbulent boundary layer. Journal Fluid Mechanics 417, 77-102.         [ Links ]

30.Ocone, R., Sundaresan, S., Jackson, R. (1993). Gas -Particle Flow in a Duct of Arbitrary Inclination with Particle-Particle Interactions. AIChE Journal 39 (8), 1261-1271.         [ Links ]

31.Ogawa, S. (1978). Multi-temperature Theory of Granular Material. Proc. US-Japan Seminar of Continuum         [ Links ]

32.Mechanical and Statistical Approaches in the Mechanics of Granular Materials, 208, Gakujutsu Bunken Fukukai, Tokyo.         [ Links ]

33.Okoli, R. (2003). Wind tunnel study Aeolian Saltation dynamics and mass flow. Journal of arid Environments 53, 569-583.         [ Links ]

34.Pasini, J.M. & Jenkins, J.T. (2005). Aeolian Transport with Collisional Suspension. Phil. Trans. R Soc. A. 363, 1625-1646.         [ Links ]

35.Pye, K. & Tsoar, H. (1990). Aeolian Sand and Sand Dunes. London: Unwin Hyman.         [ Links ]

36.Rasmussen, K. & Sorensen, M. (1999). Aeolian mass Transport near the Saltation threshold. Earth Surf. Process. Landforms 24, 413-422.         [ Links ]

37.Rasmussen, K. & Sorensen, M. (2005). Dynamics of particles in Aeolian Saltation. Powders & Grains. Balkema. Chapter report, 967-971.         [ Links ]

38.Patankar, S. V. (1980). Numerical Heat Transfer and Fluid Flow. Hemisphere,Washington,DC.         [ Links ]

39.Sauermann, G. (2001). Modeling of Wind Blown Sand and Desert Dunes. Ph.D. Thesis. Universidad de Stuttgart.         [ Links ]

40.Shaeffer, D G. (1987) Instability in the Evolution Equations Describing Incompressible Granular Flow. Journal Differ. Eqs. 66, 19.         [ Links ]

41.Schlichting, H. (1979). Boundary-Layer Theory, McGraw Hill, Seventh Edition.         [ Links ]

42.Sinclair, J. & Jackson, R. (1989). Gas Particle Flow in a Vertical Pipe with Particle-Particle Interactions. AIChE Journal 35, 1473.         [ Links ]

43.Sorensen, M. (2004). On the rate of Aeolian sand transport. Geomorphology 59, 53-62.         [ Links ]

44.Spies, P.J., Mcewan, I., Butterfiled, G.R. (1995). On wind velocity profile measurement taken in wind tunnels with saltating grains. Sedimentology 42, 515-521.         [ Links ]

45.Stout, J.E. (1998). Effect of averaging time on the apparent threshold for Aeolian transport. Journal of arid Environments 39, 395-401.         [ Links ]

46.Syamlal, M, Rogers, W., O´Brien, T. J. (1993). MFIX Documentation, Theory Guide. Technical Note DOE/ METC 94/1004. Morgantown, Wes Virginia.         [ Links ]

47.Ungar, J. E. & Haff, P. K. (1987). Steady-State Saltation in Air. Sedimentology 34, 289-299.         [ Links ]

48.Van Wachen, B.G.M., Shouten, J.C., Krishna, R., Den Bleek, Van, C. M., Sinclair, J. L. (2000). Comparative Analysis of CFD Models of Dense Gas-Solid Systems. AIChE Journal 47 (5), 1035-1051.        [ Links ]