2.0 REVIEW OF NUMERICAL LIFTING LINE MODEL AND DOWNWASH

Following Collar^{(Reference Collar1)} and McCormick^{(Reference Mccormick2)}, a wing of span b is divided into n wing segments with each segment having centre points located at p_i. Each segment has a horseshoe vortex with strength Γ_i (Fig. 1).

Figure 1. Vortex lattices.

Define a vector of vortex filament strengths

$${\bf \Gamma} = [\Gamma _1 \cdots \Gamma _N ]^T $$

Define w_i as the downwash on the wing segment at point p_i due to the streaming vortices, calculated using the Biot-Savart Law. Noting that p_i is the location of the centre of the ith wing section, write for the downwash elements:

(1) $$ w_i = {1 \over {4\pi }}\sum\limits_{j = 1}^n {\left( {{1 \over {p_i - p_j + {{\Delta x} \over 2}}} - {1 \over {p_i - p_j - {{\Delta x} \over 2}}}} \right)\Gamma _j } $$

Index i relates to the point p_i as $ p_i = i\Delta x - {{\Delta x} \over 2} $, and so $p_i - p_j = \Delta x(i - j)$. Also using $\Delta x = b/n $:

$$ w_i = {n \over {\pi b}}\sum\limits_{j = 1}^n {\left( {{1 \over {1 - 4(i - j)^2 }}} \right)\Gamma _j .} $$

Define a kernel matrix Q with elements q _i,j such that

(2) $$ q_{i,j} = {1 \over {1 - 4(i - j)^2 }}, $$

then the downwash vector w is

(3) $$ {\bf{w}} = {1 \over b}{n \over \pi }Q{\bf{\Gamma }} $$

Note that while the span b is in the definition of the downwash, the matrix Q is dimensionless. Without loss of generality for an arbitrary wing, we can choose a system of units such that span b = 1, simplifying the presentation to

(4) $$ {\bf{w}} = {n \over \pi }Q{\bf{\Gamma }} $$

Collar^{(Reference Collar1)} developed a similar but more generalised expression where he defined a matrix U as

(5) $$ U_{i,j} = {{ - 1} \over {(k + i - j)(k + i - j - 1)}} $$

which resulted in his expression for the downwash as

(6) $$ {{4\pi } \over n}{\bf{w}} = U{\bf{\Gamma }} $$

His term k could vary between 0 and 1, thus allowing for the point at which the downwash was calculated to vary between the two filaments and not just in the middle (k = 1/2). When k = 1/2, then U = 4Q. The scaling difference between Equations (3) and (6) likely hid the property related to the triangular number.

3.0 PROPERTIES OF THE MATRIX Q

Look at Q as a 4 × 4 matrix along with its inverse:

$$ \matrix{ {Q_4 = \left[ {\matrix{ 1 & { - {1 \over 3}} & { - {1 \over {15}}} & { - {1 \over {35}}} \cr { - {1 \over 3}} & 1 & { - {1 \over 3}} & { - {1 \over {15}}} \cr { - {1 \over {15}}} & { - {1 \over 3}} & 1 & { - {1 \over 3}} \cr { - {1 \over {35}}} & { - {1 \over {15}}} & { - {1 \over 3}} & 1 \cr } } \right]} & {Q_4 = \left[ {\matrix{ {{{1225} \over {1024}}} & {{{525} \over {1024}}} & {{{315} \over {1024}}} & {{{175} \over {1024}}} \cr {{{525} \over {1024}}} & {{{1425} \over {1024}}} & {{{615} \over {1024}}} & {{{315} \over {1024}}} \cr {{{315} \over {1024}}} & {{{615} \over {1024}}} & {{{1425} \over {1024}}} & {{{525} \over {1024}}} \cr {{{175} \over {1024}}} & {{{315} \over {1024}}} & {{{525} \over {1024}}} & {{{1225} \over {1024}}} \cr } } \right].} \cr } $$

The matrix Q is bi-symmetric and is a Real Hermitian Toeplitz matrix and makes for an excellent test case for ongoing research in the properties of infinite Toeplitz matrices^{(Reference Bogoya, Böttcher, Grudsky and Maximenko12,Reference Böttcher, Bogoya, Grudsky and Maximenko13)} . It is also an M-Matrix with many well-known properties that follow^{(Reference Poole and Boullion14)}. As Q is a symmetric M-matrix, it is positive definite and its inverse is a positive matrix Q ⁻¹ > 0, i.e. all of the elements of Q ⁻¹ are positive. Within the space of positive matrices, the matrix sum of elements can serve as a matrix norm, with a number of interesting properties^{(Reference Merikoski15)}. Define the function Su(X) as the sum of the elements of matrix X, which is also equal to e^T X e.

The matrix Q decomposes into either the sum or the Hadamard product of a simpler matrix H where the elements of H are

