Analysis and Application of Particle Backtracking Algorithm in Wind–Sand Two-Phase Flow Using SPH Method

Abstract

Due to the high sensitivity of grid-based micro-scale wind–sand flow models to deformation and distortion, this study employs the Smooth Particle Hydrodynamics (SPH) method for numerical simulations. The advantage of the SPH method is that it can dynamically analyze the entire trajectory of the particles, thus allowing the initial positional distribution of sand-buried particles to be traced. This study utilizes the advantages of the SPH method. It develops particle backtracking algorithms based on the SPH method using the C language. It analyses the initial location distribution, concentration, velocity, and particle size distribution of sand-buried particles to formulate targeted measures to cope with wind–sand disasters. Meanwhile, this paper improves a particle modeling algorithm to realize arbitrary mixing particle size and mixing ratio by programming in C language and combining it with pixel recognition technology. In addition, this paper will use the particle backtracking algorithm to analyze the classical embankment wind and sand flow field and then propose adequate measures for embankment wind and sand disaster management by investigating sand particle movement characteristics.

1. Introduction

Global deserts and desertification lands cover more than 7 × 10⁶ km² areas of the Earth. Since the 1990s, the desertification area has expanded at the rate of 2460 km² per year, causing a significant impact on human production and living activities and a global economic loss of more than USD 42.3 billion each year [1]. The main hazards of wind–sand flow on desert roads include sand burial and wind erosion [2]. Sand burial refers to the movement of wind and sand to the vicinity of the roadbed; the roadbed structure im-pedes the movement of wind and sand to the downwind side so that the wind and sand flow from an unsaturated state to a saturated state. This results in the accumulation of sand particles at different positions on the roadbed [3]. Many researchers have analyzed the flow field characteristics of windbreaks and the structural characteristics of wind and sand flow through on-site observation, wind tunnel experiments, and numerical simulation to obtain optimal windbreak setting parameters. Li [4] et al. investigated how windbreaks change the trajectory of wind–sand flow and found that the flow field of wind–sand flow over windbreaks varies with different heights of windbreaks, and the accumulation of sand behind the walls is also different. Avila-Sanchez [5] et al. combined wind tunnel experiments and particle image velocimetry to analyze the characteristics of the flow field in the vertical profile of a windbreak on a bridge by investigating the effect of the parameters of the windbreak on the reduction in the wind speed, such as the wind speed, the height of the enclosure, the installation spacing and the number of layers. Zhang [6] et al. determined that upright windbreaks provide good shielding by analyzing 3D flow lines and the surface pressure of trains. Sicot et al. [7] proved that the geometry of the wind barrier significantly affects the train’s aerodynamic force. Cheng Si Jin [8] et al. used wind tunnel tests to explore the critical length of wind barriers of different thicknesses and analyzed the effects of using wind barriers of various thicknesses on the aerodynamic characteristics of trains in wind tunnel tests.

According to the existing research, in exploring the process of wind–sand erosion, the existing methods are mainly used to numerically simulate the wind–sand flow by wind tunnel experimental test or Euler–Lagrange method, which analyses the parameters such as the shape, height, width, and porosity of the windbreak primarily. However, all of the above techniques need to divide the grid for analysis, and the small-scale wind–sand flow model is more sensitive to grid distortion. The grid division method seriously restricts the precision and accuracy of the wind–sand flow analysis. Therefore, this paper introduces the Smooth Particle Hydrodynamics (SPH) method to simulate the wind–sand flow numerically.

Smooth Particle Hydrodynamics (SPH) is a meshless Lagrangian particle method. In the SPH method, an object is described by a series of particles. These particles carry physical information, such as density and velocity [9,10]. The technique eliminates the sensitivity of small-scale models to large deformations or distortions of the mesh and allows the tracing of particle trajectories [11,12]. At the same time, this paper proposes a new optimization method for sand prevention based on the trajectory of particles, i.e., to explore the source of sand particles through the analysis of the concentration, velocity, and particle size distribution of sand particles to optimize the location of windbreaks and to effectively implement sand prevention measures. In this paper, a particle backtracking algorithm based on the SPH method is written in C. This algorithm can realize the capture and analysis of particles’ source, concentration, velocity, and particle size distribution and carry out sand blocking and sand treatment by investigating the movement mechanism of sand particles. In addition, considering the diversity of research needs and the fact that the particle sizes and proportions of sand particles are not fixed values, this paper combines the pixel recognition technique to improve a particle modeling algorithm that can realize arbitrary mixing of particle sizes and proportions, which further enhances the accuracy of wind–sand flow analysis.

