Comparison of CFD Simulation to UAS Measurements for Wind Flows in Complex Terrain: Application to the WINSENT Test Site

Abstract

The two modified transport equations for the turbulent kinetic energy k and the dissipation ε read: ∂(ρhk)∂t+∂(ρh ujk)∂xj=∂∂xj[(μ+μtσk)(∂k∂xj) ]+P+G+Sk−ρhε ∂(ρhε)∂t+∂(ρh ujε)∂xj=∂∂xj[(μ+μtσε)(∂ε∂xj) ]+Cε1*(P+G)+Sε−Cε2ε2k where P represents the production rate of turbulent kinetic energy due to shear and G represents the production/destruction of turbulence by buoyancy forces defined as: P=τij(∂ui∂xj)=μt(∂ui∂xj+∂uj∂xi )∂ui∂xj−23ρh ∂uk∂xkδij G=−1Tgμtσh∂θ∂xi Sk and Sε are source terms relative to the forest canopy and described later in Equations (15) and (16). σk, σε , Cε1* , Cε2 are model coefficients listed in Table 1. The thermal stratification has a large impact on the vertical wind profile and turbulence levels. [...]different thermal stratification cases will be conducted in future studies when the complete measurement set up will be running. See PDF.] Turbulence Model Constants Cµ Cε1 Cε2 σε σk Standard [29] 0.09 1.44 1.92 1.00 1.3 Adapted [26] 0.256 1.13 1.90 0.74 1.3 Closure Constants βp βd Cε4 Cε5 Value 1 5.1 0.9 0.9 Variable Level Variables Slope α (°) Height H (m) LAI −1 15 20 2 +1 30 25 5 Mesh Properties Number of Cells Maximum Cell Size Horizontal Resolution at the Ground Vertical Resolution at the Ground Basis 66.8 million 100m 25m 2.6m Nested 25.4 million 50m 6.75m 1.2m Run LAI Forest Height H (m) Run 1 (no canopy) 0 0 Run 2 2 20 Run 3 2 25 Run 4 5 20 Run 5 5 25 Author Contributions All the authors have substantial contribution to this work. Rautenberg, A.; Allgeier, J.; Jung, S.; Bange, J. Calibration Procedure and Accuracy of Wind and Turbulence Measurements with Five-Hole Probes on Fixed-Wing Unmanned Aircraft in the Atmospheric Boundary Layer and Wind Turbine Wakes.

1. Introduction

Wind power is currently one of the most promising renewable energy sources. The year 2017 was a record year for annual installations in Europe, with 16.8 GW of additional wind power capacity installed. Wind energy remains the second largest form of power generation capacity in Europe, closely approaching gas installations. In the EU, wind energy overtook nuclear energy in 2013, hydro in 2015, and coal in 2016. In 2017, offshore installed wind power capacity represented 15.8 GW against 153 GW for onshore installations [1]. Onshore installations are mainly built on flat terrain, making them easier to operate compared to those mounted on hilly terrain, where forecasts are more uncertain, wear and tear is greater, and maintenance and construction costs are higher. However, wind energy in mountainous regions has been making inroads in recent years and is of increasing interest to the wind-energy community.

WindForS, a wind energy research cluster in Southern Germany, aims to answer the question of how to optimize installations in complex terrains and extend their service life. In the framework of the project Wind Science and Engineering in Complex Terrain (WINSENT), a field-test site for research and industry, located behind an escarpment in the Swabian Alps, near the town of Geislingen an der Steige is setting up. The project working plan is divided into two stages: the first phase investigates the local wind flow without wind turbines where different measurement equipment is used to characterize the wind flow, such as towers equipped with anemometers, Lidar [2], eddy-covariance stations, and an unmanned aircraft system (UAS) (see [3,4,5]). The second stage will include two wind turbines with a nominal output of around 750 kW. The turbine-terrain interference will be studied in this second phase.

Characterizing the wind flow in complex terrain is more challenging compared to flat terrain. At these locations, wind flows are more complex, and they are influenced by changes in topography such as hills, escarpments, and roughness, leading to non-linear features such as high levels of turbulence, wind shear, unsteadiness, etc. Analytical models, mainly developed in the 1970s and 1980s such as the well-known Wind Atlas Analysis and Application Program (WAsP), are capable of predicting the mean wind field only in simple geometries and are suitable only for attached flow (see [6,7,8,9] and others). Modelling non-linear features of the flow became a feasible option with the help of computational fluid dynamics (CFD) solvers which retain the nonlinearity of the Navier-Stokes equations and simulate momentum, turbulence, and energy [10]. The commonly referenced field experiment of Bolund hill, which consists of a steep hill, has been used for a blind comparison of models with different ranges of fidelity (linearized and CFD models) as reported in [11]. One of the outcomes of this comparison was that the linearized models were not able to predict the mean flow features such as the speed-up, unlike the CFD models (see [12] for a complete review). In [13], the flow field over a large-scale model of the same hill was investigated experimentally and showed that the mean wind, wind shear, and turbulence level are extremely sensitive to the details of the terrain. Making the edge of the hill sharper resulted in a reduction of the estimated annual energy production by at least 50% and an increase in the turbulent level by a factor of five in the worst-case scenario. They concluded “the mean wind, wind shear, and turbulence level are extremely sensitive to the exact details of the terrain”. This shows that special attention should be paid not only to the model but also to the topography representation in complex terrains by using high-resolution topographic data and information on land use, such as the Corine land cover [14].

