Numerical Study of the Flow Characteristics of Downburst-like Wind over the 3D Hill Using Different Turbulence Models

Abstract

With the rapid development of computational fluid dynamics (CFD) technology, it has been widely used to study the wind field characteristics of downbursts in mountainous areas. However, there is little guidance on the selection of different turbulence models for simulating downburst wind fields over hills using CFD, and few comparative studies have been conducted. This paper used nine turbulence models to simulate the wind field of a downburst over a 3D quadratic ideal hill. The simulated values of average and transient winds were compared with wind tunnel test data, and the flow characteristics at different moments under a downburst were analyzed. The flow characteristics in the wake region of the downburst over the hill are also quantitatively analyzed using the proper orthogonal decomposition (POD) method. The results show that approximately 85% of the results from the LES and REA models fall within a 30% error range, so the large eddy simulation (LES) model and the realizable k-ε model (REA) are more accurate in simulating the mean wind field, and the transient wind field simulated by the LES model is also in good agreement with the experimental data. In addition, this paper reveals the evolution mechanism of the transient wind field structure over a hill model under a downburst and finds that the first-order mode obtained by POD may be related to the acceleration effect on the hilltop.

1. Introduction

During thunderstorms, the exchange of hot air near the surface with cold air aloft can trigger a meteorological phenomenon known as a downburst [1]. A downburst is characterized by a high-velocity downward airflow that transforms into horizontal motion upon hitting the ground (Figure 1). As the horizontal flow spreads, it can generate extremely destructive winds that may cause severe damage to engineering structures such as buildings and transmission line towers [2,3,4,5]. Downbursts frequently occur in mountainous areas where the acceleration effect of wind velocity at the hilltop can exacerbate the hazards caused by these storms. Thus, studying the wind field characteristics in mountainous areas under downbursts is of great practical significance.

The main research methods used to study the wind field characteristics of downbursts in mountainous regions include field measurements, wind tunnel tests, and numerical simulations. Although field measurements [6,7,8,9,10,11] can provide authentic wind field data, the random and short-lived nature of downbursts in time and space [12] presents challenges in obtaining adequate data. Furthermore, field measurements entail significant labor and financial investment. Hjelmfelt [13] summarized the structure and characteristics of downbursts from field measurements and noted that mean velocity profiles obtained through the impinging jet model in wind tunnel tests are comparable to field-measured data. Given these limitations, researchers have conducted wind tunnel tests based on the impinging jet model to simulate downbursts and investigate their flow characteristics [14,15,16,17]. Wind tunnel tests offer lower costs and shorter experimentation cycles while maintaining accuracy and persuasiveness compared to field measurements. However, wind tunnel tests are restricted by size and cannot provide complete information on the flow field [18], limiting their ability to replicate real-world wind fields in complex mountainous terrain.

Figure 1. Flow field structure of a downburst [13].

Figure 1. Flow field structure of a downburst [13].

In recent years, with the rapid development of computer hardware and software, the use of CFD has become increasingly popular [18,19]. As a result, many researchers have begun to use various CFD turbulence models to simulate the wind field of downbursts in flat areas. Based on the impinging jet model, Kim and Hangan [20] used the Reynolds-averaged Navier–Stokes (RANS) method to study the wind characteristics of downbursts in flat areas and obtained a reasonable simulation of the downburst wind field. Qu and Ji [21] used the Reynolds stress model (RSM) to conduct unsteady simulations of the formation and spread of downbursts. Their numerical results show that wind velocity varies greatly throughout the formation and diffusion process, with the maximum wind velocity occurring when the downward airflow hits the ground. Aboshosha et al. [22] used LES to describe the flow fields of downbursts in four different exposure areas: open air, rural areas, suburban areas, and urban areas. Their results show that the mean wind and fluctuating wind obtained by numerical simulations are consistent with field-measured data and that peak wind velocity decreases and the height corresponding to peak wind velocity increases as roughness increases. Sengupta and Sarkar [23] performed three-dimensional numerical simulations based on the impinging jet model and used different turbulence models such as standard k-ε, renormalization group k-ε (RNG), realizable k-ε, shear stress transport k-ω (SST), RSM, and LES to study the effect of turbulence models on mean wind velocity. A comparison with experimental data shows that the simulated values of LES, REA, and RSM are more consistent with experimental data than those of the other three models. Haines and Taylor [24] conducted numerical simulations using the impinging jet model and four different turbulence models, including the Improved Delayed Detached-Eddy Simulation k-ω (IDDES k-ω), SST, Scale Adaptive Simulation (SAS), and LES. Their goal was to determine the most suitable turbulence model for simulating downburst outflow. The results showed that the dynamic LES turbulence model provided simulated values closest to experimental data and accurately captured velocity at different positions. Khayrullina et al. [25] conducted a validation study of the steady-state RANS turbulence model for an impinging jet at two different Reynolds numbers. They focused on the results provided by five different turbulence models: standard k-ε, REA, RNG, SST, and RSM. Their findings indicated that the accuracy of each turbulence model varied at different locations of the jet and that the REA model predicted velocity with the best overall performance. Žužul et al. [26] performed numerical simulations using three meshing methods and five turbulence models: standard k-ε, REA, RNG, standard k-ω, and SST. Their results showed that a fine mesh and the SST model produced simulation results most consistent with experimental data.

