the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Breakdown of the velocity and turbulence in the wake of a wind turbine – Part 2: Analytical modelling
Frédéric Blondel
Valéry Masson
This work aims to develop an analytical model for the streamwise velocity and turbulence in the wake of a wind turbine where the expansion and the meandering are taken into account independently. The velocity and turbulence breakdown equations presented in the companion paper are simplified and resolved analytically, using shape functions chosen in the moving frame of reference. This methodology allows us to propose a physically based model for the added turbulence and thus to have a better interpretation of the physical phenomena at stake, in particular when it comes to wakes in a nonneutral atmosphere. Five input parameters are used: the widths (in vertical and horizontal directions) of the nonmeandering wake, the standard deviation of wake meandering (in both directions) and a modified mixing length. Two calibrations for these parameters are proposed: one if the users have access to velocity time series and the other if they do not. The results are tested on a neutral and an unstable largeeddy simulation (LES) that were both computed with MesoNH. The model shows good results for the streamwise velocity in both directions and can accurately predict modifications due to atmospheric instability. For the axial turbulence, the model misses the maximum turbulence at the top tip in the neutral case, and the proposed calibrations lead to an overestimation in the unstable case. However, the model shows encouraging behaviour as it can predict a modification of the shape function (from bimodal to unimodal) as instability and thus meandering increases.
 Article
(6288 KB)  Fulltext XML
 BibTeX
 EndNote
The CPU cost of classical computational fluid dynamic models is too high to deal with all the different cases needed to estimate and optimize the performances of a wind farm. Thus, socalled engineering models have been developed to estimate the power loss due to wakes at a low computational cost (e.g. Jensen, 1983; Larsen et al., 2008; Bastankhah and PortéAgel, 2014). These design tools are based on physical considerations and are often calibrated and validated against numerical results or measurements. Among these tools, analytical models are the quickest: they consist of a single formula that can be directly applied to the wind farm setup and atmospheric conditions, leading to fast results even for a whole farm. A very commonly used model is the one developed by Bastankhah and PortéAgel (2014), who assumed an axisymmetric and selfsimilar Gaussian velocity deficit in the wake and solved the mass and momentum conservation equations to find a relation between the amplitude and width of the Gaussian. It can be adapted for a nonaxisymmetric wake (Xie and Archer, 2014):
where $\stackrel{\mathrm{\u203e}}{U}$ is the mean velocity field; ${\stackrel{\mathrm{\u203e}}{U}}_{\mathrm{\infty}}$ is the mean velocity upstream of the turbine; C(x) is the maximum velocity deficit; C_{T} is the thrust coefficient; D is the turbine diameter; ($x,y,z$) denotes the streamwise, lateral and vertical coordinates, centred at the turbine's hub; and σ_{y,z} denotes the wake widths in the lateral and vertical directions. In the present work, the vertical and horizontal axes are centred at the hub position. Here and in the following, the Reynolds decomposition is used to write any unsteady field X(t) as a sum of a mean and a varying part: $X\left(t\right)=\stackrel{\mathrm{\u203e}}{X}+{X}^{\prime}\left(t\right)$.
The stability of the atmospheric boundary layer (ABL) influences the wake recovery (Abkar and PortéAgel, 2015), and the largescale eddies carried in this region of the atmosphere are often associated with wake meandering, i.e. oscillations in the instantaneous wake around its mean position (Larsen et al., 2008). To model the meandering, the concepts of fixed and moving frames of reference (denoted FFOR and MFOR respectively) defined in the dynamic wake meandering (DWM) model are used herein (Larsen et al., 2007). The FFOR is bound to the ground: it is the frame of reference in which we want to compute the turbulence and velocity fields. In the FFOR the effects of meandering are not differentiated from the wake expansion caused by turbulent mixing. The MFOR moves with the wake centre at each time step: in this frame of reference, only the wake expansion due to turbulent mixing is represented, making the fields in this frame of reference easier to interpret. The instantaneous streamwise velocity can be changed from one frame to another according to the following relation:
where subscripts MF and FF denote the velocity fields in the MFOR and FFOR respectively, and y_{c}(x,t) and z_{c}(x,t) are the time series of the wake centre's coordinates at the downstream position x. The concept of MFOR and FFOR can be used to write an analytical wake model for the velocity deficit as in the work of Braunbehrens and Segalini (2019):
where σ_{fy,fz}(x) denotes the standard deviations of the wake centre's coordinates in the lateral and vertical directions respectively, σ_{y,z}(x) denotes the wake widths in the MFOR, and C(x) is the maximum velocity deficit in the MFOR. Such a model allows for the independent calibration of the effects of meandering (through the variables σ_{fy,fz}) and of wake expansion due to turbulent mixing (through the variables σ_{y,z}). The former parameters are a function of atmospheric stability through lateral and vertical turbulence (Braunbehrens and Segalini, 2019; Du et al., 2021; Brugger et al., 2022), whereas the latter parameters can be a function of axial turbulence as in Eq. (1) (Fuertes et al., 2018; Niayifar and PortéAgel, 2016) or turbine operating conditions such as C_{T} and atmospheric shear (Braunbehrens and Segalini, 2019).
For the turbulent kinetic energy (TKE), it is common to model only the maximum value of added turbulence, which can be computed with the Crespo model (Crespo and Hernandez, 1996) or the Frandsen model (Frandsen, 2007) as in the IEC 614001 standard. Their approach is mainly empirical and can be extended to describe the whole profile of turbulence instead of the maximum value alone (Ishihara and Qian, 2018). This widely used model (hereafter denoted I&Q2018) is simple, since it only requires the knowledge of the thrust coefficient and the upstream turbulence intensity, but it is totally empirical and does not account for atmospheric stability. Moreover, it has been shown that the wake in an unstable ABL dissipates faster than in a neutral ABL, even at the same level of turbulence intensity (Du et al., 2021). This behaviour cannot be taken into account in the I&Q2018 model due to the limited number of inputs.
The present work aims to propose a physically based model that predicts both the mean and the variance (i.e. turbulence) of the axial velocity in the wake of a wind turbine. The advantage of basing our model on physical interpretations is that it gives more room for further improvements, as we know which assumptions were made and how it degrades the results. Moreover, the proposed model is dependent on atmospheric stability, since it influences both the velocity and the turbulence fields in the wake (see companion paper). Many models, such as the I&Q2018 model, do not take atmospheric stability into account, assuming that stable and unstable cases compensate for each other and that a calibration on neutral cases is thus sufficient. This approach is valid for monthly or yearly estimations of wind farms' performances. But some applications of the future wind industry, such as digital twins, need estimations over a day, an hour or even smaller periods. In such cases, the stability must be taken into account. Since we showed in the companion paper that stability mainly affects the wake meandering, this phenomenon must be decoupled from the wake expansion to take the ABL stability into account. To do so, the breakdowns described in the companion paper are reused and briefly explained in the following.
A field in the MFOR can be written as an unsteady translation of the same field in the FFOR through Eq. (3). To shorten this equation, the notation $\widehat{a(y,z)}=a(y{y}_{\mathrm{c}}(t),z{z}_{\mathrm{c}}(t\left)\right)$ for any field a is used. For the present work, it is also important to note that for any field a,
where $**$ denotes a 2D convolution and f_{c} is the probability density function (PDF) of the wake centre position. In the companion paper, it was shown that the velocity (Eq. 6) and turbulence (Eq. 7) in the FFOR can be expressed as a function of their counterparts in the MFOR. This is achieved by dividing these quantities into several terms, (I) and (II) in Eq. (6) and (III) to (VII) in Eq. (7).
These terms are thoroughly described and quantified in the companion paper, where they are separated into pure terms (I, III and IV) and cross terms (II, V, VI and VII).
Term (I) is the convolution of $\stackrel{\mathrm{\u203e}}{{U}_{\mathrm{MF}}}$ with f_{c}. It is a pure mean velocity term: it is 0 only if the mean velocity is 0. Conversely, term (II) is a cross term because it can be equal to 0 either if there is no meandering (operator $\widehat{.}$ has no effect) or if there is no turbulence in the MFOR (${U}_{\mathrm{MF}}^{\prime}=\mathrm{0}$). Term (III), also written as k_{m} in the following to be consistent with notation from Keck et al. (2013) and Conti et al. (2021), is the turbulence purely induced by meandering: in the case of a meandering steady wake, i.e. ${U}_{\mathrm{MF}}^{\prime}=\mathrm{0}$, Eq. (7) reduces to this term only. Term (IV) is the rotoradded turbulence, which is also written as k_{a} for consistency with other works (Conti et al., 2021). It is the turbulence purely induced by the rotor: in the absence of meandering, the equation reduces to this term only, also written as k_{a} in the following for consistency with the literature. Term (V) is the covariance of $\widehat{\stackrel{\mathrm{\u203e}}{{U}_{\mathrm{MF}}}}$ and $\widehat{{U}_{\mathrm{MF}}^{\prime}}$, term (VI) can be viewed as the varying part of the MFOR turbulence, and term (VII) is the square of term (II). It is a pure dissipation term as it is always negative. Like term (II), they are cross terms since they are equal to zero if either the turbulence in the MFOR or the meandering is null. The companion paper showed that terms (II) and (VII) are negligible in their respective equations. In the breakdown of the turbulence equation, term (V) is of lesser importance than (III) and (IV) but drives the vertical asymmetry of the turbulence profiles.
The proposed analytical model is based on the velocity and turbulence breakdowns (Eqs. 6 and 7). Similar to Eq. (4) (Braunbehrens and Segalini, 2019), the reasoning starts by writing the wake properties in the MFOR and the wake meandering with different parameters to take into account meandering due to atmospheric stability independently of the expansion due to turbulence mixing. It is common in wake modelling to assume that meandering can be entirely accounted for by increasing the wake expansion. However, it is a phenomenon of a different nature, and it leads to velocity and turbulence profiles of different shapes. In the present model, these phenomena are modelled separately, and it is assumed that they do not interact. This is equivalent to neglecting cross terms in Eqs. (6) and (7), which have been shown to take smaller values than pure terms in the companion paper. However, in the future, modelling these cross terms might be necessary to improve the results. The main added value of this work is to propose a new framework that can be used with different shape functions in the MFOR to propose other turbulence models. Nevertheless, two calibrations (one requiring the inflow time series and one that does not) are proposed for the model to demonstrate how it can be tuned and to test the model.
In the second section of this work, the datasets are presented: for the calibration of the model, a dataset from the MOMENTA project is used, and for the validation the neutral and unstable cases obtained from the largeeddy simulations (LESs) from the companion paper are reused. The third section presents the derivation of the model. The fourth section shows the chosen calibration methods, and the fifth section presents the corresponding results. All these results are discussed in the sixth section, followed by the conclusion.
2.1 Description of the LES code
The analytical model developed in this work is based on LES datasets generated with the MesoNH solver (Lac et al., 2018). It is a finitevolume and finitedifference research code for ABL simulations where the Navier–Stokes equations and the energy conservation equation are resolved on an Arakawa C grid. This solver models the stability of the ABL with a buoyancy term in the momentum equation. The Coriolis force and largescale forcing are also taken into account. The effect of the wind turbine on the surrounding flow is modelled with an actuator line method, i.e. rotating source terms in the momentum equation.
To close the set of equations, the subgrid TKE equation is resolved, and all the subgrid quantities are written as a function of this subgrid TKE, the resolved variables and a Deardorff mixing length. A grid nesting method allows for simultaneously having a vertical and horizontal mesh size of 1.5 and 0.5 m in the wake region for the two datasets and a domain large enough to compute the largest eddies of the atmosphere. The model and numerical parameters are described in more detail in the companion paper.
2.2 Simulation setup
Two different LES datasets are used in this work: the first one for creating and calibrating the model and the second one for testing the model's results. Inflow conditions of these datasets can be found in Table 1. For both datasets, only the wake mean streamwise velocity (${\stackrel{\mathrm{\u203e}}{U}}_{x}$, written as U_{x} in the following) and the streamwise turbulence (${k}_{x}=\stackrel{\mathrm{\u203e}}{{u}^{\prime}{u}^{\prime}}$) are computed. The proposed model thus only deals with the streamwise velocity and turbulence.
The calibration dataset contains six simulations, with four different ABL stabilities and three different thrust values. The simulated turbine is 92 m in diameter and has a hub height of 80 m. The turbine's data were obtained in the context of the MOMENTA project (Jézéquel, 2023).
To perform such simulations, a precursor without heat flux is first simulated in a domain of 19 km × 15 km (with a horizontal resolution of 37.5 m) for 25 h to let the turbulence establish itself and to allow the system to reach a quasisteady state. Then, a ground heat flux is applied for 4 h: 0, 30, 60 and 120 W m^{−2} for “neutral”, “weakly unstable”, “unstable” and “strongly unstable” cases respectively. This allows us to simulate three different levels of atmospheric instability, starting from the same neutral state. No stable case was simulated because of the induced veer (gradient of inflow wind direction) that leads to a deformed wake. The veer could have been modelled as in Abkar et al. (2018), but it would significantly complicate the present derivations. Moreover, meandering and meandering turbulence are negligible in a stably stratified ABL (see companion paper), and there is thus little interest in using the approach presented herein. Developing the model for veered cases is a challenge that is out of the scope of this work.
After these two steps, the coarsest computational domain (horizontal resolution of 37.5 m) is ready: two grid nestings are then applied to reach a resolution of 1.5 m in the most refined domain. Then, 10 min of dynamics is used to let the flow establish itself in the wake of the wind turbine, and the postprocessing is performed on the following 50 min of dynamics. The data are sampled at 1.2 Hz, which is the approximate limit between resolved and subgrid TKEs for these simulations (equivalent to 4 times the mesh size).
The wind turbine rotational speed and pitch are set according to the controller's database and the calculated upstream velocity. Since all the cases are computed at a similar inflow velocity, similar values of the thrust coefficient are obtained in the simulations. To include the influence of the thrust coefficient in the model, two additional cases with a degraded thrust coefficient are also computed, with the same inflow as the neutral case. To reduce the thrust, the pitch value is increased from 0 to 3 and 4.5^{∘} respectively.
The second set of simulations, hereafter called the validation dataset, is based on the neutral and unstable cases that are described in the companion paper. The simulated turbine is a modified version of the Vestas V27: it is a threebladed rotor with a diameter D=27 m and a hub height of 32.1 m. The simulation methodology is quite similar, as described in the paragraph above, except that one additional nesting is required to reach the targeted mesh size. In the validation dataset, the velocity is sampled at 1 Hz, and the simulations last for 80 and 40 min for the neutral and unstable cases respectively. This was due to benchmark requirements and computational limitations. A statistical convergence of our datasets is proposed in the Appendix of the companion paper. Overall, it concluded that increasing the duration of the simulation for the unstable case would improve the reliability of the simulations.
2.3 Wake tracking
For the validation simulations, the wake centre's y_{c}(x,t) and z_{c}(x,t) coordinates are computed at each time step and each downstream position with the constantflux waketracking algorithm, which is described in the companion paper. To facilitate the wake tracking, a reference simulation is also run. It is a simulation with the same inflow and boundary conditions but without the wind turbine. The corresponding velocity field noted as U_{ref} thus represents a developing ABL without the perturbations of a wind turbine.
Another method is proposed here to compute the unsteady wake centres in the calibration dataset. A passive scalar (similar to a pollutant) is emitted at the rotor disc with a concentration value of 1 at each time step. This new variable is only driven by the advection scheme, in accordance with the passive tracer of the DWM theory, and only marginally impairs the code's performance. By supposing that this variable follows the wake, the unsteady wake centre is deduced from the centre of mass of this pollutant at each downstream position. The results lead to a lowfrequency behaviour similar to the constantflux method used in the companion paper but with fewer outliers (see Fig. 1). Since the postprocess is more straightforward and the results seem better, this method has been used for the calibration dataset.
2.4 Inflow conditions
Table 1 shows the hub height velocity, thrust coefficients and turbulence intensities at hub height for each of the cases. The directional turbulence intensities are defined as
and the global turbulence intensity is defined as
Figures 2 and 3 show the profiles of some inflow variables for the calibration and validation cases respectively. The profiles are taken 2.5 diameters upstream of the wind turbine and are averaged along the y direction (the direction transverse to the wind turbine).
In the left panel the mean velocity is plotted. The calibration dataset (Fig. 2) was built in order to have similar hub height velocities between the cases (around 7 m s^{−1}), whereas the validation dataset comes from simulations that reproduced the SWiFT benchmark where the hub height velocities differed. The Monin–Obukhov profiles are plotted as follows using dotted lines:
where κ=0.41 is the von Kármán constant (Cheng et al., 2019),
and ${x}_{u}=(\mathrm{1}\mathrm{15}\cdot z/{L}_{\mathrm{MO}}{)}^{\mathrm{0.25}}$. Since z_{0} is known from the simulations (0.17 in the calibration dataset and 0.014 in the validation dataset), the profiles are found by fitting Eq. (10) on the corresponding velocity profile, with parameters u_{*} and L_{MO}. The results, in dotted lines, match the inflow profiles well, showing that it respects the Monin–Obukhov similarity theory around the turbine's height.
The middle panels of Figs. 2 and 3 show the inflow TKE, defined as
In the calibration dataset, the amount of TKE increases as the imposed heat flux increases. In the validation dataset, this is not the case since the neutral case is at a higher velocity at hub height, but the TI of the unstable case is indeed higher than that of the neutral case.
The right panels of Figs. 2 and 3 show the modified mixing length ${l}_{\mathrm{m},\mathrm{\infty}}^{*}$ upstream of the wind turbine. This quantity is discussed and used in Sect. 3 to compute the mixing length in the MFOR. Here, the value is computed as the ratio of turbulence and shear:
However, in the unstable cases, the velocity profile becomes nearly constant above a given height, leading to low values of $\partial {U}_{\mathrm{\infty}}/\partial z$ and thus very chaotic behaviour of ${l}_{\mathrm{m},\mathrm{\infty}}^{*}$. To have a more reliable curve, the derivative of U is resolved analytically using Eq. (10):
with L_{MO} and u_{*} fitted from the velocity profile. The resulting curve, in dotted lines, gives a more useful quantity in the turbulence to shear ratio, while still being of the order of magnitude of the directly computed ratio (in solid lines).
In this section, we derive an analytical model for the dominating terms of Eqs. (6) and (7). First, an analytical form is proposed for the velocity deficit in the MFOR $\mathrm{\Delta}{\stackrel{\mathrm{\u203e}}{U}}_{x,\mathrm{MF}}$, the turbulence in the MFOR k_{x,MF} and the meandering distribution f_{c}. Then, some terms are neglected, and the convolutions of Eqs. (6) and (7) are resolved analytically to get a model for the velocity deficit and added turbulence in the FFOR. To help the reader, the main variable notations and subscripts used in this section and afterwards are summarized in Table 2.
3.1 Independent modelling of the wake in the MFOR and meandering
3.1.1 Wake velocity deficit in the MFOR
Based on the literature (Bastankhah and PortéAgel, 2014; Xie and Archer, 2014), the mean velocity deficit is modelled with the longestablished Gaussian velocity deficit (see Eq. 1):
where subscript ._{am} stands for the analytical model; C(x) is defined in Eq. (2); and σ_{y}, σ_{z} denotes the wake widths in the MFOR. The overline is dropped because this analytical model is static. In the literature, it has been shown that doubleGaussian (Keane et al., 2016) or superGaussian (Blondel and Cathelain, 2020) shapes provide more accurate results, but here the Gaussian shape allows for a straightforward computation of the convolutions in our model and is still pertinent in the far wake. It is shown later that this approximation leads to discrepancies in the near wake.
3.1.2 Wakeadded turbulence in the MFOR
To model term (IV) or k_{x,a}, one needs an analytical form for the turbulence in the MFOR k_{x,MF}. It was first thought that it would be better to model the added turbulence in the MFOR, i.e. $\mathrm{\Delta}{k}_{x,\mathrm{MF}}={k}_{x,\mathrm{MF}}{k}_{x,\mathrm{\infty}}$, in order to separate the rotoradded turbulence Δk_{x,MF} from the ambient turbulence. This procedure was done in the companion paper; however it leads to negative values of Δk_{x,MF} (in particular near the ground), i.e. smaller turbulence in the wake compared to the turbulence upstream of the wind turbine. This is not compatible with a model that predicts only increased turbulence in the wake of a turbine (as is the case here or in I&Q2018), and this approach has thus been abandoned.
The derivation of a model for k_{x,MF} is not as straightforward as for ΔU_{MF} because turbulence comes from the unsteadiness of the flow, whereas an analytical model is by definition steady. In the DWM, the Madsen formulation (Madsen et al., 2010) is used to scale the velocity profile with an empirical function of the wakegenerated shear. One could also assume selfsimilarity of the Δk_{x,MF} function and try to derive a model as was done for the velocity in Bastankhah and PortéAgel (2014). The main issue here is that an analytical form of the model is needed in the FFOR; i.e. the convolution of f_{c,am} with the chosen shape function for $\mathrm{\Delta}{k}_{x,\mathrm{MF},\mathrm{am}}$ must have an analytical solution, which is not trivial for the aforementioned models.
It is proposed here to assume that the turbulence in the MFOR is solely driven by wakegenerated shear, as in Madsen et al. (2010). To relate the turbulence in the MFOR to mean gradients, two models for the velocity scale u_{0} are combined. In the first, it is assumed to be proportional to the square root of the TKE (Pope, 2000). However in the present work, the 3D TKE is not computed, so it is replaced with the axial turbulence k_{x}:
where C_{μ} is a constant and l_{m} is the mixing length. In the second method, the velocity scale is defined from the norm of the strainrate tensor $\mathrm{}\stackrel{\mathrm{\u203e}}{\stackrel{\mathrm{\u203e}}{S}}\mathrm{}$
From the literature (Iungo et al., 2017), it appears that in the wake of a wind turbine, the dominating term (in cylindrical coordinates) is $\frac{\partial U}{\partial r}$. It is supposed herein that these results can be transposed in Cartesian coordinates and are applicable in the MFOR. The velocity scale can thus be written as a function of the derivatives of the axial velocity.
Combining Eqs. (16) and (18) leads to
In Eq. (19), the last term $(\mathrm{1}\mathrm{\Delta}{U}_{\mathrm{\infty}})\frac{\partial {U}_{\mathrm{\infty}}\left(z\right)}{\partial z}$ represents the produced turbulence due to the interaction between wakegenerated shear and atmospheric shear. It is this term that induces a maximum of turbulence at the top tip in cases of high atmospheric shear such as neutral or stable ABLs. Even though an analytical form of this term can be found by assuming U_{∞}(z) to be a log law or a power law, the convolution product with f_{c} in Eq. (7) did not lead to any analytical solution.
It was thus decided to neglect shear in the formulation and to add the contribution of the inflow turbulence with a maximum function. This is a strong assumption that impacts the results (see Sect. 5) but allows us to compute the total added turbulence,
with
where ${K}_{\mathrm{MF}}=({U}_{\mathrm{\infty}}C{l}_{\mathrm{m}}^{*}{)}^{\mathrm{2}}$ and ${l}_{\mathrm{m}}^{*}$ is the modified mixing length ${l}_{\mathrm{m}}^{*}={l}_{\mathrm{m}}/\sqrt{\mathrm{2}}{C}_{\mathit{\mu}}^{\mathrm{1}/\mathrm{4}}$. In other words, the modified mixing length is the ratio of the axial turbulence to the quadratic sum of the vertical and horizontal gradients of the axial velocity deficit.
Figure 4 shows the profiles of the modified mixing length in the wake normalized by the modified mixing length upstream of the turbine at hub height: ${l}_{\mathrm{m}}^{*}/{l}_{\mathrm{m},\mathrm{\infty}}^{*}\left({z}_{\mathrm{hub}}\right)$, where ${l}_{\mathrm{m},\mathrm{\infty}}^{*}$ is defined in Eq. (13) and ${l}_{\mathrm{m}}^{*}$ is computed as
One can see that there are two distinct values: one inside the wake and one outside the wake. Inside the wake, the value is fairly constant (except in the bottom of the wake where it increases chaotically, probably due to the effect of the ground). It only seems to vary with the streamwise distance, and it was thus assumed that ${l}_{\mathrm{m}}^{*}/{l}_{\mathrm{m},\mathrm{\infty}}^{*}\left({z}_{\mathrm{hub}}\right)$ is only dependent on x. Theoretically, it could be possible to develop a model with two mixing lengths (one for the wake and another for the ambient turbulence), but with such an assumption, no analytical solution of Eq. (7) could be achieved.
Note that in Eq. (21), the error in the near wake due to the Gaussian shape assumption for velocity deficit in the MFOR propagates onto $\mathrm{\Delta}{k}_{x,\mathrm{MF},\mathrm{am}}$. Using a Gaussian instead of a superGaussian function leads to an underestimation of the wakegenerated shear and thus to a much weaker but more spread axial turbulence around the blade's tips. Moreover, the model does not account for the atmosphericgenerated shear. This phenomenon, which leads to a smaller value of wakegenerated turbulence at the bottom tip compared to the top tip, cannot be represented in our model. Finally, the model imposes that ${k}_{x,\mathrm{MF},\mathrm{am}}=\mathrm{0}$ at the centre of the wake, a condition that is not fulfilled in the calibration dataset. Another possible improvement would be to add the streamwise gradient $\partial {U}_{x}/\partial x$ in Eq. (18). Despite these flaws, this expression has been chosen since it has an analytical solution of its convolution with the wake centre position distribution f_{c,am}.
3.1.3 Wake meandering
For the PDF of wake meandering, the central limit theorem leads to a Gaussian distribution (Braunbehrens and Segalini, 2019). The distribution of the wake centre f_{c} is nonaxisymmetric, and its variance σ_{f} is thus defined in both dimensions:
3.2 Velocity in the FFOR
In the following, the dependency of the variables on coordinate x is omitted to shorten the equations.
In Eq. (6), the velocity in the wake is written under its dimensional form, whereas the model chosen in Eq. (15) is written under the velocity deficit form. To relate the velocity to the velocity deficit, it is needed to assume that despite its dependency on z due to the atmospheric shear, the upstream velocity U_{∞} can be considered to be a constant when applying the 2D convolution product with the wake centre distribution. For any function g(y,z), this simplification can be written as
An analytical form of term (I) can then be deduced from Eqs. (15) and (23):
Since it has been shown in the companion paper that term (II) of Eq. (6) is negligible, we approximate that ${U}_{x,\mathrm{FF},\mathrm{am}}=(\mathrm{I}{)}_{\mathrm{am}}$. The velocity deficit in the FFOR ΔU_{FF,am} is thus the convolution product of two Gaussian functions. It is known that the convolution product of two normalized Gaussian functions of variance ${\mathit{\sigma}}_{a}^{\mathrm{2}}$ and ${\mathit{\sigma}}_{b}^{\mathrm{2}}$ is a normalized Gaussian function of variance ${\mathit{\sigma}}_{a}^{\mathrm{2}}+{\mathit{\sigma}}_{b}^{\mathrm{2}}$ (Teitelbaum, 2022). Equation (25) can be written as the product of two convolution products, leading to
Even though the reasoning of Braunbehrens and Segalini (2019) is different, it is shown here that their model (Eq. 4) can be found by neglecting term (II) and assuming Eq. (24) as well as Gaussian shapes for the velocity deficit in the MFOR and the wake centre's distribution. This is still a Gaussian form, i.e. Eq. (1), with a FFOR wake width defined as ${\mathit{\sigma}}_{ty,tz}=\sqrt{{\mathit{\sigma}}_{y,z}^{\mathrm{2}}+{\mathit{\sigma}}_{fy,fz}^{\mathrm{2}}}$ and a maximum velocity deficit of
To fulfil the conservation of momentum as in Eq. (2), one would need ${C}_{\mathrm{FF}}=\mathrm{1}\sqrt{\mathrm{1}{C}_{\mathrm{T}}/(\mathrm{8}{\mathit{\sigma}}_{ty}{\mathit{\sigma}}_{tz}/{D}^{\mathrm{2}})}$, which is not the case here. In fact, with this methodology, the conservation of momentum can only be enforced in the MFOR or the FFOR. This is the consequence of neglecting term (II) in the velocity breakdown; however, the induced error is relatively low since term (II) is negligible. Combining Eqs. (25) and (26) leads to our model for the velocity in the wake of a wind turbine:
3.3 Model for the turbulence in the FFOR
For the turbulence, a model has been found for terms (III) (Eq. 30) and (IV) (Eq. 33). Even though the contribution of the three cross terms of Eq. (7) is not always negligible (see companion paper), the two modelled terms are predominant, and the result of the model limited to these two terms can be compared to the turbulence in the FFOR. The total modelled turbulence is computed here as
3.3.1 Meandering term
With the same assumptions as for term (I), it is possible to derive an analytical formulation for term (III) of Eq. (7), i.e. the turbulence induced by wake meandering. The assumption of Eq. (24) must be used again to get ${U}_{\mathrm{\infty}}^{\mathrm{2}}$ out of the convolution product, and Eq. (26) is reused to compute the righthand side of term (III): ${\stackrel{\mathrm{\u203e}}{\widehat{\stackrel{\mathrm{\u203e}}{{U}_{\mathrm{MF}}}}}}^{\mathrm{2}}$. On the lefthand side, there is a convolution of the Gaussian function f_{c,am} with $\mathrm{\Delta}{U}_{x,\mathrm{MF},\mathrm{am}}^{\mathrm{2}}$, which is also a Gaussian function of widths $\sqrt{\mathrm{0.5}}{\mathit{\sigma}}_{y}$ and $\sqrt{\mathrm{0.5}}{\mathit{\sigma}}_{z}$. It is thus possible to use the fact that the convolution of two Gaussian functions is a Gaussian function (Teitelbaum, 2022).
The shape of term (III) is thus not a doubleGaussian, as it may be interpreted in the literature (Stein and Kaltenbach, 2019; Ishihara and Qian, 2018), but rather a Gaussian of width $\sqrt{\mathrm{0.5}{\mathit{\sigma}}^{\mathrm{2}}+{\mathit{\sigma}}_{f}^{\mathrm{2}}}$ minus a thinner and less pronounced Gaussian of width $\sqrt{\mathrm{0.5}{\mathit{\sigma}}^{\mathrm{2}}+\mathrm{0.5}{\mathit{\sigma}}_{f}^{\mathrm{2}}}$. It can be verified that this expression is always larger than 0; i.e. the meandering only produces turbulence and does not dissipate it.
3.3.2 Rotoradded turbulence term
Combining the chosen models for the wake meandering distribution and the added turbulence in the MFOR (Eqs. 21 and 23) in Eq. (20) leads to an analytical form of the axial rotoradded turbulence:
At this point, the added turbulence in the FFOR is the sum of two terms that are identical if the coordinates y and z are swapped. It is the product of two convolutions: the first of $f:y\to {y}^{\mathrm{2}}\mathrm{exp}({y}^{\mathrm{2}}/{\mathit{\sigma}}_{y}^{\mathrm{2}})$ with a Gaussian function and the second of two Gaussian functions. The first convolution product has been solved with a computer algebra tool (Scherfgen, 2022), and the other has already been solved in Eq. (30). It gives
From Eq. (32), we need to add the same quantity with y←z and z←y, then factorize and simplify to deduce the model for ${k}_{x,\mathrm{a},\mathrm{am}}$ as follows:
where
It can be noted that in the absence of meandering, i.e. for ${\mathit{\sigma}}_{fy}={\mathit{\sigma}}_{fz}=\mathrm{0}$, the model retrieves its MFOR form (Eq. 21). As for terms (I) and (III), the expression of ${k}_{x,\mathrm{a},\mathrm{am}}$ is based on a Gaussian velocity deficit hypothesis, even in the near wake where the LES wake exhibits a shape closer to a tophat function. The velocity gradient that is the source of the rotoradded turbulence is thus lower and more spread in the model compared to the actual values. Another issue of the model is that it inadequately takes shear into account due to the assumptions of Eqs. (20) and (24). Indeed, the only source of vertical asymmetry in Eq. (33) is ${U}_{\mathrm{\infty}}^{\mathrm{2}}$, i.e. the velocity shear upstream of the turbine.
The model's equations are based on five variables: the wake widths in the MFORs σ_{y} and σ_{z}, the modified mixing length ${l}_{\mathrm{m}}^{*}$, and the standard deviations of the meandering distributions σ_{fy} and σ_{fz}. Each of these variables needs to be calibrated from the inflow conditions to have a usable model. To do so, the results from the calibration dataset are used. Two versions of the wake meandering calibration are used: the “base” calibration, used when the time series of the upstream velocities are known, and the “engineering” calibration, used when they are not.
4.1 Wake width in the MFOR
As described in Sect. 3, we assume that the wake in the MFOR follows a Gaussian shape function. Moreover, we assume here that the wake is axisymmetric (σ_{y}=σ_{z}), thus reducing the number of parameters in the model from five to four. The width of the wake in the MFOR is deduced from fitting the function of Eq. (35) on the velocity deficit ΔU_{MF} through a nonlinear leastsquares method.
C is fixed as a function of σ (Eq. 2 with ${\mathit{\sigma}}_{y}={\mathit{\sigma}}_{z}=\mathit{\sigma}$), and the optimization is run on parameters $\mathit{\{}{C}_{\mathrm{0}},{y}_{\mathrm{0}},{z}_{\mathrm{0}},\mathit{\sigma}\mathit{\}}$, where y_{0}, z_{0} are the mean wake centres; σ is the wake width (the parameter of interest); and C_{0} is an offset to help the algorithm.
The resulting wake widths in the MFOR as a function of the downstream distance are plotted with solid lines in Fig. 5 for the six cases of the calibration dataset. Except in the near wake, the wake width evolves linearly with the distance to the turbine. Moreover, the greater the instability (and thus the level of turbulence; see Fig. 2), the greater the slope of this linear relation. Finally, the simulations with degraded thrust seem to have the same slope as the neutral case but with a different origin.
For all these reasons, the chosen function for the calibration is the following:
where a, b and c are parameters to fit; I is the total turbulence intensity (Eq. 9); and $\mathit{\beta}=\mathrm{0.5}\left(\mathrm{1}+\sqrt{\mathrm{1}{C}_{t}}\right)/\sqrt{\mathrm{1}{C}_{t}}$ (Bastankhah and PortéAgel, 2014). A leastsquare fit method on the six different curves allowed us to compute the best values of a, b and c (see Table 3). Note that this fit is in the end very similar to what can be found in the literature (e.g. Fuertes et al., 2018), except that the slope (parameter a) is smaller because the models of the literature implicitly assume that the meandering is included in the wake expansion.
4.2 Modified mixing length
The modified mixing length ${l}_{\mathrm{m}}^{*}$ in Eq. (33) directly drives the amount of turbulence added by the turbine. In Sect. 3, it is shown that this variable in the upper part of the wake is independent of the simulation case when normalized with the upstream modified mixing length. Therefore, the evolution of ${l}_{\mathrm{m}}^{*}/{l}_{\mathrm{m},\mathrm{\infty}}^{*}$ is plotted in Fig. 6, and it shows an approximately linear behaviour with the downstream distance.
In the first approach, it is thus decided to fit the mixing length with a linear function of $x/D$:
where ${l}_{\mathrm{m},\mathrm{\infty}}^{*}$ is deduced from Eqs. (13) and (14), in which u_{*} and L_{MO} can be found from a fit of the inflow velocity profile. A leastsquare fit method on the different curves from Fig. 6 is used to fit Eq. (37). The resulting parameters d and e can be found in Table 4, and the corresponding fitted function is plotted with a dashed black line in Fig. 6. The results are quite satisfying even though all the curves are not perfectly superimposed.
4.3 Wake meandering
The widths of the wake centre's distributions, σ_{fy} and σ_{fz}, are computed as the standard deviations of the wake centre's coordinates, y_{c}(x,t) and z_{c}(x,t):
The resulting amount of meandering in the horizontal (top figure) and vertical directions (bottom figure) for the six cases of the calibration dataset can be found in Fig. 7. The LES results are plotted with solid lines. Overall, the more unstable the case, the more meandering is found. However, the meandering does not solely depend on the lateral turbulence intensity. In particular, the weakly unstable case has greater vertical meandering than the unstable case, despite having a lower I_{z} value (see Table 1). It is also worth noting that the reductions in the thrust coefficient have little to no effect on the meandering (all the neutral cases are equivalent).
To model the amount of meandering, Braunbehrens and Segalini (2019) propose the following formula:
where U_{c} is an advection velocity and 𝒜 is the autocorrelation function of the velocity (the lateral and vertical one respectively). For each case, the results of Eq. (39) are plotted with dashed lines in Fig. 7 for U_{c}=0.8U_{∞}. This model works fairly well for the amount of meandering, with the right order of magnitude in each case, and it predicts the different behaviours of the vertical and lateral directions for the unstable and weakly unstable cases. However, such calibration for σ_{fy} and σ_{fz} is not appropriate for analytical wake modelling because it requires time series of wind velocities at hub height, whereas usually only the mean values are available.
Therefore, we hereby propose (dotted lines in Fig. 7) an engineeringoriented solution to approximate the amount of meandering without access to the unsteady time series of velocities upstream of the turbine. In the first attempts to model the meandering (Ainslie, 1988), it was proposed that the wake meandering should be a linear function of the inflow wind direction's variance. However, more recent work (Doubrawa et al., 2018) has shown that the amount of meandering decreases with the rotor size. Indeed, following the theory of the DWM model, only the eddies larger than the size of the rotor are energetic enough to induce wake meandering. Thus the idea is to calibrate the amount of wake meandering only with eddies larger than this size as follows:
and similarly for σ_{fz}. In Eq. (40), ${k}_{y}^{\mathrm{D}}$ is the lateral turbulence with a size that is larger than the diameter of the turbine, i.e. the variance of the wind velocity averaged over a circle of 2 rotor diameters and centred at the hub. Note that the time variance is performed after the spatial averaging.
The issue is that ${k}_{y}^{\mathrm{D}}$ and ${k}_{z}^{\mathrm{D}}$ are not known a priori, and since the stability of the ABL modifies the lowfrequency range of the turbulence spectrum, it is expected that the share of the turbulence with a larger size than the rotor to the total turbulence is dependent on the atmospheric stability. This can be observed in Fig. 8 where the ratio between the turbulence larger than a disc of diameter d_{disc}, ${k}_{y,z}^{{\mathrm{d}}_{\mathrm{disc}}}$ to the total turbulence is computed for k_{y} and k_{z} for each case.
Figure 8 highlights two distinct behaviours, depending on the stability conditions: the unstable cases (orange and purple curves) decrease much slower than the nearneutral cases (red, grey, brown and blue), and this phenomenon is particularly marked for the lateral turbulence. It shows that the unstable cases have (in proportion) more lowfrequency or largesize eddies than the nearneutral cases.
Even though a fully physical approach would require a measure of the stability and an indepth study of the turbulence spectrum as a function of the ABL conditions, the objective here is to propose an analytical model that is easy to implement and use. It is thus proposed to model the ratios ${k}_{y}^{{\mathrm{d}}_{\mathrm{disc}}}/{k}_{y}$ and ${k}_{z}^{{d}_{\mathrm{disc}}}/{k}_{z}$ with an analytical function:
A leastsquare fit is used to determine the value of the parameter Γ. Two different fits are used in order to have one result for unstable cases and one for nearneutral cases. The results are given in Table 5, and Eq. (41) is plotted in Fig. 8 with dotted black and dashdotted black lines for the neutral and unstable values respectively.
Γ can be interpreted as a measure of the largescale eddies of the atmosphere, even though it is not defined as the integral length scale. The combination of Eqs. (40) and (41) with values from Table 5 is plotted with dotted lines in Fig. 7. Even though the model cannot predict the nonlinear behaviour in the far wake, the results remain quite good. Only the weakly unstable case gives poor results, possibly because it is at the edge between the nearneutral and unstable case, and necessitates a Γ value of its own.
In this section, we analyse the results of the new model described in the precedent sections. For the streamwise velocity, the model is described with Eq. (28) and for streamwise turbulence with the sum of Eqs. (30) and (33). This validation is done with the two validation cases (see Table 1), i.e. with the unstable and neutral SWiFT simulations. Three versions of the calibration of σ, σ_{fy}, σ_{fz} and ${l}_{\mathrm{m}}^{*}$ are shown:

The base calibration is defined with Eqs. (36), (38) and (39) as well as values for a, b, c, d and e from Tables 3 and 4. This calibration makes more sense physically but requires the time series of the inflow velocities to determine the autocorrelation necessary to compute σ_{f} from Eq. (39). It is plotted with dashed blue lines in Figs. 9–12.

The engineering calibration uses the same equations except for the wake meandering, where Eqs. (40) and (41) are used instead of Eq. (38) and parameter Γ is taken from Table 5. It is plotted with dotted red lines in Figs. 9–12.

Finally, we also propose the “best” version of the model. Knowing that the calibration produces errors, it was interesting to see what would be the results of the “best calibration possible”, i.e. with parameters σ, σ_{fy}, σ_{fz} and ${l}_{\mathrm{m}}^{*}$ directly taken from the LES of the SWiFT simulation (and not from the calibration deduced from Sect. 4). Obviously, this version of the model cannot be used, but it is helpful to determine if the discrepancies between our model and the LES come from the calibration or the construction of the model itself. It is plotted with dashdotted orange lines in Figs. 9–12.
Additionally, in Figs. 9–12 the reader will find the results directly from MesoNH (solid black line) and the result from a widely used model of the I&Q2018 model (dashdotteddotted purple line), one of the few in the literature that predicts the profiles of both mean streamwise velocity and streamwise turbulence.
5.1 Velocity field
The results for the streamwise velocity field in the FFOR can be found in Figs. 9 and 10 for the neutral and unstable cases respectively. The horizontal (top) and vertical (bottom) profiles of velocity are plotted for the reference LES, results from the literature and the three versions of the aforementioned model's calibration. The three columns are three different positions downstream of the wind turbine: $x/D=\mathrm{2}$, $x/D=\mathrm{5}$ and $x/D=\mathrm{8}$.
Our model (with any calibration) behaves very similarly to the I&Q2018 model in the neutral case (Fig. 9), and both are accurate compared to the LES data in black. The only discrepancy is in the near wake, where both models assume a Gaussian shape, whereas a superGaussian shape (Blondel and Cathelain, 2020) would be more appropriate. These overall good results confirm that the hypotheses made in Sect. 3 for the velocity in the MFOR and the wake centre distribution are good and that meandering has been correctly computed.
In the unstable case (Fig. 10), the literature model underestimates the wake dissipation, whereas the proposed model is more accurate. This is because the I&Q2018 model only uses C_{T} and I_{x} as parameters. As shown in Table 1, these values are very similar in the neutral and unstable cases of the validation case, and the I&Q2018 results are thus very similar between the neutral and unstable cases. It cannot predict the increase in meandering under unstable ABL conditions due to higher values of largescale turbulence in the lateral and vertical directions.
The proposed model is better in that regard, showing a larger wake expansion due to the higher predicted meandering compared to the neutral case. It shows that the determination of the velocity deficit in nonneutral cases necessitates more than only the total streamwise turbulence. In this case, one can note a discrepancy between the best version and the two calibrations of our model. It is due to an overestimation of the wake width in the MFOR for the unstable case (not shown here). Indeed, the neutral and unstable cases have similar wake widths in the MFOR, while they have different total turbulence intensities I (Table 1), and therefore Eq. (36) gives accurate results for the neutral case but overestimates the MFOR wake width in the unstable case. As a result, there is a compensation of error, where the calibration underestimates the velocity deficit, whereas the best version is supposed to slightly overestimate it, resulting in a very good match. Nevertheless, even without this error compensation, the best version still outperforms the literature model.
5.2 Turbulence field
With the same plotting convention as in Figs. 9 and 10, the profiles of turbulence in the horizontal and vertical directions are plotted in Figs. 11 and 12 for the neutral and unstable cases respectively.
In the neutral case (Fig. 11), the I&Q2018 model performs remarkably well. It correctly predicts the location of the double peak in the horizontal direction and the top tip peak in the vertical direction. The proposed model shows less good results: despite the order of magnitude being accurate, the shape of the function is not, and the top tip maximum is not correctly positioned. Since the calibrations do not significantly differ from the best version of the model, this is attributed to modelling errors and not to the calibration. The authors suggest that this deviation originates from the omission of shear in the modelling of rotoradded turbulence, as described in Eq. (33).
The unstable case shows the main shortcomings of the I&Q2018 model and the added value of our model. As shown previously, the I&Q2018 model gives similar results between the unstable and neutral SWiFT cases because they have similar inflow I_{x} and C_{T} values. However, the MesoNH simulations show significant differences, in particular the fact that around x=5 D, the turbulence profile is unimodal in the unstable case and bimodal in the neutral case. This difference cannot be predicted by the I&Q2018 model as it always assumes a bimodal shape with a maximum at the top tip. However, this change in shape can be predicted by our model since both Eqs. (30) and (33) are bimodal when σ≫σ_{f} and unimodal when σ≪σ_{f}.
Except for the upstream turbulence profiles, the inflow conditions used in the I&Q2018 are very similar between the neutral and unstable cases. Consequently, the purple profiles are alike in Figs. 11 and 12, whereas the stronger meandering in the unstable case leads to a Gaussianlike turbulence profile, even in the vertical direction. The maximum turbulence is thus no longer located at the top tip but rather at hub height. This property is well predicted by our model but not by the I&Q2018 model, which does not take meandering into account and which predicts quasiidentical behaviours between the neutral and unstable cases. As shown in Figs. 5 and 7, the amount of meandering starts lower but grows faster than the wake width in the MFOR, in particular in unstable conditions. Hence, one can expect a bimodal shape in the near wake and a unimodal shape in the far wake, as seen in Figs. 11 and 12.
However, the calibration of our model leads to an overestimation of the streamwise turbulence, in particular in the near wake. Since there are not many differences between the basic and engineering calibrations, it is not attributed to the meandering calibration (these two calibrations only differ by the meandering modelling) but rather to the overestimated σ in the MFOR, as well as an overestimated ${l}_{\mathrm{m}}^{*}$. When computed directly from the simulation, the values of ${l}_{\mathrm{m}}^{*}$ are very similar between the neutral and unstable cases, whereas the values of ${l}_{\mathrm{m},\mathrm{\infty}}^{*}$ are much greater in the unstable case (see Fig. 3). Therefore, the value of ${l}_{\mathrm{m}}^{*}$ is overestimated by the model, leading to an overestimation of the rotoradded turbulence and thus to the total turbulence.
The best version of the model gives interesting results, showing that if a better calibration was achieved, in particular for the modified mixing length, the results of the model would be better. This question is further detailed in the next section.
The previous section shows the results of the model developed in this paper. It is quite good for the streamwise velocity field but can be improved for the turbulence, where the fully empirical model of Ishihara and Qian (2018) shows overall better results in neutral cases but has shortcomings in unstable cases. However, since it is physically based, we know the assumptions of the present model and thus have clear opportunities for improvements. The main ones known by the authors are listed below. Moreover, this work shows that the modification of the velocity and turbulence fields when the ABL stability is modified (and not the I_{x} or C_{T}) can be predicted. This is a crucial point, as future applications of analytical models, such as digital twins, will require an estimation of the wake velocity and turbulence over small time lapses and not a yearly average like annual energy production (AEP) calculations.
The authors want to emphasize that the presented work is a first step toward a fully physically based model for turbulence profiles that depend on atmospheric stability. In the companion paper, it was shown that the turbulence in the wake of a wind turbine is the sum of several terms, and here we present a methodology to analytically model the most important of these terms. Even though a fully usable calibration is proposed for anyone who would like to test the model, the main purpose of this work is to demonstrate how the rotoradded turbulence and meandering turbulence can be modelled from simple functions.
6.1 Calibration improvement
In Figs. 9–12, there are discrepancies between the best version of the model and our proposed calibrations. This is particularly true for turbulence, and it is attributed to the calibration of ${l}_{\mathrm{m}}^{*}$. Contrarily to σ and σ_{f}, which can be computed on a wake no matter what, our computation of the modified mixing length ${l}_{\mathrm{m}}^{*}$ makes sense only if it is assumed that the rotoradded turbulence only comes from the wake shear. Additionally, the vertical velocity gradient of the ABL $\partial {U}_{\mathrm{\infty}}/\partial z$ is voluntarily omitted in Eq. (20).
On one hand, all of these assumptions make the measure of ${l}_{\mathrm{m}}^{*}$ a hardly reliable variable. On the other hand, our model is strongly dependent on this parameter. Indeed, the rotoradded turbulence is proportional to the square of ${l}_{\mathrm{m}}^{*}$. Therefore, a small over or underestimation of ${l}_{\mathrm{m}}^{*}$ is likely to happen, and it leads to large differences. Figure 13 shows the effect of multiplying parameter d of the mixing length (Table 4) by a factor of 0.8 (dotted red line), 1.2 (dashdotted orange line) and 1.5 (dashdotteddotted purple line) for the basic calibration in the neutral case. It results in large differences from one result to another, showing that even small differences in ${l}_{\mathrm{m}}^{*}$ can drastically change the conclusions.
6.2 Modellization improvements
Besides a better calibration, the model could benefit from conceptual improvements. Indeed, the best version of the model (orange curve in Figs. 11 and 12) does not match the LES results. In other words, even with a “perfect” calibration, the model still misses some features of the turbulence in the wake.
At several points of the reasoning, the atmospheric shear, i.e. the dependence of U_{∞} on z, is neglected (Eqs. 24 and 20). The first improvement that comes to mind is to model the interaction between atmospheric and wake shear. By doing so, it would be possible to have the reduction in shear near the ground and an increase in shear at the top tip, leading to a smaller value of turbulence at the bottom tip compared to the top tip, as observed in the LES datasets and modelled in the I&Q2018 model. In the current form of the model, the shear is only accounted for through ${U}_{\mathrm{\infty}}^{\mathrm{2}}$ as a factor in k_{m,am} and k_{a,am}. This small contribution is compensated for by the upstream turbulence k_{∞} that is larger at the bottom than at the top, leading to almost symmetric vertical profiles for the model, whereas LES profiles clearly show higher values at the top tip of the wake (at least for the neutral case and in the near wake of the unstable case).
A second improvement that could be done concerns the near wake. As mentioned in Sect. 5, instead of using a simple Gaussian function, a superGaussian function would be more accurate. This generic function takes a tophat form in the near wake and progressively transitions to a Gaussian function as it travels downstream. It was shown in Blondel and Cathelain (2020) that it gives more accurate results in the near wake. Such a function would improve not only the velocity model but also the meandering and rotoradded turbulence terms, which are built upon the velocity model. The latter in particular is a function of the spatial derivative of ΔU using the Gaussian function instead of the superGaussian function as done in this work, thus leading to an underestimation of the shear at the edge of the turbine.
For both of these improvements, some solutions were tried: not neglecting the $\partial {U}_{\mathrm{\infty}}/\partial z$ in the derivation of the rotoradded turbulence and using a superGaussian function instead of a Gaussian function for the velocity in the MFOR. In both cases, no analytical solution for the models was reached. If such a fully analytical resolution is indeed impossible, an approximated form (for instance, based on LES results) could be proposed in the future.
Finally, modelling the additional terms of Eq. (7), in particular the covariance term (V), could further improve the model. It was shown in the companion paper that this term can represent about 10 % of the total turbulence in the wake and redistributes the turbulence vertically. Given the order of magnitude, this is of lesser importance than the aforementioned points but would also improve the results, or at least the physical accuracy of the model.
This work is the second part of a twostep study that aims to model the turbulence in the wake of a wind turbine based on the meandering phenomenon. In the companion paper, the velocity and turbulence in the FFOR were broken down into different terms, some of which were shown to be negligible. In the present work, an analytical model is proposed for the dominating terms of the velocity and turbulence breakdowns, i.e. the meandering turbulence and the rotoradded turbulence. The originality of this work is that it allows for the independent modelling of the effects of meandering (and thus of the ABL stability) and the wake expansion and that it gives the whole turbulence profile rather than only the maximum value. For the velocity, it writes
and for the turbulence
where $C=\mathrm{1}\sqrt{\mathrm{1}{C}_{\mathrm{T}}/(\mathrm{8}{\mathit{\sigma}}_{y}{\mathit{\sigma}}_{z}/{D}^{\mathrm{2}})}$, C_{T} is the thrust coefficient, D is the turbine diameter, and k_{x,∞} and U_{∞} are the variance and mean values of the upstream axial velocity. The model's parameters are the wake widths σ_{y}, σ_{z}; the amount of meandering σ_{fy},σ_{fz}; and the modified mixing length ${l}_{\mathrm{m}}^{*}$. Two calibrations of these parameters are proposed in Table 6: the first one (base calibration) can be used if time series of the wind velocity are available and the second one (engineering calibration) if they are not. In this table, 𝒜_{ϕ} is the autocorrelation of ϕ, U_{c}=0.8U_{∞} and ${l}_{\mathrm{m},\mathrm{\infty}}^{*}$ is found by fitting the inflow velocity profile (Eq. 14). The expressions of velocity and added turbulence in the MFOR used to build Eqs. (42) and (43) can also be used as inputs to the DWM: combined with a synthetic turbulence generation, the unsteady effects of meandering can be modelled.
The model has been tested on two LESs of a single wind turbine wake in a neutral and unstable atmosphere. For the velocity, the results are satisfactory, either in the vertical or in the lateral direction. The model performs better than the model from Ishihara and Qian (2018) in the unstable case as it correctly predicts the increased dissipation due to the increase in meandering. For the turbulence profiles, however, the results are not as good. Since the atmospheric shear was neglected in several steps of the model, the maximum turbulence at the top tip in the neutral case could not be predicted. In the unstable case, the modified mixing length ${l}_{\mathrm{m}}^{*}$ was overestimated, and since the model is very sensitive to this parameter, it resulted in too large values of added turbulence. However, the model of Ishihara and Qian (2018) does not correctly predict the turbulence in the unstable case either. In particular, it still predicts a bimodal shape with a maximum at the top tip in the whole wake, whereas the proposed model successfully transitions from a bimodal to a unimodal shape, according to the LES results.
This is the first step toward a fully analytical, physically based model for turbulence and velocity profiles in the wake of a wind turbine that takes into account atmospheric stability. For future work, the treatment of shear must be improved to model vertical turbulence profiles more realistically. The MFOR velocity deficit function could be replaced by a more accurate function in the near wake to improve the model's results in this region. It would also be interesting to derive an analytical model for the other terms of the turbulence breakdown.
Finally, this model can currently only be used for one turbine, as it predicts only the streamwise velocity and turbulence but necessitates the upstream lateral and vertical turbulence. For the model to be usable for multiturbines, an expression for every term of the Reynolds stress tensor (or at least the diagonal terms to get the total TKE) would be needed, which implies a model for the lateral and vertical velocities U_{y} and U_{z}. This also implies more advanced studies of wake meandering from a turbine working in waked conditions, as most of the wake meandering studies are performed in freestream conditions.
The MesoNH code is open source and can be downloaded on the dedicated website. The authors can provide the source code of the modified version 543 that was used in this work. The data used for the plot presented here and in Part 1 (Jézéquel et al., 2024) are available through the following online deposit: https://doi.org/10.5281/zenodo.6562720 (Jézéquel, 2022). The data for the calibration can be found under the following DOI: https://doi.org/10.25326/568 (Jézéquel, 2023). The model equations were written in Python and can be found through the following online deposit: https://doi.org/10.5281/zenodo.10245174 (Jézéquel and Blondel, 2023).
EJ wrote the analytical model with FB. All the authors worked on the interpretation of the results. The paper was written by EJ with feedback from FB and VM.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This paper was edited by Raúl Bayoán Cal and reviewed by three anonymous referees.
Abkar, M. and PortéAgel, F.: Influence of atmospheric stability on windturbine wakes: A largeeddy simulation study, Phys. Fluids, 27, 035104, https://doi.org/10.1063/1.4913695, 2015. a
Abkar, M., Sørensen, J., and PortéAgel, F.: An Analytical Model for the Effect of Vertical Wind Veer on Wind Turbine Wakes, Energies, 11, 1838, https://doi.org/10.3390/en11071838, 2018. a
Ainslie, J. F.: Calculating the flowfield in the wake of wind turbines, J. Wind Eng. Indust. Aerodynam., 27, 213–224, https://doi.org/10.1016/01676105(88)900372, 1988. a
Bastankhah, M. and PortéAgel, F.: A new analytical model for windturbine wakes, Renew. Energy, 70, 116–123, https://doi.org/10.1016/j.renene.2014.01.002, 2014. a, b, c, d, e
Blondel, F. and Cathelain, M.: An alternative form of the superGaussian wind turbine wake model, Wind Energ. Sci., 5, 1225–1236, https://doi.org/10.5194/wes512252020, 2020. a, b, c
Braunbehrens, R. and Segalini, A.: A statistical model for wake meandering behind wind turbines, J. Wind Eng. Indust. Aerodynam., 193, 103954, https://doi.org/10.1016/j.jweia.2019.103954, 2019. a, b, c, d, e, f, g
Brugger, P., Markfort, C., and PortéAgel, F.: Field measurements of wake meandering at a utilityscale wind turbine with nacellemounted Doppler lidars, Wind Energ. Sci., 7, 185–199, https://doi.org/10.5194/wes71852022, 2022. a
Cheng, Y., Zhang, M., Zhang, Z., and Xu, J.: A new analytical model for wind turbine wakes based on MoninObukhov similarity theory, Appl. Energy, 239, 96–106, https://doi.org/10.1016/j.apenergy.2019.01.225, 2019. a
Conti, D., Dimitrov, N., Peña, A., and Herges, T.: Probabilistic estimation of the Dynamic Wake Meandering model parameters using SpinnerLidarderived wake characteristics, Wind Energ. Sci., 6, 1117–1142, https://doi.org/10.5194/wes611172021, 2021. a, b
Crespo, A. and Hernandez, J.: Turbulence characteristics in windturbine wakes, J. Wind Eng. Indust. Aerodynam., 61, 71–85, https://doi.org/10.1016/01676105(95)00033x, 1996. a
Doubrawa, P., MartínezTossas, L. A., Quon, E., Moriarty, P., and Churchfield, M. J.: Comparison of Mean and Dynamic Wake Characteristics between ResearchScale and FullScale Wind Turbines, J. Phys.: Conf. Ser., 1037, 072053, https://doi.org/10.1088/17426596/1037/7/072053, 2018. a
Du, B., Ge, M., Zeng, C., Cui, G., and Liu, Y.: Influence of atmospheric stability on wind turbine wakes with a certain hubheight turbulence intensity, Phys. Fluids, 33, 055111, https://doi.org/10.1063/5.0050861, 2021. a, b
Frandsen, S.: Turbulence and turbulencegenerated structural loading in wind turbine clusters, PhD thesis, risøR1188(EN), DTU, ISBN 8755034586, 2007. a
Fuertes, F. C., Markfort, C., and PortéAgel, F.: Wind Turbine Wake Characterization with NacelleMounted Wind Lidars for Analytical Wake Model Validation, Remote Sens., 10, 668, https://doi.org/10.3390/rs10050668, 2018. a, b
Ishihara, T. and Qian, G.W.: A new Gaussianbased analytical wake model for wind turbines considering ambient turbulence intensities and thrust coefficient effects, J. Wind Eng. Indust. Aerodynam., 177, 275–292, https://doi.org/10.1016/j.jweia.2018.04.010, 2018. a, b, c, d, e
Iungo, G. V., Santhanagopalan, V., Ciri, U., Viola, F., Zhan, L., Rotea, M. A., and Leonardi, S.: Parabolic RANS solver for lowcomputationalcost simulations of wind turbine wakes, Wind Energy, 21, 184–197, https://doi.org/10.1002/we.2154, 2017. a
Jensen, N.: A note on wind turbine interaction, techreport,Tech. Rep. RisøM2411, Risø National Laboratory, Denmark, 16 pp., https://orbit.dtu.dk/files/55857682/ris_m_2411.pdf (last access: 5 December 2022), 1983. a
Jézéquel, E.: Figures data from papers “Breakdown of the velocity and turbulence in the wake of a wind turbine”, parts 1 and 2 [Data set], Zenodo [data set], https://doi.org/10.5281/zenodo.6562720, 2022. a
Jézéquel, E.: SIMU ATMO, Aeris [data set], https://doi.org/10.25326/568, 2023. a, b, c
Jézéquel, E. and Blondel, F.: Implementation of the analytical model deduced from the velocity and turbulence breakdown, version of the reviewed paper, Zenodo [code], https://doi.org/10.5281/zenodo.10245174, 2023. a
Jézéquel, E., Blondel, F., and Masson, V.: Breakdown of the velocity and turbulence in the wake of a wind turbine – Part 1: Largeeddysimulation study, Wind Energ. Sci., 9, 97–117, https://doi.org/10.5194/wes9972024, 2024. a
Keane, A., Aguirre, P. E. O., Ferchland, H., Clive, P., and Gallacher, D.: An analytical model for a full wind turbine wake, J. Phys.: Conf. Ser., 753, 032039, https://doi.org/10.1088/17426596/753/3/032039, 2016. a
Keck, R.E., Maré, M. D., Churchfield, M. J., Lee, S., Larsen, G., and Madsen, H. A.: Two improvements to the dynamic wake meandering model: including the effects of atmospheric shear on wake turbulence and incorporating turbulence buildup in a row of wind turbines, Wind Energy, 18, 111–132, https://doi.org/10.1002/we.1686, 2013. a
Lac, C., Chaboureau, J.P., Masson, V., Pinty, J.P., Tulet, P., Escobar, J., Leriche, M., Barthe, C., Aouizerats, B., Augros, C., Aumond, P., Auguste, F., Bechtold, P., Berthet, S., Bielli, S., Bosseur, F., Caumont, O., Cohard, J.M., Colin, J., Couvreux, F., Cuxart, J., Delautier, G., Dauhut, T., Ducrocq, V., Filippi, J.B., Gazen, D., Geoffroy, O., Gheusi, F., Honnert, R., Lafore, J.P., Brossier, C. L., Libois, Q., Lunet, T., Mari, C., Maric, T., Mascart, P., Mogé, M., Molinié, G., Nuissier, O., Pantillon, F., Peyrillé, P., Pergaud, J., Perraud, E., Pianezze, J., Redelsperger, J.L., Ricard, D., Richard, E., Riette, S., Rodier, Q., Schoetter, R., Seyfried, L., Stein, J., Suhre, K., Taufour, M., Thouron, O., Turner, S., Verrelle, A., Vié, B., Visentin, F., Vionnet, V., and Wautelet, P.: Overview of the MesoNH model version 5.4 and its applications, Geosci. Model Dev., 11, 1929–1969, https://doi.org/10.5194/gmd1119292018, 2018. a
Larsen, G. C., Madsen Aagaard, H., Bingöl, F., Mann, J., Ott, S., Sørensen, J. N., Okulov, V., Troldborg, N., Nielsen, N. M., Thomsen, K., Larsen, T. J., and Mikkelsen, R.: Dynamic wake meandering modeling, Risø National Laboratory, ISBN 9788755036024, 2007. a
Larsen, G. C., Madsen, H. A., Thomsen, K., and Larsen, T. J.: Wake meandering: a pragmatic approach, Wind Energy, 11, 377–395, https://doi.org/10.1002/we.267, 2008. a, b
Madsen, H. A., Larsen, G. C., Larsen, T. J., Troldborg, N., and Mikkelsen, R.: Calibration and Validation of the Dynamic Wake Meandering Model for Implementation in an Aeroelastic Code, J. Sol. Energ. Eng., 132, 4, https://doi.org/10.1115/1.4002555, 2010. a, b
Niayifar, A. and PortéAgel, F.: Analytical Modeling of Wind Farms: A New Approach for Power Prediction, Energies, 625, 012039, https://doi.org/10.1088/17426596/625/1/012039, 2016. a
Pope, S. B.: Turbulent Flows, Cambridge University Press, https://doi.org/10.1017/cbo9780511840531, 2000. a
Scherfgen, D.: Integral Calculator, https://www.integralcalculator.com/ (last access: 29 April 2022), 2022. a
Stein, V. P. and Kaltenbach, H.J.: NonEquilibrium Scaling Applied to the Wake Evolution of a Model Scale Wind Turbine, Energies, 12, 2763, https://doi.org/10.3390/en12142763, 2019. a
Teitelbaum, J.: Convolution of Gaussians is Gaussian, https://jeremy9959.net/Math5800Spring2020/notebooks/convolution_of_gaussians.html (last access: 29 April 2022), 2022. a, b
Xie, S. and Archer, C.: Selfsimilarity and turbulence characteristics of wind turbine wakes via largeeddy simulation, Wind Energy, 18, 1815–1838, https://doi.org/10.1002/we.1792, 2014. a, b