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

Open
zhilongwei wants to merge 16 commits into
base: dev
Choose a base branch
from
Open

MoorDyn: Feature/syrope for polyester #3324

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 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
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
366 changes: 341 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 on lines +105 to +110
Copy link
Copy Markdown
Contributor

@RyanDavies19 RyanDavies19 May 10, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@zhilongwei if you could add in these simple checks that would be great. We have a similar check for the Line%ElasticMod==3 case in Line_GetStateDeriv:

openfast/modules/moordyn/src/MoorDyn_Line.f90

Line 1377 in 2895884

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

Copy link
Copy Markdown
Author

@zhilongwei zhilongwei May 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have done so.

https://github.com/zhilongwei/openfast/blob/6974791cb5f136be4507413026b59d6818b9c93e/modules/moordyn/src/MoorDyn_Line.f90#L1476-L1491


! 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 @@ -241,6 +255,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
Copy link
Copy Markdown
Author

@zhilongwei zhilongwei May 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Line%Tmax and Line%Tmean are already initialized when reading the SYROPE IC section in MoorDyn.f90.

                  DO il = 1, N
                     if (TempIDnums(il) >= 1 .and. TempIDnums(il) <= p%nLines) then      ! ensure line ID is in range
                        DO j = 1, m%LineList(TempIDnums(il))%N
                           m%LineList(TempIDnums(il))%Tmax(j) = Tmax0
                           m%LineList(TempIDnums(il))%Tmean(j) = Tmean0
                        END DO
                     else
                        CALL SetErrStat( ErrID_Fatal, ' Unable to parse Syrope line initial conditions '//trim(Num2LStr(l))//'. Line number '//TRIM(Int2LStr(TempIDnums(il)))//' out of bounds.', ErrStat, ErrMsg, RoutineName )
                        CALL CleanUp()
                        return
                     endif 
                  END DO

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 @@ -301,7 +327,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 @@ -412,17 +437,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 is 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 @@ -1172,6 +1223,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 @@ -1366,7 +1419,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 @@ -1420,6 +1473,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 @@ -2053,6 +2201,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