1 Introduction

Multiplane interest has been growing up in the last 20 years, as the large variety of scientific contributions show [1, 3, 5, 6, 10, 16, 17, 19]. The renewed interest in such an aircraft architecture has been boosted by the new challenges the aeronautical world is facing, because of increments in the worldwide air traffic and of stricter regulations about noise and pollution. Hence, the attention has been focusing on new aircraft architectures in order to fulfil requirements that seem impossible to be satisfied with conventional ones. In this context, the box-wing is one of the most studied aircraft architectures [4].

Multiplane theory studies go back to the 20s of the last century, in particular with respect to the minimisation of the induced drag [22, 26]. A seminal work in this sense is the 1924 Prandtl NACA Technical Note 182 [32], English translation of the original [27] of few years before. According to Prandtl, the optimal induced drag arrangement of a multiplane is the limit multiplane ideally having infinitely many wings. As a result, the Best Wing System (BWS) is a box-wing whose vertical wings are equivalent to the tip vortices of the virtual internal wings, whilst the horizontal ones generate the proper lift. Nevertheless, Prandtl did not describe the method he used to derive the BWS performance curve, as the gap between the horizontal segments of the box-wing may vary, nor he showed the resulting optimal lift distribution. Nowadays, modern computational aerodynamics shows the exactness of Prandtl prevision [2, 10, 14] and provides numerical results on the optimal lift distribution profile (see for instance [13]).

Recently, Frediani and Montanari [16] have proposed a closed-form solution to Prandtl problem. According to them, the optimal lift distribution is the superposition of a constant and of an elliptical distribution for the horizontal wings, and it is a butterfly-shaped one for the vertical wings. Their solution is reasonably close to the actual solution [2, 14] for gaps of aeronautical interest. However, for higher values of the gap, discrepancies increase.

In the following years, Demasi and co-workers have intensely studied the problem of the induced drag of multiplanes, with a focus on non-conventional configurations such as the box-wing.

They achieved many interesting aerodynamics theoretical results, in wide generality. In [14], the authors numerically show some primary equivalence properties of closed lifting systems, within a variational numerical framework. For instance, a property that will be useful in the rest of the paper is the following: a closed lifting system presents the same optimal induced drag of a multiwing having infinite number of wings (not necessarily planar), whose tips lies on the curve defining the closed system. It is noteworthy that in [26, Sec. 46, pag. 309] essentially the same claim was used to give a qualitative—yet conceptually rigorous—answer on the shape of the minimum induced drag lifting system for assigned lift and external contour:

Which is the lifting system that, for given external contour, has the minimum induced drag for a given lift? [...] Given a multiwing lifting system, if the annular lifting system that connects the extremal points of the provided wings with segments is considered, the resulting system (annular + wings) has a lower minimum induced drag than the original set of wings, for a given generated lift. However, since internal wings can be removed without any modifications of the induced drag, it follows that every annular lifting system has a minimum induced drag lower than every other lifting system embedded in its interior part, having elements which connect two points of the annular contour. The problem is hence solved.Footnote 1

In [12, 13], the authors show how the optimal lift distribution for a symmetric box-wing is not elliptic-plus-constant. They further investigate the behaviour of the box-wing as the vertical gap increases. It turns out that the lift distribution tends to be constant on both wings. In the limit case of indefinitely large vertical gap, the optimal induced drag tends to zero.

In [14], the authors show, among other properties, that in an optimal multiwing generating positive lift, each wing contributes with positive lift. Therefore, the minimum of the induced drag is achieved when the lifting elements generate positive lift.

This fact was observed in the ‘20s for the multiplane case. In fact, from Munk’s Theorem of stagger, it can be inferred that in an optimal multiplane the induced vertical velocities at each wing must be equal and constant along the spans [22, 23]. Therefore, all the induced velocities must have the same direction as well as the same intensity. In particular, if we consider a global positive lift, all the induced velocities must be directed negatively, which implies that positive lift is generated in every wing.

Moreover, in [14], the authors show that in a symmetric multiwing, the circulations are conjugate-symmetric distributed, i.e. pairs of wings symmetrically located with respect to the symmetry horizontal plane have the same circulation distribution.

Recently, many research projects have focused on box-wing applications. For example, Project IDINTOS (IDrovolante INnovativo TOScano) [21], conceived, developed and built a prototype of a box-wing (PrandtlPlane) seaplane. The design was supported by numerical simulations, wing tunnel and towing tank tests [7, 18].

The box-wing configuration is also the subject of the PARSIFAL Project (PrandtlPlane ARchitecture for the Sustainable Improvement of Future AirpLanes) [8, 9, 19], which aims at studying the applicability of the PrandtlPlane in the civil transportation market of the next years [1, 19]. For further details, see [24].

However, after about a century, it remains to understand how Prandtl achieved his results. In fact, Prandtl just reported the results, neglecting to detail the method he used. A philological study of Prandtl original notes, supposing they are preserved, could provide an unequivocal answer. In this paper, by analysing Prandtl words of the report and of other coeval sources, we attempt at recovering the method he used to evaluate the BWS performance in [27, 32]. However, some questions remain open. In particular, we will see that he used such an approach based on the potential flow only for the BWS, whilst in the same work, for biplane and triplane, he directly evaluated the optimal induced drag (say, he used a direct method). Moreover, he did not provide the optimal resulting lift distribution. In the light of these facts, three questions naturally arise, that we attempt to answer throughout the paper:

  1. 1.

    Why did Prandtl change method when passing from biplane/triplane to the BWS?

  2. 2.

    What would have happened if he continued with the direct method, the number of wing increasing?

  3. 3.

    Which is the BWS lift distribution? Why did not he provide it?

In the second part of the paper, by collecting some results that were likely to be known at Prandtl’s time, we formalise a simple model which is able to reproduce his results, with a surprising accuracy. Since many Prandtl’s works have not been translated from German, the book Aerodinamica [26], published in 1932 but based on the results achieved in the previous years, of the renowned aerodynamicist Enrico Pistolesi, coeval of Prandtl, realistically contains the basic concepts that a 1920s aerodynamicist may have known. In particular, three chapters (namely the 45, 46 and 47) are dedicated to the study of multiplanes, with a focus on the minimum induced drag problem. Indeed, Pistolesi was an attentive connoisseur of Prandtl work, as the critical review [20] confirms. Then, a fourth research question arises:

  1. 4.

    Could Prandtl extend the direct method to formalise a model similar to the one we propose?

The paper is structured as follows. Section 2 introduces the main aspect of the lifting line theory, a framework that will be assumed throughout the paper. Section 3 clarifies Prandtl approach for the Best Wing System. Section 4 introduces some preliminary concepts on the coefficients of self and mutual induction upon different hypotheses on the lift distributions. Section 5 develops the simple model for the evaluation of the induced drag of multiplanes. Section 6 presents the numerical results of our model, whilst Sect. 7 contains a critical discussion of the results and some remarks, try to answer to the four aforementioned research questions.

2 Generalities on the Lifting Line Theory: The Biplane

Hereafter, we assume the lifting line theory and planar wings of equal span. Without entering in the merit of the theory, we recall that the lift, for a single lifting segment, called wing, is expressed by

