Large Eddy Simulations of the Flow Fields over Simplified Hills with Different Roughness Conditions, Slopes, and Hill Shapes: A Systematical Study

: Turbulent ﬂow ﬁelds over topographies are important in the area of wind energy. The roughness, slope, and shape of a hill are important parameters a ﬀ ecting the ﬂow ﬁelds over topographies. However, these e ﬀ ects are always examined separately. The systematic investigations of these e ﬀ ects are limited, the coupling between these e ﬀ ects is still unrevealed, and the turbulence structures as a function of these e ﬀ ects are still unclear. Therefore, in the present study, the ﬂow ﬁelds over twelve simpliﬁed isolated hills with di ﬀ erent roughness conditions, slopes, and hill shapes are examined using large eddy simulations. The mean velocities, velocity ﬂuctuations, fractional speed-up ratios, and visualizations of the turbulent ﬂow ﬁelds are presented. It is found that as the hill slope increases, the roughness e ﬀ ects become weaker, and the roughness e ﬀ ects will further weaken as the hill changes from 3D to 2D. In addition, the fractional speed-up ratio at the summit of rough hills can even reach to three times as large as that over the corresponding smooth hills. Furthermore, the underestimation of the ratios of spanwise ﬂuctuation to the streamwise ﬂuctuation by International Electrotechnical Commission (IEC) 61400-1 is quite obvious when the hill shape is 3D. Finally, coherent turbulence structures can be identiﬁed for smooth hills, and as the hill slope increases, the coherent turbulence structures will experience clear evolutions. After introducing the ground roughness, the coherent turbulence structures break into small eddies.


Introduction
Flow structures over complex topographies are featured by the flow separations and reattachments which strongly depend on the topographic aspects, including surface roughness conditions, slopes, and shapes. On the other hand, the turbulent flow fields over complex terrain are of a great interest for many applications, such as wind turbine sittings [1], pollution diffusions [2], estimation of aerodynamic loadings on structures [3], identifications of tree damage [4], and forest fire propagation [5]. As a starting point modeling flow fields over real complex terrains, great efforts have been made to clarify the turbulent boundary layer (TBL) over simplified isolated hills. Two simplified isolated hills, i.e., two-dimensional (2D) and three-dimensional (3D) hills have been extensively examined in wind-tunnel experiments and numerical simulations.
For the roughness effects, Britter et al. experimentally found that at the lee side of the hill the roughness significantly alters the flow [6]. Then Pearse et al. examined the flow over several triangular and bell-shaped hills [7]. Importantly, it was found that increasing the surface roughness results in an increasing of amplification factors of the velocities at the hill crests. Using large eddy simulations (LES), Brown et al. concluded that the critical slope for the separation will increase as the ground roughness increases [8]. Moreover, the flow over a series of sinusoidal hills was investigated in wind tunnel experiments by Athanassiadou and Castro [9], who concluded that the effective roughness length is much greater than the length of the individual roughness blocks covering the hill. Importantly, it was pointed out by Finnigan and Belcher that even the canopy density is very limited [10], the canopy continues to play an important role in the flow above the canopy. Ross and Vosper found that the linear analytical solutions for the flow over a low and wide forested hill are still applicable [11]. From the studies by Cao and Tamura [12,13], it was determined that the speed-up ratio above the crest of a rough hill is larger than that of a smooth hill, same as the study by Pearse et al. [7], and it was also found that adding or removing ground roughness will create completely different turbulence structure in the wake, consistent with the study by Britter et al. [6]. Adopting LES, a constant mixing-length assumption was found not to be strictly valid within the canopy by Ross [14], the same conclusion reached in a recent study by Okaze et al. [15], and the structure of the turbulence over a forested hill was found to be broadly similar to that over flat ground, with sweeps and ejections dominating. Then, Loureiro et al. developed a consistent theory on the flow over rough ground using water-tank experiments [5,16]. Furthermore, the surface roughness was found to have a great influence on the flow separation point which occurs earlier with rougher surface leading to a larger recirculation area, in the studies by Takahashi et al. [17], Cao and Tamura [12], Tamura et al. [18], Loureiro et al. [5], and Cao et al. [19]. Most recently, Liu et al. concluded that with intent to capture the turbulent characteristics accurately the horizontal gird size should be at least as large as the height of the roughness canopy [20].
For the hill slope effects, Finardi et al. numerically found that the Cartesian coordinates and a brick-like terrain mesh are effective for modeling the flow over steep topographies [21]. Adopting a Reynolds stress model, Ying et al. concluded that a domain containing steeply-sloped topography requires higher horizontal resolutions for improving computational stability [3], which was also confirmed in the recent studies by Hu et al. and Ma and Liu [22,23]. Ferreira et al. found that the size of the recirculating region is strongly dependent on the hill slope and a quite evident growth of the recirculation bubble can be identified as the hill slope increases [24]. The wind tunnel experiments by Neff and Meroney and Athanassiadou and Castro as well as the numerical simulations by Griffiths and Middleton and Cao et al. showed that the flow separation occurs in the steep hill while the flow remains attached over the low-slope hill [9,19,25,26]. Most recently, the hills with different slopes were examined by Pirooz and Flay [27], who found that in comparison with the more peaked hill crests, the flat-topped hills would have a lower speed-up.
For the hill shape effects, Ishihara et al. experimentally found that the cavity zone behind a 3D hill is smaller than that behind the 2D hill attributing to the convergence of the flow in the wake of the 3D hill [28]. It was also found that the lateral velocity variances behind the 3D hill provide a secondary local maximum in the wall layer. This phenomenon was not observed in the 2D hill. Lubitz and White measured the flow past the hills with different shapes [29], and it was observed that owing to the secondary flow around the 3D hill, the speed-up and the size of the wake region of the flow over a 3D hill decrease significantly when compared with that over a 2D hill. Liu et al. adopted LES and found that the spectra of the fluids in the wake are sensitive to the oncoming turbulence condition for the 2D hill [30], which is not true for the 3D hill. Most recently, different vortices were identified numerically by Ishihara and Qi [31], in which the roller vortices were found to be significant on the lee side of the 2D hill, while horseshoe vortices appear around the 3D hill.
However, the roughness effects, the slope effects, and the effects of the shape of the isolated hills are always examined separately. The systematic examinations of these effects are limited, the coupling of these effects remains unrevealed, and turbulence structures as a function of these effects remain unclear.
In the present LES, twelve isolated hills with different ground roughness conditions, hill slopes and hill shapes were chosen. The flow fields were systematically examined, including the mean velocities, root mean square of the velocity fluctuations, speed-up ratios, fluctuation ratios; in addition, Energies 2019, 12, 3413 3 of 25 different turbulence structures in the wake were analyzed through the visualization of Q-criteria. The effects of the ground roughness, hill slopes and hill shapes, as well as the coupling between them are revealed.

