Filtered lifting line theory and application to the actuator line model

Abstract

Lifting line theory describes the cumulative effect of shed vorticity from finite span lifting surfaces. In this work, the theory is reformulated to improve the accuracy of the actuator line model (ALM). This model is a computational tool used to represent lifting surfaces, such as wind-turbine blades in computational fluid dynamics. In ALM, blade segments are represented by means of a Gaussian body force distribution with a prescribed kernel size. Prior analysis has shown that a representation of the blade using an optimal kernel width $\unicode[STIX]{x1D716}^{opt}$ of approximately one quarter of the chord size results in accurate predictions of the velocity field and loads along the blades. Also, simulations have shown that use of the optimal kernel size yields accurate representation of the tip-vortex size and the associated downwash resulting in accurate predictions of the tip losses. In this work, we address the issue of how to represent the effects of finite span wings and tip vortices when using Gaussian body forces with a kernel size larger than the optimal value. This question is relevant in the context of coarse-scale large-eddy simulations that cannot afford the fine resolutions required to resolve the optimal kernel size. For this purpose, we present a filtered lifting line theory for a Gaussian force distribution. Based on the streamwise component of the vorticity transport equation, we develop an analytical model for the induced velocity resulting from the spanwise changes in lift force for an arbitrary kernel scale. The results are used to derive a subfilter-scale velocity model that is used to correct the velocity along the blade when using kernel sizes larger than $\unicode[STIX]{x1D716}^{opt}$. Tests are performed in large-eddy simulation of flow over fixed wings with constant and elliptic chord distributions using various kernel sizes. Results show that by using the proposed subfilter velocity model, kernel-size independent predictions of lift coefficient and total lift forces agree with those obtained with the optimal kernel size.

1 Introduction

The actuator line model (ALM) has become a prominent tool to represent lifting surfaces, such as wind-turbine blades, without having to resolve the entire flow field near the blades (Sørensen & Shen 2002; Ivanell, Sørensen & Henningson 2007; Mikkelsen et al. 2007; Troldborg, Sørensen & Mikkelsen 2010; Porté-Agel et al. 2011; Churchfield et al. 2012; Jha et al. 2014; Martínez-Tossas, Churchfield & Leonardi 2015; van Kuik 2018). ALM involves distributing actuator points on a line along the span of a blade, along which the lift force is evaluated from tabulated or otherwise known lift and drag coefficients and the local incoming flow conditions. For implementation in a computational simulation code, such as in large-eddy simulations (LES), the resulting line force is filtered using a Gaussian kernel with characteristic length scale $\unicode[STIX]{x1D716}$ representing the spatial extent over which the force is spatially distributed. In recent work (Martínez-Tossas, Churchfield & Meneveau 2017), an optimal smoothing length scale $\unicode[STIX]{x1D716}^{opt}$ was determined, which minimizes the error between the actual velocity field and that modelled by ALM in two-dimensional (ideal) flow over a lifting airfoil. Specifically, when using $\unicode[STIX]{x1D716}^{opt}\approx 0.25c$ , where $c$ is the blade’s chord length, it was shown in Martínez-Tossas et al. (2017) that the flow field obtained by using the Gaussian body force provided good agreement with the potential flow solution for flow over a Joukowski airfoil. The latter was the flow solution taken to be the ground truth for the purpose of the comparison with the flow induced by the Gaussian force. Numerical tests have obtained similar values for the optimal smoothing length (Shives & Crawford 2013; Jha et al. 2014). Use of the optimal kernel scale was tested in a three-dimensional simulation of a wind turbine using the ALM (Martínez-Tossas, Churchfield & Meneveau 2016). Very good agreement was found between the simulation and blade element momentum (BEM) theory with a tip loss correction (Hansen 2007), especially near the tip. The accuracy was attributed to the capability of the simulation with $\unicode[STIX]{x1D716}^{opt}$ to generate shed vorticity and resulting tip vortices with a realistic core size.

As is known from classical lifting line theory (Prandtl & Tietjens 1957; Batchelor 2000; Anderson Jr 2010), the downwash associated with the tip vortex changes the local angle of attack along the blade, thus affecting the lift coefficient. In the ALM, the body force smooths the vorticity generation encountered by the incoming flow, resulting in a vortex with a finite core. The prior results (Martínez-Tossas et al. 2016) show that if the kernel size is chosen appropriately, good results may be obtained using the ALM without any additional tip loss correction, owing to the fact that the tip vortex core induced by the optimal Gaussian kernel size is realistic.

However, the optimal kernel size is rather small, as one quarter of the chord is often considered too expensive for practical applications of the ALM, and use of significantly larger kernel sizes is often unavoidable in practice. Prior attempts of using coarse kernel sizes (e.g. Stevens, Martínez-Tossas & Meneveau 2018) have shown good results in terms of the downstream flow (wakes in wind farm simulations), but the blade torque and power generation were generally too high, sometimes even exceeding the Betz limit. The reason is that the use of $\unicode[STIX]{x1D716}>\unicode[STIX]{x1D716}^{opt}$ leads to LES tip vortices that are too thick and too weak. Hence, the downwash is underpredicted and the lift is overpredicted, especially near the blade tip or near regions in the blade where the lift coefficient or chord (i.e. the circulation) varies significantly.

For LES using $\unicode[STIX]{x1D716}>\unicode[STIX]{x1D716}^{opt}$ , we therefore require a subfilter model for the local induced velocity that mimics the effect of a fully resolved tip vortex as well as any shed vorticity along the span. Developing such a model is the main purpose of this paper. In § 2, we briefly review elements of classic lifting line theory and show that the Gaussian filtering associated with the ALM leads to an unclosed subfilter modelling problem, even if the problem appears at first sight to be a linear one. Therefore, a lifting line theory must be derived for a Gaussian-filtered lift force distribution, which is accomplished in § 3. It is shown in § 4 that with certain simplification, a general semianalytical solution can be derived. Attention is placed in developing practically useful approximate calculation methods to facilitate implementation in LES, as described in § 5. In recent work, (Dag 2017; Dag & Sørensen 2018) researchers have proposed a similar approach based on empirical observations about velocity profiles in shed vortices with smoothed cores. As will be shown, the present theory serves as formal proof that their approach can be derived from first principles and as a starting point for detailed analysis. Sample large-eddy simulations of flow over three types of finite-length lifting surfaces are presented in § 6, testing the subfilter-scale velocity correction derived using the proposed theory. Conclusions are presented in § 7.

2 Gaussian filtering and linearity

We begin by considering the case of an infinitely long (two-dimensional (2-D)) lifting surface. In that case, the lift force (per unit fluid mass) as an actuator line may be thought of as distributed over a (singular) line, according to:

(2.1) $$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D70C}}f_{y}(x,y,z)=-\unicode[STIX]{x1D6E4}U_{\infty }\unicode[STIX]{x1D6FF}(x)\unicode[STIX]{x1D6FF}(y)=-\frac{1}{2}c_{L}cU_{\infty }^{2}\unicode[STIX]{x1D6FF}(x)\unicode[STIX]{x1D6FF}(y).\end{eqnarray}$$

Here, $\unicode[STIX]{x1D6E4}$ is the circulation, $U_{\infty }$ the inflow velocity, $c$ the chord length, $c_{L}$ the lift coefficient and $\unicode[STIX]{x1D6FF}$ are delta Dirac functions. Also, $x$ is the coordinate in the direction of the incoming flow, $y$ is the direction of the lift force and $z$ is the coordinate along the blade span. We consider the actuator line model using a (three-dimensional (3-D)) Gaussian kernel as proposed in the original introduction of the ALM (Sørensen & Shen 2002):

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D716}}(x,y,z)=\frac{1}{\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x03C0}^{3/2}}\text{e}^{-(x^{2}+y^{2}+z^{2})/\unicode[STIX]{x1D716}^{2}}.\end{eqnarray}$$

For the case of an infinitely long blade, integration along the spanwise direction (from $z=-\infty$ to $z=\infty$ ) leads to the 2-D force distribution:

(2.3) $$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D70C}}\widetilde{f}_{y}(\boldsymbol{x})=\int \frac{1}{\unicode[STIX]{x1D70C}}f_{y}(\boldsymbol{x}^{\prime })\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D716}}(\boldsymbol{x}-\boldsymbol{x}^{\prime })d^{3}\boldsymbol{x}^{\prime }=-\frac{1}{\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x03C0}}\frac{1}{2}c_{L}cU_{\infty }^{2}\text{e}^{-(x^{2}+y^{2})/\unicode[STIX]{x1D716}^{2}}.\end{eqnarray}$$

As derived in Martínez-Tossas et al. (2017), the curl of this force is a source in the transport equation for the $z$ -component vorticity. The linearized $z$ -component vorticity transport equation (linearized in the sense that the advection velocity in the $x$ -direction is set to a constant velocity $U_{\infty }$ everywhere) was integrated and the $z$ -component vorticity distribution derived. The result is a Gaussian distribution with width $\unicode[STIX]{x1D716}$ (see also Forsythe et al. 2015). Applying Biot–Savart then leads to the perturbation (tangential) velocity distribution:

(2.4) $$\begin{eqnarray}\widetilde{u}_{\unicode[STIX]{x1D703}}^{\prime }=\frac{1}{4\unicode[STIX]{x1D70B}}c_{L}U_{\infty }\frac{c}{r}(1-\text{e}^{-r^{2}/\unicode[STIX]{x1D716}^{2}}),\end{eqnarray}$$

where $r^{2}=x^{2}+y^{2}$ . This flow is a Lamb–Oseen vortex with core size $\unicode[STIX]{x1D716}$ representing the bound vorticity resulting from the Gaussian-filtered lifting surface.

