From ddb547d9f83e191d76dda3986e4f94dfca5928d6 Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Thu, 3 Apr 2025 18:53:22 -0600 Subject: [PATCH 1/3] MD: Fixes viscoelastic model initalization bug, adds error check for damping ratio, adds flag to disable calls to WrOver in initalization that caused MATLAB issues --- modules/moordyn/src/MoorDyn.f90 | 12 +++-- modules/moordyn/src/MoorDyn_Line.f90 | 60 ++++++++++++------------ modules/moordyn/src/MoorDyn_Registry.txt | 1 + modules/moordyn/src/MoorDyn_Types.f90 | 4 ++ 4 files changed, 44 insertions(+), 33 deletions(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index ca98b7e23e..e4850c9020 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -484,6 +484,8 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er read (OptValue,*) p%inertialF_rampT else if ( OptString == 'OUTSWITCH') then read (OptValue,*) p%OutSwitch + else if ( OptString == 'DISABLEOUTTIME') then + read (OptValue,*) p%disableOutTime else CALL SetErrStat( ErrID_Warn, 'Unable to interpret input '//trim(OptString)//' in OPTIONS section.', ErrStat, ErrMsg, RoutineName ) end if @@ -2789,10 +2791,12 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! provide status message - ! bjj: putting this in a string so we get blanks to cover up previous values (if current string is shorter than previous one) - Message = ' t='//trim(Num2LStr(t))//' FairTen 1: '//trim(Num2LStr(FairTensIC(1,1)))// & - ', '//trim(Num2LStr(FairTensIC(1,2)))//', '//trim(Num2LStr(FairTensIC(1,3))) - CALL WrOver( Message ) + IF (p%disableOutTime == 0) THEN ! option to turn this off if users want (helpful for matlab) + ! bjj: putting this in a string so we get blanks to cover up previous values (if current string is shorter than previous one) + Message = ' t='//trim(Num2LStr(t))//' FairTen 1: '//trim(Num2LStr(FairTensIC(1,1)))// & + ', '//trim(Num2LStr(FairTensIC(1,2)))//', '//trim(Num2LStr(FairTensIC(1,3))) + CALL WrOver( Message ) + ENDIF ! check for convergence (compare current tension at each fairlead with previous 9 values) IF (I > 9) THEN diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index a001d1597e..f40d482a75 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -121,6 +121,12 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) ! If input value is negative, it is considered to be desired damping ratio (zeta) ! from which the line's BA can be calculated based on the segment natural frequency. IF (LineProp%BA < 0) THEN + IF (Line%nEApoints > 0) THEN + ErrMsg = ' Line damping cannot be a ratio if stress-strain lookup table is used.' + ErrStat = ErrID_Fatal + RETURN + ENDIF + ! - we assume desired damping coefficient is zeta = -LineProp%BA ! - highest axial vibration mode of a segment is wn = sqrt(k/m) = 2N/UnstrLen*sqrt(EA/w) Line%BA = -LineProp%BA * Line%UnstrLen / Line%N * SQRT(LineProp%EA * LineProp%w) @@ -188,10 +194,6 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) !CALL CleanUp() RETURN END IF - ! initialize to unique values on the range 0-2pi - do I=0, Line%N - Line%phi(I) = (REAL(I)/Line%N)*2*Pi - enddo ! initialize other things to 0 Line%rdd_old = 0.0_DbKi @@ -410,6 +412,23 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg) ENDIF + ! If using the viscoelastic model, initalize the deltaL_1 to the delta L of the segment. + ! This is required here to initalize the state as non-zero, which avoids an initial + ! transient where the segment goes from unstretched to stretched in one time step. + IF (Line%ElasticMod > 1) THEN + DO I = 1, N + ! calculate current (Stretched) segment lengths and unit tangent vectors (qs) for each segment (this is used for bending calculations) + CALL UnitVector(Line%r(:,I-1), Line%r(:,I), Line%qs(:,I), Line%lstr(I)) + Line%dl_1(I) = Line%lstr(I) - Line%l(I) ! delta l of the segment + END DO + ENDIF + + ! initialize phi to unique values on the range 0-2pi if VIV model is active + IF (Line%Cl > 0) THEN + DO I=0, Line%N + Line%phi(I) = (REAL(I)/Line%N)*2*Pi + ENDDO + ENDIF CALL CleanUp() ! deallocate temporary arrays @@ -1349,7 +1368,7 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! Double check none of the assumptions were violated (this should never happen) IF (Line%alphaMBL <= 0 .OR. Line%vbeta <= 0 .OR. Line%l(I) <= 0 .OR. Line%dl_1(I) <= 0 .OR. EA_D < Line%EA) THEN ErrStat = ErrID_Warn - ErrMsg = "Viscoelastic model: Assumption for mean laod dependent dynamic stiffness violated" + ErrMsg = "Viscoelastic model: Assumption for mean load dependent dynamic stiffness violated" if (wordy > 2) then print *, "Line%alphaMBL", Line%alphaMBL print *, "Line%vbeta", Line%vbeta @@ -1535,8 +1554,13 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! Vortex Induced Vibration (VIV) cross-flow lift force Line%Lf(:,I) = 0.0_DbKi ! Zero lift force - IF ((Line%Cl > 0.0) .AND. (.NOT. Line%IC_gen)) THEN ! If non-zero lift coefficient and not during IC_gen ! Ignore the following: and internal node then VIV to be calculated .AND. (I /= 0) .AND. (I /= N) + IF ((Line%Cl > 0.0) .AND. (.NOT. Line%IC_gen)) THEN ! If non-zero lift coefficient and not during IC_gen + ! Note: This logic is slightly different than MD-C, but equivalent. MD-C runs the VIV model for only the internal + ! nodes. That means in MD-F the state vector has N+1 extra states when using the VIV model while the MD-C state + ! matrix can keep N-1 rows. That approach means in the future if we can get the end node accelerations, we can + ! easily add those into the VIV model. MD-C does not do that, and only computes the lift force for internal nodes. + ! ----- The Synchronization Model ------ ! Crossflow velocity and acceleration. rd component in the crossflow direction @@ -1563,31 +1587,11 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Line%phi_yd = 2*Pi + Line%phi_yd ! atan2 to 0-2Pi range ENDIF - ! Note: amplitude calculations and states commented out. Would be needed if a Cl vs A lookup table was ever implemented - - ! const real A_int = Misc(I)[1]; - ! const real As = Misc(I)[2]; - ! non-dimensional frequency f_hat = Line%cF + Line%dF *sin(Line%phi_yd - Line%phi(I)) ! phi is integrated from state deriv phi_dot ! frequency of lift force (rad/s) phi_dot = 2*Pi*f_hat*MagVp / d ! to be added to state - - ! ----- The rest of the model ----- - - ! ! Oscillation amplitude - ! const real A_int_dot = abs(yd); - ! ! Note: Check if this actually measures zero crossings - ! if ((yd * yd_old[i]) < 0) { ! if sign changed, i.e. a zero crossing - ! Amp[i] = A_int-A_int_old[i]; ! amplitude calculation since last zero crossing - ! A_int_old[i] = A_int; ! stores amplitude of previous zero crossing for finding Amp - ! } - ! ! Careful with integrating smoothed amplitude, as 0.1 was a calibarated value based on a very simple integration method - ! const real As_dot = (0.1/dtm)*(Amp[i]-As); ! As to be variable integrated from the state. stands for amplitude smoothed - - ! ! Lift coefficient from lookup table - ! const real C_l = cl_lookup(x = As/d); ! create function in Line.hpp that uses lookup table - + ! The Lift force IF (I==0) THEN ! Disable for end nodes for now, node acceleration needed for synch model Line%Lf(:,0) = 0.0_DbKi @@ -1610,8 +1614,6 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai else ! if only VIV state, then N+1 entries after internal node states Xd( 6*Line%N-6 + I + 1) = phi_dot endif - - ! Miscd(I)[2] = As_dot; ! unused state that could be used for future amplitude calculations ENDIF diff --git a/modules/moordyn/src/MoorDyn_Registry.txt b/modules/moordyn/src/MoorDyn_Registry.txt index 95adb6df1d..60a4ad7c8e 100644 --- a/modules/moordyn/src/MoorDyn_Registry.txt +++ b/modules/moordyn/src/MoorDyn_Registry.txt @@ -450,6 +450,7 @@ typedef ^ ^ DbKi cv - typedef ^ ^ IntKi inertialF - 0 - "Indicates MoorDyn returning inertial moments for coupled 6DOF objects. 0: no, 1: yes, 2: yes with ramp to inertialF_rampT" - typedef ^ ^ R8Ki inertialF_rampT - 30 - "Ramp time for inertial forces" - typedef ^ ^ IntKi OutSwitch - 1 - "Switch to disable outputs when running with full OF. 0: no MD main outfile, 1: write MD main outfile" "(-)" +typedef ^ ^ IntKi disableOutTime - 0 - "Disables the printing of the fairtens and timestep in init to the console, useful for the MATLAB wrapper" "(-)" # --- parameters for wave and current --- typedef ^ ^ IntKi nxWave - - - "number of x wave grid points" - typedef ^ ^ IntKi nyWave - - - "number of y wave grid points" - diff --git a/modules/moordyn/src/MoorDyn_Types.f90 b/modules/moordyn/src/MoorDyn_Types.f90 index 63eba5f6ba..cd85475cf6 100644 --- a/modules/moordyn/src/MoorDyn_Types.f90 +++ b/modules/moordyn/src/MoorDyn_Types.f90 @@ -485,6 +485,7 @@ MODULE MoorDyn_Types INTEGER(IntKi) :: inertialF = 0 !< Indicates MoorDyn returning inertial moments for coupled 6DOF objects. 0: no, 1: yes, 2: yes with ramp to inertialF_rampT [-] REAL(R8Ki) :: inertialF_rampT = 30 !< Ramp time for inertial forces [-] INTEGER(IntKi) :: OutSwitch = 1 !< Switch to disable outputs when running with full OF. 0: no MD main outfile, 1: write MD main outfile [(-)] + INTEGER(IntKi) :: disableOutTime = 0 !< Disables the printing of the fairtens and timestep in init to the console, useful for the MATLAB wrapper [(-)] INTEGER(IntKi) :: nxWave = 0_IntKi !< number of x wave grid points [-] INTEGER(IntKi) :: nyWave = 0_IntKi !< number of y wave grid points [-] INTEGER(IntKi) :: nzWave = 0_IntKi !< number of z wave grid points [-] @@ -3982,6 +3983,7 @@ subroutine MD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%inertialF = SrcParamData%inertialF DstParamData%inertialF_rampT = SrcParamData%inertialF_rampT DstParamData%OutSwitch = SrcParamData%OutSwitch + DstParamData%disableOutTime = SrcParamData%disableOutTime DstParamData%nxWave = SrcParamData%nxWave DstParamData%nyWave = SrcParamData%nyWave DstParamData%nzWave = SrcParamData%nzWave @@ -4385,6 +4387,7 @@ subroutine MD_PackParam(RF, Indata) call RegPack(RF, InData%inertialF) call RegPack(RF, InData%inertialF_rampT) call RegPack(RF, InData%OutSwitch) + call RegPack(RF, InData%disableOutTime) call RegPack(RF, InData%nxWave) call RegPack(RF, InData%nyWave) call RegPack(RF, InData%nzWave) @@ -4493,6 +4496,7 @@ subroutine MD_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%inertialF); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%inertialF_rampT); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%OutSwitch); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%disableOutTime); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nxWave); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nyWave); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nzWave); if (RegCheckErr(RF, RoutineName)) return From 0c022ff2c163af7e9ac1c86fe13b9889f76e3ff1 Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Tue, 8 Apr 2025 11:55:27 -0600 Subject: [PATCH 2/3] Update OpenFAST and FAST.Farm to give MoorDyn access to the ED Ptfm 6 DOF at initalization --- glue-codes/fast-farm/src/FAST_Farm_Subs.f90 | 6 +++--- modules/openfast-library/src/FAST_Registry.txt | 1 + modules/openfast-library/src/FAST_Subs.f90 | 5 ++++- modules/openfast-library/src/FAST_Types.f90 | 4 ++++ 4 files changed, 12 insertions(+), 4 deletions(-) diff --git a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 index 0a7f1a60b0..0e00862963 100644 --- a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 +++ b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 @@ -834,9 +834,9 @@ SUBROUTINE Farm_InitMD( farm, ErrStat, ErrMsg ) ALLOCATE( MD_InitInp%PtfmInit(6,farm%p%NumTurbines), MD_InitInp%TurbineRefPos(3,farm%p%NumTurbines), STAT = ErrStat2 ) if (Failed0("MoorDyn PtfmInit and TurbineRefPos initialization inputs in FAST.Farm.")) return; - ! gather spatial initialization inputs for Farm-level MoorDyn - DO nt = 1,farm%p%NumTurbines - MD_InitInp%PtfmInit(:,nt) = farm%FWrap(nt)%m%Turbine%MD%m%PtfmInit ! turbine PRP initial positions and rotations in their respective coordinate systems from each FAST/MD instance + ! gather spatial initialization inputs for Farm-level MoorDyn (platform locations in their respective coordinate systems and locations of the turbines in the farm global coordinate system) + DO nt = 1,farm%p%NumTurbines + MD_InitInp%PtfmInit(:,nt) = farm%FWrap(nt)%m%Turbine%p_FAST%PlatformPosInit ! platform initial positions in their respective coordinate systems from each FAST/ED instance MD_InitInp%TurbineRefPos(:,nt) = farm%p%WT_Position(:,nt) ! reference positions of each turbine in the farm global coordinate system END DO diff --git a/modules/openfast-library/src/FAST_Registry.txt b/modules/openfast-library/src/FAST_Registry.txt index 2cbaa53e41..cf20465772 100644 --- a/modules/openfast-library/src/FAST_Registry.txt +++ b/modules/openfast-library/src/FAST_Registry.txt @@ -193,6 +193,7 @@ typedef ^ FAST_ParameterType INTEGER VTK_tWidth - - - "Width of number of files typedef ^ FAST_ParameterType DbKi VTK_fps - - - "number of frames per second to output VTK data" - typedef ^ FAST_ParameterType FAST_VTK_SurfaceType VTK_surface - - - "Data for VTK surface visualization" typedef ^ FAST_ParameterType CHARACTER(4) Tdesc - - - "description of turbine ID (for FAST.Farm) screen printing" +typedef ^ FAST_ParameterType DbKi PlatformPosInit {6} - - "Platform inital 6 DOF position from ED (this is different from TurbinePos)" - # Parameters for linearization typedef ^ FAST_ParameterType LOGICAL CalcSteady - - - "Calculate a steady-state periodic operating point before linearization [unused if Linearize=False]" - diff --git a/modules/openfast-library/src/FAST_Subs.f90 b/modules/openfast-library/src/FAST_Subs.f90 index 1b0050734c..7b94ec9e78 100644 --- a/modules/openfast-library/src/FAST_Subs.f90 +++ b/modules/openfast-library/src/FAST_Subs.f90 @@ -290,6 +290,9 @@ SUBROUTINE FAST_InitializeAll( t_initial, p_FAST, y_FAST, m_FAST, ED, SED, BD, S CALL ED_Init( Init%InData_ED, ED%Input(1), ED%p, ED%x(STATE_CURR), ED%xd(STATE_CURR), ED%z(STATE_CURR), ED%OtherSt(STATE_CURR), & ED%y, ED%m, p_FAST%dt_module( MODULE_ED ), Init%OutData_ED, ErrStat2, ErrMsg2 ) CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + + ! Assign the inital positions for use by MoorDyn initalization + p_FAST%PlatformPosInit = Init%OutData_ED%PlatformPos p_FAST%ModuleInitialized(Module_ED) = .TRUE. CALL SetModuleSubstepTime(Module_ED, p_FAST, y_FAST, ErrStat2, ErrMsg2) @@ -1207,7 +1210,7 @@ SUBROUTINE FAST_InitializeAll( t_initial, p_FAST, y_FAST, m_FAST, ED, SED, BD, S Init%InData_MD%FileName = p_FAST%MooringFile ! This needs to be set according to what is in the FAST input file. Init%InData_MD%RootName = p_FAST%OutFileRoot - Init%InData_MD%PtfmInit(:,1) = Init%OutData_ED%PlatformPos ! initial position of the platform (when a FAST module, MoorDyn just takes one row in this matrix) + Init%InData_MD%PtfmInit(:,1) = p_FAST%PlatformPosInit ! initial position of the platform (when a FAST module, MoorDyn just takes one row in this matrix) Init%InData_MD%FarmSize = 0 ! 0 here indicates normal FAST module use of MoorDyn, for a single turbine Init%InData_MD%TurbineRefPos(:,1) = 0.0_DbKi ! for normal FAST use, the global reference frame is at 0,0,0 Init%InData_MD%g = p_FAST%Gravity ! This need to be according to g used in ElastoDyn diff --git a/modules/openfast-library/src/FAST_Types.f90 b/modules/openfast-library/src/FAST_Types.f90 index 125d2a0643..26c6582ae2 100644 --- a/modules/openfast-library/src/FAST_Types.f90 +++ b/modules/openfast-library/src/FAST_Types.f90 @@ -207,6 +207,7 @@ MODULE FAST_Types REAL(DbKi) :: VTK_fps = 0.0_R8Ki !< number of frames per second to output VTK data [-] TYPE(FAST_VTK_SurfaceType) :: VTK_surface !< Data for VTK surface visualization [-] CHARACTER(4) :: Tdesc !< description of turbine ID (for FAST.Farm) screen printing [-] + REAL(DbKi) , DIMENSION(1:6) :: PlatformPosInit = 0.0_R8Ki !< Platform inital 6 DOF position from ED (this is different from TurbinePos) [-] LOGICAL :: CalcSteady = .false. !< Calculate a steady-state periodic operating point before linearization [unused if Linearize=False] [-] INTEGER(IntKi) :: TrimCase = 0_IntKi !< Controller parameter to be trimmed {1:yaw; 2:torque; 3:pitch} [unused if Linearize=False; used only if CalcSteady=True] [-] REAL(ReKi) :: TrimTol = 0.0_ReKi !< Tolerance for the rotational speed convergence (>0) [unused if Linearize=False; used only if CalcSteady=True] [-] @@ -1487,6 +1488,7 @@ subroutine FAST_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return DstParamData%Tdesc = SrcParamData%Tdesc + DstParamData%PlatformPosInit = SrcParamData%PlatformPosInit DstParamData%CalcSteady = SrcParamData%CalcSteady DstParamData%TrimCase = SrcParamData%TrimCase DstParamData%TrimTol = SrcParamData%TrimTol @@ -1657,6 +1659,7 @@ subroutine FAST_PackParam(RF, Indata) call RegPack(RF, InData%VTK_fps) call FAST_PackVTK_SurfaceType(RF, InData%VTK_surface) call RegPack(RF, InData%Tdesc) + call RegPack(RF, InData%PlatformPosInit) call RegPack(RF, InData%CalcSteady) call RegPack(RF, InData%TrimCase) call RegPack(RF, InData%TrimTol) @@ -1772,6 +1775,7 @@ subroutine FAST_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%VTK_fps); if (RegCheckErr(RF, RoutineName)) return call FAST_UnpackVTK_SurfaceType(RF, OutData%VTK_surface) ! VTK_surface call RegUnpack(RF, OutData%Tdesc); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%PlatformPosInit); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%CalcSteady); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%TrimCase); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%TrimTol); if (RegCheckErr(RF, RoutineName)) return From f227e93ff27bbe92b2dccc667db630eabf8e1b45 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 14 Apr 2025 11:20:31 -0600 Subject: [PATCH 3/3] MD 2746: update r-test pointer --- reg_tests/r-test | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reg_tests/r-test b/reg_tests/r-test index 54a2998d96..8c26865157 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 54a2998d96908461eeb554f756bfc9da69055804 +Subproject commit 8c2686515773e29f102d7ae21d2c6acf103c35c5