Interaction of vortex stretching with wind power fluctuations

Some studies Politis et al. (2012) observed a slower wake recovery for wind turbines on the top of hills; however, other studies Hyvärinen and Segalini (2017) observed a faster wake recovery. Existing reports on the seasonal variability for the meteorological impact of large wind farms are inconsistent Baidya Roy, Pacala, and Walko (2004); Keith et al. (2004); Archer, Mirzaeisefat, and Lee (2013); Archer and Jacobson (2013); Pan and Archer (2018). The literature disagrees on the deflection of wakes caused by the Coriolis force Abkar and Porté-Agel (2016); van der Laan and Sørensen (2017). It is thus a challenging endeavor to conduct a parametric case study of the direct influence of vortex stretching on wind farms. Suitable characterizations of such phenomena are of great importance in designing strategies for increased penetration of renewable projects within the global energy mix Roy (2011), which can aid in wind farm layout optimization and control, as well as understanding the inherent variability of the aerodynamic power of wind farms Calaf, Meneveau, and Meyers (2010); Roy (2011); Bossuyt, Meneveau, and Meyers (2017); Xie and Archer (2017); Stevens et al. (2017); Alam, Afanassiev, and Singh (2019).

Liu and Stevens (2021) observed that surface protuberance could increase the aerodynamic power of wind farms; however, the optimal way of controlling roughness elements of complex terrain is not fully clear Bhuiyan and Alam (2020a). In a wind tunnel experiment with five wind turbines sited over hilly terrain, each with a rotor diameter of $25.4$ cm, Tian et al. (2013) observed that the efficiency of wind turbines improved significantly due to local speed-up effects. Bossuyt, Meneveau, and Meyers (2017) presented the measured velocity field over an array of $100$ micro turbines with a rotor diameter of $3$ cm. The outcome suggests that the power fluctuations follow the $-5/3$ spectrum of inertial range turbulence (see also the work of Bandi (2017)). Knowing the spectral bound helps estimate the ancillary reserves and influence of large-scale atmospheric turbulence.

This research investigates how the variability of atmospheric boundary layer flow over a complex terrain can influence power fluctuations in wind farms. Extensive experimental studies are being performed to understand the transfer of kinetic energy by atmospheric turbulence and the effects of ground roughness in the design of wind farms Tian et al. (2013); Bossuyt, Meneveau, and Meyers (2017); Chamorro and Porté-Agel (2011). However, challenges include the availability of fewer measurement data and even fewer numerical studies of the mechanism that transfers kinetic energy from the largest-scale motions, where it is produced, to the smallest-scale activities, where it dissipates. Most importantly, the dynamical characteristics of energy cascade Carbone and Bragg (2020); Hossen, Mulayath Variyath, and Alam (2021), especially for the efficient reduced-order representation of turbulence in practical applications, such as optimal design of wind farms, is of high importance.

Our recent works Alam (2011); Alam and Fitzpatrick (2018); Bhuiyan and Alam (2020a) developed immersed boundary methods for large eddy simulation (LES) of orographic influence in the atmospheric boundary layer. Unlike classical LESs using the Smagorinsky model, which considers only the ‘strain portion’ of the velocity gradient tensor, our approach considers both the ‘strain tensor’ and the ‘rotation tensor’ to model the eddy viscosity Hossen, Mulayath Variyath, and Alam (2021). This approach allows us to dynamically adapt the energy dissipation rates to the scale of the energetic eddies in the roughness sublayer. Note that the LES technique developed by Deardorff (1972, 1970) at the National Center for Atmospheric Research has become a prominent tool for simulating the effects of surface protrusions on the aerodynamic power of wind turbines Liu and Stevens (2021). For example, suppose computational resources permit many grid points to resolve wind turbines’ tip vortices. Then the choice of a more sophisticated “mixed model” over the standard “Samgorinsky model” may not be the primary criteria for determining simulation accuracy. Singh and Alam (2022) proposed an improved actuator disk model for wind farm simulations, and assessed the method’s performance using wind tunnel measurements. The current article presents a scale-adaptive LES of utility-scale wind farms considering an array of $41$ RE-power 5-MW wind turbines Alam, Afanassiev, and Singh (2019). It focuses on relating spectral features of wind turbines’ wakes to properties of atmospheric boundary layer flows over complex terrain Lundquist, Chow, and Lundquist (2010); Alam (2011); Arthur et al. (2019); Bhuiyan and Alam (2020a).

This article focuses on two approaches for the numerical simulation of atmospheric turbulence around wind farms. First, we consider the Navier-Stokes (NS) equations linearized around the mean profile of the atmospheric boundary layer flow Zare, Georgiou, and Jovanović (2020). The main idea is to leverage the underlying physics in the form of a prior model that arises from first principles, such as the NS equations Kevlahan (2005). Stochastically forced linearized dynamics transfers energy between the mean flow and velocity fluctuations. For instance, like in stochastic climate prediction models, white-in-time stochastic forcing of linearized NS equations lead to quasi-geostrophic turbulence Hirt et al. (2019). This approach is beneficial for directly observing the interaction between the array of wind turbines and the incoming atmospheric turbulence.