Configurations
The numerical models are shown in Figure 1, where x, y, and z denote the streamwise, spanwise and vertical directions, respectively. The hills and the ridges are placed at (0, 0, 0). The shapes of the 3D hills are determined by: where, L is a constant equaling to 100 mm, and h = 20 mm, 40 mm, and 80 mm. The shapes of the ridges are determined by: where h = 20 mm, 40 mm, and 80 mm. Two vertical coordinates, z denoting the absolute height and z = z − z s (x, y) denoting the height above the local surface, will be used for the convenience of discussion. Both the smooth and rough ground conditions are considered. For all of the rough cases, the height of the roughness canopy equals 5 mm. Varying the hill heights but constant radius is from the consideration that in a certain district, the higher the hill is, the steeper the hill is more likely to be in the real situation. And the constant height of roughness canopy is from the consideration that the height of the forest is nearly independent with the hill height in the real situation. The case settings are listed in Table 1 and the shape of the topographies are shown in Figure 2. Four of the cases are modelled to be in accordance with the Ishihara et al. [28] experiments for validation purposes. The wind tunnel experiments by Ishihara et al. [28] did not consider the effects of the hill slopes, which is one of the motivations of the present study.  To best reproduce the experimental results, the configuration for our numerical model is set same as that in the wind tunnel experiments by Ishihara et al. [29], except the domain size in the spanwise direction, Ly and the upstream necking zone. In the experiment of a smooth 3D hill by Ishihara et al. [29], turbulent boundary layer (TBL) was simulated using two 60 mm high cubic arrays placed downstream of the contraction exit, followed by 20 mm and 10 mm cubic roughness elements, covering 1.2 m of the test-section floor. The remaining 5.8 m of the test section floor was covered with plywood, which was as smooth as the hill surface. The 3D hills or 2D ridges were mounted 4.6 m downstream of the contraction exit. In the simulation, the domain extends over (Lx, Ly, Lz) = (9, 0.65, 0.9) m 3 = (25, 1.8, 2.5) δ 3 = (90, 6.5, 9) L 3 . Two nested domains (coarse and fine) are adopted, as illustrated by the red dashed lines in Figure 1. The fine-grid domain covers the range of (Lx′, Ly′, Lz′) = (10, 2, 3) L 3 in the x-, y-, and z-directions, respectively. Both the upstream and the downstream fine-grid regions are 5L long. To best reproduce the experimental results, the configuration for our numerical model is set same as that in the wind tunnel experiments by Ishihara et al. [28], except the domain size in the spanwise direction, L y and the upstream necking zone. In the experiment of a smooth 3D hill by Ishihara et al. [28], turbulent boundary layer (TBL) was simulated using two 60 mm high cubic arrays placed downstream of the contraction exit, followed by 20 mm and 10 mm cubic roughness elements, covering 1.2 m of the test-section floor. The remaining 5.8 m of the test section floor was covered with plywood, which was as smooth as the hill surface. The 3D hills or 2D ridges were mounted 4.6 m downstream of the contraction exit. In the simulation, the domain extends over (L x , L y , L z ) = (9, 0.65, 0.9) m 3 = (25, 1.8, 2.5) δ 3 = (90, 6.5, 9) L 3 . Two nested domains (coarse and fine) are adopted, as illustrated by the red dashed lines in Figure 1. The fine-grid domain covers the range of (L x , L y , L z ) = (10, 2, 3) L 3 in the x-, y-, and zdirections, respectively. Both the upstream and the downstream fine-grid regions are 5L long.

