Added flexibility for strip-theory hydrostatic load calculation

docs/source/user/api_change.rst

@@ -76,6 +76,7 @@ ServoDyn                                      14       PitSpr(3)            3.6E
 ServoDyn                                      15       PitDamp(1)           1.4E6         PitDamp(1)  - Blade 1 pitch damping constant
 ServoDyn                                      16       PitDamp(2)           1.4E6         PitDamp(2)  - Blade 2 pitch damping constant
 ServoDyn                                      17       PitDamp(3)           1.4E6         PitDamp(3)  - Blade 3 pitch damping constant *[unused for 2 blades]*
+HydroDyn                                      \*       HstMod               1             HstMod      - Method of computing hydrostatic loads. (0: Up to the still water level. 1: Up to the instantaneous free surface) *[overwrite to 0 when WaveStMod = 0 in SeaState]*
 ============================================= ======== ==================== ========================================================================================================================================================================================================
 
 

docs/source/user/hydrodyn/input_files.rst

@@ -536,6 +536,17 @@ accounted for if **WaveDisp** = 1. **AMMod** should only be set to 0 if wave
 stretching is causing numerical instabilities with flexible fixed-bottom support 
 structures modeled in SubDyn.
 
+**HstMod** controls the computation of distributed hydrostatic loads on 
+strip-theory members. If **HstMod** = 0, the hydrostatic pressure is always 
+integrated on the instantaneous wetted surface up to the SWL, even if wave 
+stretching is enabled. If **HstMod** = 1 and one of the wave stretching model is 
+selected with **WaveStMod** > 0, the hydrostatic pressure will be integrated up 
+to the instantaneous wave free surface, considering both wave elevation and wave 
+slope. If wave stretching is not used, **HstMod** will always be set to zero 
+internally, disregarding the user input. Irrespective of **HstMod**, the exact 
+displaced position of the structure is always used when computing hydrostatic 
+loads to obtain the correct hydrostatic restoring effect.
+
 Axial Coefficients
 ------------------
 This and the next several sections of the input file control the

modules/hydrodyn/src/HydroDyn_Input.f90

