Large Eddy Simulation of Utility-Scale Wind Farm Sited over Complex Terrain

1. Introduction

Large wind farms are growing globally to reduce the energy-related carbon footprint [1,2,3,4]. Optimal onshore locations near coastlines in heavily populated regions are often not allowed for wind farm installation [5]. The other onshore sites are in rural areas, where the meteorological effects of forests and mountains can influence the local wind [6,7,8,9]. For large wind farms in rural onshore locations, the secondary circulation associated with the land-atmosphere interaction contributes to large-scale effects in wind farm layout optimization [1,10,11]. Computational analysis of characterizing the wind and atmospheric turbulence is relatively effective because measurements over the entire site are time-consuming. Creating a complete map of the local topography of the whole wind farm is affordable with modern high-performance computing facilities. However, as detailed in several recent articles [12,13,14,15,16], the fine grid causes a computational burden and becomes highly skewed if we resolve a steep mountain, which also deteriorates the numerical error.

This article employs a computational model of the local topography and the surface roughness to simulate the effect of complex terrain. Bhuiyan and Alam [16] demonstrates a scale-adaptive large eddy simulation (LES) method to simulate the land-atmosphere exchange of complex terrain. The LES of a wind farm in complex terrain can partially resolve the turbulence kinetic energy (TKE) and the Earth’s surface topography, and the subgrid model minimizes the effects of unresolved TKE and other meteorological phenomena [3,13,16,17]. Due to wake interactions and significant variability of ground elevations, the attached eddy hypothesis may not be accurate. Thus, it can enhance the dispersive stresses of the secondary motion close to the ground. The yawed flow is the horizontal misalignment of an incoming flow to the rotor axis. Studies have shown that in wind farms, neglecting the land-atmosphere exchange, yawed flow can reduce power production and increase the load and moments at the rotor blade [1,4,18,19].

This article considers a scale-adaptive LES method to account for the land-atmosphere exchange and focuses on the effects of airflow reduction by complex terrain and forests [20]. To gain a better insight into the effects of spatial inhomogeneity of the Earth’s surface, we follow the method proposed by Bao et al. [21] and incorporate the local stress of complex terrain into the momentum equation by adapting resolved stresses at grid points adjacent to the Earth’s surface. We call this approach the ‘near-surface model’, which differs from the classical wall models by the way it formulates velocity and stress conditions on the complex terrain [22,23,24,25]. Chow et al. [26] reviewed the existing scale-adaptive LES of atmospheric boundary layer flows (see also [27,28,29]). Among them are blending the characteristic length scales [28], hybrid RANS-LES in the near-surface region [27], and the canopy stress model [29]. In the current development, we employ Helmholtz’s second vortex theorem and the vortex stretching mechanism to capture the transition of characteristic scales from one flow regime to the other, such as convergence zone, pure shear, or inflection, to dynamically adjust the dissipation rate [16,17,30,31]. We model wind turbines as actuator disks, leading to a drag source in the atmospheric boundary layer [1,4,19]. Readers interested in the scale-adaptive LES method and the actuator disk model may consult the cited reference.

1.1. Surface roughness

Wind turbines extract energy at heights between 30 m and 150 m and contribute toward the dissipation of energy by surface roughness. The atmospheric sciences community has thoroughly examined the effect of surface roughness caused by mountains, forests, vegetation, etc. Surface roughness is generally represented through the aerodynamic roughness length ($z_0$) and the Monin-Obukhov similarity theory [32,33]. The impact of an isolated hill on the turbine performance depends on the turbine size and its position relative to the mountain [6,8]. Some studies show wind farm performance can be enhanced over hilly terrain [6,8]. However, other studies have observed that overall surface drag may outweigh the benefit of the local speed-up effect resulting from the disturbances of the surface protrusion (e.g., [11] and the refs in [17]). Past studies of wind farms in complex terrain [5] have either empirically tuned the earth-atmosphere exchange [34,35,36] or ignored it [1].

Wind tunnel experiments verified that the laminar-turbulent transition occurs in the boundary layer due to the surface-mounted roughness elements [34,35,37,38,39]. Cao et al. [40] studied the effects of 2D hills on the turbulent boundary layer. Ye et al. [36] investigated the impact of various roughness elements on flow resistance in mountainous regions. The surface drag and turbulent momentum transport caused by surface roughness primarily constrain the near-surface mean wind speed, dissipate a fraction of the kinetic energy, and alter the atmospheric boundary layer profile [36,40,41,42]. The representation of subgrid-scale effects of mountains within a scale-adaptive LES of atmospheric boundary layer has improved in recent years [15,16,21]

An isolated hill with a smooth surface exerts a blockage effect on the downstream side, which may reduce the power output and enhances the fatigue load on relatively small turbines [6]. The power output of a wind turbine sited on top of a hill may be (about $82\%$) higher relative to that on a flat terrain. Vanderwende and Lundquist [2] observed that the power output of a wind farm over a vegetated land increases by 14% when the roughness length ($z_0$) changed from $0.1$ m to $0.25$ m due to the increased crop height. Tobin et al. [43] observed that optimally placed artificial windbreaks may increase the power production by at least 10% [8,11,44]. In contrast, the overall power losses due to wake effects could be up to $20\%$ of the total power for wind farms over a flat terrain [45].

Each turbine also poses an additional drag and increases near-surface kinetic energy dissipation [20,46]. Therefore, sustaining the high power output of a large wind farm depends on whether another source of kinetic energy can compensate for the extra energy dissipation by a large array of turbines. For example, 3% variations of the hub height velocity caused by a complex terrain could lead to approximately 27% fluctuations in the power output [5]. The irregular topography of complex terrain, such as the Rocky Mountains and Tibetan Plateau, coastlines, or varying land use are efficient for the vertical kinetic energy transfer from aloft, and thus, can play an essential role in wind farm optimization [2,16,46,47].

1.2. Form-induced dispersive stress

