FEM updating for damage modeling of composite cylinders under radial compression considering the winding pattern

.


Introduction
Fiber-reinforced composite structures are being increasingly used in structural applications, such as in space launch vehicles, fuselages, and energy storage systems (e.g.pressure vessels) due to their outstanding high strength-and stiffness-to-weight ratios, and high-corrosion and fatigue resistances [1][2][3][4].Another key aspect behind the systematic replacement of metallic components by fiber-reinforced composites is their high design freedom, in which placing fibers in suitable directions substantially enhances the performance of the component [5][6][7].This design freedom is, nonetheless, dependent not only on the orthotropy of the system but also on the manufacturing process [8,9].Focusing on closed shells, filament winding (FW) is widely recognized as the most adequate production technique to manufacture cylinders and pressure vessels [10,11], for instance.However, when a filament-wound structure is designed, an aspect from the manufacturing process is strength.Morozov [18] developed a finite element (FE) model to investigate the stress and strain distribution of FW cylinders tubes under internal pressure.He concluded that the WP generates a non-uniform stress field throughout the shell and the stress levels are higher for the models with WP modeled than conventional models in which the shell has nominal ± angle-ply angles.Hernández-Moreno et al. [19] analyzed composite tubes under external pressure and concluded that the WP has no influence on the buckling and collapse strength.Uddin et al. [20] developed a model similar to Ref. [18] but focusing on predicting the mechanical response of a flywheel disk under rotational loading.The authors concluded that the stress distribution around the disk is sensitive to the WP.Zakrzhevskii and Khitrov [21] studied the effect of the interweaving on torsion of thin-walled filament-wound structures and found that the load capacity can be incremented by just increasing the interweaving of the fibers, which means increasing the pattern number.Azevedo et al. [22] experimentally investigated the mechanical response of composite cylinders with WP of 1/1, 3/1, and 5/1 under axial compression.They concluded that the higher the WP, the higher the strength; also, the highest stiffness was found for cylinders with a WP of 3/1.Lisbôa et al. [14] evaluated experimentally the WP influence of composite cylinders under radial compression.The stiffness does not vary statistically among the patterns, whereas the shape of load × displacement curves, the maximum supported load, and energy at break were strongly dependent on the pattern.
Considering that the study carried out by Lisbôa et al. [14] provides comprehensive experimental results, developing a computational model able to take into account the manufacturing characteristics of the winding process becomes necessary.Moreover, the development of a computational methodology able to reproduce the winding pattern is essential to capture more accurately stress fields throughout the shell.Given the fact that simulating experimental results of large fiber-reinforced FW shells under complex loading cases (i.e.radial compression) is extremely challenging, optimization algorithms can further assist stand-alone FE models.An efficient way to approximate simulations to experiments is by ''updating'' parameters of the virtual model.This procedure is called finite element model updating (FEMU) [23] and is based on a sensitivity analysis that iteratively updates coefficients/properties (such as elastic properties, strength allowables, damage coefficients) with the aim to better correlate simulations to target values, i.e. experimental response.For a well-structured FEMU procedure, a usual outcome is a more reliable FE framework, in which the effect of uncertainties is minimized by reducing errors associated with a stand-alone FE model.
Markiewicz et al. [24] presented a comprehensive review on inverse problems for parameter identification, characterizing some methodologies, which includes FEMU.The FEMU method is a very multidisciplinary tool and it is usually applied to solve many engineering problems, such as the parameter identification for damage detection of structures [25] and plasticity [26].The procedure is also used in the dynamics of structures.Zhao et al. [27] updated both mass and stiffness distributions of a composite flying-wing aircraft by using the method and the authors considered two strategies for evaluating the fitness function: the root mean square error (RMSE) and the mean absolute error (MAE).For the latter, a particular optimization procedure was developed using sequential linear programming and it was found that the updated FEM could achieve exactly the test data for MAE while RMSE distributed the error in all responses.Kumar Bagha et al. [28] enhanced numerical predictions of a structure from its experimental modal analyses applying FEMU.In composite materials, Rahmani et al. [29] developed a procedure for accurate identification of mechanical properties of composite components using micromechanical models, where FEMU was very effective in finding the best input properties, while Lim et al. [30] developed a multiscale damage plasticity model and inverse characterization for particulate composites, where FEMU was used to find the best material properties.
As can be seen, FEMU is a very efficient approach to approximate FE models to match experimental results over a wide range of materials, loading cases, engineering scenarios, etc.Given the fact that composite cylinders present a highly non-linear response when loaded in radial compression [14], a model able to predict these non-linearities and identify damage initiation and propagation is essential.One of the main difficulties in a damage model is to properly identify input parameters.Although it is possible to do so experimentally, this is a very costly and time-consuming procedure, which is not attractive for industrial purposes.Therefore, a framework able to update input parameters of a damage model can be an interesting alternative to overcome that.
In this context, this work proposes a methodology to computationally construct the mosaic pattern formed during the filament winding process for composite cylinders.The following patterns are considered: 1/1, 2/1, and 3/1.Furthermore, the cylinders are subjected to radial compressive loading.An efficient damage model is proposed to predict the progressive failure of the cylinders.Considering that commercial FE software platforms are developed to cover a range of materials, microstructures, and geometries, this work introduces a complementary finite element model updating (FEMU) routine to calibrate simulations with respect to experiments considering the winding pattern generated by the filament winding process, aiming at evaluating the influence of the winding pattern on damage propagation.
The manuscript is further organized as follows: Section 2 describes the FE model along with the pattern generation, progressive damage model, and FEMU approaches; Section 3 briefly presents details on manufacturing and testing; Section 4 presents results and discussions; and Conclusions are drawn in Section 5. Further results are presented as the supplementary material.

