AD: updates from new bem momentum corr

modules/aerodyn/src/AeroDyn.f90

@@ -374,9 +374,14 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut
       ! set the rest of the parameters
    p%SkewMod = InputFileData%SkewMod
    do iR = 1, nRotors
-      !p%rotors(iR)%AeroProjMod = InitInp%rotors(iR)%AeroProjMod
-      p%rotors(iR)%AeroProjMod = AeroProjMod(iR)
+      p%rotors(iR)%AeroProjMod = InitInp%rotors(iR)%AeroProjMod
+      !p%rotors(iR)%AeroProjMod = AeroProjMod(iR)
       p%rotors(iR)%AeroBEM_Mod = InitInp%rotors(iR)%AeroBEM_Mod
+      call WrScr('   AeroDyn: projMod: '//trim(num2lstr(p%rotors(iR)%AeroProjMod))//', BEM_Mod:'//trim(num2lstr(p%rotors(iR)%AeroBEM_Mod)))      
+      if (AeroProjMod(iR) == APM_BEM_Polar .and. InputFileData%WakeMod==WakeMod_FVW) then
+         call WrScr('AeroProj cannot be 2 (somehow) when using OLAF - Aborting')
+         STOP
+      endif
       call SetParameters( InitInp, InputFileData, InputFileData%rotors(iR), p%rotors(iR), p, ErrStat2, ErrMsg2 )
       if (Failed()) return;
    enddo
@@ -3545,7 +3550,7 @@ subroutine SetOutputsFromBEMT( p, u, m, y )
    real(reki)                              :: c                      ! local chord length 
    real(reki)                              :: aoa                    ! local angle of attack 
    real(reki)                              :: Cl,Cd,Cm               ! local airfoil lift, drag and pitching moment coefficients 
-   real(reki)                              :: Cn,Ct                  ! local airfoil normal and tangential force coefficients 
+   real(reki)                              :: Cxa,Cya                ! local airfoil normal and tangential force coefficients 
 
 
    do k=1,p%NumBlades
@@ -3556,14 +3561,14 @@ subroutine SetOutputsFromBEMT( p, u, m, y )
          Cl  = m%BEMT_y%cl(j,k)
          Cd  = m%BEMT_y%cd(j,k)
          Cm  = m%BEMT_y%cm(j,k)
-         Cn  =  Cl*cos(aoa) + Cd*sin(aoa)
-         Ct  = -Cl*sin(aoa) + Cd*cos(aoa) ! NOTE: this is not Ct but Cy_a (y_a going towards the TE)
+         Cxa =  Cl*cos(aoa) + Cd*sin(aoa)
+         Cya = -Cl*sin(aoa) + Cd*cos(aoa)
 
          ! Dimensionalize the aero forces and moment
          q = 0.5 * p%airDens * m%BEMT_y%Vrel(j,k)**2              ! dynamic pressure of the jth node in the kth blade
          c = p%BEMT%chord(j,k)
-         forceAirfoil(1)  = Cn * q * c
-         forceAirfoil(2)  = Ct * q * c
+         forceAirfoil(1)  = Cxa * q * c
+         forceAirfoil(2)  = Cya * q * c
          forceAirfoil(3)  = 0.0_reki
          momentAirfoil(1) = 0.0_reki
          momentAirfoil(2) = 0.0_reki
@@ -4429,9 +4434,18 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x
          InitInp%BEM_Mod    = -1
          call SetErrStat(ErrID_Fatal, "AeroProjMod needs to be 1 or 2 when used with BEM", ErrStat, ErrMsg, RoutineName)   
       endif
+   else if (p%AeroBEM_Mod== BEMMod_3D_MomCorr) then
+         InitInp%BEM_Mod  = BEMMod_3D
+         InitInp%MomentumCorr = .TRUE. 
    endif
    p%AeroBEM_Mod = InitInp%BEM_Mod ! Very important, for consistency
-   !call WrScr('   AeroDyn: projMod: '//trim(num2lstr(p%AeroProjMod))//', BEM_Mod:'//trim(num2lstr(InitInp%BEM_Mod)))
+   ! 
+
+   if (InitInp%MomentumCorr) then
+      call WrScr('   AeroDyn: projMod: '//trim(num2lstr(p%AeroProjMod))//', BEM_Mod:'//trim(num2lstr(InitInp%BEM_Mod))//', Momentum Correction')
+   else 
+      call WrScr('   AeroDyn: projMod: '//trim(num2lstr(p%AeroProjMod))//', BEM_Mod:'//trim(num2lstr(InitInp%BEM_Mod)))
+   endif
       ! remove the ".AD" from the RootName
    k = len_trim(InitInp%RootName)
    if (k>3) then

modules/aerodyn/src/AeroDyn_Driver_Subs.f90

@@ -504,7 +504,7 @@ subroutine Init_ADI_ForDriver(iCase, ADI, dvr, FED, dt, errStat, errMsg)
             InitInp%AD%rotors(iWT)%AeroProjMod = wt%projMod
          endif
          InitInp%AD%rotors(iWT)%AeroBEM_Mod = wt%BEM_Mod
-         !call WrScr('   Driver:  projMod: '//trim(num2lstr(InitInp%AD%rotors(iWT)%AeroProjMod))//', BEM_Mod:'//trim(num2lstr(InitInp%AD%rotors(iWT)%AeroBEM_Mod)))
+         call WrScr('   Driver:  projMod: '//trim(num2lstr(InitInp%AD%rotors(iWT)%AeroProjMod))//', BEM_Mod:'//trim(num2lstr(InitInp%AD%rotors(iWT)%AeroBEM_Mod)))
          InitInp%AD%rotors(iWT)%HubPosition    = y_ED%HubPtMotion%Position(:,1)
          InitInp%AD%rotors(iWT)%HubOrientation = y_ED%HubPtMotion%RefOrientation(:,:,1)
          InitInp%AD%rotors(iWT)%NacellePosition    = y_ED%NacelleMotion%Position(:,1)
@@ -1012,6 +1012,9 @@ subroutine Dvr_ReadInputFile(fileName, dvr, errStat, errMsg )
          wt%BEM_Mod = -1
       else
          call ParseVar(FileInfo_In, CurLine, 'BEM_Mod'//sWT    , wt%BEM_Mod     , errStat2, errMsg2, unEc); if(Failed()) return
+        ! call WrScr('>>> Driver: ProjMod and BEM_Mod are present in AeroDyn driver input file. ProjMod: '//trim(num2lstr(wt%projMod))//' BEM_Mod: '//trim(num2lstr(wt%BEM_Mod)))
+         print*,'>>> Driver:  ProjMod, BEM_Mod:',wt%ProjMod, wt%BEM_Mod
+
       endif
       call ParseVar(FileInfo_In, CurLine, 'BasicHAWTFormat'//sWT    , wt%basicHAWTFormat       , errStat2, errMsg2, unEc); if(Failed()) return
 

modules/aerodyn/src/AeroDyn_IO.f90

@@ -2168,19 +2168,14 @@ subroutine calcCantAngle(f, xi,stencilSize,n,cantAngle)
     implicit none
     integer(IntKi), intent(in)  :: stencilSize, n 
     integer(IntKi)              :: i, j
-    integer(IntKi)              :: sortInd(n)
     integer(IntKi)              :: info
     real(ReKi),  intent(in)     :: f(n), xi(n)
     real(ReKi)                  :: cx(stencilSize), cf(stencilSize), xiIn(stencilSize)
-    real(ReKi)                  :: fIn(stencilSize), cPrime(n), fPrime(n), xiAbs(n)
+    real(ReKi)                  :: fIn(stencilSize), cPrime(n), fPrime(n)
     real(ReKi), intent(inout)   :: cantAngle(n)
-    real(ReKi), parameter       :: ThisTol = 1e-6
 
     do i = 1,size(xi)
 
-        xiAbs = abs(xi-xi(i))
-        call hpsort_eps_epw (n, xiAbs, sortInd, ThisTol)
-
         if (i.eq.1) then
             fIn = f(1:stencilSize)
             xiIn = xi(1:stencilSize)
@@ -2417,147 +2412,5 @@ subroutine r8vm_sl ( n, a, b, x, job, info )
   return
 end subroutine r8vm_sl
 
-!                                                                            
-  ! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino 
-  ! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino  
-  !                                                                            
-  ! This file is distributed under the terms of the GNU General Public         
-  ! License. See the file `LICENSE' in the root directory of the               
-  ! present distribution, or http://www.gnu.org/copyleft.gpl.txt .             
-  !                                                                            
-  ! Adapted from flib/hpsort_eps
-  !---------------------------------------------------------------------
-  subroutine hpsort_eps_epw (n, ra, ind, eps)
-  !---------------------------------------------------------------------
-  ! sort an array ra(1:n) into ascending order using heapsort algorithm,
-  ! and considering two elements being equal if their values differ
-  ! for less than "eps".
-  ! n is input, ra is replaced on output by its sorted rearrangement.
-  ! create an index table (ind) by making an exchange in the index array
-  ! whenever an exchange is made on the sorted data array (ra).
-  ! in case of equal values in the data array (ra) the values in the
-  ! index array (ind) are used to order the entries.
-  ! if on input ind(1)  = 0 then indices are initialized in the routine,
-  ! if on input ind(1) != 0 then indices are assumed to have been
-  !                initialized before entering the routine and these
-  !                indices are carried around during the sorting process
-  !
-  ! no work space needed !
-  ! free us from machine-dependent sorting-routines !
-  !
-  ! adapted from Numerical Recipes pg. 329 (new edition)
-  !
-  !use kinds, ONLY : DP
-  implicit none  
-  !-input/output variables
-  integer(IntKi), intent(in)    :: n  
-  real(ReKi), intent(in)        :: eps
-  integer(IntKi)                :: ind (n)  
-  real(ReKi)                    :: ra (n)
-  !-local variables
-  integer(IntKi)                :: i, ir, j, l, iind  
-  real(ReKi)                    :: rra  
-!
-  ! initialize index array
-  IF (ind (1) .eq.0) then  
-     DO i = 1, n  
-        ind (i) = i  
-     ENDDO
-  ENDIF
-  ! nothing to order
-  IF (n.lt.2) return  
-  ! initialize indices for hiring and retirement-promotion phase
-  l = n / 2 + 1  
-
-  ir = n  
-
-  sorting: do 
-
-    ! still in hiring phase
-    IF ( l .gt. 1 ) then  
-       l    = l - 1  
-       rra  = ra (l)  
-       iind = ind (l)  
-       ! in retirement-promotion phase.
-    ELSE  
-       ! clear a space at the end of the array
-       rra  = ra (ir)  
-       !
-       iind = ind (ir)  
-       ! retire the top of the heap into it
-       ra (ir) = ra (1)  
-       !
-       ind (ir) = ind (1)  
-       ! decrease the size of the corporation
-       ir = ir - 1  
-       ! done with the last promotion
-       IF ( ir .eq. 1 ) then  
-          ! the least competent worker at all !
-          ra (1)  = rra  
-          !
-          ind (1) = iind  
-          exit sorting  
-       ENDIF
-    ENDIF
-    ! wheter in hiring or promotion phase, we
-    i = l  
-    ! set up to place rra in its proper level
-    j = l + l  
-    !
-    DO while ( j .le. ir )  
-       IF ( j .lt. ir ) then  
-          ! compare to better underling
-          IF ( hslt( ra (j),  ra (j + 1) ) ) then  
-             j = j + 1  
-          !else if ( .not. hslt( ra (j+1),  ra (j) ) ) then
-             ! this means ra(j) == ra(j+1) within tolerance
-           !  if (ind (j) .lt.ind (j + 1) ) j = j + 1
-          ENDIF
-       ENDIF
-       ! demote rra
-       IF ( hslt( rra, ra (j) ) ) then  
-          ra (i) = ra (j)  
-          ind (i) = ind (j)  
-          i = j  
-          j = j + j  
-       !else if ( .not. hslt ( ra(j) , rra ) ) then
-          !this means rra == ra(j) within tolerance
-          ! demote rra
-         ! if (iind.lt.ind (j) ) then
-         !    ra (i) = ra (j)
-         !    ind (i) = ind (j)
-         !    i = j
-         !    j = j + j
-         ! else
-             ! set j to terminate do-while loop
-         !    j = ir + 1
-         ! endif
-          ! this is the right place for rra
-       ELSE
-          ! set j to terminate do-while loop
-          j = ir + 1  
-       ENDIF
-    ENDDO
-    ra (i) = rra  
-    ind (i) = iind  
-
-  END DO sorting    
-contains 
-
-  !  internal function 
-  !  compare two real number and return the result
-
-  logical function hslt( a, b )
-    REAL(ReKi) :: a, b
-    IF( abs(a-b) <  eps ) then
-      hslt = .false.
-    ELSE
-      hslt = ( a < b )
-    end if
-  end function hslt
-
-  !
-end subroutine hpsort_eps_epw
-
 !----------------------------------------------------------------------------------------------------------------------------------
 END MODULE AeroDyn_IO

modules/aerodyn/src/BEMT.f90

@@ -90,6 +90,21 @@ real(ReKi) function ComputePhiWithInduction( Vx, Vy, a, aprime, cantAngle, xVelC
 
 end function ComputePhiWithInduction
 
+subroutine ComputePhiFromInductions(u, p, phi, axInduction, tanInduction)
+   type(BEMT_InputType),         intent(in   ) :: u
+   type(BEMT_ParameterType),     intent(in   ) :: p
+   real(ReKi),                   intent(inout) :: phi(:,:)
+   real(ReKi),                   intent(in   ) :: axInduction(:,:)
+   real(ReKi),                   intent(in   ) :: tanInduction(:,:)
+   integer(IntKi)                              :: i, j
+   do j = 1,p%numBlades ! Loop through all blades
+      do i = 1,p%numBladeNodes ! Loop through the blade nodes / elements
+         phi(i,j) = ComputePhiWithInduction( u%Vx(i,j), u%Vy(i,j),  axInduction(i,j), tanInduction(i,j), u%cantAngle(i,j), u%xVelCorr(i,j) )
+      enddo             ! I - Blade nodes / elements
+   enddo          ! J - All blades
+end subroutine ComputePhiFromInductions
+
+
 !----------------------------------------------------------------------------------------------------------------------------------   
 subroutine BEMT_Set_UA_InitData( InitInp, interval, Init_UA_Data, errStat, errMsg )
 ! This routine is called from BEMT_Init.
@@ -175,7 +190,7 @@ subroutine BEMT_SetParameters( InitInp, p, errStat, errMsg )
    p%MomentumCorr   = InitInp%MomentumCorr
    p%BEM_Mod        = InitInp%BEM_Mod
    !call WrScr('>>>> BEM_Mod '//trim(num2lstr(p%BEM_Mod)))
-   if ((p%BEM_Mod/=BEMMod_2D .and. p%BEM_Mod/=BEMMod_3D )) then
+   if (.not.(ANY( p%BEM_Mod == (/BEMMod_2D, BEMMod_3D/)))) then
       call SetErrStat( ErrID_Fatal, 'BEM_Mod needs to be 0 or 2 for now', errStat, errMsg, RoutineName )
       return
    endif

modules/aerodyn/src/BEMTUncoupled.f90

@@ -59,6 +59,8 @@ module BEMTUnCoupled
    public :: GetEulerAnglesFromOrientation
 
    public :: VelocityIsZero
+
+   public :: BEMTU_Test_ACT_Relationship
 contains
 
 !..................................................................................................................................   
@@ -227,6 +229,20 @@ subroutine computeAirfoilOperatingAOA( BEM_Mod, phi, theta, cantAngle, toeAngle,
 
 
 end subroutine computeAirfoilOperatingAOA
+! ---
+!> Angle of attack in the airfoil reference frame
+real(ReKi) function computeAirfoilAOA(Vrel_a) result(AoA)
+   real(ReKi), intent(in   ) :: Vrel_a(3)
+   real(ReKi)                :: numer, denom, ratio
+   ! Determine angle of attack as angle between airfoil chordline (afAxialVec) and inflow
+   numer = Vrel_a(2)
+   denom = sqrt(Vrel_a(1)**2 + Vrel_a(2)**2)
+   ratio = numer / denom
+   AoA = acos( max( min( ratio, 1.0_ReKi ), -1.0_ReKi ) )
+   AoA = sign( AoA, Vrel_a(1) )
+end function computeAirfoilAOA
+
+
 !..................................................................................................................................
 !> Transform the aerodynamic coefficients (Cl,Cd,Cm) (directed based on Vrel_xy_a )
 !! from the airfoil coordinate system (a) to the without sweep pitch coordinate system (w)
@@ -476,7 +492,14 @@ subroutine ApplySkewedWakeCorrection(BEM_Mod, SkewMod, yawCorrFactor, F, azimuth
    end if
 
       !bjj: modified 22-Sep-2015: RRD recommends 32 instead of 64 in the denominator (like AD14)
+   ! TODO TODO TODO
+   ! ADLEG:
+   !if(BEM_Mod==BEMMod_2D) then
    yawCorr = ( yawCorrFactor * yawCorr_tan * (tipRatio) * sin(azimuth) ) ! bjj: note that when chi gets close to +/-pi this blows up
+   !else
+   !   ! ADENV:
+   !   yawCorr = ( yawCorrFactor * F * yawCorr_tan * (tipRatio) * cos(azimuth-azimuthOffset) ) ! bjj: note that when chi gets close to +/-pi this blows up
+   !endif
 
    a = a * (1.0 +  yawCorr)
 
@@ -1216,4 +1239,52 @@ FUNCTION GetEulerAnglesFromOrientation(EulerDCM,orientation) RESULT(theta)
 end function
 !-----------------------------------------------------------------------------------------
 
+
+
+!> Simple test for a-Ct relationship. 
+subroutine BEMTU_Test_ACT_Relationship()
+   real(R8Ki) :: chi0
+   real(R8Ki) :: delta_chi
+   real(ReKi) :: F
+   logical :: skewConvention
+   integer :: i
+   integer :: iUnit
+   real(R8Ki)              :: c2, c1, c0 ! Empirical CT = c2*a^2 + c1*a + c0 for a > a0
+   ! Get Coefficients for Empirical CT
+   iUnit = 123
+
+   ! --- No Momentum Corr, F=1
+   F=1; skewConvention=.False.
+   call parametricStudy('ACTCoeffs_F10_NoCo.csv')
+   ! --- No Momentum Corr, F=0.5
+   F=0.5; skewConvention=.False.
+   call parametricStudy('ACTCoeffs_F05_NoCo.csv')
+   ! --- Momentum Corr, F=1
+   F=1; skewConvention=.True.
+   call parametricStudy('ACTCoeffs_F10_Corr.csv')
+   ! --- Momentum Corr, F=0.5
+   F=0.5; skewConvention=.True.
+   call parametricStudy('ACTCoeffs_F05_Corr.csv')
+
+   STOP
+
+contains
+   subroutine parametricStudy(filename)
+      character(len=*) :: filename
+      chi0=-50 * D2R
+      open(unit=iUnit, file=filename)
+      write(iUnit, '(5(A15))') 'chi0', 'c0', 'c1', 'c2', 'F'
+      do i=1,21
+         call getEmpiricalCoefficients(chi0 ,F , c0, c1, c2, skewConvention)
+         write(iUnit,'(5(F15.5))') chi0*R2D, c0, c1, c2, F
+         chi0 = chi0 + 5*D2R
+      enddo
+      close(iUnit)
+   end subroutine
+
+
+
+end subroutine BEMTU_Test_ACT_Relationship
+
+
 end module BEMTUncoupled

modules/aerodyn/src/BEMT_Registry.txt

@@ -24,6 +24,7 @@ param     BEMT/BEMT                   -                              INTEGER
 
 param     BEMT/BEMT                   -                              INTEGER                 BEMMod_2D                   -      0     -   "2D BEM assuming Cx, Cy, phi, L, D are in the same plane" -
 param     BEMT/BEMT                   -                              INTEGER                 BEMMod_3D                   -      2     -   "3D BEM assuming a momentum balance system, and an airfoil system" -
+param     BEMT/BEMT                   -                              INTEGER                 BEMMod_3D_MomCorr           -      250   -   "3D BEM - Momentum Correction" 
 
 #
 #

modules/aerodyn/src/BEMT_Types.f90

@@ -43,6 +43,7 @@ MODULE BEMT_Types
     INTEGER(IntKi), PUBLIC, PARAMETER  :: SkewMod_PittPeters_Cont = 4      ! Pitt/Peters continuous formulation [-]
     INTEGER(IntKi), PUBLIC, PARAMETER  :: BEMMod_2D = 0      ! 2D BEM assuming Cx, Cy, phi, L, D are in the same plane [-]
     INTEGER(IntKi), PUBLIC, PARAMETER  :: BEMMod_3D = 2      ! 3D BEM assuming a momentum balance system, and an airfoil system [-]
+    INTEGER(IntKi), PUBLIC, PARAMETER  :: BEMMod_3D_MomCorr = 250      ! 3D BEM - Momentum Correction [-]
 ! =========  BEMT_InitInputType  =======
   TYPE, PUBLIC :: BEMT_InitInputType
     REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE  :: chord      !< Chord length at node [m]