Skip to content

Commit 055be86

Browse files
robynnemurrayNRELghaymanNREL
robynnemurrayNREL
authored and
ghaymanNREL
committed
Pull branch (#18)
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.
1 parent 2d5398b commit 055be86

10 files changed

Lines changed: 733 additions & 454 deletions

‎modules-local/aerodyn/src/AeroDyn.f90‎

Lines changed: 27 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -866,7 +866,8 @@ subroutine SetParameters( InitInp, InputFileData, p, ErrStat, ErrMsg )
866866

867867
p%AirDens = InputFileData%AirDens
868868
p%KinVisc = InputFileData%KinVisc
869-
p%SpdSound = InputFileData%SpdSound
869+
870+
870871

871872
!p%AFI ! set in call to AFI_Init() [called early because it wants to use the same echo file as AD]
872873
!p%BEMT ! set in call to BEMT_Init()
@@ -1417,7 +1418,11 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg )
14171418
if (InputFileData%AirDens <= 0.0) call SetErrStat ( ErrID_Fatal, 'The air density (AirDens) must be greater than zero.', ErrStat, ErrMsg, RoutineName )
14181419
if (InputFileData%KinVisc <= 0.0) call SetErrStat ( ErrID_Fatal, 'The kinesmatic viscosity (KinVisc) must be greater than zero.', ErrStat, ErrMsg, RoutineName )
14191420
if (InputFileData%SpdSound <= 0.0) call SetErrStat ( ErrID_Fatal, 'The speed of sound (SpdSound) must be greater than zero.', ErrStat, ErrMsg, RoutineName )
1420-
1421+
if (InputFileData%Pvap <= 0.0) call SetErrStat ( ErrID_Fatal, 'The vapour pressure (Pvap) must be greater than zero.', ErrStat, ErrMsg, RoutineName )
1422+
if (InputFileData%Patm <= 0.0) call SetErrStat ( ErrID_Fatal, 'The atmospheric pressure (Patm) must be greater than zero.', ErrStat, ErrMsg, RoutineName )
1423+
if (InputFileData%FluidDepth <= 0.0) call SetErrStat ( ErrID_Fatal, 'Fluid depth (FluidDepth) cannot be negative', ErrStat, ErrMsg, RoutineName )
1424+
1425+
14211426

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

14401445
if (.not. InputFileData%FLookUp ) call SetErrStat( ErrID_Fatal, 'FLookUp must be TRUE for this version.', ErrStat, ErrMsg, RoutineName )
14411446
end if
1442-
1447+
1448+
if ( InputFileData%CavitCheck .and. InputFileData%AFAeroMod == 2) then
1449+
call SetErrStat( ErrID_Fatal, 'Cannot use unsteady aerodynamics module with a cavitation check', ErrStat, ErrMsg, RoutineName )
1450+
end if
1451+
1452+
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 )
1453+
1454+
14431455

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

