Lightweight design of variable-angle filament-wound cylinders combining Kriging-based metamodels with particle swarm optimization

Variable-angle filament-wound (VAFW) cylinders are herein optimized for minimum mass under manufacturing constraints, and for various design loads. A design parameterization based on a second-order polynomial variation of the tow winding angle along the axial direction of the cylinders is utilized to explore the nonlinear steering-thickness dependency in VAFW structures, whereby the thickness becomes a function of the filament steering angle. Particle swarm optimization coupled with three Kriging-based metamodels is used to find the optimum designs. A single-curvature Bogner–Fox–Schmit–Castro finite element is formulated to accurately and efficiently represent the variable stiffness properties of the shells, and verifications are performed using a general purpose plate element. Alongside the main optimization studies, a vast analysis of the design space is performed using the metamodels, showing a gap in the design space for the buckling strength that is confirmed by genetic algorithm optimizations. Extreme lightweight while buckling-resistant designs are reached, along with non-conventional optimum layouts thanks to the high degree of thickness build-up tailoring.


Introduction
The development of new concepts for lightweight structures has been overwhelmingly exploited in several fields, mainly in aeronautical and aerospace structures, to comply with Green Aviation for enhancing fuel efficiency and decrease aviation emissions towards reaching carbon-neutral air transportation. A practical and direct manner to improve the energy efficiency and reduce the fuel consumption of an aircraft is by reducing the mass of its components (Zhu et al. 2018), achieving carbon footprint reduction and better flight performance. Both aeronautical and space industries have been continuously developing new concepts for lightweight structures to pursue these benefits, where fiber-reinforced composite materials play a major role. Figure 1 illustrates the evolution of the use of composites in aeronautics over the last 60 decades.
In 2017, NASA designed, optimized, manufactured and tested a variable-angle tow (VAT) wing for optimal passive load alleviation, proving that the use of fiber steering through automated fiber placement (AFP) manufacturing is more structurally efficient than conventional Responsible Editor: Qing Li. Zhihua Wang, José Humberto S. Almeida Jr. and Saullo G. P. Castro have contributed equally for this work.

3
140 Page 2 of 23 straight-fiber-reinforced composite wings (Brooks and Martins 2018). Liguori et al. (2019) optimized a wingbox section using genetic algorithms (GA) including manufacturing constraints of the AFP process. They achieved a wingbox 6.48% lighter, with higher buckling strength when compared with the initial straight-fiber baseline structure. Tailoring a particular property for laminated composites by means of in-plane fiber angle variation is proven to be efficient for variable-angle tow (VAT) composites. For example, in the first attempt to develop such VAT composites, Hyer and Lee (1991) obtained higher buckling loads for simply supported open-hole laminates under uni-axial compression loading. Thenceforth, several approaches have been developed based on the advancement of both computational approaches and manufacturing techniques for composite materials.
Modern aircraft have their fuselages made of carbon fiber-reinforced polymer (CFRP) composites, which makes cylindrical unstiffened shells (Degenhardt et al. 2014) of particular importance towards understanding and further reducing aircraft overall mass. In spacecraft, monocoque structures are commonly built of unstiffened shells, requiring conservative knock-down factors to account for the high sensitivity of the buckling strength towards geometric and load imperfections (Castro et al. 2014a). The full potential of CFRP composites in these components is still underexplored, and an efficient strategy to take full advantage of these materials consists of reducing their masses while keeping their structural performance. The state-of-the-art on mass minimization of composite shells is scarce, as discussed next, with most of the studies focusing on optimization problems that seek minimum compliance for a given mass. Pelletier and Vel (2006) developed a multi-objective optimization of composite cylinders for strength, stiffness and minimal mass via the layerwise tailoring of fiber angles and fiber volume fraction. Although some enhancements were found, it was difficult to find a layout with reasonable compromise among the three objective functions utilized. Sadeghifar et al. (2010) optimized stiffened cylindrical shells for minimum mass and maximum axial buckling load using a GA, finding that I-section stiffeners are more effective than rectangular ones to achieve panels with reduced mass. Ghasemi and Hajmohammad (2017) proposed a multi-objective optimization of composite shells under external pressure for minimum mass and maximum buckling pressure using GA. They reached a mass reduction of 12% while respecting the desired buckling pressure.
In general, optimizations based on GA require thousands of function evaluations, becoming therefore limited due to the high computational cost. An appropriate strategy to improve the efficiency of these optimizations is to employ metamodels, which are computationally efficient surrogate models that approximate the high-fidelity response using approximated functions, such as hypersurfaces. Rouhi et al. (2015) optimized VAT cylindrical shells for maximizing the buckling load using radial basis functions as surrogates and GA to find the optimum in the surrogate model. Blom et al. (2009) used design explorer tools (Booker et al. 1999) as surrogates to minimize the buckling load of VAT composite cylinders.
Recently, Wang et al. (2020) developed a new accelerated Kriging metamodeling scheme to restrict the area of interest within an optimization framework to maximize the buckling load of VAT composite cylinders manufactured by filament winding, designated as variable-angle filament-wound-VAFW.
This non-extensive state-of-the-art revealed that: -no work on mass minimization for VAT composite cylinders manufactured by filament winding (FW), the most suitable and fastest manufacturing process for composite cylinders; -lack of works using metamodels for training FE models for minimizing mass under buckling constraints of VAT cylinders; -lack of exploitation of several metamodels to find the most suitable one for VAFW cylinders; and -no reports exploiting the tow overlapping feature of the FW process.
To fulfill these gaps, the present work proposes an optimization framework to minimize the mass of VAFW cylinders constrained by design loads for uni-axial buckling strength. Three Kriging-based metamodels are developed to approximate the buckling load and mass of the VAFW cylinders, namely: classical Kriging, Co-Kriging, and Accelerated-Kriging. The finite element analyses used to accurately calculate the buckling load and mass are performed using a single-curvature Bogner-Fox-Schmit-Castro (BFSC) (Castro and Jansen 2021), consisting of an enriched rectangular four-node element with ten degrees-of-freedom per node to achieve third-order interpolation for the displacement field. The high-order approximation combined with 16 integration points makes this element a suitable and fast-converging approximation to represent shells with variable stiffness. The FE model takes into consideration the manufacturing characteristics of the FW process, and can accurately represent the variable stiffness properties according to the proposed design parameterization. The FE predictions from the single-curvature BFSC model for different design loads are verified against a general purpose shell element of a commercial FE software. The design optimizations are performed using Particle Swarm Optimization (PSO) coupled with the Kriging metamodels, and the optimum designs are verified against a conventional GA optimization coupled with the BFSC finite element models. A vast analysis on the design space is performed using Kriging-based metamodels, showing the ranges of mass and buckling strength for different number of layers; with the possible mass range for different design loads also determined.

