Python CCBlade.evaluate方法代碼示例

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
    def test_dprecurveTip2(self):

        precurve = np.linspace(1, 10, self.n)
        precurveTip = 10.1
        precone = 0.0
        rotor = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip,
            self.B, self.rho, self.mu, precone, self.tilt, self.yaw, self.shearExp,
            self.hubHt, self.nSector, derivatives=True, precurve=precurve, precurveTip=precurveTip)

        P, T, Q, dP_ds, dT_ds, dQ_ds, dP_dv, dT_dv, dQ_dv = \
            rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=False)

        dT_dprecurveTip = dT_ds[0, 5]
        dQ_dprecurveTip = dQ_ds[0, 5]
        dP_dprecurveTip = dP_ds[0, 5]

        pct = float(precurveTip)
        delta = 1e-6*pct
        pct += delta

        rotor = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip,
            self.B, self.rho, self.mu, precone, self.tilt, self.yaw, self.shearExp,
            self.hubHt, self.nSector, derivatives=False, precurve=precurve, precurveTip=pct)

        Pd, Td, Qd = rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=False)

        dT_dprecurveTip_fd = (Td - T) / delta
        dQ_dprecurveTip_fd = (Qd - Q) / delta
        dP_dprecurveTip_fd = (Pd - P) / delta

        np.testing.assert_allclose(dT_dprecurveTip_fd, dT_dprecurveTip, rtol=1e-4, atol=1e-8)
        np.testing.assert_allclose(dQ_dprecurveTip_fd, dQ_dprecurveTip, rtol=1e-4, atol=1e-8)
        np.testing.assert_allclose(dP_dprecurveTip_fd, dP_dprecurveTip, rtol=1e-4, atol=1e-8)

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
    def test_dpresweepTip3(self):

        presweep = np.linspace(1, 10, self.n)
        presweepTip = 10.1
        precone = 0.0
        rotor = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip,
            self.B, self.rho, self.mu, precone, self.tilt, self.yaw, self.shearExp,
            self.hubHt, self.nSector, derivatives=True, presweep=presweep, presweepTip=presweepTip)

        CP, CT, CQ, dCP_ds, dCT_ds, dCQ_ds, dCP_dv, dCT_dv, dCQ_dv = \
            rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=True)

        dCT_dpresweepTip = dCT_ds[0, 6]
        dCQ_dpresweepTip = dCQ_ds[0, 6]
        dCP_dpresweepTip = dCP_ds[0, 6]

        pst = float(presweepTip)
        delta = 1e-6*pst
        pst += delta

        rotor = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip,
            self.B, self.rho, self.mu, precone, self.tilt, self.yaw, self.shearExp,
            self.hubHt, self.nSector, derivatives=False, presweep=presweep, presweepTip=pst)

        CPd, CTd, CQd = rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=True)

        dCT_dpresweepTip_fd = (CTd - CT) / delta
        dCQ_dpresweepTip_fd = (CQd - CQ) / delta
        dCP_dpresweepTip_fd = (CPd - CP) / delta

        np.testing.assert_allclose(dCT_dpresweepTip_fd, dCT_dpresweepTip, rtol=1e-4, atol=1e-8)
        np.testing.assert_allclose(dCQ_dpresweepTip_fd, dCQ_dpresweepTip, rtol=1e-4, atol=1e-8)
        np.testing.assert_allclose(dCP_dpresweepTip_fd, dCP_dpresweepTip, rtol=1e-4, atol=1e-8)

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
    def test_dtheta3(self):

        dCT_dtheta = self.dCT_dv[0, 2, :]
        dCQ_dtheta = self.dCQ_dv[0, 2, :]
        dCP_dtheta = self.dCP_dv[0, 2, :]
        dCT_dtheta_fd = np.zeros(self.n)
        dCQ_dtheta_fd = np.zeros(self.n)
        dCP_dtheta_fd = np.zeros(self.n)

        for i in range(self.n):
            theta = np.array(self.theta)
            delta = 1e-6*theta[i]
            theta[i] += delta

            rotor = CCBlade(self.r, self.chord, theta, self.af, self.Rhub, self.Rtip,
                self.B, self.rho, self.mu, self.precone, self.tilt, self.yaw, self.shearExp,
                self.hubHt, self.nSector, derivatives=False)

            CPd, CTd, CQd = rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=True)

            dCT_dtheta_fd[i] = (CTd - self.CT) / delta
            dCQ_dtheta_fd[i] = (CQd - self.CQ) / delta
            dCP_dtheta_fd[i] = (CPd - self.CP) / delta

        np.testing.assert_allclose(dCT_dtheta_fd, dCT_dtheta, rtol=5e-6, atol=1e-8)
        np.testing.assert_allclose(dCQ_dtheta_fd, dCQ_dtheta, rtol=7e-5, atol=1e-8)
        np.testing.assert_allclose(dCP_dtheta_fd, dCP_dtheta, rtol=7e-5, atol=1e-8)

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
    def test_dtheta2(self):

        dT_dtheta = self.dT_dv[0, 2, :]
        dQ_dtheta = self.dQ_dv[0, 2, :]
        dP_dtheta = self.dP_dv[0, 2, :]
        dT_dtheta_fd = np.zeros(self.n)
        dQ_dtheta_fd = np.zeros(self.n)
        dP_dtheta_fd = np.zeros(self.n)

        for i in range(self.n):
            theta = np.array(self.theta)
            delta = 1e-6*theta[i]
            theta[i] += delta

            rotor = CCBlade(self.r, self.chord, theta, self.af, self.Rhub, self.Rtip,
                self.B, self.rho, self.mu, self.precone, self.tilt, self.yaw, self.shearExp,
                self.hubHt, self.nSector, derivatives=False)

            Pd, Td, Qd = rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=False)

            dT_dtheta_fd[i] = (Td - self.T) / delta
            dQ_dtheta_fd[i] = (Qd - self.Q) / delta
            dP_dtheta_fd[i] = (Pd - self.P) / delta

        np.testing.assert_allclose(dT_dtheta_fd, dT_dtheta, rtol=7e-5, atol=1e-8)
        np.testing.assert_allclose(dQ_dtheta_fd, dQ_dtheta, rtol=7e-5, atol=1e-8)
        np.testing.assert_allclose(dP_dtheta_fd, dP_dtheta, rtol=7e-5, atol=1e-8)

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
    def test_dpresweep2(self):

        presweep = np.linspace(1, 10, self.n)
        presweepTip = 10.1
        precone = 0.0
        rotor = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip,
            self.B, self.rho, self.mu, precone, self.tilt, self.yaw, self.shearExp,
            self.hubHt, self.nSector, derivatives=True, presweep=presweep, presweepTip=presweepTip)

        P, T, Q, dP_ds, dT_ds, dQ_ds, dP_dv, dT_dv, dQ_dv = \
            rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=False)

        dT_dpresweep = dT_dv[0, 4, :]
        dQ_dpresweep = dQ_dv[0, 4, :]
        dP_dpresweep = dP_dv[0, 4, :]


        dT_dpresweep_fd = np.zeros(self.n)
        dQ_dpresweep_fd = np.zeros(self.n)
        dP_dpresweep_fd = np.zeros(self.n)
        for i in range(self.n):
            ps = np.array(presweep)
            delta = 1e-6*ps[i]
            ps[i] += delta

            rotor = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip,
                self.B, self.rho, self.mu, precone, self.tilt, self.yaw, self.shearExp,
                self.hubHt, self.nSector, derivatives=False, presweep=ps, presweepTip=presweepTip)

            Pd, Td, Qd = rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=False)

            dT_dpresweep_fd[i] = (Td - T) / delta
            dQ_dpresweep_fd[i] = (Qd - Q) / delta
            dP_dpresweep_fd[i] = (Pd - P) / delta


        np.testing.assert_allclose(dT_dpresweep_fd, dT_dpresweep, rtol=3e-4, atol=1e-8)
        np.testing.assert_allclose(dQ_dpresweep_fd, dQ_dpresweep, rtol=3e-4, atol=1e-8)
        np.testing.assert_allclose(dP_dpresweep_fd, dP_dpresweep, rtol=3e-4, atol=1e-8)

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
    def test_dprecurve3(self):

        precurve = np.linspace(1, 10, self.n)
        precurveTip = 10.1
        precone = 0.0
        rotor = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip,
            self.B, self.rho, self.mu, precone, self.tilt, self.yaw, self.shearExp,
            self.hubHt, self.nSector, derivatives=True, precurve=precurve, precurveTip=precurveTip)

        CP, CT, CQ, dCP_ds, dCT_ds, dCQ_ds, dCP_dv, dCT_dv, dCQ_dv = \
            rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=True)

        dCT_dprecurve = dCT_dv[0, 3, :]
        dCQ_dprecurve = dCQ_dv[0, 3, :]
        dCP_dprecurve = dCP_dv[0, 3, :]


        dCT_dprecurve_fd = np.zeros(self.n)
        dCQ_dprecurve_fd = np.zeros(self.n)
        dCP_dprecurve_fd = np.zeros(self.n)
        for i in range(self.n):
            pc = np.array(precurve)
            delta = 1e-6*pc[i]
            pc[i] += delta

            rotor = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip,
                self.B, self.rho, self.mu, precone, self.tilt, self.yaw, self.shearExp,
                self.hubHt, self.nSector, derivatives=False, precurve=pc, precurveTip=precurveTip)

            CPd, CTd, CQd = rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=True)

            dCT_dprecurve_fd[i] = (CTd - CT) / delta
            dCQ_dprecurve_fd[i] = (CQd - CQ) / delta
            dCP_dprecurve_fd[i] = (CPd - CP) / delta

        np.testing.assert_allclose(dCT_dprecurve_fd, dCT_dprecurve, rtol=3e-4, atol=1e-8)
        np.testing.assert_allclose(dCQ_dprecurve_fd, dCQ_dprecurve, rtol=3e-4, atol=1e-8)
        np.testing.assert_allclose(dCP_dprecurve_fd, dCP_dprecurve, rtol=3e-4, atol=1e-8)

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
    def test_dhubht2(self):

        dT_dhubht = self.dT_ds[0, 2]
        dQ_dhubht = self.dQ_ds[0, 2]
        dP_dhubht = self.dP_ds[0, 2]

        hubht = float(self.hubHt)
        delta = 1e-6*hubht
        hubht += delta

        rotor = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip,
            self.B, self.rho, self.mu, self.precone, self.tilt, self.yaw, self.shearExp,
            hubht, self.nSector, derivatives=False)

        Pd, Td, Qd = rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=False)

        dT_dhubht_fd = (Td - self.T) / delta
        dQ_dhubht_fd = (Qd - self.Q) / delta
        dP_dhubht_fd = (Pd - self.P) / delta

        np.testing.assert_allclose(dT_dhubht_fd, dT_dhubht, rtol=1e-5, atol=1e-8)
        np.testing.assert_allclose(dQ_dhubht_fd, dQ_dhubht, rtol=5e-5, atol=1e-8)
        np.testing.assert_allclose(dP_dhubht_fd, dP_dhubht, rtol=5e-5, atol=1e-8)

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
    def test_dhubht3(self):

        dCT_dhubht = self.dCT_ds[0, 2]
        dCQ_dhubht = self.dCQ_ds[0, 2]
        dCP_dhubht = self.dCP_ds[0, 2]

        hubht = float(self.hubHt)
        delta = 1e-6*hubht
        hubht += delta

        rotor = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip,
            self.B, self.rho, self.mu, self.precone, self.tilt, self.yaw, self.shearExp,
            hubht, self.nSector, derivatives=False)

        CPd, CTd, CQd = rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=True)

        dCT_dhubht_fd = (CTd - self.CT) / delta
        dCQ_dhubht_fd = (CQd - self.CQ) / delta
        dCP_dhubht_fd = (CPd - self.CP) / delta

        np.testing.assert_allclose(dCT_dhubht_fd, dCT_dhubht, rtol=1e-5, atol=1e-8)
        np.testing.assert_allclose(dCQ_dhubht_fd, dCQ_dhubht, rtol=5e-5, atol=1e-8)
        np.testing.assert_allclose(dCP_dhubht_fd, dCP_dhubht, rtol=5e-5, atol=1e-8)

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
class CCBladePower(Component):
    def __init__(self, naero, npower):
        super(CCBladePower, self).__init__()
        """blade element momentum code"""

        # inputs
        self.add_param('Uhub', val=np.zeros(npower), units='m/s', desc='hub height wind speed')
        self.add_param('Omega', val=np.zeros(npower), units='rpm', desc='rotor rotation speed')
        self.add_param('pitch', val=np.zeros(npower), units='deg', desc='blade pitch setting')

        # outputs
        self.add_output('T', val=np.zeros(npower), units='N', desc='rotor aerodynamic thrust')
        self.add_output('Q', val=np.zeros(npower), units='N*m', desc='rotor aerodynamic torque')
        self.add_output('P', val=np.zeros(npower), units='W', desc='rotor aerodynamic power')

        
        # (potential) variables
        self.add_param('r', val=np.zeros(naero), units='m', desc='radial locations where blade is defined (should be increasing and not go all the way to hub or tip)')
        self.add_param('chord', val=np.zeros(naero), units='m', desc='chord length at each section')
        self.add_param('theta', val=np.zeros(naero),  units='deg', desc='twist angle at each section (positive decreases angle of attack)')
        self.add_param('Rhub', val=0.0, units='m', desc='hub radius')
        self.add_param('Rtip', val=0.0, units='m', desc='tip radius')
        self.add_param('hubHt', val=0.0, units='m', desc='hub height')
        self.add_param('precone', val=0.0, desc='precone angle', units='deg')
        self.add_param('tilt', val=0.0, desc='shaft tilt', units='deg')
        self.add_param('yaw', val=0.0, desc='yaw error', units='deg')

        # TODO: I've not hooked up the gradients for these ones yet.
        self.add_param('precurve', val=np.zeros(naero), units='m', desc='precurve at each section')
        self.add_param('precurveTip', val=0.0, units='m', desc='precurve at tip')

        # parameters
        self.add_param('airfoils', val=[0]*naero, desc='CCAirfoil instances', pass_by_obj=True)
        self.add_param('B', val=0, desc='number of blades', pass_by_obj=True)
        self.add_param('rho', val=0.0, units='kg/m**3', desc='density of air')
        self.add_param('mu', val=0.0, units='kg/(m*s)', desc='dynamic viscosity of air')
        self.add_param('shearExp', val=0.0, desc='shear exponent')
        self.add_param('nSector', val=4, desc='number of sectors to divide rotor face into in computing thrust and power', pass_by_obj=True)
        self.add_param('tiploss', val=True, desc='include Prandtl tip loss model', pass_by_obj=True)
        self.add_param('hubloss', val=True, desc='include Prandtl hub loss model', pass_by_obj=True)
        self.add_param('wakerotation', val=True, desc='include effect of wake rotation (i.e., tangential induction factor is nonzero)', pass_by_obj=True)
        self.add_param('usecd', val=True, desc='use drag coefficient in computing induction factors', pass_by_obj=True)

        self.naero = naero
        self.deriv_options['form'] = 'central'
        self.deriv_options['step_calc'] = 'relative'
        
    def solve_nonlinear(self, params, unknowns, resids):

        self.r = params['r']
        self.chord = params['chord']
        self.theta = params['theta']
        self.Rhub = params['Rhub']
        self.Rtip = params['Rtip']
        self.hubHt = params['hubHt']
        self.precone = params['precone']
        self.tilt = params['tilt']
        self.yaw = params['yaw']
        self.precurve = params['precurve']
        self.precurveTip = params['precurveTip']
        self.airfoils = params['airfoils']
        self.B = params['B']
        self.rho = params['rho']
        self.mu = params['mu']
        self.shearExp = params['shearExp']
        self.nSector = params['nSector']
        self.tiploss = params['tiploss']
        self.hubloss = params['hubloss']
        self.wakerotation = params['wakerotation']
        self.usecd = params['usecd']
        self.Uhub = params['Uhub']
        self.Omega = params['Omega']
        self.pitch = params['pitch']
        
        self.ccblade = CCBlade(self.r, self.chord, self.theta, self.airfoils, self.Rhub, self.Rtip, self.B,
            self.rho, self.mu, self.precone, self.tilt, self.yaw, self.shearExp, self.hubHt,
            self.nSector, self.precurve, self.precurveTip, tiploss=self.tiploss, hubloss=self.hubloss,
            wakerotation=self.wakerotation, usecd=self.usecd, derivatives=True)

        # power, thrust, torque
        self.P, self.T, self.Q, self.M, self.dP, self.dT, self.dQ \
            = self.ccblade.evaluate(self.Uhub, self.Omega, self.pitch, coefficients=False)
        unknowns['T'] = self.T
        unknowns['Q'] = self.Q
        unknowns['P'] = self.P
        

    def list_deriv_vars(self):

        inputs = ('precone', 'tilt', 'hubHt', 'Rhub', 'Rtip', 'yaw',
                  'Uhub', 'Omega', 'pitch', 'r', 'chord', 'theta', 'precurve', 'precurveTip')
        outputs = ('P', 'T', 'Q')
        return inputs, outputs


    def linearize(self, params, unknowns, resids):

        dP = self.dP
        dT = self.dT
        dQ = self.dQ