The use of large eddy simulation (LES) or detached eddy simulation (DES) in complex terrain has increased in recent years. These models are superior to both steady and unsteady Reynolds averaged Navier-Stokes (RANS) as they have the advantage of providing additional information on turbulent structures and non-linear features of wind flow over complex terrain. However, the requirement of high computational resources and the challenge of obtaining proper inflow boundary conditions limits their usefulness for wind engineering studies where a fast solution is required [10]. In addition, for wind-energy applications, the simulation often needs to be computed for different wind directions or stability conditions. Due to these constraints, RANS/URANS models are still appropriate for use as they provide a good balance between computational effort and model accuracy.

A common problem when using CFD models in a complex terrain is how to specify initial and boundary conditions. A standard way to proceed is to impose a standard logarithmic velocity profile at an inlet which is orientated perpendicular to the wind direction. Another approach is to use numerical weather prediction (NWP) models, also called mesoscale models, to provide more realistic boundary conditions for the CFD simulations. Mesoscale models, generally, have a low spatial resolution with a horizontal extent on the order of 2 km and they generally present a bias on the predicted wind speed and turbulent quantities in complex terrains due to unresolved topographic effects. In [15], the WRF mesoscale model for the Horns Rev wind farm was applied with a horizontal resolution of 333 m and showed that the model still underestimates the power deficit due to its coarse resolution. Nevertheless, as demonstrated in [16], much effort has been made to improve the NWP models. Coupling a mesoscale model to a microscale model (CFD), which presents a more detailed representation of the topography seems to be one of the more promising approaches for wind-energy assessment in complex terrain. Several coupling methods of mesoscale and microscale models have been developed in the last decade, and an overview of these methods can be found in [10,17]. One of the most common coupling methods is the one-way approach where the mesoscale model is coupled to the microscale model through the lateral boundaries at fixed times (time-slice). This approach has been successfully applied in complex terrain by several authors as [18,19,20] among others.

In the present paper, the WINSENT test site is studied by means of numerical simulations along with UAS data. We confine the study to a neutral stratification case. Section 2 presents the physical model with its validation, while Section 3 presents the parametric study. The parametric study is conducted with the help of a design of experiment (DoE) method applied to a two-dimensional case. This study intends to assess the sensitivity of our model results for different parameters and predict their effect on some relevant variables for wind turbines. Section 4 gives an overview of the test site and the measurement system, while Section 5 presents the model which uses the one-way coupling approach. The same section compares simulation results against UAS measurements. Section 6 draws a summary and discusses the limitation of the model and possible future improvements.

2. Methods 2.1. Numerical Model

To simulate the atmospheric boundary layer (ABL) flow, the finite volume method (FVM) on the OpenFOAM v2.4.0 (Open Source Field Operation and Manipulation) software, provided by the OpenFOAM Foundation U.K., was used [21]. The OpenFOAM toolbox includes open source C++ libraries released under the general public license (GPL). An unsteady Reynolds Averaged Navier-Stokes approach under the Boussinesq approximation, where density is only influenced by buoyancy forces, was considered [22]. The transport equations for mass, momentum, and potential temperature can be written relative to a hydrostatic state (index h ) in a Cartesian coordinate as:

∂(_ρh)∂t+∂(_{ρh uj})∂_xj=0

∂(_{ρh ui})∂t+∂(_{ρh ui uj})∂_xi=−∂p∂_xi+∂∂_xj[(μ+_μt)(∂_ui∂_xj+∂_uj∂_xi)−23k_δij]−_ρh(θ−_θhθh)g+_Fc+_Su

∂(_ρhθ)∂t+∂(_{ρh ui}θ)∂_xi=∂∂_xj[(μ+_μtσh)(∂θ∂_xi) ]

where _Ui (i = 1, 2, 3) are the velocity component in the x, y, and z directions, respectively, p and θ are the pressure and potential temperature, μ and _μt are the molecular and the turbulent eddy viscosities, g is the gravitational acceleration, and _σh is the turbulent Prandtl number. _Fc is the Coriolis force defined as _{εijk fc Uk} , where _fc=2Ωsinλ is the Coriolis parameter, a function of the Earth’s angular velocity Ω and the latitude λ . _Su is the source term representing a forest canopy and described later by Equation (14). The hydrostatic fluid density _ρh is given in a hydrostatic reference state (subscript 0) as a function of the hydrostatic pressure _ph and the temperature _Th as:

_ρh=_{phRd Th}

_Th=_T0 ²−2 A g z_Rd

_ph=_p0(−_T0A+(^_T0A2)−2 A g z_RdA)

with the constant reference pressure p₀ usually set to 1000 hPa , _T0 is the reference temperature equal to 288.5 K , A = 50 K and _Rd = 287.05 J·k^g−1·^K−1 according to [23,24].

We introduce the hydrostatic potential temperature which is, by definition, a function of the temperature and the pressure according to an adiabatic state change as follows:

_θh=_Th ^(_p0ph)RCp

To close the given set of equations above, the k−ε model was used. The standard k−ε model, when applied to ABL studies, is known to lead to a turbulent length scale l , that grows approximately linearly with the height and thus, possesses very deep boundary layers. In real ABL flows, the turbulent length scale is limited by the ABL height or the stability [25,26]. To overcome this problem, the mixing length model proposed by [25], which introduces a limiting size of turbulent eddies _lmax in the ABL is used. The two modified transport equations for the turbulent kinetic energy k and the dissipation ε read:

∂(_ρhk)∂t+∂(_{ρh uj}k)∂_xj=∂∂_xj[(μ+_μtσk)(∂k∂_xj) ]+P+G+_Sk−_ρhε

∂(_ρhε)∂t+∂(_{ρh uj}ε)∂_xj=∂∂_xj[(μ+_μtσε)(∂ε∂_xj) ]+_Cε1*(P+G)+_Sε−_Cε2^ε2k

