Reliability-based buckling optimization with an accelerated Kriging metamodel for filament-wound variable angle tow composite cylinders

A reliability ‐ based optimization framework is introduced and used to design ﬁ lament ‐ wound cylindrical shells with variable angle tow. Seven design cases are investigated to enable a comparison between constant ‐ stiffness and variable angle tow designs, also considering effects of thickness variation created due to overlapping tow paths, determined using the kinematics of the ﬁ lament winding manufacturing process. The uncertainty in the winding angle is considered in the optimization by means of metamodels constructed using the Kriging method. Moving search windows are incorporated into the Kriging metamodel to accelerate its convergence by reducing the number of training iterations. The results prove the ef ﬁ cacy of the proposed framework and clearly demonstrate the advantage of variable ‐ stiffness designs over conventional ones for achieving a maximum load carrying capacity, while keeping the robustness of the design towards manufacturing uncertainties.


Introduction
Carbon fiber reinforced polymer (CFRP) cylindrical shells are largely utilized in space, aeronautical, marine, and energy structures, essentially due to their high capacity to sustain high levels of axial and radial compressive loads, in which most of the counterpart is under a pure membrane state [1,2]. For instance, CFRP cylindrical shells are often employed as primary structural components in space launch vehicles [3,4]. Considering that these counterparts carry high axial compression load levels, buckling is one of the limiting design constraints [5][6][7].
Among the processes to manufacture such closed shells, filament winding (FW) is well-established for high-productivity rates, highquality, and fiber volume fractions as high as 70% [8], rendering cylindrical shells for axial compression [9,10,8] and hydrostatic pressure [11,12] load carrying capacity. In fact, a reliable manufacturing process is essential to generate structures with less geometrical imperfections, which will reflect in less conservative knock-down factors [13].
Exploiting the capabilities of modern manufacturing processes is essential to reach optimal designs aiming at improved mechanical performance of the final structure. In CFRP structures, an effective way to enable larger design spaces it to explore variable stiffness designs by steering the fibers along preferential directions [14][15][16]. The expected performance improvement over conventional designs is essentially attributed to the benign load redistribution through tailoring the stiffness locally [1,17]. Variable angle tow [18] (also known as variableaxial [19] and variable-stiffness [20]) laminates have been extensively investigated for flat laminates [21,22], however less attention has been given to non-flat surfaces, especially closed cylindrical shells [23]. The first work on VAT shells in axial compression has been reported by Tatting [24], who analyzes and optimizes cylindrical shells in bending, internal pressure, torsion, and axial compression. White and Weaver [16] designed bend-free VAT shells under uniform pressure and they found that with shells compatible with the membrane hypothesis, the bending which results from a variable-radius of curvature is brought to zero by varying their orthotropy with VAT designs. White et al. [1] investigated optimal designs for buckling and post-buckling of VAT cylinders under axial compression. Hao et al. [25,26] optimize VAT composite panels and shells using an isogeometric analysis. Almeida Jr. et al. [17] model and optimize composite cylinders for axial compression load based on the characteristics of the tailored fiber placement (TFP) process. Aiming at decreasing the number of finite element analyses, Blom et al. [20] optimize VAT cylindrical shells for maximum load-carrying capability under bending using design explorer to construct surrogate models. axial buckling performance using a radial basis function (RBF) technique to build up surrogate models relating the buckling load to the local ply angle. Pitton et al. [28] optimize VAT cylinders for axial compression buckling load using an artificial neural network (ANN) aiming at approximating the buckling load to the pre-buckling stiffness of cylindrical shell. Note that among the aforementioned work on CFRP shells in axial compression, the authors [20,27,28] focused on developing surrogate models, or meta-models ("the model of the model"), in order to approximate the response of the physical system using simpler and computationally cheaper models when compared to deterministic approaches.
Even in automated manufacturing process, such as filament winding, uncertainties are present and they may arise due to several factors, such as different tow tension levels, winding angle variation, tow twist at turnaround zones and when tows are placed at low angles, i.e. more parallel to the cylinder longitudinal axis [29]. Hence, the actual performance of the component always floats around a deterministic value due to the inevitable uncertainties. Traditionally, uncertainties are incorporated into the design by means of safety-factors [30][31][32]. However, in order to explore the full potential of CFRP structures, safetyfactors should be avoided and the real stochastic effect of uncertainties should be considered. Nevertheless, mapping the uncertainties will ultimately lead to a better understanding of how manufacturing and design parameters affect the structural performance and safety [32]. It is known that load and geometric imperfections play an important role in the load carrying capacity of cylindrical shells [33][34][35][36][37][38].
The probability of achieving either the desired of acceptable performance for a counterpart under uncertainties is defined as reliability [39]. The most disseminated category to consider reliability is the well-known statistical-based Monte Carlo Simulation [40] approach, but it has a serious drawback associated to low computational efficiency when used in highly nonlinear problems and large sampling. Variance reduction methods can be utilized to enhance its efficiency. A popular variance reduction method is stratified sampling using Latin Hypercube Sampling (LHS). The LHS method makes all or nearly all sample means fall within a tiny fraction of the error. Either if one or multiple simulations need to be carried out, the estimations performed via LHS, which are unbiased, have relative deviations always lower than Monte Carlo [41]; the other category is the numerical analysisbased approach, in which common methods include first order reliability method (FORM) [42], second order reliability method (SORM) [43], radial basis functions, extreme value method [44], and Kriging [45]. These methods are, essentially, metamodels or surrogate models. Among them, Kriging metamodel is often chosen due to its high computational efficiency when under low design variables [46]. Although a few works have been developed to optimize VAT cylinders using both deterministic and probabilistic approaches, no report has been found on exploiting the manufacturing characteristics of the FW technique to generate VAT cylinders using a reliability-based optimization approach.
This work focuses on the development of a framework to optimize composite cylinders for axial buckling load allowing a VAT fiber path using a reliability-based design optimization (RBDO) approach. The manufacturing characteristics of the filament winding (FW) process are considered in the computational models and uncertainties related to variation on the winding angle are stochastically evaluated using Kriging-based metamodels. Five different design strategies are compared, consisting of one constant-stiffness and four VAT designs.