In addition to studies on downburst wind fields in flat areas, there has also been research on downburst flow over hills. Selvam and Holmes [27] used the impinging jet model and the standard k-ε model to simulate the steady-state wind field of a downburst over a hill. Their results showed good agreement between simulated values and field-measured data. They also found that the wind velocity at the hilltop under a downburst was smaller than that under an atmospheric boundary layer wind field. Letchford and Illidge [15] studied the effects of radial distance and hill slope on the acceleration factor of a cosine-shaped hill based on impinging jet experiments and numerical simulations. Their results showed reasonable agreement between numerical simulations, wind tunnel tests, and field measurements. They also found that the acceleration factor at the hilltop was directly proportional to the hill slope and inversely proportional to the radial distance. Wood et al. [28] used CFD to simulate the wind field over an embankment under a downburst. Their research showed that speed-up factors were similar at the crest of the embankment but dropped off rapidly behind the crest compared with boundary layer flow. Mason et al. [29] used the impinging jet model to conduct experiments and numerical simulations on the wind field of a downburst over three topographic features: cliffs, triangular hills, and bell mountains. They studied the influence of parameters such as hill slope, radial position, and the diameter of the downburst itself. Their results showed that the acceleration factor was greatly affected by the hill slope and that its maximum acceleration effect could reach about 30%. Mason et al. [30] performed three-dimensional numerical simulations of downburst events over two topographical features (cliffs and bell mountains) using a non-hydrostatic sub-cloud model. They then examined the effect of topography on simulated wind velocity. They estimated the acceleration factor of the downburst and compared it with the acceleration factor of boundary layer wind, showing that the acceleration of downburst wind was lower than that of boundary layer wind. Abd-Elaal et al. [4] used the Detached Eddy Simulation (DES) method to numerically simulate a downburst over two real terrains. They studied the horizontal and vertical wind speed distribution of the downburst as it passed over the real terrain. Their results showed a significant change in the horizontal wind speed profile above the summit, with the corresponding downburst horizontal wind speed increasing by 25% over flatlands. Yan et al. [31] numerically simulated the downburst wind field on a quadratic curved hill using three RANS turbulence models: REA, RNG, and SST. By comparing their results with wind tunnel test data, they found that the REA model was optimal.

By reviewing the research on numerical simulations of flat wind fields and wind fields over hills under downbursts, it is evident that while there are many studies on flat wind fields simulated using different turbulence models, there are fewer studies on wind fields over hills. The main contributions and novelties of this study are shown below.

(1)

This work aims to provide a more comprehensive reference for the selection of turbulence models. The most suitable turbulence model for simulating the mean and transient wind fields when downburst over hills is determined.

(2)

The evolution mechanism of the transient wind field structure over the quadratic ideal hill under a downburst is revealed.

(3)

The flow separation characteristics and the acceleration phenomena of wind fields over hills and flatlands are analyzed, and the flow separation obtained by different turbulence models is compared.

The rest of this paper is organized as follows. Section 2 introduces the parameter settings for the numerical simulation, including the turbulence model, computational domain, boundary conditions, grid division, working condition setting, and solution algorithm. This section also provides grid independence verification and a description of the wind tunnel test layout. In Section 3, different turbulence models are used to carry out numerical simulations and their results are compared with impinging jet model tests. The accuracy of each turbulence model is then analyzed both qualitatively and quantitatively. Additionally, the turbulent flow characteristics are further analyzed using the POD method. Finally, the conclusions of this study are presented in Section 4.

2. Numerical Simulations and Wind Tunnel Tests

Since the impinging jet model can simulate the entire development process of a downburst and takes into account the wall jet region, the wind field of a downburst over a hill is studied by a numerical simulation method based on the impinging jet model. This section introduces the solved governing equations and selected turbulence model and describes the computational domain, boundary conditions, and computational settings used in the numerical simulation.

2.1. Governing Equations

In the field of structural wind engineering, the maximum wind speed is less than 0.3 Ma, and the atmospheric density is regarded as a constant and remains constant. It is generally considered as an incompressible flow, which should satisfy the law of conservation of mass, momentum, and energy. However, it is considered to be in isothermal stratiform. The air temperature is regarded as constant, and the energy change is ignored. Therefore, in the numerical simulation calculation, the energy conservation equation is generally not considered, and only the mass and momentum conservation equations need to be satisfied. The conservation equations of mass and momentum are shown in the following equation.

$$\frac{\partial u_i}{\partial x_i}=0$$

(1)

$$\rho \frac{\partial u_i}{\partial t}+\rho \frac{\partial u_i\partial u_j}{\partial x_j}=\frac{\partial }{\partial x_j}\left(\mu \frac{\partial u_i}{\partial x_j}\right)-\frac{\partial p}{\partial x}+S_i$$

(2)

where $u_i$ represents the components of fluid velocity vector in $x_i$ directions, $\rho$ is fluid density, $\mu$ is the dynamic viscosity, $p$ is the pressure, and $S_i$ represents the generalized source term of the momentum equation in the three velocity directions.

2.2. Turbulence Model

This section introduces the turbulence models used in CFD numerical simulations, which are mainly divided into three categories: RANS, LES, and DES. It also explains the differences and connections between these models.

2.2.1. RANS Method

RANS is the most widely used method for numerical simulation due to its short computational time and low computational cost. The core idea is to avoid solving the transient Navier–Stokes equations directly by time-homogenizing them using the statistical theory of turbulence, thus reducing computational effort. However, this process leads to the appearance of Reynolds stress, causing the system of equations to be unclosed. Therefore, the turbulence model is needed to enable the system of equations to be solved. Here, three RANS turbulence models, RNG, REA, and SST, are selected for the numerical simulation of wind fields over the hill based on the impinging jet model. The following is a brief description of these three turbulence models.

The RNG model was proposed by Yakhot and Orzag [32]. It extends the standard k-ε model by using a renormalization group approach and has the following main features: an additional term R_ε is added to the transport equation to improve the accuracy of fast strain flow; an analytically derived differential formulation of the effective viscosity is added to account for low Reynolds number effects. These characteristics make the results of the RNG model more accurate than those of the standard model, and RNG is more widely used.

The REA model was proposed by Shih et al. [33]. Based on the standard k-ε model, the REA model adds alternative formulas for turbulent viscosity and corrects the transport equation for the dissipation rate ε. This allows for the accurate prediction of flow field information in simulations of walls with large curvature or inverse pressure gradient and flow separation.

The SST model was improved from the standard k-ω model by Menter [34]. It revises the definition of turbulent viscosity and mainly considers the transport of the main turbulent shear stress. It can solve the standard turbulent kinetic energy k equation and the equation of a specific turbulence rate ω. The SST model combines the accuracy of the standard k-ω model in the near-wall region with the independence of the standard k-ε model in inlet-free turbulence. As a result, this model produces more accurate results than both the standard k-ω and standard k-ε models.

