Material Design for Optimal Postbuckling Behaviour of Composite Shells

Lightweight thin-walled structures are crucial for many engineering applications. Advanced manufacturing methods are enabling the realization of composite materials with spatially varying material properties. Variable angle tow fibre composites are a representative example, but also nanocomposites are opening new interesting possibilities. Taking advantage of these tunable materials requires the development of computational design methods. The failure of such structures is often dominated by buckling and can be very sensitive to material configuration and geometrical imperfections. This work is a review of the recent computational developments concerning the optimisation of the response of composite thin-walled structures prone to buckling, showing how baseline products with unstable behaviour can be transformed in stable ones operating safely in the post-buckling range. Four main aspects are discussed: mechanical and discrete models for composite shells, material parametrization and objective function definition, solution methods for tracing the load-displacement path and assessing the imperfection sensitivity, structural optimisation algorithms. A numerical example of optimal material design for a curved panel is also illustrated.


Introduction
Thin-walled composite structures are employed in a wide range of structural applications, particularly in the aerospace industry due to the high strength-to-weight ratio. Their design is dominated by buckling, which is mainly influenced by geometry and material properties. The geometry is usually constrained by the structural functionality and only little changes are possible. Conversely, the spatial distribution of the elastic properties (like for the fibre orientation) can be easily varied. Then, an optimization process of the material distribution can provide the desired mechanical behavior in terms of displacement and capacity. Different manufacturing options are also available to tailor the stiffness and reduce buckling effects: grid stiffeners [1], multi-layered and variable thickness composites [2], variable angle tows (VATs) [3]. New technologies for 3D printed products are also opening new prospective for cheaper and more tunable structures [4][5][6].
Moreover, as known, optimal properties in composites are usually sought also by controlling the orientation of the fibers in each layer, since the fibers orientation significantly affects the stiffness distribution, hence, the load-carrying capability or the elastic limits states (see, for instance, [7,8]). A promising direction in the context of material optimization is that offered by nanostructured materials which can exhibit multifunctional properties and are thus prone to more advanced multi-objective optimizations. In this field, nanocomposite materials made of thermosetting or thermoplastic polymers integrated with carbon nanotubes (CNTs) are currently subject to intense developments [9][10][11]. The mechanical behavior of shells can be simulated by different structural models. The Mindlin-Reissner model is the most common for shear flexible shells, where the kinematics is governed by displacements and rotations of the middle surface. The Kirchhoff-Love one is attractive for thin shells, where only displacements of the mid-surface are needed. Alternative shear flexible models have been proposed, like solid-shell elements, which are solid elements able to obtain the shell solution without meshing through the thickness. Such formulations use displacement degrees of freedom (DOFs) only and the number of overall DOFs can be equal to the one in Mindlin-Reissner elements [12,13], but without the rotation parametrization in large deformation problems. Beside finite element formulations [12], Isogeometric analysis (IGA) based on NURBS highly continuous shape functions [14] is an interesting alternative for the description of geometry and kinematics over the shell mid-surface.
The linearized buckling load is the objective function in many optimal design strategies proposed in the literature. However, this can lead to an elastic limit state known as buckling mode interaction, characterized by unstable post-critical behavior [15] with a deterioration of their capacity due to geometrical, load and material imperfections. Instead, a more reliable design, which takes into account the geometrically nonlinear behaviour, should be considered. In this case, the failure load of the structure can be used as objective function to maximize. This can be defined as the first limit load for the unstable load-displacement curves or as the load leading to deformation limits, taking into account the typical postbuckling stiffness reduction. Optimizing the post-buckling behaviour is a challenging task. Firstly, a suitable mechanical model and its discrete approximation are required to describe with acceptable accuracy geometry, boundary conditions and deformations. Discretization techniques are usually needed and, then, the structural problem is generally described by a high number of discrete nonlinear equations defining the equilibrium load-displacement path. Moreover, an imperfection sensitivity is generally needed for reliable estimates of the safety factor.
The static Riks method [16,17] is a standard tool for path following the solutions of a set of nonlinear equations. This approach is suitable for assigned data, but not for structural optimizations, which require a new equilibrium path for any change in the design variables, since the single run is too time-consuming with current CPUs. The same holds for an imperfection sensitivity analysis. Promising generalizations of the Riks method have been presented in [18][19][20], which are able to perform parametric analysis in a more efficient way by using the fold line concept.
In the optimal design presented in [21], the failure load is given by a nonlinear finite element (FE) buckling problem, extended in [22,23] in order to consider the worst case of the geometrical imperfection. An alternative reduced order model formulation is offered by numerical implementations of Koiter's theory of elastic stability [24], allowing to estimate the initial post-critical behaviour in terms of slope and curvature of the bifurcated branches [25].
More recently, a solution algorithm based on Koiter's theory implemented within a Finite Element environment was proposed in [26,27]. In this case, the design is able to consider general geometries, loading and boundary conditions. Moreover, a good accuracy in predicting the initial postbuckling response is given by a multi-modal asymptotic expansion which accounts also for nonlinear buckling modal interactions [28]. The strategy also provides an inexpensive sensitivity analysis with a statistical estimation of the worstcase imperfection, assumed to be a combination of the linearized buckling modes of the perfect structure. A hybrid solution strategy, referred to as the Koiter-Newton approach, was further investigated in [29,30].
Despite the difficulties associated with the prediction of the nonlinear behavior, another challenge is the solution of the optimization problem that is generally nonlinear and nonconvex. Its solution is usually computationally expensive and difficult due to the possibility of local minima. Frequently employed algorithms are the random search  [26,31], genetic algorithms [32] and gradient-based techniques such as the method of moving asymptotes [33] or sequential linear programming [34]. This paper is a review of the recent computational developments concerning the material design for optimising the post-buckling response of composite thin-walled structures, showing how baseline products with unstable behaviour can be transformed in stable ones operating safely in the post-buckling range. Four main aspects are discussed: mechanical and discrete models for composite shells, material parametrization and objective function definition, solution methods for tracing the load-displacement path and assessing the imperfection sensitivity, structural optimisation algorithms. A numerical example of optimal material design is also given, where the optimisation suppresses the snap-through instability leading a globally stable behaviour. Examples of applicability for actual engineering uses can be found in [27,35].
The paper is organized as follows: Section 2 describes the solid-shell continuum and its discrete finite element and isogeometric counterpart for elastic shell structures; Section 3 reviews the material parametrizations and introduces the optimization problem; Section 4 discusses the solution method for tracing the equilibrium path of a slender elastic structure with particular focus on the reduced order modelling provided by the Koiter method; Section 5 reviews the nonlinear optimization algorithms suitable for the problem under consideration; a numerical example of optimization for a curved panel is illustrated in Section 6; finally, conclusions are drawn in Section 7.

