Cohesive zone modelling to predict crack growth under fatigue loading

This manuscript presents the partial results in the development of numerical models based on Cohesive Zone Modelling to analyse quasi-static and fatigue crack propagation in adhesive interfaces. This work applied to adhesive interface between fibre laminates and structural epoxy adhesive of a wind turbine blade 3rd shear-web connection to the trailing edge panels. Fracture characterization experiments using the Double Cantilever Beam with Unequal Bending Moments (DCB-UBM) specimen are performed, and results are used to determine different modelling parameters, such as the mixed-mode cohesive laws and fatigue damage parameters for cohesive crack propagation. The developed numerical model is compared and validated against sub-component testing on two different scale levels, using a Modified Tear Test (MTT) specimen geometry with simplified geometry and loading, and a more complex geometry sub-component specimen with a simplified Shear Web Disbond Test (SBDT). The work is carried out as part of the CORTIR-II project, funded by the Danish Ministry of Energy.


Introduction
A wind turbine's blade comprises an outer aerodynamic shell, and shear webs that join the upwind and downwind sections of that shell.These webs are bonded to one another at the leading and trailing edges.Typically, a blade has two shear webs that run almost the entire length of the blade from the transition zone (denoted by the shear web fish mouth).Due to the increasing global demand for clean and sustainable energy, wind turbines have become more prevalent over the past few decades.This surge in interest has driven manufacturers to develop larger wind turbines that are more efficient and cost-effective than their smaller counterparts.However, the continuous increase in size has resulted in higher loads, which require additional strengthening.This challenge is addressed by adding a third shear web to the blades.Unfortunately, this leads to a failure phenomenon known as shear web disbonding, where cracks in the adhesive interface between the shear web footing and spar cap propagate due to large deflections in the aerodynamic skins near the third shear web.If the crack grows too large this could compromise the integrity of the blade, eventually causing the wind turbine to collapse.The Technical University of Denmark (DTU), Aalborg University and Bladena have partnered with wind turbine owners and operators, to initiate the CORTIR-II project.which aims to investigate shear web disbonding in large wind turbine blades.The project is divided into four stages to investigate the issue in detail.It starts with characterisation and moves on to sub-component testing with more detail, before finally conducting large-scale testing.

Figure 2 Building block approach for Cortir II project
This work presents the current state of part of this project which is focused on developing numerical models based on Cohesive Zone Modelling to analyze quasi-static and fatigue crack propagation in the adhesive interface between fibre laminates and structural epoxy adhesive.The modelling parameters are determined from results obtained in an experimental campaign using the DCB-UBM test [1], , which is a novel test method for mixed-mode fracture characterization.The developed numerical modelling will be validated using the Modified Tear Test (MTT), and subcomponent testing method.

Cohesive zone modelling
The fibre bridging phenomenon occurs when the fibres in neighbouring laminate plies bridge the debonding plane, which constitutes an additional mechanism of energy dissipation during crack propagation .As a result, the debonding resistance and inter-laminar fracture toughness are enhanced.Since no purely analytical solution is available to directly determine a bridging/cohesive law it is necessary to perform fracture experiments that allow to compute the J integral.This since the J 1293 (2023) 012012 IOP Publishing doi:10.1088/1757-899X/1293/1/0120123 integral can be used to study non-linear fracture mechanics problems involving large fracture process zones, such as fibre bridging (Sørensen [2] ).
Cohesive Zone Modeling (CZM) is a useful tool for predicting and modelling crack initiation and propagation in situations where Linear Elastic Fracture Mechanics (LEFM) cannot be assumed.For this reason and because they are widely implemented in commercial FEA packages over recent decades, researchers have increasingly shown interest in CZM applied in finite element modelling.CZM can incorporate nonlinear behaviour into the material's constitutive response, making it capable of predicting failure in ductile crack problems with significant process zones.The length of the process zone is critical in the fracture process as it can indicate if the fracture is ductile or brittle.If the process zone is much smaller than other characteristic lengths of the specimen such as crack length, specimen length, width, etc., brittle fracture occurs.CZM has previously been applied to the study of various failures in concrete, metals, composites, bi-material interfaces, and adhesive joints.In the CZM the fracture process zone (FPZ) is substituted with a cohesive region, where the exchange of stress between the two artificial crack surfaces is defined using a relationship between the traction and separation of the FPZ.In practical FEA analysis, the cohesive zone model is commonly applied by a single layer of cohesive elements, which are placed along the expected path of crack propagation.Each element follows a constitutive response which uses a strength-based criterion to predict the onset of damage and the degradation of the material stiffness.

