Modeling of the NREL 5-MW wind turbine blades

Hi,
I am modeling the NREL 5-MW wind turbine in Abaqus/CAE. The blades are generated using the properties of different sections of the airfoils that were considered the same as ref. [1]. Table 2-1 in ref. [1] lists the distributed blade structural properties of the NREL offshore 5-MW baseline wind turbine. The stiffness and orientation of the principal elastic axes can vary for different cross-sections along the length of a blade, depending on the blade geometry.
For the blade properties, 49 stations along the span of the blades were considered. A generalized section is used for different elements of each blade. The properties of this generalized section consisted of the area (A), the moment of inertia along two axes (I11, I22), and polar moment of inertia (J). The material used for blades has an elastic modulus of 13.1 GPa, and shear modulus of 8 GPa according to ref. [2]. The mass of blades is given as longitudinal mass in line with the blade elements.
The strong axis of the blades is oriented at the edgewise direction of the blades. The flapwise and edgewise section stiffness and inertia values, are given about the principal structural axes of each cross-section as oriented by the structural-twist angle (StrcTwst) [1]. Each of the 49 elements of each blade is given a change of angle based on the value of StrcTwst. In fact, elements twist along the blade span from 13.308° at the blade root to 0° at the blade tip.
The Abaqus program handling the structural simulation requires following parameters for a generalized section:
A, I11, I12, I22, J.
To my understanding, the required values from the reference are calculated as follows.
A = EAStff / E
I11 = FlpStff / E
I12 = Null*
I22 = EdgStff / E
J = GJStff / G
Density = BMassDen / Area

  • According to “Moments of Inertia” rules, the value of I12 must satisfy the inequality –(I11+I22)/2 < I12 <= (I11+I22)/2.
    where A is area (m^2), E is elastic modulus ¶, I11 is the area moment of inertia for the flapwise (m^4), I12 is the area moment of inertia with respect to flapwise and edgewise (m^4), I22 is the area moment of inertia for the edgewise (m^4), J is polar moment of inertia (m^4), Density is the density of each element (kg/m^3).
    It seems that the definition of two equations, FlpStff and EdgStff, has been mistaken in ref. [3]. So both formulas have the same x^2 expressions.
    FlpStff = ∫∫ E(x, y) x^2 dxdy
    EdgStff = ∫∫ E(x, y) x^2 dxdy
    where E(x,y) is the modulus of elasticity in N/m^2, and x and y are the flapwise and edgewise distances in meters from the blade section elastic center to the differential area element, respectively.
    Nevertheless, it seems that the FlpStff is the bending stiffness about the principal elastic x-axis, and the EdgStff is the bending stiffness about the principal elastic y-axis. According to this definition, values FlpStff and EdgStff are calculated as follows:
    FlpStff = ∫∫ E(x, y) y^2 dxdy
    EdgStff = ∫∫ E(x, y) x^2 dxdy
    where x and y are the edgewise and flapwise distances in meters from the blade section elastic center to the differential area element, respectively.

Is my interpretation of wind turbine blade modeling correct? If my interpretation is correct, do the parameters displayed in Table 2-1 in ref. [1] need to be modified?
Regards,

Dear Mohsen,

You haven’t stated what Ref. [1], [2], or [3] are. I’m assuming [1] is the NREL 5-MW baseline wind turbine specifications report: nwtc.nrel.gov/system/files/FAST.pdf. I’m not sure what [2] is, so, I can’t confirm E or G.

I agree that there is a typo in [3], but I disagree with your interpretation. The correct equations are:
FlpStff = ∫∫ E(x, y) x^2 dxdy
EdgStff = ∫∫ E(x, y) y^2 dxdy
where x and y are the flapwise and edgewise distances in meters from the blade section elastic center to the differential area element, respectively.

The values in Table 2-1 are correct.

Best regards,