Mechanical and Numerical Models for Composite Shells
The main equations of the shell model are now described. The outset is the 3D Cauchy continuum whose deformation is described by means of the Green-Lagrange strain. A linear approximation is assumed through the shell thickness for the kinematics according to the solid-shell formulation [12,25,36]. The model is rotation-free, so making simple and effective its discrete formulation. Two kinds of discretization over the shell mid-surface are presented, namely linear finite elements and isogeometric (IGA) formulations.

Constitutive Matrix for the Lamina
Deflections and buckling behavior of multi-layered composites can be modeled efficiently using a homogenized material model based on the hypothesis of the classical lamination theory.
The constitutive matrix of each lamina is usually known with respect to a local Cartesian reference system, defined by the orthogonal triad {ē 1 ,ē 2 ,ē 3 }, withē 1 the fibre direction. Then, it is necessary to express it with respect to the Cartesian reference system of the homogenised material, defined by the triad {e 1 , e 2 , e 3 }. Let assume thatē 3 ≡ e 3 and denote with ϑ the angle around e 3 betweenē 1 and e 1 . To simplify the notation, we omit to report explicitly the dependence of all quantities from the ith lamina. Moreover, it is worth noting that ϑ can vary over the mid-surface of the ply in VAT composites.
Strain vector¯ is linked to that in the global reference system bȳ where Materials 2021, 14, 1665 5 of 20 0 denotes zero matrices of suitable dimensions and The lamina elastic law linking the Green-Lagrange strain to the second Piola-Kirchoff stress is then where by standard transformation The block matrixC of the orthotropic constitutive law of the lamina with respect to the lamina reference system isC whereC p is obtained assuming a plane stress condition with a decoupling of membrane and thickness strains in order to eliminate thickness locking [12]. The coefficientC 33 , linking thickness stress and strain, is maintained in order to avoid zero energy modes (thickness stretch).

