Nodes of different colours represent the following:
Solid arrows point from a parent (sub)module to the submodule which is descended from it. Dashed arrows point from a module being used to the module or program unit using it. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | WindTurbineNum | |||
| logical, | intent(in) | :: | RunFlag |
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
SUBROUTINE CalcWindTurbine (WindTurbineNum, RunFlag)
! SUBROUTINE INFORMATION:
! AUTHOR Daeho Kang
! DATE WRITTEN Octorber 2009
! MODIFIED na
! RE-ENGINEERED na
! PURPOSE OF THIS SUBROUTINE:
!
! METHODOLOGY EMPLOYED:
! na
! REFERENCES:
! Sathyajith Mathew. 2006. Wind Energy: Fundamental, Resource Analysis and Economics. Springer,
! Chap. 2, pp. 11-15
! Mazharul Islam, David S.K. Ting, and Amir Fartaj. 2008. Aerodynamic Models for Darrieus-type sSraight-bladed
! Vertical Axis Wind Turbines. Renewable & Sustainable Energy Reviews, Volume 12, pp.1087-1109
! USE STATEMENTS:
USE ScheduleManager, ONLY: GetCurrentScheduleValue
USE Psychrometrics, ONLY: PsyWFnTdbTwbPb, PsyRhoAirFnPbTdbW
USE DataEnvironment, ONLY: OutBaroPress, WindSpeed, WindSpeedAt, OutDryBulbTempAt, OutWetBulbTempAt, OutBaroPressAt
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
! SUBROUTINE ARGUMENT DEFINITIONS:
INTEGER, INTENT(IN) :: WindTurbineNum ! System is on
LOGICAL, INTENT(IN) :: RunFlag ! System is on
! SUBROUTINE PARAMETER DEFINITIONS:
REAL(r64), PARAMETER :: MaxTheta = 90.0d0 ! Maximum of theta
REAL(r64), PARAMETER :: MaxDegree = 360.0d0 ! Maximum limit of outdoor air wind speed in m/s
REAL(r64), PARAMETER :: PitchAngle = 0.0d0 ! No pitch control, i.e. maximum rotor speed
REAL(r64), PARAMETER :: SecInMin = 60.0d0
REAL(r64), PARAMETER :: MaxTSR = 12.0d0 ! Maximum of tip speed ratio
! INTERFACE BLOCK SPECIFICATIONS
! na
! DERIVED TYPE DEFINITIONS
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
REAL(r64) :: LocalWindSpeed ! Ambient wind speed at the specific height in m/s
REAL(r64) :: RotorH ! Height of the rotor in m
REAL(r64) :: RotorD ! Diameter of the rotor in m
REAL(r64) :: LocalHumRat ! Local humidity ratio of the air at the specific height
REAL(r64) :: LocalAirDensity ! Local density at the specific height in m
REAL(r64) :: PowerCoeff ! Power coefficient
REAL(r64) :: SweptArea ! Swept area of the rotor in m2
REAL(r64) :: WTPower ! Total Power generated by the turbine in the quasi-steady state in Watts
REAL(r64) :: Power ! Actual power of the turbine in Watts
REAL(r64) :: TipSpeedRatio ! Tip speed ratio (TSR)
REAL(r64) :: TipSpeedRatioAtI ! Tip speed ratio at the ith time step
REAL(r64) :: AzimuthAng ! Azimuth angle of blades in degree
REAL(r64) :: ChordalVel ! Chordal velocity component in m/s
REAL(r64) :: NormalVel ! Normal velocity component in m/s
REAL(r64) :: AngOfAttack ! Angle of attack of a single blade in degree
REAL(r64) :: RelFlowVel ! Relative flow velocity component in m/s
REAL(r64) :: TanForce ! Tangential force in N.m
REAL(r64) :: NorForce ! Normal force in N.m
REAL(r64) :: RotorVel ! Rotor velocity in m/s
REAL(r64) :: AvgTanForce ! Average tangential force in N.m
REAL(r64) :: Constant ! Constants within integrand of tangential force
REAL(r64) :: IntRelFlowVel ! Integration of relative flow velocity
REAL(r64) :: TotTorque ! Total torque for the number of blades
REAL(r64) :: Omega ! Angular velocity of rotor in rad/s
REAL(r64) :: TanForceCoeff ! Tnagential force coefficient
REAL(r64) :: NorForceCoeff ! Normal force coefficient
REAL(r64) :: Period ! Period of sine and cosine functions
REAL(r64) :: Integrand ! Integrand of tangential force
REAL(r64) :: C1 ! Empirical power coefficient C1
REAL(r64) :: C2 ! Empirical power coefficient C2
REAL(r64) :: C3 ! Empirical power coefficient C3
REAL(r64) :: C4 ! Empirical power coefficient C4
REAL(r64) :: C5 ! Empirical power coefficient C5
REAL(r64) :: C6 ! Empirical power coefficient C6
REAL(r64) :: LocalTemp ! Ambient air temperature at the height in degree C
REAL(r64) :: LocalPress ! Local atmospheric pressure at the particular height in Pa
REAL(r64) :: InducedVel ! Induced velocity on the rotor in m/s
!unused REAL(r64) :: SysEfficiency ! Overall wind turbine system efficiency including generator and inverter
REAL(r64) :: MaxPowerCoeff ! Maximum power coefficient
REAL(r64) :: RotorSpeed ! Speed of rotors
! Estimate local velocity and density
RotorH = WindTurbineSys(WindTurbineNum)%RotorHeight
RotorD = WindTurbineSys(WindTurbineNum)%RotorDiameter
RotorSpeed = WindTurbineSys(WindTurbineNum)%RatedRotorSpeed
LocalTemp = OutDryBulbTempAt(RotorH)
LocalPress = OutBaroPressAt(RotorH)
LocalHumRat = PsyWFnTdbTwbPb(LocalTemp,OutWetBulbTempAt(RotorH),LocalPress)
LocalAirDensity = PsyRhoAirFnPbTdbW(LocalPress,LocalTemp,LocalHumRat)
LocalWindSpeed = WindSpeedAt(RotorH)
LocalWindSpeed = LocalWindSpeed / WindTurbineSys(WindTurbineNum)%WSFactor
! Flow
! Check wind conditions for system operation
IF (GetCurrentScheduleValue(WindTurbineSys(WindTurbineNum)%SchedPtr) > 0 .AND. &
LocalWindSpeed > WindTurbineSys(WindTurbineNum)%CutInSpeed .AND. &
LocalWindSpeed < WindTurbineSys(WindTurbineNum)%CutOutSpeed) THEN
! System is on
Period = 2.0d0 * pi
Omega = (RotorSpeed * Period) / SecInMin
SweptArea = (Pi * RotorD**2.0d0) / 4
TipSpeedRatio = (Omega * (RotorD/2.0d0)) / LocalWindSpeed
! Limit maximum tip speed ratio
IF (TipSpeedRatio > WindTurbineSys(WindTurbineNum)%MaxTipSpeedRatio) THEN
TipSpeedRatio = WindTurbineSys(WindTurbineNum)%MaxTipSpeedRatio
END IF
SELECT CASE (WindTurbineSys(WindTurbineNum)%RotorType)
CASE (HAWT) ! Horizontal axis wind turbine
MaxPowerCoeff = WindTurbineSys(WindTurbineNum)%MaxPowerCoeff
! Check if empirical constants are available
C1 = WindTurbineSys(WindTurbineNum)%PowerCoeffC1
C2 = WindTurbineSys(WindTurbineNum)%PowerCoeffC2
C3 = WindTurbineSys(WindTurbineNum)%PowerCoeffC3
C4 = WindTurbineSys(WindTurbineNum)%PowerCoeffC4
C5 = WindTurbineSys(WindTurbineNum)%PowerCoeffC5
C6 = WindTurbineSys(WindTurbineNum)%PowerCoeffC6
IF (C1 > 0.0d0 .AND. C2 > 0.0d0 .AND. C3 > 0.0d0 .AND. C4 >= 0.0d0 .AND. C5 > 0.0d0 .AND. C6 > 0.0d0) THEN
! Analytical approximation
! Maximum power, i.e., rotor speed is at maximum, and pitch angle is zero
TipSpeedRatioAtI = 1.0d0 / ( (1.0d0 / (TipSpeedRatio + 0.08d0 * PitchAngle)) - (0.035d0 / (PitchAngle**3 + 1.0d0)))
PowerCoeff = C1 * ((C2 / TipSpeedRatioAtI) - (C3 * PitchAngle) - (C4 * PitchAngle**1.5d0) - C5) * &
(EXP(-(C6 / TipSpeedRatioAtI)))
IF (PowerCoeff > MaxPowerCoeff) THEN
PowerCoeff = MaxPowerCoeff
END IF
WTPower = 0.5d0 * LocalAirDensity * PowerCoeff * SweptArea * LocalWindSpeed**3
ELSE ! Simple approximation
WTPower = 0.5d0 * LocalAirDensity * SweptArea * LocalWindSpeed**3 * MaxPowerCoeff
PowerCoeff = MaxPowerCoeff
END IF
! Maximum of rated power
IF (LocalWindSpeed >= WindTurbineSys(WindTurbineNum)%RatedWindSpeed .OR. &
WTPower > WindTurbineSys(WindTurbineNum)%RatedPower) THEN
WTPower = WindTurbineSys(WindTurbineNum)%RatedPower
PowerCoeff = WTPower / (0.5d0 * LocalAirDensity * SweptArea * LocalWindSpeed**3.0d0)
END IF
! Recalculated Cp at the rated power
WindTurbineSys(WindTurbineNum)%PowerCoeff = PowerCoeff
CASE (VAWT) ! Vertical axis wind turbine
RotorVel = Omega * (RotorD / 2.0d0)
! Recalculated omega, if TSR is greater than the maximum
IF (TipSpeedRatio >= WindTurbineSys(WindTurbineNum)%MaxTipSpeedRatio) THEN
RotorVel = LocalWindSpeed * WindTurbineSys(WindTurbineNum)%MaxTipSpeedRatio
Omega = RotorVel / (RotorD / 2.0d0)
END IF
AzimuthAng = MaxDegree / WindTurbineSys(WindTurbineNum)%NumOfBlade
! Azimuth angle between zero and 90 degree
IF (AzimuthAng > MaxTheta ) THEN ! Number of blades is 2 or 3
AzimuthAng = AzimuthAng - MaxTheta
IF (AzimuthAng == MaxTheta) THEN ! 2 blades
AzimuthAng = 0.0d0
END IF
ELSE IF (AzimuthAng == MaxTheta) THEN ! 4 blades
AzimuthAng = 0.0d0
END IF
InducedVel = LocalWindSpeed * 2.0d0/3.0d0
! Velocity components
ChordalVel = RotorVel + InducedVel * COS(AzimuthAng * DegToRadians)
NormalVel = InducedVel * SIN(AzimuthAng * DegToRadians)
RelFlowVel = SQRT(ChordalVel**2 + NormalVel**2)
! Angle of attack
AngOfAttack = ATAN((SIN(AzimuthAng * DegToRadians) / ((RotorVel / LocalWindSpeed) / (InducedVel / LocalWindSpeed) + &
COS(AzimuthAng * DegToRadians))))
! Force coefficients
TanForceCoeff = ABS(WindTurbineSys(WindTurbineNum)%LiftCoeff * SIN(AngOfAttack * DegToRadians) - &
WindTurbineSys(WindTurbineNum)%DragCoeff * COS(AngOfAttack * DegToRadians))
NorForceCoeff = WindTurbineSys(WindTurbineNum)%LiftCoeff * COS(AngOfAttack * DegToRadians) + &
WindTurbineSys(WindTurbineNum)%DragCoeff * SIN(AngOfAttack * DegToRadians)
! Net tangential and normal forces
TanForce = TanForceCoeff * 0.5d0 * LocalAirDensity * WindTurbineSys(WindTurbineNum)%ChordArea * RelFlowVel**2
NorForce = NorForceCoeff * 0.5d0 * LocalAirDensity * WindTurbineSys(WindTurbineNum)%ChordArea * RelFlowVel**2
Constant = (1.0d0 / Period) * (TanForce / RelFlowVel**2)
! Relative flow velocity is the only function of theta in net tangential force
! Integral of cos(theta) on zero to 2pi goes to zero
! Integrate constants only
IntRelFlowVel = RotorVel**2 * Period + InducedVel**2 * Period
! Average tangential force on a single blade
AvgTanForce = Constant * IntRelFlowVel
TotTorque = WindTurbineSys(WindTurbineNum)%NumOfBlade * AvgTanForce * (RotorD / 2.0d0)
WTPower = TotTorque * Omega
! Check if power produced is greater than maximum or rated power
IF (WTPower > WindTurbineSys(WindTurbineNum)%RatedPower) THEN
WTPower = WindTurbineSys(WindTurbineNum)%RatedPower
END IF
WindTurbineSys(WindTurbineNum)%ChordalVel = ChordalVel
WindTurbineSys(WindTurbineNum)%NormalVel = NormalVel
WindTurbineSys(WindTurbineNum)%RelFlowVel = RelFlowVel
WindTurbineSys(WindTurbineNum)%TanForce = TanForce
WindTurbineSys(WindTurbineNum)%NorForce = NorForce
WindTurbineSys(WindTurbineNum)%TotTorque = TotTorque
END SELECT
IF (WTPower > WindTurbineSys(WindTurbineNum)%RatedPower) THEN
WTPower = WindTurbineSys(WindTurbineNum)%RatedPower
END IF
! Actual power generated by the wind turbine system
Power = WTPower * WindTurbineSys(WindTurbineNum)%SysEfficiency
WindTurbineSys(WindTurbineNum)%Power = Power
WindTurbineSys(WindTurbineNum)%TotPower = WTPower
WindTurbineSys(WindTurbineNum)%LocalWindSpeed = LocalWindSpeed
WindTurbineSys(WindTurbineNum)%LocalAirDensity = LocalAirDensity
WindTurbineSys(WindTurbineNum)%TipSpeedRatio = TipSpeedRatio
ELSE ! System is off
WindTurbineSys(WindTurbineNum)%Power = 0.0d0
WindTurbineSys(WindTurbineNum)%TotPower = 0.0d0
WindTurbineSys(WindTurbineNum)%PowerCoeff = 0.0d0
WindTurbineSys(WindTurbineNum)%LocalWindSpeed = LocalWindSpeed
WindTurbineSys(WindTurbineNum)%LocalAirDensity = LocalAirDensity
WindTurbineSys(WindTurbineNum)%TipSpeedRatio = 0.0d0
WindTurbineSys(WindTurbineNum)%ChordalVel = 0.0d0
WindTurbineSys(WindTurbineNum)%NormalVel = 0.0d0
WindTurbineSys(WindTurbineNum)%RelFlowVel = 0.0d0
WindTurbineSys(WindTurbineNum)%AngOfAttack = 0.0d0
WindTurbineSys(WindTurbineNum)%TanForce = 0.0d0
WindTurbineSys(WindTurbineNum)%NorForce = 0.0d0
WindTurbineSys(WindTurbineNum)%TotTorque = 0.0d0
END IF
RETURN
END SUBROUTINE CalcWindTurbine