Dear Jason,
First, I apologize for forgetting to write the references. Ref. [2] was written by Asareh et al. Also, Asareh and Prowell [4] have written “Seismic Loading for FAST” Report. As you know, this approach allows consideration of the seismic load source within the FAST code.
In the following, the relationships I wrote are based on the assumption of Fig. 1 that represents the wind turbine blade components. As can be seen, the axis x is assumed to be horizontal, and the axis y is vertical (by default). According to this figure, the second moments area of flapwise and edgewise are used around axes x and y, respectively. If we take these assumptions, according to chapter 10 of ref. [5] and ref. [6], values FlpStff and EdgStff are calculated as follows:
FlpStff = ∫∫ E(x, y) y^2 dxdy
EdgStff = ∫∫ E(x, y) x^2 dxdy
Aren’t your hypothetical axes in line with Fig. 1?
Sincerely,


Fig. 1: Wind turbine blade components

References
[1] Jonkman J, Butterfield S, Musial W, Scott G. Definition of a 5-MW reference wind turbine for offshore system development. National Renewable Energy Lab.(NREL), Golden, CO (United States); 2009 Feb 1 (available from here: nrel.gov/docs/fy09osti/38060.pdf).
[2] Asareh MA, Schonberg W, Volz J. Fragility analysis of a 5-MW NREL wind turbine considering aero-elastic and seismic interaction using finite element method. Finite Elements in Analysis and Design. 2016 Nov 1;120:57-67. doi:10.1016/j.finel.2016.06.006.
[3] Jonkman JM, Buhl Jr ML. FAST user’s guide. National Renewable Energy Laboratory, Golden, CO, Technical Report No. NREL/EL-500-38230. 2005 Aug (available from here: nwtc.nrel.gov/system/files/FAST.pdf).
[4] Asareh MA, Prowell I. Seismic loading for FAST. No. NREL/SR-5000-53872, National Renewable Energy Laboratory, Golden, CO. 2011 Aug (available from here: nrel.gov/docs/fy12osti/53872.pdf).
[5] Hibbeler RC. Engineering mechanics: statics. 12th ed, 2010.
[6] en.wikipedia.org/wiki/Second_moment_of_area.

Dear Moshen.Motlagh,

The x and y axis are not shown correctly in your Figure (that you sent me via e-mail). Instead, in ElastoDyn, the local x-axis of the blade coordinate system is directed towards the suction side of the airfoil and y is directed towards the trailing edge. This is seen in the old FAST v7 User’s Guide: nwtc.nrel.gov/system/files/FAST.pdf. I confirm my original expressions for FlpStff and EdgStff.

Best regards,

Dear Jason,
Thank you for the clarification regarding the blade coordinate system. It has resolved my confusion.
Next, I have a question about the modeling issue I raised in the first post. As you can see, I considered E and G according to the ref. [2] to constant values and calculated the rest of the parameters for each of the 49 pieces of blades.
Could you please confirm this method? Do you recommend a better method for FE modeling of wind turbine blades?
Sincerely,

Dear Mohsen,

I don’t think I can comment on that. I would say that as long as the beam properties are equivalent, than you can partition out E, G, etc. as you want. In general, your equations assume that the beam sections are isotropic i.e., E and G are constant throughout the blade. However, this is likely not the case in the real blade.

Best regards,

Dear Jason,
Following the modeling of the NREL 5-MW wind turbine in Abaqus software and finding the natural frequencies and mode shapes of it, I considered its lumped mass in the Nacelle (CM)_x, Nacelle (CM)_y position. The coordinate axis at the center of the base of the turbine tower is assumed. The x component along the wind and the y component towards the top of the tower were considered positive. The following values were used according to ref. [1],
Nacelle (CM)_x = NacCMxn = 1.9 m
Nacelle (CM)_y = TowerHt + NacCMzn = 87.6 + 1.75 = 89.35 m
According to ref. [2], the nacelle inertia about the yaw axis was taken to be 2,607,890 kg.m^2. It is equivalent to the DOWEC turbine’s nacelle inertia about its nacelle CM, but translated to the yaw axis using the parallel-axis theorem with the nacelle mass and downwind distance to the nacelle CM.
Fig. 16 of ref. [1] shows the NcIMUxn and NcIMUzn parameters. Using these parameters, Nacelle (IMU)_x and Nacelle (IMU)_y values were obtained as follows, respectively. NcIMUxn and NcIMUzn parameters are obtained from Appendix A (FAST input files) of ref. [2].
Nacelle (IMU)_x = NcIMUxn = 3.09528 m
Nacelle (IMU)_y = TowerHt + NcIMUzn = 87.6 + 2.23336 = 89.83336 m
According to the above description, is it correct to consider the coordinates of Nacelle (IMU)_x and Nacelle (IMU)_y for nacelle inertia? Or just put it in the center of the yaw axis and at the height of Nacelle (IMU)_y?
Sincerely,

