On the multi-fidelity approach in surrogate-based multidisciplinary design optimisation of high-aspect-ratio wing aircraft

Abstract

The reduction of computational costs in the context of the Multidisciplinary Design Optimisation of a typical medium-range aircraft was investigated through an assessment of active constraints and the use of multi-fidelity models-based estimation of drag and structural stress. The results show that for this problem, from the set of considered constraints that includes flutter boundary, the active constraint is a 2.5g pull up Maximum Take Off Weight. Results show that the multi-fidelity approach reduced the required high-fidelity aerodynamic number of evaluations, for both drag assessment and stress assessment with sufficient level of accuracy for the former and conservatively for the latter. Further computational cost reduction can be achieved using a surrogate model based Multidisciplinary Design Optimisation. The best configuration attained shows an Aspect Ratio increase of 16%, a reduction of 4.5% in fuel consumption and wing structural weight increase of 2.7% relative to a predefined baseline configuration.

Nomenclature

A

real constant (k/g)

AOA

Angle-of-Attack

AR

Aspect Ratio

B

real constant (N s lb-1)

c

fuel consumption (lb/N/s)

CFD

Computational Fluid Dynamics

DS

dive speed (m/s)

EBM

Equivalent Beam Model

E[I]

expected improvement

f

generic function

FM

flutter margin

FS

flutter speed (m/s)

FSI

Fluid-Structure-Interaction

HARW

High Aspect Ratio Wing

k

correlation value

L/D

Lift-over-Drag ratio

MDA

Multidisciplinary Analysis

MDO

Multidisciplinary Design Optimisation

MTOW

Maximum Take-Off Weight

PM

Panel Method

r

correlation vector

R

correlation matrix

RANS

Reynolds-Averaged Navier-Stokes

$ \hat{s} $

mean squared error

SF

safety factor

SBW

Strut-Braced Wing

w

structural mass (kg)

x

design point

$\boldsymbol{x}_\boldsymbol{i}$

design variable

y

evaluated response

$ \hat{y} $

predicted response

Greek symbol

$\hat{\sigma}$

stress (MPa)

$\sigma_{yield}$

yield stress (MPa)

$\theta$

Kriging hyperparameter

$\phi$

distribution density function

$\Phi$

probability density function

$\hat{\mu}$

regression mean

Subscripts

i

index

1.0 Introduction

Nowadays, conceptual aircraft design requires assessment tools other than statistical data or empirical formulas to determine the best conceptual candidates to proceed with for detail design. This requirement is prompted by the emergence of new aircraft configurations that rely on design space expansion to comply with environmental impact and cost reduction requirements. As a consequence, different geometries that are intended to improved overall flight efficiency are being researched, including the High-Aspect-Ratio Wing [1, 2], Braced Wing [3, 4], Joined Wing [5, 6], Blended Wing Body [7, 8, 9] and other aircraft configurations [10, 11].

Some of these configurations reduce induced drag in cruise flight conditions based on an Aspect Ratio (AR) increase. This is the case for the High Aspect Ratio Wing (HARW), Braced Wing, such as strut-braced (SBW) or truss-braced wing, and Joined Wing configurations, among others. The downside is the structural weight penalty due to an increase in the root bending moment [12], increasing the maximum take-off weight (MTOW) for a fixed payload. Thus, these designs tend to be more prone to aeroelastic and flexibility effects caused by relatively larger deformations of the wing structures compared to conventional wing designs.

If the deformations are large enough both aerodynamic and structural non-linear effects become more prominent. A 2007 review of the non-linear aeroelastic prediction methods [13] concludes that these should be used in early design stages to identify potential problems while in Patil et al. [14] one can find a comparison of linear and non-linear aeroelastic analyses results for very high-aspect-ratio wings showing significant differences between the two.

Other studies [17] indicate that using aeroelastic instabilities as flutter as a driver for wing design may be diluted due to the possible existence of high amplitude limit cycle oscillations well below the flutter speed in very high aspect ratio wings. Furthermore, flutter boundary predictions based on high-fidelity CFD results are reported as more conservative that the traditional potential flow theory based results [15]. Non-linear effects of flexibility in aircraft trim [14, 18] and in closed loop control [19] have also been subject of research.

Since these large deformations influence the flutter boundary [15], they should be considered in an early stage of the design process, which can be done in a Multidisciplinary Design Optimisation (MDO) environment [16]. As the optimal designs tend to reduce structural weight and stiffness, deformation increases and the flutter margin is further reduced, a trend that is shown in several studies [26, 27]. A possible means of mitigation of such problem would be the incorporation of flutter suppression devices in the design [20].

When optimising a HARW or SBW configuration, for example, it is possible that not only structural stress constraints of the design, but also the flutter boundary can be an active constraint, although an assessment in an early design stage can be challenging [16]. Fluid Structure Interaction (FSI) analysis can be the first step to evaluate the effects of deformation in the flutter boundary. By solving the coupled aero-structural problem [21] and obtaining the aeroelastically deformed shape for a particular flight condition and the corresponding aerodynamic coefficients [22, 23] one can use the stiffness and aerodynamics of the deformed structure to estimate the flutter speed with current methods [24, 25]. Whether or not non-linearities are included in the deformed structure stiffness matrix depends on the computational setup available for the calculation.