The Strain Energy for the Shell Model in Generalised Quantities
Letting Ω be the mid-surface of the shell, the stored strain energy of our shell model can be conveniently rewritten, in compact notation, as where The transverse shear stiffness C t can be evaluated more accurately by means of shear correction factors as, for example, reported in the Abaqus/Standard [37] manual. Finally note that, as recently proposed in [38], thermal effects due to general temperature distributions can be easily accounted for in the model. Generally, a membrane-flexural coupling is possible. Higher order lamination theories or layer-wise interpolations [39,40] are also available for obtaining more accurate inter-laminar stresses. A recent paper [41] proposes the inter-laminar stress recovery starting from the homogenized response.

Discretization Methods
The continuum solid-shell model can be discretized using a displacement based formulation or a mixed one. In this section, a displacement based finite element [12] and isogeometric [14] formulation are reviewed, while we refer to [12,25] for a mixed one.

Geometry and Displacement Interpolation
The shell is discretized in quadrilateral elements through a mesh generation. According to the isoparametric concept, the same interpolation is used for the geometry and displacements. The geometry is described over the element as where x e = [x 0e , x ne ] collects the element discrete variables of the geometry corresponding to x 0 and x n , respectively. The matrix N u (ζ) collects the interpolation functions where N(ζ 1 , ζ 2 ) are bi-dimensional functions of the mid-surface coordinates only as proposed in [13,36]. The displacement field is interpolated using the same shape functions where d e = [d 0e , d ne ] collects the element discrete degrees of freedom (DOFs) for the displacement fields u 0 and u n . The Green-Lagrange strains in Equation (5), upon considering Equations (13) and (15), become where L(ζ 1 , ζ 2 ) := Q(ζ 1 , ζ 2 , x e ) and Q(ζ 1 , ζ 2 , d e ) has a linear dependence from d e , and its expression can be found in [13]. The stored energy of the shell can be evaluated using a numerical integration as where d is the global counterpart of d e , e denotes the element, g indicates the integration point, and w g is the corresponding weight.

Finite Element Formulation
If bilinear shape functions N(ζ 1 , ζ 2 ) are employed for the middle surface approximation, we have a hexahedron solid-shell linear element. Low order elements are however affected by shear and trapezoidal locking. In order to eliminate these undesired inaccuracies, it is possible to redefine the transverse shear strain componentsĒ ηζ ,Ē ξζ and the transverse normal strain componentĒ ζζ by the Assumed Natural Strain (ANS) technique with number and location of the sampling points as reported in [12,42]. The in-plane bending response of the element is improved by replacing the in-plane shear strainĒ ξη with its value at ξ = η = 0 that is a Selective Reduced Integration (SRI) retaining the correct matrix rank.

Isogeometric Formulation
In the solid-shell isogeometric version, NURBS of arbitrary order and continuity are employed as middle surface shape functions N(ζ 1 , ζ 2 ). Locking occurs for low order interpolations. Shear and membrane lockings are typical in small deformation problems. Furthermore, an additional locking occurs in large deformations, even for flat plates, due to the nonlinear strain measure in Equation (16). The high continuity of the NURBS functions allows the use of patch-wise numerical integrations [43,44], based on a lower number of integration points compared to Gauss rules. Moreover, well-tuned patch-wise reduced scheme can avoid locking. We refer to [13,36,45] for more details on this topic. The displacement-based IGA model represents a reliable choice from the point of view of the discrete approximation and the efficiency of the integration compared to high order FEs.
The NURBS high continuity also allows the use of a Kirchhoff-Love model for thin shells as presented in [45,46], which has the advantage of describing the kinematic using only the mid-surface displacement, thus halving the number of unknowns.

