Python CCBlade.distributedAeroLoads方法代码示例

# 需要导入模块: from ccblade import CCBlade [as 别名]
# 或者: from ccblade.CCBlade import distributedAeroLoads [as 别名]
    def test_dprecurveTip1(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)

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

        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)

        Npd, Tpd = rotor.distributedAeroLoads(self.Uinf, self.Omega, self.pitch, self.azimuth)
        dNp_dprecurveTip_fd = (Npd - Np) / delta
        dTp_dprecurveTip_fd = (Tpd - Tp) / delta

        np.testing.assert_allclose(dNp_dprecurveTip_fd, 0.0, rtol=1e-4, atol=1e-8)
        np.testing.assert_allclose(dTp_dprecurveTip_fd, 0.0, rtol=1e-4, atol=1e-8)

# 需要导入模块: from ccblade import CCBlade [as 别名]
# 或者: from ccblade.CCBlade import distributedAeroLoads [as 别名]
    def test_dprecurve1(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)

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

        dNp_dprecurve_fd = np.zeros((self.n, self.n))
        dTp_dprecurve_fd = np.zeros((self.n, 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)

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

            dNp_dprecurve_fd[i, :] = (Npd - Np) / delta
            dTp_dprecurve_fd[i, :] = (Tpd - Tp) / delta

        np.testing.assert_allclose(dNp_dprecurve_fd, dNp_dprecurve, rtol=3e-4, atol=1e-8)
        np.testing.assert_allclose(dTp_dprecurve_fd, dTp_dprecurve, rtol=3e-4, atol=1e-8)

# 需要导入模块: from ccblade import CCBlade [as 别名]
# 或者: from ccblade.CCBlade import distributedAeroLoads [as 别名]
    def test_dpresweep1(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)

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

        dNp_dpresweep = dNp_dX[5, :]
        dTp_dpresweep = dTp_dX[5, :]

        dNp_dpresweep_fd = np.zeros(self.n)
        dTp_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)

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

            dNp_dpresweep_fd[i] = (Npd[i] - Np[i]) / delta
            dTp_dpresweep_fd[i] = (Tpd[i] - Tp[i]) / delta

        np.testing.assert_allclose(dNp_dpresweep_fd, dNp_dpresweep, rtol=1e-5, atol=1e-8)
        np.testing.assert_allclose(dTp_dpresweep_fd, dTp_dpresweep, rtol=1e-5, atol=1e-8)

# 需要导入模块: from ccblade import CCBlade [as 别名]
# 或者: from ccblade.CCBlade import distributedAeroLoads [as 别名]
    def test_dhubht1(self):

        dNp_dhubht = self.dNp_dX[8, :]
        dTp_dhubht = self.dTp_dX[8, :]

        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)

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

        dNp_dhubht_fd = (Npd - self.Np) / delta
        dTp_dhubht_fd = (Tpd - self.Tp) / delta

        np.testing.assert_allclose(dNp_dhubht_fd, dNp_dhubht, rtol=1e-5, atol=1e-8)
        np.testing.assert_allclose(dTp_dhubht_fd, dTp_dhubht, rtol=1e-5, atol=1e-8)

# 需要导入模块: from ccblade import CCBlade [as 别名]
# 或者: from ccblade.CCBlade import distributedAeroLoads [as 别名]
    def test_dtheta1(self):

        dNp_dtheta = self.dNp_dX[2, :]
        dTp_dtheta = self.dTp_dX[2, :]
        dNp_dtheta_fd = np.zeros(self.n)
        dTp_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)

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

            dNp_dtheta_fd[i] = (Npd[i] - self.Np[i]) / delta
            dTp_dtheta_fd[i] = (Tpd[i] - self.Tp[i]) / delta


        np.testing.assert_allclose(dNp_dtheta_fd, dNp_dtheta, rtol=1e-6, atol=1e-8)
        np.testing.assert_allclose(dTp_dtheta_fd, dTp_dtheta, rtol=1e-4, atol=1e-8)