$$\begin{aligned} L = \rho V_\infty \int \limits _{-b/2}^{b/2} \Gamma (y) \ \text {d}y, \end{aligned}$$
(1)

where \(\Gamma (y)\) is the circulation function, \(\rho\) is the air density, \(V_\infty\) is the asymptotic flow velocity intensity, b is the wingspan, y is the coordinate in span direction. Due to the finiteness of the wingspan, a drag component, named induced drag, appears, whose expression is (see, for example, [35])

$$\begin{aligned} D = \rho \int \limits _{-b/2}^{b/2} \Gamma (y) \ w(y) \ \text {d}y, \end{aligned}$$
(2)

where w(y) is the vertical component of velocity, called downwash, induced by the vorticity elements. It is noteworthy that the downwash expression depends on the lift distribution itself.

With reference to Fig. 1, we consider a biplane whose wings, labelled 1 and 2, have a common span b and a vertical distance (gap) denoted by h. It is convenient to introduce the dimensionless geometrical parameter gap-to-span ratio

$$\begin{aligned} k:=\dfrac{h}{b}. \end{aligned}$$
(3)

Since the induced drag is invariant if lifting elements are translated parallel to the asymptotic flow direction (x direction in Fig. 1) (Munk’s Theorem of stagger [23]), we consider the wings lying on the \(y-z\) plane. For each wing, we can apply Eqs. (1) and (2) by posing, with a clear meaning of notation,

$$\begin{aligned} L_i = \rho V_\infty \int \limits _{-b/2}^{b/2} \Gamma _i(y_i) \ \text {d}y_i,\quad i\in \{1,2\}, \end{aligned}$$
(4)

and

$$\begin{aligned} D_i = \rho \int \limits _{-b/2}^{b/2} \Gamma _i(y_i) \ w_i(y_i) \ \text {d}y_i,\quad i\in \{1,2\}. \end{aligned}$$
(5)

For each wing of the biplane, the induced velocity \(w_i(y_i)\) may be decomposed into the self-induced downwash, namely \(w_{ii}(y_i)\), and the one induced by the other wing’s vorticity, namely \(w_{ij}(y_i)\). As a consequence, Eq. (5) becomes

$$\begin{aligned} \begin{aligned} D_i&= \rho \int \limits _{-b/2}^{b/2} \Gamma _i(y_i) \ w_{ii}(y_i) {\text {d}}y_i + \rho \int \limits _{-b/2}^{b/2} \Gamma _i(y_i) \ w_{ij}(y_i) {\text {d}}y_i, \\&=: D_{ii} + D_{ij} \quad i,j\in \{1,2\}, \ i\ne j. \end{aligned} \end{aligned}$$
(6)

The expressions of \(w_{ii}\) and \(w_{ij}\) are given by the following expressions (see, for example, [35]):

$$\begin{aligned} w_{ii}(y_i)= & {} \dfrac{1}{4\pi }\mathop {\int \!\!\!\!\!\!=}\limits _{-b/2}^{b/2}\dfrac{\Gamma _i(y_i)}{(t-y_i)^2} \text {d}t, \end{aligned}$$
(7)
$$\begin{aligned} w_{ij}(y_i)= & {} \dfrac{1}{4\pi }\int \limits _{-b/2}^{b/2} \Gamma _j(y_j)\dfrac{h^2 - (y_j-y_i)^2}{[h^2 + (y_j-y_i)^2]^2} \text {d}y_j, \end{aligned}$$
(8)

where \(\mathop {\int \!\!\!\!\!\!=}\) denotes the integration in the Hadamard finite-part sense.

The total lift and total induced drag for the biplane are then,

$$\begin{aligned} L = L_1 + L_2 \end{aligned}$$
(9)

and

$$\begin{aligned} D = D_1 + D_2 = D_{11}+D_{22} + 2D_{12}, \end{aligned}$$
(10)

since it can be shown that, under our assumptions, \(D_{ij}=D_{ji}\) (\(i\ne j\)) [32].

Fig. 1
figure 1

Conceptual sketch for the evaluation of the induced downwash for a biplane

Full size image

More explicitly, substituting Eq. (7) into Eq. (6), \(D_{ii} \ (i\in \{1,2\})\) terms read

$$\begin{aligned} D_{ii} = \dfrac{\rho }{4\pi } \int \limits _{-b/2}^{b/2}\mathop {\int \!\!\!\!\!\!=}\limits _{-b/2}^{b/2}\Gamma _i(y_1) \Gamma _i(y_2) \dfrac{1}{(y_2-y_1)^2} {\text {d}}y_1 {\text {d}}y_2. \end{aligned}$$
(11)

Substituting Eq. (8) into Eq. (6), the coupling term \(D_{12}\) reads [35]

$$\begin{aligned} D_{12}= & {} \dfrac{\rho }{4\pi } \int \limits _{-b/2}^{b/2}\int \limits _{-b/2}^{b/2}\Gamma _1(y_1)\Gamma _2(y_2) \nonumber \\&\times \dfrac{h^2 - (y_2-y_1)^2}{\left[ h^2+(y_2-y_1)^2\right] ^2} \ \text {d}y_1 \ \text {d}y_2. \end{aligned}$$
(12)

Condensing the contributes that depend on the circulation shapes in coefficients \(\sigma _{ii}\) and \(\sigma _{ij}\), the total induced drag for a biplane can be formally expressed as [26])

$$\begin{aligned} D = \dfrac{1}{\pi q b^2}\left( \sigma _{11}L_1^2+\sigma _{22}L_2^2+2L_1L_2\sigma _{12}\right) , \end{aligned}$$
(13)

where \(\sigma _{ii}\), \(i\in \{1,2\}\), are the self-induction coefficients and \(\sigma _{12}\) is the mutual induction coefficient.

The case hereby exposed for the biplane can be easily generalised for multiplanes, as we will see in the following.

3 Prandtl Method for the Best Wing System

In [27, 32], Prandtl evaluated the optimal biplane and triplane performances, in terms of induced drag, adopting elliptic circulations on all the wings, claiming that it is a more refined and conceptually rigorous approach than the method of Albert Betz. The latter corresponds to the evaluation of the mutual induced drag coefficients considering constant circulations. Of course, constant circulations are non-physical, since violate the zero-circulation boundary condition at wingtips. As a remarkable fact, despite the promising prediction power of assuming all elliptic circulations on the wings (the elliptic distribution is the well-known minimiser of the induced drag for the monoplane), Prandtl changed method to obtain the BWS performance. Accordingly, he claimed:

The induced drag of this wing system [the BWS] has been evaluated according to a method which I cannot here be gone into in detail.

and

For the calculation of this drag [of the BWS] I am indebted to Messrs. Grammel and Pohlhausen.

No more details are provided, apart from these laconic words.

In others works of Prandtl [30, 31], he acknowledges again the contribution of Grammel and Pohlhausen:

Dr. Grammel and Pohlhausen have treated, at my instigation, the case of the biplane made up of two straight monoplanes of the same span. [...] The calculations [...] are solved by means of elliptic integrals. I have given the formulas in my Wing Theory II [29].

This time, few more details are provided on the method that, since the repeated acknowledgements, we can claim to be the sought one. We just expose the main assumptions and the main results of this method (for further details, see [26, 28, 29, 31, 34]).

