Accuracy of high-order, discrete approximations to the lifting-line equation

Abstract

The accuracy of several numerical schemes for solving the lifting-line equation is investigated. Circulation is represented on discrete elements using polynomials of varying degree, and a novel scheme is introduced based on a discontinuous representation that permits arbitrary polynomial degrees to be used. Satisfying the Helmholtz theorems at inter-element boundaries penalises the discontinuities in the circulation distribution, which helps ensure the solution converges towards the correct, continuous behaviour as the number of elements increases. It is found that the singular vorticity at the wing tips drives the leading-order error of the solution. With constant panel widths, numerical schemes exhibit suboptimal accuracy irrespective of the basis degree; however, driving the width of the tip panel to zero at a rate faster than the domain average enables improved accuracy to be recovered for the quadratic-strength elements. In all cases considered, higher-order circulation elements exhibit higher accuracy than their lower-order counterparts for the same total degrees of freedom in the solution. It is also found that the discontinuous quadratic elements are more accurate than their continuous counterparts while also being more flexible for geometric representation.

Nomenclature

AR

aspect ratio

${a_P}$

amplitude of ${P^{th}}$ basis

b

wingspan

c

local chord

${C_D}$

drag coefficient

${C_L}$

lift coefficient

${C_{L,\alpha }}$

wing lift-curve slope

${c_{l,\alpha }}$

section lift-curve slope

e

span efficiency factor, $C_L^2/\!\left( {\pi AR{C_D}} \right)$

N

number of elements

P

polynomial degree

S

wing planform area

w

downwash velocity

y

spanwise abscissa

Subscripts

c

control point

i

induced

j, k

element indices

L, R

left- and right-side contributions

tip

wingtip

Greek Symbol

$\alpha $

angle-of-attack

${\alpha _0}$

zero-lift angle-of-attack

$\Gamma $

bound circulation strength

$\gamma $

trailing vorticity strength, $d\Gamma /dy$

$\varepsilon $

wing twist

$\eta $

isoparametric coordinate

$\xi $

non-dimensional mapping coordinate

1.0 Introduction

Surface-singularity methods are a well-established technique for predicting the inviscid characteristics of an aerodynamic configuration, including lift, pitching moment, and vortex-induced drag. The methods originate from linear potential flow, which permit exact Green’s function solutions that can be combined using linear superposition [1]. The resulting composite solution represents an exact solution to the governing equations (which are themselves an approximation of the Navier-Stokes equations). While potential flow requires the fluid domain to be irrotational, vorticity is necessary to generate lift. As such, its influence may be included by confining it to the surface, whereby it emulates the vorticity contained in the boundary layer, and in the trailing wake. Both regions are formally excluded from the potential flow domain. The singularities have a global influence on the flow field that decays with distance, meaning that the full flow field can be found by satisfying only the surface boundary condition. In other words, the governing equations do not need to be solved off the body as in Euler or Navier-Stokes computational fluid dynamics methods [1].

Notable early applications of singularity methods for lifting geometries include the well-known Joukowsky transformation, which describes two-dimensional flow about a class of aerofoils [2], and the lifting-line model for finite-span wings developed originally by Lanchester [3] and formalised by Prandtl [4] in which a bound vortex filament of varying strength, and the associated trailing vortex sheet necessary to satisfy Helmholtz’s circulation conservation theorem, is used to capture the effect of circulatory lift (per Stokes theorem) and predict the vortex-induced drag. Prandtl’s lifting-line formulation has had an undeniable impact on the field of aerodynamics, and the minimum-induced-drag solutions [5, 6] remain important benchmarks in aerodynamic design. Analogous vortex-based approaches have also been applied to screw propellers, such as by Goldstein who used it to find the minimum loss propeller [7]. The properties of a singular vortex sheet was used by Munk to construct thin-aerofoil theory by providing a strategy to support a tangential velocity discontinuity across a streamline taken to be the mean camber line [8]. Exact solutions to this theory are enabled by the Glauert integral [9]. Theodorsen [10] and Wagner [11] extended thin-aerofoil theory to include unsteady dynamics due to harmonic motion and indicial responses, respectively, with exact expressions available for small-amplitude motions.

Thin-aerofoil and finite-wing theories may be further generalised into lifting-surface theory in which the vortex sheet strength varies in both the spanwise and chordwise directions over the planform, while the strength of the trailing vortex sheet varies only in the spanwise direction [12]. Singularity methods may be further generalised to include thickness effects using various arrangements of doublets along the surface [1]. These methods are extensible to flow regimes beyond linear potential flow, such as subsonic compressible flow and linearised supersonic flow through the Prandtl-Glauert transformation. In two dimensions this manifests as a simple correction factor based on free-stream Mach number, but in three dimensions requires a detailed scaling of the geometry and application of Goethert’s rule to find local pressure coefficients.