(7) $$ H_{i,j} = {1 \over {1 + 2(i - j)}} $$

so

$$ Q = {1 \over 2}(H + H^T ) $$

or alternatively

$$ Q = H \otimes H^T $$

where ⊗ indicates the Hadamard product (an element by element multiplication). The matrix H is central to the 2D vortex lattice model, as shown below, and has its own set of interesting properties.

4.0 ELEMENT SUM OF THE INVERSE OF Q

The matrix Q has the surprising property that the sum of the elements of the inverse of $ Q_n^{ - 1} $ is always an integer, and is the Triangular Number T(n), proof at the end of the paper. The number satisfies the relations

$$ \matrix{ \eqalign{ {\rm{Su}}(Q_n^{ - 1} ) = & T(n) = 1 + 2 + \cdots + n = \left( {\matrix{ {n + 1} \cr 2 \cr } } \right) = {{n(n + 1)} \over 2} \cr = & 1,3,6,10,15,21, \ldots \cr} & . \cr } $$

While the sum of all the elements is an integer based on numerical experiments over a variety of matrix sizes, it appears that none of the sums of subsets of elements are integers (excepting special cases when choosing n/2 elements and n is even). It would be an interesting problem to confirm or deny this property.

5.0 ROW SUMS OF THE INVERSE OF Q

In the investigation of Q, the next property identified after noting the relationship with the triangular number was for the row sums of Q ⁻¹. The reason this is of interest is that, in the constant downwash case, we can compute the distribution of Γ (and therefore lift) along the wing by inverting Equation (4) to get Γ∝Q ⁻¹e where $ {\bf{e}} = \left[ {\matrix{ 1 & 1 & \cdots & 1 \cr } } \right]^T $, so for constant downwash the lift is distributed as per the row sums of Q ⁻¹. We anticipate the row sums to tend towards an ellipse.

Collar did not derive this quantity in his paper, and the author has not found the result anywhere else. The expression for the elements of the vector r = Q ⁻¹e can be written in many different forms, as is typical for combinatorial quantities. Writing it in terms of central binomial coefficient terms gives

(8) $$ {\bf{r}}_k - {{k(2n - 2k + 1)} \over {2^{2n - 1} }}\left( {\matrix{ {2k} \cr k \cr } } \right)\left( {\matrix{ {2n - 2k} \cr {n - k} \cr } } \right) $$

The benefit of writing it with the central binomial coefficients terms is that there are tight bounds that are also good approximations for the coefficients:

$$ {{4m} \over {\sqrt {\pi (m + 1/2)} }} \le \left( {\matrix{ {2m} \cr m \cr } } \right) \le {{4m} \over {\sqrt {\pi m} }}. $$

To show that r_k tends to an ellipse, look at a general term to use to approximate the central binomial coefficient,

$$ \left( {\matrix{ {2m} \cr m \cr } } \right) \simeq {{4m} \over {\sqrt {\pi (m + \beta )} }}. $$

Substituting into Equation (8) and simplifying gives

$$ {\bf{r}}_k \simeq {{2k(2n - 2k + 1)} \over {\pi \sqrt {(k + \beta )(n - k + \beta )} }} $$

Substitute k = nz where $ z \in (0,1] $ and let $ n \to \infty $, scaling by a 1/n to normalise the result:

(9) $$\matrix{ {{{\bf{r}}_z} = \mathop {\lim }\limits_{n \to \infty } {1 \over n}\left( {{{2nz(2n - 2nz + 1)} \over {\pi \sqrt {(nz + \beta )(n - nz + \beta )} }}} \right)} \hfill \cr {\quad = {4 \over \pi }\sqrt {z(1 - z)} .} \hfill \cr } $$

Equation (9) is the equation of an ellipse and corresponds to the exact solution for the constant downwash condition for a wing with a span b=1. Both the lower and upper bound on the row sums tend to the same ellipse, confirming the shape of the row sums. Figure 2 shows the correspondence between Equations (8) and (9), after scaling the exact result.

Figure 2. Comparing exact solution (solid line) to row sums (points) for inverse of Q, with every 5th point shown.

We can closely approximate the row sums with

(10) $$ {\bf r}_k \simeq {4 \over \pi }\sqrt {\left( {k - {1 \over 2}} \right)\left( {n - k + {1 \over 2}} \right)} $$

If tighter bounds for the central binomial coefficient are needed one can use the following with the upper bound particularly tight^{(Reference Johnson16)}:

$$ {{4^m } \over {\sqrt {\pi (m + 1/3)} }} \le \left( {\matrix{ {2m} \cr m \cr } } \right) \le {{4^m } \over {\sqrt {\pi (m + 1/4)} }}. $$

6.0 THE ELEMENTS OF THE INVERSE OF Q