#.........這裏部分代碼省略.........

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
plt.legend(loc='upper left')
plt.grid()
plt.show()
plt.savefig('compare_m4_m6_m66_yaw_m67_AeroForces.png')
plt.close(1)




# To get the power, thrust, and torque at the same conditions 
# (in both absolute and coefficient form), use the evaluate method. 
# This is generally used for generating power curves so it expects 
# array_like input. For this example a list of size one is used.
# Note that the outputs are numpy arrays (of length 1 for this example). 
#m1_P,  m1_T,  m1_Q  = m1_rotor.evaluate([m1_Uinf], [m1_Omega], [pitch])
m1_CP, m1_CT, m1_CQ = m1_rotor.evaluate([m1_Uinf], [m1_Omega], [pitch], coefficient=True)

#m4_P,  m4_T,  m4_Q  = m4_rotor.evaluate([m4_Uinf], [m4_Omega], [pitch])
m4_CP, m4_CT, m4_CQ = m4_rotor.evaluate([m4_Uinf], [m4_Omega], [pitch], coefficient=True)

#m6_P,  m6_T,  m6_Q  = m6_rotor.evaluate([m6_Uinf], [m6_Omega], [pitch])
m6_CP, m6_CT, m6_CQ = m6_rotor.evaluate([m6_Uinf], [m6_Omega], [pitch], coefficient=True)