Design parameterization
The pioneering work of Wang et al. (2020) on reliabilitybased design optimization of variable-angle filament-wound (VAFW) cylinders, compared different VAFW designs and found that the maximum buckling load was achieved by the design named VAFW-P, with "P" being an abbreviation for parabolic, indicating that the winding angle distribution (x) varies along the longitudinal direction following the second-order Lagrange polynomial given in Eq. 1 (de Quadros and Hernandes 2018). Each angle-ply layer consists of two plies of angles ± (x) controlled by three design variables VP 1 , VP 2 and VP 3 that are the control points for the secondorder Lagrange polynomial. In order to produce an angle distribution that is symmetric about the middle cross section at x = L∕2 , the control points are positioned in the order Figure 2 illustrates an incomplete angle-ply layer where the blue region is a group of adjacent filaments that create the ply with the positive angle distribution + (x) . The gray region represents the ply with negative angle distribution − (x) . Note how the variable-angle filaments overlap when moving from regions of lower (x) to regions of larger (x) values. In the present work, the VAFW-P parameterization is used to find optimum designs of minimum mass for different given design load levels, whereby the number of angle-ply layers and the winding angle distribution in each layer are the design variables to be determined.
where: As explained in Wang et al. (2020), for VAFW cylinders the thickness distribution of the kth ply given by h k (x) changes due to filament bending, filament shearing or any combination of both (Wang et al. 2020), and the thickness increase is necessary to keep a constant volume (Castro et al. 2019;Vertonghen and Castro 2021), leading to the following nonlinear dependence of h k (x) with the winding angle distribution (x): with h tow being the filament thickness; and (x) the filament steering angle. The steering-thickness coupling of Eq. 3 poses yet another challenge for the design and optimization, where the mass of VAFW cylinders becomes not only a function of the number of angle-ply layers, but also depends on the winding angle distribution (x) . For FW process, it is still unknown whether the variable-angle filaments are driven by either tow bending or shearing (Wang et al. 2020). However, this assumption does affect the variable thickness calculation described in Eq. 3, hence not effecting the stiffness nor the inertia models. The hoop spacing between two adjacent filaments at a given longitudinal position over the cylinder is determined by Wang et al. (2020): According to Eq. 4, when max ( (x)) is chosen to define a constant value for c to be used throughout the cylinder, gaps appear in regions where (x) < max ( (x)) , as illustrated in the gap design of Fig. 3. Conversely, if min ( (x)) is used to calculate c , it creates overlaps where (x) > min ( (x)) , as depicted in the overlap design of Fig. 3. The present study adopts the overlap design for all cases, meaning that the filament steering angle (x) of Eq. 3 is calculated with respect to the minimum filament angle as (x) = (x) − min ( (x)) , for 0 ≤ x ≤ L.
The geometry of all cylinders herein optimized is fixed, with a radius of r = 0.15 m , and length of L = 0.3 m . The orthotropic material properties of the filaments are: E 11 = 90 GPa , E 22 = 7 GPa , 12 = 0.32 , G 12 = G 13 = G 23 = 4.4 GPa ; with a filament thickness of h tow = 0.4 mm . The material density is = 1600 kg∕m 3 .