@@ -353,6 +353,10 @@ SUBROUTINE HydroDyn_ParseInput( InputFileName, OutRootName, FileInfo_In, InputFi
    call ParseVar( FileInfo_In, CurLine, 'AMMod', InputFileData%Morison%AMMod, ErrStat2, ErrMsg2, UnEc )
       if (Failed())  return;
 
+   ! HstMod - Method of computing strip-theory hydrostatic loads. {0: Up to the still water level. 1: Up to the instantaneous free surface} (switch)
+   call ParseVar( FileInfo_In, CurLine, 'HstMod', InputFileData%Morison%HstMod, ErrStat2, ErrMsg2, UnEc )
+      if (Failed())  return;
+
    !-------------------------------------------------------------------------------------------------
    !  Axial Coefficients Section
    !-------------------------------------------------------------------------------------------------
@@ -1809,7 +1813,10 @@ SUBROUTINE HydroDynInput_ProcessInitData( InitInp, Interval, InputFileData, ErrS
       CALL SetErrStat( ErrID_Fatal,'AMMod must be 0 or 1',ErrStat,ErrMsg,RoutineName)
       RETURN
    END IF
-
+   IF ( InputFileData%Morison%HstMod /= 0 .AND. InputFileData%Morison%HstMod /= 1) THEN
+      CALL SetErrStat( ErrID_Fatal,'HstMod must be 0 or 1',ErrStat,ErrMsg,RoutineName)
+      RETURN
+   END IF  
 
    !-------------------------------------------------------------------------------------------------
    ! Member Joints Section

modules/hydrodyn/src/Morison.f90

@@ -2806,6 +2806,7 @@ SUBROUTINE Morison_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, In
    p%NMOutputs  = InitInp%NMOutputs                       ! Number of members to output [ >=0 and <10]
    p%WaveDisp   = InitInp%WaveDisp
    p%AMMod      = InitInp%AMMod
+   p%HstMod     = InitInp%HstMod
    p%VisMeshes  = InitInp%VisMeshes                       ! visualization mesh for morison elements
    p%PtfmYMod   = InitInp%PtfmYMod
    p%NFillGroups = InitInp%NFillGroups
@@ -2819,6 +2820,10 @@ SUBROUTINE Morison_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, In
        p%AMMod = 0_IntKi
    END IF
 
+   ! Only compute hydrostatic loads up to the wave free surface if waves stretching is enabled
+   IF ( p%WaveField%WaveStMod .EQ. 0_IntKi ) THEN
+       p%HstMod = 0_IntKi
+   END IF
 
    ALLOCATE ( p%MOutLst(p%NMOutputs), STAT = errStat2 )
    IF ( errStat2 /= 0 ) THEN
@@ -3770,7 +3775,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat,
             ! Standard hydrostatic load calculation
             case (1)
 
-               IF ( p%WaveField%WaveStMod > 0_IntKi ) THEN ! If wave stretching is enabled, compute buoyancy up to free surface
+               IF ( p%HstMod > 0_IntKi ) THEN ! If wave stretching is enabled, compute buoyancy up to free surface
                   CALL GetTotalWaveElev(p, m, Time, pos1, Zeta1, ErrStat2, ErrMsg2 )
                   CALL GetTotalWaveElev(p, m, Time, pos2, Zeta2, ErrStat2, ErrMsg2 )
                     CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
@@ -3800,7 +3805,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat,
                ! Get free surface elevation and normal at the element midpoint (both assumed constant over the element)
                posMid = 0.5 * (pos1+pos2)
 
-               IF (p%WaveField%WaveStMod > 0) THEN
+               IF (p%HstMod > 0_IntKi) THEN
                   CALL GetTotalWaveElev(p, m, Time, posMid, ZetaMid, ErrStat2, ErrMsg2 )
                     CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
                   CALL GetFreeSurfaceNormal( p, m, Time, posMid, n_hat, ErrStat2, ErrMsg2 )
@@ -4590,7 +4595,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat,
          end if
          if (mem%i_floor == 0) then  ! both ends above or at seabed
             ! Compute loads on the end plate of node 1
-            IF (p%WaveField%WaveStMod > 0) THEN
+            IF (p%HstMod > 0_IntKi) THEN
                CALL GetTotalWaveElev(p, m, Time, pos1, Zeta1, ErrStat2, ErrMsg2)
                  CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
                CALL GetFreeSurfaceNormal(p, m, Time, pos1, n_hat, ErrStat2, ErrMsg2)
@@ -4618,7 +4623,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat,
             m%F_B_End(:, mem%NodeIndx(  1)) = m%F_B_End(:, mem%NodeIndx(  1)) + F_B_End
 
             ! Compute loads on the end plate of node N+1
-            IF (p%WaveField%WaveStMod > 0) THEN
+            IF (p%HstMod > 0_IntKi) THEN
                CALL GetTotalWaveElev(p, m, Time, pos2, Zeta2, ErrStat2, ErrMsg2)
                  CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
                CALL GetFreeSurfaceNormal(p, m, Time, pos2, n_hat, ErrStat2, ErrMsg2)
@@ -4647,7 +4652,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat,
 
          elseif ( mem%doEndBuoyancy ) then ! The member crosses the seabed line so only the upper end potentially have hydrostatic load
             ! Only compute the loads on the end plate of node N+1
-            IF (p%WaveField%WaveStMod > 0) THEN
+            IF (p%HstMod > 0_IntKi) THEN
                CALL GetTotalWaveElev(p, m, Time, pos2, Zeta2, ErrStat2, ErrMsg2)
                  CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
                CALL GetFreeSurfaceNormal(p, m, Time, pos2, n_hat, ErrStat2, ErrMsg2)
@@ -5934,7 +5939,7 @@ SUBROUTINE getElementHstLds_Mod1(p, m, mem, Time, pos1, pos2, Zeta1, Zeta2, k_ha
          ! Scaled radius of element at point where its centerline crosses the waterplane
          rh       = r1 + h0*dRdl
          ! Estimate the free-surface normal at the free-surface intersection, n_hat
-         IF ( p%WaveField%WaveStMod > 0_IntKi ) THEN ! If wave stretching is enabled, compute free surface normal
+         IF (p%HstMod > 0_IntKi) THEN ! If wave stretching is enabled, compute free surface normal
             CALL GetFreeSurfaceNormal(p, m, Time, FSInt, n_hat, ErrStat2, ErrMsg2 )
               CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
          ELSE ! Without wave stretching, use the normal of the SWL

modules/hydrodyn/src/Morison.txt

@@ -337,7 +337,8 @@ typedef   ^                            ^                             INTEGER
 #
 typedef   ^                            InitInputType                 ReKi                     Gravity                          -          -        -         "Gravity (scalar, positive-valued)"    m/s^2
 typedef   ^                            ^                             INTEGER                  WaveDisp                         -          -        -         "Method of computing Wave Kinematics. (0: use undisplaced position, 1: use displaced position, 2: use low-pass filtered displaced position) " -
-typedef   ^                            ^                             INTEGER                  AMMod                            -          -        -         "Method of computing distributed added-mass force. (0: Only and always on nodes below SWL at the undisplaced position. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveMod = 0 or 6 or when WaveStMod = 0 in SeaState]" -
+typedef   ^                            ^                             INTEGER                  AMMod                            -          -        -         "Method of computing distributed added-mass force. (0: Only and always on nodes below SWL at the undisplaced position. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveStMod = 0 in SeaState]" -
+typedef   ^                            ^                             INTEGER                  HstMod                           -          -        -         "Method of computing strip-theory hydrostatic loads. (0: Up to the still water level. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveStMod = 0 in SeaState]" -
 typedef   ^                            ^                             INTEGER                  NJoints                          -          -        -         "Number of user-specified joints"    -
 typedef   ^                            ^                             INTEGER                  NNodes                           -          -        -         "Total number of nodes in the final software model"    -
 typedef   ^                            ^                             Morison_JointType        InpJoints                        {:}        -        -         "Array of user-specified joints"    -
@@ -472,6 +473,7 @@ typedef   ^                            ParameterType                 DbKi
 typedef   ^                            ^                             ReKi                     Gravity                       -          -        -         "Gravity (scalar, positive-valued)"    m/s^2
 typedef   ^                            ^                             INTEGER                  WaveDisp                      -          -        -         "Method of computing Wave Kinematics. (0: use undisplaced position, 1: use displaced position, 2: use low-pass filtered displaced position) " -
 typedef   ^                            ^                             INTEGER                  AMMod                         -          -        -         "Method of computing distributed added-mass force. (0: Only and always on nodes below SWL at the undisplaced position. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveMod = 0 or 6 or when WaveStMod = 0 in SeaState]" -
+typedef   ^                            ^                             INTEGER                  HstMod                           -          -        -         "Method of computing strip-theory hydrostatic loads. (0: Up to the still water level. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveStMod = 0 in SeaState]" -
 typedef   ^                            ^                             INTEGER                  NMembers                      -          -        -         "number of members"        -
 typedef   ^                            ^                             Morison_MemberType       Members                       {:}        -        -         "Array of Morison members used during simulation"    -
 typedef   ^                            ^                             INTEGER                  NNodes                        -          -        -         ""        -

modules/hydrodyn/src/Morison_Types.f90

@@ -405,7 +405,8 @@ MODULE Morison_Types
   TYPE, PUBLIC :: Morison_InitInputType
     REAL(ReKi)  :: Gravity = 0.0_ReKi      !< Gravity (scalar, positive-valued) [m/s^2]
     INTEGER(IntKi)  :: WaveDisp = 0_IntKi      !< Method of computing Wave Kinematics. (0: use undisplaced position, 1: use displaced position, 2: use low-pass filtered displaced position)  [-]
-    INTEGER(IntKi)  :: AMMod = 0_IntKi      !< Method of computing distributed added-mass force. (0: Only and always on nodes below SWL at the undisplaced position. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveMod = 0 or 6 or when WaveStMod = 0 in SeaState] [-]
+    INTEGER(IntKi)  :: AMMod = 0_IntKi      !< Method of computing distributed added-mass force. (0: Only and always on nodes below SWL at the undisplaced position. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveStMod = 0 in SeaState] [-]
+    INTEGER(IntKi)  :: HstMod = 0_IntKi      !< Method of computing strip-theory hydrostatic loads. (0: Up to the still water level. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveStMod = 0 in SeaState] [-]
     INTEGER(IntKi)  :: NJoints = 0_IntKi      !< Number of user-specified joints [-]
     INTEGER(IntKi)  :: NNodes = 0_IntKi      !< Total number of nodes in the final software model [-]
     TYPE(Morison_JointType) , DIMENSION(:), ALLOCATABLE  :: InpJoints      !< Array of user-specified joints [-]
@@ -539,6 +540,7 @@ MODULE Morison_Types
     REAL(ReKi)  :: Gravity = 0.0_ReKi      !< Gravity (scalar, positive-valued) [m/s^2]
     INTEGER(IntKi)  :: WaveDisp = 0_IntKi      !< Method of computing Wave Kinematics. (0: use undisplaced position, 1: use displaced position, 2: use low-pass filtered displaced position)  [-]
     INTEGER(IntKi)  :: AMMod = 0_IntKi      !< Method of computing distributed added-mass force. (0: Only and always on nodes below SWL at the undisplaced position. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveMod = 0 or 6 or when WaveStMod = 0 in SeaState] [-]
+    INTEGER(IntKi)  :: HstMod = 0_IntKi      !< Method of computing strip-theory hydrostatic loads. (0: Up to the still water level. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveStMod = 0 in SeaState] [-]
     INTEGER(IntKi)  :: NMembers = 0_IntKi      !< number of members [-]
     TYPE(Morison_MemberType) , DIMENSION(:), ALLOCATABLE  :: Members      !< Array of Morison members used during simulation [-]
     INTEGER(IntKi)  :: NNodes = 0_IntKi      !<  [-]
@@ -3240,6 +3242,7 @@ subroutine Morison_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode, E
    DstInitInputData%Gravity = SrcInitInputData%Gravity
    DstInitInputData%WaveDisp = SrcInitInputData%WaveDisp
    DstInitInputData%AMMod = SrcInitInputData%AMMod
+   DstInitInputData%HstMod = SrcInitInputData%HstMod
    DstInitInputData%NJoints = SrcInitInputData%NJoints
    DstInitInputData%NNodes = SrcInitInputData%NNodes
    if (allocated(SrcInitInputData%InpJoints)) then
@@ -3687,6 +3690,7 @@ subroutine Morison_PackInitInput(RF, Indata)
    call RegPack(RF, InData%Gravity)
    call RegPack(RF, InData%WaveDisp)
    call RegPack(RF, InData%AMMod)
+   call RegPack(RF, InData%HstMod)
    call RegPack(RF, InData%NJoints)
    call RegPack(RF, InData%NNodes)
    call RegPack(RF, allocated(InData%InpJoints))
@@ -3892,6 +3896,7 @@ subroutine Morison_UnPackInitInput(RF, OutData)
    call RegUnpack(RF, OutData%Gravity); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%WaveDisp); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%AMMod); if (RegCheckErr(RF, RoutineName)) return
+   call RegUnpack(RF, OutData%HstMod); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%NJoints); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%NNodes); if (RegCheckErr(RF, RoutineName)) return
    if (allocated(OutData%InpJoints)) deallocate(OutData%InpJoints)
