2.1 Cantilever wing geometric description

The 3D wing geometric model is illustrated in Fig. 2. It is versatile and can be used for a cantilever wing, the upper and lower wings of a joint wing, and the wing and strut of a strut-braced wing configuration. It can also be used for any aircraft’s fin, tail and canard. The count and pitches of the wing bays and sections are variable. In addition to varying the span length ${b_{Wing}}$ , and total number of wing sections ( ${N_{Wing}}$ ), each wing section can be manipulated by changing their chord lengths ( ${c_j}$ ), thicknesses ( ${t_j}$ ), sweep angles ( ${\theta _j}$ ) and dihedrals ( ${\beta _j}$ ). Unless these design values are assigned at all the sections, they will be calculated with linear interpolation between the closest known two values or default values. To generate this geometric model in 3D axes, coordinate points of each section in $x,y$ and $z$ axes should be computed using well known geometric relations.

2.2 Wing box geometric description

Structural idealisation reduces the computational time of the structural analysis of wing boxes in the conceptual or preliminary design stage. Following the procedure in Ref. (Reference Megson43), an idealised single-cell wing box model was developed. The idealisation process is illustrated in Fig. 3. The wing box shapes are represented with simplified geometries while the stringers and the spar caps are idealised with mass booms. The following assumptions are made in the development of the new wing box model.

• • Total stringer counts, ${N_{Boom}}$ , can be changed at each wing box section, but the numbers of stringers located in the upper and lower surfaces of the wing box sections are equal.

• • The wing box is assumed symmetric around two reference axes shown in Fig. 4. This assumption provides a direct solution rather than an iterative one in calculating mass boom cross-sectional area as explained in Section 4.1 and reduces the computational time. The suitability of this assumption for the proposed ME method is presented with validation studies.

Figure 3. Illustration of the wing box idealisation and the comparison of the actual and idealised wing box sections.

Figure 4. The geometric model representation of the idealised wing box section.

The position of the front and rear spars, shown in Fig. 2, can be changed with the non-dimensional (ND) parameters of ${\varepsilon _{fs,j}},{\varepsilon _{rs,j}}$ , respectively. These local parameters are the fraction of the local wing chord length. The cross-sectional area of the spar caps can be scaled to that of the stringer area using a ratio of $R$ . Wing box edge height is a ND variable, ${\Psi _{Edge,j}}$ , as a fraction of wing thickness, ${t_{Wing}}$ . The local thickness of the idealised wing box, ${t_{wb}}$ ,can be different than that of aerofoil if ${\Psi _u}$ , shown in Fig. 4, is assigned other than zero. ${\Psi _u}$ is the ND variable as a fraction of wing thickness. The wing box shape can be transformed from a diamond to a hexagon and a rectangular shape with the change of the ${\Psi _{Edge,j}}$ as can be seen in Fig. 4. The wing box geometry is generated with these parameters and the wing planform parameters explained above.

Stringers are distributed equally, crossing their reference line with the idealised wing box skin. The stringers are assumed to be positioned parallel to the wing box reference line with a constant pitch of ${d_{stringer}}$ . The stringers are cut off when they contact with the spars; hence, the stringer counts are reduced towards the wing tip, as shown in Fig. 5. The front and rear spars are consistent along the wing from the root to the tip. The effects of gaps and overlaps on the wing mass due to the sweep angle change are not considered separately, illustrated in Fig. 5, and are assumed to balance each other to simplify the analysis.

2.3 Other components and the visualization

A simple fuselage geometric model is developed. The width and length of the fuselage can be defined by ${D_F}$ and ${l_F}$ , respectively. The nose shape and the tail cone are generated by the default value of the ND parameters as a fraction of fuselage diameter and length. For simplicity, the fuselage cross-section is assumed to be a circle. Landing gears, engines and fuel drop tanks are considered concentrated loads, and their geometry is assumed to be cylinders. For a wing-mounted engine configuration, the concentrated loads are assumed to be attached to wing sections shown in Fig. 2 for simplicity in computations.

Each component can be generated and located in a global coordinate system using the geometric models explained above to create the entire aircraft. The origin of the global coordinate system is chosen as the nose of the aircraft. Therefore, 2D and 3D visualisations of the whole aircraft can be generated with this approach and sample views are presented in the following sections.

3.1 Load distribution models

Although the weight of aircraft components can be estimated at the initial design stage, the distribution shape of these loads may not be known and should be estimated. Triangular, rectangular, trapezoidal and elliptical distribution models can provide acceptable accuracy at the initial design stage and are used in this study for initial inertia load distribution. The total magnitude of the distributed load acting along a bay of the wing shown in Fig. 2 can be assumed as a strip load and positioned at the wing local section for simplicity in shear force and bending moment calculations.