Traction separation law
Commercial FEA software packages such as Abaqus required a set of parameters that govern the traction-separation law, for example in the case of pure mode I fracture: •   , the stiffness of the cohesive element in the linear elastic regime • Traction   0 (or separation   0 ) for damage onset •    , fracture energy release rate • Damage evolution law There are multiple traction-separation laws used in literature.Some of the most common bi-linear, exponential, and trapezoidal behaviours although non-standard (tabulated) laws computed directly from experimental tests are also employed.It is generally agreed that the key parameters of traction-separation law previously mentioned should be obtained from experiments although some other studies have used calibration methods based on trial-and-error search or a combination of both.
One of the aims of the present study is to obtain traction-separation law directly from an experimental fracture test following the methodology proposed by [2] which involves the computation of the J integral and the measurement of the crack tip separation in a quasi-static DCB-UBM test.In this method, the traction-separation law is obtained by the partial differentiation of the fracture resistance   and the crack tip separation.In pure mode I (opening) for example, the traction-separation law can be obtained by:   (  ,   ) = (  (  ,   )) (  ) for model II (shear): (  ,   ) = (  (  ,   )) (  )

Modelling fatigue with cohesive zone models
A novel approach to consider the deterioration of materials caused by fatigue loading is to degrade the traction-separation response as the number of cycles increases (Khoramishad et al. [4]).This model was initially developed for aluminium joints with epoxy adhesive and it was validated against single lap joint (SLJ) and laminated doublers in bending (LDB) fatigue tests.In this approach the cohesive strength (  ) and critical fracture energy ( , ) decrease through the damage parameter (  ) which follow a fatigue damage evolution law of the following form: where   is the damage accumulated over ∆N cycles while α, β and  ℎ are material parameters which must be calibrated against experimental results.  is the maximum overall strain defined in terms of the normal and shear strains in the cohesive element: In the beginning, the damage parameter   for each element within the cohesive zone is set to zero and damage accumulation in the load cycle only occurs if the maximum overall strain is larger than the threshold value  ℎ .The damage parameter which ranges between 0 and 1 can be computed after each load cycle (or set of load cycles) and update the effective traction and energy release rate used to estimate the crack propagation.  =    0  , =    ,0 Where i denotes the current time step, and  0 and  ,0 is the initial cohesive strength and fracture energy respectively.

DCB-UBM tests (Level -0)
DCB-UBM test is the first of three different experimental tests that are carried out to study the quasistatic and fatigue debonding phenomena.This test campaign is aimed to assess the critical fracture parameters such as energy release rates and cohesive damage law under various mixed-mode loading conditions.
For validation purposes, two additional tests are proposed: First, a new test approach, referred to as the Modified Tear Test (MTT), derived from the Sandwich Tear Test (STT), is specifically designed for Sub-component level 1 testing [3].And second a setup for a simplified shear web geometry (Subcomponent II testing) which considers the knowledge acquired from conducting the aforementioned tests and includes external load conditions closer to those encountered in a full-scale blade scenario.

