```fortran 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 ```