Structural sizing and mass estimation of transport aircraft wings with distributed, hydrogen, and electric propulsions

Abstract

Current literature offers limited mass estimation methodologies and their application in the conceptual or preliminary design stages of moderate to high aspect ratio wings with electric, hydrogen or distributed propulsions. This study presents the development and application of a quasi-analytical wing mass estimation method to address this limitation. The proposed method is distinguished from the existing mass estimation methods by its expanded realistic load cases, sensitivity to several design parameters, improved accuracy with short computational time and capabilities for future applications. To achieve these features, new geometric models are introduced; 483 load cases including symmetric manoeuvre, rolling, and combined cases are covered following airworthiness requirements; the structural elements are idealised and sized with strength and buckling criteria; existing methods are evaluated and integrated cautiously for secondary structures and non-optimum masses. A computation time of 0.1s is accomplished for one load case. The developed method achieved the highest accuracy with an average error of -2.2% and a standard error of 1.8% for wing mass estimates compared with six existing methods, benchmarked against thirteen wings of different aircraft categories. The effects of engine numbers with dual- to 16-engine setups and the dry wing concepts on the wing mass are investigated. The optimised number of engines and their locations decreased the wing mass of the high aspect ratio wing significantly. In contrast, the dry wing design increased the wing masses of all baseline aircraft. The future applications and improvements of the presented method in novel configurations and multidisciplinary designed optimisation studies are explained.

Nomenclature

ATAG

Air Transport Action Group

ME

mass estimation

FEA

finite element analysis

MMI

mass moment of inertia

MDO

multidisciplinary design optimisation

EASA

European Union Aviation Safety Agency

ND

non-dimensional

1.0 Introduction

The aviation sector stands at the intersection of remarkable growth and alleviating adverse environmental impacts. European Environment Agency et al. [1] indicated an 8% surge in the total number of flights between 2014 and 2017, with 42% growth projections from 2017 to 2040. In collaboration with the Air Transport Action Group (ATAG), the Flightpath 2050 report set ambitious targets for the future of aviation [2]. The report presented the goal of a reduction of 75% in $C{O_2}$ emissions per passenger kilometer, a 90% in $N{O_x}\;$ emissions, and a 60% in perceived noise compared to a conventional aircraft in 2000. To meet these objectives, contemporary research has focused on high aspect ratio wings [3] with distributed propulsions [4] utilising hybrid electric or hydrogen fuel technologies [5, 6] which promise enhanced performances and alleviate the adverse environmental effects.

The research for effective designs of the configurations with the potential of reducing emissions necessitates advanced conceptual and preliminary design methods that are computationally efficient, accurate, and sensitive to design decisions [7, 8]. In that regard, structural sizing and mass estimation (ME) of aircraft wings emerge as paramount disciplines from the initial to detailed design stages [9, 10]. Mass estimation greatly influences the aircraft’s performance, fuel consumption and operating costs [7, 11, 12] . Hence, having a clear understanding and accurate assessment of how the mass changes in response to the design decisions is essential for meeting the design requirements of any conventional or novel aircraft.

The distributed propulsion application has gained popularity with increasing interest in electric propulsion systems as this application’s promising results have been debated in aircraft performance [13, 14]. The investigation of the impact of both the quantity and spatial arrangement of distributed propulsions on wing mass necessitates research which must be conducted concurrently with the examination of complementary sub-disciplines such as aerodynamics and performance for precise assessment and refinement.

Similarly, the dry wing concept has emerged because of growing interest in zero-carbon hydrogen or electric-powered aircraft. Although the batteries of electric-powered aircraft can be located inside the wing box for load alleviation, there has not been a feasible approach for locating the cryogenic hydrogen tanks in those spaces. Viable solutions for such cryogenic hydrogen tanks might be locating them inside the fuselage or attaching them to the wings. The effect of these design solutions on the wing mass should be studied with validated ME methods to achieve accurate aircraft conceptual and preliminary designs.

However, a broad literature review indicated a lack of wing ME methodologies and their applications for the wings from moderate to high aspect ratios using distributed, hybrid-electric, or hydrogen propulsion systems. Although some existing wing ME methods might be used in such applications with minimal modifications, there has been still a gap in the literature presenting their validity and application. Hence, this study aims to bridge this gap by developing and implementing an advanced wing ME method considering the requirements of wing configurations with advanced technologies.