Test set-up
The DCB-UBM test, as introduced by [2], involves the application of uneven bending moments to the crack flanks' tips in the specimen.By adjusting the ratio of these applied moments, crack propagation can occur under different fracture modes, ranging from pure mode I to pure mode II.Since the load applied to the specimen solely originates from bending moments, both the energy release rate and the mode-mixity angle remain unaffected by the crack length.Another advantage of employing pure bending moments for loading, particularly in cases involving large process zones, is that the applied load remains constant throughout the entire process zone.The tests were performed in a DCB-UBM servo-hydraulic testing machine developed in-house at DTU. Pure moments are applied over each laminate only on one side of the specimen by two separate servo-hydraulic torsional actuators.To maintain pure moment loading, it is required to minimize any in-plane and out-of-plane forces, so the carriage plates are mounted on raceway shafts, utilizing track rollers.
The DCB-UBM specimens have overall dimensions of 450 mm x 30 mm x 18.2 mm and are made from two GFRP laminates bonded on top of each other (Figure 6).The specimens were manufactured by Global Wind Services by hand layup, manufacturing the laminates separately and then bonding them together by secondary bonding and vacuum bag consolidation.A 50 mm pre-crack is made by inserting a Teflon film during the consolidation process.The top laminate is made from a triaxial/biaxial fabric with the following stacking sequence [45/0/-45/ ±45]8 and an average cured thickness of 10 mm.The bottom laminate is made from biaxial fabric with a stacking sequence of [±45]10 and an average cured thickness of 6.8 mm.The specimen configuration was determined through preliminary studies, utilizing Crack Surface Displacement Extrapolation (CSDE) analysis to calculate mode mixity angles and energy release rates.Due to the significant tip rotations, potential laminate failure and limited fixation options, two 6 mm thick steel reinforcement plates (doublers) were bonded to the top and bottom surfaces of the specimen using an epoxy structural adhesive.Additionally, a clamp was installed to prevent crack initiation between the specimen and plates.The moments and rotations of the loading arms are directly obtained from the sensors of the test machine through the MTS data acquisition system.Additional data such as crack length and separation of the crack tip is collected using the iMetrum Non-Contact Video Measurement live-tracking software.The displacements of the tracked crack flanks, which are utilized to determine the phase angle, should be closely monitored near the crack tip and adjacent to the crack flank surface.To ensure precise tracking of crack propagation, it was crucial to accurately define the crack path and initial crack length.Additionally, to enhance the image correlation results, all specimens' surfaces were meticulously sanded and painted white, and speckle patterns of approximately 1mm were applied to the steel plates using a sharp marker.Camera and lightening were adjusted as well as the DIC threshold values until a reliable crack propagation was achieved.Quasi-static and fatigue tests were conducted to establish the key parameters of the cohesive laws.Three moment ratios (MR = M2/M1) were tested MR=-1.2, -0.714, and -0.4 corresponding to phase angles equal to 7.4°, -6.5° and -22.4° respectively.Phase angles larger than 30° (including Mode II) were not considered due to limitations in the load capacity of the testing machine.At least three specimens were tested for each phase angle.
Fatigue tests were performed at three different load levels (M1max) corresponding to 60%, 70%, and 80% of the critical energy release rate (G/Gc) found in quasi-static tests.All tests were performed with a load ratio R = M1min/M1max= 0.2 at the same moment ratios (MR) studied in the quasi-static test.The crack extension (Δa) is measured directly from the photo frames at different numbers of cycles in such a way that the crack propagation rate (Δa/ΔN) is estimated between two adjacent measurement points.Live-tracking methods were not suitable for this fatigue test due to the relative movement of the specimens and the opening/closing of the crack flanks.Other technical challenges had to be sorted out to perform this test, including a sudden detachment of steel doubler reinforcements and high vibrations.However, the most difficulty was to maintain an accurate moment ratio between the arms when running the experiment at high frequency.In such a case the time delay in the response of both arms made the controller underestimate the moment on arm 2. Thus, the input moment ratio had to be slightly calibrated to match the desired value.