The pattern generation
The WP is an intrinsic feature of the FW process for either helical or polar winding [13].Fig. 1 depicts a scheme to generate two different patterns with the same winding angle.In Fig. 1(a), , , and  correspond to the winding angle, pitch, and tow width, respectively [13,31,32].In Fig. 1(b), the second stroke is presented and the first intercrossing regions are shown.In Fig. 1(c), the third and fourth strokes are depicted and then the pattern is defined: figures on the top and bottom correspond to the WP = 1 and WP = 2, respectively.In WP = 1, the third stroke is adjacent to the first one and, the fourth stroke is adjacent to the second one.The process continues until full coverage is reached, as shown in Fig. 1(d) on top, which results in one diamondshaped form in one pitch, .In a different manner, the third stroke of WP = 2 is deposited on a diametrically opposite position of the first one (180 • in the circumferential distance), developing then the aforementioned pattern.It is worth mentioning that the trajectories in Fig. 1 are geodesic (straight lines in a cylindrical developed surface).The most noticeable difference between the patterns is the size of the partitions (divisions created by the diamond-shaped forms in the circumferential direction).The pattern number, an integer, is defined as the number of these divisions/partitions.In greater detail, the periodic diamond-shaped form is presented in Fig. 1(e) and shows each side of the diamond with a different laminate configuration ( ).Three regions are highlighted: regular laminate, circular and helical cross-over.

The FE model
The pattern is generated by dividing the geometry into diamondshaped forms, following manufacturing observations.Then, the laminate properties (i.e.stacking) are properly considered in each region.Fig. 2(a)-(i) presents the graphical output from FW-CAM software Cadwind, the manufactured cylinders, and the modeling approach for identical cylinders except for the pattern: WP = 1 (Figs.(f), and (i)).The colors in Fig. 2(g)-(i) highlight different stacking configurations.This approach makes the FE mesh independent of the WP.
One can observe small differences among the patterns in Fig. 2, particularly to lines that construct the triangles.WP = 1 and WP = 2 have several lines inside the diamonds whereas WP = 3 has none.This refers to the algorithm for constructing such geometry and the possibility of multi-valued functions for the two former patterns.Thus, since WP = 4 encompasses both patterns [13], it is used as a base and the triangles are redefined to resemble a particular pattern (WP = 1 or WP = 2).
The algorithm to create the pattern in a cylindrical geometry is relatively simple and straightforward.Geometry generation is developed in two steps: the creation of the lines and areas and the assignment of areas to a respective stacking configuration.One can also further divide the first step into 2 parts: definition of the geodesic and circular lines and formation of the areas.The geodesic lines can be considered as straight lines in a cylindrical coordinate system.Thus, they can be directly built by just knowing an initial point and the winding angle (and, of course, the length).The number of lines required to develop the helical cross-overs of a particular pattern is twice the number of the pattern, due to the clockwise and anticlockwise lines.
Here lie the aforementioned multi-valued solutions: by choosing the initial and final points, two different trajectories are possible (clockwise and anticlockwise).That is why, for WP = 1 and WP = 2, a larger pattern is chosen to define the lines and then the algorithm that assigns a particular area is modified.The circular lines are also directly related to the WP and winding angle and they simulate the circular cross-over.For each pitch,  ℎ , which is a full turn of a tow onto the cylinder, there are 2WP + 1 circular lines, which means that the distance between two circular lines (or two circular cross-overs) is  ℎ ∕2WP.By defining both sets of lines, the generation of the areas is straightforward (regarding the boundary lines).The assignment of stacking sequence configuration to the areas considers a simple property of the diamonds: areas that share a boundary have opposite stacking configurations.
This algorithm can be slightly modified and applied to CAD and CAE environments.Here, both pattern generator and FE models are implemented as parametric codes in ANSYS Mechanical Parametric Design Language (APDL).The full model, already meshed, is depicted in Fig. 3. Information on the pattern is also observed (purple-and blue-colored elements) along with two rigid areas representing the compressive platens (in red).An 8-node quadratic element in its degenerated form (6-node triangular element) based on Mindlin's theory is used for simulating the shell and the contact behavior.The contact algorithm used is based on a surface-to-surface formulation along with the augmented Lagrangian algorithm.Furthermore, the model incorporates geometric non-linearity and the problem is solved implicitly along with the Newton-Raphson method.
The bottom platen (see Fig. 3) is fully clamped.The top platen can only move vertically and in this direction, the displacement is applied to mimic the actual testing.The reaction forces on the top-edge nodes are monitored throughout the simulation and used to compare with the experimental response.The elastic material properties used as input in the models are presented in Table 1.No direct boundary condition is applied to the cylinder since the first step of the convergence analysis is used to find the contact between the platens and the specimen.This initial evaluation can be considered a preliminary analysis of the contact, performed to avoid the need for stabilization or inertia relief techniques.The next stage of the simulation begins when the top platen finds the contact surface of the cylinder.