Consider a monoplane and the potential lifting line theory. Under optimality conditions, the downwash velocity is constant. The downstream velocity field is, approximately, a uniplanar flow around the vortex system produced by the lifting system in its motion, and this vortex system may be regarded, as a first approximation, as a solid body in the fluid (so-called rigid-wake). Under these hypotheses, it can be proved that the circulation is equal to

$$\begin{aligned} \Gamma =w\,(\varphi _2-\varphi _1), \end{aligned}$$
(14)

where w is the asymptotic downwash velocity intensity, \(\varphi _2-\varphi _1\) is the jump of the potential (evaluated for unitary downwash) when passing from the dorsal to the ventral side of the lifting line. From Eq. (1), we have

$$\begin{aligned} L=\rho V_\infty w \int \limits _{-b/2}^{b/2} (\varphi _2-\varphi _1) \ \mathrm {d}y, \end{aligned}$$
(15)

or, posing

$$\begin{aligned} \Sigma :=\int \limits _{-b/2}^{b/2} (\varphi _2-\varphi _1) \ \mathrm {d}y, \end{aligned}$$
(16)

we have

$$\begin{aligned} L=\rho V_\infty w \Sigma . \end{aligned}$$
(17)

 Substituting Eq. (14) and the expression for \(w\) from Eq. (17) into Eq. (2), we obtain

$$\begin{aligned} D = \dfrac{L^2}{\rho V_\infty ^2 \Sigma }. \end{aligned}$$
(18)

The key-point is then evaluating \(\Sigma\). It can be done through elliptic integrals, as detailed in [29]. Note that \(\Sigma\) has the dimension of a surface and depends only upon the geometrical properties of the front view of the lifting system. Moreover, the diagram of \(\varphi _2-\varphi _1\), at the due scale, represents the distribution of the lift.

For a general multiplane of N wings, the same considerations hold for each wing. Therefore, one has (under optimality conditions \(w_{i} = w\) for each \(i=1,\dots ,N\))

$$\begin{aligned} L=\rho V_\infty w \sum \limits _{i=1}^{N} \Sigma _i, \end{aligned}$$
(19)

where \(\Sigma _i\) is the quantity of Eq. (16) evaluated for the ith wing.

Figure 2 shows the fictitious surfaces \(\Sigma\) for the optimal monoplane, biplane, triplane and, finally, for the BWS.

Fig. 2
figure 2

Representation of the fictitious surface \(\Sigma\) from [30] (figures are not in scale)

Full size image

For a monoplane, it is easy to see that \(\Sigma = \pi \frac{b^2}{4} =: \Sigma _1\). To compare the performance of optimal multiplanes in terms of induced drag, for given lift and wingspan, it is sufficient to evaluate the value of the fictitious area \(\Sigma\) for the case at hand. This value can then be normalised with the monoplane counterpart, thus providing an immediate comparison of the “efficiency”. For instance, the ratio of the induced drag between the optimal N-wing multiplane and the optimal monoplane reads

$$\begin{aligned} \dfrac{L^2}{\rho V_\infty ^2 \Sigma _{N}} \dfrac{ \rho V_\infty ^2 \Sigma _1}{L^2} = \dfrac{\Sigma _1}{\Sigma _{N}}, \end{aligned}$$
(20)

where \(\Sigma _{N}\) is the fictitious surface for the generic N-wing multiplane.

As a final remark, the induced drag of Eq. (18) is inversely proportional to \(\Sigma\), i.e. the larger \(\Sigma\), the smaller the induced drag (the lift is constant in our problem). For an optimal biplane having infinitely distant wings, \(\Sigma _2\) equals \(2\Sigma _1\),Footnote 2 from which the well-known result that the optimal biplane has half of the induced drag of the optimal monoplane (of same lift and span). Conversely, for the BWS having infinitely distant horizontal wings, \(\Sigma _{BWS}\rightarrow \infty\). As a result, the corresponding induced drag tends to zero, in agreement with the results presented in [12,13,14].

4 Mutually Induced Drag Coefficients Evaluation

In the following part of the paper, we try to answer the four research questions by extending the Prandtl direct method as the number of wings increases. In fact, it is well-known the equivalence between the box-wing and a multiplane having infinite many wings, as it will be discussed more in detail. As a preliminary part to the proposed extension of the direct method, that will be practically developed in Sect. 5, some concepts and calculations must be introduced on the self and mutual induced drag coefficients. In fact, our purpose is to write the formal expression of the induced drag of a generic multiplane, by generalising Eq. (13). Once done, the first-order optimality of the induced drag is evaluated. Let us consider in the following three special circulation shapes.

4.1 Two Elliptic Circulations

It is well-known that an elliptic lift distribution is the minimiser of the induced drag for a monoplane of assigned wingspan and total lift. For a biplane, the optimal lift distribution is not elliptic, as pointed out in [11]. However, as also done by Prandtl for the biplane and triplane, we may assume elliptic distributions on all wings, which read

$$\begin{aligned} \Gamma _i(y_i) = \Gamma _{0i} \sqrt{1-\left( \dfrac{2y_i}{b}\right) ^2} \quad i\in \{1,2\}. \end{aligned}$$
(21)

Substituting Eq. (21) into Eq. (4), we obtain

$$\begin{aligned} L_i = \dfrac{\pi }{4}\rho V_\infty b \Gamma _{0i} \quad i\in \{1,2\} \end{aligned}$$
(22)

and, more expressively,

$$\begin{aligned} \Gamma _i(y_i) = \dfrac{4L_i}{\pi \rho V_\infty b} \sqrt{1-\left( \dfrac{2y_i}{b}\right) ^2} \quad i\in \{1,2\}. \end{aligned}$$
(23)

It can be shown that for elliptic lift distributions, the self-induced downwash is constant [35]

$$\begin{aligned} w_{ii} = \dfrac{\Gamma _{0i}}{2b} = \dfrac{2L_i}{\pi \rho V_\infty b^2} \quad i\in \{1,2\}. \end{aligned}$$
(24)

Substituting Eqs. (23) and (24) into Eq. (2), we obtain the well-known formula

$$\begin{aligned} D_{ii} = \dfrac{L_i^2}{\pi q b^2} \quad i\in \{1,2\}. \end{aligned}$$
(25)

Looking at Eqs. (25) and (13), coefficients \(\sigma _{ii}=1\).

Substituting Eq. (23) into Eq. (12), we obtain

$$\begin{aligned} D_{12} = \dfrac{L_1L_2}{\pi q b^2}\sigma _{12}, \end{aligned}$$
(26)

where (see [35])