Outside of simple geometries and certain asymptotic limits, closed-form solutions to surface-singularity representations are generally unavailable. Therefore, numerical approximations are required to either the governing equations or the enforcement of the boundary conditions. Spectral-like collocation methods1 are available for the lifting-line model applied to planar wings, such as the method of Multhopp [13]. More flexibility in representation is obtained by using singularity elements with finite spatial extent (but global influence) and polynomial strength distributions. The flow tangency condition on the surface is satisfied at a finite number of points on each element. In essence, the discrete surface-singularity solutions are exact solutions to the governing equation with approximate enforcement of the surface boundary condition. For lifting-line and lifting-surface methods, it is common to use constant-strength horseshoe vortex elements and/or vortex rings as the elements [1, 14]. Higher-order (quadratic) vortex elements have been used by Horstmann for fixed-wing applications [15], and this method has been extended to force-free wakes by Bramesfeld [16] and rotating systems by Basom [17].

For geometries with thickness, the surface is discretised using a water-tight set of panels (giving the common name of “panel methods”). Early and well-known panel methods are those Hess and Smith [18], Hess [19], and Woodward [20] which use low-order basic singularity elements. Later evolutions of the concept such as PANAIR [21, 22] used by Boeing and the Hess and Friedman method [23] used by Douglas Aircraft employed higher-order (i.e. non-constant) singularity elements on the panels. Owing to their low computational cost, lifting-line and lifting-surface methods remain in use for conceptual design and analysis frameworks for fixed-wing aircraft [24, 25], high-Mach number systems [26], wind turbines [27], propellers [28], rotorcraft comprehensive analyses [29–32] and real-time flight simulation [33]. Incorporating additional physics and design considerations remain an open research area, such as aerostructural optimisaton [34] and unsteady aerodynamics [35, 36].

A necessary consideration in the application of surface-singularity methods is the discretisation error. That is, the quality of the solution is dependent on the number of elements/panels that are used and the order of the singularity representation. A notable cautionary tale is found in the low-order lifting-surface studies of crescent-shaped wings by van Dam in which induced drag values were predicted to be below Munk’s theoretical minimum [37, 38]. Later studies by Mortara and Maughmer found this to be an artifact of discretisation and suggested the use of a Richardson extrapolation factor to better approximate the continuum solution [39], but their recommendations are more empirical than theoretical. The vast majority of extant lifting-line and lifting-surface methods use low-order (i.e. constant) singularity elements for the simplicity of implementation and fast computation. In contrast, two-dimensional panel methods for aerofoils that remain in widespread use, such as the XFOIL [40] and PROFIL [41] analysis and design codes, employ linear and quadratic strength panels, respectively.

By the nature of discretisation error, higher-order representations should provide superior accuracy for a given number of solution unknowns compared to lower-order methods. However, this has not been rigorously quantified in the literature to the best of the author’s knowledge. Additionally, existing higher-order methods have singularity-strength continuity requirements that make their implementation more cumbersome and, in the case of thin geometries, constrain their applicability to complex topologies. The objective of the present work is to quantify the errors and formal order of accuracy achieved by higher-order surface-singularity methods. The focus here is restricted to solutions of the lifting-line problem for planar wings. This is for two primary reasons: First, there are both exact and spectrally accurate solutions available for quantifying numerical error. Second, lifting-line methods, and the extension to lifting-surface methods using the same types of singularity elements, represent the dominant usage case, thereby providing a practical aspect to the study. A novel higher-order discretisation strategy is presented in this work that alleviates the inter-element continuity constraints while maintaining consistency to arbitrary order of accuracy, and shows improved accuracy over schemes that strictly enforce circulation continuity between elements. The concepts being investigated and associated lessons learned are readily extensible to general lifting surfaces.

2.0 Background

Lifting-line theory is the high-aspect-ratio asymptotic limit of lifting-surface theory. In this limit, chordwise variations in the bound vorticity, including those due to sweep, are assumed to occur over a length scale that is much, much smaller than the length scale associated with spanwise variations. This limit enables the chordwise and spanwise characteristics to be decoupled, allowing each spanwise station to be modeled as a two-dimensional aerofoil subject to a chordwise-constant downwash field induced by the trailing vortex sheet. By relating the local circulation strength to the induced angle-of-attack from the trailing wake, the governing equation for a planar wing,