Tests results
It was discovered that the interface failure was primarily controlled by extensive fibre bridging on a large scale, both during testing under quasi-static and fatigue loading.Fibre bridging was observed across all tested modes, with its greatest impact occurring in the proximity of mode I.The presence of voids within the adhesive, as well as variations in the thickness of the bondline, and internal preload caused by imbalanced composite face-sheets, result in variations in the fracture resistance of each specimen.In some instances, the crack deviated towards the opposite interface from the location of the pre-crack, forming a kink.Besides this, most experiments were successfully performed.The obtained data was post-processed in MATLAB to obtain the fracture resistance R-curve by computing the J-integral at each instant taking as input the applied moments (M1 and M2).The Jintegral formulation presented by [5] was modified to account for the asymmetry due to the difference in the thickness and material properties of the top and bottom laminates.The crack extension (Δa) and tip-end separation (δ n , δ t ) were also plotted.An example of the post-processed data is in Figure 8. b) The fracture resistance for crack initiation (  ) and steady-state crack growth (  ) are compared at different phase angles.Results show that   is larger than   for all studied phase angles which shows the relevance of the fibre-bridging toughening effect, particularly at low phase angles.
To obtain the cohesive traction-separation laws it is convenient to graph the fracture resistance for each test as a function of the normal crack end opening and the tangential crack end opening.Due to the high data scattering a non-linear least square surface polynomial fitting is utilized, yielding a parametric function that enables the estimation of the traction-separation cohesive laws.For the fatigue tests the crack propagation rates (da/dN) were estimated at different load levels (G/Gc) and moment ratios (M2/M1).As expected, the crack propagation rate increases with the load level following a power function similar to Pari's law [6](Figure 11 and Figure 12).However large scattering was obtained for da/dN measured despite the load level and moment ratio being the same.Through more careful analysis, it was determined that all da/dN measurements were taken at different effective crack lengths.As a result, the size of the bridging zone for each measurement was different, and this could explain the observed data scattering.In fact, crack propagation rates measured at the beginning of the test (when the crack length is small), were much larger than those da/dN measured when the crack has extensively propagated.This suggests that the fatigue crack propagation rate is not only a function of the moment ratio but also a function of the size of the bridging zone (Figure 13).From the same figure, it can be also observed that at large lengths of the bridging zone, the value da/dN asymptotically reaches a constant value which can be related to a fully developed bridging zone and the steady-state critical fracture energy (GcI = Jss) determined through the R-curve (section 3.2.) This dependency has also been previously reported for carbon/epoxy laminates (e.g.L. Yao et al. [7] and G. Murri [8]) and alternative experimental procedures and reduction methods have been proposed.In this regard, a more detailed test campaign is ongoing at DTU Construct to characterize da/dN and account for its dependency on the bridging zone length.The laminates and adhesive layer are modelled as separate geometric bodies joined at the interface using a tie constraint.The laminates are modelled as orthotropic materials while the adhesive material and the material for the doubler stiffeners an isotropic materials.Material properties for the laminates were obtained by coupon testing during Cortir I project and are summarized in Table 1.The models are mesh using 2D CPE8R, quadratic plane strain elements with reduced integration.Two concentrated moments (M1 and M2) are applied over reference points located on the left side of the specimen which are coupled to the specimen's external surface by a rigid body constraint.The specimen is restrained on the opposite side simulating the rollers Quasi-static tests are run using the standard/linear solver while the fatigue test uses an explicit solver.For the fatigue model, the damage criteria described in section 2.2. was implemented in a VUSDFLD Fortran subroutine.This subroutine enables the dynamic adjustment of field variables at a specific material point.These adjustments can be defined as functions of time or any of the material point quantities listed in the "Available output variable keys."This capability is particularly useful for introducing material properties that depend on the solution, as these properties can be conveniently expressed as functions of existing field variables.Since modelling fatigue loading cycle-by-cycle is computationally non-practical a different strategy is adopted in which bundles of multiple load cycles are represented into a single time step.

