Structural Optimization of Wind Turbine Blades for Improved Dynamic Performance

Abstract

The design of the main structure of a wind turbine blade is optimized aiming at the improvement of the overall dynamic performance. Three optimization strategies are developed and tested. The first fundamental one is based on minimizing the total structural mass of the blade spar under frequency and strength constraints. The second and third strategies are concerned with the reduction of the overall vibration level by either minimizing a frequency-placement index or maximizing the natural frequencies and placing them at their target values to avoid large amplitudes and resonance occurrence. Design variables include cross-sectional dimensions and material properties along the spanwise direction of the blade spar. The optimization problem is formulated as a nonlinear constrained problem solved by sequential quadratic programming (SQP) technique. Two specific layup configurations, namely, circumferentially asymmetric stiffness ( CAS ) and circumferentially uniform stiffness ( CUS ), are analyzed. Exact analytical methods are applied to calculate the natural modes of vibration of a composite, thin-walled, tapered blade spar. The influence of coupling on the vibration modes is identified, and the functional behavior of the frequencies with the lamination parameters is thoroughly investigated and discussed. Finite element modeling using NX Nastran solver is performed in order to validate the analytical results. As a case study, optimized blade spar designs of a 750-kW horizontal axis wind turbine are given. The attained solutions show that the approach used in this study enhances the dynamic characteristics of the optimized spar structures as compared with a known baseline design of the wind turbine blade.

1. Introduction

Among all renewable energies of different styles, wind energy is the most popularized and potentially applicable type of green energy. Because larger wind turbines have more power capture and economic advantages, the typical size of utility-scale wind turbines, as shown in Figure 1, has grown dramatically over the last three decades [1, 2]. Such large flexible configuration, operating in uncertain environments, gives rise to significant vibration problems and assesses the importance of analyzing structural dynamics in the design of successful wind turbine systems. The main supporting structures of the rotating blades are usually fabricated from thin-walled composite beams with a variety of cross-sectional types. These configurations are used extensively in the design of many aerodynamic structures because of their light weight-to-stiffness ratio and long fatigue life. To design these components, the dynamic characteristics, especially near-resonant conditions, need to be well examined to assure a safe operation.

Figure 1.

Size and power increase of commercial wind turbines.

The objective of this investigation is to optimize the structural dynamics of a thin-walled composite blade through the minimization of structural mass or reduction of the overall vibration level. The latter can be attained directly by maximizing the natural frequencies of the main blade structure under strength and mass constraints. In general conditions, however, we need a material that is as light as possible for a specified stiffness in order to satisfy the design criteria and to minimize the weight-induced fatigue loads.

The main advantages of fiber composite materials [3] are their high strength and stiffness combined with low density, their superior fatigue properties due to the prevention of crack propagation, and their ability to tailor the layup for optimum strength and stiffness. However, sharp transitions between component materials may cause stress and strain discontinuities that facilitate failure [4]. A solution that can be promising to enhance dynamics and aeroelastic stability of composite blades is the use of functionally graded materials (FGMs), in which the mechanical and physical properties vary spatially within the structure. The concept of functionally graded materials was originated in Japan in 1984 during a space project, in the form of proposed thermal barrier material capable of withstanding high-temperature gradients [5]. FGMs may also be developed using fiber-reinforced layers with a changing volume fraction of fibers, rather than constant, producing grading of the material with favorable properties [6].

Considering, next, optimization of wind turbine blades, Maalawi and Negm [7] presented an optimization model for the design of a typical blade structure of horizontal axis wind turbines. The main blade spar was represented by thin-walled tubular beam composed of uniform segments with the design variables chosen to be the cross-sectional area, radius of gyration, and length of each segment. The optimal design is pursued with respect to maximum frequency design criterion subject to mass and aeroelastic constraints. The optimization problem was solved by multidimensional search techniques, where the aeroelastic stability boundaries and steady-state response were calculated using Floquet’s transition matrix theory. Another work by Maalawi [8] developed an optimization model for placing the frequencies of a wind turbine tower/nacelle/rotor structure in free yawing motion. The mathematical formulation considered a single pole tower configuration having thin-walled circular cross section with constant taper along the tower height. The nacelle/rotor combination was modeled as a rigid mass elastically supported at the top of the tower by the torsional spring of the yawing mechanism. The resulting governing differential equation of motion was solved analytically by transforming it into a standard form of Bessel’s equation, which leads to the necessary exact solutions for the frequencies and mode shapes. Useful design charts were developed for placing the frequencies at their needed target values with no penalty of increasing the total structural weight of the system.

In the context of using the concept of material grading, Librescu and Maalawi [9] formulated an analytical approach for obtaining the optimal design of a class of solid nonuniform composite wings with improved aeroelastic stability. The objective function was measured by maximizing the wing divergence speed while maintaining the total structural mass at a value equal to that of a known baseline design. Exact solutions were obtained for different categories of unidirectionally reinforced composite wing structures: namely, the linear volume fraction (L-VF), the parabolic volume fraction (PR-VF), and the piecewise volume fraction (PW-VF) wing models. Results revealed that in general, the torsional stability of the wing can be substantially improved by using nonuniform, functionally graded composites instead of the traditional ones having uniform volume fractions of the constituent materials.

The optimization of large wind turbines was considered by Kun-Nan Chen and Pin-Yung Chen [10], who applied a two-step procedure for finding the optimum design of composite blades. The first step concerned the optimal aerodynamic shape of the blade as described by the chord and twist angle distributions in the spanwise direction. The second step yielded the optimal material distribution. A 3-MW wind turbine with blades having cross sections of NREL S818, S825, and S826 airfoil types is demonstrated as a case study. A parameterized finite element model of the aerodynamically optimized blade was created using the ANSYS software. The optimization results showed that the initial blade model is an infeasible design due to a high level of the maximum stress, exceeding the upper limit of the stress constraint, but eventually the process converges to a feasible solution with the expense of increased total mass of the blade. Another work by Maalawi and Badr [11] considered the excessive wind turbine blade vibrations induced by continuous pitching, which is necessary to limit the power output and protect the generator from damages in severe wind conditions. They utilized analytical Bessel’s functions of the first kind which yield to the exact solutions of the resulting governing differential equation. The associated optimization problem was formulated by considering two forms of the objective function. The first one was represented by a direct maximization of the fundamental frequency, while the second one considered minimization of the square of the difference between the fundamental frequency and its target or desired value. In both strategies, an equality constraint is imposed on the total structural mass in order not to violate other economic and performance requirements. Design variables encompass the blade tapering ratio, chord, and shear wall thickness distributions. Danny Sale et al. [12] developed a numerical methodology for the structural analysis and optimization of composite blades for wind and hydrokinetic turbines. They derived a structural mechanic model which is based upon a combination of classical lamination theory with a Euler-Bernoulli and shear flow theory applied to composite beams. The development of this simplified structural model was motivated by the need for an accurate and computationally efficient method that is suitable for parametric design and optimization studies of composite blades. An important characteristic of this structural model is its ability to handle complex geometric shapes and isotropic or anisotropic composite layups. For a specified design load, the objective of the structural optimization was to minimize the blade’s mass while satisfying constraints on maximum allowable stress, blade tip deflection, buckling, and placement of blade natural frequencies. Adam Chehouri et al. [13] presented an improved version of the preliminary optimization tool called CoBlade, which offers designers and engineers an accelerated design phase by providing the capabilities to rapidly evaluate alternative composite layups and study their effects on static failure and fatigue of wind turbine blades. In this study, the optimization formulations included nonlinear failure constraints, and a comparison between three formulations was made to show the importance of choosing the blade mass as the main objective function and the inclusion of failure constraints in the wind turbine blade design.

A recent work by Maalawi [14] introduced a mathematical model for optimizing dynamic performance of thin-walled functionally graded box beams with closed cross sections. The objective function was measured by maximizing the natural frequencies and places them at their target values to avoid the occurrence of large amplitudes of vibration. Variables include fiber volume fraction, fiber orientation angle, and ply thickness distributions. Various power law expressions describing the distribution of the fiber volume fraction have been implemented, where the power exponent was taken as a main optimization variable. The mass of the structure is kept at a value equal to that of a known reference beam. Side constraints were also imposed on the design variables in order to avoid having unacceptable optimal solutions. A case study, including optimization of a cantilevered, a single-cell spar beam made of carbon/epoxy composite was considered. Conspicuous design charts were developed, showing the optimum design trends for the mathematical models implemented in the study. It was shown that the developed mathematical models are adequately satisfying the required global optimization of typical composite, functionally graded, thin-walled beam structures.

This chapter focuses on the optimization of the main structure of a wind turbine blade by either minimizing structural mass under frequency and strength constraints or maximizing the natural frequencies under mass, strength, and side constraints. This model is applied to tapered, anisotropic spar beam with thin-walled closed cross section made of laminated fibrous composites with variable thickness and stiffness. The study is focused on the spar structure that represents the main load-carrying component of the wind turbine blade. Material grading concept is utilized by changing the fiber content throughout the blade structure. Design variables include the volume fraction distribution of the constituent materials of construction and geometric and cross-sectional parameters of the blade spar. The blade is assumed to have a large span-to-chord ratio, which enabled us to model the main spar as an equivalent straight beam, positioned along the elastic axis. Structural analyses are performed using simplified mathematical expressions by implementing the conventional beam and classical lamination theories. The governing differential equations of motion are derived and solved by the transfer-matrix method for the coupled extensional-torsional and flexural-torsional modes of vibration. A case study is given considering thin-walled blade spar of a 750-KW horizontal axis wind turbine. Numerical results are presented and discussed showing the success of the developed mathematical model in producing efficient blade designs with improved dynamic performance. Finally, the relevant concluding remarks and recommendations for future studies are given and discussed.

2. Governing differential equations of motion

