Prediction of crack growth at trailing edge bondlines of a wind turbine rotor blade for the assessment of remaining service life

Debonding along trailing-edge joint of a wind turbine rotor blade is one of the often reported failure modes during the rotor blade maintenance inspections. This study focuses on buckling-induced debonding along the trailing edge under in-plane compressive loads due to the gravity forces on the rotor blade. With a priori known debond crack length, an analytical model is used for predicting debond growth rate and assessing the time that takes the crack to reach a critical length. For the growth rate prediction, required input parameters are size of the initial debonding, laminate and adhesive layer thicknesses, elastic and interfacial fracture mechanics properties, gravitational load distribution and position of the initial debonding. The damage model is implemented in excel spreadsheet tool. The predictions of the analytical model were compared with the results from 3D finite element simulations and the difference between the results (energy release rate) was 8%. This study demonstrates how damage specific tools can be developed for predicting remaining useful life of rotor blades.


Introduction
Modern wind turbine rotor blades typically use load-carrying composite laminates in spar caps, leading edge and trailing edge regions of their cross-section.The laminates located at the trailing edge mainly take up the bending moments associated with the gravitational loads i.e., the mass of the blade [1].Consider point A on the trailing edge of a rotor blade, as shown in Fig. 1a.During the half rotation, due to gravity the point A is loaded in tension and for the remaining rotation, the loading type gets reversed resulting in a tension-compression cyclic loading history.With the increase in size of the rotor blades, weight significantly increases and consequently, gravitational loads become one of the key design factors.Moreover, wind turbines are designed for a service life of 20-25 years.To ensure a safe operational life of the blades, it is of great importance to be able to predict the fatigue behaviour of composite/sandwich material interfaces, and adhesive bondlines.
Adhesive joints in the rotor blade are highly susceptible to damage in the form of interfacial cracks, particularly, a joint failure in the trailing edge (see Fig. 1b) can lead to final collapse of the blade structure [2].The trailing edge can develop at least two types of bondline damage: tunnelling cracks (cracking in the adhesive along its width direction) and interfacial debonding, i.e., cracking along the interface between laminate and adhesive layer (in the blade length 1293 (2023) 012037 IOP Publishing doi:10.1088/1757-899X/1293/1/012037 2 direction).While tunnelling cracks alone do not cause failure [3], in some cases can lead to the formation of interfacial cracks which can eventually fail by buckling of the blade [4].1.1.Background Ataya and Ahmed [2] visually inspected numerous rotor blades of two different lengths 9.5 m and 14.2 m for damages at the trailing edge.The major damage type observed was the debond cracking along the laminate/adhesive interface starting from the rotor blade root section towards the tip.The damage formed was due to a combination of several factors such as extreme loading conditions, complex geometries and imperfections coming from manufacturing processes.More importantly, the length of the debond cracks varied from a few centimetres up to approx.1.5 m.Haselbach et al. [5] tested a 34 m long blade by four-point loads and found that the local buckling of the trailing edge led to the debonding of the joint.It was concluded that the edgewise loads are the primary loads driving the failure of the trailing edge.
Full-scale blade test performed by Eder and Bitsche [6] on a 34 m blade revealed that the initiation of a debond crack occurred in the vicinity of local buckled region on the trailing edge and then propagated all along the joint.It was reported that the most critical failure mode is the normal opening of the trailing edge and a method based on virtual crack closure technique was proposed to estimate the energy release rate of the debond crack-tip [7].However, the method did not considered the edgewise loading that could lead to local buckling of the trailing edge laminate which can lead to debonding.
Subcomponent testing on cut-section of a 3 m long trailing edge panel was performed by Chen et al. [8].The cut panel was loaded in compression and the sequence of damage modes were documented.With the increase in the applied load, the trailing edge panel buckled followed by adhesive debonding.Post the adhesive joint failure, damage in composite laminate was observed in the form of delamination and fibre-dominant mode.A non-linear finite element model was IOP Publishing doi:10.1088/1757-899X/1293/1/0120373 developed in the study to account for buckling, contact and matrix/fibre damages to predict the failure evolution.However, boundary conditions considered in the numerical model were based on subcomponent testing and not the full-scale blade.
In order to predict the trailing edge damage growth, detailed yet simple models are needed that takes into account the i) correct stresses/strains that act on the blade, ii) boundary conditions, iii) layers thickness, stiffness and iv) the current damage state captured by an initial debond length parameter.

