Filtered actuator disks: Theory and application to wind turbine models in large eddy simulation

TL;DR

This paper introduces a vortex cylinder correction model for filtered actuator disks in large eddy simulations, improving wind turbine power predictions without requiring high-resolution grids.

Contribution

It develops a vortex cylinder model to quantify shed vorticity and derives a correction factor that enhances power prediction accuracy in LES of wind turbines.

Findings

Correction factor collapses power curves across filter widths

Enables accurate power measurement without high grid resolution

Reduces over-prediction of power in LES simulations

Abstract

The actuator disk model (ADM) continues to be a popular wind turbine representation in large eddy simulations (LES) of large wind farms. Computational restrictions typically limit the number of grid points across the rotor of each actuator disk and require spatial filtering to smoothly distribute the applied force distribution on discrete grid points. At typical grid resolutions, simulations cannot capture all of the vorticity shed behind the disk and subsequently over-predict power by upwards of 10%. To correct these modeling errors, we propose a vortex cylinder model to quantify the shed vorticity when a filtered force distribution is applied at the actuator disk. This model is then used to derive a correction factor for numerical simulations that collapses the power curve for simulations at various filter widths and grid resolutions onto the curve obtained using axial momentum theory. The proposed correction therefore facilitates accurate power measurements in LES without resorting to highly refined numerical grids.

Figures

Equations

$$F=-\frac{1}{2}\rho \pi R^2{C}_T^{\prime }{u}_d^2.$$

$$\mathbf{f}(\mathbf{x})=F\mathcal{R}(\mathbf{x})\hat{\mathbf{x}}.$$

$$u_d=\int \mathcal{R}(\mathbf{x})\mathbf{u}(\mathbf{x})\cdot \hat{\mathbf{x}}{d}^3\mathbf{x},$$

$$\mathcal{I}(\mathbf{x})=\frac{1}{s\pi R^2}\left[H(x+s/2)-H(x-s/2)\right]H(R-r),$$

$$\mathcal{R}(\mathbf{x})=\int G(\mathbf{x}-{\mathbf{x}}^{\prime })\mathcal{I}({\mathbf{x}}^{\prime })d^3{\mathbf{x}}^{\prime }$$

$$G(\mathbf{x})={\left(\frac{6}{\pi {\Delta }^2}\right)}^{3/2}\exp \left(-\frac{6\Vert \mathbf{x}{\Vert }^2}{{\Delta }^2}\right)$$

$$\mathcal{R}(\mathbf{x})={\mathcal{R}}_1(x){\mathcal{R}}_2(r).$$

$${\mathcal{R}}_1(x)=\frac{1}{s}{\left(\frac{6}{\pi {\Delta }^2}\right)}^{1/2}\int \left[H\left(x^{\prime }+\frac{s}{2}\right)-H\left(x^{\prime }-\frac{s}{2}\right)\right]\exp \left(-6\frac{(x-x^{\prime }{)}^2}{{\Delta }^2}\right)d{x}^{\prime }$$

$${\mathcal{R}}_1(x)=\frac{1}{2s}\left[\mathrm{erf}\left(\frac{\sqrt{6}}{\Delta }\left(x+\frac{s}{2}\right)\right)-\mathrm{erf}\left(\frac{\sqrt{6}}{\Delta }\left(x-\frac{s}{2}\right)\right)\right].$$

$${\mathcal{R}}_2(r)=\frac{1}{\pi R^2}\frac{6}{\pi {\Delta }^2}\int \int H\left(R-\sqrt{y^{\prime 2}+z^{\prime 2}}\right)\exp \left(-6\frac{(y-y^{\prime }{)}^2+(z-z^{\prime }{)}^2}{{\Delta }^2}\right)d{y}^{\prime }d{z}^{\prime }$$

$$l(r)=\frac{1}{2}\rho \pi R^2{C}_T^{\prime }{u}_d^2{\mathcal{R}}_2(r)2\pi r.$$

$$\Gamma (r)=\frac{l(r)}{\rho V}=\frac{{\pi }^2{R}^2{C}_T^{\prime }{u}_d^2{\mathcal{R}}_2(r)}{\Omega }.$$

$${\ell }_{\text{adv}}=u_d\frac{2\pi }{\Omega }.$$