(1) \begin{align}\Gamma = {1 \over 2}{U_\infty }c{c_{l,\alpha }}\left( {\alpha - {\alpha _0} - \varepsilon } \right) - {{c{c_{l,\alpha }}} \over {8\pi }}\int_{ - b/2}^{b/2} {{{{d\Gamma } \over {d{y_0}}}} \over {y - {y_0}}}d{y_0}\end{align}

may be obtained in which c varies along the span and the aerofoil lift-curve is assumed to be linear. The integral in the final term of Eq. 1 is a direct result of the trailing vorticity. Through a well-known coordinate transformation and assuming a Fourier sine series expansion of $\Gamma $ , the integral may be converted to a Glauert-type integral [9] for which analytical solutions are available. It is typical for expressions of the bound circulation to enforce a zero-lift constraint on the wing tips. While physically justified, these boundary conditions are enforced naturally by the system and do not need to be strictly imposed.

For general chord distributions, exact solutions to Eq. 1 are unavailable, necessitating the use of numerical methods. The most common methods documented in the literature are all based on approximation theory in which the lifting-line equation is satisfied at discrete control points along the wingspan. Although the equation is not satisfied in between the control points, the global error distribution is the best available approximation of zero error from the discrete approximation of the circulation distribution. The method of Multhopp [13] takes advantage of this by selecting the Chebyshev zero points on $[\!-\!b/2,b/2]$ to provide the most-uniform zero approximation in conjunction with the analytic solution of the Glauert integral.

To allow greater flexibility in the representation of the circulation distribution, such as for modeling non-planar and multiple wings, it is common to express $\Gamma $ using piecewise polynomials. One widely used option is to use constant-strength horseshoe vortices placed along the quarter-chord as illustrated in Fig. 1. For each horseshoe vortex, there is a single unknown circulation strength, necessitating the lifting-line equation being satisfied at a number of points equal to the number of horseshoe vortices. The width of each horseshoe need not be constant, and it is typical for the control point to be placed at the midpoint of the bound portion of each vortex [14]. Some implementations employ clustering, such as cosine distributions (e.g. the Chebyshev points) to increase the resolution near geometric discontinuities such as taper breaks, dihedral breaks, and wingtips.

Figure 1. Linearly tapered wing approximated using horseshoe vortices.

Higher-order lifting line methods have been proposed in which the circulation is modelled as varying quadratically across each vortex element. Between adjacent elements, the strength of both the circulation, $\Gamma $ , and trailing vorticity, $d\Gamma /dy$ , are required to be continuous, thereby eliminating any singular velocities away from geometric discontinuities [15, 16]. As such a representation contains $3N$ unknowns with $2N - 2$ continuity constraints, it is well-closed by requiring $\Gamma = 0$ at both wingtips and that the lifting line equation is satisfied at a single control point on each element. This strategy, especially its lifting-surface counterpart [15, 16], have been used with great success in aerodynamic wing design, particularly for non-planar wings [42–44]. A limitation in this continuous representation is that it lacks natural handling of complex geometries such as split winglets or lifting struts prevalent on contemporary and next-generation high-efficiency aircraft as the inter-element continuity constraints become underdefined.

In this work, a novel high-order numerical strategy is employed that is based on discontinuous singularity elements. In doing so, arbitrarily high polynomial orders may be consistently used without needing to find a “sweet spot” in the total number of element control points and continuity constraints equaling the total number of unknowns as is the case with Horstmann’s method [15].

3.0 Methodology

For the following study, polynomial representations of varying order will be considered for representing the circulation on each vortex element. The maximum polynomial degree will be represented as “P”, and the number of control points as “Q”. For example, the aforementioned horseshoe vortices would be a “P0Q1” discretisation, denoting constant polynomials and a single control/quadrature point for each element.

3.1 Basic vorticity elements

For the present study, the polynomial circulation distribution on each element is represented as a Legendre series,

(2) \begin{align} \Gamma = {a_0} + {a_1}\eta + {{{a_2}} \over 2}\left( {3{\eta ^2} - 1} \right) + \ldots \end{align}

where $\eta $ is an isoparametric coordinate defined as

(3) \begin{align} \eta = {{y - {y_{0,j}}} \over {\Delta {y_j}}}\end{align}

that spans $[\!-\!1,1]$ on each element so that the discrete orthogonality relations of the Legendre polynomials hold. While this is not a necessary condition, numerical experimentation has shown it to be advantageous by improving the condition number of the resulting linear system compared to using a monomial basis. It is convenient to define the induced velocity due to isolated vortex elements in terms of influence coefficients,

