Fundamentals of Renewable Energy and Applications A Review of in-situ Loading Conditions for Mathematical Modeling of Asymmetric Wind Turbine Blades

This paper reviews generalized solutions to the classical beam moment equation for solving the deflexion and strain fields of composite wind turbine blades. A generalized moment functional is presented to effectively model the moment at any point on a blade/beam utilizing in-situ load cases. Models assume that the components are constructed from in- plane quasi-isotropic composite materials of an overall elastic modulus of 42 GPa. Exact solutions for the displacement and strains for an adjusted aerofoil to that presented in the literature and compared with another defined by the Joukowski transform. Models without stiffening ribs resulted in deflexions of the blades which exceeded the generally acceptable design code criteria. Each of the models developed were rigorously validated via numerical (Runge-Kutta) solutions of an identical differential equation used to derive the analytical models presented. The results obtained from the robust design codes, written in the open source Computer Aided Software (CAS) Maxima, are shown to be congruent with simulations using the ANSYS commercial finite element (FE) codes as well as experimental data. One major implication of the theoretical treatment is that these solutions can now be used in design codes to maximize the strength of analogues components, used in aerospace and most notably renewable energy sectors, while significantly reducing their weight and hence cost. The most realistic in-situ loading conditions for a dynamic blade and stationary blade are presented which are shown to be unique to the blade optimal tip speed ratio, blade dimensions and wind speed.


Introduction
In the wind turbine industry it is paramount that components are designed to obtain the greatest obtainable efficiency. As time has progressed selecting the correct materials and aerofoil shape has become predominant in the determination of the optimum blade tip-speed ratio, thus increasing a wind turbines power coefficient, (C p ) and overall turbine efficiency (ƞ total ). It is in the industry's best interest to produce high performance turbine components at the lowest costs, without compromising the turbines structural stability. Solving problems of this ilk drive wind turbine design, research and development forward and one such problem is the subject of this paper.
There is currently a surprising lack of computational models in literature, which accurately predict the deflexion and strain fields of wind turbine blades. However, it has proven from relatively simple physical tests on the blades that relatively straightforward and well developed underpinning theory can be used to predict the deflexion of wind turbine blades remarkably well [1]. This said, testing protocols are quite restrictive, especially for the larger expensive components. For such components therefore, effective computational predictive modeling is required to ascertain specific structural static and dynamic design limits; such models can be employed to virtually test any component in a safe and cost free environment.
The work described throughout this paper considers small scale wind turbines (5m blades) which are representative of those used on wind-turbines producing outputs of 1.5 -50 kW; additionally these size of blades are ideal for validation experimentally [1,2]. Moreover validation of larger components becomes technically and financially prohibitive [3]. Whence, modeling the performance and the actual working limits of the blades (i.e. blades angular velocity and tip speed ratio) will be taken from real case studies evident in the literature [4].

Previous work
Recently [1] researchers have used the classical beam elasticity differential equation in order to effectively model the displacement of standard aerodynamic geometries. Though this work was informative in predicting the required displacement and hence strains fields remarkably well given the simplicity modeling method, these researchers do submit that further work is required particularly with respect to the use of composite materials under salient loading conditions. Moreover, conducting such investigations better predictions can be made of blade behavior when attached to the turbine. In practice it is generally accepted that the geometry of a blade is determined by the aerodynamic properties which the designer requires. The ruling characteristics are aerodynamic coefficients: chord, twist and thickness distribution. When designed blade materials are selected with the aim of minimizing tip deflexion, which reduces the probability of a blade/ tower collision [5]. It is also important that these components last the desired lifetime, around 20 years, while under accepted loads [6].
The aforementioned Whitty et al paper [1] presented closed form solutions to the classical beam elasticity differential equation. The developed analytical method was utilized to effectively model the displacement of standard aerodynamic geometries used for wind turbine blades. These models assume that the components are constructed Page 2 of 7 aerodynamic loading conditions, tip-speed optimization for the model blades considered and salient benchmarking procedures.