Motivation
The purpose of this study is to demonstrate a methodology for predicting damage growth at trailing edge joint of the wind turbine rotor blade.A simplified one-dimensional analytical model is presented that requires input parameters such as laminate and adhesive layers thickness, size of the initial debonding, weight distribution along the blade and elastic, fracture mechanics properties of the laminate adhesive interface.Note that the model presented here is significantly simplified with the intention to be usable it in engineering practice.By identifying and measuring critical damage parameters, these simplified yet effective predictive tools enables the end-users to estimate the approximate time needed for next damage inspection or blade-repair actions.The accuracy of the analytical model is assessed by comparing it's result with results from an accurate 3D FE model.
The paper is organized as follows: the mechanics behind the trailing edge debonding, the damage geometry and the expression for energy release rate calculation is presented in section 2. Details of the finite element modelling and the boundary conditions are presented in section 3.In section 4, results from the analytical model and finite element simulations are discussed.The predictions of debond crack-tip energy release rate as a function of the size of the initial debond were provided for the stresses acting on DTU's 86 m blade.The limitations of the damage model are discussed in section 5 and the conclusions are summarized in section 6.

Analytical modelling 2.1. Stresses due to gravitational load
The wind turbine rotor blade acts like a cantilever beam with its root fixed and tip left free.Thus, the rotor blade is subjected to gravity load due its own mass resulting in edgewise stresses.Let x denote the coordinate that gives the position at any point along the blade, i.e., the distance measured from the root section.The gravity load distribution along the blade length, q(x) (N/m), can be obtained from the density of the blade material, ρ mat , gravitational constant, g (= 9.81 m/s 2 ) and the cross-sectional area of the rotor blade, A blade (x), given as subsequently, the edgewise moment per unit width can be obtained by integrating the gravity load from an arbitrary section x until the blade tip x = L, given as Following Mikkelsen [9], we assume a simplified cross-sectional area of the blade model wherein, the analysis of only load-carrying laminates located at spar caps (pressure side and suction side), leading and trailing edges, are taken into account.Since these laminates contribute to the overall stiffness of the blade, the simplified cross-sectional area can be written as A blade (x) = 2(A f lap (x) + A edge (x)) [9].Further, we assume that the laminates at spar caps take up the bending moments in the flap-wise direction while those at leading and trailing edges take up edgewise moments.From the expression of gravity load distribution, the shear force, T (x), and the edgewise moment, M 0 (x) over the blade cross section can be found by solving the two differential equations, with the boundary conditions T (L) = 0 and M 0 (L) = 0 that represent a free blade tip.The two differential equations are obtained from a force and a moment equilibrium of an infinitesimal beam element, respectively.Using the trapezoidal rule [10], the following numerical scheme is implemented with a blade split up into n elements with the nodal points: Let σ 0 (x) denote the in-plane edgewise stress acting on the trailing edge, at a given distance x from the root and with the simplified area moment of inertia [9], where h edge (x) is chord length, the edgewise stresses can be written as [9] σ 0 (x) = M 0 (x) h edge (x)A edge (x) (7)

Buckling-induced debonding
Consider a segment of length, l, along the trailing edge of a wind turbine rotor blade where there exists a local buckling driven debonding as shown in Fig. 2b.In the present study, the trailing edge consists of top and bottom laminates each with equal thickness, h 1 , and a central adhesive layer with thickness h 2 .We assume that both the laminates are made up of same material with a Young's modulus, E 1 , and Poisson's ratio, ν 1 , while the adhesive properties are denoted by E 2 and ν 2 .We define two non-dimensional parameters Σ = E 1 /E 2 and η = h 1 /h 2 .Fig. 3 shows the equilibrium of the half symmetric sandwich panel before and after buckling.If the right end of the laminate is subjected to a uniform in-plane compressive edgewise stress, σ0 (x), and a moment per unit width, M 0 .By assuming same axial strain in all the layers and from the force balance in the x 1 -direction, the compressive stresses in the top laminate is given by If this stress exceeds the critical buckling stress, σ 1,cr , then the laminate buckles locally and results in a strip of debond length 2a, see figure Fig. 2b.We assume h 1 << a so that the thin ply can be treated as a long, clamped Euler column of length 2a.Further, it is assumed that outside the debonded area the adhesive and the bottom laminate are perfectly bonded and provide a stiff support to the top laminate h 1 .This assumption allows us to use the bilayer solution developed by Hutchinson and Suo [11].The critical buckling stress based on clamped-clamped Euler column of length 2a is given by  where D is the flexural rigidity of the column per unit width, given by D = E 1 h 3 1 /(12(1 − ν 2 1 )).In case of buckling, a moment per unit width, M 1 , acts on the laminate and a change in the axial compression force ∆P per unit width near the crack-tip interfaces, see figure Fig. 3b.Here, the laminate is under axial compression and positive M 1 , meaning that the debond is open and the laminate deforms away from the adhesive layer.Using the von Karman non-linear plate theory, the maximum deflection, the buckling load and the moment are obtained as From the equilibrium condition of post buckling (see Fig. 3b), the relation between the moments can be written as In order to calculate the energy release rate per unit width at the crack-tip, we use the expression from Karihaloo and Stang [12], given as IOP Publishing doi:10.1088/1757-899X/1293/1/012037 Figure 3: Equilibrium of the half symmetric trailing edge damage geometry a) before buckling of the top laminate and b) after buckling.