In an MDO procedure where structural stress and flutter margins are constraints and flight efficiency are the, or part of, the objective function, multiple computational fidelity levels can be used to compute the constraints and the objective function. Lower fidelity models can be used to predict global stress levels and structural mass with a degree of accuracy due to their good lift predictions, namely lifting line models including non-linearities [28], expanded to swept and non-planar wings [29] or integrated in multi-fidelity analysis frameworks [30].

Low-fidelity aero-structural models are also extensively used in flutter speed calculations. Examples of two and three-dimensional low-order models used for flutter boundary calculations based on small unsteady disturbances about a known steady flow CFD solution can be found in Opgenoord et al. [31] or Latif et al. [32]. In Opgenoord et al. [33] a comparison of analytical methods for flutter and post flutter calculations in a wing section with three degrees of freedom is presented and a hybrid approach involving analytical and Finite Element models for flutter calculations is presented in Sudha et al. [34].

For accurate drag assessment Reynolds Averaged Navier-Stokes (RANS) based Computational Fluid Dynamics (CFD) computations are typically relied upon to capture viscous and compressible effects, which may also have influence in the aeroelastic behaviour. Phenomena as transonic flutter, limit cycle oscillation and the transition between the two are captured only with high-fidelity aerodynamic models [35], as is stall flutter [36] and fluid mode instabilities [37].

The MDO procedure should be as fast as possible; therefore, utilising lower fidelity and computationally cheaper models as much as possible is favoured as long as the final result is not jeopardised. Numerical solutions such as Reduced Order Models [38, 39], surrogate-models [40, 41] and multi-fidelity [42] can alleviate the computation burden. Furthermore, several methods can be used to obtain surrogate models based on samples of previously calculated results, therefore avoiding the evaluation of the constraints and objectives functions directly from expensive high-fidelity models [43]. There are some multi-fidelity models that have been applied to estimate the flutter boundary. Examples of surrogate based aeroelastic instability parametric search [44] and aerodynamic damping estimation for flutter calculations [45, 46] can be found in the literature. Multi-fidelity methods were also used to estimate flutter boundary, gust response [47] and aerodynamic optimisation [48, 49], although there is very scarce literature in multi-fidelity models applied to MDO considering flutter [50].

This study presents a quantification of the differences between FSI results using low- and high-fidelity aerodynamics models as a precursor to building a surrogate model and running a surrogate-based MDO of a HARW. It intends to provide researchers with practical information on the advantages and accuracy of multi-fidelity analyses approaches in aero-structural optimisation as well as on the design driving constraints, including flutter margin. It also reports the findings on the performance improvements of a surrogate-based optimised HARW configuration when compared to a baseline configuration for a medium-range-transport aircraft.

The paper starts with the problem description in the second chapter and following the methodology of the Multidisciplinary Design Analysis (MDA) and the surrogate modelling process in Sections 3 and 4, respectively. The investigation of the requirement for an FSI procedure for the evaluation of constraints and objective functions are concluded in Section 5. Furthermore, choices are made about which evaluations require high-fidelity and low-fidelity results for sufficient accuracy, leading to a mixed fidelity approach. Finally, the optimised results using exclusively low-fidelity or high-fidelity results are compared to the mixed-fidelity approach and between them, followed by final conclusions in Section 6.

2.0 Problem description

A baseline aircraft configuration has been provided by Bombardier Aerospace in order to develop the study of the HARW configuration. It represents a generic medium-range transport aircraft with a low wing with aspect ratio of 12, a typical fuselage and a T-tail. Both structural, fuel and system masses distribution were also provided, and a predefined jet engine database used to determine the fuel consumption given a thrust requirement and a flight condition. Figure 1 shows the mission profile. The cruise Mach number and altitude are 0.7 and 11km (36,000ft) respectively; the dive Mach number is 0.82. Throughout this study, the wing’s areofoils shape and spanwise relative position are assumed constant along with wing sweep. Hence, the optimisation result is limited in its potential, as overall changes in geometry would require a further aerodynamic optimisation. However, the computational costs would be significantly higher as the number of design variables would increase considerably. Furthermore, the benefit throughout the optimisation and the methodology of the multi-fidelity assessment can be shown with the given restraints in the detail of the optimisation. Fuselage, vertical tail geometry and horizontal tail planform are also assumed to be constant, as well as their structural components and mass distribution. The wing structure is assumed to be a wingbox-type structure made of isotropic aluminum. Variations allowed from the baseline wing configuration include parameters that affect directly both aerodynamic and structural characteristics (span, chord and twist distributions) and parameters that affect primarily the structural characteristics and indirectly the aerodynamic characteristics due to effects on deformation (spar and skin thickness distributions). As a consequence of the wing geometric changes, the horizontal tail area is changed in order to maintain the same Volume Coefficient as the baseline configuration.

