Abstract
The flow field and the sand flow field constitutive equations are analyzed at first, then the different desert highway numerical models are established by considering the crossroad and by changing the road surface height and air stream flow field, then three kinds of different models with different complex air flow fields are made for simulating the sand ripple formation process by weak coupling of air and sand flow field, and finally the numerical simulations of these models are conducted and the affect process of sand morphology under complex air flow fields are discussed. The results show that under the uniform airflow field, the straight parallel ripple formed and the flared ripple formed in the middle region of the crossroad, and the wavelength of the ripples on the desert highway is bigger than that of the ripples around the road when the road height is higher than that of the sand surface height. Under the nonuniform complex airflow field, the complex curved ripples are formed, and some of the local area, where the whirlwind exits, no ripples are formed.
1. Introduction
In China, the size of the arid areas is more than 26000 km2; occupies 27% of our country’s land area. Besides the trend of desertification gradually increasing, this is more crucial problem and cannot be ignored [1]. Here rich in natural resources and the development potential are tremendous, but because of the lack of rain and drought, the ecological environment is fragile. Desert and desertification land are quite wide in China, from west Xinjiang to east Heilongjiang, intermittent located in the north of China arid, semiarid, and part of subhumid regions. Although the desert area is not suitable for human habitation because of its sparsely populated, desolate environment, but in the heart of its vastness is richness in coal, oil, natural gas, and other mineral resources [2]. In recent years, with the exploitation of mineral resources and the development of tourism resources, the speed and scale of desert highway construction in China had entered into the phase of rapid development, to a certain extent, solved the traffic problem in desert area. At the same time, the interaction between construction of desert highway and sand morphology formation process on the desert regions became a hot research topic in recent years.
Aeolian sand ripple is one of the most common sand morphology in the desert and formed on the sand bed surface by moving sand particles. In general, the geometric patterns of aeolian ripple have a rules, however it is affected by natural conditions. Its formation and evolution processes are close to air flow field characteristics and therefore its shape is not exactly the same [3].
In this paper, due to sand flow field real characteristics, the establishing process of air flow field constitutive equations is analyzed at first; then the implication relations and independency between air flow field and the sand flow field are analyzed. Then by means of numerical simulation of the airflow field and sand flow field coupling characteristics, the two desert highway intersections that formed sand ripples processes are studied and simulated under the uniform air flow field with different heights of the desert highway. Then three kinds of different models by considering the different complex air flow field are made for simulating the sand ripple formation process. Finally the numerical simulations of these models are conducted and the affect processes of sand morphology under complex air flow fields are discussed.
2. Theoretical Basis
In this section, the theoretical basis for airflow and sand flow field constitutive equations are discussed and illustrated.
2.1. Constitutive Equations of Airflow Field
Navier Stokes equation generalizes the universal law of viscous uncompressible fluid flow and is of special significance in fluid mechanics. Due to the compression, the airflow field has little effect on the sand flow field, where the airflow field could be considered as an incompressible fluid, and also the flow is assumed to be two-dimensional, laminar, and Newtonian in the simulating area. So, the governing mathematical model of airflow field can be expressed by the following set of mass and momentum conservation equations.
Momentum equation is
and the continuity equation is where is stream velocity components on the direction, is pressure, and is the time. Re is a dimensionless parameter known as the Reynolds number and expressed using , where is density, is viscosity, is the reference length dimension, and is the reference velocity dimension.
2.1.1. Initial and Boundary Conditions
The following initial conditions are used in all airflow field as shown in Figure 1:
The no-slip boundary conditions for rigid boundary and boundary conditions for the prescribed inflow boundary are used as follows: where the superscript and denote the normal and the tangential directions to the boundary, respectively, the subscript and denote the wall and inlet respectively, and is given velocity value. Usually the fixed solid boundary has both a rotational and/or a translational velocity, and then this could equally be described by no-slip conditions as wall. In this paper, for a no-slip rigid boundary (4) is used, for the inflow boundary in cavity flow problem (5) is used, and for the other prescribed inflow boundary (6) is used.
The boundary condition for the pressure can be derived from the Navier-Stokes equation and in no-slip boundary where the positive direction is normal to wall, the boundary condition for the pressure expressed a more general form as in the following equation:
2.2. Constitutive Equations of Sand Flow Field
The blown sand by wind in the air and sand flow field mainly appear as three kinds of types, as suspension, saltation, and creep [4, 5] as shown in Figure 1. Suspension is a long-term jumping process, generally the jumping distance is from 10 meters to thousand kilometers, it is the sand moving process far away from the sand bed surface and into the airflow field. Saltation is the short-term jumping process and the jumping distance is less than 10 meters; it is the sand moving process near the sand bed surface. The creep is wriggle process as reptile, no jumping but keeps close to and moving on sand bed surface. When the sand grains are carried up into the air by the wind and move along to any direction and reach somewhere in any distance, the power and velocity of the wind always cannot maintain the sand grains continuously in the air and finally the sand grains fall back to the ground somewhere. The salting grains travel by making long jumps and the repeating grains move in hopping manner over a much smaller distance [6]. In this paper, for computational simulation of the sand morphology forming process by wind, a new mathematical model [7] which is established based on the sand ripples model proposed by Nishimori and Ouchi [4] and Miao et al. [5] is used mainly to set out the two flow fields from the theoretical analysis and introduce the air flow field and sand flow field constitutive equation, and weakly coupling methods for air and sand flow field are analyzed by modifying sand ripple forming process by considering the uniform and nonuniform complex sand flow area.
2.2.1. The Saltation Dynamics Model of Sand Grains
The sand saltation can be described by using the following Newtonian moving equation [7, 8]: where is the initial position and and are the velocity and acceleration of the sand particles and the is gravity of sand particles. Assume that is a moving distance of sand particles on the time t by wind and is expressed as where the velocity and acceleration of the sand particles are expressed as where is base value of jumping distance of sand particles, which is jump along the directions and correlated with wind power. is the height of the grain position, is wind velocity, is sand grains critical wind velocity for jumping (in this papers use 0.2 m/s), and are the average mass and diameter of the sand grains, is grains shape coefficient (take it 1.0 for spherical), and are the density of the air and sand, and is resistance coefficient [7].
Take the equations from (9), (10) and (11) into (8) and introduce the relation of jumping time step and obtain the following equation: where is gradient along the directions.
Equation (12) expresses the main movement of the sand suspension and saltation form; due to limitation, the discussion of sand creep and avalanche phenomenon are ignored and the Nishimori and Ouchi [4] and Miao et al. [5] method are directly used. Some interesting and reasonable numerical results are obtained [7–10]. In this study similar numerical methods are used for airflow and sand flow field.
3. Numerical Method for Airflow and Sand flow Field
In this paper, the airflow field constitutive equations are discrete by using finite differential method with velocity correction method and the weak coupling between airflow field and sandflow field is used for implemtation of numerical calculation.
3.1. Coupling Method for Airflow and Sand flow Field
At the air-sand interface as shown in Figure 1, the actual airflow field and sand flow field are coupled by transition conditions and kinematic conditions. These coincide with the jump conditions for a contact discontinuity that are derived in continuum mechanics. Due to the explicit time stepping, the numerical fluxes have to be computed only once in the entire domain. Then one fixed point iteration only requires updating the numerical fluxes near the air-sand interface and the data in the cells attached to the interface. This makes the fixed point iteration very costly [11]. In this study, due to the air-sand interface complexity, the interaction effect of phase transition and surface tension is neglected, and the viscous effects in the air fluid also are neglected; then the simple and weak air-sand flow interface are coupled using the following simple condition: where the and are pressure of air and sand at interface , respectively; also and are flow velocity of air and sand at interface , respectively. On the other hand, this study only focuses on the sand formation under the nonuniform airflow field near the air-sand interface, so the weak coupling method is used to save computing time by the following manner.(1)The airflow velocities are assumed to be uniform along the height in the domain and air-sand interface is flat at the initial time step; therefore, only 2D numerical approach is needed for solving the airflow field for coupling sand bed.(2)Airflow field and sand flow field can be decomposed independently, and in each of these fields discretizations are applied, respectively, by keeping the same grids number at air-sand interface and directly exchanging the related information as velocity and pressure at corresponding gird points;(3)Only one airflow field which is obtained from steady state at sufficient time step with flat sand bed surface is performed for all time steps of the sand flow process.(4)The pressure and the velocity values at the sand bed surface are directly obtained from the 2D airflow field at steady state on the corresponding grids.(5)Numerical simulations of sand flow are conducted independently by using stable and constant airflow field on every unique time step of sand flow process.
3.2. Numerical Method for Airflow Field
In this paper the staggered Marker And Cell (MAC) method is used [12, 13] to discretize the velocity and pressure fields used in the Navier-Stokes equations as described in (1) and (2). Figure 2 shows an example of a staggered grid and the arrangement of staggered grid provides stronger coupling between the pressure and velocity variables thus improving stability. The solid lines in Figure 2 represent the original primary grid, while the dashed lines means secondary grid. The pressure is assigned to the nodes on the primary, while the velocities are defined on the secondary nodes.
The merit of this method is that the pressure is defined on the primary lines crossing point, and the velocity is defined on primary and secondary lines crossing points. So the pressure can be evaluated by using the and on primary grids, and the velocity can be evaluated by using the or on the secondary grids along the direction, while the velocity evaluated by using the or on the secondary grids along the direction. As the simulation progresses, the marker particles are moved through the velocity field and then used to determine which cells contain fluid.
Consider
The Euler explicit scheme by ignoring external forces is used for solving two-dimensional airflow field problems, and the temporal flow velocity and are calculated by using Finite Differential Method (FDM) as equation (14) by adding −1/2 to the subscript for easy coding.
The MAC method actually solves for a relative quantity called velocity correction. Figure 3 shows the flow chart of MAC algorithm. At first, the initial flow field and pressure distribution in the domain are set, and set of momentum and continuity equations are coupled and its nonlinearity of the equations are solved by using iteration method as shown in Figure 3. The pressure field is assumed to be known from the previous iteration. Using this, the flow velocities and in momentum equations are solved by using differential equation (14), and at this stage the newly obtained velocities and usually do not satisfy the continuity due to the pressure field being set by assumption. So the corrections to flow velocities and pressure are proposed by using (15) to satisfy the discrete continuity equation with the appropriate error range.
Consider where , , and are expressed using the following equation: Then the corrections of velocities are solved, with the correction of pressure gradient term by using the following equation:
3.3. Validation of Numerical Method for Airflow Field
For the validation of the numerical method, the uniform grids are considered, where the width and the height of a grid cell are fixed over the entire domain. For the simplification of the coupling problem between the air and sand flow field, the size of grid cells are set to constant and equal each other in all dimensions. Simulation results for incompressible laminar flow in square cavity are modeled and the results are compared to the literature results [14] to verify the validity of the simulation solution. No-slip velocity boundary condition () is applied on all the walls, except the top lid for the cavity flow. On the top lid ( and ) is applied. The bottom boundary of the domain is modeled as wall. The Reynolds number is set as 100, 500, and 1000. The same grid cell size 0.002 along the x and y are used, and the numerical domain is divided in a total of 250000 cells. The time step is set as 0.00001 s and the errors Ew are controlled in Er = 0.0001 in every time step by using the iteration method, also the numerical calculations are conducted until the cavity flow is at steady state, where the errors Et are over the prescribed value Etr.
Figures 4(a) and 4(b) show the numerical results of Ux velocity profiles with different time steps and Reynolds number along a vertical line, respectively, passing through the geometric center of the cavity flow. It is clearly seen in Figure 4(a), the Ux velocity profiles strongly depend on the numerical processing time step, and when time steps equal 20 s, the velocity profile goes to steady state and any increase beyond these time steps would lead to an insignificant change in the resulting solution. Figure 4(b) shows the Ux velocity profiles with Reynolds number equal to 100, 500, and 1000, respectively, at steady state, and the numerical results of velocity profiles in this study close to that of the literature [14, 15] results.
(a)
(b)
4. Numerical Simulations of Sand Ripple Formation
The following simulation only considers three main movement of the sand form, namely, saltation, creep, and avalanche under the air flow field.
4.1. Sand Ripple Formation under the Uniform Airflow Field
The dimensionless parameters are used for simulation in all cases as follows. The calculation grids are 1000000 with 1000 grids both along the - and -axis, including 360000 grids for highway area just as with 200 grids widths and 1000 grids long. Initial sand bed surface roughness range is 1.5–2.0, sand ripples initial wavelength is 10, and sand grains mass range is between 1.0 and 1.5 and its diameter range is between 0.5 and 1.0. The initial roughness of sand bed surface, sand grain mass, and diameters are generated in random manner by using its range. The average clustered sand grain field density (2.65 g/cm3) and air stream field density (0.00129 g/cm3) are assumed to be constant, and the gravity acceleration of cubic clustered sand field also set a constant as −981 cm/s2. Airflow velocity is always a uniform constant quantitative 1 from left to right.
Figure 5 shows the numerical results of the highway with crossroad under the uniform air flow field. When the highway height is lower than height of sand bed surface as height , as shown in Figures 5(a) and 5(d), it is shown that with the increasing of calculation steps the ripples are formed on the sand bed and gradually into the highway. When the calculation steps are at 2000th as shown in Figure 5(d), the whole road is covered by sand and ripples are formed. The ripples shape and size nearly close to sand bed ripples shape and size.
(a)
(b)
(c)
(d)
(e)
(f)
When the highway height is the same as the sand bed surface as height = 0, it is clearly shown in Figures 5(b) and 5(e), with the increasing of calculation steps, the crossed road is also covered by sand gradually, and sand grain coverage rate in the road, where the wind direction is perpendicular to it, is faster than that of the parallel one. When the calculation steps are at 2000th as shown in Figure 5(e), the road, where the wind direction is perpendicular to it, is completely covered by sand and the sand ripples are formed, the thin sand areas are formed at the intermediate on the crossroad where the wind direction parallel to the road, and thick sand areas are formed on both sides of the crossroad. But overall the sand cover rate is slower than when highway height is −2 as shown in Figures 5(a) and 5(d). This conclusion is consistent with other scholars who study the impact of the fixed area of sand ripples. However, the interesting phenomenon is at this time it can be found that the wavelength of the ripples on the crossroad is slightly larger than that of the wavelength of the ripples in the others areas.
Figures 5(c) and 5(f) show that the highway height is higher than height of sand surface as height . With the increasing of calculation steps, the ripples are formed, but overall the sand cover rate is slower than when highway height is −2 and 0. When the calculation step is at 2000th as shown in Figure 5(f), the flared ripple area is formed clearly on the crossroad which is parallel to the wind direction and the regional erosion by sand is far less than another area. It is also very interesting that when the road height is higher than that of the sand bed surface height, the wavelength of the ripples on the crossroad get bigger than that of the sand bed surface; this is may be due to sand particles easily moving on the road surface.
4.2. Sand Ripple Formation under the Complex Airflow Field
In above simulation, we only considered uniform airflow field. However, due to the complexity of actual desert atmosphere, geomorphology, and vegetation environment, the desert airflow field is nonuniform and quiet complex. Figure 6 shows the field wind tunnel tests with semicircular obstacle, the solid arrow lines expressed wind direction, the dashed lines expressed the sand ripples ridge, and the long dash dotted line means the different sand ripples boundary. The results show different types of ripples such as arced ripples and parabolic ripples are obtained experimentally. Very flat no sand ripples area clearly appears in Figure 6(b), and in this no sand ripples area a little whirlwind can be observed one-by-one as shown in Figures 6(a) and 6(c), and this may be one reason why in this area no sand ripples appear. So, the airflow field nonuniformity for sand ripple formation process should be considered in numerical simulation.
In this paper, for more realistic consideration of sand flow field into the actual simulation of the sand ripples formation process, three kinds of special models are made by changing the boundary condition with and without obstructs as shown in Figure 7, the air flow field obtained by using the numerical solution of Navier-Stokes equation.
(a)
(b)
(c)
Figure 7 shows the initial condition of these three kinds of models, and Figure 7(a) is a cavity flow model, Figure 7(b) is a highway model with one side uniform inlet velocity and one part uniform outlet velocity, and Figure 7(c) is a square barrier or obstruction model with one part uniform inlet velocity and outlet velocity.
For the cavity flow model, no-slip velocity boundary condition () is applied on all the walls, except the top lid. On the top lid ( and ) is applied. The bottom boundary of the domain is modeled as wall. The boundary conditions describing the current simulated computational domain as well as the surface boundary layer are depicted in Figure 7(a).
For the highway and obstruction model, no-slip velocity boundary condition () is applied on all the walls, except the left upper side and right down side. On the left upper side and right down side, ( and ) is applied. The bottom boundary of the domain is modeled as wall. The boundary conditions describing the current simulated computational domain as well as the surface boundary layer is depicted in Figures 7(b) and 7(c).
The numerical results are shown from Figure 8 to Figure 10. The velocity field of cavity flow model with Re = 100 shown in Figure 8(a), Figures 8(b), and 8(c) are sand bed morphology at 20th and 130th steps, respectively. Although erosion can be observed in this cavity flow model, but it did not generate any sand ripples near the vortex area during the whole sand flow process, some other areas can be found as some heaped up sand, but cannot confirm sand ripples.
(a)
(b)
(c)
Figure 9 shows the sand ripples formation under the highway model with one part of side uniform inlet velocity and another part of opposite side outlet velocity air flow field with . It is shown that the air velocity distribution in Figure 9(a) is not so uniform but nearly smoothed from inlet to outlet, and vortex area is not so clear. The sand bed morphologies at 20th and 100th steps are shown in Figures 9(b) and 9(c), respectively. The complex sand morphology formed in this model including complex shaped ripples and small sand dune is clearly seen.
(a)
(b)
(c)
(a)
(b)
(c)
Figure 10 shows the sand ripples formation under the square barrier (obstructs) model with one part of side uniform is inlet velocity and another part of opposite side is free outlet and the . It is shown that the air velocity distribution in Figure 10(a) that is strongly complex and vortex can be found in area and , respectively. Figures 10(b) and 10(c) are the sand bed morphologies at 20th and 130th steps, respectively. It is clearly seen that the complex sand morphology were also formed in this model including complex shaped ripples and small sand dune, and in areas and some erosion and deposition process of sand can be found but cannot confirm sand ripples.
By using nonuniform airflow field, the different results with uniform fields are obtained and some phenomenons close to Figure 6 are observed. The simulation result has improved compared with simulation result of uniform and nonuniform airflow field, and there appear different extent erosions around obstruction.
5. Conclusions
The sand morphology formation processes are simulated under the uniform and nonuniform airflow condition using the weak coupling method, and the following conclusions are obtained.
Under the uniform flow field, the flared ripple area in the middle region is formed, and the horn mouth of the flared ripple area toward the wind direction on the crossroad. Through the analysis it is clear that the highway roadbed height strongly affects the sand deposition process. Wavelength of the ripples on the crossroad is bigger than that of the sand surface when the road height is higher.
Under the nonuniform complex flow field, the complex sand shapes are formed as curved ripples, small sand dune, and flat sand bed surface without ripples and the dynamic phenomenon of sand ripples formation process close to the experimental observation.
In this paper, only one-directional weak coupling of air and sand flow field is considered. However, the sand ripples shape and air flow field would ceaseless affected each other, this will be next target in the future research.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work was supported by the National Natural Science Foundation of China (11102178, 10010701, and 10062002).