Filament-wound cylinder designs
Seven different filament-wound cylinder designs were considered in this study:   Fig. 1c; VAT-L-CT: variable stiffness with linear variation of the winding angle and constant thickness, design variables θ VL 1 ; θ VL 2 - Fig. 1d; VAT-L-VT: same as VAT-L-CT with variable thickness - Fig. 1e; VAT-P-CT: variable stiffness with second-order variation of the winding angle and constant thickness, design variables θ VP 1 ; θ VP 2 ; θ VP 3 - Fig. 1f; VAT-L-VT: same as VAT-P-CT with variable thicknesssee Fig. 1g.
All cylinders had a length L and the x-axis was used to represent the longitudinal direction. All designs were made from two filamentwound layers, or plies, with winding angle functions þθðxÞ and ÀθðxÞ, generating an angle-ply balanced laminate [47]. In all cases, θðxÞ was symmetric about the plane x ¼ L=2. Twelve tows were used to illustrate each layer in Fig. 1, with the first layer represented in blue and the second layer in green. A higher transparency was used when the tow is revolving behind the cylinder. The tows were assumed to have a constant thickness in designs CS, VAT-4, VAT-8, VAT-L-CT and VAT-P-CT. In contrast, the tows had a variable thickness in designs VAT-L-VT and VAT-P-VT. Their thickness was varied according to bending-driven variable angle tow (VAT) kinematics, explained below. Note that for VAT-L-VT and VAT-P-VT the tows do not necessarily follow a parallel path, due to the VAT kinematics combined with the kinematics of the filament winding process [48,49].
In designs VAT-4 and VAT-8 the cylinders were partitioned into four and eight frames, respectively, where a frame consists of a region of constant θ. Therefore, considering the plane of symmetry, VAT-4 had two design variables: θ V4 1 and θ V4 2 ; whereas VAT-8 had four design variables: θ V8 1 ; θ V8 2 ; θ V8 3 and θ V8 4 . The VAT-L design explored a linear variation of the winding angle along its length (x-coordinate). This was expressed as: where θ VL 1 and θ VL 2 were the design variables. Two variations of the VAT-L design were considered: VAT-L-CT had a constant thickness (Fig. 1d) whereas VAT-L-VT had a variable thickness (Fig. 1e).The design VAT-P had a second-order variation of the winding angle according to: where: The three VAT-P design variables were: θ VP 1 ; θ VP 2 and θ VP 3 ; and the interpolation points were fixed at: Again, both designs with a constant (VAT-P-CT) and variable (VAT-P-VT) thickness were considered, see Fig. 1f and g.
Variable angle tows made out of CFRP prepreg materials are mainly achieved by means of bending and shearing. When the tows are sheared, as in the continuous tow shearing process [50], the tow width is kept constant when measured normal to the steering direction, as illustrated in Fig. 2b. The tow width measured perpendicularly to the tow path changes according to w tow cos ΔθðxÞ in tow shearing.
The steering direction in the present study is the axial direction of the cylinders. Note that when shearing occurs the tow thickness must change due to conservation of mass. Castro et al. [18] derived a relation for the effective thickness h e ðxÞ, which can be adapted for FW as: where h tow is the nominal tow thickness and ΔθðxÞ represents the change in filament winding angle. The higher thicknesses created by tow shearing are illustrated in Fig. 2Àb as darker regions, and the rela- (3) to calculate the local thickness value.
Alternatively, variable angle tows can be achieved through in-plane bending [14,51], which keeps the tow width constant when measured perpendicularly to the tow path, inevitably creating residual in-plane stresses on the tows to accommodate the variable-angle, as illustrated in Fig. 2a. Ultimately, this residual in-plane stresses will determine the minimum radius of curvature achievable in VAT designs [51]. In automated fiber placement (AFP), it is customary to avoid thickness variation during fiber steering by means of cut-and-restart [20], which is not an option in the filament winding (FW) process. Castro et al. [18] demonstrated that the thickness buildup due to overlaps created by adjacent tows under in-plane bending can also be represented by Eq. (3), in a smeared approach.
For FW, it remains unclear if variable angle tows (VAT) are achieved by means of tow bending or tow shearing. Consequently, in the present study we assumed that FW achieves VAT by an unknown combination of bending and shearing. Nevertheless, this assumption does affect the thickness calculation because Eq. (3) can be used to calculate the local thickness for both VAT mechanisms, i.e. tow bending or tow shearing.
For the variable thickness designs VAT-L-VT and VAT-P-VT it was assumed that bending is the main tow steering mechanism. Therefore, the circumferential spacing between two adjacent tows (see Fig. 2) was: According to Eq. (4), when maxðθðxÞÞ is chosen to define a constant Δc applied throughout the cylinder, gaps appear in regions where θðxÞ < maxðθðxÞÞ. Conversely, if minðθðxÞÞ is used to calculate Δc, it creates overlaps where θðxÞ > minðθðxÞÞ. The design options with gaps or overlaps are illustrated in Fig. 3 for the VAT-L-VT cylinder. We