Second, we investigate the scale-adaptivity of dissipation rates Kurowski and Teixeira (2018) in LES of wind farms, where we utilize the transfer of energy via vortex stretching Carbone and Bragg (2020); Johnson (2020); Hossen, Mulayath Variyath, and Alam (2021). Note that LES accurately captures the large-scale unsteady motions in wind farms Stevens and Meneveau (2017), using fine grids combined with secondary wall models Bose and Park (2018), thanks to the recent progress in computing power. However, sufficiently fine grid resolutions are often not affordable when performing the LES of wind farms with many turbines. In our scale-adaptive LES Alam and Fitzpatrick (2018); Alam, Afanassiev, and Singh (2019); Bhuiyan and Alam (2020a), stretching a vortex is the dynamic mechanism for flow structure evolution and energy cascade, highlighting that regions of enhanced strain surround the elements of enhanced vorticity Carbone and Bragg (2020); Johnson (2020). However, a detailed understanding of the mechanisms responsible for energy transfer in turbulent flows remains elusive. Thus, our knowledge of wind farm aerodynamics is primarily empirical, albeit relying on physical intuition, numerical simulation, and experimental measurements. In the scale-adaptive criteria, filtered strain-rate tensors align perfectly with stresses arising from the low pass filter.

These observations about turbulence and wind farms suggest a need to explore the potential for scale adaptivity of subgrid-scale turbulence in the LES framework across a full range of resolutions employed by wind farm simulation models. Section II focuses on developing subgrid-scale transport schemes using the vortex stretching mechanism. Section III extends the stochastic dynamical modeling approach as a tool to transfer energy between the mean flow and velocity fluctuations in idealized wind farms. Section IV investigates the spectrum of aerodynamic power fluctuations in a simulated wind farm. Finally, section V summarizes the findings of the present research.

We now briefly extend exact mathematical results toward a large-scale circulation balance in (atmospheric) turbulence, while considering low-pass filtered velocity fields Tsinober (2001). One of the best-known mathematical results of such a circulation balance is due to Foias and Temam (1989). It asserts that the regularity and uniqueness of the velocity $\bm{u}(\bm{x},t)$ are guaranteed up to a finite time if the enstrophy (i.e. $(1/2)|\bm{\nabla}\times\bm{u}|^{2}$) of the flow remains bounded Kang, Yun, and Protas (2020). The mathematical analysis of Foias and Temam (1989) is consistent with an earlier experimental result due to Taylor (1938), which depends upon an exact consequence of the NS dynamics:

$$\frac{\partial}{\partial t}\left(\frac{1}{2}|\bm{\omega}|^{2}\right)=\bm{\omega}^{T}\bm{S}\bm{\omega}-\nu|\bm{\nabla}\bm{\omega}|^{2}.$$

(1)

Now, consider the filtered NS equations governing the atmospheric boundary layer flow Garratt (1992); Pope (2000). Let $\tau$ represent the cumulative effect of the subfilter-scale stress $\tau^{s}=[\overline{u_{i}u}_{j}-\bar{u}_{i}\bar{u}_{j}]$ plus the viscous stress $2\nu\bm{S}$. An evolution equation for filtered enstrophy is (see also Eq. 51 of Meneveau (1994))

$$\frac{\partial}{\partial t}\left(\frac{1}{2}|\bar{\bm{\omega}}|^{2}\right)=\bar{\bm{\omega}}^{T}\mathcal{S}\bar{\bm{\omega}}-\left[\bm{\nabla}\cdot(\tau-\frac{1}{3}\rm{tr}~{}\tau)\right]\cdot(\bm{\nabla}\times\bar{\bm{\omega}})$$

(2)

where $\mathcal{S}$ and $\bm{S}$ represents filtered and exact (or unfiltered) rates of strain, respectively. Estimates of finite a priori bounds on the growth of enstrophy $(1/2)|\bm{\omega}|^{2}$ (a measure of the vorticity $\bm{\omega}=\bm{\nabla}\times\bm{u}$) remain an open question on the regularity problem for the 3D NS system Kang, Yun, and Protas (2020). Wind tunnel measurements suggest that vorticity is produced by vortex stretching three times as fast as it dissipates Taylor (1938). From Eq (2), in the absence of viscous effects in the inertial range, vortex stretching is expected to be highly correlated with turbulence dissipation in order to ensure a finite bound on enstrophy.

Existing dynamic modeling approaches lean indirectly on the fact that the transfer of energy from large to small eddies is a local interaction associated with vortex stretching Richardson (1922); Taylor (1938); Lumley (1992); Davidson (2004); Doan et al. (2018). In the physical space, the transfer of energy occurs through the production of vorticity Pope (2000); Tsinober (2001); Johnson (2021); Hossen, Mulayath Variyath, and Alam (2021). The vorticity of a fluid element may be reoriented or concentrated or diffused depending on the motion of that fluid element and on the torques applied to it by the surrounding fluid elements Jimenzez et al. (1993); Saffman (1993). The role of vortex stretching in the turbulence energy cascade has been thoroughly evaluated by Johnson (2020); that study also highlights the contribution of strain self-amplification. Indeed, Batchelor and Townsend (1949) (e.g. pg 253) precisely formulates that an interaction between two neighboring energetic eddies can dissipate energy only on a very small scale Tsinober (2001).

Consequently in LES, the subfilter stress tensor $\tau^{s}_{ij}=\overline{u_{i}u}_{j}-\bar{u}_{i}\bar{u}_{j}$ is typically assumed to be precisely aligned with the eigenframe of the symmetric part of the velocity gradient tensor, $\mathcal{S}_{ij}=(1/2)(\partial\bar{u}_{i}/\partial x_{j}+\partial\bar{u}_{j}/\partial x_{i})$. Following Smagorinsky (1963), a model ($\tau_{ij}$) of the subfilter scale stress $\tau^{s}_{ij}$ is defined as