The RANS method is based on the statistical averaging of flow by the turbulence model. However, it cannot capture the pulsatility of turbulence or accurately simulate high-intensity linear flow. Additionally, its turbulent kinetic energy pulsation spectrum cannot be resolved.

2.2.2. LES Method

The LES model was proposed by Smagorinsky in 1963 [35]. The fundamental concept of LES is to filter turbulent eddies using a characteristic grid. Eddies larger than the spatial filter grid are solved directly using the Navier–Stokes equation, while eddies smaller than the spatial grid are modeled using a subgrid-scale model. In this case, the standard Smagorinsky–Lilly (SM) subgrid-scale model is employed, as shown in the following equation.

$${\tau }_{ij}-\frac{1}{3}{\tau }_{kk}{\delta }_{ij}=-2{\left(C_S\mathsf{\Delta }\right)}^2\left|\overline{S}\right|{\overline{S}}_{ij}$$

(3)

where τ_ij denotes the sublattice scale stress, ${\tau }_{ij}=\overline{u_i{u}_j}-{\bar{u}}_i{\bar{u}}_j$, and v_SGS denotes the sublattice vortex viscosity coefficient, and the expression is

$$v_{SGS}=l_S\left|\overline{S}\right|={\left(C_S\mathsf{\Delta }\right)}^2\left|\overline{S}\right|$$

(4)

where ${\overline{S}}_{ij}=\frac{1}{2}\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)$; $\left|\overline{S}\right|=\sqrt{2{\overline{S}}_{ij}{\overline{S}}_{ij}}$, $\Delta ={\left({\Delta }_x{\Delta }_y{\Delta }_z\right)}^{\frac{1}{3}}$; $l_s$ is the Smagorinsky length scale; $\left|\bar{S}\right|$ is the norm of the strain rate after filtering; ${\overline{u}}_i$, ${\overline{u}}_j$ indicates the filtered velocity; ${\Delta }_i$ denotes the grid size in direction I; and C_S is the Smagorinsky constant, which takes the value of 0.1 in this paper.

2.2.3. DES Method

DES is a hybrid method that combines the RANS and LES methods to simulate turbulent flow within a flow field. RANS is utilized in the near-wall region, while LES is employed in regions of flow separation and wake. Compared to LES, DES reduces the requirement for grid resolution and enhances computational efficiency. In contrast to RANS, DES improves simulation accuracy within the separation region. Various turbulence models are employed by different DES methods and are briefly described below.

Spalart et al. [36] proposed the Detached Eddy Simulation Spalart–Allmaras (DES-SA) model, which is based on the Spalart–Allmaras equation. This model employs a hybrid length scale L_DES, instead of the wall distance d. The hybrid length scale is defined as the minimum value between the local grid spacing and the near-wall distance. The expression for the hybrid length scale is as follows:

$$L_{\mathit{DES}}=\min \left(L_{\mathit{RANS}}{,C}_{\mathit{DES}}\Delta \right)$$

(5)

where C_DES represents the scale factor, generally taken as C_DES = 0.65; Δ is the sublattice scale, taken as ∆ = max {∆x, ∆y, ∆z}; and Δx, Δy, Δz are the grid scales in the x, y, z directions, respectively.

The Detached Eddy Simulation Realizable k-ε (DES-RK) model employs the transport equations for k and ε from the REA model to compute the eddy viscosity in the RANS region and solve for the subgrid-scale viscosity in the LES region. In the DES-RK model, the dissipation term in the transport equation for turbulent kinetic energy, k, is modified to $k^{3/2}{/L}_{\mathit{DES}}$, where $L_{\mathit{DES}}$ represents the hybrid length scale. The expression for the hybrid length scale is as follows:

$$L_{\mathit{DES}}=\min \left(L_{\mathit{REA}}{,C}_{\mathit{DES}}\Delta \right)$$

(6)

where the value of the sublattice scale Δ is the same as that of Equation (5), and the value of the scale factor C_DES is 0.61.

The Detached Eddy Simulation Shear Stress Transport (DES-SST) model was proposed by Menter [34] and is based on the SST model. The DES-SST model incorporates the hybrid length scale L_DES. In the SST model, the length scale is defined as $l_{\mathit{SST}}=\sqrt{k}{/\beta }^{*}\omega$, while in the DES-SST model, the hybrid length scale is defined as follows:

$$L_{\mathit{DES}}=\min \left(L_{\mathit{SST}}{,C}_{\mathit{DES}}\Delta \right)$$

(7)

where the subgrid scale Δ and the scale factor C_DES are the same as in Equation (6).

Menter et al. [34] proposed the Delayed Detached Eddy Simulation SST k-ω (DDES-SST) model with a hybrid length scaled by

$$L_{\mathit{DES}}=\min \left(L_{\mathit{SST}},\frac{C_{\mathit{DES}}\Delta }{{1-F}_{\mathit{SST}}}\right)$$

(8)

Shur et al. [37] developed an Improved Delayed Detached Eddy Simulation SST (IDDES-SST) model, which improves the sublattice scale Δ, the expression of which is shown below.

$$\Delta =\min \left\{\max \left[C_w{d}_w{, C}_w{h}_w{, l}_w\right]{, h}_{\max }\right\}$$

(9)

where $C_w=0{.15, h}_{\max }=\max \left\{\Delta x, \Delta y, \Delta z\right\}$ and has the same meaning as the sublattice scale in Equation (5); l_w is the total length of the grid for wall normal; and d_w is the distance between the midpoint of the flow field and the nearest wall.

This section provided a brief overview of the differences between various turbulence models. Both the DES and LES methods are capable of simulating unsteady wind fields, while the RANS method is limited to simulating steady wind fields. A comparison of these three methods reveals that the RANS method has the lowest computational cost but correspondingly fewer degrees of freedom in its analysis. In contrast, the LES and DES methods require longer computational times but offer higher simulation accuracy.

2.3. Computational Domain and Boundary Condition Setting