Table 3. The parameters in the table are combined to generate rolling and combined cases

3.2 Aircraft mass moment of inertia model

A computationally fast MMI model with an error margin of ± 10% is developed considering its minimal effect on wing mass estimations. The aircraft components are split into simplified geometric shapes of cylinders and cuboids as shown in Fig. 6. Their MMIs can be calculated at their local centre of gravity (CG) and then carried to the reference lines with the parallel axes theorem. The engine, landing gears and drop tanks are modelled with cylindrical shapes. The fuselage is also represented with 50 cylindrical shapes, while the cuboids represent the wing. The MMI of the entire aircraft about $x,y$ and, $z$ axes are obtained by summing the MMI of all components.

Figure 6. Simplified geometric model for mass moment of inertia calculation of complete aircraft.

3.3 Spanwise aerodynamic loads

At the conceptual or preliminary structural design stage, semi-empirical or low-order aerodynamic methods (such as vortex lattice or panel methods) can be used to compute the aerodynamic loads. The developed ME model is compatible with both due to its discrete architecture of the lifting surface geometric models. Previously Elham et al. [Reference Elham, La Rocca and Van Tooren22] showed the minimal difference between the effect of the elliptical distribution and the distributions from low-order aerodynamic methods on wing ME. Hence, Schrenk’s lift distribution model [Reference Howe44] is used in this study considering its simplicity in application and better representation of the spanwise lift distribution compared to the elliptical distribution.

It is assumed that the wing generates a total lift ${L_{Wing,Case}}$ , which is equal in magnitude but opposite in direction to the total inertia force of the aircraft ${F_{Aircraft,\;Case}} = {W_{Aircraft,\;Case}}{n_{Loading}}$ , where ${W_{Aircraft,\;Case}}$ is the total weight of the aircraft in a particular flight case and ${n_{Loading}}$ is the loading factor. Equation (1) ensures that the lift generated by the wing balances the inertia forces acting on the aircraft during a specific flight case while also incorporating the lift generated by the tail. In Equation (1), ${E_{Tail}}$ denotes the percentage ( $\% $ ) of the total lift generated by the tail. The tail generally generates a small amount of lift, usually negative (downforce), to maintain pitch stability. At this stage, the tail’s contribution assumed as a fixed $5\% $ for simplification. This allows for easier calculations without needing detailed aerodynamic modeling, which is typically reserved for more advanced design stages. Additionally, a sensitivity analysis, as detailed in Section 5.2, confirms the minimal impact of ${E_{Tail}}$ on the total wing mass and verify the simplified approach.

(1) \begin{align}{L_{Wing,{\rm{\;}}Case}} = {\rm{\;}} - {W_{Aircraft,{\rm{\;}}Case}}{n_{Loading}}\left( {1 + \frac{{{E_{Tail}}}}{{100}}} \right)\end{align}

If the original wing planform is converted into an equivalent elliptic planform, the chord variation, ${c_{Ellip,j}}$ , along the wing can be calculated with Equation (2) and the Schrenk’s chord distribution, ${c_{Schrenk,j}}$ , is the average of the original and the elliptic chords as given in Equation (3). Once ${c_{Schrenk,j}}$ is calculated then the total lift can be redistributed along the lifting surface proportional to the value of ${c_{Schrenk,j}}$ as given in Equation (4).

(2) \begin{align}{c_{Ellip,j}} = {\rm{\;}}\frac{{4{S_{LS}}}}{{\pi b}}\sqrt {1 - {{\left( {\frac{{2{y_{LS,j}}}}{{{b_{Wing}}}}{\rm{\;}}} \right)}^2}} \end{align}

(3) \begin{align}{c_{Schrenk,j}} = \left( {{c_{Ellip,j}} + {c_j}} \right)/2\end{align}

(4) \begin{align}{P_{Schrenk,j}} = \left( {{L_{LS,{\rm{\;}}Case}}{\rm{\;}}/\left( {{S_{LS}}} \right)} \right){c_{Schrenk,j}}\end{align}

in which $j$ denotes the wing bay number, as shown in Fig. 2, ${c_j}$ is the wing section chord length, ${S_{LS}}$ is the surface area of the lifting surface, ${L_{LS,\;Case}}$ is the total lift load in a flight case, and ${P_{Schrenk,j}}$ is the lift load per unit length at wing $j$ th section.

3.3.1 Aileron loads

