From ba7c86a7e91a9e44b958ab61c98c81242aa1093d Mon Sep 17 00:00:00 2001 From: Emmanuel Branlard Date: Fri, 13 Nov 2020 17:57:56 -0700 Subject: [PATCH 1/2] Allowing one blade in ElastoDyn --- modules/elastodyn/src/ElastoDyn.f90 | 79 ++++++++++++++++---------- modules/elastodyn/src/ElastoDyn_IO.f90 | 4 +- 2 files changed, 51 insertions(+), 32 deletions(-) diff --git a/modules/elastodyn/src/ElastoDyn.f90 b/modules/elastodyn/src/ElastoDyn.f90 index 2e0c2814e2..545d3a4f58 100644 --- a/modules/elastodyn/src/ElastoDyn.f90 +++ b/modules/elastodyn/src/ElastoDyn.f90 @@ -1241,11 +1241,8 @@ SUBROUTINE ED_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg ) m%AllOuts( QD_Y ) = x%QDT( DOF_Y ) m%AllOuts( QD2_B1E1 ) = m%QD2T( DOF_BE(1,1) ) - m%AllOuts( QD2_B2E1 ) = m%QD2T( DOF_BE(2,1) ) m%AllOuts( QD2_B1F1 ) = m%QD2T( DOF_BF(1,1) ) - m%AllOuts( QD2_B2F1 ) = m%QD2T( DOF_BF(2,1) ) m%AllOuts( QD2_B1F2 ) = m%QD2T( DOF_BF(1,2) ) - m%AllOuts( QD2_B2F2 ) = m%QD2T( DOF_BF(2,2) ) m%AllOuts( QD2_DrTr ) = m%QD2T( DOF_DrTr ) m%AllOuts( QD2_GeAz ) = m%QD2T( DOF_GeAz ) m%AllOuts( QD2_RFrl ) = m%QD2T( DOF_RFrl ) @@ -1271,6 +1268,10 @@ SUBROUTINE ED_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg ) m%AllOuts( QD_B2E1 ) = x%QDT( DOF_BE(2,1) ) m%AllOuts( QD_B2F1 ) = x%QDT( DOF_BF(2,1) ) m%AllOuts( QD_B2F2 ) = x%QDT( DOF_BF(2,2) ) + + m%AllOuts( QD2_B2E1 ) = m%QD2T( DOF_BE(2,1) ) + m%AllOuts( QD2_B2F1 ) = m%QD2T( DOF_BF(2,1) ) + m%AllOuts( QD2_B2F2 ) = m%QD2T( DOF_BF(2,2) ) IF ( p%NumBl > 2 ) THEN m%AllOuts( Q_B3E1 ) = x%QT( DOF_BE(3,1) ) @@ -2092,8 +2093,9 @@ SUBROUTINE Init_DOFparameters( InputFileData, p, ErrStat, ErrMsg ) ErrStat = ErrID_None ErrMsg = '' - - IF ( p%NumBl == 2 ) THEN + IF ( p%NumBl == 1 ) THEN + p%NDOF = 18 + ELSEIF ( p%NumBl == 2 ) THEN p%NDOF = 22 ELSE p%NDOF = ED_MaxDOFs @@ -4478,21 +4480,34 @@ SUBROUTINE SetOutParam(OutList, p, ErrStat, ErrMsg ) END DO - IF ( p%NumBl < 3_IntKi ) THEN - InvalidOutput(PtchPMzc3) = .TRUE. - - InvalidOutput( Q_B3E1) = .TRUE. - InvalidOutput( Q_B3F1) = .TRUE. - InvalidOutput( Q_B3F2) = .TRUE. - - InvalidOutput( QD_B3E1) = .TRUE. - InvalidOutput( QD_B3F1) = .TRUE. - InvalidOutput( QD_B3F2) = .TRUE. - InvalidOutput( QD2_B3E1) = .TRUE. - InvalidOutput( QD2_B3F1) = .TRUE. - InvalidOutput( QD2_B3F2) = .TRUE. - ELSE IF ( p%NumBl > 2_IntKi ) THEN + ! Invalid outputs based on number of blades + IF ( p%NumBl < 3 ) THEN + InvalidOutput(PtchPMzc3) = .TRUE. + InvalidOutput( Q_B3E1) = .TRUE. + InvalidOutput( Q_B3F1) = .TRUE. + InvalidOutput( Q_B3F2) = .TRUE. + InvalidOutput( QD_B3E1) = .TRUE. + InvalidOutput( QD_B3F1) = .TRUE. + InvalidOutput( QD_B3F2) = .TRUE. + InvalidOutput( QD2_B3E1) = .TRUE. + InvalidOutput( QD2_B3F1) = .TRUE. + InvalidOutput( QD2_B3F2) = .TRUE. + ENDIF + IF ( p%NumBl < 2 ) THEN + InvalidOutput(PtchPMzc2) = .TRUE. + InvalidOutput( Q_B2E1) = .TRUE. + InvalidOutput( Q_B2F1) = .TRUE. + InvalidOutput( Q_B2F2) = .TRUE. + InvalidOutput( QD_B2E1) = .TRUE. + InvalidOutput( QD_B2F1) = .TRUE. + InvalidOutput( QD_B2F2) = .TRUE. + InvalidOutput( QD2_B2E1) = .TRUE. + InvalidOutput( QD2_B2F1) = .TRUE. + InvalidOutput( QD2_B2F2) = .TRUE. + ENDIF + ! 1-bladed or 3-bladed, no teeter + IF ( p%NumBl /= 2 ) THEN InvalidOutput( TeetPya) = .TRUE. InvalidOutput( TeetVya) = .TRUE. InvalidOutput( TeetAya) = .TRUE. @@ -8225,11 +8240,13 @@ SUBROUTINE FillAugMat( p, x, CoordSys, u, HSSBrTrq, RtHSdat, AugMat ) ! Initialize the portions of the mass matrix below the diagonal associated ! with the teeter and pure blade DOFs using the partial loads at the teeter pin; only do this if necessary: - IF ( ( p%NumBl == 2 ) .AND. ( p%DOF_Flag(DOF_Teet) ) ) THEN - DO L = 1,p%DOFs%NPSBE(K) ! Loop through all active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the blade - AugMat(DOF_Teet,p%DOFs%PSBE(K,L)) = -DOT_PRODUCT( RtHSdat%PAngVelEH(DOF_Teet,0,:), & - RtHSdat%PMomLPRot(:,p%DOFs%PSBE(K,L)) ) ! [C(q,t)]B - ENDDO ! L - All active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the blade + IF ( ( p%NumBl == 2 ) ) THEN + IF ( p%DOF_Flag(DOF_Teet) ) THEN ! NOTE: two "ifs" since DOF_Teet might be out of bound + DO L = 1,p%DOFs%NPSBE(K) ! Loop through all active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the blade + AugMat(DOF_Teet,p%DOFs%PSBE(K,L)) = -DOT_PRODUCT( RtHSdat%PAngVelEH(DOF_Teet,0,:), & + RtHSdat%PMomLPRot(:,p%DOFs%PSBE(K,L)) ) ! [C(q,t)]B + ENDDO ! L - All active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the blade + ENDIF ENDIF @@ -8512,11 +8529,13 @@ SUBROUTINE FillAugMat( p, x, CoordSys, u, HSSBrTrq, RtHSdat, AugMat ) + RtHSdat%TFrlMom ! + {-f(qd,q,t)}SpringTF + {-f(qd,q,t)}DampTF ENDIF - IF ( ( p%NumBl == 2 ) .AND. ( p%DOF_Flag(DOF_Teet) ) ) THEN - ! The teeter DOF does not affect any DOF index larger than DOF_Teet. Therefore, there is no need to perform the loop: DO I = Diag(DOF_Teet),NActvDOF + IF ( ( p%NumBl == 2 ) ) THEN + IF ( p%DOF_Flag(DOF_Teet) ) THEN ! NOTE: two "ifs" since DOF_Teet may be out of bound + ! The teeter DOF does not affect any DOF index larger than DOF_Teet. Therefore, there is no need to perform the loop: DO I = Diag(DOF_Teet),NActvDOF AugMat(DOF_Teet, DOF_Teet) = -DOT_PRODUCT( RtHSdat%PAngVelEH(DOF_Teet,0,:), RtHSdat%PMomLPRot(:,DOF_Teet) ) ! [C(q,t)]H + [C(q,t)]B AugMat(DOF_Teet, p%NAug) = DOT_PRODUCT( RtHSdat%PAngVelEH(DOF_Teet,0,:), RtHSdat%MomLPRott ) & ! {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB + RtHSdat%TeetMom ! + {-f(qd,q,t)}SpringTeet + {-f(qd,q,t)}DampTeet + ENDIF ENDIF !.................................................................................................................................. ! So far, we have only filled in the portions of the mass matrix on and below the diagonal. Because the mass matrix is symmetric @@ -10086,15 +10105,15 @@ SUBROUTINE ED_PrintSum( p, OtherState, GenerateAdamsModel, ErrStat, ErrMsg ) ELSE RotorType = 'Upwind,' ENDIF - IF ( p%NumBl == 2 ) THEN - RotorType = TRIM(RotorType)//' two-bladed rotor' - IF ( p%DOF_Flag(DOF_Teet) ) THEN + RotorType = TRIM(RotorType)//' '//trim(Num2LStr(p%NumBl))//'-bladed rotor' + IF ( p%NumBl==2 ) THEN + IF ( p%DOF_Flag(DOF_Teet) ) THEN ! NOTE: two "ifs" required since DOF_Teet might be out of bound RotorType = TRIM(RotorType)//' with teetering hub.' ELSE RotorType = TRIM(RotorType)//' with rigid hub.' ENDIF ELSE - RotorType = TRIM(RotorType)//' three-bladed rotor with rigid hub.' + RotorType = TRIM(RotorType)//' with rigid hub.' ENDIF WRITE (UnSu,'(A)') ' '//TRIM(RotorType) diff --git a/modules/elastodyn/src/ElastoDyn_IO.f90 b/modules/elastodyn/src/ElastoDyn_IO.f90 index d097d0801b..033ace77a2 100644 --- a/modules/elastodyn/src/ElastoDyn_IO.f90 +++ b/modules/elastodyn/src/ElastoDyn_IO.f90 @@ -5045,8 +5045,8 @@ SUBROUTINE ValidatePrimaryData( InputFileData, BD4Blades, Linearize, ErrStat, Er CALL SetErrStat( ErrID_Fatal, 'NumBl must be 1, 2, or 3 for BeamDyn simulations.',ErrStat,ErrMsg,RoutineName) END IF ELSE - IF ( ( InputFileData%NumBl < 2 ) .OR. ( InputFileData%NumBl > MaxBl ) ) THEN - CALL SetErrStat( ErrID_Fatal, 'NumBl must be either 2 or 3.',ErrStat,ErrMsg,RoutineName) + IF ( ( InputFileData%NumBl < 1 ) .OR. ( InputFileData%NumBl > MaxBl ) ) THEN + CALL SetErrStat( ErrID_Fatal, 'NumBl must be 1, 2, or 3.',ErrStat,ErrMsg,RoutineName) END IF END IF From fff196658a69cd8776bffa5a5b0592e793510190 Mon Sep 17 00:00:00 2001 From: Emmanuel Branlard Date: Wed, 18 Nov 2020 18:00:52 -0700 Subject: [PATCH 2/2] Additional support for servodyn with 1 blade --- modules/servodyn/src/BladedInterface.f90 | 9 ++++--- modules/servodyn/src/ServoDyn.f90 | 30 ++++++++++++++---------- 2 files changed, 24 insertions(+), 15 deletions(-) diff --git a/modules/servodyn/src/BladedInterface.f90 b/modules/servodyn/src/BladedInterface.f90 index 8fd332cdec..bfe8c41c72 100644 --- a/modules/servodyn/src/BladedInterface.f90 +++ b/modules/servodyn/src/BladedInterface.f90 @@ -882,8 +882,11 @@ SUBROUTINE Fill_CONTROL_vars( t, u, p, ErrMsgSz, dll_data ) dll_data%ShaftBrakeStatusBinaryFlag = 0 ! no brakes deployed dll_data%HSSBrDeployed = .false. - dll_data%PrevBlPitch(1:p%NumBl) = p%BlPitchInit - dll_data%BlPitchCom(1:p%NumBl) = p%BlPitchInit + dll_data%PrevBlPitch(:) = 0.0_ReKi ! Harcoded to size 3 + dll_data%BlPitchCom(:) = 0.0_ReKi ! Harcoded to size 3 + dll_data%BlAirfoilCom(:)= 0.0_ReKi ! Harcoded to size 3 + dll_data%PrevBlPitch(1:p%NumBl) = p%BlPitchInit(1:p%NumBl) + dll_data%BlPitchCom(1:p%NumBl) = p%BlPitchInit(1:p%NumBl) end if call Fill_avrSWAP( t, u, p, ErrMsgSz, dll_data ) ! we'll set the avrSWAP variable, for the legacy version of the DLL, too. @@ -972,7 +975,7 @@ SUBROUTINE Retrieve_avrSWAP( p, dll_data, ErrStat, ErrMsg ) ELSE !IF ( p%Ptch_Cntrl == GH_DISCON_PITCH_CONTROL_COLLECTIVE ) THEN ! Collective pitch control !> * Record 45: Demanded pitch angle (Collective pitch) (rad) - dll_data%BlPitchCom = dll_data%avrSWAP(45) ! Demanded pitch angle (Collective pitch) (rad) + dll_data%BlPitchCom(:) = dll_data%avrSWAP(45) ! Demanded pitch angle (Collective pitch) (rad) !> * Record 46, demanded pitch rate (Collective pitch), is ingored since record 10 is set to 0 by ServoDyn indicating pitch position actuator diff --git a/modules/servodyn/src/ServoDyn.f90 b/modules/servodyn/src/ServoDyn.f90 index ac0fd59c83..ca9206fa81 100644 --- a/modules/servodyn/src/ServoDyn.f90 +++ b/modules/servodyn/src/ServoDyn.f90 @@ -340,7 +340,7 @@ SUBROUTINE SrvD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO END IF - u%BlPitch = p%BlPitchInit + u%BlPitch = p%BlPitchInit(1:p%NumBl) u%Yaw = p%YawNeut u%YawRate = 0.0 @@ -351,7 +351,7 @@ SUBROUTINE SrvD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO u%ExternalYawPosCom = p%YawNeut u%ExternalYawRateCom = 0. - u%ExternalBlPitchCom = p%BlPitchInit + u%ExternalBlPitchCom = p%BlPitchInit(1:p%NumBl) u%ExternalGenTrq = 0. u%ExternalElecPwr = 0. u%ExternalHSSBrFrac = 0. @@ -361,11 +361,11 @@ SUBROUTINE SrvD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO u%WindDir = 0. !Inputs for the Bladed Interface: - u%RootMyc = 0. + u%RootMyc(:) = 0. ! Hardcoded to 3 u%YawBrTAxp = 0. u%YawBrTAyp = 0. u%LSSTipPxa = 0. - u%RootMxc = 0. + u%RootMxc(:)= 0. ! Hardcoded to 3 u%LSSTipMxa = 0. u%LSSTipMya = 0. u%LSSTipMza = 0. @@ -535,7 +535,7 @@ SUBROUTINE SrvD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO CALL CheckError( ErrStat2, ErrMsg2 ) IF (ErrStat >= AbortErrLev) RETURN - do i=1,size(SrvD_Indx_Y_BlPitchCom) + do i=1,size(SrvD_Indx_Y_BlPitchCom) ! NOTE: potentially limit to NumBl InitOut%LinNames_y(SrvD_Indx_Y_BlPitchCom(i)) = 'BlPitchCom('//trim(num2lstr(i))//'), rad' InitOut%RotFrame_y(SrvD_Indx_Y_BlPitchCom(i)) = .true. end do @@ -983,7 +983,7 @@ SUBROUTINE DLL_controller_call(t, u, p, x, xd, z, OtherState, m, ErrStat, ErrMsg m%FirstWarn = .FALSE. END IF ELSE - m%dll_data%PrevBlPitch(1:p%NumBl) = m%dll_data%BlPitchCom ! used for linear ramp of delayed signal + m%dll_data%PrevBlPitch(1:p%NumBl) = m%dll_data%BlPitchCom(1:p%NumBl) ! used for linear ramp of delayed signal m%LastTimeCalled = t CALL BladedInterface_CalcOutput( t, u, p, m, ErrStat2, ErrMsg2 ) @@ -1697,8 +1697,12 @@ SUBROUTINE SrvD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_o end if - do i=1,size(SrvD_Indx_Y_BlPitchCom) - y_op(SrvD_Indx_Y_BlPitchCom(i)) = y%BlPitchCom(i) + do i=1,size(SrvD_Indx_Y_BlPitchCom) ! Note: Potentially limit to NumBl + if (i<=p%NumBl) then + y_op(SrvD_Indx_Y_BlPitchCom(i)) = y%BlPitchCom(i) + else + y_op(SrvD_Indx_Y_BlPitchCom(i)) = 0.0_ReKI + endif end do y_op(SrvD_Indx_Y_YawMom) = y%YawMom y_op(SrvD_Indx_Y_GenTrq) = y%GenTrq @@ -2773,9 +2777,9 @@ SUBROUTINE SrvD_SetParameters( InputFileData, p, ErrStat, ErrMsg ) CALL AllocAry( p%PitManRat, p%NumBl, 'PitManRat', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) IF (ErrStat >= AbortErrLev) RETURN - p%TPitManS = InputFileData%TPitManS(1:p%NumBl) - p%BlPitchF = InputFileData%BlPitchF(1:p%NumBl) - p%PitManRat = InputFileData%PitManRat(1:p%NumBl) + p%TPitManS (1:p%NumBl) = InputFileData%TPitManS(1:p%NumBl) + p%BlPitchF (1:p%NumBl) = InputFileData%BlPitchF(1:p%NumBl) + p%PitManRat(1:p%NumBl) = InputFileData%PitManRat(1:p%NumBl) !............................................. ! Set generator and torque control parameters: @@ -3184,7 +3188,7 @@ SUBROUTINE Pitch_CalcOutput( t, u, p, x, xd, z, OtherState, BlPitchCom, ElecPwr, CASE ( ControlMode_EXTERN ) ! User-defined from Simulink or LabVIEW. - BlPitchCom = u%ExternalBlPitchCom ! copy entire array + BlPitchCom = u%ExternalBlPitchCom(1:p%NumBl) CASE ( ControlMode_DLL ) ! User-defined pitch control from Bladed-style DLL @@ -3372,6 +3376,8 @@ SUBROUTINE SetOutParam(OutList, p, ErrStat, ErrMsg ) ! Determine which inputs are not valid InvalidOutput(BlAirFlC3) = ( p%NumBl < 3 ) + InvalidOutput(BlAirFlC2) = ( p%NumBl < 2 ) + InvalidOutput( BlPitchC2) = ( p%NumBl < 2 ) InvalidOutput( BlPitchC3) = ( p%NumBl < 3 ) InvalidOutput( NTMD_XQ) = ( .not. p%CompNTMD ) InvalidOutput( NTMD_XQD) = ( .not. p%CompNTMD )