1.1 Ideal mass estimation method

The balance between the accuracy and complexity of a ME method is critical [15, 16] and the wing ME methods have been classified considering this aspect. Previously, Dababneh and Kipouros [7] reviewed the classifications of the ME methods in the literature and suggested classifying the wing ME methods into three main groups consisting of empirical-based, analytical-based and finite element analysis (FEA) based methods. Empirical and analytical methods can also cover semi-empirical and quasi-analytical methods, respectively. This classification approach [7] is employed in this study. Its details and other classifications can be found in the Ref. (7).

Although the empirical methods are preferred at the conceptual aircraft design stage for their ease of use, and short computational time, they might lack in sensitivity and accuracy. They might also fail in the wing ME of novel configurations as they are generally built on the historical data of conventional aircraft. Their development for novel configurations is usually not possible due to the limitation in available data.

FEA-based methods can be developed to cover novel configurations; however, resorting to high-fidelity FEA-based methods can become impractical for conceptual and preliminary design stages due to their complexity in implementation and long computational time as previously reported in Refs (17–20). Law [21] stated that ‘If an analytical solution to a mathematical model is available and is computationally efficient, it is usually desirable to study the model in this way rather than via a simulation.’

The features of an ideal wing ME method for the conceptual or early preliminary design stages of novel configurations can be listed as follows [22].

• • Sensitivity to increasing aspect ratios and other design parameters

• • Short computational time

• • High accuracy level with low average and standard errors

• • Suitability for its integration in multidisciplinary design optimisation (MDO) environments

• • Rely on physics rather than statistics to be able to cover novel configurations

• • Sensitive to the effects of emerging propulsion technologies and alternative fuels

Considering these features, the need arises for analytical or quasi-analytical methods, which are positioned between the empirical and finite-element-based methods in terms of complexity in their application [15, 16]. They might be developed to be sensitive to the studied aircraft design parameters and novel technologies, such as propulsion technologies discussed above while achieving higher accuracy levels than empirical and FEA-based methods [22].

1.2 Introduction to the advanced quasi-analytical method

The proposed quasi-analytical wing ME method blends analytical and empirical methods for optimum accuracy and computational time. Empirical methods are used in the ME of the secondary structures and non-optimum masses. On the other hand, new analytical-based ME approaches are introduced for the load-carrying structures as their structural sizes and masses change with different loading conditions caused by the novel configurations.

The author critically reviewed the existing ME methods and identified their strengths and limitations in the design studies of the novel configurations under scope. Four key features of existing ME methods are identified to improve the proposed ME method in its application area, accuracy and sensitivity compared to existing methods. These features are grouped under geometric models, broad load cases, detailed structural sizing, and architecture, which are discussed in the following sections.

1.2.1 Geometric models

The proposed ME method is developed to be sensitive to the variation of the sweep, thickness to chord ratio, dihedral, taper ratio and aerofoil types along the wings. The importance of the model’s sensitivity to these design variables in aircraft design has been discussed and shown in the literature [12, 23, 24] . Hence, a discrete geometric model approach (as shown in Fig. 2) is utilised in this study. It allows unlimited definition of wing sections in the spanwise direction, a capability inspired by earlier studies [8, 25, 26] . This model enables assigning the wing and wing box design parameters at each section or enables computation through interpolation. The model’s discrete architecture also caters to intermediate and low-fidelity level aerodynamic and structural analysis methods [27, 28]. It can also be used for the design and optimisation of dry wings and wings with distributed, hydrogen and electric propulsion systems.

In the literature, analytical or quasi-analytical methods have utilised various wing box models such as the I-beam [29], idealised [8], double plates [30], and hexagonal [26] for wing ME. The proposed method introduces an advanced idealised wing box model with a cross-section that can morph between diamond, rectangular and hexagonal shapes. This model is more representative than I-beam and double plate models, allowing for more detailed sizing of critical structural components, including stringers, skin panels, spar webs and caps.

1.2.2 Broad load cases