The damage model
A damage model based on instantaneous stiffness degradation is employed.Before the degradation process takes place, the damage must be triggered.Hashin failure criterion [34] is chosen for damage initiation.

Table 1
Material properties used in the simulations [33].The 123 material coordinate system is oriented following the fiber direction, as illustrated in Fig. 2(j Hashin criterion identifies four different failure modes: tensile fiber, compressive fiber, tensile matrix, and compressive matrix failure, as follows [35].
where   , with ,  = {1, 2, 3, 4, 5, 6}, correspond to the stresses, in which 1-axis is considered to be in the fiber direction.The parameters    ,    ,   , and   are always positive and represent the strength allowables required to trigger a particular damage when greater than 1.The description of  1 ,  1 ,  2 ,  2 ,  4 , and  6 parameters are presented in Table 1.
An intralaminar material property degradation model is used.Here, a degradation factor is multiplied by the material stiffness component after any damage initiation criterion is met.In other words, a ply discount mechanism is used to degrade the stiffness matrices.
Initially, a displacement  is applied on the top compression platen and transmitted to the composite cylinder via contact, as previously explained in Section 2.2.Throughout the analysis, this displacement incrementally increases and, for each step, the equilibrium is sought.
When any of the failure criterion (Eq.( 1)) is satisfied and damage is triggered in any finite element, the elastic properties of this element is degraded following the specified damage allowable related to that failure criterion met.Damage keeps progressing through the FE mesh until the cylinder completely loses its structural carrying capacity.This means that several failure modes may happen simultaneously and different finite elements may fail in different manners.This happens because, after damage initiates, the updated elastic properties re-distribute the stresses and consequently cause the failure of more elements at the same time increment or not.This procedure is repeated  times until the final failure is reached or the total displacement is applied.Thus, since the numerical evaluation is a simulation of the experimental test, the analyses were performed on an intact specimen up to a pre-defined maximum displacement, applied monotonically, simulating a quasi-static experiment.
The damage progression is applied by a slight modification of Hooke's Law.The nominal Cauchy stress, , written in Voigt notation, is obtained from the product of the strain,  (also written in Voigt notation), and a damaged elasticity tensor, , which is a function of the damage parameters under plane-stress stress state [35][36][37].This tensor is written as where  = ,   ,   , and   correspond to the fiber, matrix and shear damage variables defined by in which    ,    ,   , and   are damage allowables associated with the fiber tensile, fiber compression, matrix tensile, and matrix compression failure, respectively.For this particular model, the vector is written as follows This vector groups the parameters that define damage progression and, hence, need to be properly determined (see next Section) for achieving accurate predictions.Thus, once the damage is triggered (Eq.( 1)), the stiffness is reduced based on the damage allowables.

Finite element model updating (FEMU) procedure
A single-pass FEMU used here stands as a simplification of the usual methods used in vibration analysis and is adapted to find the best set of damage allowable values (presented in Section 2.3) to degrade the stiffness matrices of the structures after the damage is triggered.The main purpose of using FEMU methodology here is to fine-tune the damage allowables so that the numerical analyses are more accurate with respect to the experimental results and this will then enable the evaluation of cylinders of different winding patterns.
The use of several experimental points of the equilibrium path during an actual test as input data in the norm does not guarantee the convexity of the problem, instead, it represents the calibration of the numerical model to correct hidden flaws.It also adds to the reliability of the model to deliver good results for other deformation cases than the ones used as input.A weighted norm is used here to define the distance between experimental and numerical data as where  determines the nature of the norm,   denotes the number of evaluation points of the norm, the superscripts 'EXP' and 'NUM' correspond to experimental and numerical, respectively, and  is the weight.The right-hand side of Eq. ( 5) is the discrete part of the lefthand side equation, in which th point defines the point of evaluation of all parameters -experimental and numerical as well as their respective weight.The last is determined as in which  TD denotes the standard deviation of the experimental data for a given displacement.This approach reduces/increases the importance of points with larger/smaller standard deviation on the norm at the th point of evaluation.
It is important to mention that the experimental and numerical data sets may be different and that points, where the force and/or displacement are monitored may mismatch.So, the procedure adopted is to: (i) create a linear space using   − 1 regular intervals; (ii) select the closest points of each data set (experimental or numerical) for the specified entry; (iii) perform a linear interpolation.The norm will then be fully defined.
It is also important to notice that Eq. ( 5) is valid even for a single curve comparison, leading to a non-weighted n-norm.Since the experimental data, numerical results, and points of the norm evaluation are placed in different positions, both experimental and numerical values are linearly interpolated given the norm position evaluation.
The cost function is then set as   and the constraint minimization process is where  are the set of damage parameters (Eq.( 4)).
A gradient descent -steepest descent -algorithm is implemented to determine  that minimize   .Fig. 4 presents a flowchart of the proposed methodology.There, the norm (Eq.( 5)) is slightly modified to determine the gradient in a specific point of the cost function.First declaring that where   is the step size of the th design variable for determining the gradient of the cost function and should not be confused with the step size of the optimization procedure.Then, the norm is rewritten as

𝑁
is the norm;  ()  and  () are, respectively, the gradient and the search direction at   ( () ) while  describes the step size in  () direction with  and , which are Goldstein conditions parameters. f inal is obtained when a tolerance, , is achieved.in which  and  correspond to the loop and sensitivities indices, respectively.The gradient is approximated by finite differences as The steepest descent method defines that the search direction, , is defined by and is used to determine the optimization step size, written as where  ∈ (0, 1) and  ∈ {0} ∪ N.Both are related to the backtracking procedure of Amrijo's step size.The domain of  is a consequence of the required monotonic reduction of the step size (backtracking) and the requirement of remaining at the same search direction.These parameters ensure that the reduction in the cost function is sufficient.Another condition applied to the step size is related to the curvature of the cost function.Both correspond to the Goldstein conditions [38,39], described by where  ∈ (0, 0.5) is a constant.The left and the right inequalities regard the curvature and the Amrijo's step size [38,39], respectively.It is worth noticing that since  and  are collinear, however, with opposite directions, the difference  (+1,0)  −  (,0)  is always confined to negative values.
Furthermore, even if a design variable achieves values smaller/larger than the imposed limits, a varying restriction is adopted.Let us first consider 0 <  < 0.5 as a restriction value in which   ∈ [, 1 − ].The updating is essentially a monotonic decreasing function that depends on .The function used here is if Fig. 5 highlights the main aspects of the proposed methodology.The first step is an initial guess of the damage allowables.The points around the initial guess are also evaluated, and then the descent direction is determined, followed by the gradient of the objective function in that particular position,   ().For that, the norms between experimental and numerical responses are defined (objective function -Eq.( 8)).

Table 2
Geometric properties of manufactured and tested cylinders wound at [±60] [14].Differences in the thickness are related to the winding procedure [13].Pattern is defined by the pattern number/skip number.

Case no.
Thickness Length These norms,  (,0)  and  (,)  with  ∈ N ≤ 4, are obtained at  () and  () +   , respectively, to determine the gradient by a finite difference scheme.With the gradient vector and, consequently, the descending direction, Goldstein conditions (Eq. ( 12)) are used to derive the size of the step (backtracking).Then, another FE analysis is required, with an updated step size,  () +  () .If the Goldstein conditions are met, the last evaluation becomes the starting position of the new loop ( (+1,0)  ≐  (,0)   , see Fig. 5).If not, the step size is monotonically reduced (Amjiro's step size) until the criteria are met.These loops continue until the gradient norm is smaller than the predefined tolerance, achieving therefore an either local or global minimum.
The FEMU herein formulated searches for the optimum set of damage allowables given experimental and numerical data for radially compressed cylinders.The objective of the algorithm is to define the role of considering the pattern in wound structures and the best input parameters to predict their response under radial compression.The FEMU approach can also be used to cover unknown manufacturing parameters, such as the thickness.However, since the thickness was experimentally measured [14], it was not necessary to consider this parameter here.

Experimental details
The experimental setup and results are described in detail in a recent publication of the group [14].The cylinders have an inner diameter of 50.8 mm.The towpreg used to manufacture the cylinders is from TCR Composites, composed of T700-12K-50C carbon fiber and UF3369 epoxy resin.All cylinders have a winding angle of ±60 • .The other geometrical characteristics of each family are presented in Table 2, in which both thickness and length are shown with their average and standard deviation values.It is worth noticing that the difference in thickness comes from inserting more tows in the process in order to produce the desired pattern [13].
The experimental test follows the recommendations of the ASTM D2412 standard and the testing setup is shown in Fig. 6.Briefly, the specimen is placed between the platens and the top platen moves downwards while the bottom platen is kept stationary.The tests are performed in an Instron Universal machine model 3382 with a load cell of 100 kN under a cross-head speed of 12.5 mm/min.Five specimens free of burrs and jagged edges of each winding pattern (WP) are tested.

Results and discussion
The FEMU algorithm is applied along with the experimental data from Lisbôa et al. [14] (the supplementary material presents the raw data used) and the FE models developed here (Section 2).The meshes used in all analyses slightly differ (due to the patterns) and correspond to the mesh in Fig. 3, which is defined following a mesh sensitivity analysis.This analysis was carried out using the reaction force as the parameter (by a given displacement), in which the selected mesh (and finer meshes) produces no statistical difference in terms of reaction force.Considering that the mesh could affect the reaction force at every step of the simulation since the damage model changes element-wise the stiffness if a failure criterion is met, damage evolution was included in the convergence analysis by evaluating the norm for different meshes.For the norm determination,   = 300 and  = 2 (Euclidean norm), and the   is obtained from Ref. [14].At the beginning of the analysis, the time step is small so that the contact surface can be properly defined.In each simulation step, the displacement change is not higher than 0.5 mm.
For each pattern, five initial guesses are considered and the optimum damage parameters for each case are obtained.Due to the multi-valued solution, post-mortem analyses are performed to identify the most suitable set of damage allowables to predict the experimental response of the cylinders.Further details of the numerical results are provided as a supplementary material.

Damage analysis -numerical and experimental force × displacement curves
These initial guesses are chosen so as to consider different initial aspects of the objective function.For example,  (0)  and  (0)  correspond to a very severe and mild damage (see Eq. ( 2)), respectively.The notation on the damage allowables henceforward follows:  (0)   ,  = {, , , , }, for the initial set (Eq. ( 14));    for the optimized set; and   (without superscript) corresponds to variables of the optimization procedure, such as the evolution of the norm/damage allowables, for a particular set of variables.
It is worth presenting the predictions without considering the pattern in the model, which can be considered here as the baseline.However, since FEMU needs an experimental response to calibrate the numerical predictions, it is not feasible to do so disregarding the pattern due to manufacturing constraints of the FW process considering a helical winding [13,17].Thus, these baselines for WP = 1, WP = 2, and WP = 3 use the damage allowables obtained through FEMU (i.e.optimized sets) considering the pattern.Fig. 7(a-c) presents the experimental results (the median curve and its standard deviation -represented as an error band), the predictions for all three winding patterns with the five optimum damage allowables for each pattern, and the prediction without a pattern.Considering that the optimum five sets of damage allowables generated similar results, a simulation without a pattern is carried out using any of these optimum damage allowables from each winding pattern.Thus, the results shown in Fig. 7(a) represent the optimum design variables (damage allowables) for WP1, WP = 2 (Fig. 7(b)) and WP = 3 (Fig. 7(c)).
As can be seen in Fig. 7(a), the simulation disregarding the pattern shows satisfactory prediction in the linear part of the curve and also in the non-linear part until a displacement of 12 mm.In Fig. 7(b), the prediction for the no pattern case is worse than for WP = 1.After a displacement of 4 mm, the simulation keeps out of the standard deviation of the experimental response until the end of the simulation.The results of the simulation without pattern considering the WP = 3 experimental data are similar to the ones with WP = 2.The equilibrium path of this simulation (7(c)) is out of the experimental standard deviation area after the damage is triggered.This indicates the importance of modeling the winding pattern.Fig. 7(a)-(c) also shows the equilibrium path of the initial guesses considering the pattern (dashed lines with hollow symbols).It can be seen that, for WP = 1, the initial guesses   ,   , and   , are relatively close to the experimental results up to 12 mm displacement.Nonetheless, the same initial guesses do not yield a good fit of the experimental curve for WP = 2 and WP = 3.For those patterns,   performed better.14)); dashed line with × symbol corresponds to the modeling without the WP for a chosen converged set of damage allowables (using guess   for   = 1, guess   for   = 2 and guess   for   = 3).
Furthermore,   does not converge to the final displacement for all patterns.These results rectify the importance of modeling the pattern (i.e., the same damage allowable yields distinct behaviors depending on the pattern) and also validate the optimization procedure and the optimum damage allowable sets.Fig. 7(a) depicts the experimental curve along with standard deviation and numerical predictions for all five initial guesses for   = 1.In the linear-elastic region, all simulations are fully within the experimental range.This good agreement is not related to FEMU, but instead to a set of suitable material properties, adequate FE modeling, and especially the consideration of the winding pattern in the analysis.
Damage is triggered at a displacement level of 6 mm, i.e. any of the four different failure criteria (Eq.( 1)) is met (deeper insights into the failure criterion that triggered damage as well as damage propagation characteristics are provided in Section 4.3).A small deviation in the FEM results is observable for    > 15 mm.Nonetheless, the numerical results are inside the deviation in most of the experimental range.It is very difficult to accurately predict the whole load × displacement curve since the experiments show that the cylinder becomes slightly stiffer in the nonlinear region of the curves, and this cannot be predicted with the current damage model.Nevertheless, the numerical predictions are still within the deviation and therefore the simulations are in good agreement with experiments.It is worth noticing the high complexity in modeling these structures, since a large sample is transversely compressed, generating large displacements.Even with these drawbacks, all five guesses provide excellent predictions.
The increase in stiffness observed in the experiments (refer to the error band in Fig. 7(a)) after a displacement of 8 mm, is attributed to the reorientation around the two plastic hinges (i.e. at the two contact points with the platens) of tensile-loaded transverse yarns, which tightens the plasticized compression zones.Although this phenomenon is associated with strain, it can be applied to displacement and it happens in two phases: the elastic or elastic-plastic transition (when damage initiates), and the strain hardening zone.At large displacements, i.e. above 12 mm, the two lateral hinges (horizontal axes) might touch the compressive platens, creating additional internal supports, therefore increasing stiffness.A similar phenomenon has been reported by Calme et al. [40].
Fig. 8(a) shows the development of the norm with respect to the iterations of the minimization procedure and Table 3 presents both the  2 -norm and the design variables values at the end of the minimization procedure.The initial guesses reach distinct norm values and, hence, none of the design variables set are the same at the minimized value.This is a consequence of the dynamics of the failure process and how one failure mode can drive the others, implying a highly pathdependent problem.Thus, there is no guarantee that the value obtained by    for WP = 1 is a global minimum.The evolution of the design variables, from the initial guesses to the converged values, is presented in the supplementary material.
Furthermore, it is noticeable that for some initial guesses, the optimized set differs significantly mostly for cases with large values of damage allowables.The reason is that the objective function is relatively ill-behaved.The optimizer enters in a ''valley'' and ''travel'' through the objective function seeking for the minimum and thus the values of the design variables can have a significant change.In the supplementary material (SM2), one can observe even the oscillation that is characteristic of the steepest descent method in such problems.Hao et al. [41] proposed a multi-start gradient-based strategy using isogeometric analysis in which multiple initial designs for gradientbased (GB) optimization were determined by space tailoring method to enhance both convergence rate and efficiency.
The same procedure is then applied to the WP = 2 experimental data.The model is also updated to include the information on the pattern (Fig. 2(b)) and the results are presented in Fig. 7(b).All numerical curves are relatively distant from the experimental curve within 7 − 12 mm displacement.This takes place because in the numerical simulations the cylinders start damaging by matrix tensile and compression simultaneously earlier than in the experimental observations, therefore the degradation process starts prematurely.However, they achieve very good agreement in the following region, i.e. between 12 and 17 mm, being fully within the standard deviation.Besides, at the last region, the standard deviation is smaller than in the first one (penalizing the curves which are ''far'' in an Euclidean sense); and, in absolute values, the error is larger for greater forces, which benefits numerical responses that are closer to the experimental ones in large displacements.This indicates the importance of analyzing several samples either experimentally or numerically (through FEMU) to evaluate the numerical predictions.Fig. 8(b) shows the norm with respect to the number of FEMU iterations.It has the same characteristics of Fig. 8(a) for WP = 1.The very first steps reduce greatly the norm and then several steps are required to stabilize the damage parameters.As in WP = 1, the values obtained in FEMU for the initial guesses are different with a similar mechanical response in terms of force × displacement curve shapes, as observed in Table 3 and Fig. 7(b).Moreover, larger values of the norm, with respect to WP = 1, are expected given the aforementioned distances between experimental and numerical values.
The results obtained for WP = 3 are presented in Fig. 7(c) and, they are in excellent agreement with the experimental ones in the entire curve.Even with a smaller standard deviation region when comparing WP = 1 and WP = 2, the predictions are more accurate.The greater success of these predictions is also due to a more homogeneous evolution throughout the test.The variation of the norms regarding the number of iterations are presented in Fig. 8(c).Again, the norm decreases very quickly and then stabilizes.As one can see, in Table 3,

