Wind farm models¶

Introduction¶

ERF supports models for wind farm parametrization in which the effects of wind turbines are represented by imposing a momentum sink on the mean flow and/or turbulent kinetic energy (TKE). Currently the Fitch model (Fitch et al. 2012), Explicit Wake Parametrization (EWP) model (Volker et al. 2015), Simplified actuator disk model (See Section 3.2 in Wind Energy Handbook 2nd edition), and Generalized actuator disk model (Mirocha et. al. 2014, see Chapter 3 of Small Wind Turbines) are supported.

Instantaneous horizontal velocity magnitude at the hub-height of 89 [m] and t= 1 [hr] for: (row-wise) Fitch model, EWP model, Simple actuator disk model, and General actuator disk model.¶

Fitch model¶

The Fitch model for wind farms introduced in Fitch et al. 2012 models the effect of wind farms (See Fig. 1) as source terms in the governing equations for the horizontal components of momentum (i.e., \(x\) and \(y\) momentum) and the turbulent kinetic energy (TKE). The wind turbine is discretized only in the vertical (ie. z direction). At a given cell \((i,j,k)\), the source terms in the governing equations are

\[\begin{split}\frac{\partial u_{ijk}}{\partial t} &= \frac{u_{ijk}}{|V|_{ijk}}\frac{\partial |V|_{ijk}}{\partial t} \\ \frac{\partial v_{ijk}}{\partial t} &= \frac{v_{ijk}}{|V|_{ijk}}\frac{\partial |V|_{ijk}}{\partial t} \\ \frac{\partial \text{TKE}_{ijk}}{\partial t} &= \frac{0.5N_{ij}C_{TKE}(|V|_{ijk})|V|_{ijk}^3A_{ijk}}{\Delta x \Delta y (z_{k+1}-z_k)}\end{split}\]

where

\[\frac{\partial |V|_{ijk}}{\partial t} = \frac{0.5N_{ij}C_T(|V|_{ijk})|V|_{ijk}^2A_{ijk}}{\Delta x \Delta y (z_{k+1}-z_k)}\]

where u and v are horizontal components of velocity, |V| is the velocity magnitude, \(N_{ij}\) is the number of turbines in cell \((i,j)\) (see Fig. fig:WindFarm), \(C_T\) is the coefficient of thrust of the turbines, \(C_{TKE}\) is the fraction of energy converted to turbulent kinetic energy – both of which are functions of the velocity magnitude, and \(A_{ijk}\) is the area intersected by the swept area of the turbine between levels \(z=z_k\) and \(z= z_{k+1}\), and is explained in the next section.

Intersected area \(A_{ijk}\)¶

Consider \(A_k^{k+1}\) – the area intersected by the swept area of the wind turbine between \(z=z_k\) and \(z = z_{k+1}\). We have (see Figs. fig:WindTurbine_Fitch and fig:Fitch_Aijk below)

\[A_k = \frac{\pi R^2}{2} - A_{ks}\]

where \(A_{ks}\) is the area of the segment of the circle as shown in Fig. fig:Fitch_Aijk. We have from geometry, \(d_k = \min(|z_k - h|,R)\) is the perpendicular distance of the center of the turbine to \(z = z_k\), where \(h\) is the height of the center of the turbine from the ground. The area of the segment is

\[A_{ks} = R^2\theta_k - d_k\sqrt{R^2 - d_k^2}\]

where \(\theta_k = \cos^{-1}\left(\frac{d_k}{R}\right)\).

Hence, we have the intersected area \(A_{ijk}\equiv A_k^{k+1}\) as

\[\begin{split}A_k^{k+1} = \begin{cases} |A_k - A_{k+1}| & \text{if } (z_k - h)(z_{k+1}-h) > 0 \\ |A_k + A_{k+1}| & \text{if } (z_k - h)(z_{k+1}-h) \le 0 \\ \end{cases}\end{split}\]

An example of the Fitch model is in .Exec_dev/WindFarmTests/SingleTurbine/Fitch

Horizontal view of the wind farm in the Fitch model – shows a wind farm in cell (i,j) with 5 wind turbines. The turbines are discretized only in the vertical direction.¶

The vertical discretization of the wind turbine in the Fitch model.¶

The area terminology in the Fitch model. The circle represents the area swept by the turbine blades.¶

Explicit Wake Parametrization (EWP) model¶

