Journal cover Journal topic
Wind Energy Science The interactive open-access journal of the European Academy of Wind Energy
Journal topic
Wind Energ. Sci., 4, 369-383, 2019
https://doi.org/10.5194/wes-4-369-2019
Wind Energ. Sci., 4, 369-383, 2019
https://doi.org/10.5194/wes-4-369-2019

Research articles 28 Jun 2019

Research articles | 28 Jun 2019

# A vortex-based tip/smearing correction for the actuator line

A vortex-based smearing correction
Alexander R. Meyer Forsting, Georg Raimund Pirrung, and Néstor Ramos-García Alexander R. Meyer Forsting et al.
• DTU Wind Energy, Technical University of Denmark, Frederiksborgvej 399, 4000 Roskilde, Denmark
Abstract

The actuator line (AL) was intended as a lifting line (LL) technique for computational fluid dynamics (CFD) applications. In this paper we prove – theoretically and practically – that smearing the forces of the actuator line in the flow domain forms a viscous core in the bound and shed vorticity of the line. By combining a near-wake representation of the trailed vorticity with a viscous vortex core model, the missing induction from the smeared velocity is recovered. This novel dynamic smearing correction is verified for basic wing test cases and rotor simulations of a multimegawatt turbine. The latter cover the entire operational wind speed range as well as yaw, strong turbulence and pitch step cases. The correction is validated with lifting line simulations with and without viscous core, which are representative of an actuator line without and with smearing correction, respectively. The dynamic smearing correction makes the actuator line effectively act as a lifting line, as it was originally intended.

1 Introduction

The actuator line (AL) technique developed by is a lifting line (LL) representation of the wind turbine rotor suitable for computational fluid dynamics (CFD) simulations. It captures transient physical features like shed and trailed vorticity (including root/tip vortices), without the computational cost associated with resolving the full rotor geometry. Thus, the AL model enables large-eddy simulations (LES) of large wind farms in realistic, turbulent atmospheric boundary layers .

## 1.1 The vortex smearing hypothesis

noticed the similarity between the velocities induced across an actuator line and those predicted by a viscous vortex core model. These models include the limiting effect of viscous shear forces on the induced velocities around vortex cores. A similar comparison of the swirl velocities about an infinite vortex line is shown in Fig. 1 – here with a Lamb–Oseen vortex core model (Lamb1932; Oseen1911). Without viscosity (inviscid) the velocities approach infinity towards the vortex centre. The startling agreement between the Lamb–Oseen and AL velocities was first demonstrated by . The Gaussian body force smearing in the AL technique thus produces similar swirl velocities to a viscous vortex. Ignoring viscous effects, the AL should, in principle, induce the same velocities as a LL – equivalent to the inviscid solution. The missing induced velocity in the AL model (shaded area in Fig. 1) can be approximated following as follows:

$\begin{array}{ll}\mathrm{\Delta }{v}_{\mathit{\theta }}\left(r\right)& =\stackrel{\text{inviscid}}{\overbrace{\frac{\mathrm{\Gamma }}{\mathrm{2}\mathit{\pi }r}}}-\stackrel{\text{viscous core}}{\overbrace{\frac{\mathrm{\Gamma }}{\mathrm{2}\mathit{\pi }r}\left[\mathrm{1}-\mathrm{exp}\left(-{r}^{\mathrm{2}}/{\mathit{ϵ}}^{\mathrm{2}}\right)\right]}}\\ \text{(1)}& & =\frac{\mathrm{\Gamma }}{\mathrm{2}\mathit{\pi }r}\mathrm{exp}\left(-{r}^{\mathrm{2}}/{\mathit{ϵ}}^{\mathrm{2}}\right),\end{array}$

where Γ represents the vortex line's circulation, r is the distance from the vortex core and ϵ is the length scale used in the force smearing. This formulation can be split into an inviscid and viscous/smearing contribution:

$\begin{array}{ll}& \mathrm{\Delta }{v}_{\mathit{\theta }}\left(r\right)=\stackrel{\text{inviscid}}{\overbrace{{v}_{\mathit{\theta }}\left(r\right)}}\underset{\text{smearing}}{\underbrace{{f}_{\mathit{ϵ}}\left({r}_{\mathit{ϵ}}\right)}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{with}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{v}_{\mathit{\theta }}\left(r\right)=\frac{\mathrm{\Gamma }}{\mathrm{2}\mathit{\pi }r},\\ \text{(2)}& & {f}_{\mathit{ϵ}}\left({r}_{\mathit{ϵ}}\right)=\mathrm{exp}\left(-{r}_{\mathit{ϵ}}^{\mathrm{2}}\right),\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{r}_{\mathit{ϵ}}=\frac{r}{\mathit{ϵ}}.\end{array}$

If this viscous behaviour of the force smearing in AL simulations was limited to the bound vortex representing the blade, it would not influence the blade forces as long as the blade was straight. However, argued that the trailing vortices (in the wake) exhibit the same viscous core, as they originate from the bound vortex. Hence, the wake of an AL induces lower velocities at the blade than a LL. The missing velocity can be estimated from the viscous core equivalence and thus correct the velocities at the blade. This mostly impacts blade forces by changing the angle of attack at the blade sections.

Figure 1Distribution of the tangential velocity component in a plane orthogonal to an infinite vortex line (along x) obtained from either an inviscid or viscous (Lamb–Oseen) theoretical vortex and an actuator line (AL) CFD simulation.

## 1.2 Contributions of this paper

corrected AL simulations of a rectangular wing and two rotors with different aspect ratios by recuperating the missing induced velocity introduced by the viscous core. For all of their simulations they were able to show the beneficial effect of the correction on the blade load distribution – represented by more physical behaviour, especially towards the tip and root. However, their implementation of the correction did not fully couple the flow field with the blade forces and the induction correction.

The major contributions of this paper are as follows:

• The development of a tuning-free, dynamic and numerically robust smearing correction, which is fully coupled to the AL model.

• A theoretical proof of the force smearing – vortex core equivalence.

• Proof of the vortex core inheritance in trailed vorticity.

• The confirmation of the missing velocity assumption by comparing LL simulations with/without viscous core and AL results with/without correction.

The test cases include constantly and elliptically loaded wings as well as rotor simulations of a multimegawatt turbine covering the entire operational wind speed range. As the AL model is especially attractive for wind farm simulations, the focus here is on coarsely resolved ALs. The correct dynamic behaviour of the new correction is verified through yawed inflow and pitch step simulations.

2 Proof of force smearing – vortex core equivalence

The equivalence between the velocity field induced by an AL and a viscous vortex can be derived directly from the incompressible Navier–Stokes equations. This proof follows the approach by that successfully connected an AL's vorticity field to its force projection. Starting by taking the curl of the incompressible momentum equation, the vorticity transport equation is obtained ($\mathbit{\omega }=\mathrm{\nabla }×\mathbit{u}$):

$\begin{array}{}\text{(3)}& \underset{\begin{array}{c}=\mathrm{0}\\ \text{steady}\end{array}}{\underbrace{\frac{\partial \mathbit{\omega }}{\partial t}}}+\left(\mathbit{u}\mathrm{\nabla }\right)\mathbit{\omega }=\underset{\begin{array}{c}=\mathrm{0}\\ \mathrm{2}D\end{array}}{\underbrace{\left(\mathbit{\omega }\mathrm{\nabla }\right)\mathbit{u}}}+\underset{\begin{array}{c}=\mathrm{0}\\ \text{inviscid}\end{array}}{\underbrace{\mathit{\nu }{\mathrm{\nabla }}^{\mathrm{2}}\mathbit{\omega }}}+\mathrm{\nabla }\frac{\mathbit{f}}{\mathit{\rho }},\end{array}$

where ν is the viscosity, ρ is density and f represents the body forces from the AL. Away from the root and tip, the flow around a high-aspect-ratio blade is nearly two-dimensional, as the span-wise flow is negligible. Viscous effects are disregarded in light of the large Reynolds numbers encountered. Furthermore the relationship between the body force and flow field becomes quasi-steady assuming the flow is attached. Cancelling the respective terms and noting that in two-dimensional flow (yz plane) $\mathbit{\omega }={\mathit{\omega }}_{x}{\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}$ and $\mathrm{\nabla }=\left(\mathrm{0},\frac{\partial }{\partial y},\frac{\partial }{\partial z}\right)$:

$\begin{array}{}\text{(4)}& \left(\mathbit{u}\mathrm{\nabla }\right){\mathit{\omega }}_{x}{\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}=\mathrm{\nabla }\frac{\mathbit{f}}{\mathit{\rho }}.\end{array}$

Assuming the drag to be negligible, the force – in the form of lift – exerted by the AL on the flow in terms of its circulation Γ becomes

$\begin{array}{}\text{(5)}& & \mathbit{f}=-{\mathbit{f}}_{\text{aero}}g\left(r\right)\text{(6)}& & {\mathbit{f}}_{\text{aero}}=L=\mathit{\rho }\mathbit{u}×\mathrm{\Gamma }{\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}\phantom{\rule{2em}{0ex}}g\left(r\right)=\frac{\mathrm{1}}{\mathit{\pi }{\mathit{ϵ}}^{\mathrm{2}}}\mathrm{exp}\left(-{r}^{\mathrm{2}}/{\mathit{ϵ}}^{\mathrm{2}}\right).\end{array}$

Here g represents a two-dimensional Gaussian force projection with r indicating the distance from the AL. Inserting these expression into Eq. (4) and exploiting standard matrix transformation and mass conservation1

$\begin{array}{}\text{(7)}& \left(\mathbit{u}\mathrm{\nabla }\right){\mathit{\omega }}_{x}{\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}& =\mathrm{\nabla }\left(\mathrm{\Gamma }g{\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}×\mathbit{u}\right)=\left(\mathbit{u}\mathrm{\nabla }\right)\mathrm{\Gamma }g{\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}\text{(8)}& \left(\mathbit{u}\mathrm{\nabla }\right){\mathit{\omega }}_{x}& =\left(\mathbit{u}\mathrm{\nabla }\right)\mathrm{\Gamma }g.\end{array}$

Due to mass conservation the u term can be inverted; this gives a direct relationship between the force projection and vorticity:

$\begin{array}{}\text{(9)}& {\mathit{\omega }}_{x}=\mathrm{\Gamma }g=\frac{\mathrm{\Gamma }}{\mathit{\pi }{\mathit{ϵ}}^{\mathrm{2}}}\mathrm{exp}\left(-{r}^{\mathrm{2}}/{\mathit{ϵ}}^{\mathrm{2}}\right).\end{array}$

As the body force is axially symmetric, the vorticity only induces tangential velocities

$\begin{array}{ll}{\mathit{\omega }}_{x}\left(r\right)& =\frac{\mathrm{1}}{r}\left(\frac{\partial r{u}_{\mathit{\theta }}}{\partial r}-\underset{\begin{array}{c}=\mathrm{0}\\ \text{axisymmetry}\end{array}}{\underbrace{\frac{\partial {u}_{r}}{\partial \mathit{\theta }}}}\right)\\ \text{(10)}& & ⇒{u}_{\mathit{\theta }}=\frac{\mathrm{1}}{r}\underset{\mathrm{0}}{\overset{r}{\int }}r{\mathit{\omega }}_{x}\left(r\right)\mathrm{d}r.\end{array}$

Inserting Eq. (9) and integrating gives the swirl velocity induced by a smeared body force

$\begin{array}{}\text{(11)}& {u}_{\mathit{\theta }}=\frac{\mathrm{\Gamma }}{\mathrm{2}\mathit{\pi }r}\left[\mathrm{1}-\mathrm{exp}\left(-{r}^{\mathrm{2}}/{\mathit{ϵ}}^{\mathrm{2}}\right)\right].\end{array}$

This expression equals that of the Lamb–Oseen vortex, only with the viscous core radius replaced by the smearing coefficient2. This marks the theoretical confirmation of the observations by , which additionally indicates that a viscous core behaviour with an AL requires inviscid, two-dimensional and locally steady flow conditions.

3 Numerical methodology

## 3.1 Actuator line simulations

The discretized incompressible Navier–Stokes equations are solved using DTU's CFD code EllipSys3D . The flow is iteratively solved at each time instant by the SIMPLE algorithm . Depending on the turbulence model either the third-order accurate QUICK (Leonard1979) scheme (RANS) or a fourth-order CDS scheme (LES) discretizes the convective terms. As the flow variables are located at the cell centres, a modified Rhie and Chow algorithm avoids pressure–velocity decoupling. Further details regarding the numerical techniques are given in . For all comparisons with the LL code, the RANS equations are solved using the kω shear-stress transport turbulence closure of . Only the turbulent inflow cases in Sect. 5.2.4 are computed using the DES technique of . The AL model was implemented by in EllipSys3D. We employ a version utilizing three-dimensional Gaussian force projection, which follows the original formulation of . As the AL model is especially attractive for wind farm simulations, the focus here is on coarsely resolved ALs, with either 9 or 19 sections (Ns) along the blade. They are uniformly spaced and discretize the blade starting from the root at 1.5 m to the tip at 63 m. The smearing length scale is connected to the number of sections, such that $\mathit{ϵ}=\mathrm{2}R/\left({N}_{s}+\mathrm{1}\right)-R$ defining the rotor radius – which ensures that the forces in the domain change smoothly between sections (Nathan2018). The tower and nacelle are not modelled.

Figure 2Numerical box domain with a structured mesh and uniform spacing around the rotor at its centre. Only every eighth grid point is shown.

The numerical domain for the rotor simulations is discretized in a verified, standard manner . It consists of a box with 25R side length that contains an inner box with a uniformly spaced refined mesh of 3.2R edge length at its centre surrounding the rotor (see Fig. 2). To capture the velocity gradients around the AL correctly the mesh spacing is $\mathrm{\Delta }x=R/\mathrm{40}$. This is twice the recommended minimum ; however, it delivers more accurate angle of attack estimates at the section centres . In total, 256 cells discretize the flow domain along each dimension, resulting in 16.8×106 degrees of freedom. All variables, except pressure and its correction, which necessitate special treatment (Sørensen1995), obey symmetry conditions on the lateral boundaries, whereas at the inflow and outflow faces they follow Dirichlet and Neumann conditions, respectively. The wing test cases follow the same approach only with Ns=32 and an inner box edge length of 3b, where b is the wing's half-span. This results in 80 cells along each dimension and 5.1×105 degrees of freedom. To ensure that the blade tip remains inside a single cell during one time step $\mathrm{\Delta }t<\mathrm{\Delta }x/\left(\mathrm{\Omega }R\right)$ – with mesh spacing Δx and rotational speed Ω. Without rotation, the term ΩR is replaced by the advection speed of the wake. The kinematic viscosity and air density are kept constant at 1.789 $×{\mathrm{10}}^{-\mathrm{5}}$ kg m−1 s−1 and 1.225 kg m−3, respectively. Simulations are stopped when the thrust residual reaches $\mathrm{1}×{\mathrm{10}}^{-\mathrm{5}}$.

The sensitivity of the rotor thrust to the domain size, time step and grid size is explored in Fig. 3. The length of the domain edges is doubled to 50R, the time step is halved with respect to a set-up obeying the method described above. A simulation of the NREL 5-MW at 8 m s−1 with either 40 or 60 grid cells along the rotor depending on the smearing length scale acts as a reference. With $\mathit{ϵ}=R/\mathrm{10}$ and R∕20, this represents 4 and 3 times the recommended resolution, respectively . Although non-zero, the sensitivity of the results is acceptable in code comparison and should impact AL simulations with and without correction similarly.

Figure 3Thrust sensitivity of the NREL 5-MW AL simulations at a 8 m s−1 wind speed with respect to grid size, doubling the domain size and halving the time step at two smearing length scales (${N}_{s}=\mathit{\left\{}\mathrm{9},\mathrm{19}\mathit{\right\}}$, ${T}_{\text{ref}}=\mathit{\left\{}\mathrm{4.20},\mathrm{4.06}\mathit{\right\}}×{\mathrm{10}}^{\mathrm{5}}$ N).

## 3.2 Free-wake lifting line rotor simulations

The in-house solver MIRAS has been employed to perform the free-wake lifting line simulations. MIRAS is a multi-fidelity computational vortex model, which is mainly used for predicting the aerodynamic behaviour of wind turbines and their wakes. It has been developed at DTU over the last decade and been extensively validated for small to large size wind turbine rotors by .

The free-wake vortex method essentially models the wake of a wind turbine using a bundle of infinitely thin vortex filaments. To avoid numerical singularities, a viscous core must be introduced, which represents a more physical distribution of the velocities induced by each vortex filament, desingularizing the Biot–Savart law near the centre of the filament. The velocity induced by each one of the elements is obtained directly by evaluating the Biot–Savart law, and by summing the velocity induced by all filaments, the total wake induction is obtained as follows:

$\begin{array}{ll}& \mathbit{u}\left({\mathbit{x}}_{i}\right)=\sum _{j=\mathrm{1}}^{N}{K}_{ij}\frac{{\mathit{\gamma }}_{j}}{\mathrm{4}\mathit{\pi }}\frac{{\mathbit{t}}_{j}×{\mathbit{r}}_{ij}}{{r}_{ij}^{\mathrm{3}}},\\ \text{(12)}& & \text{where}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{K}_{ij}=\frac{{r}_{ij}^{\mathrm{2}}}{{\left({\mathit{\epsilon }}_{j}^{\mathrm{2}z}+{r}_{ij}^{\mathrm{2}z}\right)}^{\mathrm{1}/z}},\end{array}$

N is total number of filaments that form the wake, ${\mathbit{r}}_{ij}={\mathbit{x}}_{i}-{\mathbit{y}}_{j}$ is the distance vector from the vortex element yj to the evaluation point xi, γj is the circulation of the filament, tj is the unit orientation vector of the jth filament and ${r}_{ij}=|{\mathbit{r}}_{ij}|$, εj is the vortex core radius of the filament, and z defines the cut-off velocity profile where the Lamb–Oseen model (Lamb1932; Oseen1911), z=2, has been employed.

A viscous core model is applied to emulate the effect of viscosity by changing the vortex core radius as a function of time :

$\begin{array}{}\text{(13)}& {\mathit{\epsilon }}_{i}\left(t\right)=\sqrt{\mathrm{4}{\mathit{\alpha }}_{v}{\mathit{\delta }}_{v}\mathit{\nu }{t}_{i}}+{\mathit{\epsilon }}_{\mathrm{0}},\end{array}$

where αv is a constant set to 1.25643 , ν is the kinematic viscosity and ti is the time elapsed since the generation of the ith filament. In order to represent the diffusive timescales, the viscous core radius is set to change with the vortex age by adding a turbulence eddy viscosity, δv, first proposed by , and in this work set to $\mathrm{1}×{\mathrm{10}}^{-\mathrm{3}}$. To avoid the singular behaviour of newly released vortex elements, an initial core radius, ε0, is introduced. In accordance with , who found that a small core radius is necessary to have flow convergence, a core radius of 0.1 % of the local chord at the release station is used.

For the sake of the present study, two different approaches to compute the angle of attack have been followed.

• Inviscid (LL), where the non-regularized Biot–Savart law is used to compute the induction from the wake filaments at the quarter-chord location. This is the standard method used in a lifting line solvers.

• Viscous (LL+core), where the regularized Biot–Savart law is used to compute the induction at the quarter-chord location. A viscous core with a radius equal to the actuator line smearing factor is used for a direct comparison of the methods.

This enables a double validation of the models. On the one hand the corrected AL simulations can be validated against the LL calculations, and on the other hand the raw AL model, without tip correction, can be compared against the LL+core simulations which include the smearing effect in the free-wake model.

4 Tip/smearing correction for the actuator line

Applying the velocity correction methodology introduced in Sect. 1.1 in three-dimensional space yields a velocity correction vector. The viscous core behaviour of the AL bound vorticity – proven in Sect. 2 to originate from the force smearing – is inherited by the trailing vortices, as will be demonstrated in Sect. 5.1.1. Therefore, the induction from the trailed vorticity at the blade is lower than without force smearing. Figure 4 shows the path of trailed vorticity shed from in-between two sections of a blade with a strength of $\mathrm{\Delta }\mathrm{\Gamma }={\mathrm{\Gamma }}_{s}-{\mathrm{\Gamma }}_{s+\mathrm{1}}$ with

$\begin{array}{}\text{(14)}& {\mathrm{\Gamma }}_{s}=\frac{\mathrm{1}}{\mathrm{2}}\sqrt{{v}_{s}^{\mathrm{2}}+{w}_{s}^{\mathrm{2}}}{C}_{\mathrm{L}}\left(\mathit{\alpha }\right)c.\end{array}$

Here s defines the blade section index, CL is the sectional lift coefficient, c is the section chord and α is the angle of attack, which depends on the inflow angle in combination with blade pitch and twist at the section. The missing induction from this single trailed vortex at a point C is obtained by integrating along the vortex line

$\begin{array}{}\text{(15)}& {\mathbit{u}}^{*}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{f}_{\mathit{ϵ}}\mathbit{\delta }\stackrel{\mathrm{̃}}{\mathbit{u}}\phantom{\rule{0.25em}{0ex}}\text{d}l,\end{array}$

where l is the vortex following coordinate. Here $\mathbit{\delta }\stackrel{\mathrm{̃}}{\mathbit{u}}$ is the velocity induced by an infinitesimal element δl of the vortex line at point C, which is given by the Biot–Savart law:

$\begin{array}{}\text{(16)}& \mathbit{\delta }\stackrel{\mathrm{̃}}{\mathbit{u}}=\frac{\mathrm{\Delta }\mathrm{\Gamma }}{\mathrm{4}\mathit{\pi }}\frac{\mathbit{\delta }\mathbit{l}×\mathbit{x}}{|\mathbit{x}{|}^{\mathrm{3}}},\end{array}$

where x is the vector pointing from the element towards C. The smearing factor for this vortex element becomes

$\begin{array}{}\text{(17)}& {f}_{\mathit{ϵ}}=\mathrm{exp}\left(-\frac{\left(\mathbit{x}{\stackrel{\mathrm{^}}{\mathbit{e}}}_{⟂}{\right)}^{\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{2}}}\right).\end{array}$

The viscous core only acts in the plane orthogonal to the vortex element δl; hence, ${\stackrel{\mathrm{^}}{\mathbit{e}}}_{⟂}$ projects x onto this plane. This is different to using the distance, $|\mathbit{x}|$, as proposed, which violates the two-dimensional nature of the viscous core.

Figure 4Trailed vorticity path. The blade rotates in the xy plane and z points downstream. The vortex element δl with a strength of ${\mathrm{\Gamma }}_{s}-{\mathrm{\Gamma }}_{s+\mathrm{1}}$ is shed at r and transported downstream by the local velocity. The distance from the shedding location r to a point C along the blade is h ($h=r-{C}_{x}$), where δl induces tangential and axial velocities.

The total missing induction at a blade section s is obtained by summing the contribution from all trailed vortices. The number of trailed vortices Nv is directly related to the number of blade sections ${N}_{v}={N}_{s}+\mathrm{1}$. Discretizing the vortices in time, the missing induction at a certain blade section becomes

$\begin{array}{}\text{(18)}& {\mathbit{u}}_{s}^{*}=\sum _{v}^{{N}_{v}}\sum _{n}^{{N}_{t}}{f}_{s,v}^{n}\mathbf{\Delta }{\stackrel{\mathrm{̃}}{\mathbit{u}}}_{s,v}^{n}.\end{array}$

Here v denotes the trailed vortex index, n is the time index and Nt is the number of time steps. Note that n=1 is the most recently shed vortex element. As a tip/smearing correction should remain computationally inexpensive, numerically solving the Biot–Savart law in Eq. (16) to obtain $\mathbf{\Delta }{\stackrel{\mathrm{̃}}{\mathbit{u}}}_{s,v}^{n}$ is unfeasible. This would necessitate NtNvNs or Nt(Ns+1)Ns evaluations. An accurate, yet fast, alternative to solving the Biot–Savart law directly is the near-wake model (NWM) for trailed vorticity by , which also includes downwind convection. It performs well for dynamic flow cases and exhibits great numerical stability, as it was originally developed to enhance the aerodynamic accuracy of blade element momentum (BEM) models. Its formulation is based on a lifting line representation of the blade's trailed vorticity as depicted in Fig. 4 and approximates the induced velocities from a single trailed vortex line by two indicial functions. The velocity induced by a vortex element is given in the NWM as

$\begin{array}{}\text{(19)}& \mathbf{\Delta }{\stackrel{\mathrm{̃}}{\mathbit{u}}}_{s,v}^{n}=\left({X}_{s,v}^{n}+{Y}_{s,v}^{n}\right)\left[\begin{array}{c}\mathrm{0}\\ \mathrm{sin}\left({\mathit{\varphi }}^{n}\right)\\ -\mathrm{cos}\left({\mathit{\varphi }}^{n}\right)\end{array}\right],\end{array}$

with ϕ representing the helix angle of the vorticity shed in the CFD domain (see Fig. 4). The indicial functions take the following form:

$\begin{array}{ll}\left\{{X}_{s,v}^{n},{Y}_{s,v}^{n}\right\}& ={a}_{\left\{X,Y\right\}}\frac{{r}_{v}}{\mathrm{4}\mathit{\pi }{h}_{s}|{h}_{s}|}\mathrm{\Delta }{\mathrm{\Gamma }}_{v}^{n}{\mathit{\varphi }}_{s,v}^{{*}^{n}}\\ & \left[\mathrm{1}-\mathrm{exp}\left(-{b}_{\left\{X,Y\right\}}\frac{\mathrm{\Delta }{\mathit{\beta }}^{{*}^{n}}}{{\mathit{\varphi }}_{s,v}^{{*}^{n}}}\right)\right]\\ \text{(20)}& & \mathrm{exp}\left(-{b}_{\left\{X,Y\right\}}\sum _{i}^{n-\mathrm{1}}\frac{\mathrm{\Delta }{\mathit{\beta }}^{{*}^{i}}}{{\mathit{\varphi }}_{s,v}^{{*}^{i}}}\right).\end{array}$

The definitions of ${a}_{\left\{X,Y\right\}},{b}_{\left\{X,Y\right\}},{\mathit{\beta }}^{*}$ and ϕ* are those of . The indicial functions allow for the solution to be time-advanced by a mere multiplication, considerably reducing the model evaluations to NvNs+Nv. In the original formulation this removes the need for bookkeeping; however, as the smearing factor also changes with the position of the vortex element, all previously shed elements are advanced individually in this specific implementation. This is only an experimental feature for testing the smearing correction and should be simple to remove in a future, more practical implementation.

Following the lifting line formulation of the NWM shown in Fig. 4 (refer to Appendix A for a detailed mathematical description) the perpendicular distance from the vortex element to C becomes

$\begin{array}{ll}& {\mathbit{x}}_{⟂}=\frac{\mathbit{\delta }\mathbit{l}}{|\mathbit{\delta }\mathbit{l}|}×\mathbit{x}=\\ \text{(21)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}r\mathrm{cos}\mathit{\varphi }\left(\begin{array}{c}\mathrm{tan}\mathit{\varphi }\left(\mathit{\beta }\mathrm{cos}\mathit{\beta }-\mathrm{sin}\mathit{\beta }\right)\\ -\mathrm{tan}\mathit{\varphi }\left(-\mathrm{1}+h/r+\mathrm{cos}\mathit{\beta }+\mathit{\beta }\mathrm{sin}\mathit{\beta }\right)\\ -\mathrm{1}+\left(\mathrm{1}-h/r\right)\mathrm{cos}\mathit{\beta }\end{array}\right).\end{array}$

Thus, the smearing factor becomes

$\begin{array}{}\text{(22)}& {f}_{s,v}^{n}=\mathrm{exp}\left(-\frac{|{\mathbit{x}}_{\mathbf{⟂}}\left({r}_{v},{\mathit{\beta }}^{n},{h}_{s},{\mathit{\varphi }}^{n}\right){|}^{\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{2}}}\right).\end{array}$

When discretizing in time, βn is taken to be at the mid-point of the vortex element.

Finally the missing velocities computed in Eq. (18) correct the original velocities from the CFD simulations

$\begin{array}{}\text{(23)}& {\mathbit{u}}_{s}={\mathbit{u}}_{s}^{\text{CFD}}+{\mathbit{u}}_{s}^{*}.\end{array}$

Therefore, the correction influences the blade forces through the angle of attack and the velocity magnitude, although it is the former that dominates. It also changes the circulation at each blade section through Eq. (14) and, thus, the shed vorticity and its induction. Hence determining the correction velocity is an iterative procedure. The correction algorithm is executed after the flow field is solved and takes the following form:

1. Interpolate the velocity vector uCFD at the section centres from the CFD flow field.

2. Compute the helix angles ϕ, where ${\mathit{\varphi }}_{v}=-{\mathrm{tan}}^{-\mathrm{1}}\left(\frac{{w}_{v-\mathrm{1}}^{\text{CFD}}+{w}_{v}^{\text{CFD}}}{{v}_{v-\mathrm{1}}^{\text{CFD}}+{v}_{v}^{\text{CFD}}}\right)$ and ${\mathit{\varphi }}_{\left\{\mathrm{1},{N}_{v}\right\}}=-{\mathrm{tan}}^{-\mathrm{1}}\left(\frac{{w}_{\left\{\mathrm{1},{N}_{s}\right\}}^{\text{CFD}}}{{v}_{\left\{\mathrm{1},{N}_{s}\right\}}^{\text{CFD}}}\right)$.

3. Combine the CFD velocities with the respective correction from the previous time step n−1, such that ${\mathbit{u}}_{n}={\mathbit{u}}_{n}^{\text{CFD}}+{\mathbit{u}}_{n-\mathrm{1}}^{*}$.

4. Compute the smearing factor fϵ for all time steps, sections and elements (Eq. 22).

5. Determine the angle of attack and velocity magnitude from un to determine Γs (Eq. 14)3,

6. Compute the velocities from the newly released vortex element $\mathbf{\Delta }{\stackrel{\mathrm{̃}}{\mathbit{u}}}_{n}$ (Eq. 19).

7. At the first iteration of each time step, advance the previous elements in time.

8. Compute the velocity correction at the current time step ${\mathbit{u}}_{n}^{*}$ (Eq. 18).

9. Update the velocity at the sections with some form of relaxation ${\mathbit{u}}_{n}={\mathbit{u}}_{n}^{\text{CFD}}+{\mathbit{u}}_{n}^{*}$.

10. Repeat steps 5–9 until convergence is reached.

We use the technique by established for the NWM to accelerate and ensure its convergence. Furthermore, the activation of the correction is delayed until the starting vorticity of the rotor has been transported at least one blade length away from the rotor plane. This enhances its numerical stability, as induction has already built up at the blades by its time of activation.

5 Results

## 5.1 Basic wing test cases

To verify the smearing hypothesis (Sect. 1.1) and the novel smearing correction (Sect. 4) two basic wing flow cases with known theoretical solutions are modelled using CFD. Either a rectangular or an elliptic wing is represented by an AL as shown in Fig. 5, where the coordinate system is unchanged from the definition in Fig. 4. The AL is discretized in uniformly spaced sections in-line with the underlying flow grid, and the smearing parameter is twice the section width, which ensures a continuous force distribution along the wing. The common simulation parameters are given in Table 1. Unless specifically stated the sectional lift coefficient CL=1 and drag is zero along the wing, independent of the angle of attack. The chord of the rectangular wing is set to 1 m and the elliptical chord distribution is

$\begin{array}{}\text{(24)}& c\left(x\right)={c}_{\mathrm{0}}\sqrt{\mathrm{1}-{\left(\frac{x-\left(b+{r}_{v=\mathrm{1}}\right)}{b}\right)}^{\mathrm{2}}},\end{array}$

with the root chord c0=4m. All simulations are performed within the same computational domain, defined in Sect. 3.1.

Figure 5Definition of the wing test cases with either a rectangular or elliptic planform. Vortices are trailed in-between sections and the actuator line forces are computed and exerted at the sections' centres.

Table 1Input parameters common to both rectangular and elliptic wing simulations.

The theoretical predictions of the velocity field are achieved by representing the vortex system of Fig. 5 with vortex filaments. The velocity induced by a filament with a viscous vortex core at an arbitrary point C is

$\begin{array}{}\text{(25)}& \mathbit{u}& ={f}_{\mathit{ϵ}}\left({x}_{⟂}\right)\frac{\mathrm{\Gamma }}{\mathrm{4}\mathit{\pi }}\frac{\left({x}_{\mathrm{1}}+{x}_{\mathrm{2}}\right)\left({\mathbit{x}}_{\mathbf{1}}×{\mathbit{x}}_{\mathbf{2}}\right)}{{x}_{\mathrm{1}}{x}_{\mathrm{2}}+{\mathbit{x}}_{\mathbf{1}}\cdot {\mathbit{x}}_{\mathbf{2}}}\text{(26)}& {x}_{⟂}& =\frac{|{\mathbit{x}}_{\mathbf{1}}×{\mathbit{x}}_{\mathbf{2}}|}{|{x}_{\mathrm{2}}-{x}_{\mathrm{1}}|}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{and}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{x}_{i}=|{\mathbit{x}}_{i}|,\end{array}$

where x1 points from the start of the filament to C and x2 from its end. For a definition of fϵ(x) refer to Eq. (17); without viscous core ${f}_{\mathit{ϵ}}\left({x}_{⟂}\right)=\mathrm{1}$. The contribution from different segments is summed to give the overall velocity field. As the wing is lightly loaded, we assume that all vortex segments remain in the xz plane.

### 5.1.1 Trailed vorticity smearing

The vortex smearing hypothesis assumes the trailed vorticity inheriting the smeared velocity field from the bound vortex, which was confirmed theoretically by for straight wings. We test this further by simulating a rectangular wing without any correction. All vorticity is shed from the wing tips, creating the well-known horseshoe vortex. Hence, for the hypothesis to be valid, the velocity distribution in the plane orthogonal to the trailed vortices should be identical to that of the bound vortex .

Figure 6 compares the velocities induced by a rectangular wing predicted by an AL and three vortex segments (one bound, two trailed) for five different smearing parameters. Only half of the wing is presented, due to symmetry. Velocity distributions are shown for lines cutting the vortex segments at right angles for y=0. Clearly the velocity smearing is identical between trailed and bound vorticity, confirming the smearing hypothesis. Slight differences are linked to the numerical discretization of the Gaussian force projection and numerical diffusion.

Figure 6Velocities induced perpendicular to a rectangular wing predicted by an actuator line (AL) without correction and by vortex segments with a viscous core (Vortex) with different smearing parameters. Velocities are shown along lines cutting the bound and trailed vortices at right angles and y=0. Only half of the horseshoe vortex is depicted as ${x}^{\prime }=x-\left(b+{r}_{v=\mathrm{1}}\right)$.

### 5.1.2 Smearing correction verification

As mentioned in Sect. 4, the new smearing correction uses a lifting line representation of the trailed vorticity. Thus, the prediction of the velocity correction with our model or vortex segments should be identical. To simultaneously verify its numerical implementation, our model only receives the sampled velocities from the flow domain to compute the circulation at the sections. Furthermore, the body forces are not applied inside the domain to avoid influencing the trailed vortex paths, and the correction velocities are not added to the CFD velocities to keep the circulation unchanged. This holds the trailed vortices in the xz plane, simplifying the representation of the wake with vortex segments. The segments' circulation is exactly the same as in the smearing correction to avoid any numerical effects influencing the comparison. Figure 7 compares the velocity correction predicted by the analytical vortex segments and our model at each section along a rectangular and elliptic wing for different smearing factors (i.e. $\mathit{ϵ}=\mathrm{2}b/\left({N}_{v}-\mathrm{1}\right)$). With decreasing force smearing the velocity correction concentrates towards the tips, as the induced velocity gradients increase. Therefore, even at higher resolution the smearing correction remains significant, although more localized. Generally the model slightly over-predicts the missing induction at the wing, becoming more prominent with increasing resolution. At Nv=64 the difference reaches a maximum of 6.7 % (rectangular) and 1.7 % (elliptic) with respect to the inflow velocity. The average error does not breach 0.5 % in any case. The velocity jump towards the tip sections of the elliptical wing is related to the equidistant discretization of the wing .

Figure 7Analytical (An.) and corresponding model prediction of the velocity correction for varying force smearing ($\mathit{ϵ}=\mathrm{2}b/\left({N}_{v}-\mathrm{1}\right)$) along a rectangular and elliptic wing, where ${x}^{\prime }=x-\left(b+{r}_{v=\mathrm{1}}\right)$.

### 5.1.3 Coupled AL-smearing correction verification

The coupling between velocity correction and the flow domain is verified by comparing the corrected downwash at an elliptical wing to the theoretical expectation. The downwash should be constant along the wing and is given by

$\begin{array}{}\text{(27)}& {v}_{\text{th}}=-\frac{{\mathrm{\Gamma }}_{\mathrm{0}}}{\mathrm{4}b}=-\frac{{w}_{\mathrm{\infty }}{c}_{\mathrm{0}}{C}_{\mathrm{L}}}{\mathrm{8}b},\end{array}$

where Γ0 is the circulation at the wing root. Similar to the CL was not fixed, but instead followed the theoretical lift curve slope for thin airfoils CL=2π. For the wing to operate at a constant lift coefficient CL=1, its angle of attack needed to include the effect of the induced velocities:

$\begin{array}{}\text{(28)}& \mathit{\alpha }={\mathit{\alpha }}_{\text{eff}}+{\mathit{\alpha }}_{i}=\frac{{C}_{\mathrm{L}}}{\mathrm{2}\mathit{\pi }}+{\mathrm{tan}}^{-\mathrm{1}}\left(\frac{{c}_{\mathrm{0}}{C}_{\mathrm{L}}}{\mathrm{8}b}\right).\end{array}$

This represents a more rigorous test of the coupled system than prescribing the loading along the wing, as only the correct downwash leads to the desired, constant sectional lift coefficients.

Figure 8 shows the downwash predicted by AL simulations with different smearing parameters and active correction. The CFD components of the velocities are shown (vCFD) separately to emphasize the contribution of the correction to arrive at the correct, constant downwash of 1 m s−1. Clearly without the correction, the induced velocities are a function of the smearing factor and only arrive at the theoretically expected value for Nv=32. Including the correction greatly reduces the dependence of the downwash on the force smearing. The insufficient correction towards the tips feeds back to the equidistant discretization of the AL , which is linked to the uniform spacing of the underlying flow grid.

Figure 8Downwash at an elliptical wing predicted by AL simulations with different smearing factors and smearing correction. The CFD components of the velocities are shown (dashed) as well as the total downwash incorporating the correction (solid). The theoretical value acts as reference.

## 5.2 Rotor simulations – NREL 5-MW

The validity of the smearing hypothesis and its correction in rotor applications is demonstrated with simulations of the NREL 5-MW turbine using actuator line (AL) and lifting line (LL) models. The input parameters for these simulations are given in Table 2.

Table 2Input parameters for the NREL 5-MW simulations.

### 5.2.1 Uniform inflow

Figure 9 compares the AL results with and without the novel smearing correction to the LL with and without viscous core. At this wind speed of 8 m s−1 the thrust coefficient is highest (CT=0.84) – and hence induction is highest – thus lending itself as a strong verification case. Clearly, there is an equivalence between the original AL and the LL with a viscous core and the corrected AL with the LL. Therefore, the smearing correction makes the AL effectively act as a LL, as originally intended by . The impact of the viscous core is most prominent toward the blade root and tip. The sudden drop in the forces predicted by the AL/LL+core for the tip section of the blade – located at $r/R=\mathrm{0.97}$ – is not triggered by any aerodynamic tip effects, but relates to a pronounced reduction in chord. Figure 10 shows the diminishing effect of the correction on the angle of attack towards the root and tip. As the correction velocities are negligible with respect to the rotational velocity – they impact the velocity magnitude by less than 0.1 % in the lifting region of the blade – it is ultimately the change in the angle of attack that explains the observations in the force distributions. While not greatly affecting the magnitude of the forces in the mid-section of the blade, the viscous core does introduce greater fluctuations in the force distribution. Hence, the missing induction introduced by the viscous core reduces the coupling between neighbouring blade sections. The smearing correction also recovers this behaviour of the LL. Surpassing rated wind speed, forces increase inboard until cut-out. Thus, just before cut-out at 25 ms−1 loading reaches a maximum towards the root, causing an equally pronounced influence of the smearing correction in this region, as demonstrated in Fig. 11. Again, the equivalence of the AL and LL implementations is remarkable. This high wind speed case also demonstrates our correction is not only a tip correction.

Figure 9Normal and tangential forces on the NREL 5-MW blades at 8 m s−1 predicted by AL simulations with/without smearing correction and LL with/without viscous core (Ns=19, ϵ=0.1R).

Figure 10Angle of attack with/without smearing correction on the NREL 5-MW blades at 8 m s−1 (Ns=19, ϵ=0.1R).

The comparison of AL and LL is summarized in Fig. 12 in the form of local thrust and power distributions at different wind speeds. Note that for the wind speeds below rated (<11.4m s−1), the coefficients are identical. The results are only presented for simulations with Ns=9 for visibility, but compare equally well at higher resolution. As mentioned earlier, the smearing correction predominantly acts towards the tip and root. An additional overview of all results is given in Table 3. Here the total rotor thrust T and power P predicted by the corrected actuator line (AL*) and the lifting line (LL) are listed as well as the influence of adding the viscous core relative to AL* and LL, respectively. The AL and LL solutions are not directly compared to avoid including any mean bias in the comparison. The impact of the correction on the AL forces is nearly identical to removing the viscous core in the LL simulations at any wind speed, which further supports our correction methodology. In light of the large errors incurred without any correction, unsurprisingly, some form of tip correction is usually applied in AL simulations.

Figure 11Normal and tangential forces on the NREL 5-MW blades at 25 m s−1 predicted by AL simulations with/without smearing correction and LL with/without viscous core (Ns=9, ϵ=0.2R).

Figure 12Local thrust and power coefficients along the NREL 5-MW blades at different wind speeds predicted by AL simulations with/without smearing correction and LL with/without core (Ns=9, ϵ=0.2R).

Table 3An overview of the simulation inputs and results for the NREL 5-MW in uniform inflow. Results are grouped by blade/grid resolution. For the actuator line (AL) and lifting line (LL) the simulation time step Δt, the total thrust T and power P as well as the relative change in these quantities caused by the correction/removing the viscous core are listed. Note that AL* represents the corrected AL results, and the change is expressed relative to the rotor thrust and power.

### 5.2.2 Yawed inflow

As the smearing correction does not include yaw effects – the wake is assumed to advect normal to the rotor plane – we tested its influence at yaw angles χ of 15,30 and 45 at 8 m s−1. Again the LL with and without viscous core acted as a reference. The time steps remained the same as in uniform inflow. Here only the results for the most extreme case at 45 yaw are shown, as the differences are most severe in this case. Figure 13 presents the normal and tangential force variation during one rotation, averaged over three distinct regions of the blade, at a wind speed of 8 m s−1. Whilst the agreement is best towards the blade tip, the force variation with azimuthal position is similar between AL and LL simulations across all sections. AL results are shifted downwards with respect to the LL predictions at the inner sections, hinting at the AL experiencing higher induction in this region. However, for the verification of the smearing correction this shift is irrelevant, instead its impact on the AL forces needs to be assessed relative to the difference between LL with and without core. In this respect the smearing correction behaves correctly, increasing forces in a similar fashion as a LL without core in the mid-section of the blade and reducing them towards the root and tip.

Figure 13Normal (a, b, c) and tangential (d, e, f) force variation during one rotation as a function of azimuthal position of the NREL 5-MW blades at 8 m s−1 and 45 yaw. The forces are averaged over three sections (inner: a, d; middle: b, e; outer: c, f) of the blade and are predicted by AL simulations with/without smearing correction and LL with/without viscous core. The blade is facing upwind at χ=0 and is pointing vertically up at χ=90 (Ns=9, ϵ=0.2R).

### 5.2.3 Pitch step

The pitch step is defined as

$\begin{array}{}\text{(29)}& \mathit{\psi }={\mathit{\psi }}_{\mathrm{0}}+\frac{\mathrm{\Delta }\mathit{\psi }}{\mathrm{2}}\left[\mathrm{1}+\mathrm{tanh}\left(k\left(t-{t}_{\mathrm{0}}\right)\right)\right],\end{array}$

with ψ0 defining the pitch angle before the step, t0 representing the time instant of the step and Δψ denoting the pitch change. Here an extremely violent step is chosen – determined by k – to encourage an equally pronounced blade force response and test the numerical stability of our correction. The parameters governing this comparison are given in Table 4, which realize a pitch step of ±2 in 0.14 s (10 % to 90 % pitch). To capture the swift change in pitch, the time step is adjusted in both AL and LL simulations to $\mathrm{3.94}×{\mathrm{10}}^{-\mathrm{2}}$s. The blade force response is normalized as

$\begin{array}{}\text{(30)}& \stackrel{\mathrm{^}}{F}\left(t\right)=\frac{F\left(t\right)-{F}_{\mathrm{0}}}{{F}_{\mathrm{\infty }}-{F}_{\mathrm{0}}},\end{array}$

with F0 and F denoting the steady-state values before and after the pitch step, respectively.

Here only the tangential force response after a +2 step for the mid and tip blade sections are shown in Fig. 14, as they capture the main features of the response. As the definition here is positive pitch to feather, the force decreases along the blade for positive pitch changes. The AL simulations exhibit a faster response with a kink at 0.14 s, coinciding with the pitch change reaching 99 % of the step. Therefore, the AL seems to capture the pitch rate lift. The LL does not show this feature so – as in yaw – the correct behaviour of the smearing correction on the AL force response should be assessed relative to the influence of removing the viscous core in the LL model. Overall, the correction has limited effect on the dynamic response, which is also confirmed by the LL simulations, but the correction essentially acts on the forces in the same fashion as removing the core in the LL. In the mid-section it reduces the forces by a maximum of 1 % during the first 2s, dropping to 0.5 % afterwards. At the tip section it intensifies the response by a maximum of 1 %, diminishing to 0.3 % at 4s.

Table 4Inputs defining the pitch step.

Figure 14Normalized tangential force response at two blade sections (middle and tip) of the NREL 5-MW following a pitch step of +2 in 0.14 s at 14 m s−1 predicted by AL simulations with/without smearing correction and LL with/without viscous core (Ns=9, ϵ=0.2R).

### 5.2.4 Turbulent inflow

Highly turbulent inflow should challenge the numerical stability of the new smearing correction by introducing strong and abrupt changes in the angle of attack. Comparing simulations with and without inflow turbulence should also reveal whether turbulence alters the nature of the correction. Figure 15 shows the impact of the smearing correction on the time-averaged normal and tangential blade forces at an 8 m s−1 mean wind speed for AL simulations with uniform and turbulent inflow. At a turbulence intensity (TI) of 15 %, the forces are unsurprisingly slightly larger (≈20N m−1) than in uniform inflow. However, the change in forces introduced by the correction is nearly identical (<2N m−1). When comparing the standard deviation of the forces with and without correction in Fig. 16 the smoothing and dampening effect of the smearing correction on the forces in highly turbulent inflow also becomes apparent. observed a corresponding reduction of the force variations on the whole rotor blade when comparing near-wake model results against BEM results for the NM 80 rotor in turbulent inflow. This illustrates that the smearing correction leads to the same dynamic coupling between neighbouring blade sections as a lifting line model.

Figure 15Time-averaged normal and tangential forces on the NREL 5-MW blades at an 8 m s−1 mean wind speed and changing inflow turbulence predicted by the AL model with and without smearing correction (Ns=35, $\mathit{ϵ}=R/\mathrm{8}$).

Figure 16Variation of normal and tangential forces on the NREL 5-MW blades at an 8 m s−1 mean wind speed and turbulence intensity of 15 % predicted by the AL model with and without smearing correction (Ns=35, $\mathit{ϵ}=R/\mathrm{8}$).

6 Conclusions

The actuator line was intended as a lifting line technique for CFD applications. In this paper we prove – theoretically and practically – that smearing the forces of the actuator line in the flow domain leads to smeared velocity fields. For the typical Gaussian force projection, the widely known Lamb–Oseen (Lamb1932; Oseen1911) viscous core appears in both bound and trailed vorticity. This core reduces the velocities approaching the vortex centre compared with the inviscid solution of the lifting line. Thus, the trailed vorticity of an actuator line induces lower velocities at the blade owing to the force projection. We recover this missing induction by combining a near-wake model of the trailed vorticity with the Lamb–Oseen viscous core model and coupling it with the actuator line model. Basic wing test cases with theoretical solutions verify the correction, as it recovers nearly all induction independent of the severity of the force smearing. Furthermore, rotor simulations show the applicability and strength of the correction over the entire operational wind speed range as well as in yaw, strong turbulence and undergoing pitch steps. Here the correction is validated with lifting line simulations with and without viscous core, which are representative of an actuator line with and without smearing correction, respectively. The agreement between the respective actuator line and lifting line results is remarkable.

The current implementation of the smearing correction relies on heavy bookkeeping. In future versions the latter will be removed without jeopardizing stability or accuracy, making it suitable for wind farm simulations in realistic atmospheric flows. Potentially, the correction might also enable accurate rotor simulations at lower discretization.

Code and data availability
Code and data availability.

All data and parts of the code covering the smearing correction are available upon request. Commercial and research licenses for EllipSys3D can be purchased from DTU.

Appendix A: Fixed-wake equations

The equations governing the fixed-wake approach underlying the smearing correction (see Fig. 4) are summarized for completeness.

Direction vector from the vortex element to a control point along the blade:

$\begin{array}{}\text{(A1)}& \mathbit{x}=\left(\begin{array}{c}-r\mathrm{cos}\mathit{\beta }+r-h\\ r\mathrm{sin}\mathit{\beta }\\ -{v}_{h}\mathit{\beta }/\mathrm{\Omega }\end{array}\right).\end{array}$

Definition of the vortex element:

$\begin{array}{}\text{(A2)}& \mathbit{\delta }\mathbit{l}=\mathit{\delta }l\mathrm{cos}\mathit{\varphi }\left(\begin{array}{c}-\mathrm{sin}\mathit{\beta }\\ -\mathrm{cos}\mathit{\beta }\\ {v}_{h}/\left(\mathrm{\Omega }r\right)\end{array}\right)\end{array}$

with $\mathit{\delta }l=\frac{r\mathit{\delta }\mathit{\beta }}{\mathrm{cos}\mathit{\varphi }}$. Note that

$\begin{array}{}\text{(A3)}& \mathrm{cos}\mathit{\varphi }=\frac{\mathrm{\Omega }r}{\sqrt{{\left(\mathrm{\Omega }r\right)}^{\mathrm{2}}+{v}_{h}^{\mathrm{2}}}}=\frac{\mathrm{1}}{\sqrt{\mathrm{1}+{\left(\frac{{v}_{h}}{\mathrm{\Omega }r}\right)}^{\mathrm{2}}}}\end{array}$

with $\mathrm{tan}\mathit{\varphi }=\frac{{v}_{h}}{\mathrm{\Omega }r}$. Incremental velocity induced by the vortex element

$\begin{array}{ll}& \mathbit{\delta }\stackrel{\mathrm{̃}}{\mathbit{u}}=\frac{\mathrm{\Delta }\mathrm{\Gamma }}{\mathrm{4}\mathit{\pi }}\frac{\mathbit{\delta }\mathbit{l}×\mathbit{x}}{|\mathbit{x}{|}^{\mathrm{3}}}=\\ \text{(A4)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}Ar\left(\begin{array}{c}\left({v}_{h}/\mathrm{\Omega }r\right)\left(\mathit{\beta }\mathrm{cos}\mathit{\beta }-\mathrm{sin}\mathit{\beta }\right)\\ -\left({v}_{h}/\mathrm{\Omega }r\right)\left(-\mathrm{1}+{h}_{r}+\mathrm{cos}\mathit{\beta }+\mathit{\beta }\mathrm{sin}\mathit{\beta }\right)\\ \left[-\mathrm{1}+\left(\mathrm{1}-{h}_{r}\right)\mathrm{cos}\mathit{\beta }\right]\end{array}\right),\end{array}$

with

$\begin{array}{ll}& A=\frac{\mathrm{\Delta }\mathrm{\Gamma }r\mathit{\delta }\mathit{\beta }}{\mathrm{4}\mathit{\pi }}\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{\underset{|\mathbit{x}{|}^{\mathrm{2}}}{\underbrace{\left[{r}^{\mathrm{2}}\left(\mathrm{1}+{\left(\mathrm{1}-{h}_{r}\right)}^{\mathrm{2}}-\mathrm{2}\left(\mathrm{1}-{h}_{r}\right)\mathrm{cos}\mathit{\beta }+{\left(\frac{{v}_{h}\mathit{\beta }}{\mathrm{\Omega }r}\right)}^{\mathrm{2}}\right)\right]}}}^{-\frac{\mathrm{3}}{\mathrm{2}}}\\ \text{(A5)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{h}_{r}=\frac{h}{r}.\end{array}$
Author contributions
Author contributions.

ARMF was responsible for writing the paper, except for Sect. 3.2 which was authored by NRG. The original numerical implementation of the near-wake model by GRP was adapted by ARMF to incorporate the smearing correction with the help of GRP. NRG performed all simulations with the lifting line, and ARMF performed all simulations with the actuator line. All authors decided on the flow cases to simulate and commented on the paper.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

We would like to acknowledge DTU Wind Energy's internal project “Virtual Atmosphere” for partially funding this research. Furthermore many thanks to senior researcher Mac Gaunaa for his insights on vortex aerodynamics and senior researcher Niels Troldborg for his input on actuator line simulations/modelling (both from DTU Wind Energy). Thanks also to Ang Li for his help regarding the near-wake model.

Review statement
Review statement.

This paper was edited by Alessandro Bianchini and reviewed by two anonymous referees.

References

Ananthan, S. and Leishman, J. G.: Role of Filament Strain in the Free-Vortex Modeling of Rotor Wakes, J. Am. Helicopter Soc., 9, 176–191, 2004. a

Dag, K., Sørensen, J., Sørensen, N., and Shen, W.: Combined pseudo-spectral/actuator line model for wind turbine applications, PhD thesis, DTU Wind Energy, Denmark, 2017. a, b, c, d, e, f

Forsythe, J. R., Lynch, E., Polsky, S., and Spalart, P.: Coupled Flight Simulator and CFD Calculations of Ship Airwake using Kestrel, in: 53rd AIAA Aerospace Sciences Meeting, AIAA, Kissimmee, Florida, USA, https://doi.org/10.2514/6.2015-0556, 2015. a

Glauert, H.: Airplane Propellers, Springer Berlin Heidelberg, Berlin, Heidelberg, 169–360, https://doi.org/10.1007/978-3-642-91487-4_3, 1935. a

Jha, P. K., Churchfield, M. J., Moriarty, P. J., and Schmitz, S.: Guidelines for Volume Force Distributions Within Actuator Line Modeling of Wind Turbines on Large-Eddy Simulation-Type Grids, J. Sol. Energ.-T. ASME, 136, 031003, https://doi.org/10.1115/1.4026252, 2014. a

Jonkman, J., Butterfield, S., Musial, W., and Scott, G.: Definition of a 5-MW reference wind turbine for offshore system development, Tech. rep., NREL/TP-500-38060, National Renewable Energy Laboratory (NREL), Colorado, USA, 2009. a

Lamb, H.: Hydrodynamics, C.U.P, 6th Edn., Cambridge University Press, Cambridge, 1932. a, b, c

Leishman, J. G., Bhagwat, M. J., and Bagai, A.: Free-Vortex Filament Methods for the Analysis of Helicopter Rotor Wakes, J. Aircraft, 39, 759–775, 2002. a

Leonard, B.: A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comput. Method. Appl. M., 19, 59–98, 1979. a

Madsen, H. A., Sørensen, N. N., Bak, C., Troldborg, N., and Pirrung, G.: Measured aerodynamic forces on a full scale 2 MW turbine in comparison with EllipSys3D and HAWC2 simulations, J. Phys. Conf. Ser., 1037, 022011, https://doi.org/10.1088/1742-6596/1037/2/022011, 2018. a

Martínez-Tossas, L. A. and Meneveau, C.: Filtered lifting line theory and application to the actuator line model, J. Fluid Mech., 863, 269–292, https://doi.org/10.1017/jfm.2018.994, 2019. a, b

Menter, F. R.: Zonal two equation kω turbulence models for aerodynamic flows, in: 23rd Fluid Dynamics, Plasmadynamics, and Lasers Conference, Fluid Dynamics and Co-located Conferences, Orlando,FL, https://doi.org/10.2514/6.1993-2906, 1993. a

Meyer Forsting, A., Troldborg, N., Bechmann, A., and Réthoré, P.-E.: Modelling Wind Turbine Inflow: The Induction Zone, PhD thesis, DTU Wind Energy, Denmark, https://doi.org/10.11581/DTU:00000022, 2017. a, b

Michelsen, J.: Basis3D – a platform for development of multiblock PDE solvers, Tech. rep., Dept. of Fluid Mechanics, Technical University of Denmark, DTU, 1994a. a

Michelsen, J.: Block structured multigrid solution of 2D and 3D elliptic PDE's, Tech. rep., Dept. of Fluid Mechanics, Technical University of Denmark, DTU, 1994b. a

Mikkelsen, R.: Actuator Disc Methods Applied to Wind Turbines, PhD thesis, Technical University of Denmark, 2003. a, b

Nathan, J.: Application of Actuator Surface Concept in LES Simulations of the Near Wake of Wind Turbines, PhD thesis, École de Technologie Supérieure, Montreal, Canada, 2018. a

Oseen, C.: Über Wirbelbewegung in einer reibenden Flüssigkeit, Arkiv för matematik, astronomi och fysik, Ark. Mat. Astron. Fys., 7, 14–21, 1911. a, b, c

Patanker, S. and Spalding, D.: A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, Int. J. Heat Mass Tran., 15, 59–98, 1972. a

Pirrung, G. R., Hansen, M. H., and Madsen, H. A.: Improvement of a near wake model for trailing vorticity, J. Phys. Conf. Ser., 555, 012083, https://doi.org/10.1088/1742-6596/555/1/012083, 2014. a, b

Pirrung, G., Madsen, H. A., Kim, T., and Heinz, J.: A coupled near and far wake model for wind turbine aerodynamics, Wind Energy, 19, 2053–2069, https://doi.org/10.1002/we.1969, 2016. a, b

Pirrung, G., Riziotis, V., Madsen, H., Hansen, M., and Kim, T.: Comparison of a coupled near- and far-wake model with a free-wake vortex code, Wind Energ. Sci., 2, 15–33, https://doi.org/10.5194/wes-2-15-2017, 2017a. a

Pirrung, G. R., Madsen, H. A., and Schreck, S.: Trailed vorticity modeling for aeroelastic wind turbine simulations in standstill, Wind Energ. Sci., 2, 521–532, https://doi.org/10.5194/wes-2-521-2017, 2017b. a, b

Ramos-García, N., Shen, W. Z., and Sørensen, J. N.: Validation of a three-dimensional viscous-inviscid interactive solver for wind turbine rotors, Renew. Energ., 70, 78–92, https://doi.org/10.1016/j.renene.2014.04.001, 2014a. a

Ramos-García, N., Shen, W. Z., and Sørensen, J. N.: Three-dimensional viscous-inviscid coupling method for wind turbine computations, Wind Energy, 19, 67–93, https://doi.org/10.1002/we.1821, 2014b. a

Ramos-García, N., M. Mølholm, J. N. S., and Walther, J. H.: Hybrid vortex simulations of wind turbines using a three-dimensional viscous-inviscid panel method, Wind Energy, 20, 1187–1889, 2017. a, b

Réthoré, P.-E. and Sørensen, N.: A discrete force allocation algorithm for modelling wind turbines in computational fluid dynamics, Wind Energy, 15, 915–926, https://doi.org/10.1002/we.525, 2012. a

Shen, W., Sørensen, J., and Mikkelsen, R.: Tip loss correction for actuator/Navier Stokes computations, J. Sol. Energ.-T. ASME, 59, 209–213, 2005. a

Shives, M. and Crawford, C.: Mesh and load distribution requirements for actuator line CFD simulations, Wind Energy, 16, 657–669, https://doi.org/10.1002/we.1546, 2013. a, b, c, d, e

Sørensen, J. N. and Shen, W. Z.: Numerical modelling of wind turbine wakes, J. Fluid. Eng.-T. ASME, 124, 393–399, https://doi.org/10.1115/1.1471361, 2002.  a, b, c

Sørensen, N.: General purpose flow solver applied to flow over hills, PhD thesis, Risø National Laboratory, 1995. a, b

Squire, H. B.: The growth of a vortex in turbulent flow, Aeronaut. Quart., 16, 302–306, 1965. a

Strelets, M.: Detached eddy simulation of massively separated flows, in: 39th AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper 2001-0879, Reno, NV, 2001. a

Troldborg, N., Sørensen, J., and Mikkelsen, R.: Actuator Line Modeling of Wind Turbine Wakes, PhD thesis, Technical University of Denmark, 2009. a, b, c, d

Vollmer, L., Steinfeld, G., and Kühn, M.: Transient LES of an offshore wind turbine, Wind Energ. Sci., 2, 603–614, https://doi.org/10.5194/wes-2-603-2017, 2017. a

$\mathrm{\nabla }×\left({\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}×\mathbit{u}\right)=\left(\mathbit{u}\mathrm{\nabla }\right){\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}-\left({\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}\mathrm{\nabla }\right)\mathbit{u}+{\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}\left(\mathrm{\nabla }\mathbit{u}\right)-\mathbit{u}\left(\mathrm{\nabla }{\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}\right)=\left(\mathbit{u}\mathrm{\nabla }\right){\stackrel{\mathrm{^}}{\mathbit{e}}}_{x}+\mathrm{0}+\mathrm{0}+\mathrm{0}$

Note that in the xy plane the circulation would be −Γ.

Strictly, the influence of the shed vorticity on the velocity at the AL should be removed as remarked by , however its influence is negligible.