A Comprehensive Investigation of Linear and Nonlinear Beam Models on Flexible Wind Turbine Blade Load Calculations

: This study was performed to investigate the effects of structural nonlinearity and large deformations on the aeroelastic loads of flexible wind turbine blades. First, a blade structural analysis model was established using the geometrically exact beam (GEB) theory. Subsequently, the blade element momentum (BEM) theory was corrected using the geometrically exact method leading to the development of a geometrically exact blade element momentum (GE-BEM) model. The results from the GE-BEM model indicated that flapwise deformations always reduce blade fatigue loads, while torsional deformations decrease fatigue loads under low wind speeds but increase them under high wind speeds. Finally, the linear Euler–Bernoulli beam and the GEB were compared to explore the influence of geometric nonlinearity on the blade aeroelastic loads, which revealed that the Euler beam model underestimates the blade loads. The simulations that used the GEB model produced torsional root twist fatigue loads that were 57.49% greater than those generated when the Euler beam model was used. Furthermore, the flapwise bending moment fatigue loads at the root were 8.24% greater than those obtained by the Euler beam model. The smallest discrepancy between the results of the two models was 7.26%, and it corresponded to the edgewise fatigue load.


Introduction
The trend toward the upscaling of offshore wind turbines by increasing their sizes and capacities is intricately linked with the widespread development of global offshore wind farms.Technological advancements, particularly those regarding the utilization of maritime transport, have effectively addressed the logistical challenges associated with transporting large wind turbine components.These advancements have enabled the design and manufacture of turbines on a larger and more flexible scale.Driven by economic considerations, the global offshore wind power industry is progressively adopting highercapacity turbines to enhance the unit power output instead of adding more small units.These colossal offshore wind turbines feature larger, longer, and more flexible towers and blades than previously used turbines [1].For example, the mean rated power capacity of the offshore wind turbines installed in Europe in 2023 was 9.7 MW, which represents a 1.7 MW increase from those installed in 2022 [2].Chinese wind energy enterprises have successfully installed offshore wind turbines with blade lengths that reach 123 m and that have nominal capacities of 16 MW [3].On 10 November 2023, the Dongfang Electric Corporation completed an 18 MW offshore wind turbine with a blade length of 126 m [4].
Projections for the year 2030 predict the emergence of offshore wind turbines with rotor diameters of 250 m and power capacities of 20 MW [5,6].
As the development of larger wind turbines has progressed, turbine blades have become increasingly slender, and the turbine structures have grown increasingly intricate.In addition, to enhance their aerodynamic performance, blades are frequently designed with pre-bend, pre-twist, and pre-sweep characteristics.The blades are designed to withstand substantial deformations, which can passively reduce the loads and enhance the aerodynamic performance.The traditional aeroelastic methods, which assume that the blades are either rigid or undergo small deformations, are no longer reasonable.Modern blades typically have tip deflections that range from 3% to 10% of the blade length [7], though deformations equal to 15% to 20% of the blade length can occur under extreme wind conditions [8].Under normal operating loads, the blade structural deformations of modern large wind turbines have surpassed the limits under which small deformation assumptions are applicable and have exhibited pronounced nonlinearity.Traditional linear beam models that are based on small deformation assumptions are inadequate for accurate blade deformation analyses [9].Moreover, significant blade deformations influence the aerodynamic shapes and load distribution.The vibrations that result from substantial blade oscillations affect the inflow velocity of the wind, thereby significantly impacting the aerodynamic loads; this phenomenon is known as the aeroelastic effect [10].Blade aeroelasticity analyses cannot disregard the aeroelastic loads induced by these effects [11,12].The challenges associated with these aeroelastic effects, which originate from nonlinear blade deformations, have spurred extensive research.Scholars have explored the aeroelasticity problems associated with wind turbine blades by approaching the issue from the following two perspectives: unsteady aerodynamics and nonlinear structural dynamics.
Several methods are commonly used to perform aerodynamic calculations for wind turbines.These methods, which are listed according to increasing computational precision but decreasing computational efficiency, include the blade element momentum (BEM) model, the vortex model, the actuator-type model, and the computational fluid dynamics (CFD) method [13].Despite the pursuit of advanced methodologies to refine wind turbine aerodynamics calculations, engineering models that use the BEM method persist as the dominant aerodynamic models utilized by mainstream wind turbine simulation software, such as DNV Bladed, FAST, CP-Lambda, hGAST, and HAWC2 [13,14].Consequently, researchers have persistently focused on enhancing BEM models to augment computational precision.Zhong et al. (2019) [15] used the CFD method to investigate the flow fields near rotating turbine blades and proposed an innovative methodology for ascertaining the effective angle of attack for individual blade elements within the BEM model.Sun et al. (2016) [16] refined the generalized momentum theory by incorporating considerations for radial flow and improved the wind speed accuracy of the BEM model.Fritz et al. (2022) [17] investigated the influence of aft-swept blades on aerodynamic performance and proposed refinements to the BEM methodology.They found that aft-swept geometries result in a reduced axial induction, particularly around the tip.Sayed et al. (2019) [18] conducted a comparative analysis of aerodynamic loads calculated by both the engineering BEM and CFD methods; their results revealed a tendency for engineering models to underestimate power and thrust.Rezaei et al. (2015) [19] employed a nonlinear geometrically exact beam (GEB) model to represent a turbine blade; the results of this model then served as feedback to an aerodynamic model, which fully incorporated the effects of rotor deformation into the aerodynamic load calculations.Liu et al. (2017) [12] examined unsteady aerodynamic blade loads under dynamic stall conditions by incorporating vibration velocity feedback into the blade aerodynamic model; in their study, they utilized a combination of the BEM model and the Beddoes-Leishman dynamic stall model.In summary, increasing computational accuracy by continually refining the BEM method persists as a pivotal means of addressing engineering and technical challenges present in wind turbine technology.
Common methods of approaching wind turbine structural dynamics calculations include the 3D finite element method and the 1D equivalent beam element method [20].
The existing 1D equivalent beam element models can be broadly divided into linear and nonlinear models.The Euler-Bernoulli and Timoshenko beam models are typical examples of linear beam models.Both the Euler-Bernoulli and Timoshenko beam models rely on the small deformation assumption and are not applicable to wind turbine blades that experience substantial deformations.A frequently employed alternative that is used when there are large deformations is the GEB model [13].The 3D finite element method is typically investigated using commercial finite element analysis software, such as ANSYS or ABAQUS [21].Many scholars are dedicated to enhancing both the aerodynamic and structural models of wind turbines while simultaneously conducting coupled analyses, known as aeroelastic analyses.The goal of these studies is a comprehensive understanding of the aeroelastic characteristics of wind turbines, which can enable the design of more efficient and reliable wind turbines.Early efforts primarily focused on achieving structural-aerodynamic coupling and an understanding of how structural vibration feedback influences aerodynamic models.Wang et al. (2014) [9] and Li et al. (2020) [7] conducted aeroelastic analyses of flexible blades using the GEB method coupled with the BEM method.They demonstrated that the GEB method can more accurately describe large vibration deformations of flexible blades.Utilizing a large eddy simulation solver and a linear beam model, Della Posta et al. (2022) [22] conducted a coupled analysis of wind turbine aerodynamics and structural dynamics, specifically examining the influence of blade torsional vibrations on aeroelastic loads.Mo et al. (2015) [23] proposed an aeroelastic simulation program based on flexible multibody dynamics and combined the BEM method with a dynamic stall model.Sabale and Gopal (2019) [24] utilized the GEB and BEM methods, while accounting for dynamic stall, to establish an aeroelastic analysis model.Using both multibody dynamics and GEB models in hGAST software, Panteli et al. (2022) [25] conducted a comparative aeroelastic analysis that revealed a significant disparity between the aeroelastic load calculation results of linear and nonlinear structural models.Finnegan et al. (2021) [26] compared three FEM numerical models and found that models based on "layered solid elements" give a higher accuracy in FE prediction.Manolas et al. (2015) [7] conducted a study that compared a second-order nonlinear beam model, a multibody model, and a first-order linear beam model for wind turbine blades; they found that the bending-torsional coupling effect was the primary nonlinear effect that caused differences between the loads.
Figure 1 presents calculation results from three wind turbine simulation software packages (GAST, HAWC2, and CP-Lambda) and a computational structural dynamics model (FLOWer-EMPIRE-Carat++) for blade tip deflection in the flapwise direction for a 10 MW reference wind turbine (RWT) [27].The maximum flapwise tip displacement reached 0.125 of the blade length, which is nearly 11.1 m, at a wind speed of 11.4 m/s.As calculated by Panteli et al. (2022) [25], when the tip displacement reached 11.1 m, the tip deflection angle in the flapwise direction was −0.61 rad, while it was 0.031 rad in the torsional direction.The flapwise rotation caused the blade elements at the blade tip to experience horizontal wind speed reductions of nearly 18%.According to the momentum theorem, the thrust of a wind turbine is proportional to the square of the horizontal wind speed.Consequently, this horizontal speed reduction produced a decrease in thrust of approximately 32%.Furthermore, the 0.025 rad blade twist deflection reduced the angle of attack and significantly decreased the lift coefficient.These estimates are preliminary, and when other factors, such as the wind turbine diameter reduction caused by blade bending and the unsteady aerodynamic characteristics of the blades, are included, the blade aeroelastic load calculations become more complex.The effects of blade vibration and deformation on the aeroelastic loads cannot be overlooked.In practical engineering applications, more accurate load calculations performed under a variety of operating conditions would contribute to an improvement in blade design and control precision [10,14,28].Even though meaningful research has been performed by scholars regarding the aerodynamics, structural dynamics, and aeroelastic coupling involved in wind turbine calculations, the mechanisms through which wind turbine deformations affect aerodynamic loads must be further explored.The process of using structural deformation information to modify the BEM model in aeroelastic coupling analyses also requires further investigation.Both qualitative and quantitative analyses are also needed to reveal the effects of structural nonlinearity on the aeroelastic loads of the blades.It was previously stated that the BEM model is characterized by high computational efficiency and precision that meets engineering requirements and that its usage is widespread.Since this study involved extensive simulations of wind turbines operating under a variety of wind conditions for extended periods of time, the BEM model was chosen for use because of these advantages.The GEB model establishes an inherent relationship between displacement and strain, thereby enabling the calculation of blade deformations across any range as well as an accurate representation of the geometric nonlinearity of the blades.Therefore, the BEM and GEB models were coupled for the aeroelastic time-domain simulations of wind turbine blades conducted during this study.
the aerodynamics, structural dynamics, and aeroelastic coupling involved in wind turbine calculations, the mechanisms through which wind turbine deformations affect aerody namic loads must be further explored.The process of using structural deformation infor mation to modify the BEM model in aeroelastic coupling analyses also requires furthe investigation.Both qualitative and quantitative analyses are also needed to reveal the ef fects of structural nonlinearity on the aeroelastic loads of the blades.It was previously stated that the BEM model is characterized by high computational efficiency and precision that meets engineering requirements and that its usage is widespread.Since this study involved extensive simulations of wind turbines operating under a variety of wind condi tions for extended periods of time, the BEM model was chosen for use because of these advantages.The GEB model establishes an inherent relationship between displacemen and strain, thereby enabling the calculation of blade deformations across any range as wel as an accurate representation of the geometric nonlinearity of the blades.Therefore, the BEM and GEB models were coupled for the aeroelastic time-domain simulations of wind turbine blades conducted during this study.In this study, to investigate the effects of significant deformations on the aerodynamic loads of wind turbine blades, a GEB model was first employed to characterize the blades This model achieved precise structural modeling of the position, velocity, orientation, and torsional deformation of each blade element.An accurate description of the wind speed at each blade element was obtained through precise coordinate transformations.Blade deformation and vibration speed corrections were introduced into the BEM model and a local planar assumption was adopted for the deforming rotor plane.A precise geometric correction method for the BEM model, referred to as the GE-BEM model, was also devel oped.This GE-BEM model accounts for variations in the wind speed at the blade elements due to blade vibrations as well as for changes in the attack angle caused by torsional os cillations.To investigate the effects of structural deformations and unsteady aerodynamic characteristics on wind turbine blade loads, different correction conditions were applied in an aeroelastic simulation model.Comparative analyses of results from aeroelastic sim ulation models with different correction conditions were conducted.Finally, extensive aeroelastic wind turbine simulations were performed using both the Euler beam and GEB methods coupled with the GE-BEM model, after which the simulation results were sub jected to fatigue analyses.The impacts of structural nonlinearity on the blade loads were studied by comparing damage-equivalent load results from simulations performed with the Euler beam and GEB models.
This paper is organized into three additional sections.The blade structural model which is based on the GEB method, a corrected BEM model (the GE-BEM model), and an aeroelastic simulation model, are described in Section 2. Verifications of the structural aerodynamic, and aeroelastic models, as well as the effects of deformations and structura nonlinearity on aeroelastic loads, are summarized in Section 3.This study's conclusions are presented in Section 4. In this study, to investigate the effects of significant deformations on the aerodynamic loads of wind turbine blades, a GEB model was first employed to characterize the blades.This model achieved precise structural modeling of the position, velocity, orientation, and torsional deformation of each blade element.An accurate description of the wind speed at each blade element was obtained through precise coordinate transformations.Blade deformation and vibration speed corrections were introduced into the BEM model and a local planar assumption was adopted for the deforming rotor plane.A precise geometric correction method for the BEM model, referred to as the GE-BEM model, was also developed.This GE-BEM model accounts for variations in the wind speed at the blade elements due to blade vibrations as well as for changes in the attack angle caused by torsional oscillations.To investigate the effects of structural deformations and unsteady aerodynamic characteristics on wind turbine blade loads, different correction conditions were applied in an aeroelastic simulation model.Comparative analyses of results from aeroelastic simulation models with different correction conditions were conducted.Finally, extensive aeroelastic wind turbine simulations were performed using both the Euler beam and GEB methods coupled with the GE-BEM model, after which the simulation results were subjected to fatigue analyses.The impacts of structural nonlinearity on the blade loads were studied by comparing damage-equivalent load results from simulations performed with the Euler beam and GEB models.
This paper is organized into three additional sections.The blade structural model, which is based on the GEB method, a corrected BEM model (the GE-BEM model), and an aeroelastic simulation model, are described in Section 2. Verifications of the structural, aerodynamic, and aeroelastic models, as well as the effects of deformations and structural nonlinearity on aeroelastic loads, are summarized in Section 3.This study's conclusions are presented in Section 4.