Design optimization
The VAFW designs are defined by the number of layers N L and the angle distribution for each layer, controlled by VP i,L with i = 1, 2, 3 and L = 1, … , N L . The angles at the control points are limited by a lower-bound value L , being the minimum winding angle; and by a upper-bound value U , being the maximum winding angle. The designs are optimized for minimum mass M( VP i,L , N L ) subject to different design load levels d = 50 kN , 100 kN , 200 kN , 500 kN , and 1000 kN , and constrained by linear buckling cr ( VP i,L , N L ) ; in the following optimization problem: The proposed optimization scheme is illustrated in the flowchart of Fig. 4. Before the optimizations for a given design load, the metamodels are trained for the range of desired design loads to map the boundaries of linear buckling load and mass for each number of layers. The optimization starts with the definition of the design load. Next, the metamodels are trained to approximate M( VP i,L , N L ) and cr ( VP i,L , N L ) , as detailed in Section 5. The training procedure is performed by adding new samples in every successive iteration based on the predicted maximum mean squared error (MSE). With the trained metamodels, the boundaries of linear buckling load for a given number of layers N L are determined by means of Latin Hypercube Sampling with 5000 samples.
The optimization to find the designs with minimum mass for a given design load is summarized as: Part i: Find candidate cylinders ( N L = 1, 2, 3, 4) for a given buckling load: -Exploring the lower and upper boundaries in terms of buckling loads for cylinders with different number of layers ( N L = 1, 2, 3, 4); -Find and keep candidate cylinders whose buckling loads are higher than the design load.
Part ii: Find the optimum design with minimum mass while respecting the buckling load boundaries defined in Part i: -build two surface responses: buckling load and mass; -run the PSO optimizer using the fitness function defined in Eq. 5, where the buckling load and mass metamodels are called at each iteration; -repeat the previous step until the global minimum for the mass is reached.

Linear buckling constraint modeling
The linear buckling constraint cr ( VP i,L , N L ) is calculated by means of finite element modeling based on a single-curvature Bogner-Fox-Schmit-Castro element (SC-BFSC). This element is herein developed to combine the high-order interpolation of the BFSC with cylindrical shell kinematics, aiming to achieve a computationally efficient and yet accurate representation of the variable stiffness shells under investigation. The BFSC finite element (Bogner et al. 1966;Castro and Jansen 2021) is a C1 contiguous confirming plate element obtained from tensor products of cubic Hermite splines. With 4 nodes per element and 10 degrees-of-freedom per node, the BFSC approximates the in-plane and outof-plane displacements using third-order polynomials, being a fast-converging method for an accurate representation of the linear buckling constraint. Figure 5   coordinate y is curvilinear following the circumferential perimeter, such that at y = 2 r the path closes on itself. The nodal connectivity is also indicated, and the mesh is built to keep the element edges parallel to the x, y coordinates, such that x = x 2 − x 1 = x 3 − x 4 , and y = y 4 − y 1 = y 3 − y 2 . Figure 6 shows the mesh closing on itself at the intersection of elements e in with elements e i1 . The linear buckling problem consists of finding the value of that leads to: where K K K and K G0 K G0 K G0 are respectively the constitutive and geometric stiffness matrices of the system, described in detail in Appendix 1. For a system with n unknown degrees-offreedom, there are n eigenvalues that are roots of the characteristic polynomial obtained with Eq. 6. In practice, the eigenvalues and corresponding buckling modes are solved using generalized eigenvalue solvers that are able to efficiently extract only a desired number of eigenvalues and buckling modes, and in the present work the locally optimal block preconditioned conjugate gradient method (Knyazev 2001) implemented in SciPy (Virtanen et al. 2020) is used.
The integration of the constitutive and geometric stiffness matrices over the finite element domain are performed numerically using standard Gauss-quadrature with 4 × 4 integration points per element. The authors verified that this amount of integration points leads to a converged behavior that is capable to accurately represent the variable stiffness of the VAFW cylinders. For each integration point, the local shell constitutive properties given by matrices A A A, B B B, D D D are calculated based on the winding angle distribution (x) of each angle-ply layer. While computing the constitutive matrices, the local thickness of each ply is consistently calculated according to Eq. 3.
The elements of the constitutive matrices are functions of , and are calculated with: where n is the number of plies; Q ij (x) is the ply stiffness rotated by the local winding angle (x) (Reddy 2003); z k (x) defines the position of the outward face of the kth ply. For a correct representation of the filament winding process, the inner face of the filament-wound cylinder must coincide with the mandrel radius. Therefore, the values of z k (x) are offset by a quantity d(x) such that z k−1 = 0 for the first ply, as illustrated in Fig. 7. The correct d(x) value is half of the total thickness: where h k (x) is the local thickness of each ply as per Eq. 3. Appendix 3 presents a verification of the proposed finite element formulation compared with the multipurpose shell element S4R from Abaqus Smith (2019).