2. Theory of the SPH Wind and Sand Flow Modeling Approach

2.1. Funembankmententals of the SPH Methods

The basic idea of the SPH method is (1) the computational domain of the desired solution is discretized into a series of arbitrarily distributed particles with independent properties that do not need to be connected by grids. (2) The kernel approximation method is adopted; the field function can be represented by integral and then the step-by-step integration of any field function can be completed. (3) In the calculation process, the superposition summation of the values of the active particles in the target particle support domain is used to replace the integral representation of the field function and its derivatives.

2.2. Basic Equations of SPH for Wind–Sand Flow

2.2.1. Kernel Function Interpolation

The construction of SPH equations is usually carried out in two main steps, the first is integral representation, and the second is particle approximation. The arbitrary function and the smooth function are integrated step by step, and then the summation approximation is performed by the values of the nearest neighboring particles.

In the SPH method, the core is the difference; for a continuous smooth function f(x), the value of the function at a point on the domain of definition is

$$\begin{array}{c}\begin{array}{c}f\left(x\right)=\underset{\Omega }{\int }f\left(x^{\prime }\right)\delta \left(x-x^{\prime }\right)d{x}^{\prime }\end{array}\end{array}$$

(1)

Replacing the Dirac function in (1) with a smooth function w(x − x′), the integral representation of the function can be written as

$$\begin{array}{c}⟨f\left(x\right)⟩=\underset{\Omega }{\int }f\left(x^{\prime }\right)W\left(x-x^{\prime },h\right)d{x}^{\prime }\end{array}$$

(2)

where h in Wx − x′, h is the smooth length used to define the range of influence of the function, and the formula is the kernel approximation formula.

For the spatial derivative ∇⋅f(x) of the function f(x), it follows that

$$\begin{array}{c}⟨\nabla \cdot f\left(x\right)⟩=\underset{\Omega }{\int }\left[\nabla \cdot f\left(x^{\prime }\right)\right]w\left(x-x^{\prime },h\right)d{x}^d\end{array}$$

(3)

2.2.2. Particle Approximation

In this paper, the particle approximation method is used to interpolate the kernel function discretization in order to obtain the discrete control equations.

If the volume $V_j$ of the particle is used to approximate the microelement dx′ at particle j in the integral, the mass of the particle can be expressed as

$$\begin{array}{c}{m}_j=∆V_{j{\rho }_j}\end{array}$$

(4)

where ${\rho }_j$ is the density of particle j

Particle discretization of Equation (2) yields

$$\begin{array}{c}⟨f\left(x\right)⟩={\displaystyle \sum _{\dot{j}=1}^N}{\displaystyle \frac{m_i}{{\rho }_j}}f\left(x_j\right)W_{ij}\end{array}$$

(5)

where $W_{ij}=W\left(r_i-r_j,h\right)$.

Therefore, the particle estimated of the spatial derivative of the field function at particle $i$ is

$$\begin{array}{c}⟨\nabla \cdot f\left(x_i\right)⟩={\displaystyle \frac{1}{{\rho }_i}}{\displaystyle \sum _{j=1}^N}{m}_j\left[f\left(x_j\right)-f\left(x_i\right)\right]\cdot {\nabla }_i{W}_{ij}\end{array}$$

(6)

2.2.3. Selection of the Kernel Function

Pressure instability problems are triggered during the movement of wind–sand flow due to the complex energy exchange between the two types of particles [13]. To cope with this problem, the five times spline kernel function is chosen as the kernel function in this paper. The advantage of the quintic spline kernel function over other kernel functions is that its second-order continuous smooth derivative is inversely related to the interparticle distance. The smaller the distance, the larger the derivative; the more significant the distance, the smaller the derivative. This characteristic makes it highly applicable in solving pressure instability problems. The five times spline kernel function expression is given below.