Complex terrain acts as a large-scale obstruction to the flow, creating dispersive stress, of which the subgrid scale effects are essential in scale-adaptive LES of wind farms. Past studies have addressed drag reduction, empirical formulations of roughness effects in the secondary turbulent motions due to rough walls [48], forests, canopies, vegetation, coral reefs, etc. [49], and arrays of wind turbines [10,47]. The double-averaging method, introduced by Raupach and Shaw [50] for canopy turbulence [51] provides the dispersive stress of turbulent flow over wind turbines [10]. For dense canopies, the form-induced drag could significantly affect the lower canopy layer [52]. Poggi et al. [53] suggest dispersive stresses could reach up to 30% of the Reynolds stresses [49,54]. Moltchanov and Shavit [49] observed that the normal dispersive stresses are relatively more sensitive to the wake than the recirculation zone. The dispersive stress can originate from the vorticity generated by surface irregularities [55]. In the wind tunnel experiment of urban-like roughness elements, Cheng and Castro [56] observed that dispersive stresses are negligible compared to Reynolds stresses. In turbulent flow over an array of cubes, Santiago et al. [57] observed that Reynolds and dispersive stresses contribute equally [58,59]. However, recent studies suggest dispersive stresses account for 30% to 40% of the Reynolds stresses [60,61].

For wind farms in complex terrain, the secondary motion of the sparse canopy formed by wind turbines leads to enhanced dispersive stresses (e.g. [50]). However, it remains unclear whether a greater downward transport of kinetic energy may be sustained when dispersive stresses are due to the simultaneous effects of terrain complexity and wind turbines. In particular, no estimate for a sustained rate of the downward transport of kinetic energy is available, which may compensate the loss within large wind farms. This work provides the contributions of the dispersive, and Reynolds stresses to the total drag generated by wind farms in complex terrain.

1.3. Outline

This article considers a scale-adaptive LES of a wind farm consisting of $12\times 5$ turbines (with a rotor diameter of 100 m) sited over a complex terrain. Instead of explicitly resolving the terrain roughness elements, we employ the geostrophic drag law within a near-surface model to simulate the effects of the terrain. We set the vertical inflow profile of the streamwise velocity with a fixed wind speed at the hub height and employ a stochastic forcing method to inject turbulence intensity into the spanwise and wall-normal components of the inflow velocity. This approach does not recirculate nonlinear eddy dynamics on embedded inflow domains.

We assess the resolution requirements for utility-scale wind farm simulations considering standard criteria. We evaluate the influence of spatial heterogeneity and drag coefficient on the performance of wind farms. Based on the utility-scale wind farm flow fields over complex terrain, we assess the relative contributions of the dispersive and Reynolds stress to the total drag in wind farms. Then, we compare the effects of spatial heterogeneity on the atmospheric turbulence around wind farms. Here, we show that the power production of a wind farm in complex terrain outweighs the associated drag of the roughness elements. Finally, we report the relative contribution of dispersive stresses to the overall drag, where dispersive stresses are due to both the roughness elements and the turbines.

The article is organized as follows, First, Section 2 provides the mathematical description of numerical methods and the theory and implementation of the near-surface model used for the present study. Section 3 presents the validation and the analysis of the result. Finally, Section 4 summarizes the main results indicating potential future research directions.

2. Numerical Methods

2.1. Large eddy simulation (LES)

This study employs an in-house LES code [16,31,62] to simulate the high Reynolds number neutral atmospheric boundary layer flow over an array of wind turbines. The LES method solves the filtered continuity equation

$$\frac{\partial {\tilde{u}}_i}{\partial x_i}=0,$$

(1)

and the momentum equations

$$\frac{\partial {\tilde{u}}_i}{\partial t}+\widetilde{u_j}\frac{\partial {\tilde{u}}_i}{\partial x_j}=-\frac{\partial \tilde{p}}{\partial x_i}-\frac{\partial {\tau }_{ij}}{\partial x_j}+f_i\chi (\mathit{x},t),$$

(2)

where ${\tilde{u}}_i(\mathit{x},t)$ denotes the spatially filtered velocity characterized by a length scale of ${\Delta }_{\mathrm{les}}=2\sqrt[3]{{\Delta }_x{\Delta }_y{\Delta }_z}$. The subscripts $i=1,2,3$ represent streamwise, spanwise and surface-normal directions, respectively.

In Eq 2, a subgrid model should correctly represent the second last term so that $-\overline{{\tau }_{ij}{\mathcal{S}}_{ij}}$ accounts for the production of small-scale turbulence kinetic energy, where $\overline{(·)}$ represents the ensemble average, ${\mathcal{S}}_{ij}=(1/2)(\partial u_i/\partial x_j+\partial u_j/\partial x_i)$, and ${\tau }_{ij}={\widetilde{u_{i}u}}_j-{\tilde{u}}_i{\tilde{u}}_j$. Borue and Orszag [63] indicate that a model of the subgrid stress ${\tau }_{ij}$ would account for the energy cascade through the process of vortex stretching. Several recent investigations have clearly outlined the role of vortex stretching and strain self-amplication to formulate a cornerstone principle for the energy cascade in LES [30,64]. In this article, we have considered the scale-adaptive LES model [31] to represent the residual stress tensor ${\tau }_{ij}$. The scale-adaptive LES is based on the idea of energy cascade by vortex stretching [16,31], and its extension to wind farm simulation is detailed by Alam [17].

2.2. Near-surface model

Close to the Earth’s surface, resolving the energy-containing eddies with an isotropic grid is very costly [29], and an LES of the atmospheric boundary layer typically uses finer grid spacing in the vertical direction than in the horizontal [21,65]. Atmospheric LES codes treat horizontal and vertical dissipation separately, often using a damping function to adjust the eddy viscosity [15,66]. Such LES method also utilizes the conventional Monin-Obukhov similarity theory to model the effects of viscous layer.