Governing Equations
In the LES strategy, large eddies are explicitly resolved, while the small eddies are parameterized by SGS models. The governing equations are usually obtained by filtering the time-dependent Navier-Stokes equations in Cartesian coordinates (x, y, z):

Governing Equations
In the LES strategy, large eddies are explicitly resolved, while the small eddies are parameterized by SGS models. The governing equations are usually obtained by filtering the time-dependent Navier-Stokes equations in Cartesian coordinates (x, y, z): whereũ i andp are the filtered velocity and pressure, respectively; µ is the viscosity; ρ is the density; and τ ij is the SGS stress. To close the equations for the filtered velocities, a model for the anisotropic residual stress tensor τ ij is needed, which is obtained as: where µ t denotes the SGS turbulent viscosity;S ij refers to the rate-of-strain tensor for the resolved scale, and δ ij is the Kronecker delta. The Smagorinsky-Lilly model is used to parameterize the SGS turbulent viscosity as: where L s stands for the SGS mixing length; κ represents the von Kármán constant (= 0.42); d is the distance to the closest wall, and Λ is the volume of a computational cell. Here, C s denoting the Smagorinsky constant is set to a value of 0.1 as studied by Liu et al. [20]. Since z + is below 2 in all of the cases, the shear stresses are obtained from the viscous stress-strain relation: where u * is the friction velocity and δ n is the distance between the cell centre and the wall.

Method Modeling Roughness
Accurate modeling of the inflow is essential for the success of simulating the flow in the 3D hill wake. In the present study, the arrangement of the upstream roughness blocks in the simulations is exactly same as that in the wind-tunnel experiments by Ishihara et al. [32]. In addition, the distance between the roughness blocks and the location of 3D hills placed is also set exactly same as that in the wind-tunnel experiment. In the volumes occupied by the roughness blocks, the drag force term is added: in which the drag force term f u,i = C f A f ρũ i |ũ i |; since the velocity in the volume occupied by the roughness blocks should be zero, the drag coefficient C f is set as 100 to model the drag effect from the solid roughness blocks; A f is the frontal area of the roughness block. In order to model the canopy, the drag force term, f u,i is added to the momentum equations in the same form as Equation (10). The drag force term is determined by is the leaf-area density, and |ũ i | is the velocity magnitude. C f and A f are determined following the study by Liu et al. [20].

Locations Boundary Type Expression
Outlet of the domain Outflow In terms of the grid size, in the vertical direction, the grid is stretched starting from a vertical grid spacing of 0.1 mm at the surface in both the fine-and coarse-grid domains in the smooth cases. When the ground is covered by vegetation canopy, 10 grids are adopted to divide the vertical space in the canopy. The resulted vertical grid resolution at the surface z + < 1.0 covers most of the domain except for the windward part of the hill where the maximum z + does not exceed 2. A horizontal grid size of 2.0 mm is used in the fine-grid region. In the rough-grid region, the horizontal grid shape is square, and a uniform grid size of 10 mm is applied. Downstream of the fine-grid domain and upstream of the roughness blocks, ∆x is stretched at a ratio of 1.2. Total grid number is 2.2 × 10 7 for the smooth cases and 2.45 × 10 7 for the rough cases. The parameters determining the grid system are listed in Table 3.