The Explicit Wake Parametrization (EWP) model [Volker et al. 2015] is very similar to the Fitch model, and models the effect of wind farms as source terms in the governing equations for the horizontal components of momentum (i.e., \(x\) and \(y\) momentum) and the turbulent kinetic energy (TKE). At a given cell \((i,j,k)\), the source terms in the governing equations are:

\[\frac{\partial u_{ijk}}{\partial t} = -\sqrt{\frac{\pi}{8}}\frac{N_{ij}c_tr_0^2\overline{u}_0^2}{\Delta x \Delta y \sigma_e} \exp\left\{-\frac{1}{2}\left(\frac{z-h}{\sigma_e}\right)^2\right\}\cos(\phi(k))\]

\[\frac{\partial v_{ijk}}{\partial t} = -\sqrt{\frac{\pi}{8}}\frac{N_{ij}c_tr_0^2\overline{u}_0^2}{\Delta x \Delta y \sigma_e} \exp\left\{-\frac{1}{2}\left(\frac{z-h}{\sigma_e}\right)^2\right\}\sin(\phi(k))\]

\[\frac{\partial \text{TKE}_{ijk}}{\partial t} = -N_{ij}\rho A_rc_t\langle \overline{u}_{i,h}\overline{u'^2}_{i,h}\rangle\]

with

\[\sigma_e = \frac{\overline{u}_0}{3KL}\left[\left(\frac{2KL}{\overline{u}_0} + \sigma_0^2\right)^{\frac{3}{2}} - \sigma_0^3\right]\]

where \(N_{ij}\) is the number of turbines in cell \((i,j)\), \(c_t\) is the thrust coefficient, \(r_0\) is the rotor radius, \(\overline{u}_0\) is the mean advection velocity at hub height, \(h\) is the hub height, \(\sigma_0 \approx 1.7 r_0\) [Volker et al. 2017] is a length scale that accounts for near-wake expansion, \(L\) is the downstream distance that the wake travels within the cell approximated as a fraction of the cell size, \(K\) is the turbulence eddy diffusivity (\(m^2/s\)), \(\Delta x\) and \(\Delta y\) are the mesh sizes in the horizontal directions, and \(\phi(k)\) is the wind direction with respect to the x-axis. \(\overline{u}_{i,h}\) and \(\overline{u'}_{i,h}\) are the mean and the fluctuating values of the velocity components (subscript \(i\) is the direction index) at the hub height \(h\), \(A_r = \pi r^2\) is the swept area of the rotor and \(\rho\) is the density of air.

The EWP model does not have a concept of intersected area by the turbine rotor like the Fitch model (see Fitch model). The exponential factor in the source terms for the velocities models the effect of the rotor for the momentum sinks (see Fig. fig:WindTurbine_EWP), and the turbulent kinetic energy source term only depends on the density, hub-height mean velocities and fluctuations, and the total swept area of the rotor \(A_r\).

In the EWP model, the exponential factor in the source terms for the velocities models the effect of the rotor for the momentum sinks unlike the Fitch model which computes the intersected area (see Fig. fig:WindTurbine_Fitch).¶

Simplified actuator disk model¶

A simplified actuator disk model based on one-dimensional momentum theory is implemented (See Section 3.2 in Wind Energy Handbook 2nd edition). A schematic of the actuator disk is shown in Fig. fig:ActuatorDisk_Schematic. The model is implemented as source terms in the equations for the horizontal velocity components (ie. x and y directions). The thrust force from the one-dimensional momentum theory is given by

\[\begin{split}F = 2 \rho \pi R^2 (\mathbf{U}_\infty \cdot \mathbf{n})^2 a (1-a) \\ = \int_0^{2\pi}\int_0^R 2 \rho (\mathbf{U}_\infty \cdot \mathbf{n})^2 a (1 - a) r\,dr\,d\theta,\end{split}\]

where \(\rho\) is the density of incoming air, \(\mathbf{U}_\infty\) is the velocity vector of incoming air at some distance (say \(d=2.5\) times the turbine diameter) upstream of the turbine (see Fig. fig:ActuatorDisk_Sampling), \(\mathbf{n}\) is the surface normal vector of the actuator disk, and \(a = 1 - \cfrac{C_P}{C_T}\), is the axial induction factor for the turbine, and \(R\) is the radius of the wind turbine swept area. The integration is performed over the swept area of the disk. Hence, the force on an elemental annular disk of thickness \(dr\) is

\[dF = 2 \rho (\mathbf{U}_\infty \cdot \mathbf{n})^2 a (1 - a) r\,dr\,d\theta,\]

where \(r~dr~d\theta\) is the elemental area of the actuator disk. In general, this can be written as

\[\begin{split}dF = 2 \rho (\mathbf{U}_\infty \cdot \mathbf{n})^2 a (1 - a) dA \\\end{split}\]

where \(dA\) is the area of the actuator disk in the mesh cell (see Fig. fig:ActuatorDisk_Schematic). In the context of the simplified actuator disk model, the source term is imposed only on a single cell, and hence the volume over which the force \(dF\) is acting is the volume of the cell — \(dV \equiv \Delta x \Delta y \Delta z\). Hence, the source terms in the horizontal velocity equations are the acceleration (or deceleration) due to the thrust force \(dF\) and is given by

\[\begin{split}\frac{\partial u}{\partial t} &= -\frac{dF \cos{\phi}}{\rho dV} \\ &= -\frac{dF \cos{\phi}}{\rho \times \Delta x \Delta y \Delta z} \\ &= -2(\mathbf{U}_\infty \cdot \mathbf{n})^2 a (1 - a)\frac{\Delta A}{\Delta x\Delta y\Delta z} \cos{\phi},\end{split}\]

\[\frac{\partial v}{\partial t} = -2(\mathbf{U}_\infty \cdot \mathbf{n})^2 a (1 - a)\frac{\Delta A}{\Delta x\Delta y\Delta z} \sin{\phi}.\]

Schematic of the simplified actuator disk model.¶

Top view showing the freestream velocity sampling disk at a distance of \(d\) from the turbine actuator disk.¶

Generalized actuator disk model¶

The generalized actuator model (GAD) based on blade element theory (Mirocha et. al. 2014, see Chapter 3 of Small Wind Turbines) is implemented. Similar to the simplified actuator disk model, GAD also models the wind turbine as a disk, but takes into account the details of the blade geometry (see Fig. fig:GAD_Schematic). The forces on the blades in the x, y, z directions are computed, and that contributes to the source terms for the fluid momentum equations. The source terms in a mesh cell inside the actuator disk are given as:

()¶\[\begin{split}\frac{\partial u}{\partial t} &= -\frac{F_x}{\rho \Delta x\Delta y\Delta z} \\ \frac{\partial v}{\partial t} &= -\frac{F_y}{\rho \Delta x\Delta y\Delta z} \\ \frac{\partial w}{\partial t} &= -\frac{F_z}{\rho \Delta x\Delta y\Delta z},\end{split}\]

where \(\rho\) is the density of air in the cell, and \(\Delta x, \Delta y, \Delta z\) are the mesh spacing in the x, y, and z directions. The forces on the GAD are given by:

()¶\[\begin{split}F_x &= F_n \cos{\Phi} + F_t \sin\zeta \sin\Phi \\ F_y &= F_n \sin{\Phi} - F_t \sin\zeta \cos\Phi \\ F_z &= -F_t \cos\zeta,\end{split}\]

where \(F_n\) and \(F_t\) are the normal and tangential forces, and the angles are as shown in Figure fig:GAD_Schematic. The normal and tangential forces are:

()¶\[\begin{split}\begin{bmatrix} F_n \\ F_t \end{bmatrix} = \begin{bmatrix} \cos\Psi & \sin\Psi \\ \sin\Psi & -\cos\Psi \end{bmatrix} \begin{bmatrix} L \\ D \end{bmatrix},\end{split}\]

where

\[\Psi = \tan^{-1}\left(\frac{V_n}{V_t}\right),\]

and

()¶\[\begin{split}V_n &= V_0(1-a_n) \\ V_t &= \Omega(1+a_t)r,\end{split}\]

where \(V_0\) is the freestream velocity at a user specified distance upstream from the disk plane as described in Section Simplified actuator disk model (also see Fig. fig:ActuatorDisk_Sampling), \(\Omega\) is the rotational speed of the turbine, \(r\) is the radial location along the blade span, and \(a_n\) and \(a_t\) are the normal and tangential induction factors. The lift and drag forces are given by:

()¶\[\begin{split}L &= \frac{1}{2} \rho V_r^2 c C_l \\ D &= \frac{1}{2} \rho V_r^2 c C_d,\end{split}\]

where \(\rho\) is the density of air, \(c\) is the chord length of the airfoil cross-section, \(C_l\) and \(C_d\) are the sectional lift and drag coefficients on the airfoil cross-section (which is a function of the incoming angle \(\psi\), blade twist \(\xi\), and blade pitch \(\phi\). See Fig. fig:GAD_Schematic), and the relative wind velocity is \(V_r = \sqrt{V_n^2 + V_t^2}\). The normal and tangential sectional coefficients are computed as:

()¶\[\begin{split}\begin{bmatrix} C_n \\ C_t \end{bmatrix} = \begin{bmatrix} \cos\Psi & \sin\Psi \\ \sin\Psi & -\cos\Psi \end{bmatrix} \begin{bmatrix} C_l \\ C_d \end{bmatrix},\end{split}\]

and the normal and tangential induction factors are given by:

()¶\[\begin{split}a_n &= \left[1 + \frac{4F \sin^2\psi}{s C_n}\right]^{-1} \\ a_t &= \left[\frac{4F \sin\psi \cos\psi}{s C_t} - 1\right]^{-1},\end{split}\]

where

\[F = F_\text{tip} + F_\text{hub} = \frac{2}{\pi} \left[\cos^{-1}\left(\exp(-f_\text{tip})\right) + \cos^{-1}\left(\exp(-f_\text{hub})\right)\right],\]

and

\[\begin{split}f_\text{tip} &= B \frac{(r_\text{tip}-r)}{2r \sin\psi} \\ f_\text{hub} &= B \frac{(r-r_\text{hub})}{2r \sin\psi},\end{split}\]

where \(r_\text{hub}\) and \(r_\text{tip}\) are the radius of the hub and the blade tip from the center of rotation of the disk, \(r\) is the radial location along the blade span, and the solidity factor is \(s=\frac{cB}{2\pi r}\), where \(B\) is the number of blades.

An iterative procedure is needed to compute the source terms, and is as follows:

An initial value is assumed for the normal and tangential induction factors \(a_n\) and \(a_t\).
Compute the normal and tangential velocities from Eqn. (). .
From the tables for the turbine specifications, details of the blade geometry and the sectional coefficients of the airfoil cross sections, compute the values of \(C_l\) and \(C_d\) corresponding to the radial location \(r\) along the blade span and the angle of attack \(\alpha = \psi - \xi + \phi\).
Compute the normal and tangential sectional coefficients \(C_n\) and \(C_t\) from Eqn. ().
Compute the normal and tangential induction factors \(a_n\) and \(a_t\) using Eqn. ().
Repeat steps 2 to 5 until the error in the normal and tangential induction factors, \(a_n\) and \(a_t\), are less than \(1 \times 10^{-5}\).
Once the iterations converge, compute the sectional lift and drag forces, \(L\) and \(D\), using Eqn. ().
Compute the normal and tangential forces, \(F_n\) and \(F_t\), using Eqn. ().
Compute the forces on the disk using Eqn. ().
Compute the source terms in the momentum equation using Eqn. ().

Different views of the GAD showing the forces and angles involved: Blade cross section showing the normal (\(V_n\)) and tangential (\(V_t\)) components of velocity with the normal (\(a_n\)) and tangential (\(a_t\)) induction factors, relative wind velocity \(V_r\), blade twist angle \(\xi\), angle of relative wind \(\psi\), blade pitch angle \(\phi\), lift (\(L\)) and drag (\(D\)) forces, and normal (\(F_n\)) and tangential (\(F_t\)) forces; top view showing the flow direction and inclination angle \(\Phi\); and front view showing the actuator disk rotating clockwise.¶

Inputs for wind farm parametrization models¶

The following are the inputs required for wind farm simulations.

// The parametrization model - Fitch, EWP, SimpleActuatorDisk
erf.windfarm_type = "Fitch"

// How are the turbine locations specified? - using latitude-longitude
// format or x-y format? lat_lon or x_y
erf.windfarm_loc_type = "lat_lon"

// If using lat_lon, then the shift of the bounding box of the wind farm
// from the x and y axes should be given. This is to avoid boundary
// effects from the inflow boundaries. For example for 2 km shift from the
// x and y axes, it should be
erf.windfarm_x_shift     = 2000.0
erf.windfarm_y_shift     = 2000.0

// Table having the wind turbine locations
erf.windfarm_loc_table = "windturbines_1WT.txt"

// Table having the specifications of the wind turbines. All turbines are assumed to
// have the same specifications
erf.windfarm_spec_table = "wind-turbine_1WT.tbl"

// In addition to the above, for simplified actuator disk model the following parameters are needed

// The distance of the freestream velocity sampling disk from the turbine actuator
// disk
erf.sampling_distance_by_D = 2.5

// The angle of the turbine actuator disk from the x axis
erf.turb_disk_angle_from_x = 135.0

// In addition to the above, for generalized actuator disk model the following parameters are needed

// Table containing additional specification information of the wind turbine.
// See Note below
erf.windfarm_spec_table_extra = "NREL-2.82-127_performance.csv"

// Table containing details of blade geometry
// See Note below
erf.windfarm_blade_table = "NREL-2p8-127_AeroDyn15_blade.dat"

// Tables containing the sectional lift and drag coefficients for the
// blade airfoil cross-sections.
// See Note below.
erf.windfarm_airfoil_tables = "Airfoils"

Note

The format for the files required for the generalized actuator disk model are

erf.windfarm_spec_table_extra = NREL-2.82-127_performance.csv
erf.windfarm_blade_table = NREL-2p8-127_AeroDyn15_blade.dat
erf.windfarm_airfoil_tables = Airfoils.

erf.windfarm_type has to be one of the supported models - Fitch, EWP, SimpleActuatorDisk.
erf.windfarm_loc_type is a variable to specify how the wind turbine locations in the wind farm is specified. If using the latitude and longitude of the turbine location, this has to be lat_lon or if using x and y coordinates to specify the turbine locations, this input is x_y.
- If using lat_lon format, erf.windfarm_x_shift and erf.windfarm_y_shfit specifies the shift of the bounding box of the wind farm from the x and y axes, so as to place the wind farm sufficiently inside the domain to avoid inflow boundary effects.
- If using x_y format, there is no need to specify the erf.windfarm_x_shift and erf.windfarm_y_shift.
The erf.windfarm_loc_table specifies the locations of the wind turbines in the wind farm.
- For the latitude-longitude format, an example is as below. Each line specifies the latitude and longitude of the wind turbine location. The third entry simply has to be always 1 (WRF requires the third entry to be always 1, so maintaining same format here). The first entry means that the turbine is located at 35.7828 deg N, 99.0168 deg W (note the negative sign in the entry corresponding to West).
```
35.7828828829 -99.0168 1
35.8219219219 -99.0168 1
35.860960961 -99.0168 1
35.9 -99.0168 1
35.939039039 -99.0168 1
35.9780780781 -99.0168 1
36.0171171171 -99.0168 1
35.7828828829 -98.9705333333 1
```
- For the x-y format, an example is as below. Each line specifies the x and y coordinates of the wind turbine location in metres
```
89264.99080053 91233.3333309002
89259.1966417755 95566.6666710007
89254.1277665419 99900.0000000001
89249.7842982733 104233.333329
89246.1663427532 108566.6666691
89243.2739881117 112899.9999981
93458.6633652711 86900.0000019001
93450.4377452458 91233.3333309002
93442.9032518779 95566.6666710007
```
The erf.windfarm_spec_table gives the specifications of the wind turbines in the wind farm. All wind turbines are assumed to have the same specifications. Here is a sample specifications table.

4
0 178.0 0.130 2.0
 0.805    50.0
 0.805    50.0
 0.805    50.0
 0.805    50.0

The first line is the number of pairs of entries for the power curve and thrust coefficient (there are 4 entries in this table which are in the last four lines of the table). The second line gives the height in meters of the turbine hub, the diameter in meters of the rotor, the standing thrust coefficient, and the nominal power of the turbine in MW. The remaining lines (four in this case) contain the three values of: wind speed (m/s), thrust coefficient, and power production in kW.

Outputs¶

Turbine locations are written into turbine_locations.vtk.
If using an actuator disk model, all the actuator disks are written out to actuator_disks_all.vtk. The actuator disks which are enclosed by the computational domain are written out to actuator_disks_in_dom.vtk.

These vtk files can be visualized in both VisIt and ParaView. The turbine_locations.vtk can be visualized using the Points Gaussian feature in ParaView or the Mesh feature in VisIt. The actuator_disks_in_dom.vtk and actuator_disks_all.vtk files can be visualized using the Wireframe feature in ParaView or Mesh feature in VisIt.