The scale ratio of numerical simulation is 1:1000, and ANSYS FLUENT 19.0 is used to carry out the numerical simulation of the downburst over a three-dimensional hill. To reduce computational cost, a one-fourth circle is used to model the computational domain, as shown in Figure 2. The outlet diameter is D = 600 mm, and the distance from the outlet nozzle to the ground is H = 1.0D = 600 mm. The center of the jet to the farthest radial distance position is 7.0D, and the outlet wind velocity, u_j, is 10 m/s. In the numerical simulation, the hill model is placed at r/h = 0 (where r represents the radial coordinate, and h represents the height of the hill model) and at a circular position of 45°. The hill model is quadratic, as shown in Equation (10).

$$z=h\left(-\frac{\left(x^2{+y}^2\right)}{L^2}+1\right)$$

(10)

where L represents the horizontal distance from the foot of the hill to the hilltop and x and y represent the distances to the center of the hill. In this study, two hill models of different sizes are employed. The Quad-D300-H075 hill model has a diameter d_h of 300 mm and a height h of 75 mm. The Quad-D176-H075 hill model has a diameter d_h of 176 mm and a height h of 75 mm.

As shown in

Table 1. Boundary conditions settings.

Position	Boundary Conditions
Inlet	Velocity inlet	Inlet wind velocity uj = 10 m/s Turbulence intensity = 1% Inlet diameter D = 600 mm
Outlet and top	Pressure outlet	Turbulence intensity = 1%
Sides	Symmetry	$v=0, \partial \left(u, v, w, k, \mathsf{\epsilon }\right)/\partial z=0$
Bottom	No-slip wall	$u=v=w=0, \partial \left(u, v, w, k, \mathsf{\epsilon }\right)/\partial z=0$
Wall	No-slip wall	$u=v=w=0, \partial \left(u, v, w, k, \mathsf{\epsilon }\right)/\partial z=0$

, the velocity inlet condition is used for the inlet, the pressure outlet condition for the outlet and top, the symmetric boundary condition for the sides, and the no-slip wall condition for the bottom and wall of the computational domain.

2.4. Grid Division and Working Condition Setting

As shown in Figure 3, this study employs a structured grid method. The near-wall region exhibits a low Reynolds number flow. To accurately capture the complex turbulent characteristics of the flow field around the structure, the grid in the near-wall region of the hill model is refined. In regions with high wind velocity gradients, such as near-wall and wake regions of the hill model, adjacent grid sizes should not differ significantly. The grid growth rate is set at 1.05. Different turbulence models correspond to different grid numbers, and the grid settings for different working conditions are shown in

Table 2. Boundary conditions settings.

Turbulence Model	Number of Grids (×106)	First Layer Grid (m)	Y+
SST	5.6	0.002	40
RNG	5.6	0.002	40
REA	5.6	0.002	40
LES	8.3	0.0005	10
DES-SA	8.3	0.0005	10
DES-RK	8.3	0.0005	10
DES-SST	8.3	0.0005	10
DDES-SST	8.3	0.0005	10
IDDES-SST	8.3	0.0005	10

. Y⁺ is a dimensionless value commonly used in boundary layer theory and is defined by the following equation:

$$Y^{+}=\frac{u^{*}{}_{}y}{v}$$

(11)

where $u^{*}$ is the wall shear velocity, y is the height of the first grid from the wall, and v is the kinematic viscosity of the fluid. The size of Y⁺ reflects the height of the first grid from the wall. A large Y⁺ indicates that the first grid node is further from the wall.

2.5. Solution Algorithms

In the numerical simulation, the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm is used to solve the velocity–pressure coupled equations. The nonlinear convection term is discretized using the Second Order Upwind format, while the momentum equation is discretized using the Bounded Central Differencing format. Time discretization is performed using the Second Order Fully Hidden format. The time it takes for the wind to cross the computational domain from inlet to outlet is referred to as a full-basin time. A total of 30 full-basin times are simulated, with results from the latter 24 full-basin times extracted for statistical analysis to ensure the stability of the statistical results.

2.6. Grid Independence Verification

To verify the effect of grid size on the results of the hill-based impinging jet model, this section uses the REA model as an example for grid size verification. On the medium grid of REA, two sets of grids with different sizes were recreated: a coarse grid with a first-level height of 0.005 m and a total of 3 million grids, and a fine grid with a first-level height of 0.0005 m and a total of 8.8 million grids. Figure 4 shows the near-ground comparison of different grid sizes.

Figure 5a depicts the vertical wind profile of horizontal velocity at various measurement points, using three distinct grid numbers. The figure reveals a high degree of conformity in wind velocity at each position under the three grid numbers. Notably, in the near-ground portions of the hillside on the windward side, findings derived from the fine and medium grids were in agreement, but those of the coarse grid deviated from those of the medium grid. Figure 5b compares the error analysis of the coarse and fine grids with that of the medium grid. Both the coarse and medium grid computations exhibited errors, notably when the wind velocity was below 4 m/s. Error rates exceeding 5% were observed at certain measurement points. On the other hand, the fine grid yielded results aligned with those of the medium grid. Similar validation outcomes were obtained using other turbulence models. Based on these results, selecting the medium grid for numerical simulation in subsequent working conditions ensures accuracy, improved computational efficiency, and reduced computational costs.

2.7. Wind Tunnel Tests

To verify the accuracy of the numerical simulations, wind tunnel tests based on the impinging jet model were performed. The impinging jet test was conducted at the Structural Wind Engineering and Urban Wind Environment Laboratory of Beijing Jiaotong University. Figure 6 shows the impinging jet device simulating the wind field of a downburst. The outlet diameter of the downburst simulator is D = 600 mm, and the height of the nozzle to the floor can be adjusted using a hydraulic device in the range of 400 mm–800 mm. The maximum wind velocity at the outlet, u_jmax, is 12 m/s.