𝜌 𝑇
is bound to this function in iterations 4 and 5.However, with the aid of Eq. ( 13), the value of the minimum is not in the boundary of a feasible domain.
As happened for the other two patterns, FEMU has obtained five different optimum solutions (Table 3).The norm values are larger compared to the other two patterns, which is a consequence of two effects: a very small standard deviation in experimental data (which strongly penalizes any deviation from the curve); and differences between experimental and numerical curves at large displacement regions.This shows that for this pattern, the norm oscillates considerably and this is also an explanation for a large number of iterations in most evaluated cases.Nevertheless, the quality of the predictions is as good as for the other cases.In some cases, shown in detail in Figs.9(b) and 9(d), there is a ''jump'' of the ''line of damage'' at the helical-cross over.This is a very particular behavior that cannot be predicted by conventional FE models without any information on the winding pattern.These ''lines of damage'' are related to the resultant stresses that act on the cylinder during the testing: two at the region of the contact between platescylinder, with compression and tension on the external and internal surface, respectively; and other two 90 • from the contact region, with tension and compression on the external and internal surface, respectively.These lines are in regions where the maximum bending moment acts [42] and are much clearer in the internal surface.This ''jump'' is due either to local strength variations that made the crack change direction or to an initial damage that might have begun there and propagated to the line of damage.For WP = 1 the line slightly diverges and then returns, whereas, for WP = 2, it goes as a straight line up to the free edge.For WP = 3, some oscillation is observable but not ''jumps'', as the helical cross-overs are closer to each other when compared to the other patterns.Although the information on how damage initiated cannot be obtained (therefore the importance of the simulations -see Section 4.3), the final failure of all three families of cylinders present a combination of fiber-dominated and matrix-dominated mechanisms, since cracks along the fiber direction and transverse cracks are found.Furthermore, delaminations are also present especially at the plate-cylinder contacts [43].Delaminations and transverse cracks can be confirmed not only by the fracture specimens but also based on the successive load drops in the force × displacement curves (see Figs. 7(a)-(c)).

