Pull branch

modules-local/aerodyn/src/AeroDyn.f90

@@ -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()
@@ -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...
@@ -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 )   
@@ -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
@@ -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()
@@ -3752,6 +3770,7 @@ FUNCTION CheckBEMTInputPerturbations( p, m ) RESULT(ValidPerturb)
 
       end if
 
+
    else ! not UseInduction
 
          do k=1,p%NumBlades      

modules-local/aerodyn/src/AeroDyn_Driver.f90

@@ -197,4 +197,4 @@ subroutine Dvr_End()
    end subroutine Dvr_End
 !................................   
 end program AeroDyn_Driver
-
+

modules-local/aerodyn/src/AeroDyn_Registry.txt

@@ -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

modules-local/aerodyn/src/AirfoilInfo.f90

@@ -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
@@ -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 )
@@ -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
 
@@ -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.
@@ -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
 

modules-local/aerodyn/src/AirfoilInfo_Registry.txt

@@ -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):

modules-local/aerodyn/src/BEMT.f90

@@ -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
@@ -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 
@@ -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
 
@@ -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
 
@@ -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
 
@@ -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
 
@@ -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
@@ -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
 
@@ -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()
@@ -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
@@ -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
@@ -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
 
 
@@ -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  = ""
@@ -1931,5 +1979,8 @@ end subroutine BEMT_UnCoupledSolve
 
 
 
-end module BEMT
-
+    end module BEMT
+
+
+
+