Only after the relation for the row sums had been identified did I become aware of Collar’s previous work. Collar⁽¹⁾ calculated the inverse of U = 4Q and proved that the elements of U ⁻¹ can be written as

$$ U_{i,j}^{ - 1} = \left( {{1 \over 2}} \right)^{2n} ( - 1)^{n - j} {{(2j - 1) \cdots (2j - 2n - 1)} \over {(j - 1)!(n - j)!}}\sum\limits_{r = 1}^i {{( - 1)^r (2n + 1 - 2r) \cdots (3 - 2r)} \over {(r - 1)!(n + 1 - r)!(2j + 1 - 2r)}} $$

To gain the elements of Q ⁻¹, scale by a factor of 4,

$$ Q_{i,j}^{ - 1} = \left( {{1 \over 2}} \right)^{2n - 2} ( - 1)^{n - j} {{(2j - 1) \cdots (2j - 2n - 1)} \over {(j - 1)!(n - j)!}}\sum\limits_{r = 1}^i {{( - 1)^r (2n + 1 - 2r) \cdots (3 - 2r)} \over {(r - 1)!(n + 1 - r)!(2j + 1 - 2r)}} $$

This expression can be manipulated into one involving central binomial coefficients. An intermediate form is

$$ Q_{i,j}^{ - 1} = {{j(2n - 2j + 1)} \over {2^{4n - 2} }}\left( {\matrix{ {2n - 2j} \cr {n - j} \cr } } \right)\left( {\matrix{ {2j} \cr j \cr } } \right)\sum\limits_{r = 0}^{i - 1} {{\left( {\matrix{ {2r} \cr r \cr } } \right)\left( {\matrix{ {2n - 2r} \cr {n - r} \cr } } \right)} \over {2j - 2r - 1}} $$

which further simplifies to one with a generalised hypergeometric function 4F ₃.

$$ \eqalign{ Q^{ - 1} = & {{\left( {\matrix{ {2i} \cr i \cr } } \right)\left( {\matrix{ {2n - 2i} \cr {n - i} \cr } } \right)\left( {\matrix{ {2j} \cr j \cr } } \right)\left( {\matrix{ {2n - 2j} \cr {n - j} \cr } } \right)} \over {2^{4n - 2} }} \cr & \times {{j(2n - 2j + 1)} \over {2i - 2j + 1}}_4 F_3 \left( {1,i + {1 \over 2},i - j + {1 \over 2},i - n;i + 1,i - j + {3 \over 2},i - n + {1 \over 2};1} \right) \cr} $$

This expression is broken into two parts to highlight the symmetry. Since both Q ⁻¹ and the left collection of central binomial coefficient terms are symmetrical about the diagonal, so must be the term to the right.

7.0 COMBINATORIAL PROPERTIES OF THE INFLUENCE MATRIX FOR 2D VORTEX LATTICE METHOD

James^{(Reference James10)} describes the analysis of a thin 2D wing section divided into n segments, with each segment having a location of vortex μ, location of sensing point λ, and distance between them Δ. The standard values in the 2D vortex lattice method are μ = 1/4, λ=3/4, and Δ = 1/2. The matrix A in his paper relating vortex strength (and so wing lift) has elements A_i,j = 1/(i − j + 1/2), which is just a scaling of H introduced above, so that A = 2H. A has an inverse with the unique property, previously identified by Collar^{(Reference Collar8)}, that the sum of elements is

$$ {\bf e}^T A^{ - 1} {\bf e} = {1 \over 2}n $$

or writing in terms of H,

$$ {\bf e}^T H^{ - 1} {\bf e} = n. $$

8.0 THE ELEMENTS OF THE INVERSE OF A

James and Collar both provide expressions, differing from each other, for the individual elements of A ⁻¹. The expression from James is

$$ A_{i,j}^{ - 1} = {{( - 1)^{n + i + j} \left( {\matrix{ {n - 1} \cr {i - 1} \cr } } \right)\left( {\matrix{ {n - 1} \cr {j - 1} \cr } } \right)\left( {\prod\limits_{r = 1}^n (n - i + 1) - r + {1 \over 2}} \right)\left( {\prod\limits_{r = 1}^n j - r + {1 \over 2}} \right)} \over {{1 \over 2} + j - i}}. $$

As per Q, the expression can be manipulated into one involving central binomial coefficients. The values for the elements of A ⁻¹ are

$$ A_{i,j}^{ - 1} = {{2^{2 - 4n} ij(2n - 2i + 1)} \over {(2i - 1)(2j - 2i + 1)}}\left( {\matrix{ {2n - 2i} \cr {n - i} \cr } } \right)\left( {\matrix{ {2i} \cr i \cr } } \right)\left( {\matrix{ {2n - 2j} \cr {n - j} \cr } } \right)\left( {\matrix{ {2j} \cr j \cr } } \right). $$

