Skip to content

MoorDyn: Feature/syrope for polyester #3324

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrew-platt merged 19 commits into from
Jun 4, 2026
+771 −37
Merged

MoorDyn: Feature/syrope for polyester #3324

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
a976072
Defined new variables of Syrope model
zhilongwei Apr 10, 2026
4244f33
1D interpolation from x and y arrays
zhilongwei Apr 10, 2026
fc7b738
Read Syrope settings
zhilongwei Apr 10, 2026
7fe99b9
Implementation for the Syrope model
zhilongwei Apr 10, 2026
38b3dc6
Allocate memory with stat
zhilongwei Apr 20, 2026
3acaa3c
Merge branch 'dev' of https://github.com/OpenFAST/openfast into featu…
zhilongwei Apr 21, 2026
f0549ea
Merge branch 'dev' of https://github.com/OpenFAST/openfast into featu…
zhilongwei May 1, 2026
25cca6d
Updated the error message in the SYROPE IC section
zhilongwei May 6, 2026
3c146e4
Line ID must be >=1 and used a robust way to read the ICs
zhilongwei May 6, 2026
94e001c
Removed unused varaibles WCK1 and WCK2
zhilongwei May 6, 2026
4dda4d2
Handle Syrope working-curve warnings without aborting MD_Init
zhilongwei May 6, 2026
d78fa63
New tests passed
zhilongwei May 6, 2026
aca2b2f
dl_1 denotes stretch instead of strain
zhilongwei May 7, 2026
494104c
Initialize Tmax and Tmean to zero
zhilongwei May 11, 2026
6974791
Check Syrope input validity
zhilongwei May 11, 2026
9951a30
Support relative path of the Syrope input files
zhilongwei May 11, 2026
5bf0452
Corrected error message when static stiffness is larger than dynamic …
zhilongwei May 18, 2026
b823678
Make sure all lines in the SYROPE IC section have Syrope line types
zhilongwei May 18, 2026
0315ef0
Merge remote-tracking branch 'OpenFAST/dev' into feature/syrope-for-p…
andrew-platt Jun 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
372 changes: 347 additions & 25 deletions modules/moordyn/src/MoorDyn.f90
View file
Open in desktop

Large diffs are not rendered by default.

18 changes: 13 additions & 5 deletions modules/moordyn/src/MoorDyn_IO.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -241,16 +241,18 @@ END SUBROUTINE setupBathymetry


! read in stiffness/damping coefficient or load nonlinear data file if applicable
SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, LineProp_Xs, LineProp_Ys, ErrStat3, ErrMsg3)
SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, LineProp_Xs, LineProp_Ys, &
ErrStat3, ErrMsg3, PriPath)

CHARACTER(40), INTENT(IN ) :: inputString
CHARACTER(*), INTENT(IN ) :: inputString
REAL(DbKi), INTENT(INOUT) :: LineProp_c
INTEGER(IntKi), INTENT( OUT) :: LineProp_nPoints
REAL(DbKi), INTENT( OUT) :: LineProp_Xs (MD_MaxNCoef) ! MD_MaxNCoef set in registry
REAL(DbKi), INTENT( OUT) :: LineProp_Ys (MD_MaxNCoef) ! MD_MaxNCoef set in registry

INTEGER(IntKi), INTENT( OUT) :: ErrStat3 ! Error status of the operation
CHARACTER(*), INTENT( OUT) :: ErrMsg3 ! Error message if ErrStat /= ErrID_None
CHARACTER(*), OPTIONAL, INTENT(IN) :: PriPath ! Path to the primary MoorDyn input file

INTEGER(IntKi) :: I
INTEGER(IntKi) :: UnCoef ! unit number for coefficient input file
Expand All @@ -259,6 +261,7 @@ SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, Line
INTEGER(IntKi) :: ErrStat4
CHARACTER(120) :: ErrMsg4
CHARACTER(120) :: Line2
CHARACTER(1024) :: FileName