It was assumed that the additional lift generated by the low- and high-speed ailerons could be represented with one equivalent aileron. The required aileron lift, $\overline{Z}_{\xi }$ , to generate enough rolling moment, $\mathop {{L_\xi }}\limits^\_ $ , can be found by the general equation of $\overline{Z}_{\xi } = \overline{L_{\xi}}/2{y_\xi }$ where ${y_\xi }$ is the moment arm length of the ailerons to the A/C centreline. If it is assumed that the lift generated by the aileron is distributed in a rectangular shape, the lift per unit length of this distributed lift, ${P_{Ail,j}}$ , can be calculated by Equation (5) where ${\rho _{Air}}$ is air density, ${V_{TAS,Case}}$ is true airspeed of the flight case, ${S_{Aileron}}$ is aileron surface area, $\xi $ is aileron deflection, ${L_\xi }$ is ND rolling moment coefficient, ${\Lambda _{Ail,in}}$ and ${\Lambda _{Ail,out}}$ are ND locations of the aileron’s inboard and outboard sections, respectively.

(5) \begin{align}{P_{Ail,j}} = \frac{{\left( {{\rho _{Air}}{V_{TAS,Case}}^2{S_{Aileron}}\xi {L_\xi }} \right)b}}{{2\left( {{\Lambda _{Ail,out}}^2 - {\Lambda _{Ail,in}}^2} \right)}}\end{align}

3.3.2 Aerodynamic damping

Due to the rolling, each wing section has a local velocity ( ${p_{Case}}{y_{Section}})$ which is perpendicular to the flight path, and the local effective angles of attack of those sections are changed by $ta{n^{ - 1}}\left( {{p_{Case}}{y_{Wing,j}}/{V_{Case}}} \right)$ [Reference Howe44]. Hence, the local lift coefficient, ${C_{L,p,j}}$ , can be written as,

(6) \begin{align}{C_{L,p,j}} = \frac{{\left( {{p_{Case}}{y_{Wing,j}}} \right)}}{{\left( {{V_{Case}}} \right)}}{a_{1,j}}\end{align}

in which, ${V_{Case}}$ , and ${a_{1,j}}$ are the true airspeed of the flight case and the lift curve slope of the wing bay, respectively. If the equation is updated by Schrenks’s hypothesis [Reference Howe44] and multiplied by $0.5{\rho _{Air}}{V_{Case}}^2$ , the spanwise aerodynamic damping load per unit length can be written as in Equation (7).

(7) \begin{align}{P_{Damping,j}} = 0.5{\rho _{Air}}{V_{Case}}{p_{Case}}{y_{Wing,j}}{a_1}\left( {\frac{{{c_j}}}{2} + \frac{{2{S_{Wing}}}}{{\pi b}}\sqrt {1 - {{\left( {2{y_{Wing,j}}/b} \right)}^2}} } \right)\end{align}

in which ${p_{Case}}$ is roll rate, ${\rho _{Air}}$ is the density of air. Finally, the total spanwise aerodynamic strip load of each wing section, ${F_{Aero,j}}$ , can be computed from the total of the lift load per unit length values using Equations (4), (5), and (7) along the wing.

3.4 Inertia loads

Inertia loads on lifting surfaces are categorised into concentrated and distributed loads. Concentrated loads primarily arise from the mass of engines, landing gears and fuel drop tanks, and the distributed inertia loads are generated by the mass of the wing itself and the fuel stored within it. The inertia load at $j$ th wing section, ${W_{Total\;\;Inertia\;,j}}$ , is the sum of the total concentrated, ${W_{Tot\;Con,j}}$ , and distributed, ${W_{Tot\;Dis,j}}$ , loads at this wing section, respectively.

The maximum volume of the wing stored fuel is calculated from the maximum available internal space of each wing box bay. The maximum internal volume, ( ${V_{Fuel\;tank,Max,j}}$ ), of local the wing box section can be computed with a truncated irregular hexagonal pyramid model, shown in Fig. 3(c). ${V_{Fuel\;tank,Max,j}}$ is reduced due to the presence of structural components and devices located inside the wing box bay. Therefore, a fuel tank efficiency factor of $\;{\eta _{tank,j}}$ is introduced to find the final fuel tank volume using ${V_{Fuel\;tank,j}} = {V_{Fuel\;tank,Max,j}}\;{\eta _{Tank}}$ . The weight of the fuel ${W_{Wing\;fuel,Case,j}}\;$ can be computed with fuel density and the remaining fuel volume in each load case.

Table 4. The capabilities of the selected mass estimation methods and their implemented alternatives

The initial wing ME and distribution shape are required to compute the total inertia loads. The initial wing mass can be estimated as 0.1TOGW, or the existing empirical ME methods, discussed in Section 5.0 can be used. The wing mass can be distributed along the wing in a trapezoidal shape.

