Predicting Fatigue Damage of Composites Using Strength Degradation and Cumulative Damage Model

This paper discusses the prediction of fatigue response of composites using an empirical strength and stiffness degradation scheme coupled to a cumulative damage accumulation approach. The cumulative damage accumulation approach is needed to account for the non-constant stress levels that arise due to stress distributions from stiffness degradation during the fatigue loading. Degradation of strength and stiffness during fatigue loading of the composite was implemented by following the empirical model presented by Shokrieh and Lessard with some modification and correction to the non-dimensional load parameter definition. The fatigue analysis was performed using ABAQUSTM finite element software using a user-defined material subroutine UMAT developed for the material response. Implementation results were first verified for unidirectional laminate test cases and validated by predicting stress versus life (S-N) curves for several laminate coupons test simulations and residual strengths of Open Hole Tension (OHT) specimens subjected to constant amplitude fatigue loading.


Introduction
The Composite materials used in primary aircraft structures can experience fatigue damage and failure due to the repeated loads they experience.Tools and models that stimulate the response of composites under cyclic loads are necessary in order to design structures that can withstand such damage.Effective tools can also help to estimate remaining lifetimes and residual strength and stiffness of composites as they are subjected to repeated or cyclic loading.
The response of composite structures under fatigue loading investigated by a number of researchers that has led to development of a number of fatigue prediction models.These models can be classified in the following categories [1]: fatigue life concepts [2,3], phenomenological models with stiffness [4,5], strength degradation models [6,7], continuum damage mechanics (CDM) based models [8][9][10] and micromechanics models [11,12] and uncertainty and Bayesian based probabilistic models [13][14][15].
A commonly used approach in fatigue life predictions is to use stress versus life, known as S-N curves.These fatigue life models [2,3] use test data from constant amplitude tests at varying stress levels and fit empirical S-N curves that relate the number of cycles to failure to the applied stress level.The constant amplitude cyclic loads are characterized by the mean stress level (σ m ) and the amplitude (σ a ) of the stress variations about the mean.This is alternatively expressed in terms of the maximum stress (σ max = σ m + σ a ) and the stress ratio or R-ratio (R = σ min /σ max ).Goodman type diagrams [3] provide models that can fit the cycles to failure as a function of mean stress and varying stress amplitude.This requires extensive experimental work, as measuring the cycles to failure for all possible mean and alternating stress values states that do not lead to static stress failure is quite large.The laminates chosen for such tests should be able to capture the various modes of intra-laminar (fiber and matrix) failure within a ply.A limitation of this approach is that at low stress levels, it is expensive to perform the cyclic tests until the specimen fails.As a result, the tests are truncated at number of cycles that represents an upper limit of what an engineering component would experience in service.
Fatigue response of composite materials can be described by the stiffness and strength degradations they exhibit due to accumulation of damage.These are characterized by the residual stiffness and strength measurements of components tested at different stress levels and to varying fractions of their life (previously determined for life models).Phenomenological models based on empirical fits to the test data provide evolution laws that describe a gradual reduction in the laminate stiffness and strength at macroscopic scale due to the damage evolution at microscopic constituent length scales.A number of models have been proposed for stiffness degradation during fatigue [4,5] and strength degradation [6,7].In strength degradation models, the ply fails when the strength value reaches the ply stress level.
Fatigue behavior in composites has also been modelled using continuum damage mechanics-based approaches by using fatigue cycles to account for damage growth [8][9][10].These types of approaches suggest that degradation of the composite material is in direct relation with specific damage such as matrix cracking or delamination.The continuum damage mechanics (CDM) models correlate one or more material properties with defined damage variables which are determined quantitatively from the damage progression or failure mechanism.
One of the new category of approaches to model fatigue in composites are models based on uncertainty and Bayesian based probabilistic framework [13][14][15].Most of the probabilistic models are based on quantifying model uncertainties for different fatigue models.A progressive damage in composite material using a full Bayesian prediction and updating framework is presented by Chiachío et al. [13].They used a cumulative damage model based on the theory of Markov chains to infer a complete damage progression.Chiachío et al. [14] presented one such model based on the uncertainty for empirical fatigue model.The authors did not account for any interlaminar failure for the Bayesian approach.In the work of Shankaraman et al. [15], several different models such as finite element (FE), surrogate model analysis, crack growth model and the effect of the types of uncertainties such as physical variability, data uncertainty and errors were incorporated in both the inverse problem and the forward problem.Diao et al. [16] incorporated a statistical approach for reinforced composite laminates subjected to multiaxial loading for the static and fatigue properties.They obtained the statistical distribution function of the material strength parameters by fitting experimental data to a two-parameter Weibull distribution.
After a review of material models proposed in the literature, we selected an empirical model proposed by Shokrieh and Lessard [6,[17][18][19].Our model tests six basic unidirectional laminates under fatigue load to produce fatigue data for AS4 3501-6 graphite epoxy material.A material user subroutine (UMAT) for modeling the material degradation during fatigue loading was developed for use within the ABAQUS™ finite element software.A cumulative damage approach was used to account for the change in stress amplitude resulting from the stress state redistribution after failure.The material degradation model is based on life prediction models that use empirically fitted constant life curves [20] to determine the number of cycles to failure and a degradation model where the reductions in stiffness and strength occur as a function of the number of cycles of loading experienced as a fraction of the total life time.

Constant Life Diagrams
Constant life diagrams (CLD) are obtained from empirical fits to fatigue test data that relate cycles to failure to the mean and alternating stress values of the constant amplitude fatigue loading.The constant life curves are simply the contour curves corresponding to a given number of cycles of failure, plotted in the space of mean and alternating stress.Many different models and methods have been proposed to parametrically describe these curves [3]. Figure 1a shows a typical constant life diagram (CLD).The main axes on a CLD graph [20] are the mean cyclic stress (σ m ) and the cyclic stress amplitude (σ a ).In addition, radial lines of constant R-ratio are shown on the graph to easily recognize how the stress ratio affects the shape of these curves.Goodman diagrams use linear fits for constant life curves.However, for composite materials researchers find that constant life curves are nonlinear and non-symmetric with respect to the mean stress axis.These curves allow interpolation of test data to estimate lifetimes for given stress loadings.A number of analytical models for parametrizing the curves have been published in literature [3,[20][21][22].Adam et al. [20] presented an analytical model that not only parametrizes the CLD but also allows unifying the different curves into a single two-parameter fatigue curve which in this work is referred to as the normalized fatigue life model (Figure 1b).failure, plotted in the space of mean and alternating stress.Many different models and methods have been proposed to parametrically describe these curves [3]. Figure 1a shows a typical constant life diagram (CLD).The main axes on a CLD graph [20] are the mean cyclic stress ( ) and the cyclic stress amplitude ( ).In addition, radial lines of constant R-ratio are shown on the graph to easily recognize how the stress ratio affects the shape of these curves.Goodman diagrams use linear fits for constant life curves.However, for composite materials researchers find that constant life curves are nonlinear and non-symmetric with respect to the mean stress axis.These curves allow interpolation of test data to estimate lifetimes for given stress loadings.A number of analytical models for parametrizing the curves have been published in literature [3,[20][21][22].Adam et al. [20] presented an analytical model that not only parametrizes the CLD but also allows unifying the different curves into a single two-parameter fatigue curve which in this work is referred to as the normalized fatigue life model (Figure 1b).A two-parameter model of a constant life diagram was fitted to the non-dimensional mean stress q and stress amplitude a values and the ratio c of the static strengths Sref1 and Sref2 of the orthotropic ply material as shown in Equation (1).
For the composite lamina in the fiber directions, the Sref1 and Sref2 terms in Equation ( 1) take on values of the fiber tension and fiber compression strengths XT and XC, respectively.For the transverse to fiber direction, the Sref1 and Sref2 terms in Equation ( 1) take on values of the transverse to fiber compression and transverse to fiber tension strengths YC and YT, respectively.Normalization allows the definition of a unified load parameter u to fit constant life diagrams.The parametric model proposed by Adam et al. [16] is given by: where f, u and v are model fit constants.The exponents u and v determine the shapes of the left and right wings of a bell-shaped CLD curve.However, Gathercole et al. [22] showed that the degree of curve-shape asymmetry was not very great for carbon-epoxy composites; therefore, the assumption that u and v are equal is valid.Thus, u is given by: In Equation (3), f is a curve fitting parameter with a value of 1.06.The unified load parameter u is used to relate the loading to the number of cycles to failure in the empirical model.Shokrieh used a linear fit model to the logarithm of the number for cycles to failure log(Nf) as a function of the unified load parameter u, A two-parameter model of a constant life diagram was fitted to the non-dimensional mean stress q and stress amplitude a values and the ratio c of the static strengths S ref1 and S ref2 of the orthotropic ply material as shown in Equation (1).
For the composite lamina in the fiber directions, the S ref1 and S ref2 terms in Equation ( 1) take on values of the fiber tension and fiber compression strengths X T and X C , respectively.For the transverse to fiber direction, the S ref1 and S ref2 terms in Equation (1) take on values of the transverse to fiber compression and transverse to fiber tension strengths Y C and Y T , respectively.Normalization allows the definition of a unified load parameter u to fit constant life diagrams.The parametric model proposed by Adam et al. [16] is given by: where f, u and v are model fit constants.The exponents u and v determine the shapes of the left and right wings of a bell-shaped CLD curve.However, Gathercole et al. [22] showed that the degree of curve-shape asymmetry was not very great for carbon-epoxy composites; therefore, the assumption that u and v are equal is valid.Thus, u is given by: In Equation (3), f is a curve fitting parameter with a value of 1.06.The unified load parameter u is used to relate the loading to the number of cycles to failure in the empirical model.Shokrieh used a linear fit model to the logarithm of the number for cycles to failure log(N f ) as a function of the unified load parameter u, where A and B are the curve fitting constants.
The number of cycles to failure (N f ) is computed for the state of axial fiber direction loading using stress σ 11 and for the transverse to fiber direction loading using stress σ 22 by choosing the appropriate normalization stress in the definition of the load parameter u, as described above.Shokrieh and Lessard also showed that the unified load parameter for the pure shear state of stress (σ 12 , σ 13 , or σ 23 ) requires modification.There is no equivalent loading direction strength tension or compression when loaded in shear.The unified load parameter for shear loading as proposed by Shokrieh [6] is In the above expression z = σ a /σ shear and q = σ shear /S, where σ a is the stress amplitude, σ shear is the applied shear stress and S is the static shear strength.
The unified parameter u, in Equation ( 3) is also modified slightly for loading in the transverse to fiber direction (matrix failure dominated state of stress), σ 22 .Our implementation of the model found that the estimation of the number of cycles to the failure using the general form of Equation (3) used by Shokrieh [6] does not provides a realistic prediction for the loading transverse to the fiber.The model is required to be normalized by the compression strength, where the general form of Equation (3) normalized using the tensile strength in the matrix direction.Implementations of Shokrieh's model by other authors found in the literature does not address this normalization approach for the unified load parameter in the matrix direction [23][24][25][26].The unified parameter in the state of stress σ 22 is normalized by the strength parameter Y C and the unified parameter, u can be expressed as: The experimental value of the transverse compressive strength Y C is significantly higher than the tensile strength Y T , where the longitudinal tensile strength X T is generally larger than the compressive strength X C .Therefore, for the σ 22 state of stress, the normalization was performed by the compressive strength of Y C , as shown in Figure 2b.
where A and B are the curve fitting constants.
The number of cycles to failure ( ) is computed for the state of axial fiber direction loading using stress σ11 and for the transverse to fiber direction loading using stress σ22 by choosing the appropriate normalization stress in the definition of the load parameter u, as described above.Shokrieh and Lessard also showed that the unified load parameter for the pure shear state of stress (σ12, σ13, or σ23) requires modification.There is no equivalent loading direction strength tension or compression when loaded in shear.The unified load parameter for shear loading as proposed by Shokrieh [6] is In the above expression z = σa/σshear and q = σshear/S, where σa is the stress amplitude, σshear is the applied shear stress and S is the static shear strength.
The unified parameter u, in Equation ( 3) is also modified slightly for loading in the transverse to fiber direction (matrix failure dominated state of stress), σ22.Our implementation of the model found that the estimation of the number of cycles to the failure using the general form of Equation ( 3) used by Shokrieh [6] does not provides a realistic prediction for the loading transverse to the fiber.The model is required to be normalized by the compression strength, where the general form of Equation (3) normalized using the tensile strength in the matrix direction.Implementations of Shokrieh's model by other authors found in the literature does not address this normalization approach for the unified load parameter in the matrix direction [23][24][25][26].The unified parameter in the state of stress σ22 is normalized by the strength parameter YC and the unified parameter, can be expressed as: The experimental value of the transverse compressive strength YC is significantly higher than the tensile strength YT, where the longitudinal tensile strength XT is generally larger than the compressive strength XC.Therefore, for the σ22 state of stress, the normalization was performed by the compressive strength of YC, as shown in Figure 2b.Table 1 shows a comparison of the computed values of the unified load parameters, u and using the model presented by Shokrieh [6] given Equation ( 5) and the modification that was made in this paper which is given by Equation ( 6).The accurate prediction of the unified parameter, u modification is significant to predict fatigue life along the matrix direction as failure generally initiates in the matrix.Table 1 shows a comparison of the computed values of the unified load parameters, u and u using the model presented by Shokrieh [6] given Equation ( 5) and the modification that was made in this paper which is given by Equation ( 6).The accurate prediction of the unified parameter, u modification is significant to predict fatigue life along the matrix direction as failure generally initiates in the matrix.

Shokrieh's Material Degradation
Shokrieh and Lessard [17,18] proposed a power law model (Figure 3) to represent the relationship between the number of cycles and the residual strength or stiffness of a laminate.This model is capable of capturing both the "sudden death" [27] and "wear out" [28] modes for failure in laminates.The "sudden death" mode is observed in fatigue tests at high stress levels whereas the "wear out" mode is observed in fatigue tests at low stress levels.

Shokrieh's Material Degradation
Shokrieh and Lessard [17,18] proposed a power law model (Figure 3) to represent the relationship between the number of cycles and the residual strength or stiffness of a laminate.This model is capable of capturing both the "sudden death" [27] and "wear out" [28] modes for failure in laminates.The "sudden death" mode is observed in fatigue tests at high stress levels whereas the "wear out" mode is observed in fatigue tests at low stress levels.The power law fit for the degradation of the strength and the stiffness follows the form: The number of cycles to failure is normalized as follows: The residual strength and stiffness are normalized using: Thus, the residual strength and stiffness as a function of the number of cycles, stress level and stress ratio are given by: where σr is the residual strength, σ0 is static strength, σ is the stress level, n is the number of applied cycles, Nf is the number of cycles to failure at the current stress level (fatigue life), α and β are the experimental curve fitting parameters, E is the residual stiffness, σ is the stress level, n is the number of applied cycles, ES is the initial stiffness, εf is the strain at failure at applied stress σ and γ and λ are The power law fit for the degradation of the strength and the stiffness follows the form: The number of cycles to failure is normalized as follows: The residual strength and stiffness are normalized using: Thus, the residual strength and stiffness as a function of the number of cycles, stress level and stress ratio are given by: where σ r is the residual strength, σ 0 is static strength, σ is the stress level, n is the number of applied cycles, N f is the number of cycles to failure at the current stress level (fatigue life), α and β are the experimental curve fitting parameters, E is the residual stiffness, σ is the stress level, n is the number of applied cycles, E S is the initial stiffness, ε f is the strain at failure at applied stress σ and γ and λ are the experimental curve fitting parameters [17,18].The curve fitting parameters α, β, γ and λ are given in Table 2.These fitting parameters provided by Shokrieh and Lessard [6,[17][18][19] contain typographical errors where the α and β parameters are switched.The material constants determined for the residual strength and the stiffness model reported by Shokrieh and Lessard contains typographical error in all the reported publications by the authors.Figure 4 shows the use of the parameters α and β for the strength degradation model of the uniaxial 0 • test specimen as given by Equation (10).It is evident that the values of the parameters are switched [17,18] and causes significant discrepancies in residual strength prediction.the experimental curve fitting parameters [17,18].The curve fitting parameters α, β, γ and λ are given in Table 2.These fitting parameters provided by Shokrieh and Lessard [6,[17][18][19] contain typographical errors where the α and β parameters are switched.The material constants determined for the residual strength and the stiffness model reported by Shokrieh and Lessard contains typographical error in all the reported publications by the authors.Figure 4 shows the use of the parameters α and β for the strength degradation model of the uniaxial 0° test specimen as given by Equation (10).It is evident that the values of the parameters are switched [17,18] and causes significant discrepancies in residual strength prediction.

Cumulative Damage Approach
Number of authors have implemented Shokrieh and Lessards's empirical model to predict fatigue behavior of composites [23][24][25][26].The authors [23][24][25][26] did not account for any damage accumulation when the state of stress redistributes due to failure in a material point.The above models were developed using constant amplitude tests.To apply them to structural fatigue simulations under varying loads, a cumulative damage approach [29] is needed.Variations in stress amplitudes for applied stress may be due either to the applied load being time varying or to stress distributions arising from stiffness changes caused by stiffness degradation or local material failure.The cumulative damage approach introduces a damage variable to track damage accumulations and

Cumulative Damage Approach
Number of authors have implemented Shokrieh and Lessards's empirical model to predict fatigue behavior of composites [23][24][25][26].The authors [23][24][25][26] did not account for any damage accumulation when the state of stress redistributes due to failure in a material point.The above models were developed using constant amplitude tests.To apply them to structural fatigue simulations under varying loads, a cumulative damage approach [29] is needed.Variations in stress amplitudes for applied stress may be due either to the applied load being time varying or to stress distributions arising from stiffness changes caused by stiffness degradation or local material failure.The cumulative damage approach introduces a damage variable to track damage accumulations and maintains continuity of the damage variable when stress changes occur.In this way damage accumulation at a new stress level starts with the initial damage level being what was sustained already due to previous history at different stress levels.The cumulative damage not only account for the damage due to the stress state redistribution but also accounts for the existing damage.
If a material that degrades due to cyclic loading is subjected to uniaxial fatigue load with maximum stress σ, failure occurs when the residual strength σ r , reaches the applied stress level σ.Using the strength degradation model (Equation ( 10)), the ratio of the stress margin due to the degradation of strength (σ r − σ) after n cycles, to the maximum stress margin at start of the fatigue loading (σ 0 − σ) is expressed as: where σ r , is the residual strength, (n) is the number of cycles at applied load and (N f ) is the number of cycles to failure.The above ratio varies from unity at start to zero at failure.In this work, we adopt a definition of damage based on the reduction in material strength.The damage variable D therefore takes the form: Figure 5 provides a schematic of the cumulative damage approach.The curves represent the residual strength degradation as a function of the number of cycles for different constant amplitude fatigue loading stress levels (σ 1 < σ 2 < σ 3 ).Connecting the end points of each curve also gives the S-N curve.In the cumulative damage approach, we maintain equivalency of damage between fatigue cycles when there is a stress level change by equating the residual strength.Let us assume D 1e is the damage accumulated at the end of the elapsed loading block of n 1 fatigue cycles at stress level σ 1 .If the stress level changes to σ 2 after the first block of n 1 fatigue cycles at stress level σ 1 , the damage accumulation starts of the new block D 2s must be equal to D 1e .To use a degradation model for constant amplitude loading, with the new stress level σ 2 we calculate first the equivalent number of cycles n s 2 of loading at stress level σ 2 that would have created the damage level D 1e .Using the definition of damage in Equation ( 13), the equivalent cycles n s 2 is calculated as shown below in Equation (14).
The remaining life at stress level If this is greater than zero then the fatigue loading can continue and the material degrades following the degradation law corresponding to stress level σ 2 .The damage accumulation in the new block n 2 cycles at stress level σ 2 is calculated based on total cycle count equal to (n s 2 + n 2 ).For any kth new block of fatigue load, we calculate the equivalent elapsed cycles as n s k and the remaining life n rk using Equation (15).
corresponding to stress level σ2.The damage accumulation in the new block n2 cycles at stress level σ2 is calculated based on total cycle count equal to ( + ).For any kth new block of fatigue load, we calculate the equivalent elapsed cycles as and the remaining life using Equation (15).

Failure Criteria and Instant Material Degradation
In our work, we use Hashin's [30] failure criteria to determine the fiber and matrix failures in a unidirectional lamina.Equations ( 16)-( 19) summarize equation for the failure envelopes for fiber and matrix failure from Hashin's criteria.
Fiber failure in tension Fiber failure in compression In Equations ( 16)-( 19), X T and X C are the longitudinal tensile and compressive strengths, Y T and Y C are the transverse tensile and compressive strengths, S 12 is the in-plane shear strength and S 23 is the out of plane shear strength.An instantaneous ply stiffness degradation scheme is used for the progressive failure when ply failure is detected.Camanho [31] showed that a small but finite residual stiffness after failure is needed to accurately predict the evolution of progressive failure.Table 3 shows the residual stiffness of the ply following failure in each mode.

Finite Element Implementation
Finite element analysis software ABAQUS [32] was used to perform the fatigue analysis of composite laminates and laminates structures.Fatigue material degradation was implemented as a user-defined material subroutine (UMAT) for use in ABAQUS non-linear implicit solution.The laminates were modelled as a three-dimensional solid with ply by ply representation.Each ply was modeled using the eight-node bilinear reduced integration C3D8R solid elements, with enhanced hourglass control.Figure 6 shows a flow chart of the material degradation model implementation.The Jacobian update implemented in the UMAT subroutine was consistent for use with geometric nonlinear analyses.Failure indices, residual strengths and damage parameters at each integration point are stored as Solution Dependent Variables (SDVs) in ABAQUS.
The fatigue loading is simulated by two load steps: a static load step under load control up to the maximum loading during each cycle, followed by a constant load over a specified time.The degradation with cycles is performed with a forward time marching scheme in with an increment of 1000 cycles for each time increment in the second step.To determine residual stiffness, additional load steps are added to unload the structure and then load the structure under displacement control until failure.

Verification and Validation
The UMAT subroutine implemented in ABAQUS was first verified using experimental data for the six unidirectional laminate cases presented by other authors for stress ratios of R = 0.1 and 10.This was followed by simulations of laminate coupons tested under homogenous stress states.Figure 7 shows the comparison of the model predicted S-N curve with experimental data and simulations results obtained from the literature for [30°]16 laminate from Ref.

Verification and Validation
The UMAT subroutine implemented in ABAQUS was first verified using experimental data for the six unidirectional laminate cases presented by other authors for stress ratios of R = 0.1 and 10.This was followed by simulations of laminate coupons tested under homogenous stress states.Figure 7 shows the comparison of the model predicted S-N curve with experimental data and simulations results obtained from the literature for [30   ] S from Ref. [23] and [45 • /−45 • ] 2s from Ref. [33].The test data exhibits significant scatter.The empirical model predictions lie within the scatter of the test data, except for the fatigue test at loading level at 80% of ultimate strength.All tests and simulations correspond to AS4/3501-6 graphite epoxy composite from Shokrieh and Lessard [6] and listed in Tables 2 and 4.

Open-Hole Tension Specimen
Han et al. [34] published experimental data for open-hole specimen fatigue tests performed for AS4 3501-6 composite material laminates.They tested seventy open-hole-laminate coupon specimens under constant amplitude tension and compression fatigue loading.The laminate tested was a 24-ply quasi-isotropic [0/ ± 45/90]3S laminate.The tension-tension fatigue tests at R = 0.1 were performed for maximum stress levels of 40%, 50%, 60%, 70% and 80% of the Notched Tensile Strength (NTS) of the laminate.NTS is the ultimate stress at failure of the laminate under quasi-static tension loading.After one million fatigue cycles at 40%, 50%, 60%, 70% and 80% of the NTS, the specimens were loaded to static failure to obtain the Residual Tensile Strength (RTS).
The instantaneous degradation scheme is known to be mesh sensitive.In fatigue only a small region around the hole sustains damage due to the stress concentration.However, to perform progressive failure to get NTS or RTS, where damage leads to two part failure of the specimen, the entire width of the specimen at the hole was finely meshed.The region away from the damage region where the mesh was not refined was modeled as elastic and did not include the UMAT provided damage degradation.This resulted in a mesh with total 13,600 elements in each ply and

Open-Hole Tension Specimen
Han et al. [34] published experimental data for open-hole specimen fatigue tests performed for AS4 3501-6 composite material laminates.They tested seventy open-hole-laminate coupon specimens under constant amplitude tension and compression fatigue loading.The laminate tested was a 24-ply quasi-isotropic [0/ ± 45/90] 3S laminate.The tension-tension fatigue tests at R = 0.1 were performed for maximum stress levels of 40%, 50%, 60%, 70% and 80% of the Notched Tensile Strength (NTS) of the laminate.NTS is the ultimate stress at failure of the laminate under quasi-static tension loading.After one million fatigue cycles at 40%, 50%, 60%, 70% and 80% of the NTS, the specimens were loaded to static failure to obtain the Residual Tensile Strength (RTS).
The instantaneous degradation scheme is known to be mesh sensitive.In fatigue only a small region around the hole sustains damage due to the stress concentration.However, to perform progressive failure to get NTS or RTS, where damage leads to two part failure of the specimen, the entire width of the specimen at the hole was finely meshed.The region away from the damage region where the mesh was not refined was modeled as elastic and did not include the UMAT provided damage degradation.This resulted in a mesh with total 13,600 elements in each ply and 11,600 elements in the refined mesh region (Figure 8).
Figure 9 shows the stress strain response of the laminate under quasi-static loading in pristine state and post fatigue tested state.The inset shows the ungripped portion of the specimen model for FE simulations.Table 5 summarizes the comparisons of the predicted and experimental values of NTS and RTS of the notched laminate specimen.
The experimentally measured NTS of the specimen [34] was 413 MPa.Progressive failure analysis prediction of NTS using the Hashin's failure criteria and instantaneous stiffness degradation describer earlier was 405 MPa.It is to be noted the material properties for AS4 3501-6 used for analyses were from Shokrieh [6] as this was not provided in the reports of Han et al. [34].No material property calibrations were done to match test results.Progressive failure analysis of the pristine specimen shows that the specimen experiences matrix damage in the 90 • and ±45 • plies at 40% of the NTS.This means for fatigue loading above this level there is damage from the very first cycle.This damage is higher when the specimen is loaded at 50 to 80% of the NTS.Fatigue damage for these cases represents growth of the initial damage.
pristine specimen shows that the specimen experiences matrix damage in the 90° and ±45° plies at 40% of the NTS.This means for fatigue loading above this level there is damage from the very first cycle.This damage is higher when the specimen is loaded at 50 to 80% of the NTS.Fatigue damage for these cases represents growth of the initial damage.pristine specimen shows that the specimen experiences matrix damage in the 90° and ±45° plies at 40% of the NTS.This means for fatigue loading above this level there is damage from the very first cycle.This damage is higher when the specimen is loaded at 50 to 80% of the NTS.Fatigue damage for these cases represents growth of the initial damage.The prediction the fatigue residual strength RTS notched specimens subjected to constant amplitude cyclic loading at 40%, 50%, 60% and 70% for 10 6 cycles, showed an increase in the residual strength.The local damage around the hole leads to stress distributions that minimize the stress concentration at the hole.The slope of the axial stress-strain curve (Figure 9) from progressive failure analyses of the specimens in post-fatigued state with simulated fatigue damage shows that the damage leads to a relatively small stiffness loss for cases when the fatigue loading maximum stress levels were 40% to the 70% of NTS.This is due to the fatigue damage being primarily due to matrix failure and damages being localized to the region around the hole.The stiffness loss due to fatigue damage was 12% for the specimen loaded at 80% of the NTS.None of the specimens failed after one million fatigue cycles at stress levels 40-80% of NTS.Post inspection of the tested specimens tested at 80% of the NTS showed significant delamination and edge peeling of the plies [34].The delamination fatigue damage mode is not included in our fatigue analysis model.
During a fatigue loading simulation, as the strength degrades the failure index increases until it reaches the critical value of one, at which point the ply is considered to be fully failed.Figures 10-13 show contour plots of the failure indices in each ply for all four failure modes after one million fatigue cycles at 40-80% of the NTS.Failed elements are indicated by failure index of 1.0.The contour plots of the failure indices show the extent of failure and the margin to failure at the maximum stress level of the fatigue load cycle.Examination of these failure indices indicate that there is significant fatigue damage in matrix tension mode in 45 • and 90 • plies and matrix compression mode in the 0 • plies and some fatigue damage due to fiber tension mode in the 0 • plies at the edge of the hole.
The fatigue model used here is based on stiffness and strength reduction directly applied to the engineering stiffness constants and strengths that are ply properties.To quantify and visualize the level of damage, a measure of the relative reduction in the stiffness/strength parameter due to damage D P , is calculated using Equation (20).
The fatigue model used here is based on stiffness and strength reduction directly applied to the engineering stiffness constants and strengths that are ply properties.To quantify and visualize the level of damage, a measure of the relative reduction in the stiffness/strength parameter due to damage , is calculated using Equation (20).Figures [14][15][16][17][18] show the relative reduction of stiffness and strength properties due to fatigue damage that exceeded 10%.The significant reductions (80-100%) are observed in the transverse to fiber (matrix) direction moduli E2 and tension strength YT and the plane shear moduli G12 in the 90° ply along the symmetry plane perpendicular to the loading direction.A small (less than 20%) reduction in in-plane shear strength is observed in all plies.Although the laminate is loaded in its plane, reduction of in-plane properties leads to transverse shear stresses (captured by the 3D Figures [14][15][16][17][18] show the relative reduction of stiffness and strength properties due to fatigue damage that exceeded 10%.The significant reductions (80-100%) are observed in the transverse to fiber (matrix) direction moduli E 2 and tension strength Y T and the plane shear moduli G 12 in the 90 • ply along the symmetry plane perpendicular to the loading direction.A small (less than 20%) reduction in in-plane shear strength is observed in all plies.Although the laminate is loaded in its plane, reduction of in-plane properties leads to transverse shear stresses (captured by the 3D modeling) that induce damage and reduction in transverse shear strengths S 23 .Although none of the specimen experienced complete failure under tension-tension loading after a million cycles, radiograph images show significant damage.The Figures 14-18 show the X-ray radiograph images of the damage accumulated in the specimen from Han et al. [34] at load level of 40-80% of the NTS.Significant damage is visible at 50%, 60% and 70% of NTS, even though the change in residual strength is not significantly high.At 80% of the NTS the damage state is significantly visible in the specimen as several modes of failure can cause the damage.In addition to the matrix damage the specimen also experienced delamination mode of failure.
Figures 19 and 20 show the evolution of damage quantified by the relative reductions in the transverse tensile strength and in-plane shear stiffness of the ply.The progression shows that the damage starts early (10 3 cycles) and continues to grow in a stable fashion.It is to be noted that while stiffness degradations redistribute loads to other plies leading to small decreases in stress.This has the potential to arrest the progression damage.The degradation in matrix properties do not seem to play a significant role on the overall stiffness or strength of the notched specimen.This is because the transverse (matrix) stiffness is one order of magnitude lower than fiber direction stiffness in our laminate.However, quantifying matrix damage is useful as matrix damage is undesirable as it can lead to delamination.The matrix damage also makes the composite prone to moisture intrusion or leakage.The effect of the interlaminar failure and damage patterns are shown by researchers for open hole composite specimens under tensile loading.Archard et al. [35] presented a discrete ply model approach to experimentally quantify damage under static loading for a composite plate with a hole under tensile loading.Nixon-Pearson et al. [36] presented an extensive experimental framework to show the sequence of damage development throughout the life of open-hole quasi-isotropic IM7/8552 carbon fiber/epoxy laminates loaded in tension-tension fatigue.Work of Nixon-Pearson et al. [36] showed that there is delamination between 0/45 and 45/90 layers during fatigue cycles.For the open hole tension specimen, we have depicted the reduction in the strength and stiffness     The effect of the interlaminar failure and damage patterns are shown by researchers for open hole composite specimens under tensile loading.Archard et al. [35] presented a discrete ply model approach to experimentally quantify damage under static loading for a composite plate with a hole under tensile loading.Nixon-Pearson et al. [36] presented an extensive experimental framework to show the sequence of damage development throughout the life of open-hole quasi-isotropic IM7/8552 carbon fiber/epoxy laminates loaded in tension-tension fatigue.Work of Nixon-Pearson et al. [36] showed that there is delamination between 0/45 and 45/90 layers during fatigue cycles.For the open hole tension specimen, we have depicted the reduction in the strength and stiffness properties in each ply and determined the residual strength of the specimen after a million fatigue cycles.The ultimate strength was determined using Hashin's failure criteria where right after fiber failure the stiffness E 1 was reduced to 7% of its initial value (Table 3).Thus, the effect of delamination to predict the ultimate residual strength of the specimen can be neglected mostly due the nature of the intralaminar failure of the specimen for predicting residual strength.
The methodology presented in this paper was specific to the model and material parameters presented by Shokrieh and Lessard [17,18].The primary goal of this paper was to develop and present a FE framework by using a UMAT subroutine to predict residual strength of the open hole specimen and show intralaminar damage plots in each ply.The number of parameters associated with the model can be significant sources of uncertainty.The model contains forty two parameters that is required to be calibrated from unidirectional tests.Thus, a separate study for the uncertainty quantification is required for the model.Using high fidelity model such as a continuum damage mechanics (CDM) based model has shown that fewer number of such fatigue parameters [37].The methodology (FE framework) can be extrapolated to any sets of material provided the material parameters (Table 2) for static and fatigue tests are available from experiment.Shokrieh's [6,16,17] model develops the empirical S-N curves, stiffness and strength degradation using uniaxial tests along each material direction and uses these models for multiaxial fatigue.However, the model has some limitations.The model does not capture the correlation or interaction between the degradation in stiffness and the residual strength parameters along different material directions of the ply.Matrix cracking failure should cause degradation of both the transverse normal and shear properties.In this model, the transverse to fiber property degradation is only due to transverse stress and the in-plane shear property degradation is only due to shear loads.The matrix cracking failures in the plies is induced by both transverse to fiber normal stresses and shear stresses.The model can be improved fitting the empirical models to the evolution of damage variables with these damage variables used to compute the constitutive stiffness and strength properties following a continuum damage approach.The models used for fatigue analysis of composites can be categorized by the fidelity of the model.In most of the cases empirical models are considered as lower fidelity model with less computational time but often having large number of fatigue parameters.Higher fidelity model includes the Continuum Damage Mechanics based model and micro-mechanics based models that are computationally expensive.Computationally cheap lower fidelity models can significantly aid in multi-level and multifidelity modeling approaches for uncertainty quantification [38].Our revisit of the Shokrieh and Lessard model [17,18] was for its use as a low fidelity model in a multifidelity framework.However, the inherent typographical errors in the original publication and ensuing citations to the paper led to large discrepancies.In this paper, we present the corrected version for these fit parameters reported by Shokrieh and Lessard [17,18].As a result of the corrections the low fidelity model provides an improved prediction accuracy that can be further refined and calibrated with other fidelity models.

Conclusions
This paper presents the empirical stiffness and strength degradation model of Shokrieh's [6], its implementation as a user material subroutine UMAT in ABAQUS FE software and its validation using test data from literature.The results indicate this empirical model can provide a good engineering estimate of fatigue life and residual strength for preliminary design analysis.The modification in the definition of the unified load parameter for the transverse to the fiber direction improved the prediction accuracy.
We found that residual strength for open-hole tension specimen is mostly related to fatigue failure in the fiber failure in 0 degree plies.Our developed user defined material subroutine (UMAT) enabled us to predict residual strength of the open hole tension specimen without the inclusion of any delamination model.The translaminar failure is captured by the degradation of the transverse normal and shear moduli of the lamina.The inclusion of the cumulative damage model adds to the capabilities to account for the stress level changes due to the failure during fatigue loading.
FE based progressive failure analysis with instantaneous degradation is highly mesh size sensitive and requires mesh convergence studies.In uncertainty quantification, the discretization errors must therefore be included as a source of model uncertainty.Further, the value of the residual stiffness used after ply failure has a significant influence on the peak load.For example, changing the residual stiffness from the values listed in Table 1 to a uniform 4% value reduces the NTS and RTS values by 25%.The lack of correlation between degradations to transverse to fiber direction normal modulus and in plane shear due to matrix damage is also a source of model error.These two are the most significant source of epistemic uncertainty that needs higher fidelity models and test data to correct.

Figure 1 .
Figure 1.(a) Schematic of a typical constant life diagram; (b) Normalized fatigue life model derived from constant life diagram.

Figure 1 .
Figure 1.(a) Schematic of a typical constant life diagram; (b) Normalized fatigue life model derived from constant life diagram.

Figure 3 .
Figure 3.The sudden death and wear out model for the power law fit used in Shokrieh's model.

Figure 3 .
Figure 3.The sudden death and wear out model for the power law fit used in Shokrieh's model.

Figure 4 .
Figure 4.The normalized residual strength model for the 0° specimen tested by Shokrieh and Lessard showing the effect of the parameters α and β obtained from the test data.

Figure 4 .
Figure 4.The normalized residual strength model for the 0 • specimen tested by Shokrieh and Lessard showing the effect of the parameters α and β obtained from the test data.

Figure 5 .
Figure 5.A schematic showing the cumulative damage approach.Figure 5.A schematic showing the cumulative damage approach.

Figure 5 .
Figure 5.A schematic showing the cumulative damage approach.Figure 5.A schematic showing the cumulative damage approach.

20 Figure 6 .
Figure 6.Flow chart of the FEM implementation using UMAT subroutine.

Figure 6 .
Figure 6.Flow chart of the FEM implementation using UMAT subroutine.

Figure 8 .
Figure 8. Mesh used for the notched tensile test specimen analyses.

Figure 9 .
Figure 9. Stress-strain response of the quasi-isotopic OHT specimen in pristine state and post-fatigue tested state after 10 6 fatigue load cycles (R = 0.1) at different specified maximum stress levels.

Figure 8 .
Figure 8. Mesh used for the notched tensile test specimen analyses.

Figure 8 .
Figure 8. Mesh used for the notched tensile test specimen analyses.

Figure 9 .Table 5 .FatigueFigure 9 .
Figure 9. Stress-strain response of the quasi-isotopic OHT specimen in pristine state and post-fatigue tested state after 10 6 fatigue load cycles (R = 0.1) at different specified maximum stress levels.

Figure 10 .
Figure 10.Fiber tension failure index after 10 6 cycles in the quasi-isotropic Open Hole Tension (OHT) specimen.

Figure 10 .
Figure 10.Fiber tension failure index after 10 6 cycles in the quasi-isotropic Open Hole Tension (OHT) specimen.

Figure 11 .
Figure 11.Fiber compression failure index after 10 6 cycles in the quasi-isotropic OHT specimen.

Figure 11 .
Figure 11.Fiber compression failure index after 10 6 cycles in the quasi-isotropic OHT specimen.

Figure 11 .
Figure 11.Fiber compression failure index after 10 6 cycles in the quasi-isotropic OHT specimen.

Figure 12 .
Figure 12.Matrix tension failure index after 10 6 cycles in the quasi-isotropic OHT specimen.

Figure 13 .
Figure 13.Matrix compression failure index after 10 6 cycles in the quasi-isotropic OHT specimen.

Figure 14 .
Figure 14.Reduction in stiffness E2 in each ply of the OHT specimen after 10 6 cycles.Figure 14.Reduction in stiffness E 2 in each ply of the OHT specimen after 10 6 cycles.

Figure 14 .
Figure 14.Reduction in stiffness E2 in each ply of the OHT specimen after 10 6 cycles.Figure 14.Reduction in stiffness E 2 in each ply of the OHT specimen after 10 6 cycles.

Figure 14 .
Figure 14.Reduction in stiffness E2 in each ply of the OHT specimen after 10 6 cycles.

Figure 15 .
Figure 15.Reduction in stiffness G12 in each ply of the OHT specimen after 10 6 cycles.

Figure 16 .
Figure 16.Reduction in strength Y T in each ply of the OHT specimen after 10 6 cycles.

Figure 16 .
Figure 16.Reduction in strength YT in each ply of the OHT specimen after 10 6 cycles.

Figure 17 .
Figure 17.Reduction in strength S12 in each ply of the OHT specimen after 10 6 cycles.

Figure 18 .
Figure 18.Reduction in strength S 23 in each ply of the OHT specimen after 10 6 cycles.

Figure 18 .
Figure 18.Reduction in strength S23 in each ply of the OHT specimen after 10 6 cycles.

Figure 19 .
Figure 19.Reduction in strength YT in 90° plies of the OHT specimen for increasing fatigue cycles.

Figure 20 .
Figure 20.Reduction in stiffness G12 in 90° plies of the OHT specimen for increasing fatigue cycles.

Figure 19 .
Figure 19.Reduction in strength Y T in 90 • plies of the OHT specimen for increasing fatigue cycles.

Figure 18 .
Figure 18.Reduction in strength S23 in each ply of the OHT specimen after 10 6 cycles.

Figure 19 .
Figure 19.Reduction in strength YT in 90° plies of the OHT specimen for increasing fatigue cycles.

Figure 20 .
Figure 20.Reduction in stiffness G12 in 90° plies of the OHT specimen for increasing fatigue cycles.

Figure 20 .
Figure 20.Reduction in stiffness G 12 in 90 • plies of the OHT specimen for increasing fatigue cycles.

Table 1 .
A comparison to the values of the unified load parameters.

Table 1 .
A comparison to the values of the unified load parameters.

Table 3 .
Material properties to be degraded based on the mode of failure.

Table 5 .
[34]arison of the residual strength of the Open Hole Tension (OHT) specimen with test data[34].

Table 5 .
[34]arison of the residual strength of the Open Hole Tension (OHT) specimen with test data[34].