FEA model results of DCB-UBM test
The quasi-static cohesive FEA model is primarily used to calibrate and validate the cohesive law estimated analytically from the experimental test results.Due to the large data scattering in the experimental results, the calibration procedure requires a clear measure of the error between experimental and FEA results.Since the primary objective of a cohesive law approach is to be able to predict both crack initiation and crack growth, both measurements are used to define a combined mean squared error of the following form: Where  1 and are the moment in arm 1 and crack growth rate at crack initiation onset and the superscripts FEA and EXP mean FEA results and experimental results respectively.Similarly, the fatigue FEA model is used to calibrate the fatigue damage parameters α, β,  ℎ against the DCB-UBM test results (section 3.2.).Since the definition of those parameters in the material model is purely phenomenological and cannot be analytically related to any physical measurable quantity it is then necessary to obtain these values through a search method.Numerous attempts were carried out to try to fit these material parameters manually using a trial-and-error technique, however, no successful combination has been found yet.For the sake of computational efficiency, more advanced search algorithms for global optimization are proposed for future analysis including pattern search, multi-start search, genetic algorithms, or surrogate optimization.Despite this, results from the analysed parameter combinations show that the model can capture the dependency of the crack growth rate (da/dN) with the load level (G/Gc) as well as its increase following a power law.This indicates that the implemented formulation might be suitable to analyse fatigue in fully developed bridging zones.Comparing the accuracy of fatigue models with different ∆N jump values is not straightforward because ∆N serves as a prerequisite value in determining other fatigue parameters.Consequently, the fatigue parameters that yield a good fit with the data for a particular ∆N are only applicable to that chosen ∆N.Nonetheless, when comparing models based on time rather than load cycles and using the same fatigue parameters, it can be demonstrated that ∆N has a limited impact on the results.method, which is utilized to examine interface failure between the adhesive and laminates while considering different mixed mode conditions.In this test, the upper laminate is pulled upwards, so the in-plane membrane forces escalate significantly as the test proceeds.This escalation leads to an increase in the mode-mixity angle, transitioning from an initial pure mode I opening condition to mixed fracture modes (opening and shear).A new test set-up was developed and produced, incorporating some components from a previous STT setup while redesigning and manufacturing the majority of them specifically for this test.

Conclusions
A building block methodology has been implemented to study crack initiation and propagation in bonded joints of glass fibre laminates and epoxy adhesives, typically encountered in shear webs in wind turbine blades.The phenomena have been studied from an experimental and numerical perspective at three different levels of complexity.This section discusses the partial results of the firstlevel DCB-UBM characterization.

Experimental study:
Quasi-static DCB-UBM tests show that the crack propagation behaviour is dominated by a large fibre bridging zone in the bondline which acts as a toughening mechanism, particularly at small mode mixites phase angles (mode I, and opening-dominated fracture).Measurements of the fracture resistance (R) through the J integral show a large increase from the crack onset (Jo) to steady-state propagation with fully developed bridging (Jss).The measurements of the fracture energy release rate were highly scattered probably due to differences in the bondline thickness, void content and curvature (due to the asymmetry of top and bottom laminates).Since the manufacturing method, materials and sizes used in the specimen fabrication are representative of those employed in the wind industry for similar applications similar deviations in the fracture resistance could be expected in operating wind turbine blades.It was also possible to determine the traction-separation law directly from the experiments for a range of phase angles by measuring crack tip displacements and applying the data reduction proposed by Sorensen together with nonlinear surface extrapolation.However, due to limitations in the test machine, it was not possible to test phase angles near mode II and therefore the validity of the predicted fracture behaviour in this region is still unknown.
Fatigue DCB-UBM tests were found to be challenging due to large machine vibrations, and difficulties to set the PDI controller, for this reason, a limited number of tests have been performed.Besides this, it was possible to measure the crack growth rate at different loading conditions and results suggest that the crack growth rate (da/dN) follows a power law function which depends on the mode mixity and load level (G/Gc).Due to the large data scattering it was not possible to estimate this dependency quantitatively with enough reliability so additional tests are carried out at this moment to reduce the statistical uncertainty of the curve fitting estimation.Additionally, it was observed that the crack growth rate (da/dN) is highly affected by the length of the bridging zone.A similar dependency of this property to the R-curve has been found by other authors in the delamination of CFRP composites with large fibre bridging zones.