Using the same approximation as in Equation (10), the individual terms of A ⁻¹ can be approximated as

$$ A_{i,j}^{ - 1} \simeq {{16j} \over {\pi ^2 (2i - 1)(2j - 2i + 1)}}\sqrt {{{i(2n - 2i + 1)} \over {j(2n - 2j + 1)}}} . $$

9.0 THE ROW SUMS OF THE INVERSE OF A

Similar to the row sums of Q ⁻¹, the distribution of the row sums of A ⁻¹ or H ⁻¹ mirrors the distribution of lift for the 2D constant unity downwash case. Collar^{(Reference Collar8)} derived an expression for the row sums of what would be A ⁻¹, defining y(k) as the kth element of the row sum vector y,

$$ y(k) = {{( - 1)^{k - 1} \prod\limits_{j = 0}^{n - 1} \left( { - j + n - k + {1 \over 2}} \right)} \over {(k - 1)!(n - k)!}}. $$

This can be manipulated into a form using central binomial coefficients. The vector of row sums of a = A ⁻¹e, is

(11) $$ {\bf a}_k = 2^{ - 2n} k{{2n - 2k + 1} \over {2k - 1}}\left( {\matrix{ {2n - 2k} \cr {n - k} \cr } } \right)\left( {\matrix{ {2k} \cr k \cr } } \right) $$

Using the same approximations as in Equation (10) for the central binomial terms gives

(12) $$a_k \simeq {2 \over \pi }{{\sqrt {k(n - k + 1/2)} } \over {2k - 1}}$$

Letting k = nz, z ∈ (0, 1] and taking the limit as n → ∞ of Equation (12) gives

$$ {\bf a}_z = {1 \over \pi }\sqrt {{{1 - z} \over z}} . $$

This corresponds to the exact solution for the 2D flow with constant downwash (Fig. 3).

Figure 3. Comparing exact solution (solid line) to row sums (points) for inverse of A.

The ratio between the lifting line and the 2D vortex lattice exact solutions for constant downwash is

$$ 4{{\sqrt {z(1 - z)} } \over {\sqrt {{{1 - z} \over z}} }} = 4z. $$

The Row Sums of the Inverse of A has a similarly clean relation to the row sums of the inverse of Q. Dividing Equation (8) by Equation (11) gives the exact ratio of the elements of the row sum of Q ⁻¹ to the elements of the row sums of A ⁻¹,

$${\bf{r}}_k = (4k - 2){\bf{a}}_k.$$

10.0 MODEL CONVERGENCE FOR THE DISCRETE LIFTING LINE MODEL

In Collar^{(Reference Collar8)}, he compared the discrete model of the wing with the ideal closed form solution available for elliptical wings, determining the rate of convergence of the discrete model to the ideal. The analysis here follows Collar’s but is greatly simplified.

Consider an elliptical wing with a peak vortex filament strength at its centre of Γ₀. The downwash everywhere along the wing is w = Γ₀/2. Recalling that the lift L = ρVΓ and assuming that ρ = V = 1, the lift of the wing is L = (π/4)Γ₀.

For the discretised model of the wing with n segments, define a vector of constant downwash w = (Γ₀/2)e, and from Equation (3),

(13) $$ \eqalign{ {\bf{\Gamma }} = & \frac{\pi } {n}Q^{ - 1} {\bf{w}} \cr = & \frac{\pi } {n}\frac{{\Gamma _0 }} {2}Q^{ - 1} {\bf{e}} .\cr} $$

The total lift generated by the system of vortices is L _n = (1/n)e^TΓ, so

(14) $$L_n = \frac{\pi } {{2n^2 }}\Gamma _0 {\bf{e}}^T Q^{ - 1} {\bf{e}}$$

Taking the ratio of R_n = L_n/L gives

(15) $$ R_n = \frac{{L_n }} {L} = \frac{2} {{n^2 }}{\bf{e}}^T Q^{ - 1} {\bf{e}} $$

In Collar’s paper, a significant combinatorial manipulation follows to gain the final result. Here, since we know e^TQ ⁻¹e = n(n + 1)/2, this immediately resolves to

(16) $$R_n = \frac{2} {{n^2 }}\left( {\frac{{n(n + 1)}} {2}} \right) = 1 + \frac{1} {n}$$

As noted in the introduction, the triangular number property of the Q ⁻¹ matrix is what leads to the exact, clean form of convergence for the constant downwash wing. This rate of convergence likely generalises to all wing planforms.

11.0 ALTERNATE PROOF OF MINIMUM INDUCED DRAG OF THE CONSTANT DOWNWASH CONDITION

It is well-known that an elliptical lift distribution results in minimum induced drag along with a resulting constant downwash along the wing^{(Reference Munk17, Reference Jones18)}. This section provides an alternate treatment using quadratic programming to gain a direct solution for Γ and w. Again, assuming that ρ = V = 1, the induced drag and lift are