(4) \begin{align} {w_i} = {I_0}{a_0} + {I_1}{a_1} + {I_2}{a_2} + \ldots \end{align}

For a planar wing where all vortex elements and control points are coplanar, the first three influence coefficients are defined as

(5) \begin{align} {I_0} = {1 \over {4\pi \Delta {y_j}}}\left( {{2 \over {1 - {\eta ^2}}}} \right) \end{align}

(6) \begin{align} {I_1} = {1 \over {4\pi \Delta {y_j}}}\left( {{{2\eta } \over {1 - {\eta ^2}}} - \ln \left| {{{1 - \eta } \over {1 + \eta }}} \right|} \right) \end{align}

and

(7) \begin{align} {I_2} = {I_0} + {3 \over {4\pi \Delta {y_j}}}\left( {\eta \ln \left| {{{1 + \eta } \over {1 - \eta }}} \right| - 2} \right) \end{align}

This study considers up to the P2 representation as the results can be generalised to lifting-surface methods, but for lifting-line applications it can be extended to arbitrarily high order. Expressions for non-planar wings may also be derived if desired, and will include both downwash and spanwise velocity components.

3.2 Enforcement of kinematic conditions

3.2.1 Continuous elements

For P2 (and higher degree) elements, a piecewise $C1$ representation (i.e. continuous $\Gamma $ and $\gamma = d\Gamma /dy$ ) may be obtained by satisfying the lifting-line equation at a number of points on each element equal to one fewer than the polynomial degree. Quadratic elements, as considered here, would constitute a P2Q1-C1 scheme. The system of equations for each $j{\rm{th}}$ interior vortex element is written as the lifting line equation defined at the control point of element j,

(8) \begin{align} {a_{0,j}} + {a_{1,j}}{\eta _{c,j}} + {{{a_{2,j}}} \over 2}\left( {3\eta _{c,j}^2 - 1} \right) = {1 \over 2}{U_\infty }{c_j}{c_{l,\alpha }}\left( {\alpha - {\alpha _0} - {\varepsilon _j}} \right) - {{{c_j}{c_{l,\alpha }}} \over 2}\sum\limits_{k = 1}^N \sum\limits_{p = 1}^3 {I_{p,j,k}}{a_{p,k}} \end{align}

in which ${I_{p,j,k}}$ denotes the influence coefficient of the $p{\rm{th}}$ basis of element k on element j, ${\eta _{c,j}}$ is the non-dimensional position of the control point on element j, and ${c_j}$ is the wing chord at the control point, along with the continuity constraints for bound circulation,

(9) \begin{align} {a_{0,j - 1}} + {a_{1,j - 1}} + {a_{2,j - 1}} = {a_{0,j}} - {a_{1,j}} + {a_{2,j}} \end{align}

and wake strength,

(10) \begin{align} {1 \over {\Delta {y_j}}}\left( {{a_{1,j}} + 3{a_{2,j}}} \right) = {1 \over {\Delta {y_{j + 1}}}}\left( {{a_{1,j + 1}} - 3{a_{2,j + 1}}} \right) \end{align}

At the left wingtip, Eq. 9 is replaced with,

(11) \begin{align} {a_{0,1}} - {a_{1,1}} + {a_{2,1}} = 0 \end{align}

and at the right wingtip, Eq. 10 is replaced with,

(12) \begin{align} {a_{0,1}} + {a_{1,1}} + {a_{2,1}} = 0 \end{align}

Unless otherwise specified, the control point for each continuous quadratic panel is taken to the be its midpoint, ${\eta _{c,j}} = 0$ . This scheme is qualitatively illustrated in Fig. 2.

Figure 2. Approximation of a continuous function (dashed line) using piecewise quadratic elements satisfying C1 continuity.

3.2.2 Discontinuous elements

Piecewise polynomials of any degree p may be used if discontinuities in both circulation and wake strength are permitted between elements, as qualitatively illustrated in Fig. 3. Without any continuity constraints, the system of equation is closed by evaluating the discrete lifting-line equation at $p + 1$ control points on each panel, which are chosen to be the Gauss-Legendre quadrature points. In turn, the leading error term is the next higher Legendre polynomial, and thus the representation on each element is optimal in a Galerkin least squares sense. Additionally, the Gauss-Legendre points provide the optimal order of accuracy for integration, which for this case is order $2p + 2$ , thereby improving the accuracy of the calculated induced drag.