Kriging metamodels
Three Kriging metamodels are considered, namely: classical Kriging, co-Kriging, and accelerated-Kriging; aimed to approximate the buckling load and mass of the VAFW cylinders possessing different number of layers N L and different winding angles ± . After an initial Latin Hypercube Sampling (LHS), the first metamodel is assembled and further trained until all convergence requirements are satisfied, and this convergence depends on the number of variables and the degree of nonlinearity of the response surface. For all metamodels of the present study, two convergence criteria, CC1 and CC2, are adopted: CC1: The predicted response error of both lower and upper boundaries must be less than a given threshold c. Then, CC1 is determined as: where Error i lower and Error i upper are the errors of the lower and upper boundaries in the ith training process, respectively. The Error is then defined as: where ŷ( c ) and y( c ) are the predicted response and the true response of the candidate sample c in design region, respectively. Three thresholds c are set to 5%, 3%, 1% and applied to find the buckling load and mass limits of all cylinders. Each threshold generates particular limits and they will be useful to select the threshold with the widest range in terms of both buckling load and mass to enlarge the design space, therefore, increasing the chances of reaching higher improvements. CC2: The error for both lower and upper boundaries in two consecutive iterations should be less than the given threshold c, which is determined as: The training process stops if the two converge criteria are met. Then, the lower and upper boundaries of the buckling load and mass can be estimated as: In the following sub-sections, each metamodel is explained in detail.

Classical Kriging
For the classical Kriging metamodel, the spatial association between the design variables and the responses is expressed as: is a vector of regression function which is also known as the trend of the prediction based on available samples, v 1 , v 2 , … , v n is a vector of coefficients, and (x x x) is a Gaussian process with zero mean and covariance cov x x x i , x x x j .
The initial random sample space generated through LHS has a population size of ten times the number of design variables, as recommended by Forrester et al. (2007). After that, the true responses for buckling load and mass are calculated by the FE analysis (FEA) described in Sect. 4. Then, the initial samples and responses are used to construct the initial Kriging metamodel (Lophaven et al. 2002). A training process follows to locate and add new samples into the metamodel, with new training samples selected as the design space location with maximum estimated MSE, being the estimated MSE calculation detailed in Appendix 2. The training is repeated until the convergence criteria of Eqs. 9 and 11 are satisfied. In addition, the lower and upper boundaries of the response in the design region are determined by Eq. 12. A detailed description of the procedure is given in Table 1.

Co-Kriging
In co-Kriging, the metamodel is built with two sets of data, one set named expensive data ( e , y e ) that are more computationally costly, while the other set is called cheap data ( c , y c ) (Zhou et al. 2020). These two sets of data are independent of each other. The response y e at e is supposed to be more accurate than the cheap data y c at c . The required total data are combined as: where r are the required samples of random variables; e and c are random samples from expensive and cheap data, respectively; n e and n c are the number of expensive and cheap data, respectively. In this work, the number of data n e and n c are defined as where n is the number of random variables.
The responses corresponding to the random samples r are shown in Eq. 16, and the combined data ( r , y r ) are the required data to build the co-Kriging metamodel.
n e = 2 × n;n c = 10 × n, Similarly to the classical Kriging approach, co-Kriging realizes the response surface with the required data through Eqs. 14,16. With the cheap data ( c , y c ) , the response of the cheap function is described as a Gaussian process f c ( ) in the total sample space (Lophaven et al. 2002). Then, the expensive response surface is defined as: where f e ( ) is the prediction at , is a scaling factor estimated by the data, and f d ( ) is the difference between expensive and cheap responses. The initial independent cheap c and expensive e samples of random variables are generated via LHS. The buckling loads c of samples c are obtained from the classical Kriging metamodel (Sect. 5.1), and the buckling loads e of samples e are calculated by FEA. Then, the metamodel can be built with the initial required data: ( c , c ) and ( e , e ).
In each iteration of the co-Kriging training, three new training samples are added to the current metamodel. During the ith iteration ( i = 1, 2, ..., n ), the first training sample is the point at where the predicted MSE is the highest in the whole design region, expressed as i , with the corresponding estimated buckling load i . The second and third new training samples are determined as the samples at the current lower and upper boundaries of the ith metamodel, which are expressed as Once the new training samples are available, they are added to the ith metamodel to improve the accuracy of the predictions. The first sample ( i , i ) in the ith iteration is added to the expensive data since it is likely to be accurate. The expensive data in the (i + 1)th iteration is then updated as: where i+1 e , i+1 e is the expensive data in the (i + 1)th iteration, and i e , i e is the expensive data in the ith iteration. Then, the other new training samples ( L i , L i ) and ( U i , U i ) are added to the cheap data. The cheap data in the (i + 1)th iteration is then updated as: c is the cheap data in the (i + 1)th iteration, and i c , i c is the cheap data in the (i)th iteration. After the updated expensive and cheap data are available, the metamodel is then trained. The lower and upper boundaries ( L i+1 , L i+1 ) and ( U i+1 , U i+1 ) are estimated by the extreme value optimization algorithm. In this work, the Particle Swarm Optimization (PSO) algorithm is employed, as explained in section 6. The training process evolves until both convergence criteria CC1 and CC2 are met. Finally, the lower and upper boundaries of the buckling load in the design space can be calculated by Eq. 12. Table 1 Detailed procedure of the classical Kriging metamodel for an arbitrary cylinder with one ± (x) layer Step Procedure 1: Design variables [ 1 , 2 , 3 ] and N = 3 , set i = 1 2: Generate 10 initial samples ( initial ) 3: Obtain the buckling load of initial samples 4: Construct initial metamodel with DACE tool, M K Obtain the boundaries with Eq. 12. Table 2 summarizes the procedure to find the buckling load boundaries through co-Kriging approach, in which the VAFW-P design with one ± (x) layer is taken as a mere example.