if (SCAN(inputString, "abcdfghijklmnopqrstuvwxyzABCDFGHIJKLMNOPQRSTUVWXYZ") == 0) then ! "eE" are exluded as they're used for scientific notation!
Expand All @@ -273,9 +276,14 @@ SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, Line
LineProp_c = 0.0

! load lookup table data from file
IF ( PRESENT(PriPath) .AND. PathIsRelative( inputString ) ) THEN
FileName = TRIM(PriPath)//TRIM(inputString)
ELSE
FileName = TRIM(inputString)
END IF

CALL GetNewUnit( UnCoef )
CALL OpenFInpFile( UnCoef, TRIM(inputString), ErrStat4, ErrMsg4 ) ! add error handling?
CALL OpenFInpFile( UnCoef, TRIM(FileName), ErrStat4, ErrMsg4 ) ! add error handling?
IF (ErrStat4 == ErrID_Fatal) then
ErrStat3 = ErrStat4
ErrMsg3 = ErrMsg4
Expand Down Expand Up @@ -303,7 +311,7 @@ SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, Line

if (I < 2) then
ErrStat3 = ErrID_Fatal
ErrMsg3 = "Less than the minimum of 2 data lines found in file "//TRIM(inputString)//" (first 3 lines are headers)."
ErrMsg3 = "Less than the minimum of 2 data lines found in file "//TRIM(FileName)//" (first 3 lines are headers)."
LineProp_npoints = 0
Close (UnCoef)
RETURN
Expand All @@ -322,7 +330,7 @@ SUBROUTINE SplitByBars(instring, n, outstrings)

CHARACTER(*), INTENT(INOUT) :: instring
INTEGER(IntKi), INTENT( OUT) :: n
CHARACTER(40), INTENT(INOUT) :: outstrings(6) ! array of output strings. Up to 6 strings can be read
CHARACTER(*), INTENT(INOUT) :: outstrings(6) ! array of output strings. Up to 6 strings can be read

INTEGER :: pos1, pos2

Expand Down
251 changes: 245 additions & 6 deletions modules/moordyn/src/MoorDyn_Line.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,9 +90,9 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg)
! copy over elasticity data
Line%ElasticMod = LineProp%ElasticMod

if (Line%ElasticMod > 3) then
if (Line%ElasticMod > 4) then
ErrStat = ErrID_Fatal
ErrMsg = "Line ElasticMod > 3. This is not possible."
ErrMsg = "Line ElasticMod > 4. This is not possible."
RETURN
endif

Expand All @@ -101,6 +101,20 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg)
Line%stiffXs(I) = LineProp%stiffXs(I)
Line%stiffYs(I) = LineProp%stiffYs(I) ! note: this does not convert to E (not EA) like done in C version
END DO

! find static strain of the slow spring on the original working curve if Syrope model is used
if (Line%ElasticMod == 4) then
DO I = 1,Line%nEApoints
Line%stiffZs(I) = Line%stiffXs(I) - 1.0 / Line%vbeta * log(1.0 + Line%vbeta / Line%alphaMBL * Line%stiffYs(I))
END DO
end if
Comment thread
andrew-platt marked this conversation as resolved.

! set working curve formula
if (Line%ElasticMod == 4) then
Line%syropeWCForm = LineProp%syropeWCForm
line%syropeWCK1 = LineProp%syropeWCK1
line%syropeWCK2 = LineProp%syropeWCK2
endif

Line%nBApoints = LineProp%nBApoints
DO I = 1,Line%nBApoints
Expand Down Expand Up @@ -250,6 +264,18 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg)
RETURN
END IF

! allocate Syrope state variables if Syrope model is being used
if (Line%ElasticMod == 4) then
ALLOCATE ( Line%Tmax(N), Line%Tmean(N), STAT = ErrStat )
IF ( ErrStat /= ErrID_None ) THEN
ErrMsg = ' Error allocating Syrope arrays, Tmax and Tmean.'
!CALL CleanUp()
RETURN
END IF
Comment thread
andrew-platt marked this conversation as resolved.
Line%Tmax = 0.0_DbKi
Line%Tmean = 0.0_DbKi
end if