Solution Schemes
Finite volume method (FVM) is used in the present LES. The second-order central difference scheme is used for the convective and viscous terms, and the second-order implicit scheme is employed for the unsteady term [33]. Time-step size ∆t is set as 0.0001 s, and in convective time units, ∆t * = ∆tU ∞ /L = 0.0058. The Courant Friedrichs Lewy (CFL) number (Courant et al. 1928) is based on the time step size (∆t), velocity (ũ i ), and grid size (∆x i ), expressed as C = ∆tΣũ i /∆x i . Here, the CFL number is limited to be less than 2 (i.e., C max = 2) in the whole computational domain. The SIMPLE (semi-implicit pressure linked equations) algorithm, which was introduced by Ferziger & Peric [33], is applied to solve the discretized equations. The calculations are performed on 2PCs in parallel (Intel core i9-7980XE, 18 cores, 64G memory, Dell, Beijing, China). The simulation cost 2521 hours in total and the initial transient effects are found to disappear after 400 time units. Statistical convergence for the mean velocities is achieved when | ũ i−y − ũ i+y /[( ũ i−y + ũ i+y )]/2 | < 1% in the near-wake of the 3D hill, which is over 350 time units, where means the time-averaging process, −y and +y mean the velocity at two points which are symmetrical position in lateral direction. For the fluctuations, u, v, and w, 350 time unit statistics still cannot show clear convergence, , especially in the region z < 1.0L for the streamwise component, which should be due to the much large turbulence in this region. However, the simulations of these 12 cases have cost us over 5 months based on the computational resources available in our group. Further increasing the statistical time can hardly be afforded by us. In the future, more detailed examination should be carried out using larger statistical time. Table 4 summarizes the numerical schemes in the present LES.

Upcoming TBL
The resulted profiles of mean streamwise velocity, U, root mean square (r.m.s.) of streamwise velocity fluctuation, u, spanwise velocity fluctuation, v, and vertical velocity fluctuation, w, in the absence of the hills, are compared with those in the experiment by Ishihara et al. [28], as shown in Figure 3. U ref determined as the mean streamwise velocity at the height of 0.16 m is applied to normalize the flow fields, which equals to 5.3 m·s −1 in the present LES and 5.2 m·s −1 in the wind-tunnel experiment. For the mean velocity profiles, there are some differences between the experimental data and those from the simulations. The difference is mainly below z = 0.1 and L = 10 mm, where is very close to the ground. The introduction of the probe for measurement in the experiment will disturb the flow. That should be the reason of the clear differences between the experiment and the simulation. The difference between the simulation and the experiment above 0.1L is below 5% which is in the acceptable range of engineering applications. For the velocity fluctuations, the components in spanwise and the vertical directions from the simulations are comparable with those from the wind tunnel experiments. The maximum discrepancy is below 5%. However, for the streamwise component, relatively large discrepancies can be found in the region 0.6L < z < 1.0L and 0.1L < z < 0.4L, and the maximum discrepancy reaches 10% at z = 0.2L. These discrepancies may result from the insufficient statistical time in the present simulations. In the wind tunnel experiments, the statistical time is10 s, which is over 2000 time units which is about 5 times as large as that in the simulation.
The boundary layer thickness in absence of the hills, δ, was about 0.36 m (both smooth and rough cases). The scale of the simulated boundary layer, λ, was about 1:1000, on the basis of the power spectra of the longitudinal velocity component. The wind speed outside the boundary layer was measured as U ∞ = 5.8 m·s −1 . As a result, the simulated boundary layer had a bulk Reynolds number: which equals 1.4 × 10 5 , where υ was the kinematic viscosity of the air equaling 14.8 × 10 -6 m 2 ·s −1 . The roughness height z 0 was 0.01 mm in the smooth cases and 0.3 mm in the rough cases. z 0 was determined by comparing the resulted mean streamwise velocity profile with the logarithmic law U (z) = U τ κ ln z z 0 , where κ = 0.4 was the Von-Kármán constant, U τ = τ/ρ = 0.212 m·s −1 , τ denotes the time averaged wall shear stress of the streamwise component. The main TBL parameters are summarized in Table 5.

Numerical Results
In the following figures about the mean velocities profiles (Figures 4,5) and the fluctuations ( Figures 6-8), the separation boundary defined as the connection of the reversion points on U profile is illustrated by solid red lines for the smooth hills and dashed red lines for the rough hills from LES. The corresponding pink lines are from the experiment by Ishihara et al. [29]. The shear layer center determined by the connection of the peak r.m.s. streamwise fluctuations, u, is illustrated by solid yellow lines for the smooth hills and dashed yellow lines for the rough hills from LES. The corresponding brown ones are from the experiment by Ishihara et al. [29]. The discussions about the results are mainly from five aspects, i.e., the roughness effects, the hill slope effects, the hill shape effects, the coupling between these effects, and the LES performance.