Accelerated Kriging
The accelerated Kriging metamodel was first introduced by Wang et al. (2020) and has shown a higher computational efficiency, accuracy and better convergence properties than the classical Kriging. With an increasing number of variables or an increase in the nonlinearity of the response, there is a larger probability that the sample with largest estimated MSE is located at a design space region that is outside from the region of interest, negatively affecting the convergence of the classical Kriging and Co-Kriging methods during the training phase. The proposed accelerated Kriging metamodel is an extension of the classical Kriging approach, whereby a moving search window, or hyper-cube, is used during the training phase. The addition of a new sample near the current optimum can accelerate the convergence, whereas keeping the robustness of the original method.
The training process starts similarly to the classical Kriging metamodeling explained in Sect. 5.1. Next, for each training iteration, two points are added being the first training point the one with largest estimated MSE, similarly to the other two Kriging metamodels previously described. The second training point of the accelerated Kriging approach is added based on the largest estimated MSE within an adaptive moving search window centered at the current optimum. The adaptive moving search window is defined as: where L i+1 and U i+1 are the lower and upper boundaries of the (i + 1)th adaptive window, respectively; The i opt is the ith optimal design in the ith iteration of the training procedure; W is the range of the window, which is determined as: where L and U are respectively the lower and upper boundaries of the design space, while q is a scale factor, which ranges from 0 to 1.
Considering that the center of the window is the previous optimal result (Eq. 20), the search window adapts itself by moving during the training iterations. In this work, a constant window width of q = 0.25 is used. The lower and upper boundaries of the ith adaptive moving window is replaced by the corresponding lower and upper boundaries of the design space if the ith window moves out of the design space, as follows: With the first and second training samples determined, the corresponding responses are calculated by FEA. The iterations are repeated until the convergence criteria CC1 and CC2 are met. Table 3 describes a summary of the proposed accelerated Kriging approach for a VAFW-P cylinder with one ± layer. Table 2 Detailed procedure of the co-Kriging metamodel for an arbitrary cylinder with one ± (x) layer Step Procedure 1: Design variables [ 1 , 2 , 3 ] and N = 3 , set i = 1 2: Generate 10 initial samples ( initial ) 3: Obtain the buckling load of initial samples 4: Construct initial co-Kriging metamodel M K Obtain the boundaries with Eq. 12.

Table 3
Detailed procedure of the accelerated-Kriging metamodel for an arbitrary cylinder with one ± (x) layer Step Procedure Obtain the boundaries with Eq. 12.

Design space analysis for buckling loads
After the convergence criteria (CC1, CC2) and thresholds ( 5% , 3% , 1% ) are met for all three Kriging-based metamodels, the lower and upper boundaries for the buckling loads are estimated. Table 4 shows the obtained boundaries and the required number of samples N s used in the training process for the VAFW-P design with one layer ( N L = 1 ). It is observed that the accelerated-Kriging metamodel provides the widest ranges in terms of buckling loads for all thresholds applied. Focusing specifically on the threshold c = 1% , which provides wider boundaries, the more the number of layers, the higher the number of new samples ( N S ) required to the training process, as shown in Table 5. Therefore, given the wider boundaries mapped by the accelerated-Kriging metamodel, which is an indicative of a higher accuracy, its results are used to estimate the buckling loads and mass for N L = 2, 3, 4. Table 6 summarizes the estimated lower and upper boundaries of linear buckling load obtained with the accelerated-Kriging metamodel for the VAFW-P cylinders with N L = 1, 2, 3, 4 . Note the gap in the design space between N L = 1 and N L = 2 , with the upper boundary buckling load of 68, 768N for N L = 1 jumping to the lower boundary buckling load of 146, 190N for N L = 2.

Mass metamodeling
The mass of the cylinders is calculated as: where m is the mass, is the density of the carbon/epoxy composite material ( 1600 kg∕m 3 ), x is the circumferential coordinate, z is the longitudinal coordinate, r is the radius of the cylinder, and h(x, z) is the thickness.
To calculate the mass for the cylinders with several layers, a simple function to sum of the mass of every layer is considered: where m i is the mass of the ith layer, and i is the winding angles in the ith layer. Once the boundaries of the mass for two layers is computed, the boundaries for the other cylinders can be calculated by Eq. 24. Then, the lower and upper limits for all cylinders are listed in Table 7. Figure 8 illustrates the hyper-surfaces of both buckling loads and mass for the cylinders trained with the accelerated-Kriging metamodel using Monte Carlo sampling (MCS). Considering that there are three design variables, one of them is kept constant at 45 • to allow a convenient graphic representation. Then, after both metamodels for buckling loads and mass are trained, the optimization follows, as detailed next. Figure 9 illustrates the outcome of the trained metamodels by plotting the buckling load versus the mass for several sample sizes using MCS approach. Throughout the optimization procedure the samples are dynamically generated, as explained in Section 6. From Fig. 9 there      is a clear trend of higher buckling loads for higher mass values, and it can be observed a complex mass distribution created by the coupled steering-thickness coupling of the VAFW-P cylinder design. From the sampling with 2000 samples of Fig. 9f, the minimum and maximum intervals of mass and buckling loads were identified for N L = 1, 2, 3, 4 , as shown in Fig. 10. Note that for design loads below ≈ 69 kN only designs with N L = 1 should be considered, whereas for a design load of ≈ 700 kN , all designs with N L = 2, 3, 4 might lead to optimum solutions and should therefore be considered by the optimization algorithms. Figure 10 also illustrates the gap in the design space between ≈ 69 kN and ≈ 146 kN , previously observed in Table 6. 6 Optimization framework