For wind farms in complex terrain, Monin–Obukhov similarity-based boundary conditions are not appropriate if the lowest grid level is located well within the roughness sublayer [67]. The Monin-Obukhov similarity theory states that the vertical distribution of the wind shear is $\partial \overline{u}/\partial z\propto u_{*}/\left(kz\right)$, which leads to the wind profile

$$\overline{u}\left(z\right)=\frac{u_{*}}{k}ln\left(\frac{z}{z_0}\right),$$

(3)

where $z_0$ is the constant of integration and is known as the aerodynamic roughness length and $u_{*}$ is the characteristic velocity scale, also known as friction velocity. Thus, in the atmospheric boundary layer [15,22], a parametrization of the horizontal components of the turbulent stresses at $z=0$ is

$${\tau }_{i3}=-c_f{({\overline{u}}_1^2+{\overline{u}}_2^2)}^{1/2}{\overline{u}}_i,i=1,2,$$

where the frictional coefficient $c_f$ is diagnosed from the velocity $u\left(z_1\right)$ given by Eq 3 at $z=z_1$ [15,21,32], and more specifically, one would not set $z_1<z_{*}$ [67]. Here, $z_{*}$ is the height of the roughness layer.

The upper limit of the roughness sublayer ($z_{*}/z_0$) varies approximately between 10 and 150. Thus, for $z_{*}/z_0=100$ with $z_0=1$ m (e.g. for forests), the depth of the roughness sublayer ($z_{*}$) may be up to 100 m. The specification of ground shear stress via Eq 3 is not consistent with boundary layer scaling laws if $\Delta z<100$ m for wind farms in forested complex terrain.

An assumption of LES is that the subgrid-scale dissipation accounts for the subgrid-scale turbulence production [68,69]. The exact viscous dissipation $2\nu {\mathit{S}}_{\mathit{ij}}{\mathit{S}}_{\mathit{ij}}$ (note the bold face) must be equal to the sum of filtered energy flux $\Pi ={\tau }_{ij}{\mathcal{S}}_{ij}$ and resolved viscous dissipation $2\nu {\mathcal{S}}_{ij}{\mathcal{S}}_{ij}$ such that [70]

$$\Pi -2\nu {\mathcal{S}}_{ij}{\mathcal{S}}_{ij}=-2\nu {\mathit{S}}_{\mathit{ij}}{\mathit{S}}_{\mathit{ij}}.$$

A classical subgrid model based on the filtered strain rate, yielding

$$\Pi (\mathit{x},t)=c_s(\mathit{x},t){\Delta }_{\mathrm{les}}^2{\left(2{\mathcal{S}}_{ij}{\mathcal{S}}_{ij}\right)}^{3/2}$$

violates the local isotropy hypothesis in LES of wind farms. The filtered energy flux $\Pi$ will not vanish when all the scales are resolved in a wind farm [70]. Within the roughness sublayer [66], we can modify the eddy viscosity ${\nu }_{\tau }(\mathit{x},t)$ to account for the roughness effects into the filtered energy equation while evaluating Eq 3 at $z=z_1$ to provide boundary conditions. Consider the subgrid scale stress averaged over a horizontal plane in the following form [66]

$$\langle {\tau }_{i3}\rangle \approx 2[\langle \gamma \left(z\right){\nu }_{\tau }+{\nu }_T]\frac{\partial \langle u_i\rangle }{\partial z},i=1,2,$$

where $\gamma \left(z\right)$ is the damping factor that accounts for the shear-driven nature of the near-surface dynamics [29] and ${\nu }_T$ is the mean-field eddy viscosity needed to force the vertical velocity gradient to match with that from the similarity theory at $z=z_1$ [66]. The mean-field eddy viscosity model accounts for the Gabor-Heisenberg uncertainty, where the subgrid model is local in physical space and nonlocal in wave number space.

In scale-adaptive LES, we assume that i) vortex stretching drives the energy cascade and maintains the Kolmogorov energy spectrum, and ii) that turbulence production, kinetic energy cascade, and viscous dissipation are in local equilibrium [16,17]. Thus, we can relax the horizontal averaging in the mean-field eddy viscosity model proposed by Sullivan et al. [66], and consider that

$${[{\tau }_{13}^2+{\tau }_{23}^2]}^{1/4}=C_D^{1/2}u(x,y,z_1),$$

(4)

where the drag coefficient is

$$C_D={\left[k/ln(z/z_0)-{\psi }_M\left(\xi \right)\right]}^2$$

(5)

and for neutral atmospheric boundary layer flow ${\psi }_M\left(\xi \right)=0$. Combining Eq 3 with Eq 4 we estimate ${\tau }_{i3}(x,y,z_1)$ for $i=1,2$, and finally, we estimate the mean-field eddy viscosity ${\nu }_T$ from

$${\tau }_{i3}=2[{\nu }_{\tau }(x,y,z_1)+{\nu }_T]\frac{\partial u_i}{\partial z},i=1,2$$

In the above discussion, we follow Sullivan et al. [66] and Basu and Lacser [67] to account for the near-surface effects into the scale-adaptive LES method employed in this work. Recent studies, such as [16,71], have considered a similar formulation in LES of the atmospheric boundary layer flow over hilly terrain.

3. Results, and Discussions

The following analysis studies the interaction between a wind farm and complex terrain. We do not explicitly resolve irregular topography, such as mountains, coastlines, or urban variations in land use. Wind energy applications usually characterize complex terrain by the roughness class associated with a roughness length scale. This work focuses on wind farms’ flow structure and performance, considering each wind turbine as a Gaussian actuator disk [72]. First, we discuss the effect of grid resolution on ABL profiles, wake profiles, turbulence, and wind farm performance. Second, we investigate the effect of surface roughness on the performance of wind farms. We investigate various factors contributing to the increased power production in wind farms. Third, we examine the relative contribution of dispersive stresses and Reynolds stress for wind farms in complex terrain.