Finite element modeling
All composite cylinders considered in study had a length of 300 mm, diameter of 136 mm, layup consisting of an angle-ply layer, AEθ, nominal thickness of 0:8 mm, resulting in cylinders with a diameter-to-thickness ratio of 170. Note that the thickness for designs VAT-L-VT and VAT-P-VT ( Fig. 1e and g) deviated from this nominal value according to Eq. (3). The material properties used are listed in Table 1; these are representative of towpregs with Toray T700-12K-50C carbon fibers and a UF3369 epoxy resin.
Initial FE models were generated in Abaqus CAE [52] finite element package and the models are parametrized through scripts written in Python language. All cylinders were meshed using S4R shell elements, which is a finite-membrane-strain shell element with four nodes and reduced integration. Three integration points through the thickness were used. A reference point was created at the center of each edge and connected to the related edge through multi-point constraint (MPC). All degrees-of-freedom, that is, translations (u r , u θ , u z = 0) and rotations (ur r , ur θ , ur z = 0), were constrained to zero for the nodes at the bottom of the cylinder (x ¼ 0). The top nodes (x ¼ L) were allowed to move in the axial direction, but all other degrees-offreedom, i.e., translations (u r , u θ = 0) and rotations (ur r , ur θ , ur z = 0), were constrained to zero. A uniaxial buckling load (F z ) was applied to the reference point at the top. The converged FE mesh is shown in Fig. 4; it has 152 elements along the length and 213 around the circumference, generating a mesh with 32,376 elements and 32,589 nodes.
The FE predictions for linear buckling were based on the eigenvalue analysis using the Lanczos Eigensolver. The general buckling problem was based on the neutral equilibrium criterion of the total energy potential Π, given by: Following the derivation in Refs. [6,53], the general form can be obtained and expressed as: where K is the constitutive stiffness matrix which depends on the geometry and filament winding fiber path configuration; and K g the geometric stiffness matrix, mainly dependent on the initial stress. Only the first buckling load was extracted for the reliability-based optimization process, which is described in Section 4.