In this test, a three-dimensional wind velocity tester, the Cobra Probe produced by Turbulent Flow Instrumentation (TFI), and a three-dimensional side-shifting frame were used to measure wind velocity at different radial positions. The Cobra Probe is a four-hole pressure probe with a wind velocity measurement range of 2–100 m/s and a typical measurement accuracy of ±0.5 m/s. The sampling frequency was set to 256 Hz, and the sampling time was 30 s. The mean of three measurements was processed to reduce error. All other relevant settings of the impinging jet test remained the same as in the numerical simulation.

3. Analysis of Wind Field Results over the Hill

In this section, the simulated mean wind velocity of different turbulence models is compared to the experimental wind velocity to verify the accuracy of the numerical simulation through qualitative and quantitative analysis. In addition, the characteristics of the transient wind field and its flow separation for a hill under a downburst are analyzed. Unless otherwise specified, the hill model used in this section is the Quad-D300-H075 hill model.

3.1. Comparison of Mean Wind Fields with Different Turbulence Models

Figure 7 shows the mean wind velocity near the hill model computed by different turbulence models and the corresponding experimental data. From the figure, it can be seen that the simulated wind velocity obtained from each turbulence model follows the same trend as the experimental wind velocity distribution. At the front of the hill on the windward side, the simulated wind velocity obtained from each turbulence model is consistent with the experimental wind velocity, because they are not influenced by the hill model. However, the simulated wind velocity obtained from the DES-SA model decays faster than the experimental wind velocity and those obtained from other turbulence models. At the foot of the hill on the windward side, each turbulence model underestimates the experimental wind velocity near the wall. On the hillside on the windward side, each turbulence model overestimates the wind velocity near the wall, with those from the DES-SA, SST, and RNG models decaying faster at higher positions. At the hilltop, there is a significant increase in wind velocity near the wall, which tends to exacerbate its disaster-causing properties. To better compare the wind velocity at the hilltop, the mean wind velocity profile is shown in Figure 7b. It can be seen that the simulated wind velocity obtained from the LES and DES-RK models is generally consistent with the experimental wind velocity overall. In the near-wall region, the simulated wind velocity of the DES-SST and DDES-SST models is lower than the experimental wind velocity, while those of the DES-SA, SST, and RNG models overestimate the wind velocity and decay faster with increasing vertical height. The simulated wind velocity of the IDEES-SST model at the near-wall position matches well with the experimental wind velocity but underestimates it at higher vertical positions. On the hillside on the leeward side, the simulated wind velocity differs significantly from the experimental wind velocity, likely due to Cobra’s inability to measure reverse flow values. At the foot of the hill on the leeward side, the simulated wind velocity obtained from the SST, RNG, and REA models agrees well with the experimental wind velocity near the wall, while those obtained from other turbulence models underestimate it. Behind the hill on the leeward side, the simulated wind velocity obtained from the DES-SA, DES-SST, DDES-SST, and IDDES-SST models all underestimate the experimental wind velocity, while those obtained from RNG and SST overestimate it. In general, the simulated wind velocity obtained from LES and DES-RK are in good agreement with the experimental wind velocity.

To further analyze the accuracy of the simulated wind field of each turbulence model, Figure 8 shows a scatterplot between the experimental and simulated wind velocities. From the figure, it can be seen that the error of some simulation results obtained by each turbulence model is greater than 30%, with the experimental wind velocity is more underestimated. This is consistent with the phenomenon observed in Figure 7a, where the simulated wind velocity away from the wall region is less than the experimental wind velocity. Comparing wind velocity values of different SST turbulence models in Figure 8b–d,i, it can be seen that overall, the results of the IDDES-SST model are more accurate than those of DES-SST, DDES-SST, and SST. At positions with a large wind velocity (u/u_j > 0.8), the results of IDDES-SST are better than those of DES-SST, DDES-SST, and SST, with an error range concentrated within 10%. In the wind velocity range of 0.4 < u/u_j ≤ 0.8, there are more points distributed when the error range of the DES-SST and DDES-SST simulated wind velocity is greater than 30%. Comparing different DES turbulence models in Figure 8b,e,f, it can be seen that the results of DES-RK are within an error range of 10%, with more distribution points and a higher accuracy. Comparing different RANS turbulence models in Figure 8g–i, it can be seen that the error distribution of wind velocity values obtained by the SST model is wide, while that obtained by the REA model is distributed within a 10% error range, indicating that the REA model has better results than other RANS models. In summary, for different SST turbulence models, their accuracy is ranked as IDDES-SST > DES-SST > DDES-SST > SST; for different DES turbulence models, their accuracy is ranked as DES-RK > DES-SST > DES-SA; and for different RANS turbulence models, their accuracy is ranked as REA > RNG > SST.

To evaluate the accuracy of different turbulence models,

Table 3. Analysis of velocity error between different turbulence models and the experiment.

Error	LES	DES-SST	DDES-SST	IDDES-SST	DES-RK	DES-SA	REA	SST	RNG
10%	50%	31.43%	28.57%	44.29%	41.43%	38.57%	57.14%	38.57%	40%
20%	71.43%	50%	45.71%	74.29%	70%	51.43%	75.71%	54.29%	70%
30%	84.29%	64.29%	61.43%	81.43%	81.43%	67.14%	85.71%	71.43%	81.43%
MNB	−3%	−15.94%	−15.86%	−2.83%	−3.85%	−13.47%	2.04%	2.95%	5.94%

presents the percentage of data points, where the difference between the simulated wind velocity (obtained from different turbulence models) and the experimental wind velocity falls within error ranges of 10%, 20%, and 30%. Additionally, the table includes the Mean Normalized Bias (MNB) for both models. The MNB represents the mean relative error between the simulated and experimental wind velocities and is defined as follows.

$$\mathrm{MNB}=\frac{1}{N}\sum _{i=1}^N\left(\frac{V_{i,\mathit{fitting}}-V_{i,\mathit{exp}}}{V_{i,\mathit{exp}}}\right)$$

(12)

where N is the number of measurement points in the experiment and V_i,exp and V_i,fitting are the wind velocity values obtained from the experiment and turbulence models, respectively. MNB is the metric that can synthetically characterize the degree of change in the data. When the MNB is negative, the simulated wind velocity is smaller than the experimental wind velocity as a whole, and when the MNB is positive, its wind velocity is larger than the experimental wind velocity.

