Skip to content

ElastoDyn: Allow one-blade turbine model #589

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrew-platt merged 3 commits into from
Mar 19, 2021
+71 −44
Merged

ElastoDyn: Allow one-blade turbine model #589

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ba7c86a
Allowing one blade in ElastoDyn
ebranlard Nov 14, 2020
fff1966
Additional support for servodyn with 1 blade
ebranlard Nov 19, 2020
07483af
Merge branch 'dev' into f/oneblade
ebranlard Mar 19, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 49 additions & 30 deletions modules/elastodyn/src/ElastoDyn.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1255,11 +1255,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 )
Expand All @@ -1285,6 +1282,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) )
Expand Down Expand Up @@ -2111,8 +2112,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
Expand Down Expand Up @@ -4452,21 +4454,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.
Expand Down Expand Up @@ -8209,11 +8224,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


Expand Down Expand Up @@ -8496,11 +8513,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
Expand Down Expand Up @@ -10070,15 +10089,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)
Expand Down
4 changes: 2 additions & 2 deletions modules/elastodyn/src/ElastoDyn_IO.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5047,8 +5047,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

Expand Down
9 changes: 6 additions & 3 deletions modules/servodyn/src/BladedInterface.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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

Expand Down
23 changes: 14 additions & 9 deletions modules/servodyn/src/ServoDyn.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -273,7 +273,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
Expand All @@ -284,7 +284,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.
Expand All @@ -294,11 +294,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.
Expand Down Expand Up @@ -430,7 +430,7 @@ SUBROUTINE SrvD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO
CALL AllocAry( InitOut%LinNames_y, SrvD_Indx_Y_WrOutput+p%NumOuts, 'LinNames_y', ErrStat2, ErrMsg2 )
if (Failed()) 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
Expand Down Expand Up @@ -1176,7 +1176,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 )
Expand Down Expand Up @@ -1903,8 +1903,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
Expand Down Expand Up @@ -2263,6 +2267,7 @@ SUBROUTINE SrvD_SetParameters( InputFileData, p, ErrStat, ErrMsg )
CALL AllocAry( p%PitManRat, p%NumBl, 'PitManRat', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName); p%PitManRat=0.0_ReKi
IF (ErrStat >= AbortErrLev) RETURN


p%TPitManS = InputFileData%TPitManS( 1:min(p%NumBl,size(InputFileData%TPitManS)))
p%BlPitchF = InputFileData%BlPitchF( 1:min(p%NumBl,size(InputFileData%BlPitchF)))
p%PitManRat = InputFileData%PitManRat(1:min(p%NumBl,size(InputFileData%PitManRat)))
Expand Down Expand Up @@ -2676,7 +2681,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

Expand Down