A thorough literature review highlighted a lack of clear procedures for defining the number and types of critical load cases necessary for the accurate structural sizing and ME of conventional and novel wings. As shown in Table 1, the range of load cases varies across studies. To comprehensively address this, the new method includes symmetric, rolling and combined load cases, following the airworthiness guidelines [31, 32]. This approach resulted in 483 total load cases, one of the broadest ranges in the literature in structural sizing and ME at the initial design stage. This comprehensive approach to load cases enhances the accuracy of the proposed method by reducing uncertainties and increasing the likelihood of capturing the most critical load cases, which directly influence structural sizing and ME for both conventional and novel configurations. Future studies will explore critical load cases using a systematic approach, integrating the methods proposed in this study to identify a reduced set of critical load cases that do not compromise the accuracy of the method, thereby further decreasing overall computational time.

Table 1. The total number of load cases considered in the previous studies in the literature

1.2.3 Detailed structural sizing

Following the outcomes of previous research [25, 26, 41] , an idealised hexagonal wing box model is developed. It is sensitive to the effect of the internal structural design of the wing box, such as spar webs, caps, stringers and wing box cross-section, on the wing box mass. All these structural elements of the wing box are sized to resist buckling, shear and normal stresses under all load cases considered. This detailed sizing procedure also enables designing and optimising conventional and dry-wing box structures and wing planforms for conventional and distributed propulsion systems and studying the effect of those structural elements on the wing mass.

1.2.4 Model architecture

The model’s sub-modules are designed and developed considering their compatibilities with each other (as illustrated in Fig. 1), the requirements of an ideal ME method discussed above especially focusing on the accuracy and computational time, its direct application for cantilever wings with distributed propulsion, emerging fuel systems, its future applications in the design, optimisation and sizing of the composite structural elements of a wing box, its improvement and application for strut-braced wings employing novel technologies and finally the application in multi-disciplinary design optimisation environments.

Figure 1. N2 Diagram of the complete wing mass estimation model.

Figure 2. 2D top and front views of the wing geometric model.

The model’s suitability and initial applications for strut-braced wings are previously presented in Ref. (42). In this paper, the model’s development and performance evaluation in accuracy and computational time are presented and discussed. Its direct applications with distributed propulsions, hybrid-electric and hydrogen fuel systems are shown. The model is designed to ensure its suitability for future applications. Its improvement and application for novel wing configurations with advanced technologies are aimed at the following studies.

1.3 Method development and implementation

The architecture of the developed ME method is illustrated in Fig. 1 with an N2 chart. The method consists of five main modules: initialiser, geometric module, load calculator module, structural sizing module and ME module.

Initialiser module: It provides a multi-fidelity feature and reduces the complexity by providing the initial estimates or design parameters. The user can surpass the initialiser module if data is provided.

Geometric module: The geometric module is required to generate the 2D and 3D coordinates of the complete aircraft by using the input design parameters. The geometric module can be connected to a visualisation tool to generate aircraft illustrations. The module is also used to compute the volumes and surface areas of the aircraft components. The geometric models are explained in Section 2.

Load calculator module: The load calculator module consists of the mass moment of inertia (MMI), flight case generator, aerodynamic and loading actions modules. The input parameters of each of these submodules are indicated in Fig. 1. This module computes the shear force and moment distributions used by the structural sizing module. The module is created by implementing the methods presented in Section 3.

Structural sizing module: This module sizes the structural elements of the wing using the procedures presented in Section 4. The dimensions of the sized structures are transferred to the ME module.

Mass estimation module: The ME module consists of the ME formulas developed in this study and the existing methods from the literature to compute the total wing mass. The detail of this module is explained in Section 5.

This paper consists of seven main sections. Sections 2–5 are explained above, and the results and discussion with validation and preliminary results of the developed ME model are presented in Sections 5 and 6. Finally, Section 7 provides the concluding remarks and future works.

2.0 Geometric model

A complete aircraft geometric model is developed considering the requirements of the structural sizing and ME module for the wings with advanced technologies as discussed in Section 1. Hence higher fidelity level is chosen for the wing compared to the other geometric models.

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. (43), 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.0 Load analysis

Although airworthiness requirements are subject to change for novel aircraft design, they can still provide an essential guideline for load analysis and structural design. A sub-module is developed to generate load cases and flight envelopes following CS.25 requirements [31].

CS.25 [31] requires designing airplanes under different load cases generated by the combination of different manoeuvre directions, altitudes, fuel and payload capacities and speed variations, as shown in Table 2. As a result, 195 symmetric load cases are created, including limit manoeuvring and discretised gust load factors. The gust cases are generated for all combinations of manoeuvre directions, altitudes, fuel and payload capacities and speed variations listed in Table 2. The limit manoeuvring load factors are generated for all combinations of manoeuvre directions and speed variations with the design maximum take-off weight.