It turns out that the same result is obtained if one simply applies the Gaussian filter to an ideal vortex velocity distribution. That is to say, recall that if one uses Biot–Savart on the initially unfiltered line vortex, one obtains the ideal bound vortex with tangential perturbation velocity $u_{\unicode[STIX]{x1D703}}^{\prime }=(1/4\unicode[STIX]{x1D70B})c_{L}U_{\infty }(c/r)$ . When this velocity is filtered (in two dimensions) using the kernel $\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D716}}(x,y)$ , the same Lamb–Oseen vortex with tangential velocity $\tilde{u} _{\unicode[STIX]{x1D703}}^{\prime }$ and core size $\unicode[STIX]{x1D716}$ as given in (2.4) is obtained. This equivalence occurs because the relationship between perturbation velocity and applied force, under the assumption of linearized perturbation vorticity transport, is linear. That is to say, we can write $\boldsymbol{u}^{\prime }={\mathcal{L}}(\,\boldsymbol{f})$ , where ${\mathcal{L}}$ represents linear operations (taking the curl of the force, integrating the linearized vorticity transport equation and applying Biot–Savart). Applying the Gaussian filter to this expression leads to $\widetilde{\boldsymbol{u}}^{\prime }=\widetilde{{\mathcal{L}}(\,\boldsymbol{f})}={\mathcal{L}}(\,\widetilde{\boldsymbol{f}})$ , allowing us to obtain the result of Martínez-Tossas et al. (2017) by simply filtering the ideal vortex tangential velocity distribution.

Next, we consider the case of a finite-length lifting surface represented as a singular line as in (2.1) but with $\unicode[STIX]{x1D6E4}(z)=0$ for $z<0$ . Figure 1(a) shows the coordinate system for which $z=0$ is at the tip of the blade and the blade extends toward $z>0$ .

Figure 1. (a) Sketch of finite lifting surface, where the origin of the coordinate system is placed at the blade tip, where $z=0$ . (b) Illustration of the Gaussian-filtered actuator line model with a force distribution across the blade’s cross-section and $z$ -dependent lift force (or circulation) distribution. Also shown (in green) are some ‘thick’ shed vortices and (in blue along the $z$ -axis) their induced velocity, $u_{y}^{\prime }$ . The distribution of the total induced $y$ -velocity as a function of the span is obtained by an integration over induced velocity from the shed vorticity.

In this case, classical lifting line theory states that because $\unicode[STIX]{x1D6E4}$ in (2.1) depends on  $z$ , vorticity is shed, generating an induced velocity at the blade that changes the angle of attack and affects the lift coefficient and thus $\unicode[STIX]{x1D6E4}$ itself. A self-consistency condition leads to the classical integro-differential equation for the circulation (Prandtl & Tietjens 1957) along the blade:

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}(z)=\frac{1}{2}\left(c_{Lb}+\frac{\text{d}c_{Lb}}{\text{d}\unicode[STIX]{x1D6FC}}\frac{u_{y}^{\prime }(z)}{U_{\infty }}\right)cU_{\infty }=\frac{1}{2}c_{Lb}cU_{\infty }+\frac{1}{2}\frac{\text{d}c_{Lb}}{\text{d}\unicode[STIX]{x1D6FC}}c\int _{0}^{L}\frac{\unicode[STIX]{x1D6E4}^{\prime }(z^{\prime })}{4\unicode[STIX]{x03C0}(z^{\prime }-z)}\,\text{d}z^{\prime },\end{eqnarray}$$

where $c_{Lb}$ is the lift coefficient at the baseline angle of attack (i.e. corresponding to the local free-stream velocity  $U_{\infty }$ ), $\text{d}c_{Lb}/\text{d}\unicode[STIX]{x1D6FC}$ is the slope of the lift versus angle of attack curve (typically close to $2\unicode[STIX]{x03C0}$ ) and $u_{y}^{\prime }(z)$ is the induced downwash velocity at position $z$ because of changes in circulation everywhere on the blade of span  $L$ . The typical approach solves this equation by expanding $\unicode[STIX]{x1D6E4}(z)$ into a Fourier series, and using the condition that $\unicode[STIX]{x1D6E4}(z)=0$ at the blade tips leads to the choice of a sine series. In general, the Fourier coefficients must be found numerically, even for the simple case of a constant chord and baseline lift parameters, $c_{Lb}$ and $\text{d}c_{Lb}/\text{d}\unicode[STIX]{x1D6FC}$ . Stewartson (1960) provides an analytical solution for the case of a semi-infinite constant chord wing.

At this stage, consider the case when the corresponding lift force distribution

(2.6) $$\begin{eqnarray}\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}f_{y}(x,y,z) & = & \displaystyle -\unicode[STIX]{x1D6E4}(z)U_{\infty }\unicode[STIX]{x1D6FF}(x)\unicode[STIX]{x1D6FF}(y)H(z)\nonumber\\ \displaystyle & = & \displaystyle -\frac{1}{2}\left(c_{Lb}+\frac{\text{d}c_{Lb}}{\text{d}\unicode[STIX]{x1D6FC}}\frac{u_{y}^{\prime }(x,y,z)}{U_{\infty }}\right)cU_{\infty }^{2}\unicode[STIX]{x1D6FF}(x)\unicode[STIX]{x1D6FF}(y)H(z),\end{eqnarray}$$

is filtered with the Gaussian kernel, $\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D716}}$ ( $H(z)$ is the Heaviside function because there is no force at $z<0$ even if $u_{y}^{\prime }$ is non-zero there). Even though in this case we can still write a linear equation, $\boldsymbol{u}^{\prime }={\mathcal{L}}(\,\boldsymbol{f})={\mathcal{L}}(\,\boldsymbol{f}_{b}+c_{1}\boldsymbol{u}^{\prime })$ , the spatially filtered velocity no longer solves the same equation as the original lifting line theory:

(2.7) $$\begin{eqnarray}\widetilde{\boldsymbol{u}}^{\prime }={\mathcal{L}}(\,\widetilde{\boldsymbol{f}})={\mathcal{L}}(\,\widetilde{\boldsymbol{f}}_{b}+\widetilde{c_{1}\boldsymbol{u}^{\prime }})\neq {\mathcal{L}}(\tilde{\boldsymbol{f}}_{b}+c_{1}\widetilde{\boldsymbol{u}}^{\prime }).\end{eqnarray}$$

The reason is that the prefactor, $c_{1}$ , in this case contains the Heaviside function, $H(z)$ , and clearly $\widetilde{H\boldsymbol{u}_{y}^{\prime }}\neq H\widetilde{\boldsymbol{u}^{\prime }}_{y}$ . Further, the velocity field induced by a force filtered at scale $\unicode[STIX]{x1D716}$ is given by:

(2.8) $$\begin{eqnarray}\boldsymbol{u}_{\unicode[STIX]{x1D716}}^{\prime }={\mathcal{L}}(\,\widetilde{\boldsymbol{f}})={\mathcal{L}}(\,\widetilde{\boldsymbol{f}}_{b}+\widetilde{c_{1}\boldsymbol{u}_{\unicode[STIX]{x1D716}}^{\prime }}).\end{eqnarray}$$

Clearly, ${\boldsymbol{u}^{\prime }}_{\unicode[STIX]{x1D716}}\neq \widetilde{\boldsymbol{u}^{\prime }}$ because $\widetilde{c_{1}\boldsymbol{u}^{\prime }}\neq \widetilde{c_{1}\widetilde{\boldsymbol{u}}^{\prime }}$ . As a result, we cannot simply use the solution of the classical lifting line theory, $\boldsymbol{u}^{\prime }$ , for an ideal line force and spatially filter the velocity to obtain the correct result corresponding to a filtered force distribution along a finite-length lifting blade. A full derivation is therefore needed to determine $\boldsymbol{u}_{\unicode[STIX]{x1D716}}^{\prime }$ as solution to (2.8), to which the next section is devoted.

3 Gaussian-filtered lifting line theory

For generality, we assume a span-dependent chord length, $c(z)$ , inflow velocity, $U_{\infty }(z)$ and lift coefficient, $c_{L}(z;\unicode[STIX]{x1D716})$ and we note that the effective lift coefficient depends upon the kernel size, $\unicode[STIX]{x1D716}(z)$ . We define $c_{Lb}(z)$ to be the (known) lift coefficient at the baseline angle of attack (i.e. corresponding to the local free-stream velocity, $U_{\infty }(z)$ ). For simplicity, we assume that the chord and kernel size do not vary rapidly along the span (i.e. at various stages we will assume $\text{d}c/\text{d}z\ll 1$ and $\text{d}\unicode[STIX]{x1D716}/\text{d}z\ll 1$ ). Later, applications with variable $\unicode[STIX]{x1D716}(z)$ and $c(z)$ will show that accurate results may also be obtained in such cases.

As in the prior section, we consider the unfiltered actuator line to be a lifting surface represented as a ‘delta function’ distribution in the $xy$ plane and a $z$ -dependent force (per unit length) along a line (i.e. the $z$ -dependent lift distribution) is given by:

(3.1) $$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D70C}}f_{y}(x,y,z)=-G(z)\unicode[STIX]{x1D6FF}(x)\unicode[STIX]{x1D6FF}(y),\end{eqnarray}$$

where

(3.2) $$\begin{eqnarray}G(z)=\unicode[STIX]{x1D6E4}(z)U_{\infty }(z)={\textstyle \frac{1}{2}}c_{L}(z)c(z)U_{\infty }^{2}(z)H(z).\end{eqnarray}$$

To explicitly take into account modifications of the local lift coefficient due to changes in incident velocity direction, we recall that the effective lift coefficient can be written as:

