Optimal design of sand blown wind tunnel

This work investigates the airflow driven by dual axial-flow fans in an atmospheric boundary layer (ABL) wind tunnel and the expected entrainment of sand movement together. The present study is conducted via 3D numerical simulation based on modelling the entire wind tunnel, including the power fan sections. Three configurations of dual fans in the tunnel are proposed. Simulation results show that the airflow in the tunnel with dual-fan configuration can satisfy the logarithmic distribution law for ABL flows. The airflow driven by the dual fans placed together at the tunnel outlet is highly similar to that in the tunnel with single fans. Although the boundary layer thickness is reduced, the maximum airflow velocity (53.393 m/s) and turbulence intensity (12.02%), which are respectively 1.75 and 1.49 times higher than those under the single-fan configuration, can be reached when dual fans are separately placed at the tunnel inlet and outlet. The simulation and experiment manifest that the separated arrangement of dual fans in the tunnel should be suitable for the experimental study of aeolian sand transport. Some measures, such as wind tunnel construction adjustment and optimal roughness element arrangement, are necessary to guarantee the required boundary layer thickness in the wind tunnel.


Introduction
The atmospheric boundary layer (ABL) is the lowest part of the atmosphere. In this layer, airflow always displays turbulent characteristics, which can be attributed to the friction exerted by the wind against the ground surface and thermal processes [1,2]. When the ABL is neutral, that is, the thermal processes are absent, a logarithmic velocity profile u(z) can be characterized by the friction velocity u * and the roughness height z 0 [3]. The airflow in the neutral ABL can be simulated in a wind tunnel, which was initially applied to the pollutant diffusion study [4,5]. ABL flow simulation in the wind tunnel is realized using the following two methods: active [6][7][8][9][10][11][12] and passive [13][14][15][16][17][18][19] simulations. Passive simulation is more popular than the active one because the spires, barriers, and roughness elements of the former are used to facilitate airflow development into an ABL flow in the wind tunnel. ABL wind tunnel experiments are currently conducted in the study of various fields, such as environment, building, transportation, automobile, ecology, agriculture, and chemical engineering [20][21][22][23][24][25][26][27][28].
In the 1940s, Bagnold [29] and Chepil [30] processed wind tunnel experiments to investigate aeolian sediment and soil erosion. Since then, the ABL wind tunnel experiment is regarded as an important methodology in the research of environmental fluid dynamics, such as aeolian sand transport [31][32][33] and wind-blown snow [20,21,34,35]. The sand used in the wind tunnel experiment is also natural sand sampled from the field desert. This notion indicates that full-scale natural wind field and sand bed conditions cannot be simulated in a size-limited wind tunnel. Nevertheless, 80% of sand is transported within a height of 30 cm on the bed surface. Accordingly, approximate reproduction (not simulation) of field wind-sand flow near the ground surface was achieved in the wind tunnel. Thus, the airflow velocity in the wind tunnel should be sufficiently high to entrain the sand particles, that is, the flow rate must be large. The main flow velocity should be at least larger than 25 m/s to cover the entire condition of natural aeolian sand transport completely [36]. Airflow in the wind tunnel must demonstrate high turbulence intensity due to the existence of strong turbulence near the ground surface in the actual aeolian sand transport process; this phenomenon is in contrast with the requirement in conventional wind tunnels [37,38]. A large fan would also be an optimal choice to meet the two above-mentioned requirements. However, the natural wind is gusty, which is an actual problem that results in the deviation between wind tunnel experiments and field observations. This phenomenon simultaneously presents large-scale unsteady features and small-scale high-frequency pulses. The sand transport is considerably influenced by the unsteady turbulent characteristics of field gust, such as transport rate fluctuations and intermittent transports [39]. However, stimulating unsteady field gust in the wind tunnel is difficult by using only one large fan, which must possess considerable inertia. The theory of inverse Fourier transform indicates that one unsteady continuous wind speed series could be regarded as a superposition of numerous different sine curves of wind speed. The energetic frequency range of turbulence in ABL is mostly from 0.01 Hz to 1 Hz [40,41]. The field observation result at the south edge of the Taklimakan desert shows that 85% of the field gust energy corresponds to the frequency below 0.25 Hz [42]. Therefore, the low-frequency energetic part of the field gust could be approximately reproduced in the wind tunnel when the power part of the wind tunnel is constituted by multiple axial fans. Each fan would drive the airflow changes with sinusoidal variation by individually controlling the frequency converter. The high-frequency random pulse part of airflow could then be self-developed while the airflow passes over the spires and roughness elements.
The designed environmental wind tunnel must meet the following three requirements to realize the simulations of continuous migration of aeolian sand and sandstorm process: (1) the airflow velocity in the wind tunnel should be sufficiently strong; (2) the power fan inertia should be as small as possible to facilitate variable speed adjustment and ensure high flow rate; and (3) the airflow field in the wind tunnel should have high turbulence intensity to reproduce the near-surface strong turbulence atmosphere. The numerical simulation of the entire structure and corresponding flow fields of the wind tunnel is an undeniably advanced and effective technique for the optimization of wind tunnel design [43]. This work focuses on 3D computational fluid dynamics (CFD) of the wind tunnel technology. The optimal power system configuration is determined through comparative simulation and analysis of the boundary layer flow characteristics in the wind tunnel, such as airflow velocity, turbulence intensity, and boundary layer thickness, under different power system configurations (single/dual fans). The present study aims to provide simulation data references for the wind tunnel upgrading and related experiments in wind blowing sand transport.

