1. Introduction
The actuator line method (ALM) was developed to represent blades or wings as lifting lines in numerical simulations of the Navier–Stokes equations (Sorensen & Shen Reference Sorensen and Shen2002). It allows the reduction of computational costs by replacing the geometry of the blades by lines carrying body forces calculated using the local velocity and aerofoil data (Sørensen et al. Reference Sørensen, Mikkelsen, Henningson, Ivanell, Sarmast and Andersen2015).
In order to calculate the forces accurately, the local velocity needs to be corrected near the tip of the blades. This correction has been a mildly controversial issue in the past, with different proposed models (Shen, Sørensen & Mikkelsen Reference Shen, Sørensen and Mikkelsen2005; Sørensen Reference Sørensen, Peinke and van Bussel2016) and even arguments that it should not be needed (Martinez et al. Reference Martinez, Leonardi, Churchfield and Moriarty2012; Sørensen Reference Sørensen, Peinke and van Bussel2016) due to the fact that the ALM creates tip vortices that should lead to lower loads near the tip. However, the overestimation of loads near the tip on the numerical results indicates that a correction is needed (Sørensen Reference Sørensen, Peinke and van Bussel2016; Dağ Reference Dağ2017).
Some recent discoveries have shed some light on this apparent controversy. Dağ & Sørensen (Reference Dağ and Sørensen2020) (originally in Dağ Reference Dağ2017) observed that the bound vortex created by the actuator line showed a Gaussian vorticity distribution equal to a Lamb–Oseen vortex model (Oseen Reference Oseen1911; Lamb Reference Lamb1932; Saffman Reference Saffman1992). The authors then conjectured that this pattern would be extended to the vortex sheet, and proposed a correction based on this model to approximate the velocity in the ALM to the velocity induced by the singular vorticity distribution predicted by a discretized Prandtl lifting line. The simulations of rotors showed clear improvements compared to the uncorrected ALM. The method brings the forces smoothly to zero near the tip and the hub of the blades. The results for a translating planar wing show that they approximate a lifting line method for this case. Meyer Forsting, Pirrung & Ramos-García (Reference Meyer Forsting, Pirrung and Ramos-García2019a) compared the results of a similar correction with the results of a lifting line technique for rotating blades with and without a viscous core model. They showed that the forces of the ALM without correction agree with a vortex method with finite core size, while the ALM with a vortex-based smearing correction agrees with the vortex method using ideal vortices.
The mathematical connection between the vortices generated by the ALM and a Lamb–Oseen vortex model has been proven by Forsythe et al. (Reference Forsythe, Lynch, Polsky and Spalart2015) for the bound vortex and by Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) (originally in Martínez Tossas Reference Martínez Tossas2017) for the vortex sheet of a translating wing. They showed that it is a consequence of the convolution of the discrete forces with a Gaussian kernel, necessary to distribute the forces and avoid numerical instabilities. The correction of Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019), termed ‘subfilter-scale velocity correction’ or ‘filtered actuator line model’, is, in essence, a variant of a vortex-based smearing correction. It was applied to a wind turbine by Stanly et al. (Reference Stanly, Martínez-Tossas, Frankel and Delorme2022).
Caprace, Chatelain & Winckelmans (Reference Caprace, Chatelain and Winckelmans2019) modified the Prandtl lifting line by distributing the vorticity using a Gaussian kernel. Even though this work does not concern the ALM directly, the connections between the mollified (smeared) lifting line and the ALM were stated clearly by the authors. They studied the effect of three-dimensional, two-dimensional (normal to the line) and one-dimensional (normal to the line and streamwise directions) regularizations. Both the three-dimensional and two-dimensional regularizations are used in the ALM; see Mikkelsen (Reference Mikkelsen2003). In that work, a variable smoothing parameter approaching zero or a very low finite value near the tip leads to improved results, a result also observed in ALM works (Shives & Crawford Reference Shives and Crawford2013; Jha et al. Reference Jha, Churchfield, Moriarty and Schmitz2014; Jha & Schmitz Reference Jha and Schmitz2018). However, in the ALM, a low value of smoothing parameter 
$\varepsilon$ may lead to numerical instabilities. Since the minimum 
$\varepsilon$ is usually around 2–3 times the grid spacing (Troldborg Reference Troldborg2009), a variable smoothing parameter may impose a very refined grid near the tips.
Other methods may allow a more compact representation of the effect of the forces. For example, one strategy in finite-volume methods to represent an actuator disk or surface is to model the effects of the forces as pressure jumps, allowing it to be restricted to one grid element (Réthoré & Sørensen Reference Réthoré and Sørensen2012; Troldborg et al. Reference Troldborg, Sørensen, Réthoré and van der Laan2015). Also, for vortex particle methods, which shed vortex elements in the wake, the smoothing parameter may be equal to the grid spacing (Caprace, Winckelmans & Chatelain Reference Caprace, Winckelmans and Chatelain2020). Nevertheless, in numerical methods, the compactness of the representation of the effects of the forces is limited by the grid spacing and is guided by numerical reasons, not by the physical properties of the flow.
Depending on the numerical method or application, a larger 
$\varepsilon$ may be necessary or desirable in the ALM. For example, Kleusberg, Schlatter & Henningson (Reference Kleusberg, Schlatter and Henningson2019b), using a spectral-element method, observed spatially growing oscillations for 
$\varepsilon =2 \Delta x$ (where 
$\Delta x$ is the average grid spacing), while the amplitude of oscillations was bounded for 
$3 \Delta x \leq \varepsilon \leq 4 \Delta x$. The oscillations for 
$\varepsilon =2 \Delta x$ were not sufficient to cause numerical instabilities on force calculations; however, they may affect applications that require low numerical disturbances in the wake. For this reason, Kleine et al. (Reference Kleine, Franceschini, Carmo, Hanifi and Henningson2022) used 
$\varepsilon =3.5 \Delta x$ for a vortex stability study. Also, Shives & Crawford (Reference Shives and Crawford2013) recommended a smearing parameter at approximately 
$\varepsilon =4 \Delta x$ if lower errors in the angle of attack are desirable (see also Meyer Forsting & Troldborg Reference Meyer Forsting and Troldborg2020).
For this reason, corrections that take into account the difference between the vorticity created by the convoluted forces and the vorticity created by discrete singular forces, such as proposed by Dağ & Sørensen (Reference Dağ and Sørensen2020), Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) and Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a), seem to be more cost-effective than reducing 
$\varepsilon$. Some techniques to further reduce the computational cost of the smearing correction are explored by Meyer Forsting, Pirrung & Ramos-García (Reference Meyer Forsting, Pirrung and Ramos-García2020).
Meyer Forsting, Pirrung & Ramos-García (Reference Meyer Forsting, Pirrung and Ramos-García2019b) investigated the wake created by the actuator line with smearing correction, confirming that the smearing parameter continues to have an influence on the wake, despite its influence on the forces at the blades being greatly reduced.
The corrections of Dağ & Sørensen (Reference Dağ and Sørensen2020) and Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a) apply an iterative process at each time step. The circulation is dependent on the local induced velocity, which depends on the circulation. For this reason, at each iteration, the velocity or circulation is updated using a relaxation iterative process. From our experience, the choice of a relaxation factor close to unity can lead to numerical instability, while a low relaxation factor increases the number of iterations and can increase the run time of the simulation.
The correction of Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) avoids an iterative procedure at each time step by using the circulation of the previous time step to calculate the induced velocity. Then a weighted average of the correction velocity of the previous time step and the correction velocity of the current time step is used as the correction term. Hence the induced velocities are calculated from the values of circulation of the previous two time steps, without taking into consideration the current value of circulation. Conceptually, it is equivalent to a first iteration of the iterative method of other works. For steady simulations, the values converge to the steady solution. However, for unsteady problems, this procedure does not guarantee compatibility between the current circulation and local velocity at each time step. This approach can be justified due to the small difference between the circulations between time steps, for most simulations. Nevertheless, an error is introduced, which is dependent on the difference of the circulation between time steps and the weighting factor (in that work, called the ‘relaxation factor’).
In this work, a method is introduced that avoids the iterative method while maintaining the compatibility between the current circulation and velocities. We propose a direct way of computing the smearing correction, based on a linearized version of the lifting line. Using this method, the correction at each time step is found by solving a linear system of equations of the order 
$N$, where 
$N$ is the number of actuator line control points. In order to achieve this, two formulations of the discretized lifting line method based on the actuator line are presented: the nonlinear formulation and its linearized version.
Additionally, some other aspects of the smearing correction are discussed, such as the choice of correction function and formation of the vortex sheet. Regarding these aspects, we aim to keep the method as general as possible. Previous works (Martínez-Tossas & Meneveau Reference Martínez-Tossas and Meneveau2019; Meyer Forsting et al. Reference Meyer Forsting, Pirrung and Ramos-García2019a; Dağ & Sørensen Reference Dağ and Sørensen2020) already showed that the smearing correction can make the forces calculated by the ALM match closely the results of the lifting line method. By introducing fewer assumptions, avoiding especially assumptions related to rotating blades, the method could be applied to other problems beyond rotor simulations, such as simulations of fixed-wing aircrafts.
This work is structured as follows. First, we present the general idea of the actuator line method with smearing correction in § 2. Then, in § 3, we draw on the theoretical work of Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) to develop the specific correction for velocities induced by vortex segments, and describe the formation of the vortex sheet. The iterative lifting line method based on the actuator line is presented in § 4. The linearization of this lifting line method (Appendix A) is the basis for the non-iterative smearing correction detailed in § 5. The numerical solver employed in the simulations and the parameters of the lifting line method are described in § 6. Results of the simulations of a translating wing and the NREL 5-MW wind turbine are shown in § 7. Finally, the main conclusions are summarized in § 8.
2. The actuator line method
The incompressible Navier–Stokes equation written in primitive variables (pressure 
$p$ and velocity 
$\boldsymbol {u}$) is
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} = - \frac{1}{\rho}\,\boldsymbol{\nabla} p + \nu\,\nabla^2 \boldsymbol{u} + \boldsymbol{f}, \end{equation}
where 
$\rho$ and 
$\nu$ are the density and kinematic viscosity of the fluid, respectively. The body force term 
$\boldsymbol {f}$, in the case of the actuator line method, is used to model the turbine (Sorensen & Shen Reference Sorensen and Shen2002; Mikkelsen Reference Mikkelsen2003; Troldborg Reference Troldborg2009). The body forces are based on the two-dimensional force per spanwise unit length, 
$\boldsymbol {F}_{2D}$, given by
\begin{equation} \boldsymbol{F}_{2D} = (F_l,F_d) = ( \tfrac{1}{2}\rho u_r^2 c C_l,\tfrac{1}{2}\rho u_r^2 c C_d ), \end{equation}
where 
$F_l$ and 
$F_d$ are the lift and drag forces (lift is perpendicular to the relative velocity, and drag is parallel to the relative velocity; see figure 1), calculated from the relative velocity 
$u_r=\sqrt {u_y^2+u_z^2}$, the local chord 
$c$, and the two-dimensional lift and drag coefficients 
$C_l$ and 
$C_d$, which are obtained from the aerofoil data at the local Reynolds number and angle of attack 
$\alpha$ (calculated using the local relative velocity).
Figure 1. Local reference system defined by a cross-section of the blade.
The dimensions of 
$\boldsymbol {F}_{2D}$ are not consistent with the dimensions of the body force 
$\boldsymbol {f}$. As an intermediate step, an idealized three-dimensional body force 
$\boldsymbol {f}^i$ is defined based on 
$\boldsymbol {F}_{2D}$, where at each cross-section, the body force is given by
\begin{equation} \boldsymbol{f}^i = \frac{1}{\rho}\,\boldsymbol{F}_{2D}\,\delta(y)\,\delta(z) , \end{equation}
where 
$\delta$ is the Dirac delta function, which can be interpreted as the limit of the Gaussian function
\begin{equation} \delta(y) = \lim_{\varepsilon \to 0} \frac{1}{{\rm \pi}^{1/2}\varepsilon} \exp\left({-\frac{y^2}{\varepsilon^2}}\right). \end{equation}
This body force 
$\boldsymbol {f}^i$, where the superscript 
$i$ indicates the idealized case, is concentrated at the origin of the local reference system (figure 1). Considering all cross-sections of a wing, the space of non-zero 
$\boldsymbol {f}^i$ defines a line, which can be interpreted as the limit of the actuator line method, when the smearing parameter 
$\varepsilon$ goes to zero. It is easy to see that 
$\boldsymbol {f}^i$ is singular at the origin, but its integration in the two-dimensional plane gives 
$\boldsymbol {F}_{2D}/\rho$.
In the ALM, to avoid numerical problems related to singularities, the forces need to be distributed smoothly on several mesh points. This is usually performed by convolving the force with a regularization kernel. A common choice for the regularization kernel is a three-dimensional Gaussian kernel. Considering the same constant Gaussian width in the three directions,
\begin{equation} \eta_3(x,y,z) := \frac{1}{{\rm \pi}^{3/2}\varepsilon^3} \exp{\left(-\frac{x^2+y^2+z^2}{\varepsilon^2}\right)}=\eta(x)\,\eta(y)\,\eta(z), \end{equation}where
\begin{equation} \eta(x) := \frac{1}{{\rm \pi}^{1/2}\varepsilon} \exp{\left(-\frac{x^2}{\varepsilon^2}\right)}, \end{equation}
and symbol 
$:=$ denotes equal to by definition. The smeared force is
\begin{equation} \boldsymbol{f} = \boldsymbol{f}^i * \eta_3 . \end{equation}Non-uniform and anisotropic kernels have also been proposed (Mikkelsen Reference Mikkelsen2003; Shives & Crawford Reference Shives and Crawford2013; Churchfield et al. Reference Churchfield, Schreck, Martinez, Meneveau and Spalart2017; Martínez-Tossas, Churchfield & Meneveau Reference Martínez-Tossas, Churchfield and Meneveau2017; Jha & Schmitz Reference Jha and Schmitz2018; Cormier, Weihing & Lutz Reference Cormier, Weihing and Lutz2021), but these are not investigated in the present work.
This presentation of the convolution operation in the actuator line method using an idealized body force formed by Dirac delta functions as an intermediate step may seem like an unusual approach, but it is mathematically equivalent to the method of classical works (Sorensen & Shen Reference Sorensen and Shen2002; Mikkelsen Reference Mikkelsen2003; save for possible typos regarding the density term, which is not relevant for incompressible flows), which employ one-dimensional integrals. This procedure, which was already applied by Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019), formalizes the connection between the two-dimensional forces and the convolution in three dimensions, which is relevant for § 3.
Similarly to previous studies (Martínez-Tossas & Meneveau Reference Martínez-Tossas and Meneveau2019; Meyer Forsting et al. Reference Meyer Forsting, Pirrung and Ramos-García2019a; Dağ & Sørensen Reference Dağ and Sørensen2020), we focus on the lift force, which has a greater influence on the forces of the turbine, and leave the effects related to the smearing of the drag force for future studies. (A two-dimensional correction for drag has been proposed by Martínez-Tossas et al. (Reference Martínez-Tossas, Churchfield and Meneveau2017), and an outline for a three-dimensional drag force correction can be found in Kleine (Reference Kleine2022).)
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
\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
\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
\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
\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$):
\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:
\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
\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.
(i) Start from a guess of the circulation distribution, usually the circulation from the previous time step.
(ii) Form the vortex sheet by:
(a) prescribing the vortex sheet, for example, by assuming helical or horseshoe vortices; or
(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.
(iii) Calculate the missing velocity
$\boldsymbol {u}^m(\boldsymbol {x}_j)$ at every control point. Find the local corrected velocities according to (2.11).(iv) Calculate the effective angle of attack using (2.9).
(v) Find the local lift coefficient
$C_l$, interpolating from the aerofoil data table.(vi) Calculate the new value of circulation
$\varGamma _j^{new}$ using (2.10).(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}(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. Vorticity created by body forces and the missing velocity
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
\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
\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
\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}
\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
\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,
\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
\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
\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
\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
\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))
\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}
\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
\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}
\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
\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
\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
\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
\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
\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
\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
\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
\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
\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
\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.
4. Iterative lifting line consistent with the actuator line method
The formulation of the actuator line presented in §§ 2.1 and 2.3 is very similar to a nonlinear lifting line method (Anderson Reference Anderson1991; see also Phillips & Snyder Reference Phillips and Snyder2000). Hence we can develop a numerical nonlinear lifting line method consistent with the actuator line method with vortex-based smearing correction. While Anderson (Reference Anderson1991) uses the undisturbed uniform inflow velocity (
$\boldsymbol {U}_{\infty }$) to calculate the circulation, we use the local velocity, basing our method on (2.10), in an approach similar to that of Phillips & Snyder (Reference Phillips and Snyder2000). Also, in this formulation, the undisturbed velocity 
$\boldsymbol {U}=(U_x,U_y,U_z)$ (velocity as if the wing were not present in the flow) can change along the span, a feature paramount for the application of the actuator line to rotating blades. The local velocity 
$\boldsymbol {u}$, which is previously unknown, can be calculated from the local undisturbed velocity 
$\boldsymbol {U}$ and the velocity induced by the vortex system, 
$\boldsymbol {u}^{vi}$, at control point 
$\boldsymbol {x}_j$, given by
\begin{equation} \boldsymbol{u}(\boldsymbol{x}_j) = \boldsymbol{U}(\boldsymbol{x}_j) + \boldsymbol{u}^{vi}(\boldsymbol{x}_j) , \end{equation}The steps of the nonlinear lifting line are as follows.
(i) Start from a guess of the circulation distribution, for example, applying (2.10) with
$\boldsymbol {u}=\boldsymbol {U}$.(ii) Form the vortex sheet by:
(a) prescribing the vortex sheet, for example, by assuming helical or horseshoe vortices; or
(b) employing a free-vortex wake method.
(iii) Calculate the velocity induced by the bound vortex and the vortex sheet at every control point,
$\boldsymbol {u}^{vi}(\boldsymbol {x}_j)$. Find the local velocity according to (4.1).(iv) Calculate the effective angle of attack using (2.9).
(v) Find the local lift coefficient
$C_l$, interpolating from the aerofoil data table.(vi) Calculate the new value of circulation
$\varGamma _j^{new}$ using (2.10).(vii) Update the current value of circulation using a relaxation factor
$r$:
(4.2)\begin{equation} \varGamma_j=r\varGamma_j^{new} + (1-r) \varGamma_j^{old}. \end{equation}(viii) Restart from step (ii) if the local velocity affects the formation of the vortex sheet or from step (iii) if a prescribed vortex sheet is considered, using the value of circulation calculated in the previous step. Iterate until a chosen convergence criterion is reached.
This iterative lifting line method was used as the reference to compare the results of the ALM in § 7.1, with the parameters detailed in § 6.2. This formulation is able to deal with configurations in which the induced velocity can no longer be considered small compared to the undisturbed velocity, as observed in Kleine et al. (Reference Kleine, Hanifi and Henningson2023). Based on the same equations, a linearized lifting line method is developed in Appendix A, which is used to derive the non-iterative smearing correction.
5. Non-iterative smearing correction based on the linearized lifting line method
The steps to apply the smearing correction (§ 2.3) are basically equal to the steps of the nonlinear lifting line described in § 4, except for step (iii), which becomes the calculation of the missing velocities, and step (ii), which allows the formation of the vortex sheet using data from the CFD simulation. Therefore, a non-iterative smearing correction based on the formulas of the linearized lifting line described in Appendix A can be devised.
A minor difference in the process of the iterative lifting line and the smearing correction is the possibility to start from a previous time step, which provides a better guess for the circulation. This minor difference, however, is what allows the linearized, non-iterative formulation of the smearing correction to have good results even if the aerofoil lift coefficient is not linear with respect to the angle of attack.
The first phase of the method is to find an approximate solution to linearize around. This is achieved by applying the first three iterative steps of § 2.3 only once. The corrected velocity calculated from the circulation distribution 
$\boldsymbol {\varGamma }^{n-1}$ from the previous time step 
$n-1$ is
\begin{equation} \boldsymbol{u}^{c}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n-1}) = \boldsymbol{u}^{s}(\boldsymbol{x}_j) + \boldsymbol{u}^{m}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n-1}) . \end{equation}
Representing the variables around which the functions are linearized by superscript 
${\dagger}$, we have
\begin{gather} u_y^{{{\dagger}}} := u_{y}^{c}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n-1}) = u_{y}^{s}(\boldsymbol{x}_j) + u_{y}^{m}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n-1}), \end{gather}
\begin{gather} u_z^{{{\dagger}}} := u_{z}^{c}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n-1}) = u_{z}^{s}(\boldsymbol{x}_j) + u_{z}^{m}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n-1}) . \end{gather}Analogously to Appendix A,
\begin{equation} u_r^{\dagger} := \sqrt{(u_z^{{\dagger}2}+u_y^{{\dagger}2})}, \quad \alpha_{j}^{\dagger} := \alpha_{gj} + \arctan{\left(\frac{u_y^{\dagger}}{u_z^{\dagger}}\right)} \quad \text{and} \quad \varGamma_{j}^{\dagger} := \frac{1}{2}\,u_r^{\dagger} \, c_j \, C_l(\alpha_{j}^{\dagger}) . \end{equation}Phase two of the linearization process involves finding another linear relation between the corrected velocity and the circulation. It is assumed that the missing velocity is not relevant for the formation of the vortex sheet, as discussed in § 3.4, and the CFD velocity is used directly to advect the vortices. This implies that the current value of circulation does not affect the geometry of the vortex sheet, which simplifies the method and is consistent with a first-order approximation.
The unknown missing velocity for the value of circulation at the current time step, 
$\boldsymbol {\varGamma }^{n}$, can be divided into two components:
\begin{equation} \boldsymbol{u}^{m}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n}) = \boldsymbol{u}^{m_{p}}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n}) + \boldsymbol{u}^{m_{c}}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n}), \end{equation}
where 
$\boldsymbol {u}^{m_{p}}$ is the velocity induced by the vortex sheet emitted in the previous time steps, and 
$\boldsymbol {u}^{m_{c}}$ is the velocity induced by the vortex sheet emitted in the current time step and by the bound vortices. The value of circulation at the current time step does not affect the vortex wake emitted in the previous time steps, hence the velocity induced by its vortices, 
$\boldsymbol {u}^{m_{p}}$, is already computed. Considering the linearity of the relationship between velocity and circulation, (5.5) can be rewritten as
\begin{equation} \boldsymbol{u}^{m}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n}) = \boldsymbol{u}^{m_{p}}(\boldsymbol{x}_j) + \boldsymbol{u}^{m_{c}}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n-1}) + \Delta \boldsymbol{u}^{m_{c}}(\boldsymbol{x}_j,\Delta \boldsymbol{\varGamma}) = \boldsymbol{u}^{m}(\boldsymbol{x}_j,\boldsymbol{\varGamma}^{n-1}) + \Delta \boldsymbol{u}^{m_{c}}(\boldsymbol{x}_j,\Delta \boldsymbol{\varGamma}), \end{equation}
where 
$\Delta \boldsymbol {u}^{m_{c}} = \boldsymbol {u}^{m_{c}}(\boldsymbol {x}_j,\boldsymbol {\varGamma }^{n})-\boldsymbol {u}^{m_{c}}(\boldsymbol {x}_j,\boldsymbol {\varGamma }^{n-1})$ and 
$\Delta \boldsymbol {\varGamma }=\boldsymbol {\varGamma }^{n}-\boldsymbol {\varGamma }^{n-1}$. A schematic representation is shown in figure 3.
Figure 3. Schematic representation of (5.6). The missing velocity is the sum of the missing velocity due to the vortex sheet emitted in the previous time steps (leftmost figure of the second row) and the missing velocity due to the bound vortex and the vortex sheet emitted in the current time step, which can also be split into two terms (central and rightmost figures of the second row).
We write the velocities induced at control point 
$j$ by the vortex system 
$k$, which includes the bound vortex and the vortex sheet created by the segment 
$k$ at the current time step, as
\begin{gather} u_{y}^{m_{ck}}(\boldsymbol{x}_j,\varGamma_k^{n}) = a_{y,jk}^{m_c} \varGamma_k^{n}, \end{gather}
\begin{gather} u_{z}^{m_{ck}}(\boldsymbol{x}_j,\varGamma_k^{n}) = a_{z,jk}^{m_c} \varGamma_k^{n}, \end{gather}
where 
$a_{y,jk}^{m_c}$ and 
$a_{y,jk}^{m_c}$ are the influence coefficients obtained considering only the bound vortex and the vortex sheet emitted at the current time step (see figure 3). It should be noted that the coefficients 
$a_{y,jk}^{m_c}$ and 
$a_{y,jk}^{m_c}$ should already be known to perform step (iii) of the usual iterative version of the smearing correction (§ 2.3). Equations (5.7) and (5.8) written in matrix form are
\begin{gather} \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{y}}}^{m_{c}}(\boldsymbol{\varGamma}^{n}) = \boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{y}}}^{m_c} \boldsymbol{\varGamma}^{n}, \end{gather}
\begin{gather} \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{z}}}^{m_{c}}(\boldsymbol{\varGamma}^{n}) = \boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{z}}}^{m_c} \boldsymbol{\varGamma}^{n} . \end{gather} The missing velocity 
$\boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{y}}}^{m}(\boldsymbol {\varGamma }^{n})$ is then written as
\begin{equation} \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{y}}}^{m}(\boldsymbol{\varGamma}^{n}) = \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{y}}}^{m}(\boldsymbol{\varGamma}^{n-1}) + \boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{y}}}^{m_c}\,\Delta\boldsymbol{\varGamma} . \end{equation}
Combining (5.2) and (5.11), the corrected velocity 
$\boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{y}}}^{c}$ can be written as
\begin{equation} \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{y}}}^{c}(\boldsymbol{\varGamma}^{n}) = \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{y}}}^{\dagger} + \boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{y}}}^{m_c}\, \Delta\boldsymbol{\varGamma}. \end{equation}
Analogously for 
$\boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{z}}}^{c}$:
\begin{equation} \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{z}}}^{c}(\boldsymbol{\varGamma}^{n}) = \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{z}}}^{\dagger} + \boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{z}}}^{m_c}\, \Delta\boldsymbol{\varGamma} . \end{equation}Finally, the linear equation can be found following the steps of Appendix A. Similar to (A1), the linearization of (2.10) becomes
\begin{align} \varGamma_j^{n} &=
\varGamma_{j}^{\dagger} + \frac{1}{2}\,c_j \left[ \left(
C_l(\alpha_{j}^{\dagger})\, \frac{u_y^{\dagger}}{u_r^{\dagger}} +
\frac{\partial C_l}{\partial \alpha}(\alpha_{j}^{\dagger})\,
\frac{u_z^{\dagger}}{u_r^{\dagger}} \right) \Delta
u_{y}^{m_{c}}\right.\nonumber\\ &\left.\quad +\, \left( C_l(\alpha_{j}^{\dagger})\,
\frac{u_z^{\dagger}}{u_r^{\dagger}} - \frac{\partial
C_l}{\partial \alpha}(\alpha_{j}^{\dagger})\,
\frac{u_y^{\dagger}}{u_r^{\dagger}} \right) \Delta
u_{z}^{m_{c}} \right] ,
\end{align}which can be written in matrix form as
\begin{equation} (\boldsymbol{\mathsf{I}} - ({\rm diag}(\boldsymbol{\mathsf{b}}_{\boldsymbol{\mathsf{y}}}^{\dagger})\,\boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{y}}}^{m_c}) - ({\rm diag}(\boldsymbol{\mathsf{b}}_{\boldsymbol{\mathsf{z}}}^{\dagger})\,\boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{z}}}^{m_c} ) ) \Delta\boldsymbol{\varGamma} = \boldsymbol{\varGamma}^{\dagger} - \boldsymbol{\varGamma}^{n-1}, \end{equation}
where 
$\boldsymbol{\mathsf{b}}_{\boldsymbol{\mathsf{y}}}^{\dagger}$ and 
$\boldsymbol{\mathsf{b}}_{\boldsymbol{\mathsf{z}}}^{\dagger}$ can be found using the corresponding versions of (A5) and (A6). Here, 
$\Delta \boldsymbol {\varGamma }$ is found by solving the linear system, and the corrected velocities are calculated using (5.12) and (5.13).
If the lift coefficient has a linear relation with the angle of attack, then this method is expected to give good results, even if the differences in the circulation are high, just like a conventional lifting line method has good results despite starting from a completely wrong circulation distribution. However, the strength of the method is in the fact that the functions are linearized around the velocities calculated at the first iteration of a conventional smearing correction. The differences in the velocities and circulation are expected to be small in most simulations, justifying this linearization, even if the lift coefficient is not linear. The method is expected to perform poorly only if the lift coefficient slope changes with the angle of attack and the flow conditions vary greatly between time steps. For these extreme cases, the linear method presented here can be used to accelerate the convergence of an iterative method. The linear system is relatively small, of size 
$N$, usually several orders of magnitude smaller than the CFD solver, which makes its solution computationally inexpensive.
For multiple turbines, each turbine can be solved in isolation, since the distance between turbines is much greater than 
$\varepsilon$. In fact, if the distance between the blades is much larger than 
$\varepsilon$ or if forces at the hub are not relevant, then the system of equations can be divided and each blade solved independently. This last simplification was not implemented for the results of § 7.
All the aspects regarding the formation of the vortex sheet and calculation of the induced velocity are accounted for in matrices 
$\boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{y}}}^{m_c}$ and 
$\boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{z}}}^{m_c}$. This means that the method of previous works (Martínez-Tossas & Meneveau Reference Martínez-Tossas and Meneveau2019; Meyer Forsting et al. Reference Meyer Forsting, Pirrung and Ramos-García2019a, Reference Meyer Forsting, Pirrung and Ramos-García2020; Dağ & Sørensen Reference Dağ and Sørensen2020; Stanly et al. Reference Stanly, Martínez-Tossas, Frankel and Delorme2022) could benefit from the procedure outlined in this section, just by forming matrices 
$\boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{y}}}^{m_c}$ and 
$\boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{z}}}^{m_c}$ compatible with their respective methods. In particular, if desirable, all of the methods for increasing the speed of calculations detailed in Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2020) could be implemented easily and would affect only the coefficients of these matrices.
For the methods that prescribe the shape of the vortex sheet, the separation of the missing velocity in the components of (5.5) could be skipped; however, it would reduce errors relative to the linearization if implemented. This means that the non-iterative strategy should also work even without bookkeeping of previous values of circulation. A mixed strategy that could be used to reduce the cost is to use a prescribed wake to calculate 
$\boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{y}}}^{m}(\boldsymbol {\varGamma }^{n-1})$, considering the circulation along the whole wake as 
$\boldsymbol {\varGamma }^{n-1}$ (substituting older values of 
$\boldsymbol {\varGamma }$ in figure 3 by 
$\boldsymbol {\varGamma }^{n-1}$) and another prescribed near wake to define 
$\Delta \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{y}}}^{m_{c}}$, reducing the magnitude of the influence coefficients (compared to the case without any bookkeeping).
6. Numerical method
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:
\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}7. Results
7.1. Comparison with the lifting line method
The results of the simulation of the ALM with smearing correction were compared with the results of a nonlinear iterative lifting line for a planar straight wing with constant chord in uniform flow. Because the flow has reached the steady state, no difference was expected between the iterative and direct methods for this case, which is confirmed by the results. Nevertheless, we present the results for both strategies for completeness.
The parameters of the geometry and simulation are presented in table 1. Except where explicitly stated, all values are non-dimensionalized by the span (
$R$), the velocity at infinity (
$U_{ref}=U_z$) and the density (
$\rho$). An ideal aerofoil without drag and with lift coefficient 
$C_l = 2 {\rm \pi}\alpha$ is chosen (considering the standard non-dimensionalization of 
$C_l$ using the chord). The angle of attack 
$\alpha _g=1/(2{\rm \pi} )$ is such that a lift coefficient 
$C_l(\alpha _g)=1$ and circulation 
$\varGamma _0/(R U_z)=0.05$ would be expected for a two-dimensional simulation (case with zero induced velocity).
Table 1. Parameters of the simulation of the wing with constant chord in uniform flow.
For the ALM simulations, the Reynolds number based on the chord is 
$Re_c = c U_z/\nu = 10^4$ (where 
$\nu$ is the kinematic viscosity). The Reynolds number is not applicable to the lifting line method.
The average grid spacing in the region of the actuator line is 
$\Delta x=R/56$. The results for two values of smearing parameter are compared: 
$\varepsilon = 3.5 \Delta x = R/16$ and 
$\varepsilon = 7 \Delta x = R/8$. The lower chosen value is higher than the usual minimum (
$2 \Delta x \leq \varepsilon \leq 3 \Delta x$), but it provides a better reference for validation because it allows better resolution of the vortex core, a choice also adopted by Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a). It is also of practical interest: 
$\varepsilon = 3.5 \Delta x$ has been used in vortex stability studies (Kleine et al. Reference Kleine, Franceschini, Carmo, Hanifi and Henningson2022), based on parametric studies that showed that a larger smearing parameter creates lower numerical oscillations (Kleusberg et al. Reference Kleusberg, Schlatter and Henningson2019b). The last value, 
$\varepsilon = 7 \Delta x$, is included only as a reference, in order to evaluate the correction for a very large kernel, but large kernels are not usually employed in practice (as far as we are aware). The case 
$\varepsilon = 2 \Delta x = R/28$ is discussed in Appendix B, where the effect of grid resolution is analysed.
Because of symmetry, just the results for half of the wing are shown. As can be seen in figure 5, the agreement is excellent for all cases. The difference in the induced velocity 
$u_y$ is of the order of 
$10^{-4}$ (
$0.01\,\%$ of the undisturbed velocity), an agreement much better than those reported by Dağ & Sørensen (Reference Dağ and Sørensen2020) and Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a). The difference is of the order of the square of 
$u_y$, which is consistent with a first-order method, and therefore a better agreement was not expected. The other works are mentioned not as a criticism, but to bring context to the current values, showing the strength of the method. For most practical applications of the actuator line, the errors of the other methods can also be considered negligible.
Figure 5. Comparison of the results for the straight wing in uniform flow: LL indicates nonlinear iterative lifting line method of § 4; ALM
$^d$ indicates ALM with non-iterative correction of § 5; ALM
$^i$ indicates ALM with iterative correction of § 2.3. (a) Corrected velocity in the 
$y$-direction. (b) Difference between 
$y$-velocities obtained by the ALM and LL. (c) Corrected velocity in the 
$z$-direction. (d) Difference between 
$z$-velocities obtained by the ALM and LL. (e) Circulation. (f) Relative difference between circulation obtained by the ALM and LL.
We believe that the remarkable agreement can be explained mainly by three factors. First, the convolution of the force with the Gaussian kernel was performed analytically using (3.9), considering constant force in a segment, which is the case that is mathematically equivalent to the lifting line method. Second, the lifting line method was carefully constructed to be consistent with the actuator line method. Third, the use of the velocity induced by a smeared vortex segment of § 3.3 has lower associated errors when compared to other strategies where the vortices are also discretized.
The difference between the results for the two values of the smearing parameter is negligible. This is exactly the aim of the smearing correction: to make the forces on the blades insensitive to variations of the smearing correction and approach the lifting line method, which is the limit 
$\varepsilon \to 0$.
The velocity in the 
$z$-direction is also shown in figure 5, because it influences the circulation. The effect of the smearing correction on 
$u^c_z$ is minimal. For most points, the missing velocity 
$u^m_z$ is at least one order of magnitude lower than the difference observed in figure 5(d). Hence the differences observed are due mainly to the velocity sampled from the numerical simulation. The reason for the differences is probably the discrepancies in the vortex sheet created in the CFD domain and the vortex sheet of the lifting line method, as mentioned in the first paragraph of § 3.4. The prescribed horseshoe vortices of the lifting line method do not induce velocity in the 
$z$-direction. In the CFD simulation, the vortices are free. As a consequence, the downwash from the wing inclines the vortex sheet, which induces velocity in the negative 
$z$-direction. A more advanced lifting line method might be able to capture those effects. Hence the differences in 
$u_z$ (and its consequence in 
$\varGamma$) can be, at least partially, attributed to a limitation of the current implementation of the lifting line method.
The relative differences in the values of circulation are of the order of 
$10^{-3}$ (differences of the order of 
$0.1\,\%$ of 
$\varGamma _0$), which is a consequence of the errors in 
$u_y$ and 
$u_z$. This difference is negligible for practical purposes. The circulation for the higher value of the smearing parameter agrees better with the lifting line. This apparently contradictory result has two simple explanations. First, the vortex core is larger, so the representation of the vorticity is better, reducing numerical errors (see Appendix B). Second, for higher values of 
$\varepsilon$, the results depend less on the data from the simulation and more on the smearing correction. In the limit 
$\varepsilon \to \infty$, no information from the numerical simulation is used, and the force is determined only by the correction, which regresses identically to the lifting line method. Nevertheless, increasing the smearing parameter is not usually beneficial, because the value of 
$\varepsilon$ has other effects on the simulation beyond its influence on the forces, as discussed in § 7.2. It should be noted, however, that the lifting line method has its own limitations on reproducing the results of wings and blades. Notwithstanding, it is a relevant tool for middle- or low-fidelity aerodynamic calculations.
Before the development of the smearing correction, the use of the ALM was restricted almost exclusively to simulations of rotating blades, probably due to its inability to reproduce accurately the induced velocities, which affected the lift of translating wings. The remarkable agreement achieved here, which further confirms the good agreement of previous works (Martínez-Tossas & Meneveau Reference Martínez-Tossas and Meneveau2019; Meyer Forsting et al. Reference Meyer Forsting, Pirrung and Ramos-García2019a; Dağ & Sørensen Reference Dağ and Sørensen2020), might encourage the adoption of the ALM for other applications beyond rotating blades. For example, the aeronautics community could use it to simulate wings or the lifting surfaces of aircrafts. The ALM enables the integration of lifting lines with a CFD solver, allowing more complex configurations and flow conditions than a classical lifting line method. The low grid requirement is a clear advantage over traditional CFD methods that employ wall boundary conditions for the wing. Also, the free-vortex method described here maintains the generality of the method, so diverse configurations are allowed, and it does not suffer from the instability problems of traditional vortex filament methods (Leishman Reference Leishman2000).
7.2. Unsteady rotor simulations
The NREL 5-MW wind turbine subjected to a sheared inflow was simulated using the iterative and direct smearing corrections. The undisturbed wind speed varies in the vertical direction (
$y$), in order to simulate a condition in which the circulation of the blade changes every time step. Hence for every time step of this simulation, the initial guess of circulation is different from the final value of circulation, for both the iterative and non-iterative corrections.
The domain was reduced to 
$L_x/R=L_y/R=8$ to make the inflow velocity 
$U_{z}=y/5+1$ always positive and avoid problems with the outflow boundary condition, which is not stable if the velocity is negative. The distance from the inflow to the blades, 
$L_{zin}/R=4$, was also reduced, maintaining the distance from the blades to the outflow, 
$L_{zout}/R=6$. Since this is a comparative study, the reduction of the domain does not affect any analysis.
The parameters of the simulation are presented in table 2. Except where stated explicitly, all values are non-dimensionalized by the radius (
$R$), the velocity at the centre of the turbine (
$U_{ref}=U(0,0,0)$) and the density (
$\rho$). Hence the geometry, detailed in Jonkman et al. (Reference Jonkman, Butterfield, Musial and Scott2009), was non-dimensionalized by the radius. For this choice of reference values, the tip speed ratio (which is the ratio of 
$\varOmega R$ and the velocity at the centre of the turbine) is 
$\varOmega R/U_{ref}=7.55$. The average grid spacing in the region of the actuator line is the same as the simulation of the straight wing, 
$\Delta x=R/56$. The Reynolds number based on the radius is 
$Re_R = R U_{ref}/\nu = 5 \times 10^4$. A time step 
$\Delta t = T/400$ was used, where 
$T$ is the period of rotation.
Table 2. Parameters of the simulation of the NREL 5-MW wind turbine in sheared flow.
In figure 6, the circulation and forces at time 
$t=12T$ are compared for the blade that has null azimuthal angle at this time (blade aligned with the 
$x$-axis). The evolution in time of the circulation of the control point closest to the tip of this blade is shown in figure 7.
Figure 6. Circulation and forces along the radial direction (
$r$) at time 
$t/T=12$, for the blade positioned at null azimuthal angle (blade aligned with the 
$x$-axis). (a) Lift force. (b) Drag force. (c) Circulation. (d) Difference of circulation, taking as reference the circulation of ALM
$^d$ with 
$\varepsilon = 3.5 \Delta x$, normalized by the maximum of the circulation.
Figure 7. Evolution in time of the circulation at the point closest to the tip. The data are saved with a period of 
$T/20$, which is higher than the actual time step in the simulation. The correction was turned off for the first two points.
As can be seen in figures 6 and 7, the agreement is very good for all cases. The iterative and direct methods have results that are practically identical. The differences in the value of circulation are lower than 
$10^{-5}$ if the same smearing parameter 
$\varepsilon$ is used. This confirms that the non-iterative procedure does not introduce errors for this unsteady flow.
The differences between the circulation for the two values of 
$\varepsilon$ are of the order of 
$1\,\%$ of the value of the circulation. The differences in the forces are of the same order of magnitude. Similar differences have already been observed in previous studies (Meyer Forsting et al. Reference Meyer Forsting, Pirrung and Ramos-García2019b). These differences are believed to be of the order of magnitude of errors due to the actuator line approximation. For example, in the code comparison of Martinez-Tossas et al. (Reference Martinez-Tossas, Churchfield, Yilmaz, Sarlak, Johnson, Sørensen, Meyers and Meneveau2018), larger differences were observed between different implementations of the actuator line. Hence, for practical purposes, such differences can be considered negligible.
The bookkeeping of position and circulation of the free wake (see § 3.4) begins at the first time step. But the application of the correction is turned on after 
$t=T/20$. The same orders of magnitude of differences were observed in figure 7, even for the points right after the correction is turned on (after the first two points).
It should be noted that the effects of the smearing parameter can never be eliminated completely. The smearing correction has the objective of correcting the forces at the blades. Possibly, a better understanding of the role of drag, viscous effects, unsteadiness or other phenomena might reduce even further the errors caused by the Gaussian smearing. However, some errors cannot be eliminated, because the discretization used in actuator line simulations does not usually support vortices with very low vortex core size. Hence the vorticity that is shed is still smeared. The smeared vortices create a wake that is different depending on the choice of 
$\varepsilon$, as can be seen in Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019b), especially in the near wake. As a consequence, differences in the wake can have an effect on the forces. Errors in the velocity have a first-order effect in the circulation and forces, hence if the modified near wake has a slightly different blockage profile, then the velocities at the turbine, and consequently the forces, would be affected.
8. Conclusions
The conception of the smearing correction by Dağ & Sørensen (Reference Dağ and Sørensen2020) (originally in Dağ Reference Dağ2017) and by Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) (originally in Martínez Tossas Reference Martínez Tossas2017), and the formal analysis provided by Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) and its improvement and analogy to a viscous lifting line by Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a, Reference Meyer Forsting, Pirrung and Ramos-García2020), were great advancements to the ALM. These steps provided not only a more accurate, general and reliable method but also an understanding of the mathematical and physical reasons for the errors of the uncorrected actuator line. However, these methods resort to iterative procedures with a relaxation factor, which are generally slower, less stable and less deterministic than direct methods. The previous method that did not employ an iterative method (Martínez-Tossas & Meneveau Reference Martínez-Tossas and Meneveau2019) had to waive the compatibility between circulation and induced velocity at each time step, which may cause errors in unsteady simulations.
In the present work, a non-iterative vortex-based smearing correction for the actuator line method is proposed and validated. Based on the linearization of a lifting line method, the iterative procedures of previous works have been substituted by the direct solution of a small linear system. For the cases tested, no significant difference is observed in the results of the iterative method and the non-iterative method. All differences are presumed to be orders of magnitude lower than the accuracy of the actuator line method.
The smearing correction reduces the effect of the smearing parameter 
$\varepsilon$ on the forces. However, differences, which are of the order of the errors of the actuator line approximation, were observed in the forces and circulation for different smearing parameters. This indicates that the smearing parameter still affects the simulation, which agrees with the observation that smearing has an effect on the near wake (Meyer Forsting et al. Reference Meyer Forsting, Pirrung and Ramos-García2019b). Nevertheless, for practical applications, the differences observed in the forces can be considered negligible.
Another contribution of the present work includes the use of a correction function based on the velocity induced by a smeared vortex segment, which eliminates the need for approximations of the correction function. Also, we implement a free-vortex wake model to define the vortex sheet based on the CFD velocities. This keeps the generality of the method and makes it applicable to several configurations without the need for ad hoc assumptions. If computational cost is prioritized, then one could use approximate functions, prescribe the wakes, or apply one of the methods for increasing the speed from Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2020). The focus of this study is on presenting and validating the technique. Many options to reduce the computational cost even further could be implemented. Additionally, the optimal choice of parameters for a desired accuracy is a matter of investigation. For these reasons, the analysis of the computational cost of different choices of options and parameters of the correction is left for future studies.
In order to further understand the limitations of the smearing correction, the comparison to blade-resolved simulations is suggested for future studies. Also, the implementation of three-dimensional corrections for drag and the understanding of unsteady effects are current areas of active research (see Kleine Reference Kleine2022).
Additionally, by carefully constructing a nonlinear lifting line method and an actuator line method that are consistent with each other, we arrive at differences in the induced velocity that are considerably better than differences reported in the literature. The good agreement serves as an a posteriori justification for the linearization of the smearing correction, in addition to the proof that the iterative lifting line method and the ALM with vortex-based smearing correction are mathematically identical if the circulation in each segment of the actuator line is assumed constant.
So far, the actuator line method has been used primarily by the rotor aerodynamics community, possibly because of its known limitations to replicate the forces on non-rotating wings. The use of the smearing correction with a general formulation removes this limitation and reproduces the results of a lifting line method, as shown here and in previous works. This could motivate other communities, such as the aeronautical community, to also adopt this technique to simulate wings, taking advantage of the lower grid requirements allowed by the ALM, in flow conditions more complex than those allowed by the lifting line method.
Acknowledgements
The authors appreciate the contributions of the anonymous reviewers, especially the reviewer who contributed to the development of the analytical formula for the velocity induced by a smeared vortex segment.
Funding
The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the High Performance Computing Center North (HPC2N). This work was conducted within StandUp for Wind. V.G.K. thanks KTH Engineering Mechanics for partially funding this work.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Linearized lifting line
Equation (2.10) can be linearized around the undisturbed velocities using the Taylor series expansion
\begin{equation} \varGamma_j = \varGamma_{0j} \!+\! \frac{1}{2}\,c_j \left[ \left( C_l(\alpha_{0j})\,\frac{U_y}{U_0} \!+\! \frac{\partial C_l}{\partial \alpha}(\alpha_{0j})\,\frac{U_z}{U_0} \right) u_{y}^{vi} \!+\! \left( C_l(\alpha_{0j})\,\frac{U_z}{U_0} - \frac{\partial C_l}{\partial \alpha}(\alpha_{0j})\, \frac{U_y}{U_0} \right) u_{z}^{vi} \right] , \end{equation}where
\begin{gather} U_{0} := \sqrt{(U_z^2+U_y^2)}, \end{gather}
\begin{gather} \alpha_{0j} := \alpha_{gj} + \arctan{\left(\frac{U_y}{U_z}\right)}, \end{gather}
\begin{gather} \varGamma_{0j} := \tfrac{1}{2}\,U_{0} c_j\,C_l(\alpha_{0j}) , \end{gather}
and 
$({\partial C_l}/{\partial \alpha })(\alpha _{0j})$ is the slope of the lift coefficient evaluated at angle 
$\alpha _{0j}$.
Defining
\begin{gather} b_{yj} := \frac{1}{2}\,c_j \left( C_l(\alpha_{0j})\,\frac{U_y}{U_0} + \frac{\partial C_l}{\partial \alpha}(\alpha_{0j})\,\frac{U_z}{U_0} \right), \end{gather}
\begin{gather} b_{zj} := \frac{1}{2}\,c_j \left( C_l(\alpha_{0j})\,\frac{U_z}{U_0} - \frac{\partial C_l}{\partial \alpha}(\alpha_{0j})\,\frac{U_y}{U_0} \right) ,\end{gather}we arrive at
\begin{equation} \varGamma_j = \varGamma_{0j} + b_{yj} u_{y}^{vi} + b_{zj} u_{z}^{vi} , \end{equation}
where we see that coefficients 
$b_{yj}$ and 
$b_{zj}$ are a measure of the sensitivity of the circulation to a change in the velocities in the 
$y$- and 
$z$-directions.
To avoid an iterative process, another relationship between the induced velocities and circulation must be used. This relationship is the one described in step (iii) of the iterative process of § 4. The velocities induced at control point 
$j$ by the vortex system 
$k$ (which includes the bound vortex and the vortex sheet created by the segment 
$k$) can be written as
\begin{gather} u_{y}^{vi_k}(\boldsymbol{x}_j) = a_{y,jk} \varGamma_k, \end{gather}
\begin{gather} u_{z}^{vi_k}(\boldsymbol{x}_j) = a_{z,jk} \varGamma_k, \end{gather}
where 
$a_{y,jk}$ and 
$a_{y,jk}$ are called the influence coefficients. The values of the influence coefficients are obtained from the Biot–Savart law and depend only on the shape of the vortex filament and its relative position to the control point. The formula for the influence coefficient for some common vortex filament shapes can be found in Katz & Plotkin (Reference Katz and Plotkin1991). Then the total induced velocities are found by considering all 
$N$ vortex systems created by the 
$N$ segments:
\begin{equation} \begin{bmatrix} u_{y}^{vi}(\boldsymbol{x}_1) \\ u_{y}^{vi}(\boldsymbol{x}_2) \\ \vdots \\ u_{y}^{vi}(\boldsymbol{x}_N) \\ \end{bmatrix} = \begin{bmatrix} a_{y,11} & a_{y,12} & \cdots & a_{y,1N} \\ a_{y,21} & a_{y,22} & \cdots & a_{y,2N} \\ \vdots & & \ddots & \vdots \\ a_{y,N1} & a_{y,N2} & \cdots & a_{y,NN} \\ \end{bmatrix} \begin{bmatrix} \varGamma_1 \\ \varGamma_2 \\ \vdots \\ \varGamma_N \end{bmatrix}, \end{equation}which can be written in matrix form as
\begin{equation} \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{y}}}^{vi} = \boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{y}}} \boldsymbol{\varGamma}. \end{equation}
Analogously for 
$u_{z}^{v}$:
\begin{equation} \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{z}}}^{vi} = \boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{z}}} \boldsymbol{\varGamma} . \end{equation}Equation (A7) can also be written in matrix form:
\begin{equation} \boldsymbol{\varGamma} = \boldsymbol{\varGamma}_0 + \boldsymbol{\mathsf{b}}_{\boldsymbol{\mathsf{y}}} \circ \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{y}}}^{vi} + \boldsymbol{\mathsf{b}}_{\boldsymbol{\mathsf{z}}} \circ \boldsymbol{\mathsf{u}}_{\boldsymbol{\mathsf{z}}}^{vi} = \boldsymbol{\varGamma}_0 + ({\rm diag}(\boldsymbol{\mathsf{b}}_{\boldsymbol{\mathsf{y}}})\, \boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{y}}}) \boldsymbol{\varGamma} + ({\rm diag}(\boldsymbol{\mathsf{b}}_{\boldsymbol{\mathsf{z}}})\, \boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{z}}}) \boldsymbol{\varGamma}, \end{equation}
where 
$\circ$ denotes the elementwise product, and 
${\rm diag}(\boldsymbol{\mathsf{b}})$ is the square diagonal matrix formed by the elements of column vector 
$\boldsymbol{\mathsf{b}}$. The circulation is finally found by solving the linear system
\begin{equation} (\boldsymbol{\mathsf{I}} - {\rm diag}(\boldsymbol{\mathsf{b}}_{\boldsymbol{\mathsf{y}}})\,\boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{y}}} - {\rm diag}(\boldsymbol{\mathsf{b}}_{\boldsymbol{\mathsf{z}}})\, \boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{z}}} ) \boldsymbol{\varGamma} = \boldsymbol{\varGamma}_0 , \end{equation}
where 
$\boldsymbol{\mathsf{I}}$ is the identity matrix. Therefore, (A14) provides a direct way to calculate the circulation along the wing, instead of the iteration of § 4. If needed, the induced velocities can be found using relations (A11) and (A12).
It is worth remembering that the current reference system is related to the local reference system of figure 1. Equation (A1), for example, depends on only the local variables. The only terms that take into account the contributions of other sections of the blades are the influence coefficients, hence these are defined first in an absolute system of reference and then transformed into the local system of reference. Also, there is no constraint regarding the values or direction of the velocity (as opposed to the classical linearization of the lifting line method, which requires constant 
$U_z$ and 
$U_z \gg U_y$). Therefore, this method is applicable directly to any system of reference, including the case of rotating blades, in both a fixed frame of reference and a rotating frame of reference.
This formulation is reduced to one of the classical implementations of the numerical lifting line for 
$U_y=0$ and constant 
$U_z=U_{\infty }$, when an ideal aerofoil and horseshoe vortices are employed. In this case, (A14) becomes
\begin{equation} - \frac{1}{2 {\rm \pi}} \left({\rm diag}(\boldsymbol{\mathsf{c}}/\boldsymbol{\mathsf{2}})\right)^{ - 1} \boldsymbol{\varGamma} + \boldsymbol{\mathsf{A}}_{\boldsymbol{\mathsf{y}}} \boldsymbol{\varGamma} = - U_{\infty} \boldsymbol{\alpha_{g}} , \end{equation}in which each term of (A15) can be compared directly to a term in (8.11) of Katz & Plotkin (Reference Katz and Plotkin1991):
\begin{equation} - \frac{\varGamma (x)}{2 {\rm \pi}\left(c(x)/2\right)} - \frac{1}{4 {\rm \pi}}\int_{x_{min}}^{x_{max}}\frac{-\partial \varGamma(x')/{\rm d}\kern 0.06em x}{x-x'}\,{{\rm d}\kern 0.06em x}' = - U_{\infty} \alpha_g, \end{equation}where the notation of Katz & Plotkin (Reference Katz and Plotkin1991) was modified to be consistent with our notation.
Appendix B. Resolution of the Gaussian vortex core
Usually, it is desirable to employ a low smearing parameter. The minimum smearing parameter is commonly taken as 
$\varepsilon \approx 2 \Delta x$ (Troldborg Reference Troldborg2009; Martínez-Tossas, Churchfield & Leonardi Reference Martínez-Tossas, Churchfield and Leonardi2015). However, for 
$\varepsilon = 2 \Delta x$, the Gaussian vortex core may not be well represented in the numerical grid. As mentioned in § 7.1, the low resolution of the vortex core may introduce errors, because the theory of § 3 is based on a perfect representation of the Gaussian distribution of vorticity in the vortex core.
In order to exemplify this error, the case of § 7.1 was simulated with the standard grid (
$\Delta x = R/56$ in the region of the actuator line) for a low smearing parameter 
$\varepsilon =2 \Delta x = R/28$. This case is compared to a simulation with the same smearing parameter, but a higher grid resolution, in figure 8.
Figure 8. Influence of grid resolution for the same value of the smearing parameter. (a) Difference between 
$y$-velocities obtained by the ALM and LL. (b) Difference between 
$z$-velocities obtained by the ALM and LL. (c) Relative difference between circulation obtained by the ALM and LL.
For 
$\varepsilon =2 \Delta x$, the difference in the induced velocity 
$u_y$ is of the order of 
$10^{-3}$ (
$0.1\,\%$ of the undisturbed velocity). It is still considerably better than the results found in the literature, and a generally acceptable error for an ALM method. However, it is a clear deterioration from the agreement found for larger smearing parameters in § 7.1. By increasing the grid resolution to 
$\Delta x=R/112$, so that 
$\varepsilon =R/28=4 \Delta x$, the differences are back to the order of magnitude discussed in § 7.1.
Observing the velocity 
$u_y$ and the vorticity in the region of the bound vortex in the plane of symmetry (
$x=0$), in figure 9, it is possible to note that the poorer grid does not have adequate resolution to represent the Gaussian vortex core. We conclude that this is the main cause of the high values of error in 
$u_y$ for the coarser grid. This error has already been identified by Shives & Crawford (Reference Shives and Crawford2013), where a smearing parameter 
$\varepsilon =4 \Delta x$ was recommended to reduce errors in the angle of attack, which are discussed further by Meyer Forsting & Troldborg (Reference Meyer Forsting and Troldborg2020). Meyer Forsting et al. (Reference Meyer Forsting, Pirrung and Ramos-García2019a) also employed a larger smearing parameter than the recommended minimum in order to avoid this error.
Figure 9. Comparison of the influence of grid resolution for the same smearing parameter (
$\varepsilon =R/28$). Pseudo-colour of 
$u_y$ and contours of total vorticity in the vortex core: (a) 
$\Delta x = R/56 = \varepsilon /2$, (b) 
$\Delta x = R/112 = \varepsilon /4$.
Hence for 
$\varepsilon =2 \Delta x$, it can be concluded that the error in 
$u_y$ is dependent mostly on the grid resolution, not the smearing parameter. The differences in 
$u_z$, however, have a component related to the discretization error, but also a relevant component related to the smearing parameter, especially near the tip. This is due to the differences in the formulation of the ALM and LL, as discussed in § 7.1.
Open access