The sensitivity of the model accuracy to this initial wing ME is assessed using direct and convergent solvers. The direct solver uses the ME result of the first iteration step. In contrast, the convergent solver iterates until a 0.05% difference is achieved in wing mass estimations (MEs) of consecutive iteration steps. Both solvers were initialised by an initial wing ME value of $0.1 \times TOGW$ . The differences between the standard and average errors of the direct and convergent solvers in wing mass estimates of 13 aircraft (listed in Section 5.0 ) were within $ \pm 1\% $ . Hence, the direct solver approach is chosen as the default of the model, considering the computational time and convergent solver is used for verifications.

The total inertia loads at $j$ th wing section during the symmetric, $\;{F_{Man,j}}$ , and the rolling, ${F_{Roll,\;j}}$ , manoeuver cases are computed with Equations (8) and (9) where $\mathop p\limits^\_ $ is the angular rolling acceleration ( ${\rm{rad}}/{{\rm{s}}^2}$ ) of the lifting surface. And following the procedure in CS.25 [31], the total inertia loads in symmetric, ${F_{Total\;Inertia\;Sym,j}}$ , and combined load cases, ${F_{Total\;Inertia\;Comb,j}}$ , are computed with Equation (10) for each load case. As explained in CS25.349 [31] the rolling loads are combined with zero and two-third of the positive manoeuvring factor (see Table 3) used in design.

(8) \begin{align}{F_{Man,j}} = {W_{Total{\rm{\;\;}}Inertia{\rm{\;}},j}}{n_{Sym{\rm{\;}}Loading}}\end{align}

(9) \begin{align}{F_{Roll,j}} = \frac{{{W_{Total\;\;Inertia\;,j}}\left( {{y_{Section,j}} + {y_{Section,j + 1}}} \right)\mathop p\limits^\_ }}{{2g}}\end{align}

For just symmetric manoeuvre cases,

\begin{align*}{F_{Total\;Inertia\;Sym,j}} = {F_{Man,j}}\end{align*}

For combined cases,

(10) \begin{align}{F_{Total\;Inertia\;Comb,j}} = \frac{2}{3}{F_{Man,j}} + {F_{Roll,j}}\end{align}

3.5 Shear force and moment calculations

Spreadsheet approaches are used to compute shear force and moments along the wing. The total shear force at $j$ th wing section, ${V_{Total,j}}$ , can be obtained with ${V_{Total,j}} = {F_{Total,j}} + {V_{Total,j + 1}}$ , where ${F_{Total,j}}$ is ${F_{Total,j}} = {F_{Total\;Inertia,j}} + {F_{Aero,j}}$ . If the shear load variation between the wing sections is assumed to be linear, then the moment increment, $\Delta {M_{Tota,j}}$ , can be obtained with $\Delta {M_{Total,j}} = \left( {{V_{Total,j}} + {V_{Total,j + 1}}} \right){b_{Bay,j}}/2$ . And the total bending moment, ${M_{Total,j}}$ , at j $th$ wing section is formulated as ${M_{Total,j}} = \Delta {M_{Total,j}} + {M_{Total,j + 1}}$ .

Torsion can be caused by loads which are not located at the shear centre of the wing box such as engine, landing gears, drop tanks, engine trust and distributed loads, including aerodynamic and inertia. The CG of the wing stored fuel is assumed to be located at the shear centre of the wing box; hence, it does not cause torsion. The torsion increment, $\Delta {M_{y,Cm,j}}$ , caused by the pitching moment of the $j$ th wing bay, ${C_{m,Bay,j}}$ , acting at the quarter chord can be calculated with $\Delta {M_{y,Cm,j}} = \left( {MA{C_{Bay,j}}{S_{Bay,j}}{C_{m,Bay,j}}q} \right)\;$ where, $q$ $,{S_{Bay,j}},MA{C_{Bay,j}}$ are dynamic pressure, surface area of the lifting surface bay, and mean aerodynamic chord of the bay, respectively. The total torsion increment, $\Delta {M_{y,Tot,j}}$ is $\Delta {M_{y,Tot,j}} = {F_{Total\;Inertia,j}}{d_{T1,Arm,j}} + {F_{Aero,j}}{d_{T2,Arm,j}}$ , where ${d_{T1,Arm,j}}$ , ${d_{T2,Arm,j}}$ are the moment arm lengths between the shear centre of the lifting surface bay and the forces acting about that centre and caused by the total inertia, and aerodynamic loads, respectively. Hence, the total torsion at j $th$ section can be computed with ${M_{y,Tot,j}} = \Delta {M_{y,Tot,j}} + \Delta {M_{y,Cm,j}} + {M_{y,Tot,j + 1}}$ .

