موفق ومحمد وفاروق Al-Khwarizmi Engineering Journal Al-Khwarizmi Engineering Journal, Vol. 9, No. 3, P.P. 15- 25 (2013) Aeroelastic Behavior of a Wind Turbine Blade by a Fluid -Structure Interaction Analysis Mauwafak A. Tawfik* Mohammed I. Abu-Tabikh** Farouk O. Hamdoon*** *,** Department of Machines and Equipment Engineering / University of Technology *** Department of Mechanical Engineering / College of Engineering / Wassit University (Received 10 March 2013; Accepted 29 May 2013) Abstract In this paper, a numerical model for fluid-structure interaction (FSI) analysis is developed for investigating the aeroelastic response of a single wind turbine blade. The Blade Element Momentum (BEM) theory was adopted to calculate the aerodynamic forces considering the effects of wind shear and tower shadow. The wind turbine blade was modeled as a rotating cantilever beam discretized using Finite Element Method (FEM) to analyze the deformation and vibration of the blade. The aeroelastic response of the blade was obtained by coupling these aerodynamic and structural models using a coupled BEM-FEM program written in MATLAB. The governing FSI equations of motion are iteratively calculated at each time step, through exchanging data between the structure and fluid by using a Newmark’s implicit time integration scheme. The results obtained from this paper show that the proposed modeling can be used for a quick assessment of the wind turbine blades taking the fluid-structure interaction into account. This modeling can also be a useful tool for the analysis of airplane propeller blades. Keywords: Wind turbine, Aeroelasticity, FSI model, FEM, Blade element momentum. 1. Introduction The interaction between fluids and structures play an important role in a number of fields. Important applications can be found in airplane wings, blade machines, tall buildings, suspension bridges and biomechanics (e.g. elastic artery modeling and simulation of flow in blood vessels). The interaction between flowing fluids and vibrating structure is the main subject of the aeroelastisity. The flow induced vibration may affect negatively the operation and the response of the system. Aeroelasticty of wind turbine blade is one of main considerations for the design of wind turbine, because aerodynamic loading causes blade to bend mostly in flapwise direction, and causes blade section to twist to create new fluid fields surrounding the blade. This interaction between aerodynamics and deformation of wind turbine blade may lead to twist to create new fluid fields surrounding the blade. This interaction between aerodynamics and deformation of wind turbine blade may lead to aeroelastic problem. In fact, a lot of studies related to the aeroelasticity of wind turbine, widely used Computational Fluid Dynamics (CFD) codes to determine the aerodynamic loading [1, 2 and 3]. These codes require a large numerical time and powerful computers because the solution using these codes may take days to completely converge. On the other hand, using a full 3D model for the blade to compute the structural deformation is computationally too expensive, and hence the blades are treated as series of 1D beam elements in many aeroelastic codes [4]. The aim of this paper is to develop an alternative computationally cheaper and thus fast aeroelastic code for investigating the aeroelastic response of a single wind turbine blade. For this purpose, a simplified numerical modeling for fluid Mauwafak A. Tawfik Al-Khwarizmi Engineering Journal, Vol. 9, No. 3, P.P. 15- 25 (2013) 16 – structure interaction analysis is presented. The BEM method is adopted to calculate the aerodynamic forces considering wind shear and tower shadow. Structural dynamic modal is developed based on Hamilton’s principle combining with FEM. The rotating blade is treated as 1D elastic beam subject to flapwise bending and torsion. A computational flow chart for the fully fluid-structure coupling is presented where the governing aeroelastic equations are solved iteratively within each time step. 2. Aerodynamic Loading Model Aerodynamic loads acting on a wind turbine blade is modeled based on unsteady BEM theory and Prandtl’s tip loss factor. This method is very fast and provides accurate results provided that reliable airfoil data are available for the lift, drag and moment coefficients as a function of the angle of attack. 2.1. BEM Theory The algorithm of unsteady BEM theory given in [5,6] is adopted in this work to compute the aerodynamic loads. Once unsteady BEM algorithm is applied, the tangential Py and normal Pz distribution of the loads acting on each blade at each radial position, at each time step is computed from lift L and drag D forces resulting from the inflow on the blade using the relative velocity and angle of attack values as follows: = | | ( ) …(1) = | | ( ) …(2) = (∅) − cos (∅) ...(3) = cos(∅) + sin (∅) …(4) Also aerodynamic pitching moment about the aerodynamic center is calculated as follows: = | | ( ) …(5) Where Cl, Cd and Cm are lift, drag and pitching moment coefficients respectively, obtained from a lookup airfoil data table depending on the angle of attack.Figure (1) shows the velocity triangle seen locally on a blade element. The vector representation of relative velocity can be written as: ⃑ = ⃑ + ⃑ + ⃑ − ⃑ …(6) Fig. 1. Velocity Triangle Seen Locally on a Blade Element [5]. From the velocity triangle, the angle of attack, α, can be computed if the induced velocity, W and the velocity of the blade section are known: = ∅ − ( + ) …(7) Where (β + θp) is the pitch angle plus the twist of the blade at each station, and the flow angle, , is defined as: ∅ = , , …(8) Knowing α, the lift, drag and moment coefficients can be looked up in airfoil data table. The equations that are used in unsteady BEM model of the present work for normal and tangential induced velocities can be derived from Glauert’s relation between thrust and induced velocity [5, 6]. = ∅ | . ( . )| …(9) = ∅ | . ( . )| …(10) Where fg refers to as the Glauert correction, is an empirical relation between the thrust coefficient and axial induction factor. It is assumed that the induction factor does not exceed 0.2 and hence fg equals 1 [6]. Structural deformation for each blade station is used to calculate the velocity of the blade section and hence the relative velocity for each element of the blade of the wind turbine is computed. 2.2. Wind Profile Besides reliable airfoil data for the different blade section, also realistic varying wind field are required as input to an aeroelastic calculation of a wind turbine [5]. Wind speed generally increases Mauwafak A. Tawfik Al-Khwarizmi Engineering Journal, Vol. 9, No. 3, P.P. 15- 25 (2013) 17 with height, see Figure (2). Therefore, a blade pointing upwards would encounter wind speeds greater than a blade pointing downwards. The wind shear is modeled in this work by the widely used power law as [6]: V (X) = V (H) ( ) …(11) Exponent v increases with the surface roughness. It is in the range between 0.1 and 0.25 [6]. The distribution of wind is also altered due to the presence of tower. For an upwind turbine, when the blade is directly in front of the tower, it experiences minimum wind. This effect is called tower shadow. A simple model for the influence of the tower is to assume a potential flow around a cylinder of a diameter 2a , see Figure (3). The radial and tangential velocity components Vθ and Vr around the tower are determined by the following equations [6]: = 1 − cos …(12) = − 1 + sin …(13) Fig. 2. Wind Shear Model [6]. Fig. 3. Tower Effect Model. 3. Blade Structural Model The rotating blade of the wind turbine is modeled as a cantilever elastic Bernoulli-Euler beam with flapwise and torsion degrees of freedom based on Stevens [8]. A pretwist angle θ is adopted in the blade model varying linearly through the span. The corresponding structural dynamic equations are derived based on Hamilton’s principle and are then discretized using the finite element with seven degree of freedom beam element. The Main coordinate systems which are chosen in the analysis of the blade model are shown in Figures (4) and (5). Fig. 4. Blade Coordinate System and Elastic Deformation [7]. Mauwafak A. Tawfik Al-Khwarizmi Engineering Journal, Vol. 9, No. 3, P.P. 15- 25 (2013) 18 Fig. 5. Cross Section Coordinate System Figure (4) shows the main coordinate system x, y, z that is fixed in the blade root with its origin in the intersection of blade root cross-section and elastic axis. When the blade is not deformed, the x axis is exactly coincident with elastic axis. The figure also shows the deformed blade and elastic displacement u, v and w in the x, y and z direction, respectively. Figure (5) shows an arbitrary blade cross section and its local coordinate system η and ζ. The torsional deflection (Ø), due to the blade deformation can also be seen. 3.1. Hamilton’s Principle The equations of motion for a rotating cantilever, isotropic blade are obtained using Hamilton’s principle [8]: ∫ [ ( − ) − ] = 0 …(14) Where U is the strain energy, T is the kinetic energy and W is the virtual work of the external forces, where their variation can be expressed as: = ∫ ∬ ( + + ) …(15) = ∫ ∬ . …(16) = ∫ ( + ) …(17) Where M is the moment about the elastic axis in the deformed coordinate system and is calculated as follows: = − …(18) Where ea is the distance from the elastic axis to the aerodynamic center. 3.2. Finite Element Method The blade is divided into a number of beam elements. Each element is a seven degree of freedom with three nodes, two external and one internal, see Figure (6) [8]. Within each element, deformations between the nodes are described using interpolating polynomials called shape functions. Fig. 6. Seven DOF Beam Element [8]. A cubic polynomial shape functions is assumed for the flap deformation (w) and a quadratic polynomial shape functions for the torsion deformation (Ø). The continuous deflections over a beam element can be expressed in term of nodal displacements: ( )∅( ) = 00 ∅ { } …(19) Where = / is the length of ith beam element, and {u}, the element nodal displacement vector is defined as: { } = , ′ , , ′ , ∅ , ∅ , ∅ …(20) The interpolating functions are defined as: = = ⎩⎪⎨ ⎪⎧ 2 − 3 + 1 ( − 2 + )−2 3 + 3 2 ( 3 − 2) ⎭⎪⎬ ⎪⎫ …(21) ∅ = ∅ ∅ ∅ = 2 − 3 + 1−4 + 4 2 2 − …(22) Mauwafak A. Tawfik Al-Khwarizmi Engineering Journal, Vol. 9, No. 3, P.P. 15- 25 (2013) 19 Application of Hamilton’s principle in the discritized form and performing the necessary calculus and assembling, the following nonlinear equation of motion for rotating blade can be obtained: [ ]{ ̈ } + [C] { ̇ }+[K]{u}={Q} …(23) Where [M], [C], and [K] are global matrices of inertia, damping and stiffness respectively, {Q} is the aerodynamic load vector. The structural damping effect was introduced into the model using the Rayleigh approach. The formula for the Rayleigh damping matrix is characterized by the following equation: [C]= a [M] + b [K] …(24) Where a and b coefficients are computed using two characteristic system eigen- frequencies ω1, ω2 and their desired damping ratio ζ1, ζ2. 4. Fully Coupled BEM-FEM Approach The system of FSI equations of motion (23) of wind turbine blade are solved by a numerical Newmark’s time integration scheme. Assuming a linear variation of the acceleration between two instants of time, the expressions for velocity and displacements are given by: ̇ ∆ = ̇ + [(1 − ) ̈ + ̈ ∆ ]∆ …(25) ∆ = + ̇ ∆ + [ − ̈ + ̈ ∆ ]∆ …(26) In which γ and β are Newmark’s parameters that are chosen based on desired stability and accuracy. The task is to evaluate the displacement, velocity and acceleration at the time t+Δt . After the displacement increment Δu is introduced to solve the FSI equations of motion, the following expressions for ̈ ∆ , ̇ ∆ , ∆ are formulated as a function of ̈ , ̇ , . ̈ ∆ = ∆ ∆ − ∆ ̇ − ( − 1) ̈ …(27) ̇ ∆ = ∆ ∆ − − 1 ̇ ∆ − 1 ̈ …(28) ∆ = + ∆ …(29) In this work, the aero elastic response of the wind turbine blade was conducted using a coupled BEM-FEM program written in MATLAB. Assuming initial values for displacement, velocity and acceleration vectors; the program first calculates the aerodynamic force vector {Q} for the blade based on BEM theory. Then at each time step, the aerodynamic forces are applied on the corresponding nodes of the blade and the program solves the structural equation based on FEM to determine the new blade deformations at time t+Δt using equations (27)-(29). Lastly, due to these new deformations the program updates the aerodynamic forces for the next time step t+Δt . The flow chart of the fully fluid – structure interaction coupling approach is shown in figure (7). 5. Numerical Results The present work is to study a single wind turbine blade’s aeroelastic response in time domain. The combined flap - torsion oscillation is presented here. A 2MW Tjaerborg wind turbine blade was chosen to perform the numerical simulation using a fully coupled FSI analysis. The detailed structural blade properties and airfoil data are well documented in [6] and [9] respectively. The blade is originally twisted and thus structural properties were modified into the main blade coordinate system [6]. Some basic parameters and geometry of this turbine are listed in Table (1). The 29 m blade is composed a series of NACA 4412-4443 profile and is discretized by 29 elements. The simulation is performed at a constant wind speed at hub Vh=12 m/s, air density is 1.225 kg/m3, wind speed exponent ν is 0.2 and the diameter of the tower is assumed to be constant with 4.5 m. The rotating speed of rotor is fairly constant speed of 20 rpm. Based on this information the blade will make one complete revolution and passes the tower every 3 seconds. The desired damping ratio ζ1=ζ2=0.02 and the time step is chosen to be Δt =0.0003. The simulation first was performed assuming the blade is rigid. Figure (8) shows the time history of angle of attack and lift coefficient of the blade at three different radial nodes. The graphs show several dips. Each dip corresponds to the instant that the blade passes the tower. Then, FSI simulation considering flexible blade has been performed. Figure (9) shows the time history of the flapwise displacement at three different radial nodes. Maximum deformation is at free end of the blade, is 2.3 m (7.93 % of blade Mauwafak A. Tawfik Al-Khwarizmi Engineering Journal, Vol. 9, No. 3, P.P. 15- 25 (2013) 20 length). This result agrees with the hypothesis introduced in Spera [10] which states that the linear formulation is adequate with small deformation where the deflection ratio is less than 10%. Figure (10) shows the time history of torsion angle at three different cross sections. The torsion angle increases with distance from the root to tip and the torsion angle undergoes high frequency oscillations. This behavior agrees qualitatively with that obtained from [1]. Figures (9) and (10) show that once the initial transient decays, there is a sudden decrease in the flapwise displacement and torsion angle at each time corresponds to the blade passes the tower. Figure (11) shows the time history of angle of attack and lift coefficient at three different nodes. The transient dynamic effect and unsmooth curves can be seen compared with rigid blade simulation result in Figure (8). It is obvious that the outer section of the blade is highly affected by the interaction between wind load and structural vibration due to its low stiffness, in comparison with the root and middle sections which are designed to have high stiffness. Elastic deformation of blade Initial Values Start Next time step Structural dynamic calculation Aerodynamic load calculation Wind load End Yes No Fig. 7. Flow Chart of the Fully FSI Coupling Algorithm. Mauwafak A. Tawfik Al-Khwarizmi Engineering Journal, Vol. 9, No. 3, P.P. 15- 25 (2013) 21 Table 1, Basic Parameters and Geometry of the 2MW Tjareborg Turbine [9]. Rotor Blades Tower Number of blades 3 Blade length 29 m Tower height 56 m Rotor diameter 61.1 m Root chord 3.3 m Shape, upper half conical Hub height 61 m Tip chord 0.9 m Diameter at base 7.25 m Rotor speed at rated power 22.36 rpm Twist (linear) 0.333 deg./m Diameter at h=28 m 4.75 m Tilt 3 deg. Blade profiles NACA 4412 - 4443 Diameter at h=56 m 4.25 m (a) Angle of Attacks. (b) Lift Coefficients. Fig. 8. Time History of (a) Angle of Attacks and (b) Lift Coefficients for Rigid Blade Computations. 1 2 3 4 5 6 7 8 9 10 11 12 4 5 6 7 8 9 10 Time (sec.) A ng le o f a tta ck (d eg .) r =27m r =21m r =12m 1 2 3 4 5 6 7 8 9 10 11 12 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 Time (sec.) Li ft co ef fic ie nt r = 27m r = 21m r = 12m Mauwafak A. Tawfik Al-Khwarizmi Engineering Journal, Vol. 9, No. 3, P.P. 15- 25 (2013) 22 Fig. 9. Time History of Flapwise Displacement at Three Different Nodes along the Blade Radial Axis. Fig. 10. Time History of Torsional Angle for Three Different Nodes along the Blade Radial Axis. (a) Angle of Attacks. 0 2 4 6 8 10 12 -0.5 0 0.5 1 1.5 2 2.5 Interaction time (sec.) F la pw is e di sp la ce m en t ( m ) r =29m r =21m r =12m 0 2 4 6 8 10 12 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 Interaction time (sec.) T or si on al a ng le ( de g. ) r = 29 m r = 21m r =12 m 1 2 3 4 5 6 7 8 9 10 11 12 4 5 6 7 8 9 10 Interaction time (sec.) A ng le o f a tta ck ( de g. ) r =27m r = 21m r = 12m Mauwafak A. Tawfik Al-Khwarizmi Engineering Journal, Vol. 9, No. 3, P.P. 15- 25 (2013) 23 (b) Lift Coefficients. Fig. 11. Time History of (a) Angle of Attacks and (b) Lift Coefficients for Flexible Blade , FSI Computations. 6. Conclusion This paper presents a fully fluid-structure interaction (FSI) modeling to solve the aeroelastic problem of a single wind turbine blade. Effects of wind shear and tower shadow are considered in this modeling. An effective and fast coupled BEM - FEM program written in MATLAB is developed for studying the aeroelastic behavior of the blade in time domain using Newmark’s time integration scheme. The main conclusions of this paper are listed below: • The results clearly show that the blade section exhibits rotation as well as flapwise deformation and stabilize after several time interactions. • The results show that the aerodynamic loads are highly affected by the flexibility of the blade which will have some effects upon the output power. • The results indicate that the blade passing the tower produces an appreciable drop in the aeroelastic response. • It is concluded that the present modal is computationally more time efficient than other CFD aeroelastic model. • The present FSI modeling is useful and can be utilized in modeling airplane propeller blades Nomenclature a radius of tower [m] B number of blades Cl lift coefficient Cd drag coefficient Cm pitching moment coefficient c length of chord [m] D drag force [N] E modulus of elasticity [N/m2] F Prandtl’s tip loss correction fg Glauert correction G modulus of rigidity [N/m2] H hub height [m], shape function L lift force [N] M moment [N.m] n normal vector; number of blades R blade length [m] r radial position [m] ⃑ blade velocity vector ⃑ undisturbed wind velocity vector ⃑ rotational velocity vector ⃑ relative velocity vector 1 2 3 4 5 6 7 8 9 10 11 12 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 Interaction time (sec.) Li ft co ef fic ie nt r = 27m r = 21m r = 12m Mauwafak A. Tawfik Al-Khwarizmi Engineering Journal, Vol. 9, No. 3, P.P. 15- 25 (2013) 24 ⃑ induced velocity vector w flap deformation [m] X distance from surface [m] α angle of attack [deg.] β Newmark’s parameter γ Newmark’s parameter strain ζ damping ratio θ pretwist blade angle [deg.] v wind profile exponent ρ density of air [Kg/m3] flow angle [deg.] , torsional deformation [deg.] 7. References [1] Y.Bazilevs, M.C.Hsu,J .Kiendl, R. Wuchner and K. U. Bletzinger,”3D Simulation of wind turbine rotors at full scale. Part II: Fluid- Structure interaction modeling with composite blades” International Journal for numerical methods in fluids, Vol.65, pp.236- 253, 2011. [2] Dong H. K. and Yoo H. K., ”Performance prediction of a 5MW wind turbine considering aeroelastic effect” World academy of science, Engineering of technology, Vol. (81), pp.771-775, 2011. [3] Shuhe W., Liaojun Z. and Yizhi Y. “Research on aeroelasticity of horizontal axis wind turbines by a fluid-structure coupling numerical method” IEEE Conference Publications, pp.1-5,2009. [4] Macolm D. J.,Laird D. L., ”Modeling of blades as equivalent beams for aeroelastic analysis,” Proceedings of AIAA/ASME. Wind Energy Symposium. Reno, USA, pp.293-303, 2003. [5] M. O .L. Hanson, J. N. Sorensen, S. Voutsinas, N. Sorensen and H .Aa. Madsen, ”State of art in wind turbine aerodynamics and aero-elasticity” Progress in aerospace sciences,Vol.42,pp.285-330,2006. [6] Martin O. L. Hanson, ” Aerodynamics of wind turbines,” Second edition, Earthscan, UK, 2008. [7] Hodges D. H. and Dowel E.H., ”Nonlinear equation of motion for elastic bending and torsion of twisted non-uniform blades” NASA TND-7818, 1974. [8] Patricia L. S. ,”Active interrogation of helicopter main rotor faults using trailing edge flap actuation” Ph.D. thesis, The Pennsylvania State University,2001. [9] Tjaerborg wind turbine ‘Esbjerg’, “Geometric and operational data http://130.226.17.201/extra/web_docs/tjare/V K-184-901130.pdf. [10] Spera D.A.,”Wind turbine technology: fundamental concepts of wind turbine engineering”, 2nd edition, ASME press, New York, 1995. http://130.226.17.201/extra/web_docs/tjare/V )2013( 15- 25 ، صفحة3، العدد9مجلة الخوارزمي الھندسیة المجلد موفق علي توفیق 25 بین المائع التام تحلیل سلوك المرونھ الھوائیھ لشفرة توربین ریاح بأستخدام نموذج للتداخل والھیكل **فاروق عمر حمدون **محمد أدریس أبو طبیخ * موفق علي توفیق ةالجامعھ التكنولوجی/ قسم ھندسة المكائن والمعدات ** ،* جامعة واسط/ ة كلیة الھندس/ قسم الھندسھ المیكانیكیة ** الخالصة .ئع والھیكلیتناول البحث الحالي تحلیال لسلوك المرونھ الھوائیھ لشفرة توربین ریاح من خالل بناء نموذج عددي یعتمد مبدأ التداخل التام بین الما كما تم تمثیل . الشفره في حساب األحمال األیرودینامیھ المؤثره على شفرة توربین الریاح تحت تأثیر القص ووجود البرج أستخدمت طریقة زخم عنصر تمت دراسة أستجابھ المرونھ الھوائیھ للشفره من خالل بناء . في حساب االھتزازات الناتجھ الشفره بعتبھ كابولیھ دواره بأستخدام طریقة العناصر المحدده كل بأستخدام الماتالب یربط بین النموذج الریاضي للقوى األیرودینامیھ ونموذج أھتزاز الھیكل من خالل تبادل البیانات بین النموذجین عند برنامج حاسوبي یة أستخدام النموذج أظھرت النتائج المستحصلھ أمكان. لمعادالت الحركھ للتكامل الزمني) Newmark(خطوه زمنیھ بصوره تكراریھ و بأستخدام طریقة األستفاده من النموذج الحالي في تحلیل المرونھ كما یمكن. المقترح في أجراء تقییم سریع لشفرات توربین الریاح تحت تأثیر التداخل بین المائع والھیكل . شفرات الطائرات المروحیھالھوائیھ ل