An analysis of the errors reveals that approximately 85% of the results from the LES and REA models fall within a 30% error range. The largest discrepancy between the simulated wind velocity (obtained from the DES-SST and DDES-SST models) and the experimental wind velocity is observed. In general, all turbulence models, except for the RANS model, underestimate the experimental wind velocity. Focusing on the hilltop location, the error values of the wind velocity at this location are analyzed for different turbulence models. Figure 9 shows that when the wind velocity values at near-surface locations are large due to the acceleration effect at the hilltop, the results of each turbulence model fall within a 10% error range. When the wind velocity values are small, the wind velocity obtained from the RNG, SST, and DES-SA models is smaller than the experimental wind velocity, with error confidence intervals within 20%. In summary, the LES and REA turbulence models provide the most accurate simulation results for the mean wind field on a hill, followed by the DES turbulence model.

Based on the above research, in the numerical simulation of the mean wind field of the downburst over the hill, a more detailed simulation effect of the turbulence model is given. In different SST turbulence models, the accuracy effect is IDDES-SST > DES-SST > DDES-SST > SST; in different DES turbulence models, the accuracy effect is DES-RK > DES-SST > DES-SA; and in different RANS turbulence models, the accuracy effect is REA > RNG > SST. On the whole, the effect of the turbulence model on the simulation results of the mean wind field on hills is investigated. It is found that the REA turbulence model is comparable to the LES simulation results, followed by the DES turbulence model. Our results are roughly consistent with the results of the mean wind field simulated by Sengupta and Sarkar [23] with the downburst flowing over flatlands. Sengupta and Sarkar [23] compared LES with various RANS turbulence models and found that LES and REA had better results. In this paper, the simulation effects of LES, DES, and RANS are given, namely, LES is approximately equal to REA and greater than DES, which provides a more comprehensive reference for turbulence model selection. It also shows that in the numerical simulation of the mean wind field of downbursts over flatlands and hills, the effects of the LES and REA turbulence models are relatively consistent. They are not affected by conditions such as flat lands or hills. In addition, the conclusion that REA is better than RNG and SST in the RANS simulation method is consistent with the results of Khayrullina et al.’s [25] numerical simulation of a downburst flowing over a flatland and Yan et al.’s [29] numerical simulation of a downburst flowing over a hill. It further shows that the simulation effect of the turbulence model is not affected by terrain conditions.

Figure 10 presents the dimensionless mean wind velocity contours (u/u_j) obtained from each turbulence model. The figure shows a consistent overall distribution trend of the mean wind velocity among the models. An acceleration effect near the ground is noticeable from the hillside on the windward side to the hilltop. At the hilltop, the simulation results of the DES-SST, DDES-SST, and IDDES-SST turbulence models exhibit essentially the same acceleration regions, with ranges slightly smaller than those of other models. Notably, the SST model portrays a larger acceleration region than all the other models. Meanwhile, at the foot of the hill on the leeward side, the computed results of the DES-SST, DDES-SST, and DES-SA models reveal conspicuous zones of low wind velocity. Further, the computed distributions using the LES, SST, RNG, and REA models at locations farther away from the hill on the leeward side are similar, indicating a greater wind velocity than that at locations corresponding to other models. These findings suggest that the DES-SST, DDES-SST, and IDDES-SST turbulence models are preferable for the simulation of hill terrain flows in situations involving an acceleration effect near the ground, whereas the SST model may overestimate the acceleration region. Conversely, the DES-SST, DDES-SST, and DES-SA models may better predict low wind velocity regions.

3.2. Structural Characteristics of Transient Wind Fields with Different Turbulence Models

Figure 11 illustrates the pulsation component std(u)/u_j of the dimensionless horizontal wind velocity computed from each turbulence model, where std(·) denotes the standard deviation. Overall, the computed results of each turbulence model are in good agreement with the experimental results. The pulsating wind velocity computed by the DES-SA model is smaller than that of other turbulence models in the windward region and at the foot of the hill on the leeward side. The computed results of the DES-SST and DDES-SST models at the hilltop are higher than the experimental data, which may be due to an inaccurate dissipation prediction by the turbulence model. In general, the LES turbulence model can simulate a pulsating wind velocity with higher accuracy.

To obtain more intuitive results, the differences in flow separation for each turbulence model under the hill are explored by comparing the vorticity contours at different moments under each turbulence model in Figure 12. Since this paper only studies the effect of the hill on the wind field of a downburst, the impinging phase of the impinging jet (i.e., when it is not influenced by the hill model) and the fully developed phase of the flow field after completely overtopping the hill are not of concern.

Figure 12 shows the vorticity contours of the wind field obtained from simulations using different turbulence models. At T = 0.2 s, as shown in Figure 12a, the downburst impacts the hill and the results of each model have the same impinging position. The LES and DES-RK models separate the secondary vortex before the other models. At T = 0.24 s, as shown in Figure 12b, the structure and size of the separation vortex obtained by the DES-SST, DDES-SST, and IDDES-SST models are essentially the same, while the vortex structures obtained by the LES model are larger. At T = 0.28 s, as shown in Figure 12c, except for the DES-RK model, the flow obtained by other models has a secondary vortex separation and reverse movement. The vortex structure at the hilltop obtained by the LES model is more discrete. At T = 0.32 s, as shown in Figure 12d, the separated vortex structures obtained by each turbulence model develop along the hill to the leeward side. The separated vortex structure obtained by the DES-RK and DES-SA models is concentrated in the near-wall region of the hill, while that obtained by the LES model is more discrete and has a larger distribution range. At T = 0.36 s and T = 0.4 s, as shown in Figure 12e,f, the vortexes by each turbulence model reach the foot of the hill on the leeward side. Most of the separated vortexes obtained by DDES-SST, DES-RK, and DES-SA are distributed at higher positions and show a trend of further improvement. The separated vortexes obtained by the LES, DES-SST, and IDDES-SST models are mostly distributed near the wall and have a large distribution range. By comparing the vorticity obtained by each turbulence model at different times, it can be found that the vortex structure flow separation of the DES-SST, DDES-SST, and IDDES-ST models is similar throughout the development process. The vortex structure obtained by the LES model occurs earlier than other models, and when flowing over the hill model, the vortex structure is more discrete, and its size is smaller. However, the vortex shedding position of DES-RK is close to the near-wall position during development, and the interaction between the main vortex and secondary vortex does not develop along the vertical height, and the shear stress is small.