Table 2. The parameters in the table are combined to generate the symmetric load cases

Figure 5. Location of the stringers, spars and wing sections.

In addition to 195 cases, more cases arise from rolling conditions and the combination of rolling and symmetric manoeuvring. The rolling cases cover arrest, initiation and steady cases, as shown in Table 3. Hence, a total of 483 load cases are generated following CS.25 requirements [31], including discretised gusts, symmetric manoeuvre, rolling, combined cases, one taxi bump case (i.e., inertial loads only) and one landing case (i.e., aerodynamic and inertial loads) as the default setting of the tool.

The external load transferred from each of the wing mounted landing gears to the wing ( ${F_{LG,Taxi}}$ ) is computed with $ - 0.9{W_{TOGW}}{n_{Taxi}}/{N_{MLG}}$ and $ - {W_{TOGW}}{n_{Landing}}/{N_{MLG}}$ for taxi and landing case, respectively. ${N_{MLG}}$ is the total number of main landing gears and ${n_{Taxi}}$ is taken as +2. For the fuselage mounted landing gears, there is no load transferred from landing gears to the wing structure. In the landing case, the lift generated by the wing is assumed be equal to the aircraft mass in landing and ${n_{Landing}}$ is computed with ${n_{Landing}} = 1 + v/\left( {g \cdot t} \right)$ where $v$ is taken as the maximum deceleration speed of 3.05m/s [31], $g\;$ is gravity (∼9.81 $\;{\rm{m}}/{{\rm{s}}^2})$ , and $t$ is the deceleration time to take the landing gear absorptions into account simply. For $t$ of 0.3 $s$ , ${n_{Landing}}$ is calculated as 2.04 as the default of the developed tools.

Rolling cases are generated following CS.25.349 and AMC 25.147 is utilised for computing roll timing [31]. Rolling load calculation procedure is listed as follows.

1. 1. Required angular rate of roll (roll velocity, $\;p$ ) is computed.

2. 2. The required equivalent aileron (explained in Section 3.3.1) deflection ( $\xi $ ) is computed to achieve the required $p$ .

3. 3. Unless aileron deflections meet the roll timing requirement, $\xi $ is updated.

4. 4. Roll acceleration ( $\dot p$ ) and the aileron aerodynamic loads are computed (see Section 3.3.1).

5. 5. Inertia loads are computed using $\dot p$ and aerodynamic damping is computed as explained in Section 3.3.2.

6. 6. Following the above steps, the total rolling load is computed for each roll phase of arrest, initiation and steady.

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. [22] 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 [44] 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)$ [44]. 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 [44] 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.0 Structural sizing

Using the load analysis method explained in Section 3 and the geometric model in Section 2, stresses on each structural element of the wing can be computed and sized for strength and buckling constraints.

As prescribed by CS 25.305(a) [31], the aluminium alloy wing components are sized for yield stress criteria under the limit loads, ensuring that harmful permanent deformation does not occur. The composite materials are sized under the ultimate loads for the failure strength of the material following the criteria in CS 25.305(a) [31]. The ultimate loads are computed as 1.5 times of the limit load, considering the safety factor as prescribed by the CS-25.

4.1 Mass boom sizing

For preliminary structural sizing, the structural idealisation approach is used from Ref. (43). 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. (43). 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. (51) 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. (51).

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. (43). 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.0 Mass estimation

In this study, analytical methods are developed to compute the shear and bending-resistant material masses of the wing boxes. The empirical and semi-empirical methods are selected from the literature for the ME of the ribs, secondary structures and wing box other masses considering their compatibilities with the proposed analytical methods. The grouping of the component masses employed in this study is presented in Fig. 8. The secondary structure masses cover the leading edge, trailing edge, landing gear doors mechanism, ailerons, elevons, slats, leading and trailing edge flaps, spoilers, wing speed breaks, winglets, wing folding and pivot [45]. Moreover, the wing box’s other masses cover aeroelastic penalty, mounting and non-optimum masses [10, 47]. The non-optimum masses include the mass of the mounting components, and penalties due to sheet taper, joints in the panels and cut-outs [47].

Figure 8. The components of the wing mass and the type of methods used in their estimation.