Particle Swarm Optimization (PSO)
PSO is a powerful meta-heuristic technique inspired on the swarm behavior of flocks of birds and shoals of fishes, and developed by Kennedy andEberhart in 1995 (Kennedy andEberhart 1995). The method is also related to evolutionary computation and it has some ties to GA. The key advantages of PSO are that it needs only primitive mathematical operations and it is computationally cheap in terms of memory and speed requirements. In PSO, a swarm with particles P is placed in the design search space. Each ith particle in the swarm has a position vector X t i = x i1 , x i2 , x i3 , … , x in T and velocity vector , v in T at iteration t. Each particle is updated according to the dimension j, as per Eq. 25. Every particle is a potential candidate solution and they move continuously around the design space until the optimization meets the convergence criteria. The new particle locations are updated at every iteration towards finding the global optimum. The velocity of the particles can be defined as: whereas the new position is calculated as: in which V d i X d i are the velocity and position of the ith particle, respectively; c 1 and c 2 are acceleration constants; r t 1 and r 2 2 represent two random numbers distributed in the range (0,1); pbest and gbest are the best position so far of the ith particle and the best global position, respectively; w is the inertia weight, which keeps a balance between local and global searches. Note that, from the sampling previously performed and shown in Fig. 9, it is possible to determine all possible values of N L for a given design load, as illustrated in Fig. 10, and for the region with a gap in the design space from ≈ 69 kN to ≈ 146 kN it was assumed N L = 2.
The parameters utilized in the PSO-metamodeling (PSO-MM) optimization framework are presented in Table 8.
The overall mass minimization procedure to meet desired buckling loads BL d can be summarized as follows: S1: For any given desired buckling load, the number of layers is determined with the boundary analysis of the buckling load shown in Table 6. The candidate number of layers N is determined as in which, The optimization flowchart is illustrated in Fig. 4, starting by calling both buckling load metamodels (28) i, j ∈ [1, … , 4], where i ≤ j in each iteration. Then, the PSO algorithm finds the minimum mass, M min n , within the buckling load boundaries, as detailed in the following steps: S4: If n < j , n = n + 1 , go back to Step3. Else, go to S5; S5: Obtain all the candidate optimal designs and corresponding masses: S6: Select the sample with the minimum mass of all candidates using Eq. 29 as the optimum for a given design load BL d . Then, the final output, opt final , is generated;

Verification of the PSO-MM framework
Genetic algorithm (GA) coupled with the developed SC-BFSC finite element model is used to verify the PSO-MM optimization framework. The GA is implemented in Open-MDAO, an open-source optimization framework (Gray et al. 2019). It is worth mentioning that the number of layers allowed is the same as in the PSO-MM optimizations, that is, N L = 1, 2, 3, 44 . The parameters utilized in the optimizations are presented in Table 9. The objective function is similar to Eq. 5, whereby the function is changed in accordance with the GA problem definition as shown in Eq. 30: where M is the mass calculated according to Eq. 24; is the buckling load of the individual; and d is the design buckling load.