$$\begin{array}{c}W\left(R,h\right)={\alpha }_d\times \left\{\begin{array}{c}(3-R{)}^5-6(2-R{)}^5+15(1-R{)}^5,0\le R<1\\ {\displaystyle \genfrac{}{}{0pt}{}{(3-R{)}^5-6(2-R{)}^5,1\le R<2}{(3-R{)}^5,2\le R<3}}\\ 0,R>3\end{array}\right.\end{array}$$

(7)

The values of ${\alpha }_d$ in one, two and three dimensions are ${\displaystyle \frac{120}{h}},{\displaystyle \frac{7}{478\pi h^2}} \mathrm{a}\mathrm{n}\mathrm{d} {\displaystyle \frac{3}{359\pi h^3}}$.

2.3. Construction of SPH Control Equations for Wind–Sand Flow

Wind–sand flow is a complex, nonlinear, self-organizing two-phase flow [14].

In the SPH approach, we simplify the wind–sand flow into the following features:

(1)

In terms of fluid viscosity, the wind–sand flow is considered a Newtonian fluid.

(2)

In terms of flow characteristics, the wind–sand flow is regarded as incompressible flow.

(3)

In the near-surface layer flow, the gas flow that can drive the sand particles to jump and move is basically in a turbulent state [15].

Wind–sand flow is the constant movement of sand particles entrained by air currents. Since the velocity variation of the wind field in the longitudinal direction in the near-surface plane is very small, the near-surface layer is simplified to a two-dimensional space in this study [2]. The equations governing the flow in the near-surface layer are mainly described as conservation of mass and momentum, to which the SPH particle approximation is applied, and the discretization yields [16]:

Continuous density equation:

$$\begin{array}{c}{\displaystyle \frac{\partial {\rho }_i}{\partial t}}={\rho }_i{\displaystyle \sum _{j=1}^N}{\displaystyle \frac{m_j}{{\rho }_j}}{u}_{ij}^{\beta }\cdot {\displaystyle \frac{\partial w_{ij}}{\partial x_i^{\beta }}}\end{array}$$

(8)

where $m_j$ is the mass of particle j, ${\rho }_i$ is the density of particle j, β denotes the direction of the coordinate axis.

Momentum equation:

$$\begin{array}{c}{\displaystyle \frac{d{\overrightarrow{v}}_i^{\alpha }}{dt}}={\displaystyle \frac{1}{{\rho }_i}}{\displaystyle \sum _{j=1}^N}{\displaystyle \frac{m_i}{p_j}}{\sigma }^{\alpha \beta }{\displaystyle \frac{{\partial W}_{ij}}{\partial {\overrightarrow{x}}_i^{\beta }}}\end{array}$$

(9)

where α denotes the direction of the coordinate axis.

2.4. Artificial Viscosity

To solve the unphysical penetration problem generated when the particles are close together during the motion process, as well as to take into account the unavoidable dissipation problem, this paper introduces the Monaghan-type artificial viscosity ${\Pi }_{ij}$, which is expressed as follows:

$$\begin{array}{c}{\Pi }_{ij}=\left\{\begin{array}{c}{\displaystyle \frac{-{\alpha }_{\Pi }\overline{c_{ij}}{\varphi }_{ij}+{\beta }_{\Pi }{\varphi }_{ij}^2}{\overline{{\rho }_{ij}}}}\\ 0\end{array}\right.\begin{array}{cc},& v_{ij}\cdot x_{ij}<0\\ ,& v_{ij}\cdot x_{ij}\ge 0\end{array}\end{array}$$

(10)

where ${\varphi }_{ij}={\displaystyle \frac{h_{ij}{v}_{ij}\cdot x_{ij}}{{\left|x_{ij}\right|}^2+{\phi }^2}}$, $\overline{c_{ij}}={\displaystyle \frac{1}{2}}(c_i+c_j)$, $\overline{{\rho }_{ij}}={\displaystyle \frac{1}{2}}({\rho }_i+{\rho }_j)$, $h_{ij}={\displaystyle \frac{1}{2}}(h_i+h_j)$, $v_{ij}=v_i-v_j$, $x_{ij}=x_i-x_j$.

2.5. Coupled Treatment of Wind–Sand Two-Phase Flow

The wind–sand movement is a typical air–solid two-phase flow. The mutual coupling effect between air and sand particles is one of the factors that must be considered when dealing with wind–sand motion. Smooth length and smooth factor are commonly used in the SPH method to deal with the coupling between two phases of wind and sand flow. In this paper, the flow field composed of air and sand particles is discretized into a series of particles, and the particles are processed as a uniform diameter of the circle, so its support domain is also a circle, and the support domain radius formula is [17]

$$\begin{array}{c}{R}_i=k_i\cdot l_i\end{array}$$

(11)

where $R_i$ is the radius of the support domain of the $i$ particle, $l_i$ is the smooth length of the $i$ particle, which is often taken to be the diameter h of the $i$ particle, and $k_i$ is the smoothness factor of the $i$ particle

As shown in Figure 1, according to the different cases of coupling between the target particle and the acting particle, the values of the smoothing factor are categorized into the following three cases. (1) When the target particle and the acting particle are both gas-phase particles, take $k_1$ = 3.0, the interaction between the particles is more obvious, and the radius of the support domain $R_1$ = 3 h. (2) When the target particle and the acting particle are the same sand-phase particle, take $k_2$ = 1.1; they can interact only when they collide, and the support domain radius $R_2$ = 1.1 h. (3) When the target particle and the acting particle are two different particles, take $k_3$ = 2.5, and the radius of the support domain $R_3$ = 2.5 h.

2.6. Time Integral

This paper uses the leap-frog time integration format, which has the advantage of efficient computational power and low storage requirements. The particle’s density, velocity, internal energy, and position can be expressed using the formula [18].

$$\begin{array}{c}t=t+\Delta t\end{array}$$

(12)

$$\begin{array}{c}{\rho }_i\left(t+\mathsf{\Delta }t/2\right)={\rho }_i\left(t-\mathsf{\Delta }t/2\right)+\Delta t\cdot D{\rho }_i\left(t\right)\end{array}$$

(13)

$$\begin{array}{c}{v}_i\left(t+\mathsf{\Delta }t/2\right)=v_i\left(t-\mathsf{\Delta }t/2\right)+\Delta t\cdot D{v}_i\left(t\right)\end{array}$$

(14)

$$\begin{array}{c}{u}_i\left(t+\mathsf{\Delta }t/2\right)=u_i\left(t-\mathsf{\Delta }t/2\right)+\Delta t\cdot D{u}_i\left(t\right)\end{array}$$

(15)

$$\begin{array}{c}{r}_i\left(t+\mathsf{\Delta }t\right)=r_i\left(t\right)+\Delta t\cdot v_i\left(t+{\displaystyle \frac{\Delta t}{2}}\right)\end{array}$$

(16)

The leap-frog format is required to satisfy the CFL condition’s stability criterion. It guarantees that the time step is proportional to the smooth length. The time steps selected in this paper are as follows:

$$\begin{array}{c}△t=min\left(\xi h_i/\left[h_i△.v_j+c_j+1.2\left(\alpha c_j+\beta ma{x}_j\left(v_i-v_j\right)\right)\right]\right)\end{array}$$

(17)

3. Particle Modeling and Backtracking Algorithm

3.1. Particle Modeling

To provide a systematic method for generating sand particles of a single size and sand particles of mixed sizes. First, geometric modeling was performed and exported in PNG format through standard drawing software [19]. Then, the color filling is performed, and the pixels on the picture are identified using the pixel recognition algorithm to generate particles with the center of each pigment block as the center of a circle and pre-set different lengths as the diameters. The blue line is the demarcation line to generate air particles at the top and sand particles at the bottom. The final model can be transformed into the corresponding particles in the SPH numerical simulation calculations by the pre-processing program, and the corresponding particle position information can be generated.

In the sand particle generation process program, the particle generation algorithm prepared in this study can arbitrarily choose to generate sand particles of several particle sizes and the percentage of each particle. In Figure 2b,e show the generation of two-phase wind–sand flow with a single particle size, and c and f show the generation of two-phase wind–sand flow with mixed particle sizes. The radii of the sand grains in yellow are 0.050 mm, 0.100 mm, and 0.150 mm, with proportions of 50%, 30%, and 20%, respectively.

In this study, the mass of the smallest grain-size sand particles is set as the basic unit, and the mass of sand particles of different grain sizes will produce corresponding changes with changes in diameter. This will make the force analysis of sand grains more precise and more consistent with the trajectory of sand grains in nature.

In particle generation, to avoid significant gaps between the sand particles, this paper first carries out the generation of sand particles of larger particle size, sand particles, in the role of gravity for the natural accumulation. After reaching the initial stabilization stage, the smaller particles fill in the gaps. Eventually, under the interaction between the particles, a stable state of natural accumulation is achieved, which is also more in line with the existence of sand particles in the natural environment. At the end of the sand particle generation, air particles and dummy particles are generated, and the physics model is constructed using the dummy particles.

3.2. Particle Backtracking Algorithm

One of the advantages of the SPH method is the ability to observe the particle trajectory and capture the particle’s position. Therefore, this paper utilizes C language to write a particle backtracking algorithm to capture and analyze sand-buried particles’ source, concentration, velocity, and size distribution to take appropriate measures to cope with the wind–sand disaster. In addition, this paper will utilize the particle backtracking algorithm to analyze the classical embankment wind–sand flow field, to analyze the embankment wind–sand flow field with different wind speeds and mixed particle sizes, and then propose practical solutions for the management of embankment wind–sand disasters.

3.3. Parameterization

According to the wind–sand two-phase flow dynamics model, the air, sand particles, and embankment models were discretized into particles and then used for numerical simulation. By analyzing the simulation results, the particle movement trajectories are explored to trace the initial distribution positions of the deposited sand grains and their hopping motion, thereby enabling an investigation of the wind and sand flow field around the embankment. The geometric model uses the dimensions and proportions in Section 3.1, and the parameters are calculated in

Table 1. Parameter setting.

Gas-Phase Particles	Sand Particles	Calculation Parameters
Diameter of gas particle ds = 0.1 mm	Diameter of sand particle ds = 0.1 mm	Karman constant Κ = 0.4
Density of gas particle ρs = 1.293 kg·m−3	Density of sand particle ρg = 2650 kg·m−3	Friction wind velocity μ* = 0.13 m·s−1
Dynamic viscosity V = 1.8 × 10−5	Coefficient of friction μ = 0.4	Time step Δt = 1.0 × 10−5 s
Mass of gas particle m = 0.523 kg	Coefficient of restitution e = 0.85	Total time step N = 50,000

below.

3.4. Boundary Condition

Both the entrance and exit of the computational domain are set as periodic boundaries to realize infinite boundaries to facilitate the observation of particle trajectories throughout the simulation process to analyze the corresponding structural characteristics of the wind–sand flow under the influence of the embankment. Virtual particles are placed along both the upper and lower boundaries, and the repulsive forces between these virtual particles and adjacent real particles are utilized to construct a solid wall boundary, preventing gas-phase particles from penetrating sand-phase particles. As shown in Figure 2, the repulsion formula between virtual and real particles is defined as follows [20]:

$$\begin{array}{c}{f}_{\left(\mathit{r}\right)}=A\left[{\left(r_0/r\right)}^n-1\right]\mathit{r}/r\end{array}$$

(18)

where $r_0$ is the initial spacing of the particles, $r$ is the spacing between dummy and fluid particles, A is 600, and n is 4.

The velocity boundary is set at the entrance because the sand particles start to jump and move under the action of airflow. The airflow velocity follows the law of exponential distribution, so the wind speed contour line formula is [17]

$$u_z={\displaystyle \frac{u_{*}}{k}}\mathit{ln}\left({\displaystyle \frac{z}{z_0}}\right)$$

(19)

4. Validation and Analysis of Results

4.1. Model Validation and Sand Barrier Principle Analysis

The jump trajectories of sand particles with different grain sizes are shown in Figure 3. When the wind–sand flow reaches the stabilization stage, the trajectories of the leapfrog motion of the five grain sizes of sand particles show typical characteristics. These trajectories show that the sand particles rise more slowly and fall more rapidly, showing an asymmetric parabolic shape. These observations are highly consistent with the sand grain leapfrog motion recorded by Zhiqiang Li et al. [21] through wind tunnel experiments and high-speed photography. In addition, the jump height gradually decreases as the radius of the sand grains increases. This phenomenon can be attributed to the increase in particle mass due to the increase in particle size of the sand particles, which decreases the sand particles’ jump height under the same airflow conditions.

As shown in Figure 4, the concentration of sand particles of different types of mixed grain sizes was compared. The concentration of sand particles of different grain sizes in all three pictures showed an exponential decay trend along the height variation, which is consistent with the experimental results of Wang Hongtao et al. [22]. In addition, it can be observed from Figure 4a,b that the larger the diameter of the sand grains, the slower the rate of concentration decay. This is due to the fact that larger-sized sand particles acquire more kinetic energy upon collision, causing their concentration to decay more slowly along the altitude.

As shown in Figure 5, the mixed grain size sand transport rate varies along the height. The mixed grain size sand particles have a grain size of 0.05 mm, 0.10 mm, 0.15 mm, 0.20 mm, 0.25 mm, and the proportion of each is 45%, 25%, 15%, 10%, 5%. The figure reveals that both the mixture of particle sizes and the sediment transport rate for each particle size exhibit a trend of initially increasing and then decreasing. This pattern aligns with the “overshoot” phenomenon observed in the experiments by Shao and Raupach [23], where the distribution of sediment transport rate along the flow path exhibits characteristics of “increase-saturation-decline-stabilization”. The figure also indicates that the sediment transport rate of sand particles decreases with increasing particle size. This is primarily due to smaller particles being more easily mobilized under wind forces, requiring less energy to reach greater heights and achieve higher velocities, in contrast to larger particles, which exhibit the opposite behavior.

Figure 3, Figure 4 and Figure 5 verify the correctness of the numerical simulation calculation method of wind–sand two-phase flow based on particle modeling of the SPH method.

The concentration distribution of single-grain size sand particles is shown in Figure 6. At Step = 100, the sand is in the just-started phase, and the sand concentration hardly changes. At Step = 25,000, the wind–sand flow reaches a stabilization stage where the peak of sand concentration is concentrated before the windward slope of the embankment, but some sand particles still fall on the embankment surface. This is due to the height of the embankment altering the distribution characteristics of the flow field, leading to the formation of a wind speed reduction zone in front of the windward slope, which decreases the kinetic energy of the sand particles and effectively obstructs their movement. However, due to the limitations imposed by the embankment’s height, some sand particles still possess sufficient kinetic energy to land on its surface. We can find that increasing the embankment height can effectively prevent sand particles from falling on the embankment surface, so windbreaks can be installed to prevent particles of sand from falling on the embankment surface.

Figure 7 shows the source analysis of sand-buried particles at different time steps. In this study, a sand-buried particle backtracking algorithm was written to capture sand-buried particles at three different time steps by utilizing the SPH method to observe the entire motion process of the particles. At 15,000 steps, most of the sand particles have begun to move; at 25,000 steps, the sand flow has already started to enter the stabilization phase; and at 35,000 steps, the sand flow has been in the stabilization phase for a certain period. Capturing the sand-buried particles at three different time-steps, it can be roughly found that the sand-buried particles are concentrated in the range of 0.005 m–0.009 m, and the approximate location of the peak of the sand-buried particles source at 0.008 m to set up a windbreak, as shown in Figure 8.

As shown in Figure 9, Figure 10 and Figure 11, demonstrating the concentration of sand particles after the installation of the windbreak shows that there are almost no sand particles on the roadway embankment surface and the two embankment surfaces with the foot of the slopes at both locations. As shown in Figure 11, the overall decrease in particle concentration with a windbreak is 97%, with a 96% decrease on the windward slope and a 98% decrease on the leeward slope.

As shown in Figure 12, comparing the above two pictures, it can be found that the vortex location of the flow field appears near the windbreak after setting up the windbreak. This indicates that the velocity decay region of the flow field is shifted forward, and the sand particles are blocked by the windbreak wall in advance so that the sand particles do not have enough kinetic energy to run onto the embankment, thus effectively solving the sand burial problem of the embankment. This also demonstrates the validity of analyzing the initial positions of buried particles through a particle backtracking algorithm, thereby informing the implementation of effective sand control measures.

As shown in Figure 13, the sand velocity distribution is compared between the two working conditions. In the embankment wind–sand flow field, the airflow field near the embankment is disturbed due to the disturbance of the embankment. Wind speeds decrease dramatically at the foot of windward slopes, creating a wind speed attenuation zone that significantly prevents the transport of sand particles. However, some sand particles still cross the embankment, and at the leeward slope, the velocity decreases abruptly due to a sudden change in velocity in the turbulent airflow field. After setting up the windbreak, the flow field around the windbreak is disturbed, and the particle velocity changes abruptly, forming a wind speed attenuation zone, which prevents the transport of sand particles.

4.2. Kinetic Analysis of the Source Distribution of Sand-Embedded Particles in Road Embankments

4.2.1. Analysis of the Source of Sand-Embedded Particles in Road Embankments with Different Wind Speeds

As shown in Figure 14, the sand concentration on the embankment varies under different wind conditions. The higher the wind speed, the higher the concentration of sand particles on the embankment surface and at the foot of the leeward slope as the wind speed increases. This illustrates that the higher the wind speed, the more sand particles collide, the more kinetic energy the sand particles gain, the more likely they are to leap over the embankment, and the more severe the sand buildup.

As shown in Figure 15, the analysis of the sources of sand-buried particles shows that the sources of sand-buried particles are mainly concentrated in the range of 0.004 m–0.0009 m. However, as the moored wind speed increases, the source range remains almost unchanged, but the peak point shifts left and rises. The greater the wind speed, the greater the kinetic energy gained by the sand particles, and the more jumps occur in the sand particles, the shorter the time for the wind–sand flow to reach equilibrium and the peak of the initial distribution of sand-buried particles is shifted to the left. This is the exact reason for the different sand concentrations under different wind speed conditions, as in Figure 14, and verifies the correctness of the particle backtracking algorithm.

As shown in Figure 16, the windbreak embankment sand concentration distribution. Almost no sand particles were present on the embankment surface for the three wind speed conditions. Sand concentrations on both windward and leeward slopes were reduced to less than 5%. Sand concentrations in front of the windbreak gradually increased with increasing wind speed, while random fluctuations in sand concentrations were observed in the vicinity of the windbreak and behind it. After crossing the leeward slope, the effect of concentration on the variation of wind speed decreases and gradually converges.

As shown in Figure 17, a comparison of the sand velocity distribution before and after the installation of the windbreak. As shown in Figure 17a, due to the embankment’s disturbance of the flow field, the particle velocity at the windward slope of the embankment decreases substantially, preventing the transport of some sand particles. Then, the sand particles move to the leeward slope; due to the turbulence of the airflow field, the sand particles constantly collide and the kinetic energy of the particles increases. It is worth noting that the higher the wind speed, the more pronounced the turbulence of the airflow field and the greater the kinetic energy the sand particles gain. The windbreak-embankment sand velocity distribution is shown in Figure 17b, and the sand velocity around the windbreak and the embankment is significantly reduced. This is due to the obstruction of the retaining wall and barrier; the flow field around the retaining wall and embankment is disturbed to form a wind speed attenuation zone, and then effective sand blocking can be carried out. From the comparative analysis of (a) and (b), it can be seen that after setting up the windbreak, the windbreak interferes with the airflow field, and the flow field around the windbreak is disturbed, which reduces the kinetic energy of the sand particles.

4.2.2. Source Analysis of Mixed Grain Size Sand Buried Particles

Due to the microscopic scale, this paper aims to investigate the effect of different particle size compositions on the source distribution of sand-buried particles in a mixture of particle sizes. Therefore, the parameters for modeling sand grains using the particle modeling algorithm based on the SPH method in this study are shown in

Table 2. Mixed grain size sand parameters.

Computational Domain	Particle Diameter	Percentage of Particles
190 mm × 90 mm	R = 0.00005 m	100%
190 mm × 90 mm	R = 0.00005 m, R = 0.00010 m, R = 0.00015 m	50%, 30%, 20%
190 mm × 90 mm	R = 0.00005 m, R = 0.00010 m, R = 0.00015 m, R = 0.00020 m, R = 0.00025 m	45%, 25%, 15%, 10%, 5%

below:

As shown in Figure 18, the particle source distribution of sand-embedded particles from the embankment. We utilize a backtracking algorithm to track the source of sand-buried particles in mixed-grain wind–sand streams. It can be found that the peak point of the sand-buried particle source gradually shifts to the right as the number of mixed particle size species increases. This is because the more significant the diameter of the sand particles, the more kinetic energy is required for the particles to make the jump and the longer it takes for the wind–sand flow to reach equilibrium.

According to Figure 18, the peak of the source of sand-buried particles can be made out, and a windbreak wall is provided at the peak thereof. The windbreak embankment sand concentration distribution is shown in Figure 19. We can find that no sand particles are present on the embankment surface, and the concentration of sand particles at the foot of the leeward slope in the windward slope map is reduced to less than 10% below. The correctness of utilizing a particle backtracking algorithm to track the source of sand-buried particles and then installing a retaining wall for sand-blocking measures is illustrated. It is worth noting that the concentration of sand particles at the leeward slopes increases with the increase in the particle size species.

4.3. Kinetic Analysis of Sand Particle Size Distribution in Embankment Flow Field

The purpose of this study is to investigate the particle size distribution changes in the wind and flow field after setting up the windbreak at the microscopic scale, so the particle radius is set to be 0.05 mm, 0.10 mm, 0.15 mm with three kinds of particle sizes, and their proportions are 50%, 30%, and 20%, respectively.

As shown in Figure 20, sand particle size distribution around the windbreak. After installing the windbreak, the proportion of 0.05 mm sand decreased in front of the windbreak, and the proportions of 0.10 mm and 0.15 mm sand increased, but the proportions fluctuated less. However, the percentage of particles fluctuated significantly behind the windbreak, with a significant decrease in the rate of 0.05 mm sand particles and a substantial increase in the rate of 0.10 mm versus 0.15 mm sand particles. This is because the airflow in the windbreak around the formation of a vortex and movement to the windbreak before the speed began to change. The larger the size of the particles transported by the kinetic energy required to transport, the larger the vortex area, and the more pronounced the influence of the larger size of the sand particles is more likely to be here to form the accumulation of sand.

As shown in Figure 21, sand particle size distribution at the foot of the leeward slope. Before the installation of the windbreak, the airflow was disturbed by the embankment, and the sand particle size distribution was slightly altered at the foot of the leeward slope, but the overall change was not significant. After setting up the windbreak, the disturbance of the flow field becomes more extensive, the proportion of 0.05 mm particles decreases, and the proportion of 0.10 mm and 0.15 mm particles increases significantly. This indicates that the vortex region formed after the installation of the windbreak has a more significant impact on larger grain-size sand particles and that large-grain-size sand particles are more likely to accumulate at the foot of the slope.

5. Conclusions

In this paper, a particle backtracking algorithm was written to analyze the source, concentration, velocity, and size distribution of sand-buried particles by taking advantage of the fact that the SPH method can track the trajectory of the particles in order to take appropriate measures to cope with the wind–sand disaster. At the same time, a particle modeling algorithm that enables arbitrary mixing of particle sizes is also improved. In addition, this paper will analyze the classical embankment wind–sand flow field using a particle backtracking algorithm and draw the following conclusions:

(1)

The concentration of sand particles on the embankment increases as the moored wind speed increases. The concentration of sand particles at the leeward slopes increases as the number of sand particle size species increases.

(2)

As the drag wind speed increases, the range of sand-buried particle sources remains almost unchanged, but the peak point then shifts left and rises. The peak point of the sand-buried particle source gradually shifted to the right and increased with the increase of sand particle size types.

(3)

Vortices appear around the wind wall, sand particles are affected by the vortex, and deposition occurs; the larger the particle size, the more obvious the particles are affected by the vortex.