Skip to content

Commit

Permalink
Pull branch (#18)
Browse files Browse the repository at this point in the history 
AeroDyn V15.04- Updated to include cavitation check

- Updated the BEMT_CalcOutputs to include a check for cavitation. Also included some checks for the cavitation input parameters here so that we don't waste time checking these if the user isn't doing a cavitation check.
  • Loading branch information
@robynnemurrayNREL
robynnemurrayNREL authored and ghaymanNREL committed Apr 20, 2017
1 parent 2d5398b commit 055be86
Show file tree
Hide file tree
Showing 10 changed files with 733 additions and 454 deletions.
35 changes: 27 additions & 8 deletions modules-local/aerodyn/src/AeroDyn.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -866,7 +866,8 @@ subroutine SetParameters( InitInp, InputFileData, p, ErrStat, ErrMsg )

p%AirDens = InputFileData%AirDens
p%KinVisc = InputFileData%KinVisc
p%SpdSound = InputFileData%SpdSound



!p%AFI ! set in call to AFI_Init() [called early because it wants to use the same echo file as AD]
!p%BEMT ! set in call to BEMT_Init()
Expand Down Expand Up @@ -1417,7 +1418,11 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg )
if (InputFileData%AirDens <= 0.0) call SetErrStat ( ErrID_Fatal, 'The air density (AirDens) must be greater than zero.', ErrStat, ErrMsg, RoutineName )
if (InputFileData%KinVisc <= 0.0) call SetErrStat ( ErrID_Fatal, 'The kinesmatic viscosity (KinVisc) must be greater than zero.', ErrStat, ErrMsg, RoutineName )
if (InputFileData%SpdSound <= 0.0) call SetErrStat ( ErrID_Fatal, 'The speed of sound (SpdSound) must be greater than zero.', ErrStat, ErrMsg, RoutineName )

if (InputFileData%Pvap <= 0.0) call SetErrStat ( ErrID_Fatal, 'The vapour pressure (Pvap) must be greater than zero.', ErrStat, ErrMsg, RoutineName )
if (InputFileData%Patm <= 0.0) call SetErrStat ( ErrID_Fatal, 'The atmospheric pressure (Patm) must be greater than zero.', ErrStat, ErrMsg, RoutineName )
if (InputFileData%FluidDepth <= 0.0) call SetErrStat ( ErrID_Fatal, 'Fluid depth (FluidDepth) cannot be negative', ErrStat, ErrMsg, RoutineName )



! BEMT inputs
! bjj: these checks should probably go into BEMT where they are used...
Expand All @@ -1439,7 +1444,14 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg )

if (.not. InputFileData%FLookUp ) call SetErrStat( ErrID_Fatal, 'FLookUp must be TRUE for this version.', ErrStat, ErrMsg, RoutineName )
end if


if ( InputFileData%CavitCheck .and. InputFileData%AFAeroMod == 2) then
call SetErrStat( ErrID_Fatal, 'Cannot use unsteady aerodynamics module with a cavitation check', ErrStat, ErrMsg, RoutineName )
end if

if (InputFileData%InCol_Cpmin == 0 .and. InputFileData%CavitCheck) call SetErrStat( ErrID_Fatal, 'InCol_Cpmin must not be 0 to do a cavitation check.', ErrStat, ErrMsg, RoutineName )