$$\tau_{ij}-\frac{1}{3}\delta_{ij}\tau_{kk}=2\nu_{\tau}\mathcal{S}_{ij},$$

(3)

where the eddy viscosity $\nu_{\tau}(\bm{x},t)$ is prescribed in a way to generate the appropriate flux of the energy, $\tau_{ij}\mathcal{S}_{ij}$. The dynamic procedure Germano et al. (1991) proposes to compute the Leonard (or resolved) stress

$$\tau^{L}_{ij}=\widetilde{\bar{u}_{i}\bar{u}}_{j}-\tilde{\bar{u}}_{i}\tilde{\bar{u}}_{j}.$$

(4)

Eq (4) filters the resolved velocity $\bar{u}_{i}(\bm{x},t)$ over a larger box of size $\alpha\Delta_{\hbox{\tiny les}}$ than the kernel of low pass filtering. The eddy viscosity is then computed as

$$\nu_{\tau}(\bm{x},t)=2c_{s}(\bm{x},t)\Delta^{2}_{\hbox{\tiny les}}|\mathcal{S}|.$$

Assuming that the flow is homogeneous (e.g. channel flow) on a plane over which an average $\langle\cdot\rangle_{h}$ can be justified, we have

$$c_{s}(\bm{x},t)=\frac{\langle\tau^{L}_{kl}\mathcal{S}_{kl}\rangle_{h}}{\alpha^{2}\Delta^{2}_{\hbox{\tiny les}}\langle|\tilde{\mathcal{S}}|\tilde{\mathcal{S}}_{mn}\mathcal{S}_{mj}\rangle_{h}-\langle|\mathcal{S}|\mathcal{S}_{pq}\mathcal{S}_{pq}\rangle_{h}}.$$

(5)

Here, $\tilde{\mathcal{S}}_{ij}=(1/2)(\partial\tilde{\bar{u}}_{i}/\partial x_{j}+\partial\tilde{\bar{u}}_{j}/\partial x_{i})$. Modern machine learning algorithms may also gradually improve the estimates of the dynamic parameter $c_{s}(\bm{x},t)$ Novati, de Laroussilhe, and Koumoutsakos (2021). The eddy viscosity may become negative on some grid points due to negative value of $c_{s}(\bm{x},t)$. A technical details of this classical dynamic procedure is given by Germano et al. (1991).

If the velocity gradient field were to characterize the small-scale eddies, it is reasonable to assume that the eddy viscosity $\nu_{\tau}(\bm{x},t)$ depends on the velocity gradient tensor via a highly nonlinear relation $\tau^{s}_{ij}=\pi_{\nu}(\nu_{\tau},\mathcal{G}_{ij})$, where the resolved velocity gradient tensor $\mathcal{G}_{ij}={\partial\bar{u}_{i}}/{\partial x_{j}}$ is the input, and residual stress tensor $\tau^{s}_{ij}$ is the output Brunton, Noack, and Koumoutsakos (2020); Novati, de Laroussilhe, and Koumoutsakos (2021). Before discussing the velocity gradient interactions for a scale-adaptive LES, it may help readers if we refer to some of the recent developments in machine learning approaches to turbulence modeling. Novati, de Laroussilhe, and Koumoutsakos (2021) propose a reinforcement learning framework for the automated discovery of $\pi_{\nu}$ using the available flow field data from experiments, field measurements, and direct numerical simulations (DNS). For instance, $\tau^{\hbox{\tiny DNS}}_{ij}=\overline{u_{i}u}_{j}-u_{i}u_{j}$ (computed from DNS field $u_{i}$) becomes the target in a reinforcement learning framework that predicts $\nu_{\tau}(\bm{x},t)$ from a snapshot of the tensor field $\mathcal{G}_{ij}$. The goal of such data-oriented methods is to fill in the serious need to calculate turbulent flows for the design of all sorts of devices, where we need relatively cheap and reliably accurate ways for the computation of turbulent flows.

The scale-similarity hypothesis, due to Leonard (1974), states that the subfilter-scale stress $\tau^{s}_{ij}$ equals the resolved stresses at scales between $\Delta_{\hbox{\tiny les}}$ and $\alpha\Delta_{\hbox{\tiny les}}$, where Eq (4) defines the Leonard (or resolved) stress. Past investigators compared the expression in Eq (4) with the true stress from DNS data Borue and Orszag (1998). The resulting correlation coefficient (of order $90$% or higher) indicates that the Leonard stress $\tau^{L}_{ij}$ correlates very well with the true subfilter scale stress $\tau^{s}_{ij}$; however, the usage of Eq (4) as a subgrid model leads to energy build-up at small scales Moser, Haering, and Yalla (2021).

For the resolved velocity $\bar{u}_{i}(\bm{x},t)$ (that is filtered over a box of size $\Delta_{\hbox{\tiny les}}$), take an average over a larger box of size $\alpha\Delta_{\hbox{\tiny les}}$. Applying the Taylor expansion Borue and Orszag (1998); Meneveau and Katz (2000),

$$\bar{u}_{i}(\bm{x},t)\approx\tilde{\bar{u}}_{i}(\bm{r},t)+(\bm{x}-\bm{r})\mathcal{G}_{ij},$$