Numerical Results
In the following figures about the mean velocities profiles (Figures 4 and 5) and the fluctuations (Figures 6-8), the separation boundary defined as the connection of the reversion points on U profile is illustrated by solid red lines for the smooth hills and dashed red lines for the rough hills from LES. The corresponding pink lines are from the experiment by Ishihara et al. [28]. The shear layer center determined by the connection of the peak r.m.s. streamwise fluctuations, u, is illustrated by solid yellow lines for the smooth hills and dashed yellow lines for the rough hills from LES. The corresponding brown ones are from the experiment by Ishihara et al. [28]. The discussions about the results are mainly from five aspects, i.e., the roughness effects, the hill slope effects, the hill shape effects, the coupling between these effects, and the LES performance.

Mean Streamwise Velocity
As shown in Figure 4, the mean streamwise velocity, U, is found to decrease to nearly zero in the roughness region, but accelerate rapidly from the roughness top at the summit. It can be also identified that the vertical distance of the reversion point with reference to the local ground becomes higher as introducing the ground roughness, implying a larger recirculation bubble. This should be attributed to the fact that the ground roughness can produce more turbulence in the wake and in return prevent the flow from reattaching. However, this trend is not preserved for 80 mm high 2D hills, see Figure 4f, in which the recirculation bubble gets even smaller after introducing the ground roughness. It should result from a balance between two effects of the roughness, i.e., preventing the flow from reattaching, and decelerating the negative U in the bubble. The former effect may enlarge the bubble while the later one may shrink the bubble. For 80 mm high 2D hills, the deceleration effects of the roughness may be stronger than those preventing the flow from reattaching.
As the hill slope increases, the flow in the wake shows larger deceleration, while the acceleration of U in the vertical direction at the summit increases. In addition, at the lee-side, the recirculation bubble is significantly stretched and the speed-up of U above the recirculation bubble seems to extend further downstream for steeper hills, which is also an indication of a more stable recirculation bubble.
U seems insensitive to the change of the hill shape from 3D to 2D at the windward side and the summit. However, downstream the summit, the hill shape effect gets significant. Firstly, it is obvious that the blockage effects of 2D hills are more energetic than those of 3D hills, which is the result that the flow approaching the 3D hills can move around the lateral sides, whereas there is no such pass for the flow over 2D hills. Secondly, the near-ground acceleration of negative U is more obvious as the hill shape turns to 2D.
For the coupling between these effects, with increasing the slope, the roughness effects on U at the locations above z = h becomes weaker. In addition, for the hills with low slope, the change of the hill shape from 3D to 2D does not show significant evolution of the flow fields (Figure 4a,d); however, as the hill slope increases, the hill shape effects turn evident (Figure 4c,f). Furthermore, for the hills with low and moderate slopes (Figure 4a,d, and Figure 4b,e), the roughness effects and the hill shape effects seem to be independent, whereas if the hill is steeply-sloped (Figure 4c,f), the roughness effects will be weakened as the hill alters from 3D to 2D.
For the performance of LES, importantly, it is obvious that the rougher the ground, the worse LES results will be. And the present LES shows larger discrepancies for 2D hills as compared with 3D hills. This should be due to the smaller eddies caused by the roughness or the stronger disturbance from 2D hills, which requires finer space and time resolutions to reproduce the flow structures. The discrepancies concentrate at the lee-side of the hills, while the separation points and the reattachment points are well predicted. The overall comparison between the numerical and the experimental results is fairly good. Figure 5 shows that the introduction of ground roughness will weaken the acceleration of the mean vertical velocity, W, at the windward side. However, at the summit, the maximum W over the rough hills becomes even larger than that over the smooth hills, which is more obvious for 3D hills. And this maximum W occurs just at the roughness top. It is also clear that after introducing ground roughness, the recovery of W in the wake turns quicker, which should be the result of the stronger mixing effects from the roughness-generated turbulence. In addition, above the recirculation bubble, the negative W is weakened when the ground is rough.

Mean Vertical Velocity
As the hill is more steeply-sloped, the near-ground W acceleration at the windward side is obviously increased, mainly because the wind at the windward side almost flows following the shape of the hill. Importantly, the steeper the hill is, the sharper the near-ground W profiles at the windward side will be, which is the most obvious as the flow moves to the summit. In the wake region, the trend of stronger negative W with increasing the slope is evident due to the more energetic recirculation bubble. Importantly, at some locations, such as x = 1.0L, the downward W even almost equals the peak upward W at the summit for both 3D and 2D hills with low and moderate slopes (Figure 5a,b,d,e), which is however not observed for the steepest 2D hills ( Figure 5). even almost equals the peak upward W at the summit for both 3D and 2D hills with low and moderate slopes (Figure 5 a,b,d,e), which is however not observed for the steepest 2D hills ( Figure 5).      When the hill shape alters to 3D, W acceleration is found to extend to lower elevation at the windward side; however, W profiles become sharper. It is interesting that the downward W in the wake vanishes as the hill shape changes to 2D and the hill turns steeper (Figure 5f). This should result from the strong turbulence in the wake of the steep 2D hills.
As for the coupling between these effects, the distortion of W profiles owing to the ground roughness seems to be strengthened for the steeper hills, while it seems to be weakened as the hill changes from 3D to 2D. Importantly, at the summit of the hills, the acceleration of W gets quicker for 2D hills as the slope is increasing.
Similar to U, the discrepancies of W with the experimental data is concentrated in the near wake. However, different from U discrepancies, the large errors of W are found in the flow over the smooth hills. For 40 mm smooth 3D hill, the sharp near-wake deceleration of W cannot be predicted well and the negative W only reaches to about half of that measured in the experiment. Except this location, the predicted W shows satisfactory agreements with the experiment.