$$\begin{aligned} \begin{aligned} \sigma _{12}&:= \!\dfrac{2}{\pi ^2} \int \limits _{-b/2}^{b/2}\int \limits _{-b/2}^{b/2} \prod \limits _{i=1}^2\sqrt{1-\left( \frac{2y_i}{b}\right) ^2} \\&\quad \times \dfrac{h^2-\left( y_2-y_1\right) ^2}{\left[ h^2+\left( y_2-y_1\right) ^2 \right] ^2} \ \mathrm {d}y_1\mathrm {d}y_2 \\&\quad \text {(posing } \eta _i := 2y_i/b, i\in \{1,2\}) \\&= \dfrac{1}{2\pi ^2 k^2} \int \limits _{-1}^{1}\int \limits _{-1}^{1} \prod \limits _{i=1}^2 \sqrt{1-\eta _i^2} \\&\quad \times \dfrac{1-\dfrac{1}{4k^2}\left( \eta _2-\eta _1\right) ^2}{\left[ 1+\dfrac{1}{4k^2}\left( \eta _2-\eta _1\right) ^2 \right] ^2} \ \mathrm {d}\eta _1\mathrm {d}\eta _2 \end{aligned} \end{aligned}$$
(27)

is the mutual influence coefficient of the induced drag of a biplane having elliptic lift distributions on both wings. The integral in Eq. (27) has been solved through numerical integration, for some values of k. Table 1 reports the numerical values compared to Prandtl’s [32].

Table 1 \(\sigma _{12}\) coefficient for different values of k, in the case of elliptic lift distributions
Full size table

Prandtl proposes to approximate the \(\sigma _{12}\) coefficient with the formula

$$\begin{aligned} \sigma _{12} = \dfrac{1-0.66k}{1.05 + 3.7k} \end{aligned}$$
(28)

in the range \(1/15\le k\le 1/2\). However, Eq. (28) does not equal 1 when \(k=0\) (see [12, 13] for more details). Indeed, when two wings, generating a lift force \(L_1 = L_2 = L/2\), tend to coincide, the resulting induced drag can be seen as the induced drag of a single wing generating a lift force equal to L. Considering the values of k up to 0.5 obtained through numerical integration of Eq. (27), the following approximation, obtained via MATLAB® Curve Fitting Toolbox, is proposed (coefficient of determination of the interpolation \(R^2=0.9999\)):

$$\begin{aligned} \sigma _{12} = \dfrac{p_1\,k + p_0}{k^3 + q_2\,k^2 + q_1\,k + q_0}, \end{aligned}$$
(29)

where

$$\begin{aligned} \begin{aligned} p_1&= 0.1919&p_0&=0.007075\\ q_2&= 0.7529&q_1&= 0.2437&q_0&= 0.007075. \end{aligned} \end{aligned}$$

Equation (29) has been assumed in the present work for the numerical part. It can be noted that for \(k=0\), Eq. (29) provides the theoretical value of 1, thus leading to more accurate results for small k.

4.2 Two Constant Circulations

In the case of constant circulations, we assume

$$\begin{aligned} \Gamma _i(y_i) = \Gamma _i, \quad i\in \{1,2\}. \end{aligned}$$
(30)

From Eq. (4), we obtain

$$\begin{aligned} \Gamma _i = \dfrac{L_i}{\rho V_\infty b}, \end{aligned}$$
(31)

and, from Eq. (12),

$$\begin{aligned} D_{12} = \dfrac{L_1 L_2}{\pi q b^2} \sigma _{12}, \end{aligned}$$
(32)

where (see [26])

$$\begin{aligned} \begin{aligned} \sigma _{12}&:= \dfrac{1}{8 } \int \limits _{-b/2}^{b/2}\int \limits _{-b/2}^{b/2}\dfrac{h^2 - (y_2-y_1)^2}{\left[ h^2+(y_2-y_1)^2\right] ^2} \text {d}y_1 \text {d}y_2\\&(\text {posing } \eta _i:= y_i/h, i\in \{1,2\})\\&= \dfrac{1}{8 } \int \limits _{-\frac{1}{2k}}^{\frac{1}{2k}}\int \limits _{-\frac{1}{2k}}^{\frac{1}{2k}} \dfrac{1 - (\eta _2-\eta _1)^2}{\left[ 1+(\eta _2-\eta _1)^2\right] ^2} \text {d}\eta _1 \text {d}\eta _2\\&= \dfrac{1}{8 } \int \limits _{-\frac{1}{2k}}^{\frac{1}{2k}} \left[ \dfrac{\eta _2-\eta _1}{\left( \eta _2-\eta _1\right) ^2+1} \right] _{-\frac{1}{2k}}^{\frac{1}{2k}} \text {d}\eta _1\\&=\dfrac{1}{8} \ln \left( 1+\dfrac{1}{k^2}\right) . \end{aligned} \end{aligned}$$
(33)

As for \(\sigma _{ii}\) coefficients, they are not finite.

Indeed, the constant circulation violates the boundary condition of zero circulation at the tips. However, this assumption provides a good approximated result for the value of the mutual induction coefficient [26] and has been used in some preliminary methods, such as Betz’s (briefly detailed in [32]). In the following, this hypothesis will be useful to model the fact that the innermost wings of a optimal N-plane assume an almost constant circulation (apart from the tip regions) [14]. The approximation is as good as \(N\rightarrow \infty\).

Again Pistolesi [26]:

It can be proved that one can indefinitely approach the [induced] drag of the box-wing the number of planes. Moreover, to achieve the minimum induced drag condition, the outermost planes must carry the most of the lift, whilst the innermost ones carry only a small amount, as small as larger is their number.Footnote 3

This fact may suggest to asymptotically approximate the lifting distribution of the innermost wings as constant (flat) ones.

It can be noted that the so obtained coefficient tends to diverge as \(k\rightarrow 0\) due to the tacitly assumed modelling hypothesis that free vortices depart only from wing tips [26]. Since the unrealistic values of the coefficient for small k (Pistolesi suggests the range of validity as \(k>1/54\)), it is possible that the numerical results based on this coefficient may diverge as well as \(k\rightarrow 0\).

4.3 A Constant and an Elliptic Circulation

It is convenient to evaluate the mutual induced drag coefficient also for two lifting segments, the former having an elliptic lift distribution, the latter having a constant one. In analogy with the previous subsections, we have:

$$\begin{aligned} \begin{aligned} \sigma _{12}&:= \dfrac{1}{2\pi } \int \limits _{-b/2}^{b/2}\int \limits _{-b/2}^{b/2} \sqrt{1-\left( \dfrac{2y_1}{b}\right) ^2} \\&\quad \times \dfrac{h^2 - (y_2-y_1)^2}{\left[ h^2+(y_2-y_1)^2\right] ^2} \ \text {d}y_1 \ \text {d}y_2 \\&\quad (\text {posing } \eta _i:= y_i/h, i\in \{1,2\})\\&= \dfrac{1}{2\pi } \int \limits _{-\frac{1}{2k}}^{\frac{1}{2k}}\int \limits _{-\frac{1}{2k}}^{\frac{1}{2k}}\sqrt{1-\left( 2k\eta _1\right) ^2} \\&\quad \times \dfrac{1 - (\eta _2-\eta _1)^2}{\left[ 1+(\eta _2-\eta _1)^2\right] ^2} \text {d}\eta _1 \text {d}\eta _2 \\&= \dfrac{1}{2\pi } \int \limits _{-\frac{1}{2k}}^{\frac{1}{2k}} \sqrt{1-\left( 2k\eta _1\right) ^2}\left[ \dfrac{\eta _2-\eta _1}{\left( \eta _2-\eta _1\right) ^2+1} \right] _{-\frac{1}{2k}}^{\frac{1}{2k}} \text {d}\eta _1. \end{aligned} \end{aligned}$$
(34)

