A continuum damage mechanics model for fatigue and degradation of fiber reinforced materials

Abstract

Objective of the present study is the definition of a continuum damage mechanics material model describing the degradation of fiber reinforced materials under fatigue loads up to final failure. Based on the linear elastic framework, a brittle damage model for fatigue conditions is derived, where the damage constitutes the only nonlinearity. The model accounts for damage effects by successive degradation of the elastic moduli. Assuming that material damage is driven by microplastic work, a stress-driven damage evolution equation is defined. For generality, a fully three-dimensional formulation on single ply level is employed. The model is implemented into a finite element program. In a validation against experimental data on filament-wound carbon fiber reinforced material, the model proves to provide a good numerical approximation of the damage during the cyclic loading history up to final material failure.

Introduction

Fiber reinforced plastics of the carbon or glass fiber reinforced (CFRP or GFRP) types, as well as other types, are important structural materials in all fields of lightweight construction. From the classical field of the aerospace industry, their area of application today ranges over all kinds of road and rail transportation, the naval industry and even civil engineering applications to the field of wind energy and numerous types of special applications where high stiffness and strength of the material at low specific weight is required. Especially, in transport applications and for all kinds of machinery, the assessment of the fatigue limits of the material in addition to the static load limits is an essential part of the integrity assessment.

Due to their distinct microstructure, the fatigue response of fiber reinforced plastics is more complex compared to metals and other structural materials. Furthermore, it depends on the type of the fiber and matrix microstructure, whether of the continuous unidirectional, woven or nonwoven, random or aligned short fiber reinforced type. It has been shown by Talreja ¹ in his now classical contribution that the failure of fiber reinforced plastics under cyclic loads within the fiber direction is governed by static fiber failure, matrix and interface cracking. In this context, the matrix is subject to fatigue effects whereas the fibers exhibit only minor fatigue effects or even no fatigue degradation at all. For high-modulus fibers, static fiber failure is the leading failure mechanism since the stiff fibers shield the matrix from deformation such that no fatigue response in the classical sense is obtained. For softer fibers, fatigue degradation of the matrix becomes the leading mechanism. As a result, a classical Wöhler-type ² S-N-curve response is obtained. Experimental investigations on GFRP materials, by Talreja, ³ yield that up to 200 load cycles, fiber failure is the dominant mechanism whereas matrix cracking normal to the loading direction becomes the dominant mechanism at the later stages. Under off-axis loads, the leading failure mechanism turns to matrix and interface degradation, already for small load and fiber angle misalignments. Similar mechanisms are active in short fiber composites, especially in case of randomly oriented fibers. In multidirectional laminates, matrix failure in the form of inter-ply delamination constitutes an additional active failure mechanism.^3,4

A specific effect in the fatigue of fiber reinforced (and other) composite materials is the possible development of a significant stiffness degradation in the initial range of the cyclic loading history. Usually, this initial degradation comes to a stand-still after several hundreds of load cycles such that the initial rapid degradation zone is followed by a long range area with (nearly) constant stiffness, lasting through most of the fatigue lifetime. When the fatigue limit is approached, a rapidly increasing development of stiffness degradation up to a complete loss of stiffness at material failure is observed. Using computed tomography, this effect has been studied on the microscopic level by Arif et al. ⁵ for injection molded short fiber reinforced PA66GF30 as well as by Garcea et al. ⁶ on $[{\hspace{0.17em}}^0\circ ,{\hspace{0.17em}}^{90}\circ ]{\hspace{0.17em}}_n$ cross-ply laminates. In both cases, an early development of micro defects at stress concentrations due to interaction of fibers or initial defect triggers has been observed. Since the formation of these micro defects is triggered by a local deformation constraint, which decreases due to the crack formation, the growth of the micro defects is strongly decelerated once they are formed. The macroscopic stiffness degradation effect is significant for all composite types with extensive fiber interactions such as short fiber composites.^5–7 Nevertheless, in weaker form, it is also observed in continuously unidirectional materials and laminates thereof.^8–10

In engineering application, the assessment of structures under fatigue load is typically performed in terms of the classical Wöhler-Miner^2,11 S-N-curve type approach. For this purpose, Bond ¹² proposed a modified Miner rule in conjunction with rainflow counting for GFRP materials. An alternative nonlinear Miner rule for short fiber composites has been proposed by Dreißig et al. ¹³ Brunbauer et al. ¹⁴ as well as Gude et al. ¹⁵ analyzed the fatigue of fiber reinforced materials using a load cycle dependent Puck's criterion. A similar approach based on the Tsai-Hill criterion has been provided by Jen and Lee. ¹⁶ Hartmann et al. ¹⁷ as well as Zago and Springer ¹⁸ investigated the precision of fatigue predictions based on damage accumulation rules of the generalized Miner type for different structural components experimentally and thereby proved their general applicability.

The advantage of generalized Wöhler-Miner type models is their high numerical efficiency since they can be applied during the postprocessing of static structural analyses. On the other hand, they inherently do not account for possible stiffness degradation during the fatigue process and thus do not account for possible stress redistribution due to inhomogeneous stiffness degradation. Furthermore, their application to load histories with strongly variable nonharmonic load cycles is difficult. For this purpose, continuum damage models tracking the local fatigue degradation as part of the constitutive model are advantageous, including these effects in a natural manner. An early continuum damage model for two-dimensional anisotropic composites has been presented by Matzenmiller et al. ¹⁹ Nevertheless, in a similar manner as a recent model proposed by Hund et al., ²⁰ the concern of the model is static damage development due to overloading rather than fatigue damage. A continuum damage mechanics model for fatigue of composites has been proposed by Kennedy et al. ²¹ Their model employs cycle number driven damage accumulation relations based on Puck's criterion in conjunction with the Miner rule, also accounting for the S-shaped stiffness degradation during the fatigue lifetime. A similar, but more simple model has recently been proposed by Mohammadi et al. ²²

Due to their load cycle driven formulation of the damage accumulation law, the models proposed by Kennedy et al. ²¹ and Mohammadi et al. ²² are restricted to at least blockwise harmonic fatigue loads in the same manner as the model proposed by Shokrieh and Taheri-Behrooz. ²³ The present study is concerned with an alternative model, applicable to all kinds of fatigue load histories. For this purpose, a three-dimensional generalization of Matzenmiller et al. ¹⁹ brittle damage model, coupled with a damage accumulation equation driven by microplastic dissipation is proposed. Assuming a Ramberg and Osgood ²⁴ type law for the microplastic deformation, all dissipative effects can be estimated from the elastic deformation, resulting in an efficient formulation for the constitutive equation where the damage accumulation provides the only nonlinear effect, thereby generalizing an earlier approach by the authors.^25,26 The model is implemented into the finite element method and validated against an experimental data base on a filament wound carbon epoxy composite material.

Damage degradation model

Basic assumptions and hypotheses

The objective of the present contribution is the definition of a continuum damage mechanics material model for fiber reinforced composites under harmonic or nonharmonic fatigue loads. The model should be able to describe the final material failure as well as a possible continuous stiffness degradation of the material, reported in the literature for some kinds of fiber reinforced plastics. The definition of the model is based on the following basic assumptions and hypotheses:

The basic material response is linear elastic as it is the case for classical unidirectionally CFRP epoxy matrix composites or similar materials. The only relevant nonlinearity to be considered is the degradation of the elastic moduli in the anisotropic formulation of Hooke's law.

The dominant damage mechanism is the formation, growth, and coalescence of microcracks in the matrix between the fibers or breakage of individual fibers and subsequent formation of microcracks in the neighboring matrix.

The evolution of damage is assumed to be driven by dissipation of microplastic work. In this context, microplasticity is defined as a contained plasticity below the macroscopic yield limit of the material with plastic strains, which are much smaller than their elastic counterparts. Hence, the microplastic deformation is negligible compared to the elastic deformation. Nevertheless, despite the rather small plastic strains and thus the rather small plastic dissipation in individual load cycles, repeated loading and unloading at larger numbers of load cycles may result in a cumulative dissipation of a considerable amount of plastic work, resulting in considerable material damage.

Stiffness degradation may occur from the beginning of the cyclic loading. The degradation of the material stiffness develops in a nonlinear manner with a possible distinct drop in the initial range, caused by the instantaneous growth of pre-existing micro-defects under the first loading sequences. The initial drop is followed by a long range with minor additional degradation. Prior to the final material failure, the second range is followed by a range with rapidly increasing stiffness degradation.

In coincidence with the anisotropy of the material's microstructure and effective elastic properties, all damage and degradation effects are suspected to be strongly anisotropic.

These hypotheses are assumptions deriving from experimental observations from literature. Their approximation quality needs to be checked by application of the material model to be derived against experimental data. Under consideration of the mentioned hypotheses, the continuum damage model is formulated on single ply level. By this means, interaction of neighboring plies and thus also fatigue driven delaminations are accounted for in a natural manner.

Damage variables

The damage effects on the macroscopic stress–strain response of the material are modeled in the Kachanov ²⁷ and Lemaitre ²⁸ sense. For this purpose, the macroscopic “effective” stresses are defined by

$$\overline{\sigma }=(1-D)\sigma$$

(1)

as a function of the microscopic stresses σ, which are interrelated with the strains $\varepsilon (=\overline{\varepsilon })$ by the constitutive relation for the undamaged material. The damage variable D forms a history variable as an additional field quantity defining the degree of damage, with

$$\begin{array}{c}0=D\text{initialundamagedstate}\\ 0<D<1\text{partiallydamagedstate}\\ D=1\text{materialfailure}\end{array}$$

(2)

Assuming that microcrack formation, growth and coalescence is the relevant damage mechanism, the damage parameter D is related to the ratio of the crack area oriented normally to the stress direction to the total area of the cross section. However, D is not identical to this relation due to the inhomogeneity of the material on microstructural level as well as the possible development of stress and strain concentrations at the crack fronts of the microcracks.

For the present problem of an anisotropic fiber reinforced material, three independent damage variables D₁, D₂, and D₃ for crack orientations normal to the three coordinate axes as sketched in Figure 1 are required. Using the standard definition of the coordinate frame with the x₁-direction coinciding with the fiber direction, D₂ and D₃ are related to the formation of interfiber microcracks (possibly also involving the fiber and matrix interfaces) whereas D₁ describes the effect of microcracks normal to the fiber orientation thus also accounting for fiber damage and breakage.

Figure 1.

Definition of damage variables.

Apart from the damage effects, the material is assumed to exhibit a linear elastic response and thus can be modeled by Hooke's law in its orthotropic form. Introducing the damage variable according to Figure 1 and using equation (1), the orthotropic Hooke's law can be written as

$$(\begin{array}{c}{\overline{\varepsilon }}_{11}\\ {\overline{\varepsilon }}_{22}\\ {\overline{\varepsilon }}_{33}\\ 2{\overline{\varepsilon }}_{23}\\ 2{\overline{\varepsilon }}_{13}\\ 2{\overline{\varepsilon }}_{12}\end{array})=(\begin{array}{cccccc}\frac{1}{(1-D_1)E_1}& \hfill -\frac{{\nu }_{21}}{E_2}\hfill & \hfill -\frac{{\nu }_{31}}{E_3}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ -\frac{{\nu }_{12}}{E_1}& \hfill \frac{1}{(1-D_2)E_2}\hfill & \hfill -\frac{{\nu }_{32}}{E_3}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ -\frac{{\nu }_{13}}{E_1}& \hfill -\frac{{\nu }_{23}}{E_2}\hfill & \hfill \frac{1}{(1-D_3)E_3}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ 0& \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{1}{(1-D_2)(1-D_3)G_{23}}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ 0& \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{1}{(1-D_1)(1-D_3)G_{13}}\hfill & \hfill 0\hfill \\ 0& \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{1}{(1-D_1)(1-D_2)G_{12}}\hfill \end{array})(\begin{array}{c}{\overline{\sigma }}_{11}\\ {\overline{\sigma }}_{22}\\ {\overline{\sigma }}_{33}\\ {\overline{\sigma }}_{23}\\ {\overline{\sigma }}_{13}\\ {\overline{\sigma }}_{12}\end{array})$$

(3)

with the Young's and shear moduli E _i and G _ij as well as the Poisson's ratios ν _ij as the elasticity parameters for the undamaged material.^25,26 Retaining the symmetry of the compliance matrix, equation (3) can be considered as Hooke's law, where the elastic constants are degraded according to

$${\overline{E}}_i=(1-D_i)E_i$$

(4)

$${\overline{G}}_{\mathit{ij}}=(1-D_i)(1-D_j)G_{\mathit{ij}}$$

(5)

$${\overline{\nu }}_{\mathit{ij}}=(1-D_i){\nu }_{\mathit{ij}}$$

(6)

by damage effects. In these equations as well as in all subsequent equations, Einstein's sum convention with respect to repeated indices is not applied. In equation (3), it is assumed that any of the effective shear moduli ${\overline{G}}_{\mathit{ij}}$ is degraded by the presence of microcracks in both planes normal to the respective shear plane and thus must vanish if either $D_i\to 1$ or $D_j\to 1$ . By this means, the introduction of independent shear damage variables can be avoided. Equation (6) follows from the requirement of symmetry of the compliance matrix in the undamaged and damaged material state and the substitution of the initial elastic moduli E _i with their degraded counterparts ${\overline{E}}_i$ in the latter case. Notice that even if the material initially is transverse isotropic (i.e. $E_3=E_2$ , $G_{13}=G_{12}$ , ${\nu }_{31}={\nu }_{21}$ , and ${\nu }_{13}={\nu }_{12}$ ), the transverse isotropy in general is lost as soon as damage effects with nonzero D₂ and D₃ develop. A similar model has been proposed by Matzenmiller et al. ¹⁹ Nevertheless, this model is restricted to the two-dimensional case and introduces independent damage variables to the shear moduli.

Damage evolution

One-dimensional considerations

The three damage variables D _i in equation (3) form additional field variables for which additional equations are required. These equations are provided in terms of evolution equations, describing the evolution of damage from a prescribed initial value to the actual value in the computed time history. According to the considerations in Section “Basic assumptions and hypotheses,” it is assumed that the increase $\text{d}D$ in damage is controlled by the amount $\text{d}{w}^{\text{mp}}$ of microplastic work density dissipated during a time increment $\text{d}t$