$$\eta (r)=\frac{d\Gamma }{dr}\frac{1}{{\ell }_{\text{adv}}}=\frac{\pi R^2{C}_T^{\prime }{u}_d}{2}\frac{d{\mathcal{R}}_2}{dr}.$$

$$\frac{\delta u(r)}{U_{\infty }}=\frac{1}{2U_{\infty }}{\int }_r^{\infty }\eta (r)dr=-\pi R^2\frac{C_T^{\prime }}{4}\frac{u_d}{U_{\infty }}{\mathcal{R}}_2(r).$$

$$\frac{u_d}{U_{\infty }}$$

$$\underset{\Delta \to 0}{\lim }\frac{u_d}{U_{\infty }}={\left(1+\frac{C_T^{\prime }}{4}\right)}^{-1}={\left(1+\frac{a}{1-a}\right)}^{-1}=(1-a),$$

$$u_d^{\text{corr}}=M\int \mathcal{R}(\mathbf{x})\mathbf{u}(\mathbf{x})\cdot \hat{\mathbf{x}}{d}^3\mathbf{x}.$$

$$\frac{u_d^{\text{corr}}}{U_{\infty }}=(1-a)=\frac{4}{4+C_T^{\prime }}.$$

$$\frac{u_d^{\text{corr}}}{U_{\infty }}=M\left(1-\pi R^2\frac{C_T^{\prime }}{4}\frac{u_d^{\text{corr}}}{U_{\infty }}{\int }_0^{\infty }{\mathcal{R}}_2^2(r)2\pi rdr\right),$$

$$M^{-1}=1+\frac{C_T^{\prime }}{4}\left(1-\pi R^2{\int }_0^{\infty }{\mathcal{R}}_2^2(r)2\pi rdr\right).$$

$$\pi R^2{\int }_0^{\infty }{\mathcal{R}}_2^2(r)2\pi rdr=1+\Delta \frac{d}{d\Delta }{\left[\pi R^2{\int }_0^{\infty }{\mathcal{R}}_2^2(r)2\pi rdr\right]}_{\Delta =0}+\mathcal{O}({\Delta }^2).$$

$$\frac{d}{d\Delta }\pi R^2{\int }_0^{\infty }{\mathcal{R}}_2^2(r)2\pi rdr=\sqrt{\frac{2}{3\pi }}\frac{1}{R}{\int }_{-\frac{R\sqrt{6}}{\Delta }}^{\infty }\left[1-\mathrm{erf}(\xi )\right]\xi e^{-{\xi }^2}\left(1+\frac{\Delta }{\sqrt{6}R}\xi \right)d\xi ,$$

$$M={\left(1+\frac{C_T^{\prime }}{4}\frac{1}{\sqrt{3\pi }}\frac{\Delta }{R}\right)}^{-1}.$$

Full text

1 Introduction

Blade-resolved simulations of wind turbines [Vijayakumar2016a], which can directly calculate the lift, drag, loading, and deflection of the blades, are computationally costly and thus mostly limited to simulations of a single wind turbine. In wind farm simulations, turbines are instead represented by actuator models, such as the actuator line model (ALM) [Sorensen2002a, Martinez2018a] and the actuator disk model (ADM) [Burton2011a]. Detailed simulations of wake characteristics and loading are provided by ALM, which models the lift and drag along the blades as a distribution of forces [Sorensen2002a, Martinez2018a]. However, the coarse grid resolution needed when simulating very large wind farms, which are becoming more prevalent, [Allaerts2018a] or applying optimal-control methods [Goit2015a] are often too costly for ALM. As a result, ADM, which distributes the thrust force of the turbine across the swept area of the rotor, remains a popular wind turbine representation [Iungo2017a, Stevens2018a, Lignarolo2016a].

