diff --git a/docs/source/user/hydrodyn/input_files.rst b/docs/source/user/hydrodyn/input_files.rst index b93aaa8035..4dde04b2c1 100644 --- a/docs/source/user/hydrodyn/input_files.rst +++ b/docs/source/user/hydrodyn/input_files.rst @@ -715,6 +715,27 @@ Members Each member in your model will have hydrodynamic coefficients, which are specified using one of the three models (**MCoefMod**). Model 1 uses a single set of coefficients found in the SIMPLE HYDRODYNAMIC COEFFICIENTS sections. Model 2 is depth-based, and is determined via the table found in the DEPTH-BASED HYDRODYNAMIC COEFFICIENTS sections. Model 3 specifies coefficients for a particular member, by referring to the MEMBER-BASED HYDRODYNAMIC COEFFICIENTS sections. Depending on **MSecGeom**, HydroDyn will either use the hydrodynamic coefficients from the input tables for circular sections or rectangular sections as appropriate. The **MHstLMod** switch controls the computation of hydrostatic loads on strip-theory members when **PropPot** = FALSE. Setting **MHstLMod** to 0 disables hydrostatic load. If set to 1, hydrostatic loads will be computed analytically. This approach is efficient, but it only works for fully submerged or surface-piercing members that are far from horizontal without partially wetted endplates. For nearly horizontal members close to the free surface or members that experience partially wetted endplates, a semi-numerical approach for hydrostatic load can be selected by setting **MHstLMod** to 2. This approach works with any member positioning in relation to the free surface at the cost of slightly longer computing time. Note that for members with rectangular cross sections, **MHstLMod** must be either 0 or 2. The analytical approach with **MHstLMod** set to 1 is not available for rectangular members. The **PropPot** flag indicates whether the corresponding member is included in the potential-flow solution. When **PropPot** = TRUE, only viscous-drag loads and ballasting loads will be computed for that member, with the assumption that all other load components are already included in the potential-flow solution. +In addition to the required inputs for each member explained above, we have several optional inputs for more advanced transverse drag modeling for rectangular members. These are **FDMod**, **VnCOffA**, **VnCOffB**, **FDLoFScA**, and **FDLoFScB**, as shown in the example below. For each rectangular member (**MSecGeom** = 2), users can either omit all these optional inputs together, in which case the simple default drag formulation will be used, or provide all five optional inputs to enable more advanced drag formulation. It is not allowed to provide only a subset of these optional inputs, while omitting others. These inputs are not relevant to members with circular cross sections (**MSecGeom** = 1) and should be omitted for these members. HydroDyn will simply ignore these optional inputs if provided for a circular member, and a warning will be displayed at the start of the simulation. + +The first optional input **FDMod** can be 0, 1, or 2. Setting **FDMod** to 0 leads to the default transverse drag formulation for rectangular members. The two components of the distributed transverse drag force normal to Side A and normal to Side B are both computed using the flow and structure velocities along the axis/centerline of the member based on **CdA** and **CdB** defined previously. This is referred to as a centerline-based formulation. Consistent with the transverse drag along cylindrical members, a lead coefficient of 1/2 is used in the expression for the drag force in both directions. + +When **FDMod** is set to 1, a face-based drag formulation is used. In this case, four side-normal drag components are computed using the flow and structure velocities at the midpoints of the four sides of the rectangular section instead of at the center. On each side, only the normal components of the flow and structure velocities are used to compute the drag force. The components parallel to the side are ignored. This is analogous to the modeling of drag force on member endplates at the joints with **AxFDMod** set to 0, if we interpret each of the four side walls of a rectangular member as an array of endplates placed along the length of the member. Also analogous to endplate drag force, a lead coefficient of 1/4 instead of 1/2 is used when computing the drag force on each side. This is to account for the fact that we now have two drag components on the pair of sides normal to Side A and two drag components on the other pair of sides normal to Side B. This change in the lead coefficient ensures a degree of consistency with the default model of **FDMod** = 0. + +Setting **FDMod** to 2 also leads to a face-based drag formulation being used. However, drag force will only be applied on the sides where the relative flow is directed away from the structure, i.e., the suction side, where flow separation is expected, not on the sides the relative flow is impinging on where flow separation is unlikely. This is again analogous to the modeling of endplate drag force at the joints with **AxFDMod** set to 1. This suction-side-only formulation is expected to work better for hybrid members modeled by potential-flow theory. Note that with **FDMod** set to 2, a lead coefficient of 1/2 instead of 1/4 is again used for each side. This is because for each pair of parallel sides, typically only one side, the suction side, will have a non-zero drag force, not both. This change in the lead coefficient is again included to ensure a degree of consistency with the other two drag formulations. + +With **FDMod** set to either 1 or 2, users can optional enable high-pass filtering of the relative flow velocity when computing the transverse drag force. The cutoff frequency of the high-pass filter can be set through **VnCOffA** for the drag components normal to Side A and **VnCOffB** for the drag components normal to Side B. These inputs are analogous to **AxVnCOff** for the endplate drag force at the joints. The two weighting factors **FDLoFScA** for the drag components normal to Side A and **FDLoFScB** for the drag components normal to Side B control the proportion of drag force computed using the unfiltered relative flow velocity. Both inputs should be between 0 and 1. When set to 1, the drag force is computed purely based on the unfiltered relative flow velocity, equivalent to no velocity filtering. When set to 0, the drag force is computed purely based on the high-pass filtered relative flow velocity. Again, these inputs are analogous to **AxFDLoFSc** for the member endplate drag forces at the joints. The velocity filtering allows the user to attenuate the drag force in response to low-frequency motion. When applied to the drag force on rectangular pontoons, this option can help improve the predictions of low-frequency heave and pitch resonance motions of a floater if properly tuned. To disable velocity filtering, users can either set **FDLoFScA** and **FDLoFScB** to 1 or set **VnCOffA** and **VnCOffB** to a number equal to or less than zero. Velocity filtering is not supported with the centerline-based approach, and the relevant inputs will be ignored if **FDMod** = 0. + +If the optional inputs are omitted for a rectangular member, the default centerline-based drag formulation with **FDMod** = 0 will be used, and velocity filtering will not be applied. Irrespective of the choice of **FDMod** and velocity filtering, the member drag coefficients can be set using any of the three **MCoefMod** options without any additional restrictions. Finally, note that these optional inputs only affect the computation of distributed transverse drag forces perpendicular to the member axis. They do not affect the distributed axial drag force along the member axis. + +.. code:: + + -------------------- MEMBERS ------------------------------------------------- + 2 NMembers - Number of members (-) + MemberID MJointID1 MJointID2 MPropSetID1 MPropSetID2 MSecGeom MSpinOrient MDivSize MCoefMod MHstLMod PropPot FDMod VnCOffA VnCOffB FDLoFScA FDLoFScB [MCoefMod=1: use simple coeff table, 2: use depth-based coeff table, 3: use member-based coeff table] [PropPot/=0 if member is modeled with potential-flow theory] + (-) (-) (-) (-) (-) (switch) (deg) (m) (switch) (switch) (flag) (switch) (Hz) (Hz) (-) (-) + 1 1 2 1 1 1 0 0.5 1 2 FALSE ! Circular members should not have any optional inputs. If provided, these inputs will be ignored. + 2 3 4 1 2 2 60 0.5 2 2 FALSE 2 0.05 0.05 0.5 0.5 ! A rectangular member with optional inputs for drag. + .. TODO 7.5.2 is the theory section which does not yet exist. .. Section 7.5.2 discusses the difference between the user-supplied discretization and the simulation discretization. diff --git a/modules/hydrodyn/src/HydroDyn_Input.f90 b/modules/hydrodyn/src/HydroDyn_Input.f90 index d9960ed40b..8cf72660e5 100644 --- a/modules/hydrodyn/src/HydroDyn_Input.f90 +++ b/modules/hydrodyn/src/HydroDyn_Input.f90 @@ -915,18 +915,51 @@ SUBROUTINE HydroDyn_ParseInput( InputFileName, OutRootName, FileInfo_In, InputFi DO I = 1,InputFileData%Morison%NMembers ! We can't use the ParseAry here since PropPot is a logical Line = FileInfo_In%Lines(CurLine) + READ(Line,*,IOSTAT=ErrStat2) InputFileData%Morison%InpMembers(I)%MemberID, InputFileData%Morison%InpMembers(I)%MJointID1, & InputFileData%Morison%InpMembers(I)%MJointID2, InputFileData%Morison%InpMembers(I)%MPropSetID1, & InputFileData%Morison%InpMembers(I)%MPropSetID2, InputFileData%Morison%InpMembers(I)%MSecGeom, & InputFileData%Morison%InpMembers(I)%MSpinOrient, InputFileData%Morison%InpMembers(I)%MDivSize, & InputFileData%Morison%InpMembers(I)%MCoefMod, InputFileData%Morison%InpMembers(I)%MHstLMod, & - InputFileData%Morison%InpMembers(I)%PropPot + InputFileData%Morison%InpMembers(I)%PropPot, InputFileData%Morison%InpMembers(I)%FDMod, & + InputFileData%Morison%InpMembers(I)%VnCOffA, InputFileData%Morison%InpMembers(I)%VnCOffB, & + InputFileData%Morison%InpMembers(I)%FDLoFScA, InputFileData%Morison%InpMembers(I)%FDLoFScB IF ( ErrStat2 /= 0 ) THEN - ErrStat2 = ErrID_Fatal - ErrMsg2 = 'Error reading members table row '//trim( Int2LStr(I))//', line ' & - //trim( Int2LStr(FileInfo_In%FileLine(CurLine)))//' of file '//trim(FileInfo_In%FileList(FileInfo_In%FileIndx(CurLine))) - if (Failed()) return; + READ(Line,*,IOSTAT=ErrStat2) InputFileData%Morison%InpMembers(I)%MemberID, InputFileData%Morison%InpMembers(I)%MJointID1, & + InputFileData%Morison%InpMembers(I)%MJointID2, InputFileData%Morison%InpMembers(I)%MPropSetID1, & + InputFileData%Morison%InpMembers(I)%MPropSetID2, InputFileData%Morison%InpMembers(I)%MSecGeom, & + InputFileData%Morison%InpMembers(I)%MSpinOrient, InputFileData%Morison%InpMembers(I)%MDivSize, & + InputFileData%Morison%InpMembers(I)%MCoefMod, InputFileData%Morison%InpMembers(I)%MHstLMod, & + InputFileData%Morison%InpMembers(I)%PropPot + IF ( ErrStat2 /= 0 ) THEN + ErrStat2 = ErrID_Fatal + ErrMsg2 = 'Error reading members table row '//trim( Int2LStr(I))//', line ' & + //trim( Int2LStr(FileInfo_In%FileLine(CurLine)))//' of file '//trim(FileInfo_In%FileList(FileInfo_In%FileIndx(CurLine))) + if (Failed()) return; + ELSE + InputFileData%Morison%InpMembers(I)%FDMod = 0_IntKi + InputFileData%Morison%InpMembers(I)%VnCOffA = -1.0_ReKi + InputFileData%Morison%InpMembers(I)%VnCOffB = -1.0_ReKi + InputFileData%Morison%InpMembers(I)%FDLoFScA = 1.0_ReKi + InputFileData%Morison%InpMembers(I)%FDLoFScB = 1.0_ReKi + END IF + ELSE + IF ( InputFileData%Morison%InpMembers(I)%MSecGeom /= MSecGeom_Rec ) THEN + call WrScr('HydroDyn Warning: The optional member inputs FDMod, VnCOffA, VnCOffB, FDLoFScA, and FDLoFScB are only applicable to members with rectangular sections. These will be ignored for Member ID '//TRIM(num2Lstr(InputFileData%Morison%InpMembers(I)%MemberID))//'. ') + InputFileData%Morison%InpMembers(I)%FDMod = 0_IntKi + InputFileData%Morison%InpMembers(I)%VnCOffA = -1.0_ReKi + InputFileData%Morison%InpMembers(I)%VnCOffB = -1.0_ReKi + InputFileData%Morison%InpMembers(I)%FDLoFScA = 1.0_ReKi + InputFileData%Morison%InpMembers(I)%FDLoFScB = 1.0_ReKi + ELSE IF ( InputFileData%Morison%InpMembers(I)%FDMod == 0_IntKi ) THEN + call WrScr('HydroDyn Warning: Velocity filtering for rectangular-member transverse drag force is only available with FDMod = 1 or 2. The optional member inputs VnCOffA, VnCOffB, FDLoFScA, and FDLoFScB will be ignored for Member ID '//TRIM(num2Lstr(InputFileData%Morison%InpMembers(I)%MemberID))//'. ') + InputFileData%Morison%InpMembers(I)%VnCOffA = -1.0_ReKi + InputFileData%Morison%InpMembers(I)%VnCOffB = -1.0_ReKi + InputFileData%Morison%InpMembers(I)%FDLoFScA = 1.0_ReKi + InputFileData%Morison%InpMembers(I)%FDLoFScB = 1.0_ReKi + END IF END IF + InputFileData%Morison%InpMembers(I)%MSpinOrient = InputFileData%Morison%InpMembers(I)%MSpinOrient * D2R if ( InputFileData%Echo ) WRITE(UnEc, '(A)') trim(FileInfo_In%Lines(CurLine)) ! Echo this line @@ -2482,8 +2515,8 @@ SUBROUTINE HydroDynInput_ProcessInitData( InitInp, Interval, InputFileData, ErrS RETURN END IF - IF ( InputFileData%Morison%InpMembers(I)%MSecGeom == MSecGeom_Rec .AND. InputFileData%Morison%InpMembers(I)%MHstLMod /= 2 ) THEN - CALL SetErrStat( ErrID_Fatal,'MHstLMod must be 2 for rectangular members.',ErrStat,ErrMsg,RoutineName) + IF ( InputFileData%Morison%InpMembers(I)%MSecGeom == MSecGeom_Rec .AND. InputFileData%Morison%InpMembers(I)%MHstLMod /= 0 .AND. InputFileData%Morison%InpMembers(I)%MHstLMod /= 2 ) THEN + CALL SetErrStat( ErrID_Fatal,'MHstLMod must be 0 or 2 for rectangular members.',ErrStat,ErrMsg,RoutineName) RETURN END IF @@ -2492,6 +2525,22 @@ SUBROUTINE HydroDynInput_ProcessInitData( InitInp, Interval, InputFileData, ErrS RETURN END IF + ! Optional member inputs for rectangular members + IF ( InputFileData%Morison%InpMembers(I)%FDMod /= 0_IntKi .AND. InputFileData%Morison%InpMembers(I)%FDMod /= 1_IntKi .AND. InputFileData%Morison%InpMembers(I)%FDMod /= 2_IntKi ) THEN + CALL SetErrStat( ErrID_Fatal,'FDMod must be 0 (centerline-based drag), 1 (face-based drag), or 2 (face-based suction-side-only drag) for rectangular members.',ErrStat,ErrMsg,RoutineName) + RETURN + END IF + + IF ( InputFileData%Morison%InpMembers(I)%FDLoFScA < 0.0_ReKi .OR. InputFileData%Morison%InpMembers(I)%FDLoFScA > 1.0_ReKi ) THEN + CALL SetErrStat( ErrID_Fatal,'FDLoFScA and FDLoFScB for rectangular members must be between 0 and 1 inclusive.',ErrStat,ErrMsg,RoutineName) + RETURN + END IF + + IF ( InputFileData%Morison%InpMembers(I)%FDLoFScB < 0.0_ReKi .OR. InputFileData%Morison%InpMembers(I)%FDLoFScB > 1.0_ReKi ) THEN + CALL SetErrStat( ErrID_Fatal,'FDLoFScA and FDLoFScB for rectangular members must be between 0 and 1 inclusive.',ErrStat,ErrMsg,RoutineName) + RETURN + END IF + END DO END IF diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index a03095fae0..fee2da24bd 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -153,6 +153,80 @@ SUBROUTINE Morison_DirCosMtrx_noSpin( pos0, pos1, DirCos ) END SUBROUTINE Morison_DirCosMtrx_noSpin + +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 + REAL(ReKi), INTENT( OUT) :: pos(:,:) ! Displaced node positions + + REAL(ReKi) :: Orient(3,3) + INTEGER(IntKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + + ! Undisplaced node position + 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 + 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) + call GetPtfmRefYOrient(u%PtfmRefY, Orient, ErrStat2, ErrMsg2) + pos = matmul(transpose(Orient),pos) + END IF + +END SUBROUTINE GetDisplacedNodePosition + + +SUBROUTINE YawMember(member, PtfmRefY, ErrStat, ErrMsg) + Type(Morison_MemberType), intent(inout) :: member + Real(ReKi), intent(in ) :: PtfmRefY + 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 + + Character(*), parameter :: RoutineName = 'YawMember' + + 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 + + call hiFrameTransform(h2i,PtfmRefY,member%y_hat,y_hat,ErrStat2,ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + member%y_hat = y_hat + + END IF + +END SUBROUTINE YawMember + !==================================================================================================== SUBROUTINE GetDistance ( a, b, l ) ! This private subroutine computes the distance between points a and b. @@ -2635,7 +2709,22 @@ subroutine SetupMembers( InitInp, p, m, errStat, errMsg ) p%Members(i)%MHstLMod = InitInp%InpMembers(i)%MHstLMod p%Members(i)%MSecGeom = InitInp%InpMembers(i)%MSecGeom p%Members(i)%MSpinOrient = InitInp%InpMembers(i)%MSpinOrient - + p%Members(i)%FDMod = InitInp%InpMembers(i)%FDMod + IF (InitInp%InpMembers(i)%VnCOffA .LE. 0.0_ReKi) THEN + p%Members(i)%VRelNFiltConstA = 1.0_ReKi + p%Members(i)%DragLoFScA = 1.0_ReKi + ELSE + p%Members(i)%VRelNFiltConstA = exp(-2.0*Pi*InitInp%InpMembers(i)%VnCOffA * p%DT) + p%Members(i)%DragLoFScA = InitInp%InpMembers(i)%FDLoFScA + END IF + IF (InitInp%InpMembers(i)%VnCOffB .LE. 0.0_ReKi) THEN + p%Members(i)%VRelNFiltConstB = 1.0_ReKi + p%Members(i)%DragLoFScB = 1.0_ReKi + ELSE + p%Members(i)%VRelNFiltConstB = exp(-2.0*Pi*InitInp%InpMembers(i)%VnCOffB * p%DT) + p%Members(i)%DragLoFScB = InitInp%InpMembers(i)%FDLoFScB + END IF + call AllocateMemberDataArrays(p%Members(i), m%MemberLoads(i), errStat2, errMsg2) call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'SetupMembers') if (ErrStat >= AbortErrLev) return @@ -2865,7 +2954,7 @@ SUBROUTINE Morison_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, In ! Define initial system states here: x%DummyContState = 0 - !xd%DummyDiscState = 0 + ALLOCATE ( xd%V_rel_n_FiltStat(p%NJoints), STAT = ErrStat ) IF ( ErrStat /= ErrID_None ) THEN ErrMsg = ' Error allocating space for V_rel_n_FiltStat array.' @@ -2874,6 +2963,14 @@ SUBROUTINE Morison_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, In END IF xd%V_rel_n_FiltStat = 0.0_ReKi + ALLOCATE ( xd%MV_rel_n_FiltStat(4,p%NNodes), STAT = ErrStat ) + IF ( ErrStat /= ErrID_None ) THEN + ErrMsg = ' Error allocating space for MV_rel_n_FiltStat array.' + ErrStat = ErrID_Fatal + RETURN + END IF + xd%MV_rel_n_FiltStat = 0.0_ReKi + z%DummyConstrState = 0 OtherState%DummyOtherState = 0 @@ -3433,6 +3530,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, REAL(ReKi) :: FDynPFSInt REAL(ReKi) :: vrelFSInt(3) REAL(ReKi) :: FAMCFFSInt(3) + REAL(ReKi) :: CdFSInt, CdAFSInt, CdBFSInt, AxCdFSInt, CaFSInt, CaAFSInt, CaBFSInt, AxCaFSInt, CpFSInt, AxCpFSInt INTEGER(IntKi) :: MemSubStat, NumFSX REAL(DbKi) :: theta1, theta2 REAL(ReKi) :: x_hat(3), x_hat1(3), x_hat2(3), y_hat(3), y_hat1(3), y_hat2(3), z_hat(3), posMid(3), zetaMid, FSPt(3) @@ -3449,8 +3547,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, !=============================================================================================== ! Get displaced positions of the hydrodynamic nodes - CALL GetDisplacedNodePosition( .FALSE., m%DispNodePosHdn ) ! For hydrodynamic loads; depends on WaveDisp and WaveStMod - CALL GetDisplacedNodePosition( .TRUE. , m%DispNodePosHst ) ! For hydrostatic loads; always use actual displaced position + CALL GetDisplacedNodePosition( u, p, .FALSE., m%DispNodePosHdn ) ! For hydrodynamic loads; depends on WaveDisp and WaveStMod + 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 @@ -3941,10 +4039,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, 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 - f_hydro = 0.5*mem%CdB(i)*p%WaveField%WtrDens*mem%SbMG(i)*TwoNorm(vec)*Dot_Product(vec,mem%x_hat)*mem%x_hat + & ! local x-direction - 0.5*mem%CdA(i)*p%WaveField%WtrDens*mem%SaMG(i)*TwoNorm(vec)*Dot_Product(vec,mem%y_hat)*mem%y_hat + & ! local z-direction - 0.25*mem%AxCd(i)*p%WaveField%WtrDens * (dSadl_p*mem%SbMG(i) + dSbdl_p*mem%SaMG(i)) * & ! axial part - abs(dot_product( mem%k, m%vrel(:,mem%NodeIndx(i)) )) * matmul( mem%kkt, m%vrel(:,mem%NodeIndx(i)) ) ! axial part cont'd + Call GetDistDrag_Rec(Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END IF 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) @@ -4049,11 +4144,17 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, IF (mem%MSecGeom==MSecGeom_Cyl) THEN 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) + RMGFSInt = SubRatio * mem%RMG( FSElem+1) + (1.0-SubRatio) * mem%RMG( FSElem) + CdFSInt = SubRatio * mem%Cd( FSElem+1) + (1.0-SubRatio) * mem%Cd( FSElem) + AxCdFSInt = SubRatio * mem%AxCd(FSElem+1) + (1.0-SubRatio) * mem%AxCd(FSElem) + CaFSInt = SubRatio * mem%Ca( FSElem+1) + (1.0-SubRatio) * mem%Ca( FSElem) + AxCaFSInt = SubRatio * mem%AxCa(FSElem+1) + (1.0-SubRatio) * mem%AxCa(FSElem) + CpFSInt = SubRatio * mem%Cp( FSElem+1) + (1.0-SubRatio) * mem%Cp( FSElem) + AxCpFSInt = SubRatio * mem%AxCp(FSElem+1) + (1.0-SubRatio) * mem%AxCp(FSElem) vec = matmul( mem%Ak,vrelFSInt ) - F_DS = mem%Cd(FSElem)*p%WaveField%WtrDens*RMGFSInt*TwoNorm(vec)*vec + & - 0.5*mem%AxCd(FSElem)*p%WaveField%WtrDens*pi*RMGFSInt*dRdl_p * & + F_DS = CdFSInt*p%WaveField%WtrDens*RMGFSInt*TwoNorm(vec)*vec + & + 0.5*AxCdFSInt*p%WaveField%WtrDens*pi*RMGFSInt*dRdl_p * & abs(dot_product( mem%k, vrelFSInt )) * matmul( mem%kkt, vrelFSInt ) ! Hydrodynamic added mass and inertia loads @@ -4061,22 +4162,20 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! ------------------- hydrodynamic added mass loads: sides: Section 7.1.3 ------------------------ IF (p%AMMod > 0_IntKi) THEN - Am = mem%Ca(FSElem)*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt*mem%Ak + & - 2.0*mem%AxCa(FSElem)*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt*dRdl_p*mem%kkt - F_AS = -matmul( Am, & - SubRatio * u%Mesh%TranslationAcc(:,mem%NodeIndx(FSElem+1)) + & - (1.0-SubRatio) * u%Mesh%TranslationAcc(:,mem%NodeIndx(FSElem )) ) + 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 ! ------------------- hydrodynamic inertia loads: sides: Section 7.1.4 ------------------------ IF ( mem%PropMCF) THEN - F_IS= p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt * matmul( mem%Ak, FAMCFFSInt ) + & - 2.0*mem%AxCa(FSElem)*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt*dRdl_p * matmul( mem%kkt, FAFSInt ) + & - 2.0*mem%AxCp(FSElem) *pi*RMGFSInt *dRdl_pp * FDynPFSInt*mem%k + F_IS= p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt * matmul( mem%Ak, FAMCFFSInt ) + & + 2.0*AxCaFSInt*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt*dRdl_p * matmul( mem%kkt, FAFSInt ) + & + 2.0*AxCpFSInt *pi*RMGFSInt *dRdl_pp * FDynPFSInt*mem%k ELSE - F_IS=(mem%Ca(FSElem)+mem%Cp(FSElem))*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt * matmul( mem%Ak, FAFSInt ) + & - 2.0*mem%AxCa(FSElem)*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt*dRdl_p * matmul( mem%kkt, FAFSInt ) + & - 2.0*mem%AxCp(FSElem) *pi*RMGFSInt *dRdl_pp * FDynPFSInt*mem%k + F_IS=(CaFSInt+CpFSInt)*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt * matmul( mem%Ak, FAFSInt ) + & + 2.0*AxCaFSInt*p%WaveField%WtrDens*pi*RMGFSInt*RMGFSInt*dRdl_p * matmul( mem%kkt, FAFSInt ) + & + 2.0*AxCpFSInt *pi*RMGFSInt *dRdl_pp * FDynPFSInt*mem%k END IF END IF ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN @@ -4086,30 +4185,34 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, 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) + CpFSInt = SubRatio * mem%Cp( FSElem+1) + (1.0-SubRatio) * mem%Cp( FSElem) + AxCpFSInt = SubRatio * mem%AxCp(FSElem+1) + (1.0-SubRatio) * mem%AxCp(FSElem) - vec = matmul( mem%Ak,vrelFSInt ) - F_DS = 0.5*mem%CdB(FSElem)*p%WaveField%WtrDens*SbMGFSInt*TwoNorm(vec)*Dot_Product(vec,mem%x_hat)*mem%x_hat + & ! local x-direction - 0.5*mem%CdA(FSElem)*p%WaveField%WtrDens*SaMGFSInt*TwoNorm(vec)*Dot_Product(vec,mem%y_hat)*mem%y_hat + & ! local z-direction - 0.25*mem%AxCd(FSElem)*p%WaveField%WtrDens * (dSadl_p*SbMGFSInt + dSbdl_p*SaMGFSInt) * & ! axial part - abs(dot_product( mem%k, vrelFSInt )) * matmul( mem%kkt, vrelFSInt ) ! axial part cont'd + Call GetDistDrag_Rec(Time,mem,FSElem,dSadl_p,dSbdl_p,F_DS,ErrStat2,ErrMsg2,SubRatio,vrelFSInt); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ! 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*mem%CaB(FSElem) * 0.25*pi*SbMGFSInt*SbMGFSInt * Dot_Product(SAFSInt,mem%x_hat)*mem%x_hat & - -p%WaveField%WtrDens*mem%CaA(FSElem) * 0.25*pi*SaMGFSInt*SaMGFSInt * Dot_Product(SAFSInt,mem%y_hat)*mem%y_hat & - -0.5*p%WaveField%WtrDens*mem%AxCa(FSElem) * (dSbdl_p*SaMGFSInt+dSadl_p*SbMGFSInt)*SQRT(SaMGFSInt*SbMGFSInt) * Dot_Product(SAFSInt,mem%k)*mem%k + 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 ! ------------------- hydrodynamic inertia loads: sides: Section 7.1.4 ------------------------ - F_IS= mem%Cp(FSElem)*p%WaveField%WtrDens* SaMGFSInt*SbMGFSInt * matmul( mem%Ak, FAFSInt ) + & ! transver FK component - FDynPFSInt*mem%AxCp(FSElem)* (SaMGFSInt*dSbdl_pp+dSadl_pp*SbMGFSInt) *mem%k + & ! axial FK component - p%WaveField%WtrDens*mem%CaB(FSElem) * 0.25*pi*SbMGFSInt*SbMGFSInt * Dot_Product(FAFSInt,mem%x_hat)*mem%x_hat + & ! x-component of diffraction part - p%WaveField%WtrDens*mem%CaA(FSElem) * 0.25*pi*SaMGFSInt*SaMGFSInt * Dot_Product(FAFSInt,mem%y_hat)*mem%y_hat + & ! y-component of diffraction part - 0.5*p%WaveField%WtrDens*mem%AxCa(FSElem) * (dSbdl_p*SaMGFSInt+dSadl_p*SbMGFSInt)*SQRT(SaMGFSInt*SbMGFSInt) * & ! axial component of diffraction part - Dot_Product(FAFSInt,mem%k)*mem%k ! axial component of diffraction part cont'd + 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 + p%WaveField%WtrDens*CaBFSInt * 0.25*pi*SbMGFSInt*SbMGFSInt * Dot_Product(FAFSInt,mem%x_hat)*mem%x_hat + & ! x-component of diffraction part + p%WaveField%WtrDens*CaAFSInt * 0.25*pi*SaMGFSInt*SaMGFSInt * Dot_Product(FAFSInt,mem%y_hat)*mem%y_hat + & ! y-component of diffraction part + 0.5*p%WaveField%WtrDens*AxCaFSInt * (dSbdl_p*SaMGFSInt+dSadl_p*SbMGFSInt)*SQRT(SaMGFSInt*SbMGFSInt) * & ! axial component of diffraction part + Dot_Product(FAFSInt,mem%k)*mem%k ! axial component of diffraction part cont'd END IF END IF @@ -4209,7 +4312,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ELSE IF ( MemSubStat .NE. 3_IntKi) THEN ! Skip members with centerline completely out of water !----------------------------No load smoothing----------------------------! 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)) + z1 = m%DispNodePosHdn(3, mem%NodeIndx(i)) + pos1 = m%DispNodePosHdn(:, mem%NodeIndx(i)) !---------------------------------------------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 @@ -4301,10 +4405,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, 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 - f_hydro = 0.5*mem%CdB(i)*p%WaveField%WtrDens*mem%SbMG(i)*TwoNorm(vec)*Dot_Product(vec,mem%x_hat)*mem%x_hat + & ! local x-direction - 0.5*mem%CdA(i)*p%WaveField%WtrDens*mem%SaMG(i)*TwoNorm(vec)*Dot_Product(vec,mem%y_hat)*mem%y_hat + & ! local z-direction - 0.25*mem%AxCd(i)*p%WaveField%WtrDens * (dSadl_p*mem%SbMG(i) + dSbdl_p*mem%SaMG(i)) * & ! axial part - abs(dot_product( mem%k, m%vrel(:,mem%NodeIndx(i)) )) * matmul( mem%kkt, m%vrel(:,mem%NodeIndx(i)) ) ! axial part cont'd + Call GetDistDrag_Rec(Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END IF 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) @@ -4649,33 +4750,6 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, CONTAINS - SUBROUTINE GetDisplacedNodePosition( forceDisplaced, pos ) - LOGICAL, INTENT( IN ) :: forceDisplaced ! Set to true to return the exact displaced position no matter WaveDisp or WaveStMod - REAL(ReKi), INTENT( OUT ) :: pos(:,:) ! Displaced node positions - REAL(ReKi) :: Orient(3,3) - - INTEGER(IntKi) :: ErrStat2 - CHARACTER(ErrMsgLen) :: ErrMsg2 - - ! Undisplaced node position - 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 - 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) - call GetPtfmRefYOrient(u%PtfmRefY, Orient, ErrStat2, ErrMsg2) - pos = matmul(transpose(Orient),pos) - END IF - - END SUBROUTINE GetDisplacedNodePosition - SUBROUTINE GetTotalWaveElev( Time, pos, Zeta, ErrStat, ErrMsg ) REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(*) ! Position at which free-surface elevation is to be calculated. Third entry ignored if present. @@ -5723,49 +5797,6 @@ SUBROUTINE getElementHstLds_Mod1( Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1, r2, END IF END SUBROUTINE getElementHstLds_Mod1 - SUBROUTINE YawMember(member, PtfmRefY, ErrStat, ErrMsg) - Type(Morison_MemberType), intent(inout) :: member - Real(ReKi), intent(in ) :: PtfmRefY - 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 - - Character(*), parameter :: RoutineName = 'YawMember' - - 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 - - call hiFrameTransform(h2i,PtfmRefY,member%y_hat,y_hat,ErrStat2,ErrMsg2) - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - member%y_hat = y_hat - - END IF - - END SUBROUTINE YawMember - SUBROUTINE YawJoint(JointNo,PtfmRefY,AM_End,An_End,DP_Const_End,I_MG_End,ErrStat,ErrMsg) Integer(IntKi), intent(in ) :: JointNo Real(ReKi), intent(in ) :: PtfmRefY @@ -5901,6 +5932,177 @@ SUBROUTINE getMemBallastHiPt(member,z_hi, ErrStat, ErrMsg) END IF END SUBROUTINE getMemBallastHiPt + + SUBROUTINE GetDistDrag_Rec(Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat,ErrMsg,SubRatio,vrelFSInt) + ! Compute the distributed (axial and transverse) drag per unit length for rectangular sections + Real(DbKi) , intent(in ) :: Time !< Current simulation time in seconds + Type(Morison_MemberType), intent(in ) :: mem !< Current member + Integer(IntKi) , intent(in ) :: i !< Node number within the member (not the global node index) + Real(ReKi) , intent(in ) :: dSadl_p !< Slope of Side A due to tapering (absolute value) + Real(ReKi) , intent(in ) :: dSbdl_p !< Slope of Side B due to tapering (absolute value) + Real(ReKi), optional , intent(in ) :: SubRatio !< Optional input. If provided, drag force will be evaluated at the free-surface intersection. + ! SubRatio is the fraction of element i (between node i and node i+1) submerged in water. + Real(ReKi), optional , intent(in ) :: vrelFSInt(3) !< Optional input. Must be provided if SubRatio is specified. + ! vrelFSInt is the fluid velocity relative to the structure at the free-surface intersection. + Real(ReKi) , intent( out) :: f_hydro(3) !< Sectional drag force (per unit length) + Integer(IntKi) , intent( out) :: ErrStat + Character(*) , intent( out) :: ErrMsg + + ! Local variables for alternative rectangular member transverse drag calculation + Integer(IntKi) :: NodeID ! Global node id of the ith node of the current member + Integer(IntKi) :: NextNodeID ! Global node id of the (i+1)th node of the current member + Integer(IntKi) :: fNo ! Face number + Integer(IntKi) :: tmpNodeInWater + Real(ReKi) :: pos(3) ! Node (on member axis) position + Real(ReKi) :: rToFC(3,4) ! Vectors from node on member axis to the four face centers + Real(ReKi) :: n_hat(3,4) ! Normal vectors of the four faces + Real(ReKi) :: filtConst(4) ! Velocity high-pass filter constants for the four faces + Real(ReKi) :: DragLoFSc(4) ! Drag weighting factors for the four faces + Real(ReKi) :: Cd(4) ! Dimensional drag coefficients for the four faces + Real(ReKi) :: vec(3) ! Relative velocity vector (flow-structure) without the axial component + Real(ReKi) :: posFC(3) ! Position of face center + Real(ReKi) :: SVFC(3) ! Structure velocity at face center + Real(SiKi) :: FVFC(3) ! Flow velocity at face center + Real(ReKi) :: vrelFC ! Relative (flow-structure) velocity at face centers + Real(ReKi) :: vrelFCf ! High-pass filtered relative (flow-structure) velocity at face centers + Real(ReKi) :: STV(3) ! Structure translational velocity at the current node/free-surface intersection + Real(ReKi) :: STV1(3) ! STV at free surface intersection computed from Node i velocity + Real(ReKi) :: STV2(3) ! STV at free surface intersection computed from Node i+1 velocity + Real(ReKi) :: SRV(3) ! Structure rotational velocity at the current node/free-surface intersection + Real(ReKi) :: FiltStat(4) ! High-pass filter states for the four faces + Real(ReKi) :: SaMG, SbMG ! Section side lengths at the current node/free-surface intersection + Real(ReKi) :: CdA, CdB, AxCd ! Drag coefficients at the current node/free-surface intersection + Real(ReKi) :: vrel(3) ! Relative flow velocity at the current node/free-surface intersection + + Integer(IntKi) :: ErrStat2 + Character(ErrMsgLen) :: ErrMsg2 + Character(*), parameter :: RoutineName = 'GetDistDrag_Rec' + + ErrStat = ErrID_None + ErrMsg = '' + + NodeID = mem%NodeIndx(i) + + f_hydro = 0.0_ReKi + IF ( m%NodeInWater(NodeID) .EQ. 0_IntKi ) THEN ! Node out of water + Return + END IF + + ! Node in water + IF ( PRESENT(SubRatio) ) THEN + ! Linearly interpolate the relevant parameters at the free-surface intersection + NextNodeID = mem%NodeIndx(i+1) + SaMG = SubRatio * mem%SaMG(i+1) + (1.0-SubRatio) * mem%SaMG(i) + SbMG = SubRatio * mem%SbMG(i+1) + (1.0-SubRatio) * mem%SbMG(i) + CdA = SubRatio * mem%CdA( i+1) + (1.0-SubRatio) * mem%CdA( i) + CdB = SubRatio * mem%CdB( i+1) + (1.0-SubRatio) * mem%CdB( i) + AxCd = SubRatio * mem%AxCd(i+1) + (1.0-SubRatio) * mem%AxCd(i) + vrel = vrelFSInt + ELSE + ! Use the relevant parameters at node i + SaMG = mem%SaMG(i) + SbMG = mem%SbMG(i) + CdA = mem%CdA( i) + CdB = mem%CdB( i) + AxCd = mem%AxCd(i) + vrel = m%vrel(:,NodeID) + END IF + + ! Axial drag + f_hydro = 0.25 * AxCd * p%WaveField%WtrDens * (dSadl_p*SbMG + dSbdl_p*SaMG) * & ! axial part + abs(dot_product( mem%k, vrel )) * matmul( mem%kkt, vrel ) ! axial part cont'd + + ! Transverse drag + IF (mem%FDMod == 0_IntKi) THEN ! Centerline-based formulation + + vec = matmul( mem%Ak,vrel ) + f_hydro = f_hydro + & + 0.5*CdB*p%WaveField%WtrDens*SbMG*TwoNorm(vec)*Dot_Product(vec,mem%x_hat)*mem%x_hat + & ! local x-direction + 0.5*CdA*p%WaveField%WtrDens*SaMG*TwoNorm(vec)*Dot_Product(vec,mem%y_hat)*mem%y_hat ! local y-direction + + ELSE ! Face-based formulation + + ! Position of node on member axis + IF ( PRESENT(SubRatio) ) THEN + pos = SubRatio * m%DispNodePosHdn( :,NextNodeID) + (1.0-SubRatio) * m%DispNodePosHdn( :,NodeID) + STV1 = u%Mesh%TranslationVel(:, NodeID) + CROSS_PRODUCT( u%Mesh%RotationVel(:, NodeID) , mem%dl * SubRatio * mem%k ) + STV2 = u%Mesh%TranslationVel(:,NextNodeID) + CROSS_PRODUCT( u%Mesh%RotationVel(:,NextNodeID) , - mem%dl * (1.0-SubRatio) * mem%k ) + STV = SubRatio * STV2 + (1.0-SubRatio) * STV1 + SRV = SubRatio * u%Mesh%RotationVel( :,NextNodeID) + (1.0-SubRatio) * u%Mesh%RotationVel( :,NodeID) + FiltStat = SubRatio * xd%MV_rel_n_FiltStat(:,NextNodeID) + (1.0-SubRatio) * xd%MV_rel_n_FiltStat(:,NodeID) + ELSE + pos = m%DispNodePosHdn( :,NodeID) + STV = u%Mesh%TranslationVel(:,NodeID) + SRV = u%Mesh%RotationVel( :,NodeID) + FiltStat = xd%MV_rel_n_FiltStat( :,NodeID) + END IF + + ! Vector from node on member axis to face centers + rToFC(1:3,1) = mem%x_hat * 0.5 * SaMG ! Side B +x_hat side + rToFC(1:3,2) = - mem%x_hat * 0.5 * SaMG ! Side B -x_hat side + rToFC(1:3,3) = mem%y_hat * 0.5 * SbMG ! Side A +y_hat side + rToFC(1:3,4) = - mem%y_hat * 0.5 * SbMG ! Side A -y_hat side + + ! Face normal vectors + n_hat(1:3,1) = mem%x_hat + n_hat(1:3,2) = - mem%x_hat + n_hat(1:3,3) = mem%y_hat + n_hat(1:3,4) = - mem%y_hat + + ! High-pass filter constant for each face + filtConst(1) = mem%VRelNFiltConstB + filtConst(2) = mem%VRelNFiltConstB + filtConst(3) = mem%VRelNFiltConstA + filtConst(4) = mem%VRelNFiltConstA + + ! Weighting factor for each face + DragLoFSc(1) = mem%DragLoFScB + DragLoFSc(2) = mem%DragLoFScB + DragLoFSc(3) = mem%DragLoFScA + DragLoFSc(4) = mem%DragLoFScA + + ! Dimensional drag coefficient for each face + Cd(1) = CdB*p%WaveField%WtrDens*SbMG + Cd(2) = Cd(1) + Cd(3) = CdA*p%WaveField%WtrDens*SaMG + Cd(4) = Cd(3) + + ! Compute and sum the drag force on all four faces + DO fNo = 1,4 + + ! Positions of face center + posFC = pos + rToFC(1:3,fNo) + + ! Compute structure velocity at face center + SVFC = STV + cross_product( SRV, rToFC(1:3,fNo) ) + + ! Compute fluid velocity at face center + Call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., tmpNodeInWater, FVFC, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + ! Note: We force each face center to also be wetted if the center node is wetted. Otherwise, the load-smoothing procedure might not work + + ! Compute the face-normal component of the relative fluid velocity (fluid-structure) at face center + vrelFC = dot_product( Real(FVFC,ReKi)-SVFC , n_hat(1:3,fNo) ) + ! High-pass-filtered face-center normal relative velocity + vrelFCf = filtConst(fNo) * (vrelFC + FiltStat(fNo)) + + ! Compute drag force based on selected formulation + IF ( mem%FDMod == 1_IntKi ) THEN ! Without suction-side-only formulation + f_hydro = f_hydro + n_hat(1:3,fNo) * ( & + (1.0_ReKi - DragLoFSc(fNo)) * 0.25*Cd(fNo)*abs(vrelFCf)*vrelFCf & + + DragLoFSc(fNo) * 0.25*Cd(fNo)*abs(vrelFC )*vrelFC ) + ELSE ! mem%FDMod == 2_IntKi With suction-side-only formulation + f_hydro = f_hydro + n_hat(1:3,fNo) * ( & + (1.0_ReKi - DragLoFSc(fNo)) * 0.50*Cd(fNo)*abs(vrelFCf)*max(vrelFCf,0.0_ReKi) & + + DragLoFSc(fNo) * 0.50*Cd(fNo)*abs(vrelFC )*max(vrelFC ,0.0_ReKi) ) + END IF + + END DO + + END IF + + END SUBROUTINE GetDistDrag_Rec + + logical function Failed() CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Failed = ErrStat >= AbortErrLev @@ -5988,10 +6190,13 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat TYPE(Morison_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables INTEGER(IntKi), INTENT( OUT) :: errStat !< Error status of the operation CHARACTER(*), INTENT( OUT) :: errMsg !< Error message if errStat /= ErrID_None - INTEGER(IntKi) :: J - INTEGER(IntKi) :: nodeInWater + 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) :: posFC(3), SVFC(3), vrelFC, vrelFCf REAL(SiKi) :: FVTmp(3) + TYPE(Morison_MemberType) :: mem !< Current member + INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 CHARACTER(*), PARAMETER :: RoutineName = 'Morison_UpdateDiscState' @@ -6000,24 +6205,14 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat errStat = ErrID_None errMsg = "" + + CALL GetDisplacedNodePosition( u, p, .FALSE., m%DispNodePosHdn ) ! For hydrodynamic loads; depends on WaveDisp and WaveStMod + ! Update state of the relative normal velocity high-pass filter at each joint DO J = 1, p%NJoints ! Get joint position - IF (p%WaveDisp == 0 ) THEN - ! use the initial X,Y location - pos(1) = u%Mesh%Position(1,J) - pos(2) = u%Mesh%Position(2,J) - ELSE - ! Use current X,Y location - pos(1) = u%Mesh%TranslationDisp(1,J) + u%Mesh%Position(1,J) - pos(2) = u%Mesh%TranslationDisp(2,J) + u%Mesh%Position(2,J) - END IF - IF (p%WaveField%WaveStMod > 0 .AND. p%WaveDisp /= 0) THEN ! Wave stretching enabled - pos(3) = u%Mesh%Position(3,J) + u%Mesh%TranslationDisp(3,J) - p%WaveField%MSL2SWL ! Use the current Z location. - ELSE ! Wave stretching disabled - pos(3) = u%Mesh%Position(3,J) - p%WaveField%MSL2SWL ! We are intentionally using the undisplaced Z position of the node. - END IF + pos = m%DispNodePosHdn(:,J) ! Get fluid velocity at the joint CALL WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, pos, .FALSE., nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) @@ -6038,6 +6233,77 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat END DO ! J = 1, p%NJoints + ! Update state of the relative normal velocity high-pass filter for each rectangular member + DO im = 1, p%NMembers + IF ( (p%Members(im)%MSecGeom == MSecGeom_Rec) .and. (p%Members(im)%FDMod > 0_IntKi) ) THEN + + 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., nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + + IF (nodeInWater .EQ. 1_IntKi) THEN + + ! 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., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + 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 + + ! 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., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + 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 + + ! 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., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + 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 + + ! 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., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + 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 + + ELSE + + vrelFC = 0.0_ReKi + + vrelFCf = mem%VRelNFiltConstB * ( vrelFC + xd%MV_rel_n_FiltStat(1,mem%NodeIndx(I)) ) + xd%MV_rel_n_FiltStat(1,mem%NodeIndx(I)) = vrelFCf - vrelFC + + vrelFCf = mem%VRelNFiltConstB * ( vrelFC + xd%MV_rel_n_FiltStat(2,mem%NodeIndx(I)) ) + xd%MV_rel_n_FiltStat(2,mem%NodeIndx(I)) = vrelFCf - vrelFC + + vrelFCf = mem%VRelNFiltConstA * ( vrelFC + xd%MV_rel_n_FiltStat(3,mem%NodeIndx(I)) ) + xd%MV_rel_n_FiltStat(3,mem%NodeIndx(I)) = vrelFCf - vrelFC + + vrelFCf = mem%VRelNFiltConstA * ( vrelFC + xd%MV_rel_n_FiltStat(4,mem%NodeIndx(I)) ) + xd%MV_rel_n_FiltStat(4,mem%NodeIndx(I)) = vrelFCf - vrelFC + + END IF + END DO ! Iterate through member nodes + END IF ! If rectangular member + END DO ! Iterate through members + END SUBROUTINE Morison_UpdateDiscState !---------------------------------------------------------------------------------------------------------------------------------- END MODULE Morison diff --git a/modules/hydrodyn/src/Morison.txt b/modules/hydrodyn/src/Morison.txt index a6b4a536a3..41fb1fbdc8 100644 --- a/modules/hydrodyn/src/Morison.txt +++ b/modules/hydrodyn/src/Morison.txt @@ -101,6 +101,11 @@ typedef ^ ^ INTEGER typedef ^ ^ INTEGER MmbrFilledIDIndx - - - "Index into the filled group table if this is a filled member" - typedef ^ ^ LOGICAL PropPot - - - "Flag T/F for whether the member is modeled with potential flow theory" - typedef ^ ^ LOGICAL PropMCF - - - "Flag T/F for whether the member is modeled with the MacCamy-Fuchs diffraction model" - +typedef ^ ^ INTEGER FDMod - - - "Rectangular member transverse drag model (0: simple centerline based; 1: face based; 2: face based with suction-side-only formulation" - +typedef ^ ^ ReKi VnCOffA - - - "Rectangular member transverse drag relative velocity high-pass filter cutoff frequency - normal to Side A" Hz +typedef ^ ^ ReKi VnCOffB - - - "Rectangular member transverse drag relative velocity high-pass filter cutoff frequency - normal to Side B" Hz +typedef ^ ^ ReKi FDLoFScA - - - "Rectangular member transverse drag weighting factor - normal to Side A" - +typedef ^ ^ ReKi FDLoFScB - - - "Rectangular member transverse drag weighting factor - normal to Side B" - typedef ^ ^ INTEGER NElements - - - "number of elements in this member" - typedef ^ ^ ReKi RefLength - - - "the reference total length for this member" m typedef ^ ^ ReKi dl - - - "the reference element length for this member (may be less than MDivSize to achieve uniform element lengths)" m @@ -182,6 +187,11 @@ typedef ^ ^ ReKi typedef ^ ^ ReKi AxCa {:} - - "Member axial Ca at each node" - typedef ^ ^ ReKi AxCp {:} - - "Member axial Cp at each node" - typedef ^ ^ ReKi Cb {:} - - "Member Cb at each node" - +typedef ^ ^ IntKi FDMod - - - "Rectangular member transverse drag model (0: simple centerline based; 1: face based; 2: face based with suction-side-only formulation" - +typedef ^ ^ ReKi VRelNFiltConstA - - - "Rectangular member transverse drag relative velocity high-pass filter constant - normal to Side A" - +typedef ^ ^ ReKi VRelNFiltConstB - - - "Rectangular member transverse drag relative velocity high-pass filter constant - normal to Side B" - +typedef ^ ^ ReKi DragLoFScA - - - "Rectangular member transverse drag weighting factor - normal to Side A" - +typedef ^ ^ ReKi DragLoFScB - - - "Rectangular member transverse drag weighting factor - normal to Side B" - typedef ^ ^ ReKi m_fb_l {:} - - "mass of flooded ballast in lower portion of each element" kg typedef ^ ^ ReKi m_fb_u {:} - - "mass of flooded ballast in upper portion of each element" kg typedef ^ ^ ReKi h_cfb_l {:} - - "distance to flooded ballast centroid from node point in lower portion of each element" m @@ -416,6 +426,7 @@ typedef ^ ContinuousStateType SiKi # Define discrete (nondifferentiable) states here: # typedef ^ DiscreteStateType ReKi V_rel_n_FiltStat {:} - - "State of the high-pass filter for the joint relative normal velocity" m/s +typedef ^ DiscreteStateType ReKi MV_rel_n_FiltStat {:}{:} - - "State of the high-pass filter for the rectangular member relative normal velocity" m/s # # # Define constraint states here: diff --git a/modules/hydrodyn/src/Morison_Types.f90 b/modules/hydrodyn/src/Morison_Types.f90 index 0cbf7b76aa..44f27a819d 100644 --- a/modules/hydrodyn/src/Morison_Types.f90 +++ b/modules/hydrodyn/src/Morison_Types.f90 @@ -148,6 +148,11 @@ MODULE Morison_Types INTEGER(IntKi) :: MmbrFilledIDIndx = 0_IntKi !< Index into the filled group table if this is a filled member [-] LOGICAL :: PropPot = .false. !< Flag T/F for whether the member is modeled with potential flow theory [-] LOGICAL :: PropMCF = .false. !< Flag T/F for whether the member is modeled with the MacCamy-Fuchs diffraction model [-] + INTEGER(IntKi) :: FDMod = 0_IntKi !< Rectangular member transverse drag model (0: simple centerline based; 1: face based; 2: face based with suction-side-only formulation [-] + REAL(ReKi) :: VnCOffA = 0.0_ReKi !< Rectangular member transverse drag relative velocity high-pass filter cutoff frequency - normal to Side A [Hz] + REAL(ReKi) :: VnCOffB = 0.0_ReKi !< Rectangular member transverse drag relative velocity high-pass filter cutoff frequency - normal to Side B [Hz] + REAL(ReKi) :: FDLoFScA = 0.0_ReKi !< Rectangular member transverse drag weighting factor - normal to Side A [-] + REAL(ReKi) :: FDLoFScB = 0.0_ReKi !< Rectangular member transverse drag weighting factor - normal to Side B [-] 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) :: dl = 0.0_ReKi !< the reference element length for this member (may be less than MDivSize to achieve uniform element lengths) [m] @@ -235,6 +240,11 @@ MODULE Morison_Types REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: AxCa !< Member axial Ca at each node [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: AxCp !< Member axial Cp at each node [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: Cb !< Member Cb at each node [-] + INTEGER(IntKi) :: FDMod = 0_IntKi !< Rectangular member transverse drag model (0: simple centerline based; 1: face based; 2: face based with suction-side-only formulation [-] + REAL(ReKi) :: VRelNFiltConstA = 0.0_ReKi !< Rectangular member transverse drag relative velocity high-pass filter constant - normal to Side A [-] + REAL(ReKi) :: VRelNFiltConstB = 0.0_ReKi !< Rectangular member transverse drag relative velocity high-pass filter constant - normal to Side B [-] + REAL(ReKi) :: DragLoFScA = 0.0_ReKi !< Rectangular member transverse drag weighting factor - normal to Side A [-] + REAL(ReKi) :: DragLoFScB = 0.0_ReKi !< Rectangular member transverse drag weighting factor - normal to Side B [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: m_fb_l !< mass of flooded ballast in lower portion of each element [kg] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: m_fb_u !< mass of flooded ballast in upper portion of each element [kg] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: h_cfb_l !< distance to flooded ballast centroid from node point in lower portion of each element [m] @@ -483,6 +493,7 @@ MODULE Morison_Types ! ========= Morison_DiscreteStateType ======= TYPE, PUBLIC :: Morison_DiscreteStateType REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: V_rel_n_FiltStat !< State of the high-pass filter for the joint relative normal velocity [m/s] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: MV_rel_n_FiltStat !< State of the high-pass filter for the rectangular member relative normal velocity [m/s] END TYPE Morison_DiscreteStateType ! ======================= ! ========= Morison_ConstraintStateType ======= @@ -1064,6 +1075,11 @@ subroutine Morison_CopyMemberInputType(SrcMemberInputTypeData, DstMemberInputTyp DstMemberInputTypeData%MmbrFilledIDIndx = SrcMemberInputTypeData%MmbrFilledIDIndx DstMemberInputTypeData%PropPot = SrcMemberInputTypeData%PropPot DstMemberInputTypeData%PropMCF = SrcMemberInputTypeData%PropMCF + DstMemberInputTypeData%FDMod = SrcMemberInputTypeData%FDMod + DstMemberInputTypeData%VnCOffA = SrcMemberInputTypeData%VnCOffA + DstMemberInputTypeData%VnCOffB = SrcMemberInputTypeData%VnCOffB + DstMemberInputTypeData%FDLoFScA = SrcMemberInputTypeData%FDLoFScA + DstMemberInputTypeData%FDLoFScB = SrcMemberInputTypeData%FDLoFScB DstMemberInputTypeData%NElements = SrcMemberInputTypeData%NElements DstMemberInputTypeData%RefLength = SrcMemberInputTypeData%RefLength DstMemberInputTypeData%dl = SrcMemberInputTypeData%dl @@ -1105,6 +1121,11 @@ subroutine Morison_PackMemberInputType(RF, Indata) call RegPack(RF, InData%MmbrFilledIDIndx) call RegPack(RF, InData%PropPot) call RegPack(RF, InData%PropMCF) + call RegPack(RF, InData%FDMod) + call RegPack(RF, InData%VnCOffA) + call RegPack(RF, InData%VnCOffB) + call RegPack(RF, InData%FDLoFScA) + call RegPack(RF, InData%FDLoFScB) call RegPack(RF, InData%NElements) call RegPack(RF, InData%RefLength) call RegPack(RF, InData%dl) @@ -1138,6 +1159,11 @@ subroutine Morison_UnPackMemberInputType(RF, OutData) call RegUnpack(RF, OutData%MmbrFilledIDIndx); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%PropPot); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%PropMCF); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%FDMod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%VnCOffA); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%VnCOffB); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%FDLoFScA); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%FDLoFScB); 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%dl); if (RegCheckErr(RF, RoutineName)) return @@ -1724,6 +1750,11 @@ subroutine Morison_CopyMemberType(SrcMemberTypeData, DstMemberTypeData, CtrlCode end if DstMemberTypeData%Cb = SrcMemberTypeData%Cb end if + DstMemberTypeData%FDMod = SrcMemberTypeData%FDMod + DstMemberTypeData%VRelNFiltConstA = SrcMemberTypeData%VRelNFiltConstA + DstMemberTypeData%VRelNFiltConstB = SrcMemberTypeData%VRelNFiltConstB + DstMemberTypeData%DragLoFScA = SrcMemberTypeData%DragLoFScA + DstMemberTypeData%DragLoFScB = SrcMemberTypeData%DragLoFScB if (allocated(SrcMemberTypeData%m_fb_l)) then LB(1:1) = lbound(SrcMemberTypeData%m_fb_l) UB(1:1) = ubound(SrcMemberTypeData%m_fb_l) @@ -2291,6 +2322,11 @@ subroutine Morison_PackMemberType(RF, Indata) call RegPackAlloc(RF, InData%AxCa) call RegPackAlloc(RF, InData%AxCp) call RegPackAlloc(RF, InData%Cb) + call RegPack(RF, InData%FDMod) + call RegPack(RF, InData%VRelNFiltConstA) + call RegPack(RF, InData%VRelNFiltConstB) + call RegPack(RF, InData%DragLoFScA) + call RegPack(RF, InData%DragLoFScB) call RegPackAlloc(RF, InData%m_fb_l) call RegPackAlloc(RF, InData%m_fb_u) call RegPackAlloc(RF, InData%h_cfb_l) @@ -2400,6 +2436,11 @@ subroutine Morison_UnPackMemberType(RF, OutData) call RegUnpackAlloc(RF, OutData%AxCa); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%AxCp); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Cb); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%FDMod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%VRelNFiltConstA); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%VRelNFiltConstB); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%DragLoFScA); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%DragLoFScB); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%m_fb_l); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%m_fb_u); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%h_cfb_l); if (RegCheckErr(RF, RoutineName)) return @@ -4236,7 +4277,7 @@ subroutine Morison_CopyDiscState(SrcDiscStateData, DstDiscStateData, CtrlCode, E integer(IntKi), intent(in ) :: CtrlCode integer(IntKi), intent( out) :: ErrStat character(*), intent( out) :: ErrMsg - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: LB(2), UB(2) integer(IntKi) :: ErrStat2 character(*), parameter :: RoutineName = 'Morison_CopyDiscState' ErrStat = ErrID_None @@ -4253,6 +4294,18 @@ subroutine Morison_CopyDiscState(SrcDiscStateData, DstDiscStateData, CtrlCode, E end if DstDiscStateData%V_rel_n_FiltStat = SrcDiscStateData%V_rel_n_FiltStat end if + if (allocated(SrcDiscStateData%MV_rel_n_FiltStat)) then + LB(1:2) = lbound(SrcDiscStateData%MV_rel_n_FiltStat) + UB(1:2) = ubound(SrcDiscStateData%MV_rel_n_FiltStat) + if (.not. allocated(DstDiscStateData%MV_rel_n_FiltStat)) then + allocate(DstDiscStateData%MV_rel_n_FiltStat(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstDiscStateData%MV_rel_n_FiltStat.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstDiscStateData%MV_rel_n_FiltStat = SrcDiscStateData%MV_rel_n_FiltStat + end if end subroutine subroutine Morison_DestroyDiscState(DiscStateData, ErrStat, ErrMsg) @@ -4265,6 +4318,9 @@ subroutine Morison_DestroyDiscState(DiscStateData, ErrStat, ErrMsg) if (allocated(DiscStateData%V_rel_n_FiltStat)) then deallocate(DiscStateData%V_rel_n_FiltStat) end if + if (allocated(DiscStateData%MV_rel_n_FiltStat)) then + deallocate(DiscStateData%MV_rel_n_FiltStat) + end if end subroutine subroutine Morison_PackDiscState(RF, Indata) @@ -4273,6 +4329,7 @@ subroutine Morison_PackDiscState(RF, Indata) character(*), parameter :: RoutineName = 'Morison_PackDiscState' if (RF%ErrStat >= AbortErrLev) return call RegPackAlloc(RF, InData%V_rel_n_FiltStat) + call RegPackAlloc(RF, InData%MV_rel_n_FiltStat) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -4280,11 +4337,12 @@ subroutine Morison_UnPackDiscState(RF, OutData) type(RegFile), intent(inout) :: RF type(Morison_DiscreteStateType), intent(inout) :: OutData character(*), parameter :: RoutineName = 'Morison_UnPackDiscState' - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: LB(2), UB(2) integer(IntKi) :: stat logical :: IsAllocAssoc if (RF%ErrStat /= ErrID_None) return call RegUnpackAlloc(RF, OutData%V_rel_n_FiltStat); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%MV_rel_n_FiltStat); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine Morison_CopyConstrState(SrcConstrStateData, DstConstrStateData, CtrlCode, ErrStat, ErrMsg)