where P represents the production rate of turbulent kinetic energy due to shear and G represents the production/destruction of turbulence by buoyancy forces defined as:

P=_τij(∂_ui∂_xj)=_μt(∂_ui∂_xj+∂_uj∂_xi )∂_ui∂_xj−23_ρh ∂_uk∂_xkδij

G=−1Tg_μtσh∂θ∂_xi

_Sk and _Sε are source terms relative to the forest canopy and described later in Equations (15) and (16). _σk, _σε , _Cε1* , _Cε2 are model coefficients listed in Table 1. These constants were adapted to atmospheric conditions, as proposed by [26]. _Cε1* introduces the maximum mixing length _lmax in the following:

_Cε1*=_Cε1+(_Cε2−_Cε1)l_lmax

where the mixing length l is equal to the dissipation length defined as _lε=(_Cμ34 ^k32)/ε . Several mixing-length models in the literature enable estimation of the limiting turbulent size of turbulent eddies in the ABL _lmax , see [27] for a review. For neutral flows, this length is computed using the Blackadar equation [28] as:

_lmax=0.00027_Ug2Ωsinλ

where _Ug is the geostrophic wind velocity.

In order to account for the effects of vegetation on the wind flow, terms in the transport equation of momentum Equation (2), turbulent kinetic energy Equation (8), and turbulence dissipation rate Equation (9) were added:

_Su=−_{ρh Cd}LAD(z)|U|_ui

_Sk=−_{ρh Cd}LAD(z)(_βp ^|U|3−_βd|U|k)

_Sε=−_{ρh Cd}LAD(z)εk(_{Cε4 βp} ^|U|3−_{Cε5 βd}|U|k)

where |U| is the velocity magnitude and _Cd the leaf drag coefficient. The values of this drag coefficient vary between 0.1 and 0.3 for most of the vegetation [30]. A value of 0.2 will be considered for the rest of the study. _βp, _βd, _Cε4 , and _Cε5 are model constants proposed by [30] and reported in Table 2. Vegetation is discretized into finite volumes where the total amount of leaves per given volume is defined by the leaf area density ( LAD ) and the sum of each layer’s LAD value over the total canopy height H is called the leaf area index ( LAI ). A typical profile of the LAD for a deciduous forest can be seen in Figure 1. The relationship between the LAI and LAD can be expressed as follows:

LAI=_∫0hLAD(z) dz

The URANS equations were solved using a PIMPLE pressure-velocity coupling algorithm. The PIMPLE solver is considered as it offers the possibility to visualize the transient effects, which will be part of our future investigation. The gradient terms and Laplacian terms were discretized using a 2nd order linear interpolation, while a 2nd order upwind interpolation was applied for divergence terms. The convergence criteria for residuals were ensured to be less than ¹⁰⁻⁶ .

2.2. k-ε Model Validation

The numerical model was verified using the solution of the Leipzig wind profile. This profile was measured from a set of 28 double-theodolite balloons between 9:15 and 16:15 on the 20 October 1931 as reported in [31]. The ABL was considered as horizontally-homogeneous with a geostrophic wind speed _Ug=17.5 m·^s−1 and a Coriolis parameter _fc=1.13×^{10−4 s−1} corresponding to a latitude of 45 degrees N. Simulations were carried out using two k−ε models with the constants adapted to atmospheric conditions, as listed in Table 1, except for _Cε1 . _Cε1 , which were considered as constant in the first model and as a function of the size of turbulent eddies in the second model (Equation (12)), called respectively k−ε unlimited and k−ε limited. A maximum mixing length _lmax of 36 m was used in the k−ε limited, as reported by [25,26]. A computational domain of 1 km × 1 km × 3 km in the x, y, and z directions, spatially resolved with 50 × 50 × 200 cells was defined. A geostrophic wind velocity _Ug=(17.7, 0, 0) m·^s−1 at the top of the domain and at a constant roughness length of _z0=0.3 m was prescribed on the ground. The boundaries at the domain sides were defined as cyclic.

The results of the two models, together with the experimental data, are shown in Figure 2. Figure 2a,b shows vertical profiles of the wind component v perpendicular to the geostrophic wind and the component u , parallel to the geostrophic wind _Ug . The measured vertical profile of the velocity shows a low-level jet and shear on the wind direction. The Coriolis term in the model lead to a turning in the boundary layer which is depicted by the velocity component v . While the limited model simulated the profile reasonably well, the unlimited model showed a flat profile. Figure 2c shows a surface wind turned by 24.3° to the geostrophic wind for the limited model. This is close to the measured value of 26°. The unlimited version predicts a smaller cross-isobar angle of 12.0° due to the flat v -profile.

The limited variant leads to turbulent kinetic energy and viscosity decreasing over the depth of the boundary layer while the standard variant shows a more or less linear increase causing a very deep boundary layer (Figure 2d–f). Overall, this validation study shows that the limiting effect is significant. For the Leipzig test case, the limited k−ε model yields results that have a better agreement with the observed profiles. Due to overestimation of the eddy viscosity with increasing height, the effects captured by the unlimited k−ε model were underestimated.

3. Design of Experiments (DoE) Study