The integral in Eq. (34) has been solved through numerical integration, for some values of k. The results are listed in Table 2.

Table 2 \(\sigma _{12}\) coefficient, obtained via numerical integration of Eq. (34) for different k values, in the case of one wing loaded with elliptic distribution and one with a constant distribution
Full size table

Considering k in the range [0, 0.5], the following fitting approximation, obtained with MATLAB® Curve Fitting Toolbox, is proposed (coefficient of determination of the interpolation \(R^2=0.9999\)):

$$\begin{aligned} \sigma _{12} = \dfrac{p_2\,k^2 + p_1\,k + p_0}{k^4 + q_3\,k^3 + q_2\,k^2 + q_1\,k + q_0}, \end{aligned}$$
(35)

where

$$\begin{aligned} \begin{aligned} p_2&= 0.1998&p_1&=0.008414&p_0&=5.933\times 10^{-6}\\ q_3&= 0.8707&q_2&= 0.3067&q_1&= 0.009049 \\ q_0&= 5.933\times 10^{-6}. \end{aligned} \end{aligned}$$

Equation (35) has been assumed in the present work for the numerical part.

As for \(\sigma _{ii}\) coefficients, for the constant-loaded wing it is not finite, whilst for the elliptic-loaded wing it is equal to 1.

Figure 3 shows the mutual induced drag coefficients, for different values of k, for a pair of elliptic-loaded wings (Eq. (27)), for a pair of constant-loaded wings (Eq. (33)), and for a pair of one elliptic-loaded wing and one constant-loaded wing (Eq. (34)) (see also [26, Fig. 138]). As stated before, for \(k\rightarrow 0\), the mutual induction coefficient given by Eq. (33) diverges.

Fig. 3
figure 3

\(\sigma _{12}\) coefficient for the three cases of Sects. 4.1, 4.2, 4.3

Full size image

As a final remark of this section, the total induced drag for a biplane is given by Eq. (13). For an elliptic circulation, \(\sigma _{ii}\) (\(i=\{1,2\}\)) is exactly equal to 1, whilst in the other cases it diverges. In the following, we assume in any case \(\sigma _{ii}=1\). Conversely, \(\sigma _{12}\) will assume one of the three expressions derived in the previous subsections.

From a formal point of view, it is equivalent to consider elliptic distributions in the evaluation of the self-induced drag coefficients, regardless of the circulation distribution used to define the mutual induced drag component.

Even though this approach is not rigorous, since self-induction and mutual induction coefficients may be based on different hypotheses in the same computation, it is a simplifying hypothesis already appearing in the method of Betz [32].

5 Induced Drag of Optimal Multiplanes

The concepts exposed in the two previous sections can be extended to a generic N-wing multiplane (or N-plane), as depicted in Fig. 4, given once and for all b and the total lift L. Prandtl states in [32], at the very end of his report, that the minimum induced drag lifting system, i.e. the BWS, is the multiplane having, at the limit, an infinite number of wings. The result is asymptotically equivalent to a box-wing lifting system. Moreover, Prandtl provides the curve of the performances of the BWS in terms of BWS-to-monoplane induced drag ratio, approximated in the range \(k\in [0, 0.5]\) (or \(k\in [1/12, 0.5]\) as reported in [30]) as

$$\begin{aligned} \left( \dfrac{D}{D_\mathrm{mono}}\right) _\mathrm{opt} = \dfrac{1+0.45 k}{1.04+2.81k}. \end{aligned}$$
(36)

We have already conceptually summarised in Sect. 3 how the BWS curve, approximated by Eq. (36), has been obtained.

To answer the research question 2, our approach is the following: to generalise the induced drag expression of Eq. (13) for the case of N-planes and to pass numerically to the limit for \(N\rightarrow \infty\).

Fig. 4
figure 4

N-wing multiplane

Full size image

According to Sect. 2, the total induced drag of a N-wing multiplane can be decomposed as the sum of N self-induced drag contributes and \(N(N-1)\) drag contributes of mutual inductions. Assuming \(\sigma _{ii}=1\), and since \(\sigma _{ij} = \sigma _{ji}\), we have:

$$\begin{aligned} \begin{aligned} D&= \sum \limits _{i=1}^N D_i + 2\sum \limits _{i=1}^{N-1} \sum \limits _{\begin{array}{c} j=i+1 \end{array}}^N D_{ij}\\&= \dfrac{1}{\pi q b^2}\left( \sum \limits _{i=1}^N L_i^2 + 2 \sum \limits _{i=1}^{N-1} \sum \limits _{\begin{array}{c} j=i+1 \end{array}}^N L_i L_j \sigma _{ij}\right) . \end{aligned} \end{aligned}$$
(37)

Since the distance between two generic wings of the multiplane is a multiple of \(h/(N-1)\), as shown in Fig. 4, we can generalise the mutual induction coefficients simply by substituting k with \(k\dfrac{\left| j-i\right| }{N-1}\) in their aforementioned expressions. Hence, we introduce the generalised mutual induction coefficient between the ith and jth wing (\(i\ne j\)) by posing

$$\begin{aligned} \sigma _{(j-i)}(k) := \sigma _{ij}\left( k\dfrac{\left| j-i\right| }{N-1}\right) . \end{aligned}$$
(38)

It is easy to see that for the generalised coefficients the following properties hold:

$$\begin{aligned} \begin{aligned} \sigma _{(j-i)}&= \sigma _{(i-j)} = \sigma _{([j+\alpha ]-[i+\alpha ])}\\ \sigma _{(1-i)}&= \sigma _{(N-(N-i+1))} \end{aligned} \end{aligned}$$
(39)

for every \(\alpha \in {\mathbb {Z}}\) such that \(1\le i+\alpha ,\, j+\alpha \le N\).

The idea is now to find the stationary point of Eq. (37), i.e. the optimal induced drag value and the resulting minimiser in terms of wings lift. It can be proved [14] that wings in symmetrical position with respect to the \(z=h/2\) axis (see Fig. 4) have the same circulation under (true) optimal condition. In our approximate model, we can assume this fact valid as well. Therefore, in the following, the external wings, namely the 1st and the Nth, have the same lifting capacity \(L_1 = L_N=:l\). Furthermore, we assume that the internal wings, from the 2nd to the \(N-1\)th, have the same lifting capacity. Since the innermost wings are \(N-2\), the common value of their lift is equal to \(\dfrac{L-2l}{N-2}\). This further simplifying hypothesis is justified at the light of the numerical results presented in [14]. For instance, for a multiplane having 20 wings, assuming \(k=0.2\), the innermost wings contribute each one from the 2.1% to the 3.3% to the total lift, whilst the outermost ones contribute for the \(28.3\%\) each one.

Introducing the dimensionless parameter

$$\begin{aligned} \lambda := \dfrac{l}{L}, \end{aligned}$$
(40)

and the minimal induced drag for a monoplane for assigned span and lifting capacity

$$\begin{aligned} D_\mathrm{mono} = \dfrac{L^2}{\pi q b^2}, \end{aligned}$$
(41)

we have