3.1. Simulation setup

This work simulates a $12\times 5$ array of wind turbines in the neutrally stratified atmospheric boundary layer over a complex terrain. Each turbine has a hub height of 100 m and a rotor diameter of $D=100$ m. Turbines are separated by a distance of $S_x=7$D in the x-direction and $S_y=5$D in the y-direction. The wind farm is schematically shown in Figure 1. The presence of the wind turbine array affects the upwind flow conditions known as blockage effect [73,74]. Strickland and Stevens [75] observed that the induction region is more apparent and may affect up to 13D on upwind flow conditions if the interspacing between the wind turbine is relatively less. Therefore, we kept the first row at a distance of 14D downstream of the inlet boundary at $x=0$. At the inlet boundary, $x=0$, the streamwise component of the velocity was assigned to $U\left(z\right)=(u_{*}/k)ln(z/z_0)$, where $u_{*}$ is chosen according to a desired Reynolds number ($\mathcal{R}e=6.7\times {10}^7$) based on the undisturbed velocity at hub height. We have the initial condition ${\overline{u}}_i(\mathit{x},0)=\langle U\left(z\right),0,0\rangle$. A stochastic forcing is considered at the inlet boundary $x=0$ to provide transient perturbations into the spanwise and surface-normal velocities, where $\langle U\left(z\right),v^{\prime }(0,y,z,t),w^{\prime }(0,y,z,t)\rangle$ is the boundary condition at $x=0$. To apply a stochastic force at $x=0$, we assume an ensemble of eddies with centers randomly distributed in space at $n_{\mathrm{spot}}$ locations, while their axis of rotation is aligned in the x direction. Such eddies maintain an enhanced level of variance by extracting energy from the mean flow and passing it to the perturbation fields.

To illustrate the transitional flow development in scale-adaptive LES for $z_0=1$ m, we sampled the velocity time series at $(x,y,z)=(12D,y,1D)$ in front of the center of three turbines denoted by T2, T3, and T4 in Figure 1. These three time series of velocity would represent the phenomena that atmospheric turbulence produces chaotic motions in the atmospheric boundary layer, which inhibit turbulence mixing. The streamwise component of each of the three velocity time series is shown in Figure 2. This result indicates that episodes of turbulent bursts persist in the streamwise velocity and pass through the first row of wind turbines. The outcome of applying the moving average with a window of 512 [s] into the streamwise velocity in front of turbine T2 is shown at the bottom plot of Figure 2.

In the following section, we like to understand the grid resolution and time step necessary to capture the wind turbine wakes.

3.2. Effect of mesh resolution

The present study considers a range of resolutions, including isotropic and non-isotropic refinements, to ensure that we had both poorly-resolved and well-resolved simulations. Following the criteria outlined by Pope [69], we discuss the effects of mesh resolution on the wake profiles, turbulence statistics, and the performance of a wind farm [65,76]. Table 1 summarizes 7 representative cases of LES for a mesh resolution study. The aerodynamic roughness length ($z_0$) of the inlet wind profile was set to 1 m. In LES Eq 2, one assumes that ${\tilde{u}}_i(\mathit{x},t)$ is a spatially filtered velocity. However, the numerical solution ${\tilde{u}}_i(\mathit{x},t)$ of Eq 2 is, in principle, a resolved part of the filtered velocity [69]. Thus, to extract the statistics of the resolved velocity ${\tilde{u}}_i(\mathit{x},t)$, Pope [69] suggests to average the flow field within a time interval, which is $45T_{*}$ in the present study. The flow was sampled at a time step such that CFL = 1, where $T_{*}=D/u_{*}$ is the large-eddy turnover time unit, and CFL stands for Courant–Friedrichs–Lewy condition.

3.2.1. Mean profiles and turbulence

Figure 3 shows the vertical profile of the mean streamwise velocity on the vertical mid-plane at 2D upstream of the first row of the wind farm. The mean was obtained within a time interval of size $45T_{*}$, and the result has been compared with the log-law wind profile given by Eq 3. The vertical profile of the mean streamwise velocity follows the logarithmic profile given by Eq 3 for all resolutions listed in Table 1. When the grid spacing is 40 m in all three directions, the atmospheric boundary layer flow is poorly-resolved [65]. Wurps et al. [65] suggests that a grid spacing of about 10 m is necessary for well-resolved LES so that the vertical profile of mean streamwise velocity is fully captured.

Figure 4 shows vertical profiles of the mean streamwise velocity on the vertical mid-plain at four wind turbine rows 1-4. We compare the profiles at locations 2D, 3D, 4D, and 5D behind the turbine for each row. Similarly, Figure 5 shows vertical profiles for rows 5, 7, 10, and 12. The plots in Figure 4 and Figure 5 compare the wake profiles among 7 resolutions. The results indicate that the solution converges to the wake profiles of case M10. The profiles of cases M30 and M40 show deviations from that of case M10. A grid spacing of 20 m in all directions is necessary to capture the wakes using the scale-adaptive LES.

Turbulence statistics of a wind farm provide insight into how kinetic energy is entrained from aloft. For this analysis, we consider the Reynolds decomposition

$${\tilde{u}}_i(\mathit{x},t)={\overline{u}}_i+u_i^{\prime }(\mathit{x},t),$$

(6)

where ${\tilde{u}}_i(\mathit{x},t)$ is the numerical solution of Eq 2 and ${\overline{u}}_i$ is the ensemble average within a time interval of length $45T_{*}$. The Reynolds stress resolved by the LES method is ${\tau }_{ij}^R=\overline{u_i^{\prime }{u}_j^{\prime }}$. The fraction of the turbulence kinetic energy modelled by the scale-adaptive LES method is $k_{sgs}=(1/2){\tau }_{ii}$, where the resolved turbulence kinetic energy is $k_{res}=(1/2){\tau }_{ii}^R$. The ratio of resolved to the total turbulence kinetic energy is effective resolution $\gamma =k_{res}/(k_{res}+k_{sgs})$[68]. A value of $\gamma =0.8$ indicates that the LES method has captured about 80% of the total turbulence energy. In this context, we have analyzed the shear stress $u_{*}^2$, resolved kinetic energy $k_{res}$, subgrid-scale kinetic energy $k_{sgs}$, and effective resolution $\gamma$.

Figure 6 shows vertical profiles of $u_{*}^2$, $k_{res}$, $k_{sgs}$, and $\gamma$, where $u_{*}^2={\left({\overline{u^{\prime }{w}^{\prime }}}^2+{\overline{v^{\prime }{w}^{\prime }}}^2\right)}^{1/2}$ represents the shear stress. Vertical shear stress profiles exhibit negligible dependence on the grid resolution except for cases $M30$ and $M40$, (see Figure 6a). Figure 6b shows that in most part of the boundary layer, the resolved TKE for coarser grid cases ($M30$ & $M40$) were higher than that for finer grid cases. Wurps et al. [65] also noted a similar scenario in atmospheric boundary layer simulations without wind turbines. Celik et al. [77] observed that in wall-bounded flows, resolved strain ${\overline{\mathcal{S}}}_{ij}$ tends to be small, thereby resulting in relatively less resolved dissipation $ϵ=2{\nu }_{sgs}{\overline{\mathcal{S}}}_{ij}{\overline{\mathcal{S}}}_{ij}$, which leads to higher resolved TKE.

Figure 6d illustrates the ratio ($\gamma$) of the resolved TKE to the total TKE in the wind farm, which shows a clear dependence on the resolution. Notably, finer grids consistently exhibit larger $\gamma$ values than coarser grids across the entire boundary layer. It is worth noting that as the vertical resolution increases, such as in cases $M20$-2 & $M20$-3, the value of $\gamma$ approaches that of the finest resolution case $M10$. Furthermore, the value of $\gamma$ is not constant within the boundary layer. As discussed by Pope [68] and Celik et al. [77], $\gamma >80\%$ may lead to a well-resolved condition for LES. Figure 6d indicates that $\gamma$ exceeds $90\%$ even for the coarsest grid. These findings for scale-adaptive LES of wind farms are consistent with the results of Wurps et al. [65] for the atmospheric boundary layer flows without wind turbines.

Figure 7a-d demonstrate color-filled contour plots of $Q=-(1/2){\mathcal{G}}_{ij}{\mathcal{G}}_{ij}$ for $Q>0$, where the color represents the vertical vorticity ${\omega }_z$. The second invariant Q of the velocity gradient tensor ${\mathcal{G}}_{ij}$ represents the relative magnitude of vorticity over strain, and thus $Q>0$ identifies vortex dominated regions. Figure 7a & Figure 7c show the coherent structures and wind turbine wakes for cases $M10$ and $M20$-3, respectively. Figure 7b & Figure 7d show the coherent structures in cases $M20$-1 and $M40$, respectively. For grid spacing of 20 m and 40 m, coherent structures are poorly-resolved (Figure 7b,d). When we refine the vertical grid keeping horizontal grid spacing 20 m (Figure 7c), the scale-adaptive LES method captures the coherent structures similar to case $M10$.

Figure 8a-f present contour plots of the streamwise and surface normal components of the resolved velocity at $z=100$ m and $t=45T_{*}$ for cases $M10$, $M20$-1, and $M40$, respectively. These contour plots compare the results between one well-resolved case ($M10$) and two poorly-resolved cases ($M20$-1 and $M40$). The results indicate that the three cases equally resolve the horizontal flow structures at hub height $z=100$ m.

3.2.2. Effects of resolution on the power output of the wind farm

Although Churchfield et al. [1] coupled the actuator line model with the aerodynamic structural response model, there is still a need for a computational framework capable of capturing atmosphere-turbine interactions in a realistic atmospheric environment [18,78,79]. For example, a generalized actuator disk model represents wind turbine wake effects by applying instantaneous forces to each component of the Navier-Stokes equations. Here, we employ the actuator disk theory; however, we calculate the thrust force based on the local disk-averaged velocity using a modified thrust coefficient (see [72,78]).

Thus, the power extracted by a wind turbine is

$$P_t=2a/(1-a)\rho A{\langle u_d\rangle }^3,$$

where a is the axial induction factor, $\rho$ is the air density, A is the rotor swept area, and $\langle u_d\rangle$ is the mean velocity at the rotor [78]. Note that local velocity $\langle u_d\rangle$ partially accounts for the wake effect. Due to additional complexity in properly considering the effects of all the influencing parameters, it is not straightforward to evaluate wind farm performance using the equation above. If the inflow condition is fixed, wake effects are a major source of power production losses in wind farms. Vertical transport and turbulence mixing are primary mechanisms for wake recovery for turbines operating in the wake of another turbine. It is worth mentioning that the average power of large wind farms may be limited to around $1.5{\mathrm{Wm}}^{-2}$ (see [10]).

Figure 9 compares the average power of a turbine at each row with various mesh resolutions. Cases $M20$-1 and $M40$ over-predict the power production, indicating the need for a finer mesh resolution; however, as we refine the vertical grid, such as in cases $M20$-2 and $M20$-3, the average power converges towards the well-resolved high-resolution case.

3.3. Effect of vertical mixing on wind farms in complex terrain

Wind speed over the hills is often higher than those in the areas over flat land [6]. However, it remains unclear whether a greater downward transport of kinetic energy may be sustained in mountainous and forested terrain [5]. The wake of a single wind turbine often exhibits reduced vertical mixing [78]. In stably stratified atmospheric conditions, temperature rise below the rotor tip suggests the downward flux may persist in large wind farms (see [17] and the refs therein). In the atmospheric boundary layer, the flux of mean kinetic energy at height z is approximately $\Phi \left(z\right)=-\overline{u^{\prime }{w}^{\prime }}\overline{u}\left(z\right)$. Using Eq (3) and assuming that $u_{*}^2=-\overline{u^{\prime }{w}^{\prime }}$, the energy flux $\Phi \left(z\right)=(u_{*}^3/z_0)ln(z/z_0)$ at height z increases as $z_0$ increases.