Figure 1. Mission profile definition.

An MDO procedure is used to optimise the aircraft configuration taking into account environmental and operational costs metric (fuel consumption in cruise) and also a production cost metric (wing structural mass). The constraints for the optimisation include maximum allowable wing structural stress ( $\sigma$ ) in 3 different load cases and minimum flutter margin constraint. The flutter margin is on the flutter speed calculation in a 2.5g loading condition at dive speed defined for the aircraft mission.

The problem statement is as follows:

(1) \begin{equation} \begin{aligned} & \text{minimise} & & f(\boldsymbol{x}) = Aw(\boldsymbol{x}) + Bc(\boldsymbol{x})\\ & \text{w.r.t.} & & \boldsymbol{x} \\ & \text{subject to} & & \sigma_i(\boldsymbol{x}) \leq \frac{\sigma_{\rm yield}}{SF}; i =1..3\\ & & & FM(\boldsymbol{x}) \leq abs\left(\frac{FS(\boldsymbol{x})-DS}{DS}\right) \end{aligned}\end{equation}

where w and c are the objective function values of structural mass and fuel consumption, respectively. A and B are weighting constants for the multi-objective optimisation and were set by external reference according to an industrial weighting of $A/B=1/3$ , based on correspondence with the industrial partner and their earlier work [51, 52]. The load case index is described by i, the maximum allowed stress is yield stress of aluminum, $\sigma_{\rm yield}$ , paired with a safety factor, SF, of 1.5. The flutter margin is given as the difference of the estimated flutter speed and the dive speed, shown in the equation as FM, FS and DS, respectively.

The 14 design variables are listed in Table 1 and depicted in Fig. 2. The skin and spar thicknesses as well as the twist spanwise distributions are bilinear, i.e. are described by continuous functions composed of two linear segments with different slopes. The spanwise position where the slope changes is hereafter referred to as kink position. A decoupling of structural and aerodynamic kink position allowed more freedom for the optimisation. This gave the optimiser more freedom while the number of design variables could still be relatively low, considering the related high computational expenses for single function evaluations. In the next section, a description of the models and method used to evaluate the different figures of interest, i.e. drag, fuel consumption, structural stress, weight and flutter speed is given.

Table 1. Design variables for the wing parametrisation and optimisation

Figure 2. Representation of exemplary design variables.

3.0 Multidisciplinary analysis models

As stated before, the expectation is that increasing AR translates into higher deformations with effects on fuel consumption and flutter speed. Therefore, an FSI calculation coupling an Equivalent Beam Model (EBM) to a Panel Method (PM) is used as the low-fidelity model, using the code developed in Suleman et al. [53] and Afonso et al. [54]. The same EBM coupled with a RANS aerodynamic solver is used as the high-fidelity model.

The EBM model can provide estimates of the wingbox mass, as well as the distribution of the maximum absolute value of Von-Mises stress in the beam sections along the span. The model extracts the sectional properties of the wingbox sections and translates those into an Euler-Bernoulli beam model, as described in Afonso et al. [54]. The model is then used for linear structural analysis once the loading is defined.

The PM model includes fuselage, wing and tail components and calculates the lift and induced drag of the configuration. More detailed descriptions of these models regarding the formulation can be found in Katz and Plotkin [55] and regarding the validation in Suleman et al. [53] and Afonso et al. [54]. A friction drag correction is added to the calculations based on flat plate boundary layer models. Although fairly accurate to predict lift in the linear range of angles of attack (AOA), the model is not suitable to accurately predict transonic drag, shocks and shock-boundary layer interactions or separation.

Commercial software is used to perform RANS simulations of the aircraft configurations. This software is capable of more accurate drag predictions and accounts for the aerodynamic non-linearities usually found in the transonic regime, such as shocks, shock interactions with the boundary layer and separation [56, 57]. The unstructured polyhedral mesh is generated in an automated process. Refinements around the leading and trailing edges of the wing and tail ensure an adequate feature resolution. A mesh convergence study was performed on the baseline model, with the final model containing between 6 to 9 m cells, depending on the geometry, mainly span and chord. Menter’s Shear Stress Transport $k-\omega$ model [58] is used for the turbulence modelling. Convergence was assumed if a $\Delta C_L$ and $\Delta C_D$ criterium, 0.0001 and 0.00001 respectively, was fulfilled for at least 20 iterations in a row. The wall clock time for the 1g load case analysis was 5 hours and for the 2.5g load case 7 hours, where roughly 30% of the time was necessary for the meshing procedure for the deformed aircraft. The computation was executed on a server with 60 cores. ¹ An initial mesh convergence study was performed to establish the refinements and mesh sizing, exemplary values are listed in Table 2.

Table 2. Mesh properties