a leading order approximation to the Leonard stress Borue and Orszag (1998); Eyink (2006) is $\tau^{L}_{ij}=c_{k}\Delta^{2}_{\hbox{\tiny les}}\mathcal{G}_{ik}\mathcal{G}_{jk}.$ Note that the second invariant of the deviatoric Leonard stress

$$\tau_{ij}^{L^{dev}}=\frac{1}{2}\left(\tau^{L}_{ij}+\tau^{L}_{ji}\right)-(1/3)\tau^{L}_{kk}\delta_{ij},$$

is $\mathcal{Q}_{\mathcal{L}}=-(1/2)\tau^{L^{dev}}_{ij}\tau^{L^{dev}}_{ji}$. After some algebraic manipulation (see Davidson (2004)), we have the following expression

$$\mathcal{Q}_{\mathcal{L}}=-\frac{1}{4}\left[\mathcal{S}_{ij}\omega_{j}\mathcal{S}_{ik}\omega_{k}+\frac{1}{3}(\mathcal{G}_{ij}\mathcal{G}_{ji})^{2}\right],$$

(6)

where the last term in the above expression represents the second invariant $\mathcal{Q}_{\mathcal{G}}=-(1/2)\mathcal{G}_{ij}\mathcal{G}_{ji}$ of the velocity gradient tensor.

It is interesting to consider the asymptotic limit of the restricted Euler dynamics Borue and Orszag (1998), where $(27/4)\mathcal{R}^{2}_{\mathcal{G}}+\mathcal{Q}^{2}_{\mathcal{G}}=0$. Now, compare the third invariant of the velocity gradient tensor,

$$\mathcal{R}_{\mathcal{G}}=-\frac{1}{3}\left[\mathcal{S}_{ij}\mathcal{S}_{ik}\mathcal{S}_{ki}+\frac{3}{4}\omega_{i}\omega_{j}\mathcal{S}_{ij}\right]$$

with the energy flux associated with Leonard stress, Eq (4),

$$\Pi:=-\mathcal{S}_{ij}\tau^{L}_{ij}=c_{k}\Delta^{2}_{\hbox{\tiny les}}\left[-\mathcal{S}_{ij}\mathcal{S}_{ik}\mathcal{S}_{ki}+\frac{1}{4}\omega_{i}\omega_{j}\mathcal{S}_{ij}\right].$$

The asymptotic limit of the restricted Euler dynamics may be interpreted as the existence of a functional that depends on $\tau^{L}_{ij}$ and represents the local rate of dissipation in turbulent flows. A negative skewness of strain, $\mathcal{S}_{ij}\mathcal{S}_{ik}\mathcal{S}_{ki}$, along with a positive value of vorticity stretching, $\omega_{i}\omega_{j}\mathcal{S}_{ij}$, will decrease $\mathcal{R}_{\mathcal{G}}$ and increase $\mathcal{Q}_{\mathcal{G}}$. It is thus evident that the stretching of vorticity would extract energy as large-scale strain is enhanced Borue and Orszag (1998); Johnson (2020); Hossen, Mulayath Variyath, and Alam (2021). Hence, an appropriate functional of $\tau^{L}_{ij}$ – may be in the form of its second invariant – can account for the turbulence dissipation rate via the vortex stretching mechanism.

To reach in the aforementioned destination, we follow Deardorff (1972), where a global quantity $k_{\hbox{\tiny sgs}}(\bm{x},t)$ was proposed for the following form of the subgrid model

$$\tau_{ij}-\frac{1}{3}\tau_{kk}\delta_{ij}=c_{k}\Delta_{\hbox{\tiny les}}\sqrt{k_{\hbox{\tiny sgs}}}\mathcal{S}_{ij}.$$

(7)

This model is widely used by the atmospheric science community, where a transport equation is solved for $k_{\hbox{\tiny sgs}}$. The idea is to ensure that the eddy viscosity $\nu_{\tau}(\bm{x},t)=c_{k}\Delta_{\hbox{\tiny les}}\sqrt{k_{\hbox{\tiny sgs}}(\bm{x},t)}$ is dynamically adjusted while the parameter $c_{k}$ is fixed.

Based on dimensional reasoning, we form a functional that maps the space of velocity gradient tensor to the field of turbulence kinetic energy, which reads in the following expression of $k_{\hbox{\tiny sgs}}$ in the context of Eq (7) (see Bhuiyan and Alam (2020a); Hossen, Mulayath Variyath, and Alam (2021)):

$$k_{\hbox{\tiny sgs}}=\frac{\Delta^{2}_{\hbox{\tiny les}}(\frac{1}{2}\mathcal{S}_{ij}\omega_{j}\mathcal{S}_{ij}\omega_{k}+\frac{1}{6}(\mathcal{G}_{ij}\mathcal{G}_{ij})^{2})^{3}}{\left[(\mathcal{S}_{ij}\mathcal{S}_{ij})^{5/2}+(\frac{1}{2}\mathcal{S}_{ij}\omega_{j}\mathcal{S}_{ij}\omega_{k}+\frac{1}{6}(\mathcal{G}_{ij}\mathcal{G}_{ij})^{2})^{5/4}\right]^{2}}.$$

(8)