Fluctuations of Streamwise Velocity
It is evident that the strea·wise fluctuation, u, at x = −L very close to the ground for the rough hills are smaller than the corresponding smooth hills, as seen from Figure 6, which should be the result of the limitation of the eddy development due to the drag effects in the roughness canopy. However, the flow can enter the roughness upper region and generate wavy structures near the roughness top, thus creating more turbulence above the roughness. And interestingly, in the wake region the elevation of the shear layer is higher for the rough hills, but at the shear layer, u over the smooth hills is larger.
As the slope increases, at the summit, u seems to be restrained, which is more obvious when the hill shape is 3D or the ground is rough. Furthermore, the peak u in the wake is found to be sharper when the hill is more steeply-sloped. It is also important that the wake depth grows quickly as the flow moves downstream; however, the increasing speed of the wake depth is slowed for the steeper hills. For the 80 mm smooth 3D hill (Figure 6c), the wake depth is even surprisingly decreased as the flow moves downstream.
After changing the hill shape to 2D, the location of the peak u is found to be lower than the corresponding 3D hills. However, this trend becomes not obvious for the steepest hills (Figure 6c,f). In addition, the peak u is obviously enhanced when the hill shape changes to 2D, which should attribute to the more stable recirculation bubble in the wake of 2D hills, providing stronger wind shears. This trend is even more obvious for the smooth hills.
Except the findings about the coupling of the ground roughness effects, hill slope effects, and the hill shape effects in the above presentations about u, we can also observe that in general the difference between the height of peak u of the smooth and rough hills gets smaller as the hill slope increases, and this trend becomes quicker as changing the hill shape from 3D to 2D, which may imply that in the real situation the forest may play a limited role in affecting the flow fields if the hill shape is about a 2D ridge and the hill is steep enough.
In LES, at the upstream footage of the hills, there is a large increase of u, implying the formation of a horseshoe vortex, but this feature is not observed in the experiments by Ishihara et al. [28]. This may be attributed to the difficulty setting the probes at this location or the fact that the existence of the probes weakens the horseshoe vortex. Interestingly, u at the hill crest shows a sharp peak in the present LES. However, this sharp peak is not measured in the experiments by Ishihara et al. [28], which should result from the limited observation points in the experiments or the disturbance of the flow fields very close the ground from the probes. Additionally, comparing with the experimental data, in the near-wake region the LES shows obvious overestimations for the 40 mm smooth 3D hill, and obvious underestimation for the 40 mm rough 2D hill, whereas these discrepancies decrease to the acceptable level as soon as the flow moves to the locations x > 1L.   (a)

Fluctuations of Spanwise Velocity
For the spanwise fluctuation, v, shown in Figure 7, the similar trend as u can be observed. However, there are some unique features for v. Firstly, in the wakes, distinct v profiles can be identified for the 3D hills (Figure 7a-c). It is obvious that no matter how steep the 3D hill is, two peaks appear on the vertical v profiles with one locating near the ground and the other locating in the region close to the shear layer. The LES by Liu et al. has revealed the mechanism of the two peaks on v profiles [34], in which the secondary turbulence structures surrounding the major vortex core in the wake are believed to be the source. Two peaks of v still appear but get less obvious when the ground is rough. In addition, as increasing the hill slope, these two peaks become more obvious, indicating the more energetic secondary vortex.
Secondly, for 2D hills with h = 80 mm, it is surprising that v in the recirculation bubble is nearly a constant at each streamwise location, implying the full mixture from the hills and abruptly homogeneous turbulence, which is confirmed by plotting the instantaneous flow fields using Q-criteria in Figure 12. The performance of the present LES for predicting v is similar as u. Some errors can still be found in the near-wake region. Figure 8 shows the distributions of vertical fluctuations, w, which are found to be very similar to those of u and v. Therefore, no detailed discussion about w is provided. In general, w is smaller than u and v, and there is no additional peak very close to the ground.