17011713
InitInp%airDens = InputFileData%AirDens
1702-
InitInp%kinVisc = InputFileData%KinVisc
1714+
InitInp%kinVisc = InputFileData%KinVisc
1715+
InitInp%Patm = InputFileData%Patm
1716+
InitInp%Pvap = InputFileData%Pvap
1717+
InitInp%FluidDepth = InputFileData%FluidDepth
17031718
InitInp%skewWakeMod = InputFileData%SkewMod
17041719
InitInp%aTol = InputFileData%IndToler
17051720
InitInp%useTipLoss = InputFileData%TipLoss
@@ -1748,10 +1763,13 @@ SUBROUTINE Init_BEMTmodule( InputFileData, u_AD, u, p, x, xd, z, OtherState, y,
17481763
end do
17491764
end do
17501765

1751-
InitInp%UA_Flag = InputFileData%AFAeroMod == AFAeroMod_BL_unsteady
1752-
InitInp%UAMod = InputFileData%UAMod
1753-
InitInp%Flookup = InputFileData%Flookup
1754-
InitInp%a_s = InputFileData%SpdSound
1766+
InitInp%UA_Flag = InputFileData%AFAeroMod == AFAeroMod_BL_unsteady
1767+
InitInp%UAMod = InputFileData%UAMod
1768+
InitInp%Flookup = InputFileData%Flookup
1769+
InitInp%a_s = InputFileData%SpdSound
1770+
InitInp%CavitCheck = InputFileData%CavitCheck
1771+
1772+
17551773

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

37583776
end if
37593777

3778+
37603779
else ! not UseInduction
37613780

37623781
do k=1,p%NumBlades

‎modules-local/aerodyn/src/AeroDyn_Driver.f90‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -197,4 +197,4 @@ subroutine Dvr_End()
197197
end subroutine Dvr_End
198198
!................................
199199
end program AeroDyn_Driver
200-
200+

‎modules-local/aerodyn/src/AeroDyn_IO.f90‎

Lines changed: 595 additions & 423 deletions
Large diffs are not rendered by default.

‎modules-local/aerodyn/src/AeroDyn_Registry.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -65,8 +65,12 @@ typedef ^ AD_InputFile IntKi TwrPotent - - - "Type tower influence on wind based
6565
typedef ^ AD_InputFile LOGICAL TwrShadow - - - "Calculate tower influence on wind based on downstream tower shadow?" -
6666
typedef ^ AD_InputFile LOGICAL TwrAero - - - "Calculate tower aerodynamic loads?" flag
6767
typedef ^ AD_InputFile Logical FrozenWake - - - "Flag that tells this module it should assume a frozen wake during linearization." -
68+
typedef ^ AD_InputFile Logical CavitCheck - - - "Flag that tells us if we want to check for cavitation" -
6869
typedef ^ AD_InputFile ReKi AirDens - - - "Air density" kg/m^3
6970
typedef ^ AD_InputFile ReKi KinVisc - - - "Kinematic air viscosity" m^2/s
71+
typedef ^ AD_InputFile ReKi Patm - - - "Atmospheric pressure" Pa
72+
typedef ^ AD_InputFile ReKi Pvap - - - "Vapour pressure" Pa
73+
typedef ^ AD_InputFile ReKi FluidDepth - - - "Submerged hub depth" m
7074
typedef ^ AD_InputFile ReKi SpdSound - - - "Speed of sound" m/s
7175
typedef ^ AD_InputFile IntKi SkewMod - - - "Type of skewed-wake correction model {1=uncoupled, 2=Pitt/Peters, 3=coupled} [used only when WakeMod=1]" -
7276
typedef ^ AD_InputFile LOGICAL TipLoss - - - "Use the Prandtl tip-loss model? [used only when WakeMod=1]" flag

‎modules-local/aerodyn/src/AirfoilInfo.f90‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -110,6 +110,7 @@ SUBROUTINE AFI_Init ( InitInput, p, ErrStat, ErrMsg, UnEcho )
110110
p%ColCd = 2
111111
p%ColCm = 0 ! These may or may not be used; initialize to zero in case they aren't used
112112
p%ColCpmin = 0 ! These may or may not be used; initialize to zero in case they aren't used
113+
113114
IF ( InitInput%InCol_Cm > 0 ) THEN
114115
p%ColCm = 3
115116
IF ( InitInput%InCol_Cpmin > 0 ) THEN
@@ -119,8 +120,8 @@ SUBROUTINE AFI_Init ( InitInput, p, ErrStat, ErrMsg, UnEcho )
119120
p%ColCpmin = 3
120121
END IF
121122
NumCoefs = MAX(p%ColCd, p%ColCm,p%ColCpmin) ! number of non-zero coefficient columns
122-
123123

124+
124125
! Process the airfoil files.
125126

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

133+
p%AFInfo( :)%ColCpmin=p%ColCpmin
134+
p%AFInfo( :)%ColCm=p%ColCm
135+
132136

133137
DO File=1,InitInput%NumAFfiles
134138

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

498502

503+
499504
! Local declarations.
500505

501506
REAL(ReKi) :: Coords (2) ! An array to hold data from the airfoil-shape table.
@@ -805,6 +810,8 @@ SUBROUTINE ReadAFfile ( AFfile, NumCoefs, InCol_Alfa, InCol_Cl, InCol_Cd, InCol_
805810
CALL Cleanup()
806811
RETURN
807812
ENDIF
813+
814+
808815

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

‎modules-local/aerodyn/src/AirfoilInfo_Registry.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -94,6 +94,8 @@ typedef ^ ^ INTEGER NumCpminAoAkts - - - "The number of angle-of-attack knots fo
9494
typedef ^ ^ INTEGER NumCpminReKts - - - "The number of log(Re) knots for 2D splines of Cpmin" -
9595
typedef ^ ^ INTEGER NumTabs - - - "The number of airfoil tables in the airfoil file" -
9696
typedef ^ ^ AFI_Table_Type Table {:} - - "The tables of airfoil data for given Re and control setting" -
97+
typedef ^ ^ INTEGER ColCpmin - - - "Column number for Cpmin" -
98+
typedef ^ ^ INTEGER ColCm - - - "Column number for Cm" -
9799

98100
# ..... Initialization data .......................................................................................................
99101
# The following derived type stores information that comes from the calling module (say, AeroDyn):

‎modules-local/aerodyn/src/BEMT.f90‎

Lines changed: 62 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -162,6 +162,9 @@ subroutine BEMT_SetParameters( InitInp, p, errStat, errMsg )
162162
p%numBladeNodes = InitInp%numBladeNodes
163163
p%numBlades = InitInp%numBlades
164164
p%UA_Flag = InitInp%UA_Flag
165+
p%CavitCheck = InitInp%CavitCheck
166+
167+
165168

166169
allocate ( p%chord(p%numBladeNodes, p%numBlades), STAT = errStat2 )
167170
if ( errStat2 /= 0 ) then
@@ -205,10 +208,13 @@ subroutine BEMT_SetParameters( InitInp, p, errStat, errMsg )
205208
end do
206209
end do
207210

208-
209-
!p%DT = InitInp%DT
211+
212+
!p%DT = InitInp%DT
210213
p%airDens = InitInp%airDens
211-
p%kinVisc = InitInp%kinVisc
214+
p%kinVisc = InitInp%kinVisc
215+
p%Patm = InitInp%Patm
216+
p%Pvap = InitInp%Pvap
217+
p%FluidDepth = InitInp%FluidDepth
212218
p%skewWakeMod = InitInp%skewWakeMod
213219
p%useTipLoss = InitInp%useTipLoss
214220
p%useHubLoss = InitInp%useHubLoss
@@ -219,6 +225,7 @@ subroutine BEMT_SetParameters( InitInp, p, errStat, errMsg )
219225
p%numReIterations = InitInp%numReIterations
220226
p%maxIndIterations = InitInp%maxIndIterations
221227
p%aTol = InitInp%aTol
228+
222229

223230
end subroutine BEMT_SetParameters
224231

@@ -431,14 +438,15 @@ end subroutine BEMT_AllocInput
431438

432439

433440
!----------------------------------------------------------------------------------------------------------------------------------
434-
subroutine BEMT_AllocOutput( y, p, errStat, errMsg )
441+
subroutine BEMT_AllocOutput( y, p, m, errStat, errMsg )
435442
! This routine is called from BEMT_Init.
436443
!
437444
!
438445
!..................................................................................................................................
439446

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

@@ -465,6 +473,10 @@ subroutine BEMT_AllocOutput( y, p, errStat, errMsg )
465473
call allocAry( y%Cm, p%numBladeNodes, p%numBlades, 'y%Cm', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
466474
call allocAry( y%Cl, p%numBladeNodes, p%numBlades, 'y%Cl', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
467475
call allocAry( y%Cd, p%numBladeNodes, p%numBlades, 'y%Cd', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
476+
call allocAry( m%Cpmin, p%numBladeNodes, p%numBlades, 'm%Cpmin', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
477+
call allocAry( m%SigmaCavit, p%numBladeNodes, p%numBlades, 'm%SigmaCavit', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
478+
call allocAry( m%SigmaCavitCrit, p%numBladeNodes, p%numBlades, 'm%SigmaCavitCrit', errStat2, errMsg2); call setErrStat(errStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
479+
468480

469481
if (ErrStat >= AbortErrLev) RETURN
470482

@@ -482,7 +494,8 @@ subroutine BEMT_AllocOutput( y, p, errStat, errMsg )
482494
y%tanInduction = 0.0_ReKi
483495
y%AOA = 0.0_ReKi
484496
y%Cl = 0.0_ReKi
485-
y%Cd = 0.0_ReKi
497+
y%Cd = 0.0_ReKi
498+
486499

487500
end subroutine BEMT_AllocOutput
488501

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

707720
#endif
721+
708722

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

743+
744+
729745
end do
730746
end do
731747

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

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

10831099

1084-
real(ReKi) :: Re, fzero
1085-
real(ReKi) :: Rtip ! maximum rlocal value for node j over all blades
1100+
real(ReKi) :: Re, fzero, theta, Vx, Vy
1101+
real(ReKi) :: Rtip, SigmaCavitCrit, SigmaCavit ! maximum rlocal value for node j over all blades
10861102

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

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

1181+
! local velocities and twist angle
1182+
Vx = u%Vx(i,j)
1183+
Vy = u%Vy(i,j)
1184+
1185+
1186+
11651187
! Set the active blade element for UnsteadyAero
11661188
m%UA%iBladeNode = i
11671189
m%UA%iBlade = j
@@ -1235,9 +1257,34 @@ subroutine BEMT_CalcOutput( t, u, p, x, xd, z, OtherState, AFInfo, y, m, errStat
12351257
else
12361258
! 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.
12371259
call ComputeSteadyAirfoilCoefs( y%AOA(i,j), y%Re(i,j), AFInfo(p%AFindx(i,j)), &
1238-
y%Cl(i,j), y%Cd(i,j), y%Cm(i,j), errStat2, errMsg2 )
1260+
y%Cl(i,j), y%Cd(i,j), y%Cm(i,j), m%Cpmin(i,j), errStat2, errMsg2 )
12391261
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName//trim(NodeTxt))
12401262
if (errStat >= AbortErrLev) return
1263+
1264+
1265+
1266+
1267+
! 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
1268+
if ( p%CavitCheck ) then
1269+
SigmaCavit= -1* m%Cpmin(i,j) ! Cavitation number on blade node j
1270+
1271+
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
1272+
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 )
1273+
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
1274+
1275+
else
1276+
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
1277+
end if
1278+
1279+
1280+
if (SigmaCavitCrit < SigmaCavit) then
1281+
call WrScr( NewLine//'Cavitation occured at node # = '//trim(num2lstr(i)//'and blade # = '//trim(num2lstr(j))))
1282+
end if
1283+
1284+
m%SigmaCavit(i,j)= SigmaCavit
1285+
m%SigmaCavitCrit(i,j)=SigmaCavitCrit
1286+
1287+
end if
12411288
end if
12421289

12431290

@@ -1803,6 +1850,7 @@ subroutine BEMT_UnCoupledSolve( phi, numBlades, airDens, mu, AFInfo, rlocal, cho
18031850
integer :: i, TestRegionResult
18041851
logical :: IsValidSolution
18051852
real(ReKi) :: Re, Vrel
1853+
18061854

18071855
ErrStat = ErrID_None
18081856
ErrMsg = ""
@@ -1931,5 +1979,8 @@ end subroutine BEMT_UnCoupledSolve
19311979

19321980

19331981

1934-
end module BEMT
1935-
1982+
end module BEMT
1983+
1984+
1985+
1986+

0 commit comments

Comments
 (0)