For numerical stability of simulations using actuator models, the body force distribution (along the blades in ALM or across the rotor swept area in ADM) is typically filtered via a convolution kernel [Calaf2010a, Meyers2010a, Martinez2017b, Martinez2019a]. Unless the geometric footprint of the filtered force properly resolves the smallest scales of the shed vorticity, the induced velocity defect is under-predicted, and hence power generation is over-predicted [Martinez2017b, Martinez2019a]. A very finely resolved actuator disk will shed a vortex cylinder of infinitesimal thickness at the edge of the disk [Burton2011a], as shown in Figure 1a. At the coarser resolutions typically used in simulations, an ADM under-predicts the velocity gradient at the edge of the disk when compared to experiments [Lignarolo2016a], indicating under-resolved vorticity generation. With only a few points across the rotor, the under-prediction of the induced velocity distribution results in an over-prediction of power by 10% or more [Munters2017a]. An alternative ADM formulation that uses the Euler equations to directly model the infinitesimal vortex rings emanating from the actuator disk, mitigates these issues and has power coefficient errors of less than 1% [vanKuik2016a]. However, such simulations cannot model the turbulent breakdown of the wake or be embedded in large eddy simulations (LES). Clearly, corrections to ADM are needed to improve power coefficient predictions in LES of wind farms.

In this communication, we develop a correction factor for the power production of filtered actuator disks by examining the structure of the shed vorticity. First, we derive a vortex cylinder model for filtered actuator disks that predicts the disk-averaged velocity generated by the shed vorticity and compare these theoretical predictions to simulations. Then, we derive a correction factor for use in simulations that collapses the power coefficient measured at a range of grid resolutions and body force filter widths onto the power curve predicted by axial momentum theory. Recent work on using similar considerations to derive corrections for the ALM approach have led to the filtered actuator line model and improved results for ALM-LES [Martinez2017b, Martinez2019a].

2 The filtered actuator disk

We follow the filtered local actuator disk formulation of Meyers & Meneveau[Meyers2010a] and Calaf et al.[Calaf2010a], which writes the total thrust force $F$ in terms of the local thrust coefficient $C_{T}^{\prime}$, air density $\rho$, radius of the rotor swept area $R$, and disk averaged velocity $u_{d}$ as

[TABLE]

This thrust force is directed in the unit normal direction of the disk $\mathbf{\hat{x}}$ and is spatially distributed with a normalized (integrates to one) indicator function $\mathcal{R}(\mathbf{x})$

[TABLE]

The disk averaged velocity is expressed as a weighted average of the normal component of the velocity field $\mathbf{u}(\mathbf{x})$

[TABLE]

with $\mathcal{R}(\mathbf{x})$ used as the weighting function.

The normalized indicator function for a disk of finite thickness $s$ can be written with units of inverse volume as

[TABLE]

where $r$ is the radial distance along the disk written in terms of the transverse coordinate directions $r^{2}=y^{2}+z^{2}$ and $H(x)$ is the Heaviside function. In order to avoid numerical problems, such as unphysical oscillations and instabilities, that arise when applying sharp forces on grid points, a smoothed indicator function, obtained by convolving $\mathcal{I}(\mathbf{x})$ with a filtering kernel $G(\mathbf{x})$

[TABLE]

is typically used. In most implementations [Calaf2010a, Meyers2010a]

[TABLE]

is a Gaussian filtering kernel whose second moment is the same as a box filter of size $\Delta$ [Pope2000a]. The filter width $\Delta=\alpha h$ is the product of a factor of order unity $\alpha$ and the effective grid size111The effective grid size is expressed as the geometric mean $h=(\Delta x\Delta y\Delta z)^{1/3}$ of the grid spacings in some formulations [Lignarolo2016a]. $h=\sqrt{\Delta x^{2}+\Delta y^{2}+\Delta z^{2}}$, where $\Delta x$, $\Delta y$, and $\Delta z$ are the grid spacings.

For the Gaussian-filtered indicator function $\mathcal{R}(\mathbf{x})$, there exists a convenient decomposition that splits the indicator function into a product of normal and radial components

[TABLE]

The normal component

[TABLE]

has the analytic solution

[TABLE]

The radial component

[TABLE]

can be computed numerically or evaluated analytically when $\Delta/R<<1$ (see Section 4). This decomposition proves useful for implementing the smoothed indicator function (5) in simulations and computing the vortex cylinder model discussed in the following section.

3 A vortex cylinder theory for filtered actuator disks

To develop the vortex cylinder theory for filtered actuator disks, we use results from blade element methods and the vortex cylinder model [Burton2011a] with modifications to include the indicator function from the filtered actuator disk. For each radial element of the disk, the lift on the blade per unit radius is given by

[TABLE]

