A sequential finite element model updating routine to identify creep parameters for filament wound composite cylinders in aggressive environments José Humberto S. Almeida Jr.,Tales V. Lisbôa,Axel Spickenheuer,Luc St-Pierre BCs + model parameters INPUT OUTPUT OUTPUT INPUT FEM optimized parametersDirect problem Inverse problemFEMU: 1. GA: Partial optimum 2. GB: Global minimum Highlights A sequential finite element model updating routine to identify creep parameters for filament wound composite cylinders in aggressive environments José Humberto S. Almeida Jr.,Tales V. Lisbôa,Axel Spickenheuer,Luc St-Pierre •Identification of primary and secondary creep parameters •Heuristic algorithm to avoid local minima •Gradient–based algorithm to ensure finding the global minimum •Fine-tuning of elastic constants for cylinders in harsh environments A sequential finite element model updating routine to identify creep parameters for filament wound composite cylinders in aggressive environments José Humberto S. Almeida Jr.a,b,∗, Tales V. Lisbôac, Axel Spickenheuercand Luc St-Pierreb aAdvanced Composites Research Group, School of Mechanical and Aerospace Engineering, Queen’s University Belfast, Belfast, UK bDepartment of Mechanical Engineering, Aalto University, Espoo, Finland cLeibniz-Institut für Polymerforschung Dresden e.V., Hohe Str. 6, 01069 Dresden ARTICLE INFO Keywords : Creep modelling Finite element model updating Genetic algorithm Gradient–based algorithm Filament windingABSTRACT In this paper, a Finite Element Model Updating (FEMU) procedure is developed to find the best creep parameters for filament-wound cylinders under radial com- pression in harsh environmental conditions. Three winding angles are considered, each under three different hygrothermal conditions. The two–stage creep model captures i) primary creep through a time-hardening approach whilst ii) secondary creep is captured by Norton’s law. Given the high number of parameters in this two–stagecreepmodelandthecomplexityofdeterminingthemexperimentally,the FEMU routine utilises an optimisation scheme that sequentially couples a Genetic Algorithm (GA) with a gradient-based (GB) Levenberq-Marquardt Algorithm (LMA) to find all required creep input parameters to feed the model that best simulates experimental results. This framework finds the global optimum through aninitialscreeningoftheoptimumareawithinthedesignspacewithGA,clearing the path to allow the GB algorithm to find the global optimum, substantially reducing the chance or even avoiding falling in local minima. The global search is driven by experimental data of cylinders loaded in radial compression under aggressive environments. The numerical results show excellent agreement with experimental results with reasonably low computational efforts. 1. Introduction Filament winding (FW) is the most suited manufacturing process for axisymmetric structures, e.g., 1 tubes[1],cylinders[2],andpressurevessels[3].TheautomationofFWenablessimultaneousproduction 2 of identical parts, which also facilitates large volume productions, entailing a decrease in the overall 3 manufacturingcosts[4].TheFWprocessisalsosuitabletoproducelargecarbonfibre–reinforcedpolymer 4 (CFRP) composite components, for example, the Vega satellite launcher (European Space Agency) and 5 Boeing787fuselage.Anotherremarkablecharacteristicistheutilisationoflongcontinuousfibres,which 6 allows to design structures with a high stiffness and strength at the minimum possible weight [5]. 7 TheCFRPcomponentsusuallyoperateathightemperaturesandinhumidenvironments,whichaffect 8 their mechanical response over time [6, 7, 8]. The behaviour of the structure is then dependent on the 9 humberto.almeida@qub.ac.uk (J.H.S. Almeida Jr.); tales-lisboa@ipfdd.de (T.V. Lisbôa); spickenheuer@ipfdd.de (A. Spickenheuer); luc.st-pierre@aalto.fi (L. St-Pierre) ORCID(s):0000-0001-7511-2910 (J.H.S. Almeida Jr.); 0000-0003-1026-6243 (T.V. Lisbôa); 0000-0001-6485-824X (A. Spickenheuer); 0000-0003-3857-5398 (L. St-Pierre) Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 1 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders combinationoftemperatureandtime,whichaffecttheircreeprate[9].Thissetofparametershasamore 10 pronounced effect when the structure is laminated with off–axis layers, which can generate premature 11 micro-cracks, resin plasticization, and fibre-matrix interface weakening [10]. These micro-damages are 12 generated due to the viscoelastic nature of polymer matrices, which exhibit a time-dependent behaviour 13 and,hence,aresensitivetobothhygrothermalandtemperatureeffects[11,12,13].Therefore,thedesign 14 ofsuchCFRPstructuresmustconsidercreepandageingdegradationthatmightaffectthelifespanofthe 15 composite component [14, 15, 16]. 16 Experiments to characterise the long-term properties of composites are avoided or reduced, when 17 possible, given the high cost and time involved in such tests [16]. Therefore, computational models able 18 to determine such long-term behaviour of composites have been developed to fully or partially replace 19 experiments.Focussingspecificallyoncompositecylinders,theyusuallyoperateundertransverseloads, 20 eitherinundergroundorabove-groundapplications,causingchangesintheircross-sectionshape,where 21 perfect cylinders become elliptical ones [17]. This radial deformation increases with time due to the 22 viscoelasticnatureofthepolymer,whichmakespredictingsuchlong-termpropertiesachallenge[18,19]. 23 Several higher–order theories have been developed over the last decades to predict the short– and 24 long–termpropertiesofbeams,plates,andshells[20].Themostwelldisseminatedhigher–ordertheories 25 include Carrera’s (CUF) [21, 22, 23] and Generalised Unified Formulations (GUF) [24]. An accurate 26 calculation of internal stresses in laminated composites, in particular thick laminates, is essential to 27 properly design lightweight components. According to Carrera [25], there are three main approaches: 28 continuum base, asymptotic or axiomatic theories. Each of these theories has a different formulation to 29 predict displacement, strains and stress fields [26]. In this case, Unified Formulations is very attractive 30 sincetheunifiedsolutionmethodcanbechosenaccordingtotheproblembyeasilychangingtheexpansion 31 functions/terms [27], which might reflect in more precise failure predictions for composite structures. 32 Indeed, CFRP structures are designed to not fail by creep, however, their mechanical performances 33 canbehugelyaffectedbycreepdeformation.Hence,thedevelopmentofcomputationalmodelstopredict 34 creepstagesovertimeareofutmostimportance.Thecreepdeformationprocesscanbedividedintothree 35 stages: primary, secondary, and tertiary [28]. In the primary or transient stage, the strain rate is high 36 but decreases with time and strain since the structure has an increase in both creep strength and strain 37 hardening. Secondary creep is a steady–state stage, where the deformation rate decreases substantially 38 and strain increases slowly over time. In tertiary creep, a sudden increase in the creep rate takes place, 39 leading to the eventual failure of the system [29]. 40 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 2 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Computational models to predict creep deformation have been developed [9, 30, 31], however a 41 common issue in high–fidelity models is the large number of input parameters needed. Determining all 42 theseinputparametersexperimentallyisoftendeemedtooexpensiveand/ortime-consuming.Tothisend, 43 Lisbôaetal.[32]developedafiniteelementmodelupdating(FEMU)routinetoapproximatesimulations 44 toexperimentalobservations[33,34,35].FEMUisasensitivityapproachtoupdateinputpropertiesthat 45 best correlate to target results, for instance, experimental results [36]. The expected outcome of a well– 46 developed FEMU algorithm is a more efficient finite element (FE) framework, wherein the influence of 47 uncertainties is minimised through the iterative updating procedure of input values. 48 To the best of the our knowledge, there is no work in the literature reporting the optimal set of input 49 parameterstosimulatecreepincomposites.Previousworksoftenrelyonsimplecurve-fittingtechniques 50 usingempiricalpower–lawmodels,suchasBurger’s[37]andFindley’s[38]approaches.Inthesemodels, 51 thenumberofinputparametersisverylowanditisratherstraightforwardtofindthebestinputparameters 52 with simple algorithms. Both models are, nevertheless, empirically–based and, hence, not constitutive 53 laws and therefore, the evaluation and calculation of stresses over time are neither accurate nor feasible, 54 givingtheabsenceofconstitutiverelationstocapturesuchstresses.Inordertoovercomethislimitation,a 55 keymotivationofthecurrentworkreliesonthedevelopmentofaconsistent2–stagecreepmodelassisted 56 by a robust and computationally–efficient framework to find the best sets of input parameters that best 57 simulate experimental observations of composite cylinders in radial compression and under aggressive 58 environments. 59 In this paper, we develop an original multi–level FEMU framework to obtain the best set of input 60 parameters to predict primary and secondary creep stage for CFRP filament wound composite cylinders 61 under radial compression loading. Another key novelty is the determination of the best input values for 62 the cylinders under different hygrothermal conditions. The two creep stages are modelled here, where a 63 time–hardening approach is used to predict primary creep, whilst secondary creep is modelled through 64 a Norton–based model. Given that the cost function for both creep problems may have several local 65 minima,amulti–levelframeworkisdevelopedcombiningai)gradient–freeoptimisationroutine(genetic 66 algorithm – GA) to evaluate the objective function of the optimisation problem globally, that is, it is 67 usedtofindtheregioninwhichtheglobalminimumis;ii)followedbyagradient–based(GB)algorithm 68 using the optimum results from GA as its starting point to find the global minimum, strongly reducing 69 the chance of falling into a local minima. 70 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 3 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Figure 1: Details of the used equipment for creep tests (adapted from [38]). 2. Experimental details 71 A dedicated creep testing equipment was designed and fabricated to evaluate the behaviour of 72 compositecylindersunderradialcompressiveloadinginhightemperatureandhumidenvironments[38]. 73 The experimental setup is shown in Figure 1. 74 The cylinders were designed using the CadWind software [39] and manufactured by FW using a 75 KUKA 140 L100 robot with towpregs from Toray (T700-12K-50C fibre and UF3369 epoxy resin). The 76 cylindersproducedhadafibrevolumefraction( 𝑉𝑓)of 70%[12].Thecylindersweremanufacturedonto 77 a stainless-steel mandrel with a diameter of 50.8mm and a length of 1mm. Curing was performed in an 78 oven with air circulation for 5h at 120◦C. The cylinders were then cooled down and extracted from the 79 mandrel. Afterwards, the specimens were cut off at the desired length of 50mm. 80 The cylinders were manufactured with three different winding angles: [±60],[±75], and[±90] 81 (actually wound at [±89.6]). The first two configurations ( [±60]and[±75]) have a helical winding 82 trajectory while the last one ( [±90]) by hoop winding. The helical specimens were produced with a 83 winding pattern of 1:1 [2]. Moreover, three different environmental conditions were considered: room 84 temperature ( -), water at room temperature ( WRT) and hot water (40◦C) (HW). Room temperature is 85 considered 23◦C and relative humidity is 50%. 86 All cylinders were then subjected to long–term radial compression [40] at 25% of their ultimate 87 load (static tests were carried out in [38]). The specimens were kept under compression and in these 88 environmental conditions for 10 days (240h). Prior to the creep tests, both WRTandHWsamples 89 wereplacedinwatertodeterminetheirwateruptake(WU)untilequilibriumwasreached.Thecylinder’s 90 thickness, applied force and water uptake for all tests are summarised in Table 1. 91 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 4 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Table 1 Specimen stacking configuration, thickness and applied forces. Winding Angle Conditioning Thickness [mm] Force [N] WU ( 𝑀∞) [%] [±60]- 0.72±0.02 135.0- WRT 0.506±0.035 HW 1.784±0.016 [±75]- 0.77±0.05 181.5- WRT 0.399±0.085 HW 1.601±0.107 [±90]- 0.70±0.01 203.0- WRT 0.357±0.036 HW 1.716±0.107 Initial populationEvaluationFitness Calculation Selection Cross-over MutationNew populationBest members Final PopulationBest memberMaximum Generation Initial GuessPopulation ordered by FitnessJacobian MatrixStep Definition Evaluation of the new stepStop criterionJacobian Matrix Adjustment Damping decrease (GN)Update parametersOptimum solutionEvaluation Damping increase (SD) Valid StepStep 1: Genetic Algorithm Step 2: Levenberg-Marquardt Algorithm true false truefalsefalse true ReprodutionFEMU Framework Figure 2: Flowchart of the sequential optimisation procedure. 3. The FEMU framework 92 The FEMU is a procedure to fine-tune or identify unknown parameters by adjusting a parametrized 93 FE model and comparing the results against experimental data [32]. An optimisation procedure is used 94 to propose new sets of input parameters that will move the simulation closer to experiment. An usual 95 drawback of such procedure is a local minimum, which refers to a minimum within a neighbourhood 96 that may not be the global minimum. In order to decrease or even avoid the probability of falling into a 97 localminimum,asequentialprocedureusinggeneticalgorithm(GA)andLevenberg-MarquartAlgorithm 98 (LMA) is herein developed. A flowchart describing the sequential procedure is depicted in Figure 2. 99 Details of each step are provided next. 100 Figure 3 shows the FEMU procedure at a glance. In a directapproach, the geometry, boundary 101 conditions (BC), and material behaviour are known. An FE model is built to obtain the mechanical 102 response of the structure. With FEMU, the problem is solved in an indirectfashion. First, a numerical 103 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 5 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders boundary conditions + model parameters model optimized parametersFinite Element ModelSolution of the model Experimental evaluations Inverse (reversed) problemDirect (forward) problemINPUTOUTPUT INPUTOUTPUT Figure 3: Schematic illustration of a direct ×inverse problem highlighting the FEMU framework developed. modelisdeveloped,inwhichtheunknownsarebasedoninitialguesses,anditspredictionsarecompared 104 to experimental tests. Second, the input parameters are adjusted to minimise the difference between 105 simulationsandexperiments.Thisadjustmentisusuallydonewithanoptimisationschemecoupledwith 106 the parameterised FE model. Thus, FEMU solves the direct problem to compare the simulation output 107 with experimental data, whereas the inverse problem represents a conceptualisation of the framework 108 procedure. 109 3.1. The FE Model 110 TheFEmodelisdepictedinFigure4,whereinalldetailsoftheFEproblemareshown.Inordertoadd 111 realism in terms of modelling as–manufactured structures, the geometry is constructed considering the 112 filament winding pattern following the methodology proposed by Lisbôa et al. [32]. Since the problem 113 is symmetric, a symmetry plane is used to partition the cylinder in half. Initial simulations on full– 114 models have been carried out and the results were identical when compared to the half–size model. 115 therefore, in order to aid in computational efficiency, a half–size model is thereafter considered. A 116 vertical compressive force is applied to the top platen, and therefore the radial displacement is released, 117 whereas all other degrees-of-freedom are restricted. The bottom rigid platen is fully clamped. No BCs 118 areappliedtothecylinder,withtheexceptionofthecontactbetweenplatenandcylinder(Figure4).The 119 FE modelling is performed using Ansys APDL 2021 R1 FE platform, which is well suited to generate 120 parametricFEmodels.A4-nodeshellelementinitsdegeneratedtriangularshapewithequivalentsingle 121 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 6 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders contact regionrigid platecylinderzy xF clamped Figure 4: The developed FE model and its relevant details. layer formulation is utilised in the converged mesh (shown in Figure 4). The problem to be solved is 122 non-linear due to the cylinder–to–compressive plates contact and large displacements hypothesis, thus 123 geometric nonlinearity is considered in all simulations in which the Newton–Raphson method is used to 124 solve the FE problem. A surface–to–surface contact algorithm is utilised to model the contact between 125 the sample and the compressive plates, in which the contact formulation has an augmented Lagrangian 126 algorithm. 127 3.2. Norm definition 128 Since two different optimisation approaches are used in a sequential manner, the norm (distance 129 betweennumericalandexperimentalresponses)mustbedefinedinawaythatcanbeusedthroughoutthe 130 sequentialframework.Furthermore,theGAconsidersfitnessastheevaluationparameterwhereasLMA 131 uses the norm for the minimisation problem. The norm is then defined as: 132 1 𝑓𝑖𝑡=𝐿2=1 𝑛𝑝𝑡𝑠√ 𝒚T𝑾𝒚, (1) where𝑓𝑖𝑡isthefitness, 𝒚aretheresidualsand 𝑾istheweightingmatrix.ThefitnessisaGAparameter 133 and it defines how fit is the individual with respect to the environment and to other individuals while 134 the residuals is a vector that contains the differences between experimental and numerical data in each 135 evaluated point. The size of this vector is the number of points measuring the norm, 𝑛𝑝𝑡𝑠. The residuals 136 and the weighting matrix are written as: 137 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 7 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders 𝑦𝑖=𝑑E 𝑖−𝑑N 𝑖, (2) 𝑊𝑖𝑖=1 (𝑑E 𝑖)2, (3) where𝑑N 𝑖and𝑑E 𝑖arenumericalandexperimentalresponsesatthe 𝑖-thposition,respectively.Thismeans 138 that the outputs of the FE model are load and displacement, where FEMU evaluates every point of the 139 experimental curve pointwise. This is a key strategy to minimise the error and approximate simulations 140 to experiments. 141 According to Eqs.(1), (2) and (3), the norm is a percentile one. It is straightforward to see that the 142 closerthesimulationistotheexperiment,thesmallerthe 𝑦𝑖and𝐿2are.Likewise,acloserresponsewill 143 develop a larger fitness, which plays an important role in the GA as such individuals are to be preferred 144 in the selection for reproduction. 145 3.3. FEMU – step 1: Genetic algorithm 146 The global flowchart of the sequential FEMU procedure is depicted in Figure 2. A summary of the 147 GA is as follows: 148 •Initially, an initial population is considered, wherein each individual is defined by a set of real- 149 valueparameters( chromosomes ),𝝆,generatedbyarandomuniformdistributionoftheparameters 150 domain. 151 •Thispopulation is evaluated, i.e., the randomised parameters are inserted into the parametric FE 152 model for obtaining the mechanical response of the cylinder. 153 •Numerical and experimental are compared through the 𝐿2-norm (Eq.(1)), and the fitnessof each 154 individual is then calculated. 155 •Individuals (parents) are selected and, through crossovers , the new individuals (children) are 156 created. A mutation chance is applied to each new individual and a new population is generated. 157 The best members of the former population are also included in the new one. 158 •Then, the new members are evaluated and their fitnessis calculated. 159 •This procedure is repeated for a particular number of generations , at which the optimisation 160 procedure stops, and then a final population is reached. 161 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 8 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Afterallpopulationsareevaluated,theaveragefitness( ̄𝑓𝑖𝑡)ofeachmemberiscalculated,asfollows: 162 ̄𝑓(𝑖) 𝑖𝑡=𝑓(𝑖) 𝑖𝑡 𝑛𝑖𝑛𝑑∑ 𝑗𝑓(𝑗) 𝑖𝑡∴̄𝑓(𝑖) 𝑖𝑡∈(0,1), (4) With the average fitness of each member defined, the reproduction procedures (selection, crossover, 163 and mutation) are applied to the population. The parameters 𝑟𝑖∈ [0,1),𝑖∈ℕ, derived from a random 164 uniform distribution, are used in this procedure to add the required randomness to the optimisation 165 procedure. 166 Each reproduction selects individuals from the population. The crossover uses simulated binary 167 crossover (SBX) [41]. Essentially, the children are defined by their parents as: 168 𝜌(𝛼,𝑛+1) 𝑖=1 2[ 𝑟𝑣𝜌(𝛼,𝑛) 𝑖+(2−𝑟𝑣)𝜌(𝛽,𝑛) 𝑖] (5) 𝜌(𝛽,𝑛+1) 𝑖=1 2[(2−𝑟𝑣)𝜌(𝛼,𝑛) 𝑖+𝑟𝑣𝜌(𝛽,𝑛) 𝑖] , (6) where𝑖correspondstoaparticularchromosomein 𝝆,𝛼and𝛽identifytheparents, 𝑛definesthegeneration 169 and: 170 𝑟𝑣=⎧ ⎪ ⎪ ⎨ ⎪ ⎪⎩1+(2𝑟)1 𝜂+1,0≤𝑟<0.5 1+[ 1 2(1−𝑟)]1 𝜂+1,0.5≤𝑟<1,(7) inwhich𝜂∈[0,∞)isan“intensity”parameter.Forsmallvalues 𝜂<2,thechildren( 𝝆(𝛼,𝑛+1)and𝝆(𝛽,𝑛+1)) 171 are likely to be far away from the parents [41]. 172 After a new population is generated, the individuals may mutate. A non-uniform mutation operator 173 [42] is herein used. Initially, it is verified if the individual will mutate by generating an 𝑟number and 174 comparing it with some predefined value, normally very small, that represents the likelihood of an 175 individual being mutated. If the individual is mutated, then two other random numbers, 𝑟1and𝑟2, are 176 independently generated for each gene. With 𝑟1, one constructs 𝛿, a “jump” function [42], as follows: 177 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 9 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders 𝛿𝑖=⎧ ⎪ ⎨ ⎪⎩𝑈𝑖−𝜌𝑖,0≤𝑟1<0.5 𝐿𝑖−𝜌𝑖,0.5≤𝑟1<1,(8) where𝑈𝑖and𝐿𝑖definetheupperandbottomlimitsofthe 𝑖-thchromosomeof 𝝆and,with𝑟2,definesits 178 mutation as: 179 𝜌(𝑎,𝑛+1) 𝑖≐𝜌(𝑎,𝑛+1) 𝑖+𝛿[ 1−𝑟(1−𝑛∕𝑛𝑔𝑒𝑛)𝛾 2] , (9) in which𝑛𝑔𝑒𝑛corresponds to the maximum generation number and 𝛾represents the degree of non- 180 uniformity. Smaller values of 𝛾tend to push the value of the chromosome away from its original value 181 and the number of generations reduces the effect of the mutation in generations close to the final one. 182 After all these steps, a new population is constructed, considering the best members of the previous 183 one, and evaluated again until the maximum number of generations is achieved. Then, the final solution 184 is obtained. 185 3.4. FEMU – step 2: Levenberg-Marquardt Algorithm 186 TheLMAisagradient-based(GB)optimisationalgorithmthatisamixbetweenGauss-Newton(GN) 187 and the Steepest-Descent (SD) procedures [43]. Through a “damping factor”, 𝜇, the LMA moves from 188 GNtoSDmethodsandvice-versa.Whenalarge 𝜇isused,themethodswitchestoSD,whereasforsmall 189 𝜇, it remains in GN method [43]. 190 The main fundamental relation of LMA reads: 191 (𝐉T𝐉−𝜇𝐈)Δ𝝆=−𝐉T𝐲, (10) where 𝐉corresponds to the Jacobian matrix, Δ𝝆defines the step size, and 𝐈is the identity matrix. The 192 new step is defined as: 193 𝝆≔𝝆+Δ𝝆, (11) Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 10 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders The Jacobian matrix 𝐉is the variation of the residuals regarding the step and, numerically, are obtained 194 by forwarding finite differences, as follows: 195 𝐽𝑖𝑗=𝜕𝑦𝑖 𝜕𝜌𝑗≈𝑦𝑖(𝝆+𝛿𝑗Δ𝝆)−𝑦𝑖(𝝆) Δ𝜌𝑗, (12) noticing that 𝛿𝑗Δ𝝆refers to an evaluation only in the 𝑗−direction. 196 Revisiting Figure 2, one observes that in LMA: 197 •TheinitialguessisthebestmemberofthefinalpopulationobtainedthroughtheGAoptimisation. 198 •Throughanarbitraryincrementoftheinitialstep,theJacobianmatrixisconstructedbyevaluating 199 some auxiliary analyses. 200 •By solving Eq.(10), Δ𝝆is obtained, and a new step, 𝝆+Δ𝝆, is evaluated. 201 •Ifthenormdifferencebetweenthenewandoldanalysesisaboveatolerance, 𝜇decreases,andthe 202 newsetofdesignvariablesisupdated,andevaluated.otherwise,(inanon-validstep), 𝜇increases, 203 the Jacobian matrix is updated, and Eq.(10) is solved again for deriving a new step size. 204 •This cycle is completed when a tolerance in the step size is reached and the optimised solution is 205 obtained. 206 The validity of the step is determined by the following relation: 207 𝜆=𝐿2(𝝆)−𝐿2(𝝆+Δ𝝆) Δ𝝆T(𝜇Δ𝝆+𝐉T𝐲), (13) in which𝜆>𝜖and𝜆<𝜖correspond to a valid and invalid step, respectively, and 𝜖≤0. 208 The size of 𝜇is related to the problem in the evaluation and its value is modified during the 209 optimisation procedure in order to switch between the methods in an appropriate way. It is suggested 210 tovary𝜇dependingonthestep’svalidity.Thus, 𝜇ismodifiedthroughanauxiliarypositiveparameter 𝜏 211 as follows: 212 𝜇≔𝜏𝜇, (14) Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 11 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders where𝜏 >1when𝜆<𝜀and𝜏 <1when𝜆>𝜀. In other words, in a valid step, LMA “switches” to GN, 213 whereas for an invalid one, the method remains at SD. 214 For invalid steps, the Jacobian matrix is updated as: 215 𝐽≔𝐽+1 Δ𝝆TΔ𝝆(𝐲i+1−𝐲i−𝐉Δ𝝆)Δ𝝆T, (15) to improve the convergence where Δ𝝆is the step size that achieved the invalid step, as well as 𝐲i+1and 216 𝐲icorrespond to the previous and actual residuals, respectively. For each valid step, the Jacobian matrix 217 is recalculated through Eq.(12). 218 The𝐿∞-norm of the step size is considered as stopping criterion: the algorithm stops if this norm is 219 below some presettolerance. Thus, a mapping ofthe design variables is performeddue to their different 220 orders of magnitude. 221 For linear parameters, the transformation is a linear one, as follows: 222 𝑔=(𝐺−𝐿) 𝑈−𝐿∴ Δ𝑔=Δ𝐺 𝑈−𝐿, (16) while for logarithm ones, the mapping is defined as: 223 𝑔=ln𝐺−ln𝐿 ln𝑈−ln𝐿∴ Δ𝑔=ln(Δ𝐺+𝐺)−ln𝐺 ln𝑈−ln𝐿, (17) noticing that 𝑔and𝐺correspond to parameters in the mapped and non-mapped domains, respectively. 224 Essentially,thedomainsofnon-mappedandmappedvariablesare 𝐺∈[𝐿,𝑈]and𝑔∈[0,1],respectively. 225 This mapping is not required by the GA algorithm since the scale of the variables is insensitive to its 226 performance. 227 3.5. Elastic properties 228 Since the material properties change with temperature and water uptake, a simplified FEMU 229 procedureisappliedtofine-tunethefourelasticconstants– 𝐸11,𝐸22,𝐺12,and𝜈12oftheCFRPmaterial 230 systemhereinused.ThissimplificationmeansthatonlyGAwasusedasanoptimisationproceduresince 231 there was no need to use LMA to solve this problem. The material properties of the same towpreg were 232 tested[12]foranextendedtimeandforhighertemperatures.Theexperimentalresultswiththeirstandard 233 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 12 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Table 2 Non-aged (60 days at 23◦C, and relative humidity of 50%) and aged (60 days at 80◦C) elastic properties for the CFRP cylinders [12]. Properties Non-aged Aged 𝐸11[GPa] 129.3±3.6 119.7 ±4.1 𝐸22[GPa] 9.11±0.49 6.3 ±0.8 𝐺12[GPa] 5.44±0.023 3.89 ±0.19 𝜈120.322±0.023 0.331 ±0.011 deviations are presented in Table 2 and these values are considered as boundary limits for the elastic 234 properties to fine-tune. 235 Inordertoimprovethefine-tuning,thedefinitionsof 𝐺23and𝜈23followtheconsiderationsfromKuo 236 et al. [44]. They are defined by: 237 𝐺23=𝐸22 2(1+𝜈23), 𝜈23=−𝐸22[ 𝐸11( 1 2−𝜈12) +2𝐺12𝜈2 12] +Λ 2𝐸11𝐺12, (18) where 238 Λ=√ 𝐸2 22[ 𝐸11(1 2−𝜈12) +2𝐺12𝜈2 12]2 −4𝐸11𝐺12[ 𝐸11𝐸22(1 2−𝜈12) −𝐺12(𝐸11−2𝐸22𝜈2 12)] . (19) These two parameters ( 𝐺23and𝜈23) are derived solely by the four parameters listed above. 239 The experimental data available for the FEMU procedure is the instantaneous displacement. This 240 point was chosen as it retains the effect of temperature and water uptake and disregards any viscous 241 (creep) effects. The 𝐿2-norm considers only one point and is defined as: 242 𝐿2=1 2||||1−𝑑N 𝑑E||||, (20) Theupperlimitsoftheelasticmoduliareconsideredtosimplifytheproblemandreducethelikelihood 243 ofhavingmultiplesolutions.Lowerlimits(averageagedproperties),andthelimitsonthePoisson’sratio 244 (𝜈12∈[0.30,0.35]),arekeptfixedthroughouttheoptimisationprocess,whichisdoneinasequentialway: 245 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 13 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders firsttheWRTcase,andthentheHWone.ForcylindersunderWRTconditioning,non-agedpropertiesare 246 consideredasupperlimits.Then,theWRTconvergedpropertiesareusedasupperlimitsfortheparameter 247 identification of specimens under HW conditioning. This hypothesis is backed up by DMA experiments 248 from[12,38],whereonenoticesareductionofthetangentelasticmoduliforthesametowpreg,winding 249 angles,andunderthesameenvironmentalconditions.Moreover,theworkfromEggersetal.[38]reports 250 thatothermechanicalresponseswereevaluatedandthisreductionisindeedobserved.Hence,theelastic 251 moduli are defined as follows: 𝐸(−) 𝛼𝛽≥𝐸(WRT) 𝛼𝛽≥𝐸(HW) 𝛼𝛽;𝛼,𝛽={1,2}. For𝜈12, no hypothesis is made. 252 3.6. Creep constants identification 253 The experimental results for cases herein considered are presented in [38]. The creep behaviour for 254 HW conditioning in all configurations has both primary and secondary stages well defined. Thus, the 255 chosencreepmodelmustconsiderbothregions.Atime-hardeningcreepstrain-ratelawcanbedefinedas 256 a function of the stress, time, and temperature (Norton-Bailey law), as follows [45]: 257 ̇ 𝜀𝑐𝑟=𝑓(𝜎,𝑇,𝑡), (21) where𝜀𝑐𝑟istheequivalentcreepstrain, 𝜎thestress,𝑇and𝑡aretemperatureandtime,respectively,while 258 (⋅)denotestimederivative.Equation(21)canbefurthersimplifiedbyneglectinganycross-effectbetween 259 the parameters [45]. Thus: 260 ̇ 𝜀𝑐𝑟=𝑓1(𝜎)𝑓2(𝑡)𝑓3(𝑇). (22) Thecreepmodelhereinconsidereduses,essentially,amixtureoftime-hardeningandNortonmodels, 261 for describing both primary and secondary creep stages. The creep strain rate is then defined as: 262 ̇ 𝜀cr=̇ 𝜀N+̇ 𝜀TH, (23) where: 263 ̇ 𝜀TH=𝐶1𝜎𝐶2𝑡𝐶3exp(−𝐶4 𝑇) , (24) Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 14 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders ̇ 𝜀N=𝐶5𝜎𝐶6exp(−𝐶7 𝑇) , (25) in which𝐶𝑖,𝑖={1...7}are constants of the model and 𝑇is the temperature (Kelvin). It is important to 264 notice that𝐶1and𝐶5correspond to constants that depend on the material and the creep mechanism; 𝐶2265 and𝐶6definethecreepstressindices; 𝐶3determinesthetime-hardeningindex;and 𝐶4and𝐶7arerelated 266 to Arrhenius law approach (temperature–dependence) for each model [46]; 𝐶1and𝐶5may also vary for 267 materials that react with moisture [6]. 268 The creep strain is defined as: 269 𝜀𝑐𝑟=𝐶1𝜎𝐶2𝑡𝐶3+1 𝐶3+1exp(−𝐶4 𝑇) +𝐶5𝜎𝐶6𝑡exp(−𝐶7 𝑇) , (26) which was obtained simply by integrating over time Eqs.(24) and (25). 270 Adiscussiononthevaluesandrangesofvariablescanbefoundintheliterature[45].Thisisimportant 271 since this optimisation procedure requires closed domains. Thus: 272 •𝐶1and𝐶5: they must be positive to avoid violating strain energy principles. They also have order 273 of10−2to10−14(relative to creep strain rate) [45], with time 𝑡in seconds. 274 •𝐶2and𝐶6:𝜀𝑐𝑟must have1≤ (𝜎)≤5, which infers that the constants must not be necessarily 275 bounded by these limits. 276 •𝐶3: due to the tendency of the strain ×time curve,−1<𝐶3≤0. 277 •𝐶4and𝐶7:althoughnoinformationregardingthetheoreticallimitsoftheseconstantscanbefound 278 in the literature, these constants are related to the creep activation energy. Considering average 279 valuesfortheactivationenergy(65–85kJ/molofCFRP)[47],gasconstant,andtemperature,the 280 Arrhenius equation has an order of −11≤ (exp(−𝐶4∕𝑇),exp(−𝐶7∕𝑇))≤−7. 281 Giventhefactthat 𝐶1,𝐶4,𝐶5,and𝐶7areinfluencedbysimilareffects,themodelisfurthersimplified 282 to: 283 𝜀𝑐𝑟=̂𝐶1𝜎𝐶2𝑡𝐶3+1 𝐶3+1+̂𝐶5𝜎𝐶6𝑡, (27) hence,̂𝐶1and̂𝐶5vary with temperature, water uptake, and creep mechanism. 284 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 15 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders For the input parameters identification, these considerations are made: 285 •̂𝐶1and̂𝐶5varywithwateruptakeandtemperature.Foreachconfigurationandconditioning,these 286 constants must be evaluated. Moreover, their new domain is set as −25≤ (̂𝐶1,̂𝐶5)≤−9. 287 •𝐶2,𝐶3and𝐶6are solely dependent on the problem to be solved. Thus, they are independent of 288 temperature and water uptake. They are obtained once for HWas it is the most sensitive response 289 and then they are kept constant for other conditions. 290 Itisimportanttomentionthattheassessedpointsforthe 𝐿2-normevaluationinthecreepparameters 291 identification follow an exponential space, which is constructed as: 292 𝑡𝑖={ 1 𝑛𝑝𝑡𝑠−1[√𝑡𝑓𝑖𝑛(𝑖−1)+√ 𝑡𝑖𝑛𝑖(𝑛𝑝𝑡𝑠−𝑖)]}2 , (28) where𝑡𝑖,𝑡𝑓𝑖𝑛, and𝑡𝑖𝑛𝑖correspond to the time where the experimental and numerical data are assessed, 293 final and initial points of both datasets, respectively. Essentially, considering that 𝑡 >0, a linear space 294 is constructed over the square root of time and then passed onto the actual space. This results in 295 approximately half of the assessed points being condensed to the third part of the evaluation time. 296 4. Results and discussion 297 4.1. Fine-tuning of elastic properties 298 Thefirststeptoidentifythecreepconstantsistofine-tunetheelasticconstants.Thisstepisnecessary 299 since the well-known reduction of the properties given temperature and water uptake: the experimental 300 data shows that the instantaneous displacement – without any viscous effect – is different given the 301 conditioning. The optimisation procedure uses only the GA since the fitness values obtained were 302 excellent, not requiring further steps with LMA. The parameters related to the optimisation procedure 303 are listed in Table 3. 304 The results obtained are shown in Table 4, where 𝑢𝑠corresponds to the static displacement applied 305 in each case and 𝐺12=𝐺13and𝜈12=𝜈13for this geometry. The elastic moduli are slightly lower than 306 those of non-conditioned specimens (Table 1 – non-aged properties). Also, the post-processed material 307 parameters𝐺23and𝜈23are also presented for the optimised sets. The high fitness indicates an excellent 308 agreement between simulations and experiments. The convergence of the GA algorithm related to these 309 results is depicted in Figure 5, where the variable 𝑐𝑜𝑛𝑣is written as: 310 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 16 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Table 3 Parameters for GA procedures for the fine-tuning of the elastic properties. # Elastic variables 𝑛𝑖𝑛𝑑𝑛𝑔𝑒𝑛Best Members 𝜂Mutation Chance 𝛾 4 40 60 8 1 1% 1 Table 4 Elastic properties for all cylinders wound at different fibre angles and under several environmental conditions obtained by FEMU. Configuration 𝑢𝑠[mm]𝐸11[GPa]𝐸22[GPa]𝐺12[GPa]𝜈12𝐺23[GPa]𝜈23𝑓𝑖𝑡 [±60]WRT 4.10 127.87 9.09 5.43 0.348 2.63 0.728 1.06×106 HW 4.17 124.57 8.95 5.42 0.345 2.59 0.727 4.33×1013 [±75]WRT 2.75 126.96 7.69 4.72 0.329 2.25 0.708 >1030 HW 2.86 126.51 7.28 4.39 0.334 2.13 0.712 >1030 [±90]WRT 2.30 122.40 7.42 4.89 0.311 2.18 0.701 >1030 HW 2.32 121.76 7.03 3.99 0.303 2.14 0.642 >1030 𝑐(𝑖) 𝑜𝑛𝑣=log𝑥𝑖 log𝑥𝑐𝑜𝑛𝑣, (29) where𝑖isthegeneration, 𝑥𝑖istheaveragefitnessofthe 𝑖-thgeneration,and 𝑥𝑐𝑜𝑛𝑣isthebestmemberatthe 311 last generation. Excluding the ±60◦, all other analyses reached the maximum possible fitness, perfectly 312 matching the experimental data(𝑢𝑠). And even for the worst-case, ±60◦WRT, a value of 𝑐𝑜𝑛𝑣= 0.5 313 meansthatthenumericalvalueis 99.9%oftheexperimentalone.Now,withelasticpropertiesfine-tuned 314 foreachconditionandwindingangle,theelasticanalysisisfullydefined,andcreepparameterscanthen 315 be identified. 316 Itisnaturaltoassumethattheobjectivefunctionmighthavesomeglobalminima( 𝐿2=0).Toavoid 317 those minima, a qualitative approach, relying on dynamic mechanical analyses (DMA) from Almeida Jr 318 etal.[12],wascarriedouttounderstandtheglobalbehaviourofagedlaminates.Ingeneral,itisassumed 319 thatagedpropertiesarelowerthannon–agedones[6,7,8].Foreachcondition,GAanalyseswerecarried 320 outandthesetsofpropertiesshowninTable4representthebestsetofresultsconsideringi)thenorm,ii) 321 correlation with DMA data and iii) connection with aged properties from [12]. This series of qualitative 322 analysesarecrucialtoincreasethereliabilityofthedataobtainedandpresentedinTable4,andthishelps 323 with the prediction of creep constants, which is presented next. 324 4.2. Creep parameters identification 325 Forthesakeofrepeatability,allparametersusedintofindtheoptimumcreeppropertiesarelistedin 326 Table5,whereastheupperandlowerlimitsarepresentedinTable6.IntheGAstep,thedesignvariables 327 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 17 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Figure 5: Convergence of the GA for fine-tuning elastic properties. The best and average fitness are shown for (a) HW and (b) WRT conditions. Table 5 Parameters for GA and LMA optimisation procedures for creep parameters identification. General GA LMA 𝑛𝑝𝑡𝑜𝑠𝑡𝑖𝑛𝑖[s]𝑡𝑓𝑖𝑛[s] #𝑛𝑣𝑎𝑟𝑛𝑖𝑛𝑑𝑛𝑔𝑒𝑛Best Members 𝜂Mutation Chance 𝛾Δ𝜌𝑖𝑛𝑖𝜏𝑖𝑛𝑖 80 0 864000 5 50 35 10 1 1% 1 0.01 0.025 Table 6 Upper and lower limits for GA, which are further utilised for mapping the LMA step. Limits# Design variables ̂𝐶1𝐶2𝐶3̂𝐶4𝐶5 Upper10−90.01 −0.99 10−90.01 Lower10−255 0 10−255 areconstrainedtothesevalues,whereasintheLMAstep,theyareusedformappingtheparameters(Eqs. 328 (16) and (17)). 329 The optimum creep constants obtained by the sequential FEMU procedure are presented in Table 7. 330 The hypothesis of fixing 𝐶2,𝐶3and𝐶6for all conditions explains why these constants vary with the 331 winding angle only. This hypothesis seems adequate as the fitness of WRTand-cases is higher than 332 HW, with an exception for [±90]-WRT. 333 The simulated creep curves, using the constants in Table 7, are compared to experimental data in 334 Figure6.Thecorrelationbetweennumericalandexperimentalresponsesisexcellentforallcases,which 335 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 18 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Table 7 Creep constants obtained by FEMU. Stacking Condition ̂𝐶1𝐶2𝐶3̂𝐶5𝐶6𝑓𝑖𝑡 [±60]HW 5.030×10−15 3.802 -0.54891.884×10−12 1.38174.360 WRT 1.727×10−157.101×10−18135.047 -9.286×10−161.000×10−24201.803 [±75]HW 3.050×10−16 3.836 -0.3751.000×10−24 1.53693.877 WRT 7.389×10−202.223×10−14135.132 -2.557×10−173.635×10−14198.374 [±90]HW 3.401×10−18 4.095 -0.1187.562×10−24 1.944145.139 WRT 4.947×10−191.000×10−2459.422 -3.774×10−193.444×10−20101.629 was impossible to achieve by running FEMU with a single algorithm (either GA or LMA). In all cases, 336 bothprimaryandsecondarystagesareproperlyidentified,showingtheefficiencyoftheapproachesused 337 in the FE model. 338 The creep curves shown in Figure 6 are given for three winding angles: (a) ±60◦, (b)±75◦and (c) 339 ±90◦. In general, the higher the fibre angle, the lower the displacement is over time. This is expected 340 since off–axis layers (e.g., ±60◦and±75◦) are more dependent on the matrix, which has a viscoelastic 341 response. Likewise, when fibres are oriented along the loading direction (i.e. radial compression), the 342 deflection over time is lower since the response is controlled primarily by the fibres. 343 The environmental conditions affect all cylinders in a similar way. Regardless of the winding angle, 344 specimens under ML have the lowest displacement. Adding water at room temperature (WRT) slightly 345 increasesthedisplacementovertime.Incontrast,hotwater(HW)drasticallyincreasesthecreeprateand 346 hencehasthehighestdeflectionlevels.Thisisexpectedsincecreepistime–andtemperature–dependent, 347 and, since the time is constant for all cylinders, the temperature has a strong influence on the creep 348 responseofthestructures.EventhoughHWcylindershaveahighcreepratethroughouttheirsecondary 349 creep stage, they did not enter the tertiary creep stage. 350 The convergence characteristics of the GA–LMA for all cases are shown in Figure 7. A vertical line 351 separatesthecontributionofeachstepofFEMU.TheresultsshowthatwhileGAgoesthrough"jumps", 352 LMA has a more "steady" behaviour while navigating around the objective function. This behaviour is 353 noticed in the cases of [±60]-HW,[±60]-WRT, and[±90]-HW, with a focus on the latter. 354 Furthermore,itisimportanttoexplainwhyLMAhaslittleinfluenceonthe WRTandMLconditions. 355 In these two cases, only two design variables were optimised whereas the HWcylinders had five 356 parameters. The reduced complexity meant that the GA algorithm could achieve, with the same number 357 of generations, much better results and the global minimum is close to the LMA initial step. 358 ADD HERE DISCUSSION REV #2 - first question (MAYBE) - interpolation and extrapolation 359 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 19 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Figure 6: Computational results from FEMU compared to experimental data for (a) [ ±60], (b) [±75], and (c) [±90] cylinders. 4.3. Stress analysis 360 Although the main objective here is to introduce the sequential optimisation to identify creep 361 parameters,thesimulationswerealsousedtocomparethestressesandequivalentstrainrateforthethree 362 windinganglesandtheirdifferentenvironmentalconditions.Firstly,fourpoints,showninFigure8,were 363 selectedtopresenttheseresults.Moreover,sinceshellwereusedintheFEmodel,theresultsarereported 364 for both top and bottom surfaces. 365 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 20 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Figure 7: Convergence of the GA and LMA optimisation procedures for (a) [±60], (b)[±75], and (c)[±90] cylinders. The number of generations ( 𝑛𝑔𝑒𝑛) is also considered for the LMA steps. The stresses at point 1 are given in Table 8, and those for point 2 are listed in Table 9 (see Figure 8 366 forthelocationofbothpoints).Ineachtable,valuesaregivenfor 𝜎11,𝜎22,and𝜎12representingthestress 367 along the fibre, matrix and in–plane shear directions, respectively. Each component of stress is given at 368 four different times to capture the creep evolution. Important to mention that the shell elements used in 369 thisfinitestrainanalysisallowforchangesinthickness( 𝜀33≠0),andthisvariationistakenintoaccount 370 in the stress calculations. 371 Interestingly, without any viscous effects (0 s), the stress in the fibre direction is similar ( ≈ 350 372 MPa) for all winding angles. In contrast, 𝜎22and𝜎12vary. This behaviour is expected since creep is a 373 phenomenon that is more pronounced on the polymer instead of the fibres, that is, transverse stresses 374 aremoreaffectedbycreepthanlongitudinalones.Furthermore,atthisstage,thestressfieldseemstobe 375 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 21 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders zxy1 23winding angle Figure 8: The four points used to report stresses and strains. Points 1 and 3 are on the contact region, whereas points 2 and 4 are on the mid-section of the cylinder. The local coordinate system is aligned with the fibre direction. Table 8 Stress components at point 1 for both top and bottom surfaces of the cylinder. Values are in MPa and given at different times. TOP BOTTOM Cond.𝝈0s 172s 432s 864s 0s 172s 432s 864s 60◦HW𝜎11-326.71 -212.43 -201.43 -194.06 350.73 185.02 168.06 155.59 𝜎22-56.798 -76.719 -85.367 -93.820 8.5201 22.932 29.193 36.360 𝜎1230.368 40.264 42.675 44.213 55.543 56.866 58.797 60.017 WRT𝜎11-328.31 -248.17 -232.83 -221.94 352.37 231.57 212.46 199.21 𝜎22-56.744 -64.992 -67.488 -69.649 8.560 14.899 16.848 18.587 𝜎1229.732 35.151 36.386 37.336 54.869 54.064 54.127 54.243 –𝜎11-329.25 -269.86 -254.75 -243.19 353.44 259.31 239.52 224.97 𝜎22-55.940 -61.564 -63.484 -65.208 7.8301 12.224 13.693 15.022 𝜎1229.622 33.550 34.643 35.545 54.762 54.064 54.025 54.056 75◦HW𝜎11-355.31 -256.57 -239.13 -228.18 359.29 232.56 213.56 201.98 𝜎22-39.435 -46.167 -49.054 -51.756 -13.809 -6.686 -3.698 -0.7574 𝜎1223.194 25.494 26.194 26.689 33.658 33.538 33.717 33.874 WRT𝜎11-351.68 -350.47 -348.86 -346.28 354.08 351.62 349.35 345.71 𝜎22-39.494 -39.553 -39.768 -40.123 -13.467 -13.289 -13.144 -12.905 𝜎1223.922 23.838 23.932 24.085 34.351 34.460 34.518 34.613 –𝜎11-345.75 -322.75 -309.17 -295.81 342.80 310.45 292.94 276.42 𝜎22-41.545 -42.664 -43.632 -44.835 -12.867 -11.411 -10.395 -9.1887 𝜎1224.744 25.021 25.316 25.664 35.401 35.238 35.180 35.184 90◦HW𝜎11-334.87 -266.55 -250.42 -244.02 334.78 266.65 250.49 244.02 𝜎22-5.3425 -10.027 -13.979 -18.962 5.7682 10.468 14.441 19.451 𝜎12-0.3564 -0.3064 -0.2799 -0.2599 0.3562 0.3071 0.2813 0.2619 WRT𝜎11-335.01 -310.64 -292.50 -276.53 334.92 310.64 292.56 276.62 𝜎22-5.9802 -7.0207 -8.0585 -9.3667 6.4314 7.4714 8.5098 9.8211 𝜎12-0.4089 -0.3962 -0.3840 -0.3705 0.4088 0.3965 0.3845 0.3712 –𝜎11-334.55 -314.21 -297.78 -282.20 334.46 314.21 297.83 282.28 𝜎22-7.2396 -8.2008 -9.1965 -10.462 7.7487 8.7073 9.7013 10.966 𝜎12-0.4105 -0.4000 -0.3898 -0.3776 0.4103 0.4002 0.3902 0.3783 insensitivetotheenvironmentalconditions.Inallcases,thestressinthefibredirection 𝜎11decreaseswith 376 time.Thisreductionisparticularlyimportantforthecylinderwith 60◦andHW,anditismorepronounced 377 at point 1 rather than at point 2. In contrast, both 𝜎22and𝜎12increase sharply with time. Again, the HW 378 cylinder wound at 60◦shows the largest variations in stresses. 379 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 22 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Table 9 Stress components at point 2 for both top and bottom surfaces of the cylinder. Values are in MPa and given at different times. TOP BOTTOM Cond.𝝈0s 172s 432s 864s 0s 172s 432s 864s 60◦HW𝜎11179.69 143.46 131.68 122.59 -179.44 -144.69 -133.52 -125.07 𝜎2218.829 24.721 28.944 34.165 -21.153 -26.936 -31.129 -36.278 𝜎12-21.366 -24.807 -27.297 -30.065 -21.121 -24.956 -27.605 -30.508 WRT𝜎11180.66 165.75 160.10 155.09 -180.35 -166.12 -160.64 -155.75 𝜎2218.836 20.358 20.980 21.586 -21.172 -22.628 -23.225 -23.809 𝜎12-21.038 -21.842 -22.162 -22.477 -20.781 -21.728 -22.106 -22.476 –𝜎11181.17 172.19 168.26 164.50 -180.83 -172.33 -168.55 -164.91 𝜎2218.427 19.357 19.758 20.159 -20.770 -21.659 -22.044 -22.429 𝜎12-20.974 -21.468 -21.672 -21.876 -20.708 -21.287 -21.528 -21.769 75◦HW𝜎11186.13 169.29 161.34 154.92 -188.30 -171.24 -163.22 -156.75 𝜎227.4495 8.8926 9.7640 10.704 -10.760 -12.074 -12.859 -13.708 𝜎12-15.125 -16.010 -16.494 -17.001 -15.215 -16.157 -16.676 -17.218 WRT𝜎11183.45 183.28 182.79 181.99 -185.76 -185.62 -185.14 -184.37 𝜎227.5491 7.6006 7.6767 7.8027 -10.836 -10.895 -10.969 -11.092 𝜎12-15.495 -15.547 -15.598 -15.681 -15.595 -15.646 -15.699 -15.785 –𝜎11178.56 176.69 174.87 172.48 -181.15 -179.24 -177.40 -174.80 𝜎228.2541 8.5001 8.7469 9.0996 -11.607 -11.842 -12.072 -12.402 𝜎12-16.049 -16.205 -16.342 -16.536 -16.157 -16.316 -16.460 16.661 90◦HW𝜎11197.86 190.14 185.04 182.62 -203.51 -193.13 -186.95 -184.08 𝜎223.5615 4.6671 5.8848 7.6164 -3.5272 -4.6418 -5.8566 -7.5751 𝜎12-0.0043 -0.0046 -0.0051 -0.0059 -0.0086 -0.0086 -0.0087 -0.0092 WRT𝜎11197.88 196.58 194.99 192.77 -203.46 -201.56 -199.36 -196.43 𝜎223.8170 3.9995 4.2171 4.5344 -3.7866 -3.9709 -4.1914 -4.5117 𝜎12-0.0049 -0.0047 -0.0049 -0.0050 -0.0090 -0.0090 -0.0089 -0.0089 –𝜎11197.29 196.26 194.98 193.15 -202.85 -201.34 -199.56 -197.10 𝜎224.5638 4.7296 4.9283 5.2195 -4.5627 -4.7307 -4.9329 -5.2285 𝜎12-0.0052 -0.0051 -0.0051 -0.0053 -0.0094 -0.0093 0.0093 -0.0092 The large stress variations due to creep can have significant effects on the safety coefficients and 380 in–service margins employed in design. As mentioned, the variations are detrimental since there is a 381 stressreductioninthefibredirection( 𝜎11)butanincreaseinthetransversedirections( 𝜎22and𝜎12)[12]. 382 Consider, for example, a HW cylinder wound at 60◦. With the stress field at 𝑃2-bottom (and using non- 383 aged properties from [12]), the instantaneous (0 s) and 10–day (864s) failure index (Tsai-Wu failure 384 criterion [48]) is 0.387 and 0.892, respectively. This represents a considerable increase of ~230%. 385 These stress variations are attributed to three factors: i) creep deformation is changing the initial 386 cross-section from a circle to an ellipse [49], reducing the effective stiffness of the cylinder, ii) the load 387 redistribution due to the contact with the compressive platens, and iii) non-linear effects. As observed 388 by Lisbôa et al. [2, 32], during the initial loading stage, the force ×displacement curve is mostly 389 linear–elastic. At an intermediate stage, this behaviour changes slightly, and the cylinder becomes more 390 compliant. At the final stage, the specimen becomes slightly stiffer again, nonetheless, this effect is to 391 some extent shadowed by major cracking on the cylinder (i.e., damage initiation). In creep experiments, 392 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 23 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders 0 0.85 1.70 2.55 3.40 4.25 5.10 5.95 6.80 7.6543.3 mm56.7 mm 44.1 mm57.5 mm (a) (b) Figure 9: Comparison between (a) numerical and (b) experimental deformed shapes for the cylinder 60◦ (HW). The colour-bar represents the predicted vertical displacement. The photograph of the test was taken after unloading. the deflection increases up to the secondary creep stage, in which some softening is observable, mostly 393 for the HW condition. Furthermore, as the cross-section becomes elliptical (Figure 9), the contact area 394 between the cylinder and platen increases and, as the force is a fixed parameter, the contact stresses 395 are reduced. This can explain the larger stress variations at point 1 (contact area) compared to point 2. 396 Moreover,itisalsoworthmentioningtheexcellentpredictioncapabilityofthemodelinwhichthecross– 397 section measurements are extremely close. The difference in the major and minor axes measurements is 398 because the experimental values were taken after the sample was unloaded, in which stress relaxation is 399 unavoidable. 400 A comparison between the numerical and experimental deformation modes is shown in Figure 9 for 401 the 60◦(HW) cylinder. The deformed shape and eccentricity in the simulations are in good agreement 402 with experiments. The major and minor axes are 1-2% larger in the experiment because the photograph 403 was taken after unloading (whereas there is no unloading in the simulation). 404 The equivalent strain rate is plotted in Figure 10 as a function of time. The strain rate is given for 405 the four points identified earlier in Figure 8. All cases exhibit a similar decreasing response except for 406 Figure 10(e), where the strain rate is almost constant. This means that the time-hardening effects are 407 almostabsent;indeed,themeasureddisplacementfor [±75]-(WRT)ispracticallyconstantovertime,see 408 Figure 6(b). Consequently, ̂𝐶1is the smallest in all cases, see Table 7. For all other cases, the strain rate 409 starts at a high value and decreases with time. This decrease tends to slow down after the 3𝑟𝑑/4𝑡ℎday, 410 which corroborates with the displacement versus time curves (see Figure 6) and indicates the beginning 411 of the second stage of creep. The results also indicate that the contact areas (points 1 and 3) have higher 412 strainratesthanthemidsection(points2and4).Thisisattributedtotheelevatedstressintheseareas(see 413 Tables 8 and 9). 414 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 24 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders Figure 10: Log-log plots of the equivalent strain rate as a function of time. Results are given at the four points shown in Figure 8, and for (a) 60◦(HW), (b) 75◦(HW), (c) 90◦(HW), (d) 60◦(WRT), (e) 75◦(WRT), (f) 90◦(WRT), (g) 60◦(ML), (h) 75◦(ML), and (i) 90◦(ML). As can be seen, the optimisation problem here is well solved in a computationally–efficient way 415 and with high degree of accuracy, which shows the great capability of the proposal multi–step FEMU 416 framework. Other approaches to solve the same problem might also be suitable, such as the use of deep 417 neural networks, a class of neural networks, to approximate the functions herein utilised for training the 418 FE models [50]. 419 Almeida Jr. JHS, Lisbôa TV, Spickenheuer A, St-Pierre L: Preprint submitted to Elsevier Page 25 of 29 A sequential finite element model updating routine to identify creep parameters for composite cylinders 5. Conclusions 420 An original finite element model updating framework has been here developed as a new method to 421 identifythecreepparametersoffilamentwoundcylindersunderdifferentharshenvironmentalconditions. 422 An original FE model is developed, in which the winding pattern is taken into consideration. All 423 cylinders were subjected to radial compressive creep loading for 10 days. The framework is based on 424 a sequential optimisation that starts with a heuristic genetic algorithm and finishes with the gradient– 425 based Levenberg-Marquardt algorithm, aiming at reducing and/or avoiding the likelihood of finding a 426 local minima. The results were in excellent agreement with experimental observations, which indicates 427 that the computational model and its hypotheses were well formulated. 428 After identifying the parameters of the creep model, the stresses and equivalent strain rate were 429 exportedfromtheFEmodelandinterestingcharacteristicswereobserved:i)thesofteningofthecylinders 430 due to the circular-to-elliptic transformation of the cross-section, ii) more pronounced reduction in the 431 transverse direction than in the longitudinal one (typical of creep consequences for continuous fibre– 432 reinforced laminates with off–axis plies), iii) the natural reduction of failure index as a consequence of 433 point ii), even with the same applied load. These results strengthen the importance of creep behaviour 434 analysis in composite structures exposed to harsh environments. This framework can augment the 435 experimental data, reduce experimental costs (since creep tests are expensive and time–consuming), 436 and elucidate creep effects on FW composite cylinders exposed to harsh environments with an efficient 437 computational tool. 438 Acknowledgements 439 HAissupportedbytheRoyalAcademyofEngineeringundertheResearchFellowshipscheme[Grant 440 No. RF/201920/19/150]. 441