A Solution Method for the Filtered Lifting Line Theory

The work of Sørensen and Shen [1] introduced the actuator line model as a way to model wind turbine blades in computational fluid dynamics using Gaussian body forces [1]. This work opened the door for a new way of modeling wings in computational fluid dynamics without the need to resolve the full geometry of the blades. The filtered lifting line theory formalized the governing equations in the actuator line model and established its connection to the original lifting line theory as a filtered formulation of the flow over a lifting line [2–5].

The filtered lifting line theory provides a continuous set of equations to model the flow over a three-dimensional wing using Gaussian body forces. The idea behind the theory equations is to model the blade forcing on a wing using a continuous representation of the forces, as opposed to point forces that lead to flow singularities. The formulation in filtered lifting line theory is also consistent with the definition of an actuator line model [1,6,7]. These formulations are typically used in computational fluid dynamic simulations to represent wings as body forces without the need to resolve the full geometry of the blades.

Traditional lifting line theory, such as the one presented by Glauert [8], represents the bound and wake vorticity using a collection of singular vorticity distributions. The method, suitable for wings with a large aspect ratio and small angles of attack, requires careful treatment of the integral involved, but leads to analytical solutions in particular cases, such as the elliptical wing [8]. The extension of the method to general circulation distributions along the wing requires using numerical methods, which traditionally rely on a discrete distribution of singular vortices [9–12]. Yet, the method appears to lack convergence as the discretization along the blade span is changed or the resolution increased. This issue was often neglected or not addressed in early literature that only focused on using few discrete elements. Practitioners often resort to using ad hoc methods to choose the distribution of control points and vortex points (i.e., inset of the tip panels [13], a cosine spacing distribution [14], or approximations based on gradients of the spanwise spacing [15]). Vorticity-based methods, such as the vortex particle method, are known to only converge (as the grid is refined) when a proper regularization (also referred to as mollification) is used [16]. The work of Chattot [11] introduces an artificial viscosity in the lifting line formulation to alleviate the numerical instabilities from singular point forces. Most of the methods focus on solving the lifting line equations without modifying the underlying equations that possess point force singularities.

The filtered lifting line theory introduces the regularization in the lifting line formulation using Gaussian body forces, based on the argument that vorticity cannot be singular and it can be better represented by a continuous, nonsingular distribution of body forces. The formulation results in a nonsingular and continuous representation of the problem and is able to converge as the grid is refined. The solutions depend on the width of the Gaussian body force, ϵ, and previous work has suggested that values on the order of a quarter of the chord provide optimal solutions [17,18]. There have been several adaptations of the filtered lifting line theory for wind turbine aerodynamics and fixed wing applications [19–22]. The theory can also be used as a subfilter correction for simulations with coarse resolutions that cannot afford to resolve the optimal kernel widths [21–25]. The filtered lifting line theory equations have also been used successfully in fixed wing applications [3,20,26]. Despite the theory being successfully applied in large eddy simulations (LES), we do not yet have a robust method to solve the stand-alone equations.

In this work, we reformulate the filtered lifting line theory equations as a function of the local flow angle at every blade element along the wing (following the approach of Ref. [27]) to obtain a fast, simple, and robust solution method. The approach from Ning [27] formulates the two-dimensional (2D; normal and tangential) blade-element momentum theory equations into a one-dimensional problem using the flow angle as the main variable. This approach allows for a simplified set of equations that is guaranteed to converge [27]. We adopt the same logic using an integral equation for the flow angle at every blade section. We propose a solution method for the equations based on a root-finding algorithm, establish resolution guidelines, and compare the analytical solutions to results from LES.

where $ϕ$ is the flow angle with respect to the inflow direction and β is the twist angle (relative to the inflow, $U∞$⁠).

where $x′, y′$⁠, and $z′$ are the locations relative to the center of each 2D airfoil section. The value of ϵ represents the physical length scale over which the forces are distributed. The main difference between the filtered and standard lifting line theory formulations is the width of the vorticity, ϵ. In the filtered lifting line theory, this value is proportional to the chord and somewhere in the vicinity of $ϵ/c=0.1−0.25$ based on the argument that the velocity cannot be infinite in the core of the vortex elements and should scale with the length scales of the problem (i.e., chord) [3,18,28]. The standard lifting line theory takes the limit of $ϵ→0$ [5]. Equation (1) is derived from first principles using body forces [2,3]. We note that a similar formulation can also be obtained by replacing the singular vortices in the standard lifting line formulation using Lamb–Oseen vortices [19,22].

Equation (1) can be cast as a Fredholm integral equation of the second kind by linearizing the circulation and assuming that the solution is within the linear region of the lift curve [3]. The Fredholm integral form of Eq. (1) can then be solved numerically using an iterative algorithm [3,29]. However, the iterative algorithm is slow (⁠$O$ min/h), requires under-relaxation, and does not always converge [3]. We seek to write the equations as a root-finding problem in terms of only one independent variable. This approach has been used successfully to reformulate the blade element momentum theory equations [27].

We note that Eq. (4) is an integral equation defined at every spanwise location and can be solved using a multidimensional root-finding algorithm.

