- Email: [email protected]
Progressive damage modeling of composite materials subjected to mixed mode cyclic loading using cohesive zone model
Journal Pre-proof
Progressive damage modelling of composite materials subjected to mixed mode cyclic loading using cohesive zone model J. Ebadi-Rajoli , A. Akhavan-Safar , H. Hosseini-Toudeshky , LFM. da Silva PII: DOI: Reference:
S0167-6636(19)30900-7 https://doi.org/10.1016/j.mechmat.2020.103322 MECMAT 103322
To appear in:
Mechanics of Materials
Received date: Revised date: Accepted date:
16 October 2019 2 January 2020 13 January 2020
Please cite this article as: J. Ebadi-Rajoli , A. Akhavan-Safar , H. Hosseini-Toudeshky , LFM. da Silva , Progressive damage modelling of composite materials subjected to mixed mode cyclic loading using cohesive zone model, Mechanics of Materials (2020), doi: https://doi.org/10.1016/j.mechmat.2020.103322
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.
Progressive damage modelling of composite materials subjected to mixed mode cyclic loading using cohesive zone model
J. Ebadi-Rajolia, A. Akhavan-Safarb, H. Hosseini-Toudeshkyc, LFM. da Silvad a
Department of Mechanical Engineering, Amirkabir University of Technology, Hafez Avenue, Tehran, Iran. b Institute of Science and Innovation in Mechanical and Industrial Engineering (INEGI), Rua Dr. Roberto Frias, 4200-465 Porto, Portugal. c Department of Aerospace Engineering, Amirkabir University of Technology, Hafez Avenue, 424, Tehran, Iran. d Department of Mechanical Engineering, Faculty of Engineering of the University of Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal.
Abstract Cohesive zone modeling (CZM) has been extensively considered as a powerful method for analysis of initiation and propagation of delamination in composite laminates subjected to cyclic loading. By making a relation between the damage parameter and the loading cycles, damage accumulation can be calculated, and the fatigue life of a component can be estimated. The aim of the current research is to present a new fatigue damage accumulation (FDA) model. The proposed model is developed based on the concepts of CZM and by linking the fracture mechanics to damage mechanics approaches. For this purpose, a developed user material subroutine (USDFLD) is implemented in a finite element software to simulate the initiation and propagation of damage in an interface layer of a composite material subjected to cyclic loading. The method is validated for different loading conditions using Paris law results. Good agreement between the proposed method and the Paris law results is observed. Keywords: Fatigue, Composite, CZM, Delamination. 1
Nomenclature list:
: initial stiffness
Interface element properties
: cohesive nominal thickness
: normal stress
: fracture toughness
: in plane shear stress
,
: out of plane shear stress
,
: normal stress at damage initiation ,
: normal strain and separation ,
: in plane shear strain and separation
,
: out of plane shear strain and separation
: total energy release rate ,
: mode-I, II energy release rate : maximum energy release rate : minimum energy release rate
: in plane shear strain at damage initiation
: effective strain at damage initiation : effective strain at final damage (failure) : strain ratio ,
: Benzeggagh–Kenane exponents
: ratio of static damage increment to fatigue damage increment : cohesive zone area
: mode I, II threshold energy release rate, respectively
: crack growth rate
: normal strain at damage initiation : maximum effective strain in loading history
: mode I, II fracture toughness
: energy release rate amplitude : load ratio : Paris law constant ,
: mode-I,II Paris law constant : mixed-mode Paris law constant
: Paris law exponent ,
: mode-I, II Paris law exponent : mixed-mode Paris law exponent
Bulk Material Properties:
: fatigue damage
E11, E22: Young’s modulus along direction 1, 2 respectively
: static damage
G12, G13, G23: shear modulus
: damage variable
12: Poisson’s ratio List of acronyms: Cohesive zone modeling (CZM)
Fatigue crack growth (FCG)
linear elastic fracture mechanics (LEFM)
Double cantilever beams (DCB)
Experimental Data Fitting (EDF)
End notched flexure (ENF)
Linking Damage mechanics and Fracture mechanics (LDF)
Mixed mode bending (MMB)
2
1. Introduction Composite materials have been increasingly used in advanced industries where lightweight structures are considered. However, predicting the behavior of composite materials under complex loading conditions is still challenging. One of the main type of service load for such materials is cyclic loading. Different models/tools have been developed to predict the fatigue behavior of composites. However, CZM is more considered by researchers. CZM is suitable for damage analysis of composite materials and also is able to predict both fatigue initiation and propagation of composites by considering the degradation of the material properties due to cyclic loading. In the framework of finite element method, there are two approaches to simulate the damage accumulation at an interface subjected to fatigue loading (Harper and Hallett 2010): (i) By applying the actual loading conditions to the structure and degrading the stiffness of the material based on the cycle-by-cycle concept. One can derive that the interface element constitutive equation is history dependent (Roe and Siegmund 2003) which is useful for low cycle fatigue analysis of ductile materials. (ii) By utilizing a loading envelope strategy where the applied load remains constant at the maximum value of the cyclic load, and the interface element is degraded based on the amount of damage accumulated at each time increment. This method is useful for high cycle fatigue analysis of linear conditions with linear elastic fracture mechanics (LEFM) assumptions. This method reduces the computational cost and time tremendously compared to the former approach (Khoramishad et al. 2010; Peerlings et al. 2000; Robinson et al. 2005; Turon et al. 2007; Van Paepegem and Degrieck 2001). Damage accumulation methods using a loading envelope approach can be categorized into two different sub-models: 3
1) Damage accumulation based on Experimental Data Fitting (EDF): In some EDF models, damage propagation is presented as a function of the effective strain (Khoramishad et al. 2010; Peerlings et al. 2000; Robinson et al. 2005) or effective stress (Van Paepegem and Degrieck 2001). The main disadvantage of the EDF based methods is the trial and error step where the model parameters should be obtained. 2) Damage accumulation by Linking Damage mechanics and Fracture mechanics (LDF): This method was first presented by (Turon et al. 2007). LDF model parameters can be directly obtained using fatigue crack growth (FCG) tests of double cantilever beams (DCB), end notched flexure specimens (ENF) and mixed mode bending (MMB) samples. Some of the proposed approaches such as (Turon et al. 2007) are based on analytical solutions in order to calculate the length of the cohesive zone formed ahead of the crack tip. However, they don’t consider the effects of static damage on fatigue damage accumulation. Setting the cycle jump,
, is also another restriction of the mentioned approach. Pirondi and Moroni (2010)
presented a fatigue damage formulation based on the stiffness-method. Harper and Hallett (2010) considered the LDF method and the effect of static damage on fatigue damage growth. In their model, the energy release rate is calculated locally but for each analysis the cohesive zone length (Lcz) should be extracted manually at the last static load increment. Kawashita and Hallett (2012) proposed a fatigue damage model for delamination analysis of composite materials subjected to high-cycle fatigue loading using cohesive elements. The major advantage of their proposed formulation is its complete independency of the cohesive zone length. Recently, an approach was presented by some authors to take into account the degradation effects of cyclic loading on cohesive properties of adhesive materials (Costa et al. 2018; Monteiro et al. 2019). It
4
was applied for pure mode I and pure mode II loading conditions. It doesn’t need the damage zone size for predicting the degradation of cohesive properties. In this paper, a new model is proposed for fatigue analysis of composite materials based on the LDF method. The proposed model is implemented in a finite element software using the USDFLD subroutine. The proposed approach is validated for different mode mixity conditions using the Paris law model.
2. Cohesive element constitutive equations Various CZM shapes have been considered by authors (Cui and Wisnom 1993; Mi et al. 1998; Needleman 1987; Tvergaard and Hutchinson 1992). In this paper, a bilinear traction separation law is employed. Eq. 1 presents the constitutive equations of the considered bilinear cohesive element. { where
(
} ,
) { and
}
{
〈
〉
};
〈 〉
{
(1)
shown in Figure 1 are the normal, in plane shear and out of plane shear
stresses of the cohesive elements, respectively. K in Eq. 1 is the initial stiffness of the elements. Based on the finite element concepts, the strain of cohesive elements can be obtained by dividing the separation (displacement) by the thickness of the cohesive element as follows: (2) where
is the strain and
is the separation at the interface. Subscripts n, sn and tn are related to
the normal, in-plane and out-of-plane modes, respectively.
in Eq. 2 is the thickness of the
cohesive element. Damage parameter (d) in Eq. 1 is defined by Eq. 3. 5
Figure 1. Three-dimensional stress state of a solid-like interface element (Balzani and Wagner 2008)
(
) (
where
(3)
)
is the maximum effective strain (
) in loading history which in 2D problems can
be obtained by Eq. 4 (compressive loads don’t introduce any damage into the cohesive elements). √〈 〉
(4)
in Eq. 3 should be calculated using Eq. 5 (Ye 1988):
√ (
)
(
)
(5)
{ where: (6)
is also found using the Benzeggagh–Kenane (BK) failure criterion based on Eq. 7 (Benzeggagh and Kenane 1996): 6
{
(
*
),
- +
(7)
3. Fatigue damage model 3.1.
Fatigue crack growth
Paris law with a normalized strain energy (Eq. 8) was used for fatigue crack propagation analysis.
{ (
where
)
(8)
is variation of the energy release rate which should be defined as follows:
(
where and
and
)
√
(9)
are the threshold energy release rate and the fracture toughness, respectively.
can be found using BK failure criterion based on Eq. 10 (Benzeggagh and Kenane
1996): (
)(
) (10)
(
)(
)
7
where
and
are the mode I and mode II fracture energies, respectively.
the mode I and mode II threshold energy release rates, respectively. rate which is equal to
and
are
is the total energy release
.
However, based on some reported experimental data, it has been found that Eq. 8 is not necessarily valid for all loading conditions and materials (Monteiro et al. 2018). The effect of mode-mixity on m and C can be considered using Blanco et al. equations (Blanco et al. 2004):
(
)
(
)
(
) (11)
(
)
(
) (
)
where I, II and m are representatives of mode I, mode II and mixed mode conditions, respectively. The mode ratio (
)) should be obtained using the following relation. (12)
where
3.2.
is the ratio of the shear strain to the normal strain.
The proposed fatigue model
In this section a fatigue damage model is presented which considers the effect of the updated static damage on the total damage of the cohesive elements subjected to cyclic loading. Assuming that
is the initial stiffness of the cohesive elements, the stiffness of the -th element
ahead of the crack tip in the cohesive damage zone will be degraded to
(
(
) )
8
before the ith fatigue time increment, where
is the damage growth caused by fatigue loading.
This degradation approach was applied to the material by developing a USDFLD subroutine. As the applied load in the fatigue step is constant (based on the cycle jumping strategy), the damage level will be constant so, the life will be infinite. To simulate the degradation of the material properties due to fatigue loading and, consequently, to simulate the damage growth, the level of damage value should be increased. By increasing the damage level by the value of
,
ABAQUS solver should resolve the problem based on the degraded elements. Part of this degradation is applied to the initial stiffness of the elements. As the stiffness for the damaged cohesive elements decreases by factor of (1-d), the displacement of the elements due to the fatigue loads (load is constant during the fatigue analysis) increases to satisfy the static equilibrium conditions. This increase in displacement will cause a new damage which is referred to as the static damage ( ) in this work. We call it static damage as it refers to the static equilibrium equations. This phenomenon is shown in Figure 2 and the procedure of obtaining the damage values is given in Figure 4. By increasing the amount of damage, the
will increase
during the fatigue step.
9
Figure 2. Schematic view of the damage parameter as a function of distance ahead of the crack tip Accordingly, the total damage in cycle
is equal to the damage in cycle
plus the static
and fatigue damage which is mathematically shown in Eq. 13. (13)
where
is the static damage growth due to the redistribution of the stiffness of the interface
elements ahead of the crack tip. In damage mechanics an increase in damaged area ( equal to an increase in cracked area (
(
)
∑
(
)
) is
) (see Figure 3):
(14)
10
where ncz is the number of integration points in the cohesive zone area (where 0
is the
effective area of -th integration point and subscript is related to the ith increment.
Figure 3- Schematic representation of the equivalence between the increase of the damaged area and the crack growth Since the total damage is the summation of the static and fatigue damage, therefore: ( ∑
where
((
)
(
) )
15)
is the size of crack growth at ith time increment.
As the damage locus is equal to the cohesive zone ahead of the crack tip (
∑
),
accordingly, Eq. 15 can be rewritten as follows:
11
(
) ((̅̅̅̅̅̅̅̅)
̅̅̅̅̅̅̅ ( ))
(16)
( ) are the average damage increments for the fatigue and static parts, ( ) and ̅̅̅̅̅̅̅ where ̅̅̅̅̅̅̅̅ respectively. All these steps were implemented in the FE analysis by developing a USDFLD subroutine. In this study, the static damage is calculated by the FEM software and the fatigue damage should be obtained by the USDFLD subroutine. According to the proposed approach, the mathematical form of the above explanation can be presented as follows: (
̅̅̅̅̅̅̅̅ ( ) (
) ̅̅̅̅̅̅̅ ( ) ) ̅̅̅̅̅̅̅̅ ( )
) (
Based on Eq. 17, (
(
(
)
(
(17)
) for each integration point can be calculated approximately as below:
)
(18)
)(
)
is obtained from the (i-1)th increment as follows:
where (
)
(
)
Accordingly, (
)
(
)
( (
) )
(19)
for the first increment of fatigue step should be calculated as below:
(20)
12
where (
) is the difference of damage values between the last two increments of the static
step. The only remaining unknown parameter in Eq. 20 is the ( equal to (
)
at all integration points. By setting the value of
) , which is considered to be , the value of (
) can be
determined for the following time increments. The progressive fatigue delamination analysis flowchart is shown in Fig. 4.
Figure 4. Progressive fatigue delamination analysis flowchart For the USDFLD subroutine, the input parameters are ,
,
, ) and FCG coefficients ( ,
,
,
,
, cohesive zone properties ( , ,
,
). The effective strain (
are calculated from the finite element solution, then the mixed mode parameters (
,
) and
,
and
) are calculated based on the initiation and evolution criteria. Finally, the value of (
) is
calculated at each integration point using Eq. 18. It should be noted that the mentioned static 13
damage is related to the numerical procedure and not to the physical aspect of the fatigue test. Both damage parts (the static and dynamic damages) don’t have a physical meaning.
4. Results and discussions The proposed fatigue damage model was applied to DCB, ENF and MMB specimens for pure mode I, pure mode II and mixed mode (I/II) loading conditions, respectively. 50% mixed mode conditions was considered where
. The geometries and the properties of the specimens
are already reported in literature. The specimen geometries are shown in Fig. 5 and test configuration for mode I, mode II and 50% mixed mode conditions are shown in Fig. 6. Material and cohesive properties for HTA/6376C carbon/epoxy used for the numerical models are given in Table 1. Table 2 shows the fatigue properties of the cohesive elements by considering Eqs. 10 and 12. Table 1. Material and cohesive properties for HTA/6376C carbon/epoxy (Robinson et al. 2005) Composite E11 (GPa) E22 (GPa) G12=G13 (GPa) G23 (GPa) ν12
120 10.5 5.25 3.48 0.3
Cohesive layer GIc (KJ/m2) GIIc (KJ/m2) η σn0 (MPa) τ0 (MPa) K (MPa)
0.26 1.002 2 15 15 105
Table 2. Fatigue delamination properties for HTA/6376C carbon/epoxy (Turon et al. 2007) CI (mm2/cycle) CII (mm2/cycle) C50% (mm2/cycle) Cm (mm2/cycle)
0.0616 2.99 4.23 458087
mI mII m50% mm
5.4 4.5 6.41 4.94
GIth (KJ/m2) 0.060 GIIth (KJ/m2) 0.100 η2 2.73
14
Figure 5. Specimen geometry.
(a) Mode I: DCB test
(b) Mode II: 3 point ENF test (Reeder et al. 2004)
(c) Mixed-mode specimen. Figure 6. Loading pattern for Mode I, Mode II and Mixed-mode conditions Fatigue delamination was simulated under load-control condition. Simulation was performed using ABAQUS package. Loading configuration and boundary condition for each case are shown in Fig. 6. For computational efficiency in mode-II configuration support rollers are replaced by the boundary conditions shown in Fig. 7. This is a common method in Mode-II
15
simulation (Harper and Hallett 2010; Turon et al. 2007). A frictionless contact between the primary crack surfaces is defined to prevent the penetration of the upper and lower crack surfaces.
Figure 7. Loading and boundary conditions applied to mode-II sample
The load was applied in two static steps in ABAQUS package. In the first step the load increases linearly from zero to the maximum value and in the second step the maximum load remains constant. The fatigue simulation is performed at the load ratio of zero (R=0) based on the previously described loading envelope strategy. Damage only occurs in cohesive elements and the kinematic behavior of all materials (strain-displacement) is simulated as linear. Specimen (beam) meets plane strain condition when r<
)
(21)
16
where
is the applied maximum load,
(Reeder et al. 2004),
is the crack length,
and ⁄
is crack length correction term
, B is the specimen width, and h is the
specimen half-depth (see Fig. 5). For pure mode II condition, G is calculated as follows (Reeder et al. 2004): (
)
(22) 0.42μh.
where
For the 50% mixed mode, G is obtained as follows (Turon et al. 2007):
(
√
(23)
)
The fatigue crack growth rate can be calculated from the following equation:
(
)
∑
∑ ∑
(24)
Fig. 8 shows the state of the damage (d) along the bondline at the end of the static step and for different loading cycles during fatigue loading.
17
Figure 8. Damage distribution along the bondline at different steps of the numerical analysis Figures 9-11 show the variation of fatigue crack growth rate as a function of the normalized energy in mode I, mode II and 50% mixed mode conditions, respectively. It is shown in Figures 9-11 that the results of the proposed model are in good agreement with the Paris-law curve data. This is due to using an analytical relation for G, considering the actual values of Lcz and also taking into account the effect of static damage at each time increment.
18
Figure 9. Variation of fatigue crack growth rate versus
⁄
in mode I condition.
=15,
K=105
Figure 10. Variation of fatigue crack growth rate versus
⁄
in mode II condition.
=15,
K=105 19
Figure 11. Variation of fatigue crack growth rate as a function of condition.
⁄
for mixed mode I/II
=15, K=105
5. Conclusion In this study, a new approach for numerical fatigue analysis of composite materials was proposed by considering the effects of static damage growth during the cyclic loading. It was shown that the amount of static damage is not constant during fatigue loading. By a precise estimation of the strain energy release rate at each time increment during the fatigue analysis, the proposed fatigue analysis approach is able to accurately predict the fatigue crack growth in composite materials. The presented extrapolation method can be used for other damage and failure modes to determine the fatigue damage accumulation with continuum damage mechanics.
20
CRediT author statement J. Ebadi-Rajoli: Conceptualization, Methodology, Software, Validation, Writing - Original Draft, Writing - Review & Editing, Project administration A. Akhavan-Safar: Writing - Review & Editing, Validation. H. Hosseini-Toudeshky: Writing - Review & Editing, Project administration, Validation. LFM. da Silva: Writing - Review & Editing, Validation.
Declaration of competing interests none
References Anderson, TL., 2005. Fracture mechanics: fundamentals and applications. CRC press. Balzani, C., Wagner, W., 2008. An interface element for the simulation of delamination in unidirectional fiber-reinforced composite laminates. Engineering Fracture Mechanics 75(9):2597-2615. doi:10.1016/j.engfracmech.2007.03.013. Benzeggagh, M., Kenane, M., 1996. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Composites science and technology 56(4):439-449. Blanco, N., Gamstedt, EK., Asp, L., Costa, J., 2004. Mixed-mode delamination growth in carbon–fibre composite laminates under cyclic loading. International Journal of Solids and Structures 41(15):4219-4235. Costa, M., Viana, G., Créac’hcadec, R., da Silva, L., Campilho, R., 2018. A cohesive zone element for mode I modelling of adhesives degraded by humidity and fatigue. International Journal of Fatigue 112:173-182. Cui, W., Wisnom, M., 1993. A combined stress-based and fracture-mechanics-based model for predicting delamination in composites. Composites 24(6):467-474. Harper, PW., Hallett, SR., 2010. A fatigue degradation law for cohesive interface elements – Development and application to composite materials. International Journal of Fatigue 32(11):1774-1787 doi:10.1016/j.ijfatigue.2010.04.006. Kawashita, LF., Hallett, SR., 2012. A crack tip tracking algorithm for cohesive interface element analysis of fatigue delamination propagation in composite materials. International Journal of Solids and Structures 49(21):2898-2913 doi:10.1016/j.ijsolstr.2012.03.034. Khoramishad, H., Crocombe, AD., Katnam, KB., Ashcroft, IA., 2010. Predicting fatigue damage in adhesively bonded joints using a cohesive zone model. International Journal of Fatigue 32(7):1146-1158 doi:10.1016/j.ijfatigue.2009.12.013. 21
Mi, Y., Crisfield, M., Davies, G., Hellweg, H., 1998. Progressive delamination using interface elements. Journal of Composite Materials 32(14):1246-1272. Monteiro, J., et al. 2018. Influence of mode mixity and loading conditions on the fatigue crack growth behavior of an epoxy adhesive. Theoretical and Applied Fracture Mechanics Accepted. Monteiro, J., et al. 2019. Mode II modeling of adhesive materials degraded by fatigue loading using cohesive zone elements. Theoretical and Applied Fracture Mechanics 103:102253. Needleman, A., 1987. A continuum model for void nucleation by inclusion debonding. Journal of applied mechanics 54(3):525-531. Peerlings, R., Brekelmans, W., De Borst, R., Geers, M., 2000. Gradient-enhanced damage modelling of high-cycle fatigue. International Journal for Numerical Methods in Engineering 49(12):1547-1569. Pirondi, A., Moroni, F., 2010. A Progressive Damage Model for the Prediction of Fatigue Crack Growth in Bonded Joints. The Journal of Adhesion 86(5-6):501-521 doi:10.1080/00218464.2010.484305. Reeder, JR., Demarco, K., Whitley, KS., 2004. The use of doubler reinforcement in delamination toughness testing. Composites Part A: Applied Science and Manufacturing 35(11):13371344. Robinson, P., Galvanetto, U., Tumino, D., Bellucci, G., Violeau, D., 2005. Numerical simulation of fatigue-driven delamination using interface elements. International Journal for Numerical Methods in Engineering 63(13):1824-1848 doi:10.1002/nme.1338. Roe, K., Siegmund, T., 2003. An irreversible cohesive zone model for interface fatigue crack growth simulation. Engineering Fracture Mechanics 70(2):209-232. Turon, A., Costa, J., Camanho, PP., Dávila, CG., 2007. Simulation of delamination in composites under high-cycle fatigue. Composites Part A: Applied Science and Manufacturing 38(11):2270-2282 doi:10.1016/j.compositesa.2006.11.009. Tvergaard, V., Hutchinson, JW., 1992. The relation between crack growth resistance and fracture process parameters in elastic-plastic solids. Journal of the Mechanics and Physics of Solids 40(6):1377-1397. Van Paepegem, W., Degrieck, J., 2001. Fatigue degradation modelling of plain woven glass/epoxy composites. Composites Part A: Applied Science and Manufacturing. 32(10):1433-1441 doi:10.1016/S1359-835X(01)00042-2. Ye, L., 1988. Role of matrix resin in delamination onset and growth in composite laminates. Composites science and technology 33(4):257-277.
22
Progressive damage modelling of composite materials subjected to mixed mode cyclic loading using cohesive zone model J. Ebadi-Rajoli , A. Akhavan-Safar , H. Hosseini-Toudeshky , LFM. da Silva PII: DOI: Reference:
S0167-6636(19)30900-7 https://doi.org/10.1016/j.mechmat.2020.103322 MECMAT 103322
To appear in:
Mechanics of Materials
Received date: Revised date: Accepted date:
16 October 2019 2 January 2020 13 January 2020
Please cite this article as: J. Ebadi-Rajoli , A. Akhavan-Safar , H. Hosseini-Toudeshky , LFM. da Silva , Progressive damage modelling of composite materials subjected to mixed mode cyclic loading using cohesive zone model, Mechanics of Materials (2020), doi: https://doi.org/10.1016/j.mechmat.2020.103322
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.
Progressive damage modelling of composite materials subjected to mixed mode cyclic loading using cohesive zone model
J. Ebadi-Rajolia, A. Akhavan-Safarb, H. Hosseini-Toudeshkyc, LFM. da Silvad a
Department of Mechanical Engineering, Amirkabir University of Technology, Hafez Avenue, Tehran, Iran. b Institute of Science and Innovation in Mechanical and Industrial Engineering (INEGI), Rua Dr. Roberto Frias, 4200-465 Porto, Portugal. c Department of Aerospace Engineering, Amirkabir University of Technology, Hafez Avenue, 424, Tehran, Iran. d Department of Mechanical Engineering, Faculty of Engineering of the University of Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal.
Abstract Cohesive zone modeling (CZM) has been extensively considered as a powerful method for analysis of initiation and propagation of delamination in composite laminates subjected to cyclic loading. By making a relation between the damage parameter and the loading cycles, damage accumulation can be calculated, and the fatigue life of a component can be estimated. The aim of the current research is to present a new fatigue damage accumulation (FDA) model. The proposed model is developed based on the concepts of CZM and by linking the fracture mechanics to damage mechanics approaches. For this purpose, a developed user material subroutine (USDFLD) is implemented in a finite element software to simulate the initiation and propagation of damage in an interface layer of a composite material subjected to cyclic loading. The method is validated for different loading conditions using Paris law results. Good agreement between the proposed method and the Paris law results is observed. Keywords: Fatigue, Composite, CZM, Delamination. 1
Nomenclature list:
: initial stiffness
Interface element properties
: cohesive nominal thickness
: normal stress
: fracture toughness
: in plane shear stress
,
: out of plane shear stress
,
: normal stress at damage initiation ,
: normal strain and separation ,
: in plane shear strain and separation
,
: out of plane shear strain and separation
: total energy release rate ,
: mode-I, II energy release rate : maximum energy release rate : minimum energy release rate
: in plane shear strain at damage initiation
: effective strain at damage initiation : effective strain at final damage (failure) : strain ratio ,
: Benzeggagh–Kenane exponents
: ratio of static damage increment to fatigue damage increment : cohesive zone area
: mode I, II threshold energy release rate, respectively
: crack growth rate
: normal strain at damage initiation : maximum effective strain in loading history
: mode I, II fracture toughness
: energy release rate amplitude : load ratio : Paris law constant ,
: mode-I,II Paris law constant : mixed-mode Paris law constant
: Paris law exponent ,
: mode-I, II Paris law exponent : mixed-mode Paris law exponent
Bulk Material Properties:
: fatigue damage
E11, E22: Young’s modulus along direction 1, 2 respectively
: static damage
G12, G13, G23: shear modulus
: damage variable
12: Poisson’s ratio List of acronyms: Cohesive zone modeling (CZM)
Fatigue crack growth (FCG)
linear elastic fracture mechanics (LEFM)
Double cantilever beams (DCB)
Experimental Data Fitting (EDF)
End notched flexure (ENF)
Linking Damage mechanics and Fracture mechanics (LDF)
Mixed mode bending (MMB)
2
1. Introduction Composite materials have been increasingly used in advanced industries where lightweight structures are considered. However, predicting the behavior of composite materials under complex loading conditions is still challenging. One of the main type of service load for such materials is cyclic loading. Different models/tools have been developed to predict the fatigue behavior of composites. However, CZM is more considered by researchers. CZM is suitable for damage analysis of composite materials and also is able to predict both fatigue initiation and propagation of composites by considering the degradation of the material properties due to cyclic loading. In the framework of finite element method, there are two approaches to simulate the damage accumulation at an interface subjected to fatigue loading (Harper and Hallett 2010): (i) By applying the actual loading conditions to the structure and degrading the stiffness of the material based on the cycle-by-cycle concept. One can derive that the interface element constitutive equation is history dependent (Roe and Siegmund 2003) which is useful for low cycle fatigue analysis of ductile materials. (ii) By utilizing a loading envelope strategy where the applied load remains constant at the maximum value of the cyclic load, and the interface element is degraded based on the amount of damage accumulated at each time increment. This method is useful for high cycle fatigue analysis of linear conditions with linear elastic fracture mechanics (LEFM) assumptions. This method reduces the computational cost and time tremendously compared to the former approach (Khoramishad et al. 2010; Peerlings et al. 2000; Robinson et al. 2005; Turon et al. 2007; Van Paepegem and Degrieck 2001). Damage accumulation methods using a loading envelope approach can be categorized into two different sub-models: 3
1) Damage accumulation based on Experimental Data Fitting (EDF): In some EDF models, damage propagation is presented as a function of the effective strain (Khoramishad et al. 2010; Peerlings et al. 2000; Robinson et al. 2005) or effective stress (Van Paepegem and Degrieck 2001). The main disadvantage of the EDF based methods is the trial and error step where the model parameters should be obtained. 2) Damage accumulation by Linking Damage mechanics and Fracture mechanics (LDF): This method was first presented by (Turon et al. 2007). LDF model parameters can be directly obtained using fatigue crack growth (FCG) tests of double cantilever beams (DCB), end notched flexure specimens (ENF) and mixed mode bending (MMB) samples. Some of the proposed approaches such as (Turon et al. 2007) are based on analytical solutions in order to calculate the length of the cohesive zone formed ahead of the crack tip. However, they don’t consider the effects of static damage on fatigue damage accumulation. Setting the cycle jump,
, is also another restriction of the mentioned approach. Pirondi and Moroni (2010)
presented a fatigue damage formulation based on the stiffness-method. Harper and Hallett (2010) considered the LDF method and the effect of static damage on fatigue damage growth. In their model, the energy release rate is calculated locally but for each analysis the cohesive zone length (Lcz) should be extracted manually at the last static load increment. Kawashita and Hallett (2012) proposed a fatigue damage model for delamination analysis of composite materials subjected to high-cycle fatigue loading using cohesive elements. The major advantage of their proposed formulation is its complete independency of the cohesive zone length. Recently, an approach was presented by some authors to take into account the degradation effects of cyclic loading on cohesive properties of adhesive materials (Costa et al. 2018; Monteiro et al. 2019). It
4
was applied for pure mode I and pure mode II loading conditions. It doesn’t need the damage zone size for predicting the degradation of cohesive properties. In this paper, a new model is proposed for fatigue analysis of composite materials based on the LDF method. The proposed model is implemented in a finite element software using the USDFLD subroutine. The proposed approach is validated for different mode mixity conditions using the Paris law model.
2. Cohesive element constitutive equations Various CZM shapes have been considered by authors (Cui and Wisnom 1993; Mi et al. 1998; Needleman 1987; Tvergaard and Hutchinson 1992). In this paper, a bilinear traction separation law is employed. Eq. 1 presents the constitutive equations of the considered bilinear cohesive element. { where
(
} ,
) { and
}
{
〈
〉
};
〈 〉
{
(1)
shown in Figure 1 are the normal, in plane shear and out of plane shear
stresses of the cohesive elements, respectively. K in Eq. 1 is the initial stiffness of the elements. Based on the finite element concepts, the strain of cohesive elements can be obtained by dividing the separation (displacement) by the thickness of the cohesive element as follows: (2) where
is the strain and
is the separation at the interface. Subscripts n, sn and tn are related to
the normal, in-plane and out-of-plane modes, respectively.
in Eq. 2 is the thickness of the
cohesive element. Damage parameter (d) in Eq. 1 is defined by Eq. 3. 5
Figure 1. Three-dimensional stress state of a solid-like interface element (Balzani and Wagner 2008)
(
) (
where
(3)
)
is the maximum effective strain (
) in loading history which in 2D problems can
be obtained by Eq. 4 (compressive loads don’t introduce any damage into the cohesive elements). √〈 〉
(4)
in Eq. 3 should be calculated using Eq. 5 (Ye 1988):
√ (
)
(
)
(5)
{ where: (6)
is also found using the Benzeggagh–Kenane (BK) failure criterion based on Eq. 7 (Benzeggagh and Kenane 1996): 6
{
(
*
),
- +
(7)
3. Fatigue damage model 3.1.
Fatigue crack growth
Paris law with a normalized strain energy (Eq. 8) was used for fatigue crack propagation analysis.
{ (
where
)
(8)
is variation of the energy release rate which should be defined as follows:
(
where and
and
)
√
(9)
are the threshold energy release rate and the fracture toughness, respectively.
can be found using BK failure criterion based on Eq. 10 (Benzeggagh and Kenane
1996): (
)(
) (10)
(
)(
)
7
where
and
are the mode I and mode II fracture energies, respectively.
the mode I and mode II threshold energy release rates, respectively. rate which is equal to
and
are
is the total energy release
.
However, based on some reported experimental data, it has been found that Eq. 8 is not necessarily valid for all loading conditions and materials (Monteiro et al. 2018). The effect of mode-mixity on m and C can be considered using Blanco et al. equations (Blanco et al. 2004):
(
)
(
)
(
) (11)
(
)
(
) (
)
where I, II and m are representatives of mode I, mode II and mixed mode conditions, respectively. The mode ratio (
)) should be obtained using the following relation. (12)
where
3.2.
is the ratio of the shear strain to the normal strain.
The proposed fatigue model
In this section a fatigue damage model is presented which considers the effect of the updated static damage on the total damage of the cohesive elements subjected to cyclic loading. Assuming that
is the initial stiffness of the cohesive elements, the stiffness of the -th element
ahead of the crack tip in the cohesive damage zone will be degraded to
(
(
) )
8
before the ith fatigue time increment, where
is the damage growth caused by fatigue loading.
This degradation approach was applied to the material by developing a USDFLD subroutine. As the applied load in the fatigue step is constant (based on the cycle jumping strategy), the damage level will be constant so, the life will be infinite. To simulate the degradation of the material properties due to fatigue loading and, consequently, to simulate the damage growth, the level of damage value should be increased. By increasing the damage level by the value of
,
ABAQUS solver should resolve the problem based on the degraded elements. Part of this degradation is applied to the initial stiffness of the elements. As the stiffness for the damaged cohesive elements decreases by factor of (1-d), the displacement of the elements due to the fatigue loads (load is constant during the fatigue analysis) increases to satisfy the static equilibrium conditions. This increase in displacement will cause a new damage which is referred to as the static damage ( ) in this work. We call it static damage as it refers to the static equilibrium equations. This phenomenon is shown in Figure 2 and the procedure of obtaining the damage values is given in Figure 4. By increasing the amount of damage, the
will increase
during the fatigue step.
9
Figure 2. Schematic view of the damage parameter as a function of distance ahead of the crack tip Accordingly, the total damage in cycle
is equal to the damage in cycle
plus the static
and fatigue damage which is mathematically shown in Eq. 13. (13)
where
is the static damage growth due to the redistribution of the stiffness of the interface
elements ahead of the crack tip. In damage mechanics an increase in damaged area ( equal to an increase in cracked area (
(
)
∑
(
)
) is
) (see Figure 3):
(14)
10
where ncz is the number of integration points in the cohesive zone area (where 0
is the
effective area of -th integration point and subscript is related to the ith increment.
Figure 3- Schematic representation of the equivalence between the increase of the damaged area and the crack growth Since the total damage is the summation of the static and fatigue damage, therefore: ( ∑
where
((
)
(
) )
15)
is the size of crack growth at ith time increment.
As the damage locus is equal to the cohesive zone ahead of the crack tip (
∑
),
accordingly, Eq. 15 can be rewritten as follows:
11
(
) ((̅̅̅̅̅̅̅̅)
̅̅̅̅̅̅̅ ( ))
(16)
( ) are the average damage increments for the fatigue and static parts, ( ) and ̅̅̅̅̅̅̅ where ̅̅̅̅̅̅̅̅ respectively. All these steps were implemented in the FE analysis by developing a USDFLD subroutine. In this study, the static damage is calculated by the FEM software and the fatigue damage should be obtained by the USDFLD subroutine. According to the proposed approach, the mathematical form of the above explanation can be presented as follows: (
̅̅̅̅̅̅̅̅ ( ) (
) ̅̅̅̅̅̅̅ ( ) ) ̅̅̅̅̅̅̅̅ ( )
) (
Based on Eq. 17, (
(
(
)
(
(17)
) for each integration point can be calculated approximately as below:
)
(18)
)(
)
is obtained from the (i-1)th increment as follows:
where (
)
(
)
Accordingly, (
)
(
)
( (
) )
(19)
for the first increment of fatigue step should be calculated as below:
(20)
12
where (
) is the difference of damage values between the last two increments of the static
step. The only remaining unknown parameter in Eq. 20 is the ( equal to (
)
at all integration points. By setting the value of
) , which is considered to be , the value of (
) can be
determined for the following time increments. The progressive fatigue delamination analysis flowchart is shown in Fig. 4.
Figure 4. Progressive fatigue delamination analysis flowchart For the USDFLD subroutine, the input parameters are ,
,
, ) and FCG coefficients ( ,
,
,
,
, cohesive zone properties ( , ,
,
). The effective strain (
are calculated from the finite element solution, then the mixed mode parameters (
,
) and
,
and
) are calculated based on the initiation and evolution criteria. Finally, the value of (
) is
calculated at each integration point using Eq. 18. It should be noted that the mentioned static 13
damage is related to the numerical procedure and not to the physical aspect of the fatigue test. Both damage parts (the static and dynamic damages) don’t have a physical meaning.
4. Results and discussions The proposed fatigue damage model was applied to DCB, ENF and MMB specimens for pure mode I, pure mode II and mixed mode (I/II) loading conditions, respectively. 50% mixed mode conditions was considered where
. The geometries and the properties of the specimens
are already reported in literature. The specimen geometries are shown in Fig. 5 and test configuration for mode I, mode II and 50% mixed mode conditions are shown in Fig. 6. Material and cohesive properties for HTA/6376C carbon/epoxy used for the numerical models are given in Table 1. Table 2 shows the fatigue properties of the cohesive elements by considering Eqs. 10 and 12. Table 1. Material and cohesive properties for HTA/6376C carbon/epoxy (Robinson et al. 2005) Composite E11 (GPa) E22 (GPa) G12=G13 (GPa) G23 (GPa) ν12
120 10.5 5.25 3.48 0.3
Cohesive layer GIc (KJ/m2) GIIc (KJ/m2) η σn0 (MPa) τ0 (MPa) K (MPa)
0.26 1.002 2 15 15 105
Table 2. Fatigue delamination properties for HTA/6376C carbon/epoxy (Turon et al. 2007) CI (mm2/cycle) CII (mm2/cycle) C50% (mm2/cycle) Cm (mm2/cycle)
0.0616 2.99 4.23 458087
mI mII m50% mm
5.4 4.5 6.41 4.94
GIth (KJ/m2) 0.060 GIIth (KJ/m2) 0.100 η2 2.73
14
Figure 5. Specimen geometry.
(a) Mode I: DCB test
(b) Mode II: 3 point ENF test (Reeder et al. 2004)
(c) Mixed-mode specimen. Figure 6. Loading pattern for Mode I, Mode II and Mixed-mode conditions Fatigue delamination was simulated under load-control condition. Simulation was performed using ABAQUS package. Loading configuration and boundary condition for each case are shown in Fig. 6. For computational efficiency in mode-II configuration support rollers are replaced by the boundary conditions shown in Fig. 7. This is a common method in Mode-II
15
simulation (Harper and Hallett 2010; Turon et al. 2007). A frictionless contact between the primary crack surfaces is defined to prevent the penetration of the upper and lower crack surfaces.
Figure 7. Loading and boundary conditions applied to mode-II sample
The load was applied in two static steps in ABAQUS package. In the first step the load increases linearly from zero to the maximum value and in the second step the maximum load remains constant. The fatigue simulation is performed at the load ratio of zero (R=0) based on the previously described loading envelope strategy. Damage only occurs in cohesive elements and the kinematic behavior of all materials (strain-displacement) is simulated as linear. Specimen (beam) meets plane strain condition when r<
)
(21)
16
where
is the applied maximum load,
(Reeder et al. 2004),
is the crack length,
and ⁄
is crack length correction term
, B is the specimen width, and h is the
specimen half-depth (see Fig. 5). For pure mode II condition, G is calculated as follows (Reeder et al. 2004): (
)
(22) 0.42μh.
where
For the 50% mixed mode, G is obtained as follows (Turon et al. 2007):
(
√
(23)
)
The fatigue crack growth rate can be calculated from the following equation:
(
)
∑
∑ ∑
(24)
Fig. 8 shows the state of the damage (d) along the bondline at the end of the static step and for different loading cycles during fatigue loading.
17
Figure 8. Damage distribution along the bondline at different steps of the numerical analysis Figures 9-11 show the variation of fatigue crack growth rate as a function of the normalized energy in mode I, mode II and 50% mixed mode conditions, respectively. It is shown in Figures 9-11 that the results of the proposed model are in good agreement with the Paris-law curve data. This is due to using an analytical relation for G, considering the actual values of Lcz and also taking into account the effect of static damage at each time increment.
18
Figure 9. Variation of fatigue crack growth rate versus
⁄
in mode I condition.
=15,
K=105
Figure 10. Variation of fatigue crack growth rate versus
⁄
in mode II condition.
=15,
K=105 19
Figure 11. Variation of fatigue crack growth rate as a function of condition.
⁄
for mixed mode I/II
=15, K=105
5. Conclusion In this study, a new approach for numerical fatigue analysis of composite materials was proposed by considering the effects of static damage growth during the cyclic loading. It was shown that the amount of static damage is not constant during fatigue loading. By a precise estimation of the strain energy release rate at each time increment during the fatigue analysis, the proposed fatigue analysis approach is able to accurately predict the fatigue crack growth in composite materials. The presented extrapolation method can be used for other damage and failure modes to determine the fatigue damage accumulation with continuum damage mechanics.
20
CRediT author statement J. Ebadi-Rajoli: Conceptualization, Methodology, Software, Validation, Writing - Original Draft, Writing - Review & Editing, Project administration A. Akhavan-Safar: Writing - Review & Editing, Validation. H. Hosseini-Toudeshky: Writing - Review & Editing, Project administration, Validation. LFM. da Silva: Writing - Review & Editing, Validation.
Declaration of competing interests none
References Anderson, TL., 2005. Fracture mechanics: fundamentals and applications. CRC press. Balzani, C., Wagner, W., 2008. An interface element for the simulation of delamination in unidirectional fiber-reinforced composite laminates. Engineering Fracture Mechanics 75(9):2597-2615. doi:10.1016/j.engfracmech.2007.03.013. Benzeggagh, M., Kenane, M., 1996. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Composites science and technology 56(4):439-449. Blanco, N., Gamstedt, EK., Asp, L., Costa, J., 2004. Mixed-mode delamination growth in carbon–fibre composite laminates under cyclic loading. International Journal of Solids and Structures 41(15):4219-4235. Costa, M., Viana, G., Créac’hcadec, R., da Silva, L., Campilho, R., 2018. A cohesive zone element for mode I modelling of adhesives degraded by humidity and fatigue. International Journal of Fatigue 112:173-182. Cui, W., Wisnom, M., 1993. A combined stress-based and fracture-mechanics-based model for predicting delamination in composites. Composites 24(6):467-474. Harper, PW., Hallett, SR., 2010. A fatigue degradation law for cohesive interface elements – Development and application to composite materials. International Journal of Fatigue 32(11):1774-1787 doi:10.1016/j.ijfatigue.2010.04.006. Kawashita, LF., Hallett, SR., 2012. A crack tip tracking algorithm for cohesive interface element analysis of fatigue delamination propagation in composite materials. International Journal of Solids and Structures 49(21):2898-2913 doi:10.1016/j.ijsolstr.2012.03.034. Khoramishad, H., Crocombe, AD., Katnam, KB., Ashcroft, IA., 2010. Predicting fatigue damage in adhesively bonded joints using a cohesive zone model. International Journal of Fatigue 32(7):1146-1158 doi:10.1016/j.ijfatigue.2009.12.013. 21
Mi, Y., Crisfield, M., Davies, G., Hellweg, H., 1998. Progressive delamination using interface elements. Journal of Composite Materials 32(14):1246-1272. Monteiro, J., et al. 2018. Influence of mode mixity and loading conditions on the fatigue crack growth behavior of an epoxy adhesive. Theoretical and Applied Fracture Mechanics Accepted. Monteiro, J., et al. 2019. Mode II modeling of adhesive materials degraded by fatigue loading using cohesive zone elements. Theoretical and Applied Fracture Mechanics 103:102253. Needleman, A., 1987. A continuum model for void nucleation by inclusion debonding. Journal of applied mechanics 54(3):525-531. Peerlings, R., Brekelmans, W., De Borst, R., Geers, M., 2000. Gradient-enhanced damage modelling of high-cycle fatigue. International Journal for Numerical Methods in Engineering 49(12):1547-1569. Pirondi, A., Moroni, F., 2010. A Progressive Damage Model for the Prediction of Fatigue Crack Growth in Bonded Joints. The Journal of Adhesion 86(5-6):501-521 doi:10.1080/00218464.2010.484305. Reeder, JR., Demarco, K., Whitley, KS., 2004. The use of doubler reinforcement in delamination toughness testing. Composites Part A: Applied Science and Manufacturing 35(11):13371344. Robinson, P., Galvanetto, U., Tumino, D., Bellucci, G., Violeau, D., 2005. Numerical simulation of fatigue-driven delamination using interface elements. International Journal for Numerical Methods in Engineering 63(13):1824-1848 doi:10.1002/nme.1338. Roe, K., Siegmund, T., 2003. An irreversible cohesive zone model for interface fatigue crack growth simulation. Engineering Fracture Mechanics 70(2):209-232. Turon, A., Costa, J., Camanho, PP., Dávila, CG., 2007. Simulation of delamination in composites under high-cycle fatigue. Composites Part A: Applied Science and Manufacturing 38(11):2270-2282 doi:10.1016/j.compositesa.2006.11.009. Tvergaard, V., Hutchinson, JW., 1992. The relation between crack growth resistance and fracture process parameters in elastic-plastic solids. Journal of the Mechanics and Physics of Solids 40(6):1377-1397. Van Paepegem, W., Degrieck, J., 2001. Fatigue degradation modelling of plain woven glass/epoxy composites. Composites Part A: Applied Science and Manufacturing. 32(10):1433-1441 doi:10.1016/S1359-835X(01)00042-2. Ye, L., 1988. Role of matrix resin in delamination onset and growth in composite laminates. Composites science and technology 33(4):257-277.
22