The Objective Function
The optimization process is aimed at maximizing the collapse load of the composite shells. In buckling problems, the collapse load can be defined as the lower bound between the critical limit load λ lim and the load associated with a deformation limit λ de f . With α being the vector collecting the generic design optimization parameters, the objective function can thus be written as The evaluations λ lim and λ de f , and thus the objective function computation, require the construction of the equilibrium path of the structure for assigned design variables α.
Another possibility is to change the optimization problem in minimizing the displacements occurring at an assigned load level. This is particularly suitable for structure with an a priori known stable behavior [35].

Constraints
Some constraints should be considered in the optimization process (see [27]). Firstly, manufacturing constraints are needed to guarantee the applicability for actual engineering products. For example, the manufacturing of VAT composites requires a limit on the minimum steering radius of the fibers. Moreover, the material behavior of the composite can be assumed to be elastic up to failure. In this case, delamination and damage can be prevented by adding constraints on maximum values of the strains or stresses.

Design Variables
The design variables collected in vector α define how the material properties, and then the constitutive matrix in Equation (12), depend on the material properties at this point, that is, we now have C(α, ζ 1 , ζ 2 ).

Layer-Wise Parameters
The most simple approach consists of optimising the lamination using, as optimisation variables, the angles ϑ i of the stacking sequence, i.e., α = {α 1 , · · · α n } with each α i ≡ ϑ i constant over the patch. This means that for single shells we have a number of optimization variables at most equal to the number of layers.
For VAT composites, the only change is that the lamina orientation is controlled by more parameters, two according to the description proposed in [47]: We refer to [27,47] for more details.

Lamination Parameters
The previously discussed layer approach requires a large number of optimization parameters when a large number of layers is used. For this reason, other approaches have been proposed in literature in order to parametrise the homogenised constitutive matrix, without an explicit use of the stacking sequence. To this aim, the main approaches are polar decompositions [48] and lamination parameters [49].
According to the last one, the constitutive matrix can be described as a linear function of 12 lamination parameters [35,50]. Assuming to have the same material through the shell thickness, the matrices C ee , C eχ , C χχ and C t of the stiffness matrix in Equation (12), are parametrized in terms of lamination parameters using the material invariants Γ k : Lamination parameters are defined in [−1, 1] as where f i is the component of the vector f = [cos(2ϑ), sin(2ϑ), cos(4ϑ), sin(4ϑ)] and ϑ is the angle around e 3 of the VAT at a given point. Matrices Γ i and Γ s i can be evaluated as with U k reported in [50]. The lamination parameters, controlling the elastic matrix, are interpolated over the shell domain as ξ = N ξ (ζ 1 , ζ 2 )α (23) where α collects the lamination parameters ξ j i , i = 1 . . . 4 , j = A, B, D at the control points of a further grid used for the material description.
Following [51], the reconstruction of the stacking sequences over the surface in terms of the lamination parameters, once they are evaluated by the optimal design process, is obtained in a second stage by minimising a least-square distance between the target distribution (23) and the lamination parameters related to the unknown fibre angle distribution. To this aim, the vector ϑ collecting the orientations of all layers can be interpolated as where N ϑ [ζ 1 , ζ 2 ] are bi-dimensional shape functions, while ϑ d are the corresponding discrete variables. The least-square problem can be written as where C(ϑ d ) ≤ 0 is a set of manufacturing constraints that the fibre tow must satisfy, i = 1 . . . n p are a fine set of sample points over the shell surface and l p is the number of lamination parameters at each point. The solution of this non-convex and nonlinear Materials 2021, 14, 1665 9 of 20 problem can be obtained using a multi-start GCMMA. A complex-step method [52] is advised for the gradient evaluations.

Equilibrium Path Evaluation
The evaluation λ lim and λ de f , and thus the objective function computation, requires the construction of the equilibrium path of the structure for assigned design variables α. A common approach for path following the equilibrium curve is the Riks arc-length method [13,16,20,53]. In this case, the nonlinear equations in the kinematic unknowns are solved step-by-step using the Newton-Raphson method. However, this kind of analysis bears a significant computational cost due to the large size of the matrices associated with a high number of DOFs. Furthermore, a reliable evaluation of the equilibrium path should take into account the sensitivity of the structure to imperfections in order to detect the worst imperfection scenario [21]. For this reason, an alternative approach called Koiter's method [25] was proposed, which assembles a reduced order model based on Koiter's theory of elastic stability for the assigned material configuration. The corresponding reduced nonlinear equations, usually in a lower number of unknowns, are then solved to obtain a good estimate of the equilibrium path at a low computational cost. The most interesting feature of the method is the possibility of including imperfections a posteriori in the reduced system [28] of the perfect structure, thus enabling inexpensive sensitivity analyses.