Fractional Speed-Up Ratio
Fractional speed-up ratio, ∆S, has been extensively applied to assess the impact of topography. ∆S is defined as: where U(x, z ) is the mean streamwise velocity at z , and U 0 (z ) is the reference streamwise velocity at the same height in the absence of the topography. The horizontal distributions of ∆S at z = 10 m and 50 m are illustrated in Figure 9. The selection of these two heights takes into account that 10 m is the reference height adopted in most of the wind-resistant guidelines, and 50 m is about the hub height of 1 MW or 1.5 MW wind turbines in operation in many hilly regions of China. It can be clearly seen from Figure 9 that the introduction of ground roughness will obviously increase ∆S at z = 10 m, where ∆S at the summit of rough hills can even reach to three times as large as that over the corresponding smooth hills. However, ∆S near the summit at z = 50 m is not sensitive to the ground roughness condition. On the other hand, the deceleration of ∆S in the wake becomes stronger in general if the ground is rough even at high elevations, which should attribute to both the drag effects and the more energetic turbulence caused by the roughness.
With increasing the hill slope, ∆S at the summit increases obviously. ∆S (z = 10 m) at the summit over rough hills reaches 0.4 for the low-slope hills, 0.7 for the moderate-slope hills, and 0.9 for the large-slope hills. However, it is worthwhile to point out that the area with positive ∆S turns smaller as the hill gets steeper. Furthermore, as increasing the hill slope, the deceleration area in the hill wakes will also expand.
It is also interesting that the hill shape shows unapparent effects on the maximum ∆S at the summit, the area of positive ∆S near the hill top, and the minimum ∆S in the near wake. However, obviously 2D hills will provide larger area with negative ∆S. The experimental data of ∆S are superimposed on Figure 9, where satisfactory agreement is achieved. The largest discrepancies are found in the wake of the 40 mm high 2D rough hill, which is about −0.3 in the LES while nearly −0.1 in the experiment.

Fluctuation Ratios
The fluctuation ratios are important parameters for determining the turbulent flow fields. In the present LES of the flow over flat ground, v/u and w/u are predicted as 0.81 and 0.49 respectively,

Fluctuation Ratios
The fluctuation ratios are important parameters for determining the turbulent flow fields. In the present LES of the flow over flat ground, v/u and w/u are predicted as 0.81 and 0.49 respectively, close to the data suggested in IEC 61400-1 [35]. In case a wind turbine site is located within a complex terrain, IEC 61400-1 permits to increase the representative value by a factor C CT defined as [35]: C CT is intended to account for the distortion of turbulence structure by complex terrains and to be estimated based on the site specific data. In the absence of site specific data, IEC 61400-1 recommends the use of a correction factor of 1.15 [35]. When C CT = 1.15 and if the fluctuation ratios v/u and w/u are modified at the same rate, v/u and w/u become 1.0 and 0.7, respectively. The horizontal distributions of v/u in LES are shown in Figure 10. The symbols superimposed in Figure 10 are the data from wind-tunnel experiment by Ishihara et al. [28], and the dash-dotted lines indicate the data from IEC 61400-1 [35]. Same as the discussion about ∆S, the data on z = 10 m and z = 50 m are extracted.
Obviously, the introduction of ground roughness can hardly affect v/u, and the underestimation of v/u by IEC 61400-1 in the hill wake is quite obvious when the hill shape is 3D [35]. The results from LES can reach about 1.5 times as large as those in the guideline at z = 10 m. Considering the fact that the suggested values in the guideline is based on the assumption of isotropic turbulence, the larger v/u in LES for 3D hills in the wake may indicate that the turbulence structure should be coherent instead of isotropic. From the snapshots of the instantaneous flow fields shown in Figure 12, the coherent turbulence structure over 3D hills is characterized mainly by a spanwise sway motion of the fluid, inducing a large spanwise fluctuation. However, this spanwise sway motion can hardly change the streamwise size of the eddies, therefore v/u becomes larger than that in the fluids characterized by isotropic eddies. We can also find that after changing the hill shape from 3D to 2D, the stronger blocking effects induce more energetic mixture, yielding nearly equal values of v/u as those in the guideline. However, it is important that at the summits, where most of the wind turbines are located, the data from the guideline are nearly consistent with the LES for both 2D and 3D hills.
As the hill slope increases, the underestimation from the guidelines seems to be enhanced, which is more obvious when the hill shape is 2D. For 80 mm 2D hills, the large underestimation from the guideline can be observed over almost the entire hill area except the locations near the summit. Furthermore, the difference between v/u over smooth hill and that over rough hill is enlarged as increasing the hill slope. At the high elevation, z = 50 m, the predicted v/u is abruptly equal to that from IEC 61400-1 [35], which is also the indication of the different turbulence structures of the flow close the ground and that at high elevations.
The horizontal distributions of w/u are shown in Figure 11, which is found to be similar to the distributions of v/u. However, at the windward side, the large values found in v/u do not appear in w/u. And, different from v/u, w/u shows larger values in the far wake region for the steepest hills.  R8S and R8R