Other models are used to obtain the fuel distribution inside the wingbox, along with the payload and systems mass distribution related to the aircraft geometric parameters. Fuel is divided into small elements and modelled as point masses rigidly connected to the closest structural node. Fuel sloshing is not specifically considered as a steady cruise flight condition is assumed. Payload and system masses and their inertial components are user inputs and modelled as mass points with inertia and are also rigidly connected to the closest structural node. Again, a more detailed description of such modelling can be found in Oden [59].

The static FSI procedure couples the aerodynamic and structural models using an iterative procedure where the structural deformation affects the mass distribution and external geometry of the aircraft, which is then meshed and analysed in the aerodynamic solver to obtain a new load and then recalculate deformation [53, 54]. The procedure is considered converged when there is no significant change in the deformation after consecutive iterations. During the procedure, the fuel, payload and systems mass distribution is also affected by the structural deformation since these are rigidly linked to the structure.

During the FSI analysis the structure is clamped at the location of the centre of gravity because the aerodynamic loads will in general not trim the aircraft. In fact, the incidence of the horizontal tail required to trim the aircraft is calculated a posteriori together with the thrust and aircraft AOA after all the aerodynamic load components (lift, drag and pitching moment) are calculated, as will be explained later. As a consequence, the effect of the horizontal tail incidence on the wing loads and deformation are neglected in these analyses. This seems reasonable, given that the tail is located aft of the wing and that the tail loading acts in the fuselage. Figure 3 depicts the models described above.

Figure 3. Model representations: EBM (a); PM (b); RANS CFD (c); Payload and Systems mass distribution (d) and Fuel mass (e).

Finally, the propulsion model provided by Bombardier Aerospace is composed of a database with entries for thrust, fuel flow, altitude and Mach number, which are interpolated to obtain an approximation of the fuel consumption for a given thrust requirement and flight condition. This model was not incorporated into the FSI procedure in the current study, meaning that the thrust force is not included in the structural deformation calculations. Similar to the horizontal tail trimming effect, the thrust effect is deemed to be of negligible consequence in the wing’s deformed shape since the engines are located in the aft part of the aircraft and loading the fuselage, and the aircraft is assumed trimmed at a fixed fuselage pitch angle.

3.1 Trimming procedure

In order to trim the aircraft, a model of the horizontal tail with several different areas was analysed for different AOA using RANS. These generated results were then interpolated to obtain the lift, drag and pitching moment variation with AOA for the horizontal tail corresponding to wing area of the configuration being analysed. All configurations maintained the volume coefficient as constant to maintain similar control/stability characteristics.

The horizontal tail model, the propulsion model and the results obtained with the FSI procedure for a number of AOAs are then interpolated or extrapolated to obtain the trimmed AOA, horizontal tail incidence and thrust/fuel flow for a specific flight condition (speed, altitude, aircraft mass, load factor).

3.2 Flutter margin calculation

In order to calculate the flutter speed for a specific static flight condition, the FSI procedure described above is used to determine the (steady) deformed wing shape in that flight condition as predicted by the EBM model. The deformed EBM wing shape is then translated into NASTRAN for flutter speed calculation using the $p-k$ method [60] without introducing any CFD-based transonic corrections. As the translated structural model does not include the nominal steady-state stress in the structure, only the effects of the geometry are being accounted for in the analysis and no stress stiffening effects are included.

This constitutes the end of the MDA of a specific configuration.

4.0 Surrogate model based mdo

Given the workflow described in the previous section, the number of the more expensive FSI runs required per configuration should be minimised. The evaluation of the objective function and each constraint in the optimisation statement requires the FSI calculation at 2 different AOA for posterior trimming and flutter margin calculations. Gradient based optimisation algorithms require multiple evaluations per variable when the gradient functions are not directly available from solution results as in case of using the adjoint method [61, 62], in order to approximate the gradients using finite differences or other methods [63]. Non-gradient based optimisation algorithms, e.g. genetic algorithms, also require a large number of evaluations [64].

It is generally accepted that optimisation based on high-fidelity CFD results is costly in terms of computational resources and/or time. The FSI procedure alone using high-fidelity CFD is over 10 fold more time consuming than the proposed alternative, as it currently requires remeshing. This factor could potentially be reduced if a mesh deformation algorithm could be implemented but would still be about 5. In a preliminary design stage a lower level of fidelity is deemed sufficient to obtain comprehensive overall results. Therefore, taking into consideration this specific aircraft optimisation problem, strategies to minimise the number of FSI runs and also the high-fidelity CFD runs based on design knowledge are in order:

1. 1. Reduce the number of constraints assuming that it is known a priori which ones are active during the optimisation procedure. This would reduce significantly the number of evaluations required and the satisfaction of the remaining initial constraint can be verified after the optimised result is obtained.

2. 2. Substitute FSI results based on high-fidelity CFD by results based on low-fidelity CFD if and whenever the accuracy of results allows it. This could be the case if lift predictions by the low-fidelity CFD models are accurate enough and allow for satisfactory evaluations of the stress constraints.