The geostrophic drag law relates the surface friction velocity $u_{*}$ (i.e. surface shear-stress ${\tau }_w=u_{*}^2$) and the aerodynamic roughness length $z_0$ [32]. In meteorological applications, the geostrophic drag law assumes that the Coriolis and the pressure gradient force are in geostrophic balance, and the atmospheric background stratification is neglected. The commonly used wind profile based on geostrophic drag law for wind resource assessment is [80]

$$u\left(z\right)-U_g=\frac{u_{*}}{k}ln\left(\frac{fz}{u_{*}}\right)+\mathcal{C}.$$

(7)

Here, $f=2\Omega sin\varphi$ represents the Coriolis force. Using Eq (7), the depth of the atmospheric boundary layer is approximated as $H_G\approx 0.16u_{*}/f$ [10]. Comparing Eq (7) with the near-surface logarithmic wind profile (Eq 3), we can relate the geostrophic wind speed $U_g$ to the ground friction velocity $u_{*}$ as

$$u_{*}=\frac{k{U}_g}{\left[ln\left(\frac{U_g}{f{z}_0}\right)-{\mathcal{C}}_{*}\right]},$$

(8)

where ${\mathcal{C}}_{*}\approx 4$ [10,81]. In conditions with high wind above flexible surface protrusions, such as crops or forests, non-dimensional arguments suggest that $z_0$ may depend on $u_{*}$ as

$$z_0={\alpha }_c{u}_{*}^2/g,$$

(9)

where ${\alpha }_c$ is the Charnock’s constant and g is the gravitational constant [32]. An important observation from Eqs (8-9) is that there is a positive wind speed bias in the inertial sublayer whenever the ground roughness elements are to exert a relatively large roughness length $z_0$.

Here, we simulate 10 flow fields in a wind farm considering 10 values of the neutral drag coefficient $C_{DN}={\left[\kappa /ln(z/z_0)\right]}^2$ (as discussed in Section 2.2, Eq 5). For each of flow fields, we scatter the resolved streamwise velocity $\tilde{u}(x,y,z)$ against the vertical coordinate z and thus, estimate $u_{*}$ and $z_0$ to get the best fit to Eq (3). Table 2 shows estimated values of $z_0$ and $u_{*}$ from these 10 velocity fields.

To deal with the overall effect of hills and other roughness elements, the concept of an effective roughness and effective surface force ($u_{*}^2$, per area and normalized by density) is not new, particularly, for hills, obstacle arrays, forests, and urban canopies [17,32,50,62]. Also, several work have focused on ‘wall modelling’ [25] and ‘near surface modelling’ [29]. The present investigation’s near-surface model accounts for the complex terrain’s overall effects into the residual stress ${\tau }_{ij}$, and hence, the wall shear stress ${\tau }_w\left(t\right)$. We have analyzed 10 LES flow fields corresponding to 10 values of effective roughness length $z_0$ listed in Table 2 to estimate the instantaneous shear stress on the Earth’s surface, ${\tau }_w\left(t\right)$. Here, we compute the shear stress as ${\tau }_w\left(t\right)=\sqrt{\left|\right|{\tau }_{13}{\left|\right|}_{\infty }^2+\left|\right|{\tau }_{23}{\left|\right|}_{\infty }^2}$, where $\left|\right|{\tau }_{i3}{\left|\right|}_{\infty }={max}_{x,y}\left|{\tau }_{i3}(x,y,\Delta z/2,t)\right|$ for $i=1,2$. Figure 10 shows the time evolution of the friction velocity $U_{*}\left(t\right)=\sqrt{{\tau }_w\left(t\right)}$ associated with 4 values of the roughness height $z_0$. Note the upper case symbol $U_{*}\left(t\right)$ for transient friction velocity. We see that a relatively smooth surface (with $z_0=5\times {10}^{-4}\mathrm{m}$) produces relatively low friction velocity $U_{*}\left(t\right)$. However, increased roughness enhances fluctuations, impacting turbulent structures and momentum flux aloft. Based on the scaling analysis of the atmospheric boundary layer [32], the influence of complex terrain on the shear stress ${\tau }_w\left(t\right)$ indicates the corresponding impact on the downward flux of energy and the geophysical potential for the wind energy density.

The primary source of kinetic energy for a wind farm is the geostrophic wind in the free atmosphere, which is transferred to the wind turbine. A basic understanding of flow phenomena associated with performance loss for downstream turbines in a wind farm may lead to improvements in wind energy harvesting. For instance, consider the wake effects in the Horns Rev wind farm [79], which consists of 80 Vestas V-80 wind turbines covering an area of $20{\mathrm{km}}^2$. The layout of Horns Rev consists of an array of 10 rows and 8 columns. The rows are rotated approximately 7 degrees counterclockwise from the North-South direction. Wu and Porté-Agel [79] demonstrate wake effects in Horns Rev wind farm using LES in a domain of approximately $16{\mathrm{km}}^2$. The simulations assumed a fixed hub-height velocity ($z_h=80\mathrm{m}$) of $8\mathrm{m}/\mathrm{s}$ and an aerodynamic roughness length of $z_0=0.05\mathrm{m}$. The LES of Horns Rev wind farm [79] provides a reference for computational studies of wind farms.