! allocate node force vectors
ALLOCATE ( Line%W(3, 0:N), Line%Dp(3, 0:N), Line%Dq(3, 0:N), &
Line%Ap(3, 0:N), Line%Aq(3, 0:N), Line%B(3, 0:N), Line%Bs(3, 0:N), &
Expand Down Expand Up @@ -310,7 +336,6 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg)
REAL(DbKi), ALLOCATABLE :: LNodesZ(:)
INTEGER(IntKi) :: N


N = Line%N ! for convenience

! try to calculate initial line profile using catenary routine (from FAST v.7)
Expand Down Expand Up @@ -421,17 +446,43 @@ 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.
! If using the viscoelastic model (not Syrope), 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
IF (Line%ElasticMod > 1 .and. Line%ElasticMod < 4) 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

! If using Syrope model, initialize dl_1 to the static stretch based on the working curve and mean tension
IF (Line%ElasticMod == 4) THEN
DO I = 1, N
CALL UnitVector(Line%r(:,I-1), Line%r(:,I), Line%qs(:,I), Line%lstr(I))
CALL SetupWorkingCurve(Line, Line%Tmax(I), ErrStat2, ErrMsg2)
IF (ErrStat2 /= ErrID_None) THEN
CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '')
RETURN
END IF
CALL calculate1DinterpolationXY( Line%stiffys_(1:Line%nEApoints), &
Line%stiffzs_(1:Line%nEApoints), &
Line%Tmean(I), &
Line%dl_1(I), ErrStat2, ErrMsg2 )
Line%dl_1(I) = Line%dl_1(I) * Line%l(I)
IF (ErrStat2 /= ErrID_None) THEN
CALL SetErrStat(ErrStat2, &
'The Syrope slow spring strain is not a single-valued function of mean tension for Line '// &
trim(Num2LStr(Line%IdNum))//', segment '//trim(Num2LStr(I))//'. '// &
'The static stiffness on the working curve should be smaller than the dynamic stiffness. '// &
trim(ErrMsg2), &
ErrStat, ErrMsg, '')
RETURN
END IF
END DO
END IF

! 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
Expand Down Expand Up @@ -1181,6 +1232,8 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai
CHARACTER(ErrMsgLen) :: ErrMsg2
CHARACTER(120) :: RoutineName = 'Line_GetStateDeriv'

Real(DbKi) :: dl_s ! static stretch (both slow and fast components)

ErrStat = ErrID_None
ErrMsg = ""

Expand Down Expand Up @@ -1384,7 +1437,7 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai
MagTd = Line%BA* ( Line%lstrd(I) - Line%lstr(I)*Line%ld(I)/Line%l(I) )/Line%l(I)

! viscoelastic model from https://asmedigitalcollection.asme.org/OMAE/proceedings/IOWTC2023/87578/V001T01A029/1195018
else if (Line%ElasticMod > 1) then
else if (Line%ElasticMod > 1 .and. Line%ElasticMod < 4) then

if (Line%ElasticMod == 3) then
if (Line%dl_1(I) > 0.0) then
Expand Down Expand Up @@ -1438,6 +1491,101 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai
! update state derivative for static stiffness stretch (last N entries in the line state vector if no VIV model, otherwise N entries past the 6*N-6 entries in this vector)
Xd( 6*N-6 + I) = ld_1

else if (Line%ElasticMod == 4) then

if (Line%dl_1(I) < 0.0) then
CALL SetErrStat(ErrID_Fatal,"Syrope model: Slow spring stretch cannot be less than zero",ErrStat,ErrMsg,RoutineName)
return
end if

IF (Line%alphaMBL <= 0 .OR. Line%vbeta <= 0 .OR. Line%l(I) <= 0 .OR. Line%dl_1(I) <= 0) THEN
CALL SetErrStat(ErrID_Warn,"Syrope model: Assumptions violated",ErrStat,ErrMsg,RoutineName)
if (wordy > 2) then
print *, "Line%alphaMBL", Line%alphaMBL
print *, "Line%vbeta", Line%vbeta
print *, "Line%l(I)", Line%l(I)
print *, "Line%dl_1(I)", Line%dl_1(I)
endif
ENDIF