Numerical study:
For characterization (DCB-UBM) tests, 2D finite element analysis (FEA) models have been developed to calibrate the cohesive traction-separation laws obtained in experimental tests.Different traction separation laws have been tested and compared including bi-linear, exponential, and tabulated (based on traction separation curves obtained in experiments).For mode-I tabulated cohesive-law can capture with accuracy the crack onset and crack propagation observed in experimental testing.However, minor calibrations were required for some cohesive law parameters using trial-and-error search methods.
A novel approach to model fatigue crack growth in the bondline was successfully implemented through a material user subroutine in which a fatigue damage law is coupled with the cycle jump technique.The same 2D FEA model was used to calibrate the fatigue damage law parameters using trial-and-error search methods.However, it was not possible to find a combination of material parameters that fully validate the experimental data.Despite this, preliminary results suggest that the proposed approach can capture with some accuracy the crack growth rate of a fully developed bridging zone but further refinement of the experimental test and numerical model is required.

Figure 1 .
Figure 1.Example of a wind turbine blade section

Figure 3 .
Figure3.Different traction-separation functions evaluated in the current project.Source:[3].Here J0 represents fracture energy release rate (computed through J integral) for crack onset while Jss represents the fracture energy release rate for steady-state propagation which is also referred as Gc CZM also allows studying the fracture process when mode-mixity is present using the power law fracture criterion based on the energy required to propagate a crack under pure mode conditions.

Figure 4 .
Figure 4. Mixed mode traction-separation and criteria for damage onset and failure

Figure 5 .
Figure 5. Evolution of the traction-separation law due to fatigue cycles.Source:[3]

Figure 10
Figure 10.a) Surface polynomial fitting, b) Cohesive traction separation law obtained analytically from quasistatic DCB-UBM tests for pure Mode I,

da/dN = 2 .Figure 13
Figure 13 Crack growth rate at different bridging zone lengths.G/GcI = 0.6 and M2/M1 = -1.23.3.FEA models of DCB-UBM tests.An FEA model was created in Abaqus 2021 to replicate the experimental DCB-UBM test explained in the previous subsection.This model uses a cohesive zone model to predict damage initiation and crack propagation both in quasi-static and fatigue loading.Cohesive elements COH2D4 were used together with element deletion to track crack propagation in two alternative ways: when damage initiates (damage variable D>0) and complete failure occurs (D=1).

Figure 14 .
Figure 14.FEA model for DCB-UBM test using the cohesive zone method.

4 .
Validation tests and FEA model Additional tests at higher complexity levels of the test pyramid are currently ongoing at DTU Construct.A brief description of those tests is presented below: a) MTT tests (Level -1): The objective of conducting sub-component level I tests is to replicate the real-life occurrence of shear web disbonding using simple test setups and specimens.While it is commonly assumed that shear web disbonding occurs under the mode I condition, the presence of mode mixes can occur due to membrane forces in aerodynamic shells or the torsional shear induced by global loading on wind turbine blades.To explore how these factors might influence the problem, it was necessary to devise a test method in which mode mixity varies during crack propagation.Additionally, this test enables the validation of the cohesive law parameters across different mode mixity angles that are estimated in level-0 of the building block approach.The Modified Tear Test (MTT) is derived from the Sandwich Tear Test (STT)

Figure 18 .
Figure 18.MTT test setup.Source: [3] b) Subcomponent Level II test -Shear Web Disbond Test (SWDT): The purpose of the subcomponent level II test is to study the simplified geometry and loading conditions of a real shear web installed in a wind turbine blade paying particular attention to the bonded interface between the shear web flange and the blade outer shell.The test specimens are a representative simplified version of only half the shear web and a portion of the outer shell to which it is bonded with a pre-crack inserted in the bondline.Four specimens were manufactured by Global Wind Services, manufacturing the shear webs and laminates separately by hand layup and vacuum bag consolidation and then joining them by secondary bonding.The top laminate and shear web flange have the same composite layup as those used in level-I and level-II tests.

Table 1 .
Laminate material properties Figure 17.Example of fatigue crack growth estimations predicted by the implemented formulation with different combination of parameters.