# 需要导入模块: from ccblade import CCBlade [as 别名]
# 或者: from ccblade.CCBlade import distributedAeroLoads [as 别名]
class CCBlade(AeroBase):
    """blade element momentum code"""

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

    # TODO: I've not hooked up the gradients for these ones yet.
    precurve = Array(iotype='in', units='m', desc='precurve at each section')
    precurveTip = Float(0.0, iotype='in', units='m', desc='precurve at tip')

    # parameters
    airfoil_files = List(Str, iotype='in', desc='names of airfoil file')
    coordinate_files = List(Str, iotype='in', desc='names of airfoil file')
    B = Int(3, iotype='in', desc='number of blades')
    rho = Float(1.225, iotype='in', units='kg/m**3', desc='density of air')
    mu = Float(1.81206e-5, iotype='in', units='kg/(m*s)', desc='dynamic viscosity of air')
    shearExp = Float(0.2, iotype='in', desc='shear exponent')
    nSector = Int(4, iotype='in', desc='number of sectors to divide rotor face into in computing thrust and power')
    tiploss = Bool(True, iotype='in', desc='include Prandtl tip loss model')
    hubloss = Bool(True, iotype='in', desc='include Prandtl hub loss model')
    wakerotation = Bool(True, iotype='in', desc='include effect of wake rotation (i.e., tangential induction factor is nonzero)')
    usecd = Bool(True, iotype='in', desc='use drag coefficient in computing induction factors')
    af = Array(iotype='in', desc='CCBlade objects')

    missing_deriv_policy = 'assume_zero'

    def execute(self):

        if len(self.precurve) == 0:
            self.precurve = np.zeros_like(self.r)

        af = self.af

        self.ccblade = CCBlade_PY(self.r, self.chord, self.theta, af, 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)


        if self.run_case == 'power':

            # power, thrust, torque
            self.P, self.T, self.Q, self.dP, self.dT, self.dQ \
                = self.ccblade.evaluate(self.Uhub, self.Omega, self.pitch, coefficient=False)


        elif self.run_case == 'loads':

            # distributed loads
            Np, Tp, self.dNp, self.dTp \
                = self.ccblade.distributedAeroLoads(self.V_load, self.Omega_load, self.pitch_load, self.azimuth_load)

            # concatenate loads at root/tip
            self.loads.r = np.concatenate([[self.Rhub], self.r, [self.Rtip]])
            Np = np.concatenate([[0.0], Np, [0.0]])
            Tp = np.concatenate([[0.0], Tp, [0.0]])

            # conform to blade-aligned coordinate system
            self.loads.Px = Np
            self.loads.Py = -Tp
            self.loads.Pz = 0*Np

            # return other outputs needed
            self.loads.V = self.V_load
            self.loads.Omega = self.Omega_load
            self.loads.pitch = self.pitch_load
            self.loads.azimuth = self.azimuth_load


    def list_deriv_vars(self):

        if self.run_case == 'power':
            inputs = ('precone', 'tilt', 'hubHt', 'Rhub', 'Rtip', 'yaw',
                'Uhub', 'Omega', 'pitch', 'r', 'chord', 'theta', 'precurve', 'precurveTip')
            outputs = ('P', 'T', 'Q')

        elif self.run_case == 'loads':

            inputs = ('r', 'chord', 'theta', 'Rhub', 'Rtip', 'hubHt', 'precone',
                'tilt', 'yaw', 'V_load', 'Omega_load', 'pitch_load', 'azimuth_load', 'precurve')
            outputs = ('loads.r', 'loads.Px', 'loads.Py', 'loads.Pz', 'loads.V',
                'loads.Omega', 'loads.pitch', 'loads.azimuth')

        return inputs, outputs


    def provideJ(self):

        if self.run_case == 'power':

            dP = self.dP
            dT = self.dT
#.........这里部分代码省略.........

# 需要导入模块: from ccblade import CCBlade [as 别名]
# 或者: from ccblade.CCBlade import distributedAeroLoads [as 别名]
# set conditions
Uinf = 10.0
tsr = 7.55
pitch = 0.0
Omega = Uinf*tsr/Rtip * 30.0/pi  # convert to RPM
azimuth = 0.0

# # evaluate distributed loads
# Np, Tp = rotor.distributedAeroLoads(Uinf, Omega, pitch, azimuth)

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

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

Np, Tp, dNp, dTp \
    = rotor.distributedAeroLoads(Uinf, Omega, pitch, azimuth)

n = len(r)

# n x n (diagonal)
dNp_dr = dNp['dr']
dNp_dchord = dNp['dchord']
dNp_dtheta = dNp['dtheta']
dNp_dpresweep = dNp['dpresweep']

# n x n (tridiagonal)
dNp_dprecurve = dNp['dprecurve']

# n x 1
dNp_dRhub = dNp['dRhub']
dNp_dRtip = dNp['dRtip']

# 需要导入模块: from ccblade import CCBlade [as 别名]
# 或者: from ccblade.CCBlade import distributedAeroLoads [as 别名]