Path-Following Analysis
The system of discrete equilibrium equations is then obtained through enforcement of the stationarity of the total potential energy according to where r is the residual vector, s(u) is the vector of generalized stress resultants (i.e., restoring forces), f is the load vector per unit multiplier, u collects the discrete variables and λ is the load multiplier. Note that u collects the global displacements d in a displacement formulation, while it can contain other variables like stresses and strains in mixed formulations. The solutions of Equation (26) define the equilibrium paths of the structure in the u − λ space. The Riks approach [16] is the most popular strategy for solving Equation (26) by adding a constraint of the shape g(λ; u) − s = 0, which defines a surface in R N+1 . Assigning successive values to the control parameter s = s (k) , the solution of the nonlinear system defines a sequence of points (steps) z (k) ≡ {u (k) , λ (k) } belonging to the equilibrium path. Starting from a known equilibrium point z 0 ≡ z (k) , the new one z (k+1) is evaluated correcting a first extrapolation z 1 = {u 1 , λ 1 } by a sequence of estimates z j (loops) by a Newton iteration Jż = −R j z j+1 = z j +ż (28) where R j ≡ R(z j ) andJ is the Jacobian of the nonlinear system (27) at z j or a suitable estimate. The method is able to provide the equilibrium path for assigned data even in case of limit points in load or displacements. Its main drawback is the high computational cost, due to the large size of the system of equations. For this reason, the method is inadequate to assess the imperfection sensitivity. Generalized path-following methods are a promising alternative to the standard Riks method, where, for example, the critical point can be evaluated by changing the initial data directly, without reevaluating the whole load-displacement curve [20]. However, in the following, a completely different approach, called the Koiter method, is illustrated. It is based on a reducer order model (ROM) and is far more efficient for imperfection sensitivity analyses compared to path-following methods, and then more suitable for design purposes. It is worth citing the possibility to couple the Koiter method with the path-following approach, obtaining the so-called Koiter-Newton approach where the ROM is used as an accurate predictor [29,30].

The Mixed Integration Point Strategy
In geometrically nonlinear analysis, a mixed format of the nonlinear equations in stress and displacement variables provides superior performances in the solution methods. In the standard path-following method, the mixed iterative process assures a greater robustness also for large steps and with a reduced number of iterations, and then a reduced computational time [54]. In Koiter analysis, a mixed format is an indispensable prerequisite to obtain an accurate ROM (see [25,54]). The improved efficiency in path-following methods and accuracy in Koiter's method is much more evident when the slenderness of the structure gets higher and, concerning the last one, when the pre-critical path presents some nonlinearities. The mixed integration point (MIP) strategy proposed in [17] can be successfully used to exploit the advantages of a mixed format in Riks and Koiter analysis [13,36] without the need of a stress interpolation.
The main idea of the MIP method is to relax the constitutive equations at each integration point by rewriting the strain energy in a pseudo Hellinger-Reissner form as where the stresses at each integration point σ g become independent variables: The stationary condition with respect to σ g gives the constitutive law at g: When Equation (31) is exactly solved and substituted in Equation (29), we obtain, again, the displacement formulation. This means that the discrete approximation of the problem is the same as in the original displacement formulation: the equilibrium path is the same when a path-following scheme is adopted (see [17]).
In this way, MIP formulations extend the results already obtained for mixed (stressdisplacements) discrete approximations, avoiding the use of a stress interpolation. This is particularly convenient in IGA, where an effective mixed formulation is not an easy task. Moreover, the MIP method was recently used to solve, by means of collocations, the strong form of the problem equations [55,56] or applied to more involved constitutive laws [57].