The value of ${C_{m,Bay,j}}\;$ can be determined using low-order aerodynamic methods, such as the vortex lattice method or panel methods. As a default of the tool ${C_{m,Bay,j}}\;$ is set to -0.1 for each wing bay, alternatively, the user can assign this value directly if a quicker ME calculation is needed. By employing discretised geometric and structural sizing approaches, the impact of ${C_{m,Bay,j}}\;$ values – calculated for each wing bay in various flight cases – on the wing mass can be analysed. In this study, wing pitching moment values are computed for the cruise condition without any aileron deflection using an external tool named AVL (Athena vortex lattice). These constant values are then applied uniformly across the wing bays for the pitching moment analysis. More detailed studies are recommended when the presented methods are integrated with low-order aerodynamic methods such as the vortex lattice method or panel methods.

The moments calculated in the local axis of j ${\rm{th}}\;$ wing-bays, ${M_j}$ , should be transformed to, ${{M^{^{\!\!\!\!\prime}}}_{j}}$ , using well known 3D rotation matrix if there is an angle change in the load path between $\;j{\rm{th}}$ and ( $j + 1)st\;$ wing-bays.

4.1 Mass boom sizing

For preliminary structural sizing, the structural idealisation approach is used from Ref. (Reference Megson43). It is assumed that the axial stresses are carried by the mass booms, representing the spar caps and stringers, and the shear stresses are carried by the panels as illustrated in Figs 3 and 7. The axial stress in $i$ th mass boom of the $j$ th wing box section (depicted in Fig. 7) can be written as, ${\sigma _{j,i}} = - \left( {{M_{x,j}}{z_{j,i}}} \right)/{I_{xx,j}} + \left( {{M_{z,j}}{x_{j,i}}} \right)/{I_{zz,j}} + {F_{y,j}}/({N_{boom,j}}{B_{j,i}})$ , where ${M_{x,j}}$ and ${M_{z,j}}$ are the moments acting around the $x$ and $z$ reference axes of the wing box section, respectively; ${F_{y,j}}$ , and ${N_{Boom}}$ are the spanwise axial loads and the total number of mass booms in $j{\rm{th}}$ wing box section; ${x_{j,i}}$ and $\;{z_{j,i}}$ are the distances of the mass booms to the symmetry axes (reference line) of the wing box; ${I_{xx,j}}$ and ${I_{zz,j}}$ are the second moment of areas of $j$ th wing box section around the $x$ and $y$ axes, respectively.

Figure 7. The geometric model used to compute the shear flows between the mass booms of the idealised wing box section.

Following the symmetrical wing box assumption explained in Section 2.2, the cross-section areas of the stringers are equal in $j$ th wing box section and can be characterised with ${B_j}$ for simplification instead of ${B_{j,i}}$ . In addition to that, the cross-section area of the spar caps $({B_{C,j}}$ ) is different than that of the stringers ( ${B_j})$ as shown in Fig. 7 and to include this difference in the model, a parameter of $R$ is introduced. If the ${B_{C,j}}$ is written in terms of ${B_j}$ , the second moment of area of $j$ th wing box section around the $x$ axis can be formulated as, ${I_{xx,j}} = 2*{B_j}*\left( {\sum\nolimits_{i = 2}^{\left( {{N_{Boom}}/2} \right) - 1} {z_i}^2 + R*{z_1}^2 + R*{z_{{N_{Boom}}/2}}^2} \right)$ and similarly ${I_{zz,j}}$ can be derived. If the equations for ${I_{xx,j}}$ and ${I_{zz,j}}$ are substituted in the equation for ${\sigma _{j,i}}$ and ${B_j}$ is isolated, Equation (11) can be derived to compute the cross-sectional area of a mass boom $({B_j}$ ) in $j$ th wing box section.

(11) \begin{align}{B_j} & = \left| \left( \frac{{{M_{x,j}}}}{{2\left( {\mathop \sum \nolimits_{j = 2}^{\frac{{{N_{Boom}}}}{2} - 1} {z_j}^2 + R{z_1}^2 + R{z_{{N_{Boom}}/2}}^2} \right)}}{z_i} + \frac{{{M_{z,j}}}}{{2\left( {\mathop \sum \nolimits_{j = 2}^{\frac{{{N_{Boom}}}}{2} - 1} {x_j}^2 + R \times {x_1}^2 + R{x_{{N_{Boom}}/2\;}}^2} \right)}}{x_i}\right.\right.\nonumber\\& \left.\left. \vphantom{\int_{\int_{\int}}^{\int^{\int}}} + \frac{{{F_{y,j}}}}{{\left( {{N_{boom}} - 4 + 4R} \right)}} \right)\frac{1}{{{\sigma _{Y,i}}}} \right|\end{align}

