WISDEM / CCBlade

A blade element momentum method for analyzing wind turbine aerodynamic performance that is robust (guaranteed convergence), fast (superlinear convergence rate), and smooth (continuously differentiable).
Other
42 stars 38 forks source link

About the 5MW test case #33

Open zsacfder opened 2 months ago

zsacfder commented 2 months ago

Hello,I'm having some problems running the 5MW test code. Why are the power curve and thrust curve at high wind speeds lower than the data in the file provided by NREL? The error seems a bit large at high wind speeds?

zsacfder commented 2 months ago

111222

gbarter commented 2 months ago

Please post or attach the code you are using to generate this comparison.

zsacfder commented 2 months ago

Thank you very much for your reply. I used the code from the help file for the test. Here is the code I used:

class TestNREL5MW(unittest.TestCase):
    def setUp(self):
        # 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

        # atmosphere
        rho = 1.225
        mu = 1.81206e-5

        afinit = CCAirfoil.initFromAerodynFile  # just for shorthand
        basepath = "C:/Users/Admin/Desktop/all/blade/uncertainty_anasys/5MW_AFFiles"

        # load all airfoils
        airfoil_types = [0] * 8
        airfoil_types[0] = afinit(path.join(basepath, "Cylinder1.dat"))
        airfoil_types[1] = afinit(path.join(basepath, "Cylinder2.dat"))
        airfoil_types[2] = afinit(path.join(basepath, "DU40_A17.dat"))
        airfoil_types[3] = afinit(path.join(basepath, "DU35_A17.dat"))
        airfoil_types[4] = afinit(path.join(basepath, "DU30_A17.dat"))
        airfoil_types[5] = afinit(path.join(basepath, "DU25_A17.dat"))
        airfoil_types[6] = afinit(path.join(basepath, "DU21_A17.dat"))
        airfoil_types[7] = afinit(path.join(basepath, "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]]

        tilt = 0
        precone = 2.5
        yaw = 0.0

        # create CCBlade object
        self.rotor = CCBlade(r, chord, theta, af, Rhub, Rtip, B, rho, mu, precone, tilt, yaw, shearExp=0.2, hubHt=90.0,derivatives=True)

    def test_thrust_torque(self):

        Uinf = np.array([3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25])
        Omega = np.array(
            [
                6.972,
                7.183,
                7.506,
                7.942,
                8.469,
                9.156,
                10.296,
                11.431,
                11.890,
                12.100,
                12.100,
                12.100,
                12.100,
                12.100,
                12.100,
                12.100,
                12.100,
                12.100,
                12.100,
                12.100,
                12.100,
                12.100,
                12.100,
            ]
        )
        pitch = np.array(
            [
                0.000,
                0.000,
                0.000,
                0.000,
                0.000,
                0.000,
                0.000,
                0.000,
                0.000,
                3.823,
                6.602,
                8.668,
                10.450,
                12.055,
                13.536,
                14.920,
                16.226,
                17.473,
                18.699,
                19.941,
                21.177,
                22.347,
                23.469,
            ]
        )

        Pref = np.array(
            [
                42.9,
                188.2,
                427.9,
                781.3,
                1257.6,
                1876.2,
                2668.0,
                3653.0,
                4833.2,
                5296.6,
                5296.6,
                5296.6,
                5296.6,
                5296.6,
                5296.6,
                5296.6,
                5296.6,
                5296.7,
                5296.6,
                5296.7,
                5296.6,
                5296.6,
                5296.7,
            ]
        )
        Tref = np.array(
            [
                171.7,
                215.9,
                268.9,
                330.3,
                398.6,
                478.0,
                579.2,
                691.5,
                790.6,
                690.0,
                608.4,
                557.9,
                520.5,
                491.2,
                467.7,
                448.4,
                432.3,
                418.8,
                406.7,
                395.3,
                385.1,
                376.7,
                369.3,
            ]
        )
        Qref = np.array(
            [
                58.8,
                250.2,
                544.3,
                939.5,
                1418.1,
                1956.9,
                2474.5,
                3051.1,
                3881.3,
                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,
            ]
        )

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

        outputs, derivs = self.rotor.evaluate(Uinf, Omega, pitch,coefficients=True)
        P, T, Q = [outputs[key] for key in ("P", "T", "Q")]

        import matplotlib.pyplot as plt
        plt.plot(Uinf, P/1e6)
        plt.plot(Uinf, Pref/1e3)
        plt.figure()
        plt.plot(Uinf, T/1e6)
        plt.plot(Uinf, Tref/1e3)
        plt.show()

        idx = Uinf < 15
        np.testing.assert_allclose(Q[idx] / 1e6, Qref[idx] / 1e3, atol=0.15)
        np.testing.assert_allclose(P[idx] / 1e6, Pref[idx] / 1e3, atol=0.2)  # within 0.2 of 1MW
        np.testing.assert_allclose(T[idx] / 1e6, Tref[idx] / 1e3, atol=0.15)

def suite():
    suite = unittest.TestSuite()
    suite.addTest(unittest.makeSuite(TestNREL5MW))
    return suite

if __name__ == "__main__":
    result = unittest.TextTestRunner().run(suite())

    if result.wasSuccessful():
        exit(0)
    else:
        exit(1)
gbarter commented 2 months ago

Thank you. That test script for CCBlade is passing in hard-coded blade pitch angles (see the line outputs, derivs = self.rotor.evaluate(Uinf, Omega, pitch,coefficients=True)), which aren't necessarily the right pitch angles for true constant-torque performance in Region III of the powercurve where you are seeing the mismatch. This is why the unit test only checks for the near-rated and below wind speeds for comparison (where the pitch setting is very low or zero).

zsacfder commented 2 months ago

Thanks for the guidance, I seem to understand what is causing the problem, I have adjusted the pitch angle for high wind speeds so that the power is stable at 5 MW. power seems to be very sensitive at high wind speeds and large pitch angles. I can't thank you enough for your help.