$$\begin{aligned} \dfrac{D}{D_\mathrm{mono}}&= \! 2\lambda ^2\left( 1+\sigma _{(N-1)}\right) + \dfrac{(1-2\lambda )^2}{N-2} \nonumber \\&\quad + 4\lambda \dfrac{1-2\lambda }{N-2}\sum \limits _{i=2}^{N-1}\sigma _{(i-1)} \nonumber \\&\quad + 2\dfrac{(1-2\lambda )^2}{(N-2)^2}\sum \limits _{i=2}^{N-1} \sum \limits _{\begin{array}{c} j=i+1 \end{array}}^{N-1}\sigma _{(j-i)} \end{aligned}$$
(42)

or, more concisely,

$$\begin{aligned} \dfrac{D}{D_\mathrm{mono}} = \lambda ^2a_1 + \lambda (1-2\lambda )a_2 + (1-2\lambda )^2 a_3, \end{aligned}$$
(43)

where

$$\begin{aligned} a_1&:= 2\left( 1+\sigma _{(N-1)}\right) \nonumber \\ a_2&:= \dfrac{4}{N-2}\sum \limits _{i=2}^{N-1}\sigma _{(i-1)}\nonumber \\ a_3&:= \dfrac{1}{N-2} + \dfrac{2}{(N-2)^2}\sum \limits _{i=2}^{N-1} \sum \limits _{\begin{array}{c} j=i+1 \end{array}}^{N-1}\sigma _{(j-i)}. \end{aligned}$$
(44)

Imposing the first optimality condition \(\dfrac{\mathrm {d} (D/D_\mathrm{mono})}{\mathrm {d} \lambda } = 0\), we obtain

$$\begin{aligned} \lambda _\mathrm{opt} = \dfrac{-a_2+4\,a_3}{2\,a_1-4\,a_2+8\,a_3}. \end{aligned}$$
(45)

Clearly, values for \(\lambda _\mathrm{opt}\) physically acceptable are in the range [0, 0.5]. In fact, the generated lift of any wing cannot be negative under optimal conditions [14], as already discussed in the Introduction.

Even though this property holds for true optimal solutions, we assume the validity also for our approximate model.

Moreover, \(\lambda _\mathrm{opt}\) can be, at most, equal to 0.5: it is the case when the total lift is equally shared among the outermost wings, and the innermost ones are necessarily unloaded.

Finally, we pose the optimal N-plane performance as

$$\begin{aligned} \left( \dfrac{D}{D_\mathrm{mono}}\right) _\mathrm{opt} := \left( \dfrac{D}{D_\mathrm{mono}}\right) (\lambda _\mathrm{opt}). \end{aligned}$$
(46)

We stress the fact that Eq. (46) represents the performances of a N-plane, optimal in the sense of our statement of the problem.

6 Results

Equation (46) has been evaluated for different number of wings and choosing different expressions of the mutual influence coefficient given by Eqs. (27), (33), (34) and by the arithmetical average of the first two ones (hereafter referenced as hybrid coefficient).

Four cases are hereafter considered (in the following, \(i,j\in \{1,2\}\), \(i\ne j\)):

  1. Case EE

    All wings present an elliptic circulation distribution (\(\sigma _{ii}=1\), \(\sigma _{ij}\) given by Eq. (27) and generalised through Eq. (38)).

    It is the only case of the four which satisfies the zero-circulation condition at wingtips. This case corresponds to the direct method of Prandtl for the assessment of the optimal biplane and triplane performances.

  2. Case CC

    All wings present a constant circulation function. (\(\sigma _{ii}=1\)Footnote 4 and \(\sigma _{ij}\) given by Eq. (33) and generalised through Eq. (38)). Clearly, this case violates the boundary conditions of zero circulation at wingtips. However, as noticed in Sect. 4.2, this case is of interest. As the number of wings increases, the lifting distributions assume an almost flat profile (apart, obviously, near the tips) that can be preliminary approximated as an everywhere constant one. Despite the violation of tip conditions, this approximation has been used in some approximated methods.

  3. Case EC

    The outermost wings present an elliptic circulation distribution, whilst the innermost ones a constant distribution (\(\sigma _{ii}=1^4\) and \(\sigma _{ij}\) accordingly given by Eq. (27), (33), (34), generalised through Eq. (38)). Similarly to the CC case, the assumption of constant circulations violates the boundary conditions at wingtips. However, this case is an asymptotic approximation of the actual optimal lift distribution (almost flat in the innermost wings and almost elliptical in the outermost ones, at least for small k).

  4. Case HY

    All wings present a superposition of an elliptic and a constant circulation distribution. (\(\sigma _{ii}=1^4\) and \(\sigma _{ij}\) given for hypothesis by the arithmetic average of Eqs. (27) and (33)—hybrid coefficient—generalised through Eq. (38)). Similarly, the constant part of circulation violates the boundary conditions at tips. However, this case will be used to recover the approximated solution of Frediani and Montanari (constant-plus-elliptical circulation).

We remark the fact that the optimality of Eq. (46), for each of the four cases, is to be intended in the sense of the optimal solution in the particular class of circulations shaped accordingly to the case.

Figures 5, 6, 7, 8 show the results of the limit infty-plane (a finite value of 1000 wings has been assumed in the computation) for the four aforementioned cases EE, CC, EC, HY, respectively. The curves so obtained (in black with square markers) are compared with the BWS prevision of Prandtl (Eq. (36), light-blue dashed curve), with the exact solution numerically obtained by Demasi et al. [15] (red dashed curve) represented by

$$\begin{aligned} \left( \dfrac{D}{D_\mathrm{mono}}\right) _\mathrm{opt} = \dfrac{1}{1+1.7433\,k^{0.8230}}, \end{aligned}$$
(47)

and with the approximate solution of Frediani and Montanari [16] (magenta dashed curve), which assumes the upper and lower wings loaded with a constant added to an elliptical distributions.

Fig. 5
figure 5

Infty-plane performances. Case EE

Full size image
Fig. 6
figure 6

Infty-plane performances. Case CC

Full size image
Fig. 7
figure 7

Infty-plane performances. Case EC

Full size image
Fig. 8
figure 8

Infty-plane performances. Case HY

Full size image

It is observed that the EE-curve (Fig. 5) presents an important deviation from the BWS solutions of Prandtl and Demasi et al. It is surprising, since EE is the only case which does not violate the zero circulation at wingtips in the evaluation of the induced drag, as stated before. Case CC (Fig. 6) is very close to the BWS curves, apart from the range of small k, where the solution clearly diverges, as expected. Case HY (Fig. 8) is very close to the solution of Frediani and Montanari (apart from the range of small values of k), with a surprising fitness, all of the simplifying assumptions being considered. The divergent trend for small k is due to the same reasons exposed for case CC. Case EC (Fig. 7) is very adherent to the solution of Prandtl and Demasi et al. In particular, for \(k<0.1\) it is more adherent to Demasi et al.’ solution, whilst for the remaining values it is more adherent to Prandtl’s. It is observed that, despite the presence of constant circulations within the EC case, the limit curve does not diverge for \(k\rightarrow 0\). Since the (convergent) EC case differs from the (divergent) CC case for the circulation on the outermost wings, the convergent behaviour of the EC case must be due to this reason. In particular, in the EC case, some mutual induction coefficients converge to 1 as \(k\rightarrow 0\). Numerically, it is apparently sufficient to avoid the divergence of the EC limit curve for small k.

