Discrete Element Simulation on Sand-Bed Collision Considering Surface Moisture Content

Abstract

The process of aeolian sand transport is an important mechanism leading to the formation and evolution of local landforms in coastal areas and desert lakes. For a long time, the role of surface moisture in incipient motion of sand grains by wind stress has been extensively studied but, in fact, sand-bed collision is the main mechanism in steady aeolian sand flow. At present, the lack of understanding of surface moisture content on sand-bed collision limits the application of aeolian sand transport models in wet coastal areas. In this paper, we adopt numerical simulations to discuss and analyze the effect of cohesive forces formed by surface moisture content on the sand-bed collision process based on discrete element method. High density contact forces appear with the surface moisture increasing, and form a closed structure around the edge of crater to resist the avulsion in horizontal direction. Under high moisture condition, even though the ejected sand grains saltate away from the surface, the tension forces will prevent from leaving. The ejected number trend with incident velocity shows some nonlinear characteristics due to the unequally distributed force chains and liquid bridges in the unsaturated sand bed surface.

1. Introduction

Coastal aeolian sand transport is a first-order control on growth of local dunes and ripples, shaping the landforms and resulting in the generation and evolution of desertification [1,2,3,4]. The temporal variability of surface moisture on beaches plays a critical role in the process of aeolian sand transport, and the incorporation of cohesive forces formed by surface moisture into wind blown sediment transport models is an important issue in the field of aeolian physics [5,6,7]. The role of surface moisture in the incipient motion of sand grains by the wind has been extensively studied, and numerous empirical and physical models have been proposed to quantify this effect [8,9,10,11,12,13,14,15,16]. However, the sand-bed collision formed by saltation is the main mechanism in steady aeolian sand flow, in which sand grains move forward by nearly ballistic trajectories under the forces of gravity and wind drag, collide with the loose sand surface and eject new sand particles [17,18] (Figure 1). The almost complete lack of knowledge about the effect on sand-bed collision precludes the deep understanding of aeolian physics and improving numerical simulation accuracy.

A relatively small number of researchers have been reported to explain the saltation behaviour considering surface moisture in field and wind tunnel environments [19,20,21,22,23,24,25,26,27]. Many researchers suggest that small amounts of surface moisture content restrain the saltation in two ways: reducing the number of incipient motions of sand particles due to wind stress and decreasing the total number of particles ejected from the sand-bed collision subsequently. However, surface moisture also enhances the cohesive forces among sand particles and produces a harder surface, which results in sand particles saltate higher with increasing moisture content and ultimately increasing the aeolian sand transport rate [19,20]. Until now, none of the current aeolian sand models can express the effect of surface moisture on sand-bed collision in saltation under coastal conditions, where enhanced cohesive force due to high moisture content appears to be a critical parameter and greatly reduces the range of application of simple transport formulations [2].

The wetted coastal surface is composed of discrete sand grains which are drawn together by a large amount of interparticle cohesive forces. Recently, significant studies have been undertaken to develop closed-form formulations for various cohesive forces and other mechanical behaviours at the solid–liquid interface, such as auto-adhesion, friction, liquid binding and so on [28,29,30,31,32,33], but knowledge of their impact on the macroscopic behaviour of the wetted sand surface is still lacking. From the viewpoint of the microscopic surface structure, the movements of collided sand particles are determined by the cohesive forces among them, which can be presented as functions of their relative positions and velocities. In consideration of this multi-body system involved in the particles collision and energy dissipation, a numerical simulation scheme is required to reproduce this complex process.

As a widely used numerical simulation method of a dense multi-body system, the Discrete Element Method (DEM) can resolve the mechanical behaviour of sand grains under consideration of various relevant cohesive forces acting on individual grains [34,35,36]. Lian et al. [37] simulate the pair-wise collisions of wet agglomerates considering cohesive forces supplied by liquid bridges among neighbouring particles, and discuss the effect of impact velocity and fluid viscosity on the energy dissipation and agglomerates broken processes. Yin et al. [38] investigate the sand-bed collision under inclined bed conditions, and offer a quantified formula to describe the number and velocity of ejected sand grains as a function of the impact grain velocity.