#.........这里部分代码省略.........
        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.V_load = params['V_load']
        self.Omega_load = params['Omega_load']
        self.pitch_load = params['pitch_load']
        self.azimuth_load = params['azimuth_load']


        if len(self.precurve) == 0:
            self.precurve = np.zeros_like(self.r)

        # airfoil files
        n = len(self.airfoils)
        af = self.airfoils
        # af = [0]*n
        # for i in range(n):
        #     af[i] = CCAirfoil.initFromAerodynFile(self.airfoil_files[i])

        self.ccblade = CCBlade(self.r, self.chord, self.theta, af, 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)

        # distributed loads
        Np, Tp, self.dNp, self.dTp \
            = self.ccblade.distributedAeroLoads(self.V_load, self.Omega_load, self.pitch_load, self.azimuth_load)

        # concatenate loads at root/tip
        unknowns['loads_r'] = self.r

        # conform to blade-aligned coordinate system
        unknowns['loads_Px'] = Np
        unknowns['loads_Py'] = -Tp
        unknowns['loads_Pz'] = 0*Np

        # return other outputs needed
        unknowns['loads_V'] = self.V_load
        unknowns['loads_Omega'] = self.Omega_load
        unknowns['loads_pitch'] = self.pitch_load
        unknowns['loads_azimuth'] = self.azimuth_load


    def list_deriv_vars(self):
        inputs = ('r', 'chord', 'theta', 'Rhub', 'Rtip', 'hubHt', 'precone',
                  'tilt', 'yaw', 'V_load', 'Omega_load', 'pitch_load', 'azimuth_load', 'precurve')
        outputs = ('loads_r', 'loads_Px', 'loads_Py', 'loads_Pz', 'loads_V',
                   'loads_Omega', 'loads_pitch', 'loads_azimuth')

        return inputs, outputs


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

        dNp = self.dNp
        dTp = self.dTp
        n = len(self.r)

        dr_dr = np.vstack([np.zeros(n), np.eye(n), np.zeros(n)])

# 需要导入模块: from ccblade import CCBlade [as 别名]
# 或者: from ccblade.CCBlade import distributedAeroLoads [as 别名]
                precone, tilt, 20.0, shearExp, m66_hubHt, nSector)
m67_rotor = CCBlade(m67_r, m67_chord, m67_theta, m67_af, m67_Rhub, m67_Rtip, B, m67_rho, m67_mu,
                precone, tilt, yaw, shearExp, m67_hubHt, nSector)
                
# convert to RPM
m1_Omega = m1_Uinf*tsr/m1_Rtip * 30.0/pi  
m4_Omega = m4_Uinf*tsr/m4_Rtip * 30.0/pi
m6_Omega = m6_Uinf*tsr/m6_Rtip * 30.0/pi
m66_Omega = m66_Uinf*tsr/m66_Rtip * 30.0/pi
m67_Omega = m67_Uinf*tsr/m67_Rtip * 30.0/pi

m66yaw_Uinf  = 1.2
m66yaw_Omega = m66yaw_Uinf*tsr/m66_Rtip * 30.0/pi

# Evaluate the distributed loads at a chosen set of operating conditions.
m1_Np, m1_Tp = m1_rotor.distributedAeroLoads(m1_Uinf, m1_Omega, pitch, azimuth)
m4_Np, m4_Tp = m4_rotor.distributedAeroLoads(m4_Uinf, m4_Omega, pitch, azimuth)
m6_Np, m6_Tp = m6_rotor.distributedAeroLoads(m6_Uinf, m6_Omega, pitch, azimuth)
m66_Np, m66_Tp = m66_rotor.distributedAeroLoads(m66_Uinf, m66_Omega, pitch, azimuth)
m66yaw20_Np, m66yaw20_Tp = m66yaw20_rotor.distributedAeroLoads(m66yaw_Uinf, m66yaw_Omega, pitch, azimuth)
m67_Np, m67_Tp = m67_rotor.distributedAeroLoads(m67_Uinf, m67_Omega, pitch, azimuth)








###############################################################################
# Plot the flapwise and lead-lag aerodynamic loading

# 需要导入模块: from ccblade import CCBlade [as 别名]
# 或者: from ccblade.CCBlade import distributedAeroLoads [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 distributedAeroLoads [as 别名]

