2.1. Lift force in the basic formulation of the actuator line method

In the general formulation of the actuator line (Sorensen & Shen Reference Sorensen and Shen2002), the lift force $F_l$ at each spanwise section $j$ of the wing (force per spanwise length) is obtained from aerofoil data from

(2.8)\begin{equation} F_{lj} = \tfrac{1}{2}\rho u_r^2 c_j\,C_l(\alpha_j), \end{equation}

where $c_j$ is the local chord of the aerofoil, and $C_l(\alpha _j)$ is the local lift coefficient, obtained from aerofoil data at an angle of attack $\alpha _j$. The local reference system of figure 1 is used here. The velocities $u_z$ and $u_y$, relative to the actuator line, are obtained from the computational fluid dynamics (CFD) simulation at the control point position $\boldsymbol {x}_j$ of the actuator line segment. The lift coefficient $C_l(\alpha _j)$ is calculated by the interpolating tabulated aerofoil coefficients using the effective local angle of attack

(2.9)\begin{equation} \alpha_j=\alpha_{gj} + \arctan{\left(\frac{u_y}{u_z}\right)}, \end{equation}

where $\alpha _{gj}$ is the geometric angle of attack given by the twist and incidence of the wing (or, in the case of rotating blades, the opposite of the local pitch angle; Troldborg Reference Troldborg2009).

If needed, the circulation at each spanwise section can be calculated from the Kutta–Joukowski theorem

(2.10)\begin{equation} \varGamma_j = \frac{F_{lj}}{\rho u_r} = \frac{1}{2}\,u_r c_j\,C_l(\alpha_j) . \end{equation}

Usually, some form of tip or smearing correction is applied to this basic formulation. The smearing corrections employed in the present work are described in §§ 2.3 and 5.

2.2. Vortex-based smearing correction

In the actuator line method with smearing correction, the velocity $\boldsymbol {u}^s$, sampled from the CFD simulation at control point $\boldsymbol {x}_j$, is summed to the ‘missing velocity’ $\boldsymbol {u}^m$ to arrive at the corrected velocity

(2.11)\begin{equation} \boldsymbol{u}^{c}(\boldsymbol{x}_j) = \boldsymbol{u}^{s}(\boldsymbol{x}_j) + \boldsymbol{u}^{m}(\boldsymbol{x}_j) . \end{equation}

The ‘missing velocity’ is defined as the difference between the velocities induced by the vortices created by the actuator line and ‘reference’ vortices. The concept of what is considered a ‘reference’ vortex varies between the works. Dağ & Sørensen (Reference Dağ and Sørensen2020) and Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a) consider a vortex filament with an infinitesimal core, while Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) consider the vortex created using the optimum smearing parameter developed in Martínez-Tossas et al. (Reference Martínez-Tossas, Churchfield and Meneveau2017). A viscous core model with the radius of the vortex core evolving in time could also be implemented easily. If there is the possibility of the vortices impinging the actuator lines, then it might be useful to employ some desingularization method.

In the present work, we adopt the definition of reference vortices as vortex filaments with an infinitesimal core. In this method, the missing velocity is calculated from the difference of velocities induced by ideal, singular, vortices (superscript ${vi}$) and vortices with finite core created by the actuator lines (superscript $v$):

(2.12)\begin{equation} \boldsymbol{u}^{m}(\boldsymbol{x}_j) = \boldsymbol{u}^{vi}(\boldsymbol{x}_j) - \boldsymbol{u}^{v}(\boldsymbol{x}_j) . \end{equation}

Using this definition of missing velocity, the aim of the method is to reproduce the results of a lifting line method. The idea behind the vortex-based smearing correction is that the velocity sampled from the numerical simulation, in a linear approximation, is given by the sum of the local undisturbed velocity $\boldsymbol {U}$ and the velocity induced by the vortices created by the actuator line:

(2.13)\begin{equation} \boldsymbol{u}^{s}(\boldsymbol{x}_j) \approx \boldsymbol{U}(\boldsymbol{x}_j) + \boldsymbol{u}^{v}(\boldsymbol{x}_j) . \end{equation}

From (2.11), the corrected velocity becomes

(2.14)\begin{equation} \boldsymbol{u}^{c}(\boldsymbol{x}_j) \approx \boldsymbol{U}(\boldsymbol{x}_j) + \boldsymbol{u}^{vi}(\boldsymbol{x}_j) , \end{equation}