Introduction to RBDO
The actual performances of structural components are susceptible to random processes due to uncertainties in design, manufacturing and operating environment. These uncertainties can be quantified with the concept of reliability. Reliability is defined as the probability of a component to meet the desired requirements, and this can be expressed as: where R is the reliability, X is the random design variable, Pr is the probabilistic function, and GðXÞ is the performance function. A reliability-based design optimization (RBDO) is an approach that allows to optimize an objective function while considering random uncertainties [54]. Compared with traditional deterministic optimization methods, RBDO represents a substantial improvement because it satisfies both optimization objectives and reliability constraints. The general RBDO framework is defined as: where costðÞ is the cost function, d is the design variable, R jr is the jth required reliability, d L and d U are both lower and upper boundaries of the design variables, respectively. In this study, a RBDO approach was developed to consider manufacturing uncertainties in the optimization of filament-wound composite cylinders. The optimization objective was to maximize the first linear buckling load (first eigenvalue), and the design variable was the filament winding angle θ. Therefore, the RBDO model in Eq. (8) was rewritten as: where F eig ðθÞ is the first buckling load for a particular winding angle θ, whereas θ L and θ U are the lower and upper boundaries of θ, respectively. In all cases, the maximum variation in winding angle between two consecutive control points had to be less than 10 . This constraint was included to avoid abrupt changes on the winding angle and to ensure manufacturing feasibility. Therefore, the optimal design results had to respect the following constraint: In some RBDO frameworks [55], the loop optimization part is suitable for optimization problems in which the objective function can be described by mathematical expressions. However, the approaches presented in Refs. [55][56][57] are not applicable to the current optimization problem for FW cylinders. Here, the buckling load was obtained from a FE simulation (see Section 3) and this can significantly increase the computational costs if multiple design iterations are required. Therefore, a Kriging metamodel (also known as surrogate model) was introduced into the RBDO framework to increase the computational efficiency.

The Kriging metamodel
A metamodel is a mathematical function, computationally inexpensive, that approximates the output of high-fidelity and computationally intense models such as deterministic FE simulations [46]. In a metamodel-based optimization, the surrogate replaces FE simulations during the optimization process. Typical metamodels include artificial neural networks, support vector regression, polynomial response surface, support vector regression, and Kriging [32]. Among them, Kriging has demonstrated to be particularly efficient and accurate for optimization problems with a small number of design variables [20,58].
A Kriging metamodel was employed in the optimization process to alleviate the intense computational cost associated with deterministic FE simulations. The initial metamodel was constructed with samples selected using Latin Hypercube Sampling (LHS). Next, the metamodel was trained by adding new samples until a desired level of accuracy was reached. After training, the updated metamodel was used to find the optimum design. This process is illustrated in Fig. 5 and the main steps are detailed below.

Initial metamodel
The initial metamodel was created with samples selected by Latin Hypercube Sampling (LHS). The LHS approach was selected to ensure that the initial samples were uniformly distributed throughout the design space. For each design a vector of variables containing the filament winding angles was created. For instance, cylinder CS has θ CS ¼ fθ 1 g; VAT-4 has θ V4 ¼ fθ 1 ; θ 2 g; VAT-8 has θ V8 ¼ fθ 1 ; θ 2 θ 3 ; θ 4 g and so forth, such that a given sample can be represented as θ i ¼ fθ 1 ; . . . ; θ n g, for i ¼ 1; 2; . . .. If the number of initial samples is N 0 , the initial sample matrix can be written as: where each row represents one sample θ i . For each sample, the corresponding first eigenvalue F ðiÞ eig was calculated using FE simulations (see Section 3). These results were used to assemble the matrix: The initial Kriging metamodel was created using S initial [46]. The initial model is defined as E initial , and its prediction for a given sample θ i is denoted E initial ðθ i Þ. The initial model needs to be trained to obtain accurate predictions, and this procedure is explained next.

