G A Fiordilino

M I Izzi

M Montemurro
email: marco.montemurro@ensam.eu

A general isogeometric polar approach for the optimisation of variable stiffness composites: application to eigenvalue buckling problems

Keywords: Composite materials, Variable Stiffness Composites, Polar method, B-Spline surfaces, Optimisation, Buckling

This study presents a general approach for the multi-scale design of variable stiffness composites (VSCs). The first-level problem of the multi-scale two-level optimisation strategy (MS2LOS) is solved to determine the optimal distribution of the VSC stiffness properties at the macroscopic scale satisfying the requirements of the problem at hand. In this phase, the VSC laminate is modelled as an equivalent homogeneous anisotropic plate whose behaviour is described in terms of polar parameters (PPs), which vary locally over the structure. The First-order Shear Deformation Theory is used to take into account the influence of the transverse shear stiffness on the mechanical response of the VSC and Basis Spline (B-Spline) surfaces are employed to represent the PPs fields. In this background, the expression of the gradient of the buckling factor is determined analytically by exploiting the properties of the polar formalism and of the B-Spline surfaces. Moreover, the effect of the discrete variables, involved in the definition of the B-Spline surfaces, on the performances of the optimised solution is investigated. The effectiveness of the approach is proven on two benchmark problems dealing with the maximisation of the first buckling load of a VSC laminate, subject to feasibility and geometric requirements, taken from the literature. The results obtained by means of the MS2LOS based on the polar formalism outperform those reported in the literature, which are obtained through an optimisation strategy based on lamination parameters.

Introduction

Nowadays, composite materials are widely used in the aerospace field because of their outstanding properties in terms of lightness, stiffness, strength, durability, etc. The need of having cheap and efficient structures and the development of new manufacturing process, e.g. automated fibre placement (AFP) machines, or recent additive manufacturing (AM) technologies, like Fused Filament Fabrication (FFF) and Continuous Filament Fabrication (CFF) processes, are driving the research towards new solutions and products [START_REF] Wang | 3d printing of polymer matrix composites: A review and prospective[END_REF]. These innovative solutions, whose properties (mechanical, thermal, etc.) vary pointwise over the structure, are know as variable stiffness composites (VSC) and they show better performances than those characterising standard straight-fibre format solutions. However, the large amount of design variables characterising VSCs requires sophisticated design/optimisation methodologies, as highlighted by Ghiasi et al. [START_REF] Ghiasi | Optimum stacking sequence design of composite materials part ii: Variable stiffness design[END_REF]. VSCs can be obtained in different ways. One of the most interesting and commonly employed solutions is the variable angle tow (VAT) composite manufactured through the AFP or the FFF+CFF technologies [START_REF] Dirk | The engineering aspects of automated prepreg layup: History, present and future[END_REF]. These processes allows placing the fibre (i.e. the tow) along a curvilinear path within the constitutive lamina, thus implying a point-wise variation of the material properties. In this way the designer can take full advantage of the curvilinear format by orienting the tows according to the local stress field in the most effective way. AFP technology guarantees a great flexibility in building composite laminates with relatively complex geometries, but it is sensitive to the manufacturing defects as discussed in Kim et al. [START_REF] Kim | Continuous tow shearing for manufacturing variable angle tow composites[END_REF][START_REF] Kim | Manufacturing characteristics of the continuous tow shearing method for manufacturing of variable angle tow composites[END_REF], where the AFP technology is compared to the CTS (continuous tow shearing). In these works, the advantages and drawbacks of both AFP and CTS technologies are discussed by putting the accent on manufacturing defects. Nik et al. [START_REF] Nik | Optimization of variable stiffness composites with embedded defects induced by automated fiber placement[END_REF] studied the influence of the most common defects related to the AFP process, i.e. gaps and overlaps, on the buckling load and in-plane stiffness of VAT laminates.

If, on the one hand, VAT composites present superior performances with respect to conventional solutions, on the other hand, some critical issues exist. The design/optimisation process is rather complex due to the large number of design variables involved at each characteristic scale [START_REF] Setoodeh | Design of variable-stiffness laminates using lamination parameters[END_REF][START_REF] Setoodeh | Design of variable-stiffness composite panels for maximum buckling load[END_REF], the damage mechanisms are not well-understood, neither damage tolerance nor dedicated failure criteria are available in the literature, etc. As discussed in [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF], these issues are related to two intrinsic properties of composite materials, i.e heterogeneity and anisotropy, that occur at different scales of the structure and vary point-wise in the case of VAT composites. Up to now, no general guidelines are available in the literature for the multi-scale design of VAT composites. The work by Hyer and Lee [START_REF] Hyer | The use of curvilinear fiber format to improve buckling resistance of composite plates with central circular holes[END_REF] represents the first attempt to explore the advantages that can be achieved by using a curvilinear fibres-path instead of the conventional straight fibre format. The authors make use of a sensitivity analysis based on finite difference and a deterministic algorithm (DA) to determine the optimal fibres-path that maximise the first buckling load of a plate with a circular hole in the centre. Gurdal and Olmedo [START_REF] Gürdal | In-plane response of laminates with spatially varying fiber orientations-variable stiffness concept[END_REF] presented a preliminary work on VSCs, where only the in-plane stiffness response is considered: a linear variation of the fibre angle is assumed and closed-form solutions are derived. The effectiveness of stiffness tailoring through the use of curvilinear fibres-path is also analysed by Jegley et al. in [START_REF] Jegley | Optimization of elastically tailored tow-placed plates with holes[END_REF], where a considerable reduction of the stress concentration around the holes and a remarkable improvement of the load carrying capacity of the plate are observed.

The first buckling load is one of the most important design criteria when dealing with the design of lightweight composite structures. As far as the optimisation of VSCs is concerned, two different approaches are available in the literature [START_REF] Ghiasi | Optimum stacking sequence design of composite materials part ii: Variable stiffness design[END_REF]. The first one is the direct approach where the parameters tuning the shape of the fibres-path are considered as design variables. The second one is the multi-level approach wherein the optimum design of a VSC is split in two sub-problems. At the macroscopic scale the VSC is represented as an equivalent single layer plate and its mechanical response is described by a set of suitable parameters, related to the VSC stiffness matrices, which constitute the design variables. At the mesoscopic scale (i.e. the fibres-path level) the goal is to determine the optimum fibrespath, in each lamina, satisfying the optimum distribution of the macroscopic mechanical variables. Several authors have worked on this topic by showing the great potential behind VSCs [START_REF] Olmedo | Buckling response of laminates with spatially varying fiber orientations[END_REF][START_REF] Nagendra | Optimization of tow fiber paths for composite design[END_REF][START_REF] Waldhart | Analysis of tow placed, parallel fiber, variable stiffness laminates[END_REF]. The direct optimisation approach has been applied by Gurdal et al. [START_REF] Gürdal | Variable stiffness composite panels: effects of stiffness variation on the in-plane and buckling response[END_REF] who dealt with the effects of the stiffness variation on the in-plane and buckling response in flat rectangular VAT laminates employing basic fibres-path definitions based on 1-D functions [START_REF] Gürdal | Variable stiffness composite panels: effects of stiffness variation on the in-plane and buckling response[END_REF]. Lopes et al. [START_REF] Lopes | Variable-stiffness composite panels: Buckling and first-ply failure improvements over straight-fibre laminates[END_REF] proved that buckling and first-ply failure loads can be considerably increased with respect to standard solutions. An alternative representation of the fibres-path is proposed by Wu et al. [START_REF] Wu | Buckling analysis and optimisation of variable angle tow composite plates[END_REF]: a modified version of the Rayleigh Ritz method is implemented to study the buckling response of a VAT plate under uni-axial load with mixed boundary conditions (mechanical and geometric) and a genetic algorithm (GA) is used to find the optimal fibres-path maximising the first buckling load. A differential quadrature methodology based on the Airy's stress function was developed to study the pre-buckling and post-buckling behaviours of VAT plates subject to axial compression and mixed boundary conditions in [START_REF] Raju | Prebuckling and buckling analysis of variable angle tow plates with general boundary conditions[END_REF][START_REF] Raju | Postbuckling analysis of variable angle tow plates using differential quadrature method[END_REF]. Coburn et al. [START_REF] Coburn | Buckling analysis of stiffened variable angle tow panels[END_REF] proposed an analytical model for the buckling analysis of blade stiffened VAT panels. The results provided by this model are in good agreement with those resulting from a finite element (FE) analysis.

Regarding the multi-level approach for optimising VSCs, recently, a review on the optimisation methodologies for VSCs making use of lamination parameters (LPs) has been presented in [START_REF] Albazzan | Efficient design optimization of nonconventional laminated composites using lamination parameters: a state of the art[END_REF]. As known, the use of LPs, as design variables, reduces the computational cost of the problem essentially because the number of variables is strongly reduced. Setoodeh et al. [START_REF] Setoodeh | Design of variable-stiffness composite panels for maximum buckling load[END_REF] presented a generalised reciprocal approximation for the design of VSC panels for maximum buckling load. In this formulation, the buckling load is approximated using a first-order Taylor series expansion in terms of the point-wise compliance tensors. The related maximisation problem is reduced to a simple local optimisation at any discretisation point. Ijsselmuiden et al. [START_REF] Ijsselmuiden | Optimization of variable-stiffness panels for maximum buckling load using lamination parameters[END_REF] dealt with the maximisation of the first buckling load of VSCs through a conservative approximation method that expresses the buckling load as a linear combination of the in-plane and bending stiffness tensors and the corresponding inverse forms. The laminate stiffness matrices are written in terms of the parameters of Tsai and Pagano and of the LPs (which are the design variables of the problem) and the adjoint method is used to determine the gradient of the buckling factor. Liu et al. [START_REF] Liu | Buckling optimization of variable-angle-tow panels using the infinite-strip method[END_REF] used a DA to minimise the mass of VAT laminates under buckling and manufacturing constraints. The buckling analysis is performed by the infinite-strip method. Wu et al. [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF] presented a two-level optimisation procedure where LPs are used as design variables. B-Spline surfaces are used to represent the LPs fields over the VAT laminate. The design problem was formulated as the maximisation of the first buckling factor subject to feasibility constraints on LPs. A multilevel approach to optimise the fibres-path at each node of a FE model is presented in [START_REF] Peeters | Stacking sequence optimisation of variable stiffness laminates with manufacturing constraints[END_REF]. In this work, during the first step of the procedure, the VAT laminate response is described in the LPs space, whilst during the second step GaussNewton quadratic approximation is used to retrive the optimised fibres-path corresponding to the optimum LPs found during the first step. To ensure manufacturability, a constraint on the gradient of the fibre angle between nodes is imposed. Nevertheless, in this way there is no warranty to find optimum fibres-paths satisfying both the manufacturing constraints and the optimised LPs distribution resulting from the first step. Moreover, the use of simplifying hypotheses as done in [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF][START_REF] Peeters | Stacking sequence optimisation of variable stiffness laminates with manufacturing constraints[END_REF], like the use of symmetric stacks to obtain membrane-bending uncoupling or the use of balanced sequences to get membrane orthotropy, systematically leads to misleading solution (e.g. a symmetric balanced stack is totally anisotropic in bending) and extremely shrinks the design space.

