Variable–angle ﬁlament–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 cylin- dersisutilizedtoexplorethenonlinearsteering–thicknessdependencyinVAFWstructures,whereby thethicknessbecomesafunctionoftheﬁberangle. Particleswarmoptimizationcoupledwithseveral Kriging-based metamodels is developed to ﬁnd the optimum designs. A single–curvature Bogner- Fox-Schmit-Castro ﬁnite element is formulated to accurately and eﬃciently represent the variable stiﬀness properties of the shells, and veriﬁcations are performed using a general–purpose plate ele- ment. Alongsidethemainoptimizationstudies,avastanalysisonthedesignspaceisperformedusing themetamodels,showingagapinthedesignspaceforthebucklingstrengththatisconﬁrmedbyge- neticalgorithmoptimizations. Extremelightweightwhilstbucklingresistantdesignsarefound,along with non–conventional optimum layouts thanks to the high degree of thickness build–up tailoring. 1. Introduction The development of new concepts for lightweight struc- tures has been overwhelmingly exploited in several ﬁelds, mainly in aeronautical and aerospace structures, to comply with Green Aviation for enhancing fuel eﬃciency and de- crease aviation emissions towards reaching carbon–neutral air transportation. A practical and direct manner to im- prove the energy eﬃciency and reduce the fuel consump- tion of an aircraft is by reducing the mass of its compo- nents [1]. Achieving as a direct consequence carbon foot- printreductionandbetterﬂightperformance. Bothaeronau- tical and space industries have been continuously develop- ing new concepts for lightweight structures to pursue these beneﬁts, where ﬁber-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 pas- sive load alleviation, proving that the use of ﬁber steering through automated ﬁber placement (AFP) manufacturing is more structurally eﬃcient than conventional straight–ﬁber– reinforced composite wings [2]. Liguorietal. [3]optimizedawingboxsectionusinggenetic <Corresponding Authors zhwanglbb@gmail.com (Z. Wang); humberto.almeida@aalto.fi (J.H.S. Almeida Jr.); A.Ashok-2@student.tudelft.nl (A. Ashok); wzhonglai@uestc.edu.cn (Z. Wang); S.G.P.Castro@tudelft.nl (S.G.P. Castro) ORCID(s): 0000-0001-8930-4817 (Z. Wang); 0000-0002-9408-7674 (J.H.S. Almeida Jr.); 0000-0002-6587-3107 (A. Ashok); 0000-0002-1708-4197 (Z. Wang); 0000-0001-9711-0991 (S.G.P. Castro)algorithms, including manufacturing constraints pertaining AFP manufacturing, achieving a wingbox section that was 6.48% lighter, with higher buckling strength and reduced out-of-planedisplacementinthepost-buckledregime,when comparedwiththeinitialstraight-ﬁberbaselinestructure. In fact, tailoring a particular property for laminated compos- itesbymeansofin–planepoint–wiseﬁberanglechangingis proventobeeﬃcient,alreadysincetheveryﬁrstreportfrom Hyer and Lee (1991) [4] on variable-angle tow (VAT) con- ﬁgurations, with the aim to optimize the buckling capacity ofsimply–supportedopen–holelaminatesinuni-axialcom- pression. Thenceforth, several approaches have been devel- oped based on the advancement of both computational ap- proachesandmanufacturingtechniquesforcompositemate- rials. Modernaircrafthavetheirfuselagesmadeofcarbonﬁber reinforcedpolymer(CFRP)composites,whichmakescylin- dricalunstiﬀenedshells[5]ofparticularimportancetowards understandingandfurtherreducingaircraftoverallmass. In spacecraft,monocoquestructuresarecommonlybuiltofun- stiﬀened shells, requiring conservative knock-down factors toaccountforthehighsensitivityofthebucklingstrengthto- wardsgeometricandloadimperfections[6]. Thefullpoten- tialofCFRPcompositesinthesecomponentsisstillunder– explored, and an eﬃcient strategy to take full advantage of these materials consists of reducing their masses whilst keeping their structural performance. The state–of-the–art on mass minimization of composite shells is scarce, as dis- cussed next, with most of the studies focusing on optimiza- tion problems that seek minimum compliance for a given Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 1 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO 1950 1970 2000 2010 2020GFRP and balsawood - Sidewalls - Bag racks - GalleysAircra� interior Boeing 737/777: - Composite ﬂoor - Horizontal/ver�cal stabilizer - Empennage1stall-composite planes CFRP: - Fuselage - Wing - 50% of compositesAirbus A350: 53% of composites Phönix FS-24Long EZBusiness jets - Variable angle tow high-aspect ra�o Wing: NASA - lighter, more structurally eﬃcient and more ﬂexibleBoeing 707 Boeing 787 Figure 1: Utilization of ﬁber–reinforced composites in aircraft over the last decades. mass. PelletierandVel[7]developedamulti-objectiveopti- mization of composite cylinders for strength, stiﬀness and minimal mass via the layerwise tailoring of ﬁber angles and ﬁber volume fraction. Although some enhancements were found, it was diﬃcult to ﬁnd a layout with reasonable compromise among the three objective functions utilized. Sadeghifar et al. [8] optimized stiﬀened cylindrical shells forminimummassandmaximumaxialbucklingloadusing ageneticalgorithm(GA),ﬁndingthatI-sectionstiﬀenersare more eﬀective than rectangular ones to achieve panels with reduced mass. Ghasemi and Hajmohammad [9] proposed a multi-objective optimization of composite shells under ex- ternal 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 eﬃciency of these optimizations is to employ metamodels, which are computationally eﬃcient surrogate modelsthatapproximatethehigh-ﬁdelityresponseusingap- proximated functions, such as hyper-surfaces. Rouhi et al. [10] optimized VAT cylindrical shells for maximum buck- ling load using radial basis functions as surrogates and ge- neticalgorithmstoﬁndtheoptimuminthesurrogatemodel. Blom et al. [11] used design explorer tools [12] as sur- rogates to minimize the buckling load of VAT composite cylinders. Recently, Wang et al. [13] developed a new ac- celerated Kriging metamodeling scheme to restrict the areaof interest within an optimization framework to maximize thebucklingloadofVATcompositecylindersmanufactured byﬁlamentwinding,designatedasvariable–angleﬁlament– wound – VAFW. No work on minimizing mass for VAT composite cylin- ders has been found, especially using manufacturing con- straintsfromﬁlamentwinding(FW),whichisthemostsuit- able and fastest manufacturing process for composite cylin- ders [14]. To fulﬁl this gap, the present work proposes an optimization approach to minimize the mass of VAFW cylinders, constrained by design loads for uni-axial buck- ling strength. Three Kriging–based metamodels are devel- oped and assessed to be used in the optimizations, namely: Kriging, Co-Kriging, and Accelerated–Kriging. The ﬁ- niteelementanalysisareperformedusingasingle-curvature Bogner-Fox-Schmit-Castro (BFSC) [15], consisting of an enriched rectangular four–node element with ten degrees– of–freedompernodetoachievethird–orderinterpolationfor the displacement ﬁeld. The FE model takes into consider- ation the manufacturing characteristics of the FW process, andcanaccuratelyrepresentthevaryingstiﬀnessproperties oftheproposeddesignparameterization. TheFEpredictions from the single–curvature BFSC model for diﬀerent design loads are veriﬁed against a general purpose shell element of a commercial FE software. A vast analysis on the de- sign space is performed using Kriging–based metamodels, showingtherangesofmassandbucklingstrengthfordiﬀer- ent number of layers. The possible mass range for diﬀerent design loads is also shown, and the optimal conﬁgurations Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 2 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO 𝜃1𝑉𝑃𝐿/4 𝜃2𝑉𝑃𝜃3𝑉𝑃𝐿/4 𝐿/4𝐿/4 Figure 2: VAFW-P parameterization with three design vari- ablesVP 1,VP 2andVP 2per angle–ply layer consisting of two plies with winding angle distribution ,.x/. are veriﬁed against a conventional genetic algorithm (GA) optimization. 2. The design The pioneering work of Wang et al. [13] on reliability- based design optimization of variable–angle ﬁlament– wound (VAFW) cylinders, compared diﬀerent VAFW de- signs and found that the maximum buckling load was achievedbythedesignnamedVAFW-P,whosewindingan- gle distribution .x/varies along the longitudinal direction usingasecond-orderLagrangepolynomial[16],asgivenin Eq. 1. Eachangle–plylayerconsistingoftwopliesofangles ,.x/is controlled by three design variables VP 1,VP 2and VP 3; whose control points are ﬁxed at x1= 0,x2=L_4, x3=L_2,x4= 3L_4,x5=L, as shown in Figure 2. In the present work, the VAFW-P parameterization is used to ﬁndoptimumdesignsofminimummassfordiﬀerentdesign loadlevels,wherebythenumberofangle–plylayersandthe winding angle distribution in each layer are the design vari- ables to be determined. .x/=T N1VP 1+N2VP 2+NL 3VP 3; xfL 2 NR 3VP 3+N4VP 2+N5VP 1; x>L 2(1) where: N1=.x*x2/.x*x3/ .x1*x2/.x1*x3/;N2=.x*x1/.x*x3/ .x2*x1/.x2*x3/ NL 3=.x*x1/.x*x2/ .x3*x1/.x3*x2/;NR 3=.x*x4/.x*x5/ .x3*x4/.x3*x5/ N4=.x*x3/.x*x5/ .x4*x3/.x4*x5/;N5=.x*x3/.x*x4/ .x5*x3/.x5*x4/(2) As explained in Wang et al [13], for VAFW cylinders the thickness distribution of the kthply given by hk.x/changes due to ﬁlament bending, ﬁlament shearing or any combina- tion of both [13], and the thickness increase is necessary to keep a constant volume [17, 18], leading to the followingnonlinear dependence of hk.x/with the winding angle dis- tribution.x/: hk.x/=htow cos.x/(3) withhtowbeingtheﬁlamentthickness;theﬁlamentsteering angleiscalculatedsimplyby .x/=.x/*min..x//,for 0fxfL. The steering-thickness coupling of Eq. 3 poses yetanotherchallengeforthedesignandoptimization,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/. The geometry of all cylin- ders herein optimized is ﬁxed, with a radius of r=0:15m, and length of L= 0:3m. The orthotropic material proper- ties of the ﬁlaments are: E11= 90GPa,E22= 7GPa, 12= 0:32,G12=G13=G23= 4:4GPa; with a ﬁla- ment thickness of htow= 0:4mm. The material density is =1600kg_m3. 3. Bogner-Fox-Schmit-Castro Finite Element with Single Curvature The Bogner-Fox-Schmit-Castro (BFSC) ﬁnite element [19, 15] is a C1 contiguous conﬁrming plate element ob- tained by taking tensor products of cubic Hermite splines. With 4 nodes per element and 10 degrees-of-freedom per node,theBFSCapproximatesthein-planeandout-of-plane displacements using 3rd–order polynomials. In the present work, a single-curvature BFSC element (SC–BFSC) is pro- posed to combine the high-order interpolation of the BFSC with cylindrical shell kinematics, with the aim of achiev- ing a computationally eﬃcient and accurate representation of the variable stiﬀness shells under investigation. Figure 3 illustrates a SC–BFSC element and a global coordinate systemxyz, where coordinate yis curvilinear following the circumferential perimeter, such that at y= 2rthe path closes on itself. The nodal connectivity is also indicated, and the mesh is built to keep the element edges parallel to thex;ycoordinates, such that lx=x2*x1=x3*x4, and ly=y4*y1=y3=y2. Figure 4 shows the mesh closing onitselfattheintersectionofelements einwithelements ei1. The natural coordinates are deﬁned as = 2x_lx*1and = 2y_ly*1, and with the proposed nodal connectivity, the values of i;iat each node, for all elements, are: Node ii 1 *1 *1 2 +1 *1 3 +1 +1 4 *1 +1(4) The displacements along x;y;zare respectively u;v;w and can be approximated within each element as: ue;ve;we=³4 i=1SSSu;v;w iueueuei(5) Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 3 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO Figure 3: Single-curvature BFSC element and the global co- ordinate system xyz. Figure 4: SC–BFSC mesh for the VAFW cylinders. whereSSSu;v;w iandueueueiare the shape functions and the 10 degrees-of-freedom of the ithnode of the SC–BFSC ele- ment, being in the following order: u,)u_)x,)u_)y,v, )v_)x,)v_)y,w,)w_)x,)w_)y,)2w_)x)y. FortheSC– BFSCelement,thesameshapefunctionsoftheplateBFSC element [15] can be used: SSSu i=HiHx iHy i0 0 0 0 0 0 0 SSSv i=0 0 0HiHx iHy i0 0 0 0 SSSw i=0 0 0 0 0 0 HiHx iHy iHxy i(6) with the cubic Hermite functions Hi,Hx i,Hy i,Hxy icalcu-lated using natural coordinates [20, 21, 22]: Hi=1 16.+i/2.i*2/.+i/2.i*2/ Hx i=*lx 32i.+i/2.i*1/.+i/2.i*2/ Hy i=*ly 32.+i/2.i*2/i.+i/2.i*1/ Hxy i=lxly 64i.+i/2.i*1/i.+i/2.i*1/(7) wherelx;lyare respectively the ﬁnite element dimensions alongx;y, as shown in Figure 3. Using the proposed nodal connectivity for the SC–BFSC element, the nodal degrees- of-freedomueueueiandtherespectiveshapefunctions SSSu;v;w iare concatenated as: ueueue= ueueue1ueueue2ueueue3ueueue4 ñ SSSu=SSSu 1SSSu 2SSSu 3SSSu 4 SSSv=SSSv 1SSSv 2SSSv 3SSSv 4 SSSw=SSSw 1SSSw 2SSSw 3SSSw 4(8) withSSSu,SSSvandSSSwbeing matrices of shape 140. Assumingequivalentsingle-layertheory[23,24],thetotal potential energy functional for one ﬁnite element ecan be expressed as: e=1 2y4 Ê y=y1x2 Ê x=x1.NNN"""+MMM/dxdy (9) wherethemembraneforcesare NNN= Nxx;Nyy;Nxy ñand thedistributedmomentsare MMM= Mxx;Myy;Mxy ñ. The integration limits x1fxfx2andy1fyfy4deﬁne the domainofoneﬁniteelement e. Themembrane """androta- tionalstrainsareassumedtofollowvonKármánkinemat- ics,alsoreferredtointheliteratureasDonnell-type[25,26] or Kirchhoﬀ-Love non-linear equations, given by: """=h n l nj"xx "yy xyi n m nk=h n l nju;x+1 2w;2 x v;y+1 rw+1 2w;2 y u;y+v;x+w;xw;yi n m nk =h n l njxx yy xyi n m nk=h n l nj*w;;xx *w;yy *2w;xyi n m nk(10) with. /;x=). /_)xused as a compact notation for partial derivatives. At the bifurcation point, the following state of equilibrium exists, considering all neelements:  =neÉ e=1 e=neÉ e=1Ê eNNNñ"""+MMMñd e=0(11) Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 4 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO Expressingthedisplacementswithinoneelementinterms of nodal coordinates ueueueas per Eq. 5, """andcan be cal- culated as: """=b f f fdSSSu ;x+w;xSSSw ;x SSSv ;y+1 rSSSw+w;ySSSw ;y SSSu ;y+SSSv ;x+w;xSSSw ;y+w;ySSSw ;xc g g geueueue =b f f fd*SSSw ;xx *SSSw ;yy *2SSSw ;xyc g g geueueue(12) wherethepartialderivativesof SSSu;v;waredirectlycalculated from the shape functions of Eq. 6 in terms of the natural coordinates;, using the following Jacobian relations: ) )x=lx 2) ) ) )y=ly 2) )(13) The neutral equilibrium criterion also requires that 2 e=0[27], such that, from Eq. 11: 2 =neÉ e=12 e=neÉ e=1L Ê eNNNñ"""+MMMñ+d e +Ê eNNNñ2"""+MMMñ2+d eM =0 (14) TheﬁrstintegralofEq. 14becomestheconstitutivestiﬀ- nessmatrixoftheelement,calculatedusingtheconstitutive relations from classical laminated plate theory [24]: NNN=AAA"""+BBB MMM=BBB"""+DDD Note that the geometric non-linearity appears in the con- stitutive stiﬀness matrix due to w;x;w;y;w;xyin Eq. 12. Therefore, the linear constitutive stiﬀness matrix of a ﬁnite elementKeKeKeis calculated by assuming w;x;w;y;w;xy= 0,leading to a4040matrix: KeKeKe=Ë xyb f f fdSSSu ;x SSSv ;y+1 rSSSw SSSu ;y+SSSv ;xc g g geñ AAAb f f fdSSSu ;x SSSv ;y+1 rSSSw SSSu ;y+SSSv ;xc g g ge +b f f fdSSSu ;x SSSv ;y+1 rSSSw SSSu ;y+SSSv ;xc g g geñ BBBb f f fd*SSSw ;xx *SSSw ;yy *2SSSw ;xyc g g ge +b f f fd*SSSw ;xx *SSSw ;yy *2SSSw ;xyc g g geñ BBBb f f fdSSSu ;x SSSv ;y+1 rSSSw SSSu ;y+SSSv ;xc g g ge +b f f fd*SSSw ;xx *SSSw ;yy *2SSSw ;xyc g g geñ DDDb f f fd*SSSw ;xx *SSSw ;yy *2SSSw ;xyc g g gedxdy(15) The second integral of Eq. 14 becomes the geometric stiﬀness matrix of the ﬁnite element KG0eKG0eKG0e, capturing the geometricallynonlineareﬀectsofapre-bucklingmembrane stressesN0N0N0= N0xx;N0yy;N0xy ñonthemembranestiﬀ- ness. Noting that 2= 000[15, 27], the equation for KG0eKG0eKG0e becomes: KG0eKG0eKG0e=Ë xyb f f fdSSSwñ ;xN0xxSSSw ;x SSSwñ ;yN0yySSSw ;y SSSwñ ;yN0xySSSw ;x+SSSwñ ;xN0xySSSw ;yc g g gedxdy(16) The contributions all neﬁnite element are added to build theglobalconstitutivestiﬀnessmatrix KKKandgeometricstiﬀ- ness matrixKG0KG0KG0of the system: KKK=ne³ e=1KeKeKe KG0KG0KG0=ne³ e=1KG0eKG0eKG0e(17) The pre-buckling stress ﬁeld of one ﬁnite element N0xx;N0yy;N0xyis calculated from the corresponding nodal displacements u0eu0eu0eas: N0N0N0=h n l njN0xx N0yy N0xyi n m nk=AAAb f f fdSSSu ;x SSSv ;y+1 rSSSw SSSu ;y+SSSv ;xc g g geu0eu0eu0e (18) whereu0eu0eu0eisdirectlyextractedfromthefullpre-bucklingdis- placement vector u0u0u0that is obtained from a static analysis, derived from the equilibrium of Eq. 11: u0u0u0=KKK*1f0f0f0 (19) Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 5 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO withf0f0f0representsanygeneralpre-bucklingforce. Whenap- pliedtheneutralequilibriumcriterionofEq. 14oneassumes thatatthebifurcationpointthereisavalueofinternalmem- branestresses NNNbasedontheknownpre-bucklingstress N0N0N0 described by NNN=N0N0N0that leads to the condition 2 =0. Therefore,theproblemconsistsofﬁndingthevalueof that leads to: uuuñKKK+KG0KG0KG0uuu=0 (20) whichholdstrueforanyvariation uuu,suchthattherequired condition for the equality of Eq. 20is: detKKK+KG0KG0KG0=0 (21) Equation21representsalinearbucklingeigenvalueprob- lem, where a characteristic polynomial in terms of can be obtainedandsolved. Forasystemwith nunknowndegrees- of-freedom, there are neigenvalues that are roots of the characteristic polynomial. In practice, the eigenvalues and correspondingbucklingmodesaresolvedusinggeneralized eigenvalue solvers that are able to eﬃciently extract only a desired number of eigenvalues and buckling modes, and in the present work the locally optimal block preconditioned conjugate gradient method [28] implemented in SciPy [29] is used. The integration of KeKeKeandKG0eKG0eKG0eover the ﬁnite element domains are performed numerically using standard Gauss- quadrature with 44integration points per element. The authors veriﬁed that this amount of integration points leads to a converged behavior. For each integration point, the lo- cal shell constitutive properties given by matrices AAA;BBB;DDD are calculated based on the winding angle distribution .x/ of all angle–ply layers. While computing the constitutive matrices, the local thickness of each ply is consistently cal- culatedaccordingtoEq. 3. Theelementsoftheconstitutive matricesarefunctionsof x,givenbyAij.x/,Bij.x/,Dij.x/, and are calculated with: Aij.x/=n³ k=1Qij.x/zk.x/*zk*1.x/ Bij.x/=n³ k=1Qij.x/1 2zk.x/2*zk*1.x/2 Dij.x/=n³ k=1Qij.x/1 3zk.x/3*zk*1.x/3(22) wherenisthenumberofplies; Qij.x/istheplystiﬀnessro- tatedbythelocalwindingangle .x/[24];zk.x/deﬁnesthe positionoftheoutwardfaceofthe kthply. Foracorrectrep- resentationoftheﬁlamentwindingprocess,theinnerfaceof theﬁlament-woundcylindermustcoincidewiththemandrel radius. Therefore,thevaluesof zk.x/areoﬀsetbyaquantity d.x/suchthatzk*1=0fortheﬁrstply,asillustratedinFig- ure 5. The correct d.x/value is half of the total thickness: d.x/=1 2nÉ k=1hk.x/ (23) a) Incorrect ply distribution b) Correct ply distributionFigure 5: a) Incorrect ply distribution with a constant mid- surface atz=0; and b) Correct ply distribution with a varying oﬀsetd.x/applied to the mid-surface. The outward face posi- tion of all plies zk.x/fork=1;2;5;6are shown. wherehk.x/is the local thickness of each ply as per Eq. 3. AppendixApresentsaveriﬁcationoftheproposedﬁniteel- ementformulationcomparedwiththemultipurposeshellel- ement S4R from Abaqus [30]. 4. The metamodels The metamodels are herein employed to decrease the computational costs of the design optimizations and to pro- vide a better understanding of the design space, and three metamodelingapproachesareconsidered,namely: Kriging, co-Kriging, and accelerated–Kriging. The metamodels are aimed at describing the hyper-surface of the buckling load and mass of ﬁlament–wound cylinders possessing diﬀerent number of layers NLand diﬀerent winding angles ,. Initially,LatinHypercubeSampling(LHS)isusedtogen- erate the initial sample space with random variables. Af- ter this initial random sampling is generated, the ﬁrst meta- model can be assembled in form of hyper-surfaces to de- scribe the surface response of both design load and mass of the cylinders. The metamodel is then trained until all con- vergence requirements are satisﬁed, which depends on the number of variables and the degree of non–linearity of the response surface. For all metamodels two convergence cri- teria,CC1andCC2, are adopted: CC1: Thepredictedresponseerrorofbothlowerandupper boundariesmustbelessthanagiventhreshold c. Then,CC1 is determined as: Errori lower<candErrori upper<c (24) whereErrori lowerandErrori upperare the errors of the lower Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 6 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO and upper boundaries in the ithtraining process, respec- tively. TheErroris then deﬁned as: Error=óó y.c/*y.c/óó y.c/(25) where y.c/andy.c/arethepredictedresponseandthetrue responseofthecandidatesample cindesignregion,respec- tively. Three thresholds care set to5~;3~;1~and applied to ﬁnd 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 ofbothbucklingloadandmassinordertoenlargethedesign space, therefore, increasing the chances of reaching higher improvements. CC2: Theerrorforbothlowerandupperboundariesintwo consecutivecyclesshouldbelessthanthegiventhreshold c, which is determined as: T Errori lower<cand Errori upper<c Errori*1 lower<cand Errori*1 upper<c(26) Thetrainingprocessstopsifthetwoconvergecriteriaare met. Then, the lower and upper boundaries of the buckling load and mass can be estimated as: TBlower=minpi*1 lower;pi lower Bupper=max pi*1 upper;pi upper (27) In the following sub-sections, each metamodel is ex- plained in detail. 4.1. Kriging For the Kriging metamodel, the spatial association be- tween the design variables and the responses is expressed as: f.xxx/=hhh.xxx/Tvvv+.xxx/ (28) where xis the predicted response,h1.xxx/;h2.xxx/;§;hn.xxx/ is a vector of regression function which is also known as the trend of the prediction based on available sam- ples,v1;v2;§;vnis a vector of coeﬃcients, and .xxx/ is a Gaussian process with zero mean and covariance covxxxi;xxxj. The initial random sample space, which is generated throughLHS,hasapopulationsizeoftentimesthenumber ofdesignvariables,asrecommendedbyForresteretal. [31]. Afterthat,thetrueresponsesforbucklingloadandmassare calculated by the FE analysis (FEA) described in section 3. Then,theinitialsamplesandresponsesareusedtoconstruct the initial Kriging metamodel [32]. A training process fol- lowstolocateandaddnewsamplesintothemetamodel,with newtrainingsamplesselectedasthesamplewithmaximumTable 1 Detailed procedure of the Kriging metamodel for an arbitrary cylinder with one ,.x/layer. Step Procedure 1: Design variables [ 1;2;3] andN=3, seti=1 2: Generate 10 initial samples ( initial) 3: Obtain the buckling load of initial samples 4: Construct initial metamodel with DACE tool, MK i*1 5: Train the metamodel 5.1: Locate sample iwith maximum MSE 5.2: Obtain the buckling load iof i 5.3: Add iandito modelMK i*1 6: Estimate the boundaries i lowerandi upper 6.1: If CC1andCC2are satisﬁed, go to Step 7 6.2: Else, i=i+1, go back to Step 5 7: Obtain the boundaries with Eq. 27. mean square error (MSE) of the prediction, deﬁned as fol- lows: MSE=1 nnÉ k=1. y.k/*y.k//2(29) where y.k/andy.k/are the prediction and true responses for the sample k, respectively. Thetrainingprocessisrepeateduntiltheconvergencecri- teria (Eqs. 24, 26) are satisﬁed. In addition, the lower and upperboundariesoftheresponseinthedesignregionarede- terminedbyEq. 27. Adetaileddescriptionoftheprocedure is given in Table 1. 4.2. Co-Kriging Thesecondmetamodelexploitedhereusestheco-Kriging principles, which is also employed to predict the buckling loadandmassoftheVAFW-Pcylinders. Theco-Krigingis, essentially, an the extension of the Kriging approach, com- monly used to approximate a surface or hyper-surface of a system response that is computationally expensive to esti- mate while the response of the system is dependent on a set oftraineddata,becomingcomputationallyinexpensive[33]. In co-Kriging, the metamodel is built with two sets of data, one set named expensive data .e;ye/that are more computationally costly, while the other set is called cheap data.c;yc/. Thesetwosetsofdataareindependentofeach other. The response yeateis supposed to be more accu- rate than the cheap data ycatc. The required total data are combined as: r=0 e c1 = .1/ e;§;.ne/ e;.1/ c;§;.nc/ cT (30) whererare the required samples of random variables; e andcare random samples from expensive and cheap data, respectively; neandncare the number of expensive and Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 7 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO cheap data, respectively. In this work, the number of data neandncare deﬁned as ne=2n;nc=10n (31) wherenis the number of random variables. The responses corresponding to the random samples r are shown in Eq. 32, and the combined data .r;yr/are the required data to build the co-Kriging metamodel. yr= ye .1/ e ;§;ye .ne/ e ;yc .1/ c ;§;yc .nc/ cT (32) SimilarlytoKrigingapproach,co-Krigingrealizesthere- sponse surface with the required data through Eqs. 30,32. Withthecheapdata .c;yc/,theresponseofthecheapfunc- tionisdescribedasaGaussianprocess fc./inthetotalsam- ple space [32]. Then, the expensive response surface is de- ﬁned as: fe./=fc./+fd./ (33) wherefe./is the prediction at ,is a scaling factor esti- mated by the data, and fd./is the diﬀerence between ex- pensive and cheap responses. The initial independent cheap cand expensive esam- ples of random variables are generated via LHS. The buck- ling loadscof samples care obtained from the Kriging metamodel (Section4.1), and the buckling loads eof sam- ples eare calculated by FEA. Then, the metamodel can be built with the initial required data: . c;c/and. e;e/. Now that the core of the co-Kriging metamodel is set, the training process can take place. In each iteration of the training,threenewtrainingsamplesareaddedtothecurrent metamodel. During the ithiteration (i=1;2;:::;n), the ﬁrst trainingsampleisthepointatwheretheMSEisthehighest in the whole design region, expressed as i. The buckling loadiof the sample iis estimated FEA. The second and thirdnewtrainingsamplesaredeterminedasthesamplesat thecurrentlowerandupperboundariesofthe ithmetamodel, whichareexpressedas . L i;L i/and. U i;U i/,respectively. Oncethenewtrainingsamplesareavailable,theyareadded to theithmetamodel to improve the accuracy of the predic- tions. The ﬁrst sample . i;i/in theithiteration is added to the expensive data since it is likely to be accurate. The expensive data in the .i+1/thiteration is then updated as:  i+1 e;i+1 e= i e; i;i e;i(34) where i+1 e;i+1 eistheexpensivedatainthe .i+1/thiter- ation, and i e;i eis the expensive data in the ithiteration. 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/thiteration is then updated as:  i+1 c;i+1 c= i c; L i; U i;i C;L i;U i(35)Table 2 Detailed procedure of the co-Kriging metamodel for an arbi- trary cylinder with one ,.x/layer. Step Procedure 1: Design variables [ 1;2;3] andN=3, seti=1 2: Generate 10 initial samples ( initial) 3: Obtain the buckling load of initial samples 4: Construct initial co-Kriging metamodel MK i*1 5: Train the metamodel 5.1: Locate sample iwith maximum MSE 5.2: Locate sample L iand U iof the boundaries of MC i*1 5.3: Obtain buckling loads i,L i,U iof i, L i, U i 5.4: Add [ i,i] to samples [ e,e] 5.5: Add [ L i,L i] and [ U i,U i] to samples [ c,c] 5.6: Obtain new model MC i 6: Estimate the boundaries i lowerandi upperofMC i 6.1: If CC1andCC2are satisﬁed, go to Step 7 6.2: Else, i=i+1, go back to Step 5 7: Obtain the boundaries with Eq. 27. where i+1 c;i+1 cisthecheapdatainthe .i+1/thiteration, and i c;i cis the cheap data in the .i/thiteration. After the updated expensive and cheap data are avail- able, the metamodel is then trained. The lower and up- perboundaries . L i+1;L i+1/and. U i+1;U i+1/areestimatedby the extreme value optimization algorithm. In this work, the ParticleSwarmOptimization(PSO)algorithmisemployed, whichisexplainedinsection5. Thetrainingprocessevolves until both convergence criteria CC1andCC2are met. Fi- nally, the lower and upper boundaries of the buckling load in the design space can be calculated by Eq. 27. Table 2 summarizes the procedure to ﬁnd the buckling load boundaries through co-Kriging approach, in which the VAFW-Pdesignwithone ,.x/layeristakenasamereex- ample. 4.3. Accelerated Kriging During the approximation of the hyper-surface for the buckling load and mass carried out through both Kriging andCo-Krigingmetamodels,considerableeﬃciencyandac- curacy have been reached. However, their performance is signiﬁcantly aﬀected by the number of variables in the sys- tem and the level of non-linearity of the system response, especiallywithregardstothenumberofiterationsneededto train the metamodels for the required convergence criteria. With the aim to overcome this issue, an accelerated Krig- ing metamodel is introduced, which has been proven to be computationally eﬃcient, accurate and improved in terms of convergence to optimize the buckling load of VAFW-P cylinders[13]. Withanincreasingnumberofvariablesoran increaseinthenon-linearityoftheresponse,thereisalarge possibility that the sample with largest MSE is located at a design space region that is too far away from the optimum region, thus reﬂecting in a slow convergence of the conven- tional Kriging and Co-Kriging methods during the training phase. The proposed accelerated Kriging metamodel is an Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 8 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO extension of the conventional Kriging approach, whereby a moving search window, or hyper-cube, is used during the training phase, allowing the addition of a new sample near the current optimum, speeding up the convergence whereas keeping the robustness of the original method. The training process starts with initial random samples generated via LHS, in which the number of samples is the same used in the conventional Kriging metamodeling ex- plained in Section 4.1. The responses of these samples are estimated by FEA. After the initial samples and the re- sponsesarecomputed,theinitialmetamodeliscreatedusing theDACEtoolboxinMatlab[32]. Next,foreachtrainingit- eration two points are added, being the ﬁrst training point theonewithlargestMSE,similarlytotheothertwoKriging metamodelspreviouslydescribed. Thesecondtrainingpoint of the accelerated Kriging approach is added based on the largestMSEwithinanadaptivemovingsearchwindowcen- tered at the current optimum. The adaptive moving search window is deﬁned as: L i+1;U i+1= i opt*0:5W;i opt+0:5W (36) whereL i+1andU i+1are the lower and upper boundaries of the.i+1/thadaptive window, respectively; The i optis the ithoptimal design in the ithiteration of the training proce- dure;Wis the range of the window, which is determined as: W=q L+U(37) whereLandUare respectively the lower and upper boundaries of the design space, while qis a scale factor, which ranges from 0 to 1. Considering that the center of the window is the previ- ous optimal result (Eq. 36), the search window adapts itself bymovingduringthetrainingiterations. Inthiswork,acon- stantwindowwidthof q=0:25isused. Thelowerandupper boundariesofthe ithadaptivemovingwindowisreplacedby thecorrespondinglowerandupperboundariesofthedesign space if the ithwindow moves out of the design space, as follows: <L i=L;ifL i<L U i=U;ifU i>U (38) With the ﬁrst and second training samples determined, thecorrespondingresponsesareestimatedbyFEA.Thenew samplesandresponsesarethenaddedtothemetamodel. The iterations are repeated until the convergence criteria CC1 andCC2are met. Table 3 describes a summary of the pro- posed approach for a VAFW-P cylinder with one ,layer. 4.4. Design space analysis for buckling loads After the convergence criteria ( CC1,CC2) and thresh- olds(5~,3~,1~)aremetforallthreeKriging-basedmeta- models, the lower and upper boundaries for the bucklingTable 3 Detailed procedure of the accelerated-Kriging metamodel for an arbitrary cylinder with one ,.x/layer. Step Procedure 1: Design variables [ 1;2;3] andN=3, seti=1 2: Generate 10 initial samples ( initial) 3: Obtain buckling load of initial samples 4: Construct initial Kriging metamodel MA i*1 5: Train the metamodel 5.1: Locate sample iwith maximum MSE 5.2: Locate sample W iwith the adaptive moving window 5.3: Obtain buckling loads i,W iof i, W i 5.4: Add [ i,i] and [ w i,w i] to the training samples 5.5: Obtain new model MA i 6: Estimate the boundaries i lowerandi upperofMA i 6.1: If CC1andCC2are satisﬁed, go to Step 7 6.2: Else, i=i+1, go back to Step 5 7: Obtain the boundaries with Eq. 27. loads are estimated. Table 4 shows the obtained boundaries and the required number of samples Nsused in the training process for the VAFW-P design with one layer ( NL= 1). It is observed that the accelerated-Kriging metamodel pro- vides the widest ranges in terms of buckling loads for all thresholds applied. Focusing speciﬁcally on the threshold c= 1~, which provides wider boundaries, the more the number of layers, the higher the number of new samples (NS) 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 higheraccuracy,itsresultsareusedtoestimatethebuckling loads and mass for NL=2;3;4. Table6summarizestheestimatedlowerandupperbound- ariesofbucklingloadfortheVAFW-Pcylinderswith NL= 1;2;3;4,obtainedwiththeaccelerated-Krigingmetamodel. Notethatthisisakeytaskaccomplishedbythemetamodel, resulting in a better understanding of the design space. 4.5. Mass metamodeling The mass of the cylinders is calculated as: m=Ëxzh.x;z/rdxdz (39) wheremis the mass, is the density of the carbon/epoxy composite material ( 1600kg_m3),xis the circumferential coordinate,zis the longitudinal coordinate, ris the radius of the cylinder, and h.x;z/is the thickness. To calculate the mass for the cylinders with several lay- ers, a simple function to sum of the mass of every layer is considered: M.;NL/=NLÉ i=1mi.i/;NL=1;2;3;4 (40) Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 9 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO Table 4 Converged lower and upper boundaries of buckling loads for VAFW-P with NL=1for the three Kriging–based metamodels. Kriging Co-Kriging Accelerated-Kriging c=5~c=3~c=1~c=5~c=3~c=1~c=5~c=3~c=1~ 2 L[N]34,534 34,534 34,304 33,955 33,955 33,615 34,842 34,842 33,729 2 U[N]65,914 65,914 66,160 66,733 66,733 67,901 67,946 67,946 68,768 Ns2 2 4 6 6 21 6 6 30 Table 5 Number of required samples to train all three Kriging–based metamodels for a threshold of c=1~. MetamodelNumber of layers ( NL). 1 2 3 4 NS: Kriging 2 25 85 102 NS: Co-Kriging 7 51 63 126 NS: Accelerated-Kriging 21 81 108 183 Table 6 Converged boundaries of buckling loads for the VAFW-P de- sign using the accelerated-Kriging metamodel. BoundariesNumber of layers ( NL) 1 2 3 4 NL L[N]33,729 146,190 329,563 694,931 NL U[N]68,768 1,852,980 2,909,814 3,496,871 Table 7 Lower and upper mass limits for all cylinder designs. BoundariesNumber of layers ( NL) 1 2 3 4 MNL L[kg]0.3644 0.7288 1.0932 1.4576 MNL U[kg]6.9796 13.9592 20.9388 27.9184 wheremiis the mass of the ithlayer, andiis the winding angles in the ithlayer. Once the boundaries of the mass for two layers is computed, the boundaries for the other cylin- derscanbecalculatedbyEq. 40. Then,thelowerandupper limits for all cylinders are listed in in Table 7. Figure 6 illustrates the hyper–surfaces of both buck- ling loads and mass for the cylinders trained with the accelerated–Kriging metamodel using Monte Carlo sam- pling (MCS). Considering that there are three design vari- ables, one of them is kept constant at 45ýto allow a con- venient graphic representation. Then, after both metamod- elsforbucklingloadsandmassaretrained,theoptimization follows, as detailed next. Figure7illustratestheoutcomeofthetrainedmetamodels by plotting the buckling load versus the mass for several samplesizesbyusingMCSapproach. Throughouttheopti- mization procedure the samples are dynamically generated, as explained in section 5. From Figure 7 there is a clear trendofhigherbucklingloadsforhighermassvalues,andit can be observed a complex mass distribution created by the coupled steering-thickness coupling of the VAFW-P cylin- (a) Buckling load surface, kN (b) Mass surface, kg Figure 6: Contour plots using the accelerated-Kriging meta- model with NL= 1(one,layer). For plotting purposes, assuming3=45ý. der design. From the sampling with 2000 samples of Fig- ure 7(f), the minimum and maximum intervals of mass and buckling loads were identiﬁed for NL=1;2;3;4, as shown inFigure8. Notethatfordesignloadsbelow ù69kNonly designswith NL=1shouldbeconsidered,whereasforade- sign load ofù700kN, all designs with NL=2;3;4might Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 10 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO lead to optimum solutions and should therefore be consid- ered by the optimization algorithms. Figure 8 also shows a gap in the design space between ù69kNandù146kN. 5. The optimization 5.1. Particle Swarm Optimization (PSO) The PSO is a powerful meta-heuristic technique inspired ontheswarmbehaviorofﬂocksofbirdsandshoalsofﬁshes, and developed by Kennedy and Eberhart in 1995 [34]. The methodisalsorelatedtoevolutionarycomputationandithas sometiestogeneticalgorithm(GA).Thekeyadvantagesof PSO are that it needs only primitive mathematical opera- tions and it is computationally cheap in terms of memory and speed requirements. In PSO, a swarm with particles Pis placed in the design search space. Each ithparticle in the swarm has a posi- tionvectorXt i=xi1;xi2;xi3;§;xinTandvelocityvector Vt i=vi1;vi2;vi3;§;vinTat iterationt. Each particle is updated according to the dimension j, as per Eq. 41. Ev- ery particle is a potential candidate solution and they move continuouslyaroundthedesignspaceuntiltheoptimization meets the convergence criteria. The new particle locations areupdatedateveryiterationtowardsﬁndingtheglobalop- timum. The velocity of the particles can be deﬁned as: Vt+1 ij=wVt ij+c1rt 1 pbestij*Xt ij +c2rt 2 gbestj*Xt ij (41) whereas the new position is calculated as: Xt+1 ij=Xt ij+Vt+1 ij(42) inwhichVd iXd iarethevelocityandpositionofthe ithparti- cle,respectively; c1andc2areaccelerationconstants; rt 1and r2 2represent two random numbers distributed in the range (0,1);pbestandgbestare the best position so far of the ith particleandthebestglobalposition,respectively; wisthein- ertiaweight,whichkeepsabalancebetweenlocalandglobal searches. Once the buckling load and mass hyper-surfaces are con- structed (Section 4), the PSO procedure takes place, whose optimization problem is formulated as: Find:VP i;L;NL;fori=1;2;3;and:L=1;5;NL Minimize:M.VP i;L;NL/ Subjectto:.VP i;L;NL/>d LfVP i;LfU;forL=3:3ýandL=87:7ý (43) whereVP i;Lare the angle at the control points for each layer L;LandUarethelowerandupperboundariesfortheﬁl- ament angle; and NLthe total number of layers to be deter- mined. Note that, from the sampling previously performedTable 8 Parameters used in the PSO–MM optimization. Parameter Value Number of variables ( NL) 1, 2, 3, 4 Minimum inertia weight ( Wmin) 0.4 Maximum inertia weight ( Wmax) 0.9 Acceleration factors ( c1=c2) 2 Population size 20 Maximum iterations 100 Number of swarms 1 and shown in Figure 7, it is possible to determine all pos- sible values of NLfor a given design load, as illustrated in Figure 8, and for the region with a gap in the design space fromù69kNtoù146kNit was assumed NL=2. TheparametersutilizedinthePSO–metamodeling(PSO– MM) optimization framework are presented in Table 8. In each optimization run shown in Eq. 43, the particle withminimummassissearchedfor5desiredbucklingloads: 50kN,100kN,200kN,500kN,and1000kN.Then,foreach maximumdesiredload,theoptimizationﬁndsbothoptimum anglesoptand number of layers N. Theoverallmassminimizationproceduretomeetdesired buckling loads BLdcan be summarized as follows: S1:For any given desired buckling load, the number of layers is determined with the boundary analysis of the bucklingloadshowninTable6. Thecandidatenumber of layersNis determined as N=[2i;§;2j] (44) in which, i;jË[1;§;4];whereifj (45) S2:Setn=i; S3:The number of layers is then N=2nand the number of design variables ( opt n) is3n. The optimization ﬂowchart is illustrated in Figure 9, starting by calling both buckling load metamodels in eachiteration. Then,thePSOalgorithmﬁndsthemini- mummass,Mmin n,withinthebucklingloadboundaries, as detailed in the following steps: S4:Ifn<j,n=n+1, go back to Step3. Else, go to S5; S5:Obtain all the candidate optimal designs and corre- sponding masses: opt ;Mmin ; =[1;§;n] (46) S6:Selectthesamplewiththeminimummassofallcandi- dates using Eq. 46 as the optimum for a given design loadBLd. Then, the ﬁnal output, opt final, is generated; Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 11 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO (a) 200 samples (b) 400 samples (c) 600 samples (d) 800 samples (e) 1000 samples (f) 2000 samples Figure 7: Design space illustrated using the buckling load versus mass for diﬀerent samples generated using the developed metamodels. Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 12 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO Mass (kg) 0.36446.9796 0.728813.9592 1.093220.938827.9184 1.4576 NL = 1 NL = 4 N L = 3 N L = 2 33,72968,768146,1901,852,980 329,5632,909,8143,496,871 694,931Buckling load (N) NL = 1 NL = 4 N L = 3 N L = 2(a) Mass boundaries (b) Buckling load boundaries Figure 8: Mass and buckling load boundaries deﬁned by the metamodels using 5000 samples. 5.2. Optimization veriﬁcation: genetic algorithm (GA) Optimizationsusingageneticalgorithm(GA)alongwith thedevelopedSC–BFSCFEmodelareusedtoverifytheop- timization results from the PSO–MM framework. The GA is implemented in OpenMDAO, an open–source optimiza- tionframework[35]. Itisworthmentioningthatthenumber oflayersallowedarethesameasinthePSO–MMoptimiza- tions, that is, NL= 1;NL= 2;NL= 3;NL= 4. The parametersutilizedintheoptimizationsarepresentedinTa- ble 9. The objective function is similar to Eq. 43, whereby the function is changed in accordance with the GA problem deﬁnition as shown in Eq. 47: Find:;NL Minimize:M.;NL/PF Penalty Function .PF/:T 1 ifgd .1+d_/2if<d Design Space: [0;90] (47)Table 9 Parameters used in the GA optimization. Parameter Value Bit size (Bs) 6 Maximum generations 100 Population size 8Bs Mutation rate 0.02 whereMisthemasscalculatedaccordingtoEq. 40; isthe bucklingloadoftheindividual;and disthedesignbuckling load. 5.3. Optimization results The objective function of the optimizations is to mini- mizethecylindermasssubjecttothefollowingdesignloads: 50kN,100kN,200kN,500kN,and1000kN.ThePSO-MM resultsconsistofoptimumdesignsobtainedperdesignload, where diﬀerent optima are obtained for diﬀerent number of ,layers in the ﬁlament winding process. The GA 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 opti- mumdesigns,andthediﬀerencesbetweenthePSO-MMand GA results are identiﬁed: •50 kN: this design load is covered by candidates with NL= 1, as illustrated in Figure 8. The GA results is slightly heavier than the PSO-MM result; •100 kN: from Figure 8 there is no design covering this design load, and this region represents part of a gap in the design space from ù 69kNtoù 146kN. The PSO-MMoptimizationfoundacandidatecylinderwith NL= 2, with a buckling load of 165.33kN, expect- edlyhigherthanthe100kNdesignload,duetothegap inthedesignspacepreviouslyidentiﬁed. Interestingly, theGAresultconﬁrmedtheidentiﬁedgapbyﬁndingan optimum solution outside the gap region, with a buck- lingloadof288.15kN,andacomparablemasswithre- specttotheoptimumdesignobtainedbythePSO-MM; •200 kN: the only candidate for this design load level are those with NL= 2. The combination of angles leads to shell with a mass of 0.7309kg and buckling loadof207.61kN.TheGA,nevertheless,doesnotﬁnd anoptimumwithbucklingloadaroundthedesiredone. However,theoptimumfoundbyGAhasasimilarmass withasubstantiallyhigherbucklingloadof242.83kN; •500 kN: as per Figure 8 this design load is covered by cylinders with NL= 2andNL= 3, and the PSO- MM optima have respective masses of 1.2715kg and 1.0951kg; with buckling loads of 498.29kN and 504.28kN. The GA ﬁnds a unique optimum candi- date withNL= 3, with optimum mass comparable to the PSO-MM result, and a higher buckling load of 659.80kN; Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 13 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO Obtain the responses of samples Start End Locate new training samples and add them to previous metamodel Generate the initial samples with Latin Hypercube Sampling Build the initial metamodel Update the metamodel Criterion satisfied? YesNo Metamodel is trained Build metamodels for buckling load and mass Start End Desired buckling load with Ni layers Boundaries of buckling load and mass with Ni layers Candidate cylinders with Nj layers Obtain Ropt with minimum mass Optimization of each candidate cylinder - PSO Results: R N1,..., R Nj Optimization Metamodeling Figure 9: Flowchart of the proposed PSO–MM framework. •1000 kN: three candidates are found by the PSO–MM to respect this desired load, with 2, 3, and 4 ,layers. Theyhavemassesbetween1.46and1.61kgwithbuck- ling loads within a range of 1021.5 and 1092.2kN, re- spectively. TheGAﬁndsonlyonecandidatewith4 , layers, which has similar mass to the candidate with 4 layersfoundbyPSO–MMframework,with1.46kgbut with a higher buckling load of 1301.7kN. The GA optima reached considerably higher buckling loads without resulting in additional mass for the design loads of 200kN, 500kN and 1000kN, shown in Table 10. ThisisnotduetotheobjectivefunctionofEq. 47utilizedin the GA optimizations, where the buckling load only par- ticipatesinthepenaltyfunctionwhen <d,withdbeing the design load.6. Discussion The PSO-MM optimum designs listed in Table 10 are il- lustrated in Figures 10 – 17. The illustrated layers kiare stacked in an increasing zdirection while constructing the ﬁlament–wound structure, as illustrated in Figure 5. The correspondingﬁrstlinearbucklingmodesareshowninFig- ure 18. When VP 1;VP 2;VP 3are diﬀerent within the same ,layer, a variable–angle ﬁlament 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. [13], the variable–angle path creates overlapping re- gionswithhigherthicknessbetweentheﬁlaments,hererep- resentedbydarkerregions. Thiscouplingbetweenthethick- nessandthesteeringangleleadstointerestingoptimumde- signs,suchastheonesillustratedinFigure13,Figure15and Figure 16. Fordesignloadsthatallowanoptimumconﬁgurationfor Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 14 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO Table 10 Optimization results from both PSO–MM and GA approaches. Approach Design Load (kN) knVP 1.ý/VP 2.ý/VP 3.ý/Mass (kg) (kN) PSO–MM50k1,46:3 ,57:8 ,57:80.3704 50.31 GA k1,41:8 ,68:0 ,49:90.3888 49.99 PSO–MM 100k1,60:9 ,60:0 ,64:80.7291 165.33k2,41:3 ,41:6 ,39:7 GAk1,22:9 ,22:9 ,21:40.7238 288.15k2,40:0 ,40:0 ,40:0 PSO–MM 200k1,6:8 ,5:2 ,6:40.7309 207.61k2,31:0 ,38:2 ,37:0 GAk1,28:3 ,28:3 ,27:60.7238 242.83k2,75:1 ,75:1 ,75:1 PSO–MM 500k1,10:5 ,56:5 ,10:61.2715 498.29k2,9:6 ,80:7 ,58:3 k1,54:2 ,53:9 ,53:4 1.0951 504.28 k2,59:8 ,65:7 ,63:7 k3,45:9 ,46:9 ,43:0 GAk1,49:6 ,48:9 ,48:2 1.0861 659.80 k2,83:6 ,85:0 ,85:0 k3,39:0 ,39:0 ,45:4 PSO–MM 1000k1,72:3 ,67:5 ,6:61.6112 1012.51k2,7:0 ,63:3 ,87:5 k1,62:3 ,36:0 ,52:2 1.5877 1008.79 k2,18:0 ,74:4 ,87:1 k3,76:5 ,85:6 ,34:4 k1,37:3 ,37:0 ,36:2 1.4586 1092.18k2,62:1 ,65:7 ,65:9 k3,5:9 ,7:0 ,5:8 k4,31:3 ,33:0 ,35:7 GAk1,21:3 ,22:7 ,19:1 1.4487 1301.76k2,53:9 ,48:9 ,51:0 k3,79:4 ,78:0 ,78:0 k4,22:7 ,19:8 ,17:7 𝐿/4𝐿 𝐿/23L/4 Figure 10: PSO-MM tow path for design load of 50 kN. diﬀerent values of the total number of layers NL, a general trendisobserved,inwhichlower NLvaluescanachievethe design load by means of more thickness build–up created by the variable–angle ﬁlaments. On the other hand, higher NLvalues tend to show less thickness build–up, thus with less variability in the ﬁlament angles. Another trend that is observed is that the buckling modes for the conﬁgurations with more thickness build–up tend to show the well-known diamond patterns, whereas the optima with more NLhave more axisymmetric buckling modes. 𝐿/4𝐿 𝐿/23L/4 Figure 11: PSO-MM tow path for design load of 100 kN with 2,layers. 7. Conclusion Thisworkfocusedondevelopanoptimizationframework coupledwithmetamodelstominimizethemassofvariable– angle ﬁlament–wound (VAFW) cylinders. The newly pro- Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 15 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO 𝐿/4𝐿 𝐿/23L/4 Figure 12: PSO-MM tow path for design load of 200 kN with 2,layers. 𝐿/4𝐿 𝐿/23L/4 Figure 13: PSO-MM tow path for design load of 500 kN with 2,layers. posedsingle–curvatureBogner-Fox-Schmit-Castroﬁniteel- ement 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 con- verges 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 ﬁnite element evaluations. This is remarkable considering the steering–thickness coupling observed in this variable– stiﬀness structure. The moving search window, or hyper- cube,usedduringthetrainingprocessenabledtheproposed acceleratedKrigingtooutperformtheconventionalKriging and the co-Kriging approaches. Particle swarm optimization was successfully used to drive the trained Kriging models towards diﬀerent optima forbucklingdesignloadsof50kN,100kN,200kN,500kN and 1000kN. All optimum results were veriﬁed against ge- netic algorithm optimizations coupled with ﬁnite elements, showing a very good agreement in the obtained minimum 𝐿/4𝐿 𝐿/23L/4 Figure 14: PSO-MM tow path for design load of 500 kN with 3,layers. 𝐿/4𝐿 𝐿/23L/4 Figure 15: PSO-MM tow path for design load of 1000 kN with 2,layers. mass. The GA runs required above 60 generations and a to- tal of 3841 ﬁnite element evaluations to converge. Future studies will focus on manufacture the optimum cylinderconﬁgurationshereinachievedthroughtheﬁlament winding process. In addition, the cylinders will be tested in axialcompressionandnon–linearmodelswillbedeveloped to simulate their buckling performance. Acknowledgements Z. Wang is ﬁnancially supported by the Chinese Scholar- ship Council (CSC); J.H.S. Almeida Jr. is supported by the Royal Academy of Engineering under the Research Fellow- ship scheme [Grant No. RF/201920/19/150]. Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 16 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO 𝐿/4𝐿 𝐿/23L/4 Figure 16: PSO-MM tow path for design load of 1000 kN with 3,layers. Conﬂict of interest The authors declare that they have no conﬂict of interest. A. FE model: Convergence analysis and veriﬁcation The proposed SC–BFSC element is veriﬁed against the general shell element with reduced integration and linear displacement interpolation (S4R) from Abaqus FEA soft- ware [30]. To determine a converged mesh size in Abaqus, aVAFWcylinderdesignwith6pliesisselected,whoseori- entations are: plies 1 and 2 VP 1=,62:3ý,VP 2=,36:0ý, VP 3=,52:2ý;plies3and4 VP 1=,18:0ý,VP 2=,74:4ý, VP 3= ,87:1ý; and plies 5 and 6 VP 1= ,76:5ý,VP 2= ,85:6ý,VP 3=,34:4ý. Table 11 and Figure 19 show the convergence for both Abaqus S4R and the SC–BFSC element. In all cases an as- pectratioof1:1waskeptfortheﬁniteelementmeshes. Note that3rd-order interpolation of the SC–BFSC allows con- vergence with approximately two order of magnitude less degrees-of-freedom (DOF) compared to Abaqus S4R ele- ment. As expected, the ﬁrst-order shear deformation theory utilizedintheS4Rformulationleadstomorerotationalﬂex- ibility and ultimately to a slightly lower buckling load. The slowconvergenceofS4Rcanbeattributedtothelinearinter- polationofthiselementthatrequiresmoreelementstoaccu- ratelydescribethebucklingmodes,andtothefactthatonly oneintegrationpointatthecenteroftheelementisused,re- quiring a larger number of elements when variable-stiﬀness propertiesaremodelled. FortheSC–BFSC,athird-orderin- terpolation is used, with 16 integration points evaluated for each element, thus leading to a fast convergence. Assum- 𝐿/4𝐿 𝐿/23L/4 Figure 17: PSO-MM tow path for design load of 1000 kN with 4,layers. ingny= 50for the SC–BFSC model, the buckling load is 1008.79kN,comparedto920.38kNfortheS4Rmodelwith ny=400, and these are the mesh reﬁnement levels adopted for the following discussions. After the convergence analysis, results obtained with a meshofAbaqusS4Relementswith ny=400areusedtover- ify all optimum results obtained with a mesh of SC–BFSC elements with ny=50, as shown in Figure 12. Replication of Results TheAbaqusmodelsusedforveriﬁcationaregeneratedus- ing Python scripts. All single–curvature BFSC models are also created using Python scripts. The Kriging metamodels andparticleswarmoptimizationisrunusingMatlabscripts. The full datasets, as well as the developed codes, can be madeavailablewithareasonablerequesttoSaulloG.P.Cas- tro. CRediT authorship contribution statement Zhihua Wang: Methodology, formal analysis, investi- gation, software, writing - original draft preparation. José Humberto S. Almeida Jr.: Methodology, formal analy- sis, validation, investigation, writing - original draft prepa- ration.Aravind Ashok: Formal analysis, software, vali- dation.Zhonglai Wang: Supervision, funding acquisition, Z Wang, JHS Almeida Jr, A Ashok, Z Wang, SGP Castro: Preprint Page 17 of 20 Lightweight design of VAFW cylinders combining Kriging–based metamodels with PSO (a) 50 kN, NL=1 (b) 100 kN, NL=2 (c) 200 kN, NL=2 (d) 500 kN, NL=2 (e) 500 kN, NL=3 (f) 1000 kN, NL=2 (g) 1000 kN, NL=3 (h) 1000 kN, NL=4 Figure 18: Buckling modes for all cylinders optimized by the PSO–MM framework. Figure 19: Convergence of the critical buckling load CRfor Abaqus S4R and the single-curvature BFSC element. The proposed element converges with two order of magnitude less degrees-of-freedom than the S4R. resources. SaulloG.P.Castro: Conceptualization,supervi- sion, methodology, resources, software, investigation, writ- ing - original draft preparation.