which reproduces the results of a lifting line method. It should be noted that the local undisturbed velocity $\boldsymbol {U}$ is not known in ALM simulations (except for simple cases, used mostly as validation). Only the velocity sampled from the simulations $\boldsymbol {u}^{s}$ is known. From this comes the need to model the velocity induced by the vorticity created by the actuator lines, $\boldsymbol {u}^{v}$, described in § 3.

On top of the results from the vortex-based smearing correction, other corrections can be introduced, based on the known limitations of the lifting line method. For example, Dağ (Reference Dağ2017) combined the smearing correction with the decambering correction of Sørensen, Dag & Ramos-García (Reference Sørensen, Dag and Ramos-García2016). Limitations and further corrections of the lifting line method are out of the scope of the present work, and we consider it as the reference for validation of the vortex-based smearing correction.

2.3. Iterative smearing correction

The missing velocity is calculated from the circulation of the bound and wake vortices. At the same time, the circulation is calculated from the corrected velocity, using (2.10). For this reason, an iterative procedure was used by Dağ & Sørensen (Reference Dağ and Sørensen2020) and Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a). The approach of Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a) applies the relaxation factor to the velocity. In this work, we apply the relaxation factor to the circulation, because we solve for the circulation in § 5.

The steps of the iterative smearing correction, for each time step, are as follows.

1. (i) Start from a guess of the circulation distribution, usually the circulation from the previous time step.

2. (ii) Form the vortex sheet by:

   1. (a) prescribing the vortex sheet, for example, by assuming helical or horseshoe vortices; or

   2. (b) employing a free-vortex wake method, for example, by advecting the vortices with the CFD velocity or by a combination of the CFD velocities and the velocities induced by the free vortices.

3. (iii) Calculate the missing velocity $\boldsymbol {u}^m(\boldsymbol {x}_j)$ at every control point. Find the local corrected velocities according to (2.11).

4. (iv) Calculate the effective angle of attack using (2.9).

5. (v) Find the local lift coefficient $C_l$, interpolating from the aerofoil data table.

6. (vi) Calculate the new value of circulation $\varGamma _j^{new}$ using (2.10).

7. (vii) Update the current value of circulation using a relaxation factor $r$:

   (2.15)\begin{equation} \varGamma_j=r\varGamma_j^{new} + (1-r) \varGamma_j^{old}. \end{equation}

8. (viii) Restart from step (ii) if the local velocity affects the formation of the vortex sheet, or from step (iii) otherwise, using the value of circulation calculated in the previous step. Iterate until a chosen convergence criterion is reached.

The forces are calculated after convergence. There is an ambiguity regarding the choice of velocity used to calculate the forces. This ambiguity is resolved for the lift force in Kleine, Hanifi & Henningson (Reference Kleine, Hanifi and Henningson2023). For the present simulations, the difference is negligible, because the error is of second order (see Kleine et al. (Reference Kleine, Hanifi and Henningson2023) for further discussion on this topic). In the present work, all forces are calculated using the corrected velocities.

3.1. Vorticity generated by body forces

Some of the development of Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) is reproduced in this section for completeness. Neglecting viscous effects in a steady flow, the vorticity equation in steady flow becomes

(3.1)\begin{equation} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\omega} = \boldsymbol{\omega} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \times \boldsymbol{f} . \end{equation}

Equation (3.1) can be linearized considering a uniform base flow $\boldsymbol {U} = const.$, arriving at

(3.2)\begin{equation} \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\omega} = \boldsymbol{\nabla} \times \boldsymbol{f} . \end{equation}

Considering uniform flow aligned with the $z$-axis, $\boldsymbol {U} = (0,0,U_z)$, for a straight wing with only lift forces, that can be considered to act aligned with the $y$-axis, $\boldsymbol {f}=(0,f_y,0)$, we have

(3.3)\begin{gather} U_z\,\frac{\partial\omega_x}{\partial z} = - \frac{\partial f_y}{\partial z} \implies \omega_x= - \frac{f_y}{U_z}, \end{gather}

(3.4)\begin{gather} U_z\,\frac{\partial\omega_z}{\partial z} = \frac{\partial f_y}{\partial x} \implies \omega_z=\int_{-\infty}^{z} \frac{1}{U_z}\,\frac{\partial f_y}{\partial x} \,{\rm d}z, \end{gather}

while $\omega _x,\omega _z \gg \omega _y \approx 0$. (For a more detailed derivation, the reader is referred to Kleine Reference Kleine2022.)

The bound vorticity $\omega _x^i$ generated by an ideal concentrated force $f_y^i = - F_l(x)/\rho \, \delta (y)\,\delta (z)$ is

(3.5)\begin{equation} \omega_x^i = - \frac{f_y^i}{U_z} = \frac{1}{U_z \rho}\,F_l(x)\,\delta(y)\,\delta(z) , \end{equation}

where it is relevant to note that a positive lift force $F_l$ applies a negative force to the fluid. The integral of the vorticity in the spanwise direction is equal to the circulation $\varGamma$ at the section. For the case with concentrated force,

(3.6)\begin{equation} \varGamma (x) = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \omega_x^i \,{{\rm d} y} \,{\rm d}z = \frac{1}{U_z \rho}\,F_l(x)\,\delta(y)\,\delta(z) = \frac{1}{U_z \rho}\,F_l(x) , \end{equation}

which agrees with the Kutta–Joukowski theorem.

3.2. Vorticity generated by a discontinuous distribution of force

In numerical methods, usually the force is discretized and represented by interpolating functions inside segments. The theory for the semi-infinite wing of Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) can be used to construct a discretized method. In the current implementation of the actuator line method, for each segment $j$, the force is calculated at control point $x_{j}$ and considered constant for $x$-positions between the boundaries of the segment, $x_{j_-} \le x \le x_{j_+}$. The two-dimensional force can be defined based on

(3.7)\begin{equation} \boldsymbol{F}_{2D}(x)=\sum_{j=1}^{N}(H(x-x_{j_-})-H(x-x_{j_+}))\,\boldsymbol{F}_{2D}(x_{j}) , \end{equation}

where $H(z)$ is the Heaviside step function with the half-maximum convention. From (2.3), the ideal body force is then given by

(3.8)\begin{equation} \boldsymbol{f}^i(x,y,z) = \sum_{j=1}^{N}(H(x-x_{j_-})-H(x-x_{j_+}))\,\delta(y)\,\delta(z)\, \frac{\boldsymbol{F}_{2D}(x_{j})}{\rho} , \end{equation}

which, in most cases, will be discontinuous at the boundary of the segments (the average is considered at the boundary of the segments). The convolution with a Gaussian function (2.7) transforms the body forces to

(3.9)\begin{equation} \boldsymbol{f}(x,y,z) = \sum_{j=1}^{N}(H_{\varepsilon}(x-x_{j_-})-H_{\varepsilon}(x-x_{j_+}))\, \eta(y)\,\eta(z)\,\frac{\boldsymbol{F}_{2D}(x_{j})}{\rho} , \end{equation}

where the function

(3.10)\begin{equation} H_{\varepsilon}(x)=\frac{\operatorname{erf}\left(\dfrac{x}{\varepsilon}\right)+1}{2} \end{equation}

is defined as a smeared (mollified) Heaviside step function (Caprace et al. Reference Caprace, Chatelain and Winckelmans2019).

Equation (3.9) is employed directly in the current simulations. As far as the authors are aware, this may be the first implementation of the analytical convolution in an ALM code. However, it should be noted that most works do not describe how the convolution operation is performed in practice. We are aware, however, that many implementations, including the previous version of our code (Kleusberg Reference Kleusberg2019), perform the convolution operation numerically.

Considering only lift force, the bound vortex is (from (3.3))

(3.11)\begin{gather} \omega_x^i= \sum_{j=1}^{N}(H(x-x_{j_-})-H(x-x_{j_+}))\,\delta(y)\,\delta(z)\,\varGamma_j, \end{gather}

(3.12)\begin{gather} \omega_x = \sum_{j=1}^{N} (H_{\varepsilon}(x-x_{j_-})-H_{\varepsilon}(x-x_{j_+}))\,\eta(y)\,\eta(z)\, \varGamma_j , \end{gather}

where $\varGamma _j=F_l(x_j)/(U_z \rho )$. We can notice that the vorticity not only spreads as a Gaussian function in the directions normal to the actuator line, but also spread outside the boundaries of the line in the spanwise direction (terms $H_{\varepsilon }$). From (3.4), the shed vorticity is