For a large rotational speed $\Omega$, the relative velocity over the blade is $V=\Omega r$. Applying the Kuta-Joukowsky theorem [Milne-Thomson1973a] with this relative velocity, the circulation distribution is

[TABLE]

The annular lift and circulation distribution will shed concentric semi-infinite vortex cylinders, as shown in Figure 1b. The circulation strength per unit radius $d\Gamma/dr$ is shed once per rotation of the rotor. During one revolution of the rotor, the wake will travel a distance $\ell_{\text{adv}}$. Assuming that the wake moves uniformly at the disk velocity $u_{d}$, this advection distance is

[TABLE]

Therefore, the shed circulation strength per unit area is

[TABLE]

For a single semi-infinite vortex cylinder with infinitesimal thickness and circulation strength per unit area of $\eta(r)=\eta_{0}\delta(r-R)$, shown in Figure 1a, the only non-vanishing component of the axial induced velocity is within the cylinder itself [Gibson1974a]. At the start of the disk, the axial induced velocity within the cylinder is uniformly $\eta_{0}/2$ [Gibson1974a]. For the circulation distribution in (14), we can superimpose the induced velocities for each of the shed concentric semi-infinite vortex cylinders. At radial location $r$, the only contribution to the induced velocity comes from cylinders with a radius $r^{\prime}$ greater than $r$. Integrating the contributions and noting that $\lim_{r\rightarrow\infty}\mathcal{R}_{2}(r)=0$, the induced velocity $\delta u(r)$ at the rotor plane and along the radius is given by

[TABLE]

As in  (3), the disk-averaged velocity can now be found by integrating the velocity at the disk $U_{\infty}+\delta u(r)$ weighted by the radial indicator function

[TABLE]

In the limit of a vanishing filter width $\Delta\rightarrow 0$, the disk-averaged velocity of a filtered actuator disk (17) tends to the solution of standard axial momentum theory, where the disk averaged velocity $u_{d}/U_{\infty}=(1-a)$ [Burton2011a] and local thrust coefficient $C_{T}^{\prime}=4a(1-a)$ [Calaf2010a, Meyers2010a] are a function of the induction factor $a$. In this limit, the filtered indicator function tends to a Heaviside function, i.e. $\lim_{\Delta\rightarrow 0}\mathcal{R}_{2}(r)=\pi^{-1}R^{-2}H(R-r)$. Therefore, the disk-averaged velocity goes to

[TABLE]

as predicted by axial momentum theory [Burton2011a].

The disk-averaged velocity in (17) is compared to simulations of actuator disks under uniform inflow using the pseudo-spectral code LESGO [Stevens2018a]. The inflow velocity $U_{\infty}$ is prescribed using a fringe region forcing that spans 25% of the the domain [Stevens2014a]. The actuator disk with diameter $D=2R$ is placed in a domain $L_{x}=15.36D$ long with a cross-section with dimensions $L_{y}=L_{z}=5.76D$. The center of the disk is located $3D$ from the inlet of the domain and in the center of the cross-section. All simulations have the same grid cell aspect ratio of $\Delta x=\Delta y=4\Delta z$. The Smagorinsky subgrid scale model with $C_{s}=0.16$ is used for numerical stability; however, the results are not affected by the value of $C_{s}$. Molecular viscosity is neglected. A number of simulations at various $C_{T}^{\prime}$, $\Delta$, and $h$ are run and the power coefficient $C_{P}=C_{T}^{\prime}(u_{d}/U_{\infty})^{3}$ [Meyers2010a] is computed from the simulation.

Despite the large domain, the simulations experience some blockage effect with the velocity slightly accelerating at the edges of the domain. To correct for this blockage, we assume that the streamwise pressure gradient will distribute the total mass flow rate lost at the rotor disk $\rho\pi R^{2}(U_{\infty}-u_{d})$ equally across the cross section of the domain; i.e. the velocity far from the rotor is $U_{\infty}+(U_{\infty}-u_{d})\pi R^{2}/(L_{z}L_{y})$, very close to what we observe in simulations. To obtain $u_{d}$ at the rotor, we need to subtract the extra velocity, and so the resulting corrected disk-averaged velocity is $u_{d}-(U_{\infty}-u_{d})\pi R^{2}/(L_{y}L_{z})$. This correction is small; even for the largest local thrust coefficient $C_{T}^{\prime}=2$, the blockage terms correction on the power coefficient is approximately 3%, which is much smaller than the largest filter-related power coefficient error.