#.........这里部分代码省略.........
            Cp[i]           = Cp_aero[i]*eff
            P[i]            = params['control_ratedPower']
            
            
        else:
            # Rated conditions
            def get_Uhub_rated_noII12(pitch, Uhub):
                Uhub_i  = Uhub
                Omega_i = min([Uhub_i * params['control_tsr'] / params['Rtip'], Omega_max])
                pitch_i = pitch           
                P_aero_i, _, _, _ = self.ccblade.evaluate([Uhub_i], [Omega_i * 30. / np.pi], [pitch_i], coefficients=False)
                P_i, eff          = CSMDrivetrain(P_aero_i, params['control_ratedPower'], params['drivetrainType'])
                return abs(P_i - params['control_ratedPower'])
            
            bnds     = [Uhub[i-1], Uhub[i+1]]
            U_rated  = minimize_scalar(lambda x: get_Uhub_rated_noII12(pitch[i], x), bounds=bnds, method='bounded', options=options)['x']
            
            if not np.isnan(U_rated):
                Uhub[i]         = U_rated
            else:
                print('Regulation trajectory is struggling to find a solution for rated wind speed. Check rotor_aeropower.py')
            
            
            
            
            
            Omega[i] = min([Uhub[i] * params['control_tsr'] / params['Rtip'], Omega_max])
            pitch0   = pitch[i]
            
            P_aero[i], T[i], Q[i], M[i], Cp_aero[i], _, _, _ = self.ccblade.evaluate([Uhub[i]], [Omega[i] * 30. / np.pi], [pitch0], coefficients=True)
            P[i], eff    = CSMDrivetrain(P_aero[i], params['control_ratedPower'], params['drivetrainType'])
            Cp[i]        = Cp_aero[i]*eff
        
        
        for j in range(i + 1,len(Uhub)):
            Omega[j] = Omega[i]
            if self.regulation_reg_III == True:
                
                pitch0   = pitch[j-1]
                bnds     = [pitch0, pitch0 + 15.]
                pitch_regionIII = minimize_scalar(lambda x: constantPregionIII(x, Uhub[j], Omega[j]), bounds=bnds, method='bounded', options=options)['x']
                pitch[j]        = pitch_regionIII
                
                P_aero[j], T[j], Q[j], M[j], Cp_aero[j], _, _, _ = self.ccblade.evaluate([Uhub[j]], [Omega[j] * 30. / np.pi], [pitch[j]], coefficients=True)
                P[j], eff       = CSMDrivetrain(P_aero[j], params['control_ratedPower'], params['drivetrainType'])
                Cp[j]           = Cp_aero[j]*eff


                if abs(P[j] - params['control_ratedPower']) > 1e+4:
                    print('The pitch in region III is not being determined correctly at wind speed ' + str(Uhub[j]) + ' m/s')
                    P[j]        = params['control_ratedPower']
                    T[j]        = T[j-1]
                    Q[j]        = P[j] / Omega[j]
                    M[j]        = M[j-1]
                    pitch[j]    = pitch[j-1]
                    Cp[j]       = P[j] / (0.5 * params['rho'] * np.pi * params['Rtip']**2 * Uhub[i]**3)

                P[j] = params['control_ratedPower']
                
            else:
                P[j]        = params['control_ratedPower']
                T[j]        = 0
                Q[j]        = Q[i]
                M[j]        = 0
                pitch[j]    = 0
                Cp[j]       = P[j] / (0.5 * params['rho'] * np.pi * params['Rtip']**2 * Uhub[i]**3)

        
        unknowns['T']       = T
        unknowns['Q']       = Q
        unknowns['Omega']   = Omega * 30. / np.pi


        unknowns['P']       = P  
        unknowns['Cp']      = Cp  
        unknowns['V']       = Uhub
        unknowns['M']       = M
        unknowns['pitch']   = pitch
        
        
        self.ccblade.induction_inflow = True
        a_regII, ap_regII, alpha_regII = self.ccblade.distributedAeroLoads(Uhub[0], Omega[0] * 30. / np.pi, pitch[0], 0.0)
        
        # Fit spline to powercurve for higher grid density
        spline   = PchipInterpolator(Uhub, P)
        V_spline = np.linspace(params['control_Vin'],params['control_Vout'], num=self.n_pc_spline)
        P_spline = spline(V_spline)
        
        # outputs
        idx_rated = list(Uhub).index(U_rated)
        unknowns['rated_V']     = U_rated
        unknowns['rated_Omega'] = Omega[idx_rated] * 30. / np.pi
        unknowns['rated_pitch'] = pitch[idx_rated]
        unknowns['rated_T']     = T[idx_rated]
        unknowns['rated_Q']     = Q[idx_rated]
        unknowns['V_spline']    = V_spline
        unknowns['P_spline']    = P_spline
        unknowns['ax_induct_cutin']   = a_regII
        unknowns['tang_induct_cutin'] = ap_regII
        unknowns['aoa_cutin']         = alpha_regII

# 需要导入模块: from ccblade import CCBlade [as 别名]
# 或者: from ccblade.CCBlade import distributedAeroLoads [as 别名]
rotor = CCBlade(r, chord, theta, af, Rhub, Rtip, B, rho, mu,
                precone, tilt, yaw, shearExp, hubHt, nSector)

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

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

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

# evaluate distributed loads
Np, Tp = rotor.distributedAeroLoads(Uinf, Omega, pitch, azimuth)

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

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

# plot
rstar = (r - Rhub) / (Rtip - Rhub)

# append zero at root and tip
rstar = np.concatenate([[0.0], rstar, [1.0]])
Np = np.concatenate([[0.0], Np, [0.0]])
Tp = np.concatenate([[0.0], Tp, [0.0]])

plt.plot(rstar, Tp/1e3, label='lead-lag')
plt.plot(rstar, Np/1e3, label='flapwise')