$$\eqalign{ D = & \frac{1} {n}{\bf{w}}^T {\bf{\Gamma }} = \frac{1} {\pi }{\bf{\Gamma }}^T Q{\bf{\Gamma }} \cr L = & \frac{1} {n}{\bf e}^T {\bf{\Gamma }}. \cr}$$

The goal is to minimise the induced drag while meeting a required lift constraint. Casting as a quadratic programming problem gives

$$\eqalign{ & \mathop {{\text{minimize}}}\limits_{\bf \Gamma} \,\frac{1} {2}{\bf{\Gamma }}^T Q{\bf{\Gamma }} \cr & {\text{subject}}\,{\text{to}}\,{\bf{e}}^T {\bf{\Gamma }} = nL. \cr}$$

with the 1/2 in the objective function used to match the standard form for quadratic programming. Since Q is symmetric positive definite and we have equality constraints, we can use the method of Lagrange multipliers^{(Reference Fletcher19)} to immediately gain the solution by solving the linear problem

(17) $$\left[ {\matrix{ Q & {\bf e} \cr {{{\bf e}^T}} & 0 \cr } } \right]\left[ {\matrix{ {\bf \Gamma} \cr \lambda \cr } } \right] = \left[ {\matrix{ 0 \cr {nL} \cr } } \right]$$

where λ is the Lagrange multiplier. Solving for Γ gives

$${\bf{\Gamma }} = Q^{ - 1} {\bf{e}}({\bf{e}}^T Q^{ - 1} {\bf{e}})^{ - 1} nL.$$

Substituting for e^TQ ⁻¹e and noting that w = (n/π)Q Γ, we show that minimum drag implies a constant downwash

$${\bf{w}} = \frac{n} {{n + 1}}\frac{{2L}} {\pi }{\bf{e}}.$$

Note that the proof of optimality of constant downwash only depends on the symmetric positive definite nature of Q and not particularly on the Biot-Savart law or the elliptical distribution of lift. In this sense, the constant downwash is the more fundamental property.

12.0 DERIVING WING PROPERTIES FOR THE LIFTING LINE FORMULATION

This section presents the lifting line approach that was under development when the discovery about the combinatorial properties of Q were made. It is easy to implement in modern mathematical software with direct support for linear algebra and offers some interesting insights to the problem. The method calculates the coefficients of lift and induced drag for the 3D wing from the properties of the 2D aerofoil sections and the wing planform shape. To simplify presentation, assume the wing has a constant 2D aerofoil shape along the span. The assumption is easily relaxed. Choose units of measurement such that the wing has span b = 1 and is operating at some density ρ and velocity V. We will show that for calculating C_L and C_Di the velocity and density terms will cancel.

Referring back to Fig. 1, for each segment the Kutta-Joukowski theorem defines the 2D lift per unit length as

$$L_i = \rho V\Gamma _i$$

Lift is also described as

$$L_i = \frac{1} {2}\rho V^2 c_i c_{li} = \frac{1} {2}\rho V^2 c_i a\alpha L_i$$

where c_i is the chord, c_li the 2D coefficient of lift for the wing section, a is the lift slope curve (assumed constant along the wing), and α_Li is the local angle-of-attack. Equating the two and cancelling common terms gives

$$\Gamma _i = \frac{1} {2}Vc_i a\alpha _{Li}.$$

Define a vector of chord values c, a corresponding diagonal matrix C, and a vector of local angles of attack α _L:

$$\eqalign{ {\bf{c}} = & [c_1 \cdots c_N ]^T \cr C = & diag({\bf{c}}) \cr {\bf{\alpha }}_L = & [\alpha _{L1} \cdots \alpha _{LN} ]^T \cr}$$

Write in vector form

(18) $${\bf{\Gamma }} = \frac{1} {2}VaC{\bf{\alpha }}_L$$

For each wing segment, the angle-of-attack is α_Li = α_i − (w_i/V), with α_i the angle-of-attack with respect to the freestream velocity (and contains a twist component if there is one) and w_i the downwash induced by the trailing vortices. Write in vector form

(19) $${\bf{\alpha }}_L = {\bf{\alpha }} - \frac{w} {V}$$

Using this and Equation (3) with b = 1, rewrite Equation (18) as

(20) $$\eqalign{ {\bf{\Gamma }} = & \frac{1} {2}VaC({\bf{\alpha }} - \frac{w} {V}) \cr = & \frac{1} {2}VaC{\bf{\alpha }} - C\frac{{an}} {{2\pi }}Q{\bf{\Gamma }} \cr}$$

To clean up the presentation, define the influence matrix

(21) $$W = \frac{{an}} {{2\pi }}Q$$

Solve for Γ using Equations (20) and (21):