The arrival at Eq (8) is a ‘heuristic’ derivation, numerically examined by our recent work Alam and Islam (2015); Walsh and Alam (2016); Bhuiyan and Alam (2020a); Hossain and Alam (2020). Several similar investigations, such as those dealing with the statistics of velocity gradient tensor (e.g., Dallas and Alexakis (2013), Buxton, Breda, and Chen (2017), Danish and Meneveau (2018)), vorticity dynamics, and vortex stretching in turbulent flows (e.g. Meneveau (1994); Borue and Orszag (1998); Nicoud and Ducros (1999); Trias et al. (2015); Carbone and Bragg (2020); Johnson (2020, 2021)) support the mathematical and the numerical analyses involved in the present research.

In order to test the scale-adaptive LES of the vortex formation, this section briefly investigates the coherent structures in the wake of a sphere. The direct forcing immersed boundary method was considered to represent the sphere. Fig 1 shows Q-isosurfaces at two Reynolds numbers: $\mathcal{R}e=3.7\times 10^{3}$ and $\mathcal{R}e=10^{5}$. Here, $\mathcal{R}e=UD/\nu$. In both the cases, the diameter $D=1$ m and the reference velocity $U=10$ m/s were fixed. The grid was uniformly refined using only $16$ grid points (i.e. $\Delta/D=0.0625$) for every length of $D$. In our LES, such a coarse grid is insufficient to resolve the smallest flow scales, the surface of the immersed solid, and the viscous boundary layer. At $\mathcal{R}e=3.7\times 10^{3}$, the Kolmogorov length scale in the vicinity of a sphere is approximately $\eta/D=0.0134$, which would decrease by up to a factor of $12$ at $\mathcal{R}e=10^{5}$.

In the sub-critical range of $\mathcal{R}e$, the dynamics of flow behind a sphere does not change significantly; i.e. the shear layer is laminar at separation from the sphere and the wake is turbulent. Moreover, the drag is due mainly to the pressure distribution on the surface of the sphere, with only a minor contribution from viscous shear stress. Fig 2 shows that the mean pressure distribution remains about the same for the two Reynolds numbers considered here. Considering that the grid is insufficient to resolve the viscous boundary layer, the pressure distribution in Fig 2 exhibits and excellent agreement with the corresponding outcome from experiments or direct numerical simulation.

Consider a statistically stationary state of turbulence in which the dynamics of velocity and pressure fluctuations $(\bm{u}^{\prime},p^{\prime})$ around the mean $(\bar{\bm{U}},\bar{p})$ satisfies the linearized NS equations. In the absence of Coriolis effects, the mean $\bar{\bm{U}}=[U(z),0,0]^{T}$ exhibits a logarithmic dependence with distance to the Earth’s surface. The linearized model is given by the following stochastically forced linear time-invariant system:

	$\displaystyle\frac{\partial u^{\prime}}{\partial t}+U(z)\frac{\partial u^{\prime}}{\partial x}$	$\displaystyle+v^{\prime}\frac{dU}{dz}$	$\displaystyle=-\frac{\partial p^{\prime}}{\partial x}+\frac{1}{\mathcal{R}e}\nabla^{2}u^{\prime}+r_{u},$
	$\displaystyle\frac{\partial v^{\prime}}{\partial t}+U(z)\frac{\partial v^{\prime}}{\partial x}$		$\displaystyle=-\frac{\partial p^{\prime}}{\partial x}+\frac{1}{\mathcal{R}e}\nabla^{2}v^{\prime}+r_{v},$		(9)
	$\displaystyle\frac{\partial w^{\prime}}{\partial t}+U(z)\frac{\partial w^{\prime}}{\partial x}$		$\displaystyle=-\frac{\partial p^{\prime}}{\partial x}+\frac{1}{\mathcal{R}e}\nabla^{2}w^{\prime}+r_{w},$
	$\displaystyle 0$		$\displaystyle=\frac{\partial u^{\prime}}{\partial x}+\frac{\partial v^{\prime}}{\partial x}+\frac{\partial w^{\prime}}{\partial x}.$

Here, $\bm{r}_{u}=[r_{u},r_{v},r_{w}]^{T}$ is a stochastic forcing term. Fourier transforms in $x$ and $y$ can parameterize the above system of partial differential equations (in $z$ and $t$) by using the horizontal wavenumbers $\bm{k}:=(k_{x},k_{y})$. It is possible to eliminate the pressure by expressing the system in which the state is given by the wall-normal velocity $w$ and the streamwise vorticity ($\omega^{\prime}_{y}=\partial_{y}u^{\prime}-\partial_{x}v^{\prime}$). Thus, the velocity fluctuations at any $\bm{k}$ is determined by the state vector $[w^{\prime},\omega^{\prime}_{y}]$. Under the assumption of scale separation between the large-scale velocity field and the fluctuating velocity $\bm{u}^{\prime}$, the random advection of the streamwise velocity fluctuation $u^{\prime}$ may be considered approximately constant in space and time. As a result, the linearized dynamics (III.1) becomes a finite-dimensional state-space representation in which the velocity fluctuations at any $x$ is given by the state vector $\Psi=[u^{\prime},\omega^{\prime}_{x}]$ such that

	$\displaystyle\frac{\partial\Psi}{\partial t}$	$\displaystyle=$	$\displaystyle\mathcal{A}\Psi+\mathcal{B}\bm{r}_{u},$		(10)
	$\displaystyle\bm{u}^{\prime}$	$\displaystyle=$	$\displaystyle\mathcal{C}\Psi.$