Training the metamodel
During training, new samples were added to the initial metamodel to increase its accuracy. The new samples to add to the metamodel were identified using the Mean Square Error (MSE) between the predicted and true responses. For samples at points where the true response is known, the MSE is defined as: whereŷðθ k Þ and yðθ k Þ are respectively the prediction and true responses at sample θ k . However, the new sample is added in the design region where the MSE is the largest θ MaxðMSEÞ . Therefore, an MSE estimator is used, detail of which is shown in Ref. [59] and in the Appendix. The strategy of inserting training points in high MSE regions is necessary to guarantee the robustness of the Kriging metamodel, but may lead to a large number of unnecessary iterations when the high MSE regions are far from the location of the design optimum. To alleviate this issue, an enhancement in the traditional Kriging-based sampling is proposed, where a second training point was added close to the current design optimum. This second sample was selected as the point with the highest MSE within a moving search window, defined as: where θ i opt is the optimal design result in ith metamodel, and W is the range of the window, which is determined as: where q is a scale factor. Two window sizes were used q ¼ 0:15 and q ¼ 0:25, representing 15% and 25% of the full design space, respectively. Cylinder design VAT-4 (see Fig. 1b) is taken as an example to explain the moving search window. This cylinder has two design variables: θ V4 1 and θ V4 2 , which makes it ideal to visualize the moving search window. Predictions of the initial metamodel are plotted in Fig. 6a; these are based on four initial samples and their position is shown as black points. The MSE surface is shown in Fig. 6b. The maximum MSE is easy to find; however, note that the sample with the maximum MSE is far from the optimum design region. This is a situation where  the moving search window can be particularly advantageous. The search window is shown in Fig. 7, which is a top view of the MSE surface from Fig. 6b. The black star represents the optimal result in the current metamodel and the black rectangular frame is the moving search window, determined by Eq. (14). In our approach two new training samples will be added in the next iteration, and their position is denoted with red symbols in Fig. 7. The first training sample corresponds to the maximum MSE inside the moving search window, whereas the second training sample is the maximum MSE of the complete design space. This approach is computationally efficient since it ensures that the accuracy of the metamodel is increased not only in regions of high MSE, but also near the optimum region. Note that the moving search window becomes a surface constraint in design spaces with three variables and a hyper-surface constraint when more than three variables are involved. The strategy of adding new training samples is used repeatedly until the metamodel has reached the desired level of accuracy. The procedure is stopped when two criteria (SC1 and SC2) are satisfied, and these are defined as:

SC1
The error on the optimal prediction for two consecutive cycles should be less than a given threshold c 1 . This is expressed as: where the prediction error is defined as: whereŷðθ i opt Þ and yðθ i opt Þ are the prediction value and the true response, respectively, for an optimal result θ i opt . A threshold c 1 ¼ 1% was used in all cases.
SC2 The difference between optimal predictions in two consecutive cycles should be less than a given threshold c 2 . This can be written as: A threshold c 2 ¼ 2:5% was used in this study.
If the metamodel has converged at iteration i, then the optimal prediction is given by:

Metamodel: The overall procedure
The procedure of creating and training the Kriging metamodel is demonstrated here with an example. The FW cylinder design with constant stiffness (CSsee Fig. 1a) is used for this purpose since it has a single design variable and is the simplest case considered in this study. The design variable had a range of θ CS 1 ¼ ½45:0 ; 86:6 , where the lower bound represents the minimum winding angle of the mandrel herein considered [8]. The mandrel would have to be modified, for example, by adding pins to achieve Winding angles lower than 45 [49]. Therefore, the design range considered guarantees to produce generate ready-to-manufacture fiber paths. The objective of the metamodel is to find the winding angle that maximizes the buckling load. The steps of the optimization procedure are: Step 1: Generate initial values for each variable using LHS. The number of initial samples for the Kriging metamodel is related to the number of design variables and the nonlinear degree of the objective function. For all cases herein explored, the number of design variables are known while the nonlinear degree of the response function is unknown. To ensure the efficiency of the proposed method, the minimum number of the initial samples is investigated and sought. For instance, for the CS design, which has one design variable, only two initial samples are required; For the VAT-4, VAT-L-CT and VAT-L-VT designs, all with two design variables, the number of initial samples is four; whereas the VAT-P-CT and VAT-P-VT designs with three design variables required six initial samples. For the VAT-8, which has four design variables, ten initial samples were selected. The next steps are explained based on the simplest RBDO case of the CS design, with a single design variable named θ 1 ; the two initial samples are referred to as ½θ 1 1 ; θ 2 1 ; Step 2: Evaluate the fitness functions ½F eig ðθ 1 1 Þ; F eig ðθ 2 1 Þ for the samples generated in Step 1. The fitness function is the buckling load predicted by FE simulations (see Section 3); Step 3: Build the initial metamodel with the initial samples in Step 1 and corresponding fitness evaluations; Step 4: Locate new training samples θ new . In this case, the candidate sample that corresponds to the largest MSE is selected as the new training sample. The PSO (Particle Swarm Optimization) method to locate the new training samples.
Step 5: Update the metamodel with new training samples θ new and corresponding fitness evaluations through FEM eigenvalue analysis F eig ðθ new Þ; Step 6: Check the stopping criterion SC1 and SC2, see Section 4.2.2. If both criterion are satisfied, the updating process stops and the procedure goes forward to Step 7. Otherwise, the process goes back to Step 4 to continue training the metamodel; Step 7: Obtain the optimal design and predictionŷ θ Ã opt using Eq. (19) with the corresponding θ Ã opt defined as the optimal winding angle which has the highest F eig .
The initial samples and updating process are detailed in Tables 2  and 3, respectively, for the CS cylinder. Note that the moving search window was not required for the CS design since it has a single design variable and the optimization procedure is straightforward. For each iteration of the metamodel, Table 3 summarizes: the optimal design θ opt along with its predicted (ŷ) and true (F eig ) buckling loads; the stan-

RBDO for FW cylinders under uncertainties
The proposed RBDO approach is introduced here. The target is to optimize the buckling loads for FW cylinders under stochastic variation of the winding angle under axial compression. In order to solve this problem, a RBDO framework is developed. The RBDO framework is composed of the metamodel along with a reliability analysis, shown in Eq. (9). In addition, it is valid to mention that the buckling loads are calculated through FEA. Therefore, the task turns to the reliability constraints. Similarly to the explanation presented in Section 4.2.3, the RBDO procedure is explained based on the CS cylinder to aid simplicity to the explanation. For the CS cylinder, the required buckling load F req eig is 36,885 N (value reached after a deterministic optimization), achieved for a winding angle of 50 and a required reliability, R r , of 0:9974 (based on the three-sigma concept). Then, the reliability constraint (from Eq. (9)) is determined by: where: Assuming that the metamodel is accurate enough at the optimal area, Eq. (9) can then be rewritten as: In the metamodel, both mean value and standard deviation (SD) of the buckling load at each candidate point are available considering the uncertainty from the design variables. The SD of the design variable θ is selected as 1 , and the SD of the response is obtained by the sampling method from the converged metamodel. Therefore, the reliability can be estimated from Eq. (20). Then, the RBDO procedure can be applied, whose results are presented in Table 3. These results reveal that the optimal design for the CS cylinder is obtained for an angle of θ opt ¼ 57:8 with a corresponding eigenvalue of 40; 607 N and standard deviation (SD) of 122 N. These outputs are achieved only after the reliability requirement is satisfied.