$$\text{d}D=\text{d}D(w^{\text{mp}})$$

(7)

In this context, the microplastic deformation is a plastic deformation induced by limited microscopic gliding of polymeric chains in the matrix material, occurring at load levels below the macroscopic yield limit. Since these gliding effects are rather limited, the resulting microplastic strains are small compared to their elastic counterparts

$${\varepsilon }^{\text{mp}}<<{\varepsilon }^{\text{el}}$$

(8)

so that—as illustrated in Figure 2—they have no significance on the macroscopic level and thus do not need to be considered in the macroscopic constitutive equation (3).

Figure 2.

Concept of micro plasticity.

Assuming that the microplastic strains in one-dimensional form can be approximated by the Ramberg and Osgood ²⁴ equation

$${\varepsilon }^{\text{mp}}=\alpha (\frac{\sigma }{E}{)}^n$$

(9)

the microplastic strain increments are defined by

$$\text{d}{w}^{\text{mp}}=\sigma 1\mathit{mu}\text{d}\varepsilon =\frac{n\alpha }{E^n}1\mathit{mu}{\sigma }^{n}1\mathit{mu}\text{d}\sigma$$

(10)

where α and n are material parameters. Hence, the damage increment

$$\text{d}D=\text{d}D(\text{d}{w}^{\text{mp}})=\text{d}D({\sigma }^{n}1\mathit{mu}\text{d}\sigma )$$

(11)

can be estimated from the stress increment $\text{d}\sigma$ together with the actual stress level σ. Recombining the material parameters, the one-dimensional damage evolution equation

$$\text{d}D=A1\mathit{mu}\omega (D)1\mathit{mu}\left|\sigma \right|{\hspace{0.17em}}^{n}1\mathit{mu}\text{d}\sigma$$

(12)

is obtained. In equation (12), A and n are material parameters. The absolute values of stress and stress increment are used to deal with negative stresses. The warping function $\omega (D)$ is introduced, since otherwise a linear increase of the damage with the number N of load cycles would be obtained for harmonic loading. In order to comply with the experimental observations in the literature, a heuristic approach with the piecewise quadratic definition

$$\omega (D)=\begin{array}{cc}b+a_{\text{I}}{(D-D_{\text{I}})}^2& \hfill \text{ for }D\le D_{\text{I}}\hfill \\ b& \hfill \text{ for }{D}_{\text{I}}<D\le D_{\text{II}}\hfill \\ b+a_{\text{II}}{(D-D_{\text{II}})}^2& \hfill \text{ for }D>D_{\text{II}}\hfill \end{array}$$

(13)

is employed (Figure 3). The material parameters $D_{\text{I}},D_{\text{II}},a_{\text{I}},a_{\text{II}}$ , and b have to be chosen such that the predicted damage evolution over the lifetime coincides with the experimental observations, considering that

$${\int }_0^1\omega (D)1\mathit{mu}\text{d}D=1$$

(14)

Figure 3.

Warping function.

The definition (13) of the warping function results in an accelerated damage evolution for small degrees D of damage (i.e. in the initial range of fatigue loading) as well as for large degrees D of damage when the instant of final failure is approached.

Discussion

For illustration, the meaning of the one-dimensional damage evolution equation (12) is investigated in more detail. Consider for the moment a simplified form where the warping of the damage evolution according to Figure 3 is discarded by assuming $\omega (D)=1$ , the total increase of damage in the case of a stress increase from a non-negative σ₀ to a positive σ₁ is obtained as

$$\Delta D={\int }_{{\sigma }_0}^{{\sigma }_1}A1\mathit{mu}{\sigma }^{n}1\mathit{mu}\text{d}\sigma$$

(15)

and thus

$$\Delta D=\frac{A}{n+1}({\sigma }_1^{n+1}-{\sigma }_0^{n+1})$$

(16)

If a cyclic load is applied at a stress ratio of R = 0, the damage accumulation per stress cycle is obtained as

$$\Delta D^{\text{cycl}}=\frac{2A}{n+1}1\mathit{mu}\Delta {\sigma }^{n+1}$$

(17)

where it is assumed that an identical accumulation of damage is experienced during the loading and unloading parts of the cycle.

For a material failure after $N^{\text{fail}}$ cycles, a damage accumulation of $\Delta D^{\text{cycl}}=1/N^{\text{fail}}$ per cycle is necessary. Thus

$$\Delta \sigma =(\frac{n+1}{2A}{)}^{\frac{1}{n+1}}(N^{\text{fail}}{)}^{-\frac{1}{n+1}}$$

(18)

is obtained by inversion of equation (17). Taking the logarithm of equation (18) yields

$$\mathit{ln}\Delta \sigma =\frac{1}{n+1}\mathit{ln}(\frac{n+1}{2A})-\frac{1}{n+1}\mathit{ln}{N}^{\text{fail}}$$

(19)

Hence, in a double logarithmic representation, a linear relation between the logarithm of the stress amplitude and the logarithm of the numbers of cycles to failure is obtained as it is postulated in the S-N-curve concept. ² By this interrelation with the well-established Wöhler-concept, the meaningful definition (12) of the one-dimensional damage evolution equation based on the microplasticity concept together with the assumption of a Ramberg–Osgood type power-law relation between stresses and microplastic strains becomes evident.

Three-dimensional formulation

For the numerical implementation of the fatigue degradation model, the one-dimensional form (12) has to be re-written to a fully three-dimensional version. For this purpose, it is assumed that each of the damage variables D₁, D₂, and D₃ related to the three independent spatial directions x₁, x₂, and x₃ is driven by all stress components σ _ij , which act on the respective surfaces of the volume element in Figure 1. Assuming an additive decomposition of the effects of the individual stress components, the three-dimensional generalizations

$$\begin{array}{c}\text{d}{D}_1=A_{11}\omega (D_1)\left|{\sigma }_{11}\right|{\hspace{0.17em}}^{n_{11}}\left|\text{d}{\sigma }_{11}\right|+A_{12}\omega (D_1)\\ \times \left|{\sigma }_{12}\right|{\hspace{0.17em}}^{n_{12}}\left|\text{d}{\sigma }_{12}\right|+A_{13}\omega (D_1)\left|{\sigma }_{13}\right|{\hspace{0.17em}}^{n_{13}}\left|\text{d}{\sigma }_{13}\right|\end{array}$$

(20)

$$\begin{array}{c}\text{d}{D}_2=A_{21}\omega (D_2)\left|{\sigma }_{12}\right|{\hspace{0.17em}}^{n_{21}}\left|\text{d}{\sigma }_{12}\right|+A_{22}\omega (D_2)\\ \times \left|{\sigma }_{22}\right|{\hspace{0.17em}}^{n_{22}}\left|\text{d}{\sigma }_{22}\right|+A_{23}\omega (D_2)\left|{\sigma }_{23}\right|{\hspace{0.17em}}^{n_{23}}\left|\text{d}{\sigma }_{23}\right|\end{array}$$

(21)

