Analyzing the causes of sand damage on desert highways from the perspective of formation mechanisms (Figure 1), the range within the roadbed and pavement of desert highways can be viewed as a system, which inherently possesses a certain sand-passing capability (Figure 2). The formation of highway sand damage occurs when the external intensity of wind-blown sand exceeds the capacity of the roadbed and pavement to transport sand. Moreover, the sand transport capacity of the roadbed and pavement is further controlled by the cross-sectional form and geometric parameters (optimizable parameters), attribute parameters (selectable parameters), and assigned parameters (controllable parameters), as shown in Table 1. This paper takes the common roadbed structures found in the construction of the Tazhong-Desert Highway and the Tumxuk-Kunyu Desert Highway in Southern Xinjiang as the engineering background, generalized into five numerical models of roadbed structures. Utilizing optimization analysis theory methods, the paper analyzes the optimal sand transport geometric parameters of each typical roadbed type, providing engineering reference data for selective wind and sand control project construction along the route.

As shown in Table 1, during the parameterized modeling process of desert highway roadbeds, various parameters are classified, including the geometric parameters, attribute parameters, and assigned parameters of the roadbed model. Attribute parameters are selected based on practical considerations, while both geometric and assigned parameters are variable, serving as targets for modeling optimization and control. The wind-blown sand flow intensity is an assigned parameter in the modeling process. Control of wind-blown sand flow in the parameterized model involves reducing the sand phase percentage or increasing the gradient wind speed value (the sand phase refers to the sand component in two-phase flow with wind and sand). In practical engineering, this corresponds to variations in the sand-capturing effects of anti-sand engineering along the highway. When utilizing the parameterized model for optimization and analysis, with a constant wind-blown sand flow intensity, higher wind speeds result in more effective natural sand transport. In real-world projects, as long as the wind speed exceeds the threshold for sand initiation, the sand-carrying capacity remains in a non-saturated state, allowing it to smoothly pass through the highway pavement.

Using the Computational Fluid Dynamics numerical modeling method, the main focus is on analyzing the transport characteristics of sand particles near the roadbed. Therefore, the Eulerian two-fluid model is chosen for solving the model. In this model, the air phase and the sand phase are calculated separately as fluid and pseudo-fluid, respectively, and are considered as coexisting and mutually permeating continuous media. The gas-solid two-phase flow follows mass conservation equations and momentum conservation equations (Xin et al., 2021; Yu, 2022), which are as follows:

Mass conservation equations such as Eqs 1, 2: ∂ φ a ρ a u x ∂ x + ∂ φ a ρ a u y ∂ y + ∂ φ a ρ a u z ∂ z = 0 (1) ∂ φ s ρ s u x ∂ x + ∂ φ s ρ s u y ∂ y + ∂ φ s ρ s u z ∂ z = 0 (2)

Conservation of momentum equations such as Eqs 3, 4: ∂ ∂ t φ a ρ a U a + ∇ ∙ φ a ρ a U a U a = − φ a ∇ p + ∇ τ a + φ a ρ a g + f s a (3) ∂ ∂ t φ s ρ s U s + ∇ ∙ φ s ρ s U s U s = − φ s ∇ p − ∇ p s + ∇ τ s + φ s ρ s g + f s a (4)

The equation for the interaction force between gas phase and sand phase, as represented by Eq. 5: f s a = 3 C D φ s φ a ρ a 4 d U r U r φ g − 2.65 (5)

In the formula, φ _a and φ _s represent the volume fractions of the air phase and sand phase, respectively, satisfying the relation φ _a + φ _s = 1; t denotes the calculation time; ρ _a and ρ _s are the densities of the air phase and sand phase, respectively; u _x , u _y , and u _z represent the velocity components in the x, y, and z directions, respectively; τ _a and τ _s are the viscous stress tensors of the air phase and sand phase, respectively; ∇ represents the Hamiltonian operator; ∇p stands for the momentum change; P is the common pressure shared by both phases; P _s is the pressure of the sand phase; g is the acceleration due to gravity; U _a and U _s denote the velocities of the gas phase and sand phase, respectively, and U _r represents the relative velocity between the two phases; C _D denotes the drag coefficient.