! validate the AFI input data because it doesn't appear to be done in AFI
if (InputFileData%NumAFfiles < 1) call SetErrStat( ErrID_Fatal, 'The number of unique airfoil tables (NumAFfiles) must be greater than zero.', ErrStat, ErrMsg, RoutineName )
Expand Down Expand Up @@ -1699,7 +1711,10 @@ SUBROUTINE Init_BEMTmodule( InputFileData, u_AD, u, p, x, xd, z, OtherState, y,
InitInp%numBlades = p%NumBlades

InitInp%airDens = InputFileData%AirDens
InitInp%kinVisc = InputFileData%KinVisc
InitInp%kinVisc = InputFileData%KinVisc
InitInp%Patm = InputFileData%Patm
InitInp%Pvap = InputFileData%Pvap
InitInp%FluidDepth = InputFileData%FluidDepth
InitInp%skewWakeMod = InputFileData%SkewMod
InitInp%aTol = InputFileData%IndToler
InitInp%useTipLoss = InputFileData%TipLoss
Expand Down Expand Up @@ -1748,10 +1763,13 @@ SUBROUTINE Init_BEMTmodule( InputFileData, u_AD, u, p, x, xd, z, OtherState, y,
end do
end do

InitInp%UA_Flag = InputFileData%AFAeroMod == AFAeroMod_BL_unsteady
InitInp%UAMod = InputFileData%UAMod
InitInp%Flookup = InputFileData%Flookup
InitInp%a_s = InputFileData%SpdSound
InitInp%UA_Flag = InputFileData%AFAeroMod == AFAeroMod_BL_unsteady
InitInp%UAMod = InputFileData%UAMod
InitInp%Flookup = InputFileData%Flookup
InitInp%a_s = InputFileData%SpdSound
InitInp%CavitCheck = InputFileData%CavitCheck



if (ErrStat >= AbortErrLev) then
call cleanup()
Expand Down Expand Up @@ -3757,6 +3775,7 @@ FUNCTION CheckBEMTInputPerturbations( p, m ) RESULT(ValidPerturb)

end if


else ! not UseInduction

do k=1,p%NumBlades
Expand Down
2 changes: 1 addition & 1 deletion modules-local/aerodyn/src/AeroDyn_Driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,4 +197,4 @@ subroutine Dvr_End()
end subroutine Dvr_End
!................................
end program AeroDyn_Driver


1,018 changes: 595 additions & 423 deletions modules-local/aerodyn/src/AeroDyn_IO.f90
View file

Large diffs are not rendered by default.

4 changes: 4 additions & 0 deletions modules-local/aerodyn/src/AeroDyn_Registry.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,12 @@ typedef ^ AD_InputFile IntKi TwrPotent - - - "Type tower influence on wind based
typedef ^ AD_InputFile LOGICAL TwrShadow - - - "Calculate tower influence on wind based on downstream tower shadow?" -
typedef ^ AD_InputFile LOGICAL TwrAero - - - "Calculate tower aerodynamic loads?" flag
typedef ^ AD_InputFile Logical FrozenWake - - - "Flag that tells this module it should assume a frozen wake during linearization." -
typedef ^ AD_InputFile Logical CavitCheck - - - "Flag that tells us if we want to check for cavitation" -
typedef ^ AD_InputFile ReKi AirDens - - - "Air density" kg/m^3
typedef ^ AD_InputFile ReKi KinVisc - - - "Kinematic air viscosity" m^2/s
typedef ^ AD_InputFile ReKi Patm - - - "Atmospheric pressure" Pa
typedef ^ AD_InputFile ReKi Pvap - - - "Vapour pressure" Pa
typedef ^ AD_InputFile ReKi FluidDepth - - - "Submerged hub depth" m
typedef ^ AD_InputFile ReKi SpdSound - - - "Speed of sound" m/s
typedef ^ AD_InputFile IntKi SkewMod - - - "Type of skewed-wake correction model {1=uncoupled, 2=Pitt/Peters, 3=coupled} [used only when WakeMod=1]" -
typedef ^ AD_InputFile LOGICAL TipLoss - - - "Use the Prandtl tip-loss model? [used only when WakeMod=1]" flag
Expand Down
9 changes: 8 additions & 1 deletion modules-local/aerodyn/src/AirfoilInfo.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,7 @@ SUBROUTINE AFI_Init ( InitInput, p, ErrStat, ErrMsg, UnEcho )
p%ColCd = 2
p%ColCm = 0 ! These may or may not be used; initialize to zero in case they aren't used
p%ColCpmin = 0 ! These may or may not be used; initialize to zero in case they aren't used

IF ( InitInput%InCol_Cm > 0 ) THEN
p%ColCm = 3
IF ( InitInput%InCol_Cpmin > 0 ) THEN
Expand All @@ -119,8 +120,8 @@ SUBROUTINE AFI_Init ( InitInput, p, ErrStat, ErrMsg, UnEcho )
p%ColCpmin = 3
END IF
NumCoefs = MAX(p%ColCd, p%ColCm,p%ColCpmin) ! number of non-zero coefficient columns



! Process the airfoil files.

ALLOCATE ( p%AFInfo( InitInput%NumAFfiles ), STAT=ErrStat2 )
Expand All @@ -129,6 +130,9 @@ SUBROUTINE AFI_Init ( InitInput, p, ErrStat, ErrMsg, UnEcho )
RETURN
ENDIF

p%AFInfo( :)%ColCpmin=p%ColCpmin
p%AFInfo( :)%ColCm=p%ColCm


DO File=1,InitInput%NumAFfiles

Expand Down Expand Up @@ -496,6 +500,7 @@ SUBROUTINE ReadAFfile ( AFfile, NumCoefs, InCol_Alfa, InCol_Cl, InCol_Cd, InCol_
TYPE (AFInfoType), INTENT(INOUT) :: AFInfo ! The derived type for holding the constant parameters for this airfoil.



! Local declarations.

REAL(ReKi) :: Coords (2) ! An array to hold data from the airfoil-shape table.
Expand Down Expand Up @@ -805,6 +810,8 @@ SUBROUTINE ReadAFfile ( AFfile, NumCoefs, InCol_Alfa, InCol_Cl, InCol_Cd, InCol_
CALL Cleanup()
RETURN
ENDIF



DO Row=1,AFInfo%Table(Table)%NumAlf

Expand Down
2 changes: 2 additions & 0 deletions modules-local/aerodyn/src/AirfoilInfo_Registry.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,8 @@ typedef ^ ^ INTEGER NumCpminAoAkts - - - "The number of angle-of-attack knots fo
typedef ^ ^ INTEGER NumCpminReKts - - - "The number of log(Re) knots for 2D splines of Cpmin" -
typedef ^ ^ INTEGER NumTabs - - - "The number of airfoil tables in the airfoil file" -
typedef ^ ^ AFI_Table_Type Table {:} - - "The tables of airfoil data for given Re and control setting" -
typedef ^ ^ INTEGER ColCpmin - - - "Column number for Cpmin" -
typedef ^ ^ INTEGER ColCm - - - "Column number for Cm" -

# ..... Initialization data .......................................................................................................
# The following derived type stores information that comes from the calling module (say, AeroDyn):
Expand Down
73 changes: 62 additions & 11 deletions modules-local/aerodyn/src/BEMT.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,9 @@ subroutine BEMT_SetParameters( InitInp, p, errStat, errMsg )
p%numBladeNodes = InitInp%numBladeNodes
p%numBlades = InitInp%numBlades
p%UA_Flag = InitInp%UA_Flag
p%CavitCheck = InitInp%CavitCheck



allocate ( p%chord(p%numBladeNodes, p%numBlades), STAT = errStat2 )
if ( errStat2 /= 0 ) then
Expand Down Expand Up @@ -205,10 +208,13 @@ subroutine BEMT_SetParameters( InitInp, p, errStat, errMsg )
end do
end do


!p%DT = InitInp%DT
!p%DT = InitInp%DT
p%airDens = InitInp%airDens
p%kinVisc = InitInp%kinVisc
p%kinVisc = InitInp%kinVisc
p%Patm = InitInp%Patm
p%Pvap = InitInp%Pvap
p%FluidDepth = InitInp%FluidDepth
p%skewWakeMod = InitInp%skewWakeMod
p%useTipLoss = InitInp%useTipLoss
p%useHubLoss = InitInp%useHubLoss
Expand All @@ -219,6 +225,7 @@ subroutine BEMT_SetParameters( InitInp, p, errStat, errMsg )
p%numReIterations = InitInp%numReIterations
p%maxIndIterations = InitInp%maxIndIterations
p%aTol = InitInp%aTol


end subroutine BEMT_SetParameters

Expand Down Expand Up @@ -431,14 +438,15 @@ end subroutine BEMT_AllocInput


!----------------------------------------------------------------------------------------------------------------------------------
subroutine BEMT_AllocOutput( y, p, errStat, errMsg )
subroutine BEMT_AllocOutput( y, p, m, errStat, errMsg )
! This routine is called from BEMT_Init.
!
!
!..................................................................................................................................

type(BEMT_OutputType), intent( out) :: y ! output data
type(BEMT_ParameterType), intent(in ) :: p ! Parameters
type(BEMT_MiscVarType), intent(inout) :: m ! Misc/optimization variables
integer(IntKi), intent( out) :: errStat ! Error status of the operation
character(*), intent( out) :: errMsg ! Error message if ErrStat /= ErrID_None

Expand All @@ -465,6 +473,10 @@ subroutine BEMT_AllocOutput( y, p, errStat, errMsg )
call allocAry( y%Cm, p%numBladeNodes, p%numBlades, 'y%Cm', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
call allocAry( y%Cl, p%numBladeNodes, p%numBlades, 'y%Cl', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
call allocAry( y%Cd, p%numBladeNodes, p%numBlades, 'y%Cd', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
call allocAry( m%Cpmin, p%numBladeNodes, p%numBlades, 'm%Cpmin', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
call allocAry( m%SigmaCavit, p%numBladeNodes, p%numBlades, 'm%SigmaCavit', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
call allocAry( m%SigmaCavitCrit, p%numBladeNodes, p%numBlades, 'm%SigmaCavitCrit', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)


if (ErrStat >= AbortErrLev) RETURN

Expand All @@ -482,7 +494,8 @@ subroutine BEMT_AllocOutput( y, p, errStat, errMsg )
y%tanInduction = 0.0_ReKi
y%AOA = 0.0_ReKi
y%Cl = 0.0_ReKi
y%Cd = 0.0_ReKi
y%Cd = 0.0_ReKi


end subroutine BEMT_AllocOutput

Expand Down Expand Up @@ -705,6 +718,7 @@ subroutine BEMT_Init( InitInp, u, p, x, xd, z, OtherState, AFInfo, y, misc, Inte
write (69,'(A)') ' '

#endif


do j = 1,p%numBlades
do i = 1,p%numBladeNodes ! Loop over blades and nodes
Expand All @@ -726,6 +740,8 @@ subroutine BEMT_Init( InitInp, u, p, x, xd, z, OtherState, AFInfo, y, misc, Inte
call WrScr( 'Warning: Turning off Unsteady Aerodynamics because C_nalpha is 0. BladeNode = '//trim(num2lstr(i))//', Blade = '//trim(num2lstr(j)) )
end if



end do
end do

Expand Down Expand Up @@ -756,7 +772,7 @@ subroutine BEMT_Init( InitInp, u, p, x, xd, z, OtherState, AFInfo, y, misc, Inte
!call BEMT_InitOut(p, InitOut, errStat2, errMsg2)
!call CheckError( errStat2, errMsg2 )

call BEMT_AllocOutput(y, p, errStat2, errMsg2) !u is sent so we can create sibling meshes
call BEMT_AllocOutput(y, p, misc, errStat2, errMsg2) !u is sent so we can create sibling meshes
call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName )
if (errStat >= AbortErrLev) then
call cleanup()
Expand Down Expand Up @@ -1081,8 +1097,8 @@ subroutine BEMT_CalcOutput( t, u, p, x, xd, z, OtherState, AFInfo, y, m, errStat
! Local variables:


real(ReKi) :: Re, fzero
real(ReKi) :: Rtip ! maximum rlocal value for node j over all blades
real(ReKi) :: Re, fzero, theta, Vx, Vy
real(ReKi) :: Rtip, SigmaCavitCrit, SigmaCavit ! maximum rlocal value for node j over all blades

integer(IntKi) :: i ! Generic index
integer(IntKi) :: j ! Loops through nodes / elements
Expand Down Expand Up @@ -1162,6 +1178,12 @@ subroutine BEMT_CalcOutput( t, u, p, x, xd, z, OtherState, AFInfo, y, m, errStat

NodeTxt = '(node '//trim(num2lstr(i))//', blade '//trim(num2lstr(j))//')'

! local velocities and twist angle
Vx = u%Vx(i,j)
Vy = u%Vy(i,j)



! Set the active blade element for UnsteadyAero
m%UA%iBladeNode = i
m%UA%iBlade = j
Expand Down Expand Up @@ -1235,9 +1257,34 @@ subroutine BEMT_CalcOutput( t, u, p, x, xd, z, OtherState, AFInfo, y, m, errStat
else
! TODO: When we start using Re, should we use the uninduced Re since we used uninduced Re to solve for the inductions!? Probably this won't change, instead create a Re loop up above.
call ComputeSteadyAirfoilCoefs( y%AOA(i,j), y%Re(i,j), AFInfo(p%AFindx(i,j)), &
y%Cl(i,j), y%Cd(i,j), y%Cm(i,j), errStat2, errMsg2 )
y%Cl(i,j), y%Cd(i,j), y%Cm(i,j), m%Cpmin(i,j), errStat2, errMsg2 )
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName//trim(NodeTxt))
if (errStat >= AbortErrLev) return




! Calculate the cavitation number for the airfoil at the node in quesiton, and compare to the critical cavitation number based on the vapour pressure and submerged depth
if ( p%CavitCheck ) then
SigmaCavit= -1* m%Cpmin(i,j) ! Cavitation number on blade node j

if ( EqualRealNos( y%Vrel(i,j), 0.0_ReKi ) ) then !if Vrel = 0 in certain cases when Prandtls tip and hub loss factors are used, use the relative verlocity without induction
if ( EqualRealNos( Vx, 0.0_ReKi ) .and. EqualRealNos( Vy, 0.0_ReKi ) ) call SetErrStat( ErrID_Fatal, 'Velocity can not be zero for cavitation check, turn off Prandtls tip loss', ErrStat, ErrMsg, RoutineName )
SigmaCavitCrit= ( ( p%Patm + ( 9.81_ReKi * (p%FluidDepth - ( u%rlocal(i,j))* cos(u%psi(j) )) * p%airDens)) - p%Pvap ) / ( 0.5_ReKi * p%airDens * (sqrt((Vx**2 + Vy**2)))**2) ! Critical value of Sigma, cavitation if we go over this

else
SigmaCavitCrit= ( ( p%Patm + ( 9.81_ReKi * (p%FluidDepth - ( u%rlocal(i,j))* cos(u%psi(j) )) * p%airDens)) - p%Pvap ) / ( 0.5_ReKi * p%airDens * y%Vrel(i,j)**2) ! Critical value of Sigma, cavitation if we go over this
end if


if (SigmaCavitCrit < SigmaCavit) then
call WrScr( NewLine//'Cavitation occured at node # = '//trim(num2lstr(i)//'and blade # = '//trim(num2lstr(j))))
end if

m%SigmaCavit(i,j)= SigmaCavit
m%SigmaCavitCrit(i,j)=SigmaCavitCrit

end if
end if


Expand Down Expand Up @@ -1803,6 +1850,7 @@ subroutine BEMT_UnCoupledSolve( phi, numBlades, airDens, mu, AFInfo, rlocal, cho
integer :: i, TestRegionResult
logical :: IsValidSolution
real(ReKi) :: Re, Vrel


ErrStat = ErrID_None
ErrMsg = ""
Expand Down Expand Up @@ -1931,5 +1979,8 @@ end subroutine BEMT_UnCoupledSolve



end module BEMT

end module BEMT




Loading

0 comments on commit 055be86

Please sign in to comment.