Results
The results of the RBDO analysis are presented in Table 4 for the seven FW cylinders introduced in Fig. 1. The moving search window was used for all VAT designs and the results are shown for both window sizes considered in this study. The results include the optimal winding angle(s) θ opt , the buckling load F eig and standard deviation SD of the optimal design, and the number of iterations N S needed to train the Kriging metamodel. In addition, the evolution of four key parameters of the optimization procedure are plotted as a function of the training iteration in Fig. 8 for a window size q ¼ 0:15 and in Fig. 9 for q ¼ 0:25. The optimal design for each cylinder is examined below.
The optimal VAT-4 cylinder has winding angles of ½58:7 ; 57:8 , a buckling load of 40; 474 N, and a standard deviation of 331 N when the window size q ¼ 0:25. Decreasing the size of the window to q ¼ 0:15 reduces the number of training iterations from 6 to 4, but returns essentially the same optimal design (within statistical margins). Therefore, the optimal VAT-4 design is practically identical to the optimal CS cylinder, see Table 4.

Cylinder VAT-8
The VAT-8 cylinder is characterized by four design variables: Fig. 1c). Here, ten initial samples were generated by LHS to create the Kriging metamodel.
With a window size q ¼ 0:25, the VAT-8 design has optimum winding angles of ½51:9 ; 58:6 ; 55:8 ; 61:0 for a buckling load of 40; 313 N    and standard deviation of 175 N. Again, the optimal design is fairly insensitive to the size of the moving search window. Nonetheless, reducing the size of the window decreases the number of training iterations from 10 to 6. In contrast with the optimal VAT-4 design, each frame in the optimal VAT-8 cylinder has a different winding angle. However, this tow steering strategy does not lead to significant performance benefits: the optimal VAT-8 cylinder has a buckling load comparable to the CS design.

Cylinders VAT-L
Two VAT-L configurations were considered: VAT-L-CT (Fig. 1d) assumed a constant ply thickness whereas VAT-L-VT (Fig. 1e) had a variable ply thickness. In both cases, the optimization procedure had two design variables and four initial samples were used to create the metamodel.
The optimal design for VAT-L-CT cylinder is θ VL 1 ¼ 57:1 ; θ VL 2 ¼ 57:9 with a buckling load of 40; 589 N and standard deviation 651 N. This optimum was obtained with a window size q ¼ 0:15, but the results are practically insensitive to q. The optimum VAT-L-CT cylinder is practically the same as the CS design. However, there is a surprising result in Table 4: the metamodel for VAT-L-CT converged more rapidly with q ¼ 0:25 than with q ¼ 0: 15. Considering the variation in the ply thickness leads to a considerably different optimal design. The optimum VAT-L-VT has ½θ VL 1 ¼ 48:9 ; θ VL 2 ¼ 86:2 with a buckling load of 42; 573 N and standard deviation of 676 N. This design, obtained with a search window of 25%, is 4:84% stronger than the optimal CS design. It is worth mentioning that although the difference between θ VL 1 and θ VL 2 is greater than 10 , there are no abrupt angle transitions in the design. There is a linear variation of the winding angle from θ VL 1 to θ VL 2 in the VAT-L-VT cylinders, see Eq. (1). A similar observation applies to the VAT-P cylinders, which are presented next.

Cylinders VAT-P
Two VAT-P designs were considered: VAT-P-CT (Fig. 1f) and VAT-P-VT (Fig. 1g). Both had three design variables and six initial samples were used to create the metamodel.
The optimized VAT-P-CT cylinder has a buckling load of 48; 538 N and standard deviation of 142 N, which is achieved for the following design variables: ½71:3 ; 55:8 ; 62:7 . This result was obtained with a window size q ¼ 0:25. In this case, reducing the windows size to q ¼ 0:15 lead to a considerably different optimal design with a lower buckling load, see Table 4.
In contrast, the optimal VAT-P-VT design was practically insensitive to the size of the search window. The maximum buckling load is 49; 576 N and standard deviation of 149 N for the following design variables: ½45:4 ; 86:5 ; 85:8 . Among all cases, the VAT-P designs offer the best performances: their buckling loads are about 20% higher than that of the CS cylinder. Finally, the results in Table 4 demonstrate the computational efficiency of the proposed framework: only four and five iterations were needed to train the Kriging metamodel for VAT-P-CT and VAT-P-VT, respectively.