Methodology
This section presents a detailed description of the aeroelastic simulation model developed during this study.The coordinate systems required for the aeroelastic simulations and the transformation relationships between them are primarily introduced in Sections 2.1 and 2.2.An orthogonal transformation matrix is introduced for the coordinate transformations.The GEB theory and how it is applied to the modeling of wind turbine blades is discussed in Section 2.3.The GE-BEM aerodynamic calculation model and its corrections are the focus of Section 2.4.The methods used to calculate the non-aerodynamic loads on the blades are detailed in Section 2.5.Finally, Section 2.6 explains how the structural and aerodynamic models are coupled in the aeroelastic coupling simulations and provides a discussion of coupling methods.

Orthogonal Transformation Matrix
When characterizing the beam section orientation when applying the GEB model, performing wind speed calculations for the blade elements in the GE-BEM model, and conducting load transformations among different frames in aeroelastic simulations, the coordinate transformations must be handled with precision.An orthogonal transformation tenser, ℜ, and a corresponding matrix, R, were employed for this task.They must satisfy the relationship presented in Equation (1): where R belongs to the Lie group of proper orthogonal linear transformations, SO (3).If E j and t j represent the material frame and the body-fixed motion frame, respectively, their relationship to one another can be expressed by Equation (2): R is the orthogonal matrix of transformation, and it can be regarded as the rotation required to move from the initial configuration to the actual configuration.R also can be understood to represent the coordinates of the rotated basis vectors in the spatial frame.R is often defined in a special frame; however, this representation is not unique and is dependent on the frames.Various parameterization methods can be used to characterize R, such as Euler angles, quaternion, Euler parameters, and rotational pseudovector methods.However, using the rotation vector to represent R is simple and intuitive; this method is also widely employed in kinematic descriptions and coordinate transformations.It was chosen for use in this study, and is expressed by Equation (3): The magnitude, ∥Ψ∥, represents the rotation angle, and its direction is the direc- tion of the rotation axis.ψ x , ψ y and ψ z its components.The choice of using a rotation vector to represent R has some disadvantages, however, such as singularity when ∥Ψ∥= π.The following revision was therefore developed and used during this study when ∥Ψ∥ > π:Ψ * = (1 − 2π ∥Ψ∥ )Ψ.According to the Rodrigues formula [29], R can be expressed by Equation (4): where Ψ is the skew-symmetric matrix formed by the components of Ψ, which are presented in Equation ( 5):

Coordinate Frames and Transformations
To accurately calculate the wind speed, rotational speed, and position of a single blade element, it is necessary to establish different frames at various locations on the wind turbine.To achieve structural-aerodynamic coupling analyses, i.e., aeroelastic analyses of the blade, it is essential to exchange vibration velocities, deformations, and load data within a unified frame.Different frames were established during this study, and transformations between these frames were conducted using the orthogonal transformation matrix discussed in Section 2.1.The definitions of the different coordinate systems, which are illustrated in Figure 2a, are provided in the following list: (1) First, a global inertial Cartesian frame, ē j o, ē 1 , ē 2 , ē 3 , was established, in which ē 1 is oriented vertically upward and ē 3 aligns with the horizontal component of the incoming wind velocity.
(2) The nacelle frame, E nj {o n ; E n1 , E n2 , E n3 }, was constructed by rotating the global frame, ē j o; ē 1 , ē 2 , ē 3 , about ē 1 by the yaw angle, θ yaw , then rotating it in the new frame about ē 2 by the tilt angle, η.The transformation relationship can be expressed by (3) The hub frame, E hj {o h ; E h1 , E h2 , E h3 }, was constructed by rotating the nacelle frame, E hj {o h ; E h1 ,E h2 ,E h3 }, about E h2 by the azimuth angle, θ azi .The transformation rela- tionship can be expressed by E hj = R azi E nj , where R azi = R([0,0,θ azi ] T ). ( 4) The blade root frame, e rj {o;e r1, e r2, e r3 }, was established by rotating the hub frame, E hj {o h ; E h1 , E h2 , E h3 }, about E h2 by the cone angle, β, then rotating it in the new co- ordinate system about E h3 by the sweep angle, σ.The transformation relationship is expressed by e rj =R σ R β E hj , where R β =R([0, β, 0] T ) and R σ = R([0, 0, σ] T ).For modeling convenience, it was assumed that the blade root frame does not change with the pitch motion.

Coordinate Frames and Transformations
To accurately calculate the wind speed, rotational speed, and position of a single blade element, it is necessary to establish different frames at various locations on the wind turbine.To achieve structural-aerodynamic coupling analyses, i.e., aeroelastic analyses o the blade, it is essential to exchange vibration velocities, deformations, and load data within a unified frame.Different frames were established during this study, and transfor mations between these frames were conducted using the orthogonal transformation ma trix discussed in Section 2.1.The definitions of the different coordinate systems, which are illustrated in Figure 2a, are provided in the following list: (1) First, a global inertial Cartesian frame, { } o, , ,

Structural Solver Based on GEB Theory
The GEB theory, which was initially formulated by Reissner (1972) [30] and then refined by 32], enables precise analysis of structural deformations under large displacements and rotations.This theory provides accurate deformation-strain relationships, maintains objectivity with regard to internal forces and moments, and remains unaffected by rigid-body motion and changes in the reference frame [33].In recent years, the GEB theory has been extensively studied and applied to nonlinear structural mechanics problems.During this study, Cardona's approach [34] was adopted, for which rotation vectors were employed to parameterize finite rotations within, and the updated Lagrangian method was used to efficiently manage the updating of the rotation vectors.