After verifying the capability of the model to reproduce the flow features in the ABL, the next step was to identify primary contributors to the wind flow at the WINSENT test site. One of the main features of the test site is an escarpment covered by heterogeneous vegetation. The two planned test wind turbines, with a hub height of 75 m and a rotor diameter of 50 m will be installed approximately 200 m downstream of the forested escarpment and will be directly influenced by the canopy. Finding the impact of the forest height and density (which varies seasonally as the foliage grows and develops) on the wind flow can be assessed using a DoE approach. Combining CFD simulations with a DoE can be used to accurately rank the importance of the design parameters in a study [32] and can significantly reduce, for a real test site, the amount of simulated case and computational time. In our study, the influence of three parameters (namely the slope of the escarpment α, the LAI, and the forest height H) on the horizontal velocity, the flow inclination angle, and the turbulent kinetic energy at different locations, corresponding to the future position of the turbines was investigated. Table 3 lists the three parameters and their assigned values.

In order to apply the DoE, a simplified test case was defined as shown in Figure 3. The test case represents a simplification of the test site in two dimensions. We consider a 150 m high escarpment and the slope α to be of 15° or 30°, corresponding to a West and West-North-West wind direction. The escarpment is covered by a forest (green block in Figure 3).

In the real case, the separation between the forest and the ground is not as sharp, but rather shows a smooth transition. Several studies in the literature on the flow over an escarpment, including wind tunnel, full-scale and simulations exist [33,34,35]. However, none of these studies presents results for a vegetated escarpment.

At the inflow, an empirical power law is used to describe the vertical wind profile as in [22]:

u(z)=u(_zr)^(z_zr)a

where μ(z) is the average streamwise velocity at height z. μ(z_r) = 6 m·^s−1 is the reference wind speed at the reference height _zr=80 m and a is the power law exponent. This exponent depends on the surface roughness and the thermal stability parameter α. A value of 0.14 is taken as we assume the neutral case in this DoE study. It is worth mentioning that such a wind profile is an idealized one and is rarely found in a hilly complex terrain. However, our DoE study aims to assess the impact of modelling assumptions for a vegetated escarpement on certain parameters of interest and is not conducted in order to obtain the exact profiles at our specific test site with complex terrains.

We define the speed-up ratio _Su as being the mean wind speed at a height z above the ground divided by the mean wind speed of the undisturbed flow at the same height:

_Su= u(z)_uundisturbed(z)

Similarly, we define _Sk as the speed-up ratio for the turbulent kinetic energy as:

_Sk= k(z)_kundisturbed(z)

Speed-up ratio profiles for the wind speed downstream of the escarpment along _L0, _L1, _L2 and _L3 are presented in Figure 4. All the possible configurations of variables listed in Table 3 were simulated, and, additionally, the case of an escarpment with no canopy was simulated. The case without forest allows one to distinguish between the effect caused by the slope of the escarpment and the canopy. Figure 4 reveals the influence of the canopy, particularly up to 80 m above ground level (agl.). At _L0 , i.e., at the crest, the speed-up ratio shows a reduction in the velocity due to the drag effect generated by the forest, whereas for the case with no canopy, a strong acceleration of the flow can be seen. This effect is still perceptible 1 km downstream of the escarpment (at position _L3 ). The same observation was made in the work of [33]. There, a wind tunnel investigation of the flow over several escarpment shapes with a slope of 2:1 (26.7°) and 4:1 (14.0°) was performed. It was found that the region of influence of the escarpment persists 10 H ( H being the escarpment height) downstream of the crest. The largest speed-up ratios were found close to the ground and decreased with increasing height for all cases. The location of the maximum moved upward as the flow proceeded downstream. The inclusion of a forest along the escarpment deflected the position of the maximum speed-up ratio S to higher altitudes. At L₁, the position of the maximum of _Su located between 10 m and 20 m agl. with no canopy model, and between 80 m and 90 m agl. for the cases with a canopy. The maximum values of S, for a slope of 30°, are higher than in the case of 15°. For example, at _L1 , a maximum value of 1.33 and 1.26 was reached with a slope of 30° and 15°, respectively. This effect was reduced 1 km downstream of the crest, where a maximum of 1.27 and 1.24 was reached. For the same forest height (blue vs. red lines), as expected, a higher deceleration of the flow, i.e., lower value of _Su was observed for an LAI of 5 instead of an LAI of 2. All the positions, except the crest position, show an acceleration of the flow at the relevant heights for a wind turbine (between 80 m and 100 m agl.).

Figure 5 presents the results for the speed-up ratio for the turbulent kinetic energy speed, _Sk . Speed-up ratio profiles for the turbulent kinetic energy are highly dependent on the canopy height and the LAI . The crest position ( _L0) shows a maximum turbulent kinetic energy occurring near the top of the canopy and decaying rapidly above the forest. At position _L3 , the high turbulence levels indicate that the wake generated behind the crest was not dissipated. However, the case without forest seems to be almost recovered and suggests that the changes in the turbulent kinetic energy were only due to the canopy. In general, the turbulent kinetic energy levels increased with the slope angle of the escarpment, but also with the forest height.

The results of the DoE on the 2D escarpment are shown in Figure 6, Figure 7 and Figure 8. The DoE method was applied at two points, _P1 and _P2 , located at different altitudes ( 50 m , 75 m, and 100 m agl.) on the lines _L1 and _L2 . The absolute effects of each factor, listed in Table 3, on the horizontal velocity, the inclination angle and the turbulent kinetic energy were evaluated. A similar behaviour for the horizontal velocity (Figure 6) and the flow inclination angle (Figure 7) can be seen: the main variability in the response is dominated by the slope of the escarpment. This becomes even more evident for the flow inclination angle. The effect of the forest height H or the leaf area index LAI remains small relative to the slope, except for the lowest positions at 50 m agl., where the forest height parameter can exceed the same absolute value as the angle (Figure 6). Generally, as we go further downstream of the escarpment ( _P1 vs. _P2 ), the effect for the three parameters becomes smaller. The response for the turbulent kinetic energy k was different. The forest height, followed by the LAI , was as important as the slope at lower levels. All the effects were then dissipated at 100 m agl.

