A computational framework for the lifetime prediction of vertical-axis wind turbines

A novel computational framework is presented for the lifetime prediction of vertical-axis wind turbines (VAWTs). The framework uses high-fidelity computational fluid dynamics (CFD) simulations for the accurate determination of the aerodynamic loading on the wind turbine, and includes these loading characteristics in a detailed 3D finite element method (FEM) model to predict fatigue cracking in the structure with a fatigue interface damage model. The fatigue interface damage model allows to simulate high-cycle fatigue cracking processes in the wind turbine in an accurate and robust fashion at manageable computational cost. The FEM analyses show that the blade-strut connection is the most critical structural part for the fatigue life of the VAWT, particularly when it is carried out as an adhesive connection (instead of a welded connection). The sensitivity of the fatigue response of the VAWT to specific static and fatigue modeling parameters and to the presence of a structural flaw is analyzed. Depending on the flaw size and flaw location, the fatigue life of the VAWT can decrease by 25%. Additionally, the decrease of the fatigue resistance of the VAWT appears to be mainly characterized by the monotonic reduction of the tensile strength of the adhesive blade-strut connection, rather than by the reduction of its mode I toughness, such that fatigue cracking develops in a brittle fashion under a relatively small crack opening. It is emphasized that the present computational framework is generic; it can also be applied for analyzing the fatigue performance of other rotating machinery subjected to fluid– structure interaction, such as horizontal-axis wind turbines, steam turbine generators and multistage pumps and compressors.


