Static and dynamic aeroelastic tailoring with composite blending and manoeuvre load alleviation

In aircraft design, proper tailoring of composite anisotropic characteristics allows to achieve weight saving while maintaining good aeroelastic performance. To further improve the design, dynamic loads and manufacturing constraints should be integrated in the design process. The objective of this paper is to evaluate how the introduction of continuous blending constraints affects the optimum design and the retrieval of the final stacking sequence for a regional aircraft wing. The effect of the blending constraints on the optimum design (1) focuses on static and dynamic loading conditions and identifies the ones driving the optimization and (2) explores the potential weight saving due to the implementation of a manoeuvre load alleviation (MLA) strategy. Results show that while dynamic gust loads can be critical for wing design, in the case of a regional aircraft, their influence is minimal. Nevertheless, MLA strategies can reduce the impact of static loads on the final design in favour of gust loads, underlining the importance of considering such load-cases in the optimisation. In both cases, blending does not strongly affect the load criticality and retrieve a slightly heavier design. Finally, blending constraints confirmed their significant influence on the final discrete design and their capability to produce more manufacturable structures.

or maximise a specific behaviour (Fukunaga et al. 1995;Diaconu and Sekine 2004;Abdalla et al. 2007;IJsselmuiden et al. 2010) and recent advancements in manufacturing technique allow for fibre steering to further enlarge the design space . Finally, while for large composite structures local laminate optimisation can produce a lighter and more efficient design, the optimum design can be characterised by a significant thickness and/or stacking sequence variations between adjacent laminates. In this case, the obtained design could be too expensive to manufacture and might also lack structural integrity (Dillinger 2014;IJsselmuiden et al. 2009). Therefore, ply continuity among adjacent plies (i.e. blending) should be considered early in the design phase.
One common strategy used to optimise composite structures relies on a bi-step approach (Herencia et al. 2007;IJsselmuiden et al. 2009;Liu et al. 2011;Dillinger 2014), where a gradient-based (continuous) optimisation of homogenised stiffness parameters (e.g. lamination parameters) is performed to obtain the optimum design. Later, a combinatorial optimisation is used to retrieve the best stacking sequences matching the continuous design while ensuring blending (Soremekun et al. 2002;Adams et al. 2004;van Campen and Gürdal 2009;Irisarri et al. 2014). This second step (discrete optimisation) is typically performed via evolutionary algorithms (e.g. genetic algorithms). However, due to the difficulties in enforcing ply continuity in the gradient-based optimisation with homogenised stiffness parameters, blending has always been implemented during stacking sequence retrieval. As a consequence, two different sets of constraints are used in the two subsequent optimisation. This reduces the chance to find in the discrete domain a design close to the optimal continuous solution. Recently, Macquart et al. (2016a) proposed employing lamination parameters combined with a set of blending constraints to be used in the continuous optimisation to achieve more realistic and manufacturable continuous designs. In Macquart et al. (2016a), the continuous blending constraints have been applied to the benchmark case of the 18-panel horseshoe to prove the effectiveness of the blending method. In a continuing effort, Macquart et al. (2016b) and Bordogna et al. (2016) have demonstrated that the application of blending constraints during aeroelastic optimisations with strain and buckling constraints results in more realistic continuous designs.
Conventional wing structural sizing of large transport aircraft usually considers symmetric static manoeuvres, like 2.5g pull-up and −1g push-down manoeuvres, as design loads together with aeroelastic phenomena like divergence, flutter and aileron effectiveness (Torenbeek 2013). However, Kenway et al. (2014) showed that a metallic large transport aircraft wing optimised for static manoeuvres loads can fail when subjected to discrete gusts, underlying the need to include dynamic load cases during optimisation. Werter (2017) obtained similar results with a composite large transport aircraft wing and showed how wing optimised with static loads and unbalanced stacking sequence are more prone to failure under dynamic loads than wing designed with a more conventional stacking sequence (e.g. [0 60% / ± 45 30% /90 10% ] s ). Finally, dynamic loads are also influenced by the flight dynamics of the aircraft as shown by Reimer et al. (2015).
In a preliminary work of Bordogna et al. (2017), the authors proposed a strategy to optimise a composite regional aircraft wing for static and dynamic aeroelastic loads, blending constraints and manoeuvre load alleviation (MLA). The proposed strategy, together with the work from other researchers, has been then adopted by DLR and integrated into the in-house tool MONA (Bramsiepe et al. 2018) with the purpose of performing a comprehensive load analysis and designing a benchmark aeroelastic model of the Airbus XRF1 (Vassberg et al. 2008) for later studies.
This paper offers a more exhaustive analysis and followup of the activities presented in Bordogna et al. (2017). The authors focus on the effect of blending constraints on the identification of the critical loads during aeroelastic tailoring of a regional aircraft wing subjected to both static and dynamic load-cases. Moreover, the effect of such constraints on manoeuvre load alleviation on the critical loads is also assessed. Finally, the influence of blending constraints over the optimal design is presented together with the "ready-to-manufacture" quality of the retrieved stacking sequence.
The paper is divided as follows. In Section 2, the concept of blending is introduced together with the chosen composite parametrisation method. Section 3 presents the wing model used in this work and the loads considered. Then, in Section 4, the optimisation problem and strategy are explained together with the concept of equivalent static load (ESL). Finally, results and conclusions are presented in Sections 5 and 6 respectively.

