Renewable Energy xxx (2014) 1e13
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method Dong Ok Yu, Oh Joon Kwon* Department of Aerospace Engineering, Korea Advanced Institute of Science and Technology, Daejeon 305-701, South Korea
a r t i c l e i n f o
a b s t r a c t
Article history: Received 7 October 2013 Accepted 11 March 2014 Available online xxx
The aeroelastic response and the airloads of horizontal-axis wind turbine rotor blades were numerically investigated using a coupled CFDeCSD method. The blade aerodynamic loads were obtained from a NaviereStokes CFD flow solver based on unstructured meshes. The blade elastic deformation was calculated using a FEM-based CSD solver which employs a nonlinear coupled flap-lag-torsion beam theory. The coupling of the CFD and CSD solvers was accomplished in a loosely coupled manner by exchanging the information between the two solvers at infrequent intervals. At first, the present coupled CFDeCSD method was applied to the NREL 5MW reference wind turbine rotor under steady axial flow conditions, and the mean rotor loads and the static blade deformation were compared with other predicted results. Then, the unsteady blade aerodynamic loads and the dynamic blade response due to rotor shaft tilt and tower interference were investigated, along with the influence of the gravitational force. It was found that due to the aeroelastic blade deformation, the blade aerodynamic loads are significantly reduced, and the unsteady dynamic load behaviors are also changed, particularly by the torsional deformation. From the observation of the tower interference, it was also found that the aerodynamic loads are abruptly reduced as the blades pass by the tower, resulting in oscillatory blade deformation and vibratory loads, particularly in the flapwise direction. Ó 2014 Elsevier Ltd. All rights reserved.
Keywords: Horizontal-axis wind turbine Coupled CFDeCSD method Rotor blade loads Aeroelastic deformation Aerodynamic performance
1. Introduction Because the power captured by wind turbines is proportional to the swept area of the rotor disc, to be a competitive energy resource over other energy generation systems, the overall size of wind turbines has been continuously increased. At the same time, to improve the cost-effective energy efficiency, the turbines are designed to have less weight. This results in more slender, lighter, and therefore more flexible rotor blades. For these slender and flexible blades, the aeroelastic deformation is unavoidable, which leads to vibratory loads, and alters the turbine power performance [1]. The flexible blades may also induce severe instability problems that shorten the operational life of the turbines [2]. These effects become increasingly more significant for offshore wind turbines, because they usually operate at higher wind speeds than those on land and are also frequently exposed to harsh operational environments such as severe storm. Therefore, it is quite important to better understand the aeroelastic behaviors of the rotor blades such
* Corresponding author. Tel.: þ82 42 350 3720; fax: þ 82 42 350 3710. E-mail address:
[email protected] (O.J. Kwon).
that large-scale wind turbines are designed to be a more efficient and reliable energy production system. Recently, a conceptual utility-scale horizontal-axis wind turbine (HAWT) was proposed by the National Renewable Energy Laboratory (NREL) [3], known as the NREL 5MW reference wind turbine. This wind turbine has been adopted extensively by numerous researchers as a reference for assessing the capabilities of their aerodynamic and aeroelastic models [4e8]. In their studies, sophisticated computational structural dynamics (CSD) models were used to accurately describe the dynamic behavior of the rotor blades, whereas relatively simple aerodynamic models mostly based on the blade-element/momentum (BEM) theory were utilized to determine the aerodynamic loads. The BEM-based method is computationally efficient, and enables to provide reasonable estimation of the aeroelastic behavior of flexible blades. However, because of the fundamental assumptions of quasi-steady flows and two-dimensional aerodynamics, the method has inevitable shortcomings for simulating realistic flows, and as a result, is not able to accurately predict the blade loads, particularly near the edge of the operating envelope, and ultimately leads to overdesigned and inefficient wind turbines. To overcome these limitations and to improve the aerodynamic prediction accuracy,
http://dx.doi.org/10.1016/j.renene.2014.03.033 0960-1481/Ó 2014 Elsevier Ltd. All rights reserved.
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
2
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
Nomenclature Cp Cn Ct Cm c Fth, Ftq R r u VN v w
as bp q rN 4
j U
2 Þ surface pressure coefficient, ðP PN Þ=ð1=2rN Vlocal 2 cÞ sectional normal force coefficient, Fn =ð1=2rN Vtip sectional tangential force coefficient, 2 cÞ Ft =ð1=2rN Vtip sectional pitching moment coefficient at 1/4-chord, 2 c2 Þ M1=4 =ð1=2rN Vtip blade chord length at 75% span sectional thrust and torque forces rotor radius blade spanwise location spanwise deflection freestream velocity lead-lag deflection flap bending deflection rotor shaft tilt, degree precone angle, degree rigid pitch angle, degree freestream density torsional deformation at elastic axis, degree rotor azimuth angle rotor rotational speed
In the present study, a coupled CFDeCSD method has been developed for predicting unsteady aerodynamic loads and timevarying aeroelastic deformation of HAWT rotor blades. The blade aerodynamic loads were calculated using an unstructured mesh NaviereStokes CFD solver, and the blade elastic deformation was obtained by a FEM-based CSD solver based on a nonlinear coupled flap-lag-torsion beam theory. The CFD and CSD solvers were calculated simultaneously by exchanging the data in a loosely coupled manner. At first, the coupled method was applied to the 5MW reference wind turbine rotor blades under axial flow conditions, and the predicted steady blade aerodynamic loads and the static structural deformation were compared with the design mean values. Then, time-accurate coupled calculations were performed, and the unsteady blade aerodynamic loads and the aeroelastic deformation due to the rotor shaft tilt and the gravitational force were studied. In the unsteady calculations, both the rotor-alone configuration and the full wind turbine geometry with tower and nacelle were considered to examine the tower interference effect. The results were compared with those of the CFD-only rigid blade calculations to investigate the effects of blade deformation on the rotor blade aerodynamic performances. 2. Coupled CFDeCSD method 2.1. CFD Solver
current state-of-the-art computational fluid dynamics (CFD) techniques based on the Euler and NaviereStokes equations can also be utilized. Methods that couple the two numerical models for tackling fluidestructure interactions are extensively described in the literature [9,10]. The coupled CFDeCSD methods have been frequently used in the rotorcraft community for predicting the vehicle trim state with better accuracy than traditional comprehensive CSD codes based on the lifting-line/BEM theory [11,12]. However, until recently, these coupled techniques have rarely been used for the analysis of wind turbine rotor blades. One of the coupled CFDeCSD analyses for the wind turbine rotor blades was made by Bazilevs et al. [13]. In their study, the 5MW reference wind turbine rotor blade was analyzed using a NaviereStokes CFD solver and a blade CSD solver based on shell formulation. The coupling was achieved by exchanging the information between the two solvers at every time step in a tightlycoupled manner. The calculations were made at the nominal operating condition when the turbine rotor is exposed to the rated wind speed only for a single blade by imposing a rotationally periodic boundary condition between the blades. However, the periodicity does not practically exist, because the rotor is in a tilt with the incoming freestream, and also is under the influence of the gravitational force. The coupled calculations were made only for less than one rotor revolution, and thus the fully-converged solutions were not obtained for the blade aerodynamic loads and the structural response. A coupled CFDeCSD analysis was also performed by Yu, Kwon, and Kwon [14] using a loosely coupled approach. In their study, the calculations were made to obtain converged static blade responses for a wide range of turbine operating conditions. The coupled methodology was utilized as the analysis module in an aerodynamic shape design optimization framework for enhancing the aerodynamic performances of the horizontal-axis wind turbine rotor blades. The coupled analyses mentioned above were made only for the single reference blade with the assumption of rotational periodicity, and the unsteady blade load behaviors and the dynamic structural responses, including the effect of the full wind turbine configuration, have not been studied yet.
For the assessment of the aerodynamic loads, an incompressible NaviereStokes CFD solver based on unstructured meshes was utilized [15]. In the present flow solver, the inviscid fluxes were evaluated using second-order Roe’s flux-difference splitting scheme [16], whereas the viscous flux terms were computed based on central differencing [17]. A dual-time implicit time integration algorithm based on the linearized Euler backward differencing was used to advance the solution in time. To consider the turbulence effect including laminareturbulent transition, the correlationbased transition ku SST model [18] was adopted. An overset mesh technique [19] was used to handle the relative motion between the rotor blades and the turbine sub-components, such as tower and nacelle. The mesh movement due to blade deformation was taken care of by using mesh deformation techniques [20,21]. The effect of the mesh movement was taken into account in evaluating the convective fluxes of both the mean flow and turbulence equations such that the geometric conservation law [22] is satisfied. On the solid wall surface, the no-slip boundary condition was imposed by setting the flow velocity equal to the wall velocity itself. The flow solver was parallelized using MeTiS and MPI libraries.
Fig. 1. Blade coordinate system and blade deformation kinematics.
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
3
Fig. 2. Schematic of CFDeCSD loose-coupling procedures.
2.2. CSD Solver The blade structure was modeled as a nonlinear Euler-Bernoulli cantilever beam undergoing spanwise, lead-lag bending, flap bending, and torsional deformations [23]. This model has been extensively used and validated for studying the aeroelastic behavior of rotor blades with aspect ratio typically larger than 15 [24,25]. In Fig. 1, the blade coordinate system and the blade
Fig. 3. Solid modeling of the NREL 5MW reference wind turbine.
deformation kinematics are presented. To obtain the discretized equations governing the blade elastic motion, a finite-element method was adopted. The blade is divided into a number of finite elements, and each element involves fifteen nodal degrees of freedom [26]. After assembling the elemental matrices and the vectors over the all beam elements, the nonlinear equations of motion for the blade can be written as
Fig. 4. Computational mesh for the CFD calculations for static coupled CFDeCSD analysis.
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
4
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
Fig. 5. Spanwise distributions of static blade deformations at the rated wind speed of 11.5 m/s.
€ þ Gq_ þ Kq ¼ F0 þ FNL þ Fg þ Faero Mq
(1)
where M, G, and K are the mass, gyroscopic, and stiffness matrices, respectively, and F0 and FNL are the constant and nonlinear vectors. Fg and Faero indicate the external loads acting on the blade due to the gravitational force and the aerodynamic loads. This nonlinear dynamic system of equations governing the blade elastic motion is solved using the generalized-a time integration [27]. When the turbine rotor operates under steady axial flow conditions, the blade aerodynamic loads remain constant. Therefore, after neglecting the gravitational force effect, the dynamic system of equations can be reduced to static equilibrium equations as
Kq ¼ F0 þ FNL þ Faero
(2)
In the above dynamic equations of motion and the static equilibrium equations, the aerodynamic loads, Faero, is provided from a separate module. For the coupled CFDeCSD analyses in the present study, the aerodynamic loads are evaluated using the CFD flow solver. The present CSD solver is also equipped with a builtin aerodynamic module based on the blade-element method [26] coupled with linear inflow models. The unsteady effects due to the shed wake are also included in the BEM aerodynamic model by adopting the indicial response method [28]. In this method, from the known indicial aerodynamic load functions that are the responses to a step change of angle-of-attack, the unsteady aerodynamic loads to arbitrary angle-of-attack change are evaluated through the superposition of the indicial responses using Duhamel’s integral. The unsteady loads caused by pitch-rate change can also be calculated in a similar manner. This in-built aerodynamic module is used to provide the aerodynamic loads for the initialization of the coupled analyses and also for a comparison purpose. In the present study, the tower was assumed to be rigid, and thus its flexibility was not considered in the CSD solver.
2.3. Coupling methodology The coupling of the CFD and CSD solvers was made by adopting two different loose-coupling methodologies. To obtain the steady blade aerodynamic loads and the static deformation under axial flow conditions, the CSD solver employs the static equilibrium equation in Eq. (2). The overall procedure of this static coupled CFDeCSD analysis is presented in Fig. 2(a). Initially, the flow simulation is performed for the undeformed blade on the initial mesh. Once a steady-state flow solution is obtained, the blade aerodynamic loads along the undeformed elastic axis are calculated
Fig. 6. Comparison of chordwise surface pressure distributions between static coupled CFDeCSD analysis and CFD-only rigid blade calculation at the rated wind speed of 11.5 m/s.
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
5
procedure is repeated until a converged static equilibrium blade position is achieved. When the turbine rotor operates at a tilt from the axial flow and/ or under the influence of the gravitational force, the incoming flow relative to the rotating blades becomes unsteady, and the blade aerodynamic loads and the deformation change periodically with the azimuthal position. To analyze this time-varying blade behavior, the blade dynamic equations of motion in Eq. (1) are solved, coupled with fully unsteady time-accurate CFD calculations, based on the delta-airload loose-coupling methodology [12]. In this methodology, the blade aerodynamic loads and the blade elastic deformation are exchanged after each rotor revolution on a periodic basis. This loose-coupling approach is regarded as the standard technique in the rotorcraft community to obtain the vehicle trim state and the blade aeroelastic response in an efficient manner by avoiding computationally intensive tightly-coupled CFDeCSD calculations inside the trim loop. The overall procedure of the loose-coupling methodology is presented in Fig. 2(b). At first, the CSD calculation is made using the built-in BEM-based aerodynamic model. At this stage, the initial blade deformations are set to zero. Once a converged periodic blade response is obtained, the deformation data as a function of spanwise location and azimuthal step (qn(r,j)) are transferred to the CFD solver, including the pitch control angles for power control, if necessary. Then, the CFD calculation is performed for the deformed blade. This initial CFD calculation usually requires three rotor revolutions to obtain a physically meaningful flow field in a periodically repeated manner. Then the integrated blade aerodynamic loads are calculated from the surface pressure and skin-friction n ðr; jÞ data for one rotor revdistributions to set-up a completeFCFD olution. This information is provided back to the CSD solver. At the subsequent coupled iterations, the aerodynamic loads are fed back to the CSD solver by taking the CFD airloads and the difference of the BEM loads between the current and previous coupled iterations as Fig. 7. Blade tip deformations and rotor torque with varying wind speed for steady axial flows.
by integrating the surface pressure and skin-friction distributions along the blade chord. The information about the blade aerodynamic loads is transferred to the CSD solver, and the blade deformation is calculated using the static equilibrium equation. Then the blade deformation data is provided back to the CFD solver as a set of three translations (u, v, w), plus a set of three rotations (v0, w0, f: Euler angles), at a number of spanwise locations along the undeformed elastic axis. Then the CFD mesh points on the blade surface are moved to the deformed blade position based on the translational and rotational blade deformation data as:
n n1 n n1 Faero ¼ FCFD þ FBEM FBEM
(4)
The blade deformation in response to this modified blade aerodynamic loads is transferred back to the CFD solver, and the CFD calculation is conducted again by following the specified blade motion. From this stage, the CFD calculation for one rotor revolution is usually sufficient to obtain the periodically-converged blade aerodynamic loads by restarting the calculation from the previously converged solution. This iterative coupling procedure is repeated until the BEM aerodynamic loads remain unchanged between the subsequent iterations. This implies that the converged blade response entirely depends on the CFD aerodynamic loads alone. 3. Results and discussion
8 9 8 9 8 9 < u þ x0 =
The coupled CFDeCSD analyses were performed for the NREL 5MW reference wind turbine rotor [3]. In Fig. 3, the geometry modeling of the full wind turbine configuration for the present CFD calculations is presented. The modeling is based on the data reported by NREL [3], except for some approximations for the hub and the nacelle. The rotor is installed with a shaft tilt angle of 5.0 , and rotates clockwise as viewed from upwind. To increase the tower clearance, the blade is set at a precone angle of 2.5 . Initially, static coupled CFDeCSD calculations were performed for the axial flow conditions for the blade-alone configuration, ignoring the rotor shaft tilt and as a result the flow asymmetry due to the wind velocity component parallel to the rotor disc plane (VNsinas, see Fig. 3). The gravity-induced fluctuating blade loads
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
6
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
Fig. 8. Unstructured overset computational meshes for unsteady time-accurate CFD calculations.
were also ignored in these calculations. The calculations were made for several wind speed cases, and the steady blade aerodynamic loads and the static blade deformation were compared with the design mean values. Next, the unsteady blade aerodynamic loads and the dynamic blade response due to the rotor shaft tilt and the gravitational force were investigated through unsteady timeaccurate coupled calculations. Both the rotor-alone configuration and the full wind turbine geometry with tower and nacelle were considered to examine the tower interference effect. At first, unsteady CFD-only calculations were performed for the rigid blade as the reference for investigating the effects of the rotor shaft tilt and the tower interference. Then, the effect of the gravitational force on the blade deformation was examined through CSD-only calculations. Finally, time-accurate coupled CFDeCSD calculations were performed for both the rotor-alone and full wind turbine configurations, and the unsteady aerodynamic loads and the blade dynamic response were investigated.
3.1. Static deformation in steady axial flows At first, the static coupled CFDeCSD analysis was made for the blade-alone configuration at the rated wind speed of 11.5 m/s. At this wind speed, the blade rotates with an angular velocity of 12.1 RPM, and the pitch control angle is set to zero degree. The flow calculation was made for a reference blade by applying the periodic boundary condition between the other two blades as shown in Fig. 4. The computational mesh contains 3,213,240 prismatic cells and 2,712,271 tetrahedral cells, and the number of nodes is 2,100,244. For the CSD calculations, 49 finite elements were used along the span to model the structure of the blade. The elastic axis was assumed to be located at 1/4-chord as reported by NREL [3]. The coupled calculation was performed until the change in the blade deformation becomes less than 0.1%. The solution converged within five coupled iterations. The elapsed calculation time was approximately 30,060 s per coupled iteration using 120 processors
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
7
Fig. 9. Comparison of azimuthal variations of sectional airloads and pitching moment for rigid blades between rotor-alone and full wind turbine configurations.
on 2.8 GHz CPUs. In Fig. 5, the converged spanwise distribution of the blade deformation is presented. It is observed that the magnitude of the blade deformation is fairly large, particularly for the torsional deformation exhibiting approximately 3.1 (nosedown) at the blade tip. The tip bending deflections are approximately 0.01R and 0.075R in the lead-lag and flapwise directions, respectively. These values are very similar to the design properties of the reference wind turbine predicted by the FAST-Aerodyn method [3]. To examine the effect of the blade deformation on the blade aerodynamic loads, the predicted surface pressure distributions by the coupled CFDeCSD analyses are compared with those of the CFD-only rigid blade calculations in Fig. 6. In the figure, the locations of the center of pressure for the both cases are also presented. It is shown that the blade deformation leads to a significant reduction in the leading-edge suction pressure peak, particularly for the blade sections at the far outboard. This loading change is directly related to the reduced effective angle-of-attack caused by the nose-down torsion indicated in Fig. 5. It is also shown that the center of pressure is located mostly behind the elastic axis (1/4chord), causing the nose-down torsion. The pressure centers for the
deformed blade sections are moved slightly further downstream due to the reduced suction pressure peak. From these results, it is obvious that for large-scale wind turbines, the blade aeroelastic deformation has non-negligible influences on the blade aerodynamic loads, and thus this effect should be considered properly for accurately predicting the rotor performance. To further investigate the effect of blade deformation on the blade aerodynamic loads, additional coupled CFDeCSD analyses and CFD-only calculations were performed at the wind speeds of 4, 6, 8, 10, 15, and 20 m/s. In Fig. 7, the predicted blade tip deformations and the rotor torque are compared with other predictions by FAST-Aerodyn [3] and BLADMODE based on BEM [29]. The BLADMODE results were obtained for an almost identical 62.6 m-long rotor blade configuration that was used as a basis model for designing the NREL 5MW wind turbine. The nose-down torsion predicted by BLADMODE is much less than the present results. This is presumably because the aerodynamic pitching moment was not accurately predicted, which is a typical shortcoming of low-order BEM-based aerodynamic models. In the case of torque, the CFD-only calculations result in higher rotor aerodynamic loads than the coupled analysis results, indicating the
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
8
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
Fig. 10. Comparison of rotor inflow velocity distributions at the cutting plane located 1.5 chord length upstream from rotor with the reference blade at 150 azimuth angle.
significance of considering the effect of blade deformation, particularly for the nose-down torsion. The aeroelastic deformation effect is magnified further at higher wind speeds. The present coupled CFDeCSD results compare very well with the FAST-Aerodyn. 3.2. Dynamic deformation with flow unsteadiness Next, the aeroelastic deformation behaviors of the rotor blades were investigated by considering the flow unsteadiness due to rotor shaft tilt, aerodynamic interference with the tower, and the gravitational force. The analyses were made when the turbine operates at the rated wind speed of 11.5 m/s. At first, the unsteady blade aerodynamic loads due to rotor shaft tilt and tower interference were calculated in a time-accurate manner for the rigid blades without blade deformation. For these CFD-only calculations, both the rotor-alone configuration and the full wind turbine model with tower and nacelle were considered. The unstructured overset meshes for the two configurations are presented in Fig. 8. It is shown that fine cells are concentrated around the rotor disc plane and near the turbine components to capture the rotor wake and the aerodynamic interference with the tower more accurately. For the rotor-alone configuration, the computational mesh consists of
approximately 8.6 M tetrahedrons and 5.3 M prismatic elements, and the total number of nodes is about 4.2 M. For the full configuration after including the nacelle and the tower, the number of elements is slightly increased to approximately 8.7 M and 5.4 M for the tetrahedral and prismatic elements, respectively. The number of nodes is also increased to about 4.3 M. In Fig. 9, the azimuthal variations of the normal and tangential force coefficients, and the pitching moment about the quarterchord are compared between the rotor-alone configuration and the full wind turbine model for the rigid blades without aeroelastic deformation. The results presented at the three selected spanwise locations (0.45, 0.64, and 0.93 R) confirm the unsteady periodic behavior due to the rotor shaft tilt. The normal force coefficient is characterized by the higher magnitude in the azimuthal region between zero degree and 180 , where the blade is advancing toward the wind component aligned tangent to the rotor disc plane (VNsinas). In the remaining region where the blade is retreating from the tangential wind component, relatively lower values are observed. This is because the advancing blade encounters higher dynamic pressure, even at relatively low angle-of-attack flows, whereas the blade in the retreating region is under the opposite flow situation. In contrast, the behavior of the tangential force
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
Fig. 11. Blade tip deformation in response to gravitational force by CSD-only calculations.
coefficient is dictated more by the angle-of-attack change, showing lower values at the advancing blade side and higher ones at the retreating side. Because the center of pressure is mostly located behind the quarter-chord at this nominal operating condition as shown in Fig. 6, negative pitching moments are predicted at all azimuthal positions. This behavior is consistent with the normal force coefficient change. It is shown that the unsteady behaviors are almost identical between the rotor-alone and full wind turbine
9
configurations. However, an abrupt change in the aerodynamic loadings is observed near the 180 azimuth angle, which is due to the interference of the blades with the tower. To better understand the effect of the aerodynamic interference with the tower, the rotor inflow velocity distribution is examined. In Fig. 10(a) and (b), the instantaneous distributions of the Vz-velocity component normal to the rotor disc at the cutting plane located 1.5 chord lengths upstream from the rotor are compared between the rotor-alone configuration and the full wind turbine model, when the reference blade is located at the azimuth angle of 150 . The difference of the Vz-velocity between the two configurations is presented in Fig. 10(c). It is shown that in the case of the full wind turbine configuration, a noticeable reduction of the Vzvelocity occurs near the tower, indicating that the tower has a nonnegligible influence by locally reducing the blade effective angleof-attack, even for the upwind operating condition. Because of this tower blockage effect, the blade pressure loading is reduced near the tower, and as a result the aerodynamic loads are abruptly decreased as indicated in Fig. 9. As a next step, the effect of the gravitational force alone on the blade elastic deformation was investigated by using the CSD-only calculations, without including the aerodynamic loads. The calculations were made for the rotor-alone configuration at the rated rotational speed. The predicted tip deformation presented in Fig. 11 shows that the blade experiences a fairly large azimuthal variation in all directions, particularly for the inplane lead-lag bending deflection, indicating that the influence of the gravitational force is non-negligible. The maximum lead bending deflection occurs near the azimuth angle of 90 , whereas the negative maximum appears near the azimuth angle of 270 . A very similar tendency is also observed for the flap bending deflection. This is because the gravitational force is aligned tangent to the blade section plane, and thus the lead-lag and flap bending degrees of freedom are directly influenced, as shown in Fig. 12. Due to the chordwise offset of the mass center (c.g.) away from the elastic axis (e.a.), the torsional deformation is also similarly influenced. In contrast, when the blade is located at the top and bottom positions, the direction of the gravitational force is aligned along the blade span, and thus the spanwise degree of freedom is more affected. However, because the spanwise tensile stiffness is relatively very high, the corresponding deformation is inherently small. Finally, the delta-airload coupled CFDeCSD analyses were conducted by simultaneously considering rotor shaft tilt, tower and the gravitational force, and their combined effects on the unsteady aerodynamic load behavior and the dynamic blade responses were investigated. The calculations were made at the rated wind speed of 11.5 m/s for both the rotor-alone configuration and the full wind turbine model. The coupled iterations were conducted until the difference of the BEM aerodynamic loads between the subsequent iterations was reduced below 0.5%. The calculations were performed on 2.93 GHz CPUs, which took approximately 81,780 s per coupled iteration using 184 processors. A converged set of the periodic blade response and the CFD-based blade aerodynamic loads was obtained after eight coupled iterations for both configurations. In Fig. 13, variations of the blade tip deformations are presented for one rotor revolution for both the rotor-alone and full wind turbine configurations. In the figure, the coupled CFDeCSD results for the rotor-alone configuration without rotor shaft tilt are also presented to identify the effect of shaft tilt, along with the static deflections without the gravitational force presented in Fig. 5. It is shown that the mean values of the predicted tip deformations are very similar to the static values. It is also shown that the periodic variations of the blade deformation induced by the gravitational force are fairly large, even for the axial flow case without the shaft
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
10
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
Fig. 12. Schematic showing the effect of gravitational force on blade deformation.
tilt. This unsteady behavior is further magnified by adding the effect of the rotor shaft tilt. Also, due to the interference with the tower, the aeroelastic deflections undergo an abruptly change near the azimuth angle of 180 , particularly for the flap bending, as shown for the full configuration. The variation of the lead-lag bending deflection is almost identical to that induced by gravity alone in Fig. 11, indicating that the lead-lag bending behavior is mostly dictated by the gravitational force. The nose-down torsion is higher between zero degree and 180 azimuth angles, and is lower in the remaining region. This variation changes the effective angles of attack at the blade sections. In Fig. 14, the predicted wake structure and the tower clearance are compared between the CFD-only and coupled CFDeCSD calculations when the reference blade is located at the azimuth angle of 180 . It is shown that in the case of the coupled prediction, because of the flap-up bending deflection, the overall positions of the roll-up tip vortices are located further downstream, compared to the CFD-only rigid-blade calculation. Also, the blade tip-to-tower clearance is significantly reduced from 0.2 R to 0.12 R, which matches well with the intended design parameters of the wind turbine with the flexible rotor blades. In Fig. 15, the azimuthal variations of the normal and tangential force coefficients and the 1/4-chord pitching moment are compared between the rigid-blade CFD-only and coupled CFDeCSD calculations at three selected spanwise locations. It is shown that in the case of the flexible blades, due to the reduced angle-of-attack caused by the nose-down torsional deformation, the overall magnitudes of the normal and tangential force coefficients are significantly reduced, particularly at the far outboard. The normal force also changes consistently with the nose-down torsion behavior, indicating that the effect of the angle-of-attack change is dominant. This is in contrast to the fact that the dynamic pressure dominates the unsteady behavior for the rigid blades. The magnitude of variation of the tangential force is also magnified by the torsion behavior. The pitching moment is slightly increased in the negative direction, which is due to the downstream movement of the center of pressure caused by the reduced suction pressure peak. It is also shown that the overall effect of the tower interference appears as the manifested oscillatory load behaviors, compared to the rigid blades.
In Fig. 16, the azimuthal variations of the predicted rotor thrust and torque are compared between the rigid-blade and the coupled CFDeCSD calculations for the full wind turbine configuration. In the figure, the thrust and torque with the means removed and their harmonic decompositions through the FFT analysis are also presented. Because of the reduced blade aerodynamic loads by the torsional deformation, the rotor thrust and torque are also significantly reduced by approximately 11% and 6%, respectively. The vibratory loads for the rigid blade are predominately at 3/rev for both the thrust and torque, which indicates the effect of the tower interference for the three-bladed rotor. In the case of the flexible blades, the amplitudes of the high-frequency loads (6/rev w 12/rev) are increased. This is because of the oscillatory blade loads caused by the tower interference as shown in Fig. 15. This vibratory rotor loads may lead to a poor power supply quality, and may reduce the operational life of the critical turbine components. In Fig. 17, the predicted blade root bending and torsion moments are compared between the rotor-alone and full wind turbine configurations for the flexible blades. It is shown that because of the rotor shaft tilt and the gravitational force, the predominant frequency of the root structural loads is 1/rev in all directions. The tower interference effect produces additional higher harmonics, particularly in the flapwise direction. 4. Conclusion In the present study, a coupled CFDeCSD method has been developed for predicting the aerodynamic loads and the aeroelastic deformation of horizontal-axis wind turbine rotor blades. The blade aerodynamic loads were obtained by a NaviereStokes CFD flow solver based on unstructured meshes. The blade elastic deformation was calculated using a FEM-based CSD solver employing a nonlinear coupled flap-lag-torsion beam theory. The coupling of the two solvers was accomplished by adopting a loose-coupling approach. The applications of the present method were made for the NREL 5MW reference wind turbine under various operating conditions. Initially, the static coupled CFDeCSD analyses were made for the blade-alone configuration in steady axial flows. It was observed that for the flexible blades, a fairly large deformation
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
11
Fig. 14. Comparison of wake structure and tower clearance between rigid blade and coupled CFDeCSD calculations when the reference blade is at the azimuth angle of 180 .
Fig. 13. Azimuthal variations of blade tip deformations.
occurs, particularly for torsion. The maximum blade deformation occurs at the rated operating condition of 11.5 m/s. The present results for the mean rotor aerodynamic performances and the static blade aeroelastic deformations compare well with the design values. It was found that the blade aerodynamic loads are significantly reduced by the nose-down torsion, and thus the aeroelastic blade deformation should be properly considered for the accurate prediction of the wind turbine aerodynamic performances. Next, the time-accurate coupled CFDeCSD analyses were made based on the delta-airload correction method. The calculations were performed for both the rotor-alone configuration and the full wind turbine geometry with tower and nacelle. It was found that the rotor shaft tilt and the gravitational force cause a fairly large unsteady variation of the blade elastic deformations, and that the
changes in the blade aerodynamic loads associated with the blade deformations are dominated mostly by the nose-down torsion. When the turbine tower was considered, an abrupt reduction of the aerodynamic loads was observed as the blades pass by the tower. The aerodynamic interference with the tower causes oscillatory behaviors of both the blade aerodynamic loads and the aeroelastic deformation, particularly for the tangential force and the flapwise bending deflection. The tower interference also produces high-frequency vibratory loads at the blade root. The lead-lag bending behavior is dominantly influenced by the gravitational force. It was found that the present coupled CFDeCSD method is well established, and is a useful tool for accurately predicting the blade unsteady aerodynamic loads and the aeroelastic deformations of horizontal-axis wind turbines under the influence of the rotor shaft tilt, the gravitational force, and the interference with other turbine components that are beyond the prediction capability of the traditional BEM-based methods.
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
12
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
Fig. 15. Comparison of azimuthal variations of normal and tangential force coefficients, and quarter-chord pitching moment at three spanwise locations between rigid blade and coupled CFDeCSD calculations.
Fig. 16. Comparison of azimuthal variations of rotor thrust and torque, and their harmonic decompositions between rigid blade and coupled CFDeCSD calculations for full wind turbine configuration.
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033
D.O. Yu, O.J. Kwon / Renewable Energy xxx (2014) 1e13
13
Fig. 17. Comparison of azimuthal variations of blade root bending moments and their harmonic decompositions between rotor-alone and full wind turbine configurations for coupled CFDeCSD calculations.
Acknowledgments This work was supported by grant No. EEWS-2014-N01140049 from EEWS Research Project of the KAIST EEWS Research Center. The authors also would like to acknowledge the support from the National Research Foundation of Korea, Grant No. 20090083510, funded by the Korean Government Ministry of Education & Science Technology (MEST) through the Multi-phenomena Computational Fluid Dynamics (CFD) Engineering Research Center.
References [1] Ahlstrom A. Influence of wind turbine flexibility on loads and power production. Wind Energy 2006;9(3):237e49. [2] Hansen MH. Aeroelastic instability problems for wind turbines. Wind Energy 2007;10(6):551e77. [3] Jonkman JM, Butterfield S, Musial W, Scott G. Definition of a 5-MW reference wind turbine for offshore system development. NREL/TP-500e38060; 2009. [4] Kallesøe BS. Equations of motion for a rotor blade, including gravity, pitch action and rotor speed variations. Wind Energy 2007;10(3):209e30. [5] Bir G, Jonkman J. Aeroelastic instabilities of large offshore and onshore wind turbines. J Phys Conf Ser 2007;75(1):1e19. [6] Riziotis VA, Voutsinas SG, Politis ES, Chaviaropoulos PK, Hansen AM, Madsen HA, et al. Identification of structural non-linearities due to large deflections on a 5MW wind turbine blade. In: European wind energy conference and exhibition; 2008. [7] Jeong MS, Yoo SJ, Lee I. Aeroelastic analysis for large wind turbine rotor blades. AIAA-Paper; 2012:1950. [8] Xudong W, Shen WJ, Zhu JN, Sørensen JN, Jin C. Shape optimization of wind turbine blades. Wind Energy 2009;12(8):781e803. [9] Viré A, Xiang J, Milthaler F, Farrell PE, Piggott MD, Latham JP, et al. Modeling of fluid-solid interactions using an adaptive mesh fluid model coupled with a combined finite-discrete element model. Ocean Dyn 2012;62(10e12):1487e 501. [10] Hou G, Wang J, Layton A. Numerical methods for fluid-structure interaction e a review. Commun Comput Phys 2012;12(2):337e77. [11] Lim JW, Strawn RC. Computational modeling of HART II blade-vortex interaction loading and wake system. J Aircr 2008;45(3):923e33. [12] Potsdam M, Yeo HS, Johnson W. Rotor airloads prediction using loose aerodynamic/structural coupling. J Aircr 2008;45(3):732e42.
[13] Bazilevs Y, Hsu MC, Kiendl J, Wuchner R, Bletzinger KU. 3D simulation of wind turbine rotors at full scale. Part II: fluid-structure interaction modeling with composite blades. Int J Numer Methods Fluids 2011;65(1e3):236e53. [14] Yu DO, Kwon HI, Kwon OJ. Performance enhancement of HAWT rotor blades by aerodynamic shape optimization. AIAA-Paper; 2012:1292. [15] Yu DO, You JY, Kwon OJ. Numerical investigation of unsteady aerodynamics of a horizontal-axis wind turbine under yawed flow conditions. Wind Energy 2013;16(5):711e27. [16] Roe PL. Approximate Riemann solver, parameter vectors and difference. J Comp Physiol 1981;43(2):357e72. [17] Mathur SR, Murthy JY. A pressure-based method for unstructured meshes. Numer Heat Transf 1997;Part B, 31(2):195e215. [18] Langtry RB, Menter FR. Correlation-based transition modeling for unstructured parallelized computational fluid dynamics codes. AIAA J 2009;47(12): 2894e906. [19] Jung MS, Kwon OJ. A parallel unstructured hybrid overset mesh technology for unsteady viscous flow simulations. In: International conference on computational fluid dynamics; 2007. [20] Kholodar DB, Morton SA, Cummings RM. Deformation of unstructured viscous grids. AIAA-Paper; 2005:0926. [21] Bottasso CL, Detomi D, Serra R. The Ball-Vertex method: a new simple spring analogy method for unstructured dynamic meshes. Comput Methods Appl Mech Eng 2005;194(39):4244e64. [22] Nkonga B, Guillard H. Godunov type method on non-structured meshes for three-dimensional moving boundary problems. Comput Methods Appl Mech Eng 1994;113(1):183e204. [23] Hodges DH, Dowell EH. Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades. NASA TN D-7818; 1974. [24] Datta A, Chopra I. CFD/CSD prediction of rotor vibration loads in high-speed flight. J Aircr 2006;43(6):1698e709. [25] Smith MJ, Lim JW, van der Wall BG, Baeder JD, Biedron RT, Boyd DD, et al. An assessment of CFD/CSD prediction state-of-the-art using the HART II International Workshop data. In: The 68th annual forum of the American Helicopter Society; 2012. [26] Bir G, Chopra I, Ganguli R, Smith EC, Vellaichamy S, Wang J, et al. University of Maryland Advanced Rotorcraft Code (UMARC). Theory Man; 1994. [27] Chung J, Hulbert GM. A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-a method. J Appl Mech 1993;60(2):371e5. [28] Leishman JG. Validation of approximate indicial aerodynamic functions for two-dimensional subsonic flow. J Aircr 1988;25(10):914e22. [29] Lindenburg C. Aeroelastic modeling of LMH63-5 blade. DOWEC-02-KL-083/0; 2002.
Please cite this article in press as: Yu DO, Kwon OJ, Predicting wind turbine blade loads and aeroelastic response using a coupled CFDeCSD method, Renewable Energy (2014), http://dx.doi.org/10.1016/j.renene.2014.03.033