In the present study, we simulate a wind farm, where the hub height ($z_h=100\mathrm{m}$) wind speed is $10\mathrm{m}/\mathrm{s}$ and the aerodynamic roughness length is $z_0=0.05$ m. The wind farm consists of an array of $12\times 5$ wind turbines, covering approximately an area of $36{\mathrm{km}}^2$. Thus, the present wind farm is different than the Horns Rev wind farm (and that of Wu and Porté-Agel [79]). However, in Figure 11, we observe that the performance of the present wind farm is very similar to that of Wu and Porté-Agel [79]. Let us briefly discuss the impact of surface roughness on wind farm performance, where we varied the surface roughness between $z_0=5\times {10}^{-4}$ m and $z_0=1.0$ m. As discussed above, an increase of the overall roughness height $z_0$ due to hills and forests increases the downward energy flux and the wind speed in the inertial sublayer. The relative power $P_n/P_1$ is a measure for the resulting impact on the wind farm performance, where $P_n$ is the average power extracted by a turbine at n-th row. Figure 12 compares $P_n/P_1$ for 4 values of the roughness height $z_0$. For $z_0=5\times {10}^{-4}$ m, each turbine after the third row indicates to have the same performance of approximately 40% relative to the first row. There are known reasons for such an observed performance. For instance, the partial wake recovery within the wind farm after a few rows (3 in the present study) is a geophysical potential due to the vertical flux of atmospheric turbulence. For $z_0=1$ m, the wake is recovered approximately 80% after the fifth row. Thus, the overall effect of complex terrain is likely to increase the available kinetic energy for wind turbines. Clearly, there is an enhanced momentum exchange between the land and atmosphere for a fully rough surface, resulting in more energy entrainment from aloft [1,2].

The present study does not provide a detailed analysis of the local production of turbulence kinetic energy and the associated load on wind turbines [1,6]. However, Figure 12 suggests an overall reduction of turbulence kinetic energy by complex terrain in wind farms. As we know, turbulence is highly intermittent, local production of turbulence and shear stress may be responsible for fatigue load.

3.4. Relation between power and drag coefficient

As the complexity (i.e. slope, distribution, etc) of the roughness elements increases, the drag coefficient also increases, thereby causing an increase in the roughness length $z_0$ [32]. Evaluating the drag coefficient from observations and identifying the variables that influence it, such as wind speed, is a costly experimental process. For a neutrally stratified atmospheric boundary layer over a fully rough surface, $C_{DN}$ is obtained by using Eq 5 and Eq 9:

$$C_{DN}=k^2/{\left[ln(zg/{\alpha }_c{u}_{*}^2)\right]}^2.$$

(10)

Formulating a similar relationship between the power production and the drag coefficient for wind farms sited over a complex terrain would be interesting. Wind turbines may also contribute to the roughness length $z_0$. In other words, the wind farm’s density and layout can affect the drag coefficient $C_{DN}$, which in turn can affect the power output. The present study unveils a relationship between heterogeneous spatial surfaces and power dependence, which can be instrumental in optimizing power generation and drag in large-scale wind farms. We obtain the drag coefficient ($C_{DN}$) from Eq 10 for each case shown in Table 2, where $z_0$ and $u_{*}$ are computed by fitting the LES data with Eq 3. We took the value of Charnock constant, ${\alpha }_c\approx 0.0144$ [32].

Figure 13a shows a correlation between the neutral drag coefficient $C_{DN}$ and the square of the friction velocity $u_{*}$.

The best least square fit indicates a linear relationship between $C_{DN}$ and wall shear stress (i.e. squared friction velocity). This finding is consistent with Charnock’s formulation illustrating the impact of wind speed on the drag coefficient [32]. Notably, as the roughness of the terrain increases, wind speed near the ground decreases while accelerating it near the hub height or aloft and thus, leading to enhanced wind power production [6]. However, increasing the roughness of the complex terrain also increases the overall drag in the wind farm [11]. We obtain the normalized power $P=P_i/P_1$ for each of the cases listed in Table 2, where $P_i$ is power obtained for Row 1 - 12 with i representing the row number 1 to 12. We regressed $P_i/P_1$ against the drag coefficient $C_{DN}$. Figure 13b shows the variation of the power production in a wind farm as a function of the neutral drag coefficient. The best fit linear equation for the wind power as a function of $C_{DN}$ is

$$P=({\beta }_1+{\beta }_2{C}_{DN})\times {10}^{-3},$$

where the parameters ${\beta }_1$ and ${\beta }_2$ are approximately $0.747$ and $0.114$, respectively. This relationship helps advance the analytical model for analyzing the power production of wind farms in complex terrains.

3.5. Dispersive stress analysis

Characterizing secondary turbulent motions is essential to quantify the degree of spatial heterogeneity in large wind farms. However, the turbulent transport of momentum in wind farms is a very complicated phenomena due to interaction of wind turbines and complex terrain. Here, we consider the Double-Averaged Navier-Stokes (DANS) equations to evaluate the dispersive (or form-induced) stresses, which represents the turbulent momentum transport due to spatial heterogeneity.

Let us consider the Reynolds decomposition of the solution of Eq 2:

$${\tilde{u}}_i(\mathit{x},t)={\overline{u}}_i(\mathit{x},t)+u_i^{\prime }(\mathit{x},t),$$

(11)

where

$${\overline{u}}_i(\mathit{x},t)\equiv \frac{1}{\mathcal{N}}\sum _{n=1}^{\mathcal{N}}{\tilde{u}}_i(\mathit{x},t;n),$$

(12)

is the time-average of $\mathcal{N}$ snapshots of the LES resolved velocity ${\tilde{u}}_i(\mathit{x},t;n)$. In the present analysis, each snapshot is the velocity at n-th time step with CFL=1. The number of snapshot ($\mathcal{N}$) is large enough to capture the flow for a duration of 45 eddy turn over time units at a sampling rate of 1 Hz, yet average velocity varies slowly in time. Next, we consider the spatial averaging of the time-averaged velocity field,

$$\langle {\overline{u}}_i\rangle =\frac{1}{A}\int \int \overline{u}(\mathit{x},t)dxdy,$$

(13)