In summary, the LES turbulence model simulates a pulsating wind velocity more accurately and can well reproduce the flow separation of the downburst flow over the hill model. Therefore, the simulation results of LES are used to compare and analyze the outflow, sinking, impinging in the wind field of the flatland and the wind field of the hill, and the process of spreading along the radial distance. Figure 13 shows the vorticity of the downburst at different moments under the flatland and its corresponding transient wind velocity fields obtained by the LES model. At T = 0.12 s, as shown in Figure 13a, during the downward impact of the airflow, it gradually expands outward due to ground obstruction, forming peak velocity and large vorticity. In addition, an obvious main ring vortex is formed near r/D = 0.75, caused by the Kelvin–Helmholtz instability (due to shear action between the impact flow and surrounding still air) [20]. At T = 0.16 s, as shown in Figure 13b, the peak velocity develops further downward and radially to near the ground, forming a high wind velocity zone in the 0.75 < r/D < 1 region. The main ring vortex also develops downward along the radial direction to the top of the high wind velocity zone. In addition, a large vortex appears near the ground in the 1.0 < r/D < 1.25 region, indicating that the peak wind velocity has a tendency to develop radially upward and flow separation occurs near the ground. At T = 0.2 s, as shown in Figure 13c, the peak wind velocity further develops in the radial direction and upward, while the flow separation phenomenon between the upper main ring vortex and near the ground is more obvious, forming a clear secondary vortex near the ground due to separation. At T = 0.24 s, as shown in Figure 13d, the peak wind velocity continues to develop to a higher position at r/D = 1.6, and the main ring vortex and secondary vortex also develop to a higher position in the radial direction. In addition, new ring vortexes and wind velocity peaks caused by the Kelvin–Helmholtz instability continue to form and develop. At T = 0.28 s, as shown in Figure 13e, as the airflow continues to develop and lift, the peak wind velocity decays, and the secondary vortex moves along the reverse direction with the rotation of the main ring vortex, also decaying. At T = 0.32 s and T = 0.36 s, as shown in Figure 13f,g, the main ring vortex continues to rotate, and the secondary vortex, resulting from flow separation, continues to move in the reverse direction and becomes progressively weaker. At the same time, main ring vortexes caused by the Kelvin–Helmholtz instability and secondary vortexes generated by flow separation are continuously generated, weakened, and disappear, representing a process of energy consumption.

Figure 14 shows the vorticity contours and their corresponding transient wind velocity contours obtained from the LES model for different moments of the downburst over the hill. At T = 0.12 s and T = 0.16 s, as shown in Figure 14a,b, since the impinging airflow is not affected by the hill model, its vorticity distribution and transient wind velocity contours are the same as those of the flatland wind field. A clear main ring vortex also appears, and flow separation occurs when the peak wind velocity develops near the ground. At T = 0.2 s, as shown in Figure 14c, as the peak wind velocity moves along the hill, the main ring vortex also moves to the hillside on the windward side and significant flow separation occurs at the contact position between the main ring vortex and the hill. At T = 0.24 s, as shown in Figure 14d, the airflow develops to the hilltop and both the peak wind velocity and main ring vortex have a tendency to develop upward. At T = 0.28 s, as shown in Figure 14e, the secondary vortex separated from the main vortex rotates in reverse at the hilltop position, similar to the occurrence and separation of secondary vortexes in flat areas. However, due to the influence of the hill model, separation and reversal of the main vortex structure occur at a higher position and its structure is more discrete. At T = 0.32 s and T = 0.36 s, as shown in Figure 14f,g, the main vortex structure advances along the hill model to the hillside and foothills on the leeward side. The secondary vortex becomes further discrete at the hilltop, its range expands, and the transient wind velocity is smaller than that of the flatland, indicating significant deceleration. As the airflow develops along the hill model, inertia causes a negative pressure zone to appear near the ground on the leeward side. In this process, new ring vortexes impact the hill, and their structure repeats the developmental spread of the main ring vortex. At T = 0.4 s and T = 0.44 s, as shown in Figure 14h,i, after the airflow passes over the hill, it continues to spread along the radial direction with small wind velocity values. The new ring vortex structure develops at the hilltop, and acceleration occurs again. The new ring vortex passes over the hill model and then tends towards a steady state.

Figure 13 and Figure 14 show the vorticity and wind velocity contours for flat land and hill conditions at different moments. The development is the same before the airflow impinges on the hill, with identical impinging and sinking phases. When the airflow impacts the hill, the wind field behaves similarly to that of the flatland, with shedding and reattachment of the main ring vortex structure and reverse rotation of the secondary vortex. However, due to the influence of the hill model, the disappearance of the ring vortex and secondary vortex along their reverse movement is less obvious. The vortex structure mainly moves to the hilltop and back of the hill and becomes more discrete. A significant acceleration effect occurs at the hilltop due to the interaction between the ring vortex structure and the rising counter-rotating secondary vortex structure.

3.3. Flow Separation of Hill

To investigate flow separation in the wake region of a hill, two hill models (Quad-D176-H075 and Quad-D300-H075) are selected for numerical simulation. The LES turbulence model is employed in this simulation, because it can more accurately simulate pulsating wind velocity and better reproduce flow separation during the entire downburst process over the hill model.

This paper takes the Quad-D176-H075 hill model as an example and presents its mean wind velocity and streamline contours in Figure 15. The figure shows that when the airflow flows over the front of the hill model, a recirculation zone appears at the foot of the hill on the windward side. This may be due to increased pressure and a large pressure gradient in this region. At the hilltop, where the pressure is lower, there is an acceleration of wind velocity. In the leeward region, due to the presence of a reverse pressure gradient, the velocity slows down near the ground, but no obvious flow separation phenomenon is observed.