(3.3) $$\begin{eqnarray}c_{L}(z;\unicode[STIX]{x1D716})=c_{Lb}(z)+\frac{\text{d}c_{L}}{\text{d}\unicode[STIX]{x1D6FC}}(z)\frac{u_{y}^{\prime }(0,0,z;\unicode[STIX]{x1D716})}{U_{\infty }(z)},\end{eqnarray}$$

where $u_{y}^{\prime }(0,0,z;\unicode[STIX]{x1D716})$ is the local perturbation velocity (evaluated at $x=y=0$ ) due to downwash from shed vorticity. For an ideal lifting surface, one may use $\text{d}c_{L}/\text{d}\unicode[STIX]{x1D6FC}=2\unicode[STIX]{x03C0}$ but, in practice, different values may be used.

As described in the preceding section, in the actuator line model the actual force per unit volume seen by the fluid is given by the convolution of the line strength with the Gaussian kernel $\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D716}}$ according to:

(3.4) $$\begin{eqnarray}\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}\widetilde{f}_{y}(\boldsymbol{x}) & = & \displaystyle \int \frac{1}{\unicode[STIX]{x1D70C}}f_{y}(\boldsymbol{x}^{\prime })\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D716}}(\boldsymbol{x}-\boldsymbol{x}^{\prime })\,\text{d}^{3}\boldsymbol{x}^{\prime }\nonumber\\ \displaystyle & = & \displaystyle -\iiint \frac{1}{\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x03C0}^{3/2}}G(z^{\prime })\unicode[STIX]{x1D6FF}(x^{\prime })\unicode[STIX]{x1D6FF}(y^{\prime })\text{e}^{-[(x-x^{\prime })^{2}+(y-y^{\prime })^{2}+(z-z^{\prime })^{2}]/\unicode[STIX]{x1D716}^{2}}\,\text{d}x^{\prime }\,\text{d}y^{\prime }\,\text{d}z^{\prime }.\quad\end{eqnarray}$$

Evaluating the $xy$ integrals and changing variables to $z^{\prime \prime }=z^{\prime }-z$ , we obtain:

(3.5) $$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D70C}}\widetilde{f}_{y}(x,y,z)=\frac{1}{\unicode[STIX]{x1D70C}}f_{y}(x,y,z;\unicode[STIX]{x1D716})=-\frac{1}{\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x03C0}^{3/2}}\text{e}^{-(x^{2}+y^{2})/\unicode[STIX]{x1D716}^{2}}\int _{-\infty }^{\infty }G(z^{\prime \prime }+z)\text{e}^{-z^{\prime \prime 2}/\unicode[STIX]{x1D716}^{2}}\,\text{d}z^{\prime \prime },\end{eqnarray}$$

where from here on we omit tildes to denote filtered variables and indicate their dependence on $\unicode[STIX]{x1D716}$ instead.

3.1 Distribution of shed vorticity

Next, we consider the steady-state vorticity transport equation for the streamwise vorticity component in the presence of the Gaussian-filtered lift body force in ideal, steady flow (neglecting viscosity and turbulence), also neglecting changes in advective velocity due to the vorticity (linearization). The only term generating vorticity is the curl of the Gaussian body force, whose streamwise component can be expressed in the $x$ -component vorticity transport equation:

(3.6) $$\begin{eqnarray}U_{\infty }(z)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{x}}{\unicode[STIX]{x2202}x}=\frac{1}{\unicode[STIX]{x1D70C}}(\unicode[STIX]{x1D735}\times \boldsymbol{f})_{x}=-\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}f_{y}}{\unicode[STIX]{x2202}z}=\frac{\text{e}^{-(x^{2}+y^{2})/\unicode[STIX]{x1D716}^{2}}}{\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x03C0}^{3/2}}\int _{-\infty }^{\infty }\frac{\text{d}}{\text{d}z}G(z^{\prime \prime }+z)\text{e}^{-z^{\prime \prime 2}/\unicode[STIX]{x1D716}^{2}}\,\text{d}z^{\prime \prime }.\end{eqnarray}$$

Integrating in $x$ between $-\infty$ and $x$ , assuming that $\unicode[STIX]{x1D714}_{x}=0$ when $x\rightarrow -\infty$ and changing variables back to $z^{\prime }=z^{\prime \prime }+z$ , yields the following streamwise vorticity distribution:

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{x}(\boldsymbol{x};\unicode[STIX]{x1D716})=\frac{1}{2U_{\infty }(z)\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x03C0}}[\text{erf}(x/\unicode[STIX]{x1D716})+1]\text{e}^{-y^{2}/\unicode[STIX]{x1D716}^{2}}\int _{-\infty }^{\infty }\frac{\text{d}G(z^{\prime })}{\text{d}z^{\prime }}\text{e}^{-(z-z^{\prime })^{2}/\unicode[STIX]{x1D716}^{2}}\,\text{d}z^{\prime }.\end{eqnarray}$$

This $\unicode[STIX]{x1D714}_{x}(\boldsymbol{x};\unicode[STIX]{x1D716})$ distribution arises from the downstream transport of the fluid rotation induced by the transverse variations of lift induced by the Gaussian-filtered ALM force. Approaching the blade near the tip with the free stream, the vorticity rises from zero to its downstream value following the error function, it has a Gaussian profile in the cross-stream $y$ -direction, while the shape of the tip vortex in the spanwise direction is non-trivial and involves the integrated effects from the possibly varying lift coefficient. If the actual lift coefficient, chord and velocity were constant for $z>0$ and zero otherwise (i.e. $G(z)=(1/2)c_{L0}cU_{\infty 0}^{2}H(z)$ so that $G^{\prime }(z)=(1/2)c_{L0}cU_{\infty 0}^{2}\unicode[STIX]{x1D6FF}(z)$ ) and if the lift did not change with angle of attack (i.e. $\text{d}c_{L}/\text{d}\unicode[STIX]{x1D6FC}=0$ , which is unphysical), equation (3.7) shows that the tip vortex for $x/\unicode[STIX]{x1D716}\gg 1$ is a Lamb–Oseen vortex with Gaussian vorticity distribution in the $zy$ plane whose core size is precisely $\unicode[STIX]{x1D716}$ (similar to the bound vortex representing the blade circulation in the $xy$ plane analysed in Martínez-Tossas et al. 2017). However, the induced velocity changes the local lift coefficient and so we seek a more general solution that takes this effect into account.

3.2 Induced velocity perturbation

To find the induced vertical velocity along the blade associated with this 3-D vorticity distribution, we use the Biot–Savart equation:

(3.8) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle u_{y}^{\prime }(0,0,z;\unicode[STIX]{x1D716})=-\frac{1}{4\unicode[STIX]{x03C0}}\iiint \frac{\unicode[STIX]{x1D714}_{x}(\boldsymbol{x}^{\prime \prime };\unicode[STIX]{x1D716})(z-z^{\prime \prime })}{(x^{\prime \prime 2}+y^{\prime \prime 2}+(z-z^{\prime \prime })^{2})^{3/2}}\,\text{d}^{3}\boldsymbol{x}^{\prime \prime }=-\frac{1}{8\unicode[STIX]{x03C0}^{2}\unicode[STIX]{x1D716}^{2}}\int _{-\infty }^{\infty }\frac{\text{d}G(z^{\prime })}{\text{d}z^{\prime }}\\ \displaystyle \iiint \frac{1}{U_{\infty }(z^{\prime \prime })}\frac{(z-z^{\prime \prime })[\text{erf}(x^{\prime \prime }/\unicode[STIX]{x1D716})+1]}{[x^{\prime \prime 2}+y^{\prime \prime 2}+(z-z^{\prime \prime })^{2}]^{3/2}}\text{e}^{-(y^{\prime \prime 2}+(z^{\prime \prime }-z^{\prime })^{2})/\unicode[STIX]{x1D716}^{2}}\,\text{d}x^{\prime \prime }\,\text{d}y^{\prime \prime }\,\text{d}z^{\prime \prime }\,\text{d}z^{\prime }.\end{array}\right\}\end{eqnarray}$$

To proceed, it is useful to neglect the $z^{\prime \prime }$ variations of $U_{\infty }(z^{\prime \prime })$ in the denominator inside the integral and approximate $U_{\infty }(z^{\prime \prime })\approx U_{\infty }(z)$ . This can be justified because the kernel $(z-z^{\prime \prime })/|\boldsymbol{x}-\boldsymbol{x}^{\prime \prime }|^{3}$ peaks at $z^{\prime \prime }=z$ . Furthermore, we recognize that the integrand in $x^{\prime \prime }$ is even in the denominator while in the numerator the $\text{erf}(..)$ term is odd. Hence, we are left with:

(3.9) $$\begin{eqnarray}u_{y}^{\prime }(0,0,z;\unicode[STIX]{x1D716})=-\frac{1}{8\unicode[STIX]{x03C0}^{2}\unicode[STIX]{x1D716}^{2}U_{\infty }(z)}\int _{-\infty }^{\infty }\frac{\text{d}G(z^{\prime })}{\text{d}z^{\prime }}\iiint \frac{(z-z^{\prime \prime })\text{e}^{-(y^{\prime \prime 2}+(z^{\prime \prime }-z^{\prime })^{2})/\unicode[STIX]{x1D716}^{2}}}{[x^{\prime \prime 2}+y^{\prime \prime 2}+(z-z^{\prime \prime })^{2}]^{3/2}}\,\text{d}x^{\prime \prime }\,\text{d}y^{\prime \prime }\,\text{d}z^{\prime \prime }\,\text{d}z^{\prime },\end{eqnarray}$$

which, with an additional change of variables, $\unicode[STIX]{x1D701}=z^{\prime \prime }-z^{\prime }$ can be rewritten as:

(3.10) $$\begin{eqnarray}u_{y}^{\prime }(0,0,z;\unicode[STIX]{x1D716})=-\frac{1}{8\unicode[STIX]{x03C0}^{2}\unicode[STIX]{x1D716}U_{\infty }(z)}\int _{-\infty }^{\infty }\frac{\text{d}G(z^{\prime })}{\text{d}z^{\prime }}F(z-z^{\prime })\,\text{d}z^{\prime },\end{eqnarray}$$

where

(3.11) $$\begin{eqnarray}F(z)=\frac{1}{\unicode[STIX]{x1D716}}\iiint \frac{(z-\unicode[STIX]{x1D701})}{[x^{\prime \prime 2}+y^{\prime \prime 2}+(z-\unicode[STIX]{x1D701})^{2}]^{3/2}}\text{e}^{-(y^{\prime \prime 2}+\unicode[STIX]{x1D701}^{2})/\unicode[STIX]{x1D716}^{2}}\,\text{d}x^{\prime \prime }\,\text{d}y^{\prime \prime }\,\text{d}\unicode[STIX]{x1D701}.\end{eqnarray}$$

The integration yields (Martínez-Tossas 2017):

(3.12) $$\begin{eqnarray}F(z)=\frac{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}}{z}\left(1-\text{e}^{-(z/\unicode[STIX]{x1D716})^{2}}\right).\end{eqnarray}$$

Hence, we find that the velocity perturbation may be written as a superposition of shifted Lamb–Oseen vortices, weighted by the gradient of the lift function  $G(z)$ :

(3.13) $$\begin{eqnarray}u_{y}^{\prime }(z;\unicode[STIX]{x1D716})=-\frac{1}{U_{\infty }(z)}\int _{-\infty }^{\infty }\frac{\text{d}G(z^{\prime })}{\text{d}z^{\prime }}\frac{1}{4\unicode[STIX]{x03C0}(z-z^{\prime })}\left(1-\text{e}^{-(z-z^{\prime })^{2}/\unicode[STIX]{x1D716}^{2}}\right)\,\text{d}z^{\prime },\end{eqnarray}$$

where $u_{y}^{\prime }(0,0,z;\unicode[STIX]{x1D716})$ on the $x=y=0$ axis has been denoted as $u_{y}^{\prime }(z;\unicode[STIX]{x1D716})$ , also indicating that the solution depends on the kernel width, $\unicode[STIX]{x1D716}$ . This is the Gaussian force field analogue of lifting line theory in which the induced velocity from thin tip vortices may be expressed as an integral of the spanwise derivative of circulation. Because of the Gaussian smoothing, a Lamb–Oseen kernel arises in this expression, again with a core size exactly equal to $\unicode[STIX]{x1D716}$ used in the Gaussian applied body force that represents the lifting surface. Note that the same expression proposed by Dag (2017) and Dag & Sørensen (2018) was based on observations that smoothed tip vortices in ALM have the shape of Lamb–Oseen vortices. The result in (3.13) serves as formal proof that their empirically motivated approach can also be derived from first principles. It also enables further analysis and solution methods as detailed below.

Since the effective lift coefficient will vary due to the induced downwash, the strength of the vortex spanwise change, $\text{d}G(z^{\prime })/\text{d}z^{\prime }$ , will be affected by the downwash as a result of a modification of  $c_{L}(z)$ :

(3.14) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle G(z)=G_{b}(z)+\frac{1}{2}\frac{\text{d}c_{L}}{\text{d}\unicode[STIX]{x1D6FC}}(z)u_{y}^{\prime }(z;\unicode[STIX]{x1D716})U_{\infty }(z)c(z),\quad \text{where}\\ \displaystyle G_{b}(z)=\frac{1}{2}c_{Lb}(z)c(z)U_{\infty }^{2}(z)\end{array}\right\}\end{eqnarray}$$

corresponds to the ‘baseline’ case (i.e. the lift coefficient appropriate for the undisturbed inflow angle of attack) without downwash. Replacing $G(z)$ in (3.13), we obtain:

(3.15) $$\begin{eqnarray}\displaystyle u_{y}^{\prime }(z;\unicode[STIX]{x1D716})U_{\infty }(z) & = & \displaystyle \int _{-\infty }^{\infty }\frac{\text{d}G_{b}(z^{\prime })}{\text{d}z^{\prime }}g(z-z^{\prime };\unicode[STIX]{x1D716})\,\text{d}z^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\int _{-\infty }^{\infty }\frac{\text{d}}{\text{d}z^{\prime }}\left[\frac{1}{2}\frac{\text{d}c_{L}}{\text{d}\unicode[STIX]{x1D6FC}}(z^{\prime })u_{y}^{\prime }(z^{\prime };\unicode[STIX]{x1D716})U_{\infty }(z^{\prime })c(z^{\prime })\right]g(z-z^{\prime };\unicode[STIX]{x1D716})\,\text{d}z^{\prime },\qquad\end{eqnarray}$$

where

(3.16) $$\begin{eqnarray}g(z,\unicode[STIX]{x1D716})=-\frac{1}{4\unicode[STIX]{x03C0}z}(1-\text{e}^{-z^{2}/\unicode[STIX]{x1D716}^{2}}).\end{eqnarray}$$

This expression is an integro-differential equation that also appears in the classical lifting line theory. In prior applications (without the Gaussian filtering), such an equation (usually formulated for the circulation) is solved for each case using a sine series expansion. Here, we propose to first recast this expression as a proper Fredholm integral equation of the second kind, using integration by parts for the second term:

(3.17) $$\begin{eqnarray}\displaystyle u_{y}^{\prime }(z;\unicode[STIX]{x1D716})U_{\infty }(z) & = & \displaystyle \int _{-\infty }^{\infty }\frac{\text{d}G_{b}(z^{\prime })}{\text{d}z^{\prime }}g(z-z^{\prime },\unicode[STIX]{x1D716})\,\text{d}z^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\int _{-\infty }^{\infty }\unicode[STIX]{x1D703}(z^{\prime })u_{y}^{\prime }(z^{\prime };\unicode[STIX]{x1D716})U_{\infty }(z^{\prime })c(z^{\prime })k(z-z^{\prime };\unicode[STIX]{x1D716})\,\text{d}z^{\prime },\end{eqnarray}$$

where

(3.18) $$\begin{eqnarray}k(z;\unicode[STIX]{x1D716})=\unicode[STIX]{x03C0}\frac{\text{d}}{\text{d}z}\left[-\frac{1}{4\unicode[STIX]{x03C0}z}(1-\text{e}^{-z^{2}/\unicode[STIX]{x1D716}^{2}})\right]=\frac{1}{4z^{2}}(1-\text{e}^{-z^{2}/\unicode[STIX]{x1D716}^{2}})-\frac{1}{2\unicode[STIX]{x1D716}^{2}}\text{e}^{-z^{2}/\unicode[STIX]{x1D716}^{2}}.\end{eqnarray}$$

Note that we have also introduced a variable that describes possible deviations from the idealized lift curve slope:

(3.19) $$\begin{eqnarray}\unicode[STIX]{x1D703}(z)=\frac{1}{2\unicode[STIX]{x03C0}}\frac{\text{d}c_{L}}{\text{d}\unicode[STIX]{x1D6FC}}(z).\end{eqnarray}$$

The equation is cast most naturally in terms of the unknown perturbation momentum flux:

(3.20) $$\begin{eqnarray}M(z;\unicode[STIX]{x1D716})=u_{y}^{\prime }(z;\unicode[STIX]{x1D716})U_{\infty }(z)c(z).\end{eqnarray}$$

Multiplying (3.17) by $c(z)$ yields:

(3.21) $$\begin{eqnarray}M(z;\unicode[STIX]{x1D716})=\int _{-\infty }^{\infty }\frac{\text{d}G_{b}(z^{\prime })}{\text{d}z^{\prime }}c(z)g(z-z^{\prime };\unicode[STIX]{x1D716})\,\text{d}z^{\prime }+\int _{-\infty }^{\infty }\unicode[STIX]{x1D703}(z^{\prime })M(z^{\prime };\unicode[STIX]{x1D716})c(z)k(z-z^{\prime };\unicode[STIX]{x1D716})\,\text{d}z^{\prime }.\end{eqnarray}$$

For general distributions of $c(z)$ , $\unicode[STIX]{x1D716}(z)$ and $\unicode[STIX]{x1D703}(z)$ along the blade, one must solve this integral equation numerically for each individual case.

4 Simplifications, solutions and approximate analytical fit

Next, we propose an approximate method that will be shown to be of sufficient accuracy for many applications and helpful as a reference solution of an interesting special case. Specifically, we consider a lifting surface extending to lengths much larger than the chord length. The blade tip is located at, say, $z=0$ . Also, we approximate $\unicode[STIX]{x1D703}(z)$ by its average distribution over the blade, $\overline{\unicode[STIX]{x1D703}}$ (i.e. $\unicode[STIX]{x1D703}(z)\approx \overline{\unicode[STIX]{x1D703}}H(z)$ ). Under these conditions, we may rewrite (3.21) as:

(4.1) $$\begin{eqnarray}M(z;\unicode[STIX]{x1D716})=\int _{-\infty }^{\infty }\frac{\text{d}G_{b}(z^{\prime })}{\text{d}z^{\prime }}c(z)g(z-z^{\prime };\unicode[STIX]{x1D716})\,\text{d}z^{\prime }+\overline{\unicode[STIX]{x1D703}}\int _{0}^{\infty }M(z^{\prime };\unicode[STIX]{x1D716})c(z)k(z-z^{\prime };\unicode[STIX]{x1D716})\,\text{d}z^{\prime }.\end{eqnarray}$$