(22) $${\bf{\Gamma }} = \frac{1} {2}V(C^{ - 1} + W)^{ - 1} a{\bf{\alpha }}.$$

The lift is L = ρV(1/n)e^TΓ, so

(23) $$L = \frac{1} {{2n}}\rho V^2 {\bf{e}}^T (C^{ - 1} + W)^{ - 1} a{\bf{\alpha }}$$

By definition, c_L = 2L/ρV ²S. The surface area S is S = e^TC e/n. Using Equation (23), the closed form expression for c_L can be written as a function of α and the geometry of the wing as captured in the matrix C,

(24) $$c_L = \frac{{{\bf{e}}^T (C^{ - 1} + W)^{ - 1} }} {{{\bf{e}}^T C{\bf{e}}}}a{\bf{\alpha }}$$

or alternatively

(25) $$c_L = \frac{{{\bf{e}}^T C(I + WC)^{ - 1} }} {{{\bf{e}}^T C{\bf{e}}}}a{\bf{\alpha }}$$

Note the cancellation of ρ and V, as expected. For the remainder of the paper, let ρ = V = 1 to simplify presentation.

The induced drag is D = (1/n)w^T Γ. With substitutions,

$$D = \frac{a} {{2n}}{\bf{\alpha }}^T (C^{ - 1} + W)^{ - 1} W(C^{ - 1} = W)^{ - 1} {\bf{\alpha }}.$$

The induced drag coefficient is c_Di = 2D/S, so

(26) $$c_{Di} = \frac{a} {n}\frac{{{\bf \alpha} ^T \left( {C^{ - 1} + W} \right)^{ - 1} W\left( {C^{ - 1} + W} \right)^{ - 1} {\bf \alpha} }} {{{\bf e}^T C{\bf e}}}$$

Noting that Oswald’s coefficient satisfies

$$C_{Di} = \frac{{c_L^2 }} {{\pi ARe}}$$

Using Equation (24) and (26), we can write

(27) $$e = \frac{a} {{\pi AR}}\frac{{{\bf \alpha} ^T \left( {C^{ - 1} + W} \right)^{ - 1} {\bf ee}^T \left( {C^{ - 1} + W} \right)^{ - 1} {\bf \alpha} }} {{\alpha ^T \left( {C^{ - 1} + W} \right)^{ - 1} W\left( {C^{ - 1} + W} \right)^{ - 1} {\bf \alpha} }}.$$

13.0 PROPERTIES AND CONVERGENCE OF CONSTANT DOWNWASH CASE – WING WITH NO TWIST

Assuming that our wing has no twist, we have α = αe for some scalar α. Constant downwash results in a distribution of Γ as per Equation (13), which we can insert into Equation (22) to get

$$Q^{ - 1} {\bf e} \propto \left( {C^{ - 1} + \frac{{an}} {{2\pi }}Q} \right)^{ - 1} {\bf e}.$$

This can be manipulated to give

$$C{\bf e} \propto Q^{ - 1} {\bf e}$$

so, the distribution of chords tends to the elliptical as expected. Also, since for some scalar β we have C e = βW ⁻¹e or WC e = β e, WC has constant row sums, and e is an eigenvector of WC with β the associated eigenvalue. This implies

$$\matrix{ {\left( {I + WC} \right){\bf e} = \left( {1 + \beta } \right){\bf e}} \hfill \cr {{{\left( {I + WC} \right)}^{ - 1}}{\bf e} = {1 \over {1 + \beta }}{\bf e}.} \hfill \cr } $$

Substituting into Equation (25) with α = α e and writing in terms of a constant c_l along the wing such that c_l = aα gives

(28) $$c_L = \frac{1} {{1 + \beta }}c_l.$$

The aspect ratio AR of a wing is b ²/S, which for this formulation is n/e^TC e. If the distribution of chords along the wing is defined as C e = βW ⁻¹e, then for such a wing the aspect ratio AR is

$${\text{AR}} = \frac{n} {{\beta \;{\bf e}^T W^{ - 1} {\bf e}}} = \frac{{2an^2 }} {{\beta \;2\pi \;n(n + 1)}} = \frac{{an}} {{\beta \pi (n + 1)}}.$$

Solving for β gives

(29) $$\beta = \frac{{an}} {{{\text{AR}}\pi (n + 1)}}$$

Using Equations (28) and (29) and writing c_L as c_L(n) to highlight the dependency, the expression for c_L(n) is

$$c_L (n) = \frac{{\left( {1 + \frac{1} {n}} \right){\text{AR}}}} {{\left( {1 + \frac{1} {n}} \right){\text{AR}} + \frac{a} {\pi }}}c_l$$

which as n grows rapidly tends to

(30) $$c_L (n) = \left( {1 + \frac{1} {n}} \right)\frac{{{\text{AR}}}} {{{\text{AR}} + \frac{a} {\pi }}}c_l$$