Damage characteristics -numerical analysis
Fig. 10 presents damage plots of the cylinders at the final load, where red and blue regions correspond to damaged and undamaged elements, respectively.The damage plots show both the outer and the inner layers of the cylinders.In all simulations, regardless of the optimum guess, failure initiated at 30% of the applied load by matrix tension at the contact with the bottom platen.Then, matrix compression failure mode is triggered.It is noticeable that the damaged regions are very similar for the different patterns, except for the contours.In these, the pattern information -helical cross-overs -can also be seen as a modifier of the damaged contour.
It is noticeable that at the top and bottom of the cylinders (the contact regions with the platens), the damage in the fiber is dominated by both tension and compression.The same occurs for the matrix damage at the sides (approximately 90 • away from the contact area).The stress distribution of WP = 3 with     is shown in Fig. 11 in each layer.Due to the nature of the element, there are differences in the stresses on the top and on the bottom of the layer (in Fig. 11 is shown as outside and inside of the deformed cylinders).One observes that, inside a layer, given the large and complex forces acting on the cylinder, there exist both tensile and compressive stresses at levels that cause/trigger damage.That is the reason for regions having both damages allowables triggered, despite the loading being monotonic.Furthermore, by the stress results, the pattern influence can be even better evaluated.Similar results are observable for other damage allowables set and patterns.
Fig. 12 depicts the information of damage of the cylinder shown in Fig. 11 and increment the information of the triggered damages.Different from other damage plots herein presented, this one considers all damage allowables by defining: no damage, partial damage, and complete damage, whose concept are no damage triggered, some damage triggered and all damage triggered.As aforementioned, due to the complex loading case, types of damage can be triggered together.This is also observed by Banat and Mania [44], where a profile is subjected to compressive loading and regions with all damage allowables triggered are observed.
An interesting result is that, for all patterns at the platen-shell contact regions, the failure is dominated by both matrix and fiber failure modes.Fig. 13 shows the progression of the fiber and matrix tension for 50%, 60%, 65%, and 70% of the final load.For the fiber (   ), the damage over the cylinder is barely noticed with 50% loading.With 60%, a line of damage appears in the inner surface with a small region on the outer surface close to the free borders due to the anticlastic curvature, an effect of free-border boundary condition.Then, with 65% loading, the damage is already in both inner and outer surfaces.The damaged area slightly increases for 70% loading.For   , the damage due to tension at the outer surface appears within 50%-60% loading and the same characteristics of the evolution of    are observed for   .Key information from Fig. 13, which is also observed for the other winding patterns, is the contribution of tensile stresses at the contact of the cylinder with the upper platen, which should be dominated by compressive stresses.At a loading level of 50%, two small areas show tensile matrix failure -which is due to edge effects since they are free.At 60% loading, the cylinder/upper platen region presents a row of failed (matrix tension) elements along the cylinder length.Above 65% up to final failure, the amount of failed elements in that area increases.Cylinder/platen friction also has a significant contribution to it.Nevertheless, the cylinder/upper platen contact is dominated by  compressive stresses, as expected.As can be seen in the Supplementary Material, matrix and fiber tension dominate the four main damaged areas, i.e. bottom and top contacts with the platen and the two sides, whereas matrix and fiber compression failure modes dominate at the upper and right-side of the cylinders.Fig. 14 shows the damage evolution for WP = 2 (guess   ) within 40%-70% loading.The same characteristics concerning the damaged regions can be observed here (tensile damage parameters in both bottom and top layers) and they appear with a similar absolute load (range between 500 N and 600 N).It is also noticeable that the damage away from the platen-shell contact starts around 50% of the applied force for   at the helical regions.This conclusion can only be drawn using the modeling approach presented here, i.e. taking the winding pattern into account in the analysis.Otherwise, the stresses would be underestimated since the helical cross-over regions may represent stress concentration areas.Therefore, failure is more likely to be triggered in those areas instead of a regular laminate area.Although the crossover regions can be interpreted as stress concentration areas, they can, nonetheless, postpone crack initiation and propagation, as observed in Ref. [14].
The damaged regions increase significantly within 50%-60% loading.However, this is more pronounced for the matrix damage failure mode, since even at a load level of 40%, damage elements are seen at the four main failure regions of the cylinder, i.e., damage propagation by matrix tension occurs faster than the other failure modes.A few elements are damaged by fiber (compression/tension) at the helical cross-over regions at this loading range.Considering that a loading level of 70% is near the final damage pattern, and therefore the damaged regions at this load step are shown, the failure pattern is mainly dominated by these modes, in the following sequence: matrix tension > fiber tension > matrix compression > fiber compression.These observations are further supported in the Supplementary Material.
In all analyzed cases (regardless of the pattern), the aforementioned damage progression occurs (see Supplementary Material).The matrix starts to fail in tension at the bottom layer close to the plate-cylinder contact.Then, at nearly the same load -around 50%, different damage mechanisms are triggered: matrix in both tension and compression at the contact areas, and, in fiber tension at the cylinder-upper platen contact region.Until this loading stage, the contact region plays the most significant role in the progression.After that, the most significant differences among the patterns arise.In WP = 1, the damage at the lateral of the cylinder appears only after a load step of 70% (at 70%, very small regions -a few number of elements at the helical region -are damaged).This behavior takes place slightly different for WP = 2 and WP = 3.For WP = 1, considerable areas are already damaged (matrix in both tension and compression) at 60% loading.Fiber damage (both compression and tension) can also be observed in some guesses and in very small regions close to the helical cross-over regions.With 70% loading, all cylinders have a similar damage pattern compared to the one at the final load step, only small variations related to the pattern lines and damage sizes are observed.The differences between the WP = 2 and WP = 3 concern, essentially, fiber damage.For WP = 2, this damage mode is considerable while for WP = 3, only regions close to the helical cross-overs are damaged.It is important noticing that these characteristics encompass all guesses, since they have similar behavior.As mentioned before, these differences in the damage  evolution characteristics for different patterns are only observed with the current approach.This also corroborates with the results found by Lisbôa et al. [14] regarding the maximum load and the absorbed energy of radially compressed filament-wound tubes.
The post-mortem analysis presented and discussed in Section 4.2, corroborates the numerical findings.One can see compressive damage at the fiber in both top and bottom layers.Regardless the optimum guess, the final numerical results are almost the same: differences are found between 50% and 70% of the loading.It is then concluded that any of the set of optimum guesses found here through FEMU approach provides a suitable prediction of the damage initiation and propagation of the cylinders.