Introduction
Vertical-axis wind turbines (VAWTs) have recently gained increasing interest for applications in offshore and urban environment (Paulsen et al., 2014;Islam et al., 2013;Gsänger and Pitteloud, 2015;Blocken, 2014).VAWTs have a number of advantages in comparison to the more conventional, horizontal-axis wind turbines (HAWTs), namely that (i) they generate less noise due to a lower peripheral velocity of the rotor blades, (ii) they are able to harness energy from any wind direction (i.e., omni-directional wind capture), (iii) their mechanical complexity is less, leading to low maintenance that primarily takes place at ground level, (iv) their installation goes quicker, and (v) they can be placed relatively close together in wind farms (Molina et al., 2017;Rezaeiha et al., 2018a;Tavallaeinejad et al., 2022).Although VAWTs and HAWTs have about the same ideal efficiency when both designs are based on the aerodynamic lift principle (Thönniszen et al., 2016) (which approximately equals 59% of the kinetic wind power in accordance with the Betz limit), the commercial success of VAWTs is still less.This is partly due to the fact that in the past two to three decades VAWTs have not received the same amount of research and development effort as HAWTs.Additionally, the lifetime of VAWTs is relatively limited (Cheng et al., 2017); for some applications blades have been observed to fail by structural fatigue as little as two to three years after installation (Paraschivoiu, 2002).
The fatigue loading experienced by VAWTs is characterized by aerodynamic fluctuations generated under a continuously varying wind angle of attack.The complex flow phenomena induced around the blades include dynamic stall under strong vortex shedding at the suction side of the blade (Ferreira et al., 2009;Yu et al., 2017;Geng et al., 2018;Rezaeiha et al., 2018a;Almohammadi et al., 2015), blade-wake interactions due to the blades encountering their own wakes (Amet et al., 2009;Rezaeiha et al., 2018a), and flow curvatures caused by the rotational motion of the blades (van der Horst et al., 2016;Bianchini et al., 2016).The combination of these flow phenomena makes the accurate prediction of aerodynamic loads on VAWTs a challenging task, https://doi.org/10.1016/j.ijsolstr.2023.112504Received 19 May 2023; Received in revised form 29 August 2023; Accepted 2 October 2023 for which various modeling strategies have been proposed.First, there are the multiple streamtube models (Paraschivoiu, 1988;Bangga et al., 2019), whereby the swept volume of the rotor is divided into independent upstream and downstream streamtubes.The induced velocity of the streamtubes is calculated by equating the drag force on the rotor to the corresponding change in fluid momentum.Despite their usefulness in the design phase, multiple streamtube algorithms fail to predict the local aerodynamic loads on the blades with high accuracy, due to difficulties in accounting for the dynamics of the wake and neglecting the blade-wake interactions (Simão Ferreira et al., 2014;Ayati et al., 2019).Second, vortex models have been used for the simulation of the flow field around VAWTs.These models realistically account for blade-wake interactions and vortex stretching and contraction (Cardona, 1984;Zanon et al., 2015), but are meant for the analysis of relatively simple (2D) aerodynamic configurations.Further, with this approach flow separation processes cannot be adequately simulated, and it is generally challenging to mimic viscous effects (Tescione et al., 2016;Dixon et al., 2008;Roy and Bandyopadhyay, 2006).Third, computational fluid dynamics (CFD) models have been applied for a detailed assessment of the aerodynamic performance of VAWTs (Rezaeiha et al., 2017a,b;He et al., 2020;Lin et al., 2019b;Lei et al., 2017a,b;Geng et al., 2018;McNaughton et al., 2014).Although from the above-mentioned modeling strategies CFD modeling is the most computationally demanding, it allows to realistically simulate all features of the unsteady airflow around the turbine, and is able to accurately determine the pressure and shear stress profiles on the turbine blades.Since the reliability of CFD simulations of VAWTs is very sensitive to the computational settings, it is imperative to carry out solution verification and validation studies in order to rigorously establish the modeling accuracy (Rezaeiha et al., 2018b).When this requirement is met, CFD modeling is a promising tool for a detailed calculation of the aerodynamic loading profiles that are underlying the fatigue response of VAWTs (Rezaeiha et al., 2017b(Rezaeiha et al., , 2018b(Rezaeiha et al., , 2019;;Geng et al., 2018).
The analysis of the fatigue behavior of VAWTs is generally considered as complex and computationally expensive, which may be the reasons that to date it has only received limited attention in the literature.In Lin et al. (2019a), the fatigue life of a laminated VAWT blade under a representative wind loading is computed from the straintime history obtained from an elastic finite element method (FEM) analysis.The rain-flow counting method is applied to the strain-time history by dividing the computed strain into distinct ranges.The overall, accumulated fatigue damage of the blade is subsequently calculated over these strain ranges using the well-known Palmgren-Miner's rule.This approach, which has also been applied for HAWTs (Rezaeiha et al., 2017d;Cheng et al., 2017), only serves as a rough approximation of the fatigue life, as it leaves the specific growth characteristics of fatigue cracks out of consideration.
An FEM analysis of fatigue crack growth in VAWTs may be performed by applying the well-known Paris law (Paris et al., 1961), whereby the crack driving force at the crack tip, which equals the stress intensity factor, determines the cyclic evolution of the crack, see Suiker and Fleck (2006), Turon et al. (2007), Pirondi and Moroni (2010) and Harper and Hallett (2010).It is emphasized, however, that this crack driving force can only be rigorously established for a single crack with a well-defined crack tip, and thus cannot be adequately determined in the case of complicated fracture patterns characterized by multiple crack bifurcation and crack coalescence events.In addition, Paris law has been originally developed for estimating the structural fatigue life from the analysis of the intermediate, crack propagation phase, thus ignoring the initial, crack nucleation phase and the final, catastrophic failure phase.
Alternatively, fatigue cracking has been simulated within an FEM setting by applying load-cycle driven cohesive zone models (Khoramishad et al., 2011(Khoramishad et al., , 2010b;;Geng and Suiker, 2019), which mimic the local deformation in a cohesive material point near the crack tip by means of a so-called traction-separation law.Cohesive zone models describe the nucleation, propagation and catastrophic failure phases of a crack in a natural fashion, both for interfacial cracks between layers (Cid Alfaro et al., 2009;Forschelen et al., 2016) and for cracks in bulk materials (Cid Alfaro et al., 2010a;Luimes et al., 2018;Eumelen et al., 2019;Scheperboer et al., 2019;Luimes and Suiker, 2021).Due to the local nature of the constitutive formulation, the approach enables to analyze complex fracture patterns -with possible crack bifurcation and crack coalescence events -in an accurate and robust fashion (Cid Alfaro et al., 2010b,a;Eumelen et al., 2019;Luimes and Suiker, 2021;Eumelen et al., 2023).
In the present communication an advanced computational framework is presented for the lifetime prediction of a VAWT.The modeling accuracy of the computational framework is guaranteed by computing the aerodynamic loading characteristics on the wind turbine with a high-fidelity CFD model, and including these loading characteristics in a detailed 3D FEM simulation that carefully mimics structural fatigue cracking by means of a load-cycle driven cohesive zone model.For this purpose, the fatigue interface damage model developed in Geng and Suiker (2019) is used, which allows to simulate high-cycle fatigue cracking processes in a robust and accurate fashion at manageable computational cost.The FEM analyses consider the sensitivity of the fatigue response of the VAWT to specific static and fatigue modeling parameters.Further, a variation study is performed on the effect of structural flaws, which, depending on their size and location, may considerably reduce the fatigue life of the wind turbine.The modeling study particularly focuses on the fatigue resistance of an adhesive bonding at the fatigue-sensitive blade-strut connection of the VAWT.Since the manufacturing of adhesive bonds occurs at relatively low temperature, strength reductions and deformations typical of welded connections can be avoided, which makes adhesive connections potentially interesting for being applied in VAWTs.The type of adhesive selected in the fatigue analyses is the toughened epoxy extensively tested in Sugiman et al. (2013a,b) and Khoramishad et al. (2010a) on its static and cyclic bonding capacity for aluminium plates, for which the parameters of the fatigue model applied in the present study have been accurately calibrated in Geng and Suiker (2019).
The paper is organized as follows.Section 2 discusses the highfidelity CFD simulation performed for determining the aerodynamic loading used as input in the FEM fatigue model of the VAWT.Section 3 reviews the interface damage model employed for simulating high-cycle fatigue cracking in an FEM setting, whereby the model formulation is specialized towards a fatigue analysis of the VAWT.In Section 4 the FEM fatigue analyses of the VAWT are discussed.The critical fatigue loading conditions are identified from the CFD results, after which the FEM results obtained for a representative, reference case are analyzed in detail.The sensitivity of the fatigue response of the VAWT to specific static and fatigue modeling parameters is assessed, and the reduction of the fatigue life on the size and location of a structural flaw is explored.Section 5 summarizes the main conclusions of the study.

CFD simulation of a rotating VAWT
The computational framework for the lifetime prediction of verticalaxis wind turbines is composed of three successive steps.First, the aerodynamic wind flow pattern around the turbine is accurately computed during one complete turbine revolution via high-fidelity CFD simulations.Second, the aerodynamic loading profiles acting on the wind turbine blades are extracted from the CFD results per degree of turbine rotation, thus resulting in 360 aerodynamic loading profiles for a complete turbine revolution.Together with the centrifugal and gravitational loadings, all 360 aerodynamic loading profiles are successively imposed on a linear, elastic FEM model of the wind turbine.From the results of the 360 FEM simulations, the CFD loading profile that is most critical for fatigue cracking is identified.This is done by selecting the FEM response corresponding to the turbine rotation angle at which the local elastic tensile stress (weighted by the corresponding tensile strength) in the turbine is maximal.The location of the maximum elastic tensile stress defines the so-called ''most fatigue-sensitive area'' of the VAWT, at which the minimum value of the elastic tensile stress generated during the turbine revolution is also identified from the simulation results, in order to compute the effective, cyclic stress amplitude characterizing the turbine revolution.Third, the above critical loading profile is used for determining the input parameters for a loadcycle driven fatigue model, as developed in a previous work (Geng and Suiker, 2019), whereby one load cycle corresponds to one turbine revolution.The fatigue interface damage model simulates progressive fracture under a large number of load cycles, i.e., high-cycle fatigue processes, and is used in a non-linear FEM fatigue simulation of the wind turbine.Accordingly, the wind turbine is subjected to the most critical loading profile observed during a turbine revolution, and its lifetime is next computed by incrementally increasing the number of load cycles (or turbine revolutions).The lifetime of the wind turbine is set by catastrophic fatigue fracture (or structural failure) under a large number of load cycles.
In correspondence with the first step explained above, this section discusses the computation of the aerodynamic loading on the VAWT during one complete turbine revolution, as following from high-fidelity CFD simulations.Accordingly, Section 2.1 discusses the geometrical and operational characteristics of the VAWT.In Section 2.2, the computational settings and parameters of the CFD simulation are presented, and Section 2.3 treats the simulation results.

Geometrical and operational characteristics of the VAWT
The VAWT considered in this study is a domestic H-type (Darrieus) turbine composed of two blades, four struts, and a turbine tower, as shown in Fig. 1.This type of VAWT has been extensively studied in Rezaeiha et al. (2018a), and the geometrical and operational characteristics listed in Table 1 have been taken from a representative case study considered in this reference.The turbine has a rotor diameter  = 1 m and a blade height  = 1 m, leading to a swept area  s =  = 1 m 2 .The blades and struts are made of 6061-T6 aluminum, which is characterized by a medium to high strength, a low weight and good corrosion resistance and weldability (Bauccio, 1993).The wind turbine tower is made of stainless steel 304.The blade profile corresponds to a hollow NACA0018 symmetrical airfoil, which has a thickness of 0.8 mm, two inner spars of thickness 1.3 mm, a chord length  b = 60 mm, and a maximum height of 0.18 b = 11 mm, see Fig. 2. Accordingly, the blade aspect ratio is ∕ b = 16.67, and the total blade area equals  b =  b  = 0.12 m 2 (with  the number of blades), in correspondence with a solidity  =  b ∕ s = 0.12.The two struts connecting each of the two blades to the turbine tower have a NACA0030 profile with a chord length  s = 23 mm.The VAWT operates at a rotational speed  = 46.5 rad∕s under an incoming freestream velocity  ∞ of 9.3 m∕s (as characteristic of urban wind conditions), resulting in a tip speed ratio  t = ∕ ∞ = 2.5, where  = ∕2 = 0.5 m is the rotation radius of the blade.The chord-based Reynolds number is computed as  c =  rel  b ∕ = 1.03 × 10 5 , with the relative velocity .1 m∕s (as evaluated at an azimuth angle  = 90 • , see Fig. 3) and the kinematic viscosity of air  = 1.46 × 10 −5 m 2 /s (measured at 13 • C).The approach-flow total turbulence intensity   = 5%, which is assumed to be representative of urban wind conditions in typical high wind-speed regions near buildings, such as corner streams or strong shear flows over roofs.Note that under spatially bounded wind conditions in aeronautical wind tunnels the turbulence intensity may lie considerably lower, in the order of one tenth of a percent (Lee and Gerontakos, 2004;Geng et al., 2018).The turbulence length scale is set as  =  = 1 m, and estimates the size of the large, energy-containing eddies in a turbulent flow at the inlet of the CFD simulation.In the case of wind flow around a turbine, the size of these eddies is found to be  on the order of the turbine diameter (Andersen et al., 2017).Finally, the reduced frequency equals  =  b ∕(2 rel ) = 0.06.
As described in Rezaeiha et al. (2018a), the tip speed ratio defines the degree of variation of the angle of attack on the blades during each turbine revolution, whereby the current, relatively low value of  t = 2.5 relates to a relatively large degree of variation of the angle of attack.Consequently, this may result in highly unsteady aerodynamic loads on the blades, especially if the critical angle of attack for dynamic stall is exceeded (Rezaeiha et al., 2018a;Geng et al., 2018), which is expected to have an adverse effect on the fatigue life of the blade.In other words, the current, relatively low value of the tip speed ratio may be considered as critical for the fatigue life of the VAWT, which is the reason it has been selected in the current modeling study.

Computational settings and parameters of CFD simulation
The aerodynamic loading on the rotor blades is calculated using a high-fidelity CFD simulation.Fig. 3 schematizes the 2D computational domain employed in the CFD analysis, which has also been used for the study in Rezaeiha et al. (2018a) and represents the cross-section  (Tescione et al., 2014;Hand and Cashman, 2017).The dimensions of the computational domain are 35 ×20, with  the rotor diameter, whereby the distances from the turbine center to the domain inlet and outlet are 10 and 25, respectively.The blockage ratio, which follows from the ratio between turbine diameter and domain height, is ∕(20) × 100% = 5%.The size of the computational domain is selected in line with the recently published guidelines for high-fidelity CFD simulations of VAWTs (Rezaeiha et al., 2019(Rezaeiha et al., , 2018b(Rezaeiha et al., , 2017a)).A comparison with the results from a comparable 2.5D computational domain shows that the selected 2D computational domain results in acceptably small overestimations of the power and thrust coefficients of less than 6% and 2%, respectively, see Rezaeiha et al. (2018a) for more details.The computational grid (or mesh) is depicted in Fig. 4, and consists of quadrilateral cells with respective maximum and average  + values of 3.8 and 1.8 on the airfoil and 1.4 and 1.0 on the shaft of the turbine tower.The total number of cells is approximately equal to 400,000.The grid discretizations at the leading edge and trailing edge of the turbine blade are shown in more detail in the insets of Fig. 4. The boundary conditions correspond to a uniform velocity inlet, zero gauge static pressure outlet, no-slip walls for the airfoil and the shaft of the turbine tower, symmetric side boundaries, and a sliding grid interface between the rotating core and the fixed, surrounding domain.
The diameter of the rotating core of the computational domain is 1.5, see Fig. 3.The CFD software package ANSYS Fluent 16.1 is employed for solving the incompressible unsteady Reynolds-averaged Navier-Stokes (URANS) equations in the computational domain, whereby a 2nd-order discretization approach is used in both space and time.The pressure-velocity coupling is warranted by the SIMPLE algorithm (AN-SYS, 2016).The flow pattern around a VAWT blade, as characterized here by a moderate chord-based Reynolds number  c = 1.03 × 10 5 , is simulated using the correlation-based transition SST turbulence model.This model combines the 2-equation SST  −  turbulence formulation (Menter, 1994) with two transport equations that respectively solve for the transition momentum thickness Reynolds number   and the intermittency  int (which determines the percentage of time the flow is turbulent, with  int = 0 for fully laminar and  int = 1 for fully turbulent), see Menter et al. (2005) for more details.Intermittencybased models are capable of providing adequate accuracy in solving the physics of laminar-to-turbulent transitional flows, and have been successfully applied for detailed calculations of the flow field around turbine blades (Geng et al., 2018;Rezaeiha et al., 2019;Suzen et al., 2003).
The CFD analysis is initialized by a steady RANS computation, after which a transient simulation of 20 turbine revolutions is carried out using an azimuth angle increment of 0.1 degree.The 20 turbine revolutions are performed to ensure that the results statistically reach a steady-state flow field.Here, 20 iterations per time step are used to warrant that scaled residuals drop below a prescribed threshold value of 1 × 10 −5 .Subsequently, starting from the 21st turbine revolution, the computational results are sampled over 100 turbine revolutions.In order to assess the accuracy of the URANS simulation, a verification study of the computational result has been performed that includes a grid-sensitivity analysis and sensitivity studies on the time-step size and convergence criterion.In addition, two sets of validation studies have been carried out whereby the CFD results were compared with the experimental measurements presented in Raciti Castelli et al. ( 2011) and Tescione et al. (2014).In general, a good agreement has been observed.More detailed information on the solution verification and validation studies can be found in Rezaeiha et al. (2017c,a) and Rezaeiha et al. (2018a).

CFD results
The aerodynamic loading acting on the wind turbine blade is composed of a static pressure  and a (skin) shear stress  w in, respectively, the normal and tangential directions of the blade surface.These two contributions can be conveniently expressed in dimensionless form by means of the so-called pressure coefficient  p and skin friction coefficient  f , i.e., where  ref = 0.5 2 ∞ is the dynamic pressure.b).Note that the coefficient  p along the vertical axes is multiplied by −1, so that positive and negative values relate to suction and compression, respectively.Clearly, the peak values of  p are much larger than those of  f , from which it may be concluded that the fatigue behavior of the blade is mainly driven by pressure loading.At the upper surface of the blade the pressure coefficient at  = 90 • strongly decreases in magnitude under an increasing chordwise position ∕ b .This negative pressure gradient is caused by the laminar-to-turbulent transition in the boundary layer, leading to turbulent flow separation towards the trailing edge of the blade (Rezaeiha et al., 2018a).Further, the magnitude of the pressure coefficient  p at  = 270 • is considerably lower than at  = 90 • , which is caused by blade-wake interactions between vortices shed from the blades in the upwind region, 45 • <  < 135 • , and the boundary layer on the blade passing downstream (Rezaeiha et al., 2018a).The resulting normal ( N ) and tangential ( T ) force coefficients can be determined by computing the integral of the coefficients  p and  f along the blade surface, i.e., (2)  where the superscripts ''l'' and ''u'' used for  p and  f refer to the lower and upper blade surfaces, respectively.From Eqs. ( 1) and ( 2) it follows that the overall normal and shear forces (per unit depth) on the blade can be computed as shows the force coefficients  N and  T as a function of the azimuth angle , as evaluated during the final (100th) turbine revolution.It becomes clear that the fatigue behavior of the blade will be governed by the relatively large loading in the normal direction of the blade chord line  b .In the upwind region 45 • <  < 135 • , the aerodynamic force coefficients obtain a maximum, followed by a strong drop in value due to dynamic stall under boundary layer separation (Rezaeiha et al., 2018a;McCroskey et al., 1976;Lee and Gerontakos, 2004;Geng et al., 2018;Wang et al., 2012).In the downwind region 225 • <  < 315 • fluctuations in the force coefficients can be observed, which, as mentioned above, are due to interactions between vortices shed from the blades in the upwind region and the boundary layer on the blade passing downstream, see Rezaeiha et al. (2018a) for more details on this phenomenon.

Fatigue interface damage model
In order to simulate high-cycle fatigue cracking processes at manageable computational cost, in Geng and Suiker (2019) an efficient fatigue interface damage model has been developed based on the net envelope of the sequence of loading-unloading cycles.In accordance with this approach, at the start of the simulation the amplitude of the fatigue load is imposed incrementally in a quasi-static fashion, after which it is kept constant and fatigue cracking is allowed to develop  further under the application of a sequence of (large) load cycle increments.These two distinctive modeling steps have been integrated in the fatigue interface damage model presented in Geng and Suiker (2019) by combining the fatigue evolution law proposed in Khoramishad et al. (2010aKhoramishad et al. ( , 2011) ) with the static interface damage model developed in Cid Alfaro et al. (2009).The static interface damage model is reviewed in Section 3.1, and its combination with the fatigue evolution law is explained in Section 3.2 within the context of a fatigue analysis of the VAWT.

Static interface damage model
The static interface damage model developed in Cid Alfaro et al. (2009) simulates the quasi-brittle fracture behavior in the small cohesive zone at the tip of a crack, whereby in a 3D geometrical setting the tractions   in a material point of the cohesive zone and the relative displacements   across the corresponding crack faces consist of three components:  ∈ {1, 2, 3}, with the numbers respectively denoting the normal direction and the two tangential directions of the crack.For convenience, the tangential directions are taken parallel (index '2') and transversely (index '3') to the direction of crack development.The tractions and relative crack face displacements are related by means of the following damage law (Cid Alfaro et al., 2009): (3) Here,  is the elastic stiffness (with a dimension of force × length −3 ),   is the Kronecker delta symbol, and  is the damage parameter, which is bounded as 0 ⩽  ⩽ 1, with  = 0 referring to the undamaged state of a material point in the cohesive zone, and  = 1 referring to the fully damaged state whereby the corresponding crack faces have completely separated.Damage values in between these two bounds, F. Geng et al.
0 <  < 1, reflect that in the cohesive zone, also regularly named the fracture process zone, the crack faces are locally connected by fibrils that are able to effectively transfer a stress lower than the fracture strength; this stress decreases towards zero when the fibrils successively break upon further crack opening, whereby the damage parameter thus grows towards unity.Hence, the damage parameter  may be interpreted as a relative measure for the amount of fibril breakage in the cohesive zone at the crack tip, as originally proposed in Dugdale (1960) and Barenblatt (1962).Further, under crack face contact the normal crack face displacement  1 is negative and the response in the interfacial material point is elastic, which is accounted for in Eq. ( 3) through the Macauley brackets, ⟨⟩ = ( + ||)∕2.
The damage development of an interfacial material point is described by (Cid Alfaro et al., 2009) which corresponds to the bilinear traction-separation law shown in Fig. 8. Here, the ultimate traction  u determines the critical stress at which the crack nucleates, and the area under the traction-separation law equals the toughness  c =  u  u ∕2 =  0  u ∕2 that characterizes the propagation characteristics of the crack, with  0 and  u the equivalent, relative displacements of the crack faces at which damage is initiated ( = 0) and completed ( = 1), respectively.The values of  0 and  u , and thus the toughness value   , generally depend on the mode-mixity of the fracture process, see Cid Alfaro et al. (2009) for the specific expressions.These expressions allow to distinguish between the mode I (tension), mode II (in-plane shear), and mode III (out-of-plane shear) contributions to the fracture process.Further,  is a deformation history variable, which monotonically increases since damage is an irreversible process.The loading and unloading conditions are governed by a damage loading function  = F (, ) based on Eq. ( 4), which has been included in a rate-dependent kinetic law to describe the evolution of the damage parameter  as (Cid Alfaro et al., 2009): where Ḋ is the damage rate,  = || = √  2 1 +  2 2 +  2 3 is the effective, relative displacement across the crack faces and  is a relaxation parameter (with dimension of time).The upper expression in Eq. ( 5) reflects the damage rate when the effective deformation  exceeds the deformation history , whereas the lower expression sets the damage rate equal to zero when (i) the deformation history has not (yet) been reached, (ii) the interfacial material point is in a state of unloading, or (iii) the damage process has completed.The accuracy and robustness of the above static interface damage model have been demonstrated for practical applications related to cracking in a variety of materials, including metals (Cid Alfaro et al., 2009;Forschelen et al., 2016), historical paints (Eumelen et al., 2019(Eumelen et al., , 2020(Eumelen et al., , 2023)), wood (Luimes et al., 2018;Scheperboer et al., 2019;Luimes and Suiker, 2021), fibrous composites (Cid Alfaro et al., 2010b,a), and cementitious materials (Scheperboer et al., 2021;Luimes et al., 2022).

Fatigue evolution law for a VAWT
The static interface damage model outlined in the above section is combined with a fatigue evolution law, in which the fatigue amplitude during a load cycle increment is kept constant and crack development in the cohesive zone is effectively driven by load cycle increments .Such a modeling strategy is suitable for fatigue processes characterized by a large number of load cycles, i.e., high-cycle fatigue processes, for which the number of load cycles  may be interpreted as a time evolution parameter and the damage generated during a load cycle increment  evolves in an almost continuous fashion.Note that this 8. Bilinear traction-separation law adopted from Cid Alfaro et al. (2009).The initial elastic stiffness is represented by , the ultimate traction equals  u , and the toughness is  c (which equals the grey area under the traction-separation law).Damage starts when the effective relative crack face displacement  0 is reached, and is complete when  u is reached.The unloading branch with a reduced stiffness (1 − ) due to damage development is indicated by the red line, whereby the corresponding deformation history equals .
modeling approach is computationally efficient, as large computational costs associated to the calculation of the complete loading-unloading hysteresis during each individual load cycle are avoided.
At the onset of the fatigue simulation, the maximum amplitude  max of the fatigue loading is incrementally applied in a quasi-static manner, during which the amount of damage  generated in the interfacial material point is computed with the above static interface damage model.The fatigue process is next activated through the stepwise application of load cycle increments , whereby the fatigue damage parameter  f (with 0 ≤  f ≤ 1) in a material point of the cohesive zone evolves incrementally, in accordance with the following fatigue evolution law (Khoramishad et al., 2010a(Khoramishad et al., , 2011;;Geng and Suiker, 2019): Here,  f is the incremental fatigue damage,  is a load correction factor -to be defined below -and ,  and  are material constants that need to be calibrated from fatigue tests.The parameter  may be associated to the brittleness of the crack, i.e., the smaller the value of , the more brittle a fatigue crack is (Khoramishad et al., 2010a;Geng and Suiker, 2019).Hence, a higher value of  leads to a smaller amount of fatigue damage (since 0 <  ≤ 1), as a result of which the fatigue life becomes larger.Further, similar as in the static interface damage model,  = || = √  2 1 +  2 2 +  2 3 is the effective relative crack face displacement, which induces fatigue damage in the cohesive zone if its value exceeds the deformation threshold  th , see Eq. ( 6).The load cycle increment  is related to the rotational speed  of the wind turbine (in rad/s) as in which  is the time increment (in seconds) and  is the downtime factor, i.e., the relative time the wind turbine is not operational, with 0 ≤  < 1.The load correction factor  in Eq. ( 6) depends on the maximum amplitude  max of the fatigue loading, the ultimate static strength   of the wind turbine, and the ratio  (with −∞ <  < 1) between the minimum and maximum load (or stress) amplitudes of the fatigue loading (Khoramishad et al., 2010a(Khoramishad et al., , 2011;;Geng and Suiker, 2019): with  a material constant.Clearly, the values of the loading parameters in Eq. ( 8) need to be determined prior to the fatigue analysis of the VAWT, which occurs as follows.Essentially, the loading parameters are determined by the part of the turbine that is most susceptible to fatigue fracture, which may be either located (i) within the blade, or (ii) at the blade-strut interface.For identifying the critical structural part, by means of a CFD analysis the aerodynamic stress distribution along the outer circumference of the blade is calculated for one complete turbine revolution in 360 steps, whereby at each step the azimuth angle is incremented by one degree.As a particular example, Figs. 5 and 6 respectively show the pressure coefficient ( p ) distribution and friction coefficient ( f ) distribution determined from a CFD analysis for selected azimuth angles  = 0 • , 90 • , 180 • and 270 • , which, after substituting these results into Eqs.(1) 1 and (1) 2 , respectively provide the corresponding pressure () and skin shear stress ( w ) distributions along the blade circumference.These aerodynamic stress profiles are next applied on an elastic FEM model of the VAWT, which for each of the 360 azimuth angles results in the elastic stress field within the blade structure and at the blade-strut interface.From the results of all 360 FEM simulations, the maximum relative stress amplitude  R,max within the VAWT is computed as where the integer  indicates the specific azimuth angle considered in the FEM simulation,  1,max and  tb respectively denote the maximum principal (tensile) stress and ultimate tensile strength in the most fatigue-sensitive area of the blade, and  1,max and  u 1 respectively designate the maximum tensile traction and the tensile strength in the most fatigue-sensitive area of the blade-strut interface.Clearly, the ''most fatigue-sensitive area'' denotes a local area in the VAWT at which the tensile stress over one turbine revolution is maximal, and fatigue cracking thus may be expected to nucleate.Since for arbitrary azimuth angles the blade is predominantly loaded in the normal direction of the blade chord line, see Fig. 7, at the blade-strut interface the normal (tensile) traction has the largest influence on the fatigue life of the wind turbine, such that the effect by the tangential (shear) traction can be ignored in Eq. ( 9).In the most fatigue-sensitive area of the VAWT the minimum value of the elastic tensile stress over one turbine revolution is also determined, which, after dividing this value by the corresponding tensile strength, leads to the minimum relative stress amplitude  R,min .Essentially, the stress amplitudes  R,max and  R,min should be interpreted as local measures characterizing the effective fatigue load amplitude during one complete turbine revolution, and, as such, define the load ratio  in Eq. ( 8) as It is emphasized that the total load acting on the rotating VAWT is dominated by the constant centrifugal loading, supplemented with smaller contributions from the constant gravitational loading and the fluctuating wind loading.Specifically, for the wind turbine configurations analyzed in the present study, in the areas of large local stresses the centrifugal load typically determines about 80 to 90% of the total stress response.Hence, the difference between  R,max and  R,min , as caused by the fluctuating wind loading, will be moderate, i.e., the actual value of the load ratio  in Eq. ( 10) lies relatively close to unity.Analogous to the derivation of Eq. ( 10), from the definition given by Eq. ( 9), in Eq. ( 8) the ratio  max ∕ s between the maximum fatigue load (or stress) amplitude and the static failure load (or stress) follows as Inserting Eqs. ( 10) and (11) into Eq.( 8) then turns the load correction factor into 9. Reduced (fatigue) traction-separation law (indicated by the solid line), as obtained by applying Eq. ( 13) to the initial (static) traction-separation law shown in Fig. 8 (indicated by the dashed line).The strength  0 , toughness  0  and effective relative crack face displacements  00 and  u0 characterizing the initial (static) separation law are distinguished from the actual, reduced strength   , reduced toughness   , and effective relative crack face displacements  0 and  u by adding the superscript ''0''.
With Eq. ( 12), the fatigue damage  f can be updated after each load cycle increment  via Eq.( 6), whereby one load cycle thus corresponds to one revolution of the turbine blade.The updated fatigue damage parameter  f in a material point is subsequently used to linearly reduce the strength and toughness parameters defining the traction-separation law of the static interface damage model, see Fig. 9, in accordance with (Khoramishad et al., 2010a(Khoramishad et al., , 2011;;Geng and Suiker, 2019): Here,   reflects the three components of the strength and toughness parameters (i.e., tension, parallel shear, transverse shear), as respectively given by  u 1 ,  u 2 ,  u 3 and  I,c ,  II,c ,  III,c .Correspondingly, the initial values  0  of these parameters are  u0 1 ,  u0 2 ,  u0 3 and  0 I,c ,  0 II,c  0 III,c , and thus characterize the static interface damage model depicted in Fig. 8.Note from Fig. 9 that the initial values of the (mode-mix dependent) crack face separations  00 and  u0 at damage initiation and damage completion follow from the initial values of the strength and toughness as  00 =  u0 ∕ and  u0 = 2 0  ∕ u0 .From the reduced tractionseparation law illustrated in Fig. 9, the damage parameter  and the corresponding traction  and tangential material stiffness ∕ are updated at the local, material point level using the static interface damage model outlined by Eqs.(3) to (5), and in turn are sent to the global level of the FEM simulation to compute structural equilibrium and the corresponding relative crack face displacements  for the next iteration.This iterative procedure is repeated until the mechanical system has converged towards global equilibrium within a predefined tolerance.Subsequently, the above procedure is repeated under the application of the next load cycle increment , until the total number of load cycles  tot characterizing the fatigue process has been reached, or the structure has failed catastrophically under the fatigue loading, see Geng and Suiker (2019) for more details.Local fatigue failure at an interfacial material point is defined by the damage parameter  reaching unity, i.e.,  = 1.Note that the condition  = 1 is realized when (i) the fatigue damage parameter reaches unity,  f = 1, as a result of which, in accordance with Eq. ( 13), the bilinear tractionseparation law fully collapses, whereby the corresponding relative crack face displacement generally is less than the ultimate, relative crack face displacement under static failure  <  u0 , or when (ii) the tractionseparation law has only partly degraded (0 <  f < 1) and the relative crack face displacement has reached the ultimate relative crack face displacement in the traction-separation relation,  =  u0 .In principle, it depends on the material characteristics of the structure and the fatigue loading amplitude whether the first (brittle) or second (ductile) fracture mechanism occurs, although in practical cases fatigue fracture will often take place in a relatively brittle fashion (Suresh, 1998;Geng and Suiker, 2019).

FEM fatigue analyses of a VAWT
This section discusses detailed FEM fatigue analyses of the VAWT.In Section 4.1 the fatigue loading parameters are determined that serve as input for the fatigue interface damage model outlined in Section 3. Section 4.2 provides the discretization aspects and material parameters of the FEM model of the VAWT.In Section 4.3 the simulation results are presented for a reference case, which in Sections 4.4 and 4.5 is followed by a discussion of results from variation studies of, respectively, specific fatigue and static modeling parameters.Finally, in Section 4.6 the influence on the VAWT fatigue life of the size and location of an initial flaw is treated.

Representative loading parameters for fatigue interface damage model
For determining the loading parameters that are required as input for the fatigue interface damage model, see Eq. ( 9), the blade-strut part of the VAWT shown in Fig. 1 is discretized into a 3D elastic FEM model, see Fig. 10.The FEM simulations are carried out using the commercial software package ABAQUS Standard.1By making use of symmetry, the blade-strut FEM model is composed of two horizontal struts that are connected to a vertical blade.Fig. 2 shows that the other end of the struts in practice is connected to the central turbine tower of the VAWT, which has been omitted in the model as it is considered to be not critical for the fatigue life of the VAWT.The connection between a strut and the shaft of the central turbine tower is mimicked by constraining the displacements at this strut end in all three directions.For all 360 azimuth angles considered, the aerodynamic pressure () and shear stress ( w ) distributions on the turbine blades, Eqs.(1) 1 and (1) 2 , are extracted from the CFD simulations.Figs. 5 and 6 show that the corresponding aerodynamic coefficients  p and  f vary along the chordwise direction.In the spanwise direction these aerodynamic coefficients, and thus the corresponding pressure and shear stress distributions, are taken as uniform, which is a realistic assumption when aerodynamic blade tip effects can be ignored (Tescione et al., 2014).In addition, the inertial forces, which include the centrifugal loading -as determined by the rotational speed of  = 46.5 rad∕s, see Table 1 -and the gravitational loading, are imposed on the blade-strut model.
The material properties used in the FEM analysis are obtained from the ASM Material Data Sheet (Bauccio, 1993), with the Young's modulus of the 6061-T6 aluminium blades and struts taken as  = 68.9GPa, the Poisson's ratio as  = 0.33, and the density as  = 2700 kg∕m 3 .The reference finite element model consists of 64,486 elements in total.The blade is modeled with 37,686 4-node shell elements equipped with a (reduced) 1-point integration scheme, which have a thickness of 1.3 mm at the shear webs and girders and a thickness of 0.8 mm at the skin, in accordance with the blade cross-section shown in Fig. 2. The struts are modeled with 26,800 8-node hexahedral elements equipped with a (reduced) 1-point integration scheme.The length-to-width aspect ratio of the shell elements is 1.3.
Within the frames of a mesh-sensitivity study, the FEM results of this reference mesh are compared to those computed with a coarser mesh of 29,669 elements and a finer mesh of 142,378 elements, whereby the shell element aspect ratio in these meshes is taken the same as for the above reference mesh (i.e., equal to 1.3), and the aerodynamic loading on the blade is taken in correspondence with a blade azimuth angle of  = 0 • .Fig. 11 illustrates the results of the mesh sensitivity study, by depicting the total strain energy  of the blade-strut system (indicated in orange) and the maximum von Mises stress  eq (indicated in blue) for the coarse mesh (M1), the reference intermediate mesh (M2) and the fine mesh (M3).Note that the total strain energy and the maximum von Mises stress respectively serve as global and local indicators of the mesh accuracy.Observe that the difference in the maximum von Mises stress   eq computed with the meshes M3 (88.9 MPa) and M2 (85.1 MPa) is only 4%, and is considerably smaller than the difference of 15% between the meshes M3 and M1 (77.5 MPa).Further, the total strain energy  only shows a minor decrease under mesh refinement, and thus may be considered to be converged.In summary, the mesh sensitivity study illustrates that the reference mesh M2 combines a high accuracy with a manageable computational time, and therefore is an adequate choice for the elastic blade-strut simulations for a complete revolution of the turbine blade.
Fig. 12 shows the results of the 360 elastic FEM simulations corresponding to one turbine revolution,  ∈ [0 • , 1 • , … 359 • ], by plotting the relative maximum stress amplitudes   1,max ∕ tb and   1,max ∕ u 1 , with  ∈ [0, 1, … 359], in, respectively, the blade and the upper bladestrut interface.The reason for selecting the upper blade-strut interface, instead of the lower blade-strut interface (see Fig. 10), is because the The tensile strength of the blade-strut interface is taken in accordance with (i) a welded connection, and (ii) an adhesive connection, using a range of representative strength values for the weld (174 MPa ⩽  u 1 ⩽ 196 MPa, dark-gray area) and adhesive (40 MPa ⩽  u 1 ⩽ 80 MPa, light-gray area), as reported in the literature (Lin et al., 2016;da Silva and Adams, 2005).The tensile strength of the 6061-T6 aluminium blade is taken as  tb = 310 MPa (black line) (Bauccio, 1993).
upper blade-strut connection experiences a somewhat larger, and thus more critical, maximum tensile stress (i.e., about 10% larger).This is due to the gravitational loading, which in the blade induces an antisymmetric stress profile with respect to the half-height of the blade that breaks the symmetry in the overall stresses generated at the upper and lower blade-strut connections.Recall from Eq. ( 9) that from the maximum stress amplitudes mentioned above the value of  R,max can be determined, which occurs in correspondence with the following procedure.For the result plotted in Fig. 12, the tensile strength of the upper blade-strut interface is successively taken in accordance with (i) a welded connection, and (ii) an adhesive connection.The tensile strength of the welded connection relates to a 6061-T6 weld, which, depending on the type of welding method, lies within the range 174 MPa ⩽  u 1 ⩽ 196 MPa (Lin et al., 2016).For the adhesive blade-strut connection, the tensile strength at room temperature for a broad selection of adhesives lies within the range 40 MPa ⩽  u 1 ⩽ 80 MPa (da Silva and Adams, 2005).Further, the tensile strength of the 6061-T6 aluminium blade is  tb = 310 MPa (Bauccio, 1993).Observe that during one turbine revolution the adhesive blade-strut connection experiences the largest, and thus most critical relative (tensile) stresswhich thus determines the value of  R,max given by Eq. ( 9) -, followed by the welded blade-strut connection, and finally the blade itself.In the upper blade-strut connection, the largest, most critical tensile stress occurs at the bottom of the connection, relatively close to the leading edge, while in the blade itself the principal tensile stress is maximal slightly below its connection with the upper strut.Since the adhesive connection is the most critical component for the fatigue life of the VAWT, in the subsections below its behavior will be simulated in more detail by means of a dedicated fatigue model.Note that the tensile strength  u 1 =  u0 1 of the adhesive hereby should be selected substantially larger than the value corresponding to the static failure limit  R,max = 1, for which local cracking in the blade-strut connection already develops during the first load cycle.Accordingly, the blade-strut connection is assumed to be made from an FM 73 M OST toughened epoxy adhesive with a relatively high tensile strength of  u 1 = 70 MPa, which is known for its good static and fatigue resistance in the bonding of aluminium plates (Sugiman et al., 2013b,a;Khoramishad et al., 2010a;Geng and Suiker, 2019).
Independent of the value of , the maximum normal (tensile) traction occurs at the bottom of the blade-strut connection, relatively close to the leading edge.Together with the results depicted in Fig. 12, this leads to the conclusion that this location is the ''most fatigue-sensitive area'' of the VAWT.
Independent of the value of the azimuth angle , the maximum tensile traction   1,max occurs at the bottom of the upper blade-strut interface, close to the leading edge, which thus may be interpreted as the ''most fatigue-sensitive area'' as defined below Eq. ( 9).From the maximum tensile traction value of 45 MPa at an azimuth angle of  = 289 • , and the minimum value of 37 MPa at an azimuth angle of  = 62 • , together with the adhesive tensile strength  u 1 =  u0 1 = 70 MPa the maximum relative stress amplitude in Eq. ( 9) becomes  R,max = 0.643, and the minimum relative stress amplitude is obtained as  R,min = 0.529.This confirms that the maximum stress amplitude is (substantially) less than the adhesive tensile strength, i.e., 0.643 < 1, so that under the present loading conditions local, static failure will not occur, and the failure response of the bladestrut connection is fully determined by its long-term fatigue behavior.From Eq. ( 10), the values for  R,min and  R,max lead to a load ratio  = 0.82.Note that this value indeed lies relatively close to unity, as caused by the relatively large stress contribution following from the constant centrifugal force.With these values of  R,max and  R,min , the load correction factor in Eq. ( 12) becomes  = 0.135 when selecting the fatigue parameter in this equation as  = 2 (Geng and Suiker, 2019).This value of  will be used in the detailed fatigue analyses presented in the subsections below.
Although not illustrated in detail here, the shape of the profile of the maximum tensile stress   1,max in the aluminium blade under a varying azimuth angle  (as used for constructing the blade trend line in Fig. 12) is comparable to that of the blade-strut interface traction depicted in Fig. 13, with somewhat higher minimum and maximum stress values of 49.6 MPa and 57.4 MPa, respectively.Note that these stress values lie considerably below the yield strength of 6061-T6 aluminium,  y = 240 MPa, so that it is realistic to model the blades (and struts) as elastic in the fatigue analysis.

FEM discretization aspects and material parameters
From the results of the elastic FEM simulation presented in Fig. 12, it can be concluded that an adhesive connection between the blade and the strut may be considered as the most critical component for the fatigue life of the VAWT.Accordingly, a detailed 3D non-linear FEM fatigue analysis of the VAWT is performed, whereby fatigue cracking in the adhesive layer is modeled by interface elements endowed with the fatigue interface damage model outlined in Section 3. The fatigue interface damage model has been implemented in ABAQUS as a usersupplied subroutine (i.e., a UMAT), see Geng and Suiker (2019) for the F. Geng et al. details of the numerical implementation.The 8-node interface elements applied are equipped with a 4-point integration scheme, whereby the thickness of the interface elements is set equal to the thickness of the adhesive layer, i.e., 0.2 mm.In the fatigue model the mesh of the elastic FEM model described in Section 4.1 is further refined towards the blade-strut connection, in order to capture the local fatigue cracking response at the connection with high accuracy.The adhesive bladestrut connection is simulated by means of 1098 interface elements, and the total number of elements equals 131,518.A preliminary mesh sensitivity study has shown that the predicted fatigue lifetime of the VAWT changes by less than 1.5% when the above mesh is refined to a mesh whereby the number of interface elements is a factor of 6.5 larger and equal to 7208, from which it may be concluded that the current mesh is sufficiently fine for generating accurate numerical results.The boundary conditions, the centrifugal loading and the gravitational loading are taken the same as for the elastic blade-strut model, see Section 4.1.The aerodynamic pressure and shear stress distributions imposed along the blade circumference correspond to the CFD result for an azimuth angle of  = 289 • , at which the normal traction at the blade-strut interface during one turbine revolution is maximal, see Fig. 13.The non-linear equilibrium equations at system level are solved stepwisely by means of a fully implicit, incremental-iterative update procedure using an automatic time-stepping method.
The material properties of the aluminum 6061-T6 blade and strut components and the adhesive blade-strut connection are listed in Table 2.As mentioned, the adhesive is an FM 73 M OST toughened epoxy, which has proven to have a good fatigue resistance in the bonding of aluminium plates (Sugiman et al., 2013b;Khoramishad et al., 2010a;Geng and Suiker, 2019).The material parameters of the blade and strut are equal to those used in the elastic FEM model described in Section 4.1.Except for the values of  th ,  and  appearing in Eq. ( 6), the parameter values of the adhesive connection are taken from Geng and Suiker (2019), in which they were calibrated from experimental static and fatigue responses of a single lap joint of two aluminium plates as reported in Khoramishad et al. (2010a), Sugiman et al. (2013b) and Sugiman et al. (2013a).The elastic stiffness  = 1.15 × 10 4 N∕mm 3 of the interface elements is obtained from the ratio between the elastic modulus of the epoxy adhesive (2.3 GPa) and its thickness (0.2 mm).Further, for simplicity reasons the initial shear strengths and toughnesses in the directions parallel and transversely to the direction of crack development are taken the same, i.e.,  u0 2 =  u0 3 and  0 II,c =  0 III,c .The value of the relaxation parameter  in Eq. ( 5) is set relatively small, such that during the initial, quasi-static application of the fatigue load amplitude possible crack advancement occurs almost rate-independently.The deformation threshold  th that determines the onset of fatigue cracking in the blade-strut model, see Eq. ( 6), is calculated from the deformation threshold  SLJ,th calibrated from the single lap joint experiment reported in Geng and Suiker (2019), by assuming the following proportionality relation: ( th ∕ 1,max ) = ( th,SLJ ∕ 1,max,SLJ ).Here,  1,max and  1,max,SLJ are the maximum relative displacements in the normal direction of the fatigue crack in the blade-strut model and the single lap joint, respectively, as obtained from the corresponding elastic FEM analyses performed for determining the fatigue loading parameters in Eq. ( 6).At the critical azimuth angle  = 289 • indicated in Fig. 13 this results in  1,max = 3.36 × 10 −3 mm.Together with the values  th,SLJ = 6.00 × 10 −3 mm and  1,max,SLJ = 14.26 × 10 −3 mm calibrated for the single lap joint configuration (Geng and Suiker, 2019), the above proportionality relation leads to a deformation threshold of  th = 1.41 × 10 −3 mm.In addition, the value of the loading parameter  = 0.135 follows from the elastic blade-strut analysis presented in Section 4.1.Finally, the value of the fatigue parameter  has been selected as  = 1.10, but its influence on the fatigue response is known to be significant (Geng and Suiker, 2019) and therefore, together with the parameter  th , will be assessed in more detail in Section 4.4 by means of a parameter variation study.The fatigue model defined by the parameter values listed in Table 2 in the sequel will be referred to as the reference case.

FEM results for the reference case
Fig. 14 shows the simulation results of the fatigue analysis for the reference case defined in Table 2, by plotting the interfacial damage  in the adhesive at the upper blade-strut connection for three different numbers of turbine revolutions, namely  = 2.74 × 10 8 (Fig. 14(a)),  = 7.63 × 10 8 (Fig. 14(b)), and  = 8.21 × 10 8 (Fig. 14(c)).Here, the colors blue and red in the contour plots respectively denote the undamaged ( = 0) and completely damaged ( = 1) areas of the connection.It can be observed from Fig. 14(a) that at  = 2.74 × 10 8 load cycles interfacial damage initiates at the lower part of the leading edge (i.e., at the ''most fatigue-sensitive area'' of the VAWT) and, almost simultaneously, at the lower part of the trailing edge.The damage developing from these two nucleation sites subsequently coalesces, and grows towards the upper part of the leading edge, see Fig. 14(b), at which it finally leads to catastrophic failure of the complete bladestrut connection, see Fig. 14(c).From the angular frequency  = 46.5 rad/s (see Table 1) and the assumption of a downtime factor of the wind turbine of  = 0.5 (i.e., the wind turbine is not operational for half the time), with Eq. ( 7) the total number of turbine revolutions  f = 8.21 × 10 8 at catastrophic failure translates into a fatigue life of 7.03 years.Since the design life of a wind turbine is around 15 to 20 years, this fatigue life is considered to be somewhat limited, and therefore will be investigated further by a more detailed assessment of the sensitivity of the computational result to specific parameters of the fatigue model.This parameter variation study is presented in the section below, whereby, for the ease of interpretation, the fatigue life is expressed in ''years'' by adopting the above downtime factor of  = 0.5 for consistency.
In order to analyze the degradation process of the adhesive more closely, the cyclic evolutions of the damage parameters  and  f are evaluated at various locations along the strut chord line of the upper blade-strut connection depicted in Fig. 10.Fig. 15 depicts the cyclic evolutions of  (solid line) and  f (dashed line) in the integration points of 8 selected interface elements labeled ''E1'' to ''E8'', as designated by the black circles in the FEM mesh displayed in the inset of the figure.Observe that  f initially develops at a relatively low rate, which steadily increases towards the completion of damage,  = 1.The fatigue damage parameter  f first starts to grow in chord line element ''E8'' located close to the trailing edge, which occurs at  = 1.10 × 10 8 turbine revolutions.Subsequently, at  = 3.10 × 10 8 turbine revolutions the damage parameter  starts to increase in element ''E8''.The development of  lags behind that of  f because for the reference case defined in Table 2 the fatigue damage threshold  th is considerably  smaller than the (mode-mix dependent) static damage threshold value  00 depicted in Fig. 9.In other words, the fatigue damage parameter  f increases in value when the interfacial response is still in the initial, elastic range, so that the damage parameter  only starts to grow once the fatigue damage parameter  f has sufficiently reduced the initial tensile strength  u0 characterizing the static traction-separation law, in accordance with Eq. ( 13) and Fig. 9.After chord line element ''E8'', fatigue damage develops in chord line element ''E1'' near the leading edge, and then grows from both the trailing edge and the leading edge towards the central chord line element ''E3'', at which damage is completed ( = 1) at  = 8.20 × 10 8 turbine revolutions.
Observe further that at damage completion,  = 1, in all interface elements ''E1'' to ''E8'' the fatigue damage parameter has reached unity,  f = 1, so that in these elements the traction-separation law displayed in Fig. 9 has fully collapsed, whereby the corresponding relative crack face displacement is (much) less than the ultimate, relative crack face displacement under static failure,  <  u0 .Accordingly, fatigue fracture occurs in a brittle fashion (see also Section 3.2), in correspondence with a relatively short fracture process zone at the crack tip, which can be identified from Fig. 14 by the narrow color transition zone from blue (undamaged material,  = 0) to red (fully damaged material,  = 1).Fig. 16.Effect of fatigue parameter  (of the adhesive blade-strut connection) on the predicted fatigue life  f .The downtime factor of the VAWT is taken as  = 0.5.The reference case corresponds to  = 1.10, see also Table 2.

Variation study of fatigue model parameters 𝑝 and 𝜆 𝑡ℎ
From the results of the parameter variation study depicted in Figure 9 of Ref. Geng and Suiker (2019), it has been concluded that from the 4 fatigue parameters listed in Table 2 the parameter  typically has the largest influence on the fatigue response.The sensitivity of the predicted fatigue life  f of the VAWT on the parameter  characterizing the fatigue behavior of the adhesive blade-strut connection is illustrated in Fig. 16 for a selection of values  ∈ [1.00, 1.05, 1.10, 1.15] (adopting a downtime factor of  = 0.5).The other material parameters are taken as listed in Table 2. Clearly, a small increase in value from  = 1.00 to  = 1.15 greatly increases the fatigue life  f by a factor of almost 26, from 0.81 years to 20.64 years.This strong sensitivity confirms the importance of a careful calibration of  for obtaining an accurate lifetime prediction of the VAWT.Note further that the lifetime prediction of a VAWT with an adhesive (FM 73 M OST toughened epoxy) blade-strut connection exceeds the design life of 15 to 20 years typically required for domestic wind turbines if the exponent  has a value larger than  = 1.15.Fig. 17 shows the sensitivity of the predicted fatigue life  f on the deformation threshold  th for a selection of values  th ∈ [0, 5, 10, 15] × 10 −4 mm.It can be observed that the fatigue life over this range of values moderately increases by a factor of 3, from 2.14 years to 7.41 years.
F. Geng et al. Fig. 17.Effect of deformation threshold  th (of the adhesive blade-strut connection) on the predicted fatigue life  f .The downtime factor of the VAWT is taken as  = 0.5.For the reference case with  th = 1.41 × 10 −3 mm (see Table 2) the fatigue life equals  f = 7.03 years, which has not been depicted in this figure.

Variation study of static model parameters 𝑡 u0
1 and  0

I,c
In order to assess how the fatigue life  f depends on the normal (tensile) strength  u0 1 and the mode I toughness  0 I,c of the adhesive blade-strut connection, these two static material parameters are successively varied as  u0 1 ∈ [50,70,90,110] MPa and  0 I,c ∈ [1, 2, 3, 4, 5, 6] N/mm.The other material parameters correspond to those of the reference case, as listed in Table 2. Fig. 18 illustrates that an increase of the normal strength from 50 MPa to 110 MPa increases the fatigue life by approximately a factor of 1.7, from 5.43 years to 9.16 years.In contrast, the variation of the mode I toughness only has a minor influence on the fatigue life, which therefore is not depicted here; specifically, the fatigue life predictions for mode I toughness values in between 1 N/mm and 6 N/mm differ by less than 14% and lie in between 6.72 years and 7.66 years.This outcome suggests that the decrease of the fatigue resistance of the blade-strut interface is mainly characterized by the monotonic reduction of the tensile strength of the adhesive, and is relatively insensitive to the reduction (i.e., the value) of the mode I toughness.Indeed, at damage completion,  = 1, the local traction-separation law shown in Fig. 9 generally fully collapses,  f = 1, whereby the crack face separation typically appears to be less than 10% of the ultimate crack face separation under static failure,  < 0.1 u0 .In other words, less than 10% of the static mode I toughness  0 I,c =  u0 1  u0 1 ∕2 is exploited during the fatigue fracture process, so that the fracture process may be typified as ''brittle'', see also the discussion at the end of Section 3.2.The results above are in agreement with the general notion that fatigue fractures often are brittle in nature and develop under relatively small crack openings, as seen in various materials (Suresh, 1998).

Effect of an initial flaw on the fatigue life
During its lifetime, the VAWT should be damage tolerant, such that it has the ability to endure flaws and defects safely until the moment these are detected and repaired.From this aspect, it is important to understand how the presence of an initial flaw influences the fatigue life of the VAWT.Hence, FEM analyses are performed on a bladestrut connection that contains an initial, ellipse-shaped flaw, for which the influence on the fatigue life is studied by varying (i) the flaw location and (ii) the flaw size.For assessing the influence by the flaw location, in both the upper and lower blade-strut connections shown in Fig. 10 an ellipse-shaped flaw with a width of 0.1 s and a height of 0.05 s is considered (in correspondence with a width-to-height aspect 1 (of the adhesive blade-strut connection) on the predicted fatigue life  f .The downtime factor of the VAWT is taken as  = 0.5.The reference case corresponds to  u0 1 = 70 MPa, see also Table 2.
ratio of 2), whereby  is the chord length of the strut.The flaw location is varied by its center along the chord line at respective distances from the strut trailing edge of 0.25 s (close to the trailing edge), 0.50 s (at the center) and 0.75 s (close to the leading edge).
The mesh density of the FEM models is comparable to that used for studying the (flawless) reference case, as presented in Section 4.3.Further, the material parameters correspond to those of the reference case, as listed in Table 2. Fig. 19 shows the predicted fatigue life  f for the three flawed blade-strut connections, together with the fatigue life of the reference case computed in Section 4.3.Observe that both the presence and the location of the flaw have a significant influence on the fatigue life of the VAWT.Specifically, a flaw located close to the trailing edge leads to the lowest fatigue life (5.81years), followed by a flaw at the center of the blade (6.47 years), and finally a flaw located close to the leading edge (6.50 years), resulting in a maximum relative difference in fatigue life of (6.50∕5.81− 1) × 100% ≈ 12%.The reason that the trailing edge of the blade-strut connection is most sensitive to the presence of an initial flaw is because fatigue cracking develops from this location, see Fig. 14(a), whereby the flaw clearly is able to substantially accelerate the damage process.In comparison to the flawless, reference configuration, the relative decrease in fatigue life of the most unfavorable configuration with a flaw close to the trailing edge is (1 − 5.81∕7.03)× 100% ≈ 17%.
The effect of the flaw size on the fatigue life is investigated by respectively increasing and decreasing the ''medium'' flaw area used in the above variation study by a factor of 2, which will be referred to as a ''large'' flaw area and a ''small'' flaw area, respectively.Here, a width-to-height aspect ratio of 2 of the flaw is taken for all the three cases.The location of the large, medium and small flaws is taken at a distance 0.25 s from the trailing edge, as this has shown to lead to the most critical fatigue life, see Fig. 19. Fig. 20 illustrates that an increase of the flaw area by a factor of 2 decreases the fatigue life by 9%, from 5.81 years to 5.27 years, while a decrease of the flaw area by a factor of 2 increases the fatigue life by 5%, from 5.81 years to 6.12 years.In comparison to the flawless reference case, the relative decrease in the fatigue life for the blade-strut connection with the large flaw is (5.27∕7.03− 1) × 100% = 25%.
In conclusion, the study on the effect of the location and size of a flaw in the blade-strut connection shows that the presence of a flaw may lead to a substantial reduction of the fatigue life of the VAWT, making it worthwhile to dedicate special attention to the prevention of flaws during the manufacturing process.

Conclusions
This paper presents an advanced computational modeling framework for the fatigue analysis of VAWTs.The modeling framework uses the aerodynamic loading computed with high-fidelity CFD simulations as input for a detailed, high-cycle fatigue FEM analysis.The FEM analysis is performed using a load-cycle driven fatigue damage model, which enables to simulate high-cycle fatigue cracking processes in an accurate and robust fashion at manageable computational cost.The main conclusions of the modeling study are summarized below.
• The blade-strut connection turns out to be the most critical structural part for the fatigue life of the VAWT, especially when it is carried out as an adhesive connection (instead of a connection).• The fatigue damage initiates at the lower of the leading edge the upper blade-strut connection (i.e., at the most ''fatiguesensitive area'' the VAWT), almost at the lower part of trailing edge.Under an increasing number of load cycles, the damage developing from these two nucleation sites coalesces, and grows towards the upper part of the leading edge, at which it finally leads to catastrophic failure of the complete blade-strut connection, see Fig. 14. • The value of the exponent  in the fatigue evolution law, Eq. ( 6), has a strong influence on the predicted fatigue life of the VAWT, see Fig. 16.The lifetime prediction of a VAWT with an adhesive (FM 73 M OST toughened epoxy) blade-strut connection shows to exceed the design life of 15 to 20 years typically required for domestic wind turbines if the exponent  in the fatigue evolution law has a sufficiently high value, i.e.,  > 1.15.These aspects emphasize the importance of a careful calibration of  for obtaining an appropriate lifetime prediction of VAWTs.
In addition to calibration/verification procedures carried out at coupon level, i.e., the single lap joint experiments considered in Geng and Suiker (2019), calibration/verification studies should be performed at the structural level by measuring, for example, the fatigue response of a large part of the VAWT, or the complete VAWT.In this way, it can be determined up to which extent the fitting of the fatigue model parameter  is influenced by the scale of the experimental configuration used in the calibration/verification procedure.Further, the fatigue tests should be performed at a load ratio  of about 0.8, see Eq. ( 10), as representative of domestic VAWTs operating under urban wind conditions.• The decrease of the fatigue resistance of the VAWT during its lifetime is mainly characterized by the monotonic reduction of the tensile strength of the adhesive blade-strut connection, rather than by the reduction of its mode I toughness.Correspondingly, fatigue fracture of the blade-strut connection develops in a relatively brittle fashion, as characterized by a narrow fracture process zone at the crack tip, see Fig. 14, and a crack opening that generally is much smaller than the crack opening under static failure.• The presence of an initial, ellipse-shaped flaw at the blade-strut connections of the VAWT may decrease the fatigue life by 25%, depending on the flaw size and flaw location.The reduction on the fatigue life is the largest if the flaw is located close to the trailing edge of the blade-strut connection.This result indicates that it is worthwhile to dedicate special attention to the prevention of flaws during the manufacturing process of VAWTs.
As a final remark, it is emphasized that the computational modeling framework presented in this paper is generic; it can also be applied for analyzing the fatigue performance of other rotating machinery subjected to fluid-structure interaction, such as HAWTs, steam turbine generators and multistage pumps and compressors.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. The blade profile is a hollow NACA0018 symmetrical airfoil, which has a thickness of 0.8 mm, two inner spars of thickness 1.3 mm, a chord length  b = 60 mm, and a maximum height of 0.18 b = 11 mm, in accordance with the photograph (upper figure) and the corresponding schematization (lower figure).
Figs. 5 and 6 show the CFD result for respectively  p and  f , as a function of the chordwise position ∕ b , with ∕ b = 0 and ∕ b = 1 respectively corresponding to the leading edge and trailing edge of the turbine blade.For the final (100th) turbine revolution considered in the URANS simulation the pressure coefficient profiles computed at the upper and lower surfaces of the blade are shown for azimuth angles  = 0 • and  = 90 • in Fig. 5(a) and for  = 180 • and  = 270 • in Fig. 5(b).Further, the skin friction coefficient profiles at the upper and lower surfaces of the blade are depicted for azimuth angles  = 0 • and  = 90 • in Fig. 6(a) and for  = 180 • and  = 270 • in Fig. 6(

Fig. 3 .
Fig.3.Schematic of the CFD computational domain and boundary conditions (not to scale).The counter-clockwise rotation of the blades is indicated by the azimuth angle , and is realized by applying a sliding grid interface between the rotating core and the fixed, surrounding domain.

Fig. 7 .
Fig. 7. Normal ( N ) and tangential ( T ) force coefficients on the blade as a function of the azimuth angle  (evaluated during the final (100th) turbine revolution).

Fig. 10 .
Fig. 10.Blade-strut configuration, together with the local FEM discretizations of the blade (upper inset) and the blade-strut connection (lower inset).

Fig. 11 .
Fig. 11.Total strain energy  of the blade-strut system (indicated by the orange bar) and the maximum von Mises stress  eq (indicated by the blue bar) for a coarse mesh (M1) of 29,669 elements, an intermediate, reference mesh (M2) of 64,486 elements, and a fine mesh (M3) of 142,378 elements.

Fig. 13 .
Fig. 13.Maximum normal (tensile) traction   1,max , with  ∈ [0, 1, … 359], in the upper blade-strut connection during one turbine revolution,  ∈[0 • , 1 • , … 359 • ].Independent of the value of , the maximum normal (tensile) traction occurs at the bottom of the blade-strut connection, relatively close to the leading edge.Together with the results depicted in Fig.12, this leads to the conclusion that this location is the ''most fatigue-sensitive area'' of the VAWT.

Fig. 14 .
Fig. 14.Contour plot of interfacial damage  in the adhesive at the upper blade-strut connection shown in Fig.10, for three different numbers of turbine revolutions .The number of years corresponding to the value of  has been calculated by adopting a downtime factor of the VAWT of  = 0.5, using Eq.(7).The undamaged and fully damaged areas are indicated in blue and red, respectively.The parameter values used in the fatigue analysis are listed in Table 2 (reference case).

Fig. 15 .
Fig. 15.The damage parameter  and fatigue damage parameter  f in 8 cohesive elements E1-E8 along the strut chord line of the upper blade-strut connection versus the number of turbine revolutions .The integration points of the 8 cohesive elements are indicated in the upper figure of the cross-section of the blade-strut connection.

Fig. 18 .
Fig. 18.Effect of normal (tensile) strength  u01 (of the adhesive blade-strut connection) on the predicted fatigue life  f .The downtime factor of the VAWT is taken as  = 0.5.The reference case corresponds to  u0 1 = 70 MPa, see also Table2.

Fig. 19 .
Fig. 19.Effect of flaw location (at the blade-strut connections) on the predicted fatigue life  f .The blade-strut connection without a flaw (= reference case) is indicated as ''Intact''.The downtime factor of the VAWT is taken as  = 0.5.

Fig. 20 .
Fig. 20.Effect of flaw size (at the blade-strut connections) on the predicted fatigue life  f .The blade-strut connection without a flaw (= reference case) is indicated as ''Intact''.The downtime factor of the VAWT is taken as  = 0.5.

Table 2
Geng and Suiker (2019)r blade-strut fatigue model (reference case).More details on the experimental calibration of the material properties of the adhesive connection can be found inGeng and Suiker (2019).