@@ -4895,6 +4900,7 @@ subroutine Morison_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrM
    DstParamData%Gravity = SrcParamData%Gravity
    DstParamData%WaveDisp = SrcParamData%WaveDisp
    DstParamData%AMMod = SrcParamData%AMMod
+   DstParamData%HstMod = SrcParamData%HstMod
    DstParamData%NMembers = SrcParamData%NMembers
    if (allocated(SrcParamData%Members)) then
       LB(1:1) = lbound(SrcParamData%Members)
@@ -5208,6 +5214,7 @@ subroutine Morison_PackParam(RF, Indata)
    call RegPack(RF, InData%Gravity)
    call RegPack(RF, InData%WaveDisp)
    call RegPack(RF, InData%AMMod)
+   call RegPack(RF, InData%HstMod)
    call RegPack(RF, InData%NMembers)
    call RegPack(RF, allocated(InData%Members))
    if (allocated(InData%Members)) then
@@ -5297,6 +5304,7 @@ subroutine Morison_UnPackParam(RF, OutData)
    call RegUnpack(RF, OutData%Gravity); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%WaveDisp); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%AMMod); if (RegCheckErr(RF, RoutineName)) return
+   call RegUnpack(RF, OutData%HstMod); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%NMembers); if (RegCheckErr(RF, RoutineName)) return
    if (allocated(OutData%Members)) deallocate(OutData%Members)
    call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return