Numerical modelling
A 5 m prismatic beam model from a 17 m long Vindeby blade is modelled through the finite element software Abaqus 2021.The geometrical model is reconstructed according to the crosssection geometries measured at approx.7 m from the root.The detailed dimensions are shown in Fig. 4a.A straight through-width pre-crack is introduced in the adhesive trailing edge joint.The pre-crack is located at the middle section of the trailing edge with a length of 0.5 m, see Fig. 4b.Solid 8-node linear brick elements (Abaqus type C3D8R) were used in the model.The typical mesh sizes are 25 × 22 × 6 mm for spar cap regions, 25 × 11 × 14 mm for sandwich panels, and 25 × 4 × 0.8 mm for trailing edge bondlines.The mesh size is chosen based on mesh sensitivity analysis performed in [8,13,14].The mesh sizes result in a total element number of 162, 600.The energy release rate at the crack tips were computed using the virtual crack closure technique (VCCT) in Abaqus.
A rigid body reference node is used to represent the node set at the end of the blade section model.The reference node is set at the elastic centre of the cross-section.The model is loaded in displacement control.A compression load is created by prescribing a displacement u x to the left end reference node (in the x-direction).The right end of the model is fixed in all directions.Frictionless contact conditions are assigned to the two surfaces of the crack plane.

Gravitational loads, edgewise stresses and the energy release rate
To demonstrate the methodology, the example turbine considered is the 10 MW DTU reference turbine [16].Table 1 shows the blade details.The blade material analysed is a glass fibre Figure 4: A 5 m prismatic beam finite element model.A straight through-width crack is introduced in the middle section of the adhesive trailing edge joint (orange area marked in the blade section model).The pre-crack is 0.5 m long and the load and the boundary conditions are applied on the two reference nodes (shown as yellow coloured dots).The model is built from one of the measured dimensions of the Vindeby blade cut-out cross-section.The data can be found in [15].composite made up of unidirectional E-glass and epoxy resin.The internal shear forces T due to gravity load variation and the bending moment per unit width values at a given position x from the root section can be obtained numerically by solving the governing differential equations of the Bernoulli-Euler beam theory.Fig. 5 shows the variation of shear force on the 10 MW DTU reference blade [16].Here, the total mass includes the mass of the load-carrying GFRP laminates and also the mass of the aerodynamic shell and the shear web.
The edgewise stresses along the blade can be found from the moment M edge (x) working on a cross-section using the eq.( 4).Fig. 6 shows the variation of edgewise stresses calculated for Figure 5: The variation in the shear force due to gravity load acting on the 10 MW DTU wind turbine rotor blade [16].
the 10 MW DTU reference blade [16].The maximum value is at the root section and zero at the tip.
Figure 6: The distribution of edgewise stresses acting on the 10 MW DTU wind turbine [16].

Damage growth prediction model
Depending on the location and the initial size of the debonding, the post-buckling response of the trailing edge laminate depends on the stress acting at the point, layers thicknesses and their Young's moduli.Fig. 7 illustrates the required input information classified into A) Structural blade model B) Damage parameters and C) Materials test data in order to predict the debond growth.The damage model we refer to is an excel spreadsheet tool that gives the output information as debond crack-tip energy release rate, for a fixed known initial size and number of cycles required to reach a predefined critical crack length.The energy release rate for several initial debond lengths in the range of 0.1 m -1 m and for various stress levels in the 10 MW DTU reference blade is shown in Fig. 8.The results for G are only valid as long as the debond buckles (σ 0 ≥ σ 1,cr ).In case of no buckling, i.e., σ 0 < σ 1,cr , G = 0. From Fig. 8, it can be noticed that for low compressive stresses, 40 − 70 MPa, the minimum size of the debond length 2a must be greater than approximately 0.2 m in order to get buckling of the top laminate.While in case of high stresses 80 − 130 MPa, 2a should be greater than approximately 0.14 m to get buckling.It is seen that for fixed σ 0 , G increases with increasing a.
After buckling of the top laminate and consequently, with the initiation of debond crack along the interface of the top laminate and the adhesive layer, it is assumed that the debond crack-tip propagates further.Additionally, assuming that the crack-tip grows in normal crack opening under the tension-compression cyclic loading, then from the interfacial properties of the bondline of E-glass fibre composites and epoxy resin, approximate crack growth time for the debond crack-tip to reach a critical length can be estimated.For illustration, a critical length of l c = 5 m is considered in the present study and the turbine rotational speed is assumed as 15 rpm.Firstly, from Goutianos and Sørensen [17]  as , where E is the Young's modulus and R is the load ratio.For the present problem, the debond only buckles when σ 1 ≥ σ 1,cr .When σ 1 < σ 1,cr , the debond crack remains closed.Then G = 0.This holds true also when the laminate is loaded in tension.Thus, effectively, R = 0. Using the Paris-Erdogan law, the crack growth rate is written as da/dN = C(∆K) m , here C and m are Paris-Erdogan law constants obtained from the materials interface characterization tests.We use the below equations to calculate the time in days/years where f is the loading frequency.Firstly, the step size ∆a is chosen with a fixed initial value and then 2a is calculated by adding ∆a iteratively.Since G (see eq.( 12)) is dependent on σ 1,cr which in turn is dependent on 2a, therefore, G is calculated in an incremental approach.Thus, for a given 2a i , the G max (2a i ) is obtained and subsequently, (da/dN ) i is estimated.Thereafter, from the current debond length, the number of cycles, N i and the time taken t i to grow the step size is calculated.For a given critical length l c , the time taken to reach this critical length, t c , is estimated.The procedure is repeated with several different ∆a values until a convergence in t c is obtained.
Table 2 provides the time estimation approximately in days and also in years for an initially detected debond crack-tip to propagate up to 1 m length.As an example to demonstrate debond growth rate predictions, two stress levels σ 0 = 40 MPa and 60 MPa are considered in Table 2.For each stress level, two different initially detected debond lengths 2a were chosen to show the time taken for the detected crack to reach 1 m length.It can be seen that for 60 MPa stress level, when 2a varies from 0.2 m to 0.225 , the time considerably reduces from 3167 days to 2 days.Thus, it can be deduced that the higher the stress level, the earlier the detection of debond crack allows sufficient time to take preventive measures or repair actions.