Lamination parameters space and composite blending constraints
Large composite structures can be divided into sections that are subsequently locally optimised to obtain lighter and better performing structures. However, this local optimisation can lead to significant discrepancies in thickness and stacking sequence among adjacent sections resulting in an optimal solution that lacks structural integrity. To ensure a certain degree of ply continuity, the definition of blending has been first introduced by Kristinsdottir et al. (2001).
Another challenge in dealing with locally optimised composite structures is the large number of design variables proportional to the number of sections and the number of plies in each section (Bettebghor 2011). To reduce the number of design variables to a constant value regardless of the stacking sequences thicknesses, homogenised stiffness parameters (i.e. lamination parameters) are used. In this section, the lamination parameters used for composite parametrisation are introduced in Section 2.1, while different definitions of blending are given in Section 2.3 and a brief introduction to the blending constraints used in this work is given in Section 2.4.

Lamination parameters
Lamination parameters (LPs) have been first introduced by Tsai and Hahn (1980) and are used to parameterise the stiffness matrix of composite laminates in a continuous space. For stacking sequence with discrete plies of constant thickness (t ply ) and ply angle (θ i ), lamination parameters are defined as in (1). In this paper, only symmetric stacking sequences with an even number of plies and constant ply thickness are considered. Therefore, only lamination parameters for membrane (A) and bending (D) stiffness matrices are taken into account.

Membrane stiffness visualisation
Lamination parameters have the advantages of describing the stiffness matrix in a continuous form and they define a convex space (Grenestedt and Gudmundson 1993) suitable for gradient-based optimisation. Moreover, mechanical quantities often have a simple dependence on lamination parameters; for example, buckling load factors is a concave function of lamination parameters (Bettebghor and Bartoli 2012). Furthermore, any symmetric stacking sequence can be reproduced with eight continuous variables plus laminate thickness. On the other hand, the use of LPs requires an additional optimisation step that retrieves a discrete stacking sequence from the continuous optimal design. This extra step is usually performed by evolutionary algorithms. Therefore, a two-step optimisation strategy is required (see Section 4.1.7).
While lamination parameters offer many advantages, it is not straight forward to reconstruct the main stiffness direction associated to a set of parameters. As introduced by Dillinger et al. (2013), to have a sense of the main inplane stiffness distribution of a given A matrix, it is possible to calculate the its normalised elastic modulus of elasticity (Ê 11 (θ )) associated to the component A 11 along an axis rotated with an angle θ with respect to the axis of the laminate as: where T is a transformation matrix: T = ⎡ ⎣ cos 2 θ sin 2 θ 2 cos θ sin θ sin 2 θ cos 2 θ −2 cos θ sin θ − cos θ sin θ cos θ sin θ cos 2 θ − sin 2 θ ⎤ ⎦ (8) Few characteristic examples of laminates and their corresponding membrane stiffness distributions are shown in Fig. 1, where the x axis represents the axis of the laminate. The stiffness distributions presented here have all been normalised with respect to the maximum stiffness associated to the single unidirectional ply.

Blending definitions
Different definitions for blending have been proposed by different authors (Fig. 2). The purpose of these blending approaches is to guarantee a certain degree of ply continuity inside a variable stiffness composite structure and therefore increasing its structural integrity.
Inner blending and outer blending have been introduced by Adams et al. (2004); in these definitions, only the innermost and the outermost plies can be dropped. These two definitions are the least flexible because plies can be drop at only one location, reducing significantly the stacking sequence design space. On the other hand, the reduced design space allows for a quicker stacking sequence retrieval but at the expense of a possible weight penalty.
Two alternative definitions, the generalised blending and relaxed generalised blending, have been formulated by van Campen et al. (2008). Generalised blending requires all plies of the thinnest section to be continuous in the whole structure. This definition allows a ply to be dropped at any location inside the stacking sequence resulting in a wider design space than Adams's definitions. Relaxed generalised blending has an even wider design space by only demanding that no discontinuous plies should be in direct physical contact with each other. By doing so, two adjacent panels can have the same thickness but different stacking sequences, and this is not possible with any of the other definitions. For the sake of clarity, throughout this paper, blending is always associated with the generalised blending definition of van Campen et al. (2008).