Beam Frame and Kinematic Assumptions
In this study, the wind turbine blade, which is illustrated in Figure 2b, was considered to be a typical cantilever beam fixed on a hub.Initially, the blade is represented by a straight beam with length L and cross-section Λ, and it has a reference configuration at B =Λ × [0, L] ⊂ ℜ 3 .To accommodate GEB modeling, a global frame, e j {o; e 1 , e 2 , e 3 }, was established at the root of the blade; this frame coincides with the root frame e rj {o;e r1 , e r2 , e r3 }.The corresponding coordinate values are represented by (x 1 , x 2 , x 3 ).It is important to note that this frame is non-inertial.A fixed Cartesian material frame, E {O′; E 1 , E 2 , E 3 } (which was also considered an attached body frame), was established at x 1 = s (where s is the length parameter).Its origin, o', is located at the centroid of the crosssection.When s = 0, o' coincides with o, E 1 coincides with the centroidal line, while E 2 and E 3 correspond to the first and second principal axes of inertia, respectively; the corresponding coordinate values are represented by (X 1 , X 2 , X 3 ).The vector t j (s, t) o ′′ ; t j (s, t), j = 1, 2, 3 is defined as the basis vector of the motion frame attached to a typical cross-section.Its origin, o ′′ , is fixed at the centroid of the cross-section with the global position vector, x 0 .At each moment, base vectors t 2 and t 3 point in the directions of the principal axes of inertia, while t 1 remains normal to the cross-section.When t = 0 (initial state), the motion frame coincides with the material frame: t j (s, 0) = E j , where j = 1,2,3.The relationship presented in Equation ( 6) represents the basic kinematic assumptions of the beam: In Equation ( 6), x 0 represents the absolute coordinate vector to the centroid of the crosssection, X I t I is the material position vector for any point on the rigid cross-section, and x is the global absolute coordinate vector for any point.In the GEB method, it is assumed that during beam deformation, the cross-section does not deform and remains a plane; in addition, shear effects are accounted for, so the cross-section does not maintain normality with respect to the centroidal line.With these assumptions, the current configuration of the beam can be completely characterized by the position of the centroidal line and the orientation of the moving frame according to Equation ( 7): In Equation (7), R represents the orthogonal transformation matrix that maps the material frame to the moving frame.Based on the kinematic assumptions, the deformation gradient of the beam in the material frame can be expressed by Equation (8): define Γ=R T •dx 0 /ds−E 1 , which represents the material measure of deformation of the neutral axis, and K = dR/ds • R T ; this tenser describes the deformations involving cross-section rotation.Using δR = RδΘ, the co-rotational variation in the deformation gradient can be expressed by Equation ( 9): where δΘ = T(Ψ)δΨ, which represents the infinitesimal rotation variation, δΨ is the rotation vector variation, and T(Ψ) denotes the linear transformation.Detailed definitions were provided by [34].The internal work of the beam is expressed by Equation ( 10): where N and M represent the resultant force and the moment vector of Piola-Lagrange stress, which is defined such that it acts on the surface of the undeformed volume element.
According to the constitutive law, N and M are defined as shown in Equation ( 11): In Equation ( 11), C represents the matrix of elastic coefficients, C .= Diag(EA, GA 2 , GA 3 , GJ, EI 2 , EI 3 ).In this matrix, EA represents the axial stiffness, GA 2 and GA 3 are the shear bending stiffnesses along the transverse axes, GJ is the torsional stiffness, and EI 2 and EI 3 represent the bending stiffnesses.Using the kinematic assumptions in Equation ( 6), the linear and angular inertial forces were derived through the application of the Hamiltonian: In Equation ( 12), .
x 0 is the acceleration vector at the centroid and A = T(Ψ) ..
Ψ, which represents the material angular velocity of the cross-section.The external work of the beam can be expressed by Equation ( 13): where n and M represent the distributed external load vector (in the global frame) and the bending moment vector (in the material frame), respectively.For wind turbine blades, the external load primarily includes the gravitational, centrifugal, and aerodynamic forces.Additionally, n(0) and n(l) are the concentrated load vectors at the nodes of the beam element in the spatial frame, while M(0) and M(l) are the concentrated moment vectors at the nodes in the material frame.By equating the internal virtual work of the beam (Equation ( 10)) with the virtual work of the inertial and external forces (Equations ( 12) and ( 13), respectively), an expression for the principle of virtual displacements was obtained, which is presented as Equation ( 14):

Implementation Details
In order to make the application of the structural model more explicit, three points must be clarified as follows: 1.
A two-node linear Lagrangian interpolation polynomial was employed within each finite element to interpolate the position vectors and rotation vectors of the nodes.To avoid shear locking in the stiffness matrix, reduced Gaussian integration (single-point integration) was utilized when computing the stiffness matrix.2.
To overcome the singularity associated with the rotation vector when ∥Ψ∥= π, which was mentioned in Section 2.1, and to ensure that the rotation vectors of various nodes are not affected by the Ψ correction during interpolation, an updated Lagrangian method was used for the interpolation of the finite rotation vectors.Thus, after completing one step of convergent iteration, the orthogonal transformation matrix, R, of the current configuration was saved.In the next iteration step, the information saved at the node was the rotation vector corresponding to the rotation increment between the previous and current configurations; this method was extended to the two nodes of the element.Thus, the total rotation of the current configuration can be expressed by Equation ( 15): where R ref represents the rotation configuration of the previous convergent configuration, R inc is the rotation matrix for the current rotation increment, and Ψ inc is the rotation increment vector, which can be linearly interpolated between two nodes.

3.
Since an energy dissipation term was introduced into the HHT-α method, which is based on the Newmark method [35,36], to enhance numerical stability and ensure second-order accuracy, it was adopted for numerical integration of the dynamic equations to guarantee stability during lengthy simulations of wind turbine systems.