Considering a turbine hub-height of 75 m and the top of the rotor blade located at 100 m agl., this parametric study shows that the primary contributor was generally the angle, i.e., for the test site, the wind direction. The forest height H and LAI play a significant role only at 50 m agl. This is in accordance with the results from Figure 4, where the main differences for all the configurations were only perceptible in the lower levels. For the simulations of the WINSENT test site, described in Section 4.1, a sensitivity analysis using the same parameter values will be conducted.

4. WINSENT Test Site Overview and Measurement Techniques 4.1. Site Description

The test site is located in the Swabian Alb in Southern Germany. The main feature of the site is an escarpment up to 200 m in height with respect to the valley (Figure 9). The escarpment is characterized by a very steep slope of around 30° from the Northwest and a more gradual slope at around 15° further west. The escarpment is covered by a deciduous forest. Downstream, it is essentially flat with an agriculturally used land. Measurements from the German Meteorological Service at Stötten (734 m a.s.l., 48.6657° latitude, 9.8646° longitude), located approximately 2 km upstream of the escarpment, indicate a westward wind direction for most of the year. This wind direction makes the test site interesting in the sense that the predominant wind direction is perpendicular to the escarpment. More details about the WINSENT test site can be found in [36].

4.2. Description of the Experiment Setup

The field measurements were carried out on the top of the plateau, over and downstream of the escarpment using the multipurpose airborne sensor carrier (MASC) operated by the Environmental Physics group at the University of Tübingen, Germany. The MASC is a UAS which can resolve turbulent fluctuations of wind and temperature with a resolution of about 30 Hz . The wind measurement was performed by measuring the flow speed and flow angles at the aircraft nose with a multi-hole flow probe. Typical absolute errors in wind speed for the sensors mounted on the MASC are of the order of ±0.5 m·^s−1 (for more information, see [4,5,37,38]). In order to measure the three-dimensional flow field at the test site, a vertical grid of racetracks with legs going back and forth over the test site was performed. The wind speed measurements were then averaged over subsections of 20 m in length along each leg. Each racetrack consisted of two horizontal straight flight legs, one in the flow direction and the second one going against the wind direction. The legs against the wind direction had a higher spatial resolution and thus were used for comparison with the simulation results.

In the present study, we consider the date 27 March 2015. The site was dominated by dry weather, almost neutral ABL, and a northwest wind (290° at the nearby weather station of Stötten during the measuring time). A vertical grid of the racetrack flight pattern with vertical steps of 25 m was performed between 13:00 and 16:00 UTC. The lowest and the highest flight heights were 75 m and 275 m agl., respectively. Two racetracks were flown at each height. In total, three flights with a duration of one hour each were carried out between 13:00 and 16:00. Figure 10 shows a three-dimensional view of the flight path performed by the MASC. The flight path covers the upstream and downstream region of the escarpment.

5. Simulation of the WINSENT Test Site 5.1. CFD Settings for the Test Site

The simulations were conducted on two nested domains with varying size and spatial resolution. A parent domain of 20 km ×20 km × 2 km and a nested domain of 5 km × 5 km × 1.5 km , centred on the escarpment, were defined. The terrain data was based on the digital terrain model DTM from the Baden-Württemberg Authorities for Spatial Information and Rural Development (LGL) with a spatial resolution of 5 m. A horizontal grid resolution of 25 m was provided for the first domain and a finer resolution of approximately 7 m was provided for the small domain (cf. Table 4). A volume mesh, finer near the ground and conforming to the site orography, was constructed for the domains using SnappyHexMesh, the mesh generator of OpenFOAM. Details of meshes can be found in Table 4. Additionally, the digitized landscape model DLM from LGL was utilized to describe the land cover on the Earth’s surface. Three classes of land uses were used for this study: urban, ground, and forest (Figure 11). The urban and the ground classes were assigned values of _z0=0.50 m and _z0=0.02 m for the aerodynamic roughness length. A no-slip boundary condition was used for the velocity. A sensitivity analysis using several LAI values and forest height was conducted (cf. Table 5). The source terms, associated with the forest, as described in Equations (14)–(16), were applied to all cells, which were located below the top of the forest canopy.

For the nested domain, boundary conditions were extracted from the solution of the parent domain. The parent domain was initialized using input data from the COSMO-DE (Consortium for Small-scale Modeling) weather model of the German Weather Service (DWD, [23]). COSMO-DE has a horizontal resolution of 2.8 km × 2.8 km and 20 vertical levels in the lower 3000 m of the computational domain. Only hourly data are delivered by the DWD. These data provide the wind components, pressure, air and surface temperature, as well as the turbulent kinetic energy. A time-slice approach was used in order to join the COSMO-DE and the CFD model. In this method, the boundary conditions for the simulation were based on one specific time slice from the COSMO-DE weather model: 27 March 2015, 15:00 UTC.

5.2. Numerical Results

In order to gain a better understanding of the wind flow at Stötten, a global view of the flow pattern at different heights agl. is presented in Figure 12. An additional black line corresponding to the legs performed by the UAS at the same altitude was added. These legs have a length of approximately 1 km . The large differences in the velocity contour plots highlight the complexity of the flow with high gradients of wind speed, especially at lower levels. A high speed-up zone along the escarpment can be seen. This speed-up extends downstream of the escarpment but not uniformly.

For the evaluation of the results, the quantities will be plotted along a horizontal line, corresponding to the performed flight path and three vertical lines (Figure 13). Error bars, corresponding to one standard deviation, were inserted for the UAS measurements.