Writing the linearized system (III.1) in the form of Eq (10) has the benefit of drawing knowledge from the area of control theory for dynamical systems when numerical models that are based on NS equations aim to capture the dynamics and statistical features of turbulent flows.

Using the Green’s identity, we get a divergence-free velocity field induced by a concentrated distribution of vortices with arbitrary orientation by solving the Poisson equation

$$\nabla^{2}\bm{u}^{\prime}=-\bm{\nabla}\times\bm{\omega}^{\prime}.$$

The linearized dynamics (III.1) with stochastic forcing provide a velocity field induced by an aligned vorticity field such that

$$\bm{u}^{\prime}(\bm{x},t)=-\frac{1}{4\pi}\int_{V}\frac{1}{|\bm{x}-\bm{y}|}\left[\bm{\nabla}^{\prime}\times\bm{\omega}(\bm{y},t)\right]d^{3}\bm{y}.$$

One may choose the volume $V$ to be aligned so that one of its boundary surfaces becomes normal to both the vorticity vector $\bm{\omega}^{\prime}(\bm{x},t)$ and the direction of background mean flow. The fluctuating velocity is then added to the mean background velocity $\bar{\bm{U}}+\bm{u}^{\prime}=[U(z),v^{\prime}(y,z),w^{\prime}(y,z)]$ in a specific region, such as the inflow boundary $x=0$.

In contrast to conventional method of generating synthetic turbulence, (i) our algorithm does not require the solution of full NS equation to re-circulate nonlinear eddy dynamics on embedded inflow domains, (ii) provides analytically known representation of unresolved small-scale instabilities, and (iii) useful for efficient representation of a much wider range of unresolved scales.

For very light wind with an equivalent roughness length $z_{0}=0.11\nu/u_{*}$, the eddy viscosity varies like $\mathcal{O}(z^{3})$ for $z<50\nu/u_{*}$. In conditions with high wind above flexible surface protrusions, such as crops or forests, observations suggest that the roughness length $z_{0}$ can be dependent on the near surface wind speed, such as

$$z_{0}=\alpha_{c}u_{*}^{2}/g,$$

(11)

where $\alpha_{c}$ is referred to as Charnock’s constant Garratt (1992) and $u_{*}$ is the friction velocity. In the roughness sublayer $z\leq z_{*}$, the wind profile also deviates from the logarithmic wind profile given by the Monin-Obukhov similarity theory Garratt (1992); Arthur et al. (2019). Typical values of $z_{*}/z_{0}$ vary between $10$ and about $150$, so that taking $z_{*}/z_{0}=100$ implies that for an open forest ($z_{0}=0.5$ m), the depth of roughness sublayer is $z_{*}\approx 50$ m.

The drag coefficient $C_{d}$ in the roughness sublayer for a neutrally stratified atmospheric boundary layer is