The geometric structural types of desert highway subgrades mainly include five categories: non-fill cut subgrade, embankment subgrade, cutting subgrade, and semi-fill cut subgrade (uphill flow, downhill flow). However, from the perspective of parametric modeling, there are only four types, as shown in Figure 3. During parametric modeling, geometric parameters such as cutting depth, embankment height, and subgrade side slope gradient are considered as optimizable parameters; attribute parameters such as road surface roughness and subgrade slope roughness (roughness under the effect of engineering measures like grass squares on cutting and embankment slopes) are considered as selectable parameters; assigned parameters such as wind-blown sand flow intensity and incoming wind speed are considered as controllable parameters. Among these, attribute and assigned parameters remain constant during optimization, while geometric parameters are variable.

Taking the embankment subgrade as an example, a three-dimensional model of the subgrade is established using the Design Modeler in ANSYS Workbench, as shown in Figure 4. With the desert highways in the hinterland of the Taklamakan Desert as the background, the subgrade is designed according to the standards of a secondary highway, with a width of 12 m, an embankment height of H, a slope angle α, and the angle between the incoming flow direction and the embankment orientation as A. To allow the wind-blown sand flow to fully develop, the distance from the inlet to the embankment is set to 40 m, and the distance from the outlet to the embankment is 80 m (Xu et al., 2018; Mehdipour and Baniamerian, 2019). To eliminate the influence of the computational domain size on the structure, the computational domain is set to have a length of (120+2H/tanα) m, a height of 30 m, and a width of 20 m. Mesh division is conducted using ICEM, with a global mesh size of 0.1 m. To improve computational accuracy and more accurately reflect the flow field changes in the near-wall region, a boundary layer is established in the near-wall area to ensure that y+ is within the logarithmic law region, which requires 11.5–30 < Y+ < 200–300 (Choi et al., 2018; Horvat et al., 2021). Therefore, the first layer mesh size is set to 0.02 m, with 10 layers in the boundary layer, resulting in 4.95×106 meshes, a mesh quality of 0.87, skewness of 0.15, indicating good mesh quality and satisfying the requirements for independence.

Wind-blown sand flow consists of air and sand particles, requiring the setting of model parameters such as material density, particle size, and volume fraction during numerical calculations. Studies have found that the particle size of sand grains in wind-blown sand flow generally ranges from 0.075 to 0.25 mm (Cheng et al., 2017; Ding et al., 2021). For the sand particle properties of the study area, based on the particle size analysis results, simulations use a single particle size of 0.15 mm. The sand particle density is set to 2600 kg/m³ (Chen et al., 2019). Regarding the selection of sand particle volume fraction, considering the research area is located in the hinterland of the desert and combining with field meteorological data, the sand particle volume fraction is set to 0.05, uniformly distributed at the inlet. Since the Mach number of the wind-blown sand flow is less than 0.3, the density variation of the wind-blown sand flow is not considered. The pressure is set at atmospheric pressure; the inlet is a typical wind speed profile flow. The equation for solving the wind speed profile, as illustrated by Eq. 6: v y = v k ln y y 0 (6)

In the formula: v represents the friction velocity; y ₀ is the roughness length; k is the von Karman constant, taken as 0.4; y is the height; v(y) denotes the wind speed value at height y.

The solution method is set to reflect the unsteady, pressure-based two-phase flow of wind and sand. The turbulence model selected is the standard k-ε model. To ensure computational accuracy and accelerate the convergence of the iterative process, standard wall functions are used for the wall region. The equation solving method employs the SIMPLEC algorithm, with convergence criteria set to 1×10⁻⁶.

To consider the adverse effects of wind speed on desert highways, the maximum values of the maximum and average wind speeds from various azimuth angles of the desert highway are taken as the computational wind speeds. As shown in Tables 2, 3, based on the onsite measured wind speed data from 2018 to 2021, the maximum values of the maximum and average wind speeds for the desert highways in the hinterland of the Taklamakan Desert are 14.8 m/s and 8.71 m/s, respectively. Therefore, in calculations, wind speeds are set to 15 m/s and 9 m/s. The sections of the desert highway in the hinterland of the Taklamakan Desert that suffer severe sand damage have their orientations intersecting at large angles with the prevailing wind direction. To fully consider the impact of the geometric structure of the sand-transporting subgrade on the sand-transporting capacity of the desert highway itself, it is advisable to avoid small-angle intersections between the incoming flow direction and the embankment orientation. Hence, the angles between the embankment orientation and the incoming flow direction are set to 60° and 90°.