Figure 14 shows the horizontal wind speed and the flow inclination angle at the same altitudes as those presented in Figure 12. The modelled and measured values are in good agreement, although the profile shape of the modelled wind speed at higher levels is slightly over predicted (up to 2 m·^s−1 ) relative to the UAS measurements. The simulation results and measurements show a horizontal profile which is becoming flat with increasing height. This can be explained by the flow becoming more homogeneous with increasing altitude and no longer being influenced by the terrain. At a height of 75 m , the UAS recorded the maximum wind speed at a distance of around 500 m . The simulation predicts a maximum at the same position. This maximum corresponds to the area over the escarpment with a speed-up as shown in the velocity contour plot in Figure 12a,b. This maximum was shifted to the east with increasing altitude, in accordance with the DoE study. At 125 m height, the maximum was located at a distance of 600 m according to the model. The position of this maximum cannot be seen in the measurement data. Figure 14 shows that using a forest height of 25 m lead to a positive velocity bias, except at the lower level of 75 m , where the velocity was reduced. By comparing the results with different LAI values (blue vs. red lines), we see that the foliage density does not have an impact on the horizontal wind speed. At a height of 275 m , the impact of the canopy and the escarpment was strongly reduced: all the canopy variants converge to a similar profile, with the differences getting smaller.

Flow angles, measured by the UAS and predicted by the simulation, indicate a complex flow near the escarpment. Figure 14 shows the flow inclination angles to have a non-constant profile. Such a change in the inclination angle is due to a local phenomenon induced by local features of the terrain. It is noticeable that, for the first two-thirds of the flight path, positive flow angles of up to 10 degrees were reached. This position corresponds to the area located upwind and over the escarpment that is dominated by updrafts. The flow inclination angle becomes smaller with increasing height, indicating that the wind was less deflected by the orography. At a height of 75 m , zero inclination angles appear at a distance of around 800 m for the simulation and 700 m for the measurements. At higher levels, the model no longer predicts large positive values but angles approaching a zero value, indicating a flow that is rather horizontal. In accordance with the simulation results, the UAS also reports zero-angle values at the highest level. Upon observation of the different parametrization of the forest, no significant differences were found. This is in accordance with the 2D DoE study, where the effect on the inclination angle was largely dominated by the slope of the escarpment.

Figure 15 presents the vertical profiles of the horizontal velocity and the turbulent kinetic energy at three positions along the flight path. The first position was located a few meters upstream the forest, the second one in the canopy and the last one approximately 300 m downstream of the escarpment (Figure 13). The plots confirm the overprediction of the model relative to the measurements at the upper levels. At position 1, i.e., upwind of the escarpment, the flow seems to already be decelerated. At position 2, the velocity profile shows an inflexion point in the lower height due to the explicit integration of the forest in our model. The effect of the escarpment and forest can be detected even at position 3; however, the profile becomes more uniform with increasing height. This feature is almost too difficult to observe in the UAS data due to the low resolution in the vertical direction. The strong shear, observed in the horizontal velocity profiles, indicated that there is a large amount of turbulence generated by vertical wind shear. Looking at the profiles of the turbulent kinetic energy, the modelled and observed values of the turbulent kinetic energy decreased with increasing height and attain their maximum values in the region of strong wind shear. The turbulent kinetic energy profiles were highly dependent on the canopy model. Upstream of the escarpment, the maximum simulated value was located at 100 m agl. At position 2, this maximum occurred at approximately 20 m agl., i.e., at the top of the forest, where a large portion of the turbulent kinetic energy was produced, due to high shear and shear-stress. Downstream of the escarpment, this maximum was carried to the upper heights with vales between 3.5 and 4.0 ^m2·^s−2 at a height of 40 m in the model. The UAS measurement also showed an increasing turbulence value of 2.4 ^m2·^s−2 at a height of 75 m .

Due to a low vertical spatial resolution in the measurements, it is difficult to determine the location and the intensity of this maxima. The simulation results show that the location of this maximum can be critical as it is approaching the bottom of the rotor plane of the future wind turbine. The turbulent kinetic energy profiles obtained at position 3 were the result of a flow which has been decelerated but, in the meantime, production of turbulent kinetic energy near the canopy top was generated.

6. Summary and Conclusions

This study aimed to assess whether or not our model can reproduce relevant features of the flow. The results show the potential of the proposed model based on the Boussinesq approximation and considering the Coriolis effect. The implementation of a limited version of the k−ε model has been successfully applied to predict the Leipzig experiment profiles. The limited k−ε model consists of reducing C_μ to 0.0256 and to make the production of ε as a function of height. The Coriolis force introduces a velocity component v perpendicular to the direction of the geostrophic wind and causes a turning in the wind. The limiting effect successfully reproduces the Leipzig profile by generating a surface wind turned value close to the experimental one, while the standard k-ε model obtains about half the limited value.

In the second step, we analyzed the flow on a 2D case of the escarpment. The focus was to understand the impact of the slope of the escarpment, the forest height, and the leaf area density on the horizontal velocity, the flow inclination angle, and the turbulent kinetic energy at turbine-relevant heights. A DoE was applied to the simulation results and showed that the primary contributor is the slope of the escarpment, i.e., for the test site, the wind direction. For the wind speed, it was found that the forest height and LAI plays a minor role in comparison to the slope of the escarpment. However, for the turbulent kinetic energy, we showed that the LAI and canopy height is as important as the escarpment slope.

In the final step, the simulations were performed on the WINSENT test site. The aim of the study was to evaluate the accuracy of a modified version of k-ε model in complex terrain. A validation against UAS measurements was performed. For the computation, the boundary conditions were derived from the mesoscale COSMO-DE model. Despite the low resolution of the COSMO-DE model, the one-way coupling method works well in complex terrain. The microscale model captures the varying pattern in the test site, such as the deceleration of the wind speed, the upward flow, and the increased turbulent kinetic energy in the lower levels. Generally, a good agreement was found in the lower levels (75 m and 125 m agl.) but discrepancies between the simulated results and the UAS measurements were found at upper levels. While the model shows enhanced absolute turbulent values in the lower altitudes, the measurements could not confirm the location of this maximum due to a low vertical resolution. Measurement campaigns with low flight heights, down to 20 m, are planned and will confirm or disprove the region of high turbulences. A remark on the stability should be mentioned. The assumption of an almost neutral ABL in our study was considered by checking the potential temperature profiles from UAS measurements (between 50 m and 300 m above the ground). This is also true at those levels, but only measurements near the ground could confirm the flow stratification. However, we have to recognize that it was likely a convective ABL and we should acknowledge that limitation in the paper. Additional measurement systems, such as Lidar, an eddy-covariance micro-meteorological station, and a tower equipped with an anemometer are going to be installed permanently on the test site, a few meters on the upwind and downwind side of the escarpment. This will also enable the attainment of values in lower levels and further valuable data for characterizing the flow and its stratification at the WINSENT test site.

The thermal stratification has a large impact on the vertical wind profile and turbulence levels. Therefore, different thermal stratification cases will be conducted in future studies when the complete measurement set up will be running. For a better comparison of UAS measurements and simulation results, the one-way coupling method with only one time step will be replaced by transient outputs from a mesoscale model.

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

Turbulence Model Constants	Cµ	Cε1	Cε2	σε	σk
Standard [29]	0.09	1.44	1.92	1.00	1.3
Adapted [26]	0.256	1.13	1.90	0.74	1.3

Closure Constants	βp	βd	Cε4	Cε5
Value	1	5.1	0.9	0.9

Variable Level	Variables
	Slope α (°)	Height H (m)	LAI
−1	15	20	2
+1	30	25	5

Mesh Properties	Number of Cells	Maximum Cell Size	Horizontal Resolution at the Ground	Vertical Resolution at the Ground
Basis	66.8 million	100m	25m	2.6m
Nested	25.4 million	50m	6.75m	1.2m

Run	LAI	Forest Height H (m)
Run 1 (no canopy)	0	0
Run 2	2	20
Run 3	2	25
Run 4	5	20
Run 5	5	25

Author Contributions

All the authors have substantial contribution to this work.

Funding

This publication was made possible through the collaboration in the WindForS wind energy research cluster. The authors acknowledge the State of Baden-Württemberg through bwHPC for providing computational resources and the German Federal Ministry for Economic Affairs and Energy for funding measurement campaigns in the framework of the German joint research projects Lidar complex (grant number 0325519) and KonTest (grant number 0325665), as well of the simulations within WINSENT (grant number 0324129).

Conflicts of Interest

The authors declare no conflict of interest.

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

References

1. Fraile, D.; Tardieu, P. Wind in Power 2017-Annual Combined Onshore and Offshore Wind Energy Statistics; WindEurope: Brussels, Belgium, 2018; p. 25.

2. Hofsäß, M.; Clifton, A.; Cheng, P.W. Reducing the Uncertainty of Lidar Measurements in Complex Terrain Using a Linear Model Approach. Remote Sens. 2018, 10, 1465.

3. Wildmann, N.; Höfsäß, M.; Weimer, F.; Joos, A.; Bange, J. MASC-A small Remotely Piloted Aircraft (RPA) for wind energy research. Adv. Sci. Res. 2014, 11, 55-61.

4. Wildmann, N.; Bernard, S.; Bange, J. Measuring the local wind field at an escarpment using small remotely-piloted aircraft. Renew. Energy 2017, 103, 613-619.

5. Rautenberg, A.; Graf, M.; Wildmann, N.; Platis, A.; Bange, J. Reviewing Wind Measurement Approaches for Fixed-Wing Unmanned Aircraft. Atmosphere 2018, 9, 422.

6. Sykes, R.I. An asymptotic theory of incompressible turbulent boundarylayer flow over a small hump. J. Fluid Mech. 1980, 101, 647-670.

7. Mason, P.J.; King, J.C. Measurements and predictions of flow and turbulence over an isolated hill of moderate slope. Q. J. R. Meteorol. Soc. 1985, 111, 617-640.

8. Taylor, P.A. Some numerical studies of surface boundary-layer flow above gentle topography. Bound.-Layer Meteorol. 1977, 11, 439-465.

9. Jackson, P.S.; Hunt, J.C.R. Turbulent wind flow over a low hill. Q. J. R. Meteorol. Soc. 1975, 101, 929-955.

10. Rodrigo, J.S.; Arroyo, R.A.C.; Moriarty, P.; Churchfield, M.; Kosovic, B.; Rethore, R.-E.; Hansen, K.S.; Hahmann, A.; Mirocha, J.D.; Rife, D. Mesoscale to microscale wind farm flow modeling and evaluation. WIREs Energy Environ. 2017, 6, e214.

11. Berg, J.; Mann, J.; Bechmann, A.; Courtney, M.S.; Jørgensen, H.E. The Bolund experiment, part I: Flow over a steep, three-dimensional hill. Bound.-Layer Meteorol. 2011, 141, 219.

12. Bechmann, A.; Sørensen, N.N.; Berg, J.; Mann, J.; Rethore, P.-E. The Bolund Experiment, Part II: Blind Comparison of Microscale Flow Models. Bound.-Layer Meteorol. 2011, 141, 245-271.

