diff --git a/docs/source/user/api_change.rst b/docs/source/user/api_change.rst index 97040f95bd..9c54b99291 100644 --- a/docs/source/user/api_change.rst +++ b/docs/source/user/api_change.rst @@ -41,6 +41,38 @@ AeroDyn \* (-) OLAF 26 RegFunctionPart 2 RegFunctionPart - Particle regularization function {0: None, 1: Exponential, 2: Compact, "default": 1} [only if VelocityMethod=2,3] (switch) ============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== +The meaning of the HydroDyn ``WaveDisp`` input has changed, although its valid values (0 and 1) are unchanged. +``WaveDisp`` now controls both where the wave kinematics and dynamic pressure are evaluated and the orientation +of the strip-theory (Morison) members and joints. With ``WaveDisp=0``, the strip-theory members are placed at a +reference configuration consistent with the potential-flow model—the reference yaw (from ``PtfmYMod``/``PtfmRefY``) +plus a horizontal x-y drift of the platform reference point governed by ``ExctnDisp`` and ``ExctnCutOff``—rather +than being held at the undisplaced position. With ``WaveDisp=1``, the exact instantaneous displaced position and +orientation are used. The previous ``WaveDisp=0`` behavior is recovered by also setting ``PtfmYMod=0`` and +``ExctnDisp=0``. + +============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== +Modified in OpenFAST `5.1.0` +-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +Module Line Flag Name Example Value +============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== +HydroDyn \* WaveDisp 0 WaveDisp - Method of displacing strip-theory members and joints {0: use potential-flow consistent displacement, 1: use exact instantaneous displacement} (switch) +============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== + +In previous versions of OpenFAST, HydroDyn's input file included the ``AMMod`` option to control the method of computing +distributed strip-theory added-mass force. Setting ``AMMod`` to 0 forces HydroDyn to always compute the added-mass +force up to the still water level (SWL). This switch was added because of numerical stability issues encountered when +evaluating the added-mass force up to the instantaneous free surface with some hydroelastic models. With the improved +solver stability of OpenFAST `5.1.0`, the ``AMMod`` option was deemed no longer necessary based on user feedback and was +removed. + +============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== +Removed in OpenFAST `5.1.0` +-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +Module Line Flag Name Example Value +============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== +HydroDyn \* AMMod 0 AMMod - Method of computing distributed strip-theory added-mass force {0: evaluate up to SWL, 1: evaluate up to instantaneous free surface if WaveStMod > 0} (switch) +============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== + OpenFAST v4.2.x to OpenFAST v5.0.0 ----------------------------------- diff --git a/docs/source/user/hydrodyn/appendix.rst b/docs/source/user/hydrodyn/appendix.rst index 9df2179bc4..521a590c3e 100644 --- a/docs/source/user/hydrodyn/appendix.rst +++ b/docs/source/user/hydrodyn/appendix.rst @@ -67,8 +67,7 @@ structure:: 0 0 0 0 0 0 0 0 0 0 0 0 ---------------------- STRIP THEORY OPTIONS -------------------------------------- - 0 WaveDisp - Method of computing Wave Kinematics {0: use undisplaced position, 1: use displaced position} (switch) - 0 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} (switch) [overwrite to 0 when WaveStMod = 0 in SeaState] + 0 WaveDisp - Method of displacing strip-theory members and joints {0: use potential-flow consistent displacement, 1: use exact instantaneous displacement} (switch) 0 HstMod - Method of computing hydrostatic loads. {0: Up to the still water level. 1: Up to the instantaneous free surface} (switch) [overwrite to 0 when WaveStMod = 0 in SeaState] ---------------------- AXIAL COEFFICIENTS -------------------------------------- 2 NAxCoef - Number of axial coefficients (-) diff --git a/docs/source/user/hydrodyn/input_files.rst b/docs/source/user/hydrodyn/input_files.rst index 36e2e6825f..848dfcc77b 100644 --- a/docs/source/user/hydrodyn/input_files.rst +++ b/docs/source/user/hydrodyn/input_files.rst @@ -244,6 +244,14 @@ motion to prevent double counting the contributions from first-order structural motion already included in the second-order potential-flow wave excitation. +Beyond the potential-flow wave excitation, **ExctnDisp** and **ExctnCutOff** also +govern the horizontal x-y drift applied to the strip-theory (Morison) members +when **WaveDisp** = 0, as described under Strip theory options below. These inputs +are honored regardless of whether potential-flow bodies are present, so a +strip-theory-only model can use the same reference-position treatment. When +**ExctnDisp** = 2, **ExctnCutOff** must be greater than zero even without +potential-flow bodies. + HydroDyn now supports large but slow (well below wave frequencies) transient platform yaw motion with both strip-theory only and hybrid potential-flow models. To enable this capability, the inputs @@ -313,13 +321,15 @@ both **DiffQTF** and **SumQTF** to be set to 0. However, mean- or slow-drift loads based on Newman's approximation can be included through the **MnDrift** or **NewmanApp** inputs explained below. -Note that the inputs **PtfmYMod** and **PtfmRefY** also affect the -strip-theory hydrodynamic load. This is because the orientation of -the strip-theory members is updated based on **PtfmRefY** instead -of the instantaneous platform yaw rotation. Behavior of previous -versions of HydroDyn can be approximately recovered by setting -**PtfmYMod** = 0 and **PtfmRefY** = 0 deg, in which case, the -inputs **PtfmYCutoff** and **NExctnHdg** are not used. +Note that **PtfmYMod** and **PtfmRefY** also affect the strip-theory +hydrodynamic load when **WaveDisp** = 0. In that mode, the displacement +and rotation of the strip-theory members and joints are based on +**PtfmRefY** (together with the x-y drift from **ExctnDisp**) rather than +the full 6DOF instantaneous structure motion, as described under Strip +theory options below. When **WaveDisp** = 1, the exact instantaneous +member and joint orientations are used instead, so **PtfmRefY** does not +affect the strip-theory load. The inputs **PtfmYCutoff** and **NExctnHdg** +are not used when **PtfmYMod** = 0. HydroDyn has two methods for calculating the radiation memory effect. Set **RdtnMod** to 1 for the convolution method, 2 for the linear @@ -707,29 +717,29 @@ apply constant hydrostatic/buoyancy loads to the generalized modes. Strip theory options -------------------- -**WaveDisp** can be set to 0 to compute the strip-theory loads using the -wave kinematics and dynamic pressure at the undisplaced position of the -structure. If set to 1, the loads will be computed using the wave kinematics -and dynamic pressure at the instantaneous displaced positions of the strip-theory -members. Note that when wave stretching is not used (\ **WaveStMod** = 0 in -SeaState), only the *X*- and *Y*-displacements of the strip-theory member -nodes are considered when **WaveDisp** = 1, while the vertical *Z*-displacement is -ignored. This is done to avoid discontinuous nodal loads that can result in -unphysical structural vibration with a SubDyn substructure model. When -**WaveStMod** > 0 and **WaveDisp** = 1, displacements of strip-theory members -in all three directions are considered when computing the wave kinematics. -A load smoothing procedure is performed to avoid discontinuous nodal loads -in this case. - -**AMMod** controls the computation of distributed strip-theory added-mass force. -If **AMMod** = 0, the strip-theory added-mass force is always evaluated up -to the SWL while neglecting the vertical displacement of the strip-theory member -nodes, even if wave stretching is enabled. With **AMMod** = 1, the strip-theory -added-mass force is evaluated up to the instantaneous free surface if -**WaveStMod** > 0. The vertical displacement of strip-theory members will also be -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. +**WaveDisp** controls both the position at which the wave kinematics and dynamic +pressure are evaluated and the orientation of the strip-theory members and joints +used when computing the hydrodynamic loads. + +Setting **WaveDisp** = 0 evaluates the loads at a reference configuration of the +strip-theory members that is consistent with the potential-flow model. Rather +than holding the members at their undisplaced position and orientation, HydroDyn +positions and orients them using the reference yaw (set through **PtfmYMod** and +**PtfmRefY**) together with a horizontal x-y drift of the HydroDyn origin (PRP) +governed by **ExctnDisp** and **ExctnCutOff** (see above). This results in 3DOF +horizontal-plane displacements only (surge, sway, and yaw). Both the member and +joint orientations used to compute the hydrodynamic loads and the node positions +at which the wave kinematics and dynamic pressure are queried follow this +reference configuration, so a mixed potential-flow/strip-theory model remains +geometrically self-consistent. The behavior of previous versions of HydroDyn with +**WaveDisp** = 0 can be recovered by setting **PtfmYMod** = 0 and **ExctnDisp** = 0. + +With **WaveDisp** = 1, the strip-theory loads are computed considering the full +6DOF instantaneous motion of the strip-theory members and joints. The wave +kinematics and dynamic pressure are evaluated at the exact instantaneous positions +of the strip-theory member nodes and joints, and member and joint orientations +consider full 3DOF rotation of the structure. A load smoothing procedure is +performed to avoid discontinuous nodal loads as members cross the free surface. **HstMod** controls the computation of distributed hydrostatic loads on strip-theory members. If **HstMod** = 0, the hydrostatic pressure is always diff --git a/modules/hydrodyn/src/HydroDyn.f90 b/modules/hydrodyn/src/HydroDyn.f90 index b956880bc8..a932887b6a 100644 --- a/modules/hydrodyn/src/HydroDyn.f90 +++ b/modules/hydrodyn/src/HydroDyn.f90 @@ -184,6 +184,11 @@ SUBROUTINE HydroDyn_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, I p%WaveField => InitInp%WaveField p%PtfmYMod = InputFileData%PtfmYMod + ! Capture the user-specified ExctnDisp/ExctnCutOff that drive the Morison WaveDisp=0 horizontal drift (PtfmRefXY) here, + ! before HydroDynInput_ProcessInitData may force ExctnDisp=0 for the WAMIT excitation grid. This lets the Morison drift + ! follow the user settings regardless of whether potential-flow bodies are present (0: none, 1: instantaneous, 2: filtered). + p%PtfmXYMod = InputFileData%WAMIT%ExctnDisp + p%CXYFilt = exp(-TwoPi*Interval*InputFileData%WAMIT%ExctnCutOff) InputFileData%Morison%WaveField => InitInp%WaveField InputFileData%WAMIT%WaveField => InitInp%WaveField @@ -248,6 +253,10 @@ SUBROUTINE HydroDyn_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, I END IF xd%PtfmRefY = InputFileData%PtfmRefY + ! Seed the filtered Morison WaveDisp=0 drift from the initial platform offset (matches WAMIT BdyPosFilt) to avoid a startup transient + xd%PtfmRefXY(1,:) = InitInp%PlatformPos(1) + xd%PtfmRefXY(2,:) = InitInp%PlatformPos(2) + ! Open a summary of the HydroDyn Initialization. Note: OutRootName must be set by the caller because there may not be an input file to obtain this rootname from. IF ( InputFileData%HDSum ) THEN @@ -636,6 +645,7 @@ SUBROUTINE HydroDyn_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, I InputFileData%Morison%VisMeshes = p%VisMeshes ! Additional Morison inputs to be initialized just in case u%Morison%PtfmRefY = 0.0_ReKi + u%Morison%PtfmRefXY = 0.0_ReKi u%Morison%PRP = [0.0_ReKi,0.0_ReKi,0.0_ReKi] ! Initialize the Morison Element Calculations CALL Morison_Init(InputFileData%Morison, u%Morison, p%Morison, x%Morison, xd%Morison, z%Morison, OtherState%Morison, & @@ -1223,24 +1233,32 @@ SUBROUTINE HydroDyn_UpdateStates( t, n, Inputs, InputTimes, p, x, xd, z, OtherSt nTime = size(Inputs) ! Update PtfmRefY - IF (p%PtfmYMod .EQ. 1) THEN + IF (p%PtfmYMod .EQ. 1 .OR. p%PtfmXYMod .EQ. 2) THEN ! Inefficient. Only need to interp PRPMesh below. Fix later. CALL HydroDyn_CopyInput(Inputs(1), u, MESH_NEWCOPY, ErrStat2, ErrMsg2) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL HydroDyn_Input_ExtrapInterp(Inputs, InputTimes, u, t, ErrStat2, ErrMsg2) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - PRPRotation = EulerExtractZYX(u%PRPMesh%Orientation(:,:,1)) - ! Yaw angle from EulerExtractZYX might not be continous in time and can contain jumps of TwoPi - ! Adjust past xd%PtfmRefY to follow. - IF ( ABS(PRPRotation(3)-xd%PtfmRefY(1)) > ABS(PRPRotation(3)-(xd%PtfmRefY(1)-TwoPi)) ) THEN - xd%PtfmRefY = xd%PtfmRefY - TwoPi - ELSE IF ( ABS(PRPRotation(3)-xd%PtfmRefY(1)) > ABS(PRPRotation(3)-(xd%PtfmRefY(1)+TwoPi)) ) THEN - xd%PtfmRefY = xd%PtfmRefY + TwoPi + IF (p%PtfmYMod .EQ. 1) THEN + PRPRotation = EulerExtractZYX(u%PRPMesh%Orientation(:,:,1)) + ! Yaw angle from EulerExtractZYX might not be continous in time and can contain jumps of TwoPi + ! Adjust past xd%PtfmRefY to follow. + IF ( ABS(PRPRotation(3)-xd%PtfmRefY(1)) > ABS(PRPRotation(3)-(xd%PtfmRefY(1)-TwoPi)) ) THEN + xd%PtfmRefY = xd%PtfmRefY - TwoPi + ELSE IF ( ABS(PRPRotation(3)-xd%PtfmRefY(1)) > ABS(PRPRotation(3)-(xd%PtfmRefY(1)+TwoPi)) ) THEN + xd%PtfmRefY = xd%PtfmRefY + TwoPi + END IF + ! Update PtfmRefY states + xd%PtfmRefY(3) = xd%PtfmRefY(2) + xd%PtfmRefY(2) = xd%PtfmRefY(1) + xd%PtfmRefY(1) = p%CYawFilt * xd%PtfmRefY(1) + (1.0-p%CYawFilt) * PRPRotation(3) + END IF + IF (p%PtfmXYMod .EQ. 2) THEN + ! Low-pass filter the PRP x,y displacement (parallel to PtfmRefY) + xd%PtfmRefXY(:,3) = xd%PtfmRefXY(:,2) + xd%PtfmRefXY(:,2) = xd%PtfmRefXY(:,1) + xd%PtfmRefXY(:,1) = p%CXYFilt * xd%PtfmRefXY(:,1) + (1.0-p%CXYFilt) * u%PRPMesh%TranslationDisp(1:2,1) END IF - ! Update PtfmRefY states - xd%PtfmRefY(3) = xd%PtfmRefY(2) - xd%PtfmRefY(2) = xd%PtfmRefY(1) - xd%PtfmRefY(1) = p%CYawFilt * xd%PtfmRefY(1) + (1.0-p%CYawFilt) * PRPRotation(3) CALL HydroDyn_DestroyInput(u, ErrStat2, ErrMsg2) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END IF @@ -1257,6 +1275,14 @@ SUBROUTINE HydroDyn_UpdateStates( t, n, Inputs, InputTimes, p, x, xd, z, OtherSt CALL Morison_CopyInput(Inputs(i)%Morison, Inputs_Morison(i), MESH_NEWCOPY, ErrStat2, ErrMsg2) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Inputs_Morison(i)%PtfmRefY = xd%PtfmRefY(i) + select case (p%PtfmXYMod) + case (1) ! instantaneous PRP x,y drift + Inputs_Morison(i)%PtfmRefXY = Inputs(i)%PRPMesh%TranslationDisp(1:2,1) + case (2) ! low-pass filtered PRP x,y drift + Inputs_Morison(i)%PtfmRefXY = xd%PtfmRefXY(:,i) + case default ! 0: no drift + Inputs_Morison(i)%PtfmRefXY = 0.0_ReKi + end select Inputs_Morison(i)%PRP = Inputs(i)%PRPMesh%Position(:,1) + Inputs(i)%PRPMesh%TranslationDisp(:,1) END DO CALL Morison_CopyInput(Inputs(1)%Morison, u_Morison, MESH_NEWCOPY, ErrStat2, ErrMsg2) @@ -1491,6 +1517,7 @@ SUBROUTINE HydroDyn_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, integer(IntKi) :: iBody, indxStart, indxEnd, AddDOFCntr ! Counters REAL(ReKi), ALLOCATABLE :: RRg2b(:,:), RRb2g(:,:) REAL(ReKi) :: PtfmRefY + REAL(ReKi) :: PtfmRefXY(2) REAL(R8Ki) :: PRPRotation(3) REAL(ReKi) :: NLFKForce(3,p%NBody) @@ -1547,6 +1574,16 @@ SUBROUTINE HydroDyn_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, if (Failed()) return END IF + ! Fresh x,y drift of the HD origin for WaveDisp=0 Morison (parallel to PtfmRefY) + SELECT CASE (p%PtfmXYMod) + CASE (1) ! instantaneous PRP x,y drift + PtfmRefXY = u%PRPMesh%TranslationDisp(1:2,1) + CASE (2) ! low-pass filtered PRP x,y drift + PtfmRefXY = p%CXYFilt * xd%PtfmRefXY(:,1) + (1.0-p%CXYFilt) * u%PRPMesh%TranslationDisp(1:2,1) + CASE DEFAULT ! 0: no drift + PtfmRefXY = 0.0_ReKi + END SELECT + !------------------------------------------------------------------- ! Additional stiffness, damping forces. These need to be placed on a point mesh which is located at the WAMIT reference point (WRP). ! This mesh will need to get mapped by the glue code for use by either ElastoDyn or SubDyn. @@ -1835,6 +1872,7 @@ SUBROUTINE HydroDyn_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, IF ( u%Morison%Mesh%Committed ) THEN ! Make sure we are using Morison / there is a valid mesh u%Morison%PtfmRefY = PtfmRefY + u%Morison%PtfmRefXY = PtfmRefXY u%Morison%PRP = u%PRPMesh%Position(:,1)+u%PRPMesh%TranslationDisp(:,1) CALL Morison_CalcOutput( Time, u%Morison, p%Morison, x%Morison, xd%Morison, & z%Morison, OtherState%Morison, y%Morison, m%Morison, & diff --git a/modules/hydrodyn/src/HydroDyn.txt b/modules/hydrodyn/src/HydroDyn.txt index 7e64865aba..bd0b2328f4 100644 --- a/modules/hydrodyn/src/HydroDyn.txt +++ b/modules/hydrodyn/src/HydroDyn.txt @@ -116,6 +116,7 @@ typedef ^ ContinuousStateType Morison_Con typedef ^ DiscreteStateType WAMIT_DiscreteStateType WAMIT {:} - - "discrete states from the wamit module" - typedef ^ DiscreteStateType Morison_DiscreteStateType Morison - - - "discrete states from the Morison module" - typedef ^ DiscreteStateType ReKi PtfmRefY {:} - - "Reference yaw position of the PRP relative to the inertial frame - Current step and two previous steps" (radians) +typedef ^ DiscreteStateType ReKi PtfmRefXY {2}{3} - - "Low-pass filtered x,y drift of the PRP/HD origin - current and two previous steps, consistent with WAMIT ExctnDisp=2" (m) # # # Define constraint states here: @@ -168,6 +169,8 @@ typedef ^ ^ logical typedef ^ ^ SeaSt_WaveFieldType *WaveField - - - "Pointer to SeaState wave field" - typedef ^ ^ INTEGER PtfmYMod - - - "Large yaw model" - typedef ^ ^ ReKi CYawFilt - - - "Low-pass filter constant for reference platform yaw position PtfmRefY" - +typedef ^ ^ INTEGER PtfmXYMod - - - "X/Y drift mode for the PRP, from the user-specified ExctnDisp applied regardless of potential-flow bodies (0 none, 1 instantaneous, 2 filtered). Used to displace the Morison members when WaveDisp=0" - +typedef ^ ^ ReKi CXYFilt - - - "Low-pass filter constant for the PRP x,y drift position (from ExctnCutOff)" - # # # ..... Inputs .................................................................................................................... diff --git a/modules/hydrodyn/src/HydroDyn_Input.f90 b/modules/hydrodyn/src/HydroDyn_Input.f90 index b091062a9f..8de78b6334 100644 --- a/modules/hydrodyn/src/HydroDyn_Input.f90 +++ b/modules/hydrodyn/src/HydroDyn_Input.f90 @@ -364,10 +364,6 @@ SUBROUTINE HydroDyn_ParseInput( InputFileName, OutRootName, FileInfo_In, InputFi call ParseVar( FileInfo_In, CurLine, 'WaveDisp', InputFileData%Morison%WaveDisp, ErrStat2, ErrMsg2, UnEc ) if (Failed()) return; - ! AMMod - Method of computing distributed added-mass force. {0: nodes below SWL when undisplaced. 1: Up to the free surface} (switch) - 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; @@ -1506,15 +1502,16 @@ SUBROUTINE HydroDynInput_ProcessInitData( InitInp, Interval, InputFileData, ErrS END IF + ! ExctnCutOff - validate before ExctnDisp may be forced to 0 below. ExctnDisp=2 also drives the Morison WaveDisp=0 drift + ! (PtfmRefXY), so a valid cutoff is required whenever ExctnDisp=2, even without potential-flow bodies. + if ( InputFileData%WAMIT%ExctnDisp == 2 .and. InputFileData%WAMIT%ExctnCutOff <= 0.0 ) then + CALL SetErrStat( ErrID_Fatal,'ExctnCutOff must be greater than zero.',ErrStat,ErrMsg,RoutineName) + end if + ! ExctnDisp - Method of computing Wave Excitation if ( InputFileData%PotMod /= 1 .or. InputFileData%WAMIT%ExctnMod == 0 .or. InitInp%WaveField%WaveMod == WaveMod_None) then InputFileData%WAMIT%ExctnDisp = 0 !Force ExctnDisp = 0, so that the Grid of Wave Excitation forces is not computed (saves time and memory) end if - - ! ExctnCutOff - if ( InputFileData%PotMod == 1 .and. InputFileData%WAMIT%ExctnMod > 0 .and. InputFileData%WAMIT%ExctnDisp == 2 .and. InputFileData%WAMIT%ExctnCutOff <= 0.0 ) then - CALL SetErrStat( ErrID_Fatal,'ExctnCutOff must be greater than zero.',ErrStat,ErrMsg,RoutineName) - end if ! PtfmVol0 - Displaced volume of water when the platform is in its undisplaced position @@ -1874,10 +1871,6 @@ SUBROUTINE HydroDynInput_ProcessInitData( InitInp, Interval, InputFileData, ErrS CALL SetErrStat( ErrID_Fatal,'WaveDisp must be 0 or 1',ErrStat,ErrMsg,RoutineName) RETURN END IF - IF ( InputFileData%Morison%AMMod /= 0 .AND. InputFileData%Morison%AMMod /= 1) THEN - 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 @@ -2745,10 +2738,6 @@ SUBROUTINE HydroDynInput_ProcessInitData( InitInp, Interval, InputFileData, ErrS if ( InputFileData%PtfmYCutOff <= 0.0_ReKi ) then CALL SetErrStat( ErrID_Fatal, 'PtfmYCutOff must be greater than 0 Hz.',ErrStat,ErrMsg,RoutineName) end if - if ( InputFileData%Morison%WaveDisp == 0 .AND. InputFileData%Morison%NMembers > 0 ) then - call SetErrStat( ErrID_Fatal,'Dynamic reference yaw offset (PtfmYMod=1) cannot be used with WaveDisp=0. Set WaveDisp=1.',ErrStat,ErrMsg,RoutineName) - return - end if if ( InputFileData%PotMod > 0 .AND. InputFileData%WAMIT%ExctnMod == 2 ) then call SetErrStat( ErrID_Fatal,'Dynamic reference yaw offset (PtfmYMod=1) cannot be used with state-space wave excitations. Set ExctnMod=0 or 1.', ErrStat, ErrMsg, RoutineName ) return diff --git a/modules/hydrodyn/src/HydroDyn_Types.f90 b/modules/hydrodyn/src/HydroDyn_Types.f90 index 7ceb823149..857d30e4b5 100644 --- a/modules/hydrodyn/src/HydroDyn_Types.f90 +++ b/modules/hydrodyn/src/HydroDyn_Types.f90 @@ -133,6 +133,7 @@ MODULE HydroDyn_Types TYPE(WAMIT_DiscreteStateType) , DIMENSION(:), ALLOCATABLE :: WAMIT !< discrete states from the wamit module [-] TYPE(Morison_DiscreteStateType) :: Morison !< discrete states from the Morison module [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: PtfmRefY !< Reference yaw position of the PRP relative to the inertial frame - Current step and two previous steps [(radians)] + REAL(ReKi) , DIMENSION(1:2,1:3) :: PtfmRefXY = 0.0_ReKi !< Low-pass filtered x,y drift of the PRP/HD origin - current and two previous steps, consistent with WAMIT ExctnDisp=2 [(m)] END TYPE HydroDyn_DiscreteStateType ! ======================= ! ========= HydroDyn_ConstraintStateType ======= @@ -184,6 +185,8 @@ MODULE HydroDyn_Types TYPE(SeaSt_WaveFieldType) , POINTER :: WaveField => NULL() !< Pointer to SeaState wave field [-] INTEGER(IntKi) :: PtfmYMod = 0_IntKi !< Large yaw model [-] REAL(ReKi) :: CYawFilt = 0.0_ReKi !< Low-pass filter constant for reference platform yaw position PtfmRefY [-] + INTEGER(IntKi) :: PtfmXYMod = 0_IntKi !< X/Y drift mode for the PRP, from the user-specified ExctnDisp applied regardless of potential-flow bodies (0 none, 1 instantaneous, 2 filtered). Used to displace the Morison members when WaveDisp=0 [-] + REAL(ReKi) :: CXYFilt = 0.0_ReKi !< Low-pass filter constant for the PRP x,y drift position (from ExctnCutOff) [-] END TYPE HydroDyn_ParameterType ! ======================= ! ========= HydroDyn_InputType ======= @@ -233,21 +236,22 @@ MODULE HydroDyn_Types integer(IntKi), public, parameter :: HydroDyn_x_Morison_DummyContState = 4 ! HydroDyn%Morison%DummyContState integer(IntKi), public, parameter :: HydroDyn_u_Morison_Mesh = 5 ! HydroDyn%Morison%Mesh integer(IntKi), public, parameter :: HydroDyn_u_Morison_PtfmRefY = 6 ! HydroDyn%Morison%PtfmRefY - integer(IntKi), public, parameter :: HydroDyn_u_Morison_PRP = 7 ! HydroDyn%Morison%PRP - integer(IntKi), public, parameter :: HydroDyn_u_WAMITMesh = 8 ! HydroDyn%WAMITMesh - integer(IntKi), public, parameter :: HydroDyn_u_PRPMesh = 9 ! HydroDyn%PRPMesh - integer(IntKi), public, parameter :: HydroDyn_u_qAddDOF = 10 ! HydroDyn%qAddDOF - integer(IntKi), public, parameter :: HydroDyn_u_qAddDOFDot = 11 ! HydroDyn%qAddDOFDot - integer(IntKi), public, parameter :: HydroDyn_u_qAddDOFDotDot = 12 ! HydroDyn%qAddDOFDotDot - integer(IntKi), public, parameter :: HydroDyn_y_WAMIT_Mesh = 13 ! HydroDyn%WAMIT(DL%i1)%Mesh - integer(IntKi), public, parameter :: HydroDyn_y_WAMIT_FAddDOF = 14 ! HydroDyn%WAMIT(DL%i1)%FAddDOF - integer(IntKi), public, parameter :: HydroDyn_y_WAMIT2_Mesh = 15 ! HydroDyn%WAMIT2(DL%i1)%Mesh - integer(IntKi), public, parameter :: HydroDyn_y_Morison_Mesh = 16 ! HydroDyn%Morison%Mesh - integer(IntKi), public, parameter :: HydroDyn_y_Morison_VisMesh = 17 ! HydroDyn%Morison%VisMesh - integer(IntKi), public, parameter :: HydroDyn_y_Morison_WriteOutput = 18 ! HydroDyn%Morison%WriteOutput - integer(IntKi), public, parameter :: HydroDyn_y_WAMITMesh = 19 ! HydroDyn%WAMITMesh - integer(IntKi), public, parameter :: HydroDyn_y_WriteOutput = 20 ! HydroDyn%WriteOutput - integer(IntKi), public, parameter :: HydroDyn_y_FAddDOF = 21 ! HydroDyn%FAddDOF + integer(IntKi), public, parameter :: HydroDyn_u_Morison_PtfmRefXY = 7 ! HydroDyn%Morison%PtfmRefXY + integer(IntKi), public, parameter :: HydroDyn_u_Morison_PRP = 8 ! HydroDyn%Morison%PRP + integer(IntKi), public, parameter :: HydroDyn_u_WAMITMesh = 9 ! HydroDyn%WAMITMesh + integer(IntKi), public, parameter :: HydroDyn_u_PRPMesh = 10 ! HydroDyn%PRPMesh + integer(IntKi), public, parameter :: HydroDyn_u_qAddDOF = 11 ! HydroDyn%qAddDOF + integer(IntKi), public, parameter :: HydroDyn_u_qAddDOFDot = 12 ! HydroDyn%qAddDOFDot + integer(IntKi), public, parameter :: HydroDyn_u_qAddDOFDotDot = 13 ! HydroDyn%qAddDOFDotDot + integer(IntKi), public, parameter :: HydroDyn_y_WAMIT_Mesh = 14 ! HydroDyn%WAMIT(DL%i1)%Mesh + integer(IntKi), public, parameter :: HydroDyn_y_WAMIT_FAddDOF = 15 ! HydroDyn%WAMIT(DL%i1)%FAddDOF + integer(IntKi), public, parameter :: HydroDyn_y_WAMIT2_Mesh = 16 ! HydroDyn%WAMIT2(DL%i1)%Mesh + integer(IntKi), public, parameter :: HydroDyn_y_Morison_Mesh = 17 ! HydroDyn%Morison%Mesh + integer(IntKi), public, parameter :: HydroDyn_y_Morison_VisMesh = 18 ! HydroDyn%Morison%VisMesh + integer(IntKi), public, parameter :: HydroDyn_y_Morison_WriteOutput = 19 ! HydroDyn%Morison%WriteOutput + integer(IntKi), public, parameter :: HydroDyn_y_WAMITMesh = 20 ! HydroDyn%WAMITMesh + integer(IntKi), public, parameter :: HydroDyn_y_WriteOutput = 21 ! HydroDyn%WriteOutput + integer(IntKi), public, parameter :: HydroDyn_y_FAddDOF = 22 ! HydroDyn%FAddDOF contains @@ -1066,8 +1070,8 @@ subroutine HydroDyn_CopyDiscState(SrcDiscStateData, DstDiscStateData, CtrlCode, integer(IntKi), intent(in ) :: CtrlCode integer(IntKi), intent( out) :: ErrStat character(*), intent( out) :: ErrMsg - integer(B4Ki) :: i1 - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: i1, i2 + integer(B4Ki) :: LB(2), UB(2) integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 character(*), parameter :: RoutineName = 'HydroDyn_CopyDiscState' @@ -1104,14 +1108,15 @@ subroutine HydroDyn_CopyDiscState(SrcDiscStateData, DstDiscStateData, CtrlCode, end if DstDiscStateData%PtfmRefY = SrcDiscStateData%PtfmRefY end if + DstDiscStateData%PtfmRefXY = SrcDiscStateData%PtfmRefXY end subroutine subroutine HydroDyn_DestroyDiscState(DiscStateData, ErrStat, ErrMsg) type(HydroDyn_DiscreteStateType), intent(inout) :: DiscStateData integer(IntKi), intent( out) :: ErrStat character(*), intent( out) :: ErrMsg - integer(B4Ki) :: i1 - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: i1, i2 + integer(B4Ki) :: LB(2), UB(2) integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 character(*), parameter :: RoutineName = 'HydroDyn_DestroyDiscState' @@ -1137,8 +1142,8 @@ subroutine HydroDyn_PackDiscState(RF, Indata) type(RegFile), intent(inout) :: RF type(HydroDyn_DiscreteStateType), intent(in) :: InData character(*), parameter :: RoutineName = 'HydroDyn_PackDiscState' - integer(B4Ki) :: i1 - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: i1, i2 + integer(B4Ki) :: LB(2), UB(2) if (RF%ErrStat >= AbortErrLev) return call RegPack(RF, allocated(InData%WAMIT)) if (allocated(InData%WAMIT)) then @@ -1151,6 +1156,7 @@ subroutine HydroDyn_PackDiscState(RF, Indata) end if call Morison_PackDiscState(RF, InData%Morison) call RegPackAlloc(RF, InData%PtfmRefY) + call RegPack(RF, InData%PtfmRefXY) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -1158,8 +1164,8 @@ subroutine HydroDyn_UnPackDiscState(RF, OutData) type(RegFile), intent(inout) :: RF type(HydroDyn_DiscreteStateType), intent(inout) :: OutData character(*), parameter :: RoutineName = 'HydroDyn_UnPackDiscState' - integer(B4Ki) :: i1 - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: i1, i2 + integer(B4Ki) :: LB(2), UB(2) integer(IntKi) :: stat logical :: IsAllocAssoc if (RF%ErrStat /= ErrID_None) return @@ -1178,6 +1184,7 @@ subroutine HydroDyn_UnPackDiscState(RF, OutData) end if call Morison_UnpackDiscState(RF, OutData%Morison) ! Morison call RegUnpackAlloc(RF, OutData%PtfmRefY); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%PtfmRefXY); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine HydroDyn_CopyConstrState(SrcConstrStateData, DstConstrStateData, CtrlCode, ErrStat, ErrMsg) @@ -1499,6 +1506,8 @@ subroutine HydroDyn_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, Err DstParamData%WaveField => SrcParamData%WaveField DstParamData%PtfmYMod = SrcParamData%PtfmYMod DstParamData%CYawFilt = SrcParamData%CYawFilt + DstParamData%PtfmXYMod = SrcParamData%PtfmXYMod + DstParamData%CXYFilt = SrcParamData%CXYFilt end subroutine subroutine HydroDyn_DestroyParam(ParamData, ErrStat, ErrMsg) @@ -1637,6 +1646,8 @@ subroutine HydroDyn_PackParam(RF, Indata) end if call RegPack(RF, InData%PtfmYMod) call RegPack(RF, InData%CYawFilt) + call RegPack(RF, InData%PtfmXYMod) + call RegPack(RF, InData%CXYFilt) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -1739,6 +1750,8 @@ subroutine HydroDyn_UnPackParam(RF, OutData) end if call RegUnpack(RF, OutData%PtfmYMod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%CYawFilt); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%PtfmXYMod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%CXYFilt); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine HydroDyn_CopyInput(SrcInputData, DstInputData, CtrlCode, ErrStat, ErrMsg) @@ -2888,6 +2901,8 @@ subroutine HydroDyn_VarPackInput(V, u, ValAry) call MV_PackMesh(V, u%Morison%Mesh, ValAry) ! Mesh case (HydroDyn_u_Morison_PtfmRefY) VarVals(1) = u%Morison%PtfmRefY ! Scalar + case (HydroDyn_u_Morison_PtfmRefXY) + VarVals = u%Morison%PtfmRefXY(V%iLB:V%iUB) ! Rank 1 Array case (HydroDyn_u_Morison_PRP) VarVals = u%Morison%PRP(V%iLB:V%iUB) ! Rank 1 Array case (HydroDyn_u_WAMITMesh) @@ -2926,6 +2941,8 @@ subroutine HydroDyn_VarUnpackInput(V, ValAry, u) call MV_UnpackMesh(V, ValAry, u%Morison%Mesh) ! Mesh case (HydroDyn_u_Morison_PtfmRefY) u%Morison%PtfmRefY = VarVals(1) ! Scalar + case (HydroDyn_u_Morison_PtfmRefXY) + u%Morison%PtfmRefXY(V%iLB:V%iUB) = VarVals ! Rank 1 Array case (HydroDyn_u_Morison_PRP) u%Morison%PRP(V%iLB:V%iUB) = VarVals ! Rank 1 Array case (HydroDyn_u_WAMITMesh) @@ -2950,6 +2967,8 @@ function HydroDyn_InputFieldName(DL) result(Name) Name = "u%Morison%Mesh" case (HydroDyn_u_Morison_PtfmRefY) Name = "u%Morison%PtfmRefY" + case (HydroDyn_u_Morison_PtfmRefXY) + Name = "u%Morison%PtfmRefXY" case (HydroDyn_u_Morison_PRP) Name = "u%Morison%PRP" case (HydroDyn_u_WAMITMesh) diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index 4dfe7e438e..b2fa48c6f9 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -24,10 +24,8 @@ MODULE Morison USE Morison_Types USE Morison_Output USE SeaSt_WaveField - ! USE HydroDyn_Output_Types USE NWTC_Library USE YawOffset - IMPLICIT NONE @@ -35,11 +33,6 @@ MODULE Morison TYPE(ProgDesc), PARAMETER :: Morison_ProgDesc = ProgDesc( 'Morison', '', '' ) - INTERFACE Morison_DirCosMtrx - MODULE PROCEDURE Morison_DirCosMtrx_Spin - MODULE PROCEDURE Morison_DirCosMtrx_noSpin - END INTERFACE - ! ..... Public Subroutines ................................................................................................... PUBLIC:: Morison_GenerateSimulationNodes @@ -49,7 +42,7 @@ MODULE Morison CONTAINS !---------------------------------------------------------------------------------------------------------------------------------- -SUBROUTINE Morison_DirCosMtrx_Spin( pos0, pos1, spin, DirCos ) +SUBROUTINE Morison_DirCosMtrx( pos0, pos1, spin, DirCos ) ! Compute the direction cosine matrix given two end points and a spin angle for rectangular members ! Left multiplying DirCos with a vector in element local sys returns vector in global sys ! Updated to match the convention in SubDyn for consistency @@ -105,59 +98,12 @@ SUBROUTINE Morison_DirCosMtrx_Spin( pos0, pos1, spin, DirCos ) DirCos = matmul(DirCos,Rspin) END IF -END SUBROUTINE Morison_DirCosMtrx_Spin - -SUBROUTINE Morison_DirCosMtrx_noSpin( pos0, pos1, DirCos ) -! Compute the direction cosine matrix given two end points without spin for cylindrical members -! Left multiplying DirCos with a vector in element local sys returns vector in global sys -! Updated to match the convention in SubDyn for consistency - - REAL(ReKi), INTENT( IN ) :: pos0(3), pos1(3) - Real(ReKi), INTENT( OUT ) :: DirCos(3,3) - Real(DbKi) :: Le, Lexy - Real(DbKi) :: dx, dy, dz - - DirCos = 0.0 - - dx = pos1(1) - pos0(1) - dy = pos1(2) - pos0(2) - dz = pos1(3) - pos0(3) - - Lexy = sqrt( dx*dx + dy*dy ) - - IF ( EqualRealNos(Lexy, 0.0_DbKi) ) THEN - IF (dz > 0) THEN - DirCos(1,1) = 1.0 - DirCos(2,2) = 1.0 - DirCos(3,3) = 1.0 - ELSE - DirCos(1,1) = 1.0 - DirCos(2,2) = -1.0 - DirCos(3,3) = -1.0 - END IF - ELSE - Le = sqrt( dx*dx + dy*dy + dz*dz ) - - DirCos(1, 1) = dy/Lexy - DirCos(1, 2) = dx*dz/(Lexy*Le) - DirCos(1, 3) = dx/Le - - DirCos(2, 1) = -dx/Lexy - DirCos(2, 2) = dy*dz/(Lexy*Le) - DirCos(2, 3) = dy/Le - - DirCos(3, 1) = 0.0 - DirCos(3, 2) = -Lexy/Le - DirCos(3, 3) = dz/Le - END IF - -END SUBROUTINE Morison_DirCosMtrx_noSpin - +END SUBROUTINE Morison_DirCosMtrx SUBROUTINE GetDisplacedNodePosition( u, p, forceDisplaced, pos ) TYPE(Morison_InputType), INTENT(IN ) :: u !< Inputs at Time TYPE(Morison_ParameterType), INTENT(IN ) :: p !< Parameters - LOGICAL, INTENT(IN ) :: forceDisplaced ! Set to true to return the exact displaced position no matter WaveDisp or WaveStMod + LOGICAL, INTENT(IN ) :: forceDisplaced ! Set to true to return the exact displaced position regardless of WaveDisp REAL(ReKi), INTENT( OUT) :: pos(:,:) ! Displaced node positions REAL(ReKi) :: Orient(3,3) @@ -168,64 +114,57 @@ SUBROUTINE GetDisplacedNodePosition( u, p, forceDisplaced, pos ) pos = u%Mesh%Position pos(3,:) = pos(3,:) - p%WaveField%MSL2SWL ! Z position measured from the SWL IF ( (p%WaveDisp /= 0) .OR. forceDisplaced ) THEN - ! Use displaced X and Y position + ! Use the full displaced position (X, Y, and Z) pos(1,:) = pos(1,:) + u%Mesh%TranslationDisp(1,:) pos(2,:) = pos(2,:) + u%Mesh%TranslationDisp(2,:) - IF ( (p%WaveField%WaveStMod > 0) .OR. forceDisplaced ) THEN - ! Use displaced Z position only when wave stretching is enabled - pos(3,:) = pos(3,:) + u%Mesh%TranslationDisp(3,:) - END IF - ELSE ! p%WaveDisp=0 implies PtfmYMod=0 - ! Rotate the structure based on PtfmRefY (constant) + pos(3,:) = pos(3,:) + u%Mesh%TranslationDisp(3,:) + ELSE ! WaveDisp=0: position the structure using the reference yaw only + ! Rotate the structure based on PtfmRefY (reference yaw, may be static or dynamic) call GetPtfmRefYOrient(u%PtfmRefY, Orient, ErrStat2, ErrMsg2) pos = matmul(transpose(Orient),pos) + ! Add the x,y drift of the HD origin so Morison tracks the potential-flow bodies + pos(1,:) = pos(1,:) + u%PtfmRefXY(1) + pos(2,:) = pos(2,:) + u%PtfmRefXY(2) END IF END SUBROUTINE GetDisplacedNodePosition - -SUBROUTINE YawMember(member, PtfmRefY, ErrStat, ErrMsg) +SUBROUTINE RotateMemberNode(p, im, member, PtfmRefY, NOrientation, ErrStat, ErrMsg) + TYPE(Morison_ParameterType), INTENT(IN ) :: p + Integer(IntKi), intent(in ) :: im ! Member index into p%Members (pristine reference) Type(Morison_MemberType), intent(inout) :: member Real(ReKi), intent(in ) :: PtfmRefY + Real(R8Ki), intent(in ) :: NOrientation(3,3) Integer(IntKi), intent( out) :: ErrStat Character(*), intent( out) :: ErrMsg - Real(ReKi) :: k(3), x_hat(3), y_hat(3) - Real(ReKi) :: kkt(3,3) - Real(ReKi) :: Ak(3,3) - Integer(IntKi) :: ErrStat2 - Character(ErrMsgLen) :: ErrMsg2 + Real(ReKi) :: Rg2b(3,3), Rb2g(3,3) - Character(*), parameter :: RoutineName = 'YawMember' + Character(*), parameter :: RoutineName = 'RotateMemberNode' ErrStat = ErrID_None ErrMsg = '' - call hiFrameTransform(h2i,PtfmRefY,member%k,k,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - member%k = k - - call hiFrameTransform(h2i,PtfmRefY,member%kkt,kkt,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - member%kkt = kkt - - call hiFrameTransform(h2i,PtfmRefY,member%Ak,Ak,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - member%Ak = Ak - - IF (member%MSecGeom == MSecGeom_Rec) THEN - - call hiFrameTransform(h2i,PtfmRefY,member%x_hat,x_hat,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - member%x_hat = x_hat + select case (p%WaveDisp) + case (0) + call GetPtfmRefYOrient(PtfmRefY, Rg2b, ErrStat, ErrMsg) + if (ErrStat /= ErrID_None) return + case (1) + Rg2b = NOrientation + end select + Rb2g = transpose(Rg2b) - call hiFrameTransform(h2i,PtfmRefY,member%y_hat,y_hat,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - member%y_hat = y_hat + ! Rotate the reference-configuration vectors to the current global orientation at this node. + ! Source is p%Members(im) so repeated per-node calls do not compound on the overwritten local copy. + member%k = matmul(Rb2g, p%Members(im)%k) + member%kkt = matmul(Rb2g, matmul(p%Members(im)%kkt, Rg2b)) + member%Ak = matmul(Rb2g, matmul(p%Members(im)%Ak, Rg2b)) + member%x_hat = matmul(Rb2g, p%Members(im)%x_hat) + member%y_hat = matmul(Rb2g, p%Members(im)%y_hat) - END IF + ! Note: Deliberately left member%CMatrix unchanged because we need the undisplaced version outside when constructing the full nodal orientation. -END SUBROUTINE YawMember +END SUBROUTINE RotateMemberNode !==================================================================================================== SUBROUTINE GetDistance ( a, b, l ) @@ -240,261 +179,30 @@ SUBROUTINE GetDistance ( a, b, l ) END SUBROUTINE GetDistance -!==================================================================================================== -SUBROUTINE ElementCentroid ( Rs, Re, p1, h, DCM, centroid ) -! This private subroutine computes the centroid of a tapered right cylinder element. -!---------------------------------------------------------------------------------------------------- - - REAL(ReKi), INTENT ( IN ) :: Rs ! starting radius - REAL(ReKi), INTENT ( IN ) :: Re ! ending radius - REAL(ReKi), INTENT ( IN ) :: p1(3) ! starting point of the element in global coordinates - REAL(ReKi), INTENT ( IN ) :: h ! height of the element - REAL(ReKi), INTENT ( IN ) :: DCM(3,3) ! direction cosine matrix to transform local element coordinates to global coordinates - REAL(ReKi), INTENT ( OUT ) :: centroid(3) ! centroid of the element in local coordinates - - centroid(1) = 0.0 - centroid(2) = 0.0 - centroid(3) = h * (Rs*Rs + 2.0*Rs*Re + 3.0*Re*Re) / (4.0*( Rs*Rs + Rs*Re + Re*Re ) ) !( 2.0*Re + Rs ) / ( 3.0 * ( Rs + Re ) ) - centroid = matmul( DCM, centroid ) + p1 - -END SUBROUTINE ElementCentroid - -!==================================================================================================== -REAL(ReKi) FUNCTION ElementVolume ( Rs, Re, h ) -! This private function computes the volume of a tapered right cylinder element. -!---------------------------------------------------------------------------------------------------- - - REAL(ReKi), INTENT ( IN ) :: Rs ! starting radius - REAL(ReKi), INTENT ( IN ) :: Re ! ending radius - REAL(ReKi), INTENT ( IN ) :: h ! height of the element - - ElementVolume = Pi*h*( Rs*Rs + Re*Re + Rs*Re ) / 3.0 - -END FUNCTION ElementVolume - -!==================================================================================================== -SUBROUTINE FindInterpFactor( p, p1, p2, s ) - - REAL(ReKi), INTENT ( IN ) :: p, p1, p2 - REAL(ReKi), INTENT ( OUT ) :: s - - REAL(ReKi) :: dp -! find normalized interpolation factor, s, such: -! p = p1*(1-s) + p2*s -! *--------------*--------------------------------* -! p1 p p2 -! -! 0-----------------------------------------------1 -! <------- s ----> - - dp = p2 - p1 - IF ( EqualRealNos(dp, 0.0_ReKi) ) THEN - s = 0 - ELSE - s = ( p - p1 ) / dp - END IF - -END SUBROUTINE FindInterpFactor -!======================================================================= -FUNCTION InterpWrappedStpInt( XValIn, XAry, YAry, Ind, AryLen ) - - - ! This function returns a y-value that corresponds to an input x-value which is wrapped back - ! into the range [1-XAry(AryLen). It finds a x-value which corresponds to a value in the XAry where XAry(Ind-1) < MOD(XValIn, XAry(AryLen)) <= XAry(Ind) - ! It is assumed that XAry is sorted in ascending order. - ! It uses the passed index as the starting point and does a stepwise interpolation from there. This is - ! especially useful when the calling routines save the value from the last time this routine was called - ! for a given case where XVal does not change much from call to call. . - ! - ! This routine assumes YAry is INTEGER. - - - ! Function declaration. - - INTEGER :: InterpWrappedStpInt ! This function. - - - ! Argument declarations. - - INTEGER, INTENT(IN) :: AryLen ! Length of the arrays. - INTEGER, INTENT(INOUT) :: Ind ! Initial and final index into the arrays. - - REAL(SiKi), INTENT(IN) :: XAry (AryLen) ! Array of X values to be interpolated. - REAL(SiKi), INTENT(IN) :: XValIn ! X value to be interpolated. - INTEGER, INTENT(IN) :: YAry (AryLen) ! Array of Y values to be interpolated. - - REAL(SiKi) :: XVal ! X value to be interpolated. - - - - ! Wrap XValIn into the range XAry(1) to XAry(AryLen) - XVal = MOD(XValIn, XAry(AryLen)) - - ! Set the Ind to the first index if we are at the beginning of XAry - IF ( XVal <= XAry(2) ) THEN - Ind = 1 - END IF - - - ! Let's check the limits first. - - IF ( XVal <= XAry(1) ) THEN - InterpWrappedStpInt = YAry(1) - Ind = 1 - RETURN - ELSE IF ( XVal >= XAry(AryLen) ) THEN - InterpWrappedStpInt = YAry(AryLen) - Ind = MAX(AryLen - 1, 1) - RETURN - END IF - - - ! Let's interpolate! - - Ind = MAX( MIN( Ind, AryLen-1 ), 1 ) - - DO - - IF ( XVal < XAry(Ind) ) THEN - - Ind = Ind - 1 - - ELSE IF ( XVal >= XAry(Ind+1) ) THEN - - Ind = Ind + 1 - - ELSE - - InterpWrappedStpInt = YAry(Ind) - RETURN - - END IF - - END DO - - RETURN -END FUNCTION InterpWrappedStpInt ! ( XVal, XAry, YAry, Ind, AryLen ) - - -!======================================================================= -FUNCTION InterpWrappedStpLogical( XValIn, XAry, YAry, Ind, AryLen ) - - - ! This function returns a y-value that corresponds to an input x-value which is wrapped back - ! into the range [0-XAry(AryLen) by interpolating into the arrays. - ! It is assumed that XAry is sorted in ascending order. - ! It uses the passed index as the starting point and does a stepwise interpolation from there. This is - ! especially useful when the calling routines save the value from the last time this routine was called - ! for a given case where XVal does not change much from call to call. When there is no correlation - ! from one interpolation to another, InterpBin() may be a better choice. - ! It returns the first or last YAry() value if XVal is outside the limits of XAry(). - ! This routine assumes YAry is REAL. - - - ! Function declaration. - - LOGICAL :: InterpWrappedStpLogical ! This function. - - - ! Argument declarations. - - INTEGER, INTENT(IN) :: AryLen ! Length of the arrays. - INTEGER, INTENT(INOUT) :: Ind ! Initial and final index into the arrays. - - REAL(SiKi), INTENT(IN) :: XAry (AryLen) ! Array of X values to be interpolated. - REAL(SiKi), INTENT(IN) :: XValIn ! X value to be interpolated. - LOGICAL, INTENT(IN) :: YAry (AryLen) ! Array of Y values to be interpolated. - - REAL(SiKi) :: XVal ! X value to be interpolated. - - - - ! Wrap XValIn into the range XAry(1) to XAry(AryLen) - XVal = MOD(XValIn, XAry(AryLen)) - - ! Set the Ind to the first index if we are at the beginning of XAry - IF ( XVal <= XAry(2) ) THEN - Ind = 1 - END IF - - - ! Let's check the limits first. - - IF ( XVal <= XAry(1) ) THEN - InterpWrappedStpLogical = YAry(1) - Ind = 1 - RETURN - ELSE IF ( XVal >= XAry(AryLen) ) THEN - InterpWrappedStpLogical = YAry(AryLen) - Ind = MAX(AryLen - 1, 1) - RETURN - END IF - - - ! Let's interpolate! - - Ind = MAX( MIN( Ind, AryLen-1 ), 1 ) - - DO - - IF ( XVal < XAry(Ind) ) THEN - - Ind = Ind - 1 - - ELSE IF ( XVal >= XAry(Ind+1) ) THEN - - Ind = Ind + 1 - - ELSE +!---------------------------------------------------------------------------------------------------------------------------------- +subroutine GetElementAxialVec(p1, p2, k_hat, errStat, errMsg) + ! Instantaneous elemental axial unit vector from node p1 to node p2. + ! Guards against zero-length elements (do not assume rigid structure). + real(ReKi), intent(in ) :: p1(3),p2(3) + real(ReKi), intent( out) :: k_hat(3) + integer, intent( out) :: errStat ! returns a non-zero value when an error occurs + character(*), intent( out) :: errMsg ! Error message if errStat /= ErrID_None + character(*), parameter :: RoutineName = 'GetElementAxialVec' - InterpWrappedStpLogical = YAry(Ind) - RETURN + real(ReKi) :: vec(3), vecLen - END IF + errStat = ErrID_None + errMsg = "" - END DO + vec = p2 - p1 + vecLen = SQRT(Dot_Product(vec,vec)) + if ( vecLen < 0.000001 ) then + call SeterrStat(ErrID_Fatal, 'An element of the Morison structure has co-located endpoints! This should never occur. Please review your model.', errStat, errMsg, RoutineName ) + return + end if + k_hat = vec / vecLen - RETURN -END FUNCTION InterpWrappedStpLogical ! ( XVal, XAry, YAry, Ind, AryLen ) -!---------------------------------------------------------------------------------------------------------------------------------- -subroutine GetOrientationAngles(p1, p2, phi, sinPhi, cosPhi, tanPhi, sinBeta, cosBeta, k_hat, errStat, errMsg) - real(ReKi), intent(in ) :: p1(3),p2(3) - real(ReKi), intent( out) :: phi, sinPhi, cosPhi, tanPhi, sinBeta, cosBeta, k_hat(3) - integer, intent( out) :: errStat ! returns a non-zero value when an error occurs - character(*), intent( out) :: errMsg ! Error message if errStat /= ErrID_None - character(*), parameter :: RoutineName = 'GetOrientationAngles' - - real(ReKi) :: vec(3), vecLen, vecLen2D, beta - - ! Initialize errStat - - errStat = ErrID_None - errMsg = "" - - ! calculate isntantaneous incline angle and heading, and related trig values - ! the first and last NodeIndx values point to the corresponding Joint nodes idices which are at the start of the Mesh - vec = p2 - p1 - vecLen = SQRT(Dot_Product(vec,vec)) - vecLen2D = SQRT(vec(1)**2+vec(2)**2) - if ( vecLen < 0.000001 ) then - call SeterrStat(ErrID_Fatal, 'An element of the Morison structure has co-located endpoints! This should never occur. Please review your model.', errStat, errMsg, RoutineName ) - return - else - k_hat = vec / vecLen - phi = atan2(vecLen2D, vec(3)) ! incline angle - end if - if ( EqualRealNos(phi, 0.0_ReKi) ) then - beta = 0.0_ReKi - else - beta = atan2(vec(2), vec(1)) ! heading of incline - endif - sinPhi = sin(phi) - cosPhi = cos(phi) - tanPhi = tan(phi) - sinBeta = sin(beta) - cosBeta = cos(beta) - -end subroutine GetOrientationAngles +end subroutine GetElementAxialVec !---------------------------------------------------------------------------------------------------------------------------------- !function to return conical taper geometry calculations (volume and center of volume) SUBROUTINE CylTaperCalc(R1, R2, H, taperV, h_c) @@ -1391,8 +1099,6 @@ SUBROUTINE SetDepthBasedCoefs_Cyl( z, tMG, NCoefDpth, CoefDpths, Cd, Ca, Cp, AxC END DO ! Linearly interpolate the coef values based on depth - !CALL FindInterpFactor( z, CoefDpths(indx1)%Dpth, CoefDpths(indx2)%Dpth, s ) - dd = CoefDpths(indx1)%Dpth - CoefDpths(indx2)%Dpth IF ( EqualRealNos(dd, 0.0_ReKi) ) THEN s = 0 @@ -1465,8 +1171,6 @@ SUBROUTINE SetDepthBasedCoefs_Rec( z, tMG, NCoefDpth, CoefDpths, CdA, CdB, CaA, END DO ! Linearly interpolate the coef values based on depth - !CALL FindInterpFactor( z, CoefDpths(indx1)%Dpth, CoefDpths(indx2)%Dpth, s ) - dd = CoefDpths(indx1)%Dpth - CoefDpths(indx2)%Dpth IF ( EqualRealNos(dd, 0.0_ReKi) ) THEN s = 0 @@ -1989,7 +1693,6 @@ subroutine SetMemberProperties_Cyl( gravity, member, MCoefMod, MmbrCoefIDIndx, M real(ReKi) :: memLength real(ReKi) :: Za real(ReKi) :: Zb - real(ReKi) :: phi real(ReKi) :: sinPhi real(ReKi) :: cosPhi real(ReKi) :: Rmid @@ -2018,11 +1721,12 @@ subroutine SetMemberProperties_Cyl( gravity, member, MCoefMod, MmbrCoefIDIndx, M member%kkt = matmul(transpose(tk),tk) call Eye(Imat,errStat,errMsg) member%Ak = Imat - member%kkt - phi = acos( max(-1.0_ReKi, min(1.0_ReKi, vec(3)/memLength) ) ) ! incline angle - sinPhi = sin(phi) - cosPhi = cos(phi) - member%cosPhi_ref = cosPhi - + cosPhi = member%k(3) ! cosine of inclination from vertical (z-component of axial unit vector) + sinPhi = sqrt( member%k(1)**2 + member%k(2)**2 ) ! sine of inclination from vertical + CALL Morison_DirCosMtrx( InitInp%Nodes(member%NodeIndx(1))%Position, InitInp%Nodes(member%NodeIndx(N+1))%Position, 0.0_ReKi, member%CMatrix ) + member%x_hat = member%CMatrix(1:3,1) + member%y_hat = member%CMatrix(1:3,2) + ! These are all per node and not done here, yet do i = 1, member%NElements+1 @@ -2064,7 +1768,7 @@ subroutine SetMemberProperties_Cyl( gravity, member, MCoefMod, MmbrCoefIDIndx, M RETURN END IF ! Check inclination - If ( ABS(phi) .GE. 0.174533 ) THEN ! If inclination from vertical is greater than 10 deg + If ( ABS(cosPhi) <= 0.9848077530_ReKi ) THEN ! If inclination from vertical is greater than 10 deg CALL SetErrStat(ErrID_Fatal, 'MacCamy-Fuchs members must be within 10 degrees from vertical. This is not true for Member ID '//trim(num2lstr(member%MemberID)), errStat, errMsg, RoutineName ) RETURN END IF @@ -2127,7 +1831,7 @@ subroutine SetMemberProperties_Cyl( gravity, member, MCoefMod, MmbrCoefIDIndx, M call SetErrStat(ErrID_Fatal, 'The lower end-plate of a member must not cross the water plane. This is not true for Member ID '//trim(num2lstr(member%MemberID)), errStat, errMsg, RoutineName ) end if end if - if ( ( Za < -InitInp%WaveField%EffWtrDpth .and. Zb >= -InitInp%WaveField%EffWtrDpth ) .and. ( phi > 10.0*d2r .or. abs((member%RMG(N+1) - member%RMG(1))/member%RefLength)>0.1 ) ) then + if ( ( Za < -InitInp%WaveField%EffWtrDpth .and. Zb >= -InitInp%WaveField%EffWtrDpth ) .and. ( ABS(cosPhi) < 0.9848077530_ReKi .or. abs((member%RMG(N+1) - member%RMG(1))/member%RefLength)>0.1 ) ) then call SetErrStat(ErrID_Fatal, 'A member which crosses the seabed must not be inclined more than 10 degrees from vertical or have a taper larger than 0.1. This is not true for Member ID '//trim(num2lstr(member%MemberID)), errStat, errMsg, RoutineName ) end if @@ -2143,7 +1847,7 @@ subroutine SetMemberProperties_Cyl( gravity, member, MCoefMod, MmbrCoefIDIndx, M Za = InitInp%Nodes(member%NodeIndx(i))%Position(3) if (Za > -InitInp%WaveField%EffWtrDpth) then ! find the lowest node above the seabed - if (cosPhi < 0.173648178 ) then ! phi > 80 degrees and member is seabed crossing + if (ABS(cosPhi) < 0.173648178 ) then ! phi > 80 degrees and member is seabed crossing call SetErrStat(ErrID_Fatal, 'A seabed crossing member must have an inclination angle of <= 80 degrees from vertical. This is not true for Member ID '//trim(num2lstr(member%MemberID)), errStat, errMsg, RoutineName ) end if @@ -2341,8 +2045,6 @@ subroutine SetMemberProperties_Rec( gravity, member, MCoefMod, MmbrCoefIDIndx, M real(ReKi) :: memLength real(ReKi) :: Za real(ReKi) :: Zb - real(ReKi) :: phi - real(ReKi) :: sinPhi real(ReKi) :: cosPhi real(ReKi) :: SaMid, SbMid real(ReKi) :: SaMidMG, SbMidMG @@ -2350,7 +2052,7 @@ subroutine SetMemberProperties_Rec( gravity, member, MCoefMod, MmbrCoefIDIndx, M real(ReKi) :: Lmid real(ReKi) :: li real(ReKi) :: Vinner_l, Vinner_u, Vouter_l, Vouter_u, Vballast_l, Vballast_u - real(ReKi) :: tk(1,3), Imat(3,3), CMatrix(3,3) + real(ReKi) :: tk(1,3), Imat(3,3) REAL(ReKi) :: h_c ! center of mass offset from first node errStat = ErrID_None @@ -2370,15 +2072,10 @@ subroutine SetMemberProperties_Rec( gravity, member, MCoefMod, MmbrCoefIDIndx, M member%kkt = matmul(transpose(tk),tk) call Eye(Imat,errStat,errMsg) member%Ak = Imat - member%kkt - IF (member%MSecGeom == MSecGeom_Rec) THEN - CALL Morison_DirCosMtrx( InitInp%Nodes(member%NodeIndx(1))%Position, InitInp%Nodes(member%NodeIndx(N+1))%Position, member%MSpinOrient, CMatrix ) - member%x_hat = CMatrix(1:3,1) - member%y_hat = CMatrix(1:3,2) - END IF - phi = acos( max(-1.0_ReKi, min(1.0_ReKi, vec(3)/memLength) ) ) ! incline angle - sinPhi = sin(phi) - cosPhi = cos(phi) - member%cosPhi_ref = cosPhi + CALL Morison_DirCosMtrx( InitInp%Nodes(member%NodeIndx(1))%Position, InitInp%Nodes(member%NodeIndx(N+1))%Position, member%MSpinOrient, member%CMatrix ) + member%x_hat = member%CMatrix(1:3,1) + member%y_hat = member%CMatrix(1:3,2) + cosPhi = member%k(3) ! cosine of inclination from vertical (z-component of axial unit vector) ! These are all per node and not done here, yet @@ -2465,16 +2162,8 @@ subroutine SetMemberProperties_Rec( gravity, member, MCoefMod, MmbrCoefIDIndx, M ! Check the member does not exhibit any of the following conditions if (.not. member%PropPot) then - ! MHstLMod=1 is not allowed for rectangular members at the moment. Skip the following check. - ! if (member%MHstLMod == 1) then - ! if ( abs(Zb) < abs(member%Rmg(N+1)*sinPhi) ) then - ! call SetErrStat(ErrID_Fatal, 'The upper end-plate of a member must not cross the water plane. This is not true for Member ID '//trim(num2lstr(member%MemberID)), errStat, errMsg, RoutineName ) - ! end if - ! if ( abs(Za) < abs(member%Rmg(1)*sinPhi) ) then - ! call SetErrStat(ErrID_Fatal, 'The lower end-plate of a member must not cross the water plane. This is not true for Member ID '//trim(num2lstr(member%MemberID)), errStat, errMsg, RoutineName ) - ! end if - ! end if - if ( ( Za < -InitInp%WaveField%EffWtrDpth .and. Zb >= -InitInp%WaveField%EffWtrDpth ) .and. ( phi > 10.0*d2r .or. abs((member%SaMG(N+1) - member%SaMG(1))/member%RefLength)>0.1 .or. abs((member%SbMG(N+1) - member%SbMG(1))/member%RefLength)>0.1 ) ) then + ! MHstLMod=1 is not allowed for rectangular members. Skip the check for partially wetted endplates. + if ( ( Za < -InitInp%WaveField%EffWtrDpth .and. Zb >= -InitInp%WaveField%EffWtrDpth ) .and. ( ABS(cosPhi) < 0.9848077530_ReKi .or. abs((member%SaMG(N+1) - member%SaMG(1))/member%RefLength)>0.1 .or. abs((member%SbMG(N+1) - member%SbMG(1))/member%RefLength)>0.1 ) ) then call SetErrStat(ErrID_Fatal, 'A member which crosses the seabed must not be inclined more than 10 degrees from vertical or have a taper larger than 0.1. This is not true for Member ID '//trim(num2lstr(member%MemberID)), errStat, errMsg, RoutineName ) end if end if @@ -2488,7 +2177,7 @@ subroutine SetMemberProperties_Rec( gravity, member, MCoefMod, MmbrCoefIDIndx, M Za = InitInp%Nodes(member%NodeIndx(i))%Position(3) if (Za > -InitInp%WaveField%EffWtrDpth) then ! find the lowest node above the seabed - if (cosPhi < 0.173648178 ) then ! phi > 80 degrees and member is seabed crossing + if (ABS(cosPhi) < 0.173648178 ) then ! phi > 80 degrees and member is seabed crossing call SetErrStat(ErrID_Fatal, 'A seabed crossing member must have an inclination angle of <= 80 degrees from vertical. This is not true for Member ID '//trim(num2lstr(member%MemberID)), errStat, errMsg, RoutineName ) end if @@ -2809,7 +2498,6 @@ 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%OutAll = InitInp%OutAll 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 @@ -2818,12 +2506,6 @@ SUBROUTINE Morison_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, In ! Pointer to SeaState WaveField p%WaveField => InitInp%WaveField - ! Only compute added-mass force up to the free surface if wave stretching is enabled - IF ( p%WaveField%WaveStMod .EQ. 0_IntKi ) THEN - ! Setting AMMod to zero just in case. Probably redundant. - 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 @@ -3496,12 +3178,6 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, INTEGER :: N ! Number of elements within a given member REAL(ReKi) :: dl ! Element length within a given member, m REAL(ReKi) :: vec(3) ! Vector pointing from a member's 1st node to its last node - REAL(ReKi) :: phi, phi1, phi2 ! member tilt angle - REAL(ReKi) :: cosPhi, cosPhi1, cosPhi2 - REAL(ReKi) :: sinPhi, sinPhi1, sinPhi2 - REAL(ReKi) :: tanPhi - 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) :: Sa1, Sa2, Sa1b, Sa2b, SaMidb, Sa1In, Sa2In, SaMidIn @@ -3519,13 +3195,12 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, REAL(ReKi) :: a_s1(3) REAL(ReKi) :: alpha_s1(3) REAL(ReKi) :: omega_s1(3) - REAL(ReKi) :: a_s2(3) + REAL(ReKi) :: a_s2(3) REAL(ReKi) :: alpha_s2(3) REAL(ReKi) :: omega_s2(3) REAL(ReKi) :: pos1(3), pos2(3) REAL(ReKi) :: Imat(3,3) REAL(ReKi) :: iArm(3), iTerm(3), h_c, dRdl_p, dRdl_pp, dSadl_p, dSadl_pp, dSbdl_p, dSbdl_pp, f_hydro(3), Am(3,3), lstar, deltal, deltalLeft, deltalRight - REAL(ReKi) :: h, h_c_AM, deltal_AM REAL(ReKi) :: F_WMG(6), F_IMG(6), F_If(6), F_B0(6), F_B1(6), F_B2(6), F_B_End(6) REAL(ReKi) :: AM_End(3,3), An_End(3), DP_Const_End(3), I_MG_End(3,3) @@ -3578,13 +3253,13 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, !=============================================================================================== ! Get displaced positions of the hydrodynamic nodes - CALL GetDisplacedNodePosition( u, p, .FALSE., m%DispNodePosHdn ) ! For hydrodynamic loads; depends on WaveDisp and WaveStMod + CALL GetDisplacedNodePosition( u, p, .FALSE., m%DispNodePosHdn ) ! For hydrodynamic loads; depends on WaveDisp CALL GetDisplacedNodePosition( u, p, .TRUE. , m%DispNodePosHst ) ! For hydrostatic loads; always use actual displaced position !=============================================================================================== ! Calculate the fluid kinematics at all mesh nodes and store for use in the equations below CALL WaveField_GetWaveKin( p%WaveField, m%WaveField_m, Time, m%DispNodePosHdn, .FALSE., .TRUE., m%nodeInWater, m%WaveElev1, m%WaveElev2, m%WaveElev, m%FDynP, m%FV, m%FA, m%FAMCF, ErrStat2, ErrMsg2 ) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + if (Failed()) return ! Compute fluid velocity relative to the structure DO j = 1, p%NNodes m%vrel(:,j) = ( m%FV(:,j) - u%Mesh%TranslationVel(:,j) ) * m%nodeInWater(j) @@ -3628,8 +3303,6 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, DO im = 1, p%NMembers mem = p%Members(im) N = mem%NElements - call YawMember(mem, u%PtfmRefY, ErrStat2, ErrMsg2) - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) !zero member loads m%memberLoads(im)%F_B = 0.0_ReKi @@ -3641,32 +3314,22 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, m%memberLoads(im)%F_IMG = 0.0_ReKi m%memberLoads(im)%F_If = 0.0_ReKi - ! Determine member submergence status - IF ( p%WaveField%WaveStMod .EQ. 0_IntKi ) THEN ! No wave stretching - Only need to check the two ends - IF ( m%nodeInWater(mem%NodeIndx(1)) .NE. m%nodeInWater(mem%NodeIndx(N+1)) ) THEN - MemSubStat = 1_IntKi ! Member centerline crosses the SWL once - ELSE IF ( m%nodeInWater(mem%NodeIndx(1)) .EQ. 0_IntKi ) THEN + ! Determine member submergence status - check every node because members can deform + NumFSX = 0_IntKi ! Number of free-surface crossing + DO i = 1, N ! loop through member elements + IF ( m%nodeInWater(mem%NodeIndx(i)) .NE. m%nodeInWater(mem%NodeIndx(i+1)) ) THEN + NumFSX = NumFSX + 1 + END IF + END DO + IF (NumFSX .EQ. 1_IntKi) THEN + MemSubStat = 1_IntKi ! Member centerline crosses the free surface once + ELSE IF (NumFSX .GT. 1_IntKi) THEN + MemSubStat = 2_IntKi ! Member centerline crosses the free surface multiple time + ELSE ! Member centerline does not cross the free surface + IF ( m%nodeInWater(mem%NodeIndx(1)) .EQ. 0_IntKi ) THEN MemSubStat = 3_IntKi ! Member centerline completely above water ELSE - MemSubStat = 0_IntKi ! Member centerline fully submerged - END IF - ELSE IF ( p%WaveField%WaveStMod > 0_IntKi ) THEN ! Has wave stretching - Need to check every node - NumFSX = 0_IntKi ! Number of free-surface crossing - DO i = 1, N ! loop through member elements - IF ( m%nodeInWater(mem%NodeIndx(i)) .NE. m%nodeInWater(mem%NodeIndx(i+1)) ) THEN - NumFSX = NumFSX + 1 - END IF - END DO - IF (NumFSX .EQ. 1_IntKi) THEN - MemSubStat = 1_IntKi ! Member centerline crosses the free surface once - ELSE IF (NumFSX .GT. 1_IntKi) THEN - MemSubStat = 2_IntKi ! Member centerline crosses the free surface multiple time - ELSE ! Member centerline does not cross the free surface - IF ( m%nodeInWater(mem%NodeIndx(1)) .EQ. 0_IntKi ) THEN - MemSubStat = 3_IntKi ! Member centerline completely above water - ELSE - MemSubStat = 0_IntKi ! Member centerline completely submerged - END IF + MemSubStat = 0_IntKi ! Member centerline completely submerged END IF END IF @@ -3674,17 +3337,14 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, IF ( .NOT. mem%PropPot ) THEN ! Member is NOT modeled with Potential Flow Theory DO i = max(mem%i_floor,1), N ! loop through member elements that are not completely buried in the seabed - ! calculate instantaneous incline angle and heading, and related trig values + ! instantaneous elemental axial unit vector (1st -> 2nd node) ! the first and last NodeIndx values point to the corresponding Joint nodes indices which are at the start of the Mesh pos1 = m%DispNodePosHst(:, mem%NodeIndx(i )) pos2 = m%DispNodePosHst(:, mem%NodeIndx(i+1)) - call GetOrientationAngles( pos1, pos2, phi, sinPhi, cosPhi, tanPhi, sinBeta, cosBeta, k_hat, errStat2, errMsg2 ) - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - ! Compute element to global DirCos matrix for undisplaced structure first - call Morison_DirCosMtrx( u%Mesh%Position(:,mem%NodeIndx(i )), u%Mesh%Position(:,mem%NodeIndx(i+1)), mem%MSpinOrient, CMatrix ) - ! Prepend body motion - Assuming the rotation of the starting node is representative of the whole element - CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,mem%NodeIndx(i))),CMatrix) + call GetElementAxialVec( pos1, pos2, k_hat, errStat2, errMsg2 ); if (Failed()) return + ! Compute total element orientation matrix assuming the rotation of the starting node is representative of the whole element + CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,mem%NodeIndx(i))),mem%CMatrix) CTrans = transpose(CMatrix) ! Note: CMatrix is element local to global displaced. CTrans is the opposite. ! save some commonly used variables @@ -3695,14 +3355,15 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, a_s2 = u%Mesh%TranslationAcc(:, mem%NodeIndx(i+1)) alpha_s2 = u%Mesh%RotationAcc (:, mem%NodeIndx(i+1)) omega_s2 = u%Mesh%RotationVel (:, mem%NodeIndx(i+1)) - IF (mem%MSecGeom == MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) r1 = mem%RMG(i ) ! outer radius at element nodes including marine growth r2 = mem%RMG(i+1) r1b = mem%RMGB(i ) ! outer radius at element nodes including marine growth scaled by sqrt(Cb) r2b = mem%RMGB(i+1) dRdl_mg = mem%dRdl_mg(i) ! Taper of element including marine growth dRdl_mg_b = mem%dRdl_mg_b(i) ! Taper of element including marine growth with radius scaling by sqrt(Cb) - ELSE IF (mem%MSecGeom == MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) Sa1 = mem%SaMG(i ) ! outer side A at element nodes including marine growth Sa2 = mem%SaMG(i+1) Sb1 = mem%SbMG(i ) ! outer side B at element nodes including marine growth @@ -3715,7 +3376,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, dSadl_mg_b = mem%dSadl_mg_b(i) ! Taper of element side A including marine growth with radius scaling by sqrt(Cb) dSbdl_mg = mem%dSbdl_mg(i) ! Taper of element side B including marine growth dSbdl_mg_b = mem%dSbdl_mg_b(i) ! Taper of element side B including marine growth with radius scaling by sqrt(Cb) - END IF + END SELECT ! ------------------ marine growth: Sides: Section 4.1.2 -------------------- ! ----- marine growth weight @@ -3723,16 +3384,16 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! lower node F_WMG(3) = - mem%m_mg_l(i)*g ! weight force : Note: this is a constant - F_WMG(4) = - mem%m_mg_l(i)*g * mem%h_cmg_l(i)* sinPhi * sinBeta! weight force - F_WMG(5) = mem%m_mg_l(i)*g * mem%h_cmg_l(i)* sinPhi * cosBeta! weight force + F_WMG(4) = - mem%m_mg_l(i)*g * mem%h_cmg_l(i)* k_hat(2) ! weight force + F_WMG(5) = mem%m_mg_l(i)*g * mem%h_cmg_l(i)* k_hat(1) ! weight force m%memberLoads(im)%F_WMG(:,i) = m%memberLoads(im)%F_WMG(:,i) + F_WMG y%Mesh%Force (:,mem%NodeIndx(i)) = y%Mesh%Force (:,mem%NodeIndx(i)) + F_WMG(1:3) y%Mesh%Moment(:,mem%NodeIndx(i)) = y%Mesh%Moment(:,mem%NodeIndx(i)) + F_WMG(4:6) ! upper node F_WMG(3) = - mem%m_mg_u(i)*g ! weight force : Note: this is a constant - F_WMG(4) = - mem%m_mg_u(i)*g * mem%h_cmg_u(i)* sinPhi * sinBeta! weight force - F_WMG(5) = mem%m_mg_u(i)*g * mem%h_cmg_u(i)* sinPhi * cosBeta! weight force + F_WMG(4) = - mem%m_mg_u(i)*g * mem%h_cmg_u(i)* k_hat(2) ! weight force + F_WMG(5) = mem%m_mg_u(i)*g * mem%h_cmg_u(i)* k_hat(1) ! weight force m%memberLoads(im)%F_WMG(:,i+1) = m%memberLoads(im)%F_WMG(:,i+1) + F_WMG y%Mesh%Force (:,mem%NodeIndx(i+1)) = y%Mesh%Force (:,mem%NodeIndx(i+1)) + F_WMG(1:3) y%Mesh%Moment(:,mem%NodeIndx(i+1)) = y%Mesh%Moment(:,mem%NodeIndx(i+1)) + F_WMG(4:6) @@ -3740,13 +3401,14 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! ----- marine growth inertial load ! lower node Imat = 0.0_ReKi - IF (mem%MSecGeom == MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) Imat(1,1) = mem%I_rmg_l(i) Imat(2,2) = mem%I_rmg_l(i) - ELSE IF (mem%MSecGeom == MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) Imat(1,1) = mem%I_xmg_l(i) Imat(2,2) = mem%I_ymg_l(i) - END IF + END SELECT Imat(3,3) = mem%I_lmg_l(i) Imat = matmul(matmul(CMatrix, Imat), CTrans) iArm = mem%h_cmg_l(i) * k_hat @@ -3760,13 +3422,14 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! upper node Imat = 0.0_ReKi - IF (mem%MSecGeom == MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) Imat(1,1) = mem%I_rmg_u(i) Imat(2,2) = mem%I_rmg_u(i) - ELSE IF (mem%MSecGeom == MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) Imat(1,1) = mem%I_xmg_u(i) Imat(2,2) = mem%I_ymg_u(i) - END IF + END SELECT Imat(3,3) = mem%I_lmg_u(i) Imat = matmul(matmul(CMatrix, Imat), CTrans) iArm = mem%h_cmg_u(i) * k_hat @@ -3790,16 +3453,15 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, case (1) 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 ) + CALL GetTotalWaveElev(p, m, Time, pos1, Zeta1, ErrStat2, ErrMsg2 ); if (Failed()) return + CALL GetTotalWaveElev(p, m, Time, pos2, Zeta2, ErrStat2, ErrMsg2 ); if (Failed()) return ELSE ! Without wave stretching, compute buoyancy based on SWL Zeta1 = 0.0_ReKi Zeta2 = 0.0_ReKi END IF Is1stElement = ( i .EQ. 1) CALL getElementHstLds_Mod1(p, m, mem, Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1b, r2b, dl, mem%alpha(i), Is1stElement, F_B0, F_B1, F_B2, ErrStat2, ErrMsg2 ) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + if (Failed()) return ! Add nodal loads to mesh IF ( .NOT. Is1stElement ) THEN m%memberLoads(im)%F_B(:, i-1) = m%memberLoads(im)%F_B(:, i-1) + F_B0 @@ -3820,25 +3482,23 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, posMid = 0.5 * (pos1+pos2) 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 ) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL GetTotalWaveElev(p, m, Time, posMid, ZetaMid, ErrStat2, ErrMsg2 ); if (Failed()) return + CALL GetFreeSurfaceNormal( p, m, Time, posMid, n_hat, ErrStat2, ErrMsg2 ); if (Failed()) return FSPt = (/posMid(1),posMid(2),ZetaMid/) ! Reference point on the free surface ELSE FSPt = (/posMid(1),posMid(2),0.0_ReKi/) n_hat = (/0.0,0.0,1.0/) END IF - IF (mem%MSecGeom == MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) CALL GetSectionUnitVectors_Cyl( k_hat, y_hat, z_hat ) CALL getElementHstLds_Mod2_Cyl( p, pos1, pos2, FSPt, k_hat, y_hat, z_hat, n_hat, r1b, r2b, dl, F_B1, F_B2, ErrStat2, ErrMsg2) - ELSE IF (mem%MSecGeom == MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) CALL GetSectionUnitVectors_Rec( CMatrix, x_hat, y_hat ) CALL getElementHstLds_Mod2_Rec( p, pos1, pos2, FSPt, k_hat, x_hat, y_hat, n_hat, Sa1b, Sa2b, Sb1b, Sb2b, dl, F_B1, F_B2, ErrStat2, ErrMsg2) - END IF - - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + END SELECT + if (Failed()) return ! Add nodal loads to mesh m%memberLoads(im)%F_B(:,i ) = m%memberLoads(im)%F_B(:,i ) + F_B1 @@ -3858,16 +3518,14 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, zFillGroup = m%zFillGroup(mem%MmbrFilledIDIndx) DO i = max(mem%i_floor,1), N ! loop through member elements that are not completely buried in the seabed IF (mem%floodstatus(i)>0) THEN - ! calculate instantaneous incline angle and heading, and related trig values + ! instantaneous elemental axial unit vector (1st -> 2nd node) ! the first and last NodeIndx values point to the corresponding Joint nodes indices which are at the start of the Mesh pos1 = m%DispNodePosHst(:,mem%NodeIndx(i )) pos2 = m%DispNodePosHst(:,mem%NodeIndx(i+1)) - call GetOrientationAngles( pos1, pos2, phi, sinPhi, cosPhi, tanPhi, sinBeta, cosBeta, k_hat, errStat2, errMsg2 ); if (Failed()) return - ! Compute element to global DirCos matrix for undisplaced structure first - call Morison_DirCosMtrx( u%Mesh%Position(:,mem%NodeIndx(i )), u%Mesh%Position(:,mem%NodeIndx(i+1)), mem%MSpinOrient, CMatrix ) - ! Prepend body motion - Assuming the rotation of the starting node is representative of the whole element - CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,mem%NodeIndx(i))),CMatrix) + call GetElementAxialVec( pos1, pos2, k_hat, errStat2, errMsg2 ); if (Failed()) return + ! Compute total element orientation matrix assuming the rotation of the starting node is representative of the whole element + CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,mem%NodeIndx(i))),mem%CMatrix) CTrans = transpose(CMatrix) ! Note: CMatrix is element local to global displaced. CTrans is the opposite. ! save some commonly used variables @@ -3880,7 +3538,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, a_s2 = u%Mesh%TranslationAcc(:, mem%NodeIndx(i+1)) alpha_s2 = u%Mesh%RotationAcc (:, mem%NodeIndx(i+1)) omega_s2 = u%Mesh%RotationVel (:, mem%NodeIndx(i+1)) - IF (mem%MSecGeom == MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) r1In = mem%Rin(i ) ! outer radius at element nodes including marine growth r2In = mem%Rin(i+1) IF ( mem%floodstatus(i) == 1 ) THEN ! Fully flooded element @@ -3891,7 +3550,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, l = mem%h_fill/mem%dl rMidIn = r1In * (1.0-l) + r2In * l END IF - ELSE IF (mem%MSecGeom == MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) Sa1In = mem%Sain(i ) ! outer side A at element nodes including marine growth Sa2In = mem%Sain(i+1) Sb1In = mem%Sbin(i ) ! outer side B at element nodes including marine growth @@ -3906,18 +3565,19 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, SaMidIn = Sa1In * (1.0-l) + Sa2In * l SbMidIn = Sb1In * (1.0-l) + Sb2In * l END IF - END IF + END SELECT ! ------------------ flooded ballast inertia: sides: Section 6.1.1 : Always compute regardless of PropPot setting --------------------- ! lower node Imat = 0.0_ReKi - IF (mem%MSecGeom == MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) Imat(1,1) = mem%I_rfb_l(i) Imat(2,2) = mem%I_rfb_l(i) - ELSE IF (mem%MSecGeom == MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) Imat(1,1) = mem%I_xfb_l(i) Imat(2,2) = mem%I_yfb_l(i) - END IF + END SELECT Imat(3,3) = mem%I_lfb_l(i) Imat = matmul(matmul(CMatrix, Imat), CTrans) iArm = mem%h_cfb_l(i) * k_hat @@ -3931,13 +3591,14 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! upper node Imat = 0.0_ReKi - IF (mem%MSecGeom == MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) Imat(1,1) = mem%I_rfb_u(i) Imat(2,2) = mem%I_rfb_u(i) - ELSE IF (mem%MSecGeom == MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) Imat(1,1) = mem%I_xfb_u(i) Imat(2,2) = mem%I_yfb_u(i) - END IF + END SELECT Imat(3,3) = mem%I_lfb_u(i) Imat = matmul(matmul(CMatrix, Imat), CTrans) iArm = mem%h_cfb_u(i) * k_hat @@ -3952,7 +3613,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! ------------------ flooded ballast weight : sides : Section 5.1.2 & 5.2.2 : Always compute regardless of PropPot setting --------------------- F_B1 = 0.0 F_B2 = 0.0 - IF (mem%MSecGeom == MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) F_B1(3) = - p%gravity * mem%m_fb_l(i) F_B1(1:3) = F_B1(1:3) + mem%FillDens * p%gravity * pi * ( rMidIn*rMidIn*(zMid-zFillGroup) - r1In*r1In*(z1-zFillGroup) ) * k_hat F_B1(4:6) = -( p%gravity * mem%m_fb_l(i) * mem%h_cfb_l(i) + mem%FillDens * p%gravity * 0.25*pi*(rMidIn**4-r1In**4) ) * Cross_Product(k_hat,(/0.0,0.0,1.0/)) @@ -3964,7 +3626,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, F_B1(1:3) = F_B1(1:3) + mem%FillDens * p%gravity * pi * rMidIn**2* (zFillGroup - zMid) * k_hat F_B1(4:6) = F_B1(4:6) + mem%FillDens * p%gravity * 0.25 * pi * rMidIn**4* Cross_Product(k_hat,(/0.0,0.0,1.0/)) END IF - ELSE IF (mem%MSecGeom == MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) CALL GetSectionUnitVectors_Rec( CMatrix, x_hat, y_hat ) F_B1(3) = - p%gravity * mem%m_fb_l(i) F_B1(1:3) = F_B1(1:3) + mem%FillDens * p%gravity * ( SaMidIn*SbMidIn*(zMid-zFillGroup) - Sa1In*Sb1In*(z1-zFillGroup) ) * k_hat @@ -3981,7 +3643,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, F_B1(1:3) = F_B1(1:3) + mem%FillDens * p%gravity * SaMidIn*SbMidIn*(zFillGroup-zMid) * k_hat F_B1(4:6) = F_B1(4:6) + mem%FillDens * p%gravity / 12.0 * (SaMidIn**3*SbMidIn*x_hat(3)*y_hat - SaMidIn*SbMidIn**3*y_hat(3)*x_hat) END IF - END IF + END SELECT m%memberLoads(im)%F_BF(:, i ) = m%memberLoads(im)%F_BF(:, i ) + F_B1 m%memberLoads(im)%F_BF(:, i+1) = m%memberLoads(im)%F_BF(:, i+1) + F_B2 @@ -3997,10 +3659,10 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, !-----------------------------------------------------------------------------------------------------! ! External Hydrodynamic Side Loads - Start ! !-----------------------------------------------------------------------------------------------------! - IF ( p%WaveField%WaveStMod > 0 .AND. MemSubStat == 1 .AND. (m%NodeInWater(mem%NodeIndx(N+1)).EQ.0_IntKi) ) THEN + IF ( (p%WaveField%WaveStMod > 0 .OR. p%WaveDisp == 1_IntKi) .AND. MemSubStat == 1 .AND. (m%NodeInWater(mem%NodeIndx(N+1)).EQ.0_IntKi) ) THEN !----------------------------Apply load smoothing----------------------------! ! only when: - ! 1. wave stretching is enabled + ! 1. wave stretching is enabled or WaveDisp=1. Either can lead to member nodes going in and out of water ! 2. member centerline crosses the free surface exactly once ! 3. the last node is out of water, which implies the first node is in water @@ -4012,9 +3674,18 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, pos1 = m%DispNodePosHdn(:,mem%NodeIndx(i )) pos2 = m%DispNodePosHdn(:,mem%NodeIndx(i+1)) - ! Free surface elevation above or below node i and i+1 - Zeta1 = m%WaveElev(mem%NodeIndx(i)) - Zeta2 = m%WaveElev(mem%NodeIndx(i+1)) + ! Free surface elevation above or below node i and i+1. + ! Without wave stretching the free surface is the SWL (z=0), consistent with how nodeInWater/MemSubStat are determined. + IF ( p%WaveField%WaveStMod > 0_IntKi ) THEN + Zeta1 = m%WaveElev(mem%NodeIndx(i )) + Zeta2 = m%WaveElev(mem%NodeIndx(i+1)) + ELSE + Zeta1 = 0.0_ReKi + Zeta2 = 0.0_ReKi + END IF + + ! Rotate member properties based on local node orientation + call RotateMemberNode(p, im, mem, u%PtfmRefY, u%Mesh%Orientation(:,:,mem%NodeIndx(i)), ErrStat2, ErrMsg2); if (Failed()) return ! Compute deltal and h_c IF ( i == 1 ) THEN ! First node @@ -4040,7 +3711,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, END IF ! Compute the slope of member radius/side length - IF (mem%MSecGeom==MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) IF (i == 1) THEN dRdl_p = abs(mem%dRdl_mg(i)) dRdl_pp = mem%dRdl_mg(i) @@ -4051,7 +3723,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, dRdl_p = abs(mem%dRdl_mg(N)) dRdl_pp = mem%dRdl_mg(N) END IF - ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) IF (i == 1) THEN dSadl_p = abs(mem%dSadl_mg(i)) dSadl_pp = mem%dSadl_mg(i) @@ -4068,17 +3740,18 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, dSbdl_p = abs(mem%dSbdl_mg(N)) dSbdl_pp = mem%dSbdl_mg(N) END IF - END IF + END SELECT !-------------------- hydrodynamic drag loads: sides: Section 7.1.2 ------------------------! vec = matmul( mem%Ak,m%vrel(:,mem%NodeIndx(i)) ) - IF (mem%MSecGeom==MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) f_hydro = mem%Cd(i)*p%WaveField%WtrDens*mem%RMG(i)*TwoNorm(vec)*vec + & ! radial part 0.5*mem%AxCd(i)*p%WaveField%WtrDens * pi*mem%RMG(i)*dRdl_p * & ! axial part abs(dot_product( mem%k, m%vrel(:,mem%NodeIndx(i)) )) * matmul( mem%kkt, m%vrel(:,mem%NodeIndx(i)) ) ! axial part cont'd - ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN - Call GetDistDrag_Rec(p, m, u, xd, Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - END IF + CASE (MSecGeom_Rec) + Call GetDistDrag_Rec(p, m, u, xd, Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat2,ErrMsg2); if (Failed()) return + END SELECT CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_D(:, i) ) y%Mesh%Force (:,mem%NodeIndx(i)) = y%Mesh%Force (:,mem%NodeIndx(i)) + m%memberLoads(im)%F_D(1:3, i) y%Mesh%Moment(:,mem%NodeIndx(i)) = y%Mesh%Moment(:,mem%NodeIndx(i)) + m%memberLoads(im)%F_D(4:6, i) @@ -4088,45 +3761,27 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, IF ( .NOT. mem%PropPot ) THEN !-------------------- hydrodynamic added mass loads: sides: Section 7.1.3 ------------------------! - IF (mem%MSecGeom==MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) Am = mem%Ca(i)*p%WaveField%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*mem%Ak + 2.0*mem%AxCa(i)*p%WaveField%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*dRdl_p*mem%kkt f_hydro = -matmul( Am, u%Mesh%TranslationAcc(:,mem%NodeIndx(i)) ) - ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) f_hydro = -p%WaveField%WtrDens*mem%CaB(i) * 0.25*pi*mem%SbMG(i)*mem%SbMG(i) * Dot_Product(u%Mesh%TranslationAcc(:,mem%NodeIndx(i)),mem%x_hat)*mem%x_hat & -p%WaveField%WtrDens*mem%CaA(i) * 0.25*pi*mem%SaMG(i)*mem%SaMG(i) * Dot_Product(u%Mesh%TranslationAcc(:,mem%NodeIndx(i)),mem%y_hat)*mem%y_hat & -0.5*p%WaveField%WtrDens*mem%AxCa(i) * (dSbdl_p*mem%SaMG(i)+dSadl_p*mem%SbMG(i))*SQRT(mem%SaMG(i)*mem%SbMG(i)) * Dot_Product(u%Mesh%TranslationAcc(:,mem%NodeIndx(i)),mem%k)*mem%k + END SELECT + ! Compute added-mass force up to the instantaneous free surface + f_hydro = f_hydro * m%nodeInWater(mem%NodeIndx(i)) ! Zero the force if node above free surface + CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_A(:, i) ) + IF (i == FSElem) THEN ! Save the distributed load at the first node below the free surface + F_A0 = f_hydro END IF - IF ( p%AMMod .EQ. 0_IntKi ) THEN ! Compute added-mass force up to the SWL - z1 = u%Mesh%Position(3, mem%NodeIndx(i)) - p%WaveField%MSL2SWL ! Undisplaced z-position of the current node - IF ( z1 > 0.0_ReKi ) THEN ! Node is above SWL undisplaced; zero added-mass force - f_hydro = 0.0_ReKi - CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_A(:, i) ) - ELSE - ! Need to compute deltal_AM and h_c_AM based on the formulation without wave stretching. - z2 = u%Mesh%Position(3, mem%NodeIndx(i+1)) - p%WaveField%MSL2SWL ! Undisplaced z-position of the next node - IF ( z2 > 0.0_ReKi ) THEN ! Element i crosses the SWL - h = -z1 / mem%cosPhi_ref ! Length of Element i between SWL and node i, h>=0 - deltal_AM = mem%dl/2.0 + h - h_c_AM = 0.5*(h-mem%dl/2.0) - ELSE - deltal_AM = deltal; - h_c_AM = h_c - END IF - ! Note: Do not overwrite deltal and h_c here. Still need them for the fluid inertia and drag forces. - CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal_AM, h_c_AM, m%memberLoads(im)%F_A(:, i) ) - END IF - ELSE ! Compute added-mass force up to the instantaneous free surface - f_hydro = f_hydro * m%nodeInWater(mem%NodeIndx(i)) ! Zero the force if node above free surface - CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_A(:, i) ) - IF (i == FSElem) THEN ! Save the distributed load at the first node below the free surface - F_A0 = f_hydro - END IF - END IF ! AMMod 0 or 1 y%Mesh%Force (:,mem%NodeIndx(i)) = y%Mesh%Force (:,mem%NodeIndx(i)) + m%memberLoads(im)%F_A(1:3, i) y%Mesh%Moment(:,mem%NodeIndx(i)) = y%Mesh%Moment(:,mem%NodeIndx(i)) + m%memberLoads(im)%F_A(4:6, i) !--------------------- hydrodynamic inertia loads: sides: Section 7.1.4 --------------------------! - IF (mem%MSecGeom==MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) IF (mem%PropMCF) THEN f_hydro= p%WaveField%WtrDens*pi*mem%RMG(i)*mem%RMG(i) * matmul( mem%Ak, m%FAMCF(:,mem%NodeIndx(i)) ) + & 2.0*mem%AxCa(i)*p%WaveField%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*dRdl_p * matmul( mem%kkt, m%FA(:,mem%NodeIndx(i)) ) + & @@ -4136,7 +3791,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, 2.0*mem%AxCa(i) *p%WaveField%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*dRdl_p * matmul( mem%kkt, m%FA(:,mem%NodeIndx(i)) ) + & 2.0*m%FDynP(mem%NodeIndx(i))*mem%AxCp(i)*pi*mem%RMG(i)*dRdl_pp*mem%k END IF - ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) ! Note: MacCamy-Fuchs correction cannot be applied to rectangular members f_hydro= mem%Cp(i)*p%WaveField%WtrDens* mem%SaMG(i)*mem%SbMG(i) * matmul( mem%Ak, m%FA(:,mem%NodeIndx(i)) ) + & ! transver FK component m%FDynP(mem%NodeIndx(i))*mem%AxCp(i)* (mem%SaMG(i)*dSbdl_pp+dSadl_pp*mem%SbMG(i)) *mem%k + & ! axial FK component @@ -4144,7 +3799,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, p%WaveField%WtrDens*mem%CaA(i) * 0.25*pi*mem%SaMG(i)*mem%SaMG(i) * Dot_Product(m%FA(:,mem%NodeIndx(i)),mem%y_hat)*mem%y_hat + & ! y-component of diffraction part 0.5*p%WaveField%WtrDens*mem%AxCa(i) * (dSbdl_p*mem%SaMG(i)+dSadl_p*mem%SbMG(i))*SQRT(mem%SaMG(i)*mem%SbMG(i)) * & ! axial component of diffraction part Dot_Product(m%FA(:,mem%NodeIndx(i)),mem%k)*mem%k ! axial component of diffraction part cont'd - END IF + END SELECT CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_I(:, i) ) y%Mesh%Force (:,mem%NodeIndx(i)) = y%Mesh%Force (:,mem%NodeIndx(i)) + m%memberLoads(im)%F_I(1:3, i) y%Mesh%Moment(:,mem%NodeIndx(i)) = y%Mesh%Moment(:,mem%NodeIndx(i)) + m%memberLoads(im)%F_I(4:6, i) @@ -4155,12 +3810,20 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, END DO ! i =1,N+1 ! loop through member nodes + IF ( FSElem < 1_IntKi ) THEN ! No free-surface-crossing element found; guard against invalid FSElem indexing below + ErrStat2 = ErrID_Fatal + ErrMsg2 = 'Failed to locate the free-surface-crossing element for Member ID '//trim(num2lstr(mem%MemberID))//'.' + if (Failed()) return + END IF + !----------------------------------------------------------------------------------------------------! ! Compute the distributed loads at the point of intersection between the member and the free surface ! !----------------------------------------------------------------------------------------------------! + ! Rotate member properties to the free-surface-crossing element (its starting node is representative) + call RotateMemberNode(p, im, mem, u%PtfmRefY, u%Mesh%Orientation(:,:,mem%NodeIndx(FSElem)), ErrStat2, ErrMsg2); if (Failed()) return ! Get wave kinematics at the free-surface intersection. Set forceNodeInWater=.TRUE. to guarantee the free-surface intersection is in water. CALL WaveField_GetNodeWaveKin( p%WaveField, m%WaveField_m, Time, FSInt, .TRUE., .TRUE., nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat2, ErrMsg2 ) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + if (Failed()) return FDynPFSInt = REAL(FDynP,ReKi) FVFSInt = REAL(FV, ReKi) FAFSInt = REAL(FA, ReKi) @@ -4179,7 +3842,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, (1.0-SubRatio) * u%Mesh%TranslationVel(:,mem%NodeIndx(FSElem )) & ) - IF (mem%MSecGeom==MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) dRdl_p = abs(mem%dRdl_mg(FSElem)) dRdl_pp = mem%dRdl_mg(FSElem) RMGFSInt = SubRatio * mem%RMG( FSElem+1) + (1.0-SubRatio) * mem%RMG( FSElem) @@ -4199,11 +3863,9 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, IF ( .NOT. mem%PropPot ) THEN ! ------------------- hydrodynamic added mass loads: sides: Section 7.1.3 ------------------------ - IF (p%AMMod > 0_IntKi) THEN - Am = CaFSInt*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt*mem%Ak + & - 2.0*AxCaFSInt*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt*dRdl_p*mem%kkt - F_AS = -matmul( Am, SAFSInt ) - END IF + Am = CaFSInt*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt*mem%Ak + & + 2.0*AxCaFSInt*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt*dRdl_p*mem%kkt + F_AS = -matmul( Am, SAFSInt ) ! ------------------- hydrodynamic inertia loads: sides: Section 7.1.4 ------------------------ IF ( mem%PropMCF) THEN @@ -4216,16 +3878,13 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, 2.0*AxCpFSInt *pi*RMGFSInt *dRdl_pp * FDynPFSInt*mem%k END IF END IF - ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) dSadl_p = abs(mem%dSadl_mg(FSElem)) dSadl_pp = mem%dSadl_mg(FSElem) dSbdl_p = abs(mem%dSbdl_mg(FSElem)) dSbdl_pp = mem%dSbdl_mg(FSElem) SaMGFSInt = SubRatio * mem%SaMG(FSElem+1) + (1.0-SubRatio) * mem%SaMG(FSElem) SbMGFSInt = SubRatio * mem%SbMG(FSElem+1) + (1.0-SubRatio) * mem%SbMG(FSElem) - ! CdAFSInt = SubRatio * mem%CdA( FSElem+1) + (1.0-SubRatio) * mem%CdA( FSElem) - ! CdBFSInt = SubRatio * mem%CdB( FSElem+1) + (1.0-SubRatio) * mem%CdB( FSElem) - ! AxCdFSInt = SubRatio * mem%AxCd(FSElem+1) + (1.0-SubRatio) * mem%AxCd(FSElem) CaAFSInt = SubRatio * mem%CaA( FSElem+1) + (1.0-SubRatio) * mem%CaA( FSElem) CaBFSInt = SubRatio * mem%CaB( FSElem+1) + (1.0-SubRatio) * mem%CaB( FSElem) AxCaFSInt = SubRatio * mem%AxCa(FSElem+1) + (1.0-SubRatio) * mem%AxCa(FSElem) @@ -4233,18 +3892,16 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, AxCpFSInt = SubRatio * mem%AxCp(FSElem+1) + (1.0-SubRatio) * mem%AxCp(FSElem) Call GetDistDrag_Rec(p, m, u, xd, Time,mem,FSElem,dSadl_p,dSbdl_p,F_DS,ErrStat2,ErrMsg2,SubRatio,vrelFSInt) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + if (Failed()) return ! Hydrodynamic added mass and inertia loads IF ( .NOT. mem%PropPot ) THEN ! ------------------- hydrodynamic added mass loads: sides: Section 7.1.3 ------------------------ - IF (p%AMMod > 0_IntKi) THEN - F_AS = -p%WaveField%WtrDens*CaBFSInt * 0.25*pi*SbMGFSInt*SbMGFSInt * Dot_Product(SAFSInt,mem%x_hat)*mem%x_hat & - -p%WaveField%WtrDens*CaAFSInt * 0.25*pi*SaMGFSInt*SaMGFSInt * Dot_Product(SAFSInt,mem%y_hat)*mem%y_hat & - -0.5*p%WaveField%WtrDens*AxCaFSInt * (dSbdl_p*SaMGFSInt+dSadl_p*SbMGFSInt)*SQRT(SaMGFSInt*SbMGFSInt) * Dot_Product(SAFSInt,mem%k)*mem%k - END IF - + F_AS = -p%WaveField%WtrDens*CaBFSInt * 0.25*pi*SbMGFSInt*SbMGFSInt * Dot_Product(SAFSInt,mem%x_hat)*mem%x_hat & + -p%WaveField%WtrDens*CaAFSInt * 0.25*pi*SaMGFSInt*SaMGFSInt * Dot_Product(SAFSInt,mem%y_hat)*mem%y_hat & + -0.5*p%WaveField%WtrDens*AxCaFSInt * (dSbdl_p*SaMGFSInt+dSadl_p*SbMGFSInt)*SQRT(SaMGFSInt*SbMGFSInt) * Dot_Product(SAFSInt,mem%k)*mem%k + ! ------------------- hydrodynamic inertia loads: sides: Section 7.1.4 ------------------------ F_IS= CpFSInt*p%WaveField%WtrDens* SaMGFSInt*SbMGFSInt * matmul( mem%Ak, FAFSInt ) + & ! transver FK component FDynPFSInt*AxCpFSInt* (SaMGFSInt*dSbdl_pp+dSadl_pp*SbMGFSInt) *mem%k + & ! axial FK component @@ -4254,7 +3911,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, Dot_Product(FAFSInt,mem%k)*mem%k ! axial component of diffraction part cont'd END IF - END IF + END SELECT !----------------------------------------------------------------------------------------------------! ! Perform the load redistribution for smooth time series ! !----------------------------------------------------------------------------------------------------! @@ -4287,26 +3944,24 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! Hydrodynamic added mass and inertia loads IF ( .NOT. mem%PropPot ) THEN - - IF ( p%AMMod > 0_IntKi ) THEN - !-------------------- hydrodynamic added mass loads: sides: Section 7.1.3 ------------------------! - ! Apply load redistribution to the first node below the free surface - Df_hydro = ((SubRatio-1.0_ReKi)/(2.0_ReKi)-f_redist)*F_A0 + SubRatio/2.0_ReKi*F_AS - CALL LumpDistrHydroLoads( Df_hydro, mem%k, deltal, h_c, Df_hydro_lumped) - m%memberLoads(im)%F_A(:, FSElem) = m%memberLoads(im)%F_A(:, FSElem) + Df_hydro_lumped - y%Mesh%Force (:,mem%NodeIndx(FSElem)) = y%Mesh%Force (:,mem%NodeIndx(FSElem)) + Df_hydro_lumped(1:3) - y%Mesh%Moment(:,mem%NodeIndx(FSElem)) = y%Mesh%Moment(:,mem%NodeIndx(FSElem)) + Df_hydro_lumped(4:6) - - ! Apply load redistribution to the second node below the free surface - IF (FSElem > 1_IntKi) THEN - Df_hydro = f_redist * F_A0 - CALL LumpDistrHydroLoads( Df_hydro, mem%k, deltal, h_c, Df_hydro_lumped) - m%memberLoads(im)%F_A(:, FSElem-1) = m%memberLoads(im)%F_A(:, FSElem-1) + Df_hydro_lumped - y%Mesh%Force (:,mem%NodeIndx(FSElem-1)) = y%Mesh%Force (:,mem%NodeIndx(FSElem-1)) + Df_hydro_lumped(1:3) - y%Mesh%Moment(:,mem%NodeIndx(FSElem-1)) = y%Mesh%Moment(:,mem%NodeIndx(FSElem-1)) + Df_hydro_lumped(4:6) - END IF + + !-------------------- hydrodynamic added mass loads: sides: Section 7.1.3 ------------------------! + ! Apply load redistribution to the first node below the free surface + Df_hydro = ((SubRatio-1.0_ReKi)/(2.0_ReKi)-f_redist)*F_A0 + SubRatio/2.0_ReKi*F_AS + CALL LumpDistrHydroLoads( Df_hydro, mem%k, deltal, h_c, Df_hydro_lumped) + m%memberLoads(im)%F_A(:, FSElem) = m%memberLoads(im)%F_A(:, FSElem) + Df_hydro_lumped + y%Mesh%Force (:,mem%NodeIndx(FSElem)) = y%Mesh%Force (:,mem%NodeIndx(FSElem)) + Df_hydro_lumped(1:3) + y%Mesh%Moment(:,mem%NodeIndx(FSElem)) = y%Mesh%Moment(:,mem%NodeIndx(FSElem)) + Df_hydro_lumped(4:6) + + ! Apply load redistribution to the second node below the free surface + IF (FSElem > 1_IntKi) THEN + Df_hydro = f_redist * F_A0 + CALL LumpDistrHydroLoads( Df_hydro, mem%k, deltal, h_c, Df_hydro_lumped) + m%memberLoads(im)%F_A(:, FSElem-1) = m%memberLoads(im)%F_A(:, FSElem-1) + Df_hydro_lumped + y%Mesh%Force (:,mem%NodeIndx(FSElem-1)) = y%Mesh%Force (:,mem%NodeIndx(FSElem-1)) + Df_hydro_lumped(1:3) + y%Mesh%Moment(:,mem%NodeIndx(FSElem-1)) = y%Mesh%Moment(:,mem%NodeIndx(FSElem-1)) + Df_hydro_lumped(4:6) END IF - + !-------------------- hydrodynamic inertia loads: sides: Section 7.1.4 --------------------------! ! Apply load redistribution to the first node below the free surface Df_hydro = ((SubRatio-1.0_ReKi)/(2.0_ReKi)-f_redist)*F_I0 + SubRatio/2.0_ReKi*F_IS @@ -4332,12 +3987,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, F_S = F_DS F_0 = F_D0 IF ( .NOT. mem%PropPot) THEN - F_S = F_S + F_IS - F_0 = F_0 + F_I0 - IF ( p%AMMod > 0_IntKi) THEN - F_S = F_S + F_AS - F_0 = F_0 + F_A0 - END IF + F_S = F_S + F_IS + F_AS + F_0 = F_0 + F_I0 + F_A0 END IF ! First node below the free surface DM_hydro = 0.5_ReKi * SubRatio**2 * deltal * cross_product(mem%k, F_S) @@ -4353,6 +4004,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, DO i = mem%i_floor+1,N+1 ! loop through member nodes starting from the first node above seabed z1 = m%DispNodePosHdn(3, mem%NodeIndx(i)) pos1 = m%DispNodePosHdn(:, mem%NodeIndx(i)) + ! Rotate member properties based on local node orientation + call RotateMemberNode(p, im, mem, u%PtfmRefY, u%Mesh%Orientation(:,:,mem%NodeIndx(i)), ErrStat2, ErrMsg2); if (Failed()) return !---------------------------------------------Compute deltal and h_c------------------------------------------! ! Cannot make any assumption about WaveStMod and member orientation IF ( m%NodeInWater(mem%NodeIndx(i)) .EQ. 0_IntKi ) THEN ! Node is out of water @@ -4407,7 +4060,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, END IF ! Compute the slope of member radius/side length - IF (mem%MSecGeom==MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) IF (i == 1) THEN dRdl_p = abs(mem%dRdl_mg(i)) dRdl_pp = mem%dRdl_mg(i) @@ -4418,7 +4072,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, dRdl_p = abs(mem%dRdl_mg(N)) dRdl_pp = mem%dRdl_mg(N) END IF - ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) IF (i == 1) THEN dSadl_p = abs(mem%dSadl_mg(i)) dSadl_pp = mem%dSadl_mg(i) @@ -4435,68 +4089,42 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, dSbdl_p = abs(mem%dSbdl_mg(N)) dSbdl_pp = mem%dSbdl_mg(N) END IF - END IF + END SELECT !--------------------- hydrodynamic drag loads: sides: Section 7.1.2 --------------------------------! vec = matmul( mem%Ak,m%vrel(:,mem%NodeIndx(i)) ) - IF (mem%MSecGeom==MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) f_hydro = mem%Cd(i)*p%WaveField%WtrDens*mem%RMG(i)*TwoNorm(vec)*vec + & ! radial part 0.5*mem%AxCd(i)*p%WaveField%WtrDens*pi*mem%RMG(i)*dRdl_p * & ! axial part abs(dot_product( mem%k, m%vrel(:,mem%NodeIndx(i)) )) * matmul( mem%kkt, m%vrel(:,mem%NodeIndx(i)) ) ! axial part cont'd - ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN - Call GetDistDrag_Rec(p, m, u, xd, Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat2,ErrMsg2) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - END IF + CASE (MSecGeom_Rec) + Call GetDistDrag_Rec(p, m, u, xd, Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat2,ErrMsg2); if (Failed()) return + END SELECT CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_D(:, i) ) y%Mesh%Force (:,mem%NodeIndx(i)) = y%Mesh%Force (:,mem%NodeIndx(i)) + m%memberLoads(im)%F_D(1:3, i) y%Mesh%Moment(:,mem%NodeIndx(i)) = y%Mesh%Moment(:,mem%NodeIndx(i)) + m%memberLoads(im)%F_D(4:6, i) IF ( .NOT. mem%PropPot ) THEN !-------------------- hydrodynamic added mass loads: sides: Section 7.1.3 ------------------------! - IF (mem%MSecGeom==MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) Am = mem%Ca(i)*p%WaveField%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*mem%Ak + 2.0*mem%AxCa(i)*p%WaveField%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*dRdl_p*mem%kkt f_hydro = -matmul( Am, u%Mesh%TranslationAcc(:,mem%NodeIndx(i)) ) - ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) f_hydro = -p%WaveField%WtrDens*mem%CaB(i) * 0.25*pi*mem%SbMG(i)*mem%SbMG(i) * Dot_Product(u%Mesh%TranslationAcc(:,mem%NodeIndx(i)),mem%x_hat)*mem%x_hat & -p%WaveField%WtrDens*mem%CaA(i) * 0.25*pi*mem%SaMG(i)*mem%SaMG(i) * Dot_Product(u%Mesh%TranslationAcc(:,mem%NodeIndx(i)),mem%y_hat)*mem%y_hat & -0.5*p%WaveField%WtrDens*mem%AxCa(i) * (dSbdl_p*mem%SaMG(i)+dSadl_p*mem%SbMG(i))*SQRT(mem%SaMG(i)*mem%SbMG(i)) * Dot_Product(u%Mesh%TranslationAcc(:,mem%NodeIndx(i)),mem%k)*mem%k - END IF - IF ( p%AMMod .EQ. 0_IntKi ) THEN ! Always compute added-mass force on nodes below SWL when undisplaced - z1 = u%Mesh%Position(3, mem%NodeIndx(i)) - p%WaveField%MSL2SWL ! Undisplaced z-position of the current node - IF ( z1 > 0.0_ReKi ) THEN ! Node is above SWL when undisplaced; zero added-mass force - f_hydro = 0.0_ReKi - CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_A(:, i) ) - ELSE ! Node at or below SWL when undisplaced - IF ( i == 1 ) THEN - deltalLeft = 0.0_ReKi - ELSE IF ( i == mem%i_floor+1 ) THEN - deltalLeft = -mem%h_floor - ELSE - deltalLeft = 0.5_ReKi * mem%dl - END IF - IF ( i == N+1 ) THEN - deltalRight = 0.0_ReKi - ELSE - z2 = u%Mesh%Position(3, mem%NodeIndx(i+1)) - p%WaveField%MSL2SWL - IF ( z2 > 0.0_ReKi ) THEN ! Element i crosses the SWL - deltalRight = -z1 / mem%cosPhi_ref - ELSE - deltalRight = 0.5_ReKi * mem%dl - END IF - END IF - deltal_AM = deltalRight + deltalLeft - h_c_AM = 0.5_ReKi * ( deltalRight - deltalLeft ) - CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal_AM, h_c_AM, m%memberLoads(im)%F_A(:, i) ) - END IF - ELSE ! Compute added-mass force on the instantaneous wetted section of the member - f_hydro = f_hydro * m%nodeInWater(mem%NodeIndx(i)) ! Zero the force if node above free surface - CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_A(:, i) ) - END IF ! AMMod 0 or 1 + END SELECT + ! Compute added-mass force on the instantaneous wetted section of the member + f_hydro = f_hydro * m%nodeInWater(mem%NodeIndx(i)) ! Zero the force if node above free surface + CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_A(:, i) ) y%Mesh%Force (:,mem%NodeIndx(i)) = y%Mesh%Force (:,mem%NodeIndx(i)) + m%memberLoads(im)%F_A(1:3, i) y%Mesh%Moment(:,mem%NodeIndx(i)) = y%Mesh%Moment(:,mem%NodeIndx(i)) + m%memberLoads(im)%F_A(4:6, i) !-------------------- hydrodynamic inertia loads: sides: Section 7.1.4 ---------------------------! - IF (mem%MSecGeom==MSecGeom_Cyl) THEN + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) IF ( mem%PropMCF ) THEN f_hydro= p%WaveField%WtrDens*pi*mem%RMG(i)*mem%RMG(i) * matmul( mem%Ak, m%FAMCF(:,mem%NodeIndx(i)) ) + & 2.0*mem%AxCa(i)*p%WaveField%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*dRdl_p * matmul( mem%kkt, m%FA(:,mem%NodeIndx(i)) ) + & @@ -4506,7 +4134,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, 2.0*mem%AxCa(i) *p%WaveField%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*dRdl_p * matmul( mem%kkt, m%FA(:,mem%NodeIndx(i)) ) + & 2.0*m%FDynP(mem%NodeIndx(i))*mem%AxCp(i)*pi*mem%RMG(i)*dRdl_pp*mem%k END IF - ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN + CASE (MSecGeom_Rec) ! Note: MacCamy-Fuchs correction cannot be applied to rectangular members f_hydro= mem%Cp(i)*p%WaveField%WtrDens* mem%SaMG(i)*mem%SbMG(i) * matmul( mem%Ak, m%FA(:,mem%NodeIndx(i)) ) + & ! transver FK component m%FDynP(mem%NodeIndx(i))*mem%AxCp(i)* (mem%SaMG(i)*dSbdl_pp+dSadl_pp*mem%SbMG(i)) *mem%k + & ! axial FK component @@ -4514,7 +4142,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, p%WaveField%WtrDens*mem%CaA(i) * 0.25*pi*mem%SaMG(i)*mem%SaMG(i) * Dot_Product(m%FA(:,mem%NodeIndx(i)),mem%y_hat)*mem%y_hat + & ! y-component of diffraction part 0.5*p%WaveField%WtrDens*mem%AxCa(i) * (dSbdl_p*mem%SaMG(i)+dSadl_p*mem%SbMG(i))*SQRT(mem%SaMG(i)*mem%SbMG(i)) * & ! axial component of diffraction part Dot_Product(m%FA(:,mem%NodeIndx(i)),mem%k)*mem%k ! axial component of diffraction part cont'd - END IF + END SELECT CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_I(:, i) ) y%Mesh%Force (:,mem%NodeIndx(i)) = y%Mesh%Force (:,mem%NodeIndx(i)) + m%memberLoads(im)%F_I(1:3, i) y%Mesh%Moment(:,mem%NodeIndx(i)) = y%Mesh%Moment(:,mem%NodeIndx(i)) + m%memberLoads(im)%F_I(4:6, i) @@ -4551,20 +4179,13 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! We need to subtract the MSL2SWL offset to place this in the SWL reference system pos1 = m%DispNodePosHst(:,mem%NodeIndx(1)) pos2 = m%DispNodePosHst(:,mem%NodeIndx(2)) - call GetOrientationAngles( pos1, pos2, phi1, sinPhi1, cosPhi1, tanPhi, sinBeta1, cosBeta1, k_hat1, errStat2, errMsg2 ) - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call GetElementAxialVec( pos1, pos2, k_hat1, errStat2, errMsg2 ); if (Failed()) return if ( N == 1 ) then ! Only one element in member - sinPhi2 = sinPhi1 - cosPhi2 = cosPhi1 - sinBeta2 = sinBeta1 - cosBeta2 = cosBeta1 - k_hat2 = k_hat1 + k_hat2 = k_hat1 else - ! We need to subtract the MSL2SWL offset to place this in the SWL reference system pos1 = m%DispNodePosHst(:, mem%NodeIndx(N )) pos2 = m%DispNodePosHst(:, mem%NodeIndx(N+1)) - call GetOrientationAngles( pos1, pos2, phi2, sinPhi2, cosPhi2, tanPhi, sinBeta2, cosBeta2, k_hat2, errStat2, errMsg2 ) - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call GetElementAxialVec( pos1, pos2, k_hat2, errStat2, errMsg2 ); if (Failed()) return end if ! z-coordinates of the two ends of the member z1 = m%DispNodePosHst(3,mem%NodeIndx( 1)) @@ -4572,10 +4193,9 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, if (mem%MSecGeom == MSecGeom_Rec) then ! Compute total orientation matrix of starting and ending joints - call Morison_DirCosMtrx( u%Mesh%Position(:,mem%NodeIndx(1)), u%Mesh%Position(:,mem%NodeIndx(N+1)), mem%MSpinOrient, CMatrix ) - CMatrix1 = matmul(transpose(u%Mesh%Orientation(:,:,mem%NodeIndx(1 ))),CMatrix) + CMatrix1 = matmul(transpose(u%Mesh%Orientation(:,:,mem%NodeIndx(1 ))),mem%CMatrix) CALL GetSectionUnitVectors_Rec( CMatrix1, x_hat1, y_hat1 ) - CMatrix2 = matmul(transpose(u%Mesh%Orientation(:,:,mem%NodeIndx(N+1))),CMatrix) + CMatrix2 = matmul(transpose(u%Mesh%Orientation(:,:,mem%NodeIndx(N+1))),mem%CMatrix) CALL GetSectionUnitVectors_Rec( CMatrix2, x_hat2, y_hat2 ) end if @@ -4584,24 +4204,26 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, if ( mem%memfloodstatus > 0 ) then if ( mem%i_floor == 0 ) then ! If the member is not buried in the seabed, compute the internal hydrostatic load on the starting endplate - if ( mem%MSecGeom == MSecGeom_Cyl ) then + SELECT CASE (mem%MSecGeom) + CASE (MSecGeom_Cyl) m%F_BF_End(1:3, mem%NodeIndx( 1)) = m%F_BF_End(1:3, mem%NodeIndx( 1)) - mem%FillDens * g * pi * mem%Rin( 1)**2* (zFillGroup - z1) * k_hat1 m%F_BF_End(4:6, mem%NodeIndx( 1)) = m%F_BF_End(4:6, mem%NodeIndx( 1)) - mem%FillDens * g * 0.25 * pi * mem%Rin( 1)**4* Cross_Product(k_hat1,(/0.0,0.0,1.0/)) - else if ( mem%MSecGeom == MSecGeom_Rec ) then + CASE (MSecGeom_Rec) m%F_BF_End(1:3, mem%NodeIndx( 1)) = m%F_BF_End(1:3, mem%NodeIndx( 1)) - mem%FillDens * g * mem%SaIn( 1) *mem%SbIn( 1)* (zFillGroup - z1) * k_hat1 m%F_BF_End(4:6, mem%NodeIndx( 1)) = m%F_BF_End(4:6, mem%NodeIndx( 1)) - mem%FillDens * g / 12.0 * (mem%SaIn( 1)**3*mem%SbIn( 1)*x_hat1(3)*y_hat1 - mem%SaIn(1)*mem%SbIn(1)**3*y_hat1(3)*x_hat1) - end if + END SELECT end if if ( (mem%i_floor 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 SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - FSPt = (/pos1(1),pos1(2),Zeta1/) ! Reference point on the free surface - ELSE - FSPt = (/pos1(1),pos1(2),0.0_ReKi/) - n_hat = (/0.0,0.0,1.0/) - END IF - - if (mem%MSecGeom==MSecGeom_Cyl) then - CALL GetSectionUnitVectors_Cyl( k_hat1, y_hat, z_hat ) - CALL GetSectionFreeSurfaceIntersects_Cyl( REAL(pos1,DbKi), REAL(FSPt,DbKi), k_hat1, y_hat, z_hat, n_hat, REAL(r1,DbKi), theta1, theta2, secStat) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetEndPlateHstLds_Cyl(p, pos1, k_hat1, y_hat, z_hat, r1, theta1, theta2, F_B_End) - IF (mem%MHstLMod == 1) THEN ! Check for partially wetted end plates - IF ( .NOT.( EqualRealNos((theta2-theta1),0.0_DbKi) .OR. EqualRealNos((theta2-theta1),2.0_DbKi*PI_D) ) ) THEN - CALL SetErrStat(ErrID_Warn, 'End plate is partially wetted with MHstLMod = 1. The buoyancy load and distribution potentially have large error. This has happened to the first node of Member ID ' //trim(num2lstr(mem%MemberID)), errStat, errMsg, RoutineName ) - END IF - END IF - else if (mem%MSecGeom==MSecGeom_Rec) then - CALL GetSectionUnitVectors_Rec( CMatrix1, x_hat, y_hat ) - CALL GetEndPlateHstLds_Rec(p, pos1, k_hat1, x_hat, y_hat, Sa1, Sb1, FSPt, n_hat, F_B_End) - end if - 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%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 SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - FSPt = (/pos2(1),pos2(2),Zeta2/) ! Reference point on the free surface - ELSE - FSPt = (/pos2(1),pos2(2),0.0_ReKi/) - n_hat = (/0.0,0.0,1.0/) - END IF - - if (mem%MSecGeom==MSecGeom_Cyl) then - CALL GetSectionUnitVectors_Cyl( k_hat2, y_hat, z_hat ) - CALL GetSectionFreeSurfaceIntersects_Cyl( REAL(pos2,DbKi), REAL(FSPt,DbKi), k_hat2, y_hat, z_hat, n_hat, REAL(r2,DbKi), theta1, theta2, secStat) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetEndPlateHstLds_Cyl(p, pos2, k_hat2, y_hat, z_hat, r2, theta1, theta2, F_B_End) - IF (mem%MHstLMod == 1) THEN ! Check for partially wetted end plates - IF ( .NOT.( EqualRealNos((theta2-theta1),0.0_DbKi) .OR. EqualRealNos((theta2-theta1),2.0_DbKi*PI_D) ) ) THEN - CALL SetErrStat(ErrID_Warn, 'End plate is partially wetted with MHstLMod = 1. The buoyancy load and distribution potentially have large error. This has happened to the last node of Member ID ' //trim(num2lstr(mem%MemberID)), errStat, errMsg, RoutineName ) - END IF - END IF - else if (mem%MSecGeom==MSecGeom_Rec) then - CALL GetSectionUnitVectors_Rec( CMatrix2, x_hat, y_hat ) - CALL GetEndPlateHstLds_Rec(p, pos2, k_hat2, x_hat, y_hat, Sa2, Sb2, FSPt, n_hat, F_B_End) - end if - m%F_B_End(:, mem%NodeIndx(N+1)) = m%F_B_End(:, mem%NodeIndx(N+1)) - F_B_End - - 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%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 SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - FSPt = (/pos2(1),pos2(2),Zeta2/) ! Reference point on the free surface - ELSE - FSPt = (/pos2(1),pos2(2),0.0_ReKi/) - n_hat = (/0.0,0.0,1.0/) - END IF - - if (mem%MSecGeom==MSecGeom_Cyl) then - CALL GetSectionUnitVectors_Cyl( k_hat2, y_hat, z_hat ) - CALL GetSectionFreeSurfaceIntersects_Cyl( REAL(pos2,DbKi), REAL(FSPt,DbKi), k_hat2, y_hat, z_hat, n_hat, REAL(r2,DbKi), theta1, theta2, secStat) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetEndPlateHstLds_Cyl(p, pos2, k_hat2, y_hat, z_hat, r2, theta1, theta2, F_B_End) - IF (mem%MHstLMod == 1) THEN ! Check for partially wetted end plates - IF ( .NOT.( EqualRealNos((theta2-theta1),0.0_DbKi) .OR. EqualRealNos((theta2-theta1),2.0_DbKi*PI_D) ) ) THEN - CALL SetErrStat(ErrID_Warn, 'End plate is partially wetted with MHstLMod = 1. The buoyancy load and distribution potentially have large error. This has happened to the last node of Member ID ' //trim(num2lstr(mem%MemberID)), errStat, errMsg, RoutineName ) - END IF - END IF - else if (mem%MSecGeom==MSecGeom_Rec) then - CALL GetSectionUnitVectors_Rec( CMatrix2, x_hat, y_hat ) - CALL GetEndPlateHstLds_Rec(p, pos2, k_hat2, x_hat, y_hat, Sa2, Sb2, FSPt, n_hat, F_B_End) - end if - m%F_B_End(:, mem%NodeIndx(N+1)) = m%F_B_End(:, mem%NodeIndx(N+1)) - F_B_End - + END SELECT + if (mem%i_floor == 0) then ! both ends at or above the seabed: load on both end plates + call AddEndBuoyancy(mem%NodeIndx( 1), pos1, k_hat1, r1, Sa1, Sb1, CMatrix1, 1.0_ReKi, 'first node') + call AddEndBuoyancy(mem%NodeIndx(N+1), pos2, k_hat2, r2, Sa2, Sb2, CMatrix2, -1.0_ReKi, 'last node') + elseif ( mem%doEndBuoyancy ) then ! seabed-crossing member: load only on the upper (N+1) end plate + call AddEndBuoyancy(mem%NodeIndx(N+1), pos2, k_hat2, r2, Sa2, Sb2, CMatrix2, -1.0_ReKi, 'last node') ! else ! entire member is buried below the seabed end if @@ -4739,9 +4281,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! Effect of wave stretching already baked into m%FDynP, m%FA, and m%vrel. No additional modification needed. - ! Joint yaw offset - call YawJoint(p, J,u%PtfmRefY,AM_End,An_End,DP_Const_End,I_MG_End,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + ! Rotate joint properties + call RotateJoint(p, J,u%PtfmRefY,u%Mesh%Orientation(:,:,J),AM_End,An_End,DP_Const_End,I_MG_End,ErrStat2,ErrMsg2); if (Failed()) return ! Lumped added mass loads qdotdot = reshape((/u%Mesh%TranslationAcc(:,J),u%Mesh%RotationAcc(:,J)/),(/6/)) @@ -4799,13 +4340,49 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! map the motion to the visulization mesh if (p%VisMeshes) then - !FIXME: error handling is incorrect here (overwrites all previous errors/warnings) - call Transfer_Point_to_Line2( u%Mesh, y%VisMesh, m%VisMeshMap, ErrStat, ErrMsg ) + call Transfer_Point_to_Line2( u%Mesh, y%VisMesh, m%VisMeshMap, ErrStat2, ErrMsg2 ); if (Failed()) return endif CONTAINS + subroutine AddEndBuoyancy(iNode, pos, k_hatE, rE, SaE, SbE, CMatrixE, sgn, nodeLabel) + ! Accumulate the external hydrostatic (buoyancy) load on one member end plate. + integer(IntKi), intent(in) :: iNode ! mesh node index of the end plate + real(ReKi), intent(in) :: pos(3), k_hatE(3) ! end node position and axial unit vector + real(ReKi), intent(in) :: rE, SaE, SbE ! scaled end radius (Cyl) or side lengths (Rec) + real(ReKi), intent(in) :: CMatrixE(3,3) ! end section orientation (Rec only) + real(ReKi), intent(in) :: sgn ! +1 at the lower end, -1 at the upper end + character(*), intent(in) :: nodeLabel ! node description used in warning messages + real(ReKi) :: FSPtL(3), n_hatL(3), ZetaL, F_B_EndL(6) + real(DbKi) :: th1, th2 + integer(IntKi) :: secStatL + + if (p%HstMod > 0_IntKi) then + call GetTotalWaveElev(p, m, Time, pos, ZetaL, ErrStat2, ErrMsg2); if (Failed()) return + call GetFreeSurfaceNormal(p, m, Time, pos, n_hatL, ErrStat2, ErrMsg2); if (Failed()) return + FSPtL = (/pos(1),pos(2),ZetaL/) ! Reference point on the free surface + else + FSPtL = (/pos(1),pos(2),0.0_ReKi/) + n_hatL = (/0.0,0.0,1.0/) + end if + + select case (mem%MSecGeom) + case (MSecGeom_Cyl) + call GetSectionUnitVectors_Cyl( k_hatE, y_hat, z_hat ) + call GetSectionFreeSurfaceIntersects_Cyl( REAL(pos,DbKi), REAL(FSPtL,DbKi), k_hatE, y_hat, z_hat, n_hatL, REAL(rE,DbKi), th1, th2, secStatL) + call GetEndPlateHstLds_Cyl(p, pos, k_hatE, y_hat, z_hat, rE, th1, th2, F_B_EndL) + if (mem%MHstLMod == 1 .and. secStatL == 1) then ! secStatL == 1: end plate crosses the free surface (partially wetted) + call SetErrStat(ErrID_Warn, 'End plate is partially wetted with MHstLMod = 1. The buoyancy load and distribution potentially have large error. This has happened to the '//trim(nodeLabel)//' of Member ID '//trim(num2lstr(mem%MemberID)), errStat, errMsg, RoutineName ) + end if + case (MSecGeom_Rec) + call GetSectionUnitVectors_Rec( CMatrixE, x_hat, y_hat ) + call GetEndPlateHstLds_Rec(p, pos, k_hatE, x_hat, y_hat, SaE, SbE, FSPtL, n_hatL, F_B_EndL) + end select + + m%F_B_End(:, iNode) = m%F_B_End(:, iNode) + sgn * F_B_EndL + end subroutine AddEndBuoyancy + logical function Failed() CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Failed = ErrStat >= AbortErrLev @@ -6061,10 +5638,11 @@ SUBROUTINE getElementHstLds_Mod1(p, m, mem, Time, pos1, pos2, Zeta1, Zeta2, k_ha END IF END SUBROUTINE getElementHstLds_Mod1 - SUBROUTINE YawJoint(p, JointNo, PtfmRefY, AM_End, An_End, DP_Const_End, I_MG_End, ErrStat, ErrMsg) + SUBROUTINE RotateJoint(p, JointNo, PtfmRefY, JOrientation, AM_End, An_End, DP_Const_End, I_MG_End, ErrStat, ErrMsg) TYPE(Morison_ParameterType), INTENT( IN ) :: p Integer(IntKi), intent(in ) :: JointNo Real(ReKi), intent(in ) :: PtfmRefY + Real(R8Ki), intent(in ) :: JOrientation(3,3) Real(ReKi), intent( out) :: AM_End(3,3) Real(ReKi), intent( out) :: An_End(3) Real(ReKi), intent( out) :: DP_Const_End(3) @@ -6072,24 +5650,29 @@ SUBROUTINE YawJoint(p, JointNo, PtfmRefY, AM_End, An_End, DP_Const_End, I_MG_End Integer(IntKi), intent( out) :: ErrStat Character(*), intent( out) :: ErrMsg - Integer(IntKi) :: ErrStat2 - Character(ErrMsgLen) :: ErrMsg2 + Real(ReKi) :: Rg2b(3,3) + Real(ReKi) :: Rb2g(3,3) - Character(*), parameter :: RoutineName = 'YawJoint' + Character(*), parameter :: RoutineName = 'RotateJoint' ErrStat = ErrID_None ErrMsg = '' - call hiFrameTransform(h2i,PtfmRefY,p%AM_End(:,:,jointNo),AM_End,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - call hiFrameTransform(h2i,PtfmRefY,p%An_End(:,jointNo),An_End,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - call hiFrameTransform(h2i,PtfmRefY,p%DP_Const_End(:,jointNo),DP_Const_End,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - call hiFrameTransform(h2i,PtfmRefY,p%I_MG_End(:,:,jointNo),I_MG_End,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + select case (p%WaveDisp) + case (0) + call GetPtfmRefYOrient(PtfmRefY, Rg2b, ErrStat, ErrMsg) + if (ErrStat /= ErrID_None) return + case (1) + Rg2b = JOrientation + end select + Rb2g = transpose(Rg2b) - END SUBROUTINE YawJoint + AM_End = matmul(matmul(Rb2g,p%AM_End(:,:,jointNo)),Rg2b) + An_End = matmul(Rb2g,p%An_End(:,jointNo)) + DP_Const_End = matmul(Rb2g,p%DP_Const_End(:,jointNo)) + I_MG_End = matmul(matmul(Rb2g,p%I_MG_End(:,:,jointNo)),Rg2b) + + END SUBROUTINE RotateJoint SUBROUTINE getMemBallastHiPt(p, m, u, member, z_hi, ErrStat, ErrMsg) ! This subroutine returns the highest point of a member's internal ballast @@ -6102,7 +5685,7 @@ SUBROUTINE getMemBallastHiPt(p, m, u, member, z_hi, ErrStat, ErrMsg) Character(*), intent( out) :: ErrMsg Integer(IntKi) :: N, elemNo - Real(ReKi) :: CMatrix0(3,3), CMatrix(3,3) + Real(ReKi) :: CMatrix(3,3) Real(ReKi) :: k_hat(3), x_hat(3), y_hat(3), z_hat(3) Real(ReKi) :: l, rIn, SaIn, SbIn, z0, pos1(3), pos2(3) @@ -6149,13 +5732,11 @@ SUBROUTINE getMemBallastHiPt(p, m, u, member, z_hi, ErrStat, ErrMsg) z_hi = MAX( pos2(3) + rIn * z_hat(3), z_hi) END IF ELSE IF (member%MSecGeom == MSecGeom_Rec) THEN - ! DirCos matrix of undisplaced member - CALL Morison_DirCosMtrx( u%Mesh%Position(:,member%NodeIndx(1 )), u%Mesh%Position(:,member%NodeIndx(N+1)), member%MSpinOrient, CMatrix0 ) ! Check the vertices of the starting section pos1 = m%DispNodePosHst(:,member%NodeIndx(1)) pos2 = m%DispNodePosHst(:,member%NodeIndx(2)) z0 = pos1(3) - CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,member%NodeIndx(1 ))),CMatrix0) + CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,member%NodeIndx(1 ))),member%CMatrix) CALL GetSectionUnitVectors_Rec( CMatrix, x_hat, y_hat ) SaIn = member%SaIn(1) SbIn = member%SbIn(1) @@ -6168,7 +5749,7 @@ SUBROUTINE getMemBallastHiPt(p, m, u, member, z_hi, ErrStat, ErrMsg) pos1 = m%DispNodePosHst(:,member%NodeIndx(N )) pos2 = m%DispNodePosHst(:,member%NodeIndx(N+1)) z0 = pos2(3) - CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,member%NodeIndx(N+1))),CMatrix0) + CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,member%NodeIndx(N+1))),member%CMatrix) CALL GetSectionUnitVectors_Rec( CMatrix, x_hat, y_hat ) SaIn = member%SaIn(N+1) SbIn = member%SbIn(N+1) @@ -6181,9 +5762,9 @@ SUBROUTINE getMemBallastHiPt(p, m, u, member, z_hi, ErrStat, ErrMsg) pos1 = m%DispNodePosHst(:,member%NodeIndx(elemNo )) pos2 = m%DispNodePosHst(:,member%NodeIndx(elemNo+1)) if ( member%h_fill>0.5*member%dl ) then - CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,member%NodeIndx(elemNo+1))),CMatrix0) + CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,member%NodeIndx(elemNo+1))),member%CMatrix) else - CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,member%NodeIndx(elemNo ))),CMatrix0) + CMatrix = matmul(transpose(u%Mesh%Orientation(:,:,member%NodeIndx(elemNo ))),member%CMatrix) end if CALL GetSectionUnitVectors_Rec( CMatrix, x_hat, y_hat ) k_hat = Cross_Product(x_hat,y_hat) @@ -6385,54 +5966,6 @@ subroutine LumpDistrHydroLoads( f_hydro, k_hat, dl, h_c, lumpedLoad ) lumpedLoad(1:3) = f_hydro*dl lumpedLoad(4:6) = cross_product(k_hat*h_c, f_hydro)*dl end subroutine LumpDistrHydroLoads -!---------------------------------------------------------------------------------------------------------------------------------- -! Takes loads on node i in element tilted frame and converts to 6DOF loads at node i and adjacent node -PURE SUBROUTINE DistributeElementLoads(Fl, Fr, M, sinPhi, cosPhi, SinBeta, cosBeta, alpha, F1, F2) - - REAL(ReKi), INTENT ( IN ) :: Fl ! (N) axial load about node i - REAL(ReKi), INTENT ( IN ) :: Fr ! (N) radial load about node i in direction of tilt - REAL(ReKi), INTENT ( IN ) :: M ! (N-m) radial moment about node i, positive in direction of tilt angle - REAL(ReKi), INTENT ( IN ) :: sinPhi ! trig functions of tilt angle - REAL(ReKi), INTENT ( IN ) :: cosPhi - REAL(ReKi), INTENT ( IN ) :: sinBeta ! trig functions of heading of tilt - REAL(ReKi), INTENT ( IN ) :: cosBeta - REAL(ReKi), INTENT ( IN ) :: alpha ! fraction of load staying with node i (1-alpha goes to other node) - - REAL(ReKi), INTENT ( OUT ) :: F1(6) ! (N, Nm) force/moment vector for node i - REAL(ReKi), INTENT ( OUT ) :: F2(6) ! (N, Nm) force/moment vector for the other node (whether i+1, or i-1) - REAL(ReKi) :: F(6) - - F(1) = cosBeta*(Fl*sinPhi + Fr*cosPhi) - F(2) = sinBeta*(Fl*sinPhi + Fr*cosPhi) - F(3) = (Fl*cosPhi - Fr*sinPhi) - F(4) = -sinBeta * M - F(5) = cosBeta * M - F(6) = 0.0 - - F1 = F*alpha - F2 = F*(1.0_ReKi-alpha) - -END SUBROUTINE DistributeElementLoads -!---------------------------------------------------------------------------------------------------------------------------------- -! Takes loads on end node i and converts to 6DOF loads, adding to the nodes existing loads -PURE SUBROUTINE AddEndLoad(Fl, M, sinPhi, cosPhi, SinBeta, cosBeta, Fi) - - REAL(ReKi), INTENT ( IN ) :: Fl ! (N) axial load about node i - REAL(ReKi), INTENT ( IN ) :: M ! (N-m) radial moment about node i, positive in direction of tilt angle - REAL(ReKi), INTENT ( IN ) :: sinPhi ! trig functions of tilt angle - REAL(ReKi), INTENT ( IN ) :: cosPhi - REAL(ReKi), INTENT ( IN ) :: sinBeta ! trig functions of heading of tilt - REAL(ReKi), INTENT ( IN ) :: cosBeta - REAL(ReKi), INTENT ( INOUT ) :: Fi(6) ! (N, Nm) force/moment vector for end node i - - Fi(1) = Fi(1) + Fl*sinPhi*cosBeta - Fi(2) = Fi(2) + Fl*sinPhi*sinBeta - Fi(3) = Fi(3) + Fl*cosPhi - Fi(4) = Fi(4) - M*sinBeta - Fi(5) = Fi(5) + M*cosBeta - -END SUBROUTINE AddEndLoad - !---------------------------------------------------------------------------------------------------------------------------------- !> Tight coupling routine for updating discrete states @@ -6453,6 +5986,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat INTEGER(IntKi) :: I, J, im, N INTEGER(IntKi) :: nodeInWater, tmpInt REAL(ReKi) :: pos(3), vrel(3), FV(3), vmag, vmagf, An_End(3) + REAL(ReKi) :: Rg2b(3,3) REAL(ReKi) :: posFC(3), SVFC(3), vrelFC, vrelFCf REAL(SiKi) :: FVTmp(3),FATmp(3) TYPE(Morison_MemberType) :: mem !< Current member @@ -6466,7 +6000,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat errMsg = "" - CALL GetDisplacedNodePosition( u, p, .FALSE., m%DispNodePosHdn ) ! For hydrodynamic loads; depends on WaveDisp and WaveStMod + CALL GetDisplacedNodePosition( u, p, .FALSE., m%DispNodePosHdn ) ! For hydrodynamic loads; depends on WaveDisp ! Update state of the relative normal velocity high-pass filter at each joint DO J = 1, p%NJoints @@ -6476,13 +6010,19 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat ! Get fluid velocity at the joint CALL WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, pos, .FALSE., .TRUE., nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + if (Failed()) return FV = REAL(FVTmp, ReKi) vrel = ( FV - u%Mesh%TranslationVel(:,J) ) * nodeInWater - ! Transform An_End based on reference yaw offset - call hiFrameTransform(h2i,u%PtfmRefY,p%An_End(:,j),An_End,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + ! Rotate An_End based on p%WaveDisp: reference yaw only (0) or full instantaneous orientation (1) + select case (p%WaveDisp) + case (0) + call GetPtfmRefYOrient(u%PtfmRefY, Rg2b, ErrStat2, ErrMsg2) + if (Failed()) return + case (1) + Rg2b = u%Mesh%Orientation(:,:,J) + end select + An_End = matmul(transpose(Rg2b),p%An_End(:,j)) ! Compute the dot product of the relative velocity vector with the directional Area of the Joint vmag = vrel(1)*An_End(1) + vrel(2)*An_End(2) + vrel(3)*An_End(3) @@ -6499,22 +6039,22 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat N = p%Members(im)%NElements mem = p%Members(im) - call YawMember(mem, u%PtfmRefY, ErrStat2, ErrMsg2) - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) DO I = mem%i_floor+1, N+1 pos = m%DispNodePosHdn(:, mem%NodeIndx(I)) CALL WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, pos, .FALSE., .TRUE., nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + if (Failed()) return IF (nodeInWater .EQ. 1_IntKi) THEN + call RotateMemberNode(p, im, mem, u%PtfmRefY, u%Mesh%Orientation(:,:,mem%NodeIndx(i)), ErrStat2, ErrMsg2); if (Failed()) return + ! Note: We force each face center to also be wetted if the center node is wetted. Otherwise, the load-smoothing procedure might not work ! Side B - +x_hat side posFC = pos + mem%x_hat * 0.5 * mem%SaMG(i) SVFC = u%Mesh%TranslationVel(:,mem%NodeIndx(i)) + cross_product( u%Mesh%RotationVel(:,mem%NodeIndx(i)), mem%x_hat * 0.5 * mem%SaMG(i) ) - call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); if (Failed()) return vrelFC = dot_product( REAL(FVTmp,ReKi) - SVFC, mem%x_hat ) vrelFCf = mem%VRelNFiltConstB * ( vrelFC + xd%MV_rel_n_FiltStat(1,mem%NodeIndx(I)) ) xd%MV_rel_n_FiltStat(1,mem%NodeIndx(I)) = vrelFCf - vrelFC @@ -6522,7 +6062,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat ! Side B - -x_hat side posFC = pos - mem%x_hat * 0.5 * mem%SaMG(i) SVFC = u%Mesh%TranslationVel(:,mem%NodeIndx(i)) + cross_product( u%Mesh%RotationVel(:,mem%NodeIndx(i)), -mem%x_hat * 0.5 * mem%SaMG(i) ) - call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); if (Failed()) return vrelFC = dot_product( REAL(FVTmp,ReKi) - SVFC, -mem%x_hat ) vrelFCf = mem%VRelNFiltConstB * ( vrelFC + xd%MV_rel_n_FiltStat(2,mem%NodeIndx(I)) ) xd%MV_rel_n_FiltStat(2,mem%NodeIndx(I)) = vrelFCf - vrelFC @@ -6530,7 +6070,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat ! Side A - +y_hat side posFC = pos + mem%y_hat * 0.5 * mem%SbMG(i) SVFC = u%Mesh%TranslationVel(:,mem%NodeIndx(i)) + cross_product( u%Mesh%RotationVel(:,mem%NodeIndx(i)), mem%y_hat * 0.5 * mem%SbMG(i) ) - call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); if (Failed()) return vrelFC = dot_product( REAL(FVTmp,ReKi) - SVFC, mem%y_hat ) vrelFCf = mem%VRelNFiltConstA * ( vrelFC + xd%MV_rel_n_FiltStat(3,mem%NodeIndx(I)) ) xd%MV_rel_n_FiltStat(3,mem%NodeIndx(I)) = vrelFCf - vrelFC @@ -6538,7 +6078,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat ! Side A - -y_hat side posFC = pos - mem%y_hat * 0.5 * mem%SbMG(i) SVFC = u%Mesh%TranslationVel(:,mem%NodeIndx(i)) + cross_product( u%Mesh%RotationVel(:,mem%NodeIndx(i)), -mem%y_hat * 0.5 * mem%SbMG(i) ) - call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); if (Failed()) return vrelFC = dot_product( REAL(FVTmp,ReKi) - SVFC, -mem%y_hat ) vrelFCf = mem%VRelNFiltConstA * ( vrelFC + xd%MV_rel_n_FiltStat(4,mem%NodeIndx(I)) ) xd%MV_rel_n_FiltStat(4,mem%NodeIndx(I)) = vrelFCf - vrelFC @@ -6564,6 +6104,12 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat END IF ! If rectangular member END DO ! Iterate through members +contains + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + end function Failed + END SUBROUTINE Morison_UpdateDiscState !---------------------------------------------------------------------------------------------------------------------------------- END MODULE Morison diff --git a/modules/hydrodyn/src/Morison.txt b/modules/hydrodyn/src/Morison.txt index f8da63bbea..be88a89e57 100644 --- a/modules/hydrodyn/src/Morison.txt +++ b/modules/hydrodyn/src/Morison.txt @@ -130,11 +130,11 @@ typedef ^ Morison_MemberType INTEGER typedef ^ ^ INTEGER MemberID - - - "User-supplied integer ID for this member" - typedef ^ ^ INTEGER NElements - - - "number of elements in this member" - typedef ^ ^ ReKi RefLength - - - "the reference total length for this member" m -typedef ^ ^ ReKi cosPhi_ref - - - "the reference cosine of the inclination angle of the member" - typedef ^ ^ ReKi dl - - - "the reference element length for this member (may be less than MDivSize to achieve uniform element lengths)" m typedef ^ ^ ReKi k {3} - - "unit vector of the member's orientation (may be changed to per-element once additional flexibility is accounted for in HydroDyn)" m typedef ^ ^ ReKi kkt {3}{3} - - "matrix of matmul(k_hat, transpose(k_hat)" - typedef ^ ^ ReKi Ak {3}{3} - - "matrix of I - kkt" - +typedef ^ ^ ReKi CMatrix {3}{3} - - "Rotation matrix from the section local system to the global system" - typedef ^ ^ ReKi x_hat {3} - - "unit vector of rectangular member local x-axis aligned with Side A" - typedef ^ ^ ReKi y_hat {3} - - "unit vector of rectangular member local y-axis aligned with Side B" - typedef ^ ^ ReKi R {:} - - "outer member radius at each node" m @@ -339,7 +339,6 @@ 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 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" - @@ -476,7 +475,6 @@ typedef ^ ^ GridInterp_ typedef ^ ParameterType DbKi DT - - - "Time step for continuous state integration & discrete state update" (sec) 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" - @@ -511,6 +509,7 @@ typedef ^ ^ Morison_Fil # typedef ^ InputType MeshType Mesh - - - "Kinematics of each node input mesh" - typedef ^ ^ ReKi PtfmRefY - - - "Reference platform yaw offset" (rad) +typedef ^ ^ ReKi PtfmRefXY {2} - - "x,y drift of the HydroDyn origin (PRP), consistent with potential-flow body ExctnDisp. Used to displace the Morison members when WaveDisp=0" (m) typedef ^ ^ ReKi PRP {3} - - "Coordinates of the principal reference point" (m) # # diff --git a/modules/hydrodyn/src/Morison_Types.f90 b/modules/hydrodyn/src/Morison_Types.f90 index 3ab02c2918..ba31b86fb5 100644 --- a/modules/hydrodyn/src/Morison_Types.f90 +++ b/modules/hydrodyn/src/Morison_Types.f90 @@ -183,11 +183,11 @@ MODULE Morison_Types INTEGER(IntKi) :: MemberID = 0_IntKi !< User-supplied integer ID for this member [-] INTEGER(IntKi) :: NElements = 0_IntKi !< number of elements in this member [-] REAL(ReKi) :: RefLength = 0.0_ReKi !< the reference total length for this member [m] - REAL(ReKi) :: cosPhi_ref = 0.0_ReKi !< the reference cosine of the inclination angle of the member [-] REAL(ReKi) :: dl = 0.0_ReKi !< the reference element length for this member (may be less than MDivSize to achieve uniform element lengths) [m] REAL(ReKi) , DIMENSION(1:3) :: k = 0.0_ReKi !< unit vector of the member's orientation (may be changed to per-element once additional flexibility is accounted for in HydroDyn) [m] REAL(ReKi) , DIMENSION(1:3,1:3) :: kkt = 0.0_ReKi !< matrix of matmul(k_hat, transpose(k_hat) [-] REAL(ReKi) , DIMENSION(1:3,1:3) :: Ak = 0.0_ReKi !< matrix of I - kkt [-] + REAL(ReKi) , DIMENSION(1:3,1:3) :: CMatrix = 0.0_ReKi !< Rotation matrix from the section local system to the global system [-] REAL(ReKi) , DIMENSION(1:3) :: x_hat = 0.0_ReKi !< unit vector of rectangular member local x-axis aligned with Side A [-] REAL(ReKi) , DIMENSION(1:3) :: y_hat = 0.0_ReKi !< unit vector of rectangular member local y-axis aligned with Side B [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: R !< outer member radius at each node [m] @@ -407,7 +407,6 @@ 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 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 [-] @@ -543,7 +542,6 @@ MODULE Morison_Types REAL(DbKi) :: DT = 0.0_R8Ki !< Time step for continuous state integration & discrete state update [(sec)] 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 [-] @@ -577,6 +575,7 @@ MODULE Morison_Types TYPE, PUBLIC :: Morison_InputType TYPE(MeshType) :: Mesh !< Kinematics of each node input mesh [-] REAL(ReKi) :: PtfmRefY = 0.0_ReKi !< Reference platform yaw offset [(rad)] + REAL(ReKi) , DIMENSION(1:2) :: PtfmRefXY = 0.0_ReKi !< x,y drift of the HydroDyn origin (PRP), consistent with potential-flow body ExctnDisp. Used to displace the Morison members when WaveDisp=0 [(m)] REAL(ReKi) , DIMENSION(1:3) :: PRP = 0.0_ReKi !< Coordinates of the principal reference point [(m)] END TYPE Morison_InputType ! ======================= @@ -590,10 +589,11 @@ MODULE Morison_Types integer(IntKi), public, parameter :: Morison_x_DummyContState = 1 ! Morison%DummyContState integer(IntKi), public, parameter :: Morison_u_Mesh = 2 ! Morison%Mesh integer(IntKi), public, parameter :: Morison_u_PtfmRefY = 3 ! Morison%PtfmRefY - integer(IntKi), public, parameter :: Morison_u_PRP = 4 ! Morison%PRP - integer(IntKi), public, parameter :: Morison_y_Mesh = 5 ! Morison%Mesh - integer(IntKi), public, parameter :: Morison_y_VisMesh = 6 ! Morison%VisMesh - integer(IntKi), public, parameter :: Morison_y_WriteOutput = 7 ! Morison%WriteOutput + integer(IntKi), public, parameter :: Morison_u_PtfmRefXY = 4 ! Morison%PtfmRefXY + integer(IntKi), public, parameter :: Morison_u_PRP = 5 ! Morison%PRP + integer(IntKi), public, parameter :: Morison_y_Mesh = 6 ! Morison%Mesh + integer(IntKi), public, parameter :: Morison_y_VisMesh = 7 ! Morison%VisMesh + integer(IntKi), public, parameter :: Morison_y_WriteOutput = 8 ! Morison%WriteOutput contains @@ -1291,11 +1291,11 @@ subroutine Morison_CopyMemberType(SrcMemberTypeData, DstMemberTypeData, CtrlCode DstMemberTypeData%MemberID = SrcMemberTypeData%MemberID DstMemberTypeData%NElements = SrcMemberTypeData%NElements DstMemberTypeData%RefLength = SrcMemberTypeData%RefLength - DstMemberTypeData%cosPhi_ref = SrcMemberTypeData%cosPhi_ref DstMemberTypeData%dl = SrcMemberTypeData%dl DstMemberTypeData%k = SrcMemberTypeData%k DstMemberTypeData%kkt = SrcMemberTypeData%kkt DstMemberTypeData%Ak = SrcMemberTypeData%Ak + DstMemberTypeData%CMatrix = SrcMemberTypeData%CMatrix DstMemberTypeData%x_hat = SrcMemberTypeData%x_hat DstMemberTypeData%y_hat = SrcMemberTypeData%y_hat if (allocated(SrcMemberTypeData%R)) then @@ -2281,11 +2281,11 @@ subroutine Morison_PackMemberType(RF, Indata) call RegPack(RF, InData%MemberID) call RegPack(RF, InData%NElements) call RegPack(RF, InData%RefLength) - call RegPack(RF, InData%cosPhi_ref) call RegPack(RF, InData%dl) call RegPack(RF, InData%k) call RegPack(RF, InData%kkt) call RegPack(RF, InData%Ak) + call RegPack(RF, InData%CMatrix) call RegPack(RF, InData%x_hat) call RegPack(RF, InData%y_hat) call RegPackAlloc(RF, InData%R) @@ -2395,11 +2395,11 @@ subroutine Morison_UnPackMemberType(RF, OutData) call RegUnpack(RF, OutData%MemberID); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%NElements); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%RefLength); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%cosPhi_ref); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%dl); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%k); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%kkt); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Ak); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%CMatrix); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%x_hat); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%y_hat); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%R); if (RegCheckErr(RF, RoutineName)) return @@ -3268,7 +3268,6 @@ subroutine Morison_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode, E ErrMsg = '' DstInitInputData%Gravity = SrcInitInputData%Gravity DstInitInputData%WaveDisp = SrcInitInputData%WaveDisp - DstInitInputData%AMMod = SrcInitInputData%AMMod DstInitInputData%HstMod = SrcInitInputData%HstMod DstInitInputData%NJoints = SrcInitInputData%NJoints DstInitInputData%NNodes = SrcInitInputData%NNodes @@ -3717,7 +3716,6 @@ subroutine Morison_PackInitInput(RF, Indata) if (RF%ErrStat >= AbortErrLev) return 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) @@ -3924,7 +3922,6 @@ subroutine Morison_UnPackInitInput(RF, OutData) if (RF%ErrStat /= ErrID_None) return 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 @@ -4946,7 +4943,6 @@ subroutine Morison_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrM DstParamData%DT = SrcParamData%DT 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 @@ -5261,7 +5257,6 @@ subroutine Morison_PackParam(RF, Indata) call RegPack(RF, InData%DT) 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)) @@ -5352,7 +5347,6 @@ subroutine Morison_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%DT); if (RegCheckErr(RF, RoutineName)) return 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) @@ -5474,6 +5468,7 @@ subroutine Morison_CopyInput(SrcInputData, DstInputData, CtrlCode, ErrStat, ErrM call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return DstInputData%PtfmRefY = SrcInputData%PtfmRefY + DstInputData%PtfmRefXY = SrcInputData%PtfmRefXY DstInputData%PRP = SrcInputData%PRP end subroutine @@ -5497,6 +5492,7 @@ subroutine Morison_PackInput(RF, Indata) if (RF%ErrStat >= AbortErrLev) return call MeshPack(RF, InData%Mesh) call RegPack(RF, InData%PtfmRefY) + call RegPack(RF, InData%PtfmRefXY) call RegPack(RF, InData%PRP) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -5508,6 +5504,7 @@ subroutine Morison_UnPackInput(RF, OutData) if (RF%ErrStat /= ErrID_None) return call MeshUnpack(RF, OutData%Mesh) ! Mesh call RegUnpack(RF, OutData%PtfmRefY); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%PtfmRefXY); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%PRP); if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -5685,6 +5682,7 @@ SUBROUTINE Morison_Input_ExtrapInterp1(u1, u2, tin, u_out, tin_out, ErrStat, Err CALL MeshExtrapInterp1(u1%Mesh, u2%Mesh, tin, u_out%Mesh, tin_out, ErrStat2, ErrMsg2) CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) u_out%PtfmRefY = a1*u1%PtfmRefY + a2*u2%PtfmRefY + u_out%PtfmRefXY = a1*u1%PtfmRefXY + a2*u2%PtfmRefXY u_out%PRP = a1*u1%PRP + a2*u2%PRP END SUBROUTINE @@ -5746,6 +5744,7 @@ SUBROUTINE Morison_Input_ExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, ErrStat, CALL MeshExtrapInterp2(u1%Mesh, u2%Mesh, u3%Mesh, tin, u_out%Mesh, tin_out, ErrStat2, ErrMsg2) CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) u_out%PtfmRefY = a1*u1%PtfmRefY + a2*u2%PtfmRefY + a3*u3%PtfmRefY + u_out%PtfmRefXY = a1*u1%PtfmRefXY + a2*u2%PtfmRefXY + a3*u3%PtfmRefXY u_out%PRP = a1*u1%PRP + a2*u2%PRP + a3*u3%PRP END SUBROUTINE @@ -6044,6 +6043,8 @@ subroutine Morison_VarPackInput(V, u, ValAry) call MV_PackMesh(V, u%Mesh, ValAry) ! Mesh case (Morison_u_PtfmRefY) VarVals(1) = u%PtfmRefY ! Scalar + case (Morison_u_PtfmRefXY) + VarVals = u%PtfmRefXY(V%iLB:V%iUB) ! Rank 1 Array case (Morison_u_PRP) VarVals = u%PRP(V%iLB:V%iUB) ! Rank 1 Array case default @@ -6072,6 +6073,8 @@ subroutine Morison_VarUnpackInput(V, ValAry, u) call MV_UnpackMesh(V, ValAry, u%Mesh) ! Mesh case (Morison_u_PtfmRefY) u%PtfmRefY = VarVals(1) ! Scalar + case (Morison_u_PtfmRefXY) + u%PtfmRefXY(V%iLB:V%iUB) = VarVals ! Rank 1 Array case (Morison_u_PRP) u%PRP(V%iLB:V%iUB) = VarVals ! Rank 1 Array end select @@ -6086,6 +6089,8 @@ function Morison_InputFieldName(DL) result(Name) Name = "u%Mesh" case (Morison_u_PtfmRefY) Name = "u%PtfmRefY" + case (Morison_u_PtfmRefXY) + Name = "u%PtfmRefXY" case (Morison_u_PRP) Name = "u%PRP" case default diff --git a/openfast_io/openfast_io/FAST_reader.py b/openfast_io/openfast_io/FAST_reader.py index 63e8c1dc32..23b629f087 100644 --- a/openfast_io/openfast_io/FAST_reader.py +++ b/openfast_io/openfast_io/FAST_reader.py @@ -2031,7 +2031,6 @@ def read_HydroDyn(self, hd_file): #STRIP THEORY OPTIONS 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 diff --git a/openfast_io/openfast_io/FAST_writer.py b/openfast_io/openfast_io/FAST_writer.py index 477593306e..012afb23ee 100644 --- a/openfast_io/openfast_io/FAST_writer.py +++ b/openfast_io/openfast_io/FAST_writer.py @@ -1719,8 +1719,7 @@ def write_HydroDyn(self): f.write(ln) 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']['WaveDisp'], 'WaveDisp', '- Method of displacing strip-theory members and joints {0: use potential-flow consistent displacement, 1: use exact instantaneous displacement} (switch)\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') diff --git a/reg_tests/r-test b/reg_tests/r-test index e773db8bbf..deded14c9c 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit e773db8bbf04208641791de8171a4607d1997a2a +Subproject commit deded14c9c4133492d4502444b09a0b82ac26d6b