Generalized moment equation derivation
The most important equation which is needed to model these loading cases is the second order differential equation for the displacement of a beam (equation (3)) [10].
The differential equation comprises of two other function which vary depending on, z, the distance along the blade. M (X, z), is the generalized moment functional, this equation will depend upon the loading case being modeled. These will be discussed and explained in the sections to follow. I (z), is the second moment of area, this function also depends on, z.
To our knowledge no derivation of a generalized moment functional for a wind-turbine blade is evident in the literature. This derivation requires the consideration of Figure 2 which shows a blade fixed at one end with general load intensity applied. Taking moments about the root renders: From the diagram (2), there are four variables, M, the moment, W, the force, q, the load intensity and ζ, the centroid of the section of load intensity to the left of z. Additional there are spacial coordinates z and ζ. These being respectively the distance from the root and an arbitrary distance to the centroid to the left of z. consider an infinitely small section of the load distribution, dz; the force at this point is given by: The total force on the component is therefore found from the integral of the distribution over the whole length: The moment equation can be derived, considering the load distribution acting over the infinitely small distance dz, whence the moment is: from in-plane quasi-isotropic composite materials, forming a blade of shell thickness with 5 mm, with or without a stiffening web. In these cases, a composite comprising of 40 layers of laminated carbon and glass substrate with an epoxy matrix was used, the effective modulus of the material being 42 GPa. The authors [1] derive an explicit solution for the strain and displacement fields for an elliptical shell formulation which is presented for the first time in the literature. This is then expanded to encompass the FX66-S-196 and NACA 63-621 symmetric aerofoil approximations. However, it has been shown that the aerofoils stated above are not symmetric ( Figure 1). Thus this work investigates the effects of the asymmetry on the displacement and strain fields using comparable methodologies.
The experimental static testing of a composite turbine blade has also been carried out by a Jordanian group [2]. While focusing mainly on the manufacturing of the blade, it promotes a link between their local academic work and local industries. The works in depth look into the design and construction of the blade as well as a thorough application of static testing protocols provided valuable data in which to verify modeling work [1] and ultimately design codes [2].

Aerofoil selection
Since it was demonstrated that small changes in the geometry had negligible effect on the structural response, the original Whitty et al [1] aerofoil ( Figure 1) used in their simulations was significantly simplified. The aerofoil is effectively half an ellipse forming the leading edge and a triangle forming the trailing edge. Following a survey of literature it was found that there is normally no symmetry between the top half and the lower half of the aerofoil, thus the Whitty et al. (adjusted) aerofoil was designed. These are then compared with the FX66-S-196 and NACA 63-621 aerofoils which are plotted using discrete points obtained from "Airfoil Tools" [7]. Finally the Continuous NACA aerofoil is aerofoils constructed from the following equations, derived from the Joukowski transform [8]: Both equations depend on the length of chord, x, the diameter of the aerofoil at the furthest point, d top and d bottom , respectively and the total chord length, lc. The constants A, B, C, D and E are all specific to the aerofoils aerodynamic properties [9]. Equations (1) and (2) are effectively the same equation manipulated to produce two connected curves at (0,0) and (0,600). These sections were plotted using the CAS Maxima ( Figure 1).

Scope
This paper provides a review of asymmetry of particular aerofoil sections for the modeling of displacement and strain fields of a wind turbine blade in service. A generalized blade moment functional which can incorporate any loading condition on the blade is also reviewed. These are employed to define the so called generalized aerodynamic load equation and different benchmarking procedures are described which should be routinely used for verification and validation of the computational models used across the industry.

Methods
In this section we review mathematical and simulation methods which can be employed to describe standard industrial testing protocols. These take the form of a generalized moment functional If higher order infinitesimals are neglected: dM = q (z). z dz Hence, the moment at the root is found by integrating over the length of the blade: To make a general solution the moment equation must consider a changing centroid as moments are taken at points further from the root. It is important therefore that the moment equation be considered in the ζ direction. As the moment equation varies with length this creates a unique problem, when any load distribution is modelled the centroid of said load distribution will continually change the further along the blade moments are take. Hence, the well known expression for the centroid of area is employed [11]. It is now possible to substitute all the known parameters into equation (4) which gives the required generalized form: Where z, is the distance along the blade length from the root, the function X(ζ) is the load distribution function at ζ, which is the distance along z, from the root, at which the centroid of the load distribution occurs. With the displacement field now solvable, the global strain field is required. This takes a much simpler algebraic relationship: Equation (7) is employed in the models of this work and it is important to note that these models will only consider cross-sectional variations in the strain with material properties remaining constant.