#m66_P,  m66_T,  m66_Q  = m66_rotor.evaluate([m66_Uinf], [m66_Omega], [pitch])
m66_CP, m66_CT, m66_CQ = m66_rotor.evaluate([m66_Uinf], [m66_Omega], [pitch], coefficient=True)

m67_CP, m67_CT, m67_CQ = m67_rotor.evaluate([m67_Uinf], [m67_Omega], [pitch], coefficient=True)

print 'model1 CP =', m1_CP
print 'model1 CT =', m1_CT
print 'model1 CQ =', m1_CQ

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
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 = path.join(path.dirname(path.realpath(__file__)), '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 = -5.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)





    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

        P, T, Q, M = self.rotor.evaluate(Uinf, Omega, pitch)

        # 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)

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
class RegulatedPowerCurve(Component): # Implicit COMPONENT

    def __init__(self, naero, n_pc, n_pc_spline, regulation_reg_II5 = True, regulation_reg_III = False):
        super(RegulatedPowerCurve, self).__init__()

        # parameters
        self.add_param('control_Vin',        val=0.0, units='m/s',  desc='cut-in wind speed')
        self.add_param('control_Vout',       val=0.0, units='m/s',  desc='cut-out wind speed')
        self.add_param('control_ratedPower', val=0.0, units='W',    desc='electrical rated power')
        self.add_param('control_minOmega',   val=0.0, units='rpm',  desc='minimum allowed rotor rotation speed')
        self.add_param('control_maxOmega',   val=0.0, units='rpm',  desc='maximum allowed rotor rotation speed')
        self.add_param('control_maxTS',      val=0.0, units='m/s',  desc='maximum allowed blade tip speed')
        self.add_param('control_tsr',        val=0.0,               desc='tip-speed ratio in Region 2 (should be optimized externally)')
        self.add_param('control_pitch',      val=0.0, units='deg',  desc='pitch angle in region 2 (and region 3 for fixed pitch machines)')
        self.add_param('drivetrainType',     val=DRIVETRAIN_TYPE['GEARED'], pass_by_obj=True)
        self.add_param('drivetrainEff',     val=0.0,               desc='overwrite drivetrain model with a given efficiency, used for FAST analysis')
        
        self.add_param('r',         val=np.zeros(naero), units='m',   desc='radial locations where blade is defined (should be increasing and not go all the way to hub or tip)')
        self.add_param('chord',     val=np.zeros(naero), units='m',   desc='chord length at each section')
        self.add_param('theta',     val=np.zeros(naero), units='deg', desc='twist angle at each section (positive decreases angle of attack)')
        self.add_param('Rhub',      val=0.0,             units='m',   desc='hub radius')
        self.add_param('Rtip',      val=0.0,             units='m',   desc='tip radius')
        self.add_param('hubHt',     val=0.0,             units='m',   desc='hub height')
        self.add_param('precone',   val=0.0,             units='deg', desc='precone angle', )
        self.add_param('tilt',      val=0.0,             units='deg', desc='shaft tilt', )
        self.add_param('yaw',       val=0.0,             units='deg', desc='yaw error', )
        self.add_param('precurve',      val=np.zeros(naero),    units='m', desc='precurve at each section')
        self.add_param('precurveTip',   val=0.0,                units='m', desc='precurve at tip')

        self.add_param('airfoils',  val=[0]*naero,                      desc='CCAirfoil instances', pass_by_obj=True)
        self.add_param('B',         val=0,                              desc='number of blades', pass_by_obj=True)
        self.add_param('rho',       val=0.0,        units='kg/m**3',    desc='density of air')
        self.add_param('mu',        val=0.0,        units='kg/(m*s)',   desc='dynamic viscosity of air')
        self.add_param('shearExp',  val=0.0,                            desc='shear exponent')
        self.add_param('nSector',   val=4,                              desc='number of sectors to divide rotor face into in computing thrust and power', pass_by_obj=True)
        self.add_param('tiploss',   val=True,                           desc='include Prandtl tip loss model', pass_by_obj=True)
        self.add_param('hubloss',   val=True,                           desc='include Prandtl hub loss model', pass_by_obj=True)
        self.add_param('wakerotation', val=True,                        desc='include effect of wake rotation (i.e., tangential induction factor is nonzero)', pass_by_obj=True)
        self.add_param('usecd',     val=True,                           desc='use drag coefficient in computing induction factors', pass_by_obj=True)

        # outputs
        self.add_output('V',        val=np.zeros(n_pc), units='m/s',    desc='wind vector')
        self.add_output('Omega',    val=np.zeros(n_pc), units='rpm',    desc='rotor rotational speed')
        self.add_output('pitch',    val=np.zeros(n_pc), units='deg',    desc='rotor pitch schedule')
        self.add_output('P',        val=np.zeros(n_pc), units='W',      desc='rotor electrical power')
        self.add_output('T',        val=np.zeros(n_pc), units='N',      desc='rotor aerodynamic thrust')
        self.add_output('Q',        val=np.zeros(n_pc), units='N*m',    desc='rotor aerodynamic torque')
        self.add_output('M',        val=np.zeros(n_pc), units='N*m',    desc='blade root moment')
        self.add_output('Cp',       val=np.zeros(n_pc),                 desc='rotor electrical power coefficient')
        self.add_output('V_spline', val=np.zeros(n_pc_spline), units='m/s',  desc='wind vector')
        self.add_output('P_spline', val=np.zeros(n_pc_spline), units='W',    desc='rotor electrical power')
        self.add_output('V_R25',       val=0.0, units='m/s', desc='region 2.5 transition wind speed')
        self.add_output('rated_V',     val=0.0, units='m/s', desc='rated wind speed')
        self.add_output('rated_Omega', val=0.0, units='rpm', desc='rotor rotation speed at rated')
        self.add_output('rated_pitch', val=0.0, units='deg', desc='pitch setting at rated')
        self.add_output('rated_T',     val=0.0, units='N',   desc='rotor aerodynamic thrust at rated')
        self.add_output('rated_Q',     val=0.0, units='N*m', desc='rotor aerodynamic torque at rated')
        self.add_output('ax_induct_cutin',   val=np.zeros(naero), desc='rotor axial induction at cut-in wind speed along blade span')
        self.add_output('tang_induct_cutin', val=np.zeros(naero), desc='rotor tangential induction at cut-in wind speed along blade span')
        self.add_output('aoa_cutin',         val=np.zeros(naero), desc='angle of attack distribution along blade span at cut-in wind speed')

        self.naero                      = naero
        self.n_pc                       = n_pc
        self.n_pc_spline                = n_pc_spline
        self.lock_pitchII               = False
        self.regulation_reg_II5         = regulation_reg_II5
        self.regulation_reg_III         = regulation_reg_III
        self.deriv_options['form']      = 'central'
        self.deriv_options['step_calc'] = 'relative'
        
    def solve_nonlinear(self, params, unknowns, resids):
                
        self.ccblade = CCBlade(params['r'], params['chord'], params['theta'], params['airfoils'], params['Rhub'], params['Rtip'], params['B'], params['rho'], params['mu'], params['precone'], params['tilt'], params['yaw'], params['shearExp'], params['hubHt'], params['nSector'])
        
        Uhub    = np.linspace(params['control_Vin'],params['control_Vout'], self.n_pc)
        
        P_aero   = np.zeros_like(Uhub)
        Cp_aero  = np.zeros_like(Uhub)
        P       = np.zeros_like(Uhub)
        Cp      = np.zeros_like(Uhub)
        T       = np.zeros_like(Uhub)
        Q       = np.zeros_like(Uhub)
        M       = np.zeros_like(Uhub)
        Omega   = np.zeros_like(Uhub)
        pitch   = np.zeros_like(Uhub) + params['control_pitch']

        Omega_max = min([params['control_maxTS'] / params['Rtip'], params['control_maxOmega']*np.pi/30.])
        
        # Region II
        for i in range(len(Uhub)):
            Omega[i] = Uhub[i] * params['control_tsr'] / params['Rtip']
        
        P_aero, T, Q, M, Cp_aero, _, _, _ = self.ccblade.evaluate(Uhub, Omega * 30. / np.pi, pitch, coefficients=True)
        P, eff  = CSMDrivetrain(P_aero, params['control_ratedPower'], params['drivetrainType'])
        Cp      = Cp_aero*eff
        
        
        
        # search for Region 2.5 bounds
        for i in range(len(Uhub)):
#.........這裏部分代碼省略.........

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
class TestGradients(unittest.TestCase):

    def setUp(self):

        # geometry
        self.Rhub = 1.5
        self.Rtip = 63.0

        self.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])
        self.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])
        self.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])
        self.B = 3  # number of blades

        # atmosphere
        self.rho = 1.225
        self.mu = 1.81206e-5

        afinit = CCAirfoil.initFromAerodynFile  # just for shorthand
        basepath = path.join(path.dirname(path.realpath(__file__)), '5MW_AFFiles') + path.sep

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

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


        self.tilt = -5.0
        self.precone = 2.5
        self.yaw = 0.0
        self.shearExp = 0.2
        self.hubHt = 80.0
        self.nSector = 8

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

        # set conditions
        self.Uinf = 10.0
        tsr = 7.55
        self.pitch = 0.0
        self.Omega = self.Uinf*tsr/self.Rtip * 30.0/pi  # convert to RPM
        self.azimuth = 90


        self.Np, self.Tp, self.dNp_dX, self.dTp_dX, self.dNp_dprecurve, self.dTp_dprecurve = \
            self.rotor.distributedAeroLoads(self.Uinf, self.Omega, self.pitch, self.azimuth)

        self.P, self.T, self.Q, self.dP_ds, self.dT_ds, self.dQ_ds, self.dP_dv, self.dT_dv, \
            self.dQ_dv = self.rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=False)

        self.CP, self.CT, self.CQ, self.dCP_ds, self.dCT_ds, self.dCQ_ds, self.dCP_dv, self.dCT_dv, \
            self.dCQ_dv = self.rotor.evaluate([self.Uinf], [self.Omega], [self.pitch], coefficient=True)

        self.rotor.derivatives = False
        self.n = len(self.r)

        # X = [r, chord, theta, Rhub, Rtip, presweep, precone, tilt, hubHt]
        # scalars = [precone, tilt, hubHt, Rhub, Rtip, precurvetip, presweeptip]
        # vectors = [r, chord, theta, precurve, presweep]

    def test_dr1(self):

        dNp_dr = self.dNp_dX[0, :]
        dTp_dr = self.dTp_dX[0, :]
        dNp_dr_fd = np.zeros(self.n)
        dTp_dr_fd = np.zeros(self.n)

        for i in range(self.n):
            r = np.array(self.r)
            delta = 1e-6*r[i]
            r[i] += delta

            rotor = CCBlade(r, self.chord, self.theta, self.af, self.Rhub, self.Rtip,
                self.B, self.rho, self.mu, self.precone, self.tilt, self.yaw, self.shearExp,
                self.hubHt, self.nSector, derivatives=False)

            Npd, Tpd = rotor.distributedAeroLoads(self.Uinf, self.Omega, self.pitch, self.azimuth)

            dNp_dr_fd[i] = (Npd[i] - self.Np[i]) / delta
            dTp_dr_fd[i] = (Tpd[i] - self.Tp[i]) / delta