where ${\sigma _{Y,i}}$ is the yield stress of the mass boom material in the $j$ th wing box section. Yield stress ( ${\sigma _{Y,i}}$ ) is used for isotropic materials under limit loads. For quasi-isotropic composite materials failure stress, ${\sigma _{U,i}}$ , is used instead of yield stress and ultimate loads are applied instead of limit loads, as discussed above.

4.2 Panel sizing

4.2.1 Shear flow

The shear flow increment between two mass booms can be calculated using Equation (12) by employing an ‘alternative method’ from Ref. (Reference Megson43). In the case of non-constant shear force over a wing bay, the method provides an approximate solution. In Equation (12), $\Delta {P_i}$ is the axial load variation along a unit length of a boom and ${q_{i - 1}}$ is the shear flow on i–1th panel. The values of $\Delta {P_i}$ and ${q_{i - 1}}$ can be calculated with Equations (13) and (14), respectively. In Equation (14), ${q_{i - 1}}$ is arranged for the first skin panel. When unknown shear flow on the first skin panel, ${q_1}$ , is calculated with Equation (14), the set of equations in Equation (12) can be solved step by step to calculate other shear flows on the rest of the skin panels using Equation (12). Hence, the material thickness of each wing box panels along the wing to resist the shear flow can be calculated with Equation (15) for yield stress criteria under limit loads. For quasi-isotropic composite materials, ultimate shear stress will be used instead of shear yield stress, and ultimate loads will be applied instead of limit loads.

(12) \begin{align}{q_{j,i}} - {q_{j,i - 1}} = - \Delta {P_i}\end{align}

(13) \begin{align} - \Delta {P_{j,i}} = \;{\sigma _{j,i}}{B_{j,i}} - \;{\sigma _{j + 1,i}}{B_{j + 1,i}}\end{align}

(14) \begin{align}{q_{j,1}} = \frac{{\mathop \sum \nolimits_{i = 2}^{{N_{Boom}}} \left( {\mathop \sum \nolimits_{a = 2}^i \Delta {P_{j,a}}} \right){d_{Arm,j,i}}{l_{Panel,j,i}} - {M_{y,j}}}}{{\mathop \sum \nolimits_{i = 1}^{{N_{Boom}}} {d_{Arm,j,i}}{l_{Panel,j,i}}}}\end{align}

(15) \begin{align}{t_{Shearj,i}} \geq \frac{{{q_{j,i}}}}{{{\tau _{Y,j,i}}}}\end{align}

in which, $j$ , $i$ denote the wing section and the panel number in that section, respectively, ${d_{Arm,j,i}}$ is the moment arm between the panel and the wing box reference line illustrated in Fig. 7, ${l_{Panel,j,i}}$ is the panel length shown in Fig. 7, and ${\tau _{Y,j,i}}$ is the shear yield stress of $ith$ panel materials in the $\;jth$ wing box section.

4.2.2 Panel buckling

Buckling is a crucial mode of failure for the wing box panels. The minimum thickness of the skin panels between the mass booms, ${t_{E.\;\;Buckling,j,i}}$ , to avoid these phenomena can be computed with Eq. (16). The equation is derived to compute panel thickness following the methods described in Ref. (Reference Grayley51) by assuming that the panels are simply supported from all edges and use isotropic or quasi-isotropic materials.

(16) \begin{align}{t_{Buckling,j,i}} \geq \;\sqrt[3]{{{q_{j,i}}{l_{Panel,j,i}}^2/\left( {\left( {3.4{{\left( {\frac{{{l_{Panel,j,i}}}}{{{b_{bay,j}}}}} \right)}^2} + 5} \right)E{\eta _B}} \right)}}\end{align}

in which, ${q_{j,i}}$ is the shear flow acting on the $ith$ panel in $jth$ wing section, ${b_{bay,j}}$ is the length of the $jth$ wing bay, ${\eta _B}$ is the plasticity factor from Ref. (Reference Grayley51).