Governing equations
The wind flow of the neutral ABL formed in the wind tunnel is incompressible. The Reynolds-averaged Navier-Stokes (RANS) approach was adopted because the study focuses on the generation of the timeaveraged steady wind field in the controlled wind tunnel. The continuity and momentum equations for pure wind flow can be respectively expressed as follows: where x i is the coordinate in the ith axis in the Cartesian coordinate system; ρ and μ are the density and dynamic viscosity coefficient of air, respectively; u i is the ith velocity component; u i is the fluctuation part of u i ; p is pressure; over bar represents the time-averaged value; and −ρu i u j is called the Reynolds stress. This stress is characterized by the Boussinesq eddy viscosity assumption, which is used to close equations.

Turbulence model
The simulation of ABL flows is usually conducted using the commercial CFD codes with steady RANS turbulence modelling, particularly the standard kmodel and standard sand-grain rough wall functions [44][45][46]. Several studies focused on certain physical flow phenomena in ABL; hence, the airflow derivation was neglected, and inlet profiles were set as the inlet condition. The power part was directly omitted from the computational domain even though the simulation would be used for comparison with the wind tunnel experiment. Such an approach generally resulted in an unexpected decline in the velocity and turbulent profiles specified at the domain inlet before reaching the observation section within the computational domain. This behaviour is a direct consequence of the inconsistency between the fully developed ABL inlet profiles and the wall function formulation. References [47][48][49][50] proposed formulations for the turbulence model constants C μ and σ ε to ensure stream-wise homogeneity when using the k profile proposed by Yang et al. [51]. Parente et al. [52] developed an improved k-ε model for the neutral ABL flow simulation with arbitrary sets of fully developed inlet conditions. The present study demonstrates the simulation of an actual ABL wind tunnel, including the rotating fans. Thus, the inconsistency between the fully developed ABL inlet profiles and the wall function formulation can be neglected. The renormalization group (RNG) k-ε model was used in this work based on various two-equation turbulence models utilized in the ABL numerical simulation [53]. The model developed by Yakhot and Orszag [54] is a two-equation turbulence model that solves the formulas for turbulent kinetic energy k and turbulence dissipation rate . This model can enhance the swirling flow accuracy by accounting for the swirl effect on turbulence. In the steady airflow, the two groups of equations for the turbulent kinetic energy k and the turbulence dissipation rate are as follows [54]: where μ eff = μ + μ t and μ t = ρC μ k 2 /ε are the effective and turbulence viscosities, respectively.

Numerical setup
The simulation was conducted by using Fluent, a commercial CFD code, to solve the governing equations. The seven coefficients of the RNG k-ε model were set as their standard values. The near-wall treatment adopted standard wall functions. The semiimplicit method for pressure-linked equation consistent (SIMPLEC) algorithm [55], which also belongs to the family of the well-known semi-implicit method for pressure-linked equation (SIMPLE) algorithm [56], was used for the pressure-velocity coupling based on the finite-volume method in the present simulations.
The SIMPLEC algorithm was more efficient in iterative computations for fine grid systems compared with the SIMPLE algorithm. Particularly, the multiple reference frame (MRF) model was employed to deal with the interaction between the rotating and the stationary regions [57].