Among the dynamic characteristics of the blade main structures, determination of the natural frequencies and the associated mode shapes is of fundamental importance. An analytical model for the free vibration of anisotropic thin-walled beams with closed cross sections was developed by Armanios and Badir [15] using a variational asymptotic approach and Hamilton’s principles. This model was applied to arbitrary closed cross sections made of laminated fibrous composites with variable thickness and stiffness. The analysis was applied to two kinds of laminated composites: circumferentially uniform stiffness (CUS) and circumferentially asymmetric stiffness (CAS). The model was also implemented in Refs. [16, 17] to investigate the influence of coupling on the free vibration of thin-walled composite beams. Shadmehri et al. [18] studied the static and dynamic characteristics of composite thin-walled beams that are constructed from a single-cell cross section. The structural model considered incorporated a number of nonclassical effects, such as material anisotropy, transverse shear, warping inhibition, nonuniform torsion, and rotary inertia. The governing equations were derived using extended Hamilton’s principle and solved using extended Galerkin’s method. Phuong and Lee [19] presented a flexural-torsional analysis of thin-walled composite box beams. A general analytical model applicable to thin-walled composite box beams subjected to vertical and torsional loads was developed. Analytical solutions for the free vibration analysis of tapered thin-walled laminated composite beams with closed cross sections are given in Ref. [20]. The exact values of frequencies were obtained by means of power series schemes. A parametric analysis was performed for different taper ratios, stacking sequences, and materials.

Considering functionally graded constructions, Kargarnovin and Hashemi [21] investigated the free vibration of a fiber composite cylinder, in which the volume fraction of fibers varies longitudinally, using a semi-analytical method. The distribution of volume fraction of fiber in base matrix was based on power law model. Another study by Liu and Shu [22] developed an analytical solution to study the free vibration of exponential functionally graded beams with a single delamination. They showed that the natural frequencies increase as Young’s modulus ratio of the constituent materials becomes bigger.

Figure 2 shows the structural model of the blade spar, which is represented by a thin-walled cantilever beam which consists of Ns uniform segments. Each segment has different dimensions and material properties that satisfy the geometrical tapering and material grading distribution. Any segment k with length L _k has a rectangular cross section with dimensions, width b _k , depth a _k , and wall thickness H _k . Each segment is a uniform laminated fibrous composite beam which consists of N _r layers, each of which has thickness h _j , fiber volume fraction V _fj , and fiber orientation angle θ _j (j = 1,2,…., Nr).

Figure 2.

Blade spar structural model.

The constitutive relationships in terms of stress resultants and kinematic variables are [15, 16]:

$$\left[\begin{array}{c}\begin{array}{c}T\\ M_x\\ M_y\end{array}\\ M_z\end{array}\right]=\left[\begin{array}{cccc}{C}_{11}& C_{12}& C_{13}& C_{14}\\ C_{12}& C_{22}& C_{23}& C_{24}\\ C_{13}& C_{23}& C_{33}& C_{34}\\ C_{14}& C_{24}& C_{34}& C_{44}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{U}_1^{\prime }\\ {\phi }^{\prime }\\ U_3^{\prime \prime }\end{array}\\ U_2^{\prime \prime }\end{array}\right]$$ E1

where $T$ is the tensile force, $M_x$ is the torsional moment, and $M_y$ and $M_z$ are the bending moments about the y and z axes, respectively. $C_{\mathit{mn}}$ are called the beam cross-sectional stiffness coefficients, and $U_1$ , $U_2$ , and $U_3$ are the average cross-sectional displacements along x, y, and z coordinates, respectively, and $\phi \left(x\right)$ is the elastic twist about the x axis. The prime denotes differentiation with respect to x. Applying Hamilton’s principle, the equations of undamped free vibration are [15]:

$$\begin{array}{l}{C}_{11}{U}_1^{\prime \prime }+C_{12}{\phi }^{\prime \prime }+C_{13}{U}_3^{\prime \prime \prime }+C_{14}{U}_2^{\prime \prime \prime }-m\acute{U_1}=0\\ C_{12}{U}_1^{\prime \prime }+C_{22}{\phi }^{\prime \prime }+C_{23}{U}_3^{\prime \prime \prime }+C_{24}{U}_2^{\prime \prime \prime }-I\acute{\phi }-S_z\acute{U_3}+S_y\acute{U_2}=0\\ \begin{array}{c}{C}_{13}{U}_1^{\prime \prime \prime }+C_{23}{\phi }^{\prime \prime \prime }+C_{33}{U}_3^{\prime \prime \prime \prime }+C_{34}{U}_2^{\prime \prime \prime \prime }+S_z\acute{\phi }+m\acute{U_3}=0\\ C_{14}{U}_1^{\prime \prime \prime }+C_{24}{\phi }^{\prime \prime \prime }+C_{34}{U}_3^{\prime \prime \prime \prime }+C_{44}{U}_2^{\prime \prime \prime \prime }-S_y\acute{\phi }+m\acute{U_2}=0\end{array}\end{array}$$ E2

where m, $I$ , and $S_z$ and $S_y$ are the mass, polar, and first moments of inertia per unit length of the beam, respectively. The dot superscript denotes differentiation with respect to time.

A closed-form solution for the most general case of the equations of motion (Eq. 2) is not available. Two particular cases of fiber layup are considered in which some of the stiffness coefficients vanish. The first case is called circumferentially uniform stiffness (CUS) and the second circumferentially asymmetric stiffness (CAS). Figure 3 shows a rectangular cross-sectional beam segment with both CUS and CAS layup configurations. CUS layup configuration is manufactured by warping the composite layup using filament winding technique, $\theta \left(-z\right)=\theta \left(z\right)$ , while CAS layup configuration is manufactured such that the beam cross section is symmetric about the OXY plane. $\theta \left(-z\right)=-\theta \left(z\right)$ .

Figure 3.

Spar segment with (a) CUS and (b) CAS layup configurations.

2.1 CUS layup configuration

For any segment k within the blade structure, a local coordinate system with $\acute{X}$ -axis is introduced such that (see Figure 2).

$$0\le \acute{x}=X-X_k\le L_k$$ E3

In the special case of CUS layup, the equations of motion reduce to (for the k^th segment):

$$\begin{array}{l}{C}_{11}^k{U}_1^{\prime \prime }+C_{12}^k{\phi }^{\prime \prime }-m_k\acute{U_1}=0\\ C_{12}^k{U}_1^{\prime \prime }+C_{22}^k{\phi }^{\prime \prime }-I_k\acute{\phi }=0\\ \begin{array}{c}{C}_{33}^k{U}_3^{\prime \prime \prime \prime }+m_k\acute{U_3}=0\\ C_{44}^k{U}_2^{\prime \prime \prime \prime }+m_k\acute{U_2}=0\end{array}\end{array}$$ E4

The first two equations express a coupled extension-twist vibration (ETV) mode, while the third and fourth equations express vertical bending vibration (VBV) and horizontal bending vibration (HBV) modes, respectively.

Applying the integral formulas given in Ref. [15], the stiffness coefficients can be determined from:

$$\begin{array}{l}{C}_{11}^k=2K_A^k\left(b_k+a_k-2H_k\right)\\ C_{12}^k=K_B^k\left(b_k-H_k\right)\left(a_k-H_k\right)\\ \begin{array}{l}{C}_{22}^k=K_C^k\frac{{\left(\left(b_k-H_k\right)\left(a_k-H_k\right)\right)}^2}{2\left(b_k+a_k-2H_k\right)}\\ C_{33}^k=\left(K_A^k-\frac{{K_B^k}^2}{K_C^k}\right)\left(\frac{{\left(a_k-H_k\right)}^3}{6}\right)\left(1+\frac{3\left(b_k-H_k\right)}{\left(a_k-H_k\right)}\right)\\ C_{44}^k=\left(K_A^k-\frac{{K_B^k}^2}{K_C^k}\right)\left(\frac{{\left(b_k-H_k\right)}^3}{6}\right)\left(1+\frac{3\left(a_k-H_k\right)}{\left(b_k-H_k\right)}\right)\end{array}\end{array}$$ E5

where $K_A$ , $K_B$ , and $K_C$ are the reduced axial, coupled axial-shear, and shear stiffness coefficients, respectively, given by the expressions:

$$\begin{array}{l}{K}_A=A_{11}-\frac{{\left(A_{12}\right)}^2}{A_{22}}\\ K_B=2\left[A_{16}-\frac{A_{12}{A}_{26}}{A_{22}}\right]\\ K_C=4\left[A_{66}-\frac{{\left(A_{26}\right)}^2}{A_{22}}\right]\end{array}$$ E6

where A_mn are called the membrane in-plane stiffnesses, which depend on the fiber orientation angle θ, volume fractions, and mechanical properties of the fiber and matrix materials [4].

The Young’s moduli in the longitudinal and lateral directions of the lamina E₁₁ and E₂₂ , the shear modulus G₁₂ , and the major Poisson’s ratio v₁₂ are calculated using the semiempirical methods by Halpin and Tsai [23]:

$$\begin{array}{l}{E}_{11}=E_m\left(1-V_f\right)+E_f{V}_f\\ E_{22}=E_m\left(1+{\eta }_E{V}_f\right)/\left(1-{\eta }_E{V}_f\right)\\ \begin{array}{l}{G}_{12}=G_m\left(1+{\eta }_G{V}_f\right)/\left(1-{\eta }_G{V}_f\right)\\ {\nu }_{12}={\nu }_m\left(1-V_f\right)+{\nu }_f{V}_f\\ \begin{array}{c}{\eta }_E=\left(E_f-E_m\right)/\left(E_f+E_m\right)\\ {\eta }_G=\left(G_f-G_m\right)/\left(G_f+G_m\right)\end{array}\end{array}\end{array}$$ E7

Subscripts “m” and “f” refer to the properties of matrix and fiber materials, respectively, and V_ƒ is the volume fraction of fibers within each lamina. Considering the equations of coupled extension-twist vibration, the assumed solution is [24]:

$$\begin{array}{l}{U}_1\left(x,t\right)={\acute{C}}_1{e}^{\mathit{\lambda x}+\mathit{i\omega t}}\\ \phi \left(x,t\right)={\acute{C}}_2{e}^{\mathit{\lambda x}+\mathit{i\omega t}}\end{array}$$ E8

where $\omega$ is the circular natural frequency of free vibration. Substituting from Eq. (8) into Eq. (4), the associated characteristic equation can be shown to be:

$$\begin{array}{l}a{\lambda }^4+b{\omega }^2{\lambda }^2+c{\omega }^4=0\\ a=C_{11}^k{C}_{22}^k-{C_{12}^k}^2\\ \begin{array}{l}b=C_{11}^k{I}_k+C_{22}^k{m}_k\\ c=m_k{I}_k\end{array}\end{array}$$ E9

The solution can be expressed as:

$$\begin{array}{l}{\lambda }_{1,2}=\pm i{\alpha }_1\\ {\lambda }_{3,4}=\pm i{\alpha }_2\\ {\alpha }_{1,2}=\omega \sqrt{\frac{b\mp \sqrt{b^2-4\mathit{ac}}}{2a}}\end{array}$$ E10