The compatibilities of four different ME methods of wing secondary structures [11, 45–47] and two methods of rib masses [10, 48], given in Table 4, with the proposed quasi-analytical wing box ME methods were tested. The Torenbeek’s method [47] for the ribs, and York and Labell’s method [45] for the secondary structures combined with the proposed analytical methods provided the lowest standard and average errors of 1.8%, and -2.2%, respectively in wing mass estimates of thirteen real aircraft compared to other empirical and semi-empirical ME method combinations shown in Table 4. Hence, the Torenbeek’s and York’s [45, 47] semi-empirical methods are chosen as a default of the developed ME approach for the ribs and the secondary structures, respectively. Furthermore, the ME methods of Torenbeek [10] are selected for mounting masses. The non-optimum masses and the aeroelastic mass penalties are covered by employing Torenbeek’ method [47]. The limitations and capabilities of the methods are discussed in following sections.

Equation (11) was given in previous sections to compute the cross-sectional area of the mass booms (bending-resistant materials), which represents the spar caps and stringers. The length of the mass booms can be obtained from the geometric model; hence, the volume and then the mass of the materials can be computed using the material density. Similarly, the total mass of the shear-resistant materials can be computed using the thicknesses of the front and rear spar webs, upper and lower skins (explained in Section 4.2), material density and geometric models. Finally, the total wing mass is computed by summing all the mass components of ribs, secondary structures, wing box other masses, and shear- and bending moment-resistant materials shown in Table 4 and Fig. 8.

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 [45] 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. (45). The mass reduction percentages ¹ for the composite structures from Ref. (47) 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 [46] for secondary structures, Torenbeek’s methods for rib mass [47] and wing box other structures [10] 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 [46] 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 [47].

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 [36, 45]. 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. [42] 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 [52]. 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. (45), 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. (42).

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. (42) 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.

6.0 Results and discussion

The application of the proposed method in ME of aircraft wings with distributed propulsion and emerging fuel systems of hydrogen and electric is presented and discussed in this section using two baseline aircraft (illustrated in Fig. 11) of Boeing 737 and CANT (a cantilever, strut-removed variant of the NASA-SUGAR aircraft without any other change in design parameters). While the wing planform of CANT is the same as the wing of NASA-SUGAR, the internal structural components of the CANT aircraft are sized with the proposed method. It should be noted that the CANT wing might not be the optimum design in terms of overall aircraft performance, and this high aspect ratio design is only selected for comparison purposes of this study. The design parameters of the baselines were previously presented in Ref. (42). The advanced material properties of the quasi-isotropic composite ² used in CANT aircraft can be accessed from Ref. (52).

Figure 11. 2D top views of the baseline Boeing 737 and CANT aircraft with moderate and high aspect ratio wings, respectively.

The calculated wing masses of the baseline aircraft and their breakdown are presented in Table 6. The higher proportion of the bending material in the higher aspect ratio wing due to higher bending moments can be observed. The computed aeroelastic mass penalties and their ratio to the total wing mass were 79.51kg - 1.74% and 474.26kg - 5.65% for Boeing 737 and CANT baseline aircraft, respectively. The sensitivity of the semi-empirical aeroelastic methods to increasing aspect ratios is verified. Future studies can also consider upgrading the aeroelastic analysis component of this method with more sophisticated approaches to achieve better sensitivity to aeroelastic constraints.

Table 6. The mass breakdown of the baseline aircraft, and their studied variations

6.1 Effect of dry wing designs

Conventionally, wing boxes, fuselages and drop tanks have been the primary repositories for fossil fuels in conventional aircraft designs. Markedly, hydrogen-fuelled aircraft are inclined towards fuselage-based fuel storage due to the inherent complexities and limitations of storing hydrogen. Therefore, the wing box spaces without fuel tanks led to the emergence of so-called ‘dry wings’.

Formerly, the effect of the fuel amount stored in the wing on the wing box masses of four baseline aircraft with different aspect ratios was studied in Ref. (42). The spar positions of a dry wing were also investigated as they have more space to be located in the absence of the internal wing fuel tanks. The current study extends the work on dry wings by covering the total wing mass and considering the effect of wing drop tanks.

