ccblade & NREL 5 MW

Wind & Water Systems Engineering Software Tools
1

Dear ccblade users/developers,

For my own education, I have modified the “example.py” file that comes with ccblade docs folder to try and compare the results of ccblade against (what I believe they are) official power curve, Cp, thrust and torque for the NREL 5MW reference turbine.

The data for the power curve/power coefficient for the NREL 5MW reference turbine come from here: github.com/NREL/turbine-models/ … 26_RWT.csv

For thrust, torque, rotor speed and pitch, I have used the data reported here: NREL5MW Rotor Speed Data

Now, by naively running the almost-unmodified code, I was kind of expecting to see the BEM code almost reproduce the reported numbers:

Image link: postimg.cc/Whb340bb

(Seems like I am unable to link images directly in the post for some reason…)

While Cp and the first parts for the power curve and torque look very good, I am a bit puzzled by the behaviour of the thrust curve and the end part of the power/torque curve. I am sure I am misunderstanding something big or I am missing something even bigger here :slight_smile: .

I am pasting below my Python script to reproduce my graphs - the NREL data is embedded in it.

I’ll appreciate any suggestion/comment you may have. Thank you in advance.

Andrea.


import numpy as np
from math import pi
import matplotlib.pyplot as plt

from ccblade import CCAirfoil, CCBlade