The algorithm discussed next allows us to compute $F(ϕ)$ for a given $ϕ$⁠, which is then passed on to a root-finding algorithm. The goal is to find $ϕ$ that satisfies Eq. (6).

To solve Eq. (6), which is an implicit equation in terms of the variable $ϕ$⁠, we use a root-finding algorithm. The following algorithm outlines the process to evaluate $F(ϕ)$⁠.

This algorithm can be used within a root-finding algorithm that will find $ϕ(z)$ that matches Eq. (6) within a certain tolerance. The flow angle, $ϕ$⁠, is initialized to zero everywhere along the blade (⁠$α=β$⁠). We use the trapezoidal rule for numerical integration, although different numerical integration schemes could also be used. For the present work, we used the root-finding algorithms “Anderson,” “krylov,” and “df-sane” that are already implemented within the SciPy library to obtain the numerical solutions [30–33]. These algorithms provided the fastest solutions (⁠$O$ s) and arrived at consistent solutions. Other root-finding algorithms can also be used. The equations are valid for any relation between angle of attack and lift coefficient, and are not restricted to the linear portion of the lift curve, as done in the original lifting line theory and other solution methods [3,5].

To solve the filtered lifting line equations (using either the solution method presented here or an actuator line implementation in computational fluid dynamics) we need to establish the number of points required along the wing span. The results presented here use constant spacing (⁠$Δz$⁠) along the span of the wing. We chose the resolution based on the ratio of ϵ to the spanwise resolution, $Δz$⁠. For cases where the chord and ϵ vary along the blade, a nonuniform distribution of points can be employed to optimize resolution. We simulate the case from Ref. [3], a fixed wing of constant chord, NACA64A17 airfoil, twist angle of $6 deg$⁠, and blade span to chord ratio of $S/c=12.5$⁠. We solve Eq. (6) numerically and obtain the quantities along the blade for a range of $ϵ/c=0.1−5.0$⁠, following the typical range of values used in actuator line simulations [18,34].

Figure 3 shows the induced velocities along the blade obtained from the numerical solution. The value of $ϵ/c$ can have a significant effect on the induced velocities, thus the resulting blade loads, with the differences in induced velocities on the order of $0.1u′y/U∞$⁠. We performed a sensitivity study in Sec. 4.2 and determined that using (⁠$ϵ/Δz=10$⁠) provided converged results within less than 0.1% of the asymptotic solution. The results in Fig. 3 are taken to be grid-independent by using grid points along the blade such that $ϵ/Δz=10$⁠. The resolution along the blade can affect the accuracy of the solution, so we study its effect on the quantities along the blade in Sec. 4.1. In the actuator line model, there is a grid resolution of the computational fluid dynamics solver, and the resolution of the actuator line, denoting the spacing between actuator points (⁠$Δr$⁠). In the context of this paper and the filtered lifting line theory, there is only one resolution along the spanwise direction that represents the grid spacing of the computational fluid dynamics solver in the actuator line model that is necessary to resolve the flow dynamics. The spacing between actuator points is expected to be enough to provide a converged body force field (⁠$ϵ/Δr>1$⁠).

Next, we studied the impact of resolution on the induced velocities and compared it to the asymptotic solution for the same fixed wing case using different $ϵ/c$ values within the recommended range [6,7,34]. We describe the grid resolution in terms of points per width of the Gaussian kernel (⁠$ϵ/Δz$⁠). Figure 4 shows the resulting induced velocity from the numerical solution to the root-finding problem using a different number of grid points. As the grid is refined (larger $ϵ/Δz$⁠), the differences in induced velocities keep converging toward an asymptotic solution. The asymptotic solution is only a function of ϵ [2,3]. The results suggest that at least $ϵ/Δz>1$ should be used to obtain reasonable solutions, although higher values are needed for full convergence.

To better understand the effect of resolution on blade loading predictions, we focus on the integral quantity of total lift coefficient, $CL=L1/2 ρ S c U∞2$⁠. Figure 5 shows the total lift coefficient (C_L) as a function of resolution (⁠$ϵ/Δz$⁠) for different values of $ϵ/c$⁠. The range of $ϵ/c$ is within 0.15–1.0, to denote the range of optimal values typically used/recommended in the literature [17,18]. Actuator line simulations of wind turbine blades often use values much larger than the range provided here, and they can converge with fewer points [6,7]. However, solutions with too coarse values of ϵ would overpredict lift because the induced velocities are smeared too much [3,7]. Each line in Fig. 5 corresponds to simulations with fixed $ϵ/c$ and each marker is from a solution with different number of grid points along the span (⁠$ϵ/Δz$⁠). The solution converges as the grid is refined for all cases, as expected.

The results from this section can be used to provide guidelines and recommended values to obtain results within a certain tolerance of the converged solutions. For grid resolutions of $ϵ/Δz>2$⁠, the differences between solutions and the asymptotic value of total wing lift coefficient are less than 0.5% of the converged value. For differences less than 0.1% of the total wing lift coefficient, values above $ϵ/Δz>4$ should be used. Previous studies of the actuator line model have suggested to use values where $ϵ/Δz>5$ to obtain grid-independent results, and the results in Fig. 5 are consistent with these findings [7,17].