dl = Line%lstr(I) - Line%l(I) ! delta l of this segment

if (EqualRealNos(Line%Tmax(I), Line%Tmean(I))) then
! on the original working curve
CALL calculate1DinterpolationXY( Line%stiffZs(1:Line%nEApoints), &
Line%stiffYs(1:Line%nEApoints), &
Line%dl_1(I) / Line%l(I), &
Line%Tmean(I), ErrStat2, ErrMsg2 )
IF (ErrStat2 /= ErrID_None) THEN
CALL SetErrStat(ErrStat2, &
'The Syrope mean tension is not a single-valued function of slow spring strain for Line '// &
trim(Num2LStr(Line%IdNum))//', segment '//trim(Num2LStr(I))//'. '// &
'The static stiffness on the original working curve is smaller than the dynamic stiffness. '// &
trim(ErrMsg2), &
ErrStat, ErrMsg, '')
RETURN
END IF

if (dl >= 0.0) then
MagT = Line%Tmean(I)
else
MagT = 0.0_DbKi ! cable can't "push"
endif

CALL calculate1DinterpolationXY( Line%stiffYs(1:Line%nEApoints), &
Line%stiffXs(1:Line%nEApoints), &
Line%Tmean(I), &
dl_s, ErrStat2, ErrMsg2 )
IF (ErrStat2 /= ErrID_None) THEN
CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '')
RETURN
END IF
else
! on the active working curve
CALL calculate1DinterpolationXY( Line%stiffzs_(1:Line%nEApoints), &
Line%stiffys_(1:Line%nEApoints), &
Line%dl_1(I) / Line%l(I), &
Line%Tmean(I), ErrStat2, ErrMsg2 )
IF (ErrStat2 /= ErrID_None) THEN
CALL SetErrStat(ErrStat2, &
'The Syrope slow spring strain is not a single-valued function of mean tension for Line '// &
trim(Num2LStr(Line%IdNum))//', segment '//trim(Num2LStr(I))//'. '// &
'The static stiffness on the working curve is smaller than the dynamic stiffness. '// &
trim(ErrMsg2), &
ErrStat, ErrMsg, '')
RETURN
END IF

if (dl >= 0.0) then
MagT = Line%Tmean(I)
else
MagT = 0.0_DbKi ! cable can't "push"
endif

CALL calculate1DinterpolationXY( Line%stiffys_(1:Line%nEApoints), &
Line%stiffxs_(1:Line%nEApoints), &
Line%Tmean(I), &
dl_s, ErrStat2, ErrMsg2 )
IF (ErrStat2 /= ErrID_None) THEN
CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '')
RETURN
END IF
endif

ld_1 = ((dl - dl_s * Line%l(I)) * (Line%alphaMBL + Line%vbeta * Line%Tmean(I)) + Line%BA_D * Line%lstrd(I)) / (Line%BA_D + Line%BA)
MagTd = Line%BA * ld_1 / Line%l(I)
Xd( 6*N-6 + I) = ld_1

! update Tmax and working curve if Tmean > Tmax
if (Line%Tmean(I) > Line%Tmax(I)) then
Line%Tmax(I) = Line%Tmean(I)
CALL SetupWorkingCurve(Line, Line%Tmax(I), ErrStat2, ErrMsg2)
IF (ErrStat2 /= ErrID_None) THEN
CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '')
RETURN
END IF
end if

end if


Expand Down Expand Up @@ -2071,6 +2219,97 @@ subroutine VisLinesMesh_Update(p,m,y,ErrStat,ErrMsg)
enddo
end subroutine VisLinesMesh_Update

SUBROUTINE SetupWorkingCurve(Line, Tmax, ErrStat, ErrMsg)

TYPE(MD_Line), INTENT(INOUT) :: Line
REAL(DbKi), INTENT(IN ) :: Tmax ! Maximum tension used to build the working curve
INTEGER(IntKi), INTENT( OUT) :: ErrStat
CHARACTER(*), INTENT( OUT) :: ErrMsg
CHARACTER(*), PARAMETER :: RoutineName = 'SetupWorkingCurve'

