MoorDyn bugfixes

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 
 

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

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
 

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"            -

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

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]"	-

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

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