Merge pull request #3307 from OpenFAST/rc-5.0.1

Cross merge rc-5.0.1 progress to dev

Author: andrew-platt

docs/OtherSupporting/HydroDyn/Breaking_Wave_Modeling_Approach_for_FAST.docx

@@ -41,7 +41,7 @@ END of input file (the word "END" must appear in the first 3 columns of this las
 "Clrnc"     - Tower clearance at each node (based on the absolute distance to the nearest point in the tower from blade node B#N# minus the local tower radius, in the deflected configuration); please note that this clearance is only approximate because the calculation assumes that the blade is a line with no volume (however, the calculation does use the local tower radius); when blade node B#N# is above the tower top (or below the tower base), the absolute distance to the tower top (or base) minus the local tower radius, in the deflected configuration, is output
 "Vx"        - Local axial velocity
 "Vy"        - Local tangential velocity
-"GeomPhi"   - Geometric phi? If phi was solved using normal BEMT equations, GeomPhi = 1; otherwise, if it was solved geometrically, GeomPhi = 0.
+"GeomPhi"   - Geometric phi? If phi was solved using normal BEMT equations, GeomPhi = 0; otherwise, if it was solved geometrically, GeomPhi = 1.
 "Chi"       - Skew angle (used in skewed wake correction) -- not available for OLAF
 "UA_Flag"   - Flag indicating if UA is turned on for this node. -- not available for OLAF
 "CpMin"     - Pressure coefficient

docs/OtherSupporting/OutListParameters.xlsx

@@ -37,7 +37,6 @@ of HydroDyn.
 - :download:`Implementation Plan - 2nd-order Forces Within HydroDyn <../../../OtherSupporting/HydroDyn/HydroDyn_2ndOrderForces_Plan.pdf>`
 - :download:`Implementation Plan - 2nd-order Wave Kinematics Within HydroDyn <../../../OtherSupporting/HydroDyn/WAVE2_document.pdf>`
 - :download:`Plan for Adding Wave Stretching to HydroDyn <../../../OtherSupporting/HydroDyn/HydroDyn_WaveStretching_Plan.docx>`
-- :download:`Breaking Wave Modeling Approach for FAST <../../../OtherSupporting/HydroDyn/Breaking_Wave_Modeling_Approach_for_FAST.docx>`
 
 
 HydroDyn allows for multiple approaches for calculating the hydrodynamic

docs/source/user/aerodyn/examples/NodalOutputs.txt

@@ -429,11 +429,12 @@ must match the modes in the 1st-order WAMIT-style input files, i.e., the .1,
 .3, and .hst files. The convention for the WAMIT-style input files requires 
 the first 6 modes to be rigid-body modes always, and the additional generalized 
 modes are numbered from 7 onward. The total number of modes in the 1st-order 
-WAMIT-style input files should be **6+NAddDOF** for each body. Currently, only 
-one potential-flow body is allowed when additional degrees of freedom are used, 
-requiring **NBody = 1**. In this case, the choice of **NBodyMod** does not 
-matter. If more than one potential-flow bodies are present, **NAddDOF** should 
-be populated with **NBody** zeros. This because, internally, HydroDyn supports 
+WAMIT-style input files should be **6+NAddDOF** for each body. 
+
+Currently, only one potential-flow body is allowed when additional degrees of 
+freedom are used, requiring **NBody = 1**. In this case, the choice of **NBodyMod** 
+does not matter. If more than one potential-flow bodies are present, **NAddDOF** 
+should be populated with **NBody** zeros. This because, internally, HydroDyn supports 
 additional degrees of freedom for an arbitrary number of bodies with any of the 
 **NBodyMod** options. The limitation of one potential-flow body only is due 
 to constraints with coupling to ExtPtfm, the only module that can make use of the 
@@ -442,6 +443,11 @@ been implemented for 2nd-order loads; therefore, the 2nd-order WAMIT-style input
 files, if included, should always contain **6NBody** modes irrespective of 
 **NAddDOF**. Finally, when **NAddDOF > 0**, **WAMITULEN** must be 1.
 
+In addition to being limited to one potential-flow body, **NAddDOF > 0** is only 
+supported when **ExctnMod** = 0 or 1 and **RdtnMod** = 0 or 1. State-space 
+wave-excitation and wave-radiation models are not supported. **NAddDOF** is not 
+used when **PotMod** = 0 or 2.
+
 Currently, the HydroDyn loads for the additional degrees of freedom can only be 
 used by the ExtPtfm module of OpenFAST with externally generated structural 
 mass (M), stiffness (K), and damping (C) matrices. This workflow is intended for

docs/source/user/api_change.rst

@@ -1976,18 +1976,27 @@ subroutine AD_CalcWind(t, u, FLowField, p, m, o, Inflow, ErrStat, ErrMsg)
    ! OLAF points
    if (allocated(o%WakeLocationPoints) .and. allocated(Inflow%InflowWakeVel)) then
       ! If rotor is MHK, add water depth to z coordinate
-      if (p%FVW%MHK > 0) then
+      if (p%FVW%MHK /= MHK_None) then
          PosOffset = [0.0_ReKi, 0.0_ReKi, p%FVW%WtrDpth]
       else
          PosOffset = 0.0_ReKi
       end if
 
-      call IfW_FlowField_GetVelAcc(FlowField, StartNode, t, &
-                                   o%WakeLocationPoints, &
-                                   Inflow%InflowWakeVel, &
-                                   NoAcc, ErrStat2, ErrMsg2, &
-                                   BoxExceedAllow=.true., PosOffset=PosOffset)
-      if(Failed()) return
+      if (p%FVW%MHK /= MHK_None .and. p%CompSeaSt) then ! MHK turbines with waves
+         call WaveField_GetWaveVelAcc_AD(p%WaveField, m%WaveField_m, &
+                                         StartNode, t, &
+                                         o%WakeLocationPoints, &
+                                         Inflow%InflowWakeVel, &
+                                         NoAcc, ErrStat2, ErrMsg2)
+         if(Failed()) return
+      else
+         call IfW_FlowField_GetVelAcc(FlowField, StartNode, t, &
+                                      o%WakeLocationPoints, &
+                                      Inflow%InflowWakeVel, &
+                                      NoAcc, ErrStat2, ErrMsg2, &
+                                      BoxExceedAllow=.true., PosOffset=PosOffset)
+         if(Failed()) return
+      end if
       StartNode = StartNode + size(o%WakeLocationPoints)
    end if
 
@@ -2034,7 +2043,7 @@ subroutine AD_CalcWind_Rotor(t, u, FlowField, p, p_AD, m, RotInflow, StartNode,
          call WaveField_GetWaveVelAcc_AD(p_AD%WaveField, m%WaveField_m, StartNode, t, &
             real(u%HubMotion%TranslationDisp + u%HubMotion%Position, ReKi), &
             RotInflow%InflowOnHub, NoAcc, ErrStat2, ErrMsg2)
-         if(Failed()) return 
+         if(Failed()) return
       else
          RotInflow%InflowOnHub = 0.0_ReKi
       end if
@@ -2085,7 +2094,7 @@ subroutine AD_CalcWind_Rotor(t, u, FlowField, p, p_AD, m, RotInflow, StartNode,
          call IfW_FlowField_GetVelAcc(FlowField, StartNode, t, &
             real(u%HubMotion%TranslationDisp + u%HubMotion%Position, ReKi), &
             RotInflow%InflowOnHub, NoAcc, ErrStat2, ErrMsg2, PosOffset=PosOffset)
-         if(Failed()) return 
+         if(Failed()) return
       else
          RotInflow%InflowOnHub = 0.0_ReKi
       end if
@@ -2469,7 +2478,7 @@ subroutine RotCavtCrit(u, p, m, errStat, errMsg)
                SigmaCavitCrit = ( p%rotors(iR)%Patm + ( p%rotors(iR)%Gravity * ( abs( u%rotors(iR)%BladeMotion(j)%Position(3,i) + u%rotors(iR)%BladeMotion(j)%TranslationDisp(3,i) ) + p%rotors(iR)%MSL2SWL ) * p%rotors(iR)%airDens ) - p%rotors(iR)%Pvap ) / ( 0.5_ReKi * p%rotors(iR)%airDens * Vreltemp**2 )  ! Critical value of Sigma, cavitation occurs if local cavitation number is greater than this
 
                if ( ( SigmaCavitCrit < SigmaCavit ) .and. ( .not. ( m%rotors(iR)%CavitWarnSet(i,j) ) ) ) then     
-                  call WrScr( NewLine//'Cavitation occurred at blade '//trim(num2lstr(j))//' and node '//trim(num2lstr(i))//'.' )
+                  call WrScr( NewLine//'Cavitation occurred at rotor '//trim(num2lstr(iR))//', blade '//trim(num2lstr(j))//', and node '//trim(num2lstr(i))//'.' )
                   m%rotors(iR)%CavitWarnSet(i,j) = .true.
                end if 
 
@@ -3208,7 +3217,7 @@ end subroutine SetDisturbedInflow
 !! Loop on blade nodes and computed a weighted sector average inflow at each node
 subroutine SetSectAvgInflow(t, p, p_AD, u, RotInflow, m, errStat, errMsg)
    real(DbKi),                   intent(in   )  :: t                      !< Current simulation time in seconds
-   type(RotParameterType),       intent(in   )  :: p                      !< AD parameters
+   type(RotParameterType),       intent(in   )  :: p                      !< Rotor parameters
    type(AD_ParameterType),       intent(in   )  :: p_AD                   !< AD parameters
    type(RotInputType),           intent(in   )  :: u                      !< AD Inputs at Time
    type(RotInflowType),          intent(in   )  :: RotInflow              !< Rotor inflow at Time
@@ -3227,6 +3236,7 @@ subroutine SetSectAvgInflow(t, p, p_AD, u, RotInflow, m, errStat, errMsg)
    real(ReKi)              :: e_t(3)      !< Polar unit vector perpendicular to rHA_perp ("e_theta")
    real(ReKi)              :: psi         !< Azimuthal offset in the current sector, runs from -psi_bwd to psi_fwd
    real(ReKi)              :: dpsi        !< Azimuthal increment
+   real(ReKi)              :: PosOffset(3)!< IfW position offset for MHK turbines
    real(ReKi), allocatable :: SectPos(:,:)!< Points used to define a given sector (for a given blade node A)
    real(ReKi), allocatable :: SectVel(:,:)!< Inflow velocity at a given sector (Undisturbed and then disturbed)
    real(ReKi), allocatable :: SectAcc(:,:)!< Inflow velocity at a given sector (Undisturbed and then disturbed)
@@ -3293,7 +3303,21 @@ subroutine SetSectAvgInflow(t, p, p_AD, u, RotInflow, m, errStat, errMsg)
 
          ! --- Inflow on sector points
          ! Undisturbed
-         call IfW_FlowField_GetVelAcc(p_AD%FlowField, 1, t, SectPos, SectVel, SectAcc, errStat=errStat2, errMsg=errMsg2); if(Failed()) return
+         if (p%MHK /= MHK_None .and. p_AD%CompSeaSt) then ! MHK turbines with waves
+            call WaveField_GetWaveVelAcc_AD(p_AD%WaveField, m%WaveField_m, 1_IntKi, t, &
+                                            SectPos, SectVel, SectAcc, ErrStat2, ErrMsg2)
+            if(Failed()) return
+         else
+            if (p%MHK /= MHK_None) then
+               PosOffset = [0.0_ReKi, 0.0_ReKi, p%WtrDpth]
+            else
+               PosOffset = 0.0_ReKi
+            end if
+            call IfW_FlowField_GetVelAcc(p_AD%FlowField, 1_IntKi, t, &
+                                         SectPos, SectVel, SectAcc, errStat=errStat2, errMsg=errMsg2, &
+                                         PosOffset=PosOffset); if(Failed()) return
+         end if
+
          ! --- Option 1 Disturbed inflow Before averaging - SectVel is modified in place
          !if (p%TwrPotent /= TwrPotent_none .or. p%TwrShadow /= TwrShadow_none) then
          !   call TwrInflArray(p, u, RotInflow, m, SectPos, SectVel, errStat2, errMsg2); if(Failed()) return

docs/source/user/hydrodyn/index.rst

@@ -517,6 +517,7 @@ typedef	^	RotMiscVarType	ReKi	TFinVrel_i 	3	-	-	"Relative velocity at the refere
 typedef	^	RotMiscVarType	ReKi	TFinSTV_i 	3	-	-	"Structural velocity at the reference point of the fin in the inertial system"
 typedef	^	RotMiscVarType	ReKi	TFinF_i 	3	-	-	"Forces at the reference point of the fin in the inertial system"
 typedef	^	RotMiscVarType	ReKi	TFinM_i 	3	-	-	"Moments at the reference point of the fin in the inertial system"
+typedef	^	RotMiscVarType	GridInterp_MiscVarType	WaveField_m	-	-	-	"misc var information from the GridInterp module"
 
 typedef	^	MiscVarType	RotMiscVarType	rotors	{:}	-	-	"MiscVars for each rotor"	-
 typedef	^	MiscVarType	FVW_InputType	FVW_u	:	-	-	"Inputs to the FVW module"	-
@@ -525,7 +526,7 @@ typedef	^	MiscVarType	FVW_MiscVarType	FVW	-	-	-	"MiscVars from the FVW module"	-
 typedef	^	MiscVarType	ReKi            WindPos	{:}{:} - -	"XYZ coordinates to query for wind velocity/acceleration"	-
 typedef	^	MiscVarType	ReKi            WindVel	{:}{:} - -	"XYZ components of wind velocity"	-
 typedef	^	MiscVarType	ReKi            WindAcc	{:}{:} - -	"XYZ components of wind acceleration"	-
-typedef ^   MiscVarType GridInterp_MiscVarType WaveField_m - - - "misc var information from the SeaState WaveField module"   -
+typedef	^	MiscVarType	GridInterp_MiscVarType	WaveField_m	-	-	-	"misc var information from the GridInterp module"	-
 typedef	^	MiscVarType	AD_InflowType	Inflow	{:}	-	-	"Inflow storage (size of u for history of inputs)"	-
 typedef	^	MiscVarType	AD_InputType	u_perturb	-	-	-	"input perturbation for linearization"	-
 typedef	^	MiscVarType	AD_OutputType	y_lin	-	-	-	"output perturbation for linearization"	-

docs/source/user/hydrodyn/input_files.rst

@@ -557,6 +557,7 @@ MODULE AeroDyn_Types
     REAL(ReKi) , DIMENSION(1:3)  :: TFinSTV_i = 0.0_ReKi      !< Structural velocity at the reference point of the fin in the inertial system [-]
     REAL(ReKi) , DIMENSION(1:3)  :: TFinF_i = 0.0_ReKi      !< Forces at the reference point of the fin in the inertial system [-]
     REAL(ReKi) , DIMENSION(1:3)  :: TFinM_i = 0.0_ReKi      !< Moments at the reference point of the fin in the inertial system [-]
+    TYPE(GridInterp_MiscVarType)  :: WaveField_m      !< misc var information from the GridInterp module [-]
   END TYPE RotMiscVarType
 ! =======================
 ! =========  AD_MiscVarType  =======
@@ -568,7 +569,7 @@ MODULE AeroDyn_Types
     REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE  :: WindPos      !< XYZ coordinates to query for wind velocity/acceleration [-]
     REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE  :: WindVel      !< XYZ components of wind velocity [-]
     REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE  :: WindAcc      !< XYZ components of wind acceleration [-]
-    TYPE(GridInterp_MiscVarType)  :: WaveField_m      !< misc var information from the SeaState WaveField module [-]
+    TYPE(GridInterp_MiscVarType)  :: WaveField_m      !< misc var information from the GridInterp module [-]
     TYPE(AD_InflowType) , DIMENSION(:), ALLOCATABLE  :: Inflow      !< Inflow storage (size of u for history of inputs) [-]
     TYPE(AD_InputType)  :: u_perturb      !< input perturbation for linearization [-]
     TYPE(AD_OutputType)  :: y_lin      !< output perturbation for linearization [-]
@@ -5553,6 +5554,9 @@ subroutine AD_CopyRotMiscVarType(SrcRotMiscVarTypeData, DstRotMiscVarTypeData, C
    DstRotMiscVarTypeData%TFinSTV_i = SrcRotMiscVarTypeData%TFinSTV_i
    DstRotMiscVarTypeData%TFinF_i = SrcRotMiscVarTypeData%TFinF_i
    DstRotMiscVarTypeData%TFinM_i = SrcRotMiscVarTypeData%TFinM_i
+   call GridInterp_CopyMisc(SrcRotMiscVarTypeData%WaveField_m, DstRotMiscVarTypeData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2)
+   call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
+   if (ErrStat >= AbortErrLev) return
 end subroutine
 
 subroutine AD_DestroyRotMiscVarType(RotMiscVarTypeData, ErrStat, ErrMsg)
@@ -5776,6 +5780,8 @@ subroutine AD_DestroyRotMiscVarType(RotMiscVarTypeData, ErrStat, ErrMsg)
    call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
    call NWTC_Library_DestroyMeshMapType(RotMiscVarTypeData%T_P_2_T_L, ErrStat2, ErrMsg2)
    call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
+   call GridInterp_DestroyMisc(RotMiscVarTypeData%WaveField_m, ErrStat2, ErrMsg2)
+   call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
 end subroutine
 
 subroutine AD_PackRotMiscVarType(RF, Indata)
@@ -5916,6 +5922,7 @@ subroutine AD_PackRotMiscVarType(RF, Indata)
    call RegPack(RF, InData%TFinSTV_i)
    call RegPack(RF, InData%TFinF_i)
    call RegPack(RF, InData%TFinM_i)
+   call GridInterp_PackMisc(RF, InData%WaveField_m) 
    if (RegCheckErr(RF, RoutineName)) return
 end subroutine
 
@@ -6083,6 +6090,7 @@ subroutine AD_UnPackRotMiscVarType(RF, OutData)
    call RegUnpack(RF, OutData%TFinSTV_i); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%TFinF_i); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%TFinM_i); if (RegCheckErr(RF, RoutineName)) return
+   call GridInterp_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m 
 end subroutine
 
 subroutine AD_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg)

modules/aerodyn/src/AeroDyn.f90

@@ -1489,6 +1489,10 @@ subroutine FVW_CalcOutput(t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg)
    character(*),                    intent(  out)  :: ErrMsg      !< Error message if ErrStat /= ErrID_None
    ! Local variables
    integer(IntKi)                :: nP
+   integer(IntKi)                :: iStepFVW
+   logical, save                 :: FirstNotice_FreeNearWake = .true.
+   logical, save                 :: FirstNotice_NearWake = .true.
+   logical, save                 :: FirstNotice_FullWake = .true.
    integer(IntKi)                :: ErrStat2
    character(ErrMsgLen)          :: ErrMsg2
    character(*), parameter       :: RoutineName = 'FVW_CalcOutput'
@@ -1514,23 +1518,24 @@ subroutine FVW_CalcOutput(t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg)
    endif
 
    ! Write some info to screen when major milestone achieved
-   if (m%iStep == p%nNWFree .and. p%nNWFree<p%nNWMax) then
-      nP = CountCPs(p, p%nNWFree, 0)
-      call WrScr(NewLine//'[INFO] OLAF free near wake is at full extent - '//trim(num2lstr(t))//'s, '//trim(num2lstr(nP))//' points.')
-   endif
-   if (m%iStep == p%nNWMax) then
-      nP = CountCPs(p, p%nNWMax, 0)
-      if (p%nFWMax==0) then
-         call WrScr(NewLine//'[INFO] OLAF wake is at full extent - '//trim(num2lstr(t))//'s, '//trim(num2lstr(nP))//' points.')
-      else
-         call WrScr(NewLine//'[INFO] OLAF near wake is at full extent - '//trim(num2lstr(t))//'s, '//trim(num2lstr(nP))//' points.')
-      endif
-   endif
-   if (p%nFWMax>0 .and. m%iStep== p%nNWMax+p%nFWMax) then
+   iStepFVW = floor(t/p%DTfvw+1.0E-6_DbKi)
+   if (FirstNotice_FullWake .and. iStepFVW == p%nNWMax+p%nFWMax) then
       nP = CountCPs(p, p%nNWMax, p%nFWMax)
       call WrScr(NewLine//'[INFO] OLAF wake is at full extent - '//trim(num2lstr(t))//'s, '//trim(num2lstr(nP))//' points.')
-   endif
-   
+      FirstNotice_FullWake     = .false.
+      FirstNotice_NearWake     = .false.
+      FirstNotice_FreeNearWake = .false.
+   elseif (FirstNotice_NearWake .and. iStepFVW == p%nNWMax) then  ! implies p%nFWMax > 0
+      nP = CountCPs(p, p%nNWMax,  0_IntKi)
+      call WrScr(NewLine//'[INFO] OLAF near wake is at full extent - '//trim(num2lstr(t))//'s, '//trim(num2lstr(nP))//' points.')
+      FirstNotice_NearWake     = .false.
+      FirstNotice_FreeNearWake = .false.
+   elseif (FirstNotice_FreeNearWake .and. iStepFVW == p%nNWFree) then ! implies p%nNWFree<p%nNWMax
+      nP = CountCPs(p, p%nNWFree, 0_IntKi)
+      call WrScr(NewLine//'[INFO] OLAF free near wake is at full extent - '//trim(num2lstr(t))//'s, '//trim(num2lstr(nP))//' points.')
+      FirstNotice_FreeNearWake = .false.
+   end if
+
    ! --- Export to VTK
    if (m%VTKstep==-1) then 
        ! Has never been called, special handling for init