The required thickness of a panel, ${t_{Panel,j,i,n}}$ , for one load case should be the maximum of ${t_{shearj,i}}$ and ${t_{Buckling,j,i}}$ . If the same process is repeated for each load case, then a total of ${N_{Case}}$ thickness values of each panel can be obtained, and the final size should be the maxima in this thickness list. This approach results in different panel thicknesses for each panel of the wing box. However, the panel thicknesses of the same surface of the same wing box section are assumed to be equal, excluding the spar webs. Considering this requirement, the skin thickness of the upper ( ${t_{Panel\;section,u,j}}$ ,), lower surfaces ( ${t_{Panel\;section,l,j}}$ ), the rear ( ${t_{,Web\;rear,j}}$ ) and, front ( ${t_{,Web\;frontj}}$ ) spar webs of the $j$ th wing box section can be computed by selecting the maximum thickness among the thickness lists generated from different load cases.

It should be noted that the method to compute the minimum panel thickness in this section is primarily developed for isotropic materials and it can also be used for quasi-isotropic composite materials to simplify computation. Unless quasi-isotropic composite materials are used, the local buckling methods presented in this section will be less applicable and other approaches are suggested.

The wing box structures are idealised using the methods from Ref. (Reference Megson43). The idealisation method is not suitable for the realistic global buckling analysis of the skin panels as it assumes that the bending stresses are only carried by the mass booms and the shear stresses are only resisted by the skin panels. The global buckling of the wing panels is assumed to be prevented by sufficient rib pitching. The rib masses are calculated with semi-empirical methods explained in Section 5. For test cases presented in the following sections, the half wings are divided into fixed 50 structural sections. The components are sized in those sections using the presented methods. The validity of the assumptions can be seen from the validation and verification studies of the method; however, the detailed effect of these assumptions together with different stiffener shapes on the developed method’s accuracy can be studied additionally in the future.

5.1 Refinements for composite materials

The selected ME methods are shown in Table 4 and explained in the previous section. The availabilities of the selected ME methods for composite structures are summarised in Table 4. York’s method [Reference York and Labell45] can be used for any type of material as it takes into account the effect of ‘construction\materials’ with coefficients which are discussed in Ref. (Reference York and Labell45). The mass reduction percentagesFootnote ¹ for the composite structures from Ref. (Reference Torenbeek47) are used together with the selected methods of the secondary structure masses, rib masses and wing box other masses.

In addition to the selected ME methods, other compatible methods given in Table 4 can also be combined with the proposed methods for isotropic and quasi-isotropic materials. Chiozzotto’s method [Reference Chiozzotto46] for secondary structures, Torenbeek’s methods for rib mass [Reference Torenbeek47] and wing box other structures [Reference Torenbeek10] combined with the proposed ME methods provided standard and average errors of 2.21% and 2.19%, respectively. This combination can also be selected for composite material, as Chiozzotto’s method [Reference Chiozzotto46] includes the effects of composite materials. The mass reductions of ribs and wing box other structures with composite materials are covered using the percentage method explained above [Reference Torenbeek47].

5.2 Validation with actual aircraft data

The errors of the wing MEs of the developed method and six existing methods for 13 cantilever wings are used for the validation of the developed method as shown in Table 5. The aircraft input data and the material properties are taken from Refs [Reference Chiozzotto36, Reference York and Labell45]. The validation study is one of the broadest in the literature in terms of the number, and types of actual aircraft data used, and comparison with existing methods. The aircraft’s takeoff gross weight and wing masses, ranging from the G-159 at 15.9 t and 1,655.6kg, to the 747 at 322.9 t and 39,916.1kg, respectively.

Table 5. The validation and comparison of the proposed wing mass estimation method against six other methods using the actual data of thirteen aircraft. The errors are presented by percentage (%)

One of the highest accuracy levels and the shortest computational time for wing MEs in the literature are achieved. The proposed method presented the best performance with a standard error of 1.74 and an average error of -2.22% as shown in Table 5. The computational time of the presented method is recorded during the validation studies with less than $0.1$ s for one load case using an Intel i7-11800H CPU and an octa-core 2.3GHz computer. The high accuracy level and short computational time meet the requirements of an ideal ME method presented in Section 1.1. These features also make the method suitable for MDO works of sustainable and environmentally friendly aircraft discussed in Section 1 to reduce the adverse environmental effects of air transport.

The effect of the tail generated lift (discussed in Section 3.3) on the total wing mass is investigated with a sensitivity study using Boeing 737 aircraft as a baseline. Figure 9 illustrates the variation in calculated total wing mass as the tail effect percentage increases from 0% to 20%. The results show a minimal impact on the wing mass, with only a 2.6% increase at the upper bound of 20% tail effect. The small changes confirm that the tail effect assumption does not significantly influence the overall wing mass, supporting its validity in preliminary design studies as explained in Section 3.3.

Figure 9. The change in the calculated total wing mass of Boeing 737 with the change of tail effect.

5.3 Evaluation and comparison of the model sensitivity