Next, we used the metric of maximum error in lift force along the span of the blade to evaluate the convergence of the method. This metric ensures that the predicted lift at every point along the blade is within a certain tolerance of the converged solution. Table 1 provides the spanwise resolution, $ϵ/Δz$⁠, needed to obtain blade loading within 5% and 1% of the converged solution. The converged solution is defined as the solution with a very fine resolution of $ϵ/Δz=30$⁠. We can generalize these values and say that for results within a few percent of the converged solution, $ϵ/Δz=1$ should suffice. However, to be within 1% of the converged solution, a finer resolution should be used up to $ϵ/Δz=3$⁠, depending on the value of $ϵ/c$⁠. These requirements can be used as guidelines to obtain converged results in solutions to the filtered lifting line theory and can also be applicable to the resolution requirements of the actuator line model. When using a nonuniform distributions of points, these guidelines can be adapted without compromising the computational requirements of the actuator line model or the solution method proposed in this paper. Future work should focus on using a nonuniform distributions of points.

These wings (Fig. 6) are meant to be illustrative examples to highlight the range of wings that can be simulated. All wings use the tabulated lift and drag coefficient data from a NACA64A17 airfoil, twist angle of $6°$⁠, and only the chord is modified to obtain the wing profiles desired. We use AMR-Wind, an LES code for wind energy applications [35,36]. AMR-Wind uses local refinement near the wing, decreasing the amount of computational resources required for the simulations. All cases used $ϵ/c=0.25$ and 600 actuator points with constant spacing to ensure that the actuator point spacing, $ϵ/Δr>10$⁠. We refine the simulations near the wing with seven levels of refinement, leading to 256 grid points along the span, and the overall domain extents are 16S in all directions, with the wing placed at the center of the domain. This level of refinement ensures that $ϵ/Δz>5$ for all cases. This resolution ensures grid-converged results within less than 1% max difference of the converged solutions. Figure 7 shows a schematic of the mesh in the middle of the domain with flow going from left to right and every refinement region highlighted in blue. This domain size ensured that the induced velocity would change within less than 0.1% of $U∞$ as the domain is extended.

The induced velocities between the theory and LES are compared in Fig. 8. There is excellent agreement between the theory and the LES, with differences in induced velocity within less than 1%. Small differences are expected because of the nonlinearity and the lack of a drag component in the filtered lifting line theory that are present in the actuator line model. This exercise serves as proof that the solution method for the filtered lifting line theory can be used to obtain accurate solutions that match the results from the actuator line model provided they use the same underlying formulations to compute the loads.

The filtered lifting line theory equations are consistent with the definition of the actuator line model. The equations are used in computational fluid dynamics to define the body forces of wings, yet solving them independently proves challenging. In this work, we reformulated the filtered lifting line theory equations as a function of the flow angle for fixed wing applications. The reformulation allows the equations to be solved numerically using fast root-finding algorithms and we performed a convergence study to understand the grid resolution requirements of the solution method.

The filtered lifting line theory equations are continuous (no singularities) and the solution converges as the number of points along the line is increased. The solutions from the filtered lifting line theory depend on the width of the Gaussian body forces ϵ and we tested a range of commonly used values. To produce converged forces whose maximum error is within 1% of the asymptotic value, a resolution of $ϵ/Δz=3$ is sufficient. When using optimal Gaussian widths of $ϵ/c≈0.25$⁠, the resolutions have a higher impact on convergence. For coarser simulations (⁠$ϵ/c>1$⁠), the integral quantities are less sensitive to resolution; however, the solutions deviate from the optimal results.

The filtered lifting line theory can be used to correct computational fluid dynamics simulations that cannot afford to resolve the fine scales required by optimal Gaussian body forces. The method presented here can be used as a standalone tool outside of computational fluid dynamics to obtain blade loading information for fixed wings. We also note that the solutions of the forces along a wing modeled using the filtered lifting line theory will depend on the width of the Gaussian body force, ϵ.

Future work should focus on studying nonuniform distribution of points, extending the solution method to account for drag forces and rotating wings, such as performed in the actuator line model for wind turbine simulations. Also, the recommended optimal values from filtered lifting line theory (⁠$ϵ/c≈0.25$⁠) should be verified and validated with results from high-fidelity blade-resolved computational fluid dynamics simulations.

The research was performed using computational resources sponsored by the U.S. Department of Energy's Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory. Funding provided by the U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Wind Energy Technologies Office. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government.

• National Renewable Energy Laboratory, operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) (Contract No. DE-AC36-08GO28308; Funder ID: 10.13039/100006233).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The code to reproduce the figures and findings in this article is available.¹

We study the computational scaling of the solution of the equations using the “df-sane” method within the SciPy library [32,33]. The computational time scales with the number of grid points. The “df-sane” algorithm is very efficient numerically, so the scaling is near linear for most of the cases and range of $ϵ/c$ values. Figure 9 shows the time scaling for different simulations with different $ϵ/c$ values. We note that the linear scaling is consistent in most of the cases, although smaller values of $ϵ/c$ require more points, and the scaling starts deteriorating.