The fact that CC and EC curves are close can be explained by looking at Fig. 3, where the mutual induction coefficients for the CC and EC cases are very close for \(k>0.1\).

For the sake of completeness, Figs. 9, 10, 11, 12 show the “convergence” of the limit solution as the number of wings increases.

Fig. 9
figure 9

Convergence as the number of wings increases. Case EE

Full size image
Fig. 10
figure 10

Convergence as the number of wings increases. Case CC

Full size image
Fig. 11
figure 11

Convergence as the number of wings increases. Case EC

Full size image
Fig. 12
figure 12

Convergence as the number of wings increases. Case HY

Full size image

For graphical purposes, the dimensionless induced drag curves are plotted as a function of 1/N, for different values of k. Whenever for \(k\rightarrow 0\) the curves diverge, they have not been plotted.

The curves generated by the present model (continuous curves) are compared with the BWS curves (same formatting of the previous figures) and with the biplane performances, from [26], given by the approximation

$$\begin{aligned} \left( \dfrac{D}{D_\mathrm{mono}}\right) _\mathrm{opt} = \dfrac{1+1.63\,k}{1.027+3.84\,k} \end{aligned}$$
(48)

for the optimal case, and by the approximation

$$\begin{aligned} \left( \dfrac{D}{D_\mathrm{mono}}\right) _\mathrm{opt} = \dfrac{1.0275+1.57\,k}{1.055+3.7\,k} \end{aligned}$$
(49)

for the case of elliptic circulations on both wings. Apart the limit cases for \(N\rightarrow \infty\), which have been already presented, Fig. 9 shows an excellent agreement, for the biplane, with the results of Pistolesi. Moreover, the same figure shows a rapid convergence to the limit curve: increasing the number of wings from 10 to 1000 does not produce significant increments. Figure 10 shows the convergence for the CC case. It is observed an oscillation of the curves (as N increases) around the BWS one. This fact is probably due to the non-physical condition of constant circulations. In any case, the case CC is meaningful only for large N, and the actual physical interpretation of the solution for a small number of wings should be disregarded. The same oscillations are observed for the cases HY and EC in Figs. 12, 11, due to the presence of constant distributions in the model. Also in these cases, the limit solution is the only one which should be interpreted.

Figures 13, 14, 15, 16 show the values of \(\lambda _\mathrm{opt}\) for the four cases, as the number of wings increases.

Fig. 13
figure 13

\(\lambda _\mathrm{opt}\) values as the number of wings increases. Case EE

Full size image
Fig. 14
figure 14

\(\lambda _\mathrm{opt}\) values as the number of wings increases. Case CC

Full size image
Fig. 15
figure 15

\(\lambda _\mathrm{opt}\) values as the number of wings increases. Case EC

Full size image
Fig. 16
figure 16

\(\lambda _\mathrm{opt}\) values as the number of wings increases. Case HY

Full size image

For the EE case (Fig. 13), the optimal value of the lift of an outer wing is equal to 0.5 for \(k\rightarrow 0\), for all numbers of the wings. As it increases, \(\lambda _\mathrm{opt}\) decreases. A saturation effect is presented for \(N>100\). More complicated trends are depicted in Figs. 14 and 16 for the CC and HY cases. For small values of k, the optimal solution lies on the boundary of the feasible range of \(\lambda _\mathrm{opt}\). In fact, it assumes the value of 0 or 0.5 for \(k<0.1\). It is due to the non-strictly convexity of the induced drag determined from our model for these two cases and for small k. However, as already stated, only the asymptotic solution is meaningful. Figure 15 shows the counterpart plot for case EC. Remarkably, the triplane violates the condition \(\lambda _\mathrm{opt}=0.5\) for \(k=0\). This fact can be due to pure numeric issues. Consider the cases EC and EE for \(N=3\) and \(k=0.0001\) (they differ only for the central wing circulation shape). For the EE case, we obtain \(a_1=3.9985\), \(a_2=3.9985\), \(a_3=1\), whilst for the EC case we obtain a slightly different \(a_2=3.9799\). Nevertheless, substituting these numbers into Eq. (45), we obtain \(\lambda _\mathrm{opt}=0.5\) and \(\lambda _\mathrm{opt}=0.2594\), respectively. However, also in the EC case, only the asymptotic solution is meaningful.

It is noteworthy to investigate also the asymptotic trend for \(k\rightarrow \infty\), as shown in Fig. 17 for the EC case only.

Fig. 17
figure 17

Asymptotic behaviour as \(k\rightarrow \infty\) for the EC case

Full size image

It is observed that the infty-plane for the EC case is able to grasp the actual behaviour of the BWS, described by the solution of Demasi at al. (Eq. (47)). As already observed in [15], the Prandtl approximation of Eq. (36) predicts a limit value of 0.1601 for \(k\rightarrow \infty\), whilst physically the limit should be zero.

Figure 18 shows the biplane and triplane curves (EE case) and the limit infty-plane curve (EC case) plotted over the original Fig. 7 of Prandtl report [32].

Fig. 18
figure 18

Plot of biplane (EE), triplane (EE) and infty-plane (EC) curves over Fig. 7 of [32]

Full size image

7 Discussion and Conclusions

It this paper we tried to throw light on the method Prandtl used for the determination of the induced drag of the Best Wing System, i.e. of the lifting system achieving the minimum induced drag among all the other ones for given lift and span. Noticing the existence in the same work of two methods, a direct method for the evaluation of the induced drag of optimal biplane and triplane, and a mysterious method for the assessment of the Best Wing System one, we wondered four research questions that we will be answered in this section. Besides, we proposed a simple method for the evaluation of the induced drag of multiplanes, using different expressions (and assumptions) for the mutual induction coefficient, extending and generalising the direct method Prandtl used for the assessment of the performance of biplanes and triplanes. In so doing, we additionally wondered if Prandtl could have achieved the Best Wing System performance with our proposed method which, we remark, is a generalisation of a method he initially used.

We first try to answer the first and second research questions, i.e. Why did Prandtl change method when passing from biplane/triplane to the Best Wing System? and What would have happened if he continued with the method used for biplane and triplanes, the number of wing increasing?

Prandtl evaluated the optimal biplane and triplane performances adopting elliptic circulation distributions on all the wings (EE case), claiming that it would be a more refined and conceptually rigorous approach than the method of Betz, which corresponds to the evaluation of the mutual induced drag coefficient as the circulation distributions were constant (CC case). Our implementation of the former case perfectly matches Prandtl results, which have been also validated by numerical simulations, obtained through Vortex Lattice Methods and Boundary Element Methods, as reported in [2, 33].

Moreover, Prandtl’s solution for the Best Wing System was extremely close to the exact solution of the minimization of induced drag for a given lift and wingspan. A remarkable result, especially if it is observed that despite the promising prediction power of assuming all elliptic circulations on the wings, Prandtl changed the method to obtain the Best Wing System curve.

Figure 9 shows that assuming all elliptic circulations on all the wings leads to grasp the optimal performances only when N is small (biplane, triplane), whilst it fails in predicting the Best Wing System performance. This answer to the second research question.

