1. 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 mechanical 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 are 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
Mechanical fatigue in composites causes complex progressive damage development that sets it apart from more conventional materials such as steel. Progressive failure or damage is defined as damage in the material that occurs over a defined time span. In a tensile test the time span is the loading time and progressive failure in the material typically occurs towards the end of that time span just before the specimen fails. In a fatigue test of a metal, the progressive damage will typically occur towards the very last few cycles as a crack is initiated and propagates. In a composite, however, progressive fatigue damage looks rather different from a metal. Instead of damage and crack propagation occurring over a very short cycle span towards end of life, mechanical 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, delamination, and fiber-matrix 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.
Traditionally, fatigue of composites has been divided into two segments, high cycle fatigue (HCF) and low cycle fatigue (LCF). The domains of the two are defined by the failure mechanism which dominates in the final rupture of the material in question. HCF is dominated by matrix damage, while LCF by fiber failure [
8,
11,
12,
13]. Notably, a relatively large volume of the material is characterized using this traditional approach. The volume is the typical size of the gauge section of a test specimen of roughly 1000 mm
3. In this study, as will be explained later, fatigue is described locally around a fiber bundle, addressing a volume on the scale of a typical element in a FEA, which can be 0.2 mm
3 or less. When developing a finite element model to describe such local fatigue damage in a component, some highly stressed material (near a defect or geometric stress concentrator) may fail after few cycles as “low cycle fatigue”, while material in the lower stressed regions may be in the “high cycle fatigue” domain. As such, a component that catastrophically fails after many cycles fails globally in the HCF domain, but it may have local material that also fails in the LCF domain. A FEA material model addressing local fatigue failure has to take local low and high cycle fatigue into account.
High and low cycle fatigue is, however, defined for the composite material, while the FEA model in this study considers fibers and matrix by themselves, though with some interaction effects. Since the matrix is much weaker and traditionally degrades faster due to fatigue than the fiber (higher slope of the S-N curve [
14]), local matrix degradation will naturally dominate in the global high cycle fatigue range. Local fiber failure will be prevalent in the low cycle fatigue range. This study focuses on modeling a global high cycle fatigue experiment.
When using finite element models to model matrix cracking, both initiation and propagation needs 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., Turon [
7] and Nixon et al. [
15] based on the method suggested by Harper et al. [
16]. Attempts to simplify 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 [
17], Flatscher et al. [
18], and Gagani et al. [
19]. 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. For the static case, Rozylo [
20] satisfactorily managed to model crack propagation without predefined crack directions using the cohesive zone modelling approach in combination with a user element subroutine (UEL) with promising results. Still, all of the above-mentioned models are relatively academic and not easy for the average design engineer to implement or to get the correct input data for. The models have been developed with lab experiments in mind and not real designs. In this study, the model was developed with the design engineer in mind and then tested on a complex lab experiment. This study therefore has a somewhat different format than most academic publications covering the topic, having a wider scope and less in detail investigations of the experimental results and modelling. This study would however not have been possible without the past academic literature going in depth in DIC and Abaqus in particular.
There are currently some commercial composite mechanical fatigue numerical frameworks available, most notably FEMFAT [
21] and Fe-safe [
22]. 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 mechanical fatigue, most notably by Koch et al. [
23]. 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. [
24]. 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 [
25,
26,
27], matrix voids affect the mechanical fatigue properties to a great degree, which was also found in the presented work here. This effect is also present in other similar materials such as concrete, where nanoparticles can be added to fill the voids and reduce microcracking [
28]. 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 are used as input for local properties in most mechanical 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 fall 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. The DIC methods used in this study have been elaborated in two articles explaining how to trace progressive failure in composites [
29] and how to estimate a local S-N curve using DIC data [
30]. DIC has recently been proven to be a valuable tool in estimating material parameters such as the fracture toughness [
31] by tracing crack propagation visually. The very direct observation method (visual) and the vast amount of data make DIC a measurement method with huge potential for more exactly estimating and monitoring material parameters. That is, provided that the user has the ability to take advantage of the data using modern data tools.
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:
- i.
Micro matrix failures were modeled using a continuum damage approach as changes in the stiffness matrix without directionality of the cracks.
- ii.
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.
- iii.
Discontinuities in the stiffness due to crack growth were modeled using an on and off loading approach in combination with simplified cycle jumping [
15,
16].
- iv.
A range of polymer matrix properties were modeled to investigate the natural variations of material properties.
7. 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 80,000 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. However, as stated above, it was concluded that the method falls short of modelling catastrophic fiber failure correctly. The reason being that the experiment consisted of many strong and weak fibers, while the model treats all fibers the same. Such a drastic event as fiber failure therefore becomes difficult to model correctly due to statistical variations in the experiment. To model it better, some failure criteria that initiate fiber failure when a region of a certain size has an exposure factor above 1.0 may be better. However, this is a computationally expensive approach. The modelling method in Model E may therefore 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 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. Considering the material model’s poor ability to correctly model the progression of fiber failure, but good capability to model matrix damage, it can be said that it is more fit for modelling components with a high cycle fatigue issue rather than low cycle.
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 127,814 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 80,000 cycles, see
Figure 23. The fluctuations are very pronounced at 127,000 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 mechanical 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, as investigated by Turon [
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
in
Figure 7. 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
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 22,
Figure 23 and
Figure 24. The strain at the 10 mm position in the curves in
Figure 21,
Figure 22,
Figure 23 and
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. This is an attribute of the modeled geometry. In other applications the macro shear crack will matter more, as investigated by Turon [
7]. For future work, it could be possible to combine the methods from the UMAT in this work into a UEL with CZM. As demonstrated by Rozylo [
20], the CZM approach is robust and able to predict crack paths in complex strain fields. With the cycle jump approach in the UMAT presented here it would be possible to make a UEL that takes into account both macro and micro cracking, making for a very robust progressive fatigue model. The runtime would very likely be an issue; however, simpler models than the one in this study would be sufficient to serve as a proof of concept and with the computers of tomorrow runtime may not be a worry in the future.
Another simplification was the use of the cycle jump method. As found in previous studies, the method seems to work well. The analysis steps were chosen here based on experimental results. It was a convenient way to identify the critical steps, but only possible if experimental data are available. To increase applicability and improve accuracy, a better estimation of cycle jumps would be needed. This could be done by convergence of strain or damage over time by running several models with different cycle jumps. These models could be more coarsely meshed or only contain parts of the full model. More elaborate approaches to estimate cycle jump size and position already exist for simpler models [
7,
15,
23] and should be possible to implement on larger models as well.
It seems a good balance of simplicity of the FEA and accuracy of the results was found for this study. Adding more complexity to the FEA, especially fiber matrix interaction effects and cycle jump iteration schemes, may get closer to the physics of the behavior of the composite, but also increases dramatically the computational effort. Whether a more complex model will improve the results remains to be seen, because usually further uncertain assumptions need to be added on a detailed level. The variation in material properties remains in all cases and requires using a worst-case approach for design calculations at the end.
Traditional design calculations would use S-N curves for fiber dominated failure obtained from coupon tests and apply them to the stress/strain concentration points in the structure. For this material an S-N curve with the origin at 22,150 microstrain [
36] would yield very conservative first fiber failure estimates. In this study the highest strain near the hole was around 1.5 to 1.8%. This would lead to 100 to 1000 cycles to failure, far below the actual catastrophic failure. The reason for the mismatch is using an S-N curve describing large catastrophic failure on the scale of a few cm to local non-critical damage development at a stress concentration.
The FEA shown here uses an S-N curve describing local fiber failure taking local material variations of the matrix into account. This curve is more realistic for use in stress/strain fields with large gradients. Since the accuracy of the FEA currently drops after first local fiber failure, a designer could use the first fiber failure for the weak material D as a design criterion. The predicted cycles to failure are then 10,000, which are much more than predicted by the traditional method. It should still be a reliable approach, because the FEA is fully capable of reproducing local strain fields. This approach is basically designing for first ply failure, a method widely used for static analysis. The approach has been applied before in the FADAS mechanical fatigue model for simple test coupons by Passipoularidis et al. [
48]. The FEA can be used further to estimate the location of final catastrophic failure and it could indicate how the strain fluctuations would look that would develop as a warning before catastrophic failure.
8. Conclusions
A finite element analysis (FEA) material model describing mechanical fatigue of composites by a progressive degradation model using the cycle jump method was developed. It is aimed at pressure vessel applications and focuses on fatigue behavior under tensile load cases, such as that found in industrial pressure vessels. It uses typical failure criteria for transversely orthotropic materials and degrades the stiffness in each material direction based on the Miner sum damage calculation with log-log S-N curves. The model is relatively simple and requires only the typical input parameters used for composite laminate analysis.
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. The strain in the four regions around the split disk hole varied with a factor of 1.5 to 2 at the most over the course of cycling. To address the variations, the FEA model was run with a parameter study on matrix properties. The most damaged regions (with strain 1.5 to 2 times higher than the least damaged) 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 there can be in a typical filament wound material. Considering the 1.5 to 2 factor difference, a reduction of 40% in strength for the matrix can be considered reasonable.
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 and the sudden nature of catastrophic fiber failure. 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.