(3.13)\begin{gather} \omega_z^i= - \sum_{j=1}^{N} \left[ \delta(x-x_{j_-})\,\delta(y)\,H(z) - \delta(x-x_{j_+}) \delta(y) H(z) \right] \varGamma_j, \end{gather}

(3.14)\begin{gather} \omega_z= - \sum_{j=1}^{N} \left[ \eta(x-x_{j_-})\,\eta(y)\,H_{\varepsilon}(z) - \eta(x-x_{j_+})\, \eta(y)\,H_{\varepsilon}(z) \right] \varGamma_j . \end{gather}

The first two terms of (3.13) are the singular semi-infinite vortices generated by the discontinuities in the circulation at the boundary of the segment, in positions $x_{j_-}$ and $x_{j_+}$. These terms in (3.14) become semi-infinite Lamb–Oseen vortices centred in $x_{j_-}$ and $x_{j_+}$. This case regresses to a discretized Prandtl lifting line in the case of (3.13), and to a lifting line with vortices with Gaussian core in the case of (3.14)

The vorticity of (3.11) and (3.13) is identical to the vorticity considered in the lifting line method, under the assumption that the velocity $U_z$ is approximately constant inside each segment. For rotating blades, this assumption is an approximation, since the local velocity changes along the blades. Nevertheless, this assumption is consistent with the linear theory and is a good approximation if the length of the segment is small compared to the radius.

The missing velocity is defined as the velocity needed to recover the velocity induced by singular vortices, as explained in § 2.2. If the implementation of the ALM considers constant circulation inside each segment, then the correction proposed by Dağ & Sørensen (Reference Dağ and Sørensen2020) is equivalent to a lifting line method, according to the linear approximation, as proposed by Dağ & Sørensen (Reference Dağ and Sørensen2020) and Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a) (putting aside other simplifications, as discussed in §§ 3.3 and 3.4). If the implementation of the ALM considers a non-constant distribution of circulation inside each segment, then the correction by Dağ & Sørensen (Reference Dağ and Sørensen2020) is not identical to the discrete lifting line method, because (3.13) and (3.14) would have extra terms. This case is not treated here, but could be a topic for further research. Interestingly, also the discretized implementation of Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) is not exactly identical (although very similar) to the classical lifting line method, because it considers implicitly the vortices located at the control points (equation (5.7) of that reference), not at the boundaries of segments, as would be usual in a classical discretized lifting line.

The bound vortex is also affected by the convolution with the Gaussian kernel, as seen in (3.12). However, it does not affect the forces in most cases, since most geometries have straight wings or blades, and a straight bound vortex does not induce velocity on itself. For multi-blade configurations, the blades are usually too distant for the correction to make any difference (distance is several times $\varepsilon$). A possible exception may be the hub; however, the circulation at the hub is low, and the forces are not as relevant for performance. As a consequence, the effect of the convolution on the vorticity generated and the implementation of corrections for bound vorticity have been, justifiably, neglected.

3.3. Velocity induced by a smeared vortex segment

In the method of Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a) and Dağ & Sørensen (Reference Dağ and Sørensen2020), the vortex wake is discretized by segments, and the missing velocity is calculated for each finite-length segment of vortex with a Gaussian distribution of vorticity. Hence it is relevant to study the velocity induced by a finite-length segment with a cross-section of Gaussian vorticity distribution.

Without loss of generality, we consider a straight segment of vortex filament with constant circulation $\varGamma _j$, located at $(x,y)=(0,0)$ and aligned with the $z$-direction, from $z_{j_-}$ to $z_{j_+}$, given by

(3.15)\begin{equation} \omega_z^i = \delta(x)\,\delta(y) (H(z-z_{j_-})-H(z-z_{j_+}) ) \varGamma_j . \end{equation}

A ‘smeared vortex segment’ is defined as the smeared equivalent of the straight segment of vortex filament formed by the convolution of (3.15) with a Gaussian function in three dimensions. In order to compute the velocity induced by a smeared vortex segment, the approach of Leonard (Reference Leonard1980) is taken (see also Leonard Reference Leonard1985; Winckelmans Reference Winckelmans1989). The vorticity of a segment of vortex filament convoluted with a three-dimensional Gaussian function is given by