3. 3. Perform FSI based on low-fidelity CFD and use the converged result to calculate the high-fidelity CFD results for the final deformed shape obtained. If the final shape is not significantly changed by the loading calculated with high-fidelity CFD the result might still be below the convergence criteria. This is more likely to occur for low-deformation flight conditions as in cruise and would allow assessment of the drag accurately, resorting to only one high-fidelity CFD run. If the accuracy allows, the same strategy can be used to evaluate and obtain the transonic correction for the flutter speed calculation at the 2.5g dive speed flight condition.

The previous possibilities to reduce the computational time aim to either reduce unnecessary constraints or reduce the number of high-fidelity CFD analyses. Another widely used strategy to reduce computational cost is to use surrogate models to evaluate the objective function and constraints with a dramatic reduction in the computational costs during the optimisation. The drawback is the generation of the results sample, which may require a significant number of evaluations to provide a sufficiently accurate model.

The next subsections explore the previously stated strategies for computational cost reduction in the MDO procedure.

4.1 Reduction of FSI and high-fidelity CFD analyses

In order to understand the validity of the strategies described above to limit the high-fidelity CFD runs in the MDA procedure an initial set of configurations was analysed. The set consists of 5 configurations with increasing AR that include the baseline configuration. All configurations have the same wing area while the span is varied between 95% and 115% and the chord distribution is scaled accordingly. The AR of the configurations therefore varies between 90% and 132% of the baseline configuration. The relative spanwise twist distribution is maintained for all configurations in this set.

While the geometric variables become defined by the statements above, the structural variables still require definition. A structural optimisation procedure was implemented in order to minimise the structural weight for each configuration. This procedure uses the same structural stress constraints as our MDO procedure, but no FSI is performed; therefore, the loading for each constraint is constant.

4.1.1 Assessment of active/inactive constraints

Figure 4 depicts the Von Mises stress distribution on the wing structural models for each structurally optimised configuration and for each constraint. The maximum stress is located around 40% of the wing, at the aerodynamic kink (Yehudi break). It is clear from Fig. 4 that the active constraint corresponds to the full fuel @ cruise speed 2.5g pull-up load case for all configurations. Both the full fuel @ cruise speed -1g load case and the half fuel @ cruise speed 2.5g pull up still have significant margin before becoming active constraints. Based on these results, the subsequent MDO procedure dropped the above-mentioned inactive constraints.

Figure 4. Von Mises stress distribution along the wingspan of the analysed set of configurations for each structural optimisation constraint (top to bottom: Load case 1, load case 2 and load case 3).

4.1.2 Deformation and loading: Low-fidelity vs high-fidelity

Once the structural variables were determined, FSI using low-fidelity CFD alone and low-fidelity followed by one FSI iteration using high-fidelity CFD were performed for the cruise flight condition and the active constraint flight condition. This is intended to quantify the differences between the loading and deformation predicted with low-fidelity CFD as compared to the high-fidelity CFD results, and determine in which cases the high-fidelity CFD requires further FSI iterations.

At the same AOA, the wing tip vertical displacement and the wing twist displacement show extrapolated relative differences between the low-fidelity and high-fidelity assessments up to 9% with respect to the undeformed wing for the cruise flight condition. These percentage differences were actually the maximum percentage differences found during the optimisation procedure, among the 108 configurations analysed. The vast majority of the configurations showed differences below 3%, with the average difference for the tip displacement 2.8% and for the twist displacement 1.8%.

These modest deflection changes are despite the differences in total lift for the cruise condition being higher in magnitude (up to 14.8% with an average of 5.5%). In this flight condition, the low-fidelity FSI result over predicts lift for the same AOA relative to the first iteration using high-fidelity CFD, which means that the difference in the load and its distribution is not enough to cause more significant differences in the wing’s deformed shape. Therefore, it is expected that an FSI scheme such as the one described above, where the low-fidelity FSI is followed by only one high-fidelity FSI iteration, would suffice to obtain accurate drag predictions since the shape of the deformed wing is nearly maintained when changing from low-fidelity to high-fidelity. Table 3 shows the results for the differences in vertical and twist displacements and lift coefficient for the analysed AOAs and for the calculated trim AOA.

Table 3. Cruise flight tip differences for vertical and twist displacements at the analysed AOAs and at the high-fidelity CFD based extrapolated trim AOA for each configuration

The shown values describe the differences between the results of the runs, performed with low-fidelity FSI and low-fidelity FSI with a single high-fidelity CFD evaluation. The column 1 and 2 refer to the two calculated angles of attack, while the column Tr is for the extrapolated trimmed case for the cruise flight condition. Differences below 10% for the low-fidelity FSI are deemed acceptable for the preliminary aircraft design study.

For the active constraint case the situation is reversed, with the low-fidelity FSI under predicting the load for the same AOA, which results in a predicted higher AOA to obtain the same 2.5g lift load. Reasons for this difference can be the more intense transonic non-linearities, for instance shocks and flow separation onset. Consequently, a higher AOA is interpolated to retrieve the same load. Table 4 shows the calculated AOAs and the differences in vertical and twist displacement at the same load based on the low-fidelity FSI results and on the first high-fidelity FSI iteration for the 2.5g lift load.