Discussion
An important finding from Fig. 8 is the determination of critical energy release rate, G c , where an initially detected debond crack can grow in unstable manner.Thus, one can estimate the critical crack length l c for that specific stress level in the blade beyond which the detected debond crack-tip propagates .In Fig. 8, dashed lines are indicated for the case σ 0 = 130 MPa to denote G c = 3520 J/m 2 and l c = 0.35 m.Therefore, for any initial debond crack detected along the blade during the inspection and for a known stress level at that position, the limiting crack length and the approximate time to reach that limit crack length can be estimated.

Limitations of the model i)
The damage framework presented is a 1D problem which is based on Euler's theory of column buckling.The analysis takes into account only the thickness of the layers and does not include width and its variations in each layer along the cross-section of the rotor blade.Therefore, an extended framework of this model would be a 2D plate buckling analysis problem.
iii) A constant thickness in layers (laminate and adhesive) is assumed all along the blade length.ii) Constraining effects from the sandwich/core panel are not included.
iii) The analysis is based on only in-plane loads coming from the self-weight/gravity of the blade, the other external loads such as out-of-plane loads are not taken into account.iv) Since G increases with increasing debond length a, static crack growth will be unstable.
When J ext reaches G c , this could be used as a criterion for determining critical length a c .

Concluding remarks
Buckling driven debonding in the trailing edge of a wind turbine rotor blade is assessed using an analytical model that takes input from the current damage state.The parameters required are a) size of the initial debonding, b) laminate and adhesive layer thicknesses, c) elastic and interfacial fracture mechanics properties and d) weight distribution of the material along the rotor blade.The proposed approach has been demonstrated on 10 MW DTU reference turbine IOP Publishing doi:10.1088/1757-899X/1293/1/012037 12 by extracting gravitational loads and edgewise stresses.The energy release rate predictions have been provided for several initial debond lengths and for various stress levels acting along the rotor blade.An approximate time estimation is provided for the debond crack-tip to reach a (predefined) critical length.The methodology demonstrated in the current study can be used in a damage specific predictive tools that helps in the prognosis of the service life of the wind turbine rotor blades.

Figure 1 :
Figure 1: A wind turbine rotor blade a) with point A on the trailing edge shown to indicate how gravity induced stresses (edgewise loading) changes during its complete revolution and b) trailing edge damage seen in field service, taken from [2].

Figure 2 :
Figure 2: Schematic illustration of the damage geometry a) wind turbine rotor blade with trailing edge debonding caused due to local buckling b) one-dimensional representation of the damage geometry with layers thickness and size of the initial debonding.

Figure 7 :
Figure 7: A block diagram showing various fields of input information that goes into the trailing edge damage growth prediction model.

Figure 8 :
Figure 8: The energy release rate, G, shown as a function of the initial debond length for various stress levels that are dependent on the position, x, along the blade length.The steps for determining critical length l c for the case of σ 0 = 130 MPa are denoted by dashed lines.

Table 1 :
Data of DTU 10 MW reference turbine , we have the relation between ∆K and G max

Table 2 :
Illustration of trailing edge damage growth from the current state of a given initial debond length and a stress level