From f20aaeb5cd28cc89b311b237aaddc6305fce896c Mon Sep 17 00:00:00 2001 From: Emmanuel Branlard Date: Thu, 30 Nov 2023 18:28:26 -0700 Subject: [PATCH] AD: adding skewModRedistr to BEMT (partial) --- modules/aerodyn/src/AeroDyn.f90 | 3 ++- modules/aerodyn/src/AeroDyn_IO.f90 | 12 ++++++------ modules/aerodyn/src/AeroDyn_Registry.txt | 2 +- modules/aerodyn/src/AeroDyn_Types.f90 | 2 +- modules/aerodyn/src/BEMT.f90 | 12 +++++++++--- modules/aerodyn/src/BEMTUncoupled.f90 | 8 ++++++-- modules/aerodyn/src/BEMT_Registry.txt | 8 +++++--- modules/aerodyn/src/BEMT_Types.f90 | 16 +++++++++++++--- 8 files changed, 43 insertions(+), 20 deletions(-) diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index 185c7f769b..019a0fd666 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -3847,7 +3847,7 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg ) if ( InputFileData%IndToler < 0.0 .or. EqualRealNos(InputFileData%IndToler, 0.0_ReKi) ) & call SetErrStat( ErrID_Fatal, 'IndToler must be greater than 0.', ErrStat, ErrMsg, RoutineName ) - if ( InputFileData%Skew_Mod /= Skew_Mod_Orthogonal .and. InputFileData%Skew_Mod /= Skew_Mod_None .and. InputFileData%Skew_Mod /= Skew_Mod_Glauert) & + if ( InputFileData%Skew_Mod /= Skew_Mod_Orthogonal .and. InputFileData%Skew_Mod /= Skew_Mod_None .and. InputFileData%Skew_Mod /= Skew_Mod_Active) & call SetErrStat( ErrID_Fatal, 'Skew_Mod must be -1, 0, or 1.', ErrStat, ErrMsg, RoutineName ) if ( InputFileData%SectAvg) then @@ -4295,6 +4295,7 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x InitInp%airDens = InputFileData%AirDens InitInp%kinVisc = InputFileData%KinVisc InitInp%skewWakeMod = InputFileData%Skew_Mod + InitInp%skewRedistrMod = InputFileData%SkewRedistrMod InitInp%yawCorrFactor = InputFileData%SkewModFactor InitInp%aTol = InputFileData%IndToler InitInp%useTipLoss = InputFileData%TipLoss diff --git a/modules/aerodyn/src/AeroDyn_IO.f90 b/modules/aerodyn/src/AeroDyn_IO.f90 index 4c9eb931eb..a18afcf2aa 100644 --- a/modules/aerodyn/src/AeroDyn_IO.f90 +++ b/modules/aerodyn/src/AeroDyn_IO.f90 @@ -771,8 +771,8 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade ! Skew_Mod- Select skew model {0: No skew model at all, -1:Throw away non-normal component for linearization, 1: Glauert skew model, } call ParseVar( FileInfo_In, CurLine, "Skew_Mod", InputFileData%Skew_Mod, ErrStat2, ErrMsg2, UnEc ) if (newInputAbsent('Skew_Mod', CurLine, errStat2, errMsg2)) then - call WrScr(' Setting Skew_Mod to 1 (skew active, Gluaert) as the input is absent (typical behavior).') - InputFileData%Skew_Mod = Skew_Mod_Glauert + call WrScr(' Setting Skew_Mod to 1 (skew active) as the input is absent (typical behavior).') + InputFileData%Skew_Mod = Skew_Mod_Active else if (skewModProvided) then call LegacyAbort('Cannot have both Skew_Mod and SkewMod in the input file'); return @@ -1065,7 +1065,7 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade else if (InputFileData%SkewMod==1) then InputFileData%Skew_Mod = Skew_Mod_None else if (InputFileData%SkewMod==2) then - InputFileData%Skew_Mod = Skew_Mod_Glauert + InputFileData%Skew_Mod = Skew_Mod_Active else call LegacyAbort('Legacy option SkewMod is not 0, 1,2 which is not supported.') endif @@ -1589,12 +1589,12 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg ) Msg = 'orthogonal' case (Skew_Mod_None) Msg = 'no correction' - case (Skew_Mod_Glauert) - Msg = 'Glauert/Pitt/Peters' + case (Skew_Mod_Active) + Msg = 'active' case default Msg = 'unknown' end select - WRITE (UnSu,Ec_IntFrmt) InputFileData%Skew_Mod, 'Skew_Mod', 'Type of skewed-wake correction model: '//TRIM(Msg) + WRITE (UnSu,Ec_IntFrmt) InputFileData%Skew_Mod, 'Skew_Mod', 'Skewed-wake correction model: '//TRIM(Msg) ! TipLoss diff --git a/modules/aerodyn/src/AeroDyn_Registry.txt b/modules/aerodyn/src/AeroDyn_Registry.txt index 88ab542493..86be5a9aeb 100644 --- a/modules/aerodyn/src/AeroDyn_Registry.txt +++ b/modules/aerodyn/src/AeroDyn_Registry.txt @@ -182,7 +182,7 @@ typedef ^ AD_InputFile ReKi SpdSound - - - "Speed of sound" m/s typedef ^ AD_InputFile IntKi SkewMod - - - "LEGACY - Skew Mod" - typedef ^ AD_InputFile IntKi Skew_Mod - - - "Select skew model {0=No skew model at all, -1=Throw away non-normal component for linearization, 1=Glauert skew model}" - typedef ^ AD_InputFile Logical SkewMomCorr - - - "Turn the skew momentum correction on or off [used only when SkewMod=1]" - -typedef ^ AD_InputFile IntKi SkewRedistrMod - - - "Type of skewed-wake correction model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1]" - +typedef ^ AD_InputFile IntKi SkewRedistrMod - - - "Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1]" - typedef ^ AD_InputFile ReKi SkewModFactor - - - "Constant used in Pitt/Peters skewed wake model (default is 15*pi/32)" - typedef ^ AD_InputFile LOGICAL TipLoss - - - "Use the Prandtl tip-loss model? [unused when WakeMod=0]" flag typedef ^ AD_InputFile LOGICAL HubLoss - - - "Use the Prandtl hub-loss model? [unused when WakeMod=0]" flag diff --git a/modules/aerodyn/src/AeroDyn_Types.f90 b/modules/aerodyn/src/AeroDyn_Types.f90 index 276b0f7dbf..a1c0bb6eb5 100644 --- a/modules/aerodyn/src/AeroDyn_Types.f90 +++ b/modules/aerodyn/src/AeroDyn_Types.f90 @@ -211,7 +211,7 @@ MODULE AeroDyn_Types INTEGER(IntKi) :: SkewMod = 0_IntKi !< LEGACY - Skew Mod [-] INTEGER(IntKi) :: Skew_Mod = 0_IntKi !< Select skew model {0=No skew model at all, -1=Throw away non-normal component for linearization, 1=Glauert skew model} [-] LOGICAL :: SkewMomCorr = .false. !< Turn the skew momentum correction on or off [used only when SkewMod=1] [-] - INTEGER(IntKi) :: SkewRedistrMod = 0_IntKi !< Type of skewed-wake correction model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1] [-] + INTEGER(IntKi) :: SkewRedistrMod = 0_IntKi !< Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1] [-] REAL(ReKi) :: SkewModFactor = 0.0_ReKi !< Constant used in Pitt/Peters skewed wake model (default is 15*pi/32) [-] LOGICAL :: TipLoss = .false. !< Use the Prandtl tip-loss model? [unused when WakeMod=0] [flag] LOGICAL :: HubLoss = .false. !< Use the Prandtl hub-loss model? [unused when WakeMod=0] [flag] diff --git a/modules/aerodyn/src/BEMT.f90 b/modules/aerodyn/src/BEMT.f90 index 7cf1c80e97..568845bf36 100644 --- a/modules/aerodyn/src/BEMT.f90 +++ b/modules/aerodyn/src/BEMT.f90 @@ -187,7 +187,6 @@ subroutine BEMT_SetParameters( InitInp, p, errStat, errMsg ) p%numBlades = InitInp%numBlades p%UA_Flag = InitInp%UA_Flag p%DBEMT_Mod = InitInp%DBEMT_Mod - p%MomentumCorr = InitInp%MomentumCorr p%BEM_Mod = InitInp%BEM_Mod !call WrScr('>>>> BEM_Mod '//trim(num2lstr(p%BEM_Mod))) if (.not.(ANY( p%BEM_Mod == (/BEMMod_2D, BEMMod_3D/)))) then @@ -250,6 +249,13 @@ subroutine BEMT_SetParameters( InitInp, p, errStat, errMsg ) p%airDens = InitInp%airDens p%kinVisc = InitInp%kinVisc p%skewWakeMod = InitInp%skewWakeMod + if (p%skewWakeMod==Skew_Mod_Active) then + p%skewRedistrMod = InitInp%skewRedistrMod + p%MomentumCorr = InitInp%MomentumCorr + else + p%skewRedistrMod = SkewRedistrMod_None + p%MomentumCorr = .false. + endif p%yawCorrFactor = InitInp%yawCorrFactor p%useTipLoss = InitInp%useTipLoss p%useHubLoss = InitInp%useHubLoss @@ -1590,7 +1596,7 @@ subroutine ApplySkewedWakeCorrection_AllNodes(p, u, m, x, phi, OtherState, axInd !............................................ ! Apply skewed wake correction to the axial induction (y%axInduction) !............................................ - if ( p%skewWakeMod == Skew_Mod_Glauert ) then + if ( p%skewWakeMod == Skew_Mod_Active ) then if (p%BEM_Mod==BEMMod_2D) then ! do nothing else @@ -1602,7 +1608,7 @@ subroutine ApplySkewedWakeCorrection_AllNodes(p, u, m, x, phi, OtherState, axInd do i = 1,p%numBladeNodes ! Loop through the blade nodes / elements if ( .not. p%FixedInductions(i,j) ) then F = getHubTipLossCorrection(p%BEM_Mod, p%useHubLoss, p%useTipLoss, p%hubLossConst(i,j), p%tipLossConst(i,j), phi(i,j), u%cantAngle(i,j) ) - call ApplySkewedWakeCorrection( p%BEM_Mod, p%skewWakeMod, p%yawCorrFactor, F, u%psi_s(j), u%psiSkewOffset, u%chi0, u%rlocal(i,j)/m%Rtip(j), axInduction(i,j), chi(i,j), m%FirstWarn_Skew ) + call ApplySkewedWakeCorrection( p%BEM_Mod, p%skewRedistrMod, p%yawCorrFactor, F, u%psi_s(j), u%psiSkewOffset, u%chi0, u%rlocal(i,j)/m%Rtip(j), axInduction(i,j), chi(i,j), m%FirstWarn_Skew ) end if ! .not. p%FixedInductions (special case for tip and/or hub loss) enddo ! I - Blade nodes / elements enddo ! J - All blades diff --git a/modules/aerodyn/src/BEMTUncoupled.f90 b/modules/aerodyn/src/BEMTUncoupled.f90 index 42238a607d..82213fa86b 100644 --- a/modules/aerodyn/src/BEMTUncoupled.f90 +++ b/modules/aerodyn/src/BEMTUncoupled.f90 @@ -449,10 +449,10 @@ real(ReKi) function BEMTU_InductionWithResidual(p, u, i, j, phi, AFInfo, IsValid end function BEMTU_InductionWithResidual !----------------------------------------------------------------------------------------- -subroutine ApplySkewedWakeCorrection(BEM_Mod, SkewMod, yawCorrFactor, F, azimuth, azimuthOffset, chi0, tipRatio, a, chi, FirstWarn ) +subroutine ApplySkewedWakeCorrection(BEM_Mod, SkewRedistrMod, yawCorrFactor, F, azimuth, azimuthOffset, chi0, tipRatio, a, chi, FirstWarn ) integer(IntKi), intent(in ) :: BEM_Mod - integer(IntKi), intent(in ) :: SkewMod + integer(IntKi), intent(in ) :: SkewRedistrMod real(ReKi), intent(in ) :: yawCorrFactor ! set to 15*pi/32 previously; now allowed to be input (to better match data) real(ReKi), intent(in ) :: F ! tip/hub loss factor real(ReKi), intent(in ) :: azimuth @@ -467,6 +467,10 @@ subroutine ApplySkewedWakeCorrection(BEM_Mod, SkewMod, yawCorrFactor, F, azimuth real(ReKi) :: yawCorr real(ReKi) :: yawCorr_tan ! magnitude of the tan(chi/2) correction term (with possible limits) + if (SkewRedistrMod==SkewRedistrMod_None) then + print*,'>>> SkewRedistrMod is None, Shouldnt happen' + STOP + endif ! Skewed wake correction if(BEM_Mod==BEMMod_2D) then diff --git a/modules/aerodyn/src/BEMT_Registry.txt b/modules/aerodyn/src/BEMT_Registry.txt index ca226ffd01..fc6540c2fe 100644 --- a/modules/aerodyn/src/BEMT_Registry.txt +++ b/modules/aerodyn/src/BEMT_Registry.txt @@ -17,8 +17,8 @@ usefrom UnsteadyAero_Registry.txt usefrom DBEMT_Registry.txt param BEMT/BEMT - INTEGER Skew_Mod_Orthogonal - -1 - "Inflow orthogonal to rotor [-]" - -param BEMT/BEMT - INTEGER Skew_Mod_None - 0 - "No skew model (previously called uncoupled)" - -param BEMT/BEMT - INTEGER Skew_Mod_Glauert - 1 - "Pitt/Peters/Glauert skew model (should be 1)" - +param BEMT/BEMT - INTEGER Skew_Mod_None - 0 - "No skew model" - +param BEMT/BEMT - INTEGER Skew_Mod_Active - 1 - "Skew model active" - param BEMT/BEMT - INTEGER Skew_Mod_PittPeters_Cont - 4 - "Pitt/Peters continuous formulation" - param BEMT/BEMT - INTEGER SkewRedistrMod_None - 0 - "No redistribution" - @@ -39,7 +39,8 @@ typedef BEMT/BEMT InitInputType ReKi typedef ^ ^ INTEGER numBlades - - - "Number of blades" - typedef ^ ^ ReKi airDens - - - "Air density" kg/m^3 typedef ^ ^ ReKi kinVisc - - - "Kinematic air viscosity" m^2/s -typedef ^ ^ INTEGER skewWakeMod - - - "Type of skewed-wake correction model [switch] {1=uncoupled, 2=Pitt/Peters, 3=coupled}" - +typedef ^ ^ INTEGER skewWakeMod - - - "Type of skewed-wake model [switch] {0=None, 1=Glauert}" - +typedef ^ ^ INTEGER skewRedistrMod - - - "Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1]" - typedef ^ ^ ReKi aTol - - - "Tolerance for the induction solution" - typedef ^ ^ LOGICAL useTipLoss - - - "Use the Prandtl tip-loss model? [flag]" - typedef ^ ^ LOGICAL useHubLoss - - - "Use the Prandtl hub-loss model? [flag]" - @@ -144,6 +145,7 @@ typedef ^ ^ INTEGER typedef ^ ^ ReKi airDens - - - "Air density" kg/m^3 typedef ^ ^ ReKi kinVisc - - - "Kinematic air viscosity" m^2/s typedef ^ ^ INTEGER skewWakeMod - - - "Type of skewed-wake correction model [switch] {0=None, 1=Glauert/Pitt/Peters}" - +typedef ^ ^ INTEGER skewRedistrMod - - - "Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1]" - typedef ^ ^ ReKi aTol - - - "Tolerance for the induction solution" - typedef ^ ^ LOGICAL useTipLoss - - - "Use the Prandtl tip-loss model? [flag]" - typedef ^ ^ LOGICAL useHubLoss - - - "Use the Prandtl hub-loss model? [flag]" - diff --git a/modules/aerodyn/src/BEMT_Types.f90 b/modules/aerodyn/src/BEMT_Types.f90 index 9202214daf..893263d342 100644 --- a/modules/aerodyn/src/BEMT_Types.f90 +++ b/modules/aerodyn/src/BEMT_Types.f90 @@ -37,8 +37,8 @@ MODULE BEMT_Types USE NWTC_Library IMPLICIT NONE INTEGER(IntKi), PUBLIC, PARAMETER :: Skew_Mod_Orthogonal = -1 ! Inflow orthogonal to rotor [-] [-] - INTEGER(IntKi), PUBLIC, PARAMETER :: Skew_Mod_None = 0 ! No skew model (previously called uncoupled) [-] - INTEGER(IntKi), PUBLIC, PARAMETER :: Skew_Mod_Glauert = 1 ! Pitt/Peters/Glauert skew model (should be 1) [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: Skew_Mod_None = 0 ! No skew model [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: Skew_Mod_Active = 1 ! Skew model active [-] INTEGER(IntKi), PUBLIC, PARAMETER :: Skew_Mod_PittPeters_Cont = 4 ! Pitt/Peters continuous formulation [-] INTEGER(IntKi), PUBLIC, PARAMETER :: SkewRedistrMod_None = 0 ! No redistribution [-] INTEGER(IntKi), PUBLIC, PARAMETER :: SkewRedistrMod_PittPeters = 1 ! Pitt/Petesr/Glauert redistribution [-] @@ -52,7 +52,8 @@ MODULE BEMT_Types INTEGER(IntKi) :: numBlades = 0_IntKi !< Number of blades [-] REAL(ReKi) :: airDens = 0.0_ReKi !< Air density [kg/m^3] REAL(ReKi) :: kinVisc = 0.0_ReKi !< Kinematic air viscosity [m^2/s] - INTEGER(IntKi) :: skewWakeMod = 0_IntKi !< Type of skewed-wake correction model [switch] {1=uncoupled, 2=Pitt/Peters, 3=coupled} [-] + INTEGER(IntKi) :: skewWakeMod = 0_IntKi !< Type of skewed-wake model [switch] {0=None, 1=Glauert} [-] + INTEGER(IntKi) :: skewRedistrMod = 0_IntKi !< Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1] [-] REAL(ReKi) :: aTol = 0.0_ReKi !< Tolerance for the induction solution [-] LOGICAL :: useTipLoss = .false. !< Use the Prandtl tip-loss model? [flag] [-] LOGICAL :: useHubLoss = .false. !< Use the Prandtl hub-loss model? [flag] [-] @@ -154,6 +155,7 @@ MODULE BEMT_Types REAL(ReKi) :: airDens = 0.0_ReKi !< Air density [kg/m^3] REAL(ReKi) :: kinVisc = 0.0_ReKi !< Kinematic air viscosity [m^2/s] INTEGER(IntKi) :: skewWakeMod = 0_IntKi !< Type of skewed-wake correction model [switch] {0=None, 1=Glauert/Pitt/Peters} [-] + INTEGER(IntKi) :: skewRedistrMod = 0_IntKi !< Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1] [-] REAL(ReKi) :: aTol = 0.0_ReKi !< Tolerance for the induction solution [-] LOGICAL :: useTipLoss = .false. !< Use the Prandtl tip-loss model? [flag] [-] LOGICAL :: useHubLoss = .false. !< Use the Prandtl hub-loss model? [flag] [-] @@ -259,6 +261,7 @@ subroutine BEMT_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode, ErrS DstInitInputData%airDens = SrcInitInputData%airDens DstInitInputData%kinVisc = SrcInitInputData%kinVisc DstInitInputData%skewWakeMod = SrcInitInputData%skewWakeMod + DstInitInputData%skewRedistrMod = SrcInitInputData%skewRedistrMod DstInitInputData%aTol = SrcInitInputData%aTol DstInitInputData%useTipLoss = SrcInitInputData%useTipLoss DstInitInputData%useHubLoss = SrcInitInputData%useHubLoss @@ -428,6 +431,7 @@ subroutine BEMT_PackInitInput(Buf, Indata) call RegPack(Buf, InData%airDens) call RegPack(Buf, InData%kinVisc) call RegPack(Buf, InData%skewWakeMod) + call RegPack(Buf, InData%skewRedistrMod) call RegPack(Buf, InData%aTol) call RegPack(Buf, InData%useTipLoss) call RegPack(Buf, InData%useHubLoss) @@ -522,6 +526,8 @@ subroutine BEMT_UnPackInitInput(Buf, OutData) if (RegCheckErr(Buf, RoutineName)) return call RegUnpack(Buf, OutData%skewWakeMod) if (RegCheckErr(Buf, RoutineName)) return + call RegUnpack(Buf, OutData%skewRedistrMod) + if (RegCheckErr(Buf, RoutineName)) return call RegUnpack(Buf, OutData%aTol) if (RegCheckErr(Buf, RoutineName)) return call RegUnpack(Buf, OutData%useTipLoss) @@ -1578,6 +1584,7 @@ subroutine BEMT_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%airDens = SrcParamData%airDens DstParamData%kinVisc = SrcParamData%kinVisc DstParamData%skewWakeMod = SrcParamData%skewWakeMod + DstParamData%skewRedistrMod = SrcParamData%skewRedistrMod DstParamData%aTol = SrcParamData%aTol DstParamData%useTipLoss = SrcParamData%useTipLoss DstParamData%useHubLoss = SrcParamData%useHubLoss @@ -1726,6 +1733,7 @@ subroutine BEMT_PackParam(Buf, Indata) call RegPack(Buf, InData%airDens) call RegPack(Buf, InData%kinVisc) call RegPack(Buf, InData%skewWakeMod) + call RegPack(Buf, InData%skewRedistrMod) call RegPack(Buf, InData%aTol) call RegPack(Buf, InData%useTipLoss) call RegPack(Buf, InData%useHubLoss) @@ -1810,6 +1818,8 @@ subroutine BEMT_UnPackParam(Buf, OutData) if (RegCheckErr(Buf, RoutineName)) return call RegUnpack(Buf, OutData%skewWakeMod) if (RegCheckErr(Buf, RoutineName)) return + call RegUnpack(Buf, OutData%skewRedistrMod) + if (RegCheckErr(Buf, RoutineName)) return call RegUnpack(Buf, OutData%aTol) if (RegCheckErr(Buf, RoutineName)) return call RegUnpack(Buf, OutData%useTipLoss)