{'fraction_non_alphanumeric': 0.035878729892961786, 'fraction_numerical': 0.01353425422073392, 'mean_word_length': 4.667645729778581, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-scale Deterministic Optimisation of Blended Composite Structures: Case Study of a Box-Wing

This work presents a multi-scale design methodology for the deterministic optimisation of thin-walled composite structures integrating a global-local approach for the assessment of the buckling strength and a dedicated strategy to recover blended stacking sequences. The methodology is based on the multi-scale two-level optimisation strategy for anisotropic materials and structures. In the first step, focused on the macroscopic scale, several design requirements are included in the problem formulation: lightness, feasibility, manufacturing, blending, buckling failure, static failure and stiffness. The second step, which focuses on the laminate mesoscopic scale, deals with the recovery of blended stacking sequences, for the structure at hand, matching the optimal geometric and elastic properties determined in the first step. As a case study, the unconventional PrandtlPlane box-wing system is used to show the effectiveness of the proposed design methodology.


Introduction
Composite materials allow for the tailoring of material properties, with several potential benefits. As known, in aerospace structures design, the most relevant design requirement is the lightness. In recent years, many works focusing on the development of design methodologies for composite structures have been carried out. However, an efficient use of composite materials to build large optimised structures remains one of the main challenges of aircraft industry [1]. In fact, if compared to isotropic metallic structures, the design of composite laminates introduces a higher level of complexity due to the increased number of design variables, the mathematical formulation of anisotropy, and the manufacturing and feasibility constraints formulation and modelling. All of these aspects must be considered since the preliminary design phase. The increase of design variables is mainly due to the laminate anisotropy, which shall be profitably exploited to obtain better performing structures. Two kinds of approaches are available in the literature for the design of composite structures: the direct approach and the multi-level one [2]. When adopting the first approach, the design problem is formulated by directly employing the plies orientations as design variables. This approach is quite straightforward and allows an easy implementation of many constraints and requirements; however, it presents some drawbacks. The optimisation problem formulated in the space of the plies orientations is highly non-convex and prone to reach just local minima. This is due to the trigonometric functions involved in the definition of the stiffness matrices of the laminate (see Sec. 2). The interested reader can find more details in [2,3]. Furthermore, the number of design variables is not generally known a priori when optimising for lightness. Moreover, it is dramatically size-dependent, which makes the approach not suitable for the design of large structures. In order to alleviate the drawbacks of the direct approach, in the multi-level one the design problem is split into two linked sub-problems, dealing with two different characteristic scales. In the first-level problem (FLP) the focus is put on the macroscopic scale and the laminates composing the structure are modelled as equivalent single-layer homogeneous anisotropic plates. The design problem is solved by optimising the parameters describing the macroscopic behaviour of the laminate (by means of a suitable representation); of course, all the requirements involved into the design problem must be taken into account at this level (through the formulation of equivalent constraints at the macroscopic scale of the laminate). In the second-level problem (SLP) the focus is put on the laminate mesoscopic scale and the aim is to find the stacking sequences (SSs) matching the optimised mechanical properties resulting from the previous level. The design variables of the SLP are the plies orientation angles. In this work, a multi-level approach is adopted.
As far as the mathematical description of the anisotropy of laminates is concerned, the vast majority of the works dealing with multi-level design strategies available in the literature makes use of the well-known lamination parameters (LPs) coupled with the parameters of Tsai and Pagano, see [4][5][6]. These parameters unquestionably provide a compact representation of the stiffness tensors of the laminate in the framework of the classic laminate theory (CLT); although, they are not all tensor invariants, as discussed in [6]. A sound alternative for describing the anisotropic behaviour of composite materials is represented by the polar formalism introduced in [7] by Verchery, and generalised to the case of higher-order equivalent single layer theories in [8][9][10] by Montemurro. Thanks to the polar formalism it is possible to represent any plane tensor by means of tensor invariants, referred to as polar parameters (PPs), which are directly related to the symmetries of the tensor. Some interesting works using LPs can be found in the literature [1,11,12]. The multi-scale two-level optimisation strategy (MS2LOS) based on the polar formalism has been originally introduced in [13,14] and has been later generalised, expanded and used in several works, such as [15][16][17][18][19][20][21].
Regarding manufacturing requirements, the so-called blending requirement, which consists of ensuring ply continuity between adjacent panels, is of paramount importance. A proper formulation of the blending requirement in multi-level optimisation procedures is one of the major challenges for composite structures design. In the framework of a multilevel approach, deriving a correct formulation of the blending constraints in either the LPs space or the PPs space is essential to ensure the possibility of finding manufacturable solutions in the SSs recovery phase. Macquart et al. [22] derived such constraints in the LPs space, although the proposed formulation is characterised by some weaknesses, as discussed in [23], where blending constraints have been formulated in the PPs space. Recently, narrower bounds for blending constraints in the PPs space have been derived in [21]. Indeed, the formulation of the blending requirement at the macroscopic scale of the laminate is beneficial also for the stacking sequences (SSs) recovery phase, which is performed in the SLP, to find fabricable solutions [24].
Due to the intrinsic multi-scale nature of composite structures, since different criticalities may appear at different scales, in the last decades dedicated global-local (GL) approaches have been developed, in order to balance the improvement in the structural description and the computational effort. In fact, more accurate structural descriptions require complex refined models, which in turn increase the total computational time. GL approaches try to find a trade-off among these two requirements. With a focus on aeronautic structural design, Arrieta and Stritz [25] developed a GL modelling strategy dedicated to damage tolerance analyses for conventional wings. Of course, this kind of analyses needs a refined model of structural components to simulate cracks growth. This represents a first example of the need of changing structural scale in optimisation procedures. A more complete approach is presented in [26]. However, the global model does not take stringers and spar-caps into account, since stiffened panels are modelled as equivalent shells. In [27], a GL approach for a high-speed wing is presented. The main issue is that local models are re-mapped to rectangular plane stiffened plates, loosing geometry effects on instability failures. Furthermore, several constraints are evaluated using analytic formulae. Liu et al. [28] presented a GL framework for optimisation of curvilinear spars and ribs (SpaRibs). The problem formulation presents a major issue: the procedure needs a significant computational effort (hundreds of cores) to find solutions in acceptable time. It is noteworthy that recently some efforts to formalise and include GL approaches into optimisation frameworks to be used for real complex structures design have been carried out [19][20][21]29]. In the last years, the optimisation of anisotropic stiffened panels including buckling requirements and manufacturing constraints, for aeronautic structures, is of renewed interest, as suggested by recent works [30][31][32]. This fact suggests the importance of the topic in the scientific literature.
To overcome the limitations related to classic design methodologies for composite laminates, a multi-scale methodology for the deterministic optimisation of thin-walled composite structures, based on the existing MS2LOS, is presented in this work. The strength of the adopted methodology rely on the mathematical formulation of the GL approach, which takes the scale transition contributions into account in the expression of the gradient of the structural response functions [20] and on an dedicated blending requirement formulation [21]. In so doing, optimal blended stacking sequences can be recovered matching the target elastic properties resulting from the FLP. Of Of Of course, this works also for moderately thick and thick laminates. As a modelling choice, the anisotropy is described by means of the polar method in the framework of the first-order shear deformation theory (FSDT), in order to consider the influence of the transverse shear stiffness of the laminate on the optimal solution.
As a meaningful case study, the PrandtlPlane (PrP) wing has been chosen (a pictorial view of the aircraft concept is shown in Fig. 1). The PrP configuration is the engineering application of the so-called "Best Wing System", introduced by L. Prandtl [33]. In particular, as a reference configuration, the PrP geometry developed within the EU funded research project PARSIFAL (Prandtlplane ARchitecture for the Sustainable Improvement of Future AirpLanes) [34] has been investigated. More details on the PrandtlPlane design within PARSIFAL project can be found in [35][36][37][38].
The paper is structured as follows. Sec. 2 briefly recalls the fundamentals of the polar method in the context of the FSDT. Sec. 3 summarises the proposed methodology. Sec. 4 describes the case study at hand, with a focus on the modelling hypotheses, loads and aircraft geometry. Sec. 5 presents the mathematical formulation and the numerical strategy of the FLP. The results of the FLP are presented in Sec. 6. Similarly, Secs. 7 and 8 present the counterpart for the SLP. Finally, Sec. 9 discusses the achieved results, whilst Sec. 10 concludes the paper with meaningful considerations and prospects.
2 Fundamentals of the polar method for composite laminates 2.1 First-order shear deformation theory In the FSDT framework [39], the expression of the stiffness matrix of a laminate (Voigt's notation) reads where A is the membrane stiffness matrix of the laminate, D is the bending stiffness matrix, H is the out-of-plane shear stiffness matrix, B the membrane/bending coupling  It is convenient to introduce also the following normalised matrices: where h is the total thickness of the laminate. In terms of geometrical and material parameters of the stacking sequence, in the special case of a laminate made of identical plies, i.e. same material and thickness for the elementary layer, the expressions of the above matrices are: In Eq. (3), N is the number of plies of the laminate, θ k is the orientation angle of the k-th ply, Q(θ k ) is the in-plane reduced stiffness matrix of the k-th ply whose material frame is turned by an angle θ k with respect to the global reference frame of the laminate. Analogously,Q(θ k ) is the out-of-plane reduced stiffness matrix of the k-th ply. Furthermore, coefficients b k , d k and c k read:

Polar formalism
The polar method allows representing any n-th order plane tensor in terms of invariants: this method was introduced, for the first time, by Verchery in 1979 [7]. For a deeper insight in the matter, the reader is addressed to [40].
In this framework, a second order symmetric plane tensor Z can be expressed in the local frame Γ = {O; x 1 , x 2 , x 3 } as: where T is the isotropic modulus, R is the deviatoric one and Φ is the polar angle. Among them, only T and R are tensor invariants, whilst Φ is needed to set the reference frame.
If L is a fourth-order plane elasticity-like tensor, i.e. with major and minor symmetries, its Cartesian components can be expressed through four moduli and two polar angles, namely T 0 , T 1 , R 0 , R 1 , Φ 0 , Φ 1 . The complete expressions read: In Eq. (6), T 0 and T 1 are the isotropic moduli, R 0 and R 1 are the anisotropic ones, Φ 0 and Φ 1 are the polar angles. Among them, only the four moduli and the difference Φ 0 − Φ 1 are tensor invariants.
A further advantage of the polar method is that for a fourth-order elasticity-like tensor, the polar invariants are related to the elastic symmetries of the tensor. Among the others, orthotropy corresponds to the condition Φ 0 − Φ 1 = K π 4 , K = 0, 1, whilst isotropy can be obtained if R 0 = R 1 = 0.
In order to properly analyse the mechanical behaviour of a laminate, it is possible to express the stiffness matrices of Eq. (3) in terms of PPs. In particular, A * , B * , D * , and thus C * , are fourth-order elasticity-like plane tensors, while H * behaves like a secondorder symmetric plane tensor. The PPs of the laminate stiffness matrices can be expressed as functions of the PPs of the lamina reduced stiffness matrices and of the geometrical properties of the stack (i.e. layer orientation, position and number). Moreover, as shown in [9,10], the deviatoric part of tensor H * can be expressed in terms of the PPs of tensor A * . Hence, the FSDT framework does not require further variables than CLT.
As a final remark of this section, it is noteworthy that membrane-bending uncoupling is often sought in many engineering applications. This property is achieved by imposing the condition B * = 0. Moreover, in order to have the same behaviour in terms of normalised membrane and bending stiffness matrices, i.e. the same elastic symmetry for A * and D * , the homogeneity property must be imposed. This condition reads C * = 0. In many practical applications, also the hypothesis of a fully orthotropic behaviour (both membrane and bending stiffness matrices) is considered. It is possible to show that, when orthotropy condition holds, the membrane PPs assume the following form [40]: where R A * 0K = (−1) K A * R A * 0 and K A * = 0, 1. Therefore, if the hypotheses of uncoupling, homogeneity and orthotropy hold, the design of a laminate is uniquely determined by only four variables: N , R A * 0K , R A * 1 and Φ A * 1 . Finally, for optimisation purposes, it is convenient to introduce the following dimensionless variables: 3 The multi-scale two-level optimisation strategy The goal of the MS2LOS is to provide a deterministic optimisation framework for complex composite structures by integrating structural responses involved at different scales. The final result of the MS2LOS is the set of blendend SSs satisfying the requirements of the problem at hand.
1. First-level problem (FLP). The aim of this phase, which focuses on the laminate macroscopic scale, is the determination of the optimal distribution of the mechanical and geometric design variables, describing the behaviour of each laminate, which minimises a given objective function by fulfilling a set of design requirements (formulated as optimisation constraints). At this level, the generic laminate is modelled as an equivalent homogeneous anisotropic plate, whose behaviour is described in terms of laminate PPs (see [40] for a primer). Thanks to this formulation, at this stage the designer can add further requirements (e.g. manufacturing constraints, strength and damage criteria, etc.) by introducing suitable constraints on the laminate PPs.
A first strong point of the FLP formulation proposed in this study is a multi-scale global-local approach for the assessment of both global and local structural responses, as proposed in [20]. Some portions of the structure, called zones of interests (ZOIs), which are likely to undergo buckling failure, are extracted from the global finite element model (GFEM) and analysed. Dedicated and refined local finite element models (LFEMs) are then generated, and buckling factors are evaluated. The GL modelling strategy used in this study is based on the sub-modelling technique [19,[41][42][43]. It is noteworthy that the assessment of the gradient of the buckling factors takes the scale transition between GFEM and LFEMs into account, which is a quite challenging task. In so doing, the mathematical formulation is suitable for deterministic optimisation.
A second strong point is the implementation of new blending requirements as a set of equivalent optimisation constraints on the laminate PPs, with the formulation recently proposed in [21].
Finally, an intermediate step, called discrete optimisation, is needed to get laminates with a discrete number of plies [21].
2. Second-level problem (SLP). The SLP of the strategy focuses on the laminate mesoscopic scale and aims at determining optimal SSs, satisfying the blending requirement between adjacent laminates, in such a way to recover the optimised PPs and thickness resulting from the FLP. The design variables of the SLP are the plies orientation angles.
The novelty of the SLP is the implementation of a numerical strategy to recover blended SSs matching the optimal elastic properties determined in the FLP, as recently presented in [21], without introducing simplifying hypotheses on the nature of the SSs. The solution phase, addressed through a metaheuristic algorithm, is suitable also for laminates with a large number of plies.
The FLP and the SLP are detailed in the following sections, through the presentation of a case study of a complex aircraft structure. For a deeper insight in the mathematical aspects related to the methodology, the interested reader is addressed to [20,21].

Case study
The PrP lifting system considered in this work is the result of a preliminary aerodynamic study presented in [44]. The PrP lifting system can be ideally split into three (semi-)wings: the front wing (FW), the rear wing (RW) and the vertical wing (VW), as shown in Fig. 2.
The global Body reference frame T B (CG; X B , Y B , Z B ), centred at the aircraft centre of mass CG is illustrated in the same figure. Due to the symmetry of the structure with respect to the aircraft longitudinal plane Y B = 0, the structural analysis is limited to the left-side part.

Geometry and material
In this study, the PrP wing-box architecture is optimised in the framework of the preliminary design phase of the aeronautic industry. During this phase, several loading conditions are considered to properly design the main components of the structure in order to comply with certification specifications [45]. Such load cases (LCs) result from the combination of basic loading conditions (BLCs) of different nature, e.g. flight loads due to symmetrical and asymmetrical manoeuvres or to gusts, ground loads, pressurization, etc. In this work, only a sub-set of these LCs is considered, as explained in Sec. 4.3.
For each wing of the PrP configuration, the external geometry is assigned in terms of the leading edge coordinates, chords and profiles of three reference sections, viz. root, kink and tip sections. As a result, the geometrical features affecting the out-plane shape of the wing, i.e. dihedral and twist angles, are automatically taken into account. All of the aforementioned quantities, which are summarised in Tab. 1, are known at the three sections (only two for the VW), and linearly vary between them. The same table adopts the reference frame T (O; x, y, z) where the origin O coincides with the orthogonal projection of the LE of the FW onto the aircraft symmetry plane Y B = 0, the x axis is parallel to the X B axis (but opposite direction) and the z axis is parallel to the Y B axis (but opposite direction). The airfoil F15−11, taken from [46,47], is considered in this study. This profile is used to describe the shape of the wing sections parallel to the free-stream direction.
The width of the wing box is obtained from front and rear spar positions, defined as chord percent according to Tab. 2. According to [46,47], a span-wise linear behaviour of both front and rear spars, for each wing, is adopted. The position of spars, with respect to the considered wing planform, is shown in Fig. 3 (in the T reference frame) and reported in Tab. 2.
Regarding the modelling of the structural components, the following simplifications are introduced: (i) only major structural components are modelled (viz. skin, stringers, ribs and spars); (ii) perfect bonding condition applies at the interfaces between the structural components; (iii) connection zones and opening/cut-out are neglected and not considered in the preliminary design phase; (iv) both GFEM and LFEMs, presented in Sec. 5.4, do not take explicit modelling of shear tie, stringer-tie and tear-strap into account.    For both FW and RW, ribs are parallel to the free stream direction between root and kink sections, whilst they are perpendicular, respectively, to rear and front spars between kink and tip sections (see Figs. 4a and 5). For the VW, ribs are perpendicular to both front Due to the lack of reference solutions for such an aircraft architecture in the literature, a preliminary genetic optimisation run has been performed to assess the values of some variables not included in the formulation of the deterministic optimisation strategy presented in Sec. 5, viz. stringers and ribs pitch, and stringers (and spar caps) geometry and stiffness. The problem formulation is conceptually similar to that of the FLP discussed in Sec. 5, but with less optimisation regions. For the solution search, the ERASMUS (EvolutionaRy Algorithm for optimiSation of ModUlar Systems) algorithm has been used (for more details on the algorithm, see [49]). Stringers and spar caps are modelled as Z-shaped stiffeners, as shown in Fig. 4b. Qualitatively, in the same figure, the manufacturing technique is shown: a basic SS (wherein plies corresponding to different orientation angles are represented through blue, red and green colours) is folded as to reproduce the Z-shaped cross section of the stringer. The preliminary optimisation calculation performed through the ERASMUS algorithm (not reported here for the sake of brevity) furnished n 0 = 0.3467, ρ 0 = 0.6165, ρ 1 = 0.1936, φ 1 = 0. These values of PPs and thickness (common for stringers and spar caps), together with stringers and ribs pitches, whose values are reported in Tab. 4, will be kept constant during the FLP of this work.
As done in [50], the VW is not subject to design in this work: the optimal properties from the genetic optimisation have been maintained. The stringers and ribs pitches are listed in Tab. 4, whilst the ribs, stringers and spar caps properties have already been presented above. The properties of the skin of the VW are n 0 = 0.5067, ρ 0K = 0.1390, ρ 1 = 0.0903, φ 1 = 0, whilst the properties of the spar webs are n 0 = 0.42, ρ 0K = 0.1617, ρ 1 = 0.0231, φ 1 = 0. Fig. 5 gives an overview of the resulting PrP wing-box geometry.

Design criteria
Only symmetric static loads are considered in this work. Certification specifications [45] identify two types of loading conditions: limit loads and ultimate loads. Limit loads are the maximum loads expected in service that the structure must withstand without detrimental permanent deformations. Ultimate loads are equal to limit loads multiplied by a prescribed factor (usually 1.5). The structure must withstand ultimate loads without failure for at least 3 seconds. For instance, for civil aircraft, limit loads in symmetrical manoeuvres occur at load factors n z = 2.5 and n z = −1. This study focuses on this class of loads.
The following set of design criteria (DCs) is integrated in the formulation of the optimisation problem.
• DC1: The global stiffness of the structure must be greater than a predefined reference • DC2: The laminate-level failure index, obtained by using the phenomenological failure criterion of Tsai-Hill expressed in the PPs space [51][52][53], for skin and spar webs, multiplied by a safety factor F S = 1.5 2 × 1.33, must be lower than or equal to the unit value. The contribution in the F S equal to 1.5 2 comes from airworthiness regulations (CS 25.303 [45]); note that the failure criterion is quadratic in the equivalent laminate strain, and that the safety factor applies to loads. The contribution equal to 1.33 considers the fact that the FI is applied at the laminate level, i.e. it is averaged over the thickness of the laminate, thus it can be too conservative with respect to the first-ply failure criterion.
• DC3: No buckling must occur in the skin panels when limit loads are applied, considering a safety factor F S = 1.5 × 1.1, where 1.5 comes from regulations and 1.1 comes from experimental considerations expressed in [18].
• DC5: Only manufacturable solutions are considered. To this end, the blending requirement between adjacent laminates is taken into account.
• DC6: Only homogeneous, uncoupled fully orthotropic laminates are considered. DC1 is expressed in terms of maximum tip displacement of the lifting system. DC2 is expressed in terms of the laminate-level failure index by using the Tsai-Hill failure criterion, evaluated in a proper subset of elements of the mesh in order to neglect the effects of local strain concentrations due to the GFEM modelling. DC3 is expressed in terms of no-buckling condition for the FW and RW dorsal and ventral skin (as discussed in Sec. 5.4.2). Of course, the evaluation of the first buckling load for such regions is done through dedicated LFEMs. DC4 is expressed by a set of inequalities defining the laminate feasible region in the PPs space. DC5 is expressed in terms of two inequalities, to be imposed to each couple of adjacent laminates (see [21]). DC6 is expressed by imposing

Load Cases
Aerodynamic loads are calculated through the Vortex Lattice Method solver AVL [54], implemented in the preliminary aircraft design tool AEROSTATE (AERodynamic optimisation with STAtic stability and Trim Evaluator) [55,56]. AEROSTATE has been used for the preliminary study of the reference configuration definition, together with medium and high fidelity corrections. For more details, the reader is addressed to [35,44,57]. The resulting lift distribution is interpolated and decomposed into a set of point forces and moments (applied to the centre of mass of the ribs of each wing) to obtain a staticallyequivalent system of forces. The aerodynamic loads are evaluated for a Mach number of 0.79, altitude of 11000 m, adopting the Standard Atmosphere model, in cruise condition. These loads define the fundamental basic load case (BLC) used in this study, which is denoted as BLC 1g . For the sake of simplicity, two load cases (LCs) have been considered. The first one (LC1) corresponds to a load factor n z = 2.5 (pull-up manoeuvre), whilst the second one (LC2) is characterised by n z = −1 (push-down manoeuvre), according to the flight envelope for civil transport aircraft. The lift distribution for all the previous cases is obtained by a simple scaling of the BLC 1g , which has been calculated for n z = 1. DCs 1, 2 and 3 are evaluated for the two LCs. Of course, DCs 4, 5 and 6 are independent from the LCs. It is noteworthy that the two considered LCs are commonly used in the preliminary design of civil aircraft wing-box structures at the boundary of the flight envelope [50,[58][59][60].

First-level problem formulation
In this study, the FLP deals with the minimisation of the mass of the PrP wing-box structure subject to feasibility, blending, static failure, stiffness and buckling requirements. The goal is to find the optimal distribution of thickness and PPs, characterising the laminates composing the main structural components, as a feasible minimiser of the problem at hand.

Continuous optimisation
The optimisation region is composed of dorsal and ventral skins of the FW and RW, and of the spar webs. For each skin, twelve regions are defined, labelled as shown in Fig. 6. Regarding the spar webs, they correspond to IDs 25 and 26 for the FW (leading and trailing edges side, respectively), and to IDs 51 and 52 for the RW (leading and trailing edges side, respectively).  Each panel is characterised by three design variables, introduced in Eq. (8). In fact, for the sake of simplicity, φ 1 has been set to zero: in this way, the main axis of orthotropy of each laminate is parallel to the stringers direction of each wing. The structure is then described by 156 variables, being 52 the optimisation regions, collected in the design variables vector ξ T := {n j 0 , ρ j 0K , ρ j 1 | j = 1, . . . , 52}.
The objective function is expressed as follows: where m 0 is the mass of the components which do not belong to the design region, m ref = 10000 kg is the reference mass, A i , n j 0 are the area and the dimensionless number of plies of the j-th panel, respectively, ρ ply , t ply are the density and the thickness of the single ply, respectively (see Tab. 3). Note that the objective function is the dimensionless mass of the whole wing-box (twice the modelled part).
In the FLP formulation, geometric and feasibility constraints, in terms of PPs, must be considered [61] to ensure that the optimal values of PPs correspond to a feasible laminate. This requirement corresponds to DC4. Under the aforementioned hypotheses, such constraints read [61]: The first two formulae in Eq. (11) can be considered as bounds for the relative design variables (see Tab. 5); therefore, the only feasibility constraint is represented by the third inequality (11), which must be evaluated for each panel. Formally: Since the Finite Element (FE) method is employed for the assessment of structural responses (see Sec. 5.3), for each LC a linear system of the form must be solved, where K is the stiffness matrix of the structure, u the nodal generalised displacements vector, f the nodal generalised external forces vector (for the particular LC). DC1 is expressed in the form of a constraint on the vertical displacement on a node of the FW tip section (as shown in Fig. 7), which must be lower than or equal to 0.15b, where b is the wing semi-span (b = 18 m). Formally: Each LFEM must be in equilibrium under the imposed BCs (Dirichlet's problem), as deeply discussed in [20]. With reference to the same work, the equilibrium reads where K is the stiffness matrix of the LFEM, K BC is the matrix obtained from the overall (singular) stiffness matrix of the structure by considering only rows and columns associated to nodes where BCs are imposed, u is the unknown displacement field of the LFEM, P is a matrix taking the mesh transition between the GFEM and the LFEM into account. Buckling constraint (DC3) involves the evaluation, for each LFEM, of the buckling factor, i.e. of the first eigenvalue of the problem where K σ is the geometrical stiffness matrix of the LFEM, λ and ψ ψ ψ are the first eigenvalue and the related eigenvector non-trivial solution of problem (16), respectively. In this study, for the sake of simplicity, DC3 is associated with LC1 for the dorsal skins of FW and RW, and with LC2 for the ventral counterparts. The solution of Eq. (15) is crucial for the assessment of the numerical strategy (see Sec. 5.3) and the evaluation of K σ , as discussed in [20]. Formally, for each LFEM, the buckling constraint reads As far as blending requirements (DC5) are concerned, for a generic couple of adjacent panels, labelled p and q, the blending constraint reads [21]: where ∆ pq (·) = · p − · q and N is the number of plies. These constraints must be imposed for each couple of adjacent panels.

Discrete optimisation
The discrete optimisation is needed to round up the continuous solution of problem (20) to discrete numbers of plies, while complying with the full set of constraints of (20). Let ξ ξ ξ c be the solution of the continuous optimisation problem (20). The discrete optimisation problem aims at finding a solution ξ ξ ξ d at a minimal distance from ξ ξ ξ c , having integer number of plies. The formulation of the mixed-integer constrained non-linear programming program (MICNLPP), presented in [21], reads: min ξ ξ ξ ξ ξ ξ c − ξ ξ ξ 2 L 2 , subject to : Problem (21) takes only blending and manufacturing constraints into account, while allowing for changes in the optimal PPs values in order to satisfy the set of constraints. This approach can affect the satisfaction of other constraints which do not enter directly in the formulation (21). The reader is addressed to [21] for more details on these aspects. Blending constraints are evaluated with the discrete number of plies N i = n 0d i N ref (i = p, q) for each generic couple of adjacent laminates p and q. The last constraint of problem (21) imposes that the difference of number of plies between two adjacent panels must be zero or greater than a predefined ∆N min (in this work, ∆N min = 4), to avoid impractical and meaningless drops of very few plies.

Numerical strategy
Problem (20) is a non-convex CNLPP in terms of both geometrical and mechanical design variables. The non-convexity is due to the set of constraints. The solution search of the continuous problem is performed via fmincon algorithm, available in the Optimisation Toolbox of MATLAB ® [62], tuned as summarised in Tab. 6. At each iteration, the MATLAB ® script invokes the Python routines which control the generation of the FE models and the evaluation of the objective function, of the constraints and of the corresponding gradients, as detailed in [20]. The optimisation procedure is coupled with ANSYS ® FE commercial software. Fig. 8 shows the conceptual work-flow of the numerical   strategy for the FLP of the MS2LOS. At each iteration of the deterministic algorithm, a new set of design variables ξ ξ ξ is passed to a Python Class which controls the GFEM generation. ANSYS ® is then invoked: an ANSYS ® Parametric Design Language (APDL) script generates the GFEM. Once problem (13) has been solved, for each LC, fundamental information such as connectivity, nodes coordinates, displacement fields, etc. are passed back to the Python Class. Then, the Python Class evaluates the objective function and all the optimisation constraints (except the buckling ones), together with their gradients. The LFEMs are generated in a similar manner. BCs are imposed by retrieving the displacement field from the GFEM results database. Matrix K σ is assembled within the Python Class, which also solves problem (16) via the scipy.sparse.linalg.eigsh routine. The use of this approach, extensively explained in [20], is mainly due to the perfect coherence between buckling factor and its gradient. The objective function value and all of the constraints (and gradients) are then passed to the optimisation algorithm. The loop is repeated until one of the convergence criteria of the fmincon algorithm is satisfied.
Once the continuous solution is obtained, the ant-colony algorithm MIDACO ® , specialised in mixed-integer programming [63][64][65][66], is used to perform the solution search for problem (21). MIDACO ® parameters have been tuned as reported in Tab. 7. A high Focus value facilitates the exploration of the neighbourhood of the starting guess, which is set to ξ ξ ξ c . The Oracle parameter, which is suggested to be the expected (or desired) objective function optimal value [67], is set to 0.

The global-local modelling approach
As stated above, the interaction between the FE models integrated in the optimisation process is based on a GL modelling approach. In particular, two different models are created. The GFEM is used to assess the macroscopic behaviour of each wing, whilst refined LFEMs are generated to properly evaluate the first buckling load of the most critical stiffened panels.
The GFEM size is of the order of 14000 elements and 50000 degrees of freedom (DOFs). These values have been obtained after a sensitivity analysis of the structural responses to the mesh size (this analysis is not reported here for the sake of brevity but it has been conducted following the same steps presented in [29,50]). Link and shell elements are connected together by node merging. Shear-tie components are not modelled, but their mechanical effect (the transfer of shear load from ribs to skin) is ensured by the direct connection between ribs and skin elements. As shown in Fig. 5, fillets connecting FW and RW to VW are not explicitly modelled; MPC184 (multi-point constraint) elements with "rigid beam" behaviour (master-slave approach) are used to link extremal ribs nodes with the central master node, as illustrated in Fig. 9a. Aerodynamic forces and moments are applied to a reference node, which, for (a) Connection between FW and VW (b) Force and moment application to ribs Figure 9: GFEM mesh particulars the sake of simplicity, is generated at the centroid of each rib; the master node is then linked to the boundary nodes of each rib via RBE3 elements, as shown in Fig. 9b.
Since AEROSTATE-AVL gives as an output the position of the pressure centre, for each of the strips in which the lifting surface is subdivided, transport moments can be easily evaluated. As far as BCs of the GFEM are concerned, the 6 DOFs of the nodes lying at the root rib of both FW and RW are set to zero.
After solving the GFEM, for each LC, the ZOIs are generated. Bays close to the root section, to the connection with the VW and to the kink transition rib are disregarded. Furthermore, the three closest areas to front and rear spars are not included within the ZOI. As an example, Fig. 10a shows the ZOI for the dorsal skin of the FW. From the ZOIs, checking zones element subsets are created, wherein DC2 is evaluated. As an example, Fig. 10b shows the checking zone for the dorsal skin of the FW.
Results provided by the GFEM are used for the evaluation of the objective function and all the constraint functions except those related to buckling failure.

The local finite element model
LFEMs are generated in order to evaluate structural phenomena that typically appear at a smaller scale with respect to the GFEM one. In this work, the LFEMs are used to assess the first buckling load of the dorsal and ventral skins of the FW and RW. The LFEM are represented by the ZOI illustrated in Fig. 10a whose mesh is refined, as shown in Fig. 11, to correctly assess the first buckling factor. The structural components are entirely modelled by using SHELL181 elements. Each LFEM has approximately 7000 elements and 43000 DOFs (as in the case of the GFEM, these values are the outcome of a sensitivity analysis of the first buckling load to the mesh size not reported here for the sake of brevity). Dirichlet-type BCs are interpolated from the GFEM, for the predefined LC, and applied to the nodes of the LFEM located in correspondence of stringers and ribs intersections with the skin, as illustrated in Fig. 11. An eigenvalue buckling analysis is performed on the To provide a clearer picture of the optimal solution, the distribution of the design variables over the structure for dorsal and ventral regions of both FW and RW is illustrated in Figs. 12-14. From the analysis of these figures and of Tab. 8, one can infer that the panels composing the skin are characterised by different types of orthotropy. Several regions (e.g. panels 1, 2, 5, etc.) are characterised by a standard orthotropy shape (the value of ρ 0K is positive and ρ 0K and ρ 1 are of the same order of magnitude), whilst other laminates (e.g. panels 3, 8, 9, 10, etc.) are characterised by the so-called "dog-bone" orthotropy shape [9] (the value of ρ 0K is negative and ρ 0K and ρ 1 are of the same order of magnitude). Furthermore, one can notice that few regions (e.g. panels 28, 29, 31, 33, 34, 35 and similar) are characterised by a square symmetry orthotropy (the value of ρ 1 is at least two order of magnitude lower than ρ 0K ), which means that the laminate show equal elastic properties along its two main axes of orthotropy. It is noteworthy that only panel 23 is characterised by an almost isotropic behaviour (because both ρ 0K and ρ 1 are almost null).
The first two columns of Tab. 9 list the values of the constraints for the optimal solution for both the continuous and the discrete optimisation. The active constraints in the continuous optimisation are the blending one and the buckling one for the ventral skin of the FW. In particular, the latter is slightly violated, but still in tolerance (see Tab. 6). The stiffness constraint, i.e. the requirement on the maximum vertical tip displacement, is verified with a wide margin. This means that the structure is very rigid and that this requirement is not so restrictive in the region where the local optimum, found by the algorithm, is located. The discrete optimisation step generally makes the constraint values to assume more negative values, with the exception of the constraint related to the buckling of the dorsal skin of the FW, which is violated (however, the corresponding  buckling factor is 1.3993).

The second-level problem formulation
The goal of the SLP is to search for at least one blended SS, for the whole structure, recovering at the same time the optimal value of thickness and of PPs resulting from the resolution of the FLP. The process of recovering the stiffness properties will be referenced as stiffness recovery (SR). The numerical strategy is the one presented in [21].

Mathematical formulation
The basic step for the mathematical formulation of the SLP is the formalisation of the SR for a single laminate. Of course, in this particular case, the design variables vector coincides with the SS of the laminate. Therefore, there are as many independent design variables as the number of plies of the laminate. In the following, target PPs are labelled with the symbol † : ρ † 0K , ρ † 1 , φ † 1 . Nevertheless, it is more convenient to consider the following equivalent set of parameters in place of ρ † 0K : As discussed in [9], it is useful to consider the following dimensionless quantities: and In Eq. (23), M is a suitable norm that can be defined in different ways. In this work, it is assumed equal to [72]: A SS for which R = 0 will be referenced to as Recovery Stacking Sequence (RSS). Equation (24) represents a suitable distance from a laminate having desired mechanical properties. Of course, its extremal value is zero, reached when all the terms R i , i = 1, . . . , 6 are null, i.e. when all the target properties are matched. The physical meanings of R i are straightforward: R 1 represents the uncoupling condition, R 2 the homogeneity condition, R 3 the orthotropy condition, R 4 , R 5 and R 6 the conditions on target PPs [9].
At the mesoscopic scale, the blending constraint is merely a decision problem whether two adjacent panels SSs are blended or not. In this work, such problem is addressed through the proper algorithm presented in [21]. Therefore, the search of blended RSSs for the whole structure can be formulated as follows: blending constraints at the meso-scale (see [21]) (26) Problem (26) is a CNLPP. As widely discussed in [21], problem (26) can be transformed in an unconstrained non-linear programming problem imposing blending by construction.
In this work, in order to have a planar and continuous top surface of each skin region, the blending scheme shown in Fig. 15 is adopted. The reason is to have no discontinuities, due to ply drop, on the outer side, which may penalise the aerodynamic performances. Therefore, the top surface of each laminate is wet by the external air flow, whilst the Figure 15: Blending scheme bottom side overlooks the wing-box inner part. Moreover, the two innermost plies (yellow and orange in Fig. 15) are shared by all the laminates as a covering [73]. The determination of the orientation values must be carried out to minimise the global residual function R for the structure.

Numerical strategy
Problem (26) has been split into eleven sub-problems: one for each of dorsal and ventral skins of FW and RW, one for each of spar webs and stringers, one for the VW skin and one for the VW spar web. For each of these eleven sub-problems, Tab. 10 lists the number of independent variables. The number of independent variables for dorsal and ventral skins of both FW and RW is determined by considering the variables propagation schemes illustrated in Figs. 16 and 17. These schemes are used for the resolution of problem (26) and are determined according to the strategy described in [21].
The solution search of problem (26), properly adapted to each of the sub-problems, is performed via MIDACO ® optimisation software, with the parameters listed in Tab            "Considering the optimal composite solution presented in this work (Sec. 6), the mass saving is approximately 30% if compared with the conventional light-alloy solution presented in [50], where it was considered the same geometry with a different approach to the design problem. This is aligned with results presented in [74], where it was shown, with low-fidelity techniques, that the composite box-wing is about 30% lighter than the metallic counterpart. By looking at Tab. 9, some constraints values suggest that the solution may be further improved. Indeed, the complexity of the structural model and of the problem formulation may have created some difficulties in the convergence of the algorithm. However, as expected, the solution lies on the feasibility domain boundary. This can be inferred from the values of the blending constraint and of the buckling constraint for the ventral skin of the FW, which take an almost null value. As discussed in Sec. 6, the discrete optimisation effect results in an increase of one buckling constraint, which assume a positive value. It is clearly due to the stress redistribution, which is not considered in the formulation of the discrete optimisation problem: this is an expected result, as deeply discussed in [21]. As far as the SLP results are concerned, it can be appreciated that blended skin laminates (Tabs. 12-15) achieve higher residuals values than stand-alone laminates (Tab. 16). The implementation proposed in this work, when dealing with single laminates (no blending), achieve very small residuals (ranging in 10 −7 − 10 −5 ), comparable for example with those achieved by a genetic algorithm in [9,14]. The blended case is indeed more complex. The residuals range in 10 −4 − 10 −1 , with an average residual of 10 −2 . The priority, in the presented approach, is the satisfaction of ply continuity between adjacent panels. This enforcement causes a shrinkage of the admissible design region, so that the number of independent orientations to be determined is not merely the sum of the number of plies appearing in each skin laminate. Accordingly, the SLP results can be interpreted in two ways. The first one is the closeness of the recovered SSs to the target PPs. Unfortunately, in the literature there is not a unified index of merit to compare different approaches. Apart some outliers, the proposed residuals are generally small, so that the proposed strategy seems to be acceptable. A second way to interpret the SLP results is to check a posteriori that the recovered SSs are actually compliant with all the constraints of the problem at hand. The third column of Tab. 9 shows the residuals of the original optimisation problem at the design point defined by the recovered SSs. The proposed solution in terms of PPs satisfies feasibility, stiffness and strength requirement, the blending one being trivially satisfied by construction. However, the buckling constraints on the dorsal skins are violated. The SLP solution impacts up to the 16% in the buckling response assessment. This may be due, again, to the stress redistribution, no more taken into account after the Continuous Optimisation due to the intrinsic nature of the multi-step approaches, and to the error introduced by the stack recovery phase. Research is ongoing in this direction, trying also to mitigate such impact.

Conclusions
This work proposes an enhanced multi-scale two-level optimisation strategy for the deterministic optimisation of composite thin-walled structures. The main two ingredients are the implementation of blending constraints within the optimisation problem formulation and a global-local approach for the assessment of buckling response. In this context, a rigorous mathematical formulation of the scale transition between global and local models makes the strategy suitable for deterministic algorithms. In addition, an approach for the explicit recovery of blended stacking sequences completes the methodology, allowing for the effective manufacturing of the structure at hand. The methodology has been applied to a case study representative of real-world industrial applications, i.e. the optimal design of the PrP lifting system structure, to show the effectiveness of the approach when dealing with the design of complex aircraft structures. A light structure design has been found, and general blended stacking sequences have been recovered also for laminates with large number of plies, thus providing an explicit manufacturable solution. In the future, the methodology could be extended to deal with other elements formulations, i.e. shell elements whose kinematics is described by higher-order theories and by including non-linear analyses (e.g. post-buckling behaviour of the most critical regions, damage mechanics, etc.) in order to generalise the proposed modelling approach and to find improved (i.e. lighter and more efficient) optimal solutions.