$$\begin{array}{c}\text{d}{D}_3=A_{31}\omega (D_3)\left|{\sigma }_{13}\right|{\hspace{0.17em}}^{n_{31}}\left|\text{d}{\sigma }_{13}\right|+A_{32}\omega (D_3)\\ \times \left|{\sigma }_{23}\right|{\hspace{0.17em}}^{n_{32}}\left|\text{d}{\sigma }_{23}\right|+A_{33}\omega (D_3)\left|{\sigma }_{33}\right|{\hspace{0.17em}}^{n_{33}}\left|\text{d}{\sigma }_{33}\right|\end{array}$$

(22)

of equation (12) are obtained for the evolution of the damage variables. The parameters A _ij and n _ij are material parameters. Notice that in general $A_{\mathit{ij}}\ne A_{\mathit{ji}}$ and $n_{\mathit{ij}}\ne n_{\mathit{ji}}$ . Since the effect of microcracks under tension and compression is inherently different, the damage evolution of the corresponding damage variables under normal tensile and compressive stresses in general will be different. In order to account for this effect, the corresponding material parameters

$$A_{1}1=\begin{array}{c}{A}_{11}^{\text{t}}\text{ for }{\sigma }_{11}\ge 0\\ A_{11}^c\text{ for }{\sigma }_{11}<0\end{array}$$

(23)

$$A_{22}=\begin{array}{c}{A}_{22}^{\text{t}}\text{ for }{\sigma }_{22}\ge 0\\ A_{22}^c\text{ for }{\sigma }_{22}<0\end{array}$$

(24)

$$A_{33}=\begin{array}{c}{A}_{33}^{\text{t}}\text{ for }{\sigma }_{33}\ge 0\\ A_{33}^{\text{c}}\text{ for }{\sigma }_{33}<0\end{array}$$

(25)

entering equations (20) to (22) are set to different values, depending on the direction of the corresponding normal stress component.

For simplicity and in order to keep the overall number of material parameters within acceptable limits, it is assumed that the warping function $\omega (D)$ is the same for all damage parameters and stress components. Since this function only controls the warping of the damage evolution, i.e. controls the amount of acceleration of the damage evolution in the initial and final stages, this assumption can be made without any loss in generality concerning the stress changes and thus the resulting numbers of cycles to failure for the different loading modes which are controlled by the parameters A _ij . As shown in equation (19) in Section “One-dimensional considerations,” the exponents n _ij control the slope of S-N-curve and thus the amounts of damage accumulated for similar stress changes $\text{d}{\sigma }_{\mathit{ij}}$ at different actual stress levels $\left|{\sigma }_{\mathit{ij}}\right|$ . For n _ij  = 0, a linear increase of the damage variable with the change in the respective stress component σ _ij is obtained whereas all $n_{\mathit{ij}}>0$ result in an increased damage accumulation $\text{d}{D}_i$ for stress changes $\left|\text{d}{\sigma }_{\mathit{ij}}\right|$ occurring at higher actual stress levels $\left|{\sigma }_{\mathit{ij}}\right|$ compared to the effect of identical stress changes at lower actual stress levels.

Numerical implementation

Implicit time integration

For application in fatigue assessment, the material model defined in Sections “Basic assumptions and hypotheses,” “Damage variables,” and “Damage evolution,” is implemented into the finite element method. Employing an incremental approach, the stress components ${\sigma }_{\mathit{ij}}^{(k+1)}$ and damage variables $D_i^{(k+1)}$ for the time step $(k+1)$ are computed from their (known) counterparts ${\sigma }_{\mathit{ij}}^{(k)}$ and $D_i^{(k)}$ at the previous time step k by solution of the nonlinear system

$$0=\mathbf{C}(D_i^{(k+1)})·{\hspace{0.17em}}^{(k+1)}-{\sigma }^{(k+1)}$$

(26)

$$0={\mathbf{D}}^{(k)}+\mathbf{D}({\sigma }_{\mathit{ij}}^{(k+1)},D_i^{(k+1)})-{\mathbf{D}}^{(k+1)}$$

(27)

where the backward Euler rule is employed with respect to the time integration of the damage increments. In equations (26) and (27), the stress and strain vectors are defined by

$$\sigma =({\sigma }_{11},{\sigma }_{22},{\sigma }_{33},{\sigma }_{23},{\sigma }_{13},{\sigma }_{12}{)}^{\text{T}}$$

(28)

$$=({\varepsilon }_{11},{\varepsilon }_{22},{\varepsilon }_{33},2{\varepsilon }_{23},2{\varepsilon }_{13},2{\varepsilon }_{12}{)}^{\text{T}}$$

(29)

whereas the stiffness matrix C is obtained by inversion of equation (3). The damage vector is defined directly by

$$\mathbf{D}=(D_1,D_2,D_3{)}^{\text{T}}$$

(30)

whereas the finite damage increments

$$\mathbf{D}=(\Delta D_1,\Delta D_2,\Delta D_3{)}^{\text{T}}$$

(31)

are determined from the infinitesimal increments (20) to (22) by substituting the infinitesimal stress increments $\text{d}{\sigma }_{\mathit{ij}}$ with their finite counterparts $\Delta {\sigma }_{\mathit{ij}}={\sigma }_{\mathit{ij}}^{(k+1)}-{\sigma }_{\mathit{ij}}^{(k)}$ . For each time step and finite element integration point, the nonlinear system (26) to (27) is solved numerically by the Newton–Raphson method.

The application of the proposed fatigue damage model in structural fatigue assessment using the numerical implementation defined by equations (26) to (31) requires the simulation of the entire cyclic loading history. By this means, the model fully accounts for load history effects, random load amplitudes or other changes in the loading situation during the fatigue lifetime in a natural manner. Possible stress re-distributions due to an initial degradation of the material (Figure 3) or during the accelerating degradation in the vicinity of the final material failure are accounted for. The case of a stable fatigue crack growth is also accessible. Temperature change effects are accounted for by defining all material parameters to be temperature dependent, assuming a piecewise linear variation between discrete temperatures.

In order to keep the numerical effort for computation of the cyclic load history within acceptable bounds, a cycle jump technique is employed by restricting the numerical simulation to a number of individual cycles, each of which is representative for an individual part of the loading history. The damage evolution computed during these reference cycles is then multiplied by the number of cycles, for which the computation cycle is representative. In the numerical implementation, this technique is easily implemented by multiplying the damage increments $\mathbf{D}$ according to equation (31) by the number of “real” cycle, for which the respective computational load cycle is assumed to be representative.

Due to the power law formulation (12) of the damage evolution equation, resulting in equations (20) to (22), damage is accumulated especially during those portions of the stress cycles featuring high stress levels. Therefore, care has to be taken to ensure that these portions of the load cycles are supplied with a sufficient number of time steps in order to ensure that the high stress parts of the cycles are resolved in a sufficiently fine manner. In the computations in Section “Numerical example,” it is observed that for harmonic oscillations usually eight equally spaced time steps per load cycle are the minimum number of time steps necessary for a convergent solution, except for cycles with rather high damage accumulation.

Explicit time integration