The limit as n → ∞ is as expected,

$$\mathop {\lim }\limits_{n \to \infty } c_L (n) = \frac{{AR}} {{AR + \frac{a} {\pi }}}c_l.$$

Using Equation (26) with substitutions, C_Di is

(31) $$c_{Di} = \frac{{n(n + 1)\pi AR}} {{\left( {(n + 1)\pi AR + na} \right)^2 }}c_l^2$$

Manipulating Equation (31) and inserting Equation (30) to show the relationship c_L to gives

(32) $$c_{Di} (n) = \left( {1 + \frac{1} {n}} \right)\left( {\frac{{c_L (n)^2 }} {{\pi {\text{AR}}}}} \right)$$

The limit as n → ∞ of Equation (32) is again the well-known

(33) $$c_{Di} = \frac{{c_L^2 }} {{\pi AR}}.$$

From Equations (32) and (33), the estimate of the Oswald coefficient for the elliptical case convergences as

$$e = 1 + \frac{1} {n}.$$

14.0 EXPERIMENTAL RESULTS ON MODEL CONVERGENCE

The convergence of the constant downwash model as per 1/n suggests that this property may exist for other planforms and can be exploited to gain accurate results from lower accuracy models. Collar contemplated this in his 1958 paper, and computational experiments suggest merit to his claim. The four planforms in Fig. 4 served as test cases for computing c_L, using Equation (24) to generate the approximations with n varying from 200 to 1,200 in steps of 100.

Figure 4. Wing planform shapes for convergence test.

The results show that for a variety of wing planforms, the convergence is by 1/n (see Fig. 5). The same behaviour is observed for computing c_Di.

Figure 5. Plot of calculated c_L versus 1/n for 4 different planforms. Aspect Ratio for each in parentheses.

Experiments show that with the 1/n convergence, when one takes a finite difference about a value of n in Fig. 5 and projects a line to the Y-axis, significant gains in accuracy are possible with only a minimal increase in computation. For example, using the dart planform from above, by doing a finite differencing of the model about n=100 (e.g., run it at n=98 and n=102) and determining the y-intercept of the line, one gets an estimate of c_L that is as accurate as the n=1,000 model but at a factor of 50 reduction in computation effort (assuming O(n ²) computational complexity). Experiments also suggest that the values determined via direct computation and those from finite differencing provides bounds above and below, respectively, for the estimate of c_L (Fig. 6).

Figure 6. Demonstration of using finite differencing to provide estimates of values. Finite difference samples were at n+2 and n-2 around the baseline n.

It would be an interesting problem to formally prove the convergence behaviour seen in Fig. 5 for all wing shapes (or determine the limits of the behaviour) and to confirm that the property of the bounding from above and below.

15.0 PROOF OF COMBINATORIAL EXPRESSIONS RELATED TO Q

Two results are proved here, the expressions for the sum of the elements of Q ⁻¹ and for the row sums of Q ⁻¹. The proofs for properties of Q and H parallel each other, so the methods used also directly apply to H. Both identities for Q were initially discovered via experimentation, generating data sets varying in dimensionality and identifying patterns in them, using mathematical symbolic software as an aid. Once having identified the answer, the problem became one of how to prove the answer in a way that was verifiable and acceptable. The proof method is based on relatively recent (since mid-1990s) research on automated proofs of sums of hypergeometric series^{(Reference Koepf20)}. Specifically, the proof uses Zeilberger’s Algorithm ^{(Reference Petkovsek, Wilf and Zeilberger21)} as implemented in the fastZeil Footnote ² Mathematica package^{(Reference Paule and Schorn22)}. It works for problems with a summand f(n, k) that is hypergeometric in both indices, meaning both of the ratios of f(n + 1, k)/f(n, k) and f(n, k + 1)/f(n, k) result in rational functions in n and k. While the algorithm’s internal workings are opaque, it generates results that are readily verifiable. The reader is steered to the references for details. The Mathematica notebooks executing the proofs are available as supplementary materials to this paper, so what follows is a summary of the proof methods.

The first identity to be proved is for the row sums. Using the identities, QQ ⁻¹e = e and r = Q ⁻¹e gives the expression to be proven, Q r = e, which when converted into a summation states

(34) $$\forall n\;and\;\forall i \in [1,n];\;\sum\limits_{k = 1}^n {{1 \over {1 - 4{{(i - k)}^2}}}} {{k(2n - 2k + 1)} \over {{2^{2n - 1}}}}\left( {\matrix{ {2k} \cr k \cr } } \right)\left( {\matrix{ {2n - 2k{\rm{ }}} \cr {n - k} \cr } } \right) = 1$$

This proves that the expression for the row sums of Q ⁻¹ is correct. The second is that