Figure 3. Discontinuous approximations of a continuous function (dashed line) with Legendre control points (symbols).

Jumps in $\Gamma $ and $d\gamma $ between elements are effectively penalised, i.e. discouraged but not explicitly eliminated in the final solution, by singularities in the induced velocity fields. Discontinuities in $\Gamma $ result in the presence of discrete trailing vortices, and the velocity penalty for a jump at ${y_{j + 1/2}}$ takes the form,

(13) \begin{align} {w_{i,\Delta \Gamma }} = {{{\Gamma _R} - {\Gamma _L}} \over {y - {y_{j + 1/2}}}} \end{align}

Discontinuities in $\gamma $ create a jump in the strength of the trailing vorticity sheet, and the penalty is of the form,

(14) \begin{align} {w_{i,\Delta \gamma }} = \left( {{\gamma _R} - {\gamma _L}} \right)\ln \!\left( {y - {y_{j + 1/2}}} \right) \end{align}

which is a weaker singularity. In the present formulation, the penalties need not be explicitly imposed as they are included in the definitions of the influence coefficients.

3.2.3 Overconstrained least-squares solution

For both of the aforementioned strategies, the number of control points on each panel is determined based on the degree of the polynomial basis so that the resulting solution is unique. Although the discrete solution approximates the lifting-line equation, and the error is orthogonal to the solution (by the definition of the Legendre basis), they are not guaranteed to be small on a given element. In fact, the downwash error may be unbounded at specific points while the global, integrated error remains bounded. To evaluate this effect in the present study, cases are evaluated in which the residual error is overconstrained on each element, i.e. a greater number of Legendre control points are used on each element than are required for a unique solution, to hopefully provide a more uniform approximation of zero error. An example would be a P2Q4 discretisation, which would have three unknown basis amplitudes, but four control points per element. While it is impossible to exactly satisfy the lifting-line equation at every control point due to there being more equations than unknowns, the integrated error over the element is minimised in a least-squares sense.

4.0 Results

Comparisons of the various schemes are performed using a planar, untwisted, elliptical wing with span $b = 10$ and root chord ${c_0} = 1$ ( $AR \approx 12.7324$ ), and a planar, untwisted rectangular wing with the same span and root chord ( $AR = 10$ ). The aerofoil sections are assumed to be ideal with a section lift-curve slope of ${c_{l,\alpha }} = 2\pi \,{\rm{ra}}{{\rm{d}}^{ - 1}} \approx 0.109662271123\,{\rm{de}}{{\rm{g}}^{ - 1}}$ . The elliptical wing possesses a well-known analytical solution to the lifting-line equation, which for this case yields a three-dimensional lift curve slope of ${C_{L,\alpha }} = 0.094775042292695\,{\rm{de}}{{\rm{g}}^{ - 1}}$ and a span efficiency factor of exactly one. A closed-form analytical solution for the rectangular wing is not available; however, a highly precise numerical solution can be constructed using a spectrally accurate method. The resulting wing lift-curve slope and span efficiency factors of ${C_{L,\alpha }} = 0.08808311706\,{\rm{de}}{{\rm{g}}^{ - 1}}$ and $e = 0.9208891958$ , respectively, and are accurate to ten significant digits. Solutions for all schemes were sought for up to $N = 2560$ elements. In some cases, however, finer resolutions resulted in lost of precision due to very large matrix condition numbers (e.g. greater than $1 \times {10^9}$ ). Nevertheless, asymptotic convergence behaviour is observed.

4.1 Uniform distribution of elements

The accuracy of the P0Q1, P1Q2, P2Q3 and P2Q1-C1 schemes for predicting ${C_{L,\alpha }}$ and e of the two planforms have been investigated by dividing the wings into N uniformly distributed elements (i.e. each element has constant width). Convergence of the numerical solution of the elliptical wing to the analytical solution with increasing N are shown in Figs. 4 and 5 for the lift-curve slope and span efficiency, respectively. For each figure, the number of DOFs is equal to the number of unknowns per elements multiplied by the number of elements (e.g., a P2 discretisation with 10 elements will have 30 DOFs). In this sense, the error as a function of computational expense, i.e. system matrix size, can be directly compared across the schemes. All discontinuous schemes approach the lift and span efficiency monotonically from above, whereas the continuous scheme approaches from below. The higher-order schemes show significantly reduced error for a given expense. Conversely, a given error tolerance may be achieved with a significantly reduced computational expense, as a P2 scheme requires $\mathcal{O}{(10^2})$ times fewer elements relative to P0. The continuous P2 scheme is observed to have greater error than its discontinuous counterpart, but the errors approach one another as the number of DOFs increases. The P0Q1 discretisation exhibits the expected ${1^{{\rm{st}}}}$ -order accurate behaviour for both functionals; however, the P1Q2 shows sublinear convergence and the two P2 schemes both show ${1^{{\rm{st}}}}$ -order convergence. For all higher-order schemes, the observed order of accuracy is below that of the expected $p + 1$ order of the schemes as determined from the polynomial basis.