A further simplification is proposed; namely, we consider the special case when the ratio $\unicode[STIX]{x1D716}(z)/c(z)$ is held constant. We thus define the $z$ -independent dimensionless parameter:

(4.2) $$\begin{eqnarray}\unicode[STIX]{x1D700}\equiv \frac{\unicode[STIX]{x1D716}(z)}{c(z)\overline{\unicode[STIX]{x1D703}}}.\end{eqnarray}$$

As mentioned before, the optimal kernel size requires $\unicode[STIX]{x1D716}^{opt}(z)\approx 0.25c(z)$ (i.e. $\unicode[STIX]{x1D700}^{opt}\overline{\unicode[STIX]{x1D703}}\approx 0.25$ (constant)). For $\unicode[STIX]{x1D716}^{les}$ values typically used in coarse-scale LES, $\unicode[STIX]{x1D716}_{les}\gg c(z)$ , in which case the solution to the integral equation is simple. As will be shown later, in the limit $\unicode[STIX]{x1D716}_{les}\gg c(z)$ , the second term on the right-hand side of (4.1) can be neglected and the induced velocity is proportional to $g(z;\unicode[STIX]{x1D716})$ .

We now proceed assuming a constant $\unicode[STIX]{x1D716}(z)/c(z)$ . As shown in appendix A, the solution to the integral equation may be written in the form:

(4.3) $$\begin{eqnarray}M(z;\unicode[STIX]{x1D716})=\frac{1}{\unicode[STIX]{x1D700}}\int _{-\infty }^{\infty }\frac{\text{d}G_{b}(z^{\prime \prime })}{\text{d}z^{\prime \prime }}S(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime \prime },\unicode[STIX]{x1D700})\,\text{d}z^{\prime \prime },\quad \text{with}~\unicode[STIX]{x1D709}=\frac{z}{\unicode[STIX]{x1D716}(z)},\,\unicode[STIX]{x1D709}^{\prime \prime }=\frac{z^{\prime \prime }}{\unicode[STIX]{x1D716}(z)},\end{eqnarray}$$

and where $S(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime \prime },\unicode[STIX]{x1D700})$ is the solution to what we term the canonical integral equation of the Gaussian-filtered lifting line theory:

(4.4) $$\begin{eqnarray}S(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime \prime },\unicode[STIX]{x1D700})=g_{\ast }(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime \prime })+\frac{1}{\unicode[STIX]{x1D700}}\int _{0}^{\infty }S(\unicode[STIX]{x1D709}^{\prime },\unicode[STIX]{x1D709}^{\prime \prime };\unicode[STIX]{x1D700})k_{\ast }(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime })\,\text{d}\unicode[STIX]{x1D709}^{\prime }.\end{eqnarray}$$

The basic kernels are given by:

(4.5a,b ) $$\begin{eqnarray}g_{\ast }(\unicode[STIX]{x1D709})=-\frac{1}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D709}}\left(1-\text{e}^{-\unicode[STIX]{x1D709}^{2}}\right)\quad \text{and}\quad k_{\ast }(\unicode[STIX]{x1D709})=\frac{1}{4\unicode[STIX]{x1D709}^{2}}\left(1-\text{e}^{-\unicode[STIX]{x1D709}^{2}}\right)-\frac{1}{2}\text{e}^{-\unicode[STIX]{x1D709}^{2}}.\end{eqnarray}$$

Recalling the definition of $M(z;\unicode[STIX]{x1D716})$ , and with (4.3), the perturbation velocity ratio is given by:

(4.6) $$\begin{eqnarray}\frac{u_{y}^{\prime }(z;\unicode[STIX]{x1D716})}{U_{\infty }(z)}=\frac{1}{\unicode[STIX]{x1D716}(z)U_{\infty }^{2}(z)}\int _{-\infty }^{\infty }\frac{\text{d}G_{b}(z^{\prime \prime })}{\text{d}z^{\prime \prime }}(z^{\prime \prime })S\left(\frac{z}{\unicode[STIX]{x1D716}(z)},\frac{z^{\prime \prime }}{\unicode[STIX]{x1D716}(z)};\unicode[STIX]{x1D700}\right)\text{d}z^{\prime \prime }.\end{eqnarray}$$

We note from (4.4) that when $\unicode[STIX]{x1D700}\gg 1$ , the solution is simply $S(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime \prime };\unicode[STIX]{x1D700})=g_{\ast }(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime \prime })$ . For finite  $\unicode[STIX]{x1D700}$ , it is possible to solve (4.4) numerically, for various $\unicode[STIX]{x1D700}$ and $\unicode[STIX]{x1D709}^{\prime \prime }$ as parameters. An iterative algorithm is used with under-relaxation. Figure 2 shows the solutions for various $\unicode[STIX]{x1D700}$ and $\unicode[STIX]{x1D709}^{\prime \prime }$ values. As shown, the solutions for $\unicode[STIX]{x1D709}^{\prime \prime }=0$ have the expected downwash perturbation velocity but include a ‘regularized’ behaviour near the tip ( $z=0$ ), where the perturbation velocity drops to zero at the ‘core’ of the tip vortex in which the thickness depends upon  $\unicode[STIX]{x1D716}$ . For $\unicode[STIX]{x1D709}^{\prime \prime }>0$ , we begin to see the other side of the shed vortex with a positive perturbation velocity but somewhat non-trivial shape toward the blade tip. As $\unicode[STIX]{x1D709}^{\prime \prime }$ grows, at distances large (compared to chord) from the tip and for large $\unicode[STIX]{x1D716}/c$ , the solution approaches a function proportional to the function  $g$ .

Figure 2. Solutions to the canonical filtered lifting line (4.4), for various $\unicode[STIX]{x1D700}$ and $\unicode[STIX]{x1D709}^{\prime \prime }$ values; namely, $\unicode[STIX]{x1D709}^{\prime \prime }=0$  (a),  $\unicode[STIX]{x1D709}^{\prime \prime }=1$  (b),  $\unicode[STIX]{x1D709}^{\prime \prime }=5$  (c) and $\unicode[STIX]{x1D709}^{\prime \prime }=10$  (d). Solid lines: $S(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime \prime };\unicode[STIX]{x1D700})/\unicode[STIX]{x1D700}$ obtained from numerical solutions to (4.4). In each panel, the various curves shown are for $\unicode[STIX]{x1D700}=0.25$ (black), 0.5 (green), 0.75 (blue), 1.0 (red), 2.0 (cyan) and 5.0 (magenta). The dot-dashed lines correspond to the proposed empirical fit, $S^{fit}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime \prime };\unicode[STIX]{x1D700})/\unicode[STIX]{x1D700}$ , given by (4.7) (same colours as the lines for corresponding  $\unicode[STIX]{x1D700}$ ).

4.1 Empirical fits

For practical applications, one may wish to avoid having to solve numerically the canonical lifting line equation. Therefore, it is useful to provide a functional form for the solution $S(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime \prime };\unicode[STIX]{x1D700})$ , even if approximate. After some experimentation, we have found that the following expression represents the numerical solution with very good accuracy for $\unicode[STIX]{x1D700}\geqslant 0.25$ :

(4.7) $$\begin{eqnarray}S^{fit}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime \prime },\unicode[STIX]{x1D700})=-\text{sgn}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime \prime })[1-0.25\text{e}^{-\unicode[STIX]{x1D700}}(1-\text{e}^{-0.2|\unicode[STIX]{x1D709}^{\prime \prime }|})]\,f(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime \prime };\unicode[STIX]{x1D700}),\end{eqnarray}$$

where

(4.8) $$\begin{eqnarray}f(\unicode[STIX]{x1D709};\unicode[STIX]{x1D700})=\left[\frac{1}{4\unicode[STIX]{x03C0}|\unicode[STIX]{x1D709}|}\left(1-\text{e}^{-\unicode[STIX]{x1D709}^{2}}\right)-\frac{0.029\unicode[STIX]{x1D700}^{-2/3}}{\unicode[STIX]{x1D709}^{2}}\left(1-\text{e}^{-0.357|\unicode[STIX]{x1D709}|^{3}}\right)\right].\end{eqnarray}$$

The dot-dashed lines in figure 2 show the proposed fitted functional form, $S^{fit}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime \prime };\unicode[STIX]{x1D700})$ , compared to the accurate full numerical solution (solid lines). As shown, for $\unicode[STIX]{x1D700}\geqslant 0.25$ , the agreement is quite good with a root-mean-square error between the fit and the numerical integration over the range of parameter values ( $0.25\leqslant \unicode[STIX]{x1D700}\leqslant 5$ , $0\leqslant \unicode[STIX]{x1D709}^{\prime \prime }\leqslant 10$ , $0\leqslant \unicode[STIX]{x1D709}\leqslant 16$ ) of 0.0023 (i.e. less than 2 % of the peak value of $|S|$ shown). The maximum difference found over this range is 0.029, which can still be considered small for practical applications. For smaller values of $\unicode[STIX]{x1D700}$ ( $\unicode[STIX]{x1D700}<0.25$ ), the comparison (not shown) begins to deteriorate and (4.7) should not be used.

In figure 2, the solutions have been plotted as a function of the variable, $\unicode[STIX]{x1D709}$ (i.e. $z$ scaled with the filter width  $\unicode[STIX]{x1D716}$ ). To provide additional insights, the results can be plotted as a function of $z/c$ , for different  $\unicode[STIX]{x1D716}$ , as presented in figure 3. This also enables us to compare with classical lifting line theory, which is the expected behaviour of the filtered lifting line solution when $\unicode[STIX]{x1D716}\rightarrow 0$ . Note that this could not be displayed as a function of $\unicode[STIX]{x1D709}$ because the limit $\unicode[STIX]{x1D716}\rightarrow 0$ is singular for our filtered line theory. The result from the classical lifting line theory for a semi-infinite wing of constant chord is shown as a dashed black line for which we have used the solution provided by Stewartson (1960), summarized in appendix B. As expected, the solution to the filtered lifting line theory tends to the classical result in the limit $\unicode[STIX]{x1D700}\rightarrow 0$ .

Figure 3. Solutions to the canonical filtered lifting line equation for $\unicode[STIX]{x1D709}^{\prime \prime }=0$ plotted as a function of $z/c=\unicode[STIX]{x1D709}\unicode[STIX]{x1D700}$ . Solid lines: $S(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime \prime }=0;\unicode[STIX]{x1D700})/\unicode[STIX]{x1D700}$ obtained from numerical solutions to (4.4) for various $\unicode[STIX]{x1D700}$ values. The curves shown are for $\unicode[STIX]{x1D700}=0.10$ (blue), $\unicode[STIX]{x1D700}=0.25$ (black), $\unicode[STIX]{x1D700}=0.5$ (green), 0.75 (blue), 1.0 (red), 2.0 (cyan) and 5.0 (magenta). The dashed line is the ideal lifting line solution (Stewartson 1960) for a semi-infinite singular line (see appendix B). The dot-dashed lines correspond to the proposed empirical fit given by (4.7) for the $\unicode[STIX]{x1D700}\geqslant 0.25$ cases.

For practical implementations, a discrete version of (4.6) can be written on $N$  discrete element points $z_{i}=(i-1)L/(N-1)$ , covering a span of length, $L$ and where $z_{1}$ and $z_{N}$ are the two ends of the blade. However, recall that the canonical case considered, as well as the approximated fit, were derived for a semi-infinite blade. We thus approximate the real case as a superposition of the solutions of two semi-infinite blades. The summation is divided into two parts: from one tip to the middle of the span and another covering the other half of the span to the other tip. That is to say, we approximate the contributions from half of the span as coming from a semi-infinite blade with one tip at $z=z_{1}=0$ and similarly the other side starting from the other side with the tip at $z=L=z_{N}$ .

Thus, we write

(4.9) $$\begin{eqnarray}\displaystyle \frac{u_{y}^{\prime }(z_{i};\unicode[STIX]{x1D716})}{U_{\infty }(z_{i})} & = & \displaystyle \frac{1}{\unicode[STIX]{x1D716}(z_{i})U_{\infty }^{2}(z_{i})}\left[\mathop{\sum }_{j=1}^{N/2}\unicode[STIX]{x0394}G_{b}(z_{j})S\left(\frac{z_{i}}{\unicode[STIX]{x1D716}(z_{i})},\frac{z_{j}}{\unicode[STIX]{x1D716}(z_{i})};\unicode[STIX]{x1D700}\right)\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\mathop{\sum }_{j=N/2+1}^{N}\unicode[STIX]{x0394}G_{b}(z_{j})S\left(\frac{L-z_{i}}{\unicode[STIX]{x1D716}(z_{i})},\frac{L-z_{j}}{\unicode[STIX]{x1D716}(z_{i})};\unicode[STIX]{x1D700}\right)\right],\end{eqnarray}$$

where $\unicode[STIX]{x0394}G_{b}(z_{j})=(1/2)[G_{b}(z_{j+1})-G_{b}(z_{j-1})]$ , for $j=2,3,\ldots ,N-1,$ and where it is understood that $\unicode[STIX]{x0394}G_{b}(z_{1})=G_{b}(z_{1})$ because $G_{b}(z_{0})=0$ as $z_{0}$ falls outside the blade. At the other end, for  $z_{N}$ , the same holds but now $\unicode[STIX]{x0394}G_{b}(z_{N})=-G_{b}(z_{N})$ because $G_{b}(z_{N+1})=0$ as $z_{N+1}$ falls outside the blade.

4.2 Numerical tests and examples

To test the proposed approximation of (4.7) and to illustrate some sample solutions to the filtered lifting line in (3.21), we compare the predicted velocity perturbation as obtained by numerical integration of the Fredholm integral in (3.21) with the approximate solution using (4.9) and the fitted function. Different cases to be tested can be characterized by various baseline lift coefficient distributions, $c_{Lb}(z)$ , blade chord distributions, $c(z)$ , and free-stream velocity, $U_{\infty }(z)$ .

The first case tested is for a constant chord wing of length, $L$ , with chord $c(z)=c_{0}$ , $L/c_{0}=20$ and a constant lift coefficient $c_{Lb}(z)=1$ and constant lift slope ( $\unicode[STIX]{x1D703}(z)=1$ ) along the span. $U_{\infty }$ is constant as well. Three cases are compared: $\unicode[STIX]{x1D716}/c_{0}=\unicode[STIX]{x1D700}=0.25$ , 0.5 and 1.5. In all cases, numerical integration of the integral equation is done on a sufficiently fine grid for fully converged results. The approximated discrete sum in (3.21) involves $N=101$ points over the same length, $L$ (i.e. $\unicode[STIX]{x0394}z/c_{0}=0.2$ is used). Results are shown in figure 4. As shown, the agreement is very good. This same case will be considered later using an implementation of a constant chord, finite airfoil in LES.

Figure 4. Comparisons of the full numerical solution of filtered lifting line theory (lines) (3.21) applied to a constant chord, $c$ , blade of constant lift coefficient, $c_{Lb}=1$ and finite length $L/c=20$ , with discrete superposition (4.9) of the fitted solution to the canonical filtered lifting line solution (symbols). $U_{\infty }$ is constant. Black line and squares: $\unicode[STIX]{x1D700}=0.25$ , blue line and circles: $\unicode[STIX]{x1D700}=0.5$ . Red line and triangles: $\unicode[STIX]{x1D700}=1.5$ .

The second case tested is again for a constant chord, $c(z)=c_{0}$ , but with three sections with different lift coefficients along the span. Specifically, $c_{Lb}=1$ for $0\leqslant z/c_{0}<2.5$ , $c_{Lb}=1.5$ for $2.5\leqslant z/c_{0}<5$ and $c_{Lb}=0.5$ for $z/c_{0}\geqslant 5$ up to $z=L=20c_{0}$ . $U_{\infty }$ is constant as well and the case $\unicode[STIX]{x1D716}/c_{0}=\unicode[STIX]{x1D700}=0.5$ is used. Similar numerical implementations as in the prior example are used. Results are shown in figure 5. As shown, the agreement between the full numerical solution and the approximation is very good, even for cases of variable  $c_{Lb}$ .

Figure 5. Comparison of numerical solution of filtered lifting line theory (line) (3.21) applied to a constant chord, $c=c_{0}$ , wing of variable lift coefficient ( $c_{Lb}=1$ for $0<z/c_{0}<2.5$ , $c_{Lb}=1.5$ for $2.5<z/c_{0}<5$ and $c_{Lb}=0.5$ for $z/c_{0}>5$ ) and finite length, $L/c_{0}=20$ , with discrete superposition (4.9) of the fitted solution to the canonical filtered lifting line solution (symbols), using $\unicode[STIX]{x1D700}=0.75$ . $U_{\infty }$ is constant.

The last test case is a wing with a chord distribution reminiscent of a wind-turbine blade, with three distinct linear segments as shown with the bottom dot-dashed line in figure 6. The tip is at $z=0$ and the root of the blade at $z=L$ . We use a varying sectional lift coefficient along the span, with $c_{Lb}=1.5$ between $0\leqslant z/c_{0}<16$ , $c_{Lb}=0.5$ between $16\leqslant z/c_{0}<18$ and $c_{Lb}=0$ between $18\leqslant z/c_{0}\leqslant 20$ , mimicking the non-lifting, typically cylindrical portion near the nacelle. Also, in this case, we vary $U_{\infty }(z)$ as if it was the relative velocity of a blade operating at a tip-speed ratio of $U_{\infty \text{-}tip}/U_{axial}=5$ , namely:

(4.10) $$\begin{eqnarray}U_{\infty }(z)=U_{\infty \text{-}tip}\,\sqrt{\left[5\left(1-\frac{z}{L}\right)\right]^{2}+1^{2}}.\end{eqnarray}$$

In this case, we use a constant kernel size; namely, $\unicode[STIX]{x1D716}/c_{0}=0.25$ (black line and squares in figure 6) and $\unicode[STIX]{x1D716}/c_{0}=2.0$ (blue line and circles). The comparison in figure 6 shows that the approximate method gives very good results, except for small $\unicode[STIX]{x1D716}$ in the near-hub region, where the accuracy of the approximate method is less critical because $c_{Lb}$ is very small and the impact on actual lift forces would be small.

Figure 6. Numerical solution (lines) of filtered lifting line theory (3.21) applied to a case with a wind-turbine-blade-like chord distribution, compared to the discrete superposition (4.9) of the fitted solution (symbols). Black line and squares, $\unicode[STIX]{x1D716}/c_{0}=0.25$ ; blue line and circles, $\unicode[STIX]{x1D716}=2.0$ . Note that in this case, $\unicode[STIX]{x1D700}=\unicode[STIX]{x1D716}/c(z)$ varies as a function of  $z$ , causing possible differences with the approximate solution that assumes that $\unicode[STIX]{x1D700}$ is constant (and also that $\text{d}\unicode[STIX]{x1D716}/\text{d}z\ll 1$ ). The dot-dashed line at the bottom (right $y$ -axis) displays the chord distribution, $c(z)/c_{0}$ , where $c_{0}=L/20$ is a reference scale and, in this case, $c_{0}=c_{tip}$ .

A test case for an elliptical wing yielded less precise agreement between the theory and numerical integrations because for such distributions the simplifying assumptions used in this section (e.g. small variations of $\unicode[STIX]{x1D716}$ along  $z$ ) are violated, especially near the tip. However, as shown in § 5, a fully explicit numerical approach applicable in LES also yields excellent results for such cases.

5 Velocity correction as a subfilter model for LES

In the ALM, we sample the fluid velocity, $\widetilde{\boldsymbol{u}}(\boldsymbol{x}_{i})$ , at the $i$ th actuator point, $\boldsymbol{x}_{i}=(x_{i},y_{i},z_{i})$ ; this velocity is used to evaluate the local ALM force, typically using tabulated lift and drag coefficients. Both the magnitude and direction of the local velocity matter. As has been argued in Martínez-Tossas et al. (2017), the bound vorticity resulting from the application of ALM forces does not affect the sampled velocity at the actuator point itself, for an infinite span wing, and when the effects of drag on the local velocity may be neglected. For a finite wing span, however, the local velocity is affected by the induced velocity, $u_{y}(z_{i};\unicode[STIX]{x1D716}^{LES})$ (henceforth we omit the prime in denoting the perturbation velocity). Thus, an approximation to the unperturbed incoming velocity, $\boldsymbol{U}_{\infty }(\boldsymbol{x}_{i})$ (excluding the effects of the velocity perturbation caused by shed vorticity) is given by:

(5.1) $$\begin{eqnarray}\boldsymbol{U}_{\infty }(\boldsymbol{x}_{i})=\widetilde{\boldsymbol{u}}(\boldsymbol{x}_{i})-u_{y}(z_{i};\unicode[STIX]{x1D716}^{LES})\boldsymbol{j},\end{eqnarray}$$

where it is understood that $\boldsymbol{i}$ , $\boldsymbol{k}$ and $\boldsymbol{j}$ are unit vectors in the directions of $\boldsymbol{U}_{\infty }(\boldsymbol{x}_{i})$ , the blade span and the direction perpendicular to both (in the lift direction), respectively.

The real blade, however, of size typically much smaller than the filter/kernel scale, $\unicode[STIX]{x1D716}^{LES}$ , will induce more concentrated shed vortices and induced velocity corrections, especially near regions of rapidly changing circulation (e.g. near the tip). As has been determined in Martínez-Tossas et al. (2017), a good representation of the flow for a real airfoil can be obtained using an optimal filter/kernel size, $\unicode[STIX]{x1D716}^{opt}\approx 0.25c$ . Therefore, we seek the velocity (denoted as $\widehat{\boldsymbol{u}}(\boldsymbol{x}_{i})$ ) that would have occurred had we used, in the LES, a kernel scale, $\unicode[STIX]{x1D716}^{opt}$ . This velocity can be obtained by adding to $\boldsymbol{U}_{\infty }(\boldsymbol{x}_{i})$ the velocity perturbation corresponding to  $\unicode[STIX]{x1D716}^{opt}$ :

(5.2) $$\begin{eqnarray}\widehat{\boldsymbol{u}}(\boldsymbol{x}_{i})=\widetilde{\boldsymbol{u}}(\boldsymbol{x}_{i})+[u_{y}(z_{i};\unicode[STIX]{x1D716}^{opt})-u_{y}(z_{i};\unicode[STIX]{x1D716}^{LES})]\,\boldsymbol{j}.\end{eqnarray}$$

Thus, the difference $\unicode[STIX]{x0394}u_{y}=u_{y}(z_{i};\unicode[STIX]{x1D716}^{opt})-u_{y}(z_{i};\unicode[STIX]{x1D716}^{LES})$ , can be considered as the subfilter velocity model to be added to the LES-sampled velocity. The corrected velocity, $\widehat{\boldsymbol{u}}(\boldsymbol{x}_{i})$ (its magnitude and, importantly, its direction with modified angle of attack) can then be used in the determination of lift and drag coefficients as well as the direction of the applied ALM forces. The lift is applied perpendicular to $\widehat{\boldsymbol{u}}(\boldsymbol{x}_{i})$ whereas the drag force is parallel.

5.1 Numerical implementation of subfilter scale velocity model

Here, we describe the procedure for implementing the subgrid-scale velocity correction in the actuator line model. There are two approaches to determining the velocity correction, $\unicode[STIX]{x0394}u_{y}$ , from the velocity perturbations evaluated at the two (LES and optimal) $\unicode[STIX]{x1D716}$ values. A first approach is to use (4.9) with the empirical fit function for the solution to the canonical filtered lifting line equation. A second approach is to implement a direct numerical solution to the original expression in (3.13) using a fully explicit approach to solve the integral equation via successive substitution and under-relaxation. In practical applications, we prefer this approach because it has similar numerical costs and requires less assumptions (e.g. $\text{d}\unicode[STIX]{x1D716}(z)/\text{d}z\ll 1$ is not required). The method consists of the following steps, where we assume the time step is denoted by superscripts $n$ (present) and $n-1$ (prior past time step):

1. (i) To compute $\unicode[STIX]{x0394}u_{y}^{n}(z_{i})$ , the (new) velocity correction at the ALM point, $z_{i}$ , determine the distribution of $G(z)=(1/2)c_{L}(z)c(z)U_{\infty }^{2}(z)$ at all other $N$ ALM points on the blade (i.e. evaluate $G^{n-1}(z_{j})$ ) from the previous time step according to:

   (5.3) $$\begin{eqnarray}G^{n-1}(z_{j})={\textstyle \frac{1}{2}}c(z_{j})U_{\infty }^{2}(z_{j})c_{L}(z_{j}).\end{eqnarray}$$

2. (ii) To evaluate numerically the integral of (3.13), compute the finite difference of $G^{n-1}(z)$ :

   (5.4) $$\begin{eqnarray}\unicode[STIX]{x0394}G^{n-1}(z_{j})={\textstyle \frac{1}{2}}[G^{n-1}(z_{j+1})-G^{n-1}(z_{j-1})],\quad \text{for}~j=2,3,\ldots ,N-1,\end{eqnarray}$$
   where the first and last points are:
   (5.5a,b ) $$\begin{eqnarray}\unicode[STIX]{x0394}G^{n-1}(z_{1})=G^{n-1}(z_{1})\quad \text{and}\quad \unicode[STIX]{x0394}G^{n-1}(z_{N})=-G^{n-1}(z_{N}).\end{eqnarray}$$

3. (iii) Compute numerically the new (not yet relaxed) induced perturbation velocities (superscript *) for both $\unicode[STIX]{x1D716}_{i}=\unicode[STIX]{x1D716}_{i}^{opt}$ and $\unicode[STIX]{x1D716}_{i}=\unicode[STIX]{x1D716}_{i}^{LES}$ :

   (5.6) $$\begin{eqnarray}\displaystyle u_{y}^{n\ast }(z_{i};\unicode[STIX]{x1D716}_{i}) & = & \displaystyle -\frac{1}{U_{\infty }(z_{i})}\int _{-\infty }^{\infty }\frac{\text{d}G^{n-1}(z^{\prime })}{\text{d}z^{\prime }}\frac{1}{4\unicode[STIX]{x03C0}(z_{i}-z^{\prime })}\left(1-\text{e}^{-(z_{i}-z^{\prime })^{2}/\unicode[STIX]{x1D716}_{i}^{2}}\right)\,\text{d}z^{\prime }\end{eqnarray}$$
   (5.7) $$\begin{eqnarray}\displaystyle & {\approx} & \displaystyle -\frac{1}{U_{\infty }(z_{i})}\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x0394}G^{n-1}(z_{j})\frac{1}{4\unicode[STIX]{x03C0}(z_{i}-z_{j})}\left(1-\text{e}^{-(z_{i}-z_{j})^{2}/\unicode[STIX]{x1D716}_{i}^{2}}\right).\end{eqnarray}$$
   Note that the perturbation velocities are perpendicular to the velocity $\boldsymbol{U}_{\infty }(\boldsymbol{x}_{i})$ .

4. (iv) Compute the updated difference between induced velocities for $\unicode[STIX]{x1D716}_{i}^{LES}$ and $\unicode[STIX]{x1D716}_{i}^{opt}$ and apply under-relaxation:

   (5.8) $$\begin{eqnarray}\unicode[STIX]{x0394}u_{y}^{n}(z_{i})=f_{u}\left[u_{y}^{n\ast }\left(z_{i};\unicode[STIX]{x1D716}_{i}^{opt}\right)-u_{y}^{n\ast }\left(z_{i};\unicode[STIX]{x1D716}_{i}^{LES}\right)\right]+(1-f_{u})\left[u_{y}^{n-1}(z_{i};\unicode[STIX]{x1D716}_{i}^{opt})-u_{y}^{n-1}(z_{i};\unicode[STIX]{x1D716}_{i}^{LES})\right].\end{eqnarray}$$
   The under-relaxation factor, $f_{u}\leqslant 1$ , is used to ensure numerical convergence of this successive substitution algorithm that can be implemented in an explicit LES code. A factor, $f_{u}=0.1$ , has been found to lead to rapidly converging and stable results and is recommended.

5. (v) When sampling the LES velocity at the actuator point $i$ for time  $n$ , add the velocity component subgrid-scale velocity to the sampled velocity:

   (5.9) $$\begin{eqnarray}\widehat{\boldsymbol{u}}(\boldsymbol{x}_{i})=\widetilde{\boldsymbol{u}}(\boldsymbol{x}_{i})+\unicode[STIX]{x0394}u_{y}^{n}(z_{i})\boldsymbol{j}.\end{eqnarray}$$

Now, with this new sampled velocity, the actuator line model can compute all of the forces and apply them to the flow field accordingly. The correction is computed at every time step. (Note that the operation count for the velocity correction is of order $N^{2}$ because convolutions are performed directly. If $N$ is large, fast-Fourier-transform-based approaches could be developed, but in our applications, $N$ has been sufficiently small so that further optimizations were not found to be necessary.)