(35) $$\forall n;\;\sum\limits_{k = 1}^n {{{k(2n - 2k + 1)} \over {{2^{2n - 1}}}}} \left( {\matrix{ {2k} \cr k \cr } } \right)\left( {\matrix{ {2n - 2k} \cr {n - k} \cr } } \right){2 \over {n(n + 1)}} = 1$$

This proves that the sum of the row sums, i.e. the sum of the elements of Q ⁻¹, is the Triangular Number T(n). Start with the proof of the sum of the elements of Q ⁻¹ as it is simpler and is used to explain the method, but be aware that it depends on the correctness of proof of the row sums terms of Q ⁻¹. Define the term to be summed:

(36) $$f(n,k) = {{k(2n - 2k + 1)} \over {{2^{2n - 1}}}}\left( {\matrix{ {2k} \cr k \cr } } \right)\left( {\matrix{ {2n - 2k} \cr {n - k} \cr } } \right){\rm{ }}{2 \over {n(n + 1)}}$$

The goal is to verify

(37) $$\forall n \in \{ 1,2, \cdots ,\infty \} ;\;\sum\limits_{k = 1}^n f(n,k) = {\text{SUM}}(n) = 1$$

Running the problem against the Zeilberger algorithm gives the recurrence

(38) $${\text{SUM}}(n + 1) - {\text{SUM}}(n) = 0$$

or alternatively

(39) $${\text{SUM}}(n + 1) = {\text{SUM}}(n)$$

Checking SUM(1) = f(1, 1) = 1 completes the induction showing that for all n, SUM(n) = 1.

To verify the recurrence generated by the solving algorithm, the algorithm also generates a proof certificate of

(40) $$r(n,k) = \frac{{(k - 1)( - 2k + 2n + 3)}} {{2(n + 2)( - k + n + 1)}}$$

The term r(n, k) when multiplied with f(n, k)gives a new term g(n, k) = r(n, k)f(n, k). The algorithm also provides a relationship that is manually verifiable,

(41) $$f(n,k) - f(n + 1,k) = g(n,k + 1) - g(n,k)$$

This defines a Wilf-Zeilberger Pair. The left and right sides of Equation (41) are easily confirmed to be equivalent (see supplementary material for this step). Now, sum both sides of Equation (41) from k = 0 to to k → ∞ to get

(42) $$\sum\limits_{k = 0}^{n + 1} f(n + 1,k) - \sum\limits_{k = 0}^n f(n,k) = \mathop {\lim }\limits_{k \to \infty } g(n,k) - g(n,0)$$

Due to the binomial coefficients in f(n, k) and g(n, k), the terms are nonzero only for a finite range of values of k, specifically 0 ≤ k ≤ n, so the summations on the left side are over finite ranges and are the sums in Equation (38). On the right side, the summations telescope into two simple terms. Checking that g(n, 0) = 0 and ∀k > n; g(n, k) = 0 shows that Equation (38) holds and thus completes the proof.

Proving the expression related to the row sums is more complex as there are two free indices, i and k. Define

(43) $$f(n,k) = {1 \over {1 - 4{{(i - k)}^2}}}{{k(2n - 2k + 1)} \over {{2^{2n - 1}}}}\left( {\matrix{ {2k} \cr k \cr } } \right)\left( {\matrix{ {2n - 2k} \cr {n - k} \cr } } \right)$$

Running the Zeilberger algorithm leads to the recurrence

(44) $$\eqalign{ & \left( {n + 1} \right)\left( {2i - 2n - 5} \right){\text{SUM}}(n + 2{\rm{)}} + \left( { - 4in - 6i + 4{n^2} + 13n + 9} \right){\text{SUM}}\left( {n + 1} \right)\; \cr & \quad + 2\left( {n + 2} \right)\left( {i - n - 1} \right){\text{SUM}}\left( n \right) = 0 \cr} $$

Verifying that SUM(1) = 1 for i=1, and SUM(2) = 1 for i=1 and 2, and checking the recurrence in Equation (44) completes the induction for all n. The proof certificate is

(45) $$r(n,k) = \frac{{\left( {k - 1} \right)\left( { - 2i + 2k + 1} \right)\left( { - 2k + 2n + 3} \right)}} {{4\left( { - k + n + 1} \right)\left( { - k + n + 2} \right)}}$$

Noting that g(n, k) = r(n, k)f(n, k) takes values g(n, 0) = g(n, n + 3) = 0 for the bounds, and confirming

(46) $$g(n,k + 1) - g(n,k) = \left({\matrix {2\left( {n + 2} \right)\left( {i - n - 1} \right)f(n,k)} \\ { + \left( { - 4in - 6i + 4n^2 + 13n + 9} \right)f\left( {n + 1,k} \right)} \\ { + \left( {n + 1} \right)\left( {2i - 2n - 5} \right)f(n + 2,k)}\matrix}\right)$$

completes the proof. The supplementary material contains the proofs for the properties related to H.