Koiter Method
The Koiter approach, described in detail in [25,28,36], is here briefly recalled its main algorithmic steps based on the MIP formulation of the solid-shell model. As shown in [54], this is necessary to improve Koiter's method accuracy because of the direct prediction of the stress and efficiency, due to the vanishing of the fourth order strain energy variations.
By collecting in vector u the global discrete displacements d and stresses σ g at each integration point g, Koiter's method is based on the following reduced order model: where ψ i are the scalar modal amplitudes,û is the linear elastic solution (path tangent to the stress-free configuration),v i denotes the ith of the n linearized buckling modes, w ij and w are quadratic correction modes. The evaluation of these vectors requires the solution of linear systems forû, w ij andŵ and a linearized buckling analysis forv i . Details can be found in [25,28,36]. According to this choice, the equilibrium path is approximated by the following nonlinear reduced system of equilibrium equations in the unknowns ψ i λ: The scalar coefficients A ijk , C ik , B ijhk and µ k [λ] are computed as the sum of elemental contributions of the stored energy variations. Their explicit expressions can be found in [25,36].
A notable feature of the method is a computationally efficient imperfection sensitivity analysis. In fact, the imperfect structure can be studied by perturbing a posteriori the same reduced system of the perfect structure by adding to it the imperfection coefficientsμ k : This means the analysis of a new geometrical imperfection simply requires only to updateμ k and solve again the small system in Equation (34). Thousands of imperfections can be analyzed in a few seconds regardless of the global number of DOFs used for the full structural model. Two strategies were proposed for the evaluation ofμ k [28]. A first solution is very quick but with a validity range restricted to small imperfection amplitudes and almost linear pre-buckling path. A second imperfection modeling is more accurate for a wider range of structural problems and just a little more expensive than the first one [28].

The Worst-Case Geometrical Imperfection
The geometrical imperfectiond can be assumed to be a linear combination of assigned shapesḋ i with combination factorsψ i , scaled in order to have an assigned maximum amplituded max chosen, for example, from experimental measurements as in [58]. Shapesḋ i can be chosen, for example, as the displacement part of the first linearized buckling modes. The worst-case imperfection can be defined as that leading to the worst value of P (α g ,ψ 1 , . . . ,ψ m ): The solution of the previous problem can be obtained by the stochastic algorithm proposed by [26,27], which exploits the reduced order modeling of Koiter's method.

Advantages of a TL Solid-Shell Formulation
Discrete models directly derived from the 3D continuum using the Green strain measure have a low order dependence on the strain energy from the discrete parameters and in particular have a 3rd order dependence for the MIP formulation. On the contrary, geometrically exact shell and beam models or those based on corotational approaches make use of the rotation tensor. This implies that the stored energy is infinitely differentiable with respect to the discrete parameters and leads to very complex expressions for the energy variations with a high computational burden for their evaluation. On the contrary, for solid-shell finite elements, the strain energy in MIP form has the lowest polynomial dependence on the corresponding discrete parameters, i.e., just one order more than in the linear elastic case, implying the zeroing of the fourth order energy variations required to build the Koiter ROM.

Postbuckling Optimisation Algorithms
The optimisation problem is based on searching for the material distribution that maximizes the structural performance, i.e., the objective function in Equation (18).

Monte Carlo Random Search with Zoom Steps
The Monte Carlo random search method proposed in [26] consists of multiple steps and requires only the computation of the objective function. During the first step, a random population of N 1 layups is generated and, for each of them, the objective function is computed. The n = n 1 elite (best) solutions, collected in α el , represent the starting points of the second step (zoom step) that improves the previous elite. For each value in α el , the objective function is evaluated N 2 times at random points defined as where j = 1 . . . n denotes the elite value and rnd is a generator of pseudo random integer values between −R and R. The radius R can decrease during the steps, e.g., R 1 during the first zoom step and R 2 for the others. At the end of a zoom step, n best solutions are selected as the new elite population to start the next step. The algorithm stops if a satisfactory convergence is obtained; otherwise, a next zoom step is started. Although very simple, the Monte Carlo search with zoom steps provides good estimates of the optimum, for practical design purposes, with a limited number of objective function evaluations. These kinds of methods become, however, more and more demanding by increasing the number of design variables.