6 Applications in LES

We now test the proposed subfilter-scale velocity model in large-eddy simulations of flow over several finite wings with uniform, non-turbulent inflow. Simulations are performed using LESGO, an LES code (Calaf, Meneveau & Meyers 2010; Graham & Meneveau 2012; Martínez-Tossas et al. 2018) that uses pseudo-spectral discretization in the spanwise and streamwise directions and second-order finite difference in the vertical direction. A fringe region is used to drive the flow to desired inflow conditions (uniform in the present application) once it reaches the end (or inflow due to periodicity) of the domain (Stevens, Graham & Meneveau 2014). Time integration is done using a second-order Adams–Bashforth scheme. The code uses a scale-dependent dynamic eddy-viscosity subgrid-scale model with Lagrangian averaging (Bou-Zeid, Meneveau & Parlange 2005), although in the present application of uniform inflow with emphasis on the aerodynamic forces induced on the blades, the effects of the subgrid model are not crucial. An eddy-viscosity model is included only to avoid numerical instabilities.

The simulations use a range of kernel sizes, from $\unicode[STIX]{x1D716}^{LES}/c=0.25$ to $\unicode[STIX]{x1D716}^{LES}/c=4$ . The length of the blades is denoted by $S$ (as opposed to $L$ used in prior sections to avoid confusion with the domain lengths). The simulation domain is $L_{x}/S=2.5$ in streamwise length, $L_{y}/S=2.5$ in vertical length and $L_{z}/S=2.5$ along the span. Figure 7(a) shows a sketch of the simulation domain. To mitigate any possible numerical effects, a very fine grid resolution is employed for all simulations, namely $\unicode[STIX]{x1D6E5}/c=0.0195$ , equal in all directions. The range of resolutions with respect to $\unicode[STIX]{x1D716}$ is from $\unicode[STIX]{x1D716}/\unicode[STIX]{x1D6E5}\approx 13$ for the case with $\unicode[STIX]{x1D716}/c=0.25$ , to $\unicode[STIX]{x1D716}/\unicode[STIX]{x1D6E5}\approx 103$ for the $\unicode[STIX]{x1D716}/c=2$ case. The length of the blade is $S/c=12.5$ . The blade has a NACA64A17 section (the same one used in the outer part of the NREL 5-MW reference turbine blade (Jonkman et al. 2009)), placed at a $6^{\circ }$ angle of attack so that the baseline lift coefficient is $C_{Lb}=1.103$ . The drag coefficient is set to zero for the test cases presented.

Figure 7. (a) Sketch of the computational domain used for a simulation of a finite-length constant chord wing. (b) Volume rendering of non-dimensional vorticity magnitude for a simulation of a stationary constant chord wing using the actuator line model.

The simulations were run with very high resolutions to avoid any effects of numerical grid resolution. However, we have verified that even much coarser resolutions ( $\unicode[STIX]{x1D716}/\unicode[STIX]{x0394}x\approx 2.8$ ) provide essentially the same results.

When simulating a fixed wing, tip vortices are formed on both sides of the blades. Figure 7(b) shows a volume rendering of the resulting vorticity magnitude. The vorticity field is reasonably uniform along the $z$ -direction along the blade until it reaches the tip. Although some vorticity is shed along the way, it is difficult to visualize in such a plot. At the tip, the vorticity takes a $90^{\circ }$ turn and travels downstream. Next, the downwash created by the shed vorticity is compared to the analytical solutions found for different cases.

First, we examine the induced velocity from LES along the blade caused by the shed vorticity in figure 8. The LES results at the various kernel sizes are shown as symbols (note that we also include the $\unicode[STIX]{x1D716}^{LES}/c=4$ case). The analytical expression ((4.9), with the function $S(..)$ evaluated using the fit of (4.7)) is plotted as solid lines. There is excellent agreement between the theory and the LES results.

Figure 8. Induced velocity along the blade for large-eddy simulations of a fixed wing under uniform inflow and angle of attack of  $6^{\circ }$ . Note that $r$ denotes the distance from the centre of the blade of length  $S$ . Symbols denote the LES results, whereas the solid lines are the theory and fit from (4.6) and (4.7).

Next, we compare the perturbation velocity, $u_{y}$ , obtained from LES before and after applying the subfilter-scale velocity correction. The proposed correction adds a subgrid-scale velocity along the actuator line perpendicular to the inflow velocity sampled at every actuator point. The desired behaviour is that the perturbation velocity after applying the correction should be independent of $\unicode[STIX]{x1D716}^{LES}$ . As shown in figure 9(a), the desired behaviour is indeed followed very accurately in the LES. Concomitant to the $\unicode[STIX]{x1D716}^{LES}$ -independent velocity perturbation, the lift coefficient obtained from the corrected velocity (figure 9 b) is also $\unicode[STIX]{x1D716}^{LES}$ -independent, whereas without the correction (dashed lines), strong dependence on $\unicode[STIX]{x1D716}^{LES}$ is observed.

Figure 9. Simulation results from a constant chord fixed wing under uniform inflow at an angle of attack of  $6^{\circ }$ . Induced velocity (a) and lift coefficient (b) along the blade for cases with and without the correction using different  $\unicode[STIX]{x1D716}$ .

Figure 10. Simulation results from an elliptic wing under uniform inflow at a nominal angle of attack of  $6^{\circ }$ . Induced velocity (a) and angle of attack (b) along the blade for cases without (dashed lines) and with (solid lines) the correction for different $\unicode[STIX]{x1D716}/c$ .

Next, a similar application and analysis are performed on an elliptic wing under uniform inflow, also at an angle of attack of  $6^{\circ }$ . The chord distribution along the blade span is given by:

(6.1) $$\begin{eqnarray}c(r)=c_{0}\sqrt{1-(2r/S)^{2}},\end{eqnarray}$$

where $c_{0}$ is the largest chord in the middle of the blade, $r$ is the coordinate along the span centred at 0 and $S$ is the total blade span. Results for the corrected and uncorrected velocities and angles of attack are shown in figure 10. Again, we observe that the curves for all the cases with the proposed correction collapse. Also note that the correction predicts a constant induced velocity throughout most of the blade, as expected for an elliptic wing. We note that near the tips, the deviations from a constant induced velocity are also as expected because of differences between the Gaussian-filtered and ideal lifting line theories. Only in the limit of $\unicode[STIX]{x1D716}\rightarrow 0$ would one expect constant induced velocity all the way to the (singular) tip. The remarkable result is that the subfilter model enables us to correct the induced velocities and angle of attack shown by the dashed lines to the $\unicode[STIX]{x1D716}$ -independent solid lines that all coincide across the span of the blade.

We can evaluate integrated quantities such as the total lift and induced drag force on the blades as a function of $\unicode[STIX]{x1D716}^{LES}$ . Figure 11 shows the total lift and induced drag forces on the blade from LES of the constant and elliptic chord wings, both with and without applying the subfilter-scale velocity correction. We see that the correction eliminates the $\unicode[STIX]{x1D716}$ -dependence of predicted aerodynamic forces. Without the correction, these important integral quantities could differ by as much as 5–10 %.

Figure 11. Total lift (a) and induced drag (b) for simulation of a constant chord wing. Squares show results without using the subfilter-scale velocity correction; circles use the correction. Forces are non-dimensionalized by the blade span, $S$ , width, $\ell$ , density, $\unicode[STIX]{x1D70C}$ , and free-stream velocity, $U_{\infty }$ .

Figure 12. Total lift (a) and induced drag (b) for simulation of a constant chord wing (top) and one with elliptic chord distribution (bottom). Squares show results without using the subfilter-scale velocity correction; circles use the correction.

7 Conclusions

A new formulation for a filtered lifting line theory has been proposed. The main difference from Prandtl’s original lifting line theory is that there is now a length scale $\unicode[STIX]{x1D716}$ for each blade section that establishes the width over which vorticity is distributed, as opposed to an infinitesimally small vorticity source. It has been found in prior work that this thickness is related to the blade chord, and a value of $\unicode[STIX]{x1D716}/c\approx 0.25$ provides optimal results. In the limit of $\unicode[STIX]{x1D716}$ going to zero, the filtered formulation yields the original lifting line theory.

The filtered lifting line theory can be used to formulate a subfilter-scale velocity correction. This correction is intended for simulations that cannot afford to resolve the recommended resolution associated with the optimal length scale, $\unicode[STIX]{x1D716}/c\approx 0.25$ . The correction provides the subgrid-scale velocities one would obtain by running a simulation with the optimal $\unicode[STIX]{x1D716}$ while using much coarser resolutions. It is worthwhile to recall that the filtered lifting line theory developed here relies on a number of crucial assumptions and simplifications, notably the neglect of viscous and turbulence effects. As justification, we note that the lift force on the wing depends mostly upon the pressure distribution, which is known to be relatively insensitive to viscous and turbulence effects. Even when operating in close-to-stall conditions, where turbulence in the wake may become substantial, we expect the model to yield accurate results, if one is using realistic tabulated lift and drag coefficients appropriate for local flow conditions (e.g. near stall).

The correction has been tested in large-eddy simulations of flow over a finite wing with constant chord and an elliptic chord distribution. There is excellent agreement between using different $\unicode[STIX]{x1D716}$ when invoking the proposed subfilter-scale velocity correction.

Further work should focus on applications of the subfilter velocity model on wind-turbine simulations, including blade rotation. Initial tests, to be presented elsewhere, yield promising results. Moreover, tests on coarse-resolution LES of entire wind farms with turbulent inflow conditions (Calaf et al. 2010; Stevens et al. 2018) are required to explore the performance of the model in the context of more practical applications.