Figure 4. Convergence of the lift-curve slope of an elliptical wing with uniform element widths.

Figure 5. Convergence of the span-efficiency factor of an elliptical wing with uniform element widths.

In an attempt to recover the design order of accuracy of the P1Q2 scheme, the P1Q3 and P1Q4 overconstrained schemes were studied. The resulting lift-curve slope error convergence is shown in Fig. 6. From these data, it is found that the overconstrained least-squares approach does not mitigate error in the domain, but rather forces the solution towards an incorrect state with non-zero residual error.

Figure 6. Error convergence of the lift-curve slope with overconstrained P1 schemes.

As an unanticipated reduction in order of accuracy was observed for the elliptical wing, the same study was repeated with a rectangular wing to evaluate whether or not there is an influence of the chord distribution and an associated polynomial representation. The results for the rectangular wing are plotted in Figs. 7 and 8. The convergence of the lift-curve slope is qualitatively similar to that for the elliptical wing in terms of the schemes approaching the the true answer from above or below, as well as the observed order of accuracy. The P2Q1-C1 scheme, however, showed improved relative accuracy for coarser discretisations, and the error was more consistent with P2Q3 across the range of DOFs. For span-efficiency factor, the discontinuous schemes again showed the same trends as observed before, but the order of accuracy for the P2Q1-C1 scheme was greater ( $ \sim\!1.2$ ) than for the discontinuous scheme.

Figure 7. Convergence of the lift-curve slope of a rectangular wing with uniform element widths.

Figure 8. Convergence of the span-efficiency factor of a rectangular wing with uniform element widths.

Although the elliptical wing has a well-behaved solution to the governing equations, the chord distribution is not conducive to being represented using a polynomial. The rectangular wing, however, removes spanwise variation in chord, allowing the solution to be fully focused on representing the circulation distribution. Nevertheless, general solutions to the lifting-line equation exhibit a singular behaviour of $d\Gamma /dy$ at the wing tips, which is poorly represented by polynomial bases. With the discontinuous schemes, the infinite slope is approximated by non-zero circulation at the tip and the associated jump penalty. The value of $\Gamma $ at the tip does approach zero as the number of elements increases (confirming some consistency of the schemes), but as is illustrated by the error convergence plotted in Fig. 9 for the elliptical planform, it converges at a rate of ${N^{ - 1/2}}$ irrespective of the element polynomial order. The continuous P2Q1-C1 scheme strictly enforces $\Gamma = 0$ at the tip, so it is expected that the discretisation error is absorbed into the solution in a much different manner. To better understand the behaviour, the errors in predicted circulation distribution of the elliptical wing with $N = 8$ is plotted in Fig. 10 for the P2Q1-C1 and P2Q3 schemes. The continuous scheme exhibits a strong underprediction of lift towards the tip due to its finite $d\Gamma /dy$ . Because of the continuity constraints, the error propagates inboard as a low-frequency oscillation. In contrast, the discontinuous, P2Q3 scheme absorbs the tip error more quickly via the jump penalties at the boundaries of the outboard-most element. The prediction errors in the tip region of the rectangular wing were also evaluated, and found to follow identical trends to the elliptical wing. Therefore, those results have been omitted for the sake of brevity.

Figure 9. Convergence of elliptial-planform wing-tip circulation using discontinuous schemes.

Figure 10. Error in circulation distribution using continuous (red) and discontinuous (black) quadratic elements with $N = 8$ .

4.2 Observations on integration error

The influence of the $d\Gamma /dy$ wingtip singularity on global order of accuracy is well-illustrated by numerically approximating the area under a semi-circle given by $f(x) = \sqrt {1 - {x^2}} $ on $[\!-\!1,1]$ using 1-, 2-, and 3-point Gauss-Legendre quadrature rules, which are ${2^{{\rm{nd}}}}$ -, ${4^{{\rm{th}}}}$ -, and ${6^{{\rm{th}}}}$ -order accurate, respectively. First considered are constant-width subintervals, which are attainable via a linear mapping function, and the error convergence is plotted below in Fig. 11. It is observed that the formal order of accuracy of the integral does not increase as the order of the quadrature rule increases. Rather, all three quadrature rules exhibit an observed order of accuracy of 1.50 rather than their design order. Investigation of the subinterval behaviour confirms that the reduction in accuracy is due to the singularities, as the contribution from the domain interior converges with to the design order.