The numerical implementation with implicit time integration as defined in Section “Implicit time integration” requires an iterative solution of the material equations (26) and (27). The iterative solution has to be computed for each global iteration in each time step at each integration point in the problem. Thus, a significant numerical effort is required, especially since typically, a larger number of cycles need to be simulated with a sufficient number of time steps each. On the other hand, the material nonlinearity requiring this local iterative procedure in the present case is only weak. The only nonlinearity is caused by the damage evolution, which in most cases is limited in the individual time steps forming the individual load cycles.

For this purpose, an alternative implementation using the forward Euler time integration scheme is proposed by substituting equations (26) to (27) by

$$0=\mathbf{C}(D_i^{(k+1)})·{\hspace{0.17em}}^{(k+1)}-{\sigma }^{(k+1)}$$

(32)

$$0={\mathbf{D}}^{(k)}+\mathbf{D}({\sigma }_{\mathit{ij}}^{(k)},D_i^{(k)})-{\mathbf{D}}^{(k+1)}$$

(33)

using explicit time integration. Since all nonlinear terms in the damage evolution equations (33) in this system are evaluated at the time step k rather than $(k+1)$ , no iterative procedure is necessary for the update of the history variable vector D . Once the damage variables are updated, the components of $D_i^{(k+1)}$ are substituted into equation (32), then providing the stress vector ${\sigma }^{(k+1)}$ at the end of the $k^{\text{th}}$ time step in a straightforward manner without the need for the solution of a nonlinear system of equations.

Due to the unavoidable requirement of small time steps during a load cycle in order to ensure that the temporal resolution of the high-stress ranges of the load cycles is sufficiently fine, the requirements imposed to the temporal discretization by the use of an explicit time integration do not impose an additional constraint. By avoidance of the need for an iterative solution of the material equations, a numerically much more efficient simulation procedure is obtained. Notice that explicit time integration is applied only on the material level, whereas a classical implicit finite element scheme is applied on the structural level.

Analytical damage accumulation rule

An alternative approach to the direct numerical implementation of the proposed damage model is a reformulation towards a generalized Miner ¹¹ approach, which then can be applied in the postprocessing of a classical linear elastic structural analysis without consideration of any damage and degradation effects during the underlying linear elastic stress analysis.

For this purpose, a harmonic fatigue load oscillating between the minimum and maximum local stress levels ${\sigma }_{\mathit{ij}}^{\text{max}}$ and ${\sigma }_{\mathit{ij}}^{\text{min}}$ , respectively, is assumed. In this context, the damage accumulation per cycle due to variations in the stress component σ _ij is given by

$$\Delta D_i=\int {\sigma }_{\mathit{ij}}^{\text{min}}{\sigma }_{\mathit{ij}}^{\text{max}}\text{d}{D}_i$$

(34)

Using this definition together with the general damage evolution equations (20) to (22) results in

$$\begin{array}{c}\Delta D_1=\frac{A_{11}}{n_{11}+1}(\begin{array}{c}({\sigma }_{11}^{\text{max}(+)}{)}^{n_{11}+1}\\ -({\sigma }_{11}^{\text{min}(+)}{)}^{n_{11}+1}\end{array})\\ +\frac{A_{12}}{n_{12}+1}(\begin{array}{c}\text{sign}({\sigma }_{12}^{\text{max}})|{\sigma }_{12}^{\text{max}(+)}{|}^{n_{12}+1}\\ -\text{sign}({\sigma }_{12}^{\text{min}})|{\sigma }_{12}^{\text{min}(+)}{|}^{n_{12}+1}\end{array})\end{array}$$

(35)

$$\begin{array}{c}+\frac{A_{13}}{n_{13}+1}(\begin{array}{c}\text{sign}({\sigma }_{13}^{\text{max}})|{\sigma }_{13}^{\text{max}(+)}{|}^{n_{13}+1}\\ -\text{sign}({\sigma }_{13}^{\text{min}})|{\sigma }_{13}^{\text{min}(+)}{|}^{n_{13}+1}\end{array})\\ \Delta D_2=\frac{A_{21}}{n_{21}+1}(\begin{array}{c}\text{sign}({\sigma }_{12}^{\text{max}})|{\sigma }_{12}^{\text{max}(+)}{|}^{n_{21}+1}\\ -\text{sign}({\sigma }_{12}^{\text{min}})|{\sigma }_{12}^{\text{min}(+)}{|}^{n_{21}+1}\end{array})\\ +\frac{A_{22}}{n_{22}+1}(({\sigma }_{22}^{\text{max}(+)})n_{22}+1-({\sigma }_{22}^{\text{min}(+)})n_{22}+1)\end{array}$$

(36)

$$\begin{array}{c}+\frac{A_{23}}{n_{23}+1}(\begin{array}{c}\text{sign}({\sigma }_{23}^{\text{max}})|{\sigma }_{23}^{\text{max}(+)}{|}^{n_{23}+1}\\ -\text{sign}({\sigma }_{23}^{\text{min}})|{\sigma }_{23}^{\text{min}(+)}{|}^{n_{23}+1}\end{array})\\ \Delta D_3=\frac{A_{31}}{n_{31}+1}(\begin{array}{c}\text{sign}({\sigma }_{13}^{\text{max}})|{\sigma }_{13}^{\text{max}(+)}{|}^{n_{31}+1}\\ -\text{sign}({\sigma }_{13}^{\text{min}})|{\sigma }_{13}^{\text{min}(+)}{|}^{n_{31}+1}\end{array})\\ +\frac{A_{32}}{n_{32}+1}(\begin{array}{c}\text{sign}({\sigma }_{23}^{\text{max}})|{\sigma }_{23}^{\text{max}(+)}{|}^{n_{32}+1}\\ -\text{sign}({\sigma }_{23}^{\text{min}})|{\sigma }_{23}^{\text{min}(+)}{|}^{n_{32}+1}\end{array})\\ +\frac{A_{33}}{n_{33}+1}(({\sigma }_{33}^{\text{max}(+)})n_{33}+1-({\sigma }_{33}^{\text{min}(+)})n_{33}+1)\end{array}$$

(37)

where the warping functions $\omega (D_i)$ have been neglected. The coefficients A₁₁, A₂₂, and A₃₃ are again defined to be dependent on the stress direction through equations (23) to (25). The minimum and maximum stress components entering equations (35) to (37) are defined as