Table 4. Active constraint tip differences for vertical and twist displacements between at the analysed AOAs and at the high-fidelity CFD based interpolated trim AOA for each configuration

Figures 5 and 6 depict the vertical and twist angle displacement spanwise distributions results for configurations AR10.9, AR13.2 and AR15.7 for both the cruise and active constraint trimmed flight conditions. These results show that the differences between the low-fidelity FSI and the low-fidelity FSI with single high-fidelity CFD tend to increase as the AR increases and as the load magnitude increases.

Figure 5. Vertical displacement vs spanwise station for configurations with AR 10.9 (Top), AR 13.2 (Middle) and AR 15.7 (Bottom) for a cruise flight condition (left) and a 2.5g pull up @ MTOW flight condition (right).

Figure 6. Twist displacement vs spanwise station for AR 10.9 (Top), AR 13.2 (Middle) and AR 15.7 (Bottom) for a cruise flight condition (left) and a 2.5g pull up @ MTOW flight condition (right).

4.1.3 Stress assessment: Low-fidelity FSI plus one high-fidelity FSI iteration vs full high-fidelity FSI

As a result of this higher magnitude difference in the high load case, the stress distribution shows differences which are significant in the context of an MDO procedure and feasibility of a design cannot be assessed based on low-fidelity FSI alone, since it is not a conservative estimate (Fig. 7). Given the higher magnitude of the negative twist for the first high-fidelity FSI iteration, it is expected that some load alleviation would cause a redistribution of lift and a relief of the bending moment in the wing for the same 2.5g load, resulting in a stress distribution closer to the low-fidelity FSI result. The one high-fidelity FSI iteration procedure would then be a conservative way to avoid the unacceptable computational cost of a full high-fidelity FSI run to assess the stress constraint. A more dominant difference in the stress results can be seen in the distribution when comparing the results with and without one high-fidelity FSI iteration. Due to the altering geometric shape and the emerging change in lift distribution, especially for higher AR, the stress distribution distinguish significantly from the rigid results (as shown in Fig. 4).

Figure 7. Von-Mises equivalent stress spanwise distribution for AR 10.9 (Top), AR 13.2 (Middle) and AR 15.7 (Bottom) for the 2.5g pull up @ MTOW flight condition.

A full high-fidelity CFD based FSI run was performed for the best configuration (AR 14) obtained during the optimisation procedure. This enabled validation of the assumptions of reduced influence of load differences between the low-fidelity FSI followed by one high-fidelity FSI iteration and the converged results of a complete high-fidelity FSI run in the lift and drag results. It also allowed the verification that the procedure of performing one high-fidelity FSI run after a converged low-fidelity FSI produces conservative stress estimates when compared to a converged full high-fidelity FSI run. Figure 8 depicts the results.

Figure 8. High-fidelity CFD-based FSI results for the optimised configuration (AR 14) compared to low-fidelity FSI and to one high-fidelity FSI run after low-fidelity FSI convergence for a constant AOA=3º and different dynamic pressures.

One can observe that both assumptions seem to be valid for the configuration analysed, with the converged high-fidelity FSI results being closer to the low-fidelity FSI results for dynamic pressures approximating the cruise and the active constraint cases, therefore reducing the shape differences as assumed for the cruise case and showing the stress results as being conservative in the active constraint case.

Table 5 shows the differences in lift and drag of the best configuration (AR 14) for the same AOA and the dynamic pressure approximating the cruise case. The results show that the magnitude of the error of the used procedure relative to the converged high-fidelity FSI is around 1% in drag and around 0.5% in both lift and L/D for this configuration. Under the assumption that these results can be generalised for the current problem, any configuration with improvements in L/D above 1% of the baseline model is deemed as having a better performance than the baseline.

Table 5. Comparison of lift and drag results between the used procedure (UP) and the converged (C) high-fidelity FSI for the optimised configuration (AR 14) for a constant $AOA=3^\circ$

4.1.4 Flutter margin assessment

Flutter boundary assessment for flutter margin calculations was performed to understand if this would be an active constraint. Figure 9 shows the flutter boundary results for the configurations with AR 10.9, AR 13.2, AR 15.7 and the optimised configuration in undeformed and deformed states, together with the flight envelope and the flutter margin boundary requirement.

Figure 9. Flight envelope and required flutter margin and flutter boundaries for configurations AR10.9, AR13.2, AR15.7 and the optimised one in undeformed and deformed states.

It is apparent that the flutter boundary is altered, with deformation reducing the flutter margin for all configurations shown. Nevertheless, the flutter speeds calculated greatly exceed the required flutter margin. Therefore, the flutter margin does not seem to be an active constraint in the optimisation, and it is confirmed that the best configuration respects the minimum flutter margin constraint.

4.1.5 Surrogate models

