Skip to content

AD: add new projection method and BEM methods #1283

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
ebranlard merged 21 commits into from
Oct 20, 2022
Merged

AD: add new projection method and BEM methods #1283

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
6f72c4f
AD: adding optional projMod to AeroDyn input file (hardcoded to 0)
ebranlard Sep 27, 2022
fc39798
AD: renaming withoutSweepPitchTwist to orientationAnnulus (ReKi not R…
ebranlard Sep 28, 2022
f3af436
AD: adding and allocating new variables
ebranlard Sep 28, 2022
59db763
AD: starting projmod 2
ebranlard Sep 28, 2022
2adf279
AD: starting to change BEM for BEMT_Mod
ebranlard Sep 28, 2022
1507f22
AD: BEMT using uTimes for UpdateStates
ebranlard Sep 28, 2022
183953a
AD: adding function for BEMTUncoupled with BEMT_Mod switch
ebranlard Sep 28, 2022
0b0c6c6
AD: introducting new induction calculation method for BEMTMod=2
ebranlard Sep 28, 2022
f890c9a
AD: automatically computing cant angle when AeroProjMod=2
ebranlard Sep 28, 2022
6d0be23
AD: introducing parameters for AeroProjMod and BEM_Mod
ebranlard Sep 30, 2022
258231f
Move quartic.f90 from AD to NWTC_Library
bjonkman Oct 10, 2022
b4e4157
AD: removing EulerDCM parameter
bjonkman Oct 10, 2022
94336e9
Merge branch 'dev' into f/ad-newproj
bjonkman Oct 10, 2022
ae4c0af
update r-test
bjonkman Oct 10, 2022
fb983f7
AD: update spacing
bjonkman Oct 10, 2022
7f4c309
fix typo in manualRegressionTest.py
bjonkman Oct 11, 2022
52e184c
BEMT: update variable name for easier searching
bjonkman Oct 11, 2022
5f1814d
Explicitly initialize AeroProjMod in OpenFAST driver
bjonkman Oct 11, 2022
896f717
Merge pull request #19 from bjonkman/f/ad-newproj
ebranlard Oct 11, 2022
4f0d335
AD: output M instead of Mz
ebranlard Oct 12, 2022
784e4b6
Merge branch 'dev' into f/ad-newproj
ebranlard Oct 20, 2022
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
520 changes: 396 additions & 124 deletions modules/aerodyn/src/AeroDyn.f90
View file
Open in desktop

Large diffs are not rendered by default.

27 changes: 13 additions & 14 deletions modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@ MODULE AeroDyn_AllBldNdOuts_IO

! Parameters related to output length (number of characters allowed in the output data headers):

! INTEGER(IntKi), PARAMETER :: OutStrLenM1_Msuffix = ChanLen - 6 ! The NREL allowed channel name length is usually 20. We are making these of the form B#N##namesuffix


! ===================================================================================================
Expand Down Expand Up @@ -253,7 +252,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx
CASE ( BldNd_VUndx )
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,p%NumBlNds
y%WriteOutput( OutIdx ) = dot_product( m%WithoutSweepPitchTwist(1,:,IdxNode,IdxBlade), u%InflowOnBlade(:,IdxNode,IdxBlade) )
y%WriteOutput( OutIdx ) = dot_product( m%orientationAnnulus(1,:,IdxNode,IdxBlade), u%InflowOnBlade(:,IdxNode,IdxBlade) )
OutIdx = OutIdx + 1
END DO
END DO
Expand All @@ -262,15 +261,15 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx
CASE ( BldNd_VUndy )
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,p%NumBlNds
y%WriteOutput( OutIdx ) = dot_product( m%WithoutSweepPitchTwist(2,:,IdxNode,IdxBlade), u%InflowOnBlade(:,IdxNode,IdxBlade) )
y%WriteOutput( OutIdx ) = dot_product( m%orientationAnnulus(2,:,IdxNode,IdxBlade), u%InflowOnBlade(:,IdxNode,IdxBlade) )
OutIdx = OutIdx + 1
END DO
END DO

CASE ( BldNd_VUndz )
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,p%NumBlNds
y%WriteOutput( OutIdx ) = dot_product( m%WithoutSweepPitchTwist(3,:,IdxNode,IdxBlade), u%InflowOnBlade(:,IdxNode,IdxBlade) )
y%WriteOutput( OutIdx ) = dot_product( m%orientationAnnulus(3,:,IdxNode,IdxBlade), u%InflowOnBlade(:,IdxNode,IdxBlade) )
OutIdx = OutIdx + 1
END DO
END DO
Expand Down Expand Up @@ -309,23 +308,23 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx
CASE ( BldNd_VDisx )
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,p%NumBlNds
y%WriteOutput( OutIdx ) = dot_product( m%WithoutSweepPitchTwist(1,:,IdxNode,IdxBlade), m%DisturbedInflow(:,IdxNode,IdxBlade) )
y%WriteOutput( OutIdx ) = dot_product( m%orientationAnnulus(1,:,IdxNode,IdxBlade), m%DisturbedInflow(:,IdxNode,IdxBlade) )
OutIdx = OutIdx + 1
END DO
END DO

CASE ( BldNd_VDisy )
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,p%NumBlNds
y%WriteOutput( OutIdx ) = dot_product( m%WithoutSweepPitchTwist(2,:,IdxNode,IdxBlade), m%DisturbedInflow(:,IdxNode,IdxBlade) )
y%WriteOutput( OutIdx ) = dot_product( m%orientationAnnulus(2,:,IdxNode,IdxBlade), m%DisturbedInflow(:,IdxNode,IdxBlade) )
OutIdx = OutIdx + 1
END DO
END DO

CASE ( BldNd_VDisz )
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,p%NumBlNds
y%WriteOutput( OutIdx ) = dot_product( m%WithoutSweepPitchTwist(3,:,IdxNode,IdxBlade), m%DisturbedInflow(:,IdxNode,IdxBlade) )
y%WriteOutput( OutIdx ) = dot_product( m%orientationAnnulus(3,:,IdxNode,IdxBlade), m%DisturbedInflow(:,IdxNode,IdxBlade) )
OutIdx = OutIdx + 1
END DO
END DO
Expand All @@ -335,23 +334,23 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx
CASE ( BldNd_STVx )
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,p%NumBlNds
y%WriteOutput( OutIdx ) = dot_product( m%WithoutSweepPitchTwist(1,:,IdxNode,IdxBlade), u%BladeMotion(IdxBlade)%TranslationVel(:,IdxNode) )
y%WriteOutput( OutIdx ) = dot_product( m%orientationAnnulus(1,:,IdxNode,IdxBlade), u%BladeMotion(IdxBlade)%TranslationVel(:,IdxNode) )
OutIdx = OutIdx + 1
END DO
END DO

CASE ( BldNd_STVy )
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,p%NumBlNds
y%WriteOutput( OutIdx ) = dot_product( m%WithoutSweepPitchTwist(2,:,IdxNode,IdxBlade), u%BladeMotion(IdxBlade)%TranslationVel(:,IdxNode) )
y%WriteOutput( OutIdx ) = dot_product( m%orientationAnnulus(2,:,IdxNode,IdxBlade), u%BladeMotion(IdxBlade)%TranslationVel(:,IdxNode) )
OutIdx = OutIdx + 1
END DO
END DO

CASE ( BldNd_STVz )
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,p%NumBlNds
y%WriteOutput( OutIdx ) = dot_product( m%WithoutSweepPitchTwist(3,:,IdxNode,IdxBlade), u%BladeMotion(IdxBlade)%TranslationVel(:,IdxNode) )
y%WriteOutput( OutIdx ) = dot_product( m%orientationAnnulus(3,:,IdxNode,IdxBlade), u%BladeMotion(IdxBlade)%TranslationVel(:,IdxNode) )
OutIdx = OutIdx + 1
END DO
END DO
Expand Down Expand Up @@ -1166,7 +1165,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,u%BladeMotion(IdxBlade)%NNodes
Vind_s = (/ -m%BEMT_u(Indx)%Vx(IdxNode,IdxBlade)*m%BEMT_y%axInduction(IdxNode,IdxBlade), m%BEMT_u(Indx)%Vy(IdxNode,IdxBlade)*m%BEMT_y%tanInduction(IdxNode,IdxBlade), 0.0_ReKi /)
Vind_g = matmul(Vind_s, m%WithoutSweepPitchTwist(:,:,IdxNode,IdxBlade))
Vind_g = matmul(Vind_s, m%orientationAnnulus(:,:,IdxNode,IdxBlade))
y%WriteOutput( OutIdx ) = dot_product(M_pg(1,1:3,IdxBlade), Vind_g(1:3) ) ! Uihn, hub normal
OutIdx = OutIdx + 1
ENDDO
Expand All @@ -1187,7 +1186,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,u%BladeMotion(IdxBlade)%NNodes
Vind_s = (/ -m%BEMT_u(Indx)%Vx(IdxNode,IdxBlade)*m%BEMT_y%axInduction(IdxNode,IdxBlade), m%BEMT_u(Indx)%Vy(IdxNode,IdxBlade)*m%BEMT_y%tanInduction(IdxNode,IdxBlade), 0.0_ReKi /)
Vind_g = matmul(Vind_s, m%WithoutSweepPitchTwist(:,:,IdxNode,IdxBlade))
Vind_g = matmul(Vind_s, m%orientationAnnulus(:,:,IdxNode,IdxBlade))
y%WriteOutput( OutIdx ) = dot_product(M_pg(2,1:3,IdxBlade), Vind_g(1:3) ) ! Uiht, hub tangential
OutIdx = OutIdx + 1
ENDDO
Expand All @@ -1208,7 +1207,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx
DO IdxBlade=1,p%BldNd_BladesOut
DO IdxNode=1,u%BladeMotion(IdxBlade)%NNodes
Vind_s = (/ -m%BEMT_u(Indx)%Vx(IdxNode,IdxBlade)*m%BEMT_y%axInduction(IdxNode,IdxBlade), m%BEMT_u(Indx)%Vy(IdxNode,IdxBlade)*m%BEMT_y%tanInduction(IdxNode,IdxBlade), 0.0_ReKi /)
Vind_g = matmul(Vind_s, m%WithoutSweepPitchTwist(:,:,IdxNode,IdxBlade))
Vind_g = matmul(Vind_s, m%orientationAnnulus(:,:,IdxNode,IdxBlade))
y%WriteOutput( OutIdx ) = dot_product(M_pg(3,1:3,IdxBlade), Vind_g(1:3) ) ! Uihr, hub radial
OutIdx = OutIdx + 1
ENDDO
Expand Down Expand Up @@ -1304,7 +1303,7 @@ SUBROUTINE AllBldNdOuts_SetParameters( InputFileData, p, p_AD, ErrStat, ErrMsg )


! Check if the requested blades exist
IF ( (InputFileData%BldNd_BladesOut < 0_IntKi) ) then
IF ( (InputFileData%BldNd_BladesOut < 0_IntKi) ) THEN
p%BldNd_BladesOut = 0_IntKi
ELSE IF ((InputFileData%BldNd_BladesOut > p%NumBlades) ) THEN
CALL SetErrStat( ErrID_Warn, " Number of blades to output data at all blade nodes (BldNd_BladesOut) must be no more than the total number of blades, "//TRIM(Num2LStr(p%NumBlades))//".", ErrStat, ErrMsg, RoutineName)
Expand Down
1 change: 1 addition & 0 deletions modules/aerodyn/src/AeroDyn_Driver_Registry.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,7 @@ typedef ^ ^ TwrData twr
typedef ^ ^ IntKi numBlades - - - "" -
typedef ^ ^ logical basicHAWTFormat - - - "If true simply input HubRad/Pitch/Overhang/Cone, otherwise all turbine inputs" -
typedef ^ ^ logical hasTower - - - "" -
typedef ^ ^ IntKi projMod - - - "If true simply input HubRad/Pitch/Overhang/Cone, otherwise all turbine inputs" -
typedef ^ ^ logical HAWTprojection - - - "" -
typedef ^ ^ IntKi motionType - - - "" -
typedef ^ ^ ReKi motion :: - - "" "-"
Expand Down
19 changes: 16 additions & 3 deletions modules/aerodyn/src/AeroDyn_Driver_Subs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -432,11 +432,17 @@ subroutine Init_AeroDyn(iCase, dvr, AD, dt, InitOutData, errStat, errMsg)
InitInData%rotors(iWT)%numBlades = wt%numBlades
call AllocAry(InitInData%rotors(iWT)%BladeRootPosition, 3, wt%numBlades, 'BladeRootPosition', errStat2, ErrMsg2 ); if (Failed()) return
call AllocAry(InitInData%rotors(iWT)%BladeRootOrientation, 3, 3, wt%numBlades, 'BladeRootOrientation', errStat2, ErrMsg2 ); if (Failed()) return
if (wt%HAWTprojection) then
InitInData%rotors(iWT)%AeroProjMod = 0 ! default, with WithoutSweepPitchTwist
if (wt%projMod==-1)then
call WrScr('>>> Using HAWTprojection to determine projMod')
if (wt%HAWTprojection) then
InitInData%rotors(iWT)%AeroProjMod = APM_BEM_NoSweepPitchTwist ! default, with WithoutSweepPitchTwist
else
InitInData%rotors(iWT)%AeroProjMod = APM_LiftingLine
endif
else
InitInData%rotors(iWT)%AeroProjMod = 1
InitInData%rotors(iWT)%AeroProjMod = wt%projMod
endif
call WrScr('>>> Using projection method '//trim(num2lstr(InitInData%rotors(iWT)%AeroProjMod)))
InitInData%rotors(iWT)%HubPosition = wt%hub%ptMesh%Position(:,1)
InitInData%rotors(iWT)%HubOrientation = wt%hub%ptMesh%RefOrientation(:,:,1)
InitInData%rotors(iWT)%NacellePosition = wt%nac%ptMesh%Position(:,1)
Expand Down Expand Up @@ -1350,6 +1356,13 @@ subroutine Dvr_ReadInputFile(fileName, dvr, errStat, errMsg )
wt => dvr%WT(iWT)
sWT = '('//trim(num2lstr(iWT))//')'
call ParseCom(FileInfo_In, CurLine, Line, errStat2, errMsg2, unEc); if(Failed()) return
! Temporary hack, look if ProjMod is present on the line
!call ParseVar(FileInfo_In, CurLine, 'ProjMod'//sWT , wt%projMod , errStat2, errMsg2, unEc); if(Failed()) return
call ParseVar(FileInfo_In, CurLine, 'ProjMod'//sWT , wt%projMod , errStat2, errMsg2, unEc);
if (errStat2==ErrID_Fatal) then
call WrScr('>>> ProjMod is not present in AeroDyn driver input file.')
wt%projMod = -1
endif
call ParseVar(FileInfo_In, CurLine, 'BasicHAWTFormat'//sWT , wt%basicHAWTFormat , errStat2, errMsg2, unEc); if(Failed()) return

! Basic init
Expand Down
7 changes: 7 additions & 0 deletions modules/aerodyn/src/AeroDyn_Driver_Types.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -190,6 +190,7 @@ MODULE AeroDyn_Driver_Types
INTEGER(IntKi) :: numBlades !< [-]
LOGICAL :: basicHAWTFormat !< If true simply input HubRad/Pitch/Overhang/Cone, otherwise all turbine inputs [-]
LOGICAL :: hasTower !< [-]
INTEGER(IntKi) :: projMod !< If true simply input HubRad/Pitch/Overhang/Cone, otherwise all turbine inputs [-]
LOGICAL :: HAWTprojection !< [-]
INTEGER(IntKi) :: motionType !< [-]
REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: motion !< [-]
Expand Down Expand Up @@ -5821,6 +5822,7 @@ SUBROUTINE AD_Dvr_CopyWTData( SrcWTDataData, DstWTDataData, CtrlCode, ErrStat, E
DstWTDataData%numBlades = SrcWTDataData%numBlades
DstWTDataData%basicHAWTFormat = SrcWTDataData%basicHAWTFormat
DstWTDataData%hasTower = SrcWTDataData%hasTower
DstWTDataData%projMod = SrcWTDataData%projMod
DstWTDataData%HAWTprojection = SrcWTDataData%HAWTprojection
DstWTDataData%motionType = SrcWTDataData%motionType
IF (ALLOCATED(SrcWTDataData%motion)) THEN
Expand Down Expand Up @@ -6070,6 +6072,7 @@ SUBROUTINE AD_Dvr_PackWTData( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMs
Int_BufSz = Int_BufSz + 1 ! numBlades
Int_BufSz = Int_BufSz + 1 ! basicHAWTFormat
Int_BufSz = Int_BufSz + 1 ! hasTower
Int_BufSz = Int_BufSz + 1 ! projMod
Int_BufSz = Int_BufSz + 1 ! HAWTprojection
Int_BufSz = Int_BufSz + 1 ! motionType
Int_BufSz = Int_BufSz + 1 ! motion allocated yes/no
Expand Down Expand Up @@ -6337,6 +6340,8 @@ SUBROUTINE AD_Dvr_PackWTData( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMs
Int_Xferred = Int_Xferred + 1
IntKiBuf(Int_Xferred) = TRANSFER(InData%hasTower, IntKiBuf(1))
Int_Xferred = Int_Xferred + 1
IntKiBuf(Int_Xferred) = InData%projMod
Int_Xferred = Int_Xferred + 1
IntKiBuf(Int_Xferred) = TRANSFER(InData%HAWTprojection, IntKiBuf(1))
Int_Xferred = Int_Xferred + 1
IntKiBuf(Int_Xferred) = InData%motionType
Expand Down Expand Up @@ -6732,6 +6737,8 @@ SUBROUTINE AD_Dvr_UnPackWTData( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, Er
Int_Xferred = Int_Xferred + 1
OutData%hasTower = TRANSFER(IntKiBuf(Int_Xferred), OutData%hasTower)
Int_Xferred = Int_Xferred + 1
OutData%projMod = IntKiBuf(Int_Xferred)
Int_Xferred = Int_Xferred + 1
OutData%HAWTprojection = TRANSFER(IntKiBuf(Int_Xferred), OutData%HAWTprojection)
Int_Xferred = Int_Xferred + 1
OutData%motionType = IntKiBuf(Int_Xferred)
Expand Down
Loading