Figure 11. Error convergence for uniformly distributed subintervals.

An improvement on the sub-optimal accuracy is available through the use of non-uniform subinterval distributions that cluster quadrature points in regions with large slopes. For example, numerous lifting-line codes, such as AVL [24], employ cosine distributions to reduce the element width at the wingtips. Such a distribution is a mapping of the form,

(15) \begin{align} y = - {b \over 2}\cos\!\left( {\xi \pi } \right) \end{align}

where $\xi = [0,1]$ , and elements are distributed uniformly in $\xi $ , i.e. $\Delta \xi = 1/N$ . Analysis of the cosine mapping near the endpoints show that the widths of the first and last elements vanish with $\Delta {\xi ^2}$ , which contrasts with the uniform distribution in y for which the endpoint widths vanish linearly in $\Delta \xi $ . The cosine mapping may be more broadly interpreted as a sigmoid function that produces symmetric distributions about the wing centreline. Polynomial-based mappings may also be generated, such as cubic,

(16) \begin{align} y = - {b \over 2} + b\!\left( {3{\xi ^2} - 2{\xi ^3}} \right)\end{align}

quintic,

(17) \begin{align} y = - {b \over 2} + b\!\left( {10{\xi ^3} - 15{\xi ^4} + 6{\xi ^5}} \right)\end{align}

and septic,

(18) \begin{align} y = - {b \over 2} + b\!\left( {35{\xi ^4} - 84{\xi ^5} + 70{\xi ^6} - 20{\xi ^7}} \right)\end{align}

for which the endpoint element widths vanish with $\Delta {\xi ^2}$ , $\Delta {\xi ^3}$ , and $\Delta {\xi ^4}$ , respectively. These distributions are plotted in Fig. 12, while the effect of the mapping on the error convergence is illustrated in Fig. 13. With 1-point Gauss-Legendre quadrature, all mappings provide the design ${2^{{\rm{nd}}}}$ -order accuracy, further confirming that the endpoints are the cause of reduced accuracy with uniform interval widths. The 2-point quadrature rule is limited to only ${3^{{\rm{rd}}}}$ -order accuracy with the cosine and cubic mappings, while it approaches ${4^{{\rm{th}}}}$ -order accuracy with the quintic mapping and fully exhibits this accuracy with the septic mapping. The 3-point rule should exhibit ${6^{{\rm{th}}}}$ -order accuracy, but it is also limited to ${3^{{\rm{rd}}}}$ -order with the cosine and cubic mappings. The quintic mapping shows a significant improvement, though the accuracy is tending towards only ${5^{{\rm{th}}}}$ -order, whereas with the septic mapping it recovers design order as the number of subintervals increases.

Figure 12. Element/subinterval distribution functions.

Figure 13. Error and order of accuracy for various subinterval distributions and quadrature rules.

It is also observed that the cosine distribution offers a slight improvement over the cubic distribution for the higher-order quadrature rules. Although small, it is attributed to the ${3^{{\rm{rd}}}}$ derivative of the cosine mapping being zero at the endpoints, allowing the asymptotic behaviour of the mapping to be achieved for smaller N than the cubic mapping.

4.3 Non-uniform distribution of elements

The lessons learned about the influence of the element distribution on the order of accuracy have been applied back to the solution of the lifting-line problem. Here, an elliptic wing is considered with cosine, quintic and septic element mapping functions. Convergence plots of ${C_{L,\alpha }}$ and e are includes as Figs. 14 - 16. For all element distributions, the P0Q1 scheme exhibits the expected ${1^{{\rm{st}}}}$ -order accuracy in both functional quantities. The P2Q3 scheme exhibits ${2^{{\rm{nd}}}}$ -order accuracy in both quantities for the cosine distribution, and even higher accuracy for quintic and septic distributions. It is found that the formal accuracy with the quintic distribution is $ \sim 2.6$ , and that the order approaches 3 with the septic distribution before numerical truncation error begins to obscure the behaviour. The P2Q1-C1 scheme shows the same order accuracy as the P2Q3 scheme for the lift coefficient; however, it is limited to ${2^{{\rm{nd}}}}$ -order accuracy in span efficiency. This is unsurprising as ${C_L}$ is a direct consequence of the polynomial basis, while span efficiency is calculated using the associated quadrature rule.

