diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index 4dfe7e438e..e95ac590af 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -3503,7 +3503,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, REAL(ReKi) :: sinBeta, sinBeta1, sinBeta2 REAL(ReKi) :: cosBeta, cosBeta1, cosBeta2 REAL(ReKi) :: CMatrix(3,3), CMatrix1(3,3), CMatrix2(3,3), CTrans(3,3) ! Direction cosine matrix for element, and its transpose - REAL(ReKi) :: l, z1, z2, zMid, r1, r2, r1b, r2b, r1In, r2In, rMidIn, z_hi, zFillGroup + REAL(ReKi) :: l, z1, z2, zMid, r1, r2, r1b, r2b, r1In, r2In, rMidIn, rn, z_hi, zFillGroup REAL(ReKi) :: Sa1, Sa2, Sa1b, Sa2b, SaMidb, Sa1In, Sa2In, SaMidIn REAL(ReKi) :: Sb1, Sb2, Sb1b, Sb2b, SbMidb, Sb1In, Sb2In, SbMidIn REAL(ReKi) :: dRdl_mg, dSadl_mg, dSbdl_mg ! shorthand for taper including marine growth of element i @@ -3822,7 +3822,12 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, 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 ) + IF (mem%MSecGeom == MSecGeom_Cyl) THEN + rn = 0.5 * (r1b + r2b) + ELSE + rn = 0.25 * MAX(Sa1b+Sa2b, Sb1b+Sb2b) + END IF + CALL GetFreeSurfaceNormal( p, m, Time, posMid, rn, n_hat, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/posMid(1),posMid(2),ZetaMid/) ! Reference point on the free surface ELSE @@ -4623,13 +4628,15 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, Sa2 = mem%SaMGB(N+1) Sb1 = mem%SbMGB( 1) Sb2 = mem%SbMGB(N+1) + r1 = 0.5 * MAX(Sa1,Sb1) ! Only used for computing the free-surface normal with finite differencing + r2 = 0.5 * MAX(Sa2,Sb2) ! Only used for computing the free-surface normal with finite differencing 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%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) + CALL GetFreeSurfaceNormal(p, m, Time, pos1, r1, n_hat, ErrStat2, ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/pos1(1),pos1(2),Zeta1/) ! Reference point on the free surface ELSE @@ -4657,7 +4664,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, 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) + CALL GetFreeSurfaceNormal(p, m, Time, pos2, r2, n_hat, ErrStat2, ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/pos2(1),pos2(2),Zeta2/) ! Reference point on the free surface ELSE @@ -4686,7 +4693,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, 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) + CALL GetFreeSurfaceNormal(p, m, Time, pos2, r2, n_hat, ErrStat2, ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/pos2(1),pos2(2),Zeta2/) ! Reference point on the free surface ELSE @@ -4835,11 +4842,12 @@ SUBROUTINE GetTotalWaveElev(p, m, Time, pos, Zeta, ErrStat, ErrMsg ) END SUBROUTINE GetTotalWaveElev - SUBROUTINE GetFreeSurfaceNormal(p, m, Time, pos, n, ErrStat, ErrMsg) + SUBROUTINE GetFreeSurfaceNormal(p, m, Time, pos, r, n, ErrStat, ErrMsg) TYPE(Morison_ParameterType), INTENT( IN ) :: p TYPE(Morison_MiscVarType), INTENT( INOUT ) :: m REAL(DbKi), INTENT( In ) :: Time REAL(ReKi), INTENT( In ) :: pos(:) ! Position at which free-surface normal is to be calculated. Third entry ignored if present. + REAL(ReKi), INTENT( In ) :: r ! Characteristic section half-width used to filter the free-surface slope REAL(ReKi), INTENT( OUT ) :: n(3) ! Free-surface normal vector INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None @@ -4849,7 +4857,7 @@ SUBROUTINE GetFreeSurfaceNormal(p, m, Time, pos, n, ErrStat, ErrMsg) ErrStat = ErrID_None ErrMsg = "" - CALL WaveField_GetNodeWaveNormal( p%WaveField, m%WaveField_m, Time, pos, n, ErrStat2, ErrMsg2 ) + CALL WaveField_GetNodeWaveNormal( p%WaveField, m%WaveField_m, Time, pos, n, ErrStat2, ErrMsg2, r ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END SUBROUTINE GetFreeSurfaceNormal @@ -5974,7 +5982,7 @@ SUBROUTINE getElementHstLds_Mod1(p, m, mem, Time, pos1, pos2, Zeta1, Zeta2, k_ha rh = r1 + h0*dRdl ! Estimate the free-surface normal at the free-surface intersection, n_hat 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 GetFreeSurfaceNormal(p, m, Time, FSInt, rh, n_hat, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ELSE ! Without wave stretching, use the normal of the SWL n_hat = (/0.0_ReKi,0.0_ReKi,1.0_ReKi/) diff --git a/modules/seastate/src/SeaSt_WaveField.f90 b/modules/seastate/src/SeaSt_WaveField.f90 index a01e2bd89e..d73f27cac3 100644 --- a/modules/seastate/src/SeaSt_WaveField.f90 +++ b/modules/seastate/src/SeaSt_WaveField.f90 @@ -43,7 +43,7 @@ FUNCTION WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, pos, ErrS ErrStat = ErrID_None IF (ALLOCATED(WaveField%WaveElev1) .or. ALLOCATED(WaveField%WaveElev2)) then - CALL WaveField_Interp_Setup3D(Time, pos, WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2) + CALL WaveField_Interp_Setup3D(Time+WaveField%WaveTimeShift, pos, WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2) CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return end if @@ -99,7 +99,7 @@ subroutine WaveField_GetMinMaxWaveElevEstimate( WaveField, MinElev, MaxElev, Err endif end subroutine WaveField_GetMinMaxWaveElevEstimate -SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, n, ErrStat, ErrMsg ) +SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, n, ErrStat, ErrMsg, r ) type(SeaSt_WaveFieldType), intent(in ) :: WaveField type(GridInterp_MiscVarType), intent(inout) :: WaveField_m real(DbKi), intent(in ) :: Time @@ -107,15 +107,37 @@ SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, n, Er real(ReKi), intent( out) :: n(3) ! Free-surface normal vector integer(IntKi), intent( out) :: ErrStat ! Error status of the operation character(*), intent( out) :: ErrMsg ! Error message if errStat /= ErrID_None + real(ReKi), optional, intent(in ) :: r ! Section half-width; if > 0, slope is central-differenced over +/- r (filters waves shorter than the section) real(SiKi) :: slope(2) + real(SiKi) :: ZetaP, ZetaM + real(DbKi) :: TimeShifted character(*), parameter :: RoutineName = 'WaveField_GetNodeWaveNormal' integer(IntKi) :: errStat2 character(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None + ErrMsg = "" + + if (present(r)) then + if (r > 0.0_ReKi) then + ZetaP = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1)+r,pos(2)/), ErrStat2, ErrMsg2 ); if (Failed()) return; + ZetaM = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1)-r,pos(2)/), ErrStat2, ErrMsg2 ); if (Failed()) return; + slope(1) = (ZetaP-ZetaM) / real(2.0_ReKi*r,SiKi) + ZetaP = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1),pos(2)+r/), ErrStat2, ErrMsg2 ); if (Failed()) return; + ZetaM = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1),pos(2)-r/), ErrStat2, ErrMsg2 ); if (Failed()) return; + slope(2) = (ZetaP-ZetaM) / real(2.0_ReKi*r,SiKi) + n = Real( (/-slope(1),-slope(2),1.0_SiKi/), ReKi) + n = n / TwoNorm(n) + return + end if + end if - call GridInterpSetupN( (/Real(Time+WaveField%WaveTimeShift,ReKi),pos(1),pos(2)/), WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2 ) + n = (/0.0_ReKi,0.0_ReKi,1.0_ReKi/) + if (.not. ALLOCATED(WaveField%WaveElev1)) return + + TimeShifted = Time + WaveField%WaveTimeShift + call GridInterpSetupN( (/Real(TimeShifted,ReKi),pos(1),pos(2)/), WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; slope = GridInterpS( WaveField%WaveElev1, WaveField%SrfGridParams, WaveField_m ) if (ALLOCATED(WaveField%WaveElev2)) then slope = slope + GridInterpS( WaveField%WaveElev2, WaveField%SrfGridParams, WaveField_m ) @@ -460,7 +482,7 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod ELSE ! Node is above SWL - need wave stretching ! Vertical wave stretching - CALL WaveField_Interp_Setup4D( Time, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + CALL WaveField_Interp_Setup4D( Time+WaveField%WaveTimeShift, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m ) ! Extrapolated wave stretching @@ -585,7 +607,7 @@ SUBROUTINE WaveField_GetNodeWaveVelAcc( WaveField, WaveField_m, Time, pos, force ELSE ! Node is above SWL - need wave stretching ! Vertical wave stretching - CALL WaveField_Interp_Setup4D( Time, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + CALL WaveField_Interp_Setup4D( Time+WaveField%WaveTimeShift, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m ) FA(:) = GridInterp4DVec( WaveField%WaveAcc, WaveField_m )