Genetic Algorithm
A genetic algorithm (GA) is a metaheuristic method inspired by the natural selection [32,59]. It is widely used to solve optimization problems when the derivatives of the objective function are difficult or costly to evaluate. Compared to a Monte Carlo random search, it relies on biologically inspired operators such as mutation, crossover and selection. The method starts from a random population and proceeds with an iterative process, called generation, providing new individuals. At each generation, the value of the objective function (called fitness) of every individual is evaluated. The best individuals are stochastically selected and recombined also with random mutations to form a new generation for the next iteration. The termination criterion is a maximum number of generations or a satisfactory fitness level. In addition, the genetic algorithm becomes, however, costly for a large number of design variables.

Globally Convergent Method of Moving Asymptotes
The optimal design problem can be also solved via a gradient based method with convex subsequent approximations of the objective function, i.e., the Global Convergent Method of Moving Asymptotes (GCMMA) [33,[60][61][62]. It is particularly suitable for the optimisation of objective functions requiring a high computational cost and depending on many variables.
The gradient evaluation with respect to the design variables is usually done numerically. In particular, the ith component is computed as where h is a small real parameter chosen defining the discrete incremental ratio, and i i is a basis vector where the unitary ith component is the only non-zero element. Generally, the gradient evaluation is extremely time-consuming, but the relative efficiency of Koiter's method reduces it to an acceptable time. Moreover, GCMMA generally converges with a few iterations, so compensating largely the gradient computation. A multi-start GCMMA is advised for complicated objective functions with multiple local optima.

Comparison of the Optimization Algorithms
According the numerical experience matured in [26,27,35] in post-buckling optimisation, the performances of the different algorithms can be compared. In particular, the main advantage of Monte Carlo with zoom steps and genetic algorithm is a likely convergence to the global optimum and no need for gradient computation. The first one usually provides better estimates for a low number of function evaluations, but it is slower to achieve a converged solution compared to the second one. However, both methods get inefficient (large number of function evaluations) when the number of design parameters increases. On the contrary, GCMMA is the most efficient choice for a large number of design variables because the cost for computing the function gradient is largely compensated by the very fast convergence. As a drawback, multiple starts are generally needed to assure convergence to global optima.  Figure 2 and imposed as usual by zeroing the corresponding DOFs. For the axial displacement, the rigid motion is prevented by maintaining the symmetry.

Post-Buckling Optimisation of a Cylindrical Panel under Compression
A postbuckling optimisation of this structure was previously carried out in [26] using straight fibre laminates and a Monte Carlo methodology. In addition, it was used as a benchmark test for the isogeometric formulation proposed in [13]. The panel is particularly suitable to assess the potentiality of the material optimisation, since its baseline Quasi-Isotropic (QI) design, that is, when all lamination parameters are zero, exhibits an unstable post-buckling behaviour with imperfection sensitivity. In the following, we will denote with n 1 × n 2 the number of isogeometric elements in the two directions of the shell domain. The NURBS isogeometric discretization is adopted for both material description and displacement field. The number of elements for the two discretizations is not necessarily the same. In particular, for the structural analysis, the displacement discretization is that giving converged results, while, for the material description, different meshes are considered. The optimisation problem at stage 1 consists of finding the distribution of lamination parameters that maximises the failure load, i.e., the minimum between the first limit load λ lim and the load leading to a limit axial displacement of the loaded section v ax = v d . In practice, the process searches for the axially-stiffest stable behaviour or the unstable configuration with the highest limit load and is less sensitive to geometrical imperfections. The loads in the objective function (18) are normalised by the first linearised buckling load of the QI case, named λ QI .
The shell thickness is t = 10 mm. The axial displacement limit is set to v d = 2 mm. The material is a widely employed E-glass/epoxy fibrecomposite with E 11 = 30.6 GPa, E 22 = 8.7 GPa, ν 12 = 0.29 , ν 23 = 0.5 , G 12 = 3.24 GPa, G 23 = 2.9 GPa. The discretization of the displacement field is based on quadratic NURBS and a 9 × 9 control grid that gives a converged solution for different material configurations. Consequently, the number of discrete displacement DOFs is 683. The maximum amplitude of the geometrical imperfection isd max = 0.5 t.
Two different grids based on quadratic NURBS, namely 6 × 6 and 9 × 9, are used to parametrize the lamination parameters. To preserve the symmetry, only the control points of a quarter of the structure are taken as independent design parameters. An orthotropic laminate with symmetric and balanced stacking sequence is considered. Under these choices, the number of independent design variables is restricted to 64 and 144, respectively, for the two grids.
The convergence of GCMMA for the two different description of the lamination parameters is illustrated in Figure 3. The optimal solution is achieved by less than 10 iterations. Moreover, a good solution is also given for the coarsest material mesh. The evolution of the equilibrium path during the iterative process is shown in Figure 4 for the 9 × 9 material grid. Each equilibrium curve corresponds to the lamination parameters obtained at the end of each GCMMA iteration, whose number is reported in the legend. For example, the equilibrium path indicated with 0 corresponds to the starting point of GCMMA, which is the QI configuration. The curve denoted with 1 is relative to the first iteration of GCMMA, and so on. Interestingly, the solution gets significantly better already after the first iteration compared to the initial QI configuration. The third iteration gives the first stable postbuckling response. Then, the computational design process finds configurations with increasing axial stiffness and the GCMMA algorithm converges soon afterwards: the equilibrium curve after the 7th iteration is practically the same as that found at the 20th one.
The results of stage 1 are summarized in Figure 5 in terms of equilibrium curves. In particular, the load-displacement path corresponding to the optimal distribution of lamination parameters (LP) is compared with the baseline QI case. Comparison is also made with the optimal solution for straight fibre laminates with a non-symmetric stacking sequence (SF) obtained in [26], that is, from the inside out, [∓54, 0 4 ] with the fibre orientations given with respect to the local reference system in Figure 2 with e 3 aligned with the surface normal vector. Figure 5 compares the equilibrium path predicted by Koiter's method with that given by the Riks method using the full model. The end shortening and the out-of-plane displacement at the centre of the panel are monitored. We can observe a satisfactory agreement in all cases. The optimal configuration of lamination parameters is illustrated in Figure 6. Finally, the worst-case imperfection detected during the optimal design process is pictured in Figure 7.  On the left, the axial displacement (v ax ) is plotted while, on the right, we have the out-of-plane displacement at the panel center (w out ). Baseline QI case, straight fibre design [26] and optimal lamination parameters are considered. v ax w out Figure 5. Comparison of the equilibrium paths obtained using Koiter's method with Riks results. On the left the axial displacement (v ax ) is plotted while on the right we have the out-of-plane displacement at the panel center (w out ). Baseline QI case, straight fibre design [26] and optimal lamination parameters are considered.

Stage 2: Recovery of the Lamination Angles
The advantage of the two-stage methodology is the possibility of designing structures made of many layers without penalizing the computational cost of the process. Actually, having many layers is an advantage, since it allows a more accurate match with the optimal solution of stage 1. This fact is shown here, where the recovered stacking sequence (SS) is symmetric and balanced with 12 independent layers. The stacking sequence is restricted to be [±θ 1 , . . . , ±θ 12 ] s , with each θ interpolated using a bivariate quadratic NURBS (see Equation (24)) over a mesh of 9 × 9 elements. The thickness of each layer is 0.208 mm.
The solution of the problem (25) is achieved by means of a multi-start GCMMA. The algorithm converges quickly, as shown in Figure 8. The corresponding equilibrium curves are reported in Figure 9, while Table 1 shows the match between the values of objective function and the first linearised buckling loads for the optimal lamination parameters and those corresponding to the recovered fibre orientations. These last ones are pictured in Figure 10.  Figure 9. Load-axial displacement v ax paths for the baseline QI case, the optimal distribution of lamination parameters (LP, 9 × 9 elements) and the fibre orientations retrieved in stage 2 (SS). Table 1. Results of the optimised material design LP for a grid of 9 × 9 control points for the lamination parameters and for the retrieved fibre orientations (SS).

Conclusions
This review collected some recent findings in the post-buckling optimisation of thin-walled composite structures. Focus was given to structural modelling, material parametrization, post-buckling analysis with imperfection sensitivity and optimisation algorithms. A computational framework merging the developments in these different aspects led to robust and efficient optimizations of composite materials with spatially varying material proprieties. As a result, numerical analysis can make it possible to fully exploit the capability of advanced manufacturing methods for the realisation of a new generation of structures able to work safely in the post-buckling regime saving materials, costs and weight.

Conflicts of Interest:
The authors declare no conflict of interest.