If $C_{11}^k{I}_k>C_{22}^k{m}_k$ , then the minus sign in Eq. (10) will generate natural frequencies with axial-mode dominated vibration, while the plus sign will generate natural frequencies with torsion-mode dominated vibration. If $C_{11}^k{I}_k<C_{22}^k{m}_k$ , the domination is reversed. The general solution of the coupled extension-twist vibration problem (Eq. 4) can be expressed as:

$$\begin{array}{l}{U}_1\left(x,t\right)=\left(c_1,\sin {\alpha }_{1}x+c_2\cos {\alpha }_{1}x+c_3\sin {\alpha }_{2}x+c_4\cos {\alpha }_{2}x\right)e^{\mathit{i\omega t}}\\ \phi \left(x,t\right)=\left(c_5,\sin {\alpha }_{1}x+c_6\cos {\alpha }_{1}x+c_7\sin {\alpha }_{2}x+c_8\cos {\alpha }_{2}x\right)e^{\mathit{i\omega t}}\end{array}$$ E11

In order to satisfy both the coupled equations of motion, there are specific relations between the constants such that:

$$\begin{array}{l}{c}_5=q_1{c}_1,c_6=q_1{c}_2,c_7=q_2{c}_3,c_8=q_2{c}_4\\ q_{1,2}=\frac{2m_{k}a}{C_{12}^k\left(b\mp \sqrt{b^2-4\mathit{ac}}\right)}-\frac{C_{11}^k}{C_{12}^k}\end{array}$$ E12

2.2 CAS layup configuration

In this case, the equations of motion, applied to segment k , reduced to:

$$\begin{array}{l}{C}_{11}^k{U}_1^{\prime \prime }-m_k\acute{U_1}=0\\ C_{22}^k{\phi }^{\prime \prime }+C_{23}^k{U}_3^{\prime \prime \prime }-I_k\acute{\phi }=0\\ \begin{array}{l}{C}_{23}^k{\phi }^{\prime \prime \prime }+C_{33}^k{U}_3^{\prime \prime \prime \prime }+m_k\acute{U_3}=0\\ C_{44}^k{U}_2^{\prime \prime \prime \prime }+m_k\acute{U_2}=0\end{array}\end{array}$$ E13

The second and third equations of motion express a coupled bending-twist vibration (BTV) mode, while the first and fourth equations of motion express extension vibration (EV) and horizontal bending vibration (HBV) modes, respectively. The non-zero stiffness coefficients are given by:

$$\begin{array}{l}{C}_{11}^k=2K_A^k\left(b_k+a_k-2H_k\right)-2\frac{{K_B^k}^2}{K_C^k}\left(b_k-H_k\right)\\ C_{22}^k=K_C^k\frac{{\left(\left(b_k-H_k\right)\left(a_k-H_k\right)\right)}^2}{2}\\ \begin{array}{l}{C}_{23}^k=K_B^k\frac{{\left(\left(b_k-H_k\right)\left(a_k-H_k\right)\right)}^2}{2}\\ C_{33}^k=K_A^k\frac{{\left(a_k-H_k\right)}^3}{6}\left[1+\frac{3\left(b_k-H_k\right)}{\left(a_k-H_k\right)}\right]-\frac{{K_B^k}^2}{2K_C^k}\\ C_{44}^k=\frac{{\left(b_k-H_k\right)}^3}{6}\left[\left(1+\frac{3\left(a_k-H_k\right)}{\left(b_k-H_k\right)}\right)K_A^k-\frac{{K_B^k}^2}{K_C^k}\right]\end{array}\end{array}$$ E14

Assuming harmonic solution similar to that given in Eq. (8), the characteristic equation is obtained by differentiating Eq. (13) to get:

$$\begin{array}{l}a{\lambda }^6+b{\omega }^2{\lambda }^4-c{\omega }^2{\lambda }^2-d{\omega }^4=0\\ a=C_{22}^k{C}_{33}^k-{C_{23}^k}^2\\ \begin{array}{c}b=C_{33}^k{I}_k\\ c=C_{22}^k{m}_k\\ d=m_k{I}_k\end{array}\end{array}$$ E15

which has the solution:

$${\lambda }_{1,2}=\pm i{\gamma }_1,{\lambda }_{3,4}=\pm i{\gamma }_2,{\lambda }_{5,6}=\pm {\gamma }_3$$ E16

The general solution takes the form:

$$\begin{array}{c}{U}_3\left(x,t\right)=\left(c_1\mathit{sin}{\gamma }_{1}x+c_2\mathit{cos}{\gamma }_{1}x+c_3\mathit{sin}{\gamma }_{2}x+c_4\mathit{cos}{\gamma }_{2}x+c_5\mathit{sinh}{\gamma }_{3}x+c_6\mathit{cosh}{\gamma }_{3}x\right)e^{\mathit{i\omega t}}\\ \phi \left(x,t\right)=\left(c_7\mathit{sin}{\gamma }_{1}x+c_8\mathit{cos}{\gamma }_{1}x+c_9\mathit{sin}{\gamma }_{2}x+c_{10}\mathit{cos}{\gamma }_{2}x+c_{11}\mathit{sinh}{\gamma }_{3}x+c_{12}\mathit{cosh}{\gamma }_{3}x\right)e^{\mathit{i\omega t}}\end{array}$$ E17

The relations between the constants are given by:

$$\begin{array}{l}{c}_7=k_1{c}_2,c_9=k_2{c}_4,c_{11}=k_3{c}_6\\ c_8=-k_1{c}_1,c_{10}=-k_2{c}_3,c_{12}=k_3{c}_5\\ \begin{array}{c}{k}_1=\frac{C_{23}^k{\gamma }_1^3}{C_{22}^k{\gamma }_1^2-I_k{\omega }^2}\\ k_2=\frac{C_{23}^k{\gamma }_2^3}{C_{22}^k{\gamma }_2^2-I_k{\omega }^2}\\ k_3=\frac{-C_{23}^k{\gamma }_3^3}{C_{22}^k{\gamma }_3^2+I_k{\omega }^2}\end{array}\end{array}$$ E18

3. Natural frequencies and mode shapes

Applying the transfer-matrix method [25], the relation between the state vectors at one end of a segment ${\left\{Z\right\}}_k$ and that at the other end ${\left\{Z\right\}}_{k+1}$ is given by:

$${\left\{Z\right\}}_{k+1}={\left[T\right]}_k{\left\{Z\right\}}_k$$ E19

where ${\left[T\right]}_k$ is called the elementary transfer matrix associated with the k^th segment. For a spar built up of Ns uniform segments, Eq. (19) can be applied at successive stations to obtain:

$$\begin{array}{l}{\left\{Z\right\}}_{\mathit{Ns}+1}={\left[T\right]}_o{\left\{Z\right\}}_1\\ {\left[T\right]}_o={\left[T\right]}_{\mathit{Ns}}{\left[T\right]}_{\mathit{Ns}-1}\dots \dots \dots {\left[T\right]}_k\dots \dots \dots {\left[T\right]}_2{\left[T\right]}_1\end{array}$$ E20

where the matrix ${\left[T\right]}_o$ is known as the overall transfer matrix, which relates the state vectors at the spar fixed end to the free end at which the boundary conditions are specified. Therefore, applying the boundary conditions at both ends and considering only the nontrivial solutions, the frequency equations can be readily obtained. Derivations of the elementary transfer matrices for different vibration modes are discussed in the following section.

3.1 CUS layup

In this case the elements of the state vector at both ends of the k ^th segment are related by the elementary transfer matrix for the coupled extension-twist vibration as follows:

$$\left[\begin{array}{c}\begin{array}{c}{U}_1\left(L_k\right)\\ \phi \left(L_k\right)\\ C_{11}{U}_1^{\prime }\left(L_k\right)+C_{12}{\phi }^{\prime }\left(L_k\right)\end{array}\\ C_{22}{\phi }^{\prime }\left(L_k\right)+C_{12}{U}_1^{\prime }\left(L_k\right)\end{array}\right]={\left[\begin{array}{cccc}{T}_{11}& T_{12}& T_{13}& T_{14}\\ T_{21}& T_{22}& T_{23}& T_{24}\\ T_{31}& T_{32}& T_{33}& T_{34}\\ T_{41}& T_{42}& T_{43}& T_{44}\end{array}\right]}_k\left[\begin{array}{c}\begin{array}{c}{U}_1\left(0\right)\\ \phi \left(0\right)\\ C_{11}{U}_1^{\prime }\left(0\right)+C_{12}{\phi }^{\prime }\left(0\right)\end{array}\\ C_{22}{\phi }^{\prime }\left(0\right)+C_{12}{U}_1^{\prime }\left(0\right)\end{array}\right]$$ E21

Imposing the cantilevered boundary conditions on the overall transfer matrix, the overall frequency equation for coupled extension-twist vibration of the blade spar is given by:

$$T_{33}^o{T}_{44}^o-T_{34}^o{T}_{43}^o=0$$ E22

3.2 CAS layup

Here the transfer matrix equation for the segment k is:

$$\left[\begin{array}{c}\begin{array}{c}-U_3\left(L_k\right)\\ \phi \left(L_k\right)\end{array}\\ \begin{array}{c}{U}_3^{\prime }\left(L_k\right)\\ C_{22}^k{\phi }^{\prime }\left(L_k\right)+C_{23}^k{U}_3^{\prime \prime }\left(L_k\right)\end{array}\\ C_{23}^k{\phi }^{\prime }\left(L_k\right)+C_{33}^k{U}_3^{\prime \prime }\left(L_k\right)\\ C_{33}^k{U}_3^{\prime \prime \prime }\left(L_k\right)\end{array}\right]={\left[\begin{array}{cc}\begin{array}{ccc}{T}_{11}& T_{12}& T_{13}\\ T_{21}& T_{22}& T_{23}\\ T_{31}& T_{32}& T_{33}\end{array}& \begin{array}{ccc}{T}_{14}& T_{15}& T_{16}\\ T_{24}& T_{25}& T_{26}\\ T_{34}& T_{35}& T_{36}\end{array}\\ \begin{array}{ccc}{T}_{41}& T_{42}& T_{43}\\ T_{51}& T_{52}& T_{53}\\ T_{61}& T_{62}& T_{63}\end{array}& \begin{array}{ccc}{T}_{44}& T_{45}& T_{46}\\ T_{54}& T_{55}& T_{56}\\ T_{64}& T_{65}& T_{66}\end{array}\end{array}\right]}_k\left[\begin{array}{c}\begin{array}{c}-U_3\left(0\right)\\ \phi \left(0\right)\end{array}\\ \begin{array}{c}{U}_3^{\prime }\left(0\right)\\ C_{22}^k{\phi }^{\prime }\left(0\right)+C_{23}^k{U}_3^{\prime \prime }\left(0\right)\end{array}\\ C_{23}^k{\phi }^{\prime }\left(0\right)+C_{33}^k{U}_3^{\prime \prime }\left(0\right)\\ C_{33}^k{U}_3^{\prime \prime \prime }\left(0\right)\end{array}\right]$$ E23

Imposing the cantilevered boundary conditions on the overall transfer matrix, the frequency determinant is given by:

$${\left|\begin{array}{ccc}{T}_{44}& T_{45}& T_{46}\\ T_{54}& T_{55}& T_{56}\\ T_{64}& T_{65}& T_{66}\end{array}\right|}_o=0$$ E24

The natural frequencies for coupled bending-twist vibration can be obtained numerically by solving the characteristic determinant of Eq. (24) for $\gamma$ .

4. Optimization model formulation

4.1 Selection of design variables and pre-assigned parameters

In order to formulate a practical and inexpensive optimization model, the large number of design variables of a wind turbine blade has to be reduced to a reasonable number that can be easily dealt with. The design variables which are not subject to change in the optimization process are called the pre-assigned parameters. They are selected to be:

1. Type of materials of construction

2. Total blade length

3. Chord distribution along the blade axis

4. Twist angle distribution along the blade axis

5. Airfoil type and dimensions

6. Dimensions and spacing of the internal supporting ribs

7. Covering skin thickness

8. Blade-to-hub attachment, which is chosen to be of hingeless type

9. Shape of the spar cross section, which is chosen to be rectangular

On the other hand, the design variables, which are subject to change during the optimization process, are chosen to be the fiber orientation angle θ_j,k , fiber volume fraction V_fj,k , thickness of each lamina h_j,k , cross-sectional dimensions (a_k , b_k , H_k ), and length L_k of each segment composing the blade spar (refer to Figure 2). The cross-sectional aspect ratio (a_k /b_k ) can be prescribed according to the airfoil cross-sectional dimensions, and only the variable a_k is taken as a design variable. In addition, to formulate a normalized optimization model with scaled variables, it is necessary to start with a known reference beam to which all the design variables, constraints, and objective function are referred. It is chosen to be made of unidirectional laminated composites with fiber volume fraction V_fo and fiber orientation angle θ_o , and its cross-sectional dimensions are denoted by b_o for width, a_o for depth, H_o for wall thickness, and L_o for length. Therefore, the normalized variables are defined as follows:

$$\begin{array}{l}{\widehat{L}}_k=L_k/L_o\\ {\widehat{b}}_k=b_k/b_o\\ \begin{array}{l}{\widehat{a}}_k=a_k/a_o\\ {\widehat{H}}_k=H_k/H_o=\sum _{j=1}^{\mathit{Nr}}{\widehat{h}}_j\\ {\widehat{h}}_j=h_j/H_o\end{array}\end{array}$$ E25

For a spar beam composed of Ns segments, the design variables are defined by the following matrix equation:

$$\left[X\right]=\left[\begin{array}{cccccc}{\widehat{L}}_1& {\widehat{L}}_2& & {\widehat{L}}_k& & {\widehat{L}}_{\mathit{Ns}}\\ {\widehat{a}}_1{\theta }_{j,1}& {\widehat{a}}_2{\theta }_{j,2}& & {\widehat{a}}_k{\theta }_{j,k}& & {\widehat{a}}_{\mathit{Ns}}{\theta }_{j,\mathit{Ns}}\\ {V_f}_{j,1}& {V_f}_{j,2}& & {V_f}_{j,k}& & {V_f}_{j,\mathit{Ns}}\\ {\widehat{h}}_{j,1}& {\widehat{h}}_{j,2}& & {\widehat{h}}_{j,k}& & {\widehat{h}}_{j,\mathit{Ns}}\\ ⋮& ⋮& \cdots \cdots \cdots & ⋮& \cdots \cdots \cdots & ⋮\\ ⋮& ⋮& & ⋮& & ⋮\\ ⋮& ⋮& & ⋮& & ⋮\\ {\theta }_{\mathit{Nr},1}& {\theta }_{\mathit{Nr},2}& & {\theta }_{\mathit{Nr},k}& & {\theta }_{\mathit{Nr},\mathit{Ns}}\\ {V_f}_{\mathit{Nr},1}& {V_f}_{\mathit{Nr},2}& & {V_f}_{\mathit{Nr},k}& & {V_f}_{\mathit{Nr},\mathit{Ns}}\\ {\widehat{h}}_{\mathit{Nr},1}& {\widehat{h}}_{\mathit{Nr},2}& & {\widehat{h}}_{\mathit{Nr},k}& & {\widehat{h}}_{\mathit{Nr},\mathit{Ns}}\end{array}\right]$$ E26

The total number of design variables in this case is equal to (3 Nr + 2)*Ns, which is a relatively large number of optimization variables.

4.2 Different optimization strategies

Several criteria are used in the dynamic optimization of wind turbine blades; some important ones are the minimal mass design, maximum frequency design, and frequency-placement criteria. In the first strategy, the mass of the blade is minimized while imposing constraints on the blade natural frequencies and strength. In the maximum frequency criterion, the reduction of vibration level is attained by maximizing the natural frequencies of the blade without regard to the complicated stiffness/mass ratio and the exciting frequencies constraints. Higher natural frequencies are favorable for reducing both of the steady-state and transient responses of any structure being excited. The last category of the objective functions to be considered is the placement of the blade natural frequencies, where the main goal is to separate the natural frequencies of the blade from the exciting frequencies to avoid large amplitudes from occurring near the resonant conditions.

4.2.1 Minimal mass design

The minimal mass design optimization problem can be stated as follows:

$$\begin{array}{cc}\text{Find}& \left[X\right]\\ \text{Which}& \text{minimizes}f\left(X\right)=\widehat{m}\\ \begin{array}{c}\text{Subject to}:\\ \text{Strength constraint}\\ \begin{array}{c}\text{Resonance avoidance}\\ \text{Side constraints}\end{array}\end{array}& \begin{array}{l}\\ {\alpha }_T-1\le 0\\ \begin{array}{c}\left(1-{∆}_n\right)\le \frac{{\widehat{\omega }}_n}{{\widehat{\omega }}_n}\le \left(1+{∆}_n\right),n=1,2,\dots \\ \left[X_L\right]\le \left[X\right]\le \left[X_U\right]\end{array}\end{array}\end{array}$$ E27

where $\widehat{m}=m/m_o$ and ${\widehat{\omega }}_n={\omega }_n/{{\omega }_n}_o$ are the dimensionless mass and frequency and $m_o$ and ${{\omega }_n}_o$ are the mass and frequency of the reference beam, respectively. $\left[X_L\right]\wedge \left[X_U\right]$ are the lower and upper bounds imposed on the design variables.

$$\begin{array}{l}\text{Structural mass}:m=2\sum _{k=1}^{\mathit{ns}}\left[\left(\sum _{j=1}^{\mathit{nr}}{\rho }_j{h}_j\right)\left(b_k+a_k-2H_k\right)L_k\right]\\ m_o=2{\rho }_o{H}_o\left(b_o+a_o-2H_o\right)L_o\end{array}$$ E28

The mass densities ${\rho }_o$ and ${\rho }_j$ are calculated according to the volume fractions V_fo and V_fj , respectively. The symbol ${\alpha }_T$ is called the Tsai-Hill rupture coefficient [26], which depends on the ratios between the principle and rupture stresses (more details are given in appendix A). ${\widehat{\omega }}_n$ , n = 1,2,… are called the target (or desired) frequencies of a known baseline design adjusted to be well separated from the exciting frequencies, and the increments ${∆}_n$ are the associated allowable tolerance of each frequency (e.g., ${∆}_n$ = 1%).

4.2.2 Maximum frequency optimization

Minimization of the overall vibration level is one of the most cost-effective solutions for a successful wind turbine design. It enhances other important design goals such as long fatigue life, high stability, and low noise level. Reduction of the overall vibration level can be attained by maximizing the stiffness-to-weight ratio of the wind turbine blade spar. It is well-known that natural frequency is a good indicator of structural stiffness-to-mass ratio. Thus, the optimization problem considered in this investigation will seek maximization of natural frequencies of the blade spar for different modes of vibrations while maintaining its total structural mass lower than or equal to that of a baseline design. Constraints are imposed on the optimization problem such that the blade spar has to be enclosed by the outer blade skin and has the sufficient strength to carry the applied loads without mass penalty.

The maximum frequency criterion, taking into considerations the different modes of vibration, may be cast in the following form:

$$\begin{array}{cc}\text{Find}& \left[X\right]\\ \text{Which}& \text{minimizes}f\left(X\right)=-\sum _n{W}_n{\widehat{\omega }}_n,n=1,2,\dots \\ \begin{array}{c}\text{Subject to}:\\ \text{Strength constraint}\\ \begin{array}{c}\text{Mass constraint}\\ \text{Side constraints}\end{array}\end{array}& \begin{array}{c}\\ {\alpha }_T-1\le 0\\ \begin{array}{c}\widehat{m}-1\le 0\\ \left[X_L\right]\le \left[X\right]\le \left[X_U\right]\end{array}\end{array}\end{array}$$ E29

The function f(X) represents a weighted sum of the nondimensional natural frequencies of successive modes of vibration with $W_n$ the weighting factor of the n ^th frequency ${\widehat{\omega }}_n$ .

4.2.3 Frequency placement

A good design philosophy for vibration reduction is to separate the natural frequencies of the structure from the harmonics of other excitation sources. Thus, natural frequencies of the blade spar should be placed near a target value. Mathematically, the frequency-placement criterion may be expressed by minimizing an objective function constructed from a weighted sum of the squares of the differences between each important frequency ${\widehat{\omega }}_n$ and its target value ${\widehat{\omega }}_n$ as follows:

$$\begin{array}{cc}\text{Find}& \left[X\right]\\ \text{Which}& \text{minimizes}f\left(X\right)=\sum _n{W}_n{\left({\widehat{\omega }}_n-{\widehat{\omega }}_n\right)}^2\\ \begin{array}{c}\text{Subject to}:\\ \text{Strength constraint}\\ \begin{array}{c}\text{Mass constraint}\\ \text{Side constraints}\end{array}\end{array}& \begin{array}{c}\\ {\alpha }_T-1\le 0\\ \begin{array}{c}\widehat{m}-1\le 0\\ \left[X_L\right]\le \left[X\right]\le \left[X_U\right]\end{array}\end{array}\end{array}$$ E30

4.3 Side constraints

Side constraints are imposed on the design variables in order to verify various geometric, manufacturing, or logical reasons. In the present optimization model, these constraints are defined in the following.

4.3.1 Length of blade spar

The total length of the blade is kept equal to that of a known baseline design:

$$\begin{array}{l}\left(\sum _{k=1}^{\mathit{Ns}}{\widehat{L}}_k\right)-1=0\\ {\widehat{L}}_k\ge 0\end{array}$$ E31

4.3.2 Spar cross-sectional dimensions

Since the blade spar is limited by a tapered configuration with a certain type of airfoil cross sections, the height of any segment k having its end at a specific position x_k+1 (refer to Figure 2) must not exceed the outermost height (a_U ) at this position. Figure 4 shows the outermost dimensions of a wind turbine blade spar of length L, where a_r and a_t represent the spar heights at the root and tip locations, respectively. The dimensionless upper limiting height ${\widehat{a}}_U$ is defined by the relation:

Figure 4.

Outermost dimensions of a wind turbine blade spar.

$$\begin{array}{l}{\widehat{a}}_U={\widehat{a}}_r\left(1-\left(1-\frac{a_t}{a_r}\right){\widehat{x}}_{k+1}\right)\\ {\widehat{x}}_{k+1}={\widehat{L}}_1+{\widehat{L}}_2+\dots +{\widehat{L}}_k\end{array}$$ E32

Lower bounds are also imposed on the height of the k ^th segment as a reasonable percentage of the airfoil height near the tip ( ${\widehat{a}}_L=0.5{\widehat{a}}_t$ ). Thus, the inequality constraints imposed on the dimensionless height of each segment of the blade spar are expressed as:

$${\widehat{a}}_L\le {\widehat{a}}_k\le {\widehat{a}}_U,k=\mathit{1},\mathit{2},\dots .\mathit{Ns}$$ E33

Upper bounds on the segment height represent the main obstacle in blade spar optimization since it restricts the spanwise variation in spar height. Note that the width of each segment (b_k ) depends on the height of that segment (a_k ) according to the pre-assigned aspect ratio of the cross section.

In addition, the inequality constraints imposed on the walls thickness of each segment are:

$${\widehat{H}}_L\le {\widehat{H}}_k\le {\widehat{H}}_U,k=\mathit{1},\mathit{2},\dots .\mathit{Ns}$$ E34

4.3.3 Fiber volume fraction and orientation angle within any segment (k)

$$\left(0.2,-\pi /2\right)\le \left(V_{\mathit{fj}},{\theta }_j\right)\le \left(0.8,\pi /2\right),j=\mathit{1},\mathit{2},\dots \mathit{Nr}$$ E35

The optimization models described by Eqs. (27–35) belong to the constrained nonlinear mathematical programming models since the objective functions and most of the constraints contain nonlinear algebraic expressions. Such problems can be solved by a variety of mathematical programming techniques [27, 28]. The sequential quadratic programming (SQP) is one of the most powerful developed and perhaps one of the best methods of optimization. The method has a theoretical basis that is related to the solution of a set of nonlinear equations using Newton’s method and the derivation of simultaneous nonlinear equations using Karush-Kuhn-Tucker (KKT) conditions to the Lagrangian of the constrained optimization problem [29, 30]. More details are given in Appendix B.

5. Behavior of the objective function

Focusing on CUS layup configurations, extensive studies have been carried out, using MATLAB, on the objective functions in order to be able to visualize the unconstrained behavior of natural frequencies and mass of specific thin-walled composite beams. Studies are performed in a two-dimensional design space such that only two design variables are allowed to change while the others are assigned to specific values. All variables are normalized with respect to a known reference beam parameters, as given in Eq. (25).

5.1 Case of two-segment cantilever beam

For a two-segment beam with each segment constructed only from one layer, the matrix of variables (Eq. 26) can be reduced to take the form:

$$\left[X\right]=\left[\begin{array}{c}{\widehat{L}}_1\\ {\widehat{a}}_1\\ \begin{array}{c}{\theta }_{1,1}\\ {V_f}_{1,1}\\ {\widehat{h}}_{1,1}\end{array}\end{array}\begin{array}{c}{\widehat{L}}_2\\ {\widehat{a}}_2\\ \begin{array}{c}{\theta }_{1,2}\\ {V_f}_{1,2}\\ {\widehat{h}}_{1,2}\end{array}\end{array}\right]$$ E36

Note that in this case ${\widehat{H}}_k={\widehat{h}}_j$ . The effects of changing two corresponding variables simultaneously in the two segments are studied in the following.

5.1.1 Effect of changing fiber orientation angles within the two segments ( ${\theta }_1$ , ${\theta }_2$ )

Consider a two-segment beam with a height of the outboard segment half that of the inboard segment but with equal length. Fiber volume fraction and wall thickness of the two segments are the same as that of the reference beam. The following matrix of pre-assigned variables describes this case:

$$\left[X\right]=\left[\begin{array}{cc}\begin{array}{c}0.5\\ 1.5\\ \begin{array}{c}{\theta }_1\\ 0.5\\ 1\end{array}\end{array}& \begin{array}{c}0.5\\ 0.75\\ \begin{array}{c}{\theta }_2\\ 0.5\\ 1\end{array}\end{array}\end{array}\right]$$ E37

Figure 5 shows the developed contours of the normalized fundamental bending frequency in the ( ${\theta }_1$ , ${\theta }_2$ )-design space. It is noticed that the maximum bending natural frequencies occurs at the point $\left({\theta }_1,{\theta }_2\right)=\left(0,0\right)$ , which is an expected natural solution. Figure 6 shows the level curves of the normalized twist-dominated natural frequency. It indicates that the maximum occurs near the design point $\left({\theta }_1,{\theta }_2\right)=\left(45,45\right)$ . On the other hand, the design points ( ${\theta }_1$ , ${\theta }_2$ ) = (0, 0), (90, 90), (0, 90), and (90, 0) give the same values of the frequency representing a local minima solution, as shown in Figure 6. The angle sequence in this case affects the frequency level curves due to the change in heights of the two segments.

Figure 5.

Effect of fiber orientation angles on the normalized fundamental bending frequency (a two-segment cantilevered beam).

Figure 6.

Effect of fiber orientation angles on normalized twist-dominated natural frequency (a two-segment cantilevered beam).

5.1.2 Effect of changing length and height of the outboard segments ( ${\widehat{L}}_2$ , ${\widehat{a}}_2$ )

Consider a two-segment beam with the following matrix of pre-assigned variables. The height and thickness of the inboard segment are equal to those of the reference beam, while thickness of the outboard segment is half that of the reference beam. Nondimensional length of the inboard segment is $1-{\widehat{L}}_2$ in order to keep the total length equal to that of the reference beam:

$$\left[X\right]=\left[\begin{array}{cc}\begin{array}{c}1-{\widehat{L}}_2\\ 1\\ \begin{array}{c}0\\ 0.5\\ 1\end{array}\end{array}& \begin{array}{c}{\widehat{L}}_2\\ {\widehat{a}}_2\\ \begin{array}{c}0\\ 0.5\\ 0.5\end{array}\end{array}\end{array}\right]$$ E38

Figure 7 shows the developed contours of the fundamental bending frequency in the ( ${\widehat{L}}_2$ , ${\widehat{a}}_2$ )-design space. It is seen that the frequency domain is divided into three regions bonded by the reference frequency value. The extreme left and right regions contain two local minima. The central region contains the global maximum. The maximum bending natural frequencies occur near the point ( ${\widehat{L}}_2$ , ${\widehat{a}}_2$ ) ≈ (0.35, 0.3). This point corresponds to a total structural mass of about 70% of the reference beam, which represents a significant gain from the optimization process. This configuration of beam dimensions gives a tapered beam with a decreasing height and length. The same value of maximum natural frequency is also repeated in the design space at the uppermost right region at ( ${\widehat{L}}_2$ , ${\widehat{a}}_2$ ) ≈ (1, 2), but the total mass becomes higher and equal to that of the reference beam.

Figure 7.

Effect of length and height of the outboard segment on normalized fundamental bending frequency (a two-segment cantilevered beam).

From the previous investigations, it is observed that the frequency functions, even though implicit in the design variables, are well behaved and continuous in the selected design spaces. A weighted sum of successive natural frequencies is a good representative for the objective function of the optimization model. Another important observation is that the length and the height of the individual segments are so effective variables in the whole optimization process.

6. Case study

As a case study, the main spar of a medium scale composite blade of a 750-kW horizontal axis wind turbine is optimized. Full description and technical data can be found in Ref. [31]. The blades cross section having NACA 63-218 airfoil is made of E-glass/epoxy composites with properties given in Table 1. The effective length of the spar is 18.33 m extending from 20–98% of the blade length. The spar cross-sectional dimensions at root are 340 mm height and 680 mm width and at the tip 100 mm and 200 mm. The wall thickness at the root and tips are 15 mm and 5 mm, respectively. The chordwise location of the spar cross-sectional center coincides with the blade pitch axis.

6.1 Baseline design

An initial baseline design of the wind turbine blade spar with total mass = 580.55 kg is selected to be composed of 10 segments of equal length. Fiber orientation angle and volume fraction within each segment are pre-assigned to be 0° and 50%, respectively. Structural and dynamic properties of the spar baseline design are given in Tables 2 and 3, respectively.

Property	Epoxy matrix	E-glass fibers
Modulus of elasticity (GPa)	Em = 4.5	E11f = E22f = 74
Modulus of rigidity (GPa)	Gm = 1.6	G12f = 30
Poisson’s ratio	νm = 0.4	ν12f = 0.25
Density (kg/m^3)	ρm = 1200	ρf = 2600
Tensile strength (MPa)	σmr = 90	σ11r ,σ22r (Eq. A-2)
Shear strength (MPa)	τmr = 52

Table 1.

Material properties of the spar structure [4, 31].

Segment no.	Length (mm)	Height (mm)	Width (mm)	Thickness (mm)
1	1833	316	632	14.4
2		292	584	13.3
3		268	536	12.2
4		244	488	11.1
5		220	440	10
6		196	392	8.9
7		172	344	7.8
8		148	296	6.7
9		124	248	5.6
10		100	200	4.5

Table 2.

Dimensions of the baseline design.

Vibration mode (dominated)	Natural frequencies (Hz)
	First mode	Second mode	Third mode
Twist	36.73	64.88	95.18
Extension	94.9	203.5	319.8
Flap bending	1.55	4.90	11.15

Table 3.

Natural frequencies of the baseline design.

The calculated normal and tangential distributed forces (P_z , P_y ) and the corresponding bending moments (M_z , M_y ) applying to the wind turbine blade are shown in Figures 8 and 9, respectively (refer to Appendix C). The wind speed of 60 m/s (survival loading condition) is taken in calculating the applied loads, while the blades were considered in the stationary horizontal position. These forces and moments are implemented in the strength constraints of the present optimization model.

Figure 8.

Shear forces applied to the blade at 60-m/s wind speed.

Figure 9.

Flapwise and chordwise bending moments applied to the blade at 60-m/s wind speed.

6.2 Definition of the reference beam

The reference beam, to which all the design variables and frequencies are normalized, is selected to be a uniform, thin-walled cantilevered spar with rectangular cross section constructed from single unidirectional E-glass/epoxy composite layer (θ_o = 0) and fiber volume fraction of V_fo = 50%. This fiber layup configuration gives maximum natural frequencies for bending and extension modes of vibration while maintaining moderate structural mass (m_o = 580.5 kg) and strength. The dimensions of the reference beam are selected according to the mean cross-sectional dimensions of the spar baseline design (see Table 4).

Dimensions (mm)	Natural frequencies: ${{\omega }_n}_o/2$ (Hz)
	Mode	First	Second	Third
Width bo = 440
Height ao = 220	Twist	17.59	52.76	87.93
Length Lo = 18,330	Extension	61.99	185.97	309.95
Thickness Ho = 10	Bending	0.65	4.07	11.39

Table 4.

Dimensions and natural frequencies of the reference beam.

6.3 Optimization results

The blade spar described in Section 6.1 is optimized according to the three different optimization models under the specified constrains as given in Section 4.2. The method of sequential quadratic programming with MATLAB optimization toolbox is applied to obtain the needed optimal solutions (refer to Appendix B). The starting, lower, and upper bounds for each design variable are given in Table 5. Determination of these values affects the optimization results; thus, several trials were performed in order to select the appropriate values.

6.3.1 Minimum mass design

In this formulation the target frequencies are set equal to the first three flapping frequencies of the baseline design, which are well separated from resonant frequencies. The attained optimal solutions are given in Table 6 with the spar cross sections constructed from balanced and symmetric laminates. The achieved mass saving reached about 25%, where the natural frequencies are found to be too close to those of the baseline design, as given in Table 7. Considerable mass reduction can be observed with much savings for the inboard segments as compared to that for the outboard segments.

Design variable	Lower boundary	Upper boundary	Starting value
Dimensionless length of each segment ${\widehat{L}}_k$	0.0	1	1/Ns
Dimensionless height of each segment ${\widehat{a}}_k$	0.2	2	1
Fiber orientation angle θj,k	−π/2	π/2	0
Fiber volume fraction Vf j,k	0.2	0.8	0.5
Dimensionless thickness of each layer ${\widehat{h}}_{j,k}$	0.01	0.2	0.1

Table 5.

Lower, upper, and starting boundaries of design variables for multiple-segment spar.

Segment no.	Segment length (mm)	Segment cross section (mm)		Fiber volume fraction	No. of layers	Layup
		Height	Thickness
1	280.5	323.8	9.9	0.31	18	(±45/03/±45/03/±45)s
2	183.3	316.7	9.9	0.33	18	(±45/03/±45/03/±45)s
3	200.8	313.8	8.9	0.32	16	(±45/03/±45/02/±45)s
4	534.1	301.4	8.9	0.33	16	(±45/02/±45/03/±45)s
5	851.0	286.3	8.9	0.32	16	(±45/03/±45/02/±45)s
6	1031.7	263.8	7.9	0.32	14	(±45/02/±45/02/±45)s
7	1112.0	235.2	7.8	0.33	14	(±45/02/±45/02/±45)s
8	1272.3	207.4	7.8	0.34	14	(±45/05/±45)s
9	1482.5	188.9	6.8	0.33	12	(±45/05)s
10	1589.4	170.4	6.5	0.33	12	(±45/05)s
11	1273.1	154.0	6.5	0.33	12	(±45/05)s
12	1132.9	138.6	6.4	0.32	12	(±45/05)s
13	1100.9	120.9	5.7	0.33	10	(±45/04)s
14	1051.8	106.3	5.6	0.33	10	(±45/04)s
15	1177.6	94.1	5.2	0.33	8	(±45/03)s
16	880.5	82.4	4.8	0.33	8	(±45/03)s
17	693.2	70.7	4.6	0.33	8	(±45/03)s
18	540.0	59.0	4.5	0.34	8	(±45/03)s
19	456.8	48.0	4.0	0.33	6	(±45/02)s
20	1485.7	44.0	3.6	0.20	6	(±45/02)s

Table 6.

Minimum mass design optimization.

Vibration mode (dominated)	Design	Natural frequencies (Hz)			$\sum _{n=1}^3{W}_n{\omega }_n$
		First mode	Second mode	Third mode
Twist	Baseline	36.73	64.88	94.90	50.49
	Min. mass	36.16	65.62	89.13	49.52
Extension	Baseline	94.90	203.50	319.80	148.00
	Min. mass	89.13	194.88	299.66	139.77
Flapping	Baseline	1.55	4.90	11.15	2.48
	Min. mass	1.53	4.88	11.17	2.45
Mass (kg)	Baseline	580.55
	Min. mass	437.76

Table 7.

Comparison of natural frequencies and mass between baseline and the minimum mass designs.

6.3.2 Maximum frequency design

In order to obtain a blade spar with higher natural frequencies, the maximum frequency optimization model given in Eq. 29 is applied under inequality mass constraint. Table 8 presents the attained optimal solutions for the blade spar with discrete grading model.

Segment no.	Segment length (mm)	Segment cross section (mm)		Fiber volume fraction	No. of layers	Layup
		Height	Thickness
1	842	329	18.7	0.8	34	(±45/07/±45/07/±45)s
2	873	318	16.3	0.8	30	(±45/06/±45/06/±45)s
3	903	306	14.4	0.8	26	(±45/05/±45/05/±45)s
4	919	294	12.3	0.8	22	(±45/04/±45/04/±45)s
5	932	281	11	0.8	20	(±45/03/±45/04/±45)s
6	937	269	8.8	0.8	16	(±45/03/±45/02/±45)s
7	930	257	7.7	0.8	14	(±45/02/±45/02/±45)s
8	908	245	6.6	0.8	12	(±45/04/±45)s
9	882	234	5.5	0.8	10	(±45/04)s
10	861	222	5.5	0.8	10	(±45/04)s
11	875	211	5.5	0.8	10	(±45/04)s
12	936	199	5.5	0.8	10	(±45/04)s
13	1035	185	5.5	0.8	10	(±45/04)s
14	1128	170	5.5	0.78	10	(±45/04)s
15	1094	156	5.5	0.65	10	(±45/04)s
16	1031	142	4	0.51	8	(±45/03)s
17	956	130	4	0.33	8	(±45/03)s
18	709	91	3	0.33	6	(±45/02)s
19	543	61	3	0.33	6	(±45/02)s
20	1037	44	2	0.2	4	(±45/0)s

Table 8.

Maximum frequency design with discrete spanwise material grading.

Figure 10 shows the optimization gains for the different natural frequencies of the blade spar with maximum frequency designs. It is noticed that gains of 48%, 23%, and 34% have been achieved in fundamental natural frequencies for twist, extension, and flapping modes of vibration, respectively.

Figure 10.

Maximum frequency optimization gains.

It can be concluded that a good blade design with maximum natural frequencies should have smaller wall thickness in the outboard portion of the blade spar. The maximum frequency optimization process recommends the segments located inboard to have higher wall thickness and higher fiber content. The spanwise variation in the height of each segment of the blade spar is always restricted by the airfoil envelope, which is the major obstacle in the optimization process.

When the wind turbine is operating, the rotating blades are the main source of vibration. The forcing frequencies are integer multiples of the rotational speed. A common way to represent natural frequencies and search for possible resonances is to plot the Campbell diagram as shown in Figure 11. The intersection of one of the radial lines with one of the system natural frequencies indicates a potential for resonance at the rotor speed corresponding to that point of intersection. The wind turbine considered in this investigation has two rotational speeds; they are 15 and 22 rpm. It is shown that for the baseline and maximum frequency designs, no resonance can occur at these operating speeds.

6.3.3 Frequency-placement design

Figure 11.

Campbell diagram for baseline and maximum frequency designs. (ـــــــBaseline, --- max. ${\sum }_n$ ).

In the minimal mass design, the natural frequencies of the optimized blade spar are constrained to be equal to those of the baseline design, while in maximum frequency optimization, the natural frequencies are entirely maximized with constraints imposed on its structural mass. If it is desired to have a blade spar with specific target frequencies, the frequency-placement optimization model has to be applied. The desired target frequencies used in the present analysis are selected to be 1.9, 6.4, and 12.65 Hz for the first three modes of flapping vibration, respectively. These frequencies are placed away from resonant frequencies at the operating speeds according to the Campbell diagram shown in Figure 12. The attained optimal solutions of the blade spar are given in Table 9, with the flapping frequencies placed at their prescribed target frequencies.

Figure 12.

Campbell diagram for baseline and frequency-placement designs. (ــــــ baseline, ---- target).

Segment no.	Segment length (mm)	Segment cross section (mm)		Fiber volume fraction	No. of layers	Layup
		Height	Thickness
1	753	330	14.5	0.80	26	(±45/05/±45/05/±45)s
2	948	318	12.9	0.80	24	(±45/04/±45/05/±45)s
3	754	308	12.1	0.80	22	(±45/04/±45/04/±45)s
4	787	298	12.3	0.80	22	(±45/04/±45/04/±45)s
5	1008	284	11.9	0.80	22	(±45/04/±45/04/±45)s
6	540	275	11.1	0.79	20	(±45/03/±45/04/±45)s
7	460	271	10.8	0.79	20	(±45/03/±45/04/±45)s
8	256	257	9.9	0.78	18	(±45/07/±45)s
9	216	251	9.9	0.78	18	(±45/08)s
10	186	241	10.9	0.71	20	(±45/09)s
11	1431	244	7.8	0.74	14	(±45/06)s
12	2563	210	6.4	0.80	12	(±45/05)s
13	1930	185	6.5	0.80	12	(±45/05)s
14	1536	165	6.4	0.80	12	(±45/05)s
15	1374	147	5.6	0.80	10	(±45/04)s
16	1310	125	6.4	0.43	12	(±45/05)s
17	585	77	5.3	0.77	10	(±45/04)s
18	185	51	10.9	0.72	20	(±45/09)s
19	431	71	5.5	0.36	10	(±45/04)s
20	1077	44	5.5	0.63	10	(±45/04)s

Table 9.

Frequency-placement design with discrete spanwise material grading.

7. Finite element analysis

The finite element method (FEM) has been demonstrated as a powerful approach which can handle dynamic analysis of laminated composite structures. Nowadays, current commercial finite element software is capable of simulating nonlinearity of many engineering problems. Commercial software also come with advanced preprocessing and post-processing abilities [32] The preprocessing is just a way used for the data input, since the finite element method requires a large amount of data, while the post-processing is another way for presenting the results in the form of deformed shapes and contour maps. The core of the analysis is what occurs in between the two processes. In selecting a finite element analysis software, it is essential to regard the pre- and post-processor, which can directly affect the analysis speed and accuracy. One of the most appropriate finite element analysis software is Femap (finite element modeling and post-processing). Femap is an engineering analysis program developed by Siemens PLM Software that is used to build finite element models of complex engineering problems and view solution results [33]. Femap uses NX Nastran solver to analyze and solve finite element problems.

Femap is an advanced engineering simulation application for creating, editing, and importing/reusing mesh-centric finite element analysis models of complex products or systems. Femap provides powerful data-driven and graphical result visualization and evaluation.

The attained optimum design of the blade spar is modeled using FEM as a tapered beam with rectangular cross section. The root and tip cross sections are drawn using Femap, and then each line in the root is ruled to the corresponding line in the tip in order to generate the required surfaces as shown in Figure 13. The geometric model is divided into 20 segments with lengths according to the maximum frequency design given in Table 8.

Figure 13.

Finite element model of a 20-segment cantilevered spar meshed with 10,000 elements.

Material properties, element type, and thickness of each segment are to be defined. According to Table 8, a number of six different orthotropic materials with different fiber volume fractions and properties of E₁₁ , E₂₂ , G₁₂ , ν₁₂ , and ρ are to be defined and input to Femap. A total of 20 different layups are to be defined due to the change in fiber orientation angles and thickness of each segment within the spar design. Furthermore, 20 element types are to be defined according to the corresponding material and layup properties.

The normal mode eigenvalues analysis using Lanczos modal analysis is used by Femap in order to define the different mode shapes and the associated natural frequencies of the optimized blade spar design. Figure 14 shows the mode shapes with the associated natural frequencies for bending vibration of the optimized blade spar. Table 10 shows a comparison between analytical and Femap natural frequencies of the optimized blade spar design.

Figure 14.

Mode shapes and natural frequencies for vertical bending vibration of the blade spar with maximum frequency design.

Analysis	Natural frequencies (Hz)
	1VBV	2VBV	3VBV
Analytical	2.09	6.87	13.30
FEM	2.17	7.04	14.49
Diff. %	4.7	2.4	8.2

Table 10.

Comparison between analytical and FEM modeling for the maximum frequency design.

8. Conclusions and recommendations

This chapter presents an optimization model for enhancing the dynamic performance of the spar beam of a wind turbine blade. Design variables include the cross-sectional dimensions and material properties variation along the spanwise direction. Three optimization strategies are developed and tested, including the minimal mass design, maximum frequency design, and frequency-placement design, by placing the frequencies at their target values to avoid large amplitudes and resonance occurrence. Side constraints are imposed on the design variables in order to avoid abnormal-shaped optimized configurations. Based on the fact that an exact dynamic analysis of uniform thin-walled beam segment is available and well established, the dynamic analysis of tapered blade spar has been obtained by applying the transfer matrix method to calculate the natural mode of vibrations. The proposed model deals with dimensionless quantities in order to be applicable to thin-walled beams with arbitrary dimensions. Results indicated that the optimization process leads to significant increase of natural frequencies of the optimized spar when compared to the reference or baseline design without mass penalty. Finite element model showed a good agreement with the analytical model developed in this study with a variation of up to 10%. The main conclusions that can be revealed from the present work are:

1. Tapered multiple-segment spar with spanwise material grading gives natural frequencies higher than that of the reference design. However, it is proved that maximization of the fundamental frequency alone does not guarantee maximization of the other higher frequencies. Higher frequencies have been found to have many local minima and maxima in the defined design space.

2. There are optimum design variables of each segment such as the length, height, wall thickness, and fiber content at which the structural dynamic performance can be enhanced. Good designs favor minimum wall thickness and higher fiber volume fraction.

Finally, the analytical model formulated in this chapter can be extended and applied to study the forced dynamic response of a wind turbine blade. Other cross-sectional types of the blade spar, such as D-shape spar, can be considered, and the effects of blade twist, shear deformation, and rotary inertia are to be considered in future studies.

Appendix

The fracture processes induced in fibrous composite materials depend upon the nature of constituents, the architecture of the laminate, and the type of mechanical loading imposed to the laminate. The rupture of fibrous composite materials is the result of one or combined effect of fiber fracture, transverse fracture in matrix, longitudinal fracture in matrix, fracture of fiber-matrix interface, and delaminations [34] (refer to Figure A.1).

Figure A.1.

Schematic of the fracture of a unidirectional fibrous composite at critical values of (a) axial, (b) transverse, and (c) shear stresses.

Figure A.2.

Lamina under loading condition of (a) principal axes and (b) general axes.

Failure theories for composites have been proposed by extending and adapting isotropic failure theories to account for the anisotropy in stiffness and strength of the composites. One of the first fracture criteria applied to an anisotropic materials was introduced by Tsai-Hill theory [26], which is applied to unidirectional lamina under principal axis in-plane loading condition as shown in Figure A.2 (a). The theory states that no fracture will occur in the lamina if ${\alpha }_T\le 1$ such that:

$${\alpha }_T={\left(\frac{{\sigma }_{11}}{{\sigma }_{11r}}\right)}^2+{\left(\frac{{\sigma }_{22}}{{\sigma }_{22r}}\right)}^2-\frac{{\sigma }_{11}{\sigma }_{22}}{{\sigma }_{11r}^2}+{\left(\frac{{\tau }_{12}}{{\tau }_{12r}}\right)}^2$$ EA.1

where the subscript r refers to rupture strength of the material.

Rupture strengths for fibrous composite materials are given by the following approximate relations [34]:

$$\begin{array}{l}{\sigma }_{11r}={\sigma }_{11\mathit{fr}}\left(V_f+V_m{E}_m/E_{11f}\right)\\ {\sigma }_{22r}={\sigma }_{\mathit{mr}}\left(1-\sqrt{4V_f/\pi }\right)\\ {\tau }_{12r}={\tau }_{\mathit{mr}}\end{array}$$ EA.2

Subscripts f and m refer to fiber and matrix materials, respectively.

For a lamina under general loading condition as shown in Figure A.2 (b), the principal stresses are given by the following equation [4, 34]:

$$\begin{array}{l}\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\tau }_{12}\end{array}\right\}=\left[\begin{array}{ccc}{c}^2& s^2& -2\mathit{cs}\\ s^2& c^2& 2\mathit{cs}\\ \mathit{cs}& -\mathit{cs}& \left(c^2-s^2\right)\end{array}\right]\left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{\mathit{xy}}\end{array}\right\}\\ c=\mathit{\text{cosθs}}=\mathit{\text{sinθ}}\end{array}$$ EA.3

The Tsai-Hill failure theory is expressed in terms of a single criterion instead of the multiple subcriteria required in the maximum stress and maximum strain theories. The theory allows for considerable interaction among the stress components. One disadvantage, however, is that it does not distinguish directly between tensile and compressive strengths. The strength must be specified and used according to the given state of stress.

Sequential quadratic programming (SQP) is one of the mainly developed and perhaps one of the best techniques of optimization [29]. The method has a theoretical base that is related to the solution of a group of nonlinear equations using Newton’s method and the derivation of concurrent nonlinear equations using Karush-Kuhn-Tucker (KKT) conditions to the Lagrangian of the constrained optimization problem. The fundamental idea of SQP is to model the optimization problem at the present iterate x_k by a quadratic programming subproblem and to employ the minimizer of this subproblem to identify a new iterate x _{k + 1} [30]. The challenge is to design the quadratic subproblem so that it yields a good step for the constrained optimization problem and so that the overall SQP algorithm has good convergence properties and good practical performance. Possibly the simplest derivation of SQP methods views them as an application of Newton’s method to the KKT optimality conditions for the optimization problem. If the problem is a so-called convex programming problem, that is, f(X) and G_i ( X ), i = 1,…,m, are convex functions, then the KKT equations are both necessary and sufficient for a global solution point. The Kuhn-Tucker equations can be stated as:

$$\begin{array}{l}\nabla f\left(x\right)+\sum _{i=1}^m{\lambda }_i.\nabla G_i\left(x\right)=0\\ {\lambda }_i.G_i\left(x\right)=0,i=1,\dots .,m_e\\ {\lambda }_i\ge 0,i=m_e+1,\dots ,m\end{array}$$ EB.1

The first equation describes a canceling of the gradients between the objective function and the active constraints at the solution point. For the gradients to be canceled, Lagrange multipliers (λ_i, i = 1,…,m) are necessary to balance the deviations in magnitude of the objective function and constraint gradients. Because only active constraints are included in this canceling operation, constraints that are not active must not be included in this operation and so are given Lagrange multipliers equal to 0. This is stated implicitly in the last two Kuhn-Tucker equations.

B.1 Optimization using MATLAB

MATLAB is popular software that is used for the solution of a variety of scientific and engineering problems. The specific toolbox of interest for solving optimization problems is called the optimization toolbox [35]. MATLAB optimization toolbox contains a library of programs or m-files, which can be used for the solution of different optimization problems. A commonly implemented function denoted by fmincon is applied to most constrained objective functions. It is the most suitable function to the optimization problem of this investigation, since it uses sequential quadratic programming (SQP) method.