The question why he abruptly moved from considering all elliptical distributions, through which he assessed the biplane and triplane performances, to the potential approach of Grammel and Pohlhausen for the Best Wing System, asymptotically equivalent to an infty-plane, remains open.

We now try to answer to the third research question: Which is the Best Wing System lift distribution? Why did not he provide it? As a matter of fact, Prandtl presented just the performance curve of the Best Wing System, not its optimal circulation. It is reasonably unlikely that Prandtl knew the circulation distribution of the Best Wing System but did not publish it in his works. As a matter of fact, in other Prandtl’s works as [28,29,30,31], in Pistolesi Aerodynamics handwritten notes for the academic year 1923–1924 [25] and book 1932 [26], no attention is paid to the shape of the lift distribution, which indeed is a fundamental part of the solution. The qualitative representation of the optimal aerodynamic load appears in a later publication by Von Karman and Burgers in 1935 ([34], Fig. 81, reproduced in our Fig. 19), showing that the load is not made by an elliptical and constant distributions over the horizontal wings, in agreement with the findings in [13].

Fig. 19
figure 19

Reproduction of Fig. 81 from [34] showing the optimal lift distribution of the BWS

Full size image

It is noteworthy that Fig. 2d shows the fictitious surface \(\Sigma _{BWS}\), directly linked to the jump of the potential function when crossing lifting elements, from which one could derive the shape of the optimal lift distribution, at least qualitatively. It is not clear in Fig. 2d by Prandtl (also present in Pistolesi’s book) whether the non-constant part of the lift is elliptic or not. For instance, Fig. 20 shows the same picture of Fig. 2d compared with green ellipses. The matching is pretty adherent.

Fig. 20
figure 20

The fictitious surface \(\Sigma\) for the optimal box-wing (from [30]) compared with elliptical distributions

Full size image

As already noticed, apart from a scaling factor, the qualitative shape of the optimal lift distribution can be obtained by summing the value of the potential on the dorsal side minus the value of the potential on the ventral side, thus obtaining qualitatively Fig. 21 if we start from Fig. 2d.

Fig. 21
figure 21

Qualitative lift distribution from the \(\Sigma\) of Fig. 2d

Full size image

We remark that the explicit sketch of the optimal lift distribution is provided only in [34], which agrees with the modern computations.

The poor attention paid to the optimal lift distribution for the box-wing can be explained by the scepticism by which Prandtl and coeval scientists regarded the box-wing architecture from the practical engineering viewpoint. In fact [26]

These considerations have only a pure theoretical interest, since, besides the induced drag, one must consider also the friction drag. Moreover, it is not possible to imagine lifting systems without structural elements such as rods, wires, etc.

Indeed, these considerations have been done in the 1920s.

Before answering the fourth research question, another point noteworthy to discuss is the actual interpretation of the asymptotic 0.1601 value of Eq. (36), in contrast with the actual limit: zero. How to interpret this fact?

We claim that 0.1601 is the result of the particular form of the approximating function chosen by Prandtl for Eq. (36), i.e. a function of the form

$$\begin{aligned} \dfrac{1+\alpha \,k}{\beta + \gamma \,k}, \end{aligned}$$
(50)

where \(\alpha\), \(\beta\) and \(\gamma\) are constants to be determined by the curve fitting. We preliminary note that the limit for \(k\rightarrow \infty\) equals \(\alpha /\gamma\) (if \(\gamma \ne 0\)), and that it is zero if and only if \(\alpha =0\). We have interpolated the values provided by Prandtl in Tab. 3 of his report for \(k\in [0.05,0.5]\) (note that for \(k=0\) Prandtl Eq. (36) fails) with the MATLAB® Curve Fitting Toolbox using a function of the form (50). The most accurate solution (\(R^2=0.9999\)) is \(\alpha =0.4913, \beta =1.041, \gamma =2.891\), which provides an asymptotic value of 0.1699. Other accurate solutions (\(R^2=0.9999\)) are \(\alpha =0.4881, \beta =1.04, \gamma =2.89\), which provides an asymptotic value of 0.1688, and \(\alpha =0.45, \beta =1.04, \gamma =2.826\), which provides an asymptotic value of 0.1592. The three proposed curve fittings are very similar to Prandtl Eq. (36). Of course, Prandtl did not have modern calculators to provide the most accurate approximation. Therefore, it can be stated that Eq. (36) has been obtained as an approximation of the Best Wing System performance, valid only in the range \(k\le 0.5\) because it has been obtained interpolating known values of the dimensionless induced drag for k up to 0.5. Asymptotic extrapolations of Eq. (36) are not legit.

Moreover, it is likely that Prandtl knew also that the asymptotic value is zero because of a simple argument available at the time. As stated before, the ratio between the induced drag of the optimal Best Wing System and the induced drag of the optimal monoplane (same lift and span) can be expressed as \(\Sigma _1/\Sigma _{BWS}\). We have already introduced \(\Sigma _1 = \frac{\pi }{4}b^2\). Similarly, we may pose \(\Sigma _{BWS} = kb^2 + 2S\), where the first term is the portion of \(\Sigma _{BWS}\) inside the optimal box-wing of dimension \(b\times h\) and the second term is the portion of \(\Sigma _{BWS}\) outside the box-wing (see Fig. 2d). Therefore, Eq. (20) simplifies to:

$$\begin{aligned} \left( \dfrac{D}{D_\mathrm{mono}}\right) _\mathrm{opt} = \dfrac{\pi b^2}{4 \left( k\,b^2 + 2S \right) }, \end{aligned}$$
(51)

which tends to zero as \(k\rightarrow \infty\) with a decay rate of \(\frac{\pi }{4\,k}\).

Finally, we answer the fourth research question: could Prandtl extend the “direct method” to formalise a model similar to the one we have proposed?

At Prandtl’s time, the assumptions we have adopted to obtain our results were likely to be known.

The fact that the Best Wing System is asymptotically equivalent to a infty-plane was known.Footnote 5 As stated before, Pistolesi provides some evidences [25, 26].

The fact that wings symmetrically placed with respect to the axis \(z=h/2\) (Fig. 4) have the same distribution under optimal conditions could have be assumed as a simplifying hypothesis. Indeed, for the triplane, Prandtl stated

it may be assumed at once that the upper and lower wings have the same lift for the minimum drag.

This assumption could have been extended for the generic multiplane, adding the simplifying hypothesis of equal distributions on the innermost wings. It allows expressing the total induced drag as a quadratic function of the dimensionless lift of the outermost wings only (\(\lambda\)), as done in Eq. (42), and to formally evaluate the minimiser Eq. (45).

The fact that the innermost wings generate less and less lift as the number of wings increases, as already stated in Sect. 4.2, was likely to be known at Prandtl’s time. However, it could have been assumed by Prandtl looking at the triplane solution. Therefore, the approximation of almost flat distributions with constant ones naturally derives.

Moreover, it has been discussed that, at the time, also the fact that under optimal conditions every wing of a multiplane must contribute with positive lift.

All this considered, an asymptotic model such as the proposed one may have been conceived by a 1920s aerodynamicist.

Our results and considerations have been made to provide insights on Prandtl’s Best Wing System, a result of about a century ago which still influences the current aeronautic scenario.