INTEGER(IntKi) :: I
INTEGER(IntKi) :: ErrStat2
CHARACTER(ErrMsgLen) :: ErrMsg2
REAL(DbKi) :: eps_max, eps_min ! maximum and minimum strains for the working curve
REAL(DbKi) :: xi ! normalized strain on the working curve
REAL(DbKi) :: denom
REAL(DbKi) :: log_arg
REAL(DbKi) :: exp_denom

ErrStat = ErrID_None
ErrMsg = ''
ErrStat2 = ErrID_None
ErrMsg2 = ''

! input validation
IF (Tmax < Line%StiffYs(1) .OR. Tmax > Line%StiffYs(Line%nEApoints)) THEN
CALL SetErrStat(ErrID_Fatal, "Tmax is outside the working curve range.", ErrStat, ErrMsg, RoutineName)
RETURN
END IF

! eps_max is computed from Tmax and the original working curve
CALL calculate1DinterpolationXY( Line%StiffYs(1:Line%nEApoints), &
Line%StiffXs(1:Line%nEApoints), &
Tmax, eps_max, ErrStat2, ErrMsg2 )
IF (ErrStat2 /= ErrID_None) THEN
CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
RETURN
END IF

! eps_min is the strain at zero tension
eps_min = Line%StiffXs(1) + Line%syropeWCK1 * (eps_max - Line%StiffXs(1))

! for linear working curve formula, if Line%syropeWCK1 > 1, treat it as slope/stiffness of the working curve
! in general, this value is much larger than 1
IF (Line%syropeWCForm == 1 .AND. Line%syropeWCK1 >= 1.0_DbKi) THEN
eps_min = eps_max - Tmax / Line%syropeWCK1
END IF

! make sure eps_min is less than eps_max and larger than zero
IF (eps_min < 0.0_DbKi .OR. eps_min >= eps_max) THEN
CALL SetErrStat(ErrID_Fatal, "Invalid working curve parameters: the zero tension strain on the working curve is out of range.", ErrStat, ErrMsg, RoutineName)
RETURN
END IF

denom = eps_max - eps_min
IF (ABS(denom) <= TINY(1.0_DbKi)) THEN
CALL SetErrStat(ErrID_Fatal, &
"Vertical Syrope working curve detected.", &
ErrStat, ErrMsg, RoutineName)
RETURN
END IF

DO I = 1, Line%nEApoints
Line%Stiffxs_(I) = eps_min + denom * REAL(I - 1, DbKi) / REAL(Line%nEApoints - 1, DbKi)
xi = (Line%Stiffxs_(I) - eps_min) / denom ! normalized strain

IF (Line%syropeWCForm == 1) THEN ! linear working curve
Line%Stiffys_(I) = Tmax * xi

ELSE IF (Line%syropeWCForm == 2) THEN ! quadratic working curve
Line%Stiffys_(I) = Tmax * xi * (Line%syropeWCK2 * xi + (1.0_DbKi - Line%syropeWCK2))

ELSE IF (Line%syropeWCForm == 3) THEN ! exponential working curve
exp_denom = 1.0_DbKi - EXP(Line%syropeWCK2)
Line%Stiffys_(I) = Tmax * (1.0_DbKi - EXP(Line%syropeWCK2 * xi)) / exp_denom

ELSE
CALL SetErrStat(ErrID_Fatal, "Invalid working curve formula specified.", ErrStat, ErrMsg, RoutineName)
RETURN
END IF

log_arg = 1.0_DbKi + (Line%vbeta / Line%alphaMBL) * Line%Stiffys_(I)
IF (log_arg <= 0.0_DbKi) THEN
CALL SetErrStat(ErrID_Fatal, &
"Non-positive argument in log() when setting working curve in Syrope model.", &
ErrStat, ErrMsg, RoutineName)
RETURN
END IF

Line%Stiffzs_(I) = Line%Stiffxs_(I) - (1.0_DbKi / Line%vbeta) * LOG(log_arg)
END DO

END SUBROUTINE SetupWorkingCurve


END MODULE MoorDyn_Line
Loading
Loading