(3.16)\begin{equation} \omega_z = \varGamma_j \int^{z_{j_+}}_{z_{j_-}} \eta_3(x,y,z-z') \,{\rm d}z' = \eta(x)\,\eta(y) ( H_{\varepsilon}(z-z_{j_-})-H_{\varepsilon}(z-z_{j_+}) ) \varGamma_j . \end{equation}

Figure 2 shows a vortex filament of length $10 \varepsilon$ with concentrated (infinite) vorticity in black and the vorticity distribution of the corresponding smeared vortex segment.

Figure 2. Smeared vortex segment of length $z_{j_+}-z_{j_-} = 10 \varepsilon$. Isocontours and volume rendering indicate normalized vorticity $\omega _z {\rm \pi}\varepsilon ^2/\varGamma _j$, given by (3.16). The black line indicates the segment of vortex filament with concentrated vorticity.

The velocity induced by a smeared vortex segment is given by Leonard (Reference Leonard1980) as

(3.17)\begin{equation} \boldsymbol{u}^{v}(\boldsymbol{x}) = \frac{\varGamma_j}{4 {\rm \pi}} \int^{z_{j_+}}_{z_{j_-}} \frac{\hat{\boldsymbol{e}}_x \times (\boldsymbol{x}-\boldsymbol{x'})}{|\boldsymbol{x}-\boldsymbol{x'}|^3}\,g( |\boldsymbol{x}-\boldsymbol{x'}| ) {\rm d}z' , \end{equation}

where function $g(s)$ is defined as

(3.18)\begin{equation} g(s) = 4 {\rm \pi}\int^s_0 \eta_3(s')\,s'^2 \,{\rm d}s' = \operatorname{erf}{\left(\frac{s}{\varepsilon}\right)} - 2 s\, \eta(s), \end{equation}

in which the three-dimensional Gaussian kernel $\eta _3(s')$ (see (2.5)) is rewritten in terms of the spherical coordinate $s'=\sqrt {x^2+y^2+z^2}$. Hence the induced velocity, in the azimuthal direction, is

(3.19)\begin{align} u_{\theta}^{v}(r,z) &= \frac{\varGamma_j}{4 {\rm \pi}}\,r \int^{z_{j_+}}_{z_{j_-}} \frac{\operatorname{erf}{\left(\dfrac{\sqrt{r^2+(z-z')^2}}{\varepsilon}\right)} - 2 \sqrt{r^2+(z-z')^2}\, \eta\left(\sqrt{r^2+(z-z')^2}\right)}{(r^2+(z-z')^2)^{3/2}} \,{\rm d} z'\nonumber\\ &= \frac{\varGamma_j}{4 {\rm \pi}} \left( \varPhi(r,z-z_{j_+}) - \varPhi(r,z-z_{j_-}) \right) , \end{align}

where

(3.20) \begin{equation} \varPhi(r,Z) = \frac{1}{r} \left( - \frac{Z}{\sqrt{r^2+Z^2}} \operatorname{erf}{\Bigg(\frac{\sqrt{r^2+Z^2}}{\varepsilon}\Bigg)} + \exp{\left(-\frac{r^2}{\varepsilon^2}\right)} \operatorname{erf}{\left(\frac{Z}{\varepsilon}\right)} \right) . \end{equation}

Analogously, $\varPhi ^{i}(r,Z)$ can be defined as

(3.21)\begin{equation} \varPhi^{i}(r,Z) := \lim_{\varepsilon \to 0} \varPhi(r,Z) = \frac{1}{r} \left( - \frac{Z}{\sqrt{r^2+Z^2}} \right) , \end{equation}

and the velocity induced by a singular vortex segment can be calculated as

(3.22)\begin{equation} u_{\theta}^{vi}(r,z) = \frac{\varGamma_j}{4 {\rm \pi}} (\varPhi^i(r,z-z_{j_+}) - \varPhi^i(r,z-z_{j_-})) , \end{equation}

which agrees with the classical results for a vortex segment (see e.g. Katz & Plotkin Reference Katz and Plotkin1991).

Therefore, in the present work, for each vortex segment, the missing velocity is then calculated as

(3.23) \begin{align} u_{\theta}^{m}(r,z) &= u_{\theta}^{vi}(r,z)-u_{\theta}^{v}(r,z) \nonumber\\ &= \frac{\varGamma_j}{4 {\rm \pi}} ((\varPhi^i(r,z-z_{j_+}) - \varPhi^i(r,z-z_{j_-})) - (\varPhi(r,z-z_{j_+}) - \varPhi(r,z-z_{j_-}))). \end{align}

The same correction is applicable to the bound vortices and wake vortices, taking into account the different orientations of the vortices.

For $z=0$, $z_{j_-}=0$ and $z_{j_+} \to +\infty$, the induced velocity is given by

(3.24)\begin{equation} u_{\theta}^{v}(r,0) = \frac{\varGamma_j}{4 {\rm \pi}r} \left( 1 - \exp{\left(-\frac{r^2}{\varepsilon^2}\right)} \right), \end{equation}

which is the formula for a semi-infinite vortex used by Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) (and consistent with the relationship mentioned by Meyer Forsting et al. Reference Meyer Forsting, Pirrung and Ramos-García2019a; Dağ & Sørensen Reference Dağ and Sørensen2020).

However, (3.23) is clearly different from the approximations employed in Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a) and Dağ & Sørensen (Reference Dağ and Sørensen2020). The approximation for $\varPhi$ used in these works for a vortex segment introduces errors. For the formula of Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a), negative and positive errors compensate each other for a semi-infinite vortex, but the same cannot be concluded for other vortex geometries. Such errors are not expected to be present in the application of the correction of Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) to straight wings in uniform flow. However, the application of the formula of Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) to the case of rotating blades without adaptations, such as performed in Stanly et al. (Reference Stanly, Martínez-Tossas, Frankel and Delorme2022), may incur other errors, since the helical vortex structure formed by the blades is being modelled by a vortex sheet formed by straight horseshoe vortices.

Nevertheless, the good results of Dağ & Sørensen (Reference Dağ and Sørensen2020) and Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a) show that their corrections improve the results when compared to ALM without correction or with other heuristic tip corrections, even using the approximate value of $\varPhi$. A possible explanation for that is based on the observation that the missing velocity is related directly to the difference $\varPhi ^i-\varPhi$. For the vortices near the blades that tend to dominate the correction, the value of $\varPhi ^i$ is much larger than $\varPhi$, so the exact value of $\varPhi$ is not as relevant, as long as the order of magnitude is the same. The dominance of the vortices near the blades is also a possible explanation for the improved results observed by Stanly et al. (Reference Stanly, Martínez-Tossas, Frankel and Delorme2022), even considering straight horseshoe vortices instead of helical vortices. The estimate of the order of magnitude of the errors caused by these different approximations is left for further studies. In the present work, the analytical integration, represented by (3.20), is used to avoid the approximation errors.

3.4. Forming the vortex sheet

Forming the vortex sheet is an intrinsic part of both the ALM with smearing correction and the lifting line method. It should be noted, however, that the methods are not going to be mathematically identical even if the vortex sheet is formed with the same procedure in both methods. In the actuator line, the induced velocity is formed partly by the smeared vortices created in the CFD simulation and partly by the vortex sheet of the smearing correction, which is formed by one of the strategies detailed below (that may be more or less dependent on the results from the simulation). Hence even if the vortex sheet from the smearing correction is identical to the vortex sheet of the lifting line method, the induced velocities can differ slightly, because the vortex sheet created in the CFD simulation may be different from the vortex sheet of the lifting line method. This phenomenon can be observed in the results of § 7.1, where its effect is noticeable but negligible.

The formation of the vortex sheet for the smearing correction varies in past works. The method of Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) does not require an explicit vortex sheet, since the correction is based on analytical semi-infinite vortices for a straight wing under uniform flow. However, the same correction was applied by Stanly et al. (Reference Stanly, Martínez-Tossas, Frankel and Delorme2022), which would mean that the helical vortex structure is modelled implicitly by horseshoe vortices. As discussed in § 3.3, errors are expected to be introduced by this approximation.

In Dağ & Sørensen (Reference Dağ and Sørensen2020), the helical vortices were imposed considering the pitch of each of the helical vortices based on the local relative flow angle at the position where the vortex is released. In Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a), a fixed helical wake was assumed, based on the near-wake model of Pirrung et al. (Reference Pirrung, Madsen, Kim and Heinz2016) and Pirrung, Madsen & Schreck (Reference Pirrung, Madsen and Schreck2017), and the bookkeeping of the position of released vortices is necessary for only one of the factors involved in the smearing correction. The need for bookkeeping was later removed by Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2020), reducing computational cost with negligible effect on the forces.

In the present work, we implement a free-vortex wake method within our simulation, a strategy already suggested by Dağ & Sørensen (Reference Dağ and Sørensen2020), assuming that the vortices are advected by the flow velocity. Tracing particles are used to follow the positions of the boundaries of the vortex segments. Also, the past values of circulation, which are present in the wake, are used directly, which is a more realistic situation than considering the whole wake as having the value of circulation that the wing currently assumes. This method has the drawback of requiring bookkeeping of the circulation and position of the vortices emitted in previous time steps. Nevertheless, it has the advantage of requiring no ad hoc assumption and therefore the method maintains its generality.

There is an ambiguity regarding using the corrected velocity, or the uncorrected velocity sampled directly from the numerical simulation, to estimate the position of vortices. However, this difference is of second order for the velocities along the actuator line. Using the corrected velocity might even introduce errors or numerical instability if used throughout the vortex sheet, because of the singularity of ideal vortices: vortices released recently might suffer the effect of the singularity of the bound vortex, and curved vortices might induce numerical instabilities.

Hence, in the current implementation, the tracing particles that define the boundaries of the vortex filaments are advected by the velocities sampled from the simulation, without correction. This reduces the number of operations needed so that the correction velocity is calculated on only the control points of the actuator line, not on other points of the vortex sheet. More relevant for our non-iterative method, the vortex sheet can be advected before the computation of the correction velocity at a certain time step. For its lower computational cost, a first-order Euler time integration method is used to advect the tracing particles. The information related to time steps older than $n-n_w$ is not stored, where $n$ is the current time step, and $n_w$ is the number of tracing particles.

In order to save computational resources, not all tracing particles are kept. The vortices released in the last $n_{nw}$ time steps are always maintained in the memory. However, a group of vortices is fused if the distance between the tracing particles is below $d_{w}$, with the circulation taken as the average of the circulation of the vortices fused together. The strategy of fusing vortices guarantees a minimum wake length of $(n_w-n_{nw}-1)d_{w}$ near the tip of the blades, the region of interest and where the advection velocity is larger. For the current simulations, the choice of parameters is $n_w=50$, $n_{nw}=10$ and $d_w=\varepsilon /2$, guaranteeing a wake length of approximately $20\varepsilon$ near the tip. Analysis of the correction of (3.23) indicates that errors involved in discarding vortices farther than this length are negligible. The wake length is lower near the hub, but this region is less relevant for the performance of the turbine. A parameter study for the simulation of the NREL 5-MW turbine of § 7.2 showed that doubling $n_w$ and $n_{nw}$ had negligible effects on the circulation.

For unsteady simulations, a vorticity wake in the spanwise direction would be shed when the circulation changes, similar to the starting vortex of Kelvin's circulation theorem. Following a quasi-steady approach, the correction for this spanwise shed vorticity is neglected, hence it is not modelled in our method. For further discussion about this shed vorticity and unsteady effects, the reader is referred to Kleine (Reference Kleine2022). The effect of the Gaussian smearing on unsteady effects is left for future studies.

6.1. Numerical solver of Navier–Stokes equations with the actuator line

The incompressible three-dimensional Navier–Stokes equations are solved in the spectral-element code Nek5000 (Fischer, Lottes & Kerkemeier Reference Fischer, Lottes and Kerkemeier2008). The system of equations is expressed in weak form, and the solution is expanded in terms of basis and test functions on each of these elements (Offermans et al. Reference Offermans, Peplinski, Marin and Schlatter2020). In this work, the basis for the velocity space consists of seventh-order Lagrange polynomial interpolants on Gauss–Lobatto–Legendre quadrature points in each element. The $\mathbb {P}_N-\mathbb {P}_{N-2}$ formulation is used (Maday & Patera Reference Maday and Patera1989). For temporal discretization, a third-order implicit/explicit scheme (BDF3/EXT3) (Fischer Reference Fischer2003) is employed. Filtering of the higher modes is applied to stabilize the simulation (Fischer & Mullen Reference Fischer and Mullen2001).

In the current simulations, an adaptive mesh refinement strategy with a spectral error indicator (Offermans Reference Offermans2019; Offermans et al. Reference Offermans, Peplinski, Marin and Schlatter2020; Tanarro et al. Reference Tanarro, Mallor, Offermans, Peplinski, Vinuesa and Schlatter2020) is used to reduce the computational cost. We choose to force maximum discretization of the grid around the actuator lines, independently of the error indicator, to guarantee adequate discretization for the chosen value of the smearing parameter.

The actuator line was implemented previously in Nek5000 (Kleusberg et al. Reference Kleusberg, Sarmast, Schlatter, Ivanell and Henningson2016, Reference Kleusberg, Mikkelsen, Schlatter, Ivanell and Henningson2017; Kleusberg Reference Kleusberg2019) with Prandtl's tip correction. The code was validated extensively, showing good agreement between numerical and experimental results (Kleusberg et al. Reference Kleusberg, Mikkelsen, Schlatter, Ivanell and Henningson2017; Mühle et al. Reference Mühle2018; Kleusberg Reference Kleusberg2019; Kleusberg, Schlatter & Henningson Reference Kleusberg, Schlatter and Henningson2020). Since the focus of those studies was on the wake and on total forces, the errors caused by the approximation of the tip correction were not considered to greatly affect the results. The low dissipation and dispersion of the spectral-element method make it well-suited for stability analysis, hence this code was used previously in several vortex stability studies (Kleusberg Reference Kleusberg2019; Kleusberg, Benard & Henningson Reference Kleusberg, Benard and Henningson2019a; Kleine et al. Reference Kleine, Kleusberg, Hanifi and Henningson2019, Reference Kleine, Franceschini, Carmo, Hanifi and Henningson2022).

The iterative and direct smearing corrections of §§ 2 and 5, with the correction velocity detailed in § 3.3, were implemented in the code. For the iterative smearing correction, the circulation is considered converged when $\| \boldsymbol {\varGamma }^{new} - \boldsymbol {\varGamma }^{old} \|/\| \boldsymbol {\varGamma }^{new} \| < 10^{-5}$ (where $\| \boldsymbol {\varGamma } \|$ represents the Euclidean norm of the circulation vector). Also implemented is the analytical form of the convolution of the Gaussian function with a force distribution formed by segments of constant force presented in § 3.2.

For the NREL 5-MW reference turbine (Jonkman et al. Reference Jonkman, Butterfield, Musial and Scott2009), the aerofoil is defined only for a few spanwise positions in the original reference. In Martinez-Tossas et al. (Reference Martinez-Tossas, Churchfield, Yilmaz, Sarlak, Johnson, Sørensen, Meyers and Meneveau2018), two different approaches have been employed by different codes, one considering abrupt changes in the aerofoil, and the other considering interpolation of coefficients between sections. In the present implementation, the aerofoil data from Jonkman et al. (Reference Jonkman, Butterfield, Musial and Scott2009) are interpolated between sections, as performed by the DTU code EllipSys3D in Martinez-Tossas et al. (Reference Martinez-Tossas, Churchfield, Yilmaz, Sarlak, Johnson, Sørensen, Meyers and Meneveau2018), making the force curves smoother and avoiding problems with the discontinuities.

The non-iterative method requires the value of the lift coefficient slope. A third-order polynomial interpolation of $C_l$ is used to avoid discontinuities in the lift coefficient slope. The tabulated aerofoil lift coefficient is modelled using a shape-preserving piecewise cubic polynomial (The MathWorks, Inc. 2022).

A schematic representation of the computational domain can be seen in figure 4. Dirichlet boundary conditions are used for the inflow, upper, lower and lateral boundary conditions. The natural outflow boundary condition is imposed at the outlet.

Figure 4. Reference system and schematic view of the computational domain. The centre of the turbine/wing is equidistant to the upper, lower and lateral boundary conditions. Distances to the inflow ($L_{zin}$) and to the outflow ($L_{zout}$) may be different. The centre of the turbine/wing is located at the origin of the coordinate system.

6.2. Lifting line method

The nonlinear iterative lifting line method of § 4 was implemented, considering horseshoe vortices aligned with the $z$-direction. The same discretization of the actuator line method is used. Since the lifting line method is much faster to run, a stricter convergence criterion is used. The iteration is considered converged when the difference of circulation at every control point is lower than $10^{-8}$ of the mean absolute value of circulation:

(6.1)\begin{equation} \frac{\max_{1 \leq j \leq N} \left(|\varGamma_j^{new} - \varGamma_j^{old}| \right) }{\frac{1}{N}\sum_{j=1}^N|\varGamma_j^{new}|} < 10^{ - 8} . \end{equation}