References
[1] Jonkman JM, Buhl Jr ML. FAST user’s guide. National Renewable Energy Laboratory, Golden, CO, Technical Report No. NREL/EL-500-38230. 2005 Aug (available from here: nwtc.nrel.gov/system/files/FAST.pdf).
[2] Jonkman J, Butterfield S, Musial W, Scott G. Definition of a 5-MW reference wind turbine for offshore system development. National Renewable Energy Lab.(NREL), Golden, CO (United States); 2009 Feb 1 (available from here: nrel.gov/docs/fy09osti/38060.pdf).

Dear Mohsen,

Actually NcIMUxn = -3.09528 m (note the negative sign).

I’m not sure I understand your question, but just to clarify, the nacelle IMU of the NREL 5-MW baseline turbine is located at the main bearing, not the nacelle CM.

Best regards,

Dear Jason,
Thank you very much for your reply. Based on your tips, I made the wind turbine FEM. I obtained the full-system natural frequencies by Abaqus/CAE. I compared them with both the FAST model and the ADAMS model, according to Table 9-1 of ref. [1]. The agreement between the results of the FEM with FAST and ADAMS results was very good. However, my model reports the 1.0236 Hz for 1st Drivetrain Torsion frequency, which is very different from the values of the two models mentioned. Do you propose a specific solution for turbine modeling so that this frequency is also well adapted to the NREL 5-MW baseline wind turbine?
Are the values of the structural twist “StrcTwst” in the values FlpStff, EdgStff, and GJStff of Table 2-1 ref. [1] considered? Or do I have to rotate all (49) parts of blades in the FEM?
Best regards,

References
[1] Jonkman J, Butterfield S, Musial W, Scott G. Definition of a 5-MW reference wind turbine for offshore system development. National Renewable Energy Lab.(NREL), Golden, CO (United States); 2009 Feb 1 (available from here: nrel.gov/docs/fy09osti/38060.pdf).

Dear Mohsen,

I don’t know why your drivetrain torsion frequency differs, but certainly this frequency depends on the boundary conditions applied. Table 9-1 from [1] reports a fixed-free boundary condition because the generator is locked by a brake. A free-free boundary condition will result in a higher natural frequency (but I would expect higher than 1.02 Hz).

I’m not sure I understand your question regarding the structural twist (StrcTwst) and the blade stiffness values, but the blade stiffness values are about the principle axes of bending, which are rotated by the structural twist.

Best regards,

Dear Jason,
Thank you so much. In the next step, I intend to simulate the rotation of the wind turbine blades at an angular velocity between Cut-In and Rated Rotor Speed for one minute. At this point, the effect of the wind on the blades and towers is ignored. Actually, the impacts of the rotating blade and the shadow effect are examined.
For this purpose, the time history of the forces and moments or accelerations created by the rotor rotation (at regular intervals) is required. Please send me if such a record is available. Otherwise, please explain how I can create such a document by the FAST program for the land-based 5 MW baseline wind turbine.
Best regards,

Dear Mohsen,

I’m not sure I understand your question. It sounds like you want to ignore aerodynamic loads, but also include the shadow effect, which does not sound consistent to me.

Regardless, to disable aerodynamic loads, set CompAero = False in the FAST primary input file. In this case I would recommend disabling the controller as well by setting CompServo = False. You can set the rotor speed (RotSpeed) in the ElastoDyn input file .

Best regards,

Dear Jason,

I am working on implementation of Trailing edge flaps(TEF) for the 5MW turbine. I have developed the TEF mechanism and I need to scale the blade mass values on addition of the TEF. I looked into the NREL technical report and found the blade mass values reported in kg/m . Can you please tell me whether it is the mass per unit chord at the point or it is some other form.

Best Regards

Dear Kenneth,

The distributed blade mass is reported as the mass per unit length along the blade. This is equivalent to the integral of the material density over the cross-sectional area, which is unique for each radial location.

Best regards,

Thank you for the clarification Jason