def Main():

    # geometry
    Rhub = 1.5
    Rtip = 63.0

    r = np.array([2.8667, 5.6000, 8.3333, 11.7500, 15.8500, 19.9500, 24.0500,
                  28.1500, 32.2500, 36.3500, 40.4500, 44.5500, 48.6500, 52.7500,
                  56.1667, 58.9000, 61.6333])
    chord = np.array([3.542, 3.854, 4.167, 4.557, 4.652, 4.458, 4.249, 4.007, 3.748,
                      3.502, 3.256, 3.010, 2.764, 2.518, 2.313, 2.086, 1.419])
    theta = np.array([13.308, 13.308, 13.308, 13.308, 11.480, 10.162, 9.011, 7.795,
                      6.544, 5.361, 4.188, 3.125, 2.319, 1.526, 0.863, 0.370, 0.106])
    B = 3  # number of blades

    tilt = 5.0
    precone = 2.5
    yaw = 0.0

    nSector = 18  # azimuthal discretization

    # atmosphere
    rho = 1.225
    mu = 1.81206e-5

    # power-law wind shear profile
    shearExp = 0.2
    hubHt = 90.0

    # turbine characteristics
    cut_out = 25.0   # m/s

    # 2 ----------

    afinit = CCAirfoil.initFromAerodynFile  # just for shorthand

    # load all airfoils
    airfoil_types = [0]*8
    airfoil_types[0] = afinit('data/Cylinder1.dat')
    airfoil_types[1] = afinit('data/Cylinder2.dat')
    airfoil_types[2] = afinit('data/DU40_A17.dat')
    airfoil_types[3] = afinit('data/DU35_A17.dat')
    airfoil_types[4] = afinit('data/DU30_A17.dat')
    airfoil_types[5] = afinit('data/DU25_A17.dat')
    airfoil_types[6] = afinit('data/DU21_A17.dat')
    airfoil_types[7] = afinit('data/NACA64_A17.dat')

    # place at appropriate radial stations
    af_idx = [0, 0, 1, 2, 3, 3, 4, 5, 5, 6, 6, 7, 7, 7, 7, 7, 7]

    af = [0]*len(r)
    for i in range(len(r)):
        af[i] = airfoil_types[af_idx[i]]

    # 2 ----------

    # 3 ----------

    # create CCBlade object
    rotor = CCBlade(r, chord, theta, af, Rhub, Rtip, B, rho, mu,
                    precone, tilt, yaw, shearExp, hubHt, nSector)

    # 4 ----------

    # Data for NREL 5 MW
    Omega = np.array([0, 0, 0, 0, 0, 6.972, 7.0775, 7.183, 7.3445, 7.506, 7.724, 7.942, 8.2055, 8.469, 8.8125,
                      9.156, 9.726, 10.296, 10.8635, 11.431, 11.6605, 11.89, 11.995, 12.1, 12.1, 12.1, 12.1,
                      12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1,
                      12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1])

    pitch = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.9115, 3.823,
                      5.2125, 6.602, 7.635, 8.668, 9.559, 10.45, 11.2525, 12.055, 12.7955, 13.536,
                      14.228, 14.92, 15.573, 16.226, 16.8495, 17.473, 18.086, 18.699, 19.32, 19.941,
                      20.559, 21.177, 21.762, 22.347, 22.908, 23.469])

    # Results for NREL 5 MW
    # Taken from http://forums.nrel.gov/t/nrel5mw-rotor-speed-data/867/1
    ROTOR_POWER = np.array([0, 0, 0, 0, 0, 42.9, 115.55, 188.2, 308.05, 427.9, 604.6, 781.3, 1019.45, 1257.6,
                            1566.9, 1876.2, 2272.1, 2668, 3160.5, 3653, 4243.1, 4833.2, 5064.9, 5296.6, 5296.6,
                            5296.6, 5296.6, 5296.6, 5296.6, 5296.6, 5296.6, 5296.6, 5296.6, 5296.6, 5296.6,
                            5296.6, 5296.6, 5296.6, 5296.65, 5296.7, 5296.65, 5296.6, 5296.65, 5296.7, 5296.65,
                            5296.6, 5296.6, 5296.6, 5296.65, 5296.7])/1e3

    POWER_CURVE = np.array([0, 0, 0, 0, 0, 40.5, 109.1, 177.7, 290.8, 403.9, 570.75, 737.6, 962.4, 1187.2,
                            1479.15, 1771.1, 2144.85, 2518.6, 2983.5, 3448.4, 4005.45, 4562.5, 4781.25, 5000, 5000,
                            5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000,
                            5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000])/1e3

    THRUST = np.array([0, 0, 0, 0, 0, 171.7, 193.8, 215.9, 242.4, 268.9, 299.6, 330.3, 364.45, 398.6, 438.3, 478,
                       528.6, 579.2, 635.35, 691.5, 741.05, 790.6, 740.3, 690, 649.2, 608.4, 583.15, 557.9, 539.2,
                       520.5, 505.85, 491.2, 479.45, 467.7, 458.05, 448.4, 440.35, 432.3, 425.55, 418.8, 412.75,
                       406.7, 401, 395.3, 390.2, 385.1, 380.9, 376.7, 373, 369.3])

    TORQUE = np.array([0, 0, 0, 0, 0, 58.8, 154.5, 250.2, 397.25, 544.3, 741.9, 939.5, 1178.8, 1418.1, 1687.5, 1956.9,
                       2215.7, 2474.5, 2762.8, 3051.1, 3466.2, 3881.3, 4030.7, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1,
                       4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1,
                       4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1, 4180.1])

    CP = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.196405748, 0.27998048600000003, 0.363555224, 0.3933196965, 0.423084169,
                    0.435105171, 0.447126173, 0.4501645485, 0.453202924, 0.45306905350000004, 0.452935183, 0.4526530135,
                    0.452370844, 0.4519479375, 0.451525031, 0.450679475, 0.44883706, 0.430468281, 0.378869113,
                    0.3384299655, 0.297990818, 0.26828947400000003, 0.23858813, 0.21628455800000002, 0.193980986,
                    0.1769081965, 0.159835407, 0.146545612, 0.133255817, 0.122756666, 0.112257515, 0.1038533405,
                    0.095449166, 0.088642447, 0.081835728, 0.07626425349999999, 0.070692779, 0.06608858649999999,
                    0.061484394, 0.057646356999999995, 0.05380832, 0.0505834795, 0.047358639, 0.044629266, 0.041899893])


    wind_speeds = Uinf = np.arange(0.5, 25.1, 0.5)

    outputs, derivs = rotor.evaluate(Uinf, Omega, pitch)
    P = outputs['P']
    T = outputs['T']
    Q = outputs['Q']

    outputs, derivs = rotor.evaluate(Uinf, Omega, pitch, coefficients=True)
    cp = outputs['CP']

    cp[wind_speeds > cut_out] = 0.0

    fig = plt.figure()
    ax = fig.add_subplot(2, 2, 1)

    ax.plot(wind_speeds, CP, lw=2, label='NREL')
    ax.plot(wind_speeds, cp, lw=2, label='ccblade')
    ax.legend(loc=0, fancybox=True, shadow=True, fontsize=14)

    ax.grid()
    ax.set_xlabel('Wind Speed $(m/s)$', fontsize=16)
    ax.set_ylabel('Power Coefficient $(-)$', fontsize=16)

    ax = fig.add_subplot(2, 2, 2)
    idx = np.where(wind_speeds >= 25.0)[0]
    P[P<0] = 0.0
    P2 = P.copy()
    P2[P2>5e6] = 5e6
    P2[idx] = 0.0

    ax.plot(wind_speeds, POWER_CURVE, lw=2, label='NREL')
    ax.plot(wind_speeds, P2/1e6, lw=2, label='ccblade')

    ax.grid()
    ax.set_xlabel('Wind Speed $(m/s)$', fontsize=16)
    ax.set_ylabel('Power $(MW)$', fontsize=16)
    ax.legend(loc=0, fancybox=True, shadow=True, fontsize=14)
    ax.set_ylim(0, 6)


    ax = fig.add_subplot(2, 2, 3)
    ax.plot(wind_speeds, THRUST, lw=2, label='NREL')
    ax.plot(wind_speeds, T/1e3, lw=2, label='ccblade')

    ax.grid()
    ax.set_xlabel('Wind Speed $(m/s)$', fontsize=16)
    ax.set_ylabel('Thrust $(kN)$', fontsize=16)
    ax.legend(loc=0, fancybox=True, shadow=True, fontsize=14)

    Q[Q<0] = 0.0

    ax = fig.add_subplot(2, 2, 4)
    ax.plot(wind_speeds, TORQUE, lw=2, label='NREL')
    ax.plot(wind_speeds, Q/1e3, lw=2, label='ccblade')

    ax.grid()
    ax.set_xlabel('Wind Speed $(m/s)$', fontsize=16)
    ax.set_ylabel('Torque $(kNm)$', fontsize=16)
    ax.legend(loc=0, fancybox=True, shadow=True, fontsize=14)

    for ax in fig.axes:
        ax.set_xlim(0, 30)

    plt.show()