Figure 14. Error convergence with elements distributed using a cosine mapping.

Figure 15. Error convergence with elements distributed using a quintic polynomial mapping.

Figure 16. Error convergence with elements distributed using a septic polynomial mapping.

An unexpected result is that the P1Q2 scheme exhibits ${1^{{\rm{st}}}}$ -order accuracy for all non-uniform distributions, rather than recovering ${2^{{\rm{nd}}}}$ -order accuracy as implied by the basis. Although the definitive mechanism behind the reduced accuracy is unknown, one hypothesis is that it is a result of the penalisation of solution discontinuities between elements. Value and slope jumps introduce singularities of strength $1/y$ and $\ln y$ , respectively, that vanish as the number of elements increases. The P1Q2 scheme has no mechanism for satisfying C1 continuity across the span, whereas the P2Q3 scheme does. The resulting slope discontinuities for P1Q2 decay slowly as N increases, and the total contribution of the penalty scales linearly with element width along the whole span. An alternative hypothesis is that this may be a consequence of even- versus odd-order bases as is known to happen in some discontinuous Galerkin finite-element schemes [45]. Unfortunately, results with a ${3^{{\rm{rd}}}}$ -order scheme (P3Q4) could not be obtained without machine precision errors at sufficiently large N to demonstrated formal order of accuracy and evaluate both hypotheses.

5.0 Conclusion

From the preceding analyses presented in this paper, some important takeaways can be made regarding the accuracy of numerical solutions to the lifting-line equation. First, the singular nature of $d\Gamma /dy$ at the wing tips is poorly modeled by polynomial bases of any order, which introduces an error proportional to $\Delta y_{tip}^{0.5}$ into the solution. For discretisations with uniform element widths, the tip error dominates the error of integrated lift and drag coefficients (measured here using the lift-curve slope and span-efficiency factors). Consequently, the observed order of accuracy is, at most, linear, irrespective of the polynomial order.

Second, improved order of accuracy may be recovered by clustering the elements near the wing tips so that the panel widths go to zero faster than ${N^{ - 1}}$ in these regions. The specific rate at which the tip elements approach zero width affects the observed order of accuracy of a given scheme. For example, the cosine distribution ( $\Delta {y_{tip}} \sim {N^{ - 2}}$ ) allows ${2^{{\rm{nd}}}}$ -order accuracy to be attained from the schemes with a quadratic basis, while the quintic ( $\Delta {y_{tip}} \sim {N^{ - 3}}$ ) and septic ( $\Delta {y_{tip}} \sim {N^{ - 4}}$ ) allow full ${3^{{\rm{rd}}}}$ -order accuracy to be recovered. Although clustering near tips and geometric discontinuities is already well known to improve accuracy and an accepted practice, this study provides new mathematical guidance for how clustering should be performed to achieve best accuracy for a given discretisation.

Third, the linear basis does not exhibit greater than ${1^{{\rm{st}}}}$ -order accuracy, regardless of the element-width mapping function. Although the numerical results are consistent (i.e. they converge to the correct continuum value, provided the system is not overconstrained), but the net effect of the penalties caused by jumps in $\Gamma $ and $d\Gamma /dy$ between elements decays more slowly than anticipated. Nevertheless, it exhibits improved accuracy over the constant-strength elements for a given number of solution degrees of freedom.

Finally, the discontinuous quadratic basis (P2Q3 scheme) exhibits both very high accuracy, but with additional flexibility of discretisation as compared to a continuous quadratic scheme (P2Q1-C1) similar to the one used by [15]. The P2Q3 scheme shows greatly reduced error compared to the P0 scheme for lift and drag with the same number of solution degrees of freedom, and with sufficient element clustering, it exhibits ${3^{{\rm{rd}}}}$ -order accuracy as opposed to the other schemes’ ${1^{{\rm{st}}}}$ -order accuracy.

The observations and conclusions drawn from this study are anticipated to be applicable to lifting-surface methods as well as lifting line. Both methods use the same singularity forms, and are also subject to the same tip singularities. The discontinuous methodology should be broadly applicable, though additional study is likely required before it can be used with relaxed/free wake analyses. The results of this study also suggest potential benefits from adapting the discretisation/paneling to minimise solution error, such as for the discretely approximated flow tangency conditions or for output functionals like lift and induced drag, which may be accomplished using an adjoint-based sensitivity analysis.