Implementation of blending constraints
Several authors took advantage of the continuous description of composite materials offered by the lamination parameters and relied on bi-step strategies to retrieve a discrete stacking sequence. However, due to the difficulties in enforcing ply continuity during the gradient-based optimisation in the lamination parameter space, blending has always been enforced during stacking sequence retrieval. As a consequence, two different sets of constraints are used in the two subsequent optimisation steps. This reduces the chance to find an equivalent of the optimal continuous design in the discrete domain. To overcome this problem, in this work, ply continuity among adjacent sections is enforced by means of the continuous blending constraints introduced by Macquart et al. (2016a) in the lamination parameter space.
The key concept in the continuous blending constraints is to quantify the change in lamination parameters ( v) due to ply drops. A comprehensive derivation of all the blending constraints can be found in Macquart et al. (2016a). Here for the sake of completeness, the derivation of the blending constraints for a single in-plane lamination parameter (13) is presented.
Let us denote v A 1(N) and v A 1(N−X) the value of the first membrane lamination parameter when the laminate has respectively N and N − X plies. The change in lamination parameter due to any number X of ply drops denoted by the set S is presented in (11).
containing the dropped plies Term containing the plies present in both sections (11) where X is the number of dropped plies, N is the total number of plies, S represents the set of dropped plies. The maximum and minimum values of (11) occur respectively for [θ j , Equation (12) implies that no blended solution can be found if the change in v A 1 in two adjacent sections is greater than 2(X/N). By applying the same approach to the remaining membrane lamination parameters, it can be shown that this limit holds. Thus, it is possible to define a blending constraint for single membrane lamination parameter change as: For the sake of brevity, the extension of the blending constraints to the higher dimensions (i.e. taking into account more LPs at once) is not covered (the interested reader is invited to check the work of Macquart et al. (2016a)). However, to provide the reader with a visual representation of the constraints, an example of the effect of a 2D blending constraints considering v A 1 and v A 2 is hereby provided. Let us take into account the situation presented in Fig. 3, where a multi-section panel is subjected to X ply drops from a section with N plies to another with N − X plies. Let us now reproduce this situation in the lamination parameter space for v A 1 and v A 2 (Fig. 4). In this example, the starting laminate section has N equal to 20 plies and it is used to generate all possible N − X plies blended sections by removing X equal to 4 plies. The blending constraint for this ply drop configuration is shown in red and it is capable of including all possible N − X plies blended panels. As shown in Macquart et al. (2016a), and visible in Fig. 4, the blending constraints are likely to be over-conservative in most case because they have been derived considering the worst possible ply drop configuration without taking into account for the effect of manufacturing constraints in the stacking sequence. For this reason, a shrinking factor α can be added to the constraints to reduce the hypersphere and tighten the constraints, shown in blue in Fig. 4. For the case presented in this paper of a single in-plane parameter, the shrinking factor is included in the derived constraint as follows:

Model description
This section introduces the structural model used in this work together with the applied loads and the notions of equivalent static load and manoeuvre load alleviation.

Structural model
The structural model used to test the developed optimisation strategy is an ONERA's internal composite wing model of a regional aircraft. In the finite element model (FEM), the wing skins, shear ribs, and spars are represented with shell elements, while stringers in the wing skins are represented with beam elements. Shell elements in the wing skins are grouped in sections in the spanwise and chordwise direction, and each section shares the same thickness and material properties. The wing spars are also grouped in sections in the spanwise direction (Fig. 5). The wing internal structure is composed of 32 shear ribs and 13 stringers. Wing dimensions and characteristics are presented in Table 1. Two different fuel mass distributions have been considered in this work, the maximum take-off weight (MTOW) and maximum landing weight (MLW) (Fig. 6). The fuel is modelled as concentrated masses and connected to the wing box FEM via multi-point connections. The amount of fuel along the spanwise direction has been estimated by considering the volume of each rib-bays, defined as the space between two consecutive ribs and the two spars, and the jet fuel density. Once the fuel weight of each rib-bays is computed, its value is added to a concentrated mass positioned at the center of gravity of the rib-bay.

Doublet-lattice method and correction
The aeroelastic loads are calculated via the MSC.Nastran (2014) aeroelastic solver that utilises doubletlattice method (DLM). Since the aeroelastic loads coming from the DLM are used to perform the trim of the wing and calculate its displacement, it is essential to correctly represent the spanwise and chordwise load distribution along the wing. This is achieved by using the concept of separation between the rigid and elastic load components, where the rigid part utilise rigid CFD results while the elastic increment is computed via DLM. This method, often referred to as hybrid static approach (HSA) (Vincenzo 2012), allows to consider in the aeroelastic load computation for airfoil camber and wing twist law.

Static manoeuvre loads
Symmetric static manoeuvres at 2.5g and −1g are considered as load cases in this work at two different flight points (i.e. cruise and sea level) and for two mass configurations (i.e. MTOW and MLW). Static manoeuvres are usually considered as the sizing load-cases for wing where M is the structural mass matrix, K is the structural stiffness matrix, Q is the aerodynamic influence coefficient matrix, Q x matrix provides force at the structural grid points due to deflection of aerodynamic control surfaces, P is the vector of applied loads (e.g. thermal, gravity) andq is the dynamic pressure. The vectors u and u x are the wing structural deflection and the aerodynamic control surface deflection respectively.

Dynamic loads and equivalent static load
Dynamic aeroelastic analysis of discrete gust is performed to ensure the wing satisfies the CS-25 certification requirements (EASA 2018) for gust and turbulence loads. The dynamic aeroelastic analysis is performed in the frequency domain; thus, the time domain gust load is first transformed using Fourier transform in the frequency domain. A frequency analysis is then performed, and the responses are subsequently converted back to the time domain. The equations-of-motion in modal coordinates are: where M is the modal mass matrix, B is the modal damping matrix, K is the modal stiffness matrix, Q is the modal aerodynamic force matrix that is function of the Mach number m and the reduced frequency k, P (ω) is modal loading vector, ω is the angular frequency and u is the modal amplitude vector. The CS-25 requires the aircraft to be subjected to discrete symmetrical vertical and lateral gusts at level flight. In this work, only discrete symmetrical vertical gust is taken into account. The gust shape must be defined as in (17).
where U is the gust velocity in equivalent air speed (EAS) at position s, s is the distance travelled inside the gust, U ds is the design gust velocity in EAS (18), H is the gust gradient in meters, U ref is the reference gust velocity in EAS, F g the flight profile alleviation factor (19) and Z mo is the maximum operating altitude in meters. The ONERA's regional wing the has a MTOW of 60,000 kg, a MLW of 55,000 kg, a maximum zero-fuel weight (MZFW) of 50,000 kg and a maximum operating altitude Z mo of 12,192 m (40,000 ft). For gust computations, 10 different gust gradients H have been used from 9 and 107 m. Once the gust shapes (17) are defined and used in the dynamic aeroelastic analysis (16), the responses to the gusts are obtained in time the domain. To build up the total response to the gusts at flight level, a separate 1g cruise static aeroelastic analysis is performed and the result is superimposed to the gust responses. Positive and negative total gust responses are obtained by adding or subtracting the gust responses from the 1g cruise configuration.
In this work, as explained with more details in Section 4, a gradient-based optimiser is used as the number of design variables is relatively large. The main issue when accounting for gust loads is their direct dependency on the design itself. As a result, gust loads and their sensitivities need to be recomputed for each design cycle, for all gusts and at all instants. Such operations are computationally expensive, and therefore it is often difficult to include transient analysis in an optimisation process (Kang et al. 2006). Nonetheless, this process has been implemented in dedicated tools like the TU Delft Proteus aeroelastic code (Rajpal et al. 2019) and Airbus Lagrange tool (Petersson 2009). The equivalent static load method formalised by Kang et al. (2001) is used in this paper to bypass this issues.
The ESL aims at computing one or more equivalent static loads f eq capable of generating the same displacement fields of the transient load at different critical time steps. These f eq , now assumed as constants with respect to the structural design, are then applied to the FEM as static loads while the design is being optimised. Once the optimised design is obtained, a new set of transient analyses are performed to update f eq and the loop is repeated until convergence. Therefore, the ESL method relies on a  Fig. 9 and explained in Section 4.1.8. The lack of sensitivities between the design variables and the transient responses constitutes one of the main drawbacks of this method. Therefore, design changes between two consecutive ESL loop need to be small enough to ease constraints satisfaction and convergence. Otherwise, relaxation strategies could also be implemented. Still, this method offers an easy implementation regardless of the different tools used in the loop and can take advantage of the already-existing gradient-based optimisation and aeroelastic analysis code. The ESL has been used in various contexts such as non-linear structures, multi-body dynamic and crash and topology optimisation for the automotive industry. Most of these applications have been summarised in Park (2006). However, the affirmation that the primal-dual Karush-Kuhn-Tucker (KKT) point of the converged ESL solution is also a primal-dual KKT point of the original dynamic response problem is still an open question (the interested reader is invited to read the works from Park and Kang (2003) and Stolpe et al. (2018)).

Manoeuvre load alleviation
Manoeuvre load alleviation is a technique that allows aeroelastic load reduction during static manoeuvres by moving the spanwise lift distribution inboard. This is usually achieved by control surfaces deflections that shift inboard the lift distribution. Thus, MLA leads to smaller wing root bending moment that translates in lighter wing structure. While MLA is used for static manoeuvres, gust load alleviation (GLA) is used to control the dynamic loads during gust or turbulence encounter. In this work, the authors look only at MLA to reduce the static aeroelastic load during symmetric 2.5g and −1g flight manoeuvres by means of aileron deflection.
Several optimisations, each of them performed with different aileron deflection settings, are performed separately. The different optimal results are then compared to assess the potential weight saving of this technique. For the sake of clarity, when talking about aileron deflection for MLA, the authors always refer to the deflection used for the 2.5g pull-up manoeuvres δ 2.5g . The corresponding deflection for the −1g push-down manoeuvre is always equal to δ −1g = δ 2.5g / − 2.5. Positive deflections result in downwards rotation of the aileron and subsequent local lift increment.
Aileron deflection is simulated in two ways. The first way is via regular DLM panel associated to a control deflection parameter inside MSC.Nastran as in Bordogna et al. (2017). The second method relies on Reynoldsaveraged Navier-Stokes (RANS) CFD simulations to compute the forces associated with aileron deflection. These forces are then included in the MLA optimisation as external forces. The advantage of using CFD aileron forces is the ability to consider the non-linear effects between the aileron deflection and the resulting change in lift. Notably, the aileron, when deflected down, can "pull" the transonic shock toward the trailing edge. This may cause the flow to separate on the upper surface of aileron, and hence create a non-linear relation between the aileron deflection and its lift increment. Such effects can be captured using RANS simulations, as shown by Fillola (2006).

Optimisation
In this section, the optimisation problem and strategy are explained. As mentioned in Section 2, the use of lamination parameters requires a two-step approach where a continuous gradient-based optimisation is followed by an inverse optimisation problem where, from the optimum LPs, the corresponding stacking sequence is retrieved. The two steps are introduced in Sections 4.1 and 4.2 respectively.

Gradient-based optimisation
The first step is composed of the continuous optimisation where the optimal LPs and section thicknesses are obtained through a gradient-based process performed via MSC.Nastran SOL200 (2014).

Design variables
As introduced in Section 3.1 the wing structural elements are grouped into sections sharing the same thickness and material properties. In total, there are 44 sections, 14 for each skin and 8 for each spar (Fig. 5). Each section has as 9 design variables, one thickness t i and 8 lamination parameters v A i or v D i for a total of 396 design variables. Shear ribs and stiffeners do not take part in the optimisation and are made of quasi-isotropic composite. Carbon fibrereinforced polymer (CFRP) is used in the structure, and its material properties are summarised in Table 2.

Compatibility constraints
Lamination parameters have the capability to virtually model any composite materials in a continuous fashion; however, they cannot assume any arbitrary values. Compatibility constraints, also called feasibility constraints in the literature, ensure lamination parameters are compatible one with the others. Well-defined compatibility constraints that consider separately the four in-plane and out-of-plane lamination parameters have been proposed by Fukunaga and Sekine (1992). However, to the authors' best knowledge, there is a lack of definition for the compatible region that considers simultaneously in-plane and bending lamination parameters. One of the latest attempts in this direction has been made by Raju et al. (2014). In this work, the two abovementioned compatibility constraints are considered during the optimisation process.

Mechanical constraints
Two mechanical constraints have been considered during the optimisation of the regional wing: strength and local buckling.
Strength constraints have been derived from the work of IJsselmuiden et al. (2008). This approach defines an analytical expression for a conservative failure envelope based on the Tsai-Wu failure criterion in strain space. This conservative failure envelope determines a region in strain space that guarantees no failure would occur within a Static and dynamic aeroelastic tailoring with composite blending and manoeuvre load alleviation laminate, irrespective of the ply orientation angles present. Strength constraint is applied to all elements of the structure.
Local buckling is applied to all regions of the wing skins delimited by two consecutive ribs and stiffeners. The closed formula for biaxial loading of simply supported plate is used (20).
where buckling occurs for 0 < λ B < 1, N X and N Y are the stresses in the longitudinal and transverse directions, D 11 , D 12 , D 22 and D 33 are bending stiffness terms of the laminate, a and b are the corresponding region dimensions and m and n are the corresponding number of half waves. This formula does not take the effects of shear fluxes N XY and non-orthotropic composite bending stiffness effects (i.e. D 13 and D 23 non zero) into account. This approximation is known to be non-conservative and should be considered carefully (Bettebghor and Bartoli 2012). A safety margin of 50% is applied to both strength and local buckling criteria.

Blending constraints
Blending constraints, introduced in Section 2, are applied to all sections of a main structural component. The four main structural components are the two wing skins and the two spars. This means, for example, that blending is applied to all sections of the upper wing skin but not to the upper and lower wing skin simultaneously. The reason behind this choice is that usually these four components are produced separately and then fastened mechanically together. Therefore, there is no need to guarantee ply continuity among them. A shrinking factor equal to 0.5 is applied to all blending constraints to help the retrieval of a blended solution.

Aileron effectiveness constraints
Aileron effectiveness is a crucial parameter for aircraft handling quality. Aileron downwards deflection results in a local lift increment and also in a higher nose-down twist moment due to the lift being generated behind the elastic axis. If this increment in twist moment is not met by sufficient wing torsional stiffness, the resulting nose-down twist can diminish the intended lift generation. As a result, the total rolling moment is reduced and, in some extreme case, it can also be reversed. To ensure no aileron reversal, the ratio between the flexible and the rigid root bending moment is constrained.
where M rigid and M flexible are the rigid and the flexible root bending moment respectively. C min is constrained such as its value remains positive at 1.2 times maximum dynamic pressure speed.

Load-cases
The regional wing model is optimised with respect to different load-cases summarized in Table 3. In total, there are 19 load-cases: 8 symmetric static manoeuvres, 1 asymmetric manoeuvre load for aileron effectiveness and 10 gust load-cases.
The 10 gust loads, as it is explained in more details in Section 4.1.8, are the 10 worst gust loading conditions resulting from a selection process. The different gust profiles, taken from the certification CS-25 (2018), are shown in Fig. 7.

Optimisation problem and strategy
The blending constraints used in the gradient-based optimisation limit the change of lamination parameters between all sections as a function of their difference in thickness. Applying those constraints while simultaneously optimising thickness and lamination parameters leads to a non-convex optimisation problem (Macquart et al. 2016a). Therefore, the following 4-point strategy is employed during step 1 of the optimisation strategy (Fig. 8).
When the blending constraints are included in the optimisation (Fig. 8a), point 1 of this strategy is the conventional weight minimisation of the wing with respect to thicknesses and lamination parameters without blending constraints (22). This step provides a feasible starting point before the introduction of the blending constraints. The unblended design (x U ) is used as a starting point for point 2, where the objective is still to minimise wing mass but in this case, the blending constraints are included (23). In point 3, a repair function rounds up the thicknesses of blended design (x B ) to an even number of plies based on the ply thickness (t ply ). After the repaired design (x R ), lamination parameters are optimised one last time in point 4 while thicknesses are kept fixed with the intent of maximise the reserve factors of all the mechanical constraints (i.e. local buckling and strength). In this step, often referred to as margin maximization (Bettebghor and Bartoli 2012; Haftka and Watson 2005;Liu 2001), the feasibility of the structure is maximised utilising a slack variable s (24). Rounding of thicknesses and maximising of reserve factors modify the stiffness of the structure, leading to internal load redistribution. Therefore, points 2 to 4 are repeated until convergence to a final continuous design (x FC ). Next, step 2 performs a stacking sequence retrieval via GA to retrieve a blended final discrete design (x FD ).
In case blending constraints are not included in the optimisation (Fig. 8b), point 2 is not performed and point 4 does not consider blending constraints during the optimisation. Points 1, 3 and 4 are repeated until convergence.
The first step utilises the gradient-based optimiser of MSC.Nastran based on the interior point optimisation algorithm (IPOPT) (Wächter and Biegler 2006). Other optimisation algorithms could also be used. The approximation of the problem is obtained with convex linearization (CON-LIN) method and all constraints are required to be smaller that −1 for ease of comparison during the optimization runs. Manoeuvre load alleviation is introduced in the optimisation by setting a fixed value of aileron deflection for the pull-up 2.5g and push-down −1g load-cases. Once the optimum design is obtained, newly fixed deflections are set, and a new optimisation process is started. Two approaches for modelling the aileron deflection are used. The first one utilises the DLM method, while the second uses results obtained with RANS CFD.

Equivalent static load implementation
The implementation of the equivalent static load (Section 3.4) inside the optimisation strategy is summarized in Fig. 9.
The ESL starts with material properties of the structural model being updated. This ensures that the latest wing design is used during the equivalent static load creation. The structural stiffness matrix K(x) is also computed and stored for later usage.
MSC.Nastran SOL146 is used to solve the governing equations-of-motion for dynamic aeroelasticity (16) for the different gust profiles shown in Fig. 7. Four separated gust load analyses are performed at different flight points (i.e. sea-level and cruise) and for different mass configurations (i.e. MTOW and MLW). For each analysis, once the time domain response is obtained, several critical time instants are identified and their associated displacement fields u gust are extracted. Not all time instances are extracted to reduce the computational time associated with the ESL process. Two criteria are used to identify the critical time steps: the maximum wing root bending moment and wingtip maximum deflections. Ten critical gust instants for each gust analysis are extracted for a total of 40 stored displacement fields for all four gust analyses. As the wing is modelled as free-flying, the flexible structural displacements are obtained by removing the rigid body translations and rotations from the displacement vector of each grid points.
In order to include the 1g flight condition in the ESL, four separated static aeroelastic analyses are performed using MSC.Nastran SOL144, one for each combination of flight point and mass configuration. The corresponding aeroelastic displacement field u 1g is computed and extracted.
The total displacement is obtained by superimposing the 1g and gust displacement fields associated to the same flight and mass configuration. The gust displacement field is added or subtracted to the cruise 1g displacement to obtain the effect of a positive or negative gust. This operation brings the total number of critical gust displacements from 40 to 80.
Once u tot are retrieved, a set of equivalent static loads f eq can be computed using (26).
f eq allows a static analysis to achieve the displacement field that occurs during the dynamic response.
In order to further reduce the number static analysis associated with the different equivalent static loads, a second selection of the most critical f eq is performed before their introduction in the optimiser. The selection is made by

Stacking sequence retrieval
The OptiBLESS (Macquart 2016) open-source stacking sequence optimisation toolbox is used to retrieve manufacturable laminates. 1 OptiBLESS uses a guide-based genetic algorithm (GA) to retrieve blended stacking sequences matching the optimised lamination parameters achieved by the gradient-based optimiser (i.e. MSC.Nastran). According to the guide-based methodology (Adams et al. 2004), the thickest laminate is defined as the guide-laminate. Other laminates from the same structure are obtained by dropping plies from the guide-laminate, therefore ensuring the final design is blended. As a result, the plies of the thinnest laminate are ensured to span the entire structure. After the continuous optimisation, each wing section is optimised in terms of laminate thickness and lamination parameters. This design outcome is used as a target solution for the discrete stacking sequence retrieval. The thickest laminate within each main structural component (i.e. skins and spars) is identified and set as the guide. Next, the ply angles describing the guide laminate stacking sequence and ply drops are used as design variables in OptiBLESS. Doing so ensure some level of structural continuity between each of the substructure laminates. The genotype used in OptiBLESS to describe composite structures is given as: Ply angles The objective function used during the discrete optimisation represents the lamination parameter matching quality between the continuous and discrete retrieved design. In other words, OptiBLESS is set to retrieve blended stacking sequences with lamination parameters matching the lamination parameters obtained at the end of the continuous optimisation. This objective function is simply expressed as the root mean square error (RMSE) between the continuous and discrete lamination parameters as shown in (28) and (29).
where N lam is the total number of laminate sections in a structural component (i.e. upper wing skin), LP i,s is the vector of input parameters for section s, w i is a weighting factor and LP i,s is the vector of lamination parameters obtained by the GA. Stacking sequences are converted into lamination parameters in order to evaluate the fitness using the following notation: According to the fitness function given in (28), the bestretrieved stacking sequence would be a manufacturable stacking sequence best matching the optimised lamination parameters obtained by MSC.Nastran.

Effect on static and dynamic load criticality and thickness/stiffness distributions
In this section, the effects of the static loads, gust loads and blending constraints over the final design of the wing are evaluated. For this purpose, four different optimisation cases are analysed and they are presented in Table 4 together with their associated normalised optimal weight. In case A, the wing structure has been optimised only with respect to static loads and serves as a reference case. In case B, both static and gust loads are considered during the optimisation; thus, the difference between cases A and B reveals the effects of the gust loads on the final design. Finally, cases C and D have respectively the same loading conditions of cases A and B plus the blending constraints. Thus, the difference between cases A and C, and B and D shows the effect of blending constraints on the optimal design.
In Table 4, two different normalised weight results are given: W RU and W NRU . W RU represents the weight results obtained with the processes introduce in Fig. 8 where thicknesses are rounded up in point 3 of to the nearest upper multiple of t ply . This thickness "rounding up" is here referred to with the subscript (·) RU . On the other hand, W NRU is obtained with the same processes used for W RU , but in this case, point 3 is skipped. All weights are normalised with respect to the NRU weight obtained Fig. 10 a-d Load-cases associated with the most active mechanical constraints at continuous optimum for cases A, B, C and D from case A. The reason behind showing both RU and NRU weights is that the rounding up does not produce consistent results. This means that two designs with the same weights but different thickness distributions can result in different rounded up weights depending on the how close the different thickness distributions are to a multiple of t ply . Therefore, NRU weights are used in all weight comparisons while final designs comparison and stacking sequence retrieval are done using the RU results.
From Table 4, it is possible to notice that, for the composite regional wing model under study, the introduction of gust loads results in a marginal NRU weight increment. This small weight increment is either absorbed or augmented up to 0.5% by the rounding up, remaining minimal. The small impact of the gusts on the final design can also be spotted in Fig. 10 where the loadcases associated with the most active mechanical constraints (i.e. strength and buckling) are plotted. By comparing cases A with B and C with D, it is possible to notice that the introduction of the gust (cases B and D) does not result in large portion of the wing structure being sized by gusts. The blending constraints do not seem to change this behaviour; however, they tend to divide more neatly the region sized by 2.5g and −1g manoeuvre loads.
A closer look at the little effect of gust loads over the final continuous optimum can also be seen by comparing Figs. 11 with 12 and 13 with 14. These figures show the thickness distribution as well as the laminate membrane stiffness distribution associated with the membrane stiffness matrix coefficient A 11 . Laminate stiffness distributions are available for all the 12 stiffness matrix coefficients associated with a symmetric laminate; however, in this work, only the stiffness distribution of the coefficient A 11 is presented because it provides a visual interpretation of the main fibre direction (Dillinger et al. 2013). By comparing Figs. 11 with 12, it is possible to see that between case A and case B, there is no change in thickness distribution and small variation in membrane stiffness distribution. This shows that the introduction of the gusts has an impact on the final continuous design, but the new load-cases have been taken care of by a marginal change in fibre direction. Things are a bit different for cases C and D (Figs. 13 and 14), and in this case, the introduction of blending constraints does not allow a change in lamination parameters if there is no significant Concerning the effect of the blending constraints, it is immediately possible to see from Table 4 that the introduction of the blending constraints results in NRU weight increment of about 5% for both cases with and without gusts. This result is expected since these constraints reduce the design space as shown in Fig. 4. Another effect of the blending constraints is a more even distribution of the thickness across the wing. This more even distribution is clearly visible in Figs. 13 and 14 where the high thickness values of Figs. 11 and 12 are no longer reached and the chordwise thickness gap is reduced.
A final effect of the blending constraints can be seen in the membrane stiffness distribution. While this distribution for cases A and B varies significantly panel by panel from wing root to tip, this change in cases C and D is more smooth and limited. This occurs again due to the connection between thickness distribution and allowable change in lamination parameters introduced by the blending constraints.

Effect on manoeuvre load alleviation
Section 5.1 showed that, for the regional composite wing under study, the effect of the gust loads on the final design exists but it is marginal if compared with that of the static loads. This section analyses the effect of a manoeuvres load alleviation strategy for static loads on the final wing design.
In Bordogna et al. (2017), a similar study performed by the authors was carried out. In that preliminary work, the aerodynamic models used to perform the MLA study were DLM with HSA correction (Vincenzo 2012) for the wing aerodynamics and regular DLM for the lift due to aileron deflection-this method is referred here as HSA+DLM. In this work, DLM with HSA correction is again used to model the wing aerodynamics but the RANS CFD is used instead for the lift generated by aileron deflections-this method is referred to as HSA+CFD.  Fig. 15a, we compare the results obtained for case A using both HSA+CFD and HSA+DLM methods. It is possible to notice the linear behaviour of HSA+DLM; this is due to the linear formulation of the DLM method used for modelling the aileron deflection. On the contrary, the HSA+CFD behaviour is non-linear since the CFD is capable of depicting the non-linear effects between the aileron deflection and the resulting change in lift. As a result, the use of the linear HSA+DLM method results in significantly higher weight reductions (up to 12.7%) when compared with the non-linear HSA+CFD (up to 6.1%).
In order to take into account more realistic aerodynamic loads, especially high aileron deflection angles, the HSA+CFD model is used for the rest of the study. Figure 15b shows the effect of the MLA for the four different cases. The results show that MLA is capable of reducing the wing structural weight between 5.4 and 6.1% for an aileron deflection of δ = −30 • .
The weight reduction comes from the lift distribution of the static load-cases being moved inboard by the action of the ailerons and therefore resulting in lower wing root bending moment. For cases A and C, where the static loads are the only load-cases applied, the weight reduction should continue until the reduction in wing root bending moment is countered by the increment in torsion created by the action of the aileron. Alternatively, the weight reduction could also reach a limit if the CFD non-linear lift generation coming from the aileron deflection reaches a ceiling. For cases B and D, where both gusts and static loads are applied, another limit in the weight reduction could come from the presence of the gust loads that are not affected at all by the MLA. Therefore, while the influence of the critical static  Fig. 19 a, b Weight and stacking sequence retrieval RMSE for cases B and D with MLA at δ = −20 • load-cases is mitigated by the action of the ailerons, the gusts get more and more significant in the wing design and this results in larger areas of the wing being sized by gusts instead of static loads (Fig. 16). Figure 17 shows the optimized wing from case C, optimized without gust, and the wing from case D under the five most critical gusts. This figure is in agreement with the paper from Kenway et al. citeKenway2014 about the importance of considering dynamic loadcases.
The gaining influence of the gust over the final design is also visible in Fig. 15b where the gap between cases A and B for what concerns the W RU grows with the aileron deflection used during MLA. A similar behaviour is also visible between cases C and D; however, this effect is less pronounced due to the effect of the non-convex blending constraints that complicate the convergence of the optimisations.
Looking at Fig. 15b, it is possible to observe that the delta in weight due to the blending constraints remains mostly constant for various MLA settings and whether or not the gust loads are applied in the optimisation. Therefore, no apparent coupling effect between the type of load-cases and the use of blending constraints is observed.
Finally, as the gusts become more critical with MLA, the reliability of the ESL to produce a feasible structure needs to be assessed. The ESL method relies on a weak coupling between the gust simulations and the optimiser with gust loads considered frozen during the optimisation and updated at fixed iteration intervals. This means that at each iteration, the design satisfies the frozen loads, but there is no guarantee that the same design can also satisfy the gust loads updated with the latest wing design. For this reason, the optimal design is tested one last time in a dynamic aeroelastic analysis where the values of strength and buckling are checked for possible failure.
The results of this test are shown in Fig. 18 where the time history of the most critical strength and buckling constraints for the last optimisation iteration are compared with the gust check step for case D with MLA with aileron deflection at δ = −30 • . The results show that the most five most critical strength and buckling constraints do not change significantly between the last optimisation iteration and the gust check step, justifying the use of the ESL for wing sizing process. Out of all the 80 gusts considered in this study, this test shows an average constraint relative error

Stacking sequence retrieval
Once the continuous optimisation of step 1 converges to an optimal wing design, the stiffness matrix of each section of the wing is defined together with the thickness distribution. However, in order to manufacture each of the sections, it is essential to retrieve the stacking sequence associated with the obtained optimal stiffness and thickness distribution. This is achieved with the stacking sequence retrieval of step 2.
In this section, two continuous optimal designs obtained with and without the use of blending constraints are considered. The GA introduced in Section 5.3 is used for both cases to retrieve the stacking sequence that best matches the continuous optimum design. The two designs considered are one from case B and the other from case D and for both designs have been obtained with MLA at δ = −20 • and the second. During the GA, manufacturing constraints has been imposed in order to retrieve only symmetric laminates with ply orientation limited to multiples of 15 • .
The GA stacking sequence retrieval has been performed 10 times for each design, and the averaged RMSE obtained after stacking sequence retrieval is shown in Fig. 19a while the weight difference for both W RU and W NRU is summarised in Fig. 19b. As seen in Section 5.1, the introduction of the blending constraints results in a weight increment of about 5% and this fact also holds in this case. However, the linking in the change in lamination parameter to the difference in thickness among adjacent plies results in significant RMSE reduction between the target and the retrieved lamination parameters. This reduction of about 60% shows that the continuous design obtained with blending constraints is much closer to the discrete solution than the continuous design obtained without the blending constraints. If at the end of the continuous optimisation all constraints are satisfied, the retrieval of a structure that does not match 100% the continuous one is likely to result in a  Fig. 20, the constraint values for each element and buckling zone are presented for both design and for all load-cases. The blending constraints clearly help the continuous optimisation to obtain a design that is closer to a manufacturable design. The total number of failed elements together with the maximum failure indices for both strength and local buckling constraints for the two retrieved discrete structures is summarised in Table 5. In particular, for the number of failed elements, the reported number consists of the failed elements in all load cases. This means that if the same element fails in two different load-cases, it is counted twice. The smaller RMSE between the target and retrieved lamination parameters obtained thanks to the blending constraints significantly reduces the number of failed elements and the failure indices. In particular, failed elements are reduced to about 80.47% and 92.52% and maximum failure indices are reduced to about 41.87% and 35.60% for local buckling and strength respectively. Figure 21 shows the spanwise twist and displacement distributions for continuous design and the discrete retrieved one. It is possible to see that the use of blending constraints results also in a discrete structure that behaves more similarly to the optimal shape obtained via continuous optimisation. Finally, the retrieved stacking sequences of the upper and lower wing skins for the case D without MLA are shown respectively in Figs. 22 and 23.

Conclusion
This paper describes a full preliminary optimization strategy for an aircraft wing composite structure from flight envelop to the blended optimal stacking sequence. Several static and dynamic loads at various flight and mass conditions are considered. Dynamic loads are included in the optimisation as equivalent static loads, and manoeuvre load alleviation is used to mitigate internal structural loads during 2.5g and −1g static manoeuvres. Implementation of MLA has been performed using both DLM and CFD aerodynamics codes to model the lift increment generated by aileron deflection, while the rest of the wing is modelled using DLM corrected with HSA. The design variables used in the optimisation are the wingbox panel thickness and composite anisotropy characteristics. Furthermore, to ensure a certain degree of ply continuity throughout the structure, continuous blending constraints have been incorporated inside the optimisation strategy. At the end of the sizing process, an optimal blended stacking sequence is retrieved.
The objective of this paper is to evaluate how the introduction of continuous blending constraints affects the optimum design and the retrieval of the final stacking sequence for a regional aircraft wing. The effect of the blending constraints on the optimum design (1) focuses on static and dynamic loading conditions and identify the ones driving the optimization and (2) explores the potential weight saving due to the implementation of a manoeuvre load alleviation (MLA) strategy.
Results obtained in this study with a regional aircraft indicate that gusts have a limited impact on the final wing design, both in terms of thickness and stiffness distribution. The use of blending constraints did not significantly affect the load criticality; however, they tend to divide more neatly the region sized by the 2.5g and −1g manoeuvre loads. The use of MLA, when applied with high deflection angle, increases the criticality of gust loads in some areas of the wing. The use of CFD for aileron lift increment results in a non-linear relation between the optimised weight and the aileron deflection leading to less optimistic weight reduction than DLM. Overall, MLA contributes to a wing structural weight reduction between 5.4 and 6.1%. With the introduction of the blending constraints, the weight saving due to the MLA remains the same as well as the weight gained from using the constraints (due to the reduced design space). Therefore, no apparent coupling effect between the load criticality and the use of blending constraints has been observed.
Finally, blending constraints confirmed their significant influence on the final design as they strongly impact both thickness and anisotropy distribution (e.g. Figs. 11,12,13 and 14). The introduction of the blending constraints results in about 5% weight increment but also leads to optimum solutions that are closer to the retrieved discrete design obtained via GA. The gap between the desired design and the retrieved one is reduced by about 60%. This mismatch reduction results in fewer failed elements with lower failure indices for all the mechanical constraints and closer wing flight shapes.

Replication of results
Unfortunately, the FEM is restricted and cannot be shared.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.