Instantaneous Flow Fields
The Q-criterion is further adopted to shed some lights on the instantaneous turbulent structures, which is defined as: where S ij and Ω ij are the antisymmetric and symmetric components of the velocity-gradient tensor, respectively, quantifying the relative amplitude of the rotation rate as well as the strain rate of the flow, expressed as: The snapshots of instantaneous Q with values of 10000 (purple) and −10000 (green) are shown in Figure 12, where the thin yellow lines indicate the boundary of the hills, the yellow thick dashed lines show the major core of the coherent structure in the wakes, and the white dashed lines indicate the locations of the secondary vortex surrounding the major core.
For smooth 3D hills, as increasing the hill slope the major core of the coherent turbulence structure in the hill wake tends to show large spanwise sway motions, shown in Figure 12e. It is obvious that after introducing the ground roughness, the clear coherent turbulence structures are broken into small eddies. However, for the rough 3D hills, when the hill slope is large enough, the periodical vortex shedding will appear again, see Figure 12k.
Different evolution process of the turbulence structure in the hill wakes can be identified for 2D hills. When the hill is smooth and the hill slope is low, the major core of the wake vortex is nearly perpendicular to the ridge line. As the hill slope increases (Figure 12d), a kind of ejection-sweep structure of large scale occurs. Further increasing the hill slope, the ejection-sweep structure turns to be a wavy structure with the major core being parallel with the ridge line (Figure 12f). In addition, when the 2D hills are covered with roughness canopy, no clear coherent structure can be identified.

Conclusions
LES are adopted to study the flow fields over simplified topographies with different shapes, hill heights and rough conditions. The findings of this study are summarized below:

1)
For the mean streamwise velocity, as the hill slope increases, the roughness effects get weaker; however, the hill shape effects become more evident. And if the hill is very steeply-sloped, the roughness effects will be further weakened as the hill changes from 3D to 2D. 2) For the mean vertical velocity, the distortion of the profiles owing to the ground roughness is strengthened for the steeper hills, while it seems to be weakened as the hill alters from 3D to 2D. Importantly, at the hill summits, the acceleration of W as increasing the slope becomes quicker for 2D hills. 3) For the fluctuations, at the summit, u seems to be restrained as increasing the slope, which is more obvious when the hill shape is 3D or the ground is rough. After changing the hill shape to 2D, the location of the peak u is lower than the corresponding 3D hills. Two peaks on v profiles still appear but become less obvious when the ground is rough. In addition, as increasing the hill slope, these two peaks turn more obvious. 4) For the fractional speed-up ratio, ∆S at the summit of rough hills can reach to three times as large as that over the corresponding smooth hills. With increasing the hill slope, ∆S at the summit increases obviously. However, the hill shape shows unapparent effects on the maximum ∆S at the summit. In addition, the area with positive ∆S becomes smaller as the hill gets steeper. 5) For the fluctuation ratios, the introduction of ground roughness can hardly affect v/u, and the underestimation of v/u by IEC 61400-1 in the hill wake is quite obvious when the hill shape is 3D [35]. After changing the hill shape from 3D to 2D, v/u becomes almost equal as those in the guideline. 6) For the turbulence structures, clear coherent turbulence structure can be identified for smooth 3D hills, and as increasing the hills slope, the major core of the coherent structure tends to show large spanwise sway motions. For 2D hills, when the ground is smooth and the slope is low, the major core of the wake vortex is nearly perpendicular to the ridge line, and as increasing the hill slope, a kind of ejection-sweep structure of large scale occurs. Further increasing the hill slope, a wavy structure with the major core being parallel with the ridge line appears. After introducing the ground roughness, the clear coherent turbulence structures are broken into small eddies.
Energies 2019, 12, x FOR PEER REVIEW 23 of 25 3D [1]. After changing the hill shape from 3D to 2D, v/u becomes almost equal as those in the guideline.
6) For the turbulence structures, clear coherent turbulence structure can be identified for smooth 3D hills, and as increasing the hills slope, the major core of the coherent structure tends to show large spanwise sway motions. For 2D hills, when the ground is smooth and the slope is low, the major core of the wake vortex is nearly perpendicular to the ridge line, and as increasing the hill slope, a kind of ejection-sweep structure of large scale occurs. Further increasing the hill slope, a wavy structure with the major core being parallel with the ridge line appears. After introducing the ground roughness, the clear coherent turbulence structures are broken into small eddies.