For this study, each of two baseline aircraft model are configured in two ways based on fuel storage, as presented in Table 7, resulting in three variations per aircraft model. The ‘Con’ designation indicates the baseline aircraft with their conventional concept where fuel is stored in the wing box and fuselage. The ‘Dry’ designation signifies that all fuel is stored within the fuselage. The ‘1Drop’ configuration represents a hybrid approach where half of the fuel normally stored in the wing box is moved to wing-mounted drop tanks, with the remaining fuel stored in the fuselage.

Table 7. The description of the investigated aircraft concepts

Table 6 shows the estimated wing masses of the six configurations (shown in Table 7) using the developed ME method. The proposed ME method conforms to the expectation that dry wings should exhibit greater mass due to the reduced inertia relief. The HAR wing of CANT, Dry, Comp configuration shows a 9.57% increase in total wing mass compared to its baseline, driven primarily by a 22.5% increase in bending-resisting material mass. This is remarkably higher than the 3.42% increase in the MAR wings of the 737, Dry, Al configuration, where the bending-resisting material mass increased by 15.7%. The HAR wing requires more reinforcement to compensate for the loss of bending relief from the removed fuel, leading to higher sensitivity.

Moreover, mounting the drop tanks to the dry wings led to a marginal reduction in the wing mass of both the moderate and high aspect ratio wing configurations, suggesting that a cryogenic drop tank could partially compensate for the absence of load relief in hydrogen-powered aircraft. For the dropped tank configuration, the HAR wing of CANT, 1Drop, Comp shows only a 2.91% increase in total mass, aided by a 12.7% decrease in bending-resisting material compared to its dry wing variant. And the MAR wing of the 737, 1Drop, Al, shows a 0.61% increase compared to the baseline, where bending material mass decreased by 10.9% compared to its dry wing variant. The drop tanks contributed to the bending relief in both contributions.

6.2 Single disciplinary optimisation of engine numbers and locations for minimum wing mass

This section presents a single disciplinary optimisation study which only covers the constraints of the developed ME method to achieve a minimum wing mass objective. It should be noted that this study is performed to investigate different numbers of engines on wing mass. Hence the wing mass is selected as an optimisation objective, and the findings do not reflect optimised wing design with the optimisation objectives for overall aircraft performance.

The effects of the total engine numbers from 2 to 18 engines and their spanwise locations on the wing mass of moderate and high aspect ratio wings are investigated using a global optimiser called differential evaluation [53] with the developed ME method. The two baseline Boeing 737 and CANT aircraft are used for this study and their optimised variants are named as Boeing 737,Dist and CANT,Dist, respectively. Their top views are illustrated in Fig. 11 with the original engine numbers and locations on the half wings. It was assumed that the total engine masses and the thrust of the baseline aircraft are distributed equally between the selected number of engines by the optimiser. The optimiser’s objective was set to converge to the minimum wing mass by varying the number of engines from 2 to 16 and their spanwise locations.

Table 8 shows the actual and optimised engine counts and non-dimensional spanwise engine positions as a fraction of the half wing length of Boeing 737 and CANT aircraft. The upper and lower bounds of engine numbers were 2 and 16, respectively. Both baseline configurations have a total of two engines. The optimiser found two and four-engine configurations for the minimum wing mass of Boeing 737 and CANaircraft, respectively. It can be observed from Table 8 that the engine locations are moved towards the wing tip in both configurations to achieve the minimum wing mass objectives. The results correlate with the previous findings in Refs [42, 54].

Table 8. The original and optimised engine counts and their ND spanwise positions as a fraction of the half span length

The comparisons of the actual and optimised wing masses of both aircraft are presented in Table 6. A significant wing mass reduction of 17.16% is achieved in CANT aircraft with a high aspect ratio wing. In contrast, the mass alleviation was minimal for the optimised Boeing 737 aircraft with a moderate aspect ratio wing. These results indicate the significant sensitivity of HAR wings to engine numbers and locations, and this sensitivity can be explained with the high bending relief and reduced bending resisting materials shown in Table 6. The results not only illustrate the trend for achieving the minimum wing mass objective by varying the engine location but also explore the impact of engine count along with other variables.

It should be noted that the optimum number of engines and locations would be different if multidisciplinary constraints and other optimisation objectives of overall aircraft performance were considered; however, the current study did not aim for such an optimisation. Additionally, the effect of the distributed propulsions on the aerodynamic load distributions is not covered in this study due to the use of simple semi-empirical aerodynamic models. Future studies are recommended to include and explore these effects. It is advisable to validate these outcomes by comparing them with higher fidelity results before their incorporation, as the current results aim to provide insights primarily into conceptual design rather than detailed design.