Conclusions
A finite element model updating (FEMU) type approach is presented to find the damage allowables that best predict the mechanical response of composite cylinders under radial compression.The typical mosaic pattern formed during the filament winding process is taken into account in the models.Experimental tests were carried out with two purposes: (i) as input for the FEMU procedure, and (ii) to compare the final damage pattern with simulations.The optimization algorithm within FEMU, whose design variables are the damage allowables, is implemented in order to minimize the norm between experimental and numerical results.This norm is weighted by the standard deviation of the experimental data.Three winding patterns are considered (1/1, 2/1, and 3/1), each of them with five pre-selected initial guesses considered in FEMU algorithm for determining the damage parameters that best predict the experimental data.
The pattern plays a key role in the radial compressive response of filament wound cylinders and, therefore, this manufacturing characteristic should be taken into account for generating realistic stress fields.The helical cross-over regions can either postpone crack initiation and propagation or represent stress concentration areas, since stress peaks are found in these regions.In general, different winding patterns and optimum guesses produce similar damage initiation and propagation mechanisms, in which failure is triggered by matrix tension at the cylinder-upper platen contact, followed by matrix compression at both top and bottom contact areas with the platens.The final damage pattern is dominated by a combination of these failure modes at these regions and the sides, regardless of the pattern and optimum guess, as expected.This is the natural response of fiber-reinforced cylinders under radial compression regardless of the winding pattern.Nonetheless, differences in maximum load, stress, strength, and damage propagation were observed for different WPs and damage allowables set.Qualitative differences in crack propagation were also seen in the post-mortem analysis of the fractured cylinders.Therefore, the experimental setting was suitable to exploit the WP effect, and the adopted numerical methodology was effective in predicting the experimental response of the cylinders.
It should be pointed out that the optimum sets of damage allowables found here are only valid for loads and BCs of the current FE model.The multiplicity of solutions herein obtained and the associated dependence on the initial values could be attenuated by using surrogate models (or metamodels).Also, the lack of physically meaningful parameters suggests caution in the extrapolation of the obtained results to conditions far from the range of values herein explored.In that case, a new set of experiments could be used for checking purposes.The proposed FEMU framework, however, can be easily adapted to any arbitrary FE model, for multiple load cases and BCs.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .Fig. 2 .
Fig. 1.Pattern generation process, where: (a) first stroke, (b) second stroke, (c) WP = 1 (top) and WP = 2 (bottom), (d) full coverage for WP = 1 (top) and WP = 2 (bottom), and (e) the detailed diamond-shaped form.The light and dark gray correspond to a forward and a backward stroke, respectively.The red rectangle, in (e), defines the circular cross-over, the green rectangle the helical cross-over, and the blue triangle the regular laminate region.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Meshed geometry highlighting three different regions: rigid compressive platens in red, +∕ −  in blue, and +∕ −  in purple.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig.4.Flowchart of the optimization procedure, where  defines the damage allowables and, in the optimization procedure, the design variables.(,)

Fig. 7 .
Fig. 7. Experimental (error band) and numerical force × displacement curves for (a)   = 1, (b)   = 2, and (c)   = 3. Continuous lines with solid symbols represent the optimized (   ) damage allowables set (through FEMU) with WP modeled; continuous lines with hollow symbols (  ) represent the initial response (damage allowables of Eq. (14)); dashed line with × symbol corresponds to the modeling without the WP for a chosen converged set of damage allowables (using guess   for   = 1, guess   for   = 2 and guess   for   = 3).

Fig. 8 .
Fig. 8.  2 norm of the optimization procedure with respect to the number of iterations,   , the initial guesses,   , and the winding pattern: (a)   = 1, (b)   = 2, and (c)   = 3.Values close to the arrows correspond to the  2 norm at the first iteration.

Fig. 9 .
Fig. 9. Images of failed specimens for: (a-b) WP = 1, (c-d) WP = 2, and (e-f) WP = 3.The first column presents the failures at the contact between the cylinders and the platen while the second column shows the inner surface of the cylinders, highlighting the lateral compressive failures.

Fig. 10 .
Fig. 10.Damage plots at the maximum compressive load for all investigated failure modes.The cylinders are plotted in the isometric view and show a similar damage response in both top and bottom regions (the contact with the platens).Final results for the patterns   = 1,   = 2, and   = 3 are associated with guesses 'a', 'c', and 'e', respectively.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 .
Fig. 12. Damage on the WP = 3 cylinder with    .Definitions: no damage (none of the failure criteria met), some damage (some of the failure criteria triggered), and full damage (all failure criteria met).

Fig. 13 .
Fig. 13.Damage propagation -   and   -from 50% to 70% of the total applied force for the winding pattern (WP) of 1.The cylinders, plotted in isometric view, have a similar mirrored damage response in both top and bottom regions -the contact with the platens -as well as in both lateral (left and right) sides.

Fig. 14 .
Fig. 14.Damage propagation of all failure modes from 40% to 70% of the total applied force (1100 N) for the WP = 2.The cylinders, plotted in isometric view, have a very similar mirrored damage response in both top and bottom regions -the contact with the compression plates -as well as in both sides.

Table 3
2 -norm and design variables values at the end of the minimization procedure for   = 1,   = 2,   = 3 and each initial guess (0)