GE-BEM Aerodynamic Model
Under normal turbulent conditions, the blade tip deforms differently in different directions.The deformations are magnified in Figure 3, which illustrates the deformations in both the flapwise and edgewise directions, to provide a more intuitive understanding.Due to large deformations, the wind turbine rotor is no longer a flat plane but can be considered to possess a symmetric, bowl-shaped surface.Under these conditions, the flat-plane assumption is no longer reasonable.A very short segment at a distance s from the root of the blade, which is denoted as δs, was taken as the object of interest.Its sweeping motion forms the II trajectory, where II is a plane within which the momentum theorem can be adopted.
In order to apply the BEM method, a deformed rotor plane frame, e ′ j o ′′ ;e ′ 1 , e ′ 2 , e ′ 3 , was established at plane II.It was specified that o ′′ be located at the centroid of the deformed blade, while the e ′ 1 direction is along the centroidal line of the blade, e ′ 3 is normal to plane II and aligned with the incoming wind direction, and e ′ 2 is tangent to plane II and perpendicular to e ′ 1 .It is important to note that, due to the beam shear deformation, e ′ 1 no longer coincides with t 1 of the motion frame.
Due to the flapwise deformation of the blade, the deformed rotor plane, plane II, deviates from the ideally formed rotor plane, which is formed by straight blade sweeping.This deviation angle is equivalent to the blade cone angle, β, which was corrected by introducing the flapwise rotation angle obtained from the structural model.Consequently, the pre-bend angle at this point on the deformed rotor plane was adjusted to β + ∆β.Within plane II, because of the edgewise vibrational deformation of the blade, the blade element also changes in the tangential direction.This change is equivalent to the sweep angle of the blade σ.The variable ∆σ was introduced as a correction for this effect; it changes the blade sweep angle to σ + ∆σ.The transformation between the deformed rotor frame and the hub frame can be expressed by Equation ( 16): Due to the flapwise deformation of the blade, the deformed rotor plane, plane II, deviates from the ideally formed rotor plane, which is formed by straight blade sweeping.This deviation angle is equivalent to the blade cone angle, β , which was corrected by introducing the flapwise rotation angle obtained from the structural model.Consequently, the pre-bend angle at this point on the deformed rotor plane was adjusted to β β + ∆ .Within plane II, because of the edgewise vibrational deformation of the blade, the blade element also changes in the tangential direction.This change is equivalent to the sweep angle of the blade σ .The variable σ ∆ was introduced as a correction for this effect; it changes the blade sweep angle to σ σ + ∆ .The transformation between the de- formed rotor frame and the hub frame can be expressed by Equation ( 16): Figure 3 shows that a coordinate transformation must be used to transform the horizontal incoming wind speed v0 into the deformed rotor frame.The transformation relationship is given by Equation ( 17): Only the second and third components of 0 ′ v participate in the calculations in the deformed rotor plane, while the first component, which is in the radial direction of the blade, can be neglected.The vibration speed is obtained in the blade root frame, so it must be transformed into the deformed rotor plane frame according to Equation (18): In Equation ( 18 Figure 3 shows that a coordinate transformation must be used to transform the horizontal incoming wind speed v 0 into the deformed rotor frame.The transformation relationship is given by Equation ( 17): Only the second and third components of v ′ 0 participate in the calculations in the deformed rotor plane, while the first component, which is in the radial direction of the blade, can be neglected.The vibration speed is obtained in the blade root frame, so it must be transformed into the deformed rotor plane frame according to Equation ( 18): In Equation ( 18), .
x 0 represents the blade vibration speed vector, while the second term represents the velocity induced by the torsional blade vibration.Additionally, d s−e denotes the position vector of the aerodynamic center of the blade relative to the torsional center and Ω represents the angular velocity of the torsional blade vibration.The position vector of the blade element, x 0 (s, t), is given by the structural model.The location change in the blade element caused by the blade vibration at s is denoted as ∆r = x 0 (s, t) − x 0 (s, 0) (when neglecting the blade location correction, ∆r = 0; x 0 (s, t) = x 0 (s, 0)).The position vector in the hub frame is expressed by x h = [r hub , 0, 0] T + R β x 0 , where r hub is the hub radius.Utilizing the rotor rotation vector, ω, the relative velocity of the wind at s caused by the turbine blade rotation can be determined (in the hub coordinate system): Subsequently, the relative velocity in Equation ( 19) can be transformed into the deformed rotor plane frame, as shown in Equation (20): The velocity of the blade element consists of the velocities of the large-scale rigid rotation of the blade about the main axis of the rotor and the vibration of the blade itself.In this paper, this velocity is referred to as the structurally synthesized motion velocity, v s .Since velocity is relative, this velocity can be regarded as that of the wind moving toward the blade.It can be expressed by Equation (21): Figure 4 shows that in the deformed rotor plane frame, the second and third components of v s were considered, while the first component, which is in the radial direction, was neglected.Thus, the airflow velocity relative to the blade element was decomposed into normal and tangential components in the deformed rotor frame.The magnitudes of these components can be calculated using Equation ( 22): ( ) where where t ∆ is the simulation time step, , which is the normal aerodynamic force coefficient, and a  represents the time derivative of the axial induction fac- tor.After iterative convergence is achieved for  and  ′ , the normal and tangential forces obtained must be transformed back into the root frame using the expression in Equation ( 30):  The derivation of the BEM model used in this study followed that presented by [37].The axial and tangential induced velocities from momentum theory are expressed as W n = av n and W t = a ′ v t , respectively, where a is the axial induction factor and a ′ is the tangential induction factor.The inflow velocity at a blade element in the deformed rotor frame can be calculated by while the inflow angle is expressed by Equation ( 23): The attack angle at the blade element is given in Equation ( 24): where θ 0 represents the section twist angle, θ p is the pitch angle, and ∆θ is the blade twist deformation angle.For the 10 MW RWT, the corresponding lift coefficient, C L , drag coefficient, C D , and moment coefficient, C M can be obtained from the tables presented by [38].Then, the normal aerodynamic force, F N , the tangential force, F T , and the aerodynamic torque, Ma, at the blade element can be calculated using Equation ( 25): where ρ is the air density and c represents the chord length.The blade tip loss correction factor is provided in Equation ( 26): where B is the number of blades and R 0 is the rotor radius.During the iterative calculation procedure that is followed when applying the BEM method, updated axial and tangential induction factors can be computed using Equations ( 27) and ( 28), respectively: where s l = Bc/2π∥x h ∥, which represents the local solidity of the rotor.A local effective rotor radius correction was introduced by the expression J e = dx/ds ≈ cos(β + ∆β).H is a correction factor, and if a old > 0.3539, where a old is the axial induction factor in the previous time step, then H = 4a old (1−a old ) 0.6+0.61aold +0.79a 2 old .Pitt and Peter's dynamic wake model [39] was adopted when accounting for unsteady aerodynamic effects, and the updated axial induction factor can be computed using Equation (29): where ∆t is the simulation time step, C N = C L cos ϕ + C D sin ϕ, which is the normal aerodynamic force coefficient, and .
a represents the time derivative of the axial induction factor.After iterative convergence is achieved for a and a ′ , the normal and tangential forces obtained must be transformed back into the root frame using the expression in Equation ( 30): (30)

Gravitational and Centrifugal Force Load Calculations
The gravitational force is defined in the global frame and is transformed to the root frame through an orthogonal transformation matrix, where g represents the gravitational acceleration.The centrifugal force can be more conveniently calculated in the hub frame, E hj , then transformed to the root frame through an orthogonal transformation matrix, F c = R T β x h ω 2 .This expression reflects the centrifugal stiffening effect, which causes larger deformations to produce greater centrifugal force components in the corresponding directions.
Because the shear center, elastic center, center of mass, and aerodynamic blade center do not coincide, the distribution of torque on the blade is quite complex.All the torques on the beam that were considered during this study are equivalent to the shear center of the blade; thus, the torques caused by the gravitational and centrifugal forces, which both act on the center of mass, can be expressed with respect to the blade shear center: In Equation (31), d g represents the position vector from the shear center to the center of mass in the root frame.

Coupling Method
In time domain wind turbine simulations, prevalent coupling methods include both strong and weak coupling, which are also referred to as tight and loose coupling, respectively.Tight coupling demands strict adherence to conservation conditions between the fluid and structural domains, which necessitates continuous iteration and results in significant increases in computational costs.Research findings have indicated that the effects of using loose coupling rather than tight coupling for aeroelastic simulations are minimal and that loose coupling can meet computational cost requirements [22].Notably, Sayed et al. (2019) [18], Sabale and Gopal (2019) [24], and Panteli et al. (2022) [25] all employed the loose coupling method.During this study, a loose coupling approach was adopted for the aeroelastic simulations, in which the structural response calculated during the preceding time step (which consists of the blade vibration velocity, the blade deformation, and the position) serves as the input for the load calculations performed in the subsequent time step.The external loads for the current step are computed using parameters from the previous step, which are assumed to remain constant during the structural dynamic response calculation iterations within the current time step.Figure 5 provides a flowchart of the aeroelastic coupling algorithm, in which the GEB model is utilized to calculate the blade structural response.When the Euler beam model is selected for calculation of the structural response, only the green portion of the flowchart must be replaced with a generic Euler beam solver.

Numerical Results
The GEB model that was established during this study was first verified using the static and dynamic portions of a benchmark case with a 45° pre-curved cantilever beam.This was performed to substantiate the analytical capability of the established GEB model

Numerical Results
The GEB model that was established during this study was first verified using the static and dynamic portions of a benchmark case with a 45 • pre-curved cantilever beam.This was performed to substantiate the analytical capability of the established GEB model to capture structural nonlinearity.A comparative assessment was also performed with the Euler beam model to explore the distinctions between linear and nonlinear models.Alongside the structural model, the aerodynamic model and the coupling method both affect the aeroelastic simulations significantly.Subsequently, a comparison was conducted with two globally acknowledged wind turbine simulation software packages, hGAST and DNV Bladed, to assess the dependability of the proposed aeroelastic model.Because these two software packages do not have open-source code and their core algorithms are proprietary, third parties cannot make corrections to them.Therefore, they were chosen to generate baseline results that were used to validate the model developed during this study.The geometry and material data, as well as the aerodynamic and control properties, were obtained from the DTU 10 MW RWT report [38].The aeroelastic simulation model developed during this study was utilized, and an aerodynamic model without any corrections served as a baseline.A comparative analysis was performed to assess the effects of blade deformation on the aerodynamic loads and root loads of the wind turbine blade.Ultimately, a comparison was drawn between the simulation results of the aeroelastic model based on the GEB method and the Euler beam model results.Lastly, spectral and fatigue analyses were conducted to explore the impacts of structural nonlinearity on the blade loads.

Benchmark Static Analysis with the 45 • Pre-Curved Beam
A 45 • pre-curved beam, which is characterized by substantial deformations and torsional nonlinearity, is frequently utilized as a benchmark to validate nonlinear structural models.Various wind turbine simulation tools have also employed this type of beam for computational validations [25,40].To validate the GEB model developed during this study, the model was used to calculate the deformation of a 45 • pre-curved beam when a steady load was applied on the free end of the beam.The geometric and material characteristics of the beam are provided in Figure 6, and the beam was discretized into eight beam elements.Steady loads of 300 N, 450 N, and 600 N were individually applied at the free end of the beam in the positive z-direction.A comparison of the position of the free end of the beam with other published results [31,34] is presented in Table 1.The results are highly consistent and therefore confirm the accuracy and reliability of the structural model developed during this study, which utilized the same internal force derivation method for the beam elements as did Cardona and Geradin (1988) [33].Additionally, an analysis of the results revealed that, due to geometric nonlinearity, rotations around the x-and z-axes occur at the end of the beam along with rotation around the y-axis.The rotation vectors for the crosssection are presented in the fourth column of Table 1.In contrast, the linear Euler beam model could only calculate the rotation around the y-axis and was unable to capture the geometric nonlinearity.analysis of the results revealed that, due to geometric nonlinearity, rotations ar x-and z-axes occur at the end of the beam along with rotation around the yrotation vectors for the cross-section are presented in the fourth column of Table trast, the linear Euler beam model could only calculate the rotation around the y was unable to capture the geometric nonlinearity.

Static Deformation Analysis of the 10 MW RWT
To validate the analytical capability of the model developed during this study for beams with varying cross-sections and to explore the impact of geometric nonlinearity on actual blade deformations, deformation calculations were performed for the 10 MW RWT under static tip loading.Doing so involved utilizing both the GEB model and a linear Euler beam model.The blade used in this study was pre-bent by 3.3 m in the flapwise direction at the blade tip according to reference [38] and had no sweep.The blade was discretely divided into 50 beam elements, and material and geometric coupling was considered in both the flapwise and twist directions.Steady loads ranging from 25 to 250 kN were applied in the flapwise (e 3 ) direction at the blade tip with a load increment of 25 kN.
The results, which are presented in Figure 7, demonstrate a strong agreement between the GEB model established during this study and hGAST software [25], which confirms the reliability and accuracy of the GEB model when it is used to analyze nonlinear beams with varying cross-sections.Figure 7a shows that the flapwise deformation is minimal when a small static load is applied at the blade tip and that Euler beam model results align with the GEB model results.As the load increases, significant deviations arise between the Euler beam and GEB model results.This deviation reaches 16% for a 250 kN blade tip load, thereby revealing the inability of the Euler beam model, which is based on the assumption of small linear deformations, to accurately represent the substantial blade deformations that exist in this scenario.Due to the aerodynamic twist and irregular crosssections within the blades, bend-bend coupling motion arises in the flapwise and edgewise directions.The non-coincidence of the aerodynamic center, centroid, and shear center of the blade cross-section introduces a bending-torsion coupling effect, which results in blade deformations that are characterized by a coupled flexural and torsional response due to external loads.Figure 7b demonstrates that because of the existence of bend-bend coupling, displacement also occurs in the edgewise direction.Figure 7c shows that for a 250 kN static load, substantial torsional rotation occurs at the blade tip because of bend-twist coupling; this rotation reaches a maximum of 2.1 • .Simultaneous displacement in the flapwise direction and negative displacement in the longitudinal direction cause the blade to shorten, as illustrated in Figure 7d.The linear Euler beam model, which is constrained by theoretical limits and the small deformation assumption, failed to capture the bend-bend and bend-twist coupling effects, and its small deformation assumption precluded the simulation of longitudinal displacement; thus, no Euler beam results could be included in Figure 7b-d.
due to external loads.Figure 7b demonstrates that because of the existence of bend-bend coupling, displacement also occurs in the edgewise direction.Figure 7c shows that for a 250 kN static load, substantial torsional rotation occurs at the blade tip because of bendtwist coupling; this rotation reaches a maximum of 2.1°.Simultaneous displacement in the flapwise direction and negative displacement in the longitudinal direction cause the blade to shorten, as illustrated in Figure 7d.The linear Euler beam model, which is constrained by theoretical limits and the small deformation assumption, failed to capture the bendbend and bend-twist coupling effects, and its small deformation assumption precluded the simulation of longitudinal displacement; thus, no Euler beam results could be included in Figure 7b-d

Benchmark Dynamic Analysis with the 45° Pre-Curved Beam
Wind turbine blades, especially those with pre-bend and pre-twist characteristics, are typical slender nonlinear structures.Their structural dynamics closely resemble those of the benchmark 45° pre-curved cantilever beam.To validate the computational capability of the new GEB model when applied to nonlinear dynamic analyses, dynamic simulations were performed for the benchmark beam when it was subjected to external impulse loads.The simulation conditions are depicted in Figure 8a, which incorporates detailed geometric and material properties from the literature [41].A comparison with the literature [41], which is presented in Figure 8b, substantiated the highly consistent and reliable dynamic computational capability of the model.

Benchmark Dynamic Analysis with the 45 • Pre-Curved Beam
Wind turbine blades, especially those with pre-bend and pre-twist characteristics, are typical slender nonlinear structures.Their structural dynamics closely resemble those of the benchmark 45 • pre-curved cantilever beam.To validate the computational capability of the new GEB model when applied to nonlinear dynamic analyses, dynamic simulations were performed for the benchmark beam when it was subjected to external impulse loads.The simulation conditions are depicted in Figure 8a, which incorporates detailed geometric and material properties from the literature [41].A comparison with the literature [41], which is presented in Figure 8b, substantiated the highly consistent and reliable dynamic computational capability of the model.
were performed for the benchmark beam when it was subjected to external impulse loads.The simulation conditions are depicted in Figure 8a, which incorporates detailed geometric and material properties from the literature [41].A comparison with the literature [41], which is presented in Figure 8b, substantiated the highly consistent and reliable dynamic computational capability of the model.

Load Calculations for the 10 MW RWT
To validate the aerodynamic and load calculation models established during this study, simulations were conducted for the 10 MW RWT using the proposed model.The simulation results were then compared with those obtained from DNV Bladed software [42] (version 4.3, which uses Timoshenko beam elements for structural modeling and the BEM method for aerodynamic modeling).This is a widely recognized wind turbine simulation software endorsed by many international wind turbine certification authorities.The simulations, which included no corrections to the BEM model, neglected the blade vibration effects on the inflow wind speed, thereby essentially treating the blade as rigid.A three-dimensional turbulent wind field with an average wind speed of 11.4 m/s and a turbulence intensity of 20%, which considered wind shear and a ground roughness coefficient, z 0 , of 0.2, was generated by DNV Bladed.Yaw misalignment and wind turbine faults were not considered.Initially, aerodynamic loads per unit length, p n and p t , were extracted for comparison at 0.8 R 0 from the rotor center, which is the primary powerproducing region of the wind turbine blade.The results were compared over a 100 s interval, as shown in Figure 9a, for clarity.The comparison revealed excellent agreement between the aerodynamic load calculation results obtained during this study and those generated by DNV Bladed.Minor discrepancies were observed, and these were attributed to differences between the correction methods for large axial induction factors, a, applied in this study and DNV Bladed [42].To further validate the capability of the aeroelastic model to compute the overall blade loads, flapwise and edgewise root moments obtained from the model for the same simulation conditions were compared with those from DNV Bladed.The comparison, which is illustrated in Figure 9b, demonstrates strong consistency between the results from this model and those from DNV Bladed.Minimal differences can be observed in the edgewise direction, and any slight variations in the flapwise direction were primarily attributed to differences between the correction methods for large axial induction factors in the BEM model.
from the model for the same simulation conditions were compared with those from DNV Bladed.The comparison, which is illustrated in Figure 9b, demonstrates strong consistency between the results from this model and those from DNV Bladed.Minimal differences can be observed in the edgewise direction, and any slight variations in the flapwise direction were primarily attributed to differences between the correction methods for large axial induction factors in the BEM model.

Aeroelastic Simulations of the 10 MW RWT
To validate the aeroelastic coupling model proposed in this paper and to investigate the impacts of structural nonlinearity on the dynamic response of the blade, the Euler beam and GEB models were independently integrated into the GE-BEM model, thereby forming an aeroelastic model.Simulations were conducted for the 10 MW RWT wind turbine under the 3D turbulent wind conditions outlined in Section 3.4.These simulations

Aeroelastic Simulations of the 10 MW RWT
To validate the aeroelastic coupling model proposed in this paper and to investigate the impacts of structural nonlinearity on the dynamic response of the blade, the Euler beam and GEB models were independently integrated into the GE-BEM model, thereby forming an aeroelastic model.Simulations were conducted for the 10 MW RWT wind turbine under the 3D turbulent wind conditions outlined in Section 3.4.These simulations accounted for aeroelastic coupling while excluding deformation corrections to the BEM model, and the blade was discretely divided into 50 beam elements.The simulation results were compared with results obtained from DNV Bladed.For clarity, the blade tip displacement and velocity were isolated within an 80 s interval; the results are depicted in Figure 10. were compared with results obtained from DNV Bladed.For clarity, the blade tip displacement and velocity were isolated within an 80 s interval; the results are depicted in Figure 10.
Figure 10a shows that both aeroelastic coupling models established during this study, which were based on the GEB and Euler beam models, produced flapwise and edgewise blade tip displacements that closely aligned with those obtained from DNV Bladed.The tip displacements calculated using the Euler beam model were more similar to those predicted by DNV Bladed than to those obtained from the GEB model.This result primarily occurred because the structural model employed by DNV Bladed is the linear Timoshenko beam model.The linear Euler beam model produced blade tip displacements that were different from those obtained by the GEB model in the flapwise direction, which was characterized by significant deformation.This overestimation of flapwise deformation is consistent with the static analysis results presented in Section 3.2.In the edgewise direction, which exhibits only minor deformation, the blade tip displacement disparities are negligible.Furthermore, due to the small deformation assumption made by the Euler beam model, the longitudinal displacement calculated by the Euler beam model significantly deviated from those produced by the GEB model and DNV Bladed.The GEB model demonstrated more pronounced longitudinal displacement variations than DNV Bladed.According to momentum theory, the flapwise vibration velocity most significantly impacts the wind speed and, consequently, the wind turbine load; thus, the flapwise vibration velocities calculated by the GEB and Euler beam models were specifically compared, as shown in Figure 10b.It is apparent that the flapwise vibration velocity of the Euler beam model surpassed that of the GEB model.This result suggests that the linear beam model also overestimates the flapwise vibration velocity.This overestimation intensifies the feedback effect of the Euler beam model in the aeroelastic coupling model, thereby resulting in a decrease in the inflow velocity for the aerodynamic model and subsequent reductions in the aerodynamic loads.This phenomenon was further verified in a subsequent analysis.

Influence of Deformation on the Blade Loads
To investigate the effects of deformation caused by vibration on the blade aerodynamic loads and the root loads, seven different aeroelastic models were established by Figure 10a shows that both aeroelastic coupling models established during this study, which were based on the GEB and Euler beam models, produced flapwise and edgewise blade tip displacements that closely aligned with those obtained from DNV Bladed.The tip displacements calculated using the Euler beam model were more similar to those predicted by DNV Bladed than to those obtained from the GEB model.This result primarily occurred because the structural model employed by DNV Bladed is the linear Timoshenko beam model.The linear Euler beam model produced blade tip displacements that were different from those obtained by the GEB model in the flapwise direction, which was characterized by significant deformation.This overestimation of flapwise deformation is consistent with the static analysis results presented in Section 3.2.In the edgewise direction, which exhibits only minor deformation, the blade tip displacement disparities are negligible.Furthermore, due to the small deformation assumption made by the Euler beam model, the longitudinal displacement calculated by the Euler beam model significantly deviated from those produced by the GEB model and DNV Bladed.The GEB model demonstrated more pronounced longitudinal displacement variations than DNV Bladed.According to momentum theory, the flapwise vibration velocity most significantly impacts the wind speed and, consequently, the wind turbine load; thus, the flapwise vibration velocities calculated by the GEB and Euler beam models were specifically compared, as shown in Figure 10b.It is apparent that the flapwise vibration velocity of the Euler beam model surpassed that of the GEB model.This result suggests that the linear beam model also overestimates the flapwise vibration velocity.This overestimation intensifies the feedback effect of the Euler beam model in the aeroelastic coupling model, resulting in a decrease in the inflow velocity for the aerodynamic model and subsequent reductions in the aerodynamic loads.This phenomenon was further verified in a subsequent analysis.

Influence of Deformation on the Blade Loads
To investigate the effects of deformation caused by vibration on the blade aerodynamic loads and the root loads, seven different aeroelastic models were established by introducing various corrections to the GE-BEM model that was previously validated.The corrections made in each model are described in the following list: 1.
A model was established with no deformation corrections or dynamic wake effects; it is referred to as GEB.

2.
A model was established with only a torsional deformation correction, ∆θ; it is referred to as GEB + ∆θ.

4.
Only an edgewise deformation correction, ∆σ, was included in the fourth model, which is referred to as GEB + ∆σ.

5.
A model was established with a blade element location change correction, ∆r; this model is referred to as GEB + ∆r.

6.
A sixth model, referred to as GEB + ∆full, was established; it included all four corrections mentioned previously.

7.
The last model included all four corrections and a dynamic wake model; this model is referred to as GEB + ∆full + DW.
Utilizing the correction scheme described in this list, seven aeroelastic models using the Euler beam model were also established and named; they are denoted as Eu, Eu + ∆θ, Eu + ∆β, Eu + ∆σ, Eu + ∆full, and Eu + ∆full + DW, where Eu is an abbreviation for Euler.
Utilizing the seven aeroelastic models based on the GEB model, simulations were conducted for the 10 MW RWT for 600 s under turbulent wind conditions.The simulations used average wind speeds of 7 m/s, 11.4 m/s, and 24 m/s, which represent typical low, rated, and high wind speed operating conditions, respectively.Using all three wind speeds allowed for a comprehensive investigation of the effects of vibration-induced blade deformations on the aerodynamic loads at different wind speeds.
To examine the effects of blade deformation on the aerodynamic loads, the effects of the deformations on the input parameters for the aerodynamic calculations were first investigated.As inferred from the previous analysis, the flapwise deformations were substantial.First, the effects of the flapwise rotation deformation correction, ∆β, on the axial inflow velocity, v n , at the blade element were explored.For clarity and a distinct comparison, the results were extracted at 0.99 R 0 for 50 s under turbulent wind conditions with an average wind speed of 11.4 m/s, as illustrated in Figure 11a.It is apparent that incorporating the ∆β correction results in a decrease in v n , although the reduction is not significant.According to momentum theory, normal aerodynamic loads are proportional to the square of v n , and as a result, a decrease in v n implies significant load reductions.As depicted in Figure 11b, including the ∆θ correction produces a distinct decrease in the attack angle of the blade element at 0.99 R 0 .Based on these two analyses, it can be tentatively concluded that blade deformations contribute to passive reductions in aerodynamic loads.To obtain a better understanding of the variations in blade loads, spectral analyses were conducted for the loads of the 10 MW RWT computed using the seven models discussed previously.They were conducted for 600 s under turbulent wind conditions with an average wind speed of 11.4 m/s. Figure 12 presents the power spectral densities (PSDs) of the aerodynamic loads at 0.8 0 R , while Figure 13 depicts the PSDs of the root bending moments.Clearly, both the normal and tangential aerodynamic force PSDs have dominant peaks at a rotational frequency of 0.16 Hz (1P).These peaks are primarily caused by the turbulence rotational sampling of the blade as it rotates and repeatedly passes through turbulent eddies.It is evident from Figure 12 that the dynamic wake effect most significantly impacts the aerodynamic loads when the mean wind speed is 11.4 m/s.This effect substantially reduces the PSDs of the aerodynamic loads, t p and n p , because the dynamic wake model acts as a dynamic filter for the axial induction factor, thereby preventing it from changing rigidly.There is therefore a gradual change in the induced velocity, which prevents sharp fluctuations in the aerodynamic loads, thereby reducing the overall load.With respect to blade deformation, the flapwise rotation angle, β ∆ , has the most significant impact on the aerodynamic loads in both the flapwise and edgewise directions.
The torsional deformation angle, θ ∆ , has the second-largest effect, while the edgewise deformation, σ ∆ , and the blade element location change, r ∆ , have minimal effects that are nearly negligible.All the deformation factors contribute to reductions in the blade aerodynamic loads.To obtain a better understanding of the variations in blade loads, spectral analyses were conducted for the loads of the 10 MW RWT computed using the seven models discussed previously.They were conducted for 600 s under turbulent wind conditions with an average wind speed of 11.4 m/s. Figure 12 presents the power spectral densities (PSDs) of the aerodynamic loads at 0.8 R 0 , while Figure 13 depicts the PSDs of the root bending moments.Clearly, both the normal and tangential aerodynamic force PSDs have dominant peaks at a rotational frequency of 0.16 Hz (1P).These peaks are primarily caused by the turbulence rotational sampling of the blade as it rotates and repeatedly passes through turbulent eddies.It is evident from Figure 12 that the dynamic wake effect most significantly impacts the aerodynamic loads when the mean wind speed is 11.4 m/s.This effect substantially reduces the PSDs of the aerodynamic loads, p t and p n , because the dynamic wake model acts as a dynamic filter for the axial induction factor, thereby preventing it from changing rigidly.There is therefore a gradual change in the induced velocity, which prevents sharp fluctuations in the aerodynamic loads, thereby reducing the overall load.With respect to blade deformation, the flapwise rotation angle, ∆β, has the most significant impact on the aerodynamic loads in both the flapwise and edgewise directions.The torsional deformation angle, ∆θ, has the second-largest effect, while the edgewise deformation, ∆σ, and the blade element location change, ∆r, have minimal effects that are nearly negligible.All the deformation factors contribute to reductions in the blade aerodynamic loads.
Figure 13 presents comparisons of the PSDs of the root moments calculated by the seven models for a mean wind speed of 11.4 m/s.As for the aerodynamic loads, the dynamic wake effect significantly reduces the flapwise loads at the blade root.However, its impacts on the edgewise root loads are relatively minor, as shown in Figure 13b.With respect to the effects of the different deformations on the root loads, the flapwise rotation deformation exerts the most significant influence, followed by the torsional deformation.The other deformations only minimally affect the flapwise root loads, as shown in Figure 13a.All the deformation corrections exhibit minor impacts on the edgewise root moment because the flapwise loads primarily result from the aerodynamic forces, while the edgewise loads are primarily influenced by the gravitational forces.Figure 13 presents comparisons of the PSDs of the root moments calculated by the seven models for a mean wind speed of 11.4 m/s.As for the aerodynamic loads, the dynamic wake effect significantly reduces the flapwise loads at the blade root.However, its impacts on the edgewise root loads are relatively minor, as shown in Figure 13b.With respect to the effects of the different deformations on the root loads, the flapwise rotation deformation exerts the most significant influence, followed by the torsional deformation.The other deformations only minimally affect the flapwise root loads, as shown in Figure 13a.All the deformation corrections exhibit minor impacts on the edgewise root moment because the flapwise loads primarily result from the aerodynamic forces, while the edgewise loads are primarily influenced by the gravitational forces.To compare the blade deformation effects at different wind speeds, spectral analyses were conducted for the aerodynamic loads at 0.8 0 R under turbulent wind conditions with average wind speeds of 7 m/s and 24 m/s, as shown in Figure 14.Comparisons of the Figure 13 presents comparisons of the PSDs of the root moments calculated by the seven models for a mean wind speed of 11.4 m/s.As for the aerodynamic loads, the dynamic wake effect significantly reduces the flapwise loads at the blade root.However, its impacts on the edgewise root loads are relatively minor, as shown in Figure 13b.With respect to the effects of the different deformations on the root loads, the flapwise rotation deformation exerts the most significant influence, followed by the torsional deformation.The other deformations only minimally affect the flapwise root loads, as shown in Figure 13a.All the deformation corrections exhibit minor impacts on the edgewise root moment because the flapwise loads primarily result from the aerodynamic forces, while the edgewise loads are primarily influenced by the gravitational forces.To compare the blade deformation effects at different wind speeds, spectral analyses were conducted for the aerodynamic loads at 0.8 0 R under turbulent wind conditions with average wind speeds of 7 m/s and 24 m/s, as shown in Figure 14.Comparisons of the To compare the blade deformation effects at different wind speeds, spectral analyses were conducted for the aerodynamic loads at 0.8 R 0 under turbulent wind conditions with average wind speeds of 7 m/s and 24 m/s, as shown in Figure 14.Comparisons of the PSDs of the normal aerodynamic force, p n , at R 0 (as depicted in Figures 11a and 13a,b) reveal that ∆β has the most significant impact on the blade aerodynamic loads at low wind speeds, while ∆θ has the greatest influence at high wind speeds.This result is primarily caused by the torsional vibration of the blade, which gradually increases with the wind speed due to the geometric coupling effect.In contrast, the flapwise deformation of the blade reaches its maximum value at the rated wind speed and then decreases as the wind speed increases.The impact of the blade torsional deformation on the aerodynamic loads also gradually intensifies as the wind speed increases.The torsional deformation weakens the loads at low wind speeds, while at high wind speeds, it enhances them.The primary reasons for this phenomenon are presented next.When the radius exceeds 41.72 m, the FFA-W3-241 airfoil is used for the blade.The aerodynamic data for this airfoil are presented in Figure 15a.Figure 15b shows that, when the wind speed at the blade is 24 m/s, the FFA-W3-241 airfoil operates in the negative attack angle range.The introduction of the torsional deformation correction (which causes the blade to rotate toward the feathering direction) results in a negative twist deformation angle, ∆θ.According to Equation ( 24), a negative twist deformation angle leads to decreases in the attack angle, α.As α continues to grow in the negative direction, as shown in Figure 15a, C l and C d increase; additionally, according to Equation ( 25), the resultant force, F N , also increases.Since the blade elements in this location represent the primary power-production region of the blade, the blade flapwise loads ultimately increase as well.This situation is exactly opposite to that when the blade operates at low wind speeds.For low wind speeds, the attack angle of the main part of the blade is positive; thus, the torsional deformation correction reduces the attack angle, which consequently reduces the resultant aerodynamic force and, therefore, reduces the loads.
wind speed increases.The impact of the blade torsional deformation on the aerodynamic loads also gradually intensifies as the wind speed increases.The torsional deformation weakens the loads at low wind speeds, while at high wind speeds, it enhances them.The primary reasons for this phenomenon are presented next.When the radius exceeds 41.72 m, the FFA-W3-241 airfoil is used for the blade.The aerodynamic data for this airfoil are presented in Figure 15a.Figure 15b shows that, when the wind speed at the blade is 24 m/s, the FFA-W3-241 airfoil operates in the negative attack angle range.The introduction of the torsional deformation correction (which causes the blade to rotate toward the feathering direction) results in a negative twist deformation angle, θ ∆ .According to Equation (24), a negative twist deformation angle leads to decreases in the attack angle, α.As α continues to grow in the negative direction, as shown in Figure 15a, Cl and Cd increase; additionally, according to Equation ( 25), the resultant force, FN, also increases.Since the blade elements in this location represent the primary power-production region of the blade, the blade flapwise loads ultimately increase as well.This situation is exactly opposite to that when the blade operates at low wind speeds.For low wind speeds, the attack angle of the main part of the blade is positive; thus, the torsional deformation correction reduces the attack angle, which consequently reduces the resultant aerodynamic force and, therefore, reduces the loads.To further validate these deductions, damage equivalent load (DEL) analyses were conducted for the aerodynamic loads calculated by the different simulation models under turbulent wind conditions.DELs were used to represent the fatigue damage data contributed by each load.Specifically, the fatigue damage data for each load obtained by rain flow counting was converted into equivalent fatigue damage caused by a single load range repeated at a single frequency.According to the Miner's rule [43], the DEL can be expressed by Equation (32): To further validate these deductions, damage equivalent load (DEL) analyses were conducted for the aerodynamic loads calculated by the different simulation models under turbulent wind conditions.DELs were used to represent the fatigue damage data contributed by each load.Specifically, the fatigue damage data for each load obtained by rain flow counting was converted into equivalent fatigue damage caused by a single load range repeated at a single frequency.According to the Miner's rule [43], the DEL can be expressed by Equation (32): where N ref = 1.0 × 10 7 , which corresponds to a 20-year lifetime damage count, and N i,k represents the number of cycles at a certain load range, S M i,k , and a certain mean load, S M i,k , determined by rainflow counting.T life,i represents the period during which the turbine will operate at the given wind speed of the simulation time series number, i, T sim,i is the simulation time, and m is the slope of the SN curve, for which a value of 10 is typically chosen for glass-reinforced plastic (GRP) [44].The mean value of the aerodynamic load in the normal direction exceeds its vibration amplitude, as shown in Figure 9a.The impact of the mean stress on the fatigue should not be overlooked; therefore, the Goodman correction for the mean stress effects on the fatigue was introduced [45], where S b represents the correction at the ultimate load of the blade (which corresponds to the load applied when the blade reaches its ultimate stress).
Fatigue analyses were conducted for the 600 s simulation effective results of the normal load, p n , at 0.8 R 0 generated by the seven simulation models with average wind speeds of 7 m/s, 11.4 m/s, and 24 m/s.To facilitate the comparison, the results from the GE model were used as a baseline.The relative differences in the DEL values of the baseline results and those calculated by the other six models are presented in Figure 16.It is evident from Figure 16 that the flapwise rotation and torsional angle corrections most significantly influence the aerodynamic fatigue loads, while the edgewise rotation angle and the blade element location change have relatively minor impacts.The dynamic wake correction has a smaller impact on the blade fatigue loads when the deformations are minimal, but it becomes more pronounced as the deformations increase.Blade deformations result in decreases in the fatigue loads at low wind speeds but increases at high wind speeds; these results provide further validation for the conclusions presented previously in this section.

Structural Nonlinearity Effects on the Aeroelastic Loads
It is evident from the analysis presented in Section 3.5 that the linear Euler beam model overestimates the blade deformations and vibration velocity.The further analyses discussed in Section 3.6 revealed that the flapwise and torsional blade deformations, along with the dynamic wake, significantly impact the blade aerodynamic loads.The factors that

Structural Nonlinearity Effects on the Aeroelastic Loads
It is evident from the analysis presented in Section 3.5 that the linear Euler beam model overestimates the blade deformations and vibration velocity.The further analyses discussed in Section 3.6 revealed that the flapwise and torsional blade deformations, along with the dynamic wake, significantly impact the blade aerodynamic loads.The factors that affect the blade aeroelastic loads are diverse, thereby rendering the situation complex.This section presents an analysis of the influence of structural nonlinearity on blade loads by contrasting the simulation results for the Euler beam model with those for the GEB model.To enhance the computational efficiency, six models were selectively employed as follows: the Eu, (Eu + ∆full), (Eu + ∆full + DW), GEB, (GEB + ∆full), and (GEB + ∆full + DW) models.Previous analyses indicated that blade deformation most significantly affects the flapwise loads.Therefore, the DELs of the normal aerodynamic loads, p n , at 0.8 R 0 and the PSDs of the flapwise root moments calculated by the six models for the 10 MW RWT blade under turbulent wind conditions with a mean wind speed of 11.4 m/s for a duration of 600 s were initially compared, as depicted in Figure 17.Subsequently, the DELs of the tangential aerodynamic loads, p t , and the PSDs of the edgewise root moments calculated by the six models were also compared, and these results are presented in Figure 18. Figure 17 shows that for turbulent wind conditions with an average wind speed of 11.4 m/s, the introduction of deformation corrections or the dynamic wake effect decreases both the aerodynamic loads and the root loads, regardless of whether the overall simulation model is based on the GEB model or the Euler beam model.Figure 17a demonstrates that the DELs of the normal aerodynamic force per unit length at 0.8 0 R generated by the simulation models based on the Euler beam model were consistently smaller than those produced by the GEB method.The DEL result produced by the Euler beam model was nearest to the GEB result when no corrections were considered, while the difference was largest when all the corrections were included.Figure 17b shows that, regardless of which corrections were considered, the Euler beam model consistently yielded smaller flapwise root moments than the GEB model.
Figure 18a shows that the comparison between the t p values calculated using the Euler beam model and the GEB model is similar to that for the n p values.The results from the Euler beam model are consistently smaller than those from the GEB model.This trend is also evident in Figure 18b, although the differences between the edgewise root moments calculated by the Euler beam and GEB models are not as pronounced as those between the flapwise root moments.The smaller differences occurred primarily because of the smaller in-plane deformations, which resulted in decreased aeroelastic coupling effects on the aerodynamic loads.Additionally, since the tangential loads are primarily Figure 17 shows that for turbulent wind conditions with an average wind speed of 11.4 m/s, the introduction of deformation corrections or the dynamic wake effect decreases both the aerodynamic loads and the root loads, regardless of whether the overall simulation model is based on the GEB model or the Euler beam model.Figure 17a demonstrates that the DELs of the normal aerodynamic force per unit length at 0.8 R 0 generated by the simulation models based on the Euler beam model were consistently smaller than those produced by the GEB method.The DEL result produced by the Euler beam model was nearest to the GEB result when no corrections were considered, while the difference was largest when all the corrections were included.Figure 17b shows that, regardless of which corrections were considered, the Euler beam model consistently yielded smaller flapwise root moments than the GEB model.
Figure 18a shows that the comparison between the p t values calculated using the Euler beam model and the GEB model is similar to that for the p n values.The results from the Euler beam model are consistently smaller than those from the GEB model.This trend is also evident in Figure 18b, although the differences between the edgewise root moments calculated by the Euler beam and GEB models are not as pronounced as those between the flapwise root moments.The smaller differences occurred primarily because of the smaller in-plane deformations, which resulted in decreased aeroelastic coupling effects on the aerodynamic loads.Additionally, since the tangential loads are primarily caused by gravitational effects, deformations have relatively minor effects.To investigate the blade deformation and structural nonlinearity effects on the fatigue loads during the 20-year operational lifetime of a wind turbine, extensive simulations and fatigue analyses were conducted in accordance with the specification of the IEC 61400-1 standard [46] (IEC 61400-1, 2003), where the design load case (DLC) 1.2 is considered.The simulations were performed under turbulent wind conditions with 11 different mean wind speeds that ranged from 5 to 25 m/s at 2 m/s intervals, and the blade was discretely divided into 50 beam elements.For each mean wind speed, three random turbulence "wind seeds" were used, and each "wind seed" had a duration of 630 s.These conditions resulted in the simulation of 33 load cases for each wind simulation model, and each load case had an effective duration of 600 s.These simulations accounted for wind shear (with a shear exponent, z0, of 0.20), while yaw misalignment and tower shadow were omitted for the sake of simplicity.The wind turbine class selected was 1A, which is the category for higher turbulence characteristics.The wind speed distribution followed a Weibull distribution [38] with a mean annual wind speed of 11.4 m/s and a shape parameter of 2. Throughout the 20-year lifetime, the total number of load cases for each average speed wind speed was scaled according to Weibull hours.
The DELs of the flapwise, edgewise, and torsional blade root moments calculated by the six different models are depicted in Figure 19.For a clearer comparison, the simulation results obtained from the GEB model were used as a reference.The results from the other five models were compared with these reference results, and the relative deviations are presented in Table 2. Furthermore, the relative deviations between the Eu + ∆full + DW and GEB + ∆full + DW models are also presented in Table 2.  To investigate the blade deformation and structural nonlinearity effects on the fatigue loads during the 20-year operational lifetime of a wind turbine, extensive simulations and fatigue analyses were conducted in accordance with the specification of the IEC 61400-1 standard [46] (IEC 61400-1, 2003), where the design load case (DLC) 1.2 is considered.The simulations were performed under turbulent wind conditions with 11 different mean wind speeds that ranged from 5 to 25 m/s at 2 m/s intervals, and the blade was discretely divided into 50 beam elements.For each mean wind speed, three random turbulence "wind seeds" were used, and each "wind seed" had a duration of 630 s.These conditions resulted in the simulation of 33 load cases for each wind simulation model, and each load case had an effective duration of 600 s.These simulations accounted for wind shear (with a shear exponent, z 0 , of 0.20), while yaw misalignment and tower shadow were omitted for the sake of simplicity.The wind turbine class selected was 1A, which is the category for higher turbulence characteristics.The wind speed distribution followed a Weibull distribution [38] with a mean annual wind speed of 11.4 m/s and a shape parameter of 2. Throughout the 20-year lifetime, the total number of load cases for each average speed wind speed was scaled according to Weibull hours.
The DELs of the flapwise, edgewise, and torsional blade root moments calculated by the six different models are depicted in Figure 19.For a clearer comparison, the simulation results obtained from the GEB model were used as a reference.The results from the other five models were compared with these reference results, and the relative deviations are presented in Table 2. Furthermore, the relative deviations between the Eu + ∆full + DW and GEB + ∆full + DW models are also presented in Table 2. Figure 19 shows that, regardless of whether the simulations utilized the Euler beam model or the GEB model, the deformations reduced the loads.A detailed comparison of the results is shown in Table 2. Additionally, an analysis and discussion of the results were conducted from four primary perspectives as follows: (1) The differences between the results of the GEB + ∆full and GEB models can characterize the wind turbine blade deformation effects on the fatigue loads of the blade during its lifetime.The relative differences between the GEB + ∆full and GEB results in Table 2 show that the torsional root load is most affected by blade deformations, which caused a 5.73% increase in the DEL.Blade deformations have relatively smaller impacts on the flapwise root fatigue load, reducing the DEL by only 1.88%, while they have almost no effect on the edgewise root fatigue load.
(2) The differences between the results generated by the GEB + ∆full + DW and GEB models reflect the dynamic wake effects on the blade fatigue loads.The relative differences between the results of these models presented in Table 2 indicate that the dynamic wake increases both the flapwise and torsional root fatigue loads but reduces the edgewise root fatigue load.(3) The differences between the Eu and GEB model results characterize the fatigue load differences solely by structural nonlinearity effects.The relative differences between the EU and GEB model results in Table 2 indicate that, throughout the lifetime of the wind turbine, the Euler beam model calculates fatigue loads that are smaller than those generated by the GEB model in all three directions.The largest difference was obtained for the DEL of the torsional root load; the value calculated by the Eu model was 51.97% less than that generated by the GEB model.The DELs of the flapwise and edgewise root bending moments obtained by the Euler beam model were 7.16% less than those calculated by the GEB model.These differences demonstrate that the linear structural model tends to underestimate fatigue loads.(4) The differences between the Eu + ∆full + DW and GEB + ∆full + DW models can characterize the variations in the fatigue loads that result from the combined wind turbine structural nonlinearity, blade deformation, and dynamic wake effects; accounting for these effects allows the simulations to align more closely with practical considerations.These relative differences are presented in the last column of Table 2, and they indicate that when accounting for the dynamic wake, the Euler beam model calculated a torsional fatigue load that was 54.49% less than that produced by the GE model, while the calculated flapwise and edgewise loads were 8.24% and 7.26%  Figure 19 shows that, regardless of whether the simulations utilized the Euler beam model or the GEB model, the deformations reduced the loads.A detailed comparison of the results is shown in Table 2. Additionally, an analysis and discussion of the results were conducted from four primary perspectives as follows: (1) The differences between the results of the GEB + ∆full and GEB models can characterize the wind turbine blade deformation effects on the fatigue loads of the blade during its lifetime.The relative differences between the GEB + ∆full and GEB results in Table 2 show that the torsional root load is most affected by blade deformations, which caused a 5.73% increase in the DEL.Blade deformations have relatively smaller impacts on the flapwise root fatigue load, reducing the DEL by only 1.88%, while they have almost no effect on the edgewise root fatigue load.
(2) The differences between the results generated by the GEB + ∆full + DW and GEB models reflect the dynamic wake effects on the blade fatigue loads.The relative differences between the results of these models presented in Table 2 indicate that the dynamic wake increases both the flapwise and torsional root fatigue loads but reduces the edgewise root fatigue load.(3) The differences between the Eu and GEB model results characterize the fatigue load differences solely by structural nonlinearity effects.The relative differences between the EU and GEB model results in Table 2 indicate that, throughout the lifetime of the wind turbine, the Euler beam model calculates fatigue loads that are smaller than those generated by the GEB model in all three directions.The largest difference was obtained for the DEL of the torsional root load; the value calculated by the Eu model was 51.97% less than that generated by the GEB model.The DELs of the flapwise and edgewise root bending moments obtained by the Euler beam model were 7.16% less than those calculated by the GEB model.These differences demonstrate that the linear structural model tends to underestimate fatigue loads.(4) The differences between the Eu + ∆full + DW and GEB + ∆full + DW models can characterize the variations in the fatigue loads that result from the combined wind turbine structural nonlinearity, blade deformation, and dynamic wake effects; accounting for these effects allows the simulations to align more closely with practical considerations.These relative differences are presented in the last column of Table 2, and they indicate that when accounting for the dynamic wake, the Euler beam model calculated a torsional fatigue load that was 54.49% less than that produced by the GE model, while the calculated flapwise and edgewise loads were 8.24% and 7.26% smaller, respectively.Furthermore, a comparison of the sixth and third columns reveals that including the dynamic wake effects amplifies the differences in the fatigue loads caused by the structural nonlinearity effects.

Conclusions
In this study, a blade structural analysis model was established using the geometrically exact beam (GEB) theory.It achieved an accurate representation of the significant deformations experienced by wind turbine blades.Subsequently, the blade element momentum (BEM) theory was corrected using the geometrically exact method to describe the wind velocity for a deformed blade, and a GE-BEM model was developed.Finally, an aeroelastic simulation model was established by coupling the structural and aerodynamic models.The validity of the developed model was confirmed by comparing its results for a pre-curved beam and a 10 MW RWT with those of commercial wind turbine simulation software packages.Simultaneously, the impacts of blade deformation on the aerodynamic performance were investigated using the GE-BEM model.Finally, a comparative analysis of aeroelastic simulation results obtained by the linear Euler-Bernoulli beam model (Euler beam model) and the GEB model was conducted to explore the influence of geometric nonlinearity on the blade aeroelastic loads.This study produced four primary conclusions: (1) A comparison with two commercial wind turbine simulation software packages demonstrated that the aeroelastic simulation model developed during this study is reliable.The GE-BEM model established during this study more accurately calculates the wind velocity at the blade elements, thereby bringing the simulation results closer to reality.(2) When a wind turbine operates near the rated wind speed, the flapwise deformation of the blades is at a maximum, and this deformation has the greatest impact on the aerodynamic loads.When a blade operates at high wind speeds, however (near the cut-out wind speed), the torsional vibration of the blade most significantly affects the aerodynamic loads.

31 Figure 3 .
Figure 3. Schematic of the blade deformation and the deformed rotor plane frame.

Figure 4 .
Figure 4. Airflow velocity vectors and aerodynamic forces in the deformed rotor frame.Figure 4. Airflow velocity vectors and aerodynamic forces in the deformed rotor frame.

Figure 4 .
Figure 4. Airflow velocity vectors and aerodynamic forces in the deformed rotor frame.Figure 4. Airflow velocity vectors and aerodynamic forces in the deformed rotor frame.

Figure 5 .
Figure 5. Flowchart of the aeroelastic module algorithm.

Figure 5 .
Figure 5. Flowchart of the aeroelastic module algorithm.

Figure 6 .
Figure 6.Deformation diagram and calculation conditions for the 45° pre-curved beam.

Figure 6 .
Figure 6.Deformation diagram and calculation conditions for the 45 • pre-curved beam.

Figure 8 .
Figure 8.(a) Loading and being of the 45° pre-curved cantilever beam [41] and (b) beam free end time-history responses compared with reference [41].Figure 8. (a) Loading and being of the 45 • pre-curved cantilever beam [41] and (b) beam free end time-history responses compared with reference [41].

Figure 8 .
Figure 8.(a) Loading and being of the 45° pre-curved cantilever beam [41] and (b) beam free end time-history responses compared with reference [41].Figure 8. (a) Loading and being of the 45 • pre-curved cantilever beam [41] and (b) beam free end time-history responses compared with reference [41].

Figure 9 . 8 0R
Figure 9. (a) Aerodynamic loads at 0.8 0 R for the 10 MW RWT and (b) blade root moments of the 10 MW RWT.

Figure 9 .
Figure 9. (a) Aerodynamic loads at 0.8 R 0 for the 10 MW RWT and (b) blade root moments of the 10 MW RWT.

J
. Mar. Sci.Eng.2024, 12, x FOR PEER REVIEW 19 of 31 accounted for aeroelastic coupling while excluding deformation corrections to the BEM model, and the blade was discretely divided into 50 beam elements.The simulation results

Figure 10 .
Figure 10.(a) Blade tip response under turbulent wind conditions and (b) blade tip flapwise velocity.

Figure 10 .
Figure 10.(a) Blade tip response under turbulent wind conditions and (b) blade tip flapwise velocity.

JFigure 11 .R
Figure 11.(a) Axial inflow wind speed at 0.99 0 R in the time domain and (b) attack angle at 0.99

Figure 11 .
Figure 11.(a) Axial inflow wind speed at 0.99 R 0 in the time domain and (b) attack angle at 0.99 R 0 in the time domain.

Figure 12 . 8 0R
Figure 12.PSDs of the aerodynamic forces per unit length at 0.8 0 R under turbulent wind conditions with an average wind speed of 11.4 m/s: (a) normal aerodynamic force PSD, n p , and (b) tangential aerodynamic force PSD, t p .

Figure 13 .
Figure 13.PSDs of the root moments under turbulent wind conditions with an average wind speed of 11.4 m/s: (a) flapwise direction and (b) edgewise direction.

Figure 12 .Figure 12 . 8 0R
Figure 12.PSDs of the aerodynamic forces per unit length at 0.8 R 0 under turbulent wind conditions with an average wind speed of 11.4 m/s: (a) normal aerodynamic force PSD, p n , and (b) tangential aerodynamic force PSD, p t .

Figure 13 .
Figure 13.PSDs of the root moments under turbulent wind conditions with an average wind speed of 11.4 m/s: (a) flapwise direction and (b) edgewise direction.

Figure 13 .
Figure 13.PSDs of the root moments under turbulent wind conditions with an average wind speed of 11.4 m/s: (a) flapwise direction and (b) edgewise direction.

Figure 14 . 8 0R
Figure 14.PSDs of the aerodynamic forces per unit length, n p , at 0.8 0 R under turbulent wind conditions with average wind speeds of (a) 7 m/s and (b) 24 m/s.

Figure 14 .
Figure 14.PSDs of the aerodynamic forces per unit length, p n , at 0.8 R 0 under turbulent wind conditions with average wind speeds of (a) 7 m/s and (b) 24 m/s.

Figure 16 .
Figure 16.Relative differences in the DELs of the aerodynamic forces per unit length at 0.8 0 R calculated by the six correction models and the baseline model under turbulent wind fields with mean wind speeds of (a) 7 m/s, (b) 11.4 m/s, and (c) 24 m/s.

Figure 16 .
Figure 16.Relative differences in the DELs of the aerodynamic forces per unit length at 0.8 R 0 calculated by the six correction models and the baseline model under turbulent wind fields with mean wind speeds of (a) 7 m/s, (b) 11.4 m/s, and (c) 24 m/s.

Figure 17 . 8 0R
Figure 17.Results from six different simulation models for a turbulent wind field with a mean wind speed of 11.4 m/s: (a) DELs of the normal aerodynamic forces per unit length at 0.8 0 R and (b) PSDs of the flapwise root bending moments.

Figure 17 .
Figure 17.Results from six different simulation models for a turbulent wind field with a mean wind speed of 11.4 m/s: (a) DELs of the normal aerodynamic forces per unit length at 0.8 R 0 and (b) PSDs of the flapwise root bending moments.

JFigure 18 .
Figure 18.Results from six simulation models for a turbulent wind field with a mean wind speed of 11.4 m/s: (a) DELs of the tangential aerodynamic forces per unit length at 0.8 and (b) PSDs of the edgewise root bending moments.

Figure 18 .
Figure 18.Results from six simulation models for a turbulent wind field with a mean wind speed of 11.4 m/s: (a) DELs of the tangential aerodynamic forces per unit length at 0.8 and (b) PSDs of the edgewise root bending moments.

Figure 19 .
Figure 19.DELs of the blade root loads calculated by the six different models: (a) flapwise and edgewise root moments and (b) torsional root moments.

Figure 19 .
Figure 19.DELs of the blade root loads calculated by the six different models: (a) flapwise and edgewise root moments and (b) torsional root moments.

( 3 )
Flapwise deformation of a turbine blade reduces the blade loads, while torsional deformation of the blade decreases the blade loads at low wind speeds but increases them at high wind speeds.(4) Due to the overestimation of both the blade deformation and the vibration speed by the Euler beam model, and consequently the overestimation of feedback to the aerodynamic model, the blade loads are underestimated by this model.The difference between DEL values calculated by the Euler beam model and the GEB model developed during this study was largest in the torsional root direction, with a maximum deviation of 53.5%.

erj{o;er1,er2,er3}, was established by rotating the hub frame
Eby the cone angle, β, then rotating it in the new coor dinate system about which represents the local solidity of the rotor.A local effective ro- tor radius correction was introduced by the expression

Table 1 .
Position of the free end of the 45 • pre-curved beam under free end loading.

Table 1 .
Position of the free end of the 45° pre-curved beam under free end loading.

Table 2 .
Comparison of the blade root moment DEL values of the six different models.

Table 2 .
Comparison of the blade root moment DEL values of the six different models.