In the proposed optimisation procedure, a Bayesian approach also known as Kriging [65] was used as a surrogate model, which assumes that the error in the objective function is a realisation of a Gaussian process. It is a spatial analysis method, considering a global low order regression, paired with local deviations based on covariances between the samples to fit and interpolate the data. The model constructs a response at any point consisting of a mean value and stochastic deviation:

(2) \begin{equation} \hat{y}(\boldsymbol{x}) = \hat{\mu} + \boldsymbol{r}^T_{\boldsymbol{xX}} \boldsymbol{R}^{-1} \boldsymbol{r}_{\boldsymbol{xX}} (\boldsymbol{y} - \textbf{1}),\end{equation}

where $\mu$ is the mean value at $\boldsymbol{x}$ , $\boldsymbol{r}_{\boldsymbol{xX}}$ is the vector of covariances between the new point and the existing ones, $\boldsymbol{R}$ is the covariance matrix of sample points, $\boldsymbol{y}$ the response of the evaluated sample points and $\textbf{1}$ is a column vector of length of the number of initial designs. The mean squared error can be estimated in a Gaussian process by:

(3) \begin{equation} \hat{s}(\boldsymbol{x}) = \hat{\sigma}^2 (\boldsymbol{1} - \boldsymbol{r}_{\boldsymbol{xX}}^T \boldsymbol{R}^{-1} \boldsymbol{r}_{\boldsymbol{xX}}),\end{equation}

with $\sigma^2$ the standard deviation. The correlation function used is the Gaussian exponential function:

(4) \begin{equation} k(\boldsymbol{x},\boldsymbol{x}^*) = \exp \left(- \sum_{j=1}^{n} \theta_j |x_j^*-x_j|^2\right)\end{equation}

where $\theta$ is the hyperparameter to be estimated in the Kriging model building process. Each $k(x_j,x_j^*)$ is an entry in the correlation matrix R, for which with the maximum likelihood estimation (MLE) the hyperparameter are determined. More detailed information on Kriging can be found for instance in Bhosekar and Ierapetritou [66]. The current work is based on the DACE (Design and Analysis of Computer Experiments) implementation [67], which is based on the work of Sacks et al. [68], with an additional extension enabling the surrogate model to consider regression rather than interpolation as noise treatment. This is necessary as inaccuracies due to meshing error, numerical rounding effects and likewise could otherwise result in a poorly built model if close points deviate too much. This additional hyperparameter is also estimated by the DACE algorithm [69].

To find new samples with the purpose of achieving an optimised result, different search criteria can be employed. In this work an Expected Improvement (E[I]) based infill criterion is used, using the uncertainty quantification provided by the Kriging model [70]. With the distributed samples in the design space the criterion can require only a small number of iterations to achieve a satisfactory optimisation result.

(5) \begin{equation} E[I(\boldsymbol{x})] = \begin{cases} \left( y_{\min} - \hat{y}(\boldsymbol{x})\right)\!\Phi\!\left(\dfrac{y_{\min} - \hat{y}(\boldsymbol{x})}{\hat{s}(\boldsymbol{x})} \right) +s \phi\!\left(\dfrac{y_{\min} - \hat{y}(\boldsymbol{x})}{\hat{s}(\boldsymbol{x})} \right), & \rm{if} s > 0, \\ 0, & \rm{if} s = 0, \end{cases}\end{equation}

where $\Phi$ and $\phi$ are cumulative distribution and probability density functions, respectively and s is the mean squared error.

The initial sample is generated by Latin Hypercube Sampling (LHS) [71] to ensure an adequate design space sampling. When the simulator cannot evaluate a specific sample, which can happen for some design variable combination, meshing errors or likewise, the sample was taken out and replaced by a different sample in a similar region (short Euclidian distance of not more than 5% of the normalised design space). After the initial sampling, the surrogate based optimisation with iterative sample addition was started and run until a maximum number of iterations was reached. In each iteration the new sample is evaluated for the objective and constraint function and added to the initial database. This improves the model locally and, as each sample has a correlated influence on the other samples, the global prediction as well. In this procedure the prediction of surrogate model becomes less erroneous and the predictions of the surrogate model therefore closer to the real response.

5.0 MDO Results

The MDO process followed the scheme pictured in Fig. 10. Starting with the initial sample and corresponding evaluated results, surrogate models were built for the objective (weighted function of wing structural mass and fuel consumption) and constraint function (stress distribution along the wing structure). The new sample, determined by the infill criterion, was then evaluated. If a convergence criterion or the maximum number of iterations is reached, the process is terminated.

Figure 10. MDO scheme.

In the present optimisation the process was terminated after 110 iterations reaching the maximum available number of iterations and no remarkable improvement could be achieved in the later iterations.

Figure 11 shows the fuel savings relative to the baseline throughout the optimisation process, where black square markers indicate feasible configurations, white circle markers indicate unfeasible points. Gray triangle markers indicate configurations for which the stress constraint was not assessed since they did not show fuel savings improvement relative to the best configuration at that point of the optimisation.