To go beyond these limitations, a general multi-scale two-level optimisation strategy (MS2LOS) for the optimum design of VSCs has been presented by Montemurro and Catapano in [START_REF] Montemurro | Variational analysis and aerospace engineering : mathematical challenges for the aerospace of the future[END_REF]. The MS2LOS relies, on the one hand, on the use of the polar formalism introduced by Verchery [START_REF] Verchery | Les invariants des tenseurs dordre 4 du type de lélasticité[END_REF] and later extended to the case of higher-order theories [START_REF] Montemurro | An extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | Corrigendum to" an extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | The polar analysis of the Third-order Shear Deformation Theory of laminates[END_REF] for the description of the anisotropic behaviour of the composite. On the other hand, the MS2LOS is based on the use of B-Spline and NURBS surfaces to represent both the laminate polar parameters (PPs) fields over the structure and the fibres-path within the generic lamina. Moreover, the MS2LOS makes the use of a hybrid optimisation tool composed by the union of a classical DA and an original GA able to deal with a special class of optimisation problems defined over a domain of variable dimension, see [START_REF] Montemurro | A contribution to the development of design strategies for the optimisation of lightweight structures[END_REF][START_REF] Costa | A General Hybrid optimization Strategy for Curve Fitting in the Non-Uniform Rational Basis Spline Framework[END_REF][START_REF] Montemurro | Design of elastic properties of laminates with minimum number of plies[END_REF][START_REF] Montemurro | A general multiscale two-level optimisation strategy for designing composite stiffened panels[END_REF][START_REF] Panettieri | Blending constraints for composite laminates in polar parameters space[END_REF][START_REF] Montemurro | Least-weight composite plates with unconventional stacking sequences: design, analysis and experiments[END_REF][START_REF] Cappelli | cation of the viscoelastic behaviour of composite materials through a non-destructive test[END_REF][START_REF] Delucia | A stochastic approach for predicting the temperature-dependent elastic properties of cork-based composites[END_REF]. Recently, the mathematical framework of the MS2LOS has been enhanced by developing a suitable formulation of the minimum admissible radius of curvature of the tow in terms of the VAT laminate PPs in [START_REF] Montemurro | On the effective integration of manufacturability constraints within the multi-scale methodology for designing variable angle-tow laminates[END_REF]. Montemurro and Catapano [9] also derived, in a closed form, the expression of the gradient of the VAT laminate strain energy and of the feasibility constraint by exploiting the properties of the B-Spline surfaces.

This study focuses on the first-level problem (FLP) of the MS2LOS, that aims at determining, at the macroscopic scale, the optimum distribution of the VAT laminate PPs in order to satisfy the requirements of the problem at hand. At the macroscopic scale, the VAT laminate is modelled as an equivalent homogeneous anisotropic plate whose mechanical behaviour is described in terms of PPs (which vary locally over the structure). The First-order Shear Deformation Theory (FSDT) is used to take into account for the influence of the transverse shear stiffness on the overall mechanical response of the VAT composite in the Finite Element (FE) model and B-Spline surfaces are employed to represent the PPs fields. The advantage of using the Vercherys polar method [START_REF] Verchery | Les invariants des tenseurs dordre 4 du type de lélasticité[END_REF] is in the fact that the elastic response of the structure at the macro-scale is described in terms of tensor invariants, i.e. the PPs, having a precise physical meaning (related to the elastic symmetries of the material).

In this study, the MS2LOS is applied to the maximisation of the first buckling load of a VAT composite subject to feasibility constraints. The main contribution consists in the derivation of the expression of the gradient of each physical response in terms of the design variables involved in the definition of the B-Spline surfaces representing the PPs distributions over the structure. In particular, the derivation of the analytical form of the gradient of the first buckling factor is anything but trivial and exploits two main properties of B-Spline blending functions: the strong convex-hull property and the local support property. Thanks to the B-Spline surfaces formalism, the optimised solution are no longer related to the mesh of the FE model, rather they are fully CAD-compatible and can be directly passed to the G-Code behind AFP and AM processes for manufacturing purposes. The effectiveness of the proposed approach is proven on two meaningful benchmarks taken from the literature.

The paper is organised as follows: the overview of the design problem and a brief reminder on the MS2LOS are reported in Section 2, while the mathematical formulation of the first-level optimisation problem is presented in Section 3. The FE model of the VAT laminate plate, for both benchmarks, is described in Section 4, while the numerical results and discussed in Section 5. Finally, Section 6 is devoted to concluding remarks and prospects.

Notation. Upper-case bold letters are used to indicate tensors (matrices), while lowercase bold letters indicate column vectors.

Multi-scale optimisation of VAT composites

Problem description

The goal of the design strategy presented in this work is the maximisation of the first buckling load of multilayer VAT plates. The effectiveness of the MS2LOS is tested on two benchmarks: a square plate (benchmark 1, taken from [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF]) and a quarter of square plate with a circular hole (benchmark 2). Both the structures are subjected to in-plane compressive loads and their geometric parameters are illustrated in Fig. 1. The main hypotheses [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF] governing the behaviour of the VAT laminate at the macroscopic scale are listed here below.

• Each ply of the VAT laminate is made of the same material and the same thickness.

• The number of plies is kept unchanged during the optimisation process.

• The material behaviour is linear elastic.

• The VAT laminate is quasi-homogeneous (i.e. uncoupled and with the same behaviour in terms of normalised membrane and bending stiffness matrices) and fully orthotropic (both in membrane and bending) point-wise, i.e. these properties apply locally at each point of the structure;

The properties of the T300/5208 carbon-epoxy pre-preg strip, which constitute the tow of the lamina, are reported in Table 1 for the sake of completeness. Of course, the technical constants of the lamina are taken from [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF], while the polar parameters are calculated accordingly.

Technical constants Polar parameters of