13. Lange, J.; Mann, J.; Berg, J.; Parvu, D.; Kilpatrick, R.; Costache, A.; Jubayer, C.; Hangan, H. For wind turbines in complex terrain, the devil is in the detail. Environ. Res. Lett. 2017, 12, 094020.

14. CORINE Land Cover-Copernicus Land Monitoring Service. Available online: https://land.copernicus.eu/pan-european/corine-land-cover (accessed on 11 February 2019).

15. Jimenez, P.A.; Navarro, J.; Palomares, A.M.; Dudhia, J. Mesoscale modeling of offshore wind turbine wakes at the wind farm resolving scale: A composite-based anlysis with the Weather Research and Forecasting model over Horns Rev. Wind Energy 2015, 18, 559-566.

16. Khani, S.; Porté-Agel, F. A Modulated-Gradient Parametrization for the Large-Eddy Simulation of the Atmospheric Boundary Layer Using the Weather Research and Forecasting Model. Bound.-Layer Meteorol. 2017, 165, 385-404.

17. Schlünzen, H.K.; Grawe, D.; Bohnenstengel, S.I.; Schlüter, I.; Koppmann, R. Joint modelling of obstacle induced and mesoscale changes-Current limits and challenges. J. Wind Eng. Ind. Aerodyn. 2011, 99, 217-225.

18. Zajaczkowski, F.; Haupt, S.; Schmehl, K. A preliminary study of assimilating numerical weather prediction data into computational fluid dynamics models for wind prediction. J. Wind Eng. Ind. Aerodyn. 2011, 99, 320-329.

19. Knaus, H.; Höfsäß, M.; Rautenberg, A.; Bange, J. Application of Different Turbulence Models Simulation Wind Flow in Complex Terrain: A Case Study for the WindForS Test Site. Computation 2018, 6, 25.

20. Knaus, H.; Rautenberg, A.; Bange, J. Model comparaison of two different non-hydrostatic formulations for the Navier-Stokes equations simulationg wind flow in complex terrain. J. Wind Eng. Ind. Aerodyn. 2017, 169, 290-307.

21. Weller, H.G.; Tabor, G.; Jasak, H.; Fureby, C. A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput. Phys. 1998, 12, 620-631.

22. Etling, D. Theoretische Meteorologie; Springer: Berlin, Germany, 2008; ISBN 978-3-540-75978-2.

23. Doms, G.; Baldauf, M. A Description of the Nonhydrostatic Regional COSMO-Model. Part I: Dynamics and Numerics; Consortium for Small-Scale Modelling (COSMO): Offenbach, Germany, 2018; p. 158.

24. Dudhia, J. A Nonhydrostatic Version of the Penn State-NCAR Mesoscale Model: Validation Tests and Simulation of an Atlantic Cyclone and Cold Front. Monthly Weather Rev. 1993, 121, 1493-1513.

25. Apsley, D.D.; Castro, I.P. A limited-lenght-scae k-ε model for the neutral and stably-stratified atmospheric boundary layer. Bound.-Layer Meteorol. 1997, 83, 75-98.

26. Detering, H.W.; Etling, D. Application of the k-ε turbulence model to the atmospheric boundary layer. Bound.-Layer Meteorol. 1985, 33, 113-133.

27. Pena, A.; Gryning, S.-E.; Mann, J.; Hasager, C.B. Lenght scales of the neutral wind profile over homogeneous terrain. J. Appl. Meteorol. Climatol. 2010, 49, 792-806.

28. Blackadar, A.K. The vertical distribution of wind and turbulent exchange in a neutral atmosphere. J. Geophys. Res. 1962, 67, 3095-3102.

29. Launder, B.E.; Spalding, D.B. The numerical computation of turbulent flows. Numer. Predict. Flow Heat Transf. Turbul. Combust. 1974, 3, 269-289.

30. Katul, G.G.; Larry, M.; Sanz, C. one- and two-equation models for canopy turbulence. Bound.-Layer Meteorol. 2004, 113, 81-109.

31. Lettau, H. A Re-examination of the "Leipzig Wind Profile" Considering some Relations between Wind and Turbulence in the Frictional Layer. Tellus 1950, 2, 125-129.

32. Elder, R.; Tourlidakis, A.; Yates, M. Advances in CFD in Fluid Machinery Design; Wiley: Hoboken, NJ, USA, 2003; ISBN 978-1-86058-353-7.

33. Bowen, A.J.; Lindley, D. A wind-tunnel investigation of the wind speed and turbulence characteristics close to the ground over various escarpment shapes. Bound.-Layer Meteorol. 1977, 12, 259-271.

34. Jensen, N.O.; Peterson, E. On the escarpment wind profile. Q. J. R. Meteorol. Soc. 1978, 719-728.

35. Emeis, S.; Frank, H.; Fiedler, F. Modification of air flow over an escarpment-Results from the Hjardemal experiment. Bound.-Layer Meteorol. 1994, 74, 131-161.

36. WINSENT-WindForS. Available online: https://www.windfors.de/en/projects/test-site/winsent/ (accessed on 5 February 2019).

37. Wildmann, N.; Rau, G.A.; Bange, J. Observations of the Early Morning Boundary-Layer Transition with Small Remotely-Piloted Aircraft. Bound.-Layer Meteorol. 2015, 157, 345-373.

38. Rautenberg, A.; Allgeier, J.; Jung, S.; Bange, J. Calibration Procedure and Accuracy of Wind and Turbulence Measurements with Five-Hole Probes on Fixed-Wing Unmanned Aircraft in the Atmospheric Boundary Layer and Wind Turbine Wakes. Atmosphere 2019, 10, 124.