Computational domains, grids, and boundary conditions
The ABL tunnel was designed as an open-type uniflow wind tunnel with an experimental section of 9000 mm L × 940 mm W × 940 mm H . Figure 1 shows the structural diagram of the entire wind tunnel. When both axial fans were placed together in part B, L A = 500 mm and L B = 2800 mm. When the dual fans were individually placed in parts A and B, L A = L B = 1800 mm. The two fans had the same arrangement in parts A and B. The flow driven by a single fan in the wind tunnel was also simulated to validate the computational model (Section 3). When the single fan was placed at part B, L A = 500 mm and L B = 1800 mm. Figures 2 and 3 show the arrangements of single/dual fans in the power parts and their geometric models. The fans have five uniform crosssectional blades with 370 mm height and 0 mm blade tip-clearance. Each fan was equipped with a circular cowling. The 3D geometric modelling and meshing were implemented by Gambit. In those studies, the X-, Y-, and Z-axes in the Cartesian coordinate system respectively represented axial (stream-wise), vertical, and span-wise directions of the wind tunnel. The origin of the coordinate system is located at the centre of the inlet cross-section of the experimental section. Symbols u, v, and w correspondingly represented the axial, vertical, and span-wise mean velocity components of the airflow. The TGrid meshing scheme was employed in the region enclosing the axial fan to generate the grids comprising tetrahedral/hybrid elements. The grids with hexahedral/wedge elements were created by the Cooper scheme in other regions of the wind tunnel model. Figure 4 shows the grids on the vertical section (X-Y plane) of the partial region, including the single fan. The grids of the model with dual fans located in one power part have similar meshing.    The rotational speed of each fan in the simulation was set as 3000 rpm. Specifically, the rotational speed of the rotating reference frame should be −3000 rpm. The rotating subdomain in MRF is also marked by a rectangular wireframe in Figure 4. The inlet and outlet conditions of the wind tunnel adopted the zero gauge pressure condition. The roughness height and constant of the wall surface were also respectively set to 300 mm and 0.7 to match the field conditions accurately.

Model validation
The wind tunnel model with only a single fan was simulated to validate the numerical model because many flow research reports in the ABL wind tunnel driven by a single fan are available. The fan was placed in part B to facilitate air entry into the wind tunnel [37].
A set of pre-simulations was conducted to check grid independence. The monitoring variable is used to select the airflow mass rate at the tunnel inlet. Figure 5 shows that the relative difference in mass rate is less than 4% and that for the last three grid meshing conditions is less than1.5%. Such a relative difference essentially satisfied the evaluation criteria setting of the engineering problems in CFD. Considering the residual errors in the simulations, the simulated result based on the third grid meshing scheme is selected in this study for analysis.
The grid convergence index (GCI) method proposed by Roache [58] was used in this study to evaluate the numerical uncertainty on the grid. This method is based on the generalized theory of Richardson extrapolation and involves the comparison of discrete solutions at two different grids of spacing (h) [59,60]. The GCI formula can be expressed as follows: where ε-relative error between solutions based on coarse and fine grids; p-formal order of accuracy of the algorithm; r-refinement factor between coarse and fine grids; F s -safety factor. Table 1 lists the GCI results. The GCI value decreased with mesh refinement. Considering the GCI evaluation criteria, Karimi et al. [61] indicated that the numerical uncertainty caused by the grid number has slight effects on the deviation between numerical and experimental results when the GCI is less than 4.5%. Thus, the solution accuracy in the cases presented in     Table 1 could be regarded as optimal when the grid number increased to 1,745,621. Figure 6 shows the contours of the time-averaged axial velocity u of the airflow on the vertical plane at the centre of the wind tunnel (X-Y plane, Z = 0). As previously mentioned, no decline can be observed in the velocity and turbulent profiles. The mean velocity u profiles of the airflow in the experimental section are uniformly distributed and demonstrate vertical gradients. Figure 7 presents the velocity profiles of airflow at three different locations in the experimental section of the tunnel. Locations A1-A3 denote the three locations whose distances away from the experimental section inlet are 1/6, 1/2, and 5/6 of the experimental section length. Among the three locations, A2 is the middle location of the experimental section, while A1 and A3 are symmetric around A2. The three locations are representative, and the typical flow field characteristics of Figure 9. Profiles of vs. z at locations A1-A3. the upstream, middle, and downstream of the experimental section can be provided. Symbol z is the height above the wind tunnel bottom. Only the flow information in the lower half part of the wind tunnel according to the airflow axial symmetry is shown in the figure. Figure 7 indicates that all profiles exhibit the boundary layer flow feature near the wind tunnel bottom. The profiles show that boundary layer thickness and mainstream velocity increased along the stream-wise airflow. This increase indicates a definite progressive development of the boundary layer flow in the wind tunnel. The boundary layer flow could be regarded as well-developed from location A2 ( Figure 6). The profiles of the turbulence kinetic energy k with height z at locations A1-A3 are highly similar (Figure 8). The difference among the k profiles could be attributed to the flow disruption effect of the air-breathing fan. The effect is remarkable when the distance away from the fan is close. Figure 9 exhibits the log-log plots of turbulence dissipation versus height z. The profiles under the boundary layer, such as those of velocity u, are overlapped together. The results in Figures 7-9 imply that the flow features of the boundary layer flow in the wind tunnel are well maintained along the stream-wise airflow.
The results in Figure 7 show that the boundary layer thickness δ at location A2 is equal to 184 mm. Figure 10(a) shows the velocity u profile of the airflow under the boundary layer at location A2 and its fitting curve based on the modified formula of the logarithmic distribution law for neutral ABL flow proposed by Richards and Hoxey [45].
where u * denotes frictional velocity, z 0 represents the aerodynamic roughness of the underlying surface, and Karman constant κ = 0.4. The results from the fitting curve equation (Figure 10(a)) indicate that frictional velocity u * = 2.782 m/s and the roughness z 0 = 2.39 mm; these values are consistent with the actual condition of the underlying surface covered vegetation over 50-60 mm in height [2]. The k and profiles are also respectively shown in Figure 10(b) and (c); these profiles present similarities to the results proposed by Juretic and Kozmar [53]. Figure 11 shows the profiles of the relative velocity u/U and turbulence intensity T u at location A2.   Herein, U denotes the velocity at the top of the boundary layer, which approximately equates to mainstream velocity. T u , a dimensionless relative variable, is used to represent the turbulent characteristic. The expression was defined herein as T u = 2/3k U. The maximum turbulence intensity T u,max equates to 8.09%, which appears at height z = 9 mm near the bottom surface of the wind tunnel. The profiles are highly similar to the wind tunnel experimental results processed by Zhang et al. [62]. Apart from the different dimensions between actual and simulated wind tunnels, the singledirectional velocity measurement and sand bed surface topography in the actual wind tunnel experiment could also be regarded as the factors that caused the specific difference between the simulated and experimental results. The analysis and comparisons indicated that the numerical model used in this study is suitable for simulating airflow in the entire ABL wind tunnel, including the blade-rotating fan.