#.........這裏部分代碼省略.........

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
    precone = 2.5
    yaw = 0.0
    shearExp = 0.2
    hubHt = high
    nSector = 8

    # create CCBlade object
    aeroanalysis = CCBlade(rad, c, alpha, af, Rhub, Rtip, blade, rho, mu, precone, tilt, yaw, shearExp, hubHt, nSector)

    tsr = np.linspace(0,20,100) # tip-speed ratio
    Uinf = Vinf*np.ones_like(tsr) # free stream wind speed
    Omega = ((Uinf*tsr)/Rtip)*(30./np.pi) # rotation rate (rad/s)
    pitch = np.ones_like(tsr)*0. # pitch (deg)

    # Calculating power coefficients at each tip-speed ratio
    CP,_,_ = aeroanalysis.evaluate(Uinf, Omega, pitch, coefficient=True)

    # Creating a power curve for the turbine (tip-speed ratio vs. power coefficient)
    powercurve = Akima1DInterpolator(tsr,CP)

    # Adjusting geometry for SPL calculations
    r_nrel = 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, 63.0]) # radial positions (m)
    rad = r_nrel*r_ratio

    # Initialize input variables
    rotorDiameter = np.ones(nturb)*rotor_diameter
    generator_efficiency = np.ones(nturb)*0.944
    yaw = np.ones((nwind,nturb))*0.
    rpm = np.ones(nwind*nturb)*rpm_max

    # Optimization

# 需要導入模塊: from ccblade import CCBlade [as 別名]
# 或者: from ccblade.CCBlade import evaluate [as 別名]
dNp_dtilt = dNp['dtilt']
dNp_dhubHt = dNp['dhubHt']
dNp_dyaw = dNp['dyaw']
dNp_dazimuth = dNp['dazimuth']
dNp_dUinf = dNp['dUinf']
dNp_dOmega = dNp['dOmega']
dNp_dpitch = dNp['dpitch']

# 5 ----------


# 6 ----------


P, T, Q, dP, dT, dQ \
    = rotor.evaluate([Uinf], [Omega], [pitch])

npts = len(P)

# npts x 1
dP_dprecone = dP['dprecone']
dP_dtilt = dP['dtilt']
dP_dhubHt = dP['dhubHt']
dP_dRhub = dP['dRhub']
dP_dRtip = dP['dRtip']
dP_dprecurveTip = dP['dprecurveTip']
dP_dpresweepTip = dP['dpresweepTip']
dP_dyaw = dP['dyaw']

# npts x npts
dP_dUinf = dP['dUinf']