The present paper investigates this sand-bed collision using the Discrete Element Method, under consideration of the sand grain size distribution associated with aeolian sand in nature and the cohesive forces formed by liquid bridges among neighbouring sand grains. Our aim is to explore the evolution of the movement behaviour of sand grains in the air and the microstructure under the sand surface, as characterized by the force chains in the dense multi-body system.

2. Model Description

2.1. Cohesive Force Equations

In this work, we consider the interstitial liquid bridges among sand particles as the main form of cohesive forces. All discussion about the inviscid liquid bridge must be limited to a certain scope, the critical separation distance $S_c$; it is a limit value that can maintain the existence of liquid bridge [39]

$$S_c=(1+0.5\alpha )V^{1/3},$$

(1)

where V is the liquid bridge volume and $\alpha$ is the contact angle of solid and liquid.

Generally, the liquid bridge exists in a form of pendular, funicular or capillary according to degree of saturation, and the attractive capillary force is the primary mechanism of cohesive force when the relative velocity between two particles is small [30,40]. The capillary forces are static and can be solved by the Laplace–Young equation using the reduced Laplace pressure $\Delta p$ and the surface tension $\gamma$ of liquid bridge [9,11]. Figure 2 shows the the liquid bridge arc radius ${\rho }_1$ and the neck radius ${\rho }_2$, the geometry formulas can be expressed as

$${\rho }_1=\frac{s_i+R_i(1-cos{\varphi }_i)}{cos({\varphi }_i+\alpha )},$$

(2)

$${\rho }_2=R_{i}sin{\varphi }_i-{\rho }_1[1-sin({\varphi }_i+\alpha )],$$

(3)

where $i=1,2$ presents the serial number of two spheres, $s_i$ is the half-separation distances, $R_i$ is the radii and ${\varphi }_i$ is the half-filling angles.

The mean curvature of the liquid bridge at the neck H can be shown as

$$H=\frac{\Delta p}{2\gamma }=\frac{({\rho }_2-{\rho }_1)}{2{\rho }_2{\rho }_1},$$

(4)

and the static capillary force can be simply given by

$$F_c=2\pi \gamma {\rho }_2(1+H{\rho }_2),$$

(5)

2.2. Particle Motion Equations

There is a broad range of inter-particle contact force models for application in DEM, e.g., the linear spring model, HertzMindlin and JKR model, and linear spring–dashpot model. The linear spring–dashpot model is adopted in this study as it has proven to accurately represent collisional behaviour in the framework of particle-based simulations of Aeolian transport [38]. Each sand particle’s motion during sand-bed collision can be solved by Newton’s equations, and the translation and rotation of each grain can be expressed as

$$m_i\frac{d^2{x}_i}{d{t}^2}=\left(\sum _{j\ne i}{\mathit{F}}_{n,ij}+\sum _{j\ne i}{\mathit{F}}_{\tau ,ij}\right)·{\mathit{e}}_1,$$

(6)

$$m_i\frac{d^2{y}_i}{d{t}^2}=\left(\sum _{j\ne i}{\mathit{F}}_{n,ij}+\sum _{j\ne i}{\mathit{F}}_{\tau ,ij}+{\mathit{G}}_i\right)·{\mathit{e}}_2,$$

(7)

$$I_i\frac{d^2{\theta }_i}{d{t}^2}=\sum _{j\ne i}\left({\mathit{r}}_{0,ij}\times {\mathit{F}}_{\mathit{\tau },ij}\right)·{\mathit{e}}_3,$$

(8)