openfast_io/openfast_io/FAST_reader.py

@@ -1966,6 +1966,7 @@ def read_HydroDyn(self, hd_file):
         f.readline()
         self.fst_vt['HydroDyn']['WaveDisp']     = int_read(f.readline().split()[0])
         self.fst_vt['HydroDyn']['AMMod']        = int_read(f.readline().split()[0])
+        self.fst_vt['HydroDyn']['HstMod']       = int_read(f.readline().split()[0])
 
         #AXIAL COEFFICIENTS
         f.readline()

openfast_io/openfast_io/FAST_writer.py

@@ -1699,6 +1699,7 @@ def write_HydroDyn(self):
         f.write('---------------------- STRIP THEORY OPTIONS --------------------------------------\n')
         f.write('{:<22d} {:<11} {:}'.format(self.fst_vt['HydroDyn']['WaveDisp'], 'WaveDisp', '- Method of computing Wave Kinematics {0: use undisplaced position, 1: use displaced position) } (switch)\n'))
         f.write('{:<22d} {:<11} {:}'.format(self.fst_vt['HydroDyn']['AMMod'], 'AMMod', '- Method of computing distributed added-mass force. (0: Only and always on nodes below SWL at the undisplaced position. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveStMod = 0 in SeaState]]\n'))
+        f.write('{:<22d} {:<11} {:}'.format(self.fst_vt['HydroDyn']['HstMod'], 'HstMod', '- Method of computing hydrostatic loads. (0: Up to the still water level. 1: Up to the instantaneous free surface) [overwrite to 0 when WaveStMod = 0 in SeaState]]\n'))
 
         f.write('---------------------- AXIAL COEFFICIENTS --------------------------------------\n')
         f.write('{:<22d} {:<11} {:}'.format(self.fst_vt['HydroDyn']['NAxCoef'], 'NAxCoef', '- Number of axial coefficients (-)\n'))