Discussion
Comparing the performances of the cylinders using the buckling load F eig can be slightly misleading since constant thickness and variable thickness designs have a different mass m. To ensure a fair comparison, the specific buckling load F S eig ¼ F eig =m is given in Table 5. In addition, the first buckling mode is shown in Fig. 10 for the optimized CS design, and in Figs. 11 and 12 for the optimized VAT designs with q ¼ 0:15 and q ¼ 0:25, respectively. The optimized fiber paths are presented in Fig. 13 for all designs considered in this study.
• CS, VAT-4, VAT-8 and VAT-L-CT: These four optimized designs have practically the same buckling load, see Table 5. Accordingly, the optimized fibre paths are very similar as well, see Fig. 13 • VAT-P-CT with windows of 15% and 25%: The RBDO of cylinder VAT-P-CT led to an interesting conclusion regarding the use of the moving search window. When the narrower window of 15% is used, the metamodel is not able to find the optimum encountered using the wider search window of 25%. Therefore, narrow search windows should be avoided to prevent undesired early convergence of the RBDO; • VAT-L-VT and VAT-P-VT: These are the two designs that enabled variable thickness coupled with variable angle tows using Eq. (3). The optimized fiber paths consist on a unique combination of helical filaments at the edges that transition to hoop-oriented filaments at the center, see Fig. 13. This provides additional circumferential stiffening at the center of the cylinder and, consequently, the buckling modes for these designs are localized at one extremity of the cylinder, see Figs. 11 and 12.

Effectiveness of the moving search window
The moving search window was proposed to accelerate the convergence of the traditional Kriging metamodel (see Section 4.2.2). To quantify its efficiency, an additional optimization was carried out for without the search window for the VAT-4 cylinder and the results are summarized in Table 6. Without the moving search window, 16 training iterations are required for the metamodel to converge. There is a significant gain in efficiency with the search window: the number of training iterations is reduced to 7 with q ¼ 0:25 and to only 5 when q ¼ 0:15. Similar gains in efficiency were observed for the other cylinder designs. We recognize that measuring the efficiency with the number of training iterations can be slightly misleading since two new samples (instead of one) are added for each training iteration when the moving search window is used. Nonetheless, the moving search window still provide a gain in efficiency: the number of samples needed was reduced from 15 to 8 by introducing a search window q ¼ 0:15. In addition, the size of the window affects the number of iterations in the training process. For larger window sizes, half of the generated samples are less focused near the optimum, which might increase the number of iterations until convergence of the training process. Conversely, for smaller window sizes, half of the generated samples are more focused near the current optimum, improving the efficiency of the training process.

Conclusions
The present study proposes a reliability-based design optimization (RBDO) approach for improving the buckling load of variable angle tow (VAT) filament-wound cylinders subject to axial compression. Stochastic variations of the winding angle are considered in the optimization of seven different designs, ranging from constant stiffness to second-order variation of the winding angle with variable thickness. For an efficient RBDO, a Kriging-based metamodel is used, for which a new approach to accelerate the training convergence is developed, based on moving search windows. Latin Hypercube Sampling is used to initially construct the metamodel and to generate sample points for the successive stochastic evaluations. In the metamodel, new training samples are determined using the Particle Swarm Optimization algorithm. The metamodels are updated until the reliability constrains were satisfied. Furthermore, both the kinematics and manufacturing constraints of the filament winding process are taken into consideration. The novel RBDO framework shows a fast convergence of the metamodel, thereby enabling a highly computational efficient optimization for all cases. For instance, cylinder VAT-4 requires only 5 iterations to find the global optimum, whereas the traditional Kriging method needs 16 iterations. In the largest optimization design herein investigated (VAT-8), 10 iterations are needed using the accelerated Kriging approach.
Designs VAT-4, VAT-8, VAT-L-CT and VAT-L-VT present buckling performance comparable to the CS design. Designs VAT-P-CT and VAT-P-VT show an improvement of ≈20% on the specific buckling load when compared to the constant-stiffness layout; consisting of a unique combination of helical filaments at the edges that transition to hoop-oriented filaments at the center, thus providing additional circumferential stiffening at the center region of the cylinder.
Next steps of this research will be to apply experimental techniques to quantify the uncertainty of the filament winding process to be incorporated in the RBDO, such as: variability on the winding angle for different manufacturing parameters; and variable thickness pattern produced by steered tows. Moreover, the RBDO approach herein developed will be applied to future designs considering nonlinear constraints related to post-buckling and damage.

Declaration of Competing Interest
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.