$$\begin{array}{c}{\sigma }_{11}^{\text{max}(+)}=\{\begin{array}{ccc}{\sigma }_{11}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{11}^{\text{max}}>0\hfill \\ 0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{11}^{\text{max}}\le 0\hfill \end{array},\\ {\sigma }_{11}^{\text{min}(+)}=\{\begin{array}{ccc}{\sigma }_{11}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{11}^{\text{min}}>0\hfill \\ 0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{11}^{\text{min}}\le 0\hfill \end{array}\end{array}$$

(38)

$$\begin{array}{c}{\sigma }_{22}^{\text{max}(+)}=\{\begin{array}{ccc}{\sigma }_{22}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{22}^{\text{max}}>0\hfill \\ 0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{22}^{\text{max}}\le 0\hfill \end{array}\\ {\sigma }_{22}^{\text{min}(+)}=\{\begin{array}{ccc}{\sigma }_{22}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{22}^{\text{min}}>0\hfill \\ 0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{22}^{\text{min}}\le 0\hfill \end{array}\end{array}$$

(39)

$$\begin{array}{c}{\sigma }_{33}^{\text{max}(+)}=\{\begin{array}{ccc}{\sigma }_{33}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{33}^{\text{max}}>0\hfill \\ 0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{33}^{\text{max}}\le 0\hfill \end{array},\\ {\sigma }_{33}^{\text{min}(+)}=\{\begin{array}{ccc}{\sigma }_{33}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{33}^{\text{min}}>0\hfill \\ 0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{33}^{\text{min}}\le 0\hfill \end{array}\end{array}$$

(40)

$$\begin{array}{c}{\sigma }_{23}^{\text{max}(+)}=\{\begin{array}{ccc}{\sigma }_{23}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{23}^{\text{max}}>0\hfill \\ -{\sigma }_{23}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{23}^{\text{max}}\le 0\hfill \end{array},\\ {\sigma }_{23}^{\text{min}(+)}=\{\begin{array}{ccc}{\sigma }_{23}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{23}^{\text{min}}>0\hfill \\ -{\sigma }_{23}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{23}^{\text{min}}\le 0\hfill \end{array}\end{array}$$

(41)

$$\begin{array}{c}{\sigma }_{13}^{\text{max}(+)}=\{\begin{array}{ccc}{\sigma }_{13}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{13}^{\text{max}}>0\hfill \\ -{\sigma }_{13}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{13}^{\text{max}}\le 0\hfill \end{array},\\ {\sigma }_{13}^{\text{min}(+)}=\{\begin{array}{ccc}{\sigma }_{13}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{13}^{\text{min}}>0\hfill \\ -{\sigma }_{13}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{13}^{\text{min}}\le 0\hfill \end{array}\end{array}$$

(42)

$$\begin{array}{c}{\sigma }_{12}^{\text{max}(+)}=\{\begin{array}{ccc}{\sigma }_{12}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{12}^{\text{max}}>0\hfill \\ -{\sigma }_{12}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{12}^{\text{max}}\le 0\hfill \end{array},\\ {\sigma }_{12}^{\text{min}(+)}=\{\begin{array}{ccc}{\sigma }_{12}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{12}^{\text{min}}>0\hfill \\ -{\sigma }_{12}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{12}^{\text{min}}\le 0\hfill \end{array}\end{array}$$

(43)

$$\begin{array}{c}{\sigma }_{11}^{\text{max}(-)}=\{\begin{array}{ccc}0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{11}^{\text{max}}>0\hfill \\ -{\sigma }_{11}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{11}^{\text{max}}\le 0\hfill \end{array},\\ {\sigma }_{11}^{\text{min}(-)}=\{\begin{array}{ccc}0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{11}^{\text{min}}>0\hfill \\ -{\sigma }_{11}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{11}^{\text{min}}\le 0\hfill \end{array}\end{array}$$

(44)

$$\begin{array}{c}{\sigma }_{22}^{\text{max}(-)}=\{\begin{array}{ccc}0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{22}^{\text{max}}>0\hfill \\ -{\sigma }_{22}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{22}^{\text{max}}\le 0\hfill \end{array},\\ {\sigma }_{11}^{\text{min}(-)}=\{\begin{array}{ccc}0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{22}^{\text{min}}>0\hfill \\ -{\sigma }_{22}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{22}^{\text{min}}\le 0\hfill \end{array}\end{array}$$

(45)

$$\begin{array}{c}{\sigma }_{33}^{\text{max}(-)}=\{\begin{array}{ccc}0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{33}^{\text{max}}>0\hfill \\ -{\sigma }_{33}^{\text{max}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{33}^{\text{max}}\le 0\hfill \end{array},\\ {\sigma }_{11}^{\text{min}(-)}=\{\begin{array}{ccc}0& \hfill \text{ for }:\hfill & \hfill {\sigma }_{33}^{\text{min}}>0\hfill \\ -{\sigma }_{33}^{\text{min}}& \hfill \text{ for }:\hfill & \hfill {\sigma }_{33}^{\text{min}}\le 0\hfill \end{array}\end{array}$$

(46)

The analytical damage accumulation rules (35) to (37) together with the minimum and maximum stress definition (38) to (46) can be applied in the same manner as the classical Miner ¹¹ rule to a numerically efficient life time analysis using the stress output of a static structural analysis of a single load cycle without consideration of damage effects as the stress input. Assuming a damage free initial state, the damage after k cycles at each position is then estimated as $D_i=k1\mathit{mu}\Delta D_i$ . Once D _i  = 1 is reached, fatigue failure at the respective material point is assumed. In the case of block-wise applied load sequences, the damage accumulation may also be computed in a block-wise manner from the output of static structural analyses for the individual load blocks.

As an alternative to a pure postprocessing analysis of the fatigue damage accumulation in the Miner ¹¹ sense, the analytical damage accumulation rules (35) to (37) may also be used in a simplified cumulative analysis as proposed by Shokrieh and Taheri-Behrooz. ²³ In this context, the stress response in the individual load cycles is computed in a simplified linear elastic analysis. Subsequently, the damage accumulation during the respective cycle is estimated from the damage accumulation rules. After completion of the load cycle, the accumulated amount of damage is applied at once and the analysis is continued with the updated stiffness matrix $\mathbf{C}(D_i)$ .

Experimental data base

Reference material and methods

For validation of the damage degradation model developed in Section “Damage degradation model,” an experimental investigation has been performed. The reference material considered is a filament wound CFRP epoxy matrix composite as it is typically used for hydrogen or other composite pressure vessels. The fiber volume fraction is approximately ${\rho }^{\text{f}}=0.64$ . Plane plates for specimen manufacture were obtained by winding the material around a thick plane plate as a winding core. Both, uniaxially reinforced plates as well as multidirectionally reinforced plates with a $[\pm 17.{\hspace{0.17em}}^3\circ ,{\hspace{0.17em}}^{90}\circ ,\pm {\hspace{0.17em}}^{70}\circ ]{\hspace{0.17em}}_{\text{s}}$ stacking sequence were processed. The unidirectionally reinforced material has been manufactured to thicknesses of approximately $11\mathit{mu}\text{mm}$ for the tests within the fiber direction and approximately 2mm for all other experiments, respectively. For the multidirectionally reinforced plates, the ply thicknesses were ∼ $0.311\mathit{mu}\text{mm}$ for the $[\pm 17.{\hspace{0.17em}}^3\circ ]$ -plies, $0.481\mathit{mu}\text{mm}$ for the $[{\hspace{0.17em}}^{90}\circ ]$ -plies and two times $0.131\mathit{mu}\text{mm}$ for the center $[\pm {\hspace{0.17em}}^{70}\circ ]$ -ply. Notice that due to the filament winding process, the $[\pm \alpha ]$ -plies form interlaced layers rather than being two separate $[+\alpha ]$ - and $[-\alpha ]$ -plies.

From the available unidirectionally fiber reinforced material with $11\mathit{mu}\text{mm}$ thickness, specimens with outer dimensions of $2401\mathit{mu}\text{mm}\times 151\mathit{mu}\text{mm}$ were manufactured for tensile testing within the fiber direction. Some of the specimens were manufactured to a reduced width of $101\mathit{mu}\text{mm}$ in order to keep the required axial forces within acceptable limits. For the tensile experiments perpendicular to the fiber direction and the tensile experiments on multidirectional laminates, specimens with the outer dimensions of $2501\mathit{mu}\text{mm}\times 251\mathit{mu}\text{mm}$ and $2501\mathit{mu}\text{mm}\times 201\mathit{mu}\text{mm}$ , respectively, were machined from the thick unidirectionally fiber reinforced plates and the multidirectionally reinforced plates. For experiments under compression, specimens with outer dimensions $1551\mathit{mu}\text{mm}\times 101\mathit{mu}\text{mm}$ were machined from the thick unidirectionally fiber reinforced plates and the multidirectional laminates respectively.

The material has been characterized under quasi-static and cyclic loading conditions. For the quasi-static tests under tension, the specimens were clamped by a mechanical clamping system and tested in an electromechanical Hegewald and Peschke testing machine till failure. The experiments were performed under displacement control with a prescribed cross head velocity of $21\mathit{mu}\text{mm}/\text{min}$ . For the compressive tests, a specially designed parallel clamping test rig according to Figure 4 has been employed. The specimens were again tested under a prescribed cross head displacement rate of $21\mathit{mu}\text{mm}/\text{min}$ . All experiments were performed at ambient temperature. During the experiments, the cross head displacement, reaction force and the axial strain were continuously recorded and evaluated for the stress–strain response.

Figure 4.

Parallel clamping test rig.

The experiments under cyclic loads were performed in a hydraulic MTS testing machine, using the same clamping systems. The tests were performed under tensile, compressive and alternating loads with stress ratios of $R={\sigma }^{\text{min}}/{\sigma }^{\text{max}}=0.1$ , 10 and –1, respectively. Test frequencies between $0.11\mathit{mu}\text{Hz}$ and 5Hz were employed, depending on the applied fatigue load level and the resulting cyclic strains. The specimens were tested in the low cycle fatigue regime until failure or till a cycle count of 100,000 cycles was reached. Experiments reaching this cycle count were considered as run-outs. The maximum number of loading cycles is derived from the requirements for composite gaseous fuel tanks, featuring up to 60,000 filling cycles during their lifetime. The cyclic experiments were performed under force control. During the experiments, the force and displacement peaks as well as the cycle count were recorded. Selected specimens were additionally equipped with strain gauges in order to determine the stiffness degradation during the fatigue process. All experiments were performed at ambient temperature.

Results

The results of the quasi-static experiments are presented in Figure 5. The top three subfigures are related to tensile experiments whereas the results of the experiments under compressive loads are presented in the bottom three subfigures. In the left hand column of subfigures, the results of the experiments within the fiber direction are presented whereas the center and right hand columns of subfigures are related to experiments on unidirectionally fiber reinforced material perpendicular to the fiber direction and on the multidirectionally reinforced laminates respectively. In most cases, a brittle response is obtained where an approximately linear stress–strain curve ends with a sudden failure of the specimen without distinct plastic deformation. The only case with a distinctively nonlinear stress–strain curve prior to failure is observed in the case of UD material tested within the fiber direction for specimen no. ZU1-1.2-01B-02. Nevertheless, even for the case of the UD material tested under compression perpendicular to the fiber direction, all other specimens exhibit a visible, but limited amount of nonlinearity prior to final failure. For all fiber orientations, similar elastic stiffnesses are obtained in tension and compression.

Figure 5.

Quasi-static stress–strain curves.

In Figure 6, the results of the fatigue experiments are presented. Again, the three columns of subfigures are related to unidirectionally fiber reinforced material tested within and perpendicular to the fiber direction and multidirectional laminates tested in the ${\hspace{0.17em}}^0\circ$ -direction, respectively. The top line of subfigures is related to tensile cyclic experiments with a stress ratio of R = 0.1, whereas the center and bottom lines of subfigures are related to fatigue experiments under compressive loads with R = 10 and alternating cyclic loads with $R=-1$ . The fatigue results are evaluated in terms of S-N-curves where the 50% quantile is determined as the linear regression line to the experimental results in the double-logarithmic representations whereas the 5%- and 95%-quantiles are determined as parallel lines to the 50%-quantile under the assumption of a Gaussian distribution which is equal for all load levels. The run-outs as well as the static results—transformed into pseudo-amplitudes according to the respective stress (R-) ratio—are added for information. In all cases, S-N-curves with small slopes are obtained. The absolute scatter of the results with respect to the stress amplitude ( ${\sigma }^{\text{a}}$ -) direction is limited. For the experiments perpendicular to the fiber direction (center row of subfigures), higher relative amounts of scatter are obtained. Nevertheless, considering the small absolute stress amplitudes ${\sigma }^{\text{a}}$ for these experiments, the scatter for these cases can also be regarded as limited.

Figure 6.

S-N-curves under cyclic loads.

The determination of the stiffness degradation from the stress and strain envelopes during the fatigue loading tends to be difficile due to heavy white noise development. Further problems were encountered due to surface damage developing during the fatigue process, which increasingly disturbed the tactile strain measurement using clip gauges on the surfaces. The results are presented in Figure 7 where the normalized stiffness $E/E_0$ is plotted as a function of the cycle count N normalized with the number $N^{\text{f}}$ of cycles to failure of the respective specimen. No unique picture is obtained for the reference material investigated. Whereas a distinct initial degradation is obtained for both experiments in cyclic tension with valid strain envelope determination (specimen no. ZU1-1-1-01A-02 experienced an early failure of the local strain measurement due to surface damage), only one out of four specimens tested in cyclic compression within the fiber direction exhibits such an effect. No initial degradation is observed for the specimens tested in tension or compression perpendicular to the fiber direction. For the multidirectional laminates tested in cyclic tension, more distinct stiffness degradation effects are observed due to interaction of the plies featuring different fiber orientations. For the multidirectional laminates tested under cyclic compression, no valid local strain measurements were obtained.

Figure 7.

Degradation under cyclic loads.

Numerical example

Parameter identification

As an exemplary application, the fatigue damage model defined in Section “Damage degradation model” is applied to the experimental data base described in Section “Experimental data base.” In a first step, the elastic moduli are determined as the secant moduli to the quasi-static stress strain curves between 0.05% and 0.25% nominal stain (Figure 5). Assuming initial transverse isotropy, the x₂- and x₃-directions are considered to be equivalent. Subsequently, the material parameters in the damage evolution equations (20) to (22) with (23) to (25) are determined by means of a reverse engineering approach. For this purpose, finite element computations are performed, where a single element is subjected to a cyclic unidirectional stress state with a sinusoidal characteristic. For acceleration of the fatigue degradation, a cycle jump technique with an acceleration factor of 1000 is used (i.e. a single computational cycle represents 1000 real loading cycles, see Section “Implicit time integration”). For reasons of computational efficiency, the explicit time integration scheme is employed. Nevertheless, competitive analyses based on both schemes show only negligible differences. The finite element analysis is continued until the element fails due to damage accumulation. In this context, the element is considered as failed, once any of the damage variables D₁, D₂, or D₃ exceeds a level of 90%. This approximation definition is required since the theoretical “kink” in the stress–strain response when reaching a damage level of 100% had to be smoothed out mathematically in order to ensure numerical stability.

In the reverse engineering approach, the damage accumulation parameters are determined. Considering two stress levels within the relevant fatigue range in Figure 6, the respective parameters for the x₁- and x₂-direction are varied from guessed initial values, until S-N-curves fitting with the experimental data are obtained. Both, the tensile and compressive cyclic modes are considered. In this context, it is found that a common value for all exponents n _ij in the damage evolution equations (20) to (22) provides a reasonable approximation. Hence, in order to reduce the number of parameters to a necessary minimum, all exponents n _ij are assumed to be equal, keeping the scaling parameters A _ij as the only different parameters for the different loading modes.

The identified parameter set is listed in Table 1. Regarding the initial elastic material parameters, transverse isotropy is assumed. However, it has to be pointed out that the transverse isotropy with respect to the stiffness in general is lost as soon as anisotropic damage with $D_2\ne D_3$ develops due to nonhydrostatic loading conditions.

Table 1.

Material parameters.

E 1	E 2	G 12	G 23	ν 12	ν 23
$(\text{GPa})$	$(\text{GPa})$	$(\text{GPa})$	$(\text{GPa})$	$(-)$	$(-)$
168.6	10.80	5.140	4.012	0.276	0.346
$a_{\text{I}}$	$D_{\text{I}}$	$a_{\text{II}}$	$D_{\text{II}}$	b	n
$(-)$	$(-)$	$(-)$	$(-)$	$(-)$	$(-)$
5.0	0.1	10.0	0.5	0.5	10
$A_{11}^{\text{t}}$	$A_{11}^{\text{c}}$	A 12	A 21	$A_{22}^{\text{t}}$	$A_{22}^{\text{c}}$	A 23
$(\text{MP}{\text{a}}^{-n+1})$	$(\text{MP}{\text{a}}^{-n+1})$	$(\text{MP}{\text{a}}^{-n+1})$	$(\text{MP}{\text{a}}^{-n+1})$	$(\text{MP}{\text{a}}^{-n+1})$	$(\text{MP}{\text{a}}^{-n+1})$	$(\text{MP}{\text{a}}^{-n+1})$
$3.25·{10}^{-40}$	$1.10·{10}^{-34}$	$0.50·{10}^{-8}$	$7.75·{10}^{-17}$	$2.75·{10}^{-17}$	$9.25·{10}^{-20}$	0

The resulting S-N-curves are presented in Figure 8. In all four cases, the numerical results for the determined data set are found in good agreement with the experimental data. Regarding the shape of the numerical S-N-curves, a linear material response in the double logarithmic representation is obtained (notice that the linear lines presented in Figure 8 are based on multiple data rather than by their end points only). The linear response is due to the assumption of the power-law definition (20) to (22) of the damage evolution equations. The fact that the well known linear shape of the S-N-curve in the double logarithmic Wöhler diagrams is obtained in the numerical simulations underlines the usefulness of the definition of the damage evolution in the present model.

Figure 8.

Parameter determination.

Validation

For a further validation and an evaluation of the transferability of the model and the parameter set presented in Table 1, the model with the previously determined parameter set is applied to the fatigue experiments from Section “Experimental data base” with alternating loads as well as to the experiments on multidirectional laminates.

The results for application of the model to the experimental data base on unidirectionally fiber reinforced specimens with alternating cyclic loads featuring a stress ratio of $R={\sigma }_{\text{min}}/{\sigma }_{\text{max}}=-1$ are presented in Figure 9. Again, a linear shape is obtained for the numerical S-N-curves in the double logarithmic Wöhler diagrams for both, the experiments within ( ${\hspace{0.17em}}^0\circ$ ) and perpendicular ( ${\hspace{0.17em}}^{90}\circ$ ) to the fiber direction. For the cyclic load applied within the fiber direction ( ${\hspace{0.17em}}^0\circ$ ), a good agreement with the experimental data is observed; however, for a limited number of data. Nevertheless, especially for the data at higher numbers $N^{\text{f}}$ of cycles to failure the quantitative agreement is good. For the experiments with cyclic loads applied perpendicular to the fiber direction ( ${\hspace{0.17em}}^{90}\circ$ ), the numerical results are found at the lower edge of the experimental scatter band.

Figure 9.

Validation–fatigue under alternating loads.

Multidirectional laminates recombining plies with different fiber orientation imply the development of more complex stress states in the individual plies even when loaded in uniaxial modes. Due to the interaction of the individual plies and the deformation constraints imposed to individual plies by their neighbors with alternative fiber and thus stiffness orientations, the stress states in all plies will be of the bi- or tri-axial type, irrespectively of the external loading direction. In order to investigate the response of the proposed damage model under more complex cyclic stress states, the model is applied to the experiments from Section “Experimental data base” on multidirectional laminates.

The results for stress ratios of R = 0.1 (tensile cyclic loads), 10 (compressive cyclic loads), and −1 (alternating cyclic loads) are presented in Figure 10. Again, perfectly linear S-N-curves in the double logarithmic representation of the Wöhler diagrams are obtained. For tensile cyclic loads, the experimentally observed fatigue strength is overestimated whereas an underestimate of the experimentally determined fatigue strength by the numerical results is observed for compressive cyclic loads. Nevertheless, in both cases, the deviation of experimental and numerical data is found in the same order of magnitude for typical scatter band widths obtained in the experimental investigation, especially when data sets consisting of a larger number of data are considered. For alternating cyclic loads with R = –1, a qualitative and quantitative agreement of experimental data and numerical results based on the proposed fiber reinforced fatigue damage model is obtained.

Figure 10

Validation–fatigue of multidirectional laminates.

Conclusions

Objective of the present study is the definition of a continuum damage mechanics model for fatigue and degradation of fiber reinforced plastics. In order to keep the model and formulation as simple as possible, a brittle formulation based on the anisotropic linear elastic Hooke's law is employed. The model is enhanced by introduction of three independent damage variables, representing damage effects within three orthogonal planes. For the damage evolution, a power-law stress formulation depending linearly on the stress increment and in power-law form on the actual stress level is defined. This definition ensures that linear S-N-curves are obtained in the double logarithmic Wöhler fatigue strength diagram.

The model is implemented as a user defined material model into a commercial finite element system. Both, implicit and explicit time integrations are considered. In a comparison of the performance and accuracy, the explicit forward Euler time integration scheme is found to provide an efficient and accurate solution, since damage effects provide the only nonlinearity in the model and the damage accumulation in individual time steps in general is small. As an alternative to the cycle by cycle type of damage accumulation during the finite element analysis, a damage accumulation rule in the Miner sense is formulated, based on the proposed continuum damage model. This damage accumulation rule allows for the determination of damage accumulation during the postprocessing of stress analyses without consideration of damage effects in a numerically efficient manner.

The model is applied to and validated against an experimental data base on a filament wound CFRP epoxy matrix composite. The model proves to provide a reasonable qualitative approximation of the experimental data for all cases of unidirectionally reinforced specimens as well as multidirectional laminates, considering tensile, compressive, and alternating cyclic loading conditions. Since the model accounts for the damage effects on the local material stiffness during the structural analysis, the material model allows a damage mechanics analysis of fiber reinforced structures in order to assess their fatigue degradation and successive failure.