Discussion
The objective function of the optimizations is to minimize the cylinder mass subject to the following design loads: 50 kN, 100 kN, 200 kN, 500 kN, and 1000 kN    results consist of only one optimum design per design load. The results are given in Table 10. From Table 10, the following characteristics of the optimum designs, and the differences between the PSO-MM and GA results are identified: -50 kN: this design load is covered by candidates with N L = 1 , as illustrated in Fig. 10. The GA results is slightly heavier than the PSO-MM result; -100 kN: from Fig. 10 there is no design covering this design load, and this region represents part of a gap in the design space from ≈ 69 kN to ≈ 146 kN . The PSO-MM optimization found a candidate cylinder with N L = 2 , with a buckling load of 165.33 kN, expectedly higher than the 100 kN design load, due to the gap in the design space previously identified. Interestingly, the GA result confirmed the identified gap by finding an optimum solution outside the gap region, with a buckling load of 288.15 kN, and a comparable mass with respect to the optimum design obtained by the PSO-MM; -200 kN: the only candidate for this design load level are those with N L = 2 . The combination of angles leads to shell with a mass of 0.7309 kg and buckling load of 207.61 kN. The GA, nevertheless, does not find an optimum with buckling load around the desired one. However, the optimum found by GA has a similar mass with a substantially higher buckling load of 242.83 kN; -500 kN: as per Fig. 10 this design load is covered by cylinders with N L = 2 and N L = 3 , and the PSO-MM optima have respective masses of 1.2715 and 1.0951 kg; with buckling loads of 498.29 and 504.28 kN. The GA finds a unique optimum candidate with N L = 3 , with optimum mass comparable to the PSO-MM result, and a higher buckling load of 659.80 kN; -1000 kN: three candidates are found by the PSO-MM to respect this desired load, with 2, 3, and 4 ± layers. They have masses between 1.46 and 1.61 kg with buckling loads within a range of 1021.5 and 1092.2 kN, respectively. The GA finds only one candidate with 4 ± layers, which has similar mass to the candidate with 4 layers found by PSO-MM framework, with 1.46 kg but with a higher buckling load of 1301.7 kN.
The GA optima reached considerably higher buckling loads without resulting in additional mass for the design loads of 200, 500, and 1000 kN, shown in Table 10. Note that this is not due to the objective function of Eq. 30 utilized in the GA optimizations, because the buckling load only participates in the penalty function when < d , with d being the design load. The PSO-MM optimum designs listed in Table 10 are illustrated in Figs. 11,12,13,14,15,16,17,and 18. The illustrated layers k i are stacked in an increasing z direction while constructing the filament-wound structure, as illustrated in Fig. 7. The corresponding first linear buckling modes are shown in Fig. 19. When VP 1 , VP 2 , VP 3 are different within the same ± layer, a variable-angle filament path is produced, which is well illustrated in some of the figures. Due to the tow bending and tow shearing mechanisms discussed in Wang et al. (2020), the variable-angle path creates overlapping regions with higher thickness between the filaments, here represented by darker regions. This coupling between the thickness and the steering angle leads to interesting optimum designs, such as the ones illustrated in Figs. 14, 16 and 17. The predicted thickness patterns showing overlapping filaments correspond to what has been experimentally observed by the authors in VAFW cylinders of smaller diameters Almeida et al. 2021), and by Labans and Bisagni in large cylinders manufactured by automated fiber placement (Labans and Bisagni 2019).
For design loads that allow an optimum configuration for different values of the total number of layers N L , a general trend is observed, in which lower N L values can achieve the design load by means of more thickness build-up created by the variable-angle filaments. On the other hand, higher N L values tend to show less thickness build-up, thus with less variability in the filament angles. Another trend that is observed is that the buckling modes for the configurations with more thickness build-up tend to show the well-known diamond patterns, whereas the optima with more N L have more axisymmetric buckling modes.
In order to illustrate the efficiency of both heuristic optimization algorithms herein used, Fig. 20 presents convergence results for both optimizers for an example design load of 50 kN. It is noted that the PSO-MM framework takes only

Conclusion
This work focused on develop an optimization framework coupled with metamodels to minimize the mass of variable-angle filament-wound (VAFW) cylinders. The newly proposed single-curvature Bogner-Fox-Schmit-Castro finite element with 4 nodes and 10 degrees-of-freedom (DOF) per node showed to be a fast-converging element to simulate the buckling behavior of VAFW cylinders. The model converges with approximately two orders of magnitude less DOF than a general purpose linear shell element, such as the S4R available in Abaqus.
The accelerated Kriging metamodel used to represent the buckling behavior and the mass of VAFW cylinders proved to be accurate and could be trained using only 183 finite element evaluations. This is remarkable considering the steering-thickness coupling observed in this variable stiffness structure. The moving search window, or hypercube, used during the training process enabled the proposed accelerated Kriging to outperform the classical Kriging and the co-Kriging approaches.
Particle swarm optimization was successfully used to drive the trained Kriging models towards different optima for buckling design loads of 50 kN, 100 kN, 200 kN, 500 kN and 1000 kN. All optimum results were verified against GA optimizations coupled with finite elements, showing a excellent agreement in the obtained minimum mass. Future studies will focus on manufacture the optimum cylinder configurations herein achieved through the filament winding process. In addition, the cylinders will be tested in axial compression and nonlinear models will be developed to simulate their buckling performance.