This paper proposes an optimization method for the geometric parameters of desert highway embankments based on DOE (Design of Experiments) and the Kriging model. Initially, a DOE is conducted to obtain initial sample points; subsequently, a surrogate model is generated using the Kriging method, with the optimization process illustrated in Figure 5.

The optimization process is completed within the Design Exploration module provided by Ansys Workbench. DOE experimental design methods include Central Composite Design (CCD), Box-Behnken Design, Latin Hypercube Sampling Design (LHS), and Optimal Space-Filling Design (OSF), among others. In this study, when generating sample points for input parameters using DOE experiments, the Optimal Space-Filling sampling method is adopted with the sample type chosen as CCD. This method enables uniform distribution of design parameters throughout the entire design space, allowing for maximum insight into the design with the fewest possible number of sample points (Xin et al., 2018). A total of 15 samples are required, and the optimization experimental sample combination is illustrated in Figure 6.

Due to its significant efficiency improvement in engineering optimization design problems, surrogate model methods are widely applied in engineering design. Commonly used surrogate models include the Standard Response Surface and Kriging (Zhang and Qiu, 2021; Zhu et al., 2022). Upon comparison, it has been found that the Kriging model demonstrates a better fit between the observed values at design points and the predicted values by the response surface. The criterion for fit assessment is: the closer to the line y=x, the better the fit, as shown in Figure 7. As an unbiased estimator model that minimizes estimation variance, Kriging exhibits ideal fitting effects for problems with high nonlinearity and maintains high accuracy both globally and locally. Therefore, in optimizing the geometric parameters of the subgrade, the Kriging model is chosen for optimization in this study.

The Kriging optimization model is a regression algorithm for spatial modeling and prediction of stochastic fields based on the theory of correlation functions (Wan et al., 2022; Liu et al., 2023). It is an optimal unbiased estimation technique, composed of a global trend function F(β, x) and a stochastic deviation Z(x). Its mathematical expression is given by Eq. 10: y x = F β , x + Z x = f T x β + Z x (10)

In the formula, F(β, x) = f ^T (x) β, where β are the coefficients of the regression function f ^T (x), and Z(x) indicates that local errors are approximately key to the model’s accuracy. The characteristics of Z(x) are shown in Eqs 11–13: E z x = 0 (11) V a r z x = σ 2 (12) C ov Z x i , Z x j = σ 2 R c , x i , x j (13)

The mathematical expectation of the random deviation is zero, and its variance σ ² is minimized, which can be calculated by Eq. 4. Here, C_ov[Z(x _i ), Z(x _j )] represents the covariance of the random deviations, and R(c, x _i , x _j ) is the correlation function, indicating the spatial correlation between any two points. The expression for R(c, x _i , x _j ) is shown in Eq. 14: R c , x i , x j = ∏ k = 1 n R k θ k , x i k − x j k = ∏ k = 1 n R k θ k , d k (14)

Where n is a variable, θ _k is a scale factor, and d _k represents the Euclidean distance between two points. It can be seen from the above definition that the key to the accuracy of the Kriging Surrogate Model lies in the selection of the correlation function. Commonly used correlation functions include the exponential function, Gaussian function, linear function, and cubic spline function.

Figure 8 presents three-dimensional response surfaces of the output parameter S with respect to the input parameters h and α for different subgrade cross-sections. From Figures 8A, B, it can be observed that the conductivity of both embankment and cut section subgrades decreases with an increase in subgrade height/cut depth; however, the change in conductivity with slope gradient varies at different subgrade heights, demonstrating specific trends: at a subgrade height/cut depth of 0.5 m, the sand transport rate increases with an increase in the slope gradient of the embankment; at a height/depth of 4 m, the conductivity decreases with an increase in the slope gradient. The primary reason for these results is that, at a height/depth of 0.5 m, sand particles can pass smoothly over the embankment, but the sand transport path is longer at lower gradients than at higher ones, leading to lower conductivity at smaller gradients; when the height/depth reaches 4 m, sand particles cannot pass smoothly over the embankment, and only when the gradient decreases can the sand pass, hence the transport rate increases as the gradient decreases. For the semi-filled and excavated uphill subgrade cross-section, conductivity decreases with an increase in cut-and-fill height and with an increase in slope gradient, as shown in Figure 8C. The main reason for these results is that the semi-filled and excavated uphill subgrade cross-section type differs from the embankment and cut section subgrades, as wind-blown sand experiences two decelerations when passing over such subgrades, making it difficult for sand particles to pass over the subgrade even at lower cut-and-fill heights, and only decreasing the slope gradient facilitates sand passage. The conductivity of the semi-filled and excavated downhill subgrade decreases with an increase in cut-and-fill height and increases with an increase in gradient, as depicted in Figure 8D.

From Figures 9A, B, it can be seen that when the incoming wind speed is 15 m/s and the embankment orientation is 60°, the embankment height being less than 1.5 m results in a transport rate greater than 1. This indicates that, for embankment heights below 1.5 m, sand particles can successfully pass over the embankment. Meanwhile, at an incoming wind speed of 15 m/s and an embankment orientation of 90°, the embankment height must be below 1.25 m to ensure the smooth passage of sand particles, as shown in Figure 9A. The higher the incoming wind speed, the greater the embankment transport rate; under the condition of an incoming wind speed of 15 m/s, an embankment height of less than 1.25 m results in a transport rate greater than 1; when the incoming wind speed is 9 m/s, an embankment height of less than 0.85 m results in a transport rate greater than 1, as illustrated in Figure 9B.

Figures 9C, D illustrate the variation in the conductivity of cut section subgrades with cut depth under different cut orientations and incoming wind speeds. As can be observed from Figure 9C, the transport rate increases with a decrease in the depth of the cut across different orientations, and it is slightly higher when the cut orientation is at 60° compared to 90°. From Figure 9D, it is known that under an incoming wind speed of 15 m/s and a cut orientation of 60°, a cut depth of less than 1.2 m results in a transport rate greater than 1. Similarly, with an incoming wind speed of 15 m/s and a cut orientation of 90°, a cut depth of less than 1.0 m yields a transport rate greater than 1.

From the analysis above, it is known that for an embankment orientation of 60° and an incoming wind speed of 15 m/s, the upper limit for the sand transport embankment height is 1.5 m; for an embankment orientation of 90° and an incoming wind speed of 15 m/s, the upper limit is 1.25 m; for an embankment orientation of 90° and an incoming wind speed of 9 m/s, the upper limit is 0.85 m. For a cut orientation of 60° and an incoming wind speed of 15 m/s, the upper limit for the sand transport cut depth is 1.2 m; for a cut orientation of 90° and an incoming wind speed of 15 m/s, the upper limit is 1.0 m; for a cut orientation of 90° and an incoming wind speed of 9 m/s, the upper limit is 0.62 m.

To determine the specific sand transport cut section parameters under different conditions, the embankment heights are set to 1.5, 1.2, and 0.8 m, respectively, and the cut depths to 1.2, 1.0, and 0.6 m. The variations in the slope gradient of the subgrade at different sand transport subgrade heights and their impact on conductivity are then studied, thereby determining the optimal slope gradient for sand transport subgrades.

From Figure 10, it can be observed that the optimal sand transport slope for an embankment height of 1.5 m is 21.7°; for an embankment height of 1.2 m, the optimal slope is 22.4°; and for an embankment height of 0.8 m, the corresponding slope for sand transport is 24.5°. Thus, the parameters of the sand transport cross-section for embankment-style road foundations under different conditions can be determined, as shown in Table 4. For a cut depth of 1.0 m, the optimal sand transport slope is 21.7°; for a cut depth of 0.5 m, the optimal slope is 26.5°. Therefore, the parameters of the sand transport cross-section for cut-style road foundations under different conditions can be derived, as indicated in Table 5.