Force Distribution Derivations
In reality there are a number of forces on a turbine blade that can cause deflexion. The purpose of designing an aerofoil for a wind turbine application is so that at relatively low wind speeds the lift force, produced by the aerodynamic properties of the blade, is large enough to overcome the drag force and the force on the blade due to gravity. This said, in this work the models will not consider the effect of gravity on the system. It is important to simulate the worst possible scenario when testing the blades after manufacturing; this ensures that the blade meets the design standards and is to for purpose. When considering tip deflexion the worst case scenario is to assert a point load on the blade at the end (i.e. furthest distance from the rotor axis). The magnitude of this force would be the maximum possible dictated by the particular standard under consideration [1,12,13], although this is usually unrealistic, though informative as it imposes a further factor of safety into the blade design process. When the blade is modeled, determining a suitable load is very important. The lift force, F L , can be determined by utilizing equation (8): Here C L , is the Lift coefficient, ρ, is the density of air, 1.225 kg/m 3 , L, is the Blade length, L C , is the chord length in meters which, due to the tapering, depends on the blade length thus: Where b r , is the breadth at the root, and b t , is the breadth at the tip of the blade. The parameter z here is the changing distance along the blade. Finally the important function V 2 R , must be explained, this is the relative velocity along the surface of the blade, depending upon the tangential wind velocity, V t. m/s, the angular velocity of the blade, ɷ, the radius of the turbine, r, and the axial velocity, V θ , equation (10): Equation (10), can be determined from the Bernoulli momentum equation for a free flow stream. From equations (8) and (10) a worst case lift force can be determined and used in the models, the worst case occurring at the highest possible tangential velocity that the turbine could experience in its lifetime [5]. As previously indicated ( Figure 2 and equation (6) an important parameter which must be determined for all the models is load intensity, q(z), that is, the force per unit length; which is relent on equation (8).
where, ω(z) is the load intensity for the particular loading case (e.g. PL, UDL, etc.), as indicated in Table 1.
These load distributions were developed after some key findings during modeling. Originally the UDL distribution was considered to be analogous to the stationary distribution on the blade, along with the RTDL which was an approximation to the dynamic load distribution. However, further investigation reveled that the most realistic loading case was that of the so called Bernoulli Dynamic Distributed Load (BDDL) and the special case 1 of this the BDDL termed the Bernoulli Stationary Distributed Load BSDL. It should be noted that the Point Load (PL) is modeled here using the theory of distributions. That is, multiplication of the Dirac distribution 2 thus concentrating the force at a point, rendering what is referred to in Table 1, as the Dirac Distributed Load (DDL). Figure 3, depicts the comparisons between each of the dierent loading cases, here application of equation (5) shows that the total load in each case corresponds to the, required Det Norske Veritas standard DNV-DS-J102, test load of 7600N.
Whence, in that case of the DDL, equation (12), the resulting distribution is that of an impulse showing a zero value until the tip  Table 2. We note here for completeness that unrealistically high speeds of 68 m/s and 110:4 m/s must be applied in order to maintain the required factor of safety as dictated by the standards [12,13].

Tip speed optimization
When the initial blade selection took place, the blade dimensions were taken directly from reference [1] and then modified to give the asymmetry. The blade is then designed for the optimization of the tip speed ratio which is the most important factor. The tip speed ratio is defined by: where ω, is the rotational speed in radians per second and U ∞ , is the free stream velocity. Typically accepted values for the tip speed ratio lie between 4 -8 [5,3]. This is the classical mean of blade design (e.g as used in industry). From reference [4], it is shown that the optimal tip speed ratio for a wind turbine depends on the number of blade n, and can be calculated: In the case utilized in this work, the turbine has 3 blades; as such application of this expression renders an optimum tip speed ratio of 5.45. This tip speed ratio should be maintained at rated wind speed, which for the vast majority of small scale turbines is 10 m/s [3]. Rearrangement of equation (18) and the substitution of this expression, the rotational speed of the turbine at 10 m/s is 10.9 rad/s. As this work considers the worst case scenarios which occurs at wind speeds of greater than 24 m/s, when the optimal tip speed ratio is maintained the rotational speed of the turbine is increased by a factor of 2.6. 1

Blade testing:
In static testing, the blade should be loaded to its most severe loading. In accordance to the Det Norske Veritas standard DNV-DS-J102, deflexion of the blade should always be less than a specified upper limit in order to avoid blade contact with the tower or other components. The design criterion for deflexion of the blade is given by the following inequality [12]: Here, F k is the characteristic load, ∆d is the largest tip deflection when passing the tower, ∆i is the smallest distance from blade tip to the tower or other obstacle in the unloaded condition, f is the load factor, ϒ m is the material factor and ϒ n is the consequence of failure factor. Typically, the load factor takes a value between 0.9 and 1.5 depending upon the desired situation. In the case where the blade could experience the most load (transportation and installation) the suggested value is 1.5 [13]; when under favorable loads this can be as low as 0.9. According to the aforementioned design code the so called material safety factor is a product of a number of other factors in the range (0.95 -1.3). Finally the consequence of failure factor is normally taken to be 1.1 [13]. Utilizing equation (19), explicit knowledge of Classical Lamination Theory (CLT) from composite materials mechanics that composites and the maximum strain failure criterion, a maximum tip deflexion of 10% of the blade length can be determined. For the blade modeled in this work therefore the design criterion was such that the blade could not deflect more that 420 mm under the test load of 7600 N.
During static testing the applied test load is greater than the design load to account for in influences from temperature, humidity, production variations and other environmental aspects during the life of the blade. The test load, F T , is determined as: . .
where, ϒsu is the blade-to-blade variation factor considered to be 1 (worse case testing), ϒ f is the load factor and F k is the characteristic load (i.e. the in-situ load). The characteristic load for this blade is ~5050N, rendering the aforementioned test load of ~7600N [12].

Model verification (benchmarking)
The benchmarking case was selected to show the confidence in the analytical model developed by comparing the results between the benchmarking cases in Whitty et.al. [1], the ANSYS-FE model data, this comparison is shown in Figure 4. Here, the benchmarking case used the original symmetrical aerofoil, under a point load of 8500 N (Figure 4). The work described in this paper slightly underestimates 1 obtained by setting ω to zero.    the deflexion when compared to the original published work. This is due to previous work allows the chord lengths and aerofoil maximum height to vary independently, however the current work varies these parameters simultaneously in an attempt to maintain the aerodynamic properties of the aerofoil.
The Whitty et.al. aerofoil [1] varies the chord length and the height 3 of the blade linearly using two separate functions which are independent of each other and therefore the deflexion is greater as the blade is modeled to be less stiff; a comparison of the maximum deflexion values is shown in Table 3.
The fact that the FE method over estimates the stiffness of any continuum system coupled with the shell element formulation employed in the model explains the lower predictions from the ANSYS software. The new analytical method predicts 13.6% greater deflexion than the FE model, whereas the method employed in the reference [1] suggests a value 22.7% greater than the FE model. These modeling methods suggest a displacement of 567± 89 mm 4 .
Using the new analytical model the maximum strain on the blade occurs not at the root of the blade, as originally expected, but it occurs at ~2250mm along the blade giving strain between about ~0.36%-0.41% depending upon the method employed. When compared to the results obtained in [1] (Figure 5) then it is clear that although the strain is expected to occur at that position, there is also an increase in the maximum strain prediction. Prima-facia this is unexpected as the maximum displacement is now less due to the increased overall stiffness (Table 3). However, this may be attributed to blade not being in a fully strained condition as reported in the reference [1]. That is, up to half the span of the blade the newly developed analytical model is most probably stiffer than that reported previously and ispo-facto is more optimally strained. Accordingly CLT dictates that de-lamination should not occur until much larger strains are realized.
When comparing the strain fields, it is important to note, for all results discussed in the proceeding section, that ANSYS-FE model maximum strain values at the root are too high due to the Saint Venant edge effect and therefore a more realistic maximum strain must be interpolated.

Result
The following results were all obtained using the parameters displayed in Table 2 and Table 4, both the developed analytical and ANSYS-FE models used these parameters. Figure 6 shows the application of equations (13) through (16); the calculated deflexion of the adjusted aerofoil under the various loading cases ( Table 1). As expected, the significant reduction in the aerofoils second moment of areas due to the asymmetry of the new design, greater deflexions are predicted by all the modeling procedures.

Displacement field results
The maximum deflexion is observed on application of DDL, where the total force is distributed such that the moment is at its greatest. The lowest deflexion is realized using the SDL, where the centroid of the load distribution is closer to the rotor axis when compared to all others. Table 5 illustrates the maximum deflexion under each of these load distributions.
As the loading cases change as does the centroid of the loading distribution this affects the deflexion. The model predicts that the further along the centroid of the load distribution from the rotor axis, the greater the deflexion.

Strain field results
Analytical model results for the strain fields are shown in Figure  7. As expected, the largest strains are observed on the application of a point load, ~0.72%. The strain due to the point load, unlike the other loading cases, shows that the maximum strain occurs 2250 mm from the root. This was expected when consulting literature on classical lamination [14] and previously developed tapered beam theory [1].
The maximum strain for the remainder of the loading cases all occur at the root (Table 6), however there is a trend, as the centroid of the load distribution moves further from the root of the blade the strain field becomes increasingly non-linear.
Tables 5 and 6 compare the strain fields which are generally increase with displacement fields, as expected. These results also hold when considering the failure over time, from [5] it is known that wind turbine blades tend to fail due to maximum strain which occur at the root. The only loading case contradicting this phenomenon at present is the point load but as demonstrated previously the point load is not a realistic in-situ loading condition.

Stiffening web
A single 5mm thick stiffening web, in the centre of the aerofoil, along the full length of the blade was added (Figure 8) in line with the literature [1,2] and current industrial practices. The analytical and ANSYS_FE model comparisons including the stiffening rib are shown in Figure 9. Application of a point load renders the maximum deflexion predicted is ~630 mm and 720mm respectively for the analytical and FE model calculations respectively.
The points in Figure 9 show the deflexion predictions from the ANSYS-FE UDL case and the PL case. Following this solution, in line with previous published work [1], the second moment of area of the stiffening web was increased using the analytical model such that the DDL deflexion field was within the design criteria of less than 10% deflexion, therefore reducing proportionally the delfexion the other loading cases considered ( Figure 10).
It was found the analytical model increases the second moment of areas from 12.5 million mm 4 to 15 million mm 4 when the stiffening rib is included (Figure 11).

Conclusions
A mathematical method for solving the deflexion and strain fields of composite wind turbine blades has been reviewed. The model makes         . This development will enable the deflexion and strain fields of any type of aerofoil constructed from quasi-isotropic materials. This being indiscriminate, whether said aerofoil is for a wind turbine, commercial aircraft or indeed other such industrial aerodynamic component. The salient conclusions gained from the work depicted in this paper are as follows: The asymmetric nature of the blade sections designs reduces the second moment of area significantly from between 20.5 -24 million mm4 for the symmetrical aerofoil to between 9.5 -12 million mm 4 , the asymmetry results in a reduction of between around 51% -55%.
Due to the reduction of the second moment of area, the stiffness of the blades is reduced, rendering a higher predicted deflexion. The       maximum deflexion of the asymmetrical blade and their corresponding maximum strain values obtained from the analytical model are seen in Table 7.
In general for this and the previous work [1,2] the solid formulation predicts greater than the shell counterparts, indicating that the shell formulations are stiffer, as expected.
The most realistic loading condition for a dynamic blade is that produced by the BDDL, equation (16). Whereas the most realistic loading condition for a stationary blade is that produced by the BSDL equation seen in equation (16). It is most important to note however, the lift force, equation (8), is unique to the blade. That is the lift force is unique to the optimal rip speed ratio, blade dimensions and wind speed.