Results and discussion
The simulation of the boundary layer flow in the wind tunnel with single fans is coded in this study as F1. Meanwhile, F2c and F2s respectively denote the simulations of the airflow driven by dual fans located at only Part B or at Parts A and B (Section 2.4). Similar to case F1, a set of pre-simulations was conducted to validate grid independence ( Figure 12 and Table 2). A shown in Figure 12 and Table 2, the error values of the mass flow rate at the third grid system and the corresponding GCI for the two different simulation cases essentially reached the set criteria during grid refinement with the increase in grid number.
Airflow in case F2c is highly similar to that in case F1. The profiles of non-dimensional velocity variable u/U of the airflow under the boundary layer at location A2 in cases F1 and F2c are drawn in Figure 13   case F2c. Herein, δ = 196 mm, u * = 4.326 m/s, z 0 = 3.65 mm, and U = 42.851 m/s. The turbulent intensity T u also increases in case F2. However, the profiles of T u in cases F1 and F2c are only slightly translated (Figure 13(b)). The Pearson correlation coefficient of both profiles of T u equates to 0.9995. Such a value also presents high similarity of the boundary layer flows in both cases.
The airflow in case F2s is different from that in cases F1 and F2c because the blade rotation of the airblowing axial fan at the wind tunnel inlet should cause the negative pressure zone behind the central region of the fan. Figure 14 shows the contours of the timeaveraged axial velocity u of the airflow in case F2s (on the X-Y plane, Z = 0). The airflow velocity u in the middle part of the tunnel decreases with heights similar to those in the boundary layer region. Although the influence of the air-blowing axial fan decreases with the flow distance, this reduction still affects the entire experimental section. Accordingly, the boundary layer thickness always increases along the airflow direction, Figure 15. Profiles of u/U in boundary layer at locations B1-B5. and the velocity U decreases along the airflow direction based on the continuity principle. Five locations were selected and labelled as B1-B5 to analyze the characteristics of the developing boundary layer flow in case F2s comprehensively. Herein, B1-B5 denote the locations whose distances away from the inlet of the experimental section are 1/10, 3/10, 1/2, 3/5, and 4/5 of the experimental section length. Locations A2 and B3 are identical (A2↔B3), that is, B3 is the middle location of the experimental section. However, locations B1, B2, B4, and B5 are asymmetric around B3. Locations B4 and B5 are purposefully moved upstream after the addition of a fan i at the inlet side of the wind tunnel to analyze the blowing flow pattern from the upstream accurately. The flow field difference under the action of single-and double-side fans is also distinguished. Figure 15 shows the profiles of the non-dimensional velocity variable u/U of the airflow in the boundary layer at the five locations (B1-B5). The airflow is not yet fully developed until the flow run passes location B3. Figure 16 and Table 3 represent the boundary layer flow characteristics of e airflow at location A2 (i.e. B3) among the three cases. Considering the effect of air-blowing axial fans, the velocity and turbulence intensity profiles of airflow in case F2s are different from those in the two other cases. Although the boundary layer thickness in case F2s is low, the values of other flow characteristic parameters, such as u * , U, and T u,max , are considerably larger than those in cases F1 and F2c.
Compared with the single-fan power configuration, the two dual-fan power arrangements could provide substantially large air quantities in the wind tunnel with the same inertia effect of the rotating fan as the single-fan power system. Considering the special requirements of the experimental research on aeolian   sand transport, dual-fan arrangements, such as that in case F2s, is appropriate for generating a boundary layer flow with high wind velocity and turbulence. A set of carefully arranged roughness elements is necessary to adjust the boundary layer flow structure and enhance the boundary layer thickness [63][64][65]. One wind tunnel has been created with a dual-fan power system similar to that in case F2s [66]. The dimension of the tunnel's experimental section is 12,000 mm L × 600 mm W × 800 mm H . Two axis fans with 3.5 kW power and 1440 rpm are individually fixed at both ends of the tunnel. The declining roof board of the tunnel is removed to offset the air-blowing fan effect. This step would generate stream-wise additional pressure in the middle part of the tunnel. Figure 17 shows the airflow velocity profile at half-length of the tunnel with/without roughness element arrangement in the tunnel while the two fans individually or simultaneously run. The figure also demonstrates that the effect of the rotating airblowing fan influences the terminal flow state of the airflow driven by the dual fans. The roughness element arrangement is more efficient in enhancing the boundary layer thickness compared with the velocity profiles shown in Figure 17. The boundary layer thickness in the tunnel with roughness elements is more than twice that in the case without these elements. The roughness element arrangement could also reduce the effect of the rotating air-blowing fan on the airflow state.