The proposed model’s sensitivity was previously demonstrated by Taflan et al. [Reference Taflan, Smith and Loughlan42] studying wing box mass variations with the change in aircraft design parameters. In that study, the validation of the model’s structural sizing and semi-empirical aeroelastic mass penalty approaches was also validated using the high aspect ratio strut-braced wing configuration from NASA’s SUGAR project [Reference Bradley, Allen and Droney52]. The SUGAR project employed higher-fidelity methods, including detailed FEA, computational fluid dynamics (CFD), and nonlinear aeroelastic analysis, providing a robust benchmark for comparison.

In this section, the presented wing ME method’s sensitivity to aspect ratio and sweep angle variations is evaluated in more detail by comparing with the existing methods. The Boeing 720 aircraft, detailed specifications of which are provided in Ref. (Reference York and Labell45), was used as a baseline aircraft for this study. The baseline aircraft’s aspect ratio and wing leading edge sweep angle are taken as 7.03 and 35.77 degrees, respectively. These variables were manipulated to construct the carpet plots while maintaining other variables constant, such as the wing taper ratio and wing surface area. Instead, the variation in aspect ratio was achieved by adjusting the wingspan and wing root chord values.

Figure 10 clearly illustrates that all methods’ wing mass estimates are similar for the aspect ratio values less than 10. This unity can be attributed to the empirical nature of existing methods, which rely on historical data of aircraft with moderate aspect ratio wings. However, a divergence in ME becomes more evident beyond an aspect ratio of 10 where the wing configurations reach beyond the range of the data on which existing methods were built. The estimations from existing methods either increase linearly or with a diminishing gradient, while the new method reflects an increasing gradient as aspect ratio values exceed 10. The carpet plots in Fig. 10 illustrate the higher sensitivity of the proposed method to alterations in the design variables and it also meets the expectations of the increasing mass growth in high aspect ratios.

Figure 10. The two surface plots illustrate the variations in wing sweep, aspect ratio and estimated wing mass.

5.4 Limitations and future improvements

The proposed method has been developed in line with the ideal ME requirements outlined in Section 1.1. It has a balance between accuracy and complexity. The method’s advanced accuracy has been demonstrated through one of the most detailed validation studies in the literature, particularly for moderate aspect ratio wings, achieving high accuracy and short computational times. The verification of the method for high aspect ratio wings has been conducted, accompanied by a sensitivity analysis. Additionally, the proposed methods’ accuracy for wing box ME of HAR strut-braced wings was previously validated in Ref. (Reference Taflan, Smith and Loughlan42).

The proposed method has been designed to address the effects of dry wings, distributed, hydrogen, and electric propulsion systems on total wing mass, targeting a significant gap in the existing literature. The method also offers valuable insights for the preliminary structural design and sizing of wings, benefiting from its detailed geometric model.

In the development of advanced quasi-analytical methods for conceptual and early design stages, the trade-off between model complexity and accuracy is critical. While higher fidelity analytical aeroelastic methods can enhance the accuracy of aeroelastic analyses, they may lead to reduced accuracy in overall wing MEs due to necessary simplifications in the geometric model aimed at reducing computational time and complexity. To address these challenges, the proposed method integrates a semi-empirical approach to compute aeroelastic mass penalties within the quasi-analytical methods presented in this study. The sensitivity of the semi-empirical aeroelastic mass penalty method to increasing aspect ratios is verified in the following section. However, more detailed evaluation and validation of the method for high to ultra-high aspect ratio wings is recommended for future research, particularly as more data for these wing configurations become available.

The method has been developed with the intention of applying it to composite materials, as discussed in Section 5.1. The methods validation for a wing box with advanced quasi-isotropic materials has been previously shown in Ref. (Reference Taflan, Smith and Loughlan42) using mass data of a HAR strut-braced wing. Due to the limited availability of data on composite wing boxes from existing transport aircraft, more detailed validation of the method with composite materials is suggested for future studies.

The detailed geometric model, coupled with load and stress analysis, facilitates the application of structural sizing procedures for various structural elements of a wing box, including each skin panel, spar web, spar cap and stringer. These capabilities render the method suitable for the detailed design and optimisation of composite structural elements. Although these applications were beyond the scope of the current study, they are planned for future research.

The presented method achieved an effective balance between accuracy and complexity, making it highly suitable for preliminary and conceptual wing design, particularly under the influences of distributed propulsion, hydrogen fuel and dry wing configurations. While its advanced accuracy for moderate aspect ratio wings has been rigorously validated, further research is necessary to refine its application to high and ultra-high aspect ratio wings, as well as to composite materials. Future work will focus on enhancing the aeroelastic analysis component and detailed applications for composite wing structures.