[Q] a Polar parameters of [ Q] b E 1 [MPa] 181000.0 T 0 [

The multi-scale two-level optimisation strategy

In the context of the MS2LOS [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF][START_REF] Montemurro | On the effective integration of manufacturability constraints within the multi-scale methodology for designing variable angle-tow laminates[END_REF], the optimum design of the VAT laminate is articulated in two sub-problems stated at different scales. The FLP aims at determining the optimum distribution of the VAT laminate mechanical design variables satisfying the requirements of the design problem. At this level, the optimisation is performed at the macroscopic scale of the VAT plate, which is modelled as an equivalent homogeneous anisotropic continuum. The design variables are the laminate PPs [START_REF] Montemurro | An extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | Corrigendum to" an extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | The polar analysis of the Third-order Shear Deformation Theory of laminates[END_REF], varying pointwise over the structure, which describe the mechanical behavior of the VAT plate. The second level problem (SLP) is devoted to the determination of a suitable lay-up which satisfies the optimum distribution of PPs resulting from the FLP.

In this work only the FLP of the MS2LOS is faced and the related mathematical framework, presented in [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF], is extended to eigenvalue buckling problems. Inasmuch as the SLP formulation is not affected by the modifications introduced in the FLP (and the main steps of the related resolution strategy remain unchanged), this part will not be detailed in the following Sections. For more details on the SLP formulation and on the related mathematical background, the reader is addressed to [START_REF] Montemurro | Variational analysis and aerospace engineering : mathematical challenges for the aerospace of the future[END_REF][START_REF] Montemurro | On the effective integration of manufacturability constraints within the multi-scale methodology for designing variable angle-tow laminates[END_REF].

Mathematical formulation of the first-level problem

The FLP mathematical background essentially relies on the FSDT and on the B-Spline surface theory. On the one hand, the FSDT allows for integrating the influence of the transverse shear stiffness on the physical responses of the VAT laminate. On the other hand, B-Spline surfaces are used to describe the PPs distribution over the structure. The use of B-Spline surfaces to describe the PPs fields has a fundamental consequence: the anisotropy field description is completely unrelated from the mesh of the FE model. Therefore, a reduction of the number of the design variables, i.e. the PPs computed only at the B-Spline surfaces control points (CPs), can be obtained with respect to classical FE-based representation [START_REF] Setoodeh | Design of variable-stiffness laminates using lamination parameters[END_REF][START_REF] Setoodeh | Design of variable-stiffness composite panels for maximum buckling load[END_REF]. Moreover, as detailed in the following of this section, the use of B-Spline entities also allows to reduce the number of optimisation constraints to be checked during the optimisation process. These aspects lead to a strong reduction of the computational cost required to carry out the solution search. In particular, the advantages related to the use of B-Spline surfaces come from two fundamental properties: the local support property and the strong convex-hull property, see [START_REF] Piegl | The NURBS book[END_REF]. The way these properties are exploited is clarified in the following.

Without loss of generality, the following subsections focus on a meaningful design problem: the maximisation of the first buckling factor of the VAT laminate subject to feasibility constraints on the PPs and to a requirement on the total thickness.

Mechanical design variables

In the framework of the FSDT [START_REF] Reddy | Mechanics of laminated composite plates and shells. Theory and Analysis[END_REF] the constitutive law of the multilayer plate (expressed within the global frame of the laminate Γ = {O; x, y, z}) can be stated as:

r = K lam ε, (1)
where r and ε are the vectors of the generalised forces per unit length and strains of the laminate middle plane, respectively, whilst K lam is the laminate stiffness matrix (Voigt's notation). In this framework, the analytical form of these arrays is:

r =          n m q          , K lam =      A B O D O sym H      , ε =          ε 0 χ 0 γ 0          . (2)
In Eq. ( 2), A, B and D are the membrane, membrane/bending coupling and bending stiffness matrices of the laminate, while H is the out-of-plane shear stiffness matrix. n, m and q are the vectors of membrane forces, bending moments and shear forces per unit length, respectively, whilst ε 0 , χ 0 and γ 0 are the vectors of in-plane strains, curvatures and out-of-plane shear strains of the laminate middle plane, respectively, [START_REF] Reddy | Mechanics of laminated composite plates and shells. Theory and Analysis[END_REF]. The laminate homogenised stiffness matrices are introduced for optimisation purposes:

A * := 1 h A, B * := 2 h 2 B, D * := 12 h 3 D, H * := 1 h H. ( 3
)
where h is the total thickness of the plate. Consider a quasi-homogeneous, fully orthotropic laminate (i.e. a laminate with B * = 0, A * = D * and characterised by the same orthotropy type and direction for membrane and bending stiffness matrices): it can be proved that, under these hypotheses, the overall number of independent mechanical design variables describing its mechanical response reduces to only three, i.e. the anisotropic moduli R A * 0K and R A * 1 and the polar angle Φ A * 1 (which represents the orientation of the main orthotropy axis) of the homogenised membrane stiffness matrix A * [START_REF] Montemurro | An extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | Corrigendum to" an extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF]. The description of the polar framework and its application in the context of the FSDT are not reported in this work for the sake of brevity. The reader is addressed to [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF][START_REF] Montemurro | An extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | Corrigendum to" an extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF] for a detailed dissertation on the aforementioned topics.

For optimisation purposes, it is useful to introduce the dimensionless PPs:

ρ 0 := R A * 0K R 0 , ρ 1 := R A * 1 R 1 , φ 1 := Φ A * 1 π/2 , (4)
where R 0 , R 1 and Φ 1 are taken from Table 1. The expression of the laminate stiffness matrices (together with their gradient) in terms of the dimensionless PPs is given in Appendix A. In the most general case, for a VAT composite, the three independent PPs vary point-wise over the structure. As stated beforehand, the variation of the generic PP ζ is expressed by means of a B-Spline scalar function:

ζ (u 1 , u 2 ) = n 1 i 1 =0 n 2 i 2 =0 N i 1 ,p 1 (u 1 )N i 2 ,p 2 (u 2 )ζ (i 1 ,i 2 ) , ζ = ρ 0 , ρ 1 , φ 1 . (5)
The dimensionless coordinates u 1 and u 2 can be arbitrarily defined: a natural choice is to relate them to the Cartesian coordinates of the laminate global frame:

u 1 = x a x , u 2 = y a y , (6)
where a j (j = x, y) is the problem characteristic length along the j axis. In Eq. ( 5),

ζ (i 1 ,i 2 )
is the value of the dimensionless PP at the generic CP, whereas N i 1 ,p 1 (u 1 ) and N i 2 ,p 2 (u 2 ) are the B-Spline blending functions computed at the dimensionless coordinates u 1 and u 2 respectively. The parameters p 1 and p 2 are the degrees of the aforementioned blending functions, while n j + 1 represents the number of CPs along the j-th parametric direction. The set of the (n 1 + 1) × (n 2 + 1) CPs is called control net. As detailed in [START_REF] Piegl | The NURBS book[END_REF], the generic blending function N i j ,p j (u j ), with j = 1, 2, can be calculated in a recursive way. The definition of the blending functions and, thus, of the B-Spline entity, is completed by the introduction of the knot-vector v (j) T := {0, . . . , 0

p j +1 , v (j)
p j +1 , . . . , v (j)
m j -p j -1 , 1, . . . , 1

p j +1 }, j = 1, 2, (7)
where the size of the knot-vector along the j-th parametric direction is m j + 1, with m j = n j + p j + 1. It is noteworthy that the B-Spline scalar function (one for each PPs) of Eq. ( 5) constitutes the third coordinate of a generic B-Spline surface [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF][START_REF] Piegl | The NURBS book[END_REF]. In particular, for each surface, the Cartesian coordinates of every CP can be calculated by means of the Greville's abscissae formulae. For further details the reader is addresses to [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF]. For insights into the B-Spline surfaces theoretical framework the reader is addressed to [START_REF] Piegl | The NURBS book[END_REF].

The integer parameters of the B-Spline scalar function of Eq. ( 5) are set a priori and do not take part in the optimisation process. As done in [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF], also the non-trivial knot-vector components appearing in Eq. ( 7) are set as evenly distributed in the interval ]0, 1[ and are kept unchanged during the optimisation analysis. Therefore, the only design variables are the PPs defined at the CPs which can be grouped into the following vector:

x T = ρ (0,0) 0 , • • • , ρ (n 1 ,n 2 ) 0 , ρ (0,0) 1 , • • • , ρ (n 1 ,n 2 ) 1 , φ (0,0) 1 , • • • , φ (n 1 ,n 2 ) 1 . (8)
In the most general case, the total number of design variables is equal to 3 × (n 1 + 1) × (n 2 + 1). Furthermore, the feasibility constraints, which arise from the combination of the layer orientations and positions within the stack [START_REF] Vannucci | A note on the elastic and geometric bounds for composite laminates[END_REF], must be integrated into the FLP formulation. Inasmuch as the laminate is quasi-homogeneous [START_REF] Vannucci | A note on the elastic and geometric bounds for composite laminates[END_REF], it is sufficient to impose these constraints only for the PPs of matrix A * , i.e.

       -1 ≤ ρ 0 ≤ 1 , 0 ≤ ρ 1 ≤ 1 , 2ρ 2 1 -1 -ρ 0 ≤ 0 . (9)
It is noteworthy that first and second inequalities of Eq. ( 9) can be taken into account as lower and upper bounds for the dimensionless anisotropic moduli, i.e. ρ 0 and ρ 1 . Therefore, the expression of the resulting feasibility constraint on the laminate PPs, which must be considered for the generic CP, is:

g ij (x) = 2 ρ (i,j) 1 2 -1 -ρ (i,j) 0 ≤ 0 , i = 0, • • • , n 1 , j = 0, • • • , n 2 . ( 10
)
The strong convex-hull property of the B-Spline blending functions ensures that, if the feasibility constraint of Eq. ( 10) is met on the CPs, it is satisfied for each point belonging to the B-Spline surface (see [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF] for more details). The feasibility constraint of Eq. ( 10) must be imposed in order to ensure that a feasible stack, satisfying the optimal PPs fields resulting from the FLP, could be found as a result of the SLP. The reader is addressed to [START_REF] Vannucci | A note on the elastic and geometric bounds for composite laminates[END_REF] for a wide discussion on feasibility constraints in the PPs space. From Eq. ( 10), it is evident that, in the most general case, the total number of feasibility constraints is equal to (n 1 + 1) × (n 2 + 1).

Mathematical statement of the problem

The FLP aims at determining the optimum distribution of the laminate PPs maximising the first buckling factor of the structure and satisfying, simultaneously, the feasibility constraints on the laminate PPs at each point of the plate. Consider the eigenvalue buckling problem:

(K + λ k K σ ) ψ k = 0, ∀ ψ k = 0, ( 11
)
where K is the stiffness matrix of the FE model, K σ is the geometric stiffness matrix, λ k is the k-th buckling factor and ψ k is the k-th eigenvector associated to λ k . In this study, only the first buckling factor is computed, thus k = 1. In the following the first buckling factor λ 1 is denoted as λ to simplify the nomenclature. Formally, the design problem can be stated as a constrained non-linear programming problem (CNLPP) as:

min x -λ(x),
subject to:

         (K + λ k K σ ) ψ k = 0, h =const. , g ij (x) ≤ 0, i = 0, • • • , n 1 , j = 0, • • • , n 2 . (12)
The design space of the FLP, together with the type of each design variable, is detailed in Table 2. In Eq. ( 12) the overall thickness of the laminate h is kept constant.

Design variable Type

Lower bound Upper bound

ρ 0 continuous -1.0 1.0 ρ 1 continuous 0.0 1.0 φ 1 continuous -1.0 1.0
Table 2: Design space of the first-level problem.

In order to solve problem ( 12) by means of a suitable deterministic algorithm, the derivatives of both the objective function and the constraints functions with respect to the design variables, i.e. the laminate dimensionless PPs at each CP, must be computed. The derivation of the gradient of the feasibility constraints of Eq. ( 10) is straightforward:

∂g st ∂ρ (i,j) 0 = -1 if s = i, t = j, 0 otherwise, ∂g st ∂ρ (i,j) 1 = 4ρ (i,j) 1 if s = i, t = j, 0 otherwise, ∂g st ∂φ (i,j) 1 = 0. (13)
Conversely, the gradient of the k-th buckling factor requires a special attention. In order to derive its analytical expression, the local support property of the B-Spline blending functions [START_REF] Piegl | The NURBS book[END_REF] as well as the adjoint method [START_REF] Errico | What is an adjoint model?[END_REF] are exploited. To this end, consider the following proposition.

Proposition 3.1. Consider a deformable anisotropic plate subject to given external loads.

Under the hypotheses of small generalised displacements and strains, the gradient of the k-th buckling factor λ k reads:

         ∂λ k ∂ζ (i,j) = λ k w k   e∈S ij ∂ζ e ∂ζ (i,j) w * ek + λ k s T ek ∂K lame ∂ζ e ε 0e + µ T ∂K ∂ζ (i,j) d 0   , ζ = ρ 0 , ρ 1 , φ 1 , Kµ = -λ k η k . ( 14
)
where µ is the solution of the adjoint system.

The proof of proposition 3.1, the meaning of the quantities appearing in Eq. ( 14) and the pseudo-code of the algorithm used to compute the gradient of λ k are provided in Appendix B.

Remark 3.1. In Eq. ( 14), ζ e represents the generic dimensionless PP of Eq. ( 5) evaluated at the element centroid, i.e. ζ e = ζ (u 1e , u 2e ), whilst, according to Eq. ( 5), its derivative reads:

∂ζ e ∂ζ (i,j) = N i,p 1 (u 1e )N j,p 2 (u 2e ). ( 15
)
Remark 3.2. In Eq. ( 14), S ij is the discretised version of the generic CP local support defined as:

S ij := e : (u 1e , u 2e ) ∈ U (1) i , U (1)
i+p 1 +1 × U (2) j , U (2)
j+p 2 +1 . ( 16
)
The local support property of B-Spline blending functions relates the control net, which affects the shape of the B-Spline surface, to the mesh of the FE model through the definition of an influence zone for each CP, i.e. the so-called local support. Therefore the variation of the generic dimensionless PP assigned to a given CP influences only the elements belonging to the local support of such a CP. As discussed in [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF], the size of the local support zone depends on the B-Spline surface discrete parameters, i.e. the degrees of the blending functions and the number of CPs.

Remark 3.3. In Eq. ( 14), d 0 is the vector of the nodal generalised displacements, solution of the preliminary linear static analysis:

Kd 0 = f 0 , ( 17
)
where f 0 is the vector of external nodal forces (of arbitrary value, generally a unit value is used). ε 0e , (e = 1, • • • , N e ) is the vector of generalised strains of Eq. ( 2) evaluated for each element composing the FE model (the number of elements is N e ). The solution of Eq. ( 17) is used to compute the pre-stress state involved in the definition of the geometric stiffness matrix K σ . All details are provided in Appendix B.

Numerical strategy

Problem ( 12) is a non-convex CNLPP in the space of the laminate PPs. Its nonlinearity is due, on the one hand, to the nature of the objective function, the buckling factor, which is a non-convex function with respect to the orthotropy orientation field. On the other hand, the complexity of such a problem is also due to the non-linear feasibility constraints imposed on the PPs of the plate, see Eq. [START_REF] Hyer | The use of curvilinear fiber format to improve buckling resistance of composite plates with central circular holes[END_REF].

For the resolution of problem ( 12) a hybrid optimisation tool called VISION (VarIable Stiffness composItes Optimisation based on NURBS ), developed by Montemurro [START_REF] Montemurro | A contribution to the development of design strategies for the optimisation of lightweight structures[END_REF] and made by the union of the ERASMUS (EvolutionaRy Algorithm for optimiSation of Mod-Ular Systems) algorithm [START_REF] Montemurro | A contribution to the development of design strategies for the optimisation of lightweight structures[END_REF] and a general deterministic algorithm, is employed. VISION is available in both MATLAB and Python versions. In this work, the Python version has been further developed through the implementation of the objective function derivatives in Appendix B in order to deal with problem [START_REF] Jegley | Optimization of elastically tailored tow-placed plates with holes[END_REF]. The Sequential Least Square Quadratic Programming (SLSQP) algorithm available in Numpy package, from the Python-based environment Scipy, has been used as a deterministic algorithm [START_REF] Nocedal | Numerical Optimization[END_REF]. VISION is interfaced with the FE model of the VAT plate (used essentially to compute the physical/geometrical responses involved into the definition of objective function and constraint functions), see Fig. 2. The VISION tool is composed of three main phases: pre-processing phase (steps 1-5), optimisation solvers (steps 6-8) and post-processing phase (step 9). A synthetic description of each step is given here below. (DRs) and non-design regions (NDRs). Once the mesh is generated, the elements belonging to DRs are selected and their data are passed to the MATLAB/PYTHON environment. For each element, the relevant information are: its identifier, the Cartesian coordinates of its centroid as well as its area (shell element) and volume (solid element). obtain dimensionless objective and constraint functions are set here. Their definition is not unique and the algorithm allows the external user to define the reference quantities. However, pre-set strategies can be selected for the most common problems (compliance minimisation, first buckling load maximisation, etc.). 6. Responses evaluation. The system responses involved into the definition of both objective function and constraint functions can be of different nature: geometrical, physical and technological. The user can implement some of them directly into the MATLAB/PYTHON environment (e.g. geometrical and technological requirements). For those responses requiring a FE analysis to be computed (e.g. buckling load, compliance, natural frequencies, etc.) an automatic interface between the optimisation algorithms (both meta-heuristics and deterministic) and the FE codes (both commercial and in-house codes) has been implemented into the VISION tool. 7. Genetic optimisation. During the first phase, solely the ERASMUS GA is interfaced with the FE model of the VAT composite: for each individual at each generation, a FE analysis is invoked for the evaluation of the physical responses. The FE model makes use of the design variables, given by the GA and elaborated by the MAT-LAB/PYTHON code which generates the B-spline surfaces representing the dimensionless PPs fields. These fields are then projected over the FE model of the VAT plate in order to calculate the desired physical response as well as the feasibility constraint at each CP. At the end of the FE analysis, the GA elaborates the results provided by the FE model in order to execute the genetic operations. These operations are repeated until the GA meets the user-defined convergence criterion. 8. Deterministic optimisation. Due to the strong non-convex nature of the VAT design problem, the aim of the genetic calculation is to provide a potential sub-optimal point in the design space, which constitutes the initial guess for the subsequent phase, i.e. the local optimisation, where a suitable deterministic algorithm is interfaced with the same FE model of the VAT composite. In this case, the optimisation calculation is speed-up by giving explicitly the expression of the gradient of both constraint functions and objective function. The convergence of the optimiser is achieved when constraints are met (or barely met, i.e. the optimum solution is on the boundary of the search domain) and one of the following conditions occurs: (a) the predicted change of the objective function is lower than a prescribed tolerance (10 -6 ); (b) the norm of the gradient of the Lagrangian functional related to the problem at hand is very close to zero (10 -6 ); (c) the predicted change of design variables is lower than a prescribed threshold value (10 -6 ); (d) a maximum value of iterations (set by the user) has been attained. The upper bound on the iterations number is 250 for the benchmarks presented in this study. 9. Results visualisation. The optimised distribution of the laminate PPs can be converted in two different standard formats: IGS file (which can be easily imported in a CAD environment) and VTK format. The latter can be exploited and manipulated into the ParaView R environment. In the benchmarks considered in this study, the VTK format has been chosen to visualise the PPs fields.

For each benchmark, two cases have been considered.

• Case 1. Only φ 1 varies over the structure, whilst ρ 0 and ρ 1 are considered uniformly distributed. In this case, the overall number of design variables for problem ( 12) is 2 + (n 1 + 1) × (n 2 + 1). Of course, in this case there is only one feasibility constraint on the laminate PPs to be checked.

• Case 2. This represents the most general case wherein ρ 0 , ρ 1 and φ 1 vary over the structure. The overall number of design variables is 3 × (n 1 + 1) × (n 2 + 1), while the number of optimisation constraints is (n

1 + 1) × (n 2 + 1)
Inasmuch as problem ( 12) is highly non-convex, the optimised solution strongly depends upon the choice of the initial guess x 0 . In this work, the initial guess x 0 is generated in three different ways. In the first case, x 0 is null, which means that an isotropic initial guess is selected. In the the second case, x 0 is the best individual resulting from the genetic optimisation performed by means of the ERASMUS algorithm [START_REF] Montemurro | A contribution to the development of design strategies for the optimisation of lightweight structures[END_REF]. In particular, the generic individual of the GA represents a potential solution for the problem at hand. In the most general case, the genotype of the individual for problem [START_REF] Jegley | Optimization of elastically tailored tow-placed plates with holes[END_REF] is composed of (n 1 + 1) × (n 2 + 1) chromosomes with three genes coding the dimensionless PPs at each CP of the B-Spline surface. In the third case, the initial guess is randomly generated.

It is noteworthy that, since problem [START_REF] Jegley | Optimization of elastically tailored tow-placed plates with holes[END_REF] involves a requirement on the first buckling factor, the well-known mode switching phenomenon [START_REF] Costa | Eigen-frequencies and harmonic responses in topology optimisation: A cad-compatible algorithm[END_REF], typical of eigenvalue problems, can occur during the iterations of the optimisation process. For example, the mode switching phenomenon (together with the spurious mode phenomenon) affects topology optimisation [START_REF] Costa | Eigen-frequencies and harmonic responses in topology optimisation: A cad-compatible algorithm[END_REF] when considering eigenvalue problems (i.e. buckling or modal analyses). As far as topology optimisation problems involving requirements on critical buckling loads (or natural frequencies) are concerned, the most common numerical strategies available in the literature to deal with the mode switching issue have been integrated in a new topology optimisation method coupling NURBS hyper-surfaces with the classical Solid Isotropic Material with Penalisation (SIMP) method developed by Montemurro and coworkers [START_REF] Costa | Eigen-frequencies and harmonic responses in topology optimisation: A cad-compatible algorithm[END_REF][START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS hyper-surfaces for 3D topology optimization problems[END_REF]. In particular, in [START_REF] Costa | Eigen-frequencies and harmonic responses in topology optimisation: A cad-compatible algorithm[END_REF], dedicated numerical artefacts and guidelines to avoid the mode switching during the topology optimisation for problems involving modal analysis are provided.

Fortunately, in the case of VSCs the situation is not as complicated as in topology optimisation problems. In particular, only the buckling mode (eigenvector) associated to the lowest buckling factor (eigenvalue) is considered at each iteration during the optimisation. Therefore, if a mode switch occurs in between two iterations, in the next one, the correct buckling mode (i.e. that associated to the lowest eigenvalue) is considered (for the expression of the objective function and the related gradient). Moreover, a further check is introduced within the SLSQP algorithm to ensure that the lowest buckling factor of the structure continuously increases (and, thus, the merit function of problem [START_REF] Jegley | Optimization of elastically tailored tow-placed plates with holes[END_REF] monotonically decreases) during iterations. Particularly, if the mode switching occurs and the value of the merit function of problem [START_REF] Jegley | Optimization of elastically tailored tow-placed plates with holes[END_REF] at the current iteration is greater than that of the previous iteration, a new descendent direction is calculated in order to get a non-increasing value of the merit function along iterations.

Of course, a sound alternative could be to consider a different objective function involving a linear combination (i.e. a weighted sum) of the first N buckling factors, or to introduce a suitable constraint on the imposed gap between two consecutive eigenvalues. Both strategies reveal effective for topology optimisation problems as discussed in [START_REF] Costa | Eigen-frequencies and harmonic responses in topology optimisation: A cad-compatible algorithm[END_REF]. However, for the test cases presented in this work the simplified strategy discussed above is sufficient to correctly take into account for the mode switching phenomenon.

The finite element model

In order to assess the buckling factor and its gradient, an eigenvalue buckling analysis has to be carried out for the FE model of the VAT composite. The FE model is generated in ANSYS R environment by using SHELL181 elements, which have four nodes and six degrees of freedom (DOFs) per node. The FE model is capable of updating the mechanical properties of each element of the mesh through the use of the PPs distributions at each iteration of the optimisation problem. The FE models of the benchmarks analysed in this work are shown in Fig. 3. For each benchmark, the BCs are set as follows.

F x x y F x A B C D (a) y x f x f y A B C D E (b)
• Benchmark 1 is a VAT composite square plate taken from the literature [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF]. In this case, the plate is subjected to an uni-axial external compressive load F x applied on nodes located at: x = 0, y = a/2 and x = a, y = a/2; see Fig. 3a. In order to get a uniform displacement u x for sides AD and BC as done in [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF], constraint equations are imposed on nodes located at x = 0, a. The mesh size is chosen after a preliminary mesh sensitivity analysis (not reported here for the sake of brevity) on the convergence of the first buckling load, for the set of BCs reported in Table 3. It has been observed that a mesh having 12005 DOFs is sufficient to properly evaluate the first buckling load of the structure.

• Benchmark 2 represents a quarter of a square plate with a hole subject to bi-axial compressive loads; The BCs on the generalised displacements are given in Table 3, while the uniform compressive forces per unit length, i.e. f x and f y , are applied on nodes located on sides BC and CD, respectively. Also in this case the mesh size is the result of a sensitivity analysis (not reported here for the sake of brevity): an overall number of 2480 DOFs is sufficient to correctly assess the first buckling factor. It is noteworthy that a module composed of appropriate routines defining a general B-Spline surface has been coded in Python (within the VISION tool) and interfaced with the FE model of the VAT plate. The following algorithm has been implemented to compute the local stiffness properties of each element of the FE model and to determine the buckling factor and its gradient.

Algorithm 1 Buckling analysis of VAT laminates through B-Spline surfaces and polar parameters.

1: For a given set of dimensionless PPs defined at each CP, build the corresponding B-Spline surfaces. 2: Discretise the plate into N e elements. 3: Retrieve the Cartesian coordinates of the e-th element centroid, i.e. (x e , y e ) and calculate the corresponding dimensionless coordinates (u 1e , u 2e ) according to Eq. ( 6). 4: Determine the dimensionless PPs (and hence the Cartesian components of the laminate stiffness matrices) according to Eq. ( 5) and assign material properties to element e; 5: If e < N e set e = e + 1 and go to step 3, otherwise stop. 6: Execute Algorithm 2 (see Appendix B) for assessing the buckling factor and its gradient.

Numerical results

Benchmark 1: square plate subject to bi-axial compressive loads

This first benchmark refers to the geometry illustrated in Fig. 1a and it is used to assess the effectiveness of the MS2LOS by comparing the optimised solutions with those available in [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF]. It is noteworthy that, since a quasi-homogeneous fully orthotropic laminate is considered in this work, the feasible region, described by the inequality of Eq. [START_REF] Hyer | The use of curvilinear fiber format to improve buckling resistance of composite plates with central circular holes[END_REF], is smaller than that considered in [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF] The plate edge length is a = 254 mm, while its overall thickness is h = 1.524 mm. The BCs on the generalised displacements are given in Table 3. The value of the axial compressive loads applied to the VAT laminate nodes is equal to the first buckling load evaluated for a quasi-homogeneous isotropic laminate subject to the same BCs, which is equal to F x = F cr x-ISO = 2700.81 N. As done in [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF], results are reported in terms of the buckling load of the VAT plate normalised with respect to this value:

K cr = F cr x-VAT /F cr x-ISO (18)
Firstly, the results concerning Case 1 of Section 3.3 are presented (φ 1 variable over the structure and ρ 0 and ρ 1 uniformly distributed). The sensitivity of the optimised normalised buckling load K cr to the size of the control net together with the number of iterations to get convergence are reported in Table 4, where a comparison with those available in the literature [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF] is also provided. The starting point (for each number of CPs) is the isotropic solution (i.e. ρ 0 = ρ 1 = 0.0), referred as ISO in the following. For each number of CPs, the optimal value of the dimensionless anisotropic moduli are reported in Table 5, while the optimal distribution of the polar angle Φ A * 1 is illustrated in Fig. 4. The following considerations can be inferred from the analysis of these results.

• From Table 4, one can infer that the results obtained by means of the MS2LOS based on the polar formalism outperform those reported in the literature. When comparing results with the same CPs number, the increase of the normalised buckling factor varies between 12.8% (in the case of 25 CPs) and 16.3% (in the case of 81 CPs) with respect to the literature solution [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF].

• Table 5 and Fig. 4 show that, even if the solutions with 49 CPs and 81 CPs are very close in terms of dimensionless anisotropic moduli values, the refined optimal distribution of the polar angle Φ A * 1 for the case with a higher number of CPs brings to an increment of 5.86% of the normalised buckling load with respect to the case of a lower number of CPs.

As a consequence, a 7 × 7 control net has been chosed for the following optimisation analyses concerning benchmark 1. This choice represents the best compromise between computational costs and performances. Two analyses are performed in the most general case (Case 2 reported in Section 3.3) wherein ρ 0 , ρ 1 , and φ 1 vary over the structure, by considering two different starting points. In this case, the CNLPP is characterised by 147 design variables and 49 optimisation constraints. In particular, Case 2a makes use of a quasi-homogeneous isotropic solution as initial guess, while Case 2b considers, as starting point, a random-generated feasible set of PPs. It must be noticed that, for Case 2b, the random starting point has been selected as the best solution (in terms of initial value of the buckling load) among twenty different (feasible) PPs distributions generated randomly.

The normalised buckling load K cr of the optimised solutions and the number of iterations to achieve convergence are reported in Table 6, where it is compared to that of the quasi-homogeneous isotropic solution and to that of the corresponding Case 1 solution of Table 4. The optimal PPs fields related to the optimised solutions Case 2a and Case 2b are illustrated in Figs. 5 and6, respectively.

By looking at the distributions of Figs. 5 and6 and by considering the results listed in Table 6, the following remarks can be done.

• As expected, the solutions obtained in Case 2a and Case 2b are different. Indeed, the buckling factor is a non-convex function in the PPs space, thus the optimised solution depends upon the choice of the initial guess.

• Table 6 shows that both optimised solutions are better than the reference one. The solution obtained in Case 2a (when the ISO solution is used as initial guess) has a normalised buckling load 2.22 times higher than that of the ISO solution. In this case, the optimised value of K cr is identical to the one obtained in case Case 1, as well as the distributions of R A * 0K , R A * 1 and Φ A * 1 (compare results reported in Fig. 5 to those in Table 5 and Fig. 4b).

• When the initial guess is randomly generated, the algorithm converges towards a different solution characterised by a normalised buckling load 2.30 times higher than the reference value. As it can been inferred from Fig. 6, the orthotropy type continuously changes over the structure because of the variation of the dimensionless polar moduli ρ 0 and ρ 1 .

As a final remark, it must be highlighted that, although problem ( 12) is formulated in a smaller design space (when compared to that used in [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF]), more efficient results (in terms of the first buckling factor) have been found in this study. How is it possible? The reason is twofold. Firstly, LPs are not as effective as PPs, which are tensor invariants, in describing the elastic symmetries of a tensor [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF][START_REF] Montemurro | Variational analysis and aerospace engineering : mathematical challenges for the aerospace of the future[END_REF][START_REF] Montemurro | An extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | Corrigendum to" an extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | The polar analysis of the Third-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | On the effective integration of manufacturability constraints within the multi-scale methodology for designing variable angle-tow laminates[END_REF]. Therefore, when PPs are used as design variables of the FLP, the optimisation algorithm has a direct control on the elastic symmetries of tensors and it is able to search for the optimum distribution of orthotropy type (by means of the polar moduli ρ 0 and ρ 1 ) and orientation (by means of the polar angle φ 1 ) maximising the first buckling factor of the VSC. Conversely, the optimisation algorithm has no control on the orthotropy type and orientation when LPs are used as design variables of the FLP as done in [START_REF] Wu | Framework for the buckling optimization of variableangle tow composite plates[END_REF]. Accordingly, in this case the optimisation algorithm is, probably, trapped in a feasible local minimum characterised by a value of the first buckling factor lower than that characterising the local minimiser in the PPs space. Secondly, as discussed in [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF][START_REF] Montemurro | Variational analysis and aerospace engineering : mathematical challenges for the aerospace of the future[END_REF][START_REF] Montemurro | An extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | Corrigendum to" an extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | The polar analysis of the Third-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | On the effective integration of manufacturability constraints within the multi-scale methodology for designing variable angle-tow laminates[END_REF], a thin (or moderately thick) VSC (but also a standard composite reinforced with straight fibres) show an "optimal" response (in terms of strain energy, buckling factor, strength, etc.) when the membrane and bending tensors tend to exhibit the same elastic behaviour. Therefore, using a quasi-homogeneous laminate with the same orthotropic behaviour in membrane and bending (i.e. A * = D * ) helps in speeding up the convergence of the algorithm towards an efficient feasible local minimiser, i.e. a VSC configuration characterised by a very good value of the first buckling factor corresponding to an optimised membrane/bending behaviour.

Benchmark 2: square plate with a hole subject to bi-axial compressive loads

This second benchmark, taken from [START_REF] Montemurro | On the effective integration of manufacturability constraints within the multi-scale methodology for designing variable angle-tow laminates[END_REF], refers to the geometry illustrated in Fig. 1b. The values of the characteristic problem sizes are: a = 90 mm, r = 30 mm and the plate thickness is h = 3 mm. The applied BCs are those listed in Table 3. The value of the bi-axial compressive loads applied to the VAT laminate middle plane corresponds to the buckling load of a quasi-homogeneous isotropic laminate (indicated as ISO in the following), which is equal to f x = f y = f cr ISO = 118.78 Nmm -1 . For benchmark 2, a preliminary study on the convergence of the normalised buckling load K cr to the control net size is done in order to determine the number of CPs to be used during optimisation. This sensitivity analysis has been conducted by considering that only φ 1 varies over the structure, whilst ρ 0 and ρ 1 are uniformly distributed. The starting point (for each number of CPs) is the ISO solution. The results of this sensitivity analysis are given in Table 7. It has been observed that the model with 6×6 CPs represents the best compromise between accuracy and computational costs. Therefore, this choice results in 38 design variables and one feasibility constraint for Case 1 and 108 design variables and 36 feasibility constraints for Case 2.

• Case 1 -Two analyses are performed in this case, by considering two different starting points. In particular, Case 1a makes use of the ISO solution as a starting point, whilst for Case 1b the initial guess is the best solution provided by the ERASMUS algorithm at the end of the genetic search. The values of the optimised dimensionless anisotropic moduli are reported in Table 9, whereas the optimal distributions of the polar angle Φ A * 1 over the structure for Case 1a and Case 1b are illustrated in Fig. 7.

ρ 0 ρ 1
Case 1a -0.70 0.39 Case 1b

1.00 1.00 From an analysis of the results reported in Tables 8 and9 and by looking at the polar angle fields of Fig. 7, the following remarks can be inferred.

-The two optimisation calculations converge towards different local minima, showing, again, the non-convex nature of the problem at hand.

-Both the optimised solutions are better than the reference one. In particular, when the solution provided by the GA is used as initial guess, the normalised buckling load of the optimised solution is 3.23 times greater than the ISO solution and almost two times greater than that obtained by using the ISO solution as initial guess. Moreover, the optimised solution of Case 1b compared has a normalised buckling load which is about 41% greater than that of the initial guess resulting from the GA calculation.

-The solution of Case 1a is characterised by a special orthotropic behaviour with K A * = 1 because the optimised value of ρ 0 is negative. The main orthotropy axis direction is aligned with the direction of the x axis of the plate almost everywhere over the VAT laminate. Conversely, the solution of Case 1b shows a standard orthotropic behaviour characterised by a unit value of the dimensionless anisotropic moduli and by an angular field which follows the boundary of the hole.

• Case 2 -This is the most general case where all the three dimensionless PPs vary over the structure. As in the first case, three different starting points have been considered corresponding to: (a) the ISO solution ( As it can be immediately inferred from the analysis of these results, regardless of the initial guess, the optimiser converges towards the same local minimum. Accordingly, the optimal PPs fields are identical for each considered case and are illustrated in Fig. 8. In this case, the normalised buckling load of the optimised solution is increased of about 25.8% with respect to the solution found by the GA. Of course, the number of iterations to achieve convergence depends on the distance between the initial guess and the local minimum. It can be remarked that, also for benchmark 2, the point-wise variation of the anisotropic polar moduli over the structure does not necessarily imply a significant increase in the K cr factor (the percentage difference between the two best solutions provided in Tables 8 and 10 is about 1.6%). This means that the influence of the point-wise variation of the main orthotropy axis orientation (related to the polar angle Φ A * 1 ) on K cr is predominant over that of the anisotropic moduli. In fact, when looking at the optimal PPs distributions shown in Fig. 8, one can notice that the anisotropic moduli are almost uniformly distributed over the plate (apart for a variation of the R A * 1 field near the upper right vertex of the plate). The optimised VAT laminate shows a standard orthotropic behaviour in each point.

Benchmark 2: effect of the loading condition on the optimised solution

The influence of the loading conditions on the optimised solutions is evaluated here. For the sake of brevity, the optimisation analyses have been performed only for the case of uniform dimensionless anisotropic moduli and point-wise variable normalised polar angle. Three values of the load ratio (LR := f y /f x ) have been considered: 0 (i.e. only f x is applied), 0.5 and 1. For each case, the isotropic solution has been considered as a starting guess.

The normalised buckling load of the optimised solutions, the optimised anisotropic moduli and the number of iterations are listed in Table 11, for each value of LR. The effect of LR on the optimal distribution of the polar angle Φ A * 1 and on the first buckling mode shape are reported in Figs. 9 and 10, respectively. Table 11: Effect of the loading conditions on the optimised solutions.

By looking at Figs. 9 and 10 and by considering the results listed in Table 11, the following remarks can be done.

• Each optimised solution is better than the ISO solution (for a given value of LR). The solution with LR = 0 presents the biggest improvement in the normalised buckling load.

• The optimised solutions with LR = 0 and LR = 0.5 have the same optimised anisotropic moduli and a quite similar distribution of the polar angle Φ A * 1 . These solutions represent a VAT laminate characterised by a square symmetry with the main orthotropy axes oriented at ±45 deg with respect to the local value of Φ A * 1 , see [START_REF] Montemurro | An extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF][START_REF] Montemurro | Corrigendum to" an extension of the polar method to the First-order Shear Deformation Theory of laminates[END_REF].

• The LR = 1 case corresponds to the Case 1a reported in Table 9. As discussed beforehand, the optimisation algorithm stops into a local minima when the ISO solution is used as a starting guess. In this case a better result can be obtained by using the best individual provided by the ERASMUS code as a starting point.

Conclusions

In this work the FLP formulation of the MS2LOS for VAT composites has been extended in order to deal with eigenvalue buckling problems. In particular, the main contribution of this study consists of the derivation of the analytical form of the gradient of the buckling factor by taking advantage of the B-Spline surface properties, which are 11.

used to represent the laminate PPs fields of the VAT composite. The considered design problem focuses on the maximisation of the first buckling load of a VAT laminate subject to feasibility constraints and geometrical requirements. The closed-form expression of the buckling factor gradient with respect to the design variables has been analytically derived by exploiting the properties of the B-Spline blending functions. The algorithms required for the buckling factor gradient assessment have also been presented in this study: they ensure the minimum number of operations for the computation of the derivatives in order to speed-up the optimisation process while guaranteeing a high accuracy. Two test cases have been analysed to prove the effectiveness of the MS2LOS based on the polar formalism. The first benchmark problem, taken from the literature, has been considered to validate the performances of the optimised solutions resulting from the FLP of the MS2LOS through a comparison with those resulting from multi-level approaches based on LPs. The numerical results found in this work show that the optimised solutions provided by the MS2LOS based on the polar formalism outperform those presented in the literature: the improvement in the normalised buckling factor varies from 12.8% (when using a B-Spline surface with 25 CPs to describe the polar angle distribution) to 16.3% (in the case of 81 CPs). In particular, the results presented in this study are very encouraging and show that a significant increase in the buckling strength of VAT composites can be 11.

obtained with respect to a reference quasi-homogeneous isotropic solution: up to 130% for benchmark 1 and up to 229% for benchmark 2. The advantages related to the use of B-Spline surfaces and to the polar formalism in the context of the MS2LOS, presented in previous publications on the topic (i.e. a strong reduction of the computational costs when compared to FE-based approaches and the continuity of the PPs fields) are once again confirmed in this work. In addition, two important aspects of the MS2LOS based on the polar formalism should be pointed out: (a) the effectiveness of the representation of the anisotropy based on PPs which allows determining locally the optimal symmetry group; (b) the absence of simplifying hypotheses either on the nature of the stack or on the kind of orthotropy of the laminate stiffness tensors during the FLP (unlike the approaches based on the use of LPs).

Regarding the prospects of this study, some aspects of theoretical, numerical and technological nature deserve a particular attention. The optimisation strategy needs to be tested on more complex benchmarks, e.g. a representative stiffened panel extracted either from the wing or from the fuselage of an aircraft made of VAT composites. In addition, the design problem requires the formalisation of further technological constraints (e.g. gap and overlap between adjacent tows, tow width, the variation of the fibre volume fraction due to imperfections, etc.), related either to the AFP process or to the FFF+CFF technologies, in the FLP theoretical framework in order to get not only an optimised solution but also a manufacturable one. Research is ongoing on these aspects.

Taking into account for the above expressions, the derivatives of matrices A * and H * read:

∂A * ∂ρ 0 = R 0 A * 1 , ∂A * ∂ρ 1 = R 1 A * 2 , ∂A * ∂φ 1 = 2πR 0 ρ 0   -s 4 s 4 c 4 -s 4 -c 4 sym s 4   + πR 1 ρ 1   -4s 2 0 2c 2 4s 2 2c 2 sym 0   , (A.7) ∂H * ∂ρ 0 = O, ∂H * ∂ρ 1 = RH * 1 , ∂H * ∂φ 1 = -πRρ 1 -s H * 2 c H * 2 sym s H * 2 . (A.8)
Moreover, inasmuch as the quasi-homogeneity hypothesis holds, the derivatives of matrix B * are null, while those of matrix D * are equal to those of matrix A * . Finally, the derivative of the laminate stiffness matrix K lam of Eq. ( 1) with respect to the generic PP ζ reads:

∂K lam ∂ζ = diag ∂A * ∂ζ , ∂D * ∂ζ , ∂H * ∂ζ , ζ = ρ 0 , ρ 1 , φ 1 . (A.9)

Appendix B. Buckling factor gradient

The details of the proof of Proposition 3.1 are given here below.

Proof. Inasmuch as external forces do not depend upon the laminate dimensionless PPs, the derivative of the right-hand side of Eq. ( 17) is

∂f 0 ∂ζ (i,j) = 0. (B.1)
Consider, now, the following functional:

:= ψ T k (K + λ k K σ ) ψ k + µ T (Kd 0 -f 0 ) = 0, with : µ = 0. (B.2)
In the above equation, is identically null because Eqs. ( 11) and ( 17) holds. Moreover, µ is the arbitrarily defined adjoint vector. The derivative of Eq. (B.2) reads:

∂ψ T k ∂ζ (i,j) (K + λ k K σ ) ψ k + ψ T k (K + λ k K σ ) ∂ψ k ∂ζ (i,j) + +ψ T k ∂ (K + λ k K σ ) ∂ζ (i,j) ψ k + µ T ∂K ∂ζ (i,j) d 0 + K ∂d 0 ∂ζ (i,j) = 0. (B.3)
The first two terms of Eq. (B.3) are identically null due to Eq. ( 11); thus Eq. (B.3) simplifies to:

ψ T k ∂K ∂ζ (i,j) ψ k +λ k ψ T k ∂K σ ∂ζ (i,j) ψ k + ∂λ k ∂ζ (i,j) ψ T k K σ ψ k +µ T ∂K ∂ζ (i,j) d 0 + K ∂d 0 ∂ζ (i,j) = 0. (B.4)
By multiplying both sides of Eq. ( 11) to ψ T k one obtains:

ψ T k K σ ψ k = - ψ T k Kψ k λ k := - w k λ k . (B.5)
Therefore, by injecting the above formula in Eq. (B.4), the derivative of the k-th buckling factor reads:

∂λ k ∂ζ (i,j) = λ k w k ψ T k ∂K ∂ζ (i,j) ψ k + λ k ψ T k ∂K σ ∂ζ (i,j) ψ k + µ T ∂K ∂ζ (i,j) d 0 + µ T K ∂d 0 ∂ζ (i,j) . (B.6)
As discussed in [START_REF] Setoodeh | Design of variable-stiffness composite panels for maximum buckling load[END_REF], the geometric stiffness matrix of the generic element can be expressed as:

K σe = 8 i=1 r 0ei K i , (B.7)
where r 0ei is the i-th component of the vector of generalised forces per unit length of Eq. ( 2), for the generic element e, whilst K i are matrices which depend only of the geometry of the element. The expression of each K i is given in Appendix C. The vector of generalised forces per unit length r 0e is obtained as solution of the system (17), i.e.

r 0e = K lame ε 0e = K lame B e L e d 0 . (B .8)
In Eq. (B.8), K lame is the laminate stiffness matrix of Eq. ( 2), computed for the generic element, while ε 0e is the vector of the generalised strains of the laminate middle plane. B e is the product of the linear differential operator and the shape function matrix (see [START_REF] Reddy | Mechanics of laminated composite plates and shells. Theory and Analysis[END_REF]) for element e, whilst L e is the connectivity matrix of element e. Therefore, the geometric stiffness matrix of the structure K σ can be expressed as:

K σ = Ne e=1 L T e K σe L e = Ne e=1 L T e 8 i=1 r 0ei K i L e . (B .9)
Consider, now, the term ψ T k ∂K σ ∂ζ (i,j) ψ k in Eq. (B.6). By taking into account for Eq. (B.9), it reads: ∂ζ (i,j) . (B.12) By injecting Eq. (B.8) into the above formula and by taking into account for the local support property of Eq. ( 16), one obtains:

ψ T k ∂K σ ∂ζ (i,j) ψ k = Ne e=1 ψ T k L T e 8 i=1 ∂r 0ei ∂ζ (i,j) K i L e ψ k . (
ψ T k ∂K σ ∂ζ (i,j) ψ k = e∈S ij
∂ζ e ∂ζ (i,j) s T ek ∂K lame ∂ζ e ε 0e + η T k ∂d 0 ∂ζ (i,j) , (B. [START_REF] Olmedo | Buckling response of laminates with spatially varying fiber orientations[END_REF] where ∂ζ e ∂ζ (i,j) represents the partial derivative of the generic dimensionless PP expressed through a B-Spline scalar function, given in Eq. ( 15), whilst ∂K lame ∂ζ e is the partial derivative of the laminate stiffness matrix with respect to the dimensionless PPs evaluated at the element centroid, whose expression is reported in Appendix A. In Eq. (B.13), the vector η k is defined as: Consider, now, the term ψ T k ∂K ∂ζ (i,j) ψ k in Eq. (B.6). As discussed in [START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF], the partial derivative of the structure stiffness matrix reads: (B.18)

The adjoint vector µ can be chosen in such a way that the term multiplying ∂d 0 ∂ζ (i,j) vanishes, i.e.

Kµ = -λ k η k . (B.19)
Finally, Eq. ( 14) can be easily got by injecting Eq. (B.19) into Eq. (B.18). This last passage concludes the proof.

It is noteworthy that the assessment of the k-th buckling factor gradient requires the resolution of three systems, i.e. Eq. ( 11), Eq. ( 17) and the adjoint system of Eq. (B. [START_REF] Raju | Prebuckling and buckling analysis of variable angle tow plates with general boundary conditions[END_REF]). Therefore, for each iteration of the optimisation process, the following algorithm is invoked to carry out all the necessary steps for computing the gradient of λ k .

Algorithm 2 Computation of the gradient of the k-th buckling factor.

1: Solve Eq. ( 17) and get d 0 , ε 0e and r 0e , ∀e. 2: Build K σ according to Eq. (B.9). 3: Solve Eq. ( 11) and get λ k , ψ k , ε ek , ∀e and ∀k. (i,j) from the first of Eq. ( 14).

Algorithm 3 Derivation of matrices K i .

1: Build a FE model made of a single element.

2: Set arbitrary material properties for the element. The material properties should be conveniently set in order to obtain a diagonal laminate stiffness matrix K lam (in order to avoid coupling effect). 3: Impose an elementary strain field (ε 0ei = 0, ε 0ej = 0, j = 1, • • • , 8 and j = i) by using suitable BCs at the four nodes. 4: Solve Eq. ( 17) and activate the pre-stress option (in this way the commercial FE code builds K σe according to the usual definition [START_REF] Ansys | ANSYS Mechanical APDL Basic Analysis Guide[END_REF]). 5: Get r 0ei and K σe from the FE software.

6: Calculate K i = K σe r 0ei 7: If i < 8 set i = i + 1 and go to step 3, otherwise stop.

The expressions of K i for a square SHELL181 element of side L are provided here below. Each matrix K i is a symmetric partitioned matrix, composed of symmetric blocks, containing several null components. Therefore, only non-null terms are provided in the following:

K 1 = 1 8        Â1 B1 BT 1 Ĉ1 -Â1 -B1
-BT Of course, the analytical form of matrices K i can be determined in the most general case of rectangular (or pseudo-rectangular) elements. These expressions are not reported here for the sake of brevity.

K 2 = 1 8        Â2 B2 BT 2 Ĉ2 -Â2 -B2 -BT 2 -Ĉ2 -Â2 -B2 -BT 2 -Ĉ2 Â2
K 7 = 1 72            - 24

Figure 1 :

 1 Figure 1: Geometric parameters of (a) the square plate and (b) the square plate with a circular hole.

Figure 2 :

 2 Figure 2: VISION flowchart[32].

1 .

 1 B-Spline surfaces parametrisation. During this phase the user must set the discrete parameters (number of CPs, degrees of the blending functions) tuning the shape of the B-Spline surfaces representing the laminate PPs fields as well as the values of the non-trivial components of each knot-vector. The above quantities are just geometric parameters of the B-Spline blending functions which are not included in the vector of design variables. If the user does not provide the knot-vector components, they are uniformly distributed in the interval [0, 1]. 2. FE model information. The FE model of the problem at hand must be properly prepared in the external FE code, in terms of both geometry and mesh of design regions

Figure 3 :

 3 Figure 3: FE model and related BCs for (a) the square plate and (b) the square plate with a hole subject to bi-axial compression loads.

  Uy = Uz = θy = θz = 0 AB Uy = θx = θz = 0 CD Ux = Uz = 0 BC, DA Uy = Uz = θx = θz = 0 BC Uy = Uz = 0 Ux(0, y) = Ux(0, a/2), Ux(a, y) = Ux(a, a/2) ∀y ∈ [0, a 2 [∪] a 2 , a] DE Ux = θy = θz = 0

ρ 0 ρ 1 VAT 5 × 5 -

 155 Case 1 0.73 0.93 VAT 7 × 7 -Case 1 0.81 0.95 VAT 9 × 9 -Case 1 0.80 0.95

Figure 4 :

 4 Figure 4: Benchmark 1 -Case 1 : optimal distribution of the polar angle Φ A * 1 .

1 Figure 5 : Benchmark 1 -

 151 Figure 5: Benchmark 1 -Case 2a: optimal distribution of the PPs over the structure.

1 Figure 6 : Benchmark 1 -

 161 Figure 6: Benchmark 1 -Case 2b: optimal distribution of the PPs over the structure.

Figure 7 :

 7 Figure 7: Benchmark 2: optimal distribution of the polar angle Φ A * 1 over the structure for Case 1a (a) and Case 1b (b).

1 Figure 8 :

 18 Figure 8: Benchmark 2 -Case 2: optimal distribution of the PPs over the structure.

Figure 9 :

 9 Figure 9: Optimal distribution of the polar angle Φ A * 1 over the structure for the optimised solution listed in Table11.

Figure 10 :

 10 Figure 10: Mode shape related to the first buckling load for the optimised solutions listed in Table11.

B. 10 )

 10 By introducing the following definitions T ek := {ψ T ek K 1 ψ ek , • • • , ψ T ek K 8 ψ ek }, with : ψ ek = L e ψ k , (B.11)Eq. (B.10) simplifies to:ψ T k ∂K σ ∂ζ (i,j) ψ k = Ne e=1 s T ek ∂r 0e

  ek K lame B e L e . (B.14)

  e B e dS L e , (B.15)where A e represents the area of the generic shell element. By injecting the above formula into the product ψ T k ∂K ∂ζ (i,j) ψ k , one obtains:ψ T k ∂K ∂ζ (i,j) ψ k = e∈S ij∂ζ e ∂ζ (i,j) w * ek , expression, ε ek represents the vector of the generalised strains of element e for the k-th mode shape. By injecting Eqs. (B.13) and (B.16) in Eq. (B.6), one obtains ∂λ k ∂ζ (i,j) (i,j) w * ek + λ k s T ek ∂K lame ∂ζ e ε 0e + µ T ∂K ∂ζ (i,j) d 0 + µ T K + λ k η T

K 4 =

 4 K 5 = K 6 = O, (C.4)

7 = L 2 , 8 , 8 = L 2 ,

 72882 Â8 = ÂT 8 , B8 = -BT 8 , Ĉ8 = ĈT 8

  Out-of-plane shear stiffness matrix of the pre-preg strip.

		MPa]	26878.6659 T [MPa]	5393.6268
	E 2 [MPa]	10273.0 T 1 [MPa]	24738.3141 R [MPa]	1776.3732
	G 12 [MPa]	7170.5 R 0 [MPa]	19708.6659 Φ [deg]	90.0
	ν 12	0.28 R 1 [MPa]	21436.2608
	ν 23	0.42 Φ 0 [deg]	0.0
		Φ 1 [deg]	0.0
	Thickness
	h ply [mm] 0.127
	a In-plane reduced stiffness matrix of the pre-preg strip.

b

Table 1 :

 1 Material properties of the T300/5208 carbon-epoxy pre-preg strip.

  [START_REF] Dirk | The engineering aspects of automated prepreg layup: History, present and future[END_REF]. Identification of DRs and NDRs. Not all the zones of the structure have to be optimised and this requirement can be fulfilled by forecasting proper NDRs within the FE model at specific locations. Therefore, including the PPs defined at all CPs (even those CPs whose local support falls within NDRs) among the design variables is useless. To this purpose, before launching the optimisation, a dedicated function checks all CPs local supports. Consider the local support of the generic CP S ij : if the dimensionless coordinates (u 1e , u 2e ) of at least one element centroid belonging to the DR are in S ij , then the corresponding PPs defined at that CP are inserted in the design variables array; otherwise, they are deleted from this array. The discarded PPs are set equal to their lamina counterparts for those CPs whose local support is the empty set. 4. Symmetries Application. Symmetries can be interpreted as variables saving from a computational viewpoint. Only independent CPs are effectively collected in the design variables array. 5. Reference Quantities. Geometrical, physical and manufacturing responses used to

Table 3 :

 3 BCs on generalised displacements for benchmarks 1 and 2.

Table 4 :

 4

		Kcr	Iterations
	ISO	1.00	-
	VAT 5 × 5, [25]	1.88	-
	VAT 7 × 7, [25]	1.96	-
	VAT 9 × 9, [25]	2.02	-
	VAT 5 × 5 -Case 1	2.12 (+12.8%)	158
	VAT 7 × 7 -Case 1	2.22 (+13.3%)	144
	VAT 9 × 9 -Case 1	2.35 (+16.3%)	118

Benchmark 1 -Case 1 : sensitivity of the optimum normalised buckling load to the number of CPs (the percentages in parentheses indicate the improvements with respect to corresponding the literature solutions).

Table 5 :

 5 Benchmark 1 -Case 1 : dimensionless anisotropic moduli for the optimised solutions.

Table 6 :

 6 Benchmark 1 -Case 2 : optimised solutions.

		Kcr	Iterations
	ISO	1.00	-
	VAT 7 × 7 -Case 1	2.22	144
	VAT 7 × 7 -Case 2a	2.22	192
	VAT 7 × 7 -Case 2b	2.30	166

Table 7 :

 7 Benchmark 2 -Sensitivity of the normalised buckling load to the number of CPs.

		Kcr	Iterations
	ISO	1.00	-
	VAT 4 × 4 1.66	28
	VAT 5 × 5 1.67	46
	VAT 6 × 6 1.79	37
	VAT 7 × 7 1.79	83
	VAT 8 × 8 1.80	90

Table 8 :

 8 Table 8 lists the normalised buckling load of the ISO solution and of the ERASMUS initial guess in the optimised configurations together with the number of iterations to achieve convergence.

		Kcr	Iterations
	ISO	1.00	-
	ERASMUS -GA 2.29	-
	Case 1a	1.79	37
	Case 1b	3.23	57

Benchmark 2 -Case 1 : optimised solutions.

Table 9 :

 9 Benchmark 2 -Case 1: dimensionless anisotropic moduli for the optimised solution.

Table 10 :

 10 Case 2a); (b) the best solution provided by the ERASMUS algorithm (Case 2b); (c) a feasible initial guess randomly generated (Case2c). The normalised buckling load relative to all initial guess and the optimised solutions are reported in Table10. Benchmark 2 -Case 2: optimised solutions.

		Kcr	Iterations
	ISO	1.00	-
	RANDOM	1.02	-
	ERASMUS-GA 2.44	-
	Case 2a	3.29	232
	Case 2b	3.29	99
	Case 2c	3.29	86

4 :

 4 Calculate w k and w * ek , ∀e and ∀k, according to Eqs. (B.5) and (B.17), respectively. 5: Calculate s ek and η k , ∀e and ∀k, from Eqs. (B.11) and (B.14), respectively. 6: Solve the adjoint system of Eq. (B.19); get µ. Calculate ∂K ∂ζ (i,j) from Eq. (B.15).

	8: Calculate	∂λ k ∂ζ

7:

  K7 -12 Â7 + 4 Ĉ7 24 K7 + 12 B7 + 2 Ĉ7 24 K7 -12 B7 + 2 Ĉ7-24 K7 + 12 Â7 + 4 Ĉ7

						12 K7 + 6 B7 + Ĉ7	-12 K7 -6 Â7 + 2 Ĉ7		
		12 K7 -6 B7 + Ĉ7	-12 K7 + 6 Â7 + 2 Ĉ7	 	 	-12 K7 + 6 Â7 + 2 Ĉ7 -24 K7 + 12 Â7 + 4 Ĉ7	12 K7 -6 B7 + Ĉ7 24 K7 -12 B7 + 2 Ĉ7	 	      	,
		-12 K7 -6 Â7 + 2 Ĉ7	12 K7 + 6 B7 + Ĉ7			24 K7 + 12 B7 + 2 Ĉ7	-24 K7 -12 Â7 + 4 Ĉ7		
	with, K7 = KT

[START_REF] Setoodeh | Design of variable-stiffness laminates using lamination parameters[END_REF]

, Â7 = ÂT 7 , B7 = -BT 7 , Ĉ7 = ĈT 7 , and K

[START_REF] Wang | 3d printing of polymer matrix composites: A review and prospective[END_REF][START_REF] Dirk | The engineering aspects of automated prepreg layup: History, present and future[END_REF]

Acknowledgements G. A. Fiordilino is grateful to the Nouvelle-Aquitaine region for its contribution to this paper through the SMARTCOMPOSITE project. This paper presents also part of the activities carried out within the research project PARSIFAL (Prandtlplane ARchitecture for the Sustainable Improvement of Future AirpLanes), which has been funded by the European Union under the Horizon 2020 Research and Innovation Program (Grant Agreement n.723149).

We understand that the Corresponding Author is the sole contact for the Editorial process.

He/she is responsible for communicating with the other authors about progress, submissions of revisions and final approval of proofs Signed by all authors as follows:

Giacinto Alberto Fiordilino Michele Iacopo Izzi Marco Montemurro (corresponding author)

Data availability

The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Appendix A. Analytical expression of the laminate stiffness matrices and of their gradient

Under the hypothesis of quasi-homogeneous laminate, i.e. B * = 0 and A * = D * , the expression of the homogenised membrane stiffness matrix in terms of the dimensionless PPs reads:

where matrices A * 0 , A * 1 , A * 2 are defined as

with c 2 = cos (πφ 1 ) , s 2 = sin (πφ 1 ) , c 4 = cos (2πφ 1 ) , s 4 = sin (2πφ 1 ) . (A.3)

Similarly, matrix H * can be decomposed as:

where matrices H * 0 and H * 1 are defined as:

It is noteworthy that the quantities T 0 , T 1 , T , R, Φ 1 and Φ, appearing in Eqs. (A.1)-(A.6), are the PPs of the pre-preg tow listed in Table 1.

Appendix C. Expression of the element geometric stiffness matrix

The analytical expression of matrices K i , (i = 1, • • • , 8), appearing in Eq. (B.7), is reported here below. In particular, the definition of the element geometric stiffness matrix, given in Eq. (B.7), is not the classical one used in commercial FE software, like ANSYS R [START_REF] Ansys | ANSYS Mechanical APDL Basic Analysis Guide[END_REF]. However, unlike the classical definition of the element geometric stiffness matrix, Eq. (B.7) is really useful for determining a straightforward expression of the buckling factor derivative.

In the following, the algorithm for retrieving the expression of each matrix K i ∈ R 24×24 for a shell element with four nodes and six DOFs per node (like the SHELL181 ANSYS R shell element), whose kinematics is described in the framework of the FSDT, is presented. Of course, this algorithm must be executed off-line, i.e. before the optimisation process, once the element type has been selected.

to feasibility and geometrical constraints is addressed in this work

• The problem is formulated and solved in the framework of the multi-scale two-level optimisation strategy (MS2LOS) • The approach is based on the polar method to describe the macroscopic response of the VSC • The polar parameters (PPs) of the VSC are described by means of B-Spline surfaces in the framework of the first-order shear deformation theory • The analytical expression of the buckling factor gradient is derived in order to speed-up the optimisation process

Author Agreement Statement

We the undersigned declare that this manuscript is original, has not been published before and is not currently being considered for publication elsewhere.

We confirm that the manuscript has been read and approved by all named authors and that there are no other persons who satisfied the criteria for authorship but are not listed. We further confirm that the order of authors listed in the manuscript has been approved by all of us.

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: