Progressive Fatigue Failure Analysis of a Filament Wound Ring Specimen with a Hole

A progressive FEA fatigue degradation model for composites was developed and implemented using a UMAT user material subroutine in Abaqus. Numerical results were compared to experimental strain field data from high frequency Digital Image Correlation (DIC) of split disk fatigue testing of pressure vessel cut outs with holes. The model correctly predicted the onset and evolution of damage in the matrix as well as the onset of fiber failure. The model use progressive failure analysis based on the maximum strain failure criterion, the cycle jump method and Miner sum damage accumulation rule. A parameter study on matrix properties was needed to capture the scatter in strain fields observed experimentally by DIC. S-N curve for the matrix material had to be lowered by 0% to 60% to capture the experimental scatter. The onset of local fiber failure had to be described by local S-N curves measured by DIC having 2.5 times greater strain than that of S-N curves found from standard coupon testing.


General
Material direction specific Material direction, with examples: Fiber direction - 11 Matrix direction - 22 In plane shear direction - 12 Through thickness direction - 33 General sign of peak/max value ̂ ̂-

Introduction
There is an urgent need for composite pressure vessels that can safely and economically transport hydrogen at 700 bar [1,2]. The technology and design standards exist, however, cost is high due to very strict testing and acceptance requirements [3,4,5] even at lower pressures. For the acceptance tests, a perfect structure is assumed. However, during a vessel´s lifetime, small damages such as a minor impact damage may occur from use. It is currently an unknown how much damage can be tolerated in the vessels due to unknown fatigue resistance. Damaged vessels are replaced by new ones, which is very costly, especially for large vessels.
Today´s pressure vessels have a static strength exceeding the design pressure of 700 bar by a factor of about 2.5 or more as required by the design standards. The factor was also identified by Berro et. al. in the OSIRHYS IV project [6]. Uncertainties of the effect of the presence of damage are largely related to mechanical fatigue. Numerical analysis in combination with well-chosen experimental data is the key to better understand how damage and fatigue may affect the mechanical performance and strength [6]. In turn, better numerical models may answer how much wear and damage can be tolerated on in use vessels, avoiding early and costly decommissioning as well as reducing costly testing requirements of new designs Fatigue in composites cause complex progressive damage development that sets it apart from more conventional materials such as steel. Instead of damage and crack propagation occurring over a very short cycle span towards end of life, fatigue damage in composites occurs steadily over the whole lifetime, gradually changing the structural behavior and redistributing loads [7]. The dominating mechanism for changing strain fields under fatigue is the development of matrix cracking over time [8,9,10]. Matrix cracking/matrix damage is seen here in its widest meaning, including cracks in the polymer part of the composite, delaminations and fibermatrix debonding. Developing matrix cracks change how forces are distributed between the load bearing fibers and cause the strain fields to change. Initiation and propagation of the various forms of matrix cracks is a complex phenomenon. When using finite element models to model matrix cracking, both initiation and propagation need to be predicted, including the propagation direction. Recent developments in such models have managed to satisfactorily predict matrix damage dominated fatigue damage propagation in laboratory test specimen having simple geometries and known direction of crack propagation, e.g. Travesa [7] and Nixon et. al. [11] based on the method suggested by Harper et. al. [12]. Attempts of simplifying the matrix crack growth by smearing matrix cracking over a larger region and modeling it by plastic behavior were reported for the static case by NASA [13], Flatscher et. al. [14] and Gagani et. al. [15]. It is however difficult to tell whether the plasticity approach matches experiments only for the particular geometry of the specimen investigated or whether it is a general way to model the material.
There are currently some commercially available composite fatigue numerical frameworks available, most notably FEMFAT [16] and Fe-safe [17]. While the models offer simple and fast fatigue evaluations, they do not include progressive damage and do fatigue analysis based on a static solution. The models are only tested on simple lab coupon specimen and lack experimental comparisons with local strain fields. Recent developments have expanded a modified smearing approach into fatigue, most notably by Koch et. al. [18]. This takes progressive fatigue damage into account. The inherent discontinuity in the stress/strain relationship upon matrix cracking and fiber failure yields challenges in finite element analysis when attempting to reduce the material stiffness at integration points during a constant load. Koch found that a cycle jump approach with constant properties for each loading cycle had to be applied. Degradation was carried out between the jumps according to the size of the cycle jump. A similar approach is used in this work. While Koch compared the model to global experimental data, this work aims to estimate the local experimental strain data as obtained from digital image correlation monitoring of the modelled test specimen. The discrepancy between local and global properties was most notably highlighted by Sevenois et. al. [19]. Sevenois argued that matrix crack initiation and propagation on the local scale happens long before catastrophic failure of typical composite fatigue test specimen. As noted by others [20,21,22], matrix voids affect the fatigue properties to a great degree, also found in the presented work here. Matrix voids induce matrix cracking on the micro level. It is essential to establish when matrix cracking is initiated on the local level to estimate fatigue life in a finite element model. However, so far this has not been taken into account and global cycles to failure for the material used as input for local properties in most fatigue models. Sevenois also highlighted the scale problem. For example, atomistic bonds break long before any typical matrix crack is initiated in the structure. In engineering terms damage causing changes to the structural behavior on a component level is important. In this work the scale of interest is that of strain field changes observable through standard scale in DIC (Digital Image Correlation). Sevenois also argued that detailed local models and sophisticated failure criteria come short of modelling anything but a perfect structure. Matrix voids and variations in fiber volume fraction throughout the structure will make the real damage development complex. In this work, all the above has been acknowledged and addressed through parameter studies on matrix material properties. The resolution of the DIC method enables comparison between model and experiment on a very detailed level, taking into account local variations in material properties. For the fiber properties, local material properties were successfully found using a damage calculation method on the DIC data. Good correlation was found between model and experiment using the local fiber properties.
This study suggests a simplified modeling approach that could be sufficient for understanding how local strain fields develop under fatigue loading in the composite material and how this may affect the global behavior. The modelling approach was implemented as a combination of cohesive surfaces and UMAT (User material subroutine) in Abaqus. A few simplifying methods were used: • Micro matrix failures were modeled using a continuum damage approach as changes in the stiffness matrix without directionality of the cracks. • Macro (visible) matrix failures were modeled as discrete cracks permitted to propagate along predefined surfaces when certain strain states are met. They were used for through-the-thickness cracks in a ply and for delamination. Only selected macro cracks were modelled. • Discontinuities in the stiffness due to crack growth were modeled using an on and off loading approach in combination with simplified cycle jumping [12,11]. • A range of polymer matrix properties were modeled to investigate the natural variations of material properties. Figure 1 Experimental setup with geometry of test sample and layup, fibers not to scale. Figure 1 outlines the experimental split disk setup with vital dimensions included. The two holes, located at both sides of the disk, were designed to simulate extreme damage in the composite, and a tension-tension load control fatigue loading with an R-ratio of 0.1 was imposed.

Experimental setup
The test rings were cut from filament wound specimen with a layup of hoop/axial/hoop/axial fibers as seen in cut A of Figure 1. The fiber orientation angles were [±89°2, ±15°1, ±89°2, ±15°1], typical for filament wound pressure vessels [6]. Figure 2 shows winding of the first axial layer. The fibers were HiPerTex W2020 glass fibers from 3B [23] and the resin was Epikote MGS RIMR 135 mixed with curing agent Epikure RIMH 137 [24]. The thickness of each layer was found through microscopy of the cross sections and is described in Figure 1. The holes were cut with a composite specific milling tool; 40200-HEMI produced by Seco Tools. Figure 3 shows the test setup in the laboratory. The cameras were timed against the load signal and images were taken at peak load with a frequency of 50 cycles between each image. The DIC data was post processed using Vic-2D from Correlated Solutions. Python scripts developed by the composite group at NTNU were used for extracting results and performing data analysis. The resolution of the processed data was 4 points per mm 2 . Adequate resolution was found through a sensitivity study trying several different resolutions.

Failure Criteria
General approach Failure criteria predict the onset of defined failure mechanisms. Once a state variable (here strain) in the material reaches the limit set by the failure criterion, the constitutive properties are degraded. The scale at which the failure criteria apply also needs to be defined. The scale on which the criteria operate can range from that of the component (load displacement curve from a test machine) and down to the atomistic level. A finite element analysis as used in this study typically represents the mesoscale.
This section describes the failure mechanisms and failure criteria used on the different scales and how the changes in properties are reflected in the constitutive relations. The following failure mechanisms were modeled: • Micro fiber failure  Modelling micro damage with UMAT Fatigue was described by a strain based S-N curve in the log-log format: where ̂1 and are the intercept and slope of the S-N curve respectively and subscript denotes the strain components as defined in Figure 5. The number of cycles to failure was defined as: If the static strain to failure ̂ is less than ̂ a cut-off was added to the S-N curve.
Since the strain field changes with the development of partial damage under cycling, the Miner sum was used, expressed in equation 3.0. Instead of only calculating the Miner sum, which is a non-physical number, it was useful to define an exposure factor . The factor describes the ratio of load to material strength at the applied strain after a certain strain history described by the Miner sum. Failure happens when The exposure factor enables the introduction of a partial degradation of the material before the particular failure mechanism has happened. The exposure factor varies for the different strain components since applicable strains and S-N curves may differ. The full set of exposure factors are given in equation 6.0. Note that the matrix dominated strains 22 and 33 are influenced by the strain 11 in the fiber direction due to the Poisson's effect as follows from equation 6.0. The opposite coupling is however neglected since this coupling gives neglectable differences in 11 due to the stiffness difference between longitudinal and transverse direction. As the failure criteria were defined such that failure in one matrix associated component or plane (22, 33, 12, 13 or 23) gives failure in all matrix associated component, the coupling between 22 and 33 was also neglected.  The Young´s and shear moduli are denoted and and the Poisson´s ratios are . How the stiffness reduction factors were changed is described in Table 1. The choices behind the factors and their interaction will be explained below. Table 1 How reduction factors are changed depending on failure criterion. As can be seen in Table 1, the factors changed the stiffness extensively, from full stiffness to 10% stiffness. To ease the change, the stiffness was gradually reduced from the intact to the failed value over an exposure factor span from 0.8 to 1.0, as expressed in equation 9.0. This is schematically shown in Figure 6 for a reduction from 1.0 to 0.1. Particularly for elements with an exposure factor fluctuating about 1.0, the softening ease the iterative scheme, avoiding distorted elements with a large internal stiffness difference. The largest exposure factor of components 22, 33, 12, 13 and 23 was used as basis for reducing the constitutive properties in these components, in accordance with Table 1.

Modeling macro damage with Abaqus contact definitions
Delaminations and macro shear cracks penetrating the entire thickness of a ply were described as macro damage, see Figure 4. Delaminations may occur between all the layers. The layers were therefor modeled separately with Abaqus' cohesive surface contact defined on the interfaces. Macro shear cracks were known to develop and extend at four locations in the test specimen. They initiated at the equator of the hole and extended in the loading directions along the shear strain bands. The locations can be seen as white gaps in the DIC strain measurements shown in Figure 17 extending out tangentially from the hole. For the 1/8 FE model the shear crack is evident in the mesh in Figure 13 in between the two mesh regions in the surface layers. Modelling the location of the cracks directly into the FE models significantly simplified the modeling approach.
Abaqus' cohesive surface definition was used with a triangular traction-separation response, illustrated in Figure 8 with the values from the normal direction ( ). The separation is defined in mm between the two surfaces. The triangle is defined by an initial elastic stiffness ( ) defined by the characteristic element length as given in equation 11.0. Upon reaching the maximum stress ( ), the contact stiffness is reduced according to the fracture energy . The elastic stiffness ( ) was defined according to suggestions by Diehl [26]. Diehl found that the contact stiffness in the elastic regime was best described using a factor of 0.05 on the characteristic element length as described in equation 11.0. The characteristic length was defined as the mean element length, d, along the shear crack in the finely meshed region in Figure 14, giving a of 0.255 mm. Based on suggestions by Perillo [27] on the used material, the Benzeggagh-Kenane (BK) mixed mode behavior was used along with an energy mixed mode ratio with a BK exponent of 1.4.
To account for fatigue in the cohesive surface definition, the contact properties were reduced. The separation stress and the fracture energy were scaled with a factor of 0.6. 0.6 was chosen based on the graph in Figure 7 showing the S-N curve for the matrix shear with the globally reduced property in red. Details on how the S-N curve was found are given in the Material properties section. The red line is set close to the convergence of the matrix S-N curve within the cycle span of interest, from 0-100000 cycles; as the specimen failed at 127814 cycles. Any strain above this threshold will lead to failure within very few cycles relative to the cycles to failure. The scaling method was done as there is no fatigue definition built into the standard delamination crack definition of Abaqus and it was not possible to the authors' knowledge of running a separate user subroutine on contact properties and micro damage UMAT simultaneously.

Cycle jump method
Overview Fiber and matrix material have little or no yielding once the failure criterion is met. Upon failure they cease to carry the load and the stress-strain curve is discontinuous as a result, shown schematically in Figure 9. The ideal fatigue-degradation material subroutine would iterate the stiffness based on the applied cycles and stress. The local stiffness for each cycle would be decided by a changing and discontinuous stress/strain curve as displayed in Figure  9. As cycles increase, ̂ and ̂ would become lower. Upon reaching ̂ and ̂ the routine would have to be able to handle a negative tangent stiffness. Negative tangent stiffness is in theory impossible using conventional iterative schemes. The fatigue material subroutine outlined above is therefore impossible. The problem was avoided in the presented work by using a cycle jump method, similar to that explored by Harper and Koch [12,18].

Implementation
The cycle jump approach has two distinct phases; (i) loading and offloading and (ii) cycle iteration, similar to that explored by Koch [18]. In the loading phase, peak exposure factors and strain in all integration points are recorded. Stiffness is kept constant and not changed as the load is increased to avoid local negative tangent stiffness. When the peak load is reached, the offloading and cycle iteration phase is initiated. Here the Miner sum is calculated on the peak strains from the loading phase and the exposure factors are changed with cycles. The stiffness is still kept constant in this phase, as the structure offloads to the strain and stress state before loading. Upon initiation of the next loading phase the stiffness is changed according to the new exposure factors from the cycle iteration and the procedure is repeated.
In case of simulations with varying fatigue load, this can be achieved by loading to different loads.
Referring to equations 6 to 8, , and are kept constant during the loading phase and changed; following the strains in the integration points. During the offloading and cycle iteration phase, is kept constant at the peak strains recorded in the loading phase. is therefore free from the strain in the integration points and the structure is left to unload by itself. Now , is changed according to the cycle in the given iteration. Upon initiation of the next loading phase is changed according to equation 9.0.
Due to the cohesive surface contact definition, the unloading phase may yield singularities in the stiffness matrix. To overcome this, damping was introduced in this phase to have the structure relax without causing singularities. Figure 10 shows the cycle jump approach explained schematically. The damping is evident in Figure 10 as the displacement curve lags behind the loading curve in the offloading and cycle iteration phase. As the iterative scheme is dependent on the structural response from the damping, the numbers of cycles where results are available in the offloading and cycle iteration phase are not fixed, however at the end of this phase/step they are. The cycle jumps were chosen based on experimental data in this work.
For complicated models the runtime is long even with high memory and high CPU capacity computers, therefor the cycle jumps has to be placed with care and at critical points in the load history. In most cases this is at the start and end of the components lifetime. At the start there will be initial matrix cracking and at the end of life there will be extensive fiber failure [28]. In between, the strain distribution will be relatively stable and cycle jumps may be relatively big as a result. Given a load to displacement history from an experiment, it is therefore advisable to place the loading phases at cycles where changes occur in the load to displacement history. Due to the long runtime it is time consuming to study several different cycle jumps. This study implemented only four well-chosen cycle jumps. Despite of this coarse approach the method gave a good indication of where damage initiated and how this damage affected the strain distribution and material behavior over time. A further expansion of the method would be to include an automatic cycle jump procedure. This could be done by assigning a maximum damaged volume and have the loading phases occur when damage extend over this specified volume. As such, the method would be completely independent of experimental data. This was done with success by Koch et. al. [18] on simple models, however it was not explored in this work due to the high computational cost. Alternatively, cycle jumps could be set at decreasing intervals until the results of a few cycle jumps converge. This procedure would be easy to implement, but would also require high computational times using a model as big as in this study.

Material properties
Material properties of an orthotropic ply with transverse isotropic behavior were used for FEA modeling. A summary of all material properties is given in Table 2.
Two main assumptions were made: • The plies have transverse isotropic behavior as per classic composite material models.
• Only material properties in the tensile direction were considered. The ring on the split disk experienced some compressive stresses. However, since these stresses were small and not in critical regions for the structural integrity, they were not considered and simply modeled with the tensile fatigue data.
Most of the material properties were measured in our laboratory during previous projects from standard coupons made of the same glass fiber and epoxy matrix. Data were obtained for unidirectional flat materials. Filament wound materials have a curvature and are strictly speaking not unidirectional. Properties were scaled to apply to filament wound material using a FEA approach by Perillo et. al. [29]. The static properties in the direction of the fiber and matrix ( 1 , 2 ,̂1 ,̂1 ,̂2 ,̂2 ,̂1 2 ) were linearly scaled to account for differences in the fiber volume fraction between the filament wound material in this study and that of Perillo. Table 2 shows the material properties and the methods used to obtain them. Only one static property was measured for this particular study; the maximum static strain to failure of the fibers. Compared to the original maximum static strain found by Perillo of 22150 microstrain [29], the value reported here of 40 000 microstrain is considerably greater. The value was found from DIC strain measurements taken from static split disk tests [30]. 40 000 microstrain was the highest strain recorded upon catastrophic failure. This maximum strain deviated by a margin of almost two from strain at the exact point of failure of 22150 microstrain, similar to what was found by Perillo from standard coupon testing. The weakest point had similar properties to data obtained by coupon testing, as coupon tests measure the weakest part of the sample. The DIC data showed however that the local strains can be much greater without causing failure. The greatest local strain (40 000 ms) was used for the fatigue analysis.
Fatigue properties were described by strain based S-N curves for the three in-plane material components of an orthotropic ply: fiber, transverse and shear components. Figure 12 shows the three individual S-N curves in a linear strain presentation. In all cases the S-N curves could be well described by the log-log presentation as given in equation 1. All testing was done for an R-ratio of 0.1 (tension-tension). This is the most relevant loading condition for pressure vessels being cycled between nearly empty and full.
Through-the-thickness shear (intralaminar shear) was measured on Short Beam Shear SBS specimen cut from the filament wound vessel with geometry according to the ASTM D2344M standard. The slope of the log-log S-N curve was found by linear regression according to equation 1.0. Figure 11 shows the data points and curve fit of the SBS testing. The slope of the curve is 0.051. Nearly the same slope (0.054) was found for the same constituents tested on flat specimens by Gagani [15].
Fatigue properties for in plane matrix cracking were not measured, but taken from the intralaminar shear data. The slope from the shear data was used for the tensile matrix curve with the origin from the static properties. Using the same slope in both tensile and shear S-N curves is controversial, however, it has been shown before that the slope in the tensile matrix direction is in general low [8], as was also found for the shear.
DIC monitoring of the SBS and split disk testing revealed a high discrepancy between local and global fatigue failure and a large scatter in matrix properties. To account for the variations in local properties a parameter study on the matrix fatigue properties was done and is reported in the Results section, changing the intercept/origin strain of the S-N curves. The parameter study made the exact knowledge of the matrix dominated S-N curves less critical, as the properties were changed in the analysis anyway to capture the experimental scatter.
Local fatigue properties of the fiber was obtained by using DIC data from testing of three split disk ring tests with a hole [30]. Strain based S-N curves were found by an iterative process. Local fiber failures were predicted using Miner sum calculations on the DIC strain data, as in equation 1.0 to 3.0. The predictions were compared with the measured failure of the samples. The S-N curve was changed until a good match between predictions and experiments was achieved. Details of the procedure are given in [30]. The S-N curve giving the best fit with the experimental data had a slope of 0.1, the same as usually measured for this type of material using standard coupon data [31,32,33,10,34,35,36,37]. However, the origin of the local S-N curve had to be greater than for typically obtained S-N curves to match the experimental data. The local S-N curve for fiber failure had to be shifted up by about a factor three compared to typically reported curves from coupon testing. This resulted in an intercept of the S-N curve for one cycle that is greater than the static strain found from DIC data, giving a cutoff of 30 cycles on the S-N curve. *See section FEA model for factor explanation

Geometry
The composite ring was modelled as a 1/8 model with the geometry in the schematic in Figure  1. The model can be seen in Figure 13. The layers (hoop/axial/hoop/axial) were modelled with five elements thickness each, as can be seen in Figure 14. Each layer was defined through the composite layup function in Abaqus. The layers were defined in this function as 10 +/-layers (+15, -15, +15…, -15 and +89, -89, +89, -89…, -89), essentially smearing the properties.
The ring was meshed using hexagonal eight node reduced integration elements, C3D8R and had 223734 elements. As explained in the Failure Criteria section, the layers had a cohesive surface contact definition defined between them to model delamination. The crack evident in the top hoop layer in Figure 13 and Figure 14 also had the cohesive surface contact definition as explained in the same section. The crack was only present in the hoop layers, while the axial layers were meshed as shown to the left in Figure 14.   Table 3 gives detailed constraint definitions with reference to Figure 15 for surface name definitions. The cylindrical coordinate system is defined with a radial and a tangential vector. Figure 15 Overview of constraints, surface names, reference point and the cylindrical coordinate system.

Experimental Results
This study compares experimentally obtained strain fields to FEA modeling of one composite ring with a hole as shown in Figure 1, tested with the split disk test method in fatigue. The strain fields in the vicinity of the hole were measured by DIC every 50 cycles. The sample failed at 127814 cycles, shortly after the last DIC frame at 127768 cycles in Figure 17. Figure 16 shows the cycles to displacement curve from the test machine. The dotted lines are the cycle jumps in the FEA model, these will be further explained in the FEA results section. Catastrophic layer failure/fiber failure in the individual layers happened over relatively short cycle spans indicated by sudden displacement jumps in the curve. The layer failures are highlighted with arrows and text in Figure 16.
There were four regions around the hole in the disk that concentrated strain. The regions are highlighted in Figure 17 over the contour plot of hoop strain close to catastrophic failure. The first fiber failure has already occurred as indicated by the white gap above the black rectangle. The gap was first visible at about 122000 cycles, roughly 6000 cycles before catastrophic failure of the surface layer. This large fiber failure led to the sudden jumps in the displacement curve visible in Figure 16. It initiated progressive fiber failure that ripped the layer across the region indicated in the black line slice in Figure 17 before it progressed to the other side of the hole.
Intuitively there should be only two regions that concentrate strain, one on each side of the hole at the equator. Macro splits (matrix cracks through the thickness of a ply) develop at these points after very few cycles and the splits grow rapidly in the load direction. The splits and bending of the material between the splits move the strain concentrations to the ends of the splits [41] at the four regions shown in Figure 17.
Provided a perfect material, the four regions around the hole should have equal strain fields throughout the test and equal damage development. Figure 18 shows the strain curves over the length of the four regions at selected numbers of cycles. Cycle 350 was the first recorded cycle by the DIC (cycle 1 was not recorded). The x-axis is in absolute values (note: no negative x values), 10 mm is at the edge of the hole and 24-26 mm is at the outer edge of the specimen. The curves are evidently not equal. The differences are due to variations in void content in the matrix, fiber density and layer thickness; these variations can easily be seen under the microscope, as shown in Figure 19.
It can be seen in Figure 18 that the strains measured by DIC have a quite pronounced strain concentration at the splits, at the 10 mm position, at a low number of cycles. With increasing number of cycles, the average strain increases but the strain concentration diminishes or moves even to the outer edge of the specimen. This effect is due to matrix damage spreading in the material changing the constitutive properties of the material [42].    Figure 17 for different cycles. Figure 19 Microscopy of the split disk specimen. As can be seen there are many voids (black) and a big variation in layer thickness over the cross section.

FEA results
The goal of the FEA modelling was to predict the correct strain field and damage development throughout the fatigue life of the specimen caused by the chosen failure mechanisms using their fatigue failure criteria with corresponding material properties. As shown in Figure 17 and Figure 18 the specimen had big variations in the shape of the experimentally measured strain field over the surface. To capture the variations, the FEA model was run for four cases with different matrix material properties and one case with the failure mechanism in the fiber direction disabled. Table 4 and Table 5 show the chosen variations of properties in the models. Model A had nominal properties as defined in the Material properties section. Model B-D had degradation on matrix properties down to 40% of the original values. In addition, the matrix damage effect on stiffness in the fiber direction was lowered down to 0.7 for the fiber for Model C and D, while Model E had the property degradation in the fiber direction disabled. The properties for the models are a representative selection of a parameter study that explored what correlated best with the range of experimental results. All models apart from Model A evidently deviate from the assumed physics of the problem, as defined in the Failure Criteria section. It was however interesting to explore what degradation was necessary to better capture the experimental results.
Variations in layer thickness were not included in this analysis as this is difficult to model correctly, but the effect was modeled indirectly as reduced matrix properties. The shear crack reduction factors were chosen based on the same evaluation as for models A to E.
FEA fatigue calculations were done using the cycle jump method as described previously. Figure 16 shows the cycle jumps over the cycle to displacement curve. The jumps were placed at interesting points in the cycle to displacement curve. Fiber/layer failure manifests itself as displacement jumps in the experimental curve. An increasing amount of matrix failure is expressed as the gradual increase in displacement between the fiber/layer failures particularly prominent over the first 10 000 cycles. The analysis cycle jumps were put before or at critical changes to the displacement curve. The cycle 1 jump is there to estimate the initial state, the cycle 10 000 jump to capture the initial matrix damage, the cycle 80 000 jump to capture the state just before first layer/fiber failure and the cycle 127 000 jump to capture the state just before catastrophic failure at 127 814 cycles.
Ideally there should be put more cycle jumps between 80 000 and 127 000 cycles considering the relatively big increase in the displacement. Due to analysis time this was however not feasible. The analysis time for one set of material parameters was six days on a fast 8 slot 64 GB RAM computer with an Intel Xeon W-2155 3.3 GHz CPU.   An example for comparing hoop strain fields calculated by the FEA with measured DIC data is shown in Figure 20 for Model D at 80 000 cycles. The most highly strained regions are qualitatively similar. The absolute values match however only at a few locations. This is to be expected since the experimental data show quite high variations in the four sectors around the hole. The FE model produce the same results in each sector due to assumed symmetries.
Running the FE analysis for each matrix cracking material model described in Table 4 and Table 5 allowed comparing the FEA against experimental data for different matrix properties. Figure 21, Figure 22, Figure 23 and Figure 24 show the lowest and highest experimental strain across the most highly strained cross sections outlined in Figure 18 compared to the experimental results for the four models A, B, C, D and E for cycle 1, 10 000, 80 000 and 127 000 respectively. The full colored curves correspond to colors of the regions in Figure 17.
There is generally good agreement between experimental data and FEA calculated strains for all but the 127000 th cycle. The models fall between the extremes from the experiment for the 1 st , 10 000 th and 80 000 th cycle. This shows that the chosen matrix material models represent the various degrees of random material variations of the matrix properties fairly well.    While the low strain curve from the experiment is relatively smooth throughout cycling ( Figure 21 -Figure 24), the high strain curve is not. The high strain curve has also got considerably greater mean strain. The strain curves when using matrix crack models C and D can be seen to have the same uneven shape and also a greater mean strain. The unevenness and greater magnitude are due to more matrix damage. It is interesting to see that the modeled strains start to fluctuate for the highly degraded matrix properties, even though the model treats the properties as the same throughout the model. The black curve in Figure 18 has the same tendency, which is the curve from the region with first observed fiber failure and also catastrophic failure. While the matrix properties are well described using the chosen envelope of degradation in models A-D, it is evident from the 127000 cycle curves in Figure 24 that strain fields after the bottom hoop layer failure are not as well described as when predominantly matrix damage is present as for cycles 1 -80000. The experimental strains are about 40% greater than the simulated strains after first experimental layer failure. Model D can also be seen to have a substantial fiber failure at 127 000 cycles, which the other models do not have. It is however not consistent with where the experimental fiber failure occurred, which was at about 22 mm along the length axis in Figure 24, as commented in the Experimental Results section. The fact that the fiber failure location deviated from the experimental results was the main motivation for disabling the fiber failure degradation in Model E. Also, the drastic stiffness change of fiber failure caused an extensive runtime for Model D. Figure 25 shows the peak exposure factor for the 22, 33, 12, 13 and 23 (matrix) components for the outer ply. The maximum was in all cases in the 22 direction. There is a region to the right of the hole at the split with a high exposure factor, corresponding with the high experimental hoop strain region in Figure 21 - Figure 24. As can be seen, models A-E predict an increasing amount of matrix damage. It can be seen that the degree of strain fluctuations in the model's hoop strain graphs corresponds to the degree of matrix damage. The FEA calculations for the onset of fiber failure (first recorded fiber exposure factor above 1.0 in any integration point) are shown in Table 6, from equation 6.0. Initial fiber failure was predicted to happen in the inner hoop layer in all cases, same as in the experiment. For Model A and B it was predicted at the equator of the hole, while for models C-E it occurred at the end of the splits. The table and cycle numbers will later be discussed and suggested as design criteria.  Figure 26 shows the FEA calculated exposure factor in the fiber direction for model A, D and E at 80 000 cycles for the bottom hoop layer. While Model A have relatively little damage in the fiber direction, Model D can be seen to have extensive damage, nearly half the load bearing cross section of the layer has an exposure factor above 1.0. Model E falls in between the two others. Model B and C are not shown, these were also in between A and D in damage extent. Catastrophic failure of the ring happened at 127 814 cycles. At that number of cycles the FEA predicted strains were lower than the experimental strains, as described above, so the accuracy of the FEA model was not too good anymore. This is evidently due to the fact that catastrophic failure of the bottom hoop layer in the model is very difficult to model correctly. Any exposure factor evaluation in top hoop layer is therefor difficult to evaluate and will not be presented. While the strain plots give a good overview of the local behavior of the ring specimen, the displacement curve gives a good indication of how changes to the local stiffness affect the global behavior. Figure 27 shows Model A-E compared to the displacement from the experiment. All models predict a rapid increase in displacement within the first 10 000 cycles. This increase is due to developing matrix damage, making the sample more compliant. The experimental curve show that this damage develops much faster, mostly within the first cycle. The discrepancy is due to the fact that the cycle 1 step was run without any damage in the FEA analysis. This is different from the experiment, where the cycle 1 loading gave initial matrix damage and evidently greater displacement than the models. The models converge with the experiment after the second step at 10 000 cycles. Putting in more cycle jumps would reduce the discrepancy, but since the focus of this study is not on the short term behavior, no further investigations of this phase was done.
Further matrix cracking created a gradual increase of displacement for all models. The increase was lowest for model A and highest for model D, as would be expected from the material properties used in the models given in Table 4. The experimental data show a mainly flat curve up to 80000 cycles and then a gradual increase in displacement. However, within an error of about 10% models A, B and C match the experimental data, models C and E being the best. Model D gives a much too compliant behavior.
The displacement curves show that the models did not properly account for spreading of fiber failure. Displacement increases due to matrix cracking were however relatively well described as evident between 1 to 10 000 cycles. The experimentally observed jumps in the displacement curve due to fiber failure are difficult to capture with such few cycle jumps. Figure 27 The displacement differential from the experiment and test compared.
The macro shear crack length is shown in Figure 28. The crack length is defined as the length of the crack where the cohesive parameters 0 / / / in Figure 8 were exceeded. The shear crack length as defined by the gaps in the DIC contour plot (see Figure 17 and Figure 20) is shorter than this, as voids first come when there is a visual shear crack. It can however be seen that shear crack modeled by FEA extended outside the frame of the DIC (15 mm) for all the models at almost all cycles. Delaminations were not measured experimentally except for some visual investigation of the edges of the ring specimen during testing, see Figure 29. The cracks highlighted with red arrows in Figure 29 are evidently macro cracks, but were not included in the model as they have much less effect on the strain distribution in the loadbearing layers compared to cracks in the hoop layers, while also being harder to predict. Figure 30 shows the delaminations in the models between layer 3 and 4, as they were the most extensive. As can be seen, the delaminations differ from the experiment when comparing Figure 29 and Figure 30. The delaminations in the experiment run through the equator between layer 3 and 4 and are more extensive than in the model. The discrepancy is likely due to that the macro cracks created free edges inside the specimen which delaminations could grow from due to more shear stress in the layer interfaces. The cracks on their own reduced the bending stiffness while also causing delaminations to grow and are the likely reason for the apparent low bending stiffness of the material in the split. The low bending stiffness gives the low strain in the center of the disk in the experiment compared to the model, as shown in Figure 20.

Discussion
The experimental strain field measurements performed by high frequency DIC have shown that the local variations in material properties significantly change the strain. Figure 17 and Figure 18 show the differing strains in the highly strained and critical regions. Using simple symmetry arguments, the strains in the four quadrants around the hole of the ring test should be the same for a component with identical material properties.
The local changes in material properties are a result of the production process and natural variations in material properties. Especially for the filament wound material investigated here local variations in fiber content, fiber placement and presence of voids are quite pronounced, amplifying these effects compared to better-controlled materials, such as prepregs. However, to some extent these variations are present in all composite materials.
Modeling the random variations of the material´s behavior is a challenge. This work used a simplified approach by modeling the ring specimens several (five) times with a constant set of material properties for each modeling run. The initial matrix properties and their degradation characteristics were changed for the different runs. The high frequency DIC measurements allowed comparing experimental fatigue strains with FEA simulations on a high level of detail.
Looking at up to 80000 cycles the comparison showed that the nominal material properties (model A) described the least damaged parts of the specimen well. Model D with degraded properties described the most damaged part of the specimen well. Model B, C and E were between the two. Modeling the entire specimen with a constant set of matrix properties is not ideal, as it deviates from the real physical conditions, but the agreement between experiments and FEA shows that this simplified approach manages to capture the strain envelope reasonably well. The approach should be sufficient for most practical purposes. Model D was the only model that had enough fiber damage to show on the strain curves, with strains in the fiber direction up to 0.8 %. However, the fiber failure occurred at a different spot from the experiment. Seeing as the modeling method did not capture fiber failure correctly, Model E was run without any degradation of fiber failure associated properties and with degradation of matrix properties in between Model D and Model C.
The FEA calculates that the first local fiber failure happens already at 10 000 -20 000 cycles for all models. This is the first point in the models with an integration point showing an exposure factor above 1.0 for the fiber direction. Due to the initial matrix cracking phase being more or less the same for all models, the first registered exposure factor above 1.0 falls within a short cycle window. As stated above it was however concluded that the method falls short of modelling catastrophic fiber failure correctly. The reason being that the experiment consist of many strong and weak fibers, while the model treats all fibers the same. Such a drastic event as fiber failure therefor becomes difficult to model correctly due to statistical variations in the experiment. To model it better, some failure criteria that initiates fiber failure when a region of a certain size have an exposure factor above 1.0 may be better, this is however a computationally expensive approach. The modelling method in Model E may therefor be the most fit for purpose as it gives the user an idea of how big a region may give fiber failure without the added computational cost of modeling the failure (relative to the other models, Model E was computationally faster with less iterations). Knowing when the exposure factor in the fiber direction reaches 1.0 in the model may be a good design input, as seen it occurs about a decade before catastrophic failure, which gives a reasonable safety factor for design purposes. Most importantly, the modelling gives a good indicator of how the strain fields will develop throughout the fatigue life. As can be seen, the general trend is the same for all models, with a flattening of the strain field with increasing cycles and damage.
After initial local fiber failure further local fiber failures develop according to the FEA. As shown in Figure 26, the region with local fiber damage spreads mainly in the loading/hoop direction along the splits and across the width of the specimen. The first global response from fiber damage between 80 000 and 90 000 cycles as seen in Figure 16 is due to an accumulation of local fiber failure that can connect via matrix damage to create a more global crack. The current FEA model is not capable of describing the accumulation of fiber damage properly, as element deletion or contact breaching has to be used in addition to stiffness reduction of the elements to properly characterize the failed regions. Such routines are computationally expensive. The FEA model should be accurate up to first fiber failure. Afterwards, the model shows reasonably well what is happening in the ring specimen, but it should be seen more as a qualitative characterization. Model D and E with the weak matrix describe much more fiber damage than models A, B and C with a stronger matrix, as would be expected.
Catastrophic failure happens at 127814 cycles; a decade after the first fiber failure was predicted. The progressive development of damage leads to fluctuations in the strain field across the width, both in FEA predictions and in experimental DIC measurements. It seems that these fluctuations indicate the onset of serious fiber damage, damage that leads to a global response of the structure. The first fluctuations were already seen for the C, D and E models at 80000 cycles, see Figure 23. The fluctuations are very pronounced at 127000 cycles, see Figure 24, even though the absolute strain values between experiment and FEA do not agree too well. These fluctuations could potentially be used as a Non Destructive Evaluation NDE method indicating imminent catastrophic failure. Qualitatively it can be observed that local fiber failures develop and spread without causing a recognizable global response. Matrix damage increases in parallel. Once a combination of local fiber damage and sufficient matrix damage exists the benign local fiber failures can rapidly combine into global fiber damage causing macroscopic/catastrophic failure.
The FEA used here addresses all failure mechanisms and degradation of material properties after failure, creating a full progressive fatigue analysis. It is based on very simple maximum strain failure criteria and easily obtained material data. Nevertheless, the set of input data needed is large, as shown in Table 2. The good agreement with experimental results up to first fiber failure is an indication that the modeling approach was successful. It is worth looking at some of the simplifications made. All micro failure mechanisms, axial and shear matrix cracking and local fiber failure, are described by simple non-interacting maximum strain criteria. The scatter in experimental data, especially the large effects of locally varying material properties dominates the result making the simple failure criteria adequate. Whether the simple criteria would also work under more complex multiaxial loading conditions is currently unclear and would need further experimental work to find the answer. In principle the cycle jump approach described here can be easily extended to more complex failure criteria if needed.
Another simplification was to prescribe in advance that the dominant shear crack would develop from the equator of the hole in the ring specimen parallel to the load direction in the hoop layers. This simplification reduced the computational effort significantly. For simple and well-defined loading conditions, the position of the shear cracks can be easily estimated in advance and possibly confirmed by simple experiments. The experimentally observed axial cracks were not modeled in advance and subsequently ignored by the FEA. It was argued here that these cracks were not critical for the ring specimens tested. In principle such cracks could be easily added. If the loading directions are completely unknown, the prescribed crack direction can be a severe limitation and the macro shear crack will matter more, such as investigated by Travesa [7].
The planes in which delaminations would develop were prescribed in the same way as for the shear cracks. This simplification should work well under most loading conditions. However, the method for defining the macro matrix crack and delamination is somewhat in contrast to other studies covering this topic, taking a highly simplified approach. Figure 28 shows the crack length defined through breaching of 0 in Figure 8. Compared to the gaps in the strain field in Figure 20 it may appear that the actual shear crack is shorter than the modeled, however, the gaps may appear long after the material has cracked and certainly long after the elastic regime defined by 0 is breached. The gaps simply indicate when the DIC is not able to pick up any displacement in the speckle pattern. The key role of the shear crack from a strain distribution perspective is to hinder transfer of shear strain across the crack line leading to a greater strain at the 10 mm position in the strain graphs in Figure 21 - Figure 24. The strain at the 10 mm position in the curves in Figure 21 - Figure 24 varies very little even though the C and D models have a shear crack twice as long as the A and B models. This indicates that for the split disk ring, the shear crack length is not critical for getting the right strain field in the peak strain regions as long as the crack is longer than a certain minimum.
Strain fields from FEA were compared with DIC strain fields from split disk testing of a composite ring specimen. Extensive variations in damage development and strain fields were measured by DIC over the specimen due to variations in void content, layer thickness and fiber volume fraction. To address the variations, the FEA model was run with a parameter study on matrix properties. The most damaged regions were best modeled by using S-N curves for matrix properties degraded by 40% compared to the original values. The original values described damage in the regions with less defects well. The experimental strain fields fell at or between the modelled material cases, showcasing how much variation can be in a typical filament wound material.
Initial fiber failure could be characterized by an S-N curve measured locally (about 0.5 mm range) by DIC. Despite the curve´s high strength values, fiber failure was predicted conservatively within a decade of the experimental failure, much better than using traditional S-N curves obtained from typical coupon specimens. The model did, however, fall short of being able to correctly describe catastrophic fiber failure (accumulation of local fiber failures), due to its relatively simple nature. The experimental results showed that regions developing fluctuations in the strain fields were the areas where catastrophic fiber failure was initiated. The weak matrix model showed the same fluctuations in the FEA. These fluctuations can be measured by DIC and can be used as a warning for eminent failure. The results indicate that fiber failure and matrix failure are linked.
It can be concluded that the developed model is sufficient to model the complex strain and damage state in a split disk test if run with weak and strong matrix properties. The model is able to show how the strain fields develop and what shape the fields will attain in regions suspect to catastrophic fiber failure.