Assessment of flow characteristics over complex terrain covered by the heterogeneous forest at slightly varying mean flow directions: (A case study

The impact of heterogeneous forest canopies and complex terrain on the horizontal distortion of the inflow is studied. Large-Eddy Simulation (LES) of the neutral Atmospheric Boundary Layer (ABL) flow is performed for a wind farm in Sweden for three cases associated with three different wind directions at the range of the static yaw misalignment ( ≃ ±6 ◦ ) where the yaw control system is not activated. The ground topography and forest properties for the numerical modeling are extracted from the Airborne Laser Scanning (ALS) 3D data. The wind turbines within the wind farm are introduced using the actuator disk model. To focus on the airflow deflection only by the complex terrain and vegetation, the study is limited to upstream wind turbines without any wake interaction. The predicted mean wind speed and turbulence intensity for the upstream wind turbines are compared against the nacelle-mounted anemometers taken from the wind farm’s turbine SCADA data. To quantify the additional load and moments induced at the rotor blades by the horizontal misalignment of the incoming flow, aero-structural simulation of the upstream wind turbines in the wind farm for all three cases is performed. The results show that the horizontal distortion of the inflow over the rotor swept area is usually kept below the range of static yaw misalignment ( ≃ 6 ◦ ) for the majority of the upstream wind turbines for all three cases. However, the impact of a large vertical shear exponent leading to misinterpretation of the results must be taken into consideration. Furthermore, the load imbalance of the rotor due to the vertical wind shear has the least direct contribution to the yaw moment. However, for a mean vertical shear exponent larger than 𝛼 = 0 . 25 , contrary to expectation, a positive mean yaw moment under the positive-yawed inflow may be observed.


Introduction
The technological development of wind farms requires in-depth knowledge of different areas of expertise.Challenges arise when predictions and calculations are not further developed as fast as wind power technology [1].Sweden have access to excellent wind resources (usually expressed by annual mean wind speed at specific height) for wind power developments.Generally, more wind resources are found at coastal and mountains which means that they are suitable locations for wind power.However, Swedish environmental variables such as complex terrain, forest canopies and cold climates increase the wind speed variability (normally defined by turbulence intensity and wind shear exponent) which in turn increases the uncertainty of wind power production and the turbine component cost breakdown.Moreover, wind turbines always operate within the Atmospheric Boundary Layer vegetation.Therefore, the airflow in a wind farm strongly depends on local topography and surface roughness [7].The local topography and surface roughness increase the uncertainty of on-site wind resource assessment by changing the airflow in and around wind farms.Additionally, the flow features significantly affect the aero-structural loading of rotor blades and annual power production for a wind farm located in a complex terrain.
Creating a full map of the wind field over the entire wind farm using the on-site measurement [8] cannot be affordable because of the great inhomogeneity of the flow features in a complex terrain.Apart from the time-consuming and expensive process of collecting measured data by means of on-site meteorological masts, the topographic complexity increases the uncertainty of local field measurement for a detailed description of a wind resource [9,10].Therefore, performing numerical models [11,12] to predict the flow field for wind farms located in hilly regions is upheld.
The early ABL modeling over complex terrain (mainly 2D/3D low hill) was performed using analytical theory [13][14][15] as a reliable approach in the region without separation.Advances in computational resources (computer hardware) and development of numerical algorithms for non-linear equations enabled to simulate ABL using Computational Fluid Dynamics (CFD).Nowadays, CFD has become a common numerical simulation tool to predict the atmospheric boundary layer turbulence structure ranging from mesoscale [16] to microscale modeling.Although advanced CFD methods may provide a more realistic wind farm flow prediction, it requires large computational resources.Moreover, the predicted wind farm flow by CFD is highly affected by the computational grid resolution, the specified boundary conditions, the heterogeneous surface roughness and the choice of turbulence models.
Two common CFD methodologies to simulate the ABL flow are the Reynolds-Averaged Navier-Stokes (RANS) and Large-Eddy Simulation (LES).The RANS models only predict the mean flow where Reynolds stresses are entirely modeled using different turbulence closures with reasonable computational costs.On the other hand, LES is able to predict unsteady flow features and chaotic motion of turbulent structures at very high-Reynolds number flows where the large-scale flow structures are dominant.Unlike RANS, in LES the large flow structures are resolved whereas the small flow structures are modeled using a subgrid model [17].Therefore, it can provide more accurate flow field than RANS.However, LES is a computationally expensive model, especially for simulations of flows with very high-Reynolds number wall turbulence such as ABL flow [18].
In complex terrain, the highly variable ground elevations increase the complexity of the airflow.Furthermore, the wake interactions between the turbines and the local ground topography make the flow complexity greater so that the attached boundary layer flow assumption is no longer valid.The separated flow and re-circulation regions increase the unsteadiness of the airflow and the RANS-based turbulence models are no longer capable of predicting the flow accurately.Hence, the LES based turbulence model may be considered as a promising approach to simulate the ABL flow over complex terrains.
Like terrain complexity, the vegetation as a part of ground topography also has a significant impact on ABL flow.Apart from the practical risks associated with wind turbines in forest regions [73], a higher turbulence level and wind shear due to the forest canopies have been reported in various numerical and experimental studies [36,67,[74][75][76][77][78].
Contrary to the horizontally homogeneous forest, in heterogeneous forest, LAD/PAD profile varies in both horizontal and vertical directions.Boudreault [85] and Arnqvist [86] presented numerical algorithms to extract the detailed distribution of PAD/PAI profiles for heterogeneous forest canopies from airborne laser (lidar) scan data for large areas with higher accuracy than the previous methods [87].Inclusion of forest heterogeneity in numerical flow modeling has a significant impact on the variability of wind field [85,88,89] and it is more pronounced for wind farms located in a complex terrain.Apart from the recent studies [39,72,89], to the best of authors' knowledge, there are few studies comparing the flow field characteristics extracted from LES using heterogeneous forest assumption to operational or met mast data at several locations.
As mentioned previously, the airflow within a wind farm is distorted by the surface roughness heterogeneity including complex topography and ground's heterogeneous vegetation.As a consequence, the flow structure over the swept rotor may rapidly change in both horizontal and vertical directions.The horizontal misalignment of an incoming flow with respect to the rotor axis, the so-called yawed flow, induce additional load and moments at the rotor blades [90].Moreover, it decreases the power production by a factor of  3 () where  is the yaw error in degrees [91].Therefore, a proper prediction of yaw loads and moments is crucial for the design of yaw control system and free-yaw motions, respectively [92].In the review article by Yang et al. [93], various destructive impacts of yaw misalignment error on aero-structural dynamics of the wind turbine such as imbalanced load and moments distributions has also been discussed.It has been reported [93] that for wind speeds below 20 m∕s, the statistical static yaw error (i.e., non-zero mean error of yaw misalignment) was about 5 degrees for a 3.6 MW wind turbine.In other words, in the static yaw misalignment (contrary to the dynamic yaw misalignment when the yaw motors are activated to adjust the nacelle direction back to a zero degree yaw error) the yaw control system is not activated to align the rotor with the incoming wind by moving the nacelle in order to make the zero-degree yaw error.Hence, wind turbines are always operating under the undesirable yaw misalignment condition where they encounter additional load and moments acting on the wind turbines' structures.
To study the impact of heterogeneous forest canopies and complex terrain on the horizontal distortion of the inflow and in the range of static yaw misalignment (≃6 • ), the numerical simulations of the neutral ABL flow for three different mean wind directions are performed.The neutral ABL flow assumption helps to study the turbulence generation mechanism only by means of local topography and surface roughness [48] where the impact of the buoyancy forces are neglected.
For this purpose, a high-fidelity CFD method -the so-called Large-Eddy Simulation (LES) -is employed to model the airflow inside and over complex terrains and around each wind turbine in a midwestern Swedish wind farm.To focus on the airflow deflection by the ground topography (complex terrain and vegetation), the wake interaction between adjacent turbines inside the wind farm must be avoided.Therefore, the study is limited to upstream wind turbines without wake interaction.To gain a better insight into the effect of the ground topography, the structural dynamics response of the upstream wind turbines in the wind farm for three different mean wind directions is also investigated.A commercial software (STAR-CCM+) is used to predict the atmospheric turbulence and wind profile over the wind farm.The predicted mean wind speed and turbulence intensity for the upstream wind turbines are compared against the SCADA (Supervisory Control And Data Acquisition) data of the wind farm.To study and compare the impact of three different mean wind directions on the dynamic response of the wind turbines installed in the wind farm, the computed flow field by STAR-CCM+ is supplied to the aeroelastic wind turbine simulator FAST [94].FAST is an open-source CAE tool for predicting the power production and simulating the structural and system response of wind turbines developed by NREL (National Renewable Energy Laboratory) in US.

Complex topography
The Röbergsfjället wind farm, located at Röbergskullen in southern part of Vansbro municipality in Sweden (60 • 16 ′ 49.8 ′′ N, 14 • 12 ′ 59.6 ′′ E), is used as the case study where it is partly covered by heterogeneous forest with some clearings (see Fig. 1).The farm, built in 2007, has the highest point at 543 m above sea level (a.s.l).Moreover, the difference between the highest and lowest points of the wind farm is   ≃ 284 m (see Fig. 2).It consists of eight Vestas V90 horizontal axis wind turbines (referred to as WT1-WT8) with a hub height of 90 m and a rotor diameter of 90 m, each with a capacity of 2 MW.
The complex topography of Röbergsfjället wind farm was extracted from Airborne Laser Scanning (ALS) 3D-data delivered on May 2017 by Swedish University of Agricultural Sciences, SLU (www.slu.se).The ALS data contains a point cloud with a density of 0.5-1 points per square meter for the terrain and 1-2 points per square meter for the vegetation (forest).A commercial software called Global Mapper [95] was used to derive the complex terrain coordinates with the specified horizontal resolution of 7 × 7 m 2 .The coordinates were imported into STAR-CCM+ in the STL format to generate the computational grid for the numerical simulations.The forest properties, Plant Area Density (PAD) and trees height, have been extracted from the ALS data [86,96] with the certain resolution of 20 × 20 m 2 and 2 m in the horizontal and vertical directions, respectively.
The terrain elevation with respect to sea level and forest properties in a local domain of 14 × 16 km 2 surrounding the wind farm can be seen in Fig. 2. The location of eight wind turbines of the wind farm and the location of the met mast are indicated as the red and black dots, respectively.

Measurement data
The wind speed and direction at Röbergsfjället site have been determined by the on-site meteorology mast data located at the clearing region of the North of the wind farm, close to the wind turbine 1.The met mast has been equipped with two heated cup and vane wind measurement stations designed for arctic conditions (Vaisala Wind Set WA25) at 40 and 60 m above ground.The measurement data include a full year mean and standard deviation of the wind speed and direction measured in 2006 (before the wind farm operation) with the frequency of an hour.These data are used to identify the dominant wind direction for the numerical simulation and to set-up the proper inlet boundary condition.
Fig. 3 displays the wind rose of the mean wind speed collected in 2006 at the met mast location (60 m above ground).The measurement data reveals that at 60 m above the ground the dominant wind direction is about 216 degrees with respect to the North.In the previous study done for the Röbergsfjället site [72], the choice of neutral atmospheric boundary layer flow assumption was justified by comparison between the full-day mean wind speed measurement at the met mast point (60 m above ground) with those measured only in the day-time and nighttime.In addition, the computed turbulence intensity from the measured data is associated with the class C of the international standard IEC 61400-1.

Simulation set-up
The numerical simulations are carried out over the complex terrains covered by the heterogeneous forest for three different wind directions.For each wind direction, a high-fidelity CFD approach, Large-Eddy Simulation (LES) is employed over a rectangular computational domain to predict the neutral atmospheric turbulence and time-varying wind profile for a period of 60 min with output sampling of 10 Hz.
The computational domain has dimensions 9Hx6HxH where H is set equal to 1.0 km i.e., approximately three and a half times larger than the difference between the highest and lowest points of the wind farm (  ≃ 300 m) in the vertical direction ().In other words, LES simulations are performed over a rectangular box of size L = 9.0 km × W = 6.0 km in horizontal plane with the variable wall-distance with an average height of 1.0 km where the maximum and minimum distance from the ground level are approximately 1600 m and 900 m, respectively.
The computational grids for all three cases (associated with three wind directions) have the same set-up.For all cases the size of the computational domain and simulation run-time are limited based on the available computational resources.Fig. 4 shows the layouts of the computational grids and the location of the wind turbines inside the farm.The computational domain enclosed by the black solid line has been chosen to be aligned with the dominant wind direction, 216 degrees w.r.t. to the North, at the site and it is called the reference case.The other two computational domains, enclosed by the white and red solid lines, are aligned with the wind directions 210 and 222 degrees, respectively.These two directions have been chosen to match within the static yaw misalignment limit [93] i.e., ≃ ±6 • with respect to the reference case direction.
The rotated domains (with respect to the three wind directions) provide the perpendicular inflow at inlet boundary condition meaning that the length and width of all three domains are associated with the local streamwise () and spanwise () directions, respectively.In addition, the wind farm is placed in the middle of each computational domain letting the turbulent structures to be fully developed within the distance between the inlet boundary and the wind farm.For each simulation case, all sides of the computational domain are parallel to each other; and they are flat except the floor which is made up of the topography of the area extracted from the ALS data.The computational grid size is constant in the horizontal directions with  =  = 17 m where  and  denote the local streamwise and spanwise directions, respectively.The horizontally and vertically varying heterogeneous PAD profiles over the wind farm (with a maximum trees height of 40 m) is used to take the impact of vegetation on the airflow into account.
In STAR-CCM+, the LES solver can only treat the ground as the smooth wall.Moreover, the first cell height is considered equal to 2 m to avoid increasing the number of grid cells and to maintain enough resolution near the wall for the ABL flow simulation.Thirteen grid cells are used to discretize the highest canopy height in the vertical direction while thirty vertically stretched grid cells are used to refine the wall region until 210 m.Above the height of 210 m, a constant grid spacing of  = 17 m is used.Because of the complex topography over the entire domain, the number of grid cells in the vertical direction vary between 78 and 91.
As previously stated in Section 2, the wind farm consists of eight Vestas V90 horizontal axis wind turbines with a hub height of 90 m and a rotor diameter of 90 m, each with a capacity of 2 MW.Despite the development of a generic 2 MW wind turbine model for a Vestas V90 machine in the previous study by the author et al. [72] and due to lack of official information about the power, thrust and rotational speed curves for the Vestas V90-2 MW machine, the NREL 5-MW reference wind turbine [97] will be used in this study, instead.The hub height ( ℎ ), the rotor diameter (D) and the hub diameter of the 5-MW reference wind turbine are equal to 90 m, 126 m and 3 m, respectively.Some physical and operational properties of the 5-MW reference wind turbine can be found in Appendix A.
All wind turbines are numerically modeled using the actuator disk approach within the wind farm.This can be done using the virtual disk model in the STAR-CCM+ solver.In this approach, instead of resolving the geometry of the rotor blades, the aerodynamic forces acting on the rotor blades are distributed over a cylindrical disk with a finite thickness to model the effects of a wind turbine on the surrounding flow field on the basis of the ideal horizontal axis wind turbine (1D momentum method) with wake rotation.For more details, see Section 3.2.In the STAR-CCM+, the physical properties of each actuator disk are defined using the disk thickness, inner and outer radii.For the NREL 5-MW reference wind turbine, the disk thickness, the inner radius and the outer radius are set to be 5.1 m (covering three computational cells in the streamwise direction considering the wake refinement of 10% with respect to the base cell size, i.e., 17 m), 3.0 m, 63.0 m, respectively.Moreover, to translate the power and the thrust coefficient into the axial and tangential forces using the 1D momentum method with wake rotation, an inflow velocity plane is required.The inflow velocity plane is located upstream of the virtual disk and is aligned with the unit normal vector of the virtual disk.In a transient simulation, the air velocity components are spatially averaged over the entire inflow velocity plane and they are projected at the rotor plane to obtain the averaged velocity vector at each time step.The averaged velocity vector is then used as the input variable to extract the power, thrust coefficient and rotational speed from the performance table of the turbine as presented in Table B.1.The extracted power and thrust coefficient are then used to compute the tangential and axial forces, respectively which are accordingly introduced as the body forces to the momentum equations i.e., Navier-Stokes equations.To prescribe the undisturbed inflow velocity vector projected at the actuator disk, the offset and radius of the inflow velocity plane (located at upstream the rotor disk) are taken as 1.0D and 1.10D, respectively.
A three-step grid mesh refinement using varying resolution and length (in the streamwise direction) rectangular boxes is employed for each wind turbine in the wind farm.The width (in the spanwise direction) and the height (in the wall-normal direction) of each refinement box are the same and they are equal to 170 m (=1.35D) and 215 m, respectively.The three-step grid mesh refinement has the resolution of 10%, 25% and 50%, respectively with respect to the grid cell base size (i.e.,  =  = 17 m).In addition, the length of boxes based on the rotor diameter (D = 126 m) are equal to 0.3D, 5D and 10D, respectively.A schematic of the refinement scheme for the reference case (216 • ) can be seen in Fig. 5. Adding eight actuator disks including the grid mesh refinement (around the rotor and the wake regions) increases the number of grid meshes making a total number of approximately 31 million grid cells for the entire of each computational domain.

Local topography and forest properties
Fig. 6 displays the local elevation (a.s.l), trees height and Plant Area Index (PAI) within the computational domain for the reference case (Dir.216 • ).As stated in Section 2, the forest properties were extracted from ALS data delivered on May 2017.To avoid including too many figures, only the reference case properties is presented in Fig. 6.A probability distribution of trees height and Plant Area Index (PAI) within the computational domain for all three cases are presented in Fig. 7 where they look like almost identical.The average of the trees height are in the range of 10-20 m and the maximum of trees height is restricted to 40 m.However, there are many scattered clearing zones in the local region violating the homogeneity assumption for the Röbergsfjället region.In addition, the wind farm has not been covered by a dense forest [80] because the mean value of the PAI distribution is about one which is below the averaged PAI values for the coniferous and deciduous forests in global temperate ecosystems [98].The mean PAI for all three cases is presented by Table 1.As seen, the mean PAI for the reference case (Dir.216 • ) is about 10% denser than the other two cases.However, all three cases are considered as sparse forests.The average PAD profile can be also computed from the vertical distribution of PAD over the entire computational domain as presented in Fig. 8.The reference case (Dir.216 • ) PAD profile is rather similar to the local profile.Moreover, all profiles except for the first two meters near the ground are not much different from each other.

Governing equations
Turbulent flow structures in Large-Eddy Simulation (LES) approach are usually divided into large, resolvable scales and small, subgrid scales (SGS) where the small scales are modeled.The distinction between large and SGS scales is done implicitly by the grid, referred to the grid filtering.The incompressible, grid-filtered, Navier-Stokes equations are expressed as where v , , p and  denote the velocity component in the   -direction (  = {, , }), the air density, the pressure and the kinematic viscosity, respectively and the overbar implies grid-filtered quantities.Inclusion of forest canopies and wind turbine in Eq. ( 1b) is done through the source terms   , and  , added to the momentum equations.In Eq. (1b), the Smagorinsky subgrid model is used to model the small scales given by   respectively.  and  are the constant model coefficients assumed to be 0.10 and 25.0, respectively.For the spatial and temporal discretization of the governing equations, a second-order bounded-central differencing scheme and a second-order time integration scheme are used, respectively.For the forest canopies, the source term in Eq. ( 1b) is given by where   = 0.15 denotes the forest drag coefficient [99],   is the vertical Plant Area Density (PAD) of the forest and | v| is the velocity magnitude.
The wind turbines are modeled as the actuator disk within the wind farm.This can be done using the virtual disk model in the STAR-CCM+ solver.In this approach, instead of resolving the geometry of the rotor blades, the aerodynamic forces acting on the rotor blades are introduced as the volumetric body forces through the source term  , in Eq. (1b).The aerodynamic forces are distributed over a cylindrical disk volume (virtual disk of a finite thickness) to model the effects of a wind turbine on the surrounding flow field.The source term  , in Eq. (1b) takes both axial and the tangential forces into account which are computed on the basis of the ideal horizontal axis wind turbine (1D momentum method) with wake rotation.In addition, combination of the ideal horizontal axis wind turbine with wake rotation gives the   induced velocities across the rotor plane including both the axial (a) and the tangential (a ′ ) components [100].In the STAR-CCM+, the axial and tangential forces for each grid cell whose cell center lies within the actuator disk are computed as where ,   ,   , ,  denote the torque, radius, grid cell volume, number of grid cells within the actuator disk and thrust, respectively.The computed axial and tangential forces make the source term vector  , in Eq. (1b).In Eq. (4b), the torque () and thrust ( ) are calculated using the rotor power ( ), rotational speed () and thrust coefficient (  ) table at various inflow velocities for a given wind turbine (see Table B.1) where  =  ∕.

Boundary conditions
A precursor simulation is performed over a flat rectangular domain (without any wind turbine model) with periodic boundary conditions in the streamwise and spanwise directions to supply the inlet boundary condition profiles for the simulations of the three cases.The top of the computational domain is considered as the symmetry boundary condition.Moreover, it is assumed that the ground is covered by homogeneous forest.The properties of the homogeneous forest is take from the averaged vertical PAD profile computed from the heterogeneous forest distribution at direction 216 • (red solid line in (Fig. 8) as the reference case.Moreover, the precursor simulation is set-up so that it yields the calculated mean wind speed and the turbulence intensity at the met mast point (located at 60 m above the ground, in the vicinity of wind turbine 1 (WT1) with 5.2 km distance from the inlet boundary condition and in the mid-span of the computational domain) similar to the on-site measurement data (in the absence of wind turbine) i.e., 10 m∕s and 0.15, respectively for the reference case.
For a fair comparison, an identical boundary conditions are specified for all three cases.The computational domain consists of four different types of boundary conditions: • Wall: No-slip wall boundary condition is set for the ground.The ground surface is assumed to be a smooth wall.No wall-function correction, which is usually imposed by specifying the momentum flux from standard similarity theory [80] based on roughness length, is used.The horizontally and vertically varying PAD profiles with a horizontal resolution of 20 m × 20 m is specified at each grid point of the computational domain using the nearest-neighbor interpolation scheme.
• Inlet: The inlet boundary condition is specified as Velocity Inlet.The inlet boundary condition consists of the sheared mean velocity profile taken from the precursor simulation which is superimposed to turbulent fluctuations required in LES approach.
The inlet turbulent fluctuations are generated in STAR-CCM+ using Synthetic Eddy Method (SEM) [101].The derived Reynolds stresses from the precursor simulation are specified as input into the SEM to provide the required correlation function by SEM.Since the flow is incompressible, the SEM scales the inflow fluctuations to maintain constant mass flow rate across the domain.Fig. 9 demonstrates the prescribed inlet boundary condition profiles for all three simulation cases.• Outlet: The Pressure Outlet boundary condition with zero gauge pressure is chosen at the outlet boundary condition.• Top and Sides: Symmetry boundary condition is specified for the top and the sides boundaries.An artificial flow acceleration [102] may lead by the symmetry boundary condition at the top boundary because of the short distance between the highest point of the complex terrain and the top boundary.Hence, to reduce the impact of the accelerated flow on the wind farm flow, the minimum height of the domain has been set 900 m which is equal to three times the difference between the highest and lowest points of the wind farm.Moreover, because of the imposed symmetry boundary conditions, the Coriolis force is not taken into account.
A constant time step of  = 0.1 s is used in the simulation to ensure that the Courant number is below one over the entire domain.However, for a few highly skewed grid meshes (due to the ground complex topography), there are instances that the Courant number becomes greater than one.The simulation is performed for 120 min.But, the atmospheric turbulence and time-varying wind profile for the last 61 min of the simulation (representing of six 10-minute sampling data) over the rotor swept area of all eight turbines are then extracted and exported to a aero-structural solver the so-called FAST for the aerodynamic and aeroelastic analyses.
As mentioned in Section 1, the focus of this study is on the impact of the complex topography and heterogeneous forest at three different mean wind directions limited to the static yaw misalignment range of ±6 • .Therefore, the results are only presented for WT5, WT7 and WT8 located upstream the other turbines excluding the wake interference.In addition, because of using NREL 5-MW turbine model in the simulations (instead of the existing Vestas 2-MW machine in the wind farm), for validation against SCADA data, WT5, WT7 and WT8 are again the most appropriate candidates.

Aeroelastic simulation of wind turbine using FAST
The aerodynamic and aeroelastic simulation of the wind turbines in the wind farm is done using FAST (Fatigue, Aerodynamics, Structures, and Turbulence) code [94].It is an open-source code developed by the NREL widely used within the wind energy community computing structure loads and energy production of a wind turbine.In FAST, the aerodynamic loads are computed based on the Blade Element Momentum (BEM) method while it employs a combined modal and multibody dynamics formulation for the dynamic response of the structure.
The aero-elastic simulation of the upstream wind turbines (WT5, WT7 and WT8) exposed to the unsteady inflow, provides more detailed information about the impact of complex topography and surface roughness heterogeneity on the air at slightly varying mean wind direction e.g., ±6 • .As mentioned in Section 3.3, the aerodynamic and aeroelastic simulations are carried out for 61 min where to avoid the start-up effect, the first 60 s is neglected from all time series.Among various aero-structural time series from FAST, five variables are chosen representing the cumulative effect of the unsteady inflow on various wind turbine components.These variables are electrical generator power (GenPwr), blade 1 flapwise moment (caused by flapwise forces) at the blade root (RootMyb1), tower-top pitching moment (YawBrMyp), low-speed shaft bending moment (LSSGagMys) and tower base pitching moment (caused by fore-aft forces).In FAST global coordinate system,  and z coordinates point out downwind and vertically upward (opposite to gravity), respectively.Hence, the pitching moment is acting on the pitch axis, i.e.,  axis.

Validation against SCADA
The LES predicted wind fields are compared to wind speed as recorded by nacelle anemometer and collected through the wind turbine's SCADA system on WT5, WT7 and WT8.For the site under consideration, 1 Hz SCADA data is available (with some losses) from 20th June 2017 to 3rd February 2019.The LES simulations are carried out under the assumption of neutral atmospheric condition with the heterogeneous PAD distribution measured in September when trees had leaves.As previously mentioned, the focus of this study is on the flow characteristics around WT5, WT7 and WT8 as the upstream turbines without wake-interaction.Hence, the SCADA data between 11:00-17:00 h from May 1st to October 1st (118 days) in 2018 (the year without missing measurement data) are extracted.In addition, the SCADA data are filtered for the mean wind speed equal to 9.72 m∕s i.e., the average of minimum and maximum mean wind speed predicted by the simulations for the three cases.The tolerance for binning the mean wind speed and mean wind direction (for each case) by the measurement are assumed to be ±1 and ±3, respectively.For the three cases, the upstream turbines (WT5, WT7 and WT8) are not disturbed by the wake of other turbines.However, the measured data by the nacelle anemometer (located downstream of the rotor) are likely disturbed by the rotating blades resulting in a reduction in the mean wind speed and an increase in the turbulence intensity to a slight extent.
Figs. 10 and 11 show the histogram of probability frequency of the mean wind speed (s) and turbulence intensity (Ti) measured by the nacelle anemometer at three cases for WT5, WT7 and WT8.The turbulence intensity of the measurement data is calculated as the standard deviation of wind speed for the given hub height divided by the mean wind speed (at hub height).For the LES simulation, the turbulence intensity is defined as Ti = ⟨ v′2 1 ⟩ 0.5 ∕ ⟨ v1 ⟩ at the hub height (for a period of 90 min) where ⟨ v′2 1 ⟩ 0.5 and ⟨ v1 ⟩ denote the time-averaged gridfiltered streamwise wind velocity and the streamwise Reynolds stress obtained.Since the nacelle anemometer is located at the downwind of the rotor, the measured wind velocity includes additional disturbances due to the rotating blades.Therefore, it can be expected that the predicted turbulence intensity by simulations would be less than the ones obtained from the nacelle anemometer.As seen, SCADA data reports a slightly higher mean wind speed for all three turbines than LES.In addition, the turbulence intensity predicted by LES underestimates the actual values, maybe for the reason mentioned above.The minimum difference between the predicted turbulence intensity by LES and the measurement data is seen for case 3 (Dir.222 • ) whereas the maximum difference is observed for the reference case i.e., case 2 (Dir.216 • ).Table 2 displays a summary of the comparison between the measurement and simulations for WT5, WT7 and WT8 at three cases.Except for some few cases, there is a fairly good agreement between the measurement data and LES for the mean wind speed.

Flow characteristics over the entire wind farm
Fig. 12 displays the iso-surface of the horizontal mean wind speed ) and turbulent kinetic energy, defined as  = 1 2 , at hub height (i.e., 90 m above the ground) for a period of 60 min.As seen, wind turbine WT1, WT7 and WT8 have greater hub-height mean wind speed than the other turbines, mostly related to their elevations level (see Fig. 6b).On the other hand, WT2, WT6 and WT4 experience higher inflow turbulence intensity at hub height than the others.This motivates the study of inflow properties at each wind turbine station in more detail.along the vertical axis (y) at 1D upstream for WT1, WT7 and WT8 for three wind directions extracted for a period of 60 min.The reason to choose 1D upstream is to avoid the impact of the induced velocity by the rotor on the upstream flow.

Inflow characteristics
As expected, because of the terrain complexity and forest canopies heterogeneity, the mean wind speed and turbulent kinetic energy profiles vary w.r.t. the turbines' location and wind directions.The main  focus here is on the mean yaw angle across the rotor recalling that the inlet inflow profile is veer and yaw-free and the Coriolis force is not taken into account.As seen in Fig. 13(c), for all cases a larger yawed flow occur at the lower part of the rotor plane (from 25 m to 125 m) while by increasing the distance from the ground, it decreases and approaches zero.
Contrary to WT8 experiencing positive mean yaw angle on the entire rotor swept area, WT7 mainly undergoes a negative inflow yaw misalignment for the lower-half of the rotor and a positive inflow yaw misalignment for the upper-half of the rotor for all three wind directions.Moreover, WT5 is only negatively yawed for the lower rotor region below 40 m for the case Dir.210 • .For all three wind directions, WT5.WT7 and WT8 are operating under yawed flow below the statistical static yaw error where the yaw control system is not activated.

Comparison of aerodynamic and aeroelastic characteristics
The yaw misalignment of the incoming flow the so-called cross-wind is one of the main sources of yaw loads/moments acting on a wind turbine.Moreover, wind turbines located in complex terrain are usually subjected to both vertical and horizontal wind shears, including the atmospheric turbulence.In addition, the vertical and horizontal wind shears have large contributions on the yaw moment.The importance of studying yaw moment is because the yaw loads are classified as the fully reversing fatigue cycles (acting on the yaw drive of the horizontal axis wind turbines' yaw system) where the peak loads are usually much larger than the mean loads [92].
To investigate the impact of the airflow distortion due to the complex terrain and heterogeneous forest canopies in the range of static yaw misalignment (≃6 • ), the dynamic response of the wind turbines without wake interference i.e., WT5, WT7 and WT8 for three different directions (210 • , 216 • and 222 • ) are computed using FAST simulation tool.The extracted 60-minute time-varying flow field (from the CFD simulation over the rotor swept area of the turbines) are fed into to FAST for aerodynamic and aeroelastic analysis.Since the main focus of this study is about the airflow distortion by the ground topography within a wind farm (in the range of static yaw misalignment), four relevant aero-structural variables i.e., ''Wind1VelX'' (axial wind velocity at hub-height), ''GenPwr'' (electrical generator power), ''YawBrMyp'' (pitch moment at the top of the tower about the lateral axis) and ''YawBrMzp'' (yaw moment at the top of the tower about the vertical axis) are investigated.As seen in Fig. 14(a), there is no significant difference in axial mean wind speed for all three cases (wind directions).Except for ''YawBrMyp'' and ''YawBrMzp'', the other aero-structural signals follow the trend of the streamwise mean wind speed to a great extent and are therefore not presented.In addition, the predicted power ''GenPwr'' (see Fig. 14(b) is proportionate to the ''Wind1VelX'' (axial wind velocity at hub-height) as expected.However, a slightly different behavior in ''GenPwr'' can be observed for WT8 with respect to the ''Wind1VelX'' values.This is because of a higher turbulence intensity at directions 210 and 222 degrees w.r.t.direction 216 degrees.Normally, for mean wind speeds lower than the turbine's rated wind speed (i.e., 11.4 m/s for the NREL 5-MW machine), a small increase of power due to turbulence (the time instants of the wind speed above the mean velocity) will recompense for the time instants of the wind speed below the mean velocity [103].
Generally, the overall yaw moment is caused by the out of plane, in-plane and tension forces together with flapwise bending moment.Here, in addition to ''YawBrMzp'' (yaw moment at the top of the tower about the vertical axis), ''YawBrMyp'' is presented because the azimuthally varying flapwise bending moment (or pitch moment) is the dominant component of the yaw moment [92].Furthermore, since the hub-height mean wind speeds (see Fig. 13(a)) are below the rated wind speed for NREL 5-MW (i.e., 11.4 m/s), there is no contribution of the cyclic pitch on the mean yaw moment.The sign of the yaw moment usually varies with respect to the inflow yaw angle.According to the coordinate system used in FAST, the negative and positive tower-top yaw moments are associated with the positive and negative inflow yaw angles, respectively excluding the sheared inflow ( = 0).This is also valid for low vertical shear exponents (e.g.,  ≤ 0.1).For higher vertical shear exponents (e.g.,  ≥ 0.25) than the normal wind profile model equal to 0.2 indicated in IEC 61400-1 standard [104], the yaw moment may be positive under the positive-yawed inflow.Furthermore, below the rated mean wind speed, the mean tower-top yaw moment for the inflow with zero yaw misalignment is usually positive.On the other hand, the tower-top pitch moment is always positive for both negative and positive yawed inflow while it is larger for inflows with negative yaw misalignments.Recalling the predicted mean yaw angle for the incoming flow at 1D upstream WT5, WT7 and WT8 illustrated in Fig. 14(c), WT7 mainly experiences a negative inflow yaw misalignment for the lower-half of the rotor and a positive inflow yaw misalignment for the upper-half of the rotor for all three wind directions.Noted that the predicted yawed inflow at 1D upstream may be deflected more while approaching the rotor because of the secondary flow induced by the rotor [105].In addition, all three cases for WT7 have the same mean vertical shear exponent across the rotor swept area ( ≃ 0.23).Hence, the lowest tower-top pitch and yaw moments (see Fig. 14(c) and (d)) in case 3 for WT7 denote a slightly positive yawed inflow contrary to the case 2 and case 3.
WT8 experiences the positive yawed inflow (because of the negative yaw moment values) upheld by Fig. 13(c).The vertical shear exponent across the rotor for WT8 as seen in Fig. 13(a) is the lowest ( ≃ 0.13) in all three wind directions and it is below the power law exponent assumed to be 0.2 as for the normal wind profile model in the international standard IEC 61400-1 [104].Hence, the impact of the vertical shear exponent on the yaw moment may be not pronounced as for WT5 and WT7.For WT5, the mean vertical shear exponent across the rotor is the largest (see Fig. 13(a)) and it is approximately equal to 0.4, i.e., three times and twice greater than the ones hold for WT7 and WT8, respectively.This may be certified by the largest predicted tower-top pitch moment demonstrated in Fig. 14(c).Recalling the role of the pitch moment as the dominant source of the yaw moment [92], a larger pitch moment may significantly increase the yaw moment.Hence, contrary to the aforementioned general statement under zero vertically sheared inflow condition (which may determine the sign of yaw moment w.r.t. the sign of the cross-wind), the positive yaw misalignments for WT5 at all three wind directions despite the positive ''YawBrMzp'' for the case 1 (Dir.210 • ) and case 2 (Dir.216 • ) may be deduced.This is fairly supported by Fig. 13(c), except for the reference case (Dir.216 • ) where it is expected a smaller yaw moment at Dir. 216 • than Dir.210 • .
The NREL 5-MW wind turbine (used in this study) is a three-bladed wind turbine encountering the cyclic yaw load (as the cumulative load of all three rotor blades) at multiples of 3P i.e., 0P, 3P, 6P, 9P and so on where 1P denotes the rotor rotational frequency.Figs. 15 and 16 show the power spectral density of tower-top pitch and yaw moments, respectively.In both figures, the peaks associated with 3P, 6P, 9P are clearly visible.For the power spectra of the yaw moment (see Fig. 16), there is another pick at f ≃ 0.33 Hz indicated by  resembling the first natural frequency of NREL 5-MW wind turbine tower [97].This peak may not be identified in case of a simulated flow field by LES with a poor resolution [81,105] resulting in underestimation of the fatigue life of a wind turbine.

Conclusion
The results from the numerical simulations display that for an onshore wind farm located at the complex terrain covered by heterogeneous forest for three different mean wind directions (in the range of static yaw misalignment) with a varying PAI about 8%-10%, the horizontal distortion of the undisturbed inflow over the entire rotor is usually kept below the range of static yaw misalignment (≃6 • ) for the majority of the upstream wind turbines for all three cases.However, the impact of a large vertical shear exponent leading to misinterpretation of the results must be taken into consideration.
Among various forces/moments giving rise to the yaw moment, the load imbalance of the rotor due to the vertical wind shear has the least direct contribution on the yaw moment.However, for a mean vertical shear exponent larger than  = 0.25, a positive yaw moment (contrary to expectation) under the positive-yawed inflow may be observed.
The power spectral density of yaw moment reveals that how a LES-based turbulent inflow with a poor grid resolution may lead to underrate of the fatigue life of a wind turbine.Therefore, to excite the first natural frequency of the wind turbine structural motions, the cut-off frequency in the predicted inflow by LES must be similar.On the other hand, obtaining a higher cut-off frequency for high-Reynolds number flows increases rapidly the computational cost.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.The Röbergsfjället wind farm, located at Röbergskullen in southern part of Vansbro municipality, Sweden.The numbers denote the wind turbines' number.

Fig. 2 .
Fig. 2. (a) Topographic map of the Röbergsfjället region (red dots indicate the location of eight wind turbines of the wind farm), (b) Trees height of the Röbergsfjället region.(The red dots and the black dot display the location of eight wind turbines and the met mast, respectively.).

Fig. 3 .
Fig. 3. Wind rose of the mean wind speed measurement in 2006 at 60 m above the ground for Röbergsfjället site.

Fig. 4 .
Fig. 4. Layout of the three computational domains aligned with the specified wind directions, i.e., 210 • (white box), 216 • (black box) and 222 • (red box) (a) Topographic elevation (above sea level) at the Röbergsfjället wind farm, (b) Canopies height covered the Röbergsfjället wind farm.(The red dots and the black dot indicate the location of eight wind turbines and the met mast, respectively in the wind farm.)The dimension of the colorbar is meter.

Fig. 5 .
Fig. 5. (a) Schematic of the wake refinement showing the three refinement steps.(b) The grid mesh at the ground surface generated by the STAR-CCM+ solver.

Fig. 6 .
Fig. 6.Ground elevation (a.s.l) and forest properties for the reference case, (a), (c) and (e): within the computational domain, (b), (d) and (f): zoom around the wind turbines in the wind power plant.

Fig. 9 .
Fig. 9.The Prescribed inlet boundary condition profiles obtained from the precursor simulation.(a) Streamwise mean velocity and (b) Turbulent kinetic energy k, normal and shear Reynolds stresses.

Fig. 10 .
Fig. 10.Relative frequency of samples of wind speed collected at WT5, WT7 and WT8, likely from neutral conditions (118 days in 2018).Red line is the mean value of the 118 samples from the measurement data and Blue line is the mean value predicted by the LES simulations.

Fig. 11 .
Fig. 11.Relative frequency of samples of turbulence intensity collected at WT5, WT7 and WT8, likely from neutral conditions (118 days in 2018).Red line is the mean value of the 118 samples from the measurement data and Blue line is the mean value predicted by the LES simulations.

Fig. 14
Fig. 14 displays the mean values of four aforementioned signals.Recalling the predicted mean yaw angle for the incoming flow at 1D upstream WT5, WT7 and WT8 illustrated in Fig.14(c), WT7 mainly experiences a negative inflow yaw misalignment for the lower-half of the rotor and a positive inflow yaw misalignment for the upper-half of the rotor for all three wind directions.Noted that the predicted yawed inflow at 1D upstream may be deflected more while approaching the rotor because of the secondary flow induced by the rotor[105].In addition, all three cases for WT7 have the same mean vertical shear exponent across the rotor swept area ( ≃ 0.23).Hence, the lowest tower-top pitch and yaw moments (see Fig.14(c) and (d)) in case 3 for WT7 denote a slightly positive yawed inflow contrary to the case 2 and case 3.WT8 experiences the positive yawed inflow (because of the negative yaw moment values) upheld by Fig.13(c).The vertical shear exponent across the rotor for WT8 as seen in Fig.13(a) is the lowest ( ≃ 0.13) in all three wind directions and it is below the power law exponent assumed to be 0.2 as for the normal wind profile model in the international standard IEC 61400-1[104].Hence, the impact of the vertical shear exponent on the yaw moment may be not pronounced as for WT5 and WT7.For WT5, the mean vertical shear exponent across the rotor is the largest (see Fig.13(a)) and it is approximately equal to 0.4, i.e., three times and twice greater than the ones hold for WT7 and WT8, respectively.This may be certified by the largest predicted tower-top pitch moment demonstrated in Fig.14(c).Recalling the role of the pitch moment as the dominant source of the yaw moment[92], a larger pitch moment may significantly increase the yaw moment.Hence, contrary to the aforementioned general statement under zero vertically sheared inflow condition (which may determine the sign of yaw moment w.r.t. the sign of the cross-wind), the positive yaw misalignments for WT5 at all three wind directions despite the positive ''YawBrMzp'' for the case 1 (Dir.210 • ) and case 2 (Dir.216 • ) may be deduced.This is fairly supported by Fig.13(c), except for the reference case (Dir.216 • ) where it is expected a smaller yaw moment at Dir. 216 • than Dir.210 • .The NREL 5-MW wind turbine (used in this study) is a three-bladed wind turbine encountering the cyclic yaw load (as the cumulative load of all three rotor blades) at multiples of 3P i.e., 0P, 3P, 6P, 9P and so on where 1P denotes the rotor rotational frequency.Figs.15 and 16show the power spectral density of tower-top pitch and yaw moments, respectively.In both figures, the peaks associated with 3P, 6P, 9P are clearly visible.For the power spectra of the yaw moment (see Fig.16), there is another pick at f ≃ 0.33 Hz indicated by  resembling the first natural frequency of NREL 5-MW wind turbine tower[97].This peak may not be identified in case of a simulated flow field by LES with a poor resolution[81,105] resulting in underestimation of the fatigue life of a wind turbine.

Table 1
Mean Plant Area Index (PAI) for three different directions.

Table 2
Comparison between measurement data and simulation.

Table A . 1
Physical and operational properties of NREL 5-MW wind turbine.

Table B . 1
Performance data for NREL 5-MW wind turbine.