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

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 [Reference Piperni, Abdo, Kafyeke and Isikveren51, Reference Priem, Gagnon, Chittick, Dufresne, Diouane and Bartoli52]. 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. [Reference Suleman, Afonso, Vale, Oliveira and Lau53] and Afonso et al. [Reference Afonso, Vale, Oliveira, Lau and Suleman54]. 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. [Reference Afonso, Vale, Oliveira, Lau and Suleman54]. 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 [Reference Katz and Plotkin55] and regarding the validation in Suleman et al. [Reference Suleman, Afonso, Vale, Oliveira and Lau53] and Afonso et al. [Reference Afonso, Vale, Oliveira, Lau and Suleman54]. 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 [Reference Blazek56, Reference Rodriguez57]. 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 [Reference Menter58] 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.Footnote ¹ 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 [Reference Oden59].

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 [Reference Suleman, Afonso, Vale, Oliveira and Lau53, Reference Afonso, Vale, Oliveira, Lau and Suleman54]. 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.

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 [Reference Martins and Hwang61, Reference Gray, Hwang, Martins, Moore and Naylor62], in order to approximate the gradients using finite differences or other methods [Reference Snyman and Wilke63]. Non-gradient based optimisation algorithms, e.g. genetic algorithms, also require a large number of evaluations [Reference Schaefer64].

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 [Reference Matheron65] 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 [Reference Bhosekar and Ierapetritou66]. The current work is based on the DACE (Design and Analysis of Computer Experiments) implementation [Reference Lophaven, Nielsen and Søndergaard67], which is based on the work of Sacks et al. [Reference Sacks, Welch, Mitchell and Wynn68], 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 [Reference Forrester, Sbester and Keane69].

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 [Reference Jones and Schonlau70]. 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) [Reference McKay, Beckman and Conover71] 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%.