Fatigue delamination damage analysis in composite materials through a rule of mixtures approach

Contents lists available at ScienceDirect Composite Structures journal homepage: www.elsevier.com/locate/compstruct Fatigue delamination damage analysis in composite materials through a rule of mixtures approach Alireza Taherzadeh-Farda,b,∗, Sergio Jiméneza,b, Alejandro Cornejoa,b, Eugenio Oñatea,b, Lucia Gratiela Barbua,b aPolytechnic University of Catalonia (UPC), Campus Nord, 08034 Barcelona, Spain bInternational Center for Numerical Methods in Engineering (CIMNE), Campus Nord UPC, 08034 Barcelona, Spain ARTICLE INFO Keywords: Rule of mixtures Homogenization Fatigue Delamination Composites Finite element method ABSTRACT The present study investigates delamination damage initiation and propagation within a homogenization theory of mixtures, using the concept of virtual layers and virtual interfaces. It eliminates spatial discretization of layers, introducing a resultant damage variable to capture structure’s bulk response under both monotonic and cyclic loads. Fatigue-induced deterioration is classified into sub-critical, critical, and over-critical stages based on interfacial stresses. Calibration is conducted employing the widely-available Wöhler curves for each loading mode independently. An advance-in-time strategy is included in the model to enhance the simulation speed. The reliability of the approach is assessed for crack initiation and propagation separately through standard test coupons, showing good correlation with experimental data in mode I, mode II, and mixedmode loading conditions. Depending on the calibration procedure adopted, the model is applicable to a wide range of stress ratios. In addition, it could be integrated into any standard finite element framework using the desired number of elements through the thickness regardless of the physical amount of layers. This allows easy modification of stacking sequences or the number of layers within the constitutive law without mesh structure changes, facilitating simulation of large-scale composite laminates with minimal accuracy loss and reduced computational costs. 1. Introduction Over the past recent years, composites have attracted much attention as the potential materials to be utilized in different industries including aerospace, automotive, marine, etc. [1,2], due to their superior mechanical properties such as corrosion and fatigue resistance, light weight, high strength and ease of transportation [3–6]. The multimaterial nature of composite materials gives more freedom to engineers for designing purposes, however, this would lead to more damage mechanisms being engaged once these materials are under structural loading as well as more uncertainty in their bulk response due to manufacturing complexities [7–9]. Therefore, more advanced methodologies are required in studying the behavior of composites. Damage mechanisms in composites are mainly comprised of intralaminar and inter-laminar modes [10]. According to the complex nature of these damages and their interactions within the structure, simulation schemes have encountered challenges to compromise between accuracy and computational cost [11]. Among all degradation modes, interface failure or delamination has been of significant importance for ∗Corresponding author at: Polytechnic University of Catalonia (UPC), Campus Nord, 08034 Barcelona, Spain. E-mail addresses: [email protected],[email protected] (A. Taherzadeh-Fard). its major effect on the strength and stiffness of the composite [12– 15]. Indeed, laminated composites are fabricated through building up the in-plane fiber-reinforced laminae and are mainly expected to bear loading in the fiber direction. However, general loading conditions usually impose some through-the-thickness stresses to the composite, and can lead to failure of the matrix and initiation of the delamination phenomenon [16,17]. Debonding in the adjacent layers may not be catastrophic within the in-plane loading, nonetheless, for composites under out-of-plane loading or bending conditions, it could be fatal [18, 19]. Prediction of this initiation and consequent propagation of the delamination crack has been associated with experimental full-scale testing procedures owing to the lack of enough knowledge and accurate numerical models, which, in turn, has made the design procedures to be highly expensive in terms of time and cost [20]. Fatigue delamination damage is mainly composed of three distinct stages of initiation, onset and propagation, and each of them has been characterized by different numerical approaches [21]. The initiation phase refers to the formation of the crack from an intact material and commonly originates from the imperfections formed during the https://doi.org/10.1016/j.compstruct.2024.118613 Received 11 June 2024; Received in revised form 8 August 2024; Accepted 23 September 2024 Composite Structures 351 (2025) 118613 Available online 27 September 2024 0263-8223/© 2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). A. Taherzadeh-Fard et al. manufacturing process [11]. A common method to account for this stage is defining an initiation zone where damage could start appearing according to the models based on S-N curves [22,23]. In this way, SN curves are formulated using experimentally-defined parameters and each mode is calibrated independently. Delamination onset, however, is defined as a visible increment in the crack through the medium from a crack starter. Although there is an unstable delamination crack propagation after crack initiation, crack onset is followed by a stable condition of crack extention [24]. Crack onset could be analyzed in a similar manner as initiation by considering fitting parameters to the experimental data [25], however this method lacks the ability to consider the mixed-mode and stress ratios. Linear elastic fracture mechanics (LEFM) is a widely-used approach in simulating the fatigue crack in the propagation phase [26,27]. These models mainly employ any forms of the Paris law to make a connection between the crack growth and the energy release rate (ERR) [11]. Therefore, other methods should also be incorporated within this formulation to calculate the energy release rate or the stress intensity factor such as Virtual Crack Closure Technique (VCCT) [8,28– 30], virtual crack extension [31,32], or crack surface displacement extrapolation methods [33]. Despite this, LEFM based models suffer from several problems such as the inability to be used in hybrid composites [34] or failure to initiate the fatigue crack within a pristine material [24]. Cyclic Cohesive Zone Method (CCZM) can also be used to simulate fatigue crack in composites which is based on defining a tractionseparation law for damage evolution within the interface [23,35–39]. Based on how the cyclic load effects are taken into account, two different methods of hysteresis loop and envelope load damage models have been considered [40]. Hysteresis loop damage models degrade the strength and the stiffness within a cycle-by-cycle algorithm and are well suited for low-cycle fatigue regime [41]. Nevertheless, calibration of the parameters is considered to be a highly challenging task [40]. Envelope load damage models, on the other side, utilize jumping strategies combined with the Paris law to evolve delamination damage within the interface. This formulation could be beneficial in the high-cycle fatigue regime and is capable of considering different factors such as mode mix and stress ratios. However, implementing Paris law within the cohesive framework can be complex [42], and linking the damage variable to the traction-separation formulations is not straightforward [20]. In addition, it is not clear if the Paris law should be considered as a material property or rather a problem dependent variable. Consequently, fatigue delamination damage models that are not based on the Paris law are highly valuable. Moreover, cohesive zone methods have several shortcomings. Mainly since they are dependent on the definition of an interface medium, crack path should be known a priori [43]. Additionally, incorporating interface elements within all interfaces would be a challenging task especially in composites with large number of layers. Since CCZMs are dependent on the definition of a traction-separation law with an initial linear part, they could manipulate the bulk response of the structure in its elastic region when there is not damage evolved within the interface [44]. Selection of the slope in CCZM linear part is also ad-hoc and requires a tuning process [45,46]. Recently, a continuum mechanics based fatigue model has been proposed by unifying different phenomena such as damage, plasticity, viscosity and temperature effects [47], and has been improved in other applications as well [48,49]. Indeed, the formulation triggers strength and stiffness deterioration once a specific stress criterion is satisfied. The degradation is accounted for by a damage-type constitutive law, while being sensitive to the cyclic loads through a fatigue reduction factor [50]. The factor also affects the fracture energy so that an energy dissipation is introduced to the model as a result of the cyclic loading. Model parameters are calibrated by using curve fitting to the common S-N curves. The model is expected to be applicable in a wide spectrum of fatigue scenarios while being capable of capturing both damage initiation and propagation over different amount of stress ratios. In the current study, a homogenization theory of mixtures has been presented to model fatigue delamination damage response in laminated composites. The work is based on a previous research by the authors [51] where loading conditions were limited to the monotonic regime. Rather than considering cohesive elements at each interface, the proposed model employs a homogenized medium where there is no need to explicitly discretize composite layers, which could significantly facilitate the pre-processing stage as well as the optimization of the design process, since no re-meshing is required when the internal configuration of the composite changes. This is even more pronounced when modeling composites with large number of layers by accounting the number of elements needed through the thickness irrespective of the actual number of layers. The cyclic loading effect is introduced to the formulation by defining a fatigue reduction function which could be calibrated according to the common S-N curves for mode I and mode II, independently. Delamination damage effect, in both contexts of crack initiation and propagation, is reflected to the adjacent virtual bulk layers so that the bulk response of the structure is accurately reproduced, while not imposing any limitations on the constitutive law used for each layer. Therefore, in addition to the inter-layer damages, there is the possibility of intera-layer mechanical phenomena such as damage and plasticity. Since the ad-hoc interface penalty stiffness concept is eliminated from the formulation, structural response will no longer be manipulated in the elastic regime. The proposed model is compatible with any finite element code and is considered to be wellsuited for simulating fatigue delamination in large-scale composites while maintaining the accuracy within an acceptable range. 2. Mathematical formulation 2.1. Preliminary concepts In this section, some introductory information regarding the underlying rule of mixtures method and interfacial stresses is provided. The reader is referred to the previous work [51] for a complete description and fundamental hypotheses. Composite laminates are mainly composed of a specific number of fiber reinforced laminae oriented at different directions. Trusdell and Toupin [52] firstly presented a method to unify all layers’ behavior through a homogenization theory. Later on, the method has been developed by other researchers as well [53,54]. It is mainly based on the assumption that strains experienced by all laminae are the same whilst the stress of the composite is obtained by taking a super-position of the stress in each laminae according to its volumetric participation, i.e.: 𝜺1=𝜺2=⋯=𝜺𝑛(1) 𝝈= 𝑛 ∑ 𝑖=1 (𝑘𝑖⋅𝝈𝑖)(2) where 𝜺𝑖,𝑘𝑖and 𝝈𝑖are the strain tensor, volumetric participation and the stress tensor of the 𝑖th layer, respectively. From the homogenization point of view, the composite is made up of virtual bulk layers and virtual interfaces, as presented in Fig. 1, in a way that the response of the structure is reproduced correctly in the macro scale. The first step in formulating the delamination is to prepare a foundation to calculate interfacial stresses. Since free edge effects are not accounted for in the homogenized medium, a straightforward averaging scheme is employed to obtain stresses at the interfaces, as follows: 𝝈𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑖= 𝝈𝑏𝑢𝑙𝑘 𝑖+𝝈𝑏𝑢𝑙𝑘 𝑖−1 2(3) where 𝝈𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑖is the interfacial stress tensor at the 𝑖th interface while 𝝈𝑏𝑢𝑙𝑘 𝑖is the bulk layer stress tensor at the 𝑖th layer. Composite Structures 351 (2025) 118613 2 A. Taherzadeh-Fard et al. Fig. 1. Virtual bulk layers and virtual interfaces through the homogenized medium [51]. According to Balzani et al. [55], delamination is driven by a normal stress component acting perpendicular to the delamination plane and two shear stress components acting on the delamination plane. Considering the coordinate system presented in Fig. 1, two equivalent stresses are calculated to be fed into the delamination damage model for mode I and mode II loadings: 𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑛,𝑖 =⟨𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑧𝑧,𝑖 ⟩(4) 𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑠,𝑖 =√(𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑥𝑧,𝑖 )2+ (𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑦𝑧,𝑖 )2(5) where 𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑛,𝑖 and 𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑠,𝑖 are the interfacial equivalent normal and shear stresses at the 𝑖th interface. ⟨∙⟩are the Macaulay brackets to prevent damage evolution in normal compression loads [51]. 2.2. Fatigue delamination damage evolution In order to accumulate the fatigue effect into the current delamination formulation, a similar scheme to that described in [49] is adopted here. To this end, the failure indicators 𝐹𝑛,𝑖 and 𝐹𝑠,𝑖 previously developed at each interface 𝑖in [51] are modified by a fatigue reduction function (𝑓𝑟𝑒𝑑 ) as below: 𝐹𝑛,𝑖 = 𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑛,𝑖 𝑓𝑟𝑒𝑑 𝑛,𝑖 (𝑁𝑐 𝑛,𝑖, 𝑅𝑛,𝑖, 𝑆𝑚𝑎𝑥 𝑛,𝑖 )−𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑛,𝑡ℎ,𝑖 >0(6) 𝐹𝑠,𝑖 = 𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑠,𝑖 𝑓𝑟𝑒𝑑 𝑠,𝑖 (𝑁𝑐 𝑠,𝑖, 𝑅𝑠,𝑖, 𝑆𝑚𝑎𝑥 𝑠,𝑖 )−𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑠,𝑡ℎ,𝑖 >0(7) where subscriptions n and s refer to normal mode I and shear mode II loadings, respectively. 𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑛,𝑡ℎ,𝑖 and 𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑠,𝑡ℎ,𝑖 are historical variables of the interface threshold stress in mode I and mode II, and their initial value would consider to be normal and shear interface strengths, respectively. These values will be updated at each step to the maximum historical normal (𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑛,𝑖 ) and shear (𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑠,𝑖 ) uni-axial stresses obtained previously. Once inequalities (6) and (7) are satisfied, the crack is permitted to grow. 𝑓𝑟𝑒𝑑 (𝑁𝑐 𝑚,𝑖, 𝑅𝑚,𝑖, 𝑆𝑚𝑎𝑥 𝑚,𝑖 )is introduced to the failure indicators to take the cyclic load effects into account by amplifying the stress state. Its value ranges from 1, indicating no cyclic load effect, down to 0 at the asymptote. According to the fatigue crack propagation through the interface, which is shown in Fig. 2, three different stages could be considered at each integration point, whilst each stage requires specific formulation for 𝑓𝑟𝑒𝑑 , as described in the following. 2.2.1. Sub-critical fatigue stage Once fatigue crack initiates and starts propagating, there are some areas in the far-field where the maximum cyclic stress at the interface (𝑆𝑚𝑎𝑥) is below the fatigue limit, 𝑆𝑡ℎ, at a specific stress ratio, 𝑅= Fig. 2. Different stages experienced by a band of elements in (a) Sub-critical (blue), (b) Critical (gray), and (c) Over-critical (red) regions as fatigue crack propagates in the direction illustrated by arrows. 𝑆𝑚𝑖𝑛∕𝑆𝑚𝑎𝑥, matching the blue area in Fig. 2(a). An element located in this area (black dashed strip) should not experience any cyclic load effects, and hence, the fatigue related quantities will be maintained in their initial values, as below: 𝑓𝑟𝑒𝑑 𝑚,𝑖 = 1 (8) 𝑁𝑐 𝑚,𝑖 = 0 (9) where 𝑁𝑐 𝑚,𝑖 is the number of cycles at the 𝑖th interface and 𝑚stands for 𝑛or 𝑠in normal and shear loadings, respectively. As the crack Composite Structures 351 (2025) 118613 3 A. Taherzadeh-Fard et al. propagates, the black dashed area becomes closer to the crack tip and the maximum cyclic stress levels will increase, which necessitates a different description for the 𝑓𝑟𝑒𝑑 . 2.2.2. Critical fatigue stage When a specific element experiences stresses greater than the fatigue limit and below the static strength, the cyclic load effects should be considered through 𝑓𝑟𝑒𝑑 . This situation is observed by the elements within the shaded area in the gray region in Fig. 2(b), where variation of 𝑓𝑟𝑒𝑑 could trigger two kinds of non-linearity. The first one is the gradual amplifying of the interfacial stresses (or literally, the gradual decrease of the strength) as a result of the 𝑓𝑟𝑒𝑑 evolution when 𝑁𝑐is growing. The second one, which is a consequence of the first one, is the rapid propagation of the crack once 𝑁𝑐reaches the critical number of cycles (𝑁𝑓) and inequalities (6) and (7) are met. The definition of the 𝑓𝑟𝑒𝑑 function can be done based on either GN (energy release rate vs. number of cycles) or S-N (Interfacial stress vs. number of cycles) diagrams. Indeed, the model in this study has originated from the high-cycle fatigue (HCF) formulation presented by Oller et al. [47], which has been developed for metals based on SN curves. Successful implementation and satisfactory results of this formulation in HCF regime [49,56] triggered the idea of applying the same approach in composite materials [57]. Therefore, although calibrating the model based on G-N curves seems to be more applicable as for more availability of these curves for composite materials, present formulation has been established based on the stress curves. Supporting this approach, there are recently-published papers where S-N curves have been employed, and a calibration process has been conducted accordingly to model the delamination propagation in laminated composites [58,59]. This may indicate the potential of the S-N curves in this context. Moreover, according to [42], some inter-connections could be established between the common Paris law (da/dN vs. G) and the S-N curve in composites. This would be an asset for the current approach in this study, since there is the potential to feed the Paris law to the model and make calibrations accordingly. This issue has not been addressed in current study, and would be a matter of future investigations. Definition of 𝑓𝑟𝑒𝑑 for the interface is conducted through a fitting procedure with respect to the common S-N surfaces of the type [49]: 𝑆𝑚,𝑖(𝑅𝑚,𝑖, 𝑁𝑐 𝑚,𝑖) =𝑆𝑡ℎ 𝑚,𝑖(𝑅𝑚,𝑖) + (𝑆𝑢 𝑚,𝑖 −𝑆𝑡ℎ 𝑚,𝑖(𝑅𝑚,𝑖))⋅exp{−𝛼𝑡 𝑚,𝑖(𝑅𝑚,𝑖) ⋅(log 𝑁𝑐 𝑚,𝑖)𝛽𝑓 𝑚,𝑖 }(10) 𝑅𝑚,𝑖 and 𝑆𝑢 𝑚,𝑖 being the stress ratio and static strength for mode 𝑚at the 𝑖th interface, respectively. 𝑆𝑡ℎ 𝑚,𝑖(𝑅𝑚,𝑖)is the fatigue limit at the 𝑖th interface: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ if |𝑅|≤1⇒𝑆𝑡ℎ 𝑚,𝑖(𝑅𝑚,𝑖) = 𝑆𝑒 𝑚,𝑖 +(𝑆𝑢 𝑚,𝑖 −𝑆𝑒 𝑚,𝑖)⋅(1 + 𝑅𝑚,𝑖 2)𝑆𝑅1 𝑚,𝑖 if |𝑅|>1⇒𝑆𝑡ℎ 𝑚,𝑖(𝑅𝑚,𝑖) = 𝑆𝑒 𝑚,𝑖 +(𝑆𝑢 𝑚,𝑖 −𝑆𝑒 𝑚,𝑖)⋅(1 + 𝑅𝑚,𝑖 2𝑅𝑚,𝑖 )𝑆𝑅2 𝑚,𝑖 (11) where 𝑆𝑒 𝑚,𝑖 is the fatigue limit at 𝑅𝑚,𝑖 = −1.𝛼𝑡 𝑚,𝑖(𝑅𝑚,𝑖)is a parameter which depends on the state of loading as described below: ⎧ ⎪ ⎨ ⎪ ⎩ if |𝑅|≤1⇒𝛼𝑡 𝑚,𝑖(𝑅𝑚,𝑖) = 𝛼𝑓 𝑚,𝑖 +(1 + 𝑅𝑚,𝑖 2)⋅𝐴𝑈𝑋𝑅1 𝑚,𝑖 if |𝑅|>1⇒𝛼𝑡 𝑚,𝑖(𝑅𝑚,𝑖) = 𝛼𝑓 𝑚,𝑖 −(1 + 𝑅𝑚,𝑖 2𝑅𝑚,𝑖 )⋅𝐴𝑈𝑋𝑅2 𝑚,𝑖 (12) 𝑆𝑅1 𝑚,𝑖 ,𝑆𝑅2 𝑚,𝑖 ,𝛼𝑓 𝑚,𝑖,𝛽𝑓 𝑚,𝑖,𝐴𝑈𝑋𝑅1 𝑚,𝑖 and 𝐴𝑈𝑋𝑅2 𝑚,𝑖 are material properties and could be calibrated at each interface 𝑖for each mode of loading 𝑚. The more Wöhler diagrams in different stress ratios are considered in the calibration, the more accurate description of the S-N surface will be obtained in a wide spectrum of stress ratios. Once S-N surface has been defined and formulated, an expression could be proposed for fatigue reduction function in the form 𝑓𝑟𝑒𝑑 𝑚,𝑖 (𝑁𝑐 𝑚,𝑖, 𝑅𝑚,𝑖, 𝑆𝑚𝑎𝑥 𝑚,𝑖 ) = exp ⎡⎢⎢⎢⎢⎣ ln (𝑆𝑚𝑎𝑥 𝑚,𝑖 ∕𝑆𝑚,𝑢,𝑖) (log 𝑁𝑓 𝑚,𝑖)(𝛽𝑓 𝑚,𝑖)2(log 𝑁𝑐 𝑚,𝑖)(𝛽𝑓 𝑚,𝑖)2⎤⎥⎥⎥⎥⎦ .(13) According to Fig. 3, when a cyclic load at a specific 𝑆𝑚𝑎𝑥 𝑚,𝑖 is applied within the critical fatigue stage, a unique 𝑓𝑟𝑒𝑑 𝑚,𝑖 function could be found which intersects the normalized S-N curve at 𝑁𝑓 𝑚,𝑖.𝑁𝑓 𝑚,𝑖 is a critical quantity which will be used in the next sessions for the jumping strategy. Its value could be obtained by setting 𝑁𝑐 𝑚,𝑖 =𝑁𝑓 𝑚,𝑖 in Eq. (13), then equating the outcome to the result from Eq. (10): 𝑁𝑓 𝑚,𝑖(𝑅𝑚,𝑖, 𝑆𝑚𝑎𝑥 𝑚,𝑖 ) = 10 ⎡⎢⎢⎢⎢⎢⎢⎣ ⎡⎢⎢⎢⎣ −1 𝛼𝑡 𝑚,𝑖(𝑅𝑚,𝑖)⋅ln⎛⎜⎜⎜⎝ 𝑆𝑚𝑎𝑥 𝑚,𝑖 −𝑆𝑡ℎ 𝑚,𝑖(𝑅𝑚,𝑖) 𝑆𝑢 𝑚,𝑖 −𝑆𝑡ℎ 𝑚,𝑖(𝑅𝑚,𝑖)⎞⎟⎟⎟⎠⎤⎥⎥⎥⎦ 1 𝛽𝑓 𝑚,𝑖 ⎤⎥⎥⎥⎥⎥⎥⎦(14) Generally, as long as 𝑁𝑐 𝑚,𝑖 < 𝑁𝑓 𝑚,𝑖, the non-linearities would be strength deterioration and fracture energy dissipation, according to the 𝑓𝑟𝑒𝑑 𝑚,𝑖 evolution. Once 𝑁𝑐 𝑚,𝑖 ≥𝑁𝑓 𝑚,𝑖, in addition to the strength, the stiffness will be degrading as well, which corresponds to a rapid failure process in the integration point level according to inequalities (6) and (7). Matlab software [60] is employed in the fitting process and to find proper values for the engaged parameters. To this end, the more experimental data is provided, the wider the reliable region will be in terms of stress ratio, 𝑅, and maximum stress, 𝑆𝑚𝑎𝑥. In other words, the model should be well-prepared prior to simulations to be able to treat as close as possible to the provided experimental Wöhler surface. In this way, the formulation could be applicable in different loading scenarios. 2.2.3. Over-critical fatigue stage Over-critical fatigue stage refers to the region where predictive stresses are greater than the static strength of the material. This situation happens as the crack propagates through the interface and the crack opening increases cycle-by-cycle, as is the case for the elements located within the red region in Fig. 2(c). The irreversible degradation process of the interface has to be reflected in the 𝑓𝑟𝑒𝑑 𝑚,𝑖 quantity. However, a different formulation than the one for the critical region should be employed. Consequently, a similar approach to [61] is adopted by considering a fatigue-induced degradation as: d𝑓𝑟𝑒𝑑 𝑚,𝑖 d𝑁𝑐 𝑚,𝑖 = − 1 𝛾𝑚,𝑖 (𝑁𝑐 𝑚,𝑖)−𝜂𝑚,𝑖 (15) where 𝛾𝑚,𝑖 and 𝜂𝑚,𝑖 are interface parameters describing the rate of degradation, and their value could be obtained according to the experimental campaign. 2.3. Fatigue delamination damage effect Once the failure criteria in Eqs. (6) and/or (7) are satisfied – either by evolution of 𝑓𝑟𝑒𝑑 𝑚,𝑖 or increasing the 𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑚,𝑖 – the damage variable for each mode 𝑚and interface 𝑖is calculated independently according to an exponential softening law [51,62]: 𝑑𝑚,𝑖 = 1 − 𝜎0,𝑡ℎ 𝑚,𝑖 𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑚,𝑖 exp [𝐴𝑚,𝑖(1 − 𝜎𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑚,𝑖 𝜎0,𝑡ℎ 𝑚,𝑖 )](16) where 𝐴𝑚,𝑖 =1 𝐶𝑚,𝑖𝐺𝑚,𝑖 (𝜎0,𝑡ℎ 𝑚,𝑖 )2𝑙𝑐 −1 2 (17) Composite Structures 351 (2025) 118613 4 A. Taherzadeh-Fard et al. Fig. 3. Normalized Wöhler and fatigue reduction function curves in a general loading state. 𝜎0,𝑡ℎ 𝑚,𝑖 ,𝐶𝑚,𝑖,𝐺𝑚,𝑖, and 𝑙𝑐being either mode I or mode II interfacial strength, modulus, fracture toughness and characteristic length, respectively. As explained in previous work [51], transverse Young’s and transverse shear moduli of the bulk layers are considered as the interfacial normal and shear moduli, whilst half of the element length is regarded as the characteristic length. Since there is not dedicated medium for the interface layer, delamination damage effects should be properly reflected in the virtual bulk layers so that the behavior of the structure at macro-scale is wellcaptured. To this end, each bulk layer is affected by two maximum normal (𝑑𝑛) and shear (𝑑𝑛) damages as below [51]: 𝜎𝑏𝑢𝑙𝑘, 𝑑 𝑧𝑧,𝑖 = (1 − max[𝑑𝑛,𝑖, 𝑑𝑛,𝑖−1])𝜎𝑏𝑢𝑙𝑘 𝑧𝑧,𝑖 (18) 𝜎𝑏𝑢𝑙𝑘, 𝑑 𝑥𝑧,𝑖 = (1 − max[𝑑𝑠,𝑖, 𝑑𝑠,𝑖−1])(1 − max[𝑑𝑛,𝑖, 𝑑𝑛,𝑖−1])𝜎𝑏𝑢𝑙𝑘 𝑥𝑧,𝑖 (19) 𝜎𝑏𝑢𝑙𝑘, 𝑑 𝑦𝑧,𝑖 = (1 − max[𝑑𝑠,𝑖, 𝑑𝑠,𝑖−1])(1 − max[𝑑𝑛,𝑖, 𝑑𝑛,𝑖−1])𝜎𝑏𝑢𝑙𝑘 𝑦𝑧,𝑖 (20) where (𝜎𝑏𝑢𝑙𝑘 𝑧𝑧,𝑖 , 𝜎𝑏𝑢𝑙𝑘 𝑥𝑧,𝑖 , 𝜎𝑏𝑢𝑙𝑘 𝑦𝑧,𝑖 )are predictive stresses acting on the delamination plane. They would be considered as the integrated stresses (𝜎𝑏𝑢𝑙𝑘, 𝑑 𝑧𝑧,𝑖 , 𝜎𝑏𝑢𝑙𝑘, 𝑑 𝑥𝑧,𝑖 and 𝜎𝑏𝑢𝑙𝑘, 𝑑 𝑦𝑧,𝑖 ) once affected by the delamination damage. All other stress components remain unaffected in terms of the delamination phenomenon [55]. To illustrate the performance of the presented formulation in different scenarios, two types of loading conditions are schematically presented in Fig. 4. In Fig. 4(a), a cyclic load is applied within the critical fatigue stage followed by a monotonic loading part. The model is expected to activate the first non-linearity as the result of the 𝑓𝑟𝑒𝑑 𝑚,𝑖 evolution and decrease the interface strength accordingly. However, in Fig. 4(b), a monotonic load is applied followed by a cyclic load within the over-critical fatigue region. As it is obvious, the response follows the same path as in the monotonic loading regime in the first cycle, whilst a progressive degradation in both strength and stiffness happens once the load is switched to a cyclic one. 2.4. Advance-in-time strategy According to Eqs. (13) and (15), the fatigue delamination method developed in the present study employs three independent quantities of 𝑁𝑐 𝑚,𝑖,𝑅𝑚,𝑖 and 𝑆𝑚𝑎𝑥 𝑚,𝑖 in the critical regime, and one quantity of 𝑁𝑐 𝑚,𝑖 in the over-critical region to identify the state of the sample within the cyclic loading. Although 𝑅𝑚,𝑖 and 𝑆𝑚𝑎𝑥 𝑚,𝑖 are computed at each step and are completely dependent on the exerted external loads, 𝑁𝑐 𝑚,𝑖 is a continuous variable in the time domain and requires progressive growth within the simulation. This could be problematic, specifically in highcycle fatigue situations due to the high computational cost associated. Therefore, a procedure should be considered to skip a number of cycles while compensating for its effects within the formulation to reduce the computational costs. To this end, an advance-in-time (AIT) strategy is considered once stable conditions are obtained in terms of the cyclic loading exerted. The method is based on the jump strategy proposed in [49,50,57,63], however it is enhanced to account for two different loading mode conditions. Indeed, previous version of the AIT method has been presented for isotropic materials where the distinction between mode I or mode II loadings was not necessary. In the case of the numerical framework proposed in this study, that is no longer the case, a method is developed to identify to most critical mode of loading, and the jump is set up accordingly. Let us consider that mode I and mode II loadings are activated in the critical fatigue region, as shown schematically in Fig. 5. To maintain generality, loading modes are referred to as modes 𝑘and 𝑙. Depending on the time when either mode has reached the critical fatigue region, the number of cycles before jump (𝑁𝑐 𝑖) could be different, since in the sub-critical fatigue stage cycles are not counted according to Eq. (9). On the other hand, the required number of cycles to reach the damage initiation conditions could be obtained using Eq. (14) for each loading mode, independently. It would be reasonable to employ the most critical loading mode in terms of the distance to the damage initiation conditions, so that the other mode remains in the allowable S-N area after the jump. As a result, the loading mode with the minimum amount of 𝑁𝑓−𝑁𝑐will be fed into the AIT strategy. Other combinations of loading modes are still possible, where the selection of the mode for the AIT strategy is quite straightforward. Within all other possibilities, loading mode in the over-critical fatigue stage would be selected for AIT strategy compared to the one in the critical fatigue stage, while critical fatigue stage will be prioritized over the sub-critical fatigue region. The jump is expected to occur once stable conditions are reached in the stress state. While the crack is propagating, the maximum cyclic stress is continuously evolving at each integration point. In this way, stability condition means a plateau in the damage and stress evolution which could be monitored employing two indicators: 𝜑𝑚,𝑖 =||||| 𝑆𝑚𝑎𝑥 𝑚,𝑖,𝑗+1 −𝑆𝑚𝑎𝑥 𝑚,𝑖,𝑗 𝑆𝑚𝑎𝑥 𝑚,𝑖,𝑗+1 |||||< 𝑡𝑜𝑙𝑒𝑟𝑎𝑛𝑐𝑒 (21) 𝜓𝑚,𝑖 =||||| 𝑅𝑚,𝑖,𝑗+1 −𝑅𝑚,𝑖,𝑗 𝑅𝑚,𝑖,𝑗+1 |||||< 𝑡𝑜𝑙𝑒𝑟𝑎𝑛𝑐𝑒 (22) Composite Structures 351 (2025) 118613 5 A. Taherzadeh-Fard et al. Fig. 4. Model performance in two different loading scenarios: (a) cyclic load followed by a monotonic tension within the critical area and (b) monotonic tension followed by a cyclic loading within the over-critical region. Fig. 5. Criticality level of two different loading modes in terms of the number of remaining cycles to failure. 𝜑𝑚,𝑖 and 𝜓𝑚,𝑖 being maximum stress and stress ratio stabilization norms, while 𝑆𝑚𝑎𝑥 𝑚,𝑖,𝑗 and 𝑅𝑚,𝑖,𝑗 being the maximum stress and the stress ratio of mode 𝑚, interface 𝑖and time step 𝑗, respectively. Once stability conditions in Eqs. (21) and (22) are satisfied, the process time as well as the number of cycles at each integration point will be updated, and the 𝑓𝑟𝑒𝑑 𝑚,𝑖 is calculated based on the newly-obtained number of cycles. 3. Numerical examples and results The potential of the proposed model is studied in this section through different loading scenarios in both damage initiation and damage propagation regimes. Unidirectional carbon fiber/epoxy prepreg (IM7/8552) properties have been employed, and the validations are conducted on this composite in all fatigue crack evolution assessments. It should be noted that, in the current study, the effect of fiber bridging has not been considered in the simulations. Indeed, the presented fatigue model acts in an integration-point-wise approach, where all the calculations are conducted at each Gauss-point within the finite element model. This implies that the basis of the fatigue model in current assessment is a local framework. However, taking into account the bridging effect necessitates additional information regarding the extent of delamination damage crack, since the R-curve phenomenon depends upon the length of the crack propagated. To this end, non-local methods should be taken into account to efficiently consider the fiber bridging effects. Ignorance of the R-curve effect would not impose any problems for the damage initiation stage, however, damage propagation would be affected by this assumption. It will be demonstrated that in the absence of the fiber bridging effects, the results still fall within an acceptable range of the experimental data. The presented fatigue formulation is implemented through the open-source code Kratos Multi-physics [64,65]. Geometry and mesh generation has been conducted utilizing pre and post processing software GiD [66]. 3.1. Constitutive law parameters calibration As it has been mentioned in Section 2.2.2, a Matlab fitting algorithm is employed to calibrate the fatigue parameters of the model. To this end, Wöhler diagrams of the interface in each loading mode is required to effectively characterize the fatigue delamination behavior. Transverse tension [18,67] and shear [22,68] fatigue tests have been conducted on the IM7/8552 carbon/epoxy composite at different stress ratios (𝑅𝑚,𝑖) and stress levels (𝑆𝑚𝑎𝑥 𝑚,𝑖 ) until the crack is initiated within the interface. The resulted experimental S-N diagrams are utilized to fit a curve in the form of Eq. (10) for each mode I and mode II. The experimental Wöhler curves along with the calibrated model responses are presented in Fig. 6. Tables 1 and 2present the material properties as well as the obtained values for the calibrated interfacial parameters, respectively. These properties will be used within the simulations in the following sections for both fatigue damage initiation and propagation regimes. Since all the simulations are conducted at 𝑅𝑚,𝑖 <1, there is no need to calibrate 𝑆𝑅2 𝑚,𝑖 and 𝐴𝑈𝑋𝑅2 𝑚,𝑖 parameters. Composite Structures 351 (2025) 118613 6 A. Taherzadeh-Fard et al. Fig. 6. Fatigue material model calibration in (a) transverse normal (mode I) and (b) shear (mode II) configurations compared to the experimental data [18,68]. Table 1 Mechanical properties of the IM7/8552 carbon/epoxy prepreg [22]. Mechanical property Value 𝐸1[GPa] 161 𝐸2=𝐸3[GPa] 11.4 𝐺12 =𝐺13 [GPa] 5.17 𝐺23 [GPa] 3.98 𝜈12 =𝜈13 0.32 𝜈23 0.43 3.2. Fatigue delamination damage initiation assessment In order to investigate the model capability to initiate fatigue crack in an undamaged medium, two standard tests have been employed in each loading mode, as illustrated in Fig. 7. Three-point bending test is utilized to initiate the fatigue crack in mode I loading. To this end, a central displacement is cycled with a frequency of 1 Hz while the two ends of the cantilever beam are set to be fixed through the whole width (𝑤= 6.35 mm), as shown in Fig. 7(a). Cyclic displacement domains are chosen in a way that they produce maximum normal stresses of 90,93,97,101, and 105.5 MPa within the central region. Carbon fibers are aligned in the loading direction and Table 2 Calibrated interfacial fatigue parameters for modes I and II [68]. Interface property Normal load Shear load (Mode I) (Mode II) Energy release rate (𝐺) [J/m2] 212 600 Interfacial strength (𝜎0 𝑡ℎ) [MPa] 111 85 Interfacial modulus (𝐶) [GPa] 11.4 5.17 𝑆𝑒[MPa] 87.7 39.1 𝑆𝑅1[MPa] 48.88 55.81 𝛼𝑓2.14e−3 1.69e−4 𝛽𝑓4.49 6.71 𝐴𝑈𝑋𝑅1−2.80e−3−2.58e−4 𝛾25 30 𝜂0.72 0.72 in 90◦angle with respect to the longitudinal direction. Hexahedron elements with the size of 0.1 mm are employed in the central region, where the crack is expected to initiate, and coarse meshes are opted for the far-field area. The cycle is defined to induce 0.1 stress ratio states according to the experimental procedure, and the results are presented in Fig. 8 compared against the experimental data. Within the experimental tests, specimens in both standard tests fail catastrophically once damage has initiated, which allows simple detection of the damage Composite Structures 351 (2025) 118613 7 A. Taherzadeh-Fard et al. Fig. 7. Schematic view of (a) three-point bending and (b) double-notched shear tests under cyclic load (dimensions in millimeters). Fig. 8. Initiation of the fatigue crack in the three-point bending configuration under mode I loading compared with the experimental data at 𝑅= 0.1[18,23]. onset. In the numerical framework, however, the necessary number of cycles to initiate the delamination damage is calculated according to Eq. (14). This value has been considered as the crack initiation point within the simulations. As it can be observed, the model captures the crack initiation cycle fairly well at each stress level, and the numerical output lies within the experimental scatter region. In the simulations, fatigue crack propagates catastrophically once it has been initiated, which aligns with the experimental findings [18,23]. In Fig. 9, the crack initiation step and its complete propagation is presented at the stress level of 97 MPa. Through the simulation, a deactivation algorithm is employed to erase fully damaged elements so that the model converges more easily. Validating the performance of the model in mode I loading, it would be valuable to check its applicability in shear mode as well. Therefore, experimental double-notched shear test results are utilized for running pure mode II simulations [68]. The configuration of the test is shown in Fig. 7(b) where carbon fibers are aligned in the longitudinal direction. A cyclic axial displacement is applied to one end through the whole width (𝑤= 12.7 mm) so that it produces shear stresses at the central part. From the numerical point of view, considering solely the central part would be enough to run the simulations. Two stress ratios of 0.1 and 0.3 are used along with a cyclic load of 1 Hz. The element size and type are the same as in the three-point bending test. Fig. 10 summarizes the numerical and experimental results. The model output for mode II loading shows a good agreement with the experimental results. This indicates that the proposed formulation would be able to deal with the crack initiation phenomenon in both mode I and mode II loading conditions in an accurate manner. The next stage is assessing the performance of the model in the propagation regime, which will be discussed in the following section. Composite Structures 351 (2025) 118613 8 A. Taherzadeh-Fard et al. Fig. 9. Three-point bending fatigue simulations in the stress level of 97 MPa in (a) initiation stage at 𝑁𝑐= 218 964 and (b) complete failure at 𝑁𝑐= 218 986. Fig. 10. Initiation of the fatigue crack in double-notch shear configuration in (a) 𝑅= 0.1and (b) 𝑅= 0.3compared with the experimental data [68]. Composite Structures 351 (2025) 118613 9 A. Taherzadeh-Fard et al. [50] Gonçalves LA, Jiménez S, Cornejo A, Barbu L, Parareda S, Casellas D. Numerical simulation of a rapid fatigue test of high Mn-TWIP steel via a high cycle fatigue constitutive law. Int J Fatigue 2023;168:107444. http://dx.doi.org/10.1016/j. ijfatigue.2022.107444, URL https://www.sciencedirect.com/science/article/pii/ S0142112322006946. [51] Taherzadeh-Fard A, Cornejo A, Jiménez S, Barbu LG. A rule of mixtures approach for delamination damage analysis in composite materials. Compos Sci Technol 2023;242:110160. http://dx.doi.org/10.1016/j.compscitech.2023.110160, URL https://www.sciencedirect.com/science/article/pii/S0266353823002531. [52] Truesdell C, Toupin R. The classical field theories. In: Principles of classical mechanics and field theory/prinzipien der klassischen mechanik und feldtheorie. Springer; 1960, p. 226–858. [53] Cornejo Velázquez A. A fully Lagrangian formulation for fluid-structure interaction between free-surface flows and multi-fracturing solids and structures (Thesis), 2020. [54] Car E, Oller S, Oñate E. An anisotropic elastoplastic constitutive model for large strain analysis of fiber reinforced composite materials. Comput Methods Appl Mech Engrg 2000;185(2):245–77. http://dx.doi.org/10.1016/ S0045-7825(99)00262-5, URL https://www.sciencedirect.com/science/article/ pii/S0045782599002625. [55] Balzani C, Wagner W. An interface element for the simulation of delamination in unidirectional fiber-reinforced composite laminates. Eng Fract Mech 2008;75(9):2597–615. [56] Gonçalves LA, Jiménez S, Cornejo A, Tedesco MM, Barbu LG. A high cycle fatigue numerical framework for component-level virtual fatigue testing: Application to a light-duty vehicle lower control arm. Eng Struct 2024;311:118198. [57] Alcayde B, Merzkirch M, Cornejo A, Jiménez S, Marklund E, Barbu L. Fatigue behaviour of glass-fibre-reinforced polymers: Numerical and experimental characterisation. Compos Struct 2024;337:118057. http://dx.doi.org/10.1016/ j.compstruct.2024.118057, URL https://www.sciencedirect.com/science/article/ pii/S0263822324001855. [58] Dávila CG, Joosten MW. A cohesive fatigue model for composite delamination based on a new material characterization procedure for the Paris law. Eng Fract Mech 2023;284:109232. [59] Dávila CG. From SN to the Paris law with a new mixed-mode cohesive fatigue model. Tech. rep., 2018. [60] The MathWorks, Inc. MATLAB R2013a (8.1.0.604). Natick, Massachusetts, USA; 2013, Software. [61] Maiti S, Geubelle PH. A cohesive model for fatigue failure of polymers. Eng Fract Mech 2005;72(5):691–708. [62] Oller S. Nonlinear dynamics of structures. Springer; 2014. [63] Barbu LG, Oller S, Martinez X, Barbat A. High cycle fatigue simulation: A new stepwise load-advancing strategy. Eng Struct 2015;97:118–29. http://dx.doi.org/ 10.1016/j.engstruct.2015.04.012, URL https://www.sciencedirect.com/science/ article/pii/S0141029615002527. [64] Ferrándiz VM, Bucher P, Rossi R, Cotela J, Carbonell J, Zorrilla R, Tosi R. KratosMultiphysics (version 8.0). 2020. [65] Dadvand P, Rossi R, Oñate E. An object-oriented environment for developing finite element codes for multi-disciplinary applications. Arch Comput Methods Eng 2010;17:253–97. [66] Ribó R, Pasenau M, Escolano E. GiD user manual. CIMNE; 2007. [67] Dávila CG, Rose CA, Murri GB, Jackson WC, Johnston WM. Evaluation of fatigue damage accumulation functions for delamination initiation and propagation. Tech. rep., 2020. [68] May M, Hallett SR. An assessment of through-thickness shear tests for initiation of fatigue failure. Composites A 2010;41(11):1570–8. [69] Murri GB. Evaluation of delamination onset and growth characterization methods under mode I fatigue loading. Tech. rep., 2013. [70] Nojavan S, Schesser D, Yang Q. An in situ fatigue-CZM for unified crack initiation and propagation in composites under cyclic loading. Compos Struct 2016;146:34–49. [71] Reeder JR, Demarco K, Whitley KS. The use of doubler reinforcement in delamination toughness testing. Composites A 2004;35(11):1337–44. [72] O’Brien TK, Johnston WM, Toland GJ. Mode II interlaminar fracture toughness and fatigue characterization of a graphite epoxy composite material. Tech. rep., 2010. Composite Structures 351 (2025) 118613 16