35 KiB
35 KiB
SUBROUTINE Coeff(p,InputFileData, ErrStat, ErrMsg)
!..................................................................................................................................
! Passed variables
TYPE(ED_ParameterType), INTENT(INOUT) :: p !< Parameters of the structural dynamics module
TYPE(ED_InputFile), INTENT(IN) :: InputFileData !< all the data in the ElastoDyn input file
INTEGER(IntKi), INTENT(OUT) :: ErrStat !< Error status
CHARACTER(*), INTENT(OUT) :: ErrMsg !< Error message when ErrStat =/ ErrID_None
! Local variables.
REAL(ReKi) :: AxRdBld (3,3) ! Temporary result holding the current addition to the p%AxRedBld() array.
REAL(ReKi) :: AxRdBldOld(3,3) ! Previous AxRdBld (i.e., AxRdBld from the previous node)
REAL(ReKi) :: AxRdTFA (2,2) ! Temporary result holding the current addition to the AxRedTFA() array.
REAL(ReKi) :: AxRdTFAOld(2,2) ! Previous AxRdTFA (i.e., AxRdTFA from the previous node)
REAL(ReKi) :: AxRdTSS (2,2) ! Temporary result holding the current addition to the AxRedTSS() array.
REAL(ReKi) :: AxRdTSSOld(2,2) ! Previous AxRdTSS (i.e., AxRdTSS from the previous node)
REAL(ReKi) :: TmpDist ! Temporary distance used in the calculation of the aero center locations.
REAL(ReKi) :: TmpDistj1 ! Temporary distance used in the calculation of the aero center locations.
REAL(ReKi) :: TmpDistj2 ! Temporary distance used in the calculation of the aero center locations.
REAL(ReKi) :: ElmntStff ! (Temporary) stiffness of an element.
REAL(ReKi) :: ElStffFA ! (Temporary) tower fore-aft stiffness of an element
REAL(ReKi) :: ElStffSS ! (Temporary) tower side-to-side stiffness of an element
REAL(ReKi) :: FMomAbvNd (p%NumBl,p%BldNodes) ! FMomAbvNd(K,J) = portion of the first moment of blade K about the rotor centerline (not root, like FirstMom(K)) associated with everything above node J (including tip brake masses).
REAL(ReKi) :: KBECent (p%NumBl,1,1) ! Centrifugal-term of generalized edgewise stiffness of the blades.
REAL(ReKi) :: KBFCent (p%NumBl,2,2) ! Centrifugal-term of generalized flapwise stiffness of the blades.
REAL(ReKi) :: KTFAGrav (2,2) ! Gravitational-term of generalized fore-aft stiffness of the tower.
REAL(ReKi) :: KTSSGrav (2,2) ! Gravitational-term of generalized side-to-side stiffness of the tower.
REAL(ReKi) :: MBE (p%NumBl,1,1) ! Generalized edgewise mass of the blades.
REAL(ReKi) :: MBF (p%NumBl,2,2) ! Generalized flapwise mass of the blades.
REAL(ReKi) :: MTFA (2,2) ! Generalized fore-aft mass of the tower.
REAL(ReKi) :: MTSS (2,2) ! Generalized side-to-side mass of the tower.
REAL(ReKi) :: Shape ! Temporary result holding a value from the SHP function
REAL(ReKi) :: Shape1 ! Temporary result holding a value from the SHP function
REAL(ReKi) :: Shape2 ! Temporary result holding a value from the SHP function
REAL(ReKi) :: TMssAbvNd (p%TwrNodes) ! Portion of the tower mass associated with everything above node J (including tower-top effects)
REAL(ReKi) :: TwstdSF (2,3,0:1) ! Temperory result holding the current addition to the TwistedSF() array.
REAL(ReKi) :: TwstdSFOld(2,3,0:1) ! Previous TwstdSF (i.e., TwstdSF from the previous node)
INTEGER(IntKi) :: I ! Generic index.
INTEGER(IntKi) :: J ! Loops through nodes / elements.
INTEGER(IntKi) :: K ! Loops through blades.
INTEGER(IntKi) :: L ! Generic index
ErrStat = ErrID_None
ErrMsg = ''
!...............................................................................................................................
! Calculate the distances from point S on a blade to the aerodynamic center in the j1 and j2 directions:
!...............................................................................................................................
DO K = 1,p%NumBl ! Loop through the blades
DO J = 1,p%BldNodes ! Loop through the blade nodes / elements
TmpDist = ( 0.25 - p%PitchAxis(K,J) )*p%Chord(J) ! Distance along the chordline from point S (25% chord) to the aerodynamic center of the blade element J--positive towards the trailing edge.
TmpDistj1 = TmpDist*p%SAeroTwst(J) ! Distance along the j1-axis from point S (25% chord) to the aerodynamic center of the blade element J
TmpDistj2 = TmpDist*p%CAeroTwst(J) ! Distance along the j2-axis from point S (25% chord) to the aerodynamic center of the blade element J
p%rSAerCenn1(K,J) = TmpDistj1*p%CThetaS(K,J) - TmpDistj2*p%SThetaS(K,J)
p%rSAerCenn2(K,J) = TmpDistj1*p%SThetaS(K,J) + TmpDistj2*p%CThetaS(K,J)
ENDDO ! J - Blade nodes / elements
ENDDO ! K - Blades
!...............................................................................................................................
! Calculate the structure that furls with the rotor inertia term:
!...............................................................................................................................
p%RrfaIner = InputFileData%RFrlIner - p%RFrlMass*( (p%rVDxn**2 )*( 1.0 - p%CRFrlSkw2*p%CRFrlTlt2 ) &
+ (p%rVDzn**2 )* p%CRFrlTlt2 &
+ (p%rVDyn**2 )*( 1.0 - p%SRFrlSkw2*p%CRFrlTlt2 ) &
- 2.0*p%rVDxn*p%rVDzn* p%CRFrlSkew*p%CSRFrlTlt &
- 2.0*p%rVDxn*p%rVDyn* p%CSRFrlSkw*p%CRFrlTlt2 &
- 2.0*p%rVDzn*p%rVDyn* p%SRFrlSkew*p%CSRFrlTlt )
IF ( p%RrfaIner < 0.0 ) THEN
ErrStat = ErrID_Fatal
ErrMsg = ' RFrlIner must not be less than RFrlMass*( perpendicular distance between rotor-furl'// &
' axis and CM of the structure that furls with the rotor [not including rotor] )^2.'
RETURN
END IF
!...............................................................................................................................
! Calculate the tail boom inertia term:
!...............................................................................................................................
p%AtfaIner = p%TFrlIner - p%BoomMass*( p%rWIxn*p%rWIxn*( 1.0 - p%CTFrlSkw2*p%CTFrlTlt2 ) &
+ p%rWIzn*p%rWIzn* p%CTFrlTlt2 &
+ p%rWIyn*p%rWIyn*( 1.0 - p%STFrlSkw2*p%CTFrlTlt2 ) &
- 2.0*p%rWIxn*p%rWIzn* p%CTFrlSkew*p%CSTFrlTlt &
- 2.0*p%rWIxn*p%rWIyn* p%CSTFrlSkw*p%CTFrlTlt2 &
- 2.0*p%rWIzn*p%rWIyn* p%STFrlSkew*p%CSTFrlTlt )
IF ( p%AtfaIner < 0.0 ) THEN
ErrStat = ErrID_Fatal
ErrMsg = ' TFrlIner must not be less than BoomMass*( perpendicular distance between tail-furl'// &
' axis and tail boom CM )^2.'
RETURN
ENDIF
!...............................................................................................................................
! Calculate the nacelle inertia terms:
!...............................................................................................................................
p%Nacd2Iner = InputFileData%NacYIner - p%NacMass*( p%NacCMxn**2 + p%NacCMyn**2 ) ! Nacelle inertia about the d2-axis
IF ( p%Nacd2Iner < 0.0 ) THEN
ErrStat = ErrID_Fatal
ErrMsg = ' NacYIner must not be less than NacMass*( NacCMxn^2 + NacCMyn^2 ).'
RETURN
END IF
! Calculate hub inertia about its centerline passing through its c.g..
! This calculation assumes that the hub for a 2-blader is essentially
! a uniform cylinder whose centerline is transverse through the cylinder
! passing through its c.g.. That is, for a 2-blader, Hubg1Iner =
! Hubg2Iner is the inertia of the hub about both the g1- and g2- axes. For
! 3-bladers, Hubg1Iner is simply equal to HubIner and Hubg2Iner is zero.
! Also, Initialize RotMass and RotIner to associated hub properties:
IF ( p%NumBl == 2 ) THEN ! 2-blader
p%Hubg1Iner = ( InputFileData%HubIner - p%HubMass*( ( p%UndSling - p%HubCM )**2 ) )/( p%CosDel3**2 )
p%Hubg2Iner = p%Hubg1Iner
IF ( p%Hubg1Iner < 0.0 ) THEN
ErrStat = ErrID_Fatal
ErrMsg = ' HubIner must not be less than HubMass*( UndSling - HubCM )^2 for 2-blader.'
RETURN
END IF
ELSE ! 3-blader
p%Hubg1Iner = InputFileData%HubIner
p%Hubg2Iner = 0.0
ENDIF
p%RotMass = p%HubMass
p%RotIner = p%Hubg1Iner
!...............................................................................................................................
! Initialize several variables to 0.0:
p%KBF = 0.0
p%KBE = 0.0
KBFCent = 0.0
KBECent = 0.0
p%TwrMass = 0.0
p%KTFA = 0.0
p%KTSS = 0.0
KTFAGrav = 0.0
KTSSGrav = 0.0
DO K = 1,p%NumBl ! Loop through the blades
! Initialize BldMass(), FirstMom(), and SecondMom() using TipMass() effects:
p%BldMass (K) = p%TipMass(K)
p%FirstMom (K) = p%TipMass(K)*p%BldFlexL
p%SecondMom(K) = p%TipMass(K)*p%BldFlexL*p%BldFlexL
DO J = p%BldNodes,1,-1 ! Loop through the blade nodes / elements in reverse
! Calculate the mass of the current element
p%BElmntMass(J,K) = p%MassB(K,J)*p%DRNodes(J) ! Mass of blade element J
! Integrate to find some blade properties which will be output in .fsm
p%BldMass (K) = p%BldMass (K) + p%BElmntMass(J,K)
p%FirstMom (K) = p%FirstMom (K) + p%BElmntMass(J,K)*p%RNodes(J)
p%SecondMom(K) = p%SecondMom(K) + p%BElmntMass(J,K)*p%RNodes(J)*p%RNodes(J)
! Integrate to find FMomAbvNd:
FMomAbvNd (K,J) = ( 0.5*p%BElmntMass(J,K) )*( p%HubRad + p%RNodes(J ) + 0.5*p%DRNodes(J ) )
IF ( J == p%BldNodes ) THEN ! Outermost blade element
! Add the TipMass() effects:
FMomAbvNd(K,J) = FmomAbvNd(K,J) + p%TipMass(K)*p%TipRad
ELSE ! All other blade elements
! Add to FMomAbvNd(K,J) the effects from the (not yet used) portion of element J+1
FMomAbvNd(K,J) = FMomAbvNd(K,J) + FMomAbvNd(K,J+1) &
+ ( 0.5*p%BElmntMass(J+1,K) )*( p%HubRad + p%RNodes(J+1) - 0.5*p%DRNodes(J+1) )
ENDIF
ENDDO ! J - Blade nodes / elements in reverse
IF (.NOT. p%BD4Blades) THEN
! Calculate BldCG() using FirstMom() and BldMass(); and calculate
! RotMass and RotIner:
p%BldCG (K) = p%FirstMom (K) / p%BldMass (K)
p%RotMass = p%RotMass + p%BldMass (K)
p%RotIner = p%RotIner + ( p%SecondMom(K) + p%BldMass (K)*p%HubRad*( 2.0*p%BldCG(K) + p%HubRad ) )*( p%CosPreC(K)**2 )
END IF
ENDDO ! K - Blades
DO K = 1,p%NumBl ! Loop through the blades
! Initialize the generalized blade masses using tip mass effects:
MBF(K,1,1) = p%TipMass(K)
MBF(K,2,2) = p%TipMass(K)
MBE(K,1,1) = p%TipMass(K)
DO J = 1,p%BldNodes ! Loop through the blade nodes / elements
! Integrate to find the generalized mass of the blade (including tip mass effects).
! Ignore the cross-correlation terms of MBF (i.e. MBF(i,j) where i /= j) since
! these terms will never be used.
Shape1 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl1Sh(:,K), 0, ErrStat, ErrMsg )
Shape2 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl2Sh(:,K), 0, ErrStat, ErrMsg )
MBF (K,1,1) = MBF (K,1,1) + p%BElmntMass(J,K)*Shape1*Shape1
MBF (K,2,2) = MBF (K,2,2) + p%BElmntMass(J,K)*Shape2*Shape2
Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldEdgSh(:,K), 0, ErrStat, ErrMsg )
MBE (K,1,1) = MBE (K,1,1) + p%BElmntMass(J,K)*Shape *Shape
! Integrate to find the generalized stiffness of the blade (not including centrifugal
! effects).
ElmntStff = p%StiffBF(K,J)*p%DRNodes(J) ! Flapwise stiffness of blade element J
Shape1 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl1Sh(:,K), 2, ErrStat, ErrMsg )
Shape2 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl2Sh(:,K), 2, ErrStat, ErrMsg )
p%KBF (K,1,1) = p%KBF (K,1,1) + ElmntStff*Shape1*Shape1
p%KBF (K,1,2) = p%KBF (K,1,2) + ElmntStff*Shape1*Shape2
p%KBF (K,2,1) = p%KBF (K,2,1) + ElmntStff*Shape2*Shape1
p%KBF (K,2,2) = p%KBF (K,2,2) + ElmntStff*Shape2*Shape2
ElmntStff = p%StiffBE(K,J)*p%DRNodes(J) ! Edgewise stiffness of blade element J
Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldEdgSh(:,K), 2, ErrStat, ErrMsg )
p%KBE (K,1,1) = p%KBE (K,1,1) + ElmntStff*Shape *Shape
! Integrate to find the centrifugal-term of the generalized flapwise and edgewise
! stiffness of the blades. Ignore the cross-correlation terms of KBFCent (i.e.
! KBFCent(i,j) where i /= j) since these terms will never be used.
ElmntStff = FMomAbvNd(K,J)*p%DRNodes(J)*p%RotSpeed**2 ! Centrifugal stiffness of blade element J
Shape1 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl1Sh(:,K), 1, ErrStat, ErrMsg )
Shape2 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl2Sh(:,K), 1, ErrStat, ErrMsg )
KBFCent(K,1,1) = KBFCent(K,1,1) + ElmntStff*Shape1*Shape1
KBFCent(K,2,2) = KBFCent(K,2,2) + ElmntStff*Shape2*Shape2
Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldEdgSh(:,K), 1, ErrStat, ErrMsg )
KBECent(K,1,1) = KBECent(K,1,1) + ElmntStff*Shape *Shape
! Calculate the 2nd derivatives of the twisted shape functions:
Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl1Sh(:,K), 2, ErrStat, ErrMsg )
p%TwistedSF(K,1,1,J,2) = Shape*p%CThetaS(K,J) ! 2nd deriv. of Phi1(J) for blade K
p%TwistedSF(K,2,1,J,2) = -Shape*p%SThetaS(K,J) ! 2nd deriv. of Psi1(J) for blade K
Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl2Sh(:,K), 2, ErrStat, ErrMsg )
p%TwistedSF(K,1,2,J,2) = Shape*p%CThetaS(K,J) ! 2nd deriv. of Phi2(J) for blade K
p%TwistedSF(K,2,2,J,2) = -Shape*p%SThetaS(K,J) ! 2nd deriv. of Psi2(J) for blade K
Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldEdgSh(:,K), 2, ErrStat, ErrMsg )
p%TwistedSF(K,1,3,J,2) = Shape*p%SThetaS(K,J) ! 2nd deriv. of Phi3(J) for blade K
p%TwistedSF(K,2,3,J,2) = Shape*p%CThetaS(K,J) ! 2nd deriv. of Psi3(J) for blade K
! Integrate to find the 1st derivatives of the twisted shape functions:
DO I = 1,2 ! Loop through Phi and Psi
DO L = 1,3 ! Loop through all blade DOFs
TwstdSF ( I,L, 1) = p%TwistedSF(K,I,L,J,2)*0.5*p%DRNodes(J)
p%TwistedSF (K,I,L,J,1) = TwstdSF ( I,L, 1)
ENDDO ! L - All blade DOFs
ENDDO ! I - Phi and Psi
IF ( J /= 1 ) THEN ! All but the innermost blade element
! Add the effects from the (not yet used) portion of element J-1
DO I = 1,2 ! Loop through Phi and Psi
DO L = 1,3 ! Loop through all blade DOFs
p%TwistedSF(K,I,L,J,1) = p%TwistedSF(K,I,L,J,1) + p%TwistedSF(K,I,L,J-1,1) &
+ TwstdSFOld( I,L, 1)
ENDDO ! L - All blade DOFs
ENDDO ! I - Phi and Psi
ENDIF
! Integrate to find the twisted shape functions themselves (i.e., their zeroeth derivative):
DO I = 1,2 ! Loop through Phi and Psi
DO L = 1,3 ! Loop through all blade DOFs
TwstdSF ( I,L, 0) = p%TwistedSF(K,I,L,J,1)*0.5*p%DRNodes(J)
p%TwistedSF (K,I,L,J,0) = TwstdSF ( I,L, 0)
ENDDO ! L - All blade DOFs
ENDDO ! I - Phi and Psi
IF ( J /= 1 ) THEN ! All but the innermost blade element
! Add the effects from the (not yet used) portion of element J-1
DO I = 1,2 ! Loop through Phi and Psi
DO L = 1,3 ! Loop through all blade DOFs
p%TwistedSF(K,I,L,J,0) = p%TwistedSF(K,I,L,J,0) + p%TwistedSF(K,I,L,J-1,0) &
+ TwstdSFOld( I,L, 0)
ENDDO ! L - All blade DOFs
ENDDO ! I - Phi and Psi
ENDIF
! Integrate to find the blade axial reduction shape functions:
DO I = 1,3 ! Loop through all blade DOFs
DO L = 1,3 ! Loop through all blade DOFs
AxRdBld ( I,L ) = 0.5*p%DRNodes(J)*( &
p%TwistedSF(K,1,I,J,1)*p%TwistedSF(K,1,L,J,1) &
+ p%TwistedSF(K,2,I,J,1)*p%TwistedSF(K,2,L,J,1) )
p%AxRedBld (K,I,L,J) = AxRdBld(I,L)
ENDDO ! L - All blade DOFs
ENDDO ! I - All blade DOFs
IF ( J /= 1 ) THEN ! All but the innermost blade element
! Add the effects from the (not yet used) portion of element J-1
DO I = 1,3 ! Loop through all blade DOFs
DO L = 1,3 ! Loop through all blade DOFs
p%AxRedBld(K,I,L,J) = p%AxRedBld(K,I,L,J) + p%AxRedBld(K,I,L,J-1) &
+ AxRdBldOld(I,L)
ENDDO ! L - All blade DOFs
ENDDO ! I - All blade DOFs
ENDIF
! Store the TwstdSF and AxRdBld terms of the current element (these will be used for the next element)
TwstdSFOld = TwstdSF
AxRdBldOld = AxRdBld
ENDDO ! J - Blade nodes / elements
IF (p%BD4Blades) THEN
!p%KBF ( K,:,: ) = 0.0_ReKi
! the 1st and zeroeth derivatives of the twisted shape functions at the blade root:
p%TwistedSF(K,:,:,:,1) = 0.0_ReKi
p%TwistedSF(K,:,:,:,0) = 0.0_ReKi
p%AxRedBld( K,:,:,: ) = 0.0_ReKi
ELSE
! Apply the flapwise modal stiffness tuners of the blades to KBF():
DO I = 1,2 ! Loop through flap DOFs
DO L = 1,2 ! Loop through flap DOFs
p%KBF(K,I,L) = SQRT( p%FStTunr(K,I)*p%FStTunr(K,L) )*p%KBF(K,I,L)
ENDDO ! L - Flap DOFs
ENDDO ! I - Flap DOFs
! Calculate the blade natural frequencies:
DO I = 1,NumBF ! Loop through flap DOFs
p%FreqBF(K,I,1) = Inv2Pi*SQRT( p%KBF(K,I,I) /( MBF(K,I,I) - p%TipMass(K) ) ) ! Natural blade I-flap frequency w/o centrifugal stiffening nor tip mass effects
p%FreqBF(K,I,2) = Inv2Pi*SQRT( p%KBF(K,I,I) / MBF(K,I,I) ) ! Natural blade I-flap frequency w/o centrifugal stiffening, but w/ tip mass effects
p%FreqBF(K,I,3) = Inv2Pi*SQRT( ( p%KBF(K,I,I) + KBFCent(K,I,I) )/ MBF(K,I,I) ) ! Natural blade I-flap frequency w/ centrifugal stiffening and tip mass effects
ENDDO ! I - Flap DOFs
p%FreqBE (K,1,1) = Inv2Pi*SQRT( p%KBE(K,1,1) /( MBE(K,1,1) - p%TipMass(K) ) ) ! Natural blade 1-edge frequency w/o centrifugal stiffening nor tip mass effects
p%FreqBE (K,1,2) = Inv2Pi*SQRT( p%KBE(K,1,1) / MBE(K,1,1) ) ! Natural Blade 1-edge frequency w/o centrifugal stiffening, but w/ tip mass effects
p%FreqBE (K,1,3) = Inv2Pi*SQRT( ( p%KBE(K,1,1) + KBECent(K,1,1) )/ MBE(K,1,1) ) ! Natural Blade 1-edge frequency w/ centrifugal stiffening and tip mass effects
! Calculate the generalized damping of the blades:
DO I = 1,NumBF ! Loop through flap DOFs
DO L = 1,NumBF ! Loop through flap DOFs
p%CBF(K,I,L) = ( 0.01*p%BldFDamp(K,L) )*p%KBF(K,I,L)/( Pi*p%FreqBF(K,L,1) )
ENDDO ! L - Flap DOFs
ENDDO ! I - Flap DOFs
p%CBE (K,1,1) = ( 0.01*p%BldEDamp(K,1) )*p%KBE(K,1,1)/( Pi*p%FreqBE(K,1,1) )
! Calculate the 2nd derivatives of the twisted shape functions at the blade root:
Shape = SHP( 0.0_ReKi, p%BldFlexL, p%BldFl1Sh(:,K), 2, ErrStat, ErrMsg )
p%TwistedSF(K,1,1,0,2) = Shape*p%CThetaS(K,0) ! 2nd deriv. of Phi1(0) for blade K
p%TwistedSF(K,2,1,0,2) = -Shape*p%SThetaS(K,0) ! 2nd deriv. of Psi1(0) for blade K
Shape = SHP( 0.0_ReKi, p%BldFlexL, p%BldFl2Sh(:,K), 2, ErrStat, ErrMsg )
p%TwistedSF(K,1,2,0,2) = Shape*p%CThetaS(K,0) ! 2nd deriv. of Phi2(0) for blade K
p%TwistedSF(K,2,2,0,2) = -Shape*p%SThetaS(K,0) ! 2nd deriv. of Psi2(0) for blade K
Shape = SHP( 0.0_ReKi, p%BldFlexL, p%BldEdgSh(:,K), 2, ErrStat, ErrMsg )
p%TwistedSF(K,1,3,0,2) = Shape*p%SThetaS(K,0) ! 2nd deriv. of Phi3(0) for blade K
p%TwistedSF(K,2,3,0,2) = Shape*p%CThetaS(K,0) ! 2nd deriv. of Psi3(0) for blade K
! Calculate the 2nd derivatives of the twisted shape functions at the tip:
Shape = SHP( 1.0_ReKi, p%BldFlexL, p%BldFl1Sh(:,K), 2, ErrStat, ErrMsg )
p%TwistedSF(K,1,1,p%TipNode,2) = Shape*p%CThetaS(K,p%TipNode) ! 2nd deriv. of Phi1(p%TipNode) for blade K
p%TwistedSF(K,2,1,p%TipNode,2) = -Shape*p%SThetaS(K,p%TipNode) ! 2nd deriv. of Psi1(p%TipNode) for blade K
Shape = SHP( 1.0_ReKi, p%BldFlexL, p%BldFl2Sh(:,K), 2, ErrStat, ErrMsg )
p%TwistedSF(K,1,2,p%TipNode,2) = Shape*p%CThetaS(K,p%TipNode) ! 2nd deriv. of Phi2(p%TipNode) for blade K
p%TwistedSF(K,2,2,p%TipNode,2) = -Shape*p%SThetaS(K,p%TipNode) ! 2nd deriv. of Psi2(p%TipNode) for blade K
Shape = SHP( 1.0_ReKi, p%BldFlexL, p%BldEdgSh(:,K), 2, ErrStat, ErrMsg )
p%TwistedSF(K,1,3,p%TipNode,2) = Shape*p%SThetaS(K,p%TipNode) ! 2nd deriv. of Phi3(p%TipNode) for blade K
p%TwistedSF(K,2,3,p%TipNode,2) = Shape*p%CThetaS(K,p%TipNode) ! 2nd deriv. of Psi3(p%TipNode) for blade K
! Integrate to find the 1st and zeroeth derivatives of the twisted shape functions
! at the tip:
DO I = 1,2 ! Loop through Phi and Psi
DO L = 1,3 ! Loop through all blade DOFs
p%TwistedSF(K,I,L,p%TipNode,1) = p%TwistedSF(K,I,L,p%BldNodes,1) + TwstdSFOld(I,L,1)
p%TwistedSF(K,I,L,p%TipNode,0) = p%TwistedSF(K,I,L,p%BldNodes,0) + TwstdSFOld(I,L,0)
ENDDO ! L - All blade DOFs
ENDDO ! I - Phi and Psi
! the 1st and zeroeth derivatives of the twisted shape functions at the blade root:
p%TwistedSF(K,:,:,0,1) = 0.0_ReKi
p%TwistedSF(K,:,:,0,0) = 0.0_ReKi
p%AxRedBld( K,:,:,0 ) = 0.0_ReKi
! Integrate to find the blade axial reduction shape functions at the tip:
DO I = 1,3 ! Loop through all blade DOFs
DO L = 1,3 ! Loop through all blade DOFs
p%AxRedBld(K,I,L,p%TipNode) = p%AxRedBld(K,I,L,p%BldNodes) + AxRdBldOld(I,L)
ENDDO ! L - All blade DOFs
ENDDO ! I - All blade DOFs
END IF ! p%BD4Blades
ENDDO ! K - Blades
! Calculate the tower-top mass:
p%TwrTpMass = p%RotMass + p%RFrlMass + p%BoomMass + p%TFinMass + p%NacMass + p%YawBrMass
DO J = p%TwrNodes,1,-1 ! Loop through the tower nodes / elements in reverse
! Calculate the mass of the current element
p%TElmntMass(J) = p%MassT(J)*abs(p%DHNodes(J)) ! Mass of tower element J
! Integrate to find the tower mass which will be output in .fsm
p%TwrMass = p%TwrMass + p%TElmntMass(J)
! Integrate to find TMssAbvNd:
TMssAbvNd (J) = 0.5*p%TElmntMass(J)
IF ( J == p%TwrNodes ) THEN ! Uppermost tower element
! Add the TwrTpMass effects:
TMssAbvNd(J) = TMssAbvNd(J) + p%TwrTpMass
ELSE ! All other tower elements
! Add to TMssAbvNd(J) the effects from the (not yet used) portion of element J+1
TMssAbvNd(J) = 0.5*p%TElmntMass(J+1) + TMssAbvNd(J) + TMssAbvNd(J+1)
ENDIF
ENDDO ! J - Tower nodes / elements in reverse
! Initialize the generalized tower masses using tower-top mass effects:
DO I = 1,2 ! Loop through all tower modes in a single direction
MTFA(I,I) = p%TwrTpMass
MTSS(I,I) = p%TwrTpMass
ENDDO ! I - All tower modes in a single direction
! set values for tower base (note that we haven't corrctly defined the values for (:,0,2) in the arrays below):
p%TwrFASF( :,0,0:1) = 0.0_ReKi
p%TwrSSSF( :,0,0:1) = 0.0_ReKi
p%AxRedTFA(:,:,0) = 0.0_ReKi
p%AxRedTSS(:,:,0) = 0.0_ReKi
DO J = 1,p%TwrNodes ! Loop through the tower nodes / elements
! Calculate the tower shape functions (all derivatives):
p%TwrFASF(1,J,2) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM1Sh(:), 2, ErrStat, ErrMsg )
p%TwrFASF(2,J,2) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM2Sh(:), 2, ErrStat, ErrMsg )
p%TwrFASF(1,J,1) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM1Sh(:), 1, ErrStat, ErrMsg )
p%TwrFASF(2,J,1) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM2Sh(:), 1, ErrStat, ErrMsg )
p%TwrFASF(1,J,0) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM1Sh(:), 0, ErrStat, ErrMsg )
p%TwrFASF(2,J,0) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM2Sh(:), 0, ErrStat, ErrMsg )
p%TwrSSSF(1,J,2) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM1Sh(:), 2, ErrStat, ErrMsg )
p%TwrSSSF(2,J,2) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM2Sh(:), 2, ErrStat, ErrMsg )
p%TwrSSSF(1,J,1) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM1Sh(:), 1, ErrStat, ErrMsg )
p%TwrSSSF(2,J,1) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM2Sh(:), 1, ErrStat, ErrMsg )
p%TwrSSSF(1,J,0) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM1Sh(:), 0, ErrStat, ErrMsg )
p%TwrSSSF(2,J,0) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM2Sh(:), 0, ErrStat, ErrMsg )
! Integrate to find the generalized mass of the tower (including tower-top mass effects).
! Ignore the cross-correlation terms of MTFA (i.e. MTFA(i,j) where i /= j) and MTSS
! since these terms will never be used.
DO I = 1,2 ! Loop through all tower DOFs in one direction
MTFA (I,I) = MTFA (I,I) + p%TElmntMass(J)*p%TwrFASF(I,J,0)**2
MTSS (I,I) = MTSS (I,I) + p%TElmntMass(J)*p%TwrSSSF(I,J,0)**2
ENDDO ! I - through all tower DOFs in one direction
! Integrate to find the generalized stiffness of the tower (not including gravitational
! effects).
ElStffFA = p%StiffTFA(J)*abs(p%DHNodes(J)) ! Fore-aft stiffness of tower element J
ElStffSS = p%StiffTSS(J)*abs(p%DHNodes(J)) ! Side-to-side stiffness of tower element J
DO I = 1,2 ! Loop through all tower DOFs in one direction
DO L = 1,2 ! Loop through all tower DOFs in one direction
p%KTFA (I,L) = p%KTFA (I,L) + ElStffFA *p%TwrFASF(I,J,2)*p%TwrFASF(L,J,2)
p%KTSS (I,L) = p%KTSS (I,L) + ElStffSS *p%TwrSSSF(I,J,2)*p%TwrSSSF(L,J,2)
ENDDO ! L - All tower DOFs in one direction
ENDDO ! I - through all tower DOFs in one direction
! Integrate to find the gravitational-term of the generalized stiffness of the tower.
! Ignore the cross-correlation terms of KTFAGrav (i.e. KTFAGrav(i,j) where i /= j)
! and KTSSGrav since these terms will never be used.
ElmntStff = -TMssAbvNd(J)*abs(p%DHNodes(J))*p%Gravity ! Gravitational stiffness of tower element J
DO I = 1,2 ! Loop through all tower DOFs in one direction
KTFAGrav(I,I) = KTFAGrav(I,I) + ElmntStff*p%TwrFASF(I,J,1)**2
KTSSGrav(I,I) = KTSSGrav(I,I) + ElmntStff*p%TwrSSSF(I,J,1)**2
ENDDO
! Integrate to find the tower axial reduction shape functions:
DO I = 1,2 ! Loop through all tower DOFs in one direction
DO L = 1,2 ! Loop through all tower DOFs in one direction
AxRdTFA (I,L) = 0.5*p%DHNodes(J)*p%TwrFASF(I,J,1)*p%TwrFASF(L,J,1)
AxRdTSS (I,L) = 0.5*p%DHNodes(J)*p%TwrSSSF(I,J,1)*p%TwrSSSF(L,J,1)
p%AxRedTFA(I,L,J) = AxRdTFA(I,L)
p%AxRedTSS(I,L,J) = AxRdTSS(I,L)
ENDDO ! L - All tower DOFs in one direction
ENDDO
IF ( J /= 1 ) THEN ! All but the lowermost tower element
! Add the effects from the (not yet used) portion of element J-1
DO I = 1,2 ! Loop through all tower DOFs in one direction
DO L = 1,2 ! Loop through all tower DOFs in one direction
p%AxRedTFA(I,L,J) = p%AxRedTFA(I,L,J) + p%AxRedTFA(I,L,J-1)+ AxRdTFAOld(I,L)
p%AxRedTSS(I,L,J) = p%AxRedTSS(I,L,J) + p%AxRedTSS(I,L,J-1)+ AxRdTSSOld(I,L)
ENDDO ! L - All tower DOFs in one direction
ENDDO
ENDIF
! Store the AxRdTFA and AxRdTSS terms of the current element (these will be used for the next element)
AxRdTFAOld = AxRdTFA
AxRdTSSOld = AxRdTSS
ENDDO ! J - Tower nodes / elements
! Apply the modal stiffness tuners of the tower to KTFA() and KTSS():
DO I = 1,2 ! Loop through all tower DOFs in one direction
DO L = 1,2 ! Loop through all tower DOFs in one direction
p%KTFA(I,L) = SQRT( InputFileData%FAStTunr(I)*InputFileData%FAStTunr(L) )*p%KTFA(I,L)
p%KTSS(I,L) = SQRT( InputFileData%SSStTunr(I)*InputFileData%SSStTunr(L) )*p%KTSS(I,L)
ENDDO ! L - All tower DOFs in one direction
ENDDO ! I - through all tower DOFs in one direction
! Calculate the tower natural frequencies:
DO I = 1,2 ! Loop through all tower DOFs in one direction
if ( EqualRealNos(( MTFA(I,I) - p%TwrTpMass ), 0.0_ReKi) ) then
p%FreqTFA(I,1) = NaN ! Avoid creating a divide by zero signal, but set p%FreqTFA(I,1) = NaN
else
p%FreqTFA(I,1) = Inv2Pi*SQRT( p%KTFA(I,I)/( MTFA(I,I) - p%TwrTpMass ) ) ! Natural tower I-fore-aft frequency w/o gravitational destiffening nor tower-top mass effects
end if
if ( EqualRealNos(( MTSS(I,I) - p%TwrTpMass ), 0.0_ReKi) ) then
p%FreqTSS(I,1) = NaN ! Avoid creating a divide by zero signal, but set p%FreqTFS(I,1) = NaN
else
p%FreqTSS(I,1) = Inv2Pi*SQRT( p%KTSS(I,I)/( MTSS(I,I) - p%TwrTpMass ) ) ! Natural tower I-side-to-side frequency w/o gravitational destiffening nor tower-top mass effects
end if
p%FreqTFA(I,2) = Inv2Pi*SQRT( ( p%KTFA(I,I) + KTFAGrav(I,I) )/MTFA(I,I) ) ! Natural tower I-fore-aft frequency w/ gravitational destiffening and tower-top mass effects
p%FreqTSS(I,2) = Inv2Pi*SQRT( ( p%KTSS(I,I) + KTSSGrav(I,I) )/MTSS(I,I) ) ! Natural tower I-side-to-side frequency w/ gravitational destiffening and tower-top mass effects
ENDDO ! I - All tower DOFs in one direction
! Calculate the generalized damping of the tower:
DO I = 1,2 ! Loop through all tower DOFs in one direction
DO L = 1,2 ! Loop through all tower DOFs in one direction
p%CTFA(I,L) = ( 0.01*InputFileData%TwrFADmp(L) )*p%KTFA(I,L)/( Pi*p%FreqTFA(L,1) )
p%CTSS(I,L) = ( 0.01*InputFileData%TwrSSDmp(L) )*p%KTSS(I,L)/( Pi*p%FreqTSS(L,1) )
ENDDO ! L - All tower DOFs in one direction
ENDDO ! I - All tower DOFs in one direction
! Calculate the tower shape functions (all derivatives) at the tower-top:
p%TwrFASF(1,p%TTopNode,2) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM1Sh(:), 2, ErrStat, ErrMsg )
p%TwrFASF(2,p%TTopNode,2) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM2Sh(:), 2, ErrStat, ErrMsg )
p%TwrFASF(1,p%TTopNode,1) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM1Sh(:), 1, ErrStat, ErrMsg )
p%TwrFASF(2,p%TTopNode,1) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM2Sh(:), 1, ErrStat, ErrMsg )
p%TwrFASF(1,p%TTopNode,0) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM1Sh(:), 0, ErrStat, ErrMsg )
p%TwrFASF(2,p%TTopNode,0) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM2Sh(:), 0, ErrStat, ErrMsg )
p%TwrSSSF(1,p%TTopNode,2) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM1Sh(:), 2, ErrStat, ErrMsg )
p%TwrSSSF(2,p%TTopNode,2) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM2Sh(:), 2, ErrStat, ErrMsg )
p%TwrSSSF(1,p%TTopNode,1) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM1Sh(:), 1, ErrStat, ErrMsg )
p%TwrSSSF(2,p%TTopNode,1) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM2Sh(:), 1, ErrStat, ErrMsg )
p%TwrSSSF(1,p%TTopNode,0) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM1Sh(:), 0, ErrStat, ErrMsg )
p%TwrSSSF(2,p%TTopNode,0) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM2Sh(:), 0, ErrStat, ErrMsg )
! Integrate to find the tower axial reduction shape functions at the tower-top:
DO I = 1,2 ! Loop through all tower DOFs in one direction
DO L = 1,2 ! Loop through all tower DOFs in one direction
p%AxRedTFA(I,L,p%TTopNode) = p%AxRedTFA(I,L,p%TwrNodes)+ AxRdTFAOld(I,L)
p%AxRedTSS(I,L,p%TTopNode) = p%AxRedTSS(I,L,p%TwrNodes)+ AxRdTSSOld(I,L)
ENDDO ! L - All tower DOFs in one direction
ENDDO
! Calculate the turbine mass:
p%TurbMass = p%TwrTpMass + p%TwrMass
RETURN
END SUBROUTINE Coeff