6.3 Impact of distributed propulsion on wing mass: a sensitivity analysis

In this section, a sensitivity analysis examines the effect of varying engine counts on wing mass using the baseline configurations from the previous study. Unlike the earlier optimisation, which adjusted both engine count and spanwise location, this study fixed the engines at equal spanwise distances and mid-chord positions. The engine masses were computed by evenly dividing the total engine mass of the baseline aircraft by the engine counts. Figure 12 shows minimal wing mass change, with a maximum variation of 1.05% in the CANT aircraft and 6.61% in the Boeing 737, contrasting with the 17.16% reduction observed when optimising both engine count and location.

Figure 12. Impact of distributed propulsion on the wing mass of the baseline Boeing 737 and the CANT aircraft.

These findings emphasise that engine location, rather than count alone, plays a pivotal role in reducing wing mass, particularly due to increased inertia relief when engines are positioned near the wing tips. This highlights the necessity of engine placement optimisation in achieving structural efficiency in the conceptual design phase.

7.0 Concluding remarks

The proposed wing mass estimation method offers a novel approach for structural sizing and mass estimation studies at the conceptual and preliminary design stage, displaying enhanced accuracy, sensitivity and computational efficiency compared to existing methods using datasets of existing transport aircraft. Its application in wing mass estimations of various configurations, including those with composite materials, hydrogen, electric or distributed propulsion systems, has been shown through modified geometric models, realistic load cases and detailed structural sizing methods.

The modified geometric models are designed considering the proposed method’s requirements in accuracy, computational speed, future improvements and its wide application area in novel configurations, parametric analysis and detailed structural optimisations. The structural sizing methods are developed considering their compatibility with the fidelity level of the geometric model. Moreover, one of the broadest load case analyses in the literature also contributes to the method’s accuracy preventing uncertainties in the determination of the dominant load cases driving the structural sizes.

The presented quasi-analytical method blends the semi-empirical and analytical mass estimation approaches cautiously to benefit from their superiority to achieve high accuracy and short computational time. Structural idealisation methods are applied to reduce the complexity of implementation while accomplishing detailed sizing of the wing structural elements including, stringers, upper and lower skin panels, front and rear spar caps and webs. The design parameters of these elements can be changed along the wing in each wing section in addition to other design parameters of the wing and the overall aircraft.

Validation against data of 13 real aircraft from various categories and companies demonstrated the method’s credibility and its highest accuracy compared to six existing methods from the literature with an average error of −2.2% and a standard error of 1.8%. This elevated level of accuracy is achieved within an exceptionally short computational time of 0.1 seconds per load case, utilising an Intel i7-11800H CPU on an octa-core 2.3GHz computer, a notable benchmark in the literature. The verification and validation studies also increase the proposed method’s credibility as a valuable preliminary structural sizing tool, aiding in early design stages.

Moreover, the model’s enhanced sensitivity to novel aircraft configurations, particularly those with high aspect ratio wings, is also verified by comparing the sensitivity results with those of five other methods using surface plots. The sensitivity analysis underscores the adaptability of the proposed method to diverse design requirements.

One of the first applications for the mass estimation of dry wing concepts for hydrogen or electric aircraft in literature is presented using the proposed method. The results highlighted the possible mass penalties of these concepts due to the absence of load relief on the wings. The potential alleviation strategies, such as wing-mounted cryogenic drop tanks, are proposed and their effects on the wing mass reduction of moderate and high aspect ratio wings are presented. The optimisation studies using the proposed method with minimum wing mass objectives reveal insights into engine placement and distributed designs from dual to 16 engine setups on moderate and high aspect ratio wings. The high sensitivity of the high aspect ratio wings to the number of engines and their locations is observed.

Looking ahead, integration into multidisciplinary design optimisation environments, more detailed design and optimisation studies with composite materials and further improvements and validations for high and ultra-high aspect ratio wings using higher fidelity aeroelastic models are promising avenues. A systematic approach will be applied to identify a reduced set of critical load cases using the proposed method, aiming to decrease overall computational time. Additionally, future studies will explore applications to strut- and truss-braced wings, expanding the method’s utility across a range of innovative aircraft configurations.