where $x_i$ and $y_i$ present the centroid coordinates of grain i, $m_i$ is its mass, ${\mathit{e}}_1$ and ${\mathit{e}}_2$ are unit vectors in the directions x and y, and ${\mathit{e}}_3={\mathit{e}}_1\times {\mathit{e}}_2$. $I=2m{R}^2/5$ is the rotation inertia, ${\theta }_i$ is the angle of rotation, ${\mathit{r}}_{0,ij}$ is the direction vector from the centroid of grain i to j.

${\mathit{F}}_{n,ij}$ and ${\mathit{F}}_{\tau ,ij}$ are the normal and tangential forces, respectively, acting on grain i due to contact with grain j. The expression for the contact forces read,

$${\mathit{F}}_{n,ij}=\left\{\begin{array}{cc}-k_n{\delta }_{ij}{\mathit{n}}_{ij}-{\beta }_n{\dot{\delta }}_{ij}{\mathit{n}}_{ij},\hfill & {\delta }_{ij}\ge 0,\hfill \\ -{\mathit{F}}_{c,ij},\hfill & -S_c\le {\delta }_{ij}<0,\hfill \\ 0,\hfill & {\delta }_{ij}<-S_c,\hfill \end{array}\right.$$

(9)

$${\mathit{F}}_{\tau ,ij}=\left\{\begin{array}{cc}-\mu |{\mathit{F}}_{n,ij}|{\mathit{t}}_{ij},\hfill & |{\mathit{F}}_{t,ij}|\ge \mu |{\mathit{F}}_{n,ij}|,\hfill \\ -k_t{\xi }_{ij}{\mathit{t}}_{ij}-{\beta }_t{\dot{\xi }}_{ij}{\mathit{t}}_{ij},\hfill & |{\mathit{F}}_{t,ij}|<\mu |{\mathit{F}}_{n,ij}|,\hfill \end{array}\right.$$

(10)

where ${\delta }_{ij}$ and ${\xi }_{ij}$ are the grain deformation in the normal and tangential directions to the contact plane associated with the inter-particle collision, respectively. $\dot{{\delta }_{ij}}$ and $\dot{{\xi }_{ij}}$ are derivatives of ${\delta }_{ij}$ and ${\xi }_{ij}$ with respect to time, ${\mathit{n}}_{ij}$ is unit vector from the mass centre of grain i to j, and ${\mathit{t}}_{ij}$ is the unit tangent vector. ${\mathit{F}}_{c,ij}$ and $S_c$ are the capillary force and critical separation distance mentioned above. $k_n=E·min(R_i,R_j)/(R_i+R_j)$ and $k_t=k_n/\left(\lambda \nu \right)$ are the moduli of elasticity of normal and tangential overlap, E is the Young’s modulus, $\nu$ is the Poisson’s ratio, and $\lambda$ is the empirical coefficient [41]. ${\beta }_n=\sqrt{4k_n\frac{m_i{m}_j}{m_i+m_j}\frac{{ln}^2\epsilon }{\sqrt{{ln}^2\epsilon +{\pi }^2}}}$ and ${\beta }_t=\sqrt{4k_t\frac{m_i{m}_j}{m_i+m_j}}$ are the normal and tangential damping coefficients, respectively. $\epsilon$ is the coefficient of restitution and $\mu$ is the coefficient of friction.

2.3. Simulation Procedures

First, a mixed sand granular bed should be generated before the simulation of a sand-bed collision process. In order to obtain a sand bed as similar to nature as possible in the simulation, we use the random number generator to generate the mixed granular system [42]. A Gaussian distribution $f=p\left(d\right)$ is chosen to character the probability density of grain diameter d [43,44], as shown in Figure 3. In total, 5000 grains are dropped into the box and the grain cooling technique to make the bed surface to become stable under gravity and cohesive forces [45]. We generate a flat sand surface of 100 <d> long and 50 <d> high until the maximum velocity of each sand particle is less than ${10}^{-5} \mathrm{m}/\mathrm{s}$.

Next, a incident sand grain is shot to sand bed at initial 10 <d> high, and the incident angle of sand grain ${\theta }_{in}$ is set as $11.5^{\circ }$, which is the most frequent in aeolian sand transport [38]. A value of 0.05 s was chosen as the final moment of calculation, because the ejected sand grains saltate to the air and the distances among them are big enough to avoid mid-air collision at that moment. Information about the position, velocity and force of each sand grain is recorded during the numerical simulation process to analyse the movement behaviour, microstructure, energy dissipation and statistical characteristics of sand grains.

Finally, we carried out 280 simulations; the parameters are incident velocity $V_{in}\left(\mathrm{m}/\mathrm{s}\right)$ (1, 2, 3, 4, 5, 6, 7) and the mass moisture content $w\left(\mathrm{g}/\mathrm{g}\right)$ (0, 0.01, 0.02, 0.03), each case being repeated 10 times for statistical analysis. The values of moisture content used in this work can ensure that the water is in an unsaturated state under the surface and form liquid bridges and attractive capillary forces among sand grains [46]. The liquid bridges are evenly distribute between neighbouring sand grains [37].

A particle flow analysis program Particle Flow Code (PFC) is used to simulate the sand-bed collision process in this work. The saltation sand grains in aeolian sand transport process move mainly downwind and vary in non-obvious way along the spanwise direction, which can be approximated as a two-dimensional motion [47]. By comparing with the experimental results, the 2D simulation method has proved to be highly accurate and are widely performed in the numerical simulation of sand-bed collision [34,35,38], so do this work.

The mechanical properties parameters of sand grains for computer simulations were set as sand grain density ${\rho }_p=2650\mathrm{kg}/{\mathrm{m}}^3$, Poisson’s ratio $\nu =0.3$, Young’s modulus $E=70\mathrm{GPa}$ and the coefficient of friction $\mu =0.3$. For properties of the liquid bridge, liquid density ${\rho }_p=1000\mathrm{kg}/{\mathrm{m}}^3$, contact angle between liquid and solid $\alpha =0^{\circ }$, surface tension $\gamma =0.025\mathrm{N}/\mathrm{m}$ and fluid viscosity $\eta =10\mathrm{m}\mathrm{Pa}\mathrm{s}$. The volume of each liquid bridge can be calculated by the volumes of two spheres as $V=0.4/3\pi (R_1^3+R_2^3)$.

3. Results and Discussion

In this section, our aim is to explore the evolution of the movement behaviour of sand grains in the air and the microstructure under sand surface. Meanwhile, we also compare the numerical simulation results with theoretical and experimental data under the dry sand surface condition.

3.1. Motion and Force Characteristics

Figure 4 presents the motion behaviours of sand grains in sand-bed collision under different surface moistures in time increments of 0.15 s, and we emphatically introduce the $V_{in}=5\mathrm{m}/\mathrm{s}$ cases because they happen most often in nature. Information about the position and velocity of each sand grain is shown in the form of circles and colour. The incident grains move from the top left to the bottom right direction in the box, collide roughly in the centre of the sand surface bed and rebound into the air on the right side. A large amount of sand grains originally belong to the surface are ejected into the air, and the total number decreases with the increasing surface moisture. This is reasonable because the combination of moisture and sand grains produces a dense and hard system, resulting in the underlying saltation cascade phenomenon observed by numerous experimental studies [20,21,23,25,26,27]. Simultaneously, largish impact craters are generated near the collision positions similar to craters on Earth or Mars. It is interesting to note that the shape of crater is wide and shallow when the surface moisture is low, and higher moisture content significantly decreases the area and slightly increases the depth of the crater. The shape and volume changes of crater not only affect the ejected grain number from there but also turn the rebound direction of incident grain, which may be an important influencing factor for the structure of drifting sand flux [20,23].

In order to investigate the deep reason for the changes of crater affected by moisture, we draw the distribution of contact forces among sand grains as shown in Figure 5, and the simulation cases are same with Figure 4. Contact forces are represented as a form of lines and the colour indicates magnitudes of forces. Although the positions of sand grains are not shown to make a clear figure, but they can be estimated by observing the endpoints of the segment. The shape and volume changes of crater also can be seen through the partial loss of contact force; however, some very important additional information is provided in this state. The large contact force points are distributed near the bottom of the box in all cases and the small contact force points beside the crater are key influence factors. The contact forces are low density and uniform distributed in low moisture cases, and seem unlikely to build up a strong resistance. High density contact forces appear with the surface moisture increasing, and form a closed structure around the edge of crater, which inevitably has a better performance to resist the avulsion in the horizontal direction.

3.2. Force Chains and Energy Dissipation Characteristics

In this part, we aim to research the effect of cohesive forces formed by liquid bridges. Figure 6 presents the distribution of force chains among sand grains under different surface moistures and the compression and tension forces are differentiated by colour. The positions of sand grains can be estimated by observing the endpoints of the segment also in this figure. It is obvious that there are only compression forces in the dry sand surface and the force chains are sparse and cluttered around the edge of crater. As the moisture content increases, the average length of tension force chains under the surface increase significantly, indicating the sand grains are pulled more tightly. Especially around the crater, the tension force chains extend outward by connecting the ejected sand grains. Hence, even though the ejected sand grains saltate away from the surface and generate a considerable distance from the adjacent particles, the tension forces will pull them from leaving until the liquid bridges are stretched to break. In others word, under high moisture condition, the ejected sand grains may be unable to consider as a liftoff state as usual, because they can be pulled back to ground soon. This is a very important mechanism in sand-bed collision, which has been neglected for so long. It explains the reason why current aeolian sand models are difficult to generalize to coastal areas. At the same time, higher moisture content supports larger liquid bridge volumes, enhances the cohesive forces and extends the critical separation distance $S_c$ described by Equation (1).

Figure 7 shows the energy dissipation in sand-bed collision under different surface moisture values in time increments of 0.15 s, and we define $P=m_{in}{V}_{in}^2/{\sum }_{i=1}^n{m}_i{V}_i^2$ as ratio of the kinetic energy of all the grains before and after the collision, where $m_{in}$ and $V_{in}$ are the mass and velocity of incident grain before collision, respectively, and the $m_i$ and $V_i$ are the masses and velocities of all grains after collision, respectively, n is the total number of sand grains. A dimensionless parameter $V_i/{(g<d>)}^{1/2}$ is used to represent the velocity of incident grain, the g is acceleration of gravity and $<d>$ is the average grain diameter. When the incident velocity is low, the increasing surface moisture intensifies the energy dissipation due to the enhanced cohesive forces. However, the trends almost reverse when the incident velocity is high; high surface moisture content resists the energy dissipation. In fact, the incident grain carries most of the kinetic energy in whole multi-body system, and its dynamic behaviour in collision determines the change of kinetic energy. The incident grain can release more energy to break the cohesive forces under low moisture condition, but it will be rebounded by wet and hard surface by keeping most of the kinetic energy.

3.3. Statistical Analysis of Sand Grains

On account of the key role of incident grains in whole multi-body system, the restitution coefficient should be the focus of the study, which affects most of the kinetic energy. Figure 8 presents the restitution coefficient $e_r$ of incident grains in sand-bed collisions under different surface moisture in times of 0.15 s. We first compare our numerical simulation results with previous experimental results under a flat surface condition. Simulation results show that the $e_r$ of incident grain distributes in the range from 0.5 to 0.6 when the surface is dry, and is consistent with the previous studies from Anderson et al. [34] and Chen et al. [47]. The reason that our $e_r$ results smaller than that of Zhou et al. [35] may be the different conditions and materials. The $e_r$ of of incident grain in multi-body system is not a constant just related to the properties of the system; in fact, this is a growth variable with the increasing incident velocity. The variation law of restitution coefficient is similar to that of energy dissipation. Interestingly, there is a slight decrease in restitution coefficient with the increasing moisture content under the low moisture condition. This can explain why the effect of the tension force between the incident grain and sand surface in the collision increases faster than the resistance effect in that case.

The average velocity of ejected grains $V_{ej}(\mathrm{m}/\mathrm{s})$ is another important value determining the structure of drifting sand flux, and it must remain constant with incident velocity to satisfy momentum and energy conservation. As shown in Figure 9, the $V_{ej}$ by the incident grain in sand-bed collision under different surface moisture values presents a relatively flat trend with the increase in incident velocity, indicating they are independent of each other. In the dry surface case, the $V_{ej}$ can reach 0.3 in almost all the incident velocity cases, and it is similar to previous experimental and simulation results [34,35,47]. With the surface moisture increases, the $V_{ej}$ drops down to a range of 0.1–0.15 $(\mathrm{m}/\mathrm{s})$, indicating that the surface moisture has a significant effect on decreasing the kinetic energy of ejected grains.

Figure 10 presents the ejected number $N_{ej}$ by incident grain in sand-bed collision under different surface moisture in times of 0.15 s, and we can see the $N_{ej}$ increases linearly with the incident velocity in the dry case, which has been sufficiently studied by theory and numerical methods [34,35]. Take the surface moisture content into consideration, the trend of $N_{ej}$ with incident velocity shows some nonlinear characteristics. It could have been caused by the unequally distributed force chains and liquid bridges in the unsaturated sand bed surface as shown in Figure 6. The average velocity and ejected numbers significantly decrease when the sand surface becomes wet, and constant with numerous experimental studies [20,21,23,25,26,27]. This work reproduces the process of surface moisture restraining the saltation, in which the liquid bridges increases kinetic energy consumption of ejected grains and decreases the total number of particles ejected from the sand-bed collision.

The information about statistical significant are also given in Figure 8, Figure 9 and Figure 10. It is clear that the standard deviations of $e_r$, $V_{ej}$ and $N_{ej}$, shown as error bars, all follow a decreasing trend generally with the increasing velocity of incident grains. In the low incident velocity cases, less kinetic energy is available for collision and the randomness due to impact location on mixed loose surface is significant. The randomness reduces with the increasing incident velocity and this change law has been confirmed by experiments [38]. In low energy environments, the surface anisotropism of dense multi-body systems enhances the randomness of the sand-bed collision and this effect should be considered in aeolian sand transport models.

4. Conclusions

In this paper, we investigate sand-bed collision using the Discrete Element Method considering the cohesive forces formed by liquid bridges among neighbouring sand grains, and explore the evolution of the movement behaviour, microstructure, energy dissipation and statistical characteristics of sand grains.

Morphologically, largish impact craters are generated near the collision positions, and higher moisture content significantly decreases the area and slightly increases the depth of the crater. High density contact forces appear with the surface moisture increasing, and form a closed structure around the edge of teh crater, which inevitably has a better performance to resist the avulsion in the horizontal direction. Those different morphologies and microstructures of craters give us a unique insight into the deformation resistability capacity of surfaces consisting of loose grains under different moisture contents.

Under high moisture conditions, even though the ejected sand grains saltate away from the surface, the tension forces will pull them from leaving until the liquid bridges are stretched to break. This is a very important mechanism in sand-bed collision, which has been neglected for so long. Higher moisture content enhances the cohesive forces, extends the critical separation distance and resists the energy dissipation. The simulation results reveal how the liquid bridges affect the evolution of movement behaviour of sand grains in the coastal aeolian sand transport process.

The effect of tension force between incident grains and sand surfaces in collision increases faster than the resistance effect in the low moisture case, resulting in the slight decreases in restitution coefficient with the increasing moisture content. The eject velocity presents a relatively flat trend and is independent of incident velocity, and the ejected number trend with incident velocity shows some nonlinear characteristics due to the unequally distributed force chains and liquid bridges in the unsaturated sand bed surface. Statistical results of moisture effects can be used in the future to perfect the aeolian sand transport model and improve calculation accuracy essentially in wet coastal areas.