Figure 11. Optimisation progress.

The attained results are seen as satisfying, considering the available resources, rather than a rigorously defensible globally optimal result. Table 6 below shows the relative changes of the significant parameters of the best optimised versus baseline configuration, as well as the maximum occurring stress.

Table 6. Optimal configuration differences relative to baseline configuration

Figure 12 shows a comparison of the outer shapes clearly showing the noticeable AR increase above the baseline. It is expected that with further continuation of the optimisation additional improvement could be achieved, based on the progress of the optimisation specifically in the end phase, where the average improvements increased significantly compared to the baseline configuration.

Figure 12. Overlapping comparison of the baseline (dark gray) and the optimised wing (light gray) configuration planforms.

6.0 Concluding remarks

In order to reduce the computational costs of MDO of a typical medium-range commercial aircraft including flutter constraints, studies were performed on the most likely active constraints in the MDO problem and on the validity and accuracy of estimating drag and stress using low-fidelity FSI followed by one iteration high-fidelity FSI versus a converged high-fidelity FSI calculation. The following remarks summarise the findings:

• Evaluation of the constraints can be minimised through a proper assessment of which of them are likely to be active. This includes the flutter margin constraint which, according to the results is also not likely to be active within the feasibility region for a HARW design. Although this seems an obvious step for an efficient optimisation, without prior knowledge of which constraints are active, it could not be ensured a priori. This work provides a basis for future researchers to non-arbitrarily disregard the flutter constraint during the optimisation of slenderer wings and perform a verification later on.

• 110 configurations where analysed using low- and high-fidelity CFD showing that differences in load distribution and stress for the same load magnitude do not produce significant differences in the displacement and stress levels, rendering the use of high-fidelity CFD for a high load case unnecessary, even for what can be considered high aspect ratio wings in commercial aircraft.

• With the high load case assessment, this work also indicates that by using low-fidelity CFD followed by one high-fidelity analysis in a FSI procedure the results are conservative. This knowledge helps future research to setup more efficient optimisation problems.

• With the cruise load case assessment, this work shows that FSI with full high-fidelity CFD is hardly justifiable if the multi-fidelity approach is available. Note that this is relative to the multidisciplinary analysis and not related to the surrogate based MDO (i.e. there is not a low-fidelity based surrogate model. The low-fidelity model is used as a speeding up process to obtain higher fidelity results in the multidisciplinary analysis).

• The methodology was applied to a representative configuration of a medium-range transport aircraft, showing that performance improvements can be expected by increasing the Aspect Ratio relative to current configurations even with structural mass increase, considering the stress constraint and also the flutter margin, although the flutter margin is degraded as the wing flexibility increases.

The results show that for this problem the active constraint is only one of the load cases, which is the 2.5g pull up at MTOW. The other load cases considered appeared less severe and therefore were left unassessed during the MDO procedure. The flutter constraint is also not active, although the effect of deformation on the flutter margin reduction is visible in the results. As novel aircraft configurations tend towards higher aspect ratios with more flexible wings, it is expected that dynamic aeroelastic behaviour will become a more prominent design driver. Slenderer wings such as of HARW or SBW are more prone to high deformations, closing into the flutter margin. Altering the number of (structural) design variables is expected to bring the stress distribution closer to the boundary along the whole wing (outer wing) and increase the authority of the flutter margin as a constraint driver.

Drag evaluation in the cruise case can be performed using the low-fidelity FSI followed by one high-fidelity FSI. The low-fidelity FSI allows speed-up of the deformed shape calculation to an acceptable level of accuracy when compared to a converged high-fidelity FSI run. Maximum displacement differences between low-fidelity FSI and after one high-fidelity FSI iteration are generally below 3%, and tend to be reduced for the converged high-fidelity FSI results. Differences in drag and L/D between the first high-fidelity FSI iteration and the converged high-fidelity FSI were estimated to be around $-1$ % and 0.5% for the best configuration in the analysed set.

Stress evaluation can also be performed using the same procedure, but this time not with the same level of accuracy although still conservative. Converged high-fidelity FSI results tend to redistribute the loads for this wing geometry and reduce the stress to values between the low-fidelity FSI result and the first high-fidelity FSI iteration results. Unfortunately, the low-fidelity FSI load predictions alone are not accurate or conservative enough to dismiss the high-fidelity evaluation.

A further reduction in computational costs was attained using a surrogate based MDO approach. Using surrogate models for both objective function and active constraint, the E[I] was used to determine the surrogate models infill configuration designs and the best feasible configuration after 110 iterations was chosen as the optimisation result. For the different used models and to evaluating parameters, the chosen approach is a viable alternative compared to more complex approaches like adjoint class methods. The outcome of the optimisation and the consideration of different levels of fidelity and complexity was manageable in a straightforward way. After all, the chosen configuration had an increase in AR of 16% and yield a reduction in fuel consumption of 4.5% while increasing the wing weight by 2.7%.