Appendix 1: Finite element model: formulation of the single-curvature BFSC
The natural coordinates are defined as = 2x∕ x − 1 and = 2y∕ y − 1 , and with the proposed nodal connectivity, the values of i , i at each node, for all elements, are: The displacements along x, y, z are respectively u, v, w and can be approximated within each element as: where S S S u,v,w i and u e u e u ei are the shape functions and the 10 degrees-of-freedom of the ith node of the SC-BFSC element, being in the following order: u, u∕ x , u∕ y , v, v∕ x , v∕ y , w, w∕ x , w∕ y , 2 w∕ x y . For the SC-BFSC element, the same shape functions of the plate BFSC element (Castro and Jansen 2021)  where x , y are respectively the finite element dimensions along x, y, as shown in Fig. 5. Using the proposed nodal connectivity for the SC-BFSC element, the nodal degreesof-freedom u e u e u ei and the respective shape functions S S S u,v,w i are concatenated as: with S S S u , S S S v and S S S w being matrices of shape 1 × 40.
Assuming equivalent single-layer theory (Barbero et al. 1990;Reddy 2003), the total potential energy functional for one finite element e can be expressed as: where gration limits x 1 ≤ x ≤ x 2 and y 1 ≤ y ≤ y 4 define the domain of one finite element e . The membrane and rotational strains are assumed to follow von Kármán kinematics, also referred to in the literature as Donnell-type (Donnell 1933, 1934 or Kirchhoff-Love nonlinear equations, given by: with (⋅), x = (⋅)∕ x used as a compact notation for partial derivatives. At the bifurcation point, the following state of equilibrium exists, considering all n e elements: Expressing the displacements within one element in terms of nodal coordinates u e u e u e as per Eq. 32, and can be calculated as: where the partial derivatives of S S S u,v,w are directly calculated from the shape functions of Eq. 33 in terms of the natural coordinates , , using the following Jacobian relations: (35) u e u e u e = { u e u e u e1 u e u e u e2 u e u e u e3 u e u e u e4 } ⊤ The neutral equilibrium criterion also requires that 2 e = 0 (Castro et al. 2014b), such that, from Eq. 38: The first integral of Eq. 41 becomes the constitutive stiffness matrix of the element, calculated using the constitutive relations from classical laminated plate theory (Reddy 2003): Note that the geometric nonlinearity appears in the constitutive stiffness matrix due to w ,x , w ,y , w ,xy in Eq. 39. Therefore, the linear constitutive stiffness matrix of a finite element K e K e K e is calculated by assuming w ,x , w ,y , w ,xy = 0 , leading to a 40 × 40 matrix: The second integral of Eq. 41 becomes the geometric stiffness matrix of the finite element K G0e K G0e K G0e , capturing the geometrically nonlinear effects of a pre-buckling membrane stresses N 0 N 0 N 0 = N 0xx , N 0yy , N 0xy ⊤ on the membrane stiffness. Noting that 2 = 0 0 0 (Castro and Jansen 2021; Castro et al. 2014b), the equation for K G0e K G0e K G0e becomes: The contributions all n e finite element are added to build the global constitutive stiffness matrix K K K and geometric stiffness matrix K G0 K G0 K G0 of the system: The pre-buckling stress field of one BFSC finite element N 0xx , N 0yy , N 0xy is calculated from the corresponding nodal displacements u 0e u 0e u 0e as: where u 0e u 0e u 0e is directly extracted from the full pre-buckling displacement vector u 0 u 0 u 0 that is obtained from a static analysis, derived from the equilibrium of Eq. 38: with f 0 f 0 f 0 representing any general pre-buckling force. When applying the neutral equilibrium criterion of Eq. 41, one assumes that at the bifurcation point there is a value of internal membrane stresses N N N based on the known pre-buckling stress N 0 N 0 N 0 described by N N N = N 0 N 0 N 0 that leads to the condition 2 = 0 . Therefore, the linear buckling problem consists of finding the value of that leads to: which holds true for any variation u u u , such that the required condition for the equality of Eq. 47 is the linear buckling equation given by:

Appendix 2: MSE estimation for the Kriging metamodels
For the unknown response function, f (x x x) , where x x x is a vector with the design variables, the response function can be approximately expressed with the Kriging metamodel (Lophaven et al. 2002) as: ] is a vector of regression function, and [v 1 , v 2 , … , v p ] is a vector of unknown coefficients.
In Eq.13, h(x x x) T v v v indicates the prediction trend, (x x x) is a Gaussian process with zero mean and covariance cov( x x x i , x x x j ) . This covariance is determined by: in which 2 is the variance of the Gaussian process and R x x x i , x x x j is the correlation function of the Gaussian process.
For the Kriging meta model with n is initial samples, x x x i , f x x x i , where i = 1, 2, … , n is , the coefficient vector v v v in Eq. 49 is calculated by where R R R is the correlation matrix whose elements are R x x x i , x x x j , i, j = 1, 2, … , n is , With the above equations, the predicted response for a new point x x x n e w can then be estimated by where: The MSE of the prediction at the new point is then calculated by

Appendix 3: Finite element model: convergence analysis and verification
The proposed SC-BFSC element is verified against the general shell element with reduced integration and linear displacement interpolation (S4R) from Abaqus FEA software (Smith 2019). To determine a converged mesh size in Abaqus, a VAFW cylinder design with 6 plies is selected, whose orientations are: plies 1 and 2  Table 11 and Fig. 21 show the convergence for both finite elements Abaqus S4R and SC-BFSC. In all cases, an aspect ratio of 1:1 was kept for the finite element meshes. Note that the third-order interpolation of the SC-BFSC allows convergence with approximately two order of magnitude less degrees-of-freedom (DOF) compared to Abaqus S4R. As expected, the first-order shear deformation theory utilized in the S4R formulation leads to more rotational flexibility and ultimately to a slightly lower buckling load. The slow convergence of S4R can be attributed to the linear interpolation of this element that requires more elements to accurately describe the buckling modes, and to the fact that only one integration point at the center of the element is used, requiring a larger number of elements when variable stiffness properties are modeled. For the SC-BFSC, a third-order interpolation is used, with 16 integration points evaluated for each element, thus leading to a fast convergence. Assuming n y = 50 for the SC-BFSC model, the buckling load is 1008.79 kN, compared to 920.38 kN for the S4R model with n y = 400 , and these are the mesh refinement levels adopted for the following discussions. After the convergence analysis, the results obtained with Abaqus and a mesh of 4R elements with n y = 400 are used to verify all optimum results obtained with a mesh of SC-BFSC elements with n y = 50 , as shown in Table 12.