fmincon starts at “x0” and attempts to find a minimizer “x” of the objective function described in the m-file named “fun” subject to the linear inequalities “ $A\times x\le b$ ” and the linear equalities “ $\mathit{Aeq}\times x=\mathit{beq}$ .” If there are no linear equalities or inequalities, A, b, Aeq, and beq are replaced with “[].” fmincon can define lower and upper bounds on the design variables in “x” so that the solution is always in the range “lb ≤ x ≤ ub.” fmincon can subject the minimization process to the nonlinear inequalities “c(x)” or equalities “ceq(x)” defined in the m-file constraint function named “nonlcon.” fmincon optimizes such that “c(x) ≤ 0” and “ceq(x) = 0.” fmincon minimizes with the optimization options involved in the structure “options.” Table B.1 defines all input and output arguments related to fmincon optimization function.

Argument	Description
fun	The function to be minimized
x0, lb., and ub	Starting, lower, and upper boundaries of design variables
nonlcon	The function that computes the nonlinear inequality constraints “c(x) ≤ 0” and the nonlinear equality constraints “ceq(x) = 0”
fval	Value of the function “fun” at x
exitflag	Integer identifying the reasons causing the algorithm terminated
grad	Gradient at x
hessian	Hessian at x
lambda	Structure containing the Lagrangian multipliers at the solution x
output	Structure containing information about the optimization process such as number of function evaluations and iterations taken and the used optimization technique

Table B.1.

Arguments related to fmincon function [35].

The distributed load vectors are expressed in the undeformed coordinates (see Figure C.1), as follows:

Figure C.1.

Deformed and undeformed configurations of a blade segment.

$$\text{Distributed forces}:\underset{¯}{\mathrm{P}}={\underset{¯}{\mathrm{P}}}_{\mathrm{A}}+{\underset{¯}{\mathrm{P}}}_{\mathrm{I}}+{\underset{¯}{\mathrm{P}}}_{\mathrm{G}}+{\underset{¯}{\mathrm{P}}}_{\mathrm{D}}$$ EC.1

$$\text{Distributed moments}:\underset{¯}{\mathrm{q}}={\underset{¯}{\mathrm{q}}}_{\mathrm{A}}+{\underset{¯}{\mathrm{q}}}_{\mathrm{I}}+{\underset{¯}{\mathrm{q}}}_{\mathrm{G}}+{\underset{¯}{\mathrm{q}}}_{\mathrm{D}}$$ EC.2

where the subscripts A, I, G, and D refer to the aerodynamic, inertial, gravitational, and damping contributions, respectively. The aerodynamic forces P _A and moments q _A can be obtained using the quasi-steady blade-element strip theory [1, 2].

Figure C.2 shows the velocity triangle and the coordinate system used for computing aerodynamic loads. V_yr and V_zr are the tangential and axial velocity components, respectively. The resultant velocity V_r can be calculated from:

Figure C.2.

Velocity triangle of an arbitrary airfoil section in the deformed state.

$${\mathrm{V}}_r\simeq \sqrt{V_{\mathrm{yr}}^2+V_{\mathrm{zr}}^2}=-V_{\mathrm{yr}}\sqrt{1+{\left(\frac{V_{\mathrm{zr}}}{V_{\mathrm{yr}}}\right)}^2}$$ EC.3

The distributed lift and drag force vectors on an arbitrary airfoil section are given by the aerodynamic formulas (see Figure C-2):

$$\begin{array}{c}\text{Lift}\underset{¯}{L}=\frac{1}{2}{\mathrm{\rho }}_a{V}_r^2\mathrm{C}{\mathrm{C}}_L\left(sin{\mathrm{\phi }}_i{\widehat{j}}^{\prime }+cos{\mathrm{\phi }}_i{\widehat{k}}^{\prime }\right)\\ \text{Drag}\underset{¯}{D}=\frac{1}{2}{\mathrm{\rho }}_a{V}_r^2\mathrm{C}{\mathrm{C}}_D\left(-cos{\mathrm{\phi }}_i{\widehat{j}}^{\prime }+sin{\mathrm{\phi }}_i{\widehat{k}}^{\prime }\right)\end{array}$$ EC.4

where ρ_a = air density; C = local chord at spanwise location r; C_L = C_L ( $\alpha$ ) = C_L ${}_{\alpha }\alpha$ , lift coefficient; C_D = C_D ( $\alpha$ ), drag coefficient; C_L ${}_{\alpha }$ = lift – curve slope; $\alpha$ = angle of attack; ${\phi }_i=\text{inflow angle}={\mathrm{\theta }}_B+\alpha$ .

Expressed in the undeformed system, the aerodynamic force vector per unit length of the blade is therefore given as:

$${\underset{¯}{P}}_A=\underset{¯}{L}+\underset{¯}{D}={\mathrm{P}}_{\mathrm{xA}}\widehat{i}+{\mathrm{P}}_{\mathrm{yA}}\widehat{j}+{\mathrm{P}}_{\mathrm{zA}}\widehat{k}$$ EC.5

The distributed inertia loads per unit blade span are obtained by applying D’Alembert’s principle as follows:

$$\begin{array}{l}\text{Inertial force vector}:{\underset{¯}{P}}_{\mathrm{IB}}=-\int \int {\rho }_B{\underset{¯}{\ddot{r}}}_B\mathrm{dy}\mathrm{dz}\\ \text{Inertial moment vector}:{\underset{¯}{q}}_{\mathrm{IB}}=-\int \int {\rho }_B\left[\left({\underset{¯}{r}}_B^{\prime }-r_{E.C}^{\prime }\right)\times {\underset{¯}{\ddot{r}}}_B\right]\mathrm{dy}\mathrm{dz}\end{array}$$ EC.6

where ρ_B is the blade material mass density, ${\underset{¯}{\ddot{r}}}_B$ is the acceleration vector, and ${\underset{¯}{r}}_B^{\prime }$ and ${\underset{¯}{r}}_{E.C}^{\prime }$ are the position vectors of an arbitrary point and the elastic center of the blade cross section, respectively.

Next, considering the equilibrium of a differential element of the deformed blade, the equilibrium equations expressed in the rotating (xyz) axes are:

$$\begin{array}{l}\text{Force equilibrium}:\frac{\partial \underset{¯}{F}}{\partial x}+\underset{¯}{P}=\underset{¯}{0}\\ \text{Moment equilibrium}:\frac{\partial \underset{¯}{M}}{\partial x}+{\widehat{i}}^{\prime }\times \underset{¯}{F}+\underset{¯}{q}=\underset{¯}{0}\end{array}$$ EC.7

$$\begin{array}{l}\text{Internal force vector}:\underset{¯}{F}={\mathrm{F}}_x\widehat{i}+{\mathrm{F}}_y\widehat{j}+{\mathrm{F}}_z\widehat{k}\\ \text{Internal moment vector}:\underset{¯}{M}={\mathrm{M}}_x\widehat{i}+{\mathrm{M}}_y\widehat{j}+{\mathrm{M}}_z\widehat{k}\end{array}$$ EC.8

A wind turbine blade cross section is composed of thin-walled, closed, single or multicellular composite beams. The periphery of each beam cross section is assumed to be constructed of flat composite laminates as shown in Figure C.5. The stress-strain relations of the composite laminates which discretize each cross section are computed using the classical laminate theory CLT. Although each laminate is actually an assembly of multiple layers with different constitutive properties, CLT is used to calculate a set of effective stiffness coefficients that allows a composite laminate to be treated as a single structural element [12]. Therefore the blade cross section may be considered to be composed of discrete sections of homogenous materials.

Generally, for a heterogeneous composite section shown in Figure C.3 and C.4, the modulus-weighted cross-sectional properties are defined as [12]:

Figure C.3.

Blade cross section with discretized composite laminated plates.

Figure C.4.

Cross section for a heterogeneous composite beam.

$$\begin{array}{l}{A}_c=\frac{1}{E_{\mathit{ref}}}\sum _{i=1}^n{E}_i{A}_i\\ y_c=\frac{1}{E_{\mathit{ref}}{A}_c}\sum _{i=1}^n{E}_i{A}_i{y}_i\\ \begin{array}{l}{z}_c=\frac{1}{E_{\mathit{ref}}{A}_c}\sum _{i=1}^n{E}_i{A}_i{z}_i\\ I_y=\frac{1}{E_{\mathit{ref}}}\sum _{i=1}^n{E}_i\left(I_{v_o,i}+A_i{z}_i^2\right)\\ \begin{array}{c}{I}_z=\frac{1}{E_{\mathit{ref}}}\sum _{i=1}^n{E}_i\left(I_{w_o,i}+A_i{y}_i^2\right)\\ I_{\mathit{yz}}=\frac{1}{E_{\mathit{ref}}}\sum _{i=1}^n{E}_i\left(I_{{\mathit{vw}}_o,i}+A_i{y}_i{z}_i\right)\end{array}\end{array}\end{array}$$ EC.9

where $E_{\mathit{ref}}$ is a reference modulus of elasticity, ( $y_i$ , $z_i$ ) denotes the geometric centroid of each discrete element of the cross section, and (v_o,i , w_o,i ) denotes the principal axes of each discrete element.

The parallel axes theorem can be applied to compute the second moments of area about the cross-sectional principal axes as follows:

$$\begin{array}{l}{I}_v=I_y-A_c{\left(z_c\right)}^2\\ I_w=I_z-A_c{\left(y_c\right)}^2\\ I_{\mathit{vw}}=I_{\mathit{yz}}-A_c{y}_c{z}_c\end{array}$$ EC.10

Once the global cross-sectional properties are computed using the method of modulus-weighted properties, the effective axial stress applied to the blade cross section can be given by:

$${\sigma }_x\left(y,z\right)=\frac{F_x}{A_c}-\frac{M_z{I}_v+M_y{I}_{\mathit{vw}}}{I_v{I}_w-{I_{\mathit{vw}}}^2}\left(y-y_c\right)+\frac{M_y{I}_w+M_z{I}_{\mathit{vw}}}{I_v{I}_w-{I_{\mathit{vw}}}^2}\left(z-z_c\right)$$ EC.11

Finally, by converting the distribution of the effective beam stresses into equivalent in-plane distributed loads on the flat laminates composing the cross section, as shown in Figure C.5, the lamina-level strains and stresses can be computed using the CLT.

Figure C.5.

Resultant loads and moments applying to a composite laminate.