where A is the area of a horizontal plane parallel to the ground surface at $z=0$. The application of the dual operation of time and spatial averaging in Eq 12 & Eq 13 results in the triple decomposition of the LES resolved field $\widetilde{u_i}(\mathit{x},t)$. Thus, we have

$${\tilde{u}}_i(\mathit{x},t)\equiv \underset{{\overline{u}}_i(\mathit{x},t)}{\underbrace{\langle {\overline{u}}_i\rangle +u_i^{\prime \prime }(\mathit{x},t)}}+u_i^{\prime }(\mathit{x},t).$$

(14)

According to the above decomposition, the dispersive stress associated with the horizontal spatial averaging is

$$D_{ij}\left(z\right)\equiv \langle ({\overline{u}}_i-\langle {\overline{u}}_i\rangle )({\overline{u}}_j-\langle {\overline{u}}_j\rangle )\rangle =\langle u_i^{\prime \prime }{u}_j^{\prime \prime }\rangle .$$

(15)

Note, however, that the Reynolds stress is due to the joint variability of the randomness between the velocity fluctuations,

$$R_{ij}\left(\mathit{x}\right)=\overline{(\widetilde{u_i}-{\overline{u}}_i)(\widetilde{u_j}-{\overline{u}}_j)}=\overline{u_i^{\prime }{u}_j^{\prime }}.$$

(16)

We have computed the dispersive, and the Reynolds stresses using Eq 15 and Eq 16, respectively, and compare their relative contributions to the total drag. We obtained the maximum value of the stresses as a function of the drag coefficient $C_{DN}$ (See Eq 10). Figure 14 shows the variations of dispersive stress components with respect to the drag coefficient. The first notable observation from Figure 14a is that the streamwise normal component $\langle u_1^{\prime \prime }{u}_1^{\prime \prime }\rangle$ of the dispersive stress is relatively insensitive to the drag coefficient. However, the spanwise $\langle u_2^{\prime \prime }{u}_2^{\prime \prime }\rangle$ and surface-normal components $\langle u_3^{\prime \prime }{u}_3^{\prime \prime }\rangle$ vary linearly with the drag coefficient (Figure 14b and Figure 14c). In contrast, the corresponding Reynolds stresses vary linearly with the drag coefficient, as shown in Figure 15. It is interesting to note that the magnitude of the streamwise component of dispersive stress $\langle u_1^{\prime \prime }{u}_1^{\prime \prime }\rangle$ is higher than its Reynolds stress $\overline{u_1^{\prime }{u}_1^{\prime }}$ counterpart. This result suggests that the dominant source of spatial heterogeneity within a wind farm is in the streamwise direction, which are the dominant "wake-occupied" regions. Additionally, the dispersive stress shows minimum variations because the present study considers an aligned arrangement of 60 wind turbines. Further research could explore this idea by quantifying the degree of spatial heterogeneity with different layouts and numbers of wind turbines and by resolving the complex terrain. Furthermore, the relative contribution of the dispersive shear stress component $\langle u_1^{\prime \prime }{u}_3^{\prime \prime }\rangle$ is almost $45\%$ that of the Reynolds shear stress to the drag.

4. Conclusion

This article presents a scale-adaptive LES of wind farms in complex terrain, incorporating a near-surface model. Using the LES data of an idealized wind farm, we have analyzed the impact of complex mountainous terrain on the performance of wind farms, where the simulation does not explicitly resolve the shape of the mountain or other roughness elements. We observe that the near-surface model helps improve the accuracy of capturing the atmospheric boundary layer turbulence without needing an extremely refined mesh. We investigated the effect of grid resolution on the ABL and wake profiles and observed that the ABL wind profile is not very sensitive to grid resolution, whereas the wake profile depends on grid resolutions. When the grid has about ten grid points across the rotor, at least in the vertical direction, the wake profiles converge according to the grid refinement analysis.

Our findings indicate a 68% increase in normalized power production in wind farms situated in mountainous and forested terrain compared to a flat or offshore location. This study treats the subgrid-scale effects of terrain complexity without resolving the shape of roughness elements. We observe a linear correlation between the power production and the drag coefficient and show that drag increases as the surface roughness of the terrain increases. These findings help us understand the consequences of atmospheric turbulence on the performance of wind farms.

We provide a quantitative assessment of the secondary motion due to spatial heterogeneity in onshore sites of wind farms. We analyze the dispersive and Reynolds stresses. We follow the double-averaged Navier-Stokes equations to formulate the dispersive stresses. We show that a dominant source impacting a wind farm through spatial heterogeneity is due to the incoming turbulent flow, where the wind turbines are perpendicular to the streamwise direction. Furthermore, we observed that the dispersive shear stress contributes nearly 45% of the drag of wind farms in complex terrain compared to the Reynolds shear stress.

The development of wind farms in complex terrain poses significant challenges. The site assessment requires assessing performance due to the interaction of wind turbine wakes and atmospheric boundary layer flow. Different from flat terrain, it is challenging to characterize the complex terrain’s wind patterns because measurements at one point may not represent the flow state at points only a few hundred meters away. The present study highlights the importance of accurately modelling the subgrid-scale effects of the secondary motion due to the spatial heterogeneity imposed by the complex terrain. Future research may focus on capturing the grid-scale variation of the topography while modelling the relevant subgrid-scale effects. There are two crucial questions. Can artificial windbreaks be considered to compensate for the blockage of mountains on the reduced performance of turbines? If trees can provide potential windbreaks, it would be interesting to characterize the optimal placement of trees on mountainous terrain. Mountain and internal waves are essential to assess the potential role of atmospheric turbulence on wind farms in complex terrain. Since the middle of the twentieth century, the atmospheric science community has confirmed that complex terrain substantially impacts large-scale circulation and stationary waves, subsequently affecting the regional climate. Suppose we aim for global net-zero emissions by 2050 in the context of 1.5 degrees change scenario. In that case, we must incorporate such a large-scale effect of mountains in optimizing the wind farm layout.