$$C_{d}=[\kappa/\ln(z_{*}/z_{0}-\Psi_{M}(z_{*}/z_{0})]^{2}$$

(12)

where $\Psi_{M}(z/z_{0})$ is an empirical function such that $\Psi_{M}=0$ for $z/z_{*}\geq 1$. In presence of flexible surface protrusions, the surface boundary condition Sullivan, McWilliams, and Moeng (1994)

$$\sqrt{\langle\tau_{13}\rangle^{2}+\langle\tau_{23}\rangle^{2}}+\sqrt{\langle u^{\prime}w^{\prime}\rangle^{2}+\langle v^{\prime}w^{\prime}\rangle^{2}}=u^{2}_{*}$$

is applied. It assumes that atmospheric boundary layer flows take on the character of a shear-driven channel flow in close proximity of a bounding surface. More specifically in the mean energy equation, the exact viscous dissipation $-2\nu\bm{\mathcal{S}}_{ij}\bm{\mathcal{S}}_{ij}$ (note the bold face) must equal the filter-scale energy flux plus the resolved viscous dissipation, $-2\nu\mathcal{S}_{ij}\mathcal{S}_{ij}$; i.e.,

$$\tau_{ij}\mathcal{S}_{ij}-2\nu\mathcal{S}_{ij}\mathcal{S}_{ij}=-2\nu\bm{\mathcal{S}}_{ij}\bm{\mathcal{S}}_{ij}.$$

Thus, a correction of the eddy viscosity in the vicinity of the surface would account for the correct energy flux and a balance of the production and the dissipation.

Following Schumann (1975), the eddy viscosity is corrected by considering the horizontally-averaged subgrid-scale stresses. Thus, in close proximity of Earth’s surface Sullivan, McWilliams, and Moeng (1994), horizontally-averaged stresses are

$$\langle\tau_{i3}\rangle\approx 2\left[\langle\gamma(z)\nu_{\tau}\rangle+\nu_{\hbox{\tiny T}}\right]\frac{\partial\langle u_{i}\rangle}{\partial z},\quad i=1,\,2,$$

where the damping factor $\gamma(z)$ switches off to Reynolds-averaged Navier-Stokes (RANS) model with a constant eddy viscosity $\nu_{\hbox{\tiny T}}$. However, as it is pointed out by Tennekes and Lumley (1976), the small-scale motions in shear layers will roughly be independent of the mean strain rate, in particular, $\partial\langle u_{i}\rangle/\partial z$. To deal with small-scale motions through vortex stretching, this work has relaxed the horizontal averaging and calculate $\nu_{\hbox{\tiny T}}$ in correcting the eddy viscosity. Thus, instead of the above formulation, we consider

$$\left[\tau_{13}^{2}+\tau_{23}^{2}\right]^{1/4}=C^{1/2}_{d}u(x,y,z_{1})$$

and find the RANS contribution $\nu_{\hbox{\tiny T}}$ from

$$\tau_{i3}=2\left[\nu_{\tau}(x,y,z_{1})+\nu_{\hbox{\tiny T}}\right]\frac{\partial u_{i}}{\partial z},\quad i=1,\,2.$$

Here, $z_{1}$ is the vertical level of grid points adjacent to the bounding surface. Recent work Bhuiyan and Alam (2018, 2020a, 2020b) studied such a profile of near-surface turbulence in atmospheric boundary layer flow over hilly terrain. They observed an excellent agreement between the results of the scale-adaptive LES and the corresponding measurements. Below, we study the $k^{-5/3}$ scaling result in the context of the large scale influence of atmospheric turbulence.

Consider a neutrally stratified atmospheric boundary layer flow in the domain $7056~{}\rm{m}\times 3024~{}\rm{m}\times 756~{}\rm{m}$, which is discretized with $896\times 384\times 96$ grid points (i.e. about $33$ millions), where the isotropic grid spacing is $\Delta x=\Delta y=\Delta z=7.875$ m ($\equiv\Delta$). Using a length scale of $D=126$ m, the computational domain is $[-8D,48D]\times[0,24D]\times[0,6D]$. The flow was initialized with an undisturbed mean profile $\bar{\bm{U}}=[(u_{*}/\kappa)\log(z/z_{0}),0,0]^{T}$.

In an accompanying domain $x\leq-8D$, we employed Eq (III.1) to simulate an ensemble of eddies – attached to the wall (i.e. Earth’s surface) – whose centers are randomly but uniformly distributed in space, and whose axes of rotation are all aligned in the $x$-direction.

In the absence of wind turbines or any other disturbances, Fig 3 shows the transition of initially quiescent flow to a turbulent state. Indeed, the visualization is consistent with the findings of Townsend (1956) that eddies that are attached to a wall can maintain an enhanced level of variance by extracting energy from the background flow and passing it to the perturbation field. The attached eddies, injected by the linear model, are stretched by the logarithmic mean profile in the domain, where the mean strain rate tensor – normalized by the angular velocity $2(u_{*}z_{0})/(\kappa z)$ – is

$$\tilde{\mathcal{S}}=\left[\begin{array}[]{ccc}0&0&1\\ 0&0&0\\ 1&0&0\\ \end{array}\right]$$

Two eigenvectors of $\tilde{\mathcal{S}}$ corresponding to nonzero eigenvalues $\lambda_{1}=1$ and $\lambda_{2}=-1$ are $[1/\sqrt{2},\,0,\,1/\sqrt{2}]^{T}$ and $[1/\sqrt{2},\,0,\,-1/\sqrt{2}]^{T}$, respectively.

To critically analyze the transitional flow development, this study has sampled velocity signals at $41$ locations on the horizontal plane $z=90$ m, as indicated in Fig 4 (with T1-T41). The velocity signals were analyzed to understand the energy content of the eddying motion. In this simulation of the atmospheric boundary layer. For this analysis, consider the auto spectral density function (ASDF) defined by the Fourier transform of the autocorrelation

$$\mathcal{R}_{uu}(t^{\prime})=\frac{1}{T}\int_{0}^{T}u(t)u(t+t^{\prime})dt$$

such as

$$\mathcal{S}_{uu}(f)=\int_{-\infty}^{\infty}\mathcal{R}_{uu}(t^{\prime})e^{-2\pi ift^{\prime}}dt^{\prime}.$$

Since $\mathcal{R}_{uu}(t^{\prime})$ is real and symmetric, setting $t^{\prime}=0$ and applying inverse Fourier transform, we get

$$\overline{u^{2}}=\int_{\infty}^{\infty}\mathcal{S}_{uu}(f)df,$$

which indicates that $\mathcal{S}_{uu}(f)df$ represents the turbulence kinetic energy in a frequency band of $df$.

Fig 5-6 show auto spectral density functions ($S_{uu},\,S_{vv},\,S_{ww}$) for each of the three components of the velocity. We have compared $3$ velocity signals (T1, T3, T5) on the row indicated by ‘T1’ in Fig 5. and ‘T41’, respectively; see in Fig 4. A time-scale $L_{x}/U$ was considered to rescale the frequency, where the length $L_{x}=56D$ accounts for $56$ rotor diameter and the velocity scale $U$ corresponds to the upstream wind at hub-height. From Townsend’s model of attached eddies, larger-than-inertial-scale coherent eddies can sense the effects of ground boundary as they are attached up to a displacement height from the surface Townsend (1956). These attached eddies are anisotropic and “active” in production of shear through the interaction between the subgrid scale momentum flux ($-\overline{u^{\prime}w^{\prime}}$) and the spanwise vorticity ($\partial U/\partial z$) of the mean flow. A comparison of the spectra with $k^{-2/3}$ indicates that attached eddies are in the range from Taylor microscale to approximately the height of the atmospheric boundary layer. A comparison of the spectra with $k^{-5/3}$ indicates a central tenet of the attached-eddy theory that the existence of a self-similar inertial subrange, where large eddies of size $D$ are detached from the ground surface (i.e. cannot sense the effect of the ground boundary).

In the atmosphere, and elsewhere, turbulent flows are characterized by apparently random and chaotic motions. The response to the nonlinear model from the linearized equation are not entirely unorganized in space and time. If the viscous terms were eliminated from the momentum principle, the resulting restricted Euler dynamics can lead to a singularity in finite time Kang, Yun, and Protas (2020). We know that the growth rate of the velocity gradient magnitude $||\bm{\nabla}\bm{u}||$ is driven by $-\rm{Tr}(\mathcal{S}\cdot\mathcal{S}\cdot\mathcal{S})+(1/4)\bm{\omega}^{T}\cdot\mathcal{S}\cdot\bm{\omega}$. The growth of velocity gradient is opposed by the coherence of vortices as they align with the strain rate eigenframe. During such interactions, rapid rotation of vortices result in a turbulence stress on the surrounding strain field, and thus, the strain field must spend energy to oppose such stresses Eyink (2006); Hossen, Mulayath Variyath, and Alam (2021); Johnson (2021). The resulting autonomous dynamics of individual vortices display many traits that are observed in turbulent flows. Other than the relative alignment of vorticity vectors with the strain rate eigenframe, dynamical mechanisms that could advert the formation of singularities due to the growth of velocity gradient tensor are subgrid scale procecess and are not accessible in the LES methods. Recall the consideration of vortex stretching and alignment in Sec II.

The coherent structure of the atmospheric surface layer is characterized by ejections (of low-momentum eddies upward) and sweeps (of high-momentum eddies downward) Garratt (1992). In this section, we briefly study such a coherent structure using the proper orthogonal decomposition (POD) method. POD of complex turbulent flow fields extracts a low-dimensional coherent structure, which is a modal analysis methodology that is extensively used in data science (albeit under different acronyms). Briefly, a truncated decomposition of $u(x,t)$ with space-time separation is

$$u(\bm{x},t)=\sum_{k=1}^{r}c_{k}(t)\phi_{k}(\bm{x}).$$

(16)

A goal of the POD analysis is to identify a data-driven basis $\{\phi_{k}\}_{k=1}^{r}$ which represents the flow kinematics, while the time evolution of the coefficients $\{c_{k}(t)\}_{k=1}^{r}$ captures the low-dimensional coherent dynamics. The functions $\phi_{k}(\bm{x})$ are orthogonal, $\langle\phi_{i},\phi_{j}\rangle=\delta_{ij}$, and they are called POD modes. The coefficients $c_{k}(t)$ are called principal directions of the time dynamics. The first coefficient vector $c_{1}(t)$ of the corresponding POD mode $\phi_{1}(\bm{x})$ – when sampled at discrete times – represents the direction along which the temporal dynamics of turbulence fluctuates the most. The projection of $u(\bm{x},t)$ along the direction of $c_{2}(t)$ to be uncorrelated with that along $c_{1}(t)$ is equivalent to constrain $c_{1}(t)$ to be orthogonal to $c_{2}(t)$. Thus, the principal components are linearly uncorrelated, $\langle c_{i}(t),c_{j}(t)\rangle=\sigma_{i}^{2}\delta_{ij}$, where $\sigma^{2}_{i}$ denotes the variability (or energy) of the signals decomposed by the POD method.

Consider a set of velocity signals $\{u(\bm{x}_{d},t)\}$ corresponding to $n$ turbines, $d=1,\ldots,n$. Here, we investigate the POD modes of the velocity fluctuation $u^{\prime}(\bm{x}_{d},t)$, thereby leading to a set of surrogate signals $\{c_{k}(t)\}$, $k=1,\dots,r\leq m$. The decomposition is carried out by solving the following optimization problem:

$$\operatorname*{\hbox{max}}_{\forall k,c_{k}(t)}\left[\frac{1}{n}\sum_{d=1}^{n}\left\langle c_{k}(t),u^{\prime}(\bm{x}_{d},t)\right\rangle^{2}\right]\quad\hbox{s. t.}\quad\langle c_{i},c_{j}\rangle=\delta_{ij},$$

which is equivalent to the singular value decomposition (SVD) of the $n\times m$ snapshot of the velocity signals.

We have projected the fluctuating part of the aerodynamic power $P^{\prime}_{d}(t)$, defined by Eq (15), along the principal components of the POD modes. Fig 9$a$ shows scattered plots of four pairs of the coefficients of POD modes: $[c_{1}(t),c_{2}(t)]$, $[c_{1}(t),c_{3}(t)]$, $[c_{3}(t),c_{4}(t)]$, $[c_{3}(t),c_{5}(t)]$. There is almost zero correlations among the POD coefficients, as expected. Fig 9$b$ shows the principal components. We see that the first $10$ coefficients capture about $50$% of the turbulence kinetic energy.

Fig 10 show the spectrum of wind power fluctuations for the five wind turbines on the last row of the wind farm. Here, we reconstruct the power fluctuations using Eq (16) for $r=1$, $r=3$, and $r=5$. Fig 10$a$ indicates that for $r=1$ the most variability of wind power fluctuations along the first principal direction is highly sensitive to self-attenuation of vortex stretching. As the number of POD modes $r$ increases, a trend in the convergence to $k^{-5/3}$ is observed. The spectrum of the wind power fluctuation does not vary in the presence of high non-linearity and vortex stretching in wind turbine wakes.