if __name__ == '__main__':
    Main()
2

Dear Andrea,

Just a couple comments:

  • The results you label as “NREL” where calculated by FAST, which includes terms not accounted for in ccblade, like weight, blade deflection, and tower deflection;
  • Regarding the offset in rotor thrust, the main reason for the offset of rotor thrust is that FAST includes the weight of the rotor projected along the tilted shaft) in its calculation of RotThrust, not just the pure aerodynamic applied thrust output by ccblade.

Best regards,

3

Dear Jason,

Thank you for your answer, I really appreciate it. Regarding your point number (2), I had indeed found today (after quite some search) that it is exactly as you wrote, so I had modified the thrust array you posted in this way:

m_rotor = 110.0  # kg
g = 9.81
tilt = 5 * pi / 180.0
thrust = thrust - m_rotor * g * sin(tilt)  # remove weight of rotor that is included in reported results

And the match on thrust looks much better indeed:

1474×930 92.7 KB

What I don’t honestly understand is the behavior at the tail end, and especially for the power curve.

Why would ccblade tell me that there is no way to produce at rated power after 22 m/s? I understand if differences (even large ones) are found in the ramp up and/or around rated speed, but the only way that I can interpret the power curve dropping off so badly at the end is that the Cp is way too low at that point for a BEM model to actually represent “reality” (whatever that means, of course).

But it also mean that, at least in this case, there is no way to model a power curve from cut-in to cut-out beside changing some of the input assumptions - which in a way invalidates the potential comparison between a BEM model and more sophisticated tools like FAST).

Andrea.

4

Dear Andrea,

Indeed that change to thrust seems to fix the main problem. Please note that in FAST, the tower deflection may also influence the tilt a bit.

Regarding the difference at high wind speeds, I’m not sure. Is the result sensitive to the blade-pitch angle at those wind speeds? You are using the optimal pitch angles from FAST, but the aerodynamic modeling is a bit different between CCBlade and FAST and the FAST solution also accounts for blade deflection, not accounted for in CCBlade. So, it is possible that CCBlade would require a bit different pitch angles to achieve rated power at high wind speeds than OpenFAST.

Best regards,