In order to obtain a more comprehensive understanding of the complex flow in the wake region, wind velocity data from the plane where the hill is located were subjected to a proper orthogonal decomposition (POD) analysis. Figure 16 presents the accumulated energy of the first 40 modes in the plane, which reveals that the first mode of the Quad-D300-H075 hill model accounts for nearly 20% of the total wind energy. Similarly, the first mode of the Quad-D176-H075 hill model accounts for over 25% of the total wind energy. Furthermore, the sum of the first four eigenvalues for both hill models amount to almost 50% of the total energy, indicating the importance of the early modes. Thus, the first four modes were selected to understand the flow characteristics of the hill models.

Figure 17 shows the contours of the mean wind velocity U component and the first four modes for both hill models. A comparison of the first four modes with the mean wind velocity contours reveals large fluctuations at the hilltop and in the wake region. The first mode represents vortex shedding and is related to airflow pulsation. The amplitude of this mode reaches its maximum at and near the hilltop, with the peak amplitude near the wall of the hilltop. This may be related to the acceleration effect at the hilltop. In the second mode, there are two regions with clearly opposite velocities in Figure 17b. The region of positive velocity indicates that the airflow is flowing downstream from here, while the region of negative velocity indicates a reverse flow. That is, there is fluctuation in the same direction at the hilltop on the windward side, while opposite fluctuation occurs at the same position on the leeward side, with two peaks. These fluctuations have an S-shape and may be caused by vortex shedding on both sides of the hill. In the third mode, there are still two opposite fluctuations whose positions differ from those of the second mode, indicating a certain phase difference. In the fourth mode, the wind velocity distributions of the two hill models differ and do not agree well, but both are similar to further diffusion of the vortex structure. In general, the first three modes show a strong similarity between the two hill models, with the first mode potentially related to an acceleration effect at the hilltop. However, significant differences appear in the fourth mode.

4. Conclusions

Compared with previous studies (as shown in

Table 4. A brief review of related literature that compare turbulence models.

	References	Turbulent Model	Research Content
Flatlands	Qu and Ji [21]	RSM	The maximum wind velocity phenomenon.
	Aboshosha et al. [22]	LES	Relation between peak wind speed and roughness; Turbulence correlations in the wide frequency range characterized by the turbulent length scales.
	Sengupta and Sarkar [23]	Standard k-ε, renormalization group k-ε (RNG), realizable k-ε, shear stress transport k-ω (SST), RSM, and LES	Comparison of simulation effects of turbulence models.
	Haines and Taylor [24]	IDDES k-ω, SST, Scale Adaptive Simulation (SAS), and LES	Comparison of simulation effects of turbulence models.
	Khayrullina et al. [25]	Standard k-ε, REA, RNG, SST, and RSM	Comparison of simulation effects of turbulence models.
	Žužul et al. [26]	Unsteady RANS (standard k-ε, REA, RNG, standard k-ω, and SST) and SAS.	Comparison of simulation effects of turbulence models; The effect of the exact geometry of the bell nozzle and the hexagonal shape of the test chamber on the simulation; The physical characteristics of downburst wind in a WindEEE dome are studied from the perspective of structural wind load.
Terrain	Wood et al. [28]	Differential Reynolds stress model (DSM) and k-ε model	Acceleration phenomenon at the crest of the embankment.
	Mason et al. [30]	/	The influence that topographic features have on the near-ground wind structure of a downburst.
	Abd-Elaal et al. [4]	DES	The profiles of downburst wind speeds as they pass over real topography, and the consequent changes in horizontal and vertical downburst wind speeds.
	Yan et al. [31]	REA, RNG, and SST	Comparison of simulation effects of turbulence models; A semi-empirical model is established for the acceleration effect of the downburst wind field over a hill.
	This paper	LES, DES-SA, DES-RK, DES-SST, DDES-SST, IDDES-SST, REA, RNG, and SST	Comparison of simulation effects of turbulence models; The evolution mechanism of the transient wind field structure over the quadratic ideal hill under a downburst is revealed; The flow separation characteristics and the acceleration phenomena of wind fields over hills and flatlands.

), this paper adopted more turbulence models to verify the effectiveness of simulating the wind field of a downburst over hills. In addition, vorticity contours were used to examine the transient wind fields obtained from each turbulence model. The separation mechanism of a three-dimensional hill model under a downthrust is further studied by using the POD method. Our analysis leads to several key conclusions.

(1)

For the mean wind field over the 3D hill model, the CFD numerical simulation method can simulate the wind field characteristics of the downburst. The results of the LES model, the DES-RK model, and the REA model are more accurate than other models, and the proportion of errors within 30% is as high as 80%. Furthermore, the simulation effect of each turbulence model is REA ≈ LES > DES.

(2)

By comparing the fluctuating wind speed, it was found that the simulation performance of each turbulence model in the transient wind field was ranked as follows: LES > DES. When comparing flow separation at the same time, it was found that the flow separation obtained by the DES-SST, DDES-SST, and IDDES-SST models showed little difference. The ring vortex structure simulated by the LES model underwent flow separation earlier than other models, and the vortex structure was more discrete when flowing over the hill model. Additionally, the vortex shedding position from the DES-RK turbulence model was concentrated near the wall, and the interaction between the main vortex and secondary vortex did not develop along the vertical height.

(3)

By comparing the development processes of downbursts in flatlands and hills, it has been found that the main stages of development are the same. The maximum wind speed is achieved through the interaction between the primary ring vortex structure and the rising counter-rotating secondary vortex structure. The key difference is that the presence of a hill causes the separation position of the secondary vortex to rise, resulting in a noticeable acceleration effect at the hilltop.

(4)

Using the LES turbulence model, flow characteristics were analyzed based on 3D hill models, Quad-D176-H075 and Quad-D300-H075. The velocity contour of the first four modes was analyzed using the POD method. The results indicate that the acceleration effect at the hilltop may be related to the first-order mode.