The vortex cylinder theory described herein is compared to axial momentum theory and simulations with various filter widths and grid resolutions in Figure 2a. All simulations over-predict the power coefficient compared to axial momentum theory. The magnitude of this over-prediction is primarily dependent on the filter width $\Delta$ and relatively insensitive to resolution; i.e. for a given filter width $\Delta=\alpha h$, all simulations collapse onto the same curve regardless of $\alpha$. Furthermore, the power curves predicted by the filtered vortex cylinder model largely agree with the simulations, with some discrepancy at coarse resolutions with large local thrust coefficients.

4 Correction factor for simulations

The theory in Section 3 provides a means for correcting the disk-averaged velocity computed in simulations. Instead of directly averaging using the filtered indicator function, as in (3), we introduce a correction factor $M$ to compute the disk-averaged velocity as

[TABLE]

The correction factor $M$ is found by ensuring that the disk-averaged velocity for the force satisfies the momentum equation

[TABLE]

Replacing $\int\!\mathcal{R}(\mathbf{x})\,\mathbf{u}(\mathbf{x})\cdot\mathbf{\hat{x}}\,d^{3}\mathbf{x}$ in (19) with (16) expressed using the corrected disk-averaged velocity $u_{d}^{\text{corr}}$

[TABLE]

and replacing $u_{d}^{\text{corr}}/U_{\infty}$ with (20), we can solve for the correction factor

[TABLE]

Figure 2b compares $C_{P}$ computed using $u_{d}$ from simulations evaluated as in (19) and with M given by (22) to the result of axial momentum theory. Except for the correction factor, the simulations are identical to those of Section 3. As already noted in Section 3, for the uncorrected simulations, the power coefficient depends on the filter width and in cases with large filter widths, greatly exceeds the Betz limit. After applying the correction factor in  (19), the dependence on the filter width is eliminated. Furthermore, the power measurements collapse onto a single curve that differs only slightly from the axial momentum theory predictions.

For small filter sizes, we can expand the integral in (22) using a Taylor series expansion in $\Delta$ around $\Delta=0$ and recalling that the integral is unity at $\Delta=0$

[TABLE]

When $\Delta/R<<1$, the curvature of the disk is negligible when evaluating the convolution (10), and the solution for a Gaussian filter is equal to the one-dimensional convolution of the Gaussian with the Heaviside function, i.e. $\mathcal{R}_{2}(r)=\frac{1}{\pi R^{2}}\frac{1}{2}\left[1-\mathrm{erf}\left(\sqrt{6}\frac{r-R}{\Delta}\right)\right]$. Replacing this equation into (23), making a change of variables $\xi=\sqrt{6}(r-R)/\Delta$

[TABLE]

and evaluating at $\Delta=0$ yields $\sqrt{\frac{2}{3\pi}}\frac{1}{R}\int_{-\infty}^{\infty}\left[1-\mathrm{erf}(\xi)\right]\xi e^{-\xi^{2}}\,d\xi=-\frac{1}{R\sqrt{3\pi}}$. Therefore, for small $\Delta/R$ the correction factor can be written as

[TABLE]

The approximate correction factor is compared to numerical integration of (22) in Figure 3, showing that the approximation of $M$ using the Taylor series expansion of the integral in (22) is valid up to $\Delta/R=1.25$ and widely applicable for most LES grids.

5 Conclusions

The filtered ADM, which is often used in simulations of large wind farms, over-predicts power production, particularly at coarse grid resolutions. This discrepancy arises because the filtered body force does not adequately resolve the semi-infinite vortex cylinder shed behind the actuator disk, which in turn under-estimates the velocity defect induced at the rotor plane. The theory derived in this communication accurately predicts the disk-averaged velocity induced by the concentric semi-infinite vortex cylinders shed behind the turbine when compared to numerical simulations. The theory then provides a means for correcting the disk-averaged velocity in LES, which return the expected power curve from axial-momentum theory for various grid resolutions and filter widths.

Acknowledgements

The authors acknowledge funding from the National Science Foundation (grant CMMI 1635430). Computations made possible by the Maryland Advanced Research Computing Center (MARCC).