Conclusions
Boundary layer flows in the wind tunnel driven by single/dual fans have been successfully simulated by modelling the entire wind tunnel construction, including the power fan sections. This modelling method in the numerical simulation can ensure that the flow features of the boundary layer flow in the wind tunnel are efficiently maintained along the stream-wise airflow. This method is also effective in numerically simulating the boundary layer flow in the ABL wind tunnel without considering the inconsistency between the fully developed ABL inlet profiles and the wall function formulation. Compared with the simulation result of the boundary layer flow driven by a single fan in the wind tunnel, the airflow in the tunnel with dual fans possesses considerably high velocity and the same inertia effect of the rotating fan. When dual fans are both placed at the tunnel outlet, the boundary layer flow in the tunnel is similar to that driven by a single fan. If dual fans are separately placed at the tunnel inlet and outlet, then the airflow can be different from the two above-mentioned flows. This occurrence can be attributed to the effect of the air-blowing fan at the tunnel inlet. The maximum airflow velocity and turbulence intensity of the boundary layer flow in such a separate dual fan configuration can correspondingly reach 53.393 m/s and 12.02%, which are 1.75 and 1.49 times higher than those in the single-fan configuration one, respectively. Nevertheless, such a fan layout scheme may cause reductions in boundary layer thickness to a certain extent. Considering the requirements of wind velocity and turbulent intensity to move a mass of sand continuously, such as field aeolian sand transport, the wind tunnel with dual fans at both ends is suitable for the experimental study of aeolian sand transport. The tunnel construction adjustment and optimal roughness element arrangement should also be adopted during wind tunnel improvement and practical application to offset the effect of air-blowing fans. Accordingly, the required boundary layer thickness in the wind tunnel can be guaranteed.

Disclosure statement
No potential conflict of interest was reported by the author(s).