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 13 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
13 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
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
334 changes: 318 additions & 16 deletions modules/moordyn/src/MoorDyn.f90
View file
Open in desktop

Large diffs are not rendered by default.

219 changes: 213 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

! 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,16 @@ 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

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 +325,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 +435,38 @@ 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, 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 +1216,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 +1412,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 +1466,76 @@ 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

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, 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, 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 +2169,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
70 changes: 70 additions & 0 deletions modules/moordyn/src/MoorDyn_Misc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@ MODULE MoorDyn_Misc
PUBLIC :: calculate4Dinterpolation
PUBLIC :: calculate3Dinterpolation
PUBLIC :: calculate2Dinterpolation
PUBLIC :: calculate1DinterpolationXY

PUBLIC :: getDepthFromBathymetry

Expand Down Expand Up @@ -824,6 +825,75 @@ SUBROUTINE calculate1Dinterpolation(f, ix0, fx, c)
c = c0*(1.0-fx) + c1*fx
END SUBROUTINE calculate1Dinterpolation

SUBROUTINE calculate1DinterpolationXY(x, y, xi, yi, ErrStat, ErrMsg)
REAL(DbKi), INTENT(IN) :: x(:)
REAL(DbKi), INTENT(IN) :: y(:)
REAL(DbKi), INTENT(IN) :: xi
REAL(DbKi), INTENT(OUT) :: yi
INTEGER(IntKi), INTENT(OUT) :: ErrStat
CHARACTER(*), INTENT(OUT) :: ErrMsg
CHARACTER(*), PARAMETER :: RoutineName = 'calculate1DinterpolationXY'

INTEGER(IntKi) :: n
INTEGER(IntKi) :: lo, hi, mid
INTEGER(IntKi) :: ix0
REAL(DbKi) :: fx, dx

ErrStat = ErrID_None
ErrMsg = ''
yi = 0.0_DbKi

n = SIZE(x)

IF (n /= SIZE(y)) THEN
CALL SetErrStat(ErrID_Fatal, 'x and y must have the same size.', ErrStat, ErrMsg, RoutineName)
RETURN
END IF

IF (n < 2) THEN
CALL SetErrStat(ErrID_Fatal, 'At least two points are required for interpolation.', ErrStat, ErrMsg, RoutineName)
RETURN
END IF

IF (ANY(x(2:n) <= x(1:n-1))) THEN
CALL SetErrStat(ErrID_Fatal, 'x array must be strictly increasing.', ErrStat, ErrMsg, RoutineName)
RETURN
END IF

! Clamp to endpoints
IF (xi <= x(1)) THEN
yi = y(1)
RETURN
ELSE IF (xi >= x(n)) THEN
yi = y(n)
RETURN
END IF

! Binary search for interval
lo = 1
hi = n
DO WHILE (hi - lo > 1)
mid = lo + (hi - lo) / 2
IF (x(mid) <= xi) THEN
lo = mid
ELSE
hi = mid
END IF
END DO

ix0 = lo
dx = x(ix0+1) - x(ix0)

IF (ABS(dx) <= TINY(1.0_DbKi)) THEN
CALL SetErrStat(ErrID_Fatal, 'Zero or too-small x interval encountered.', ErrStat, ErrMsg, RoutineName)
RETURN
END IF

fx = (xi - x(ix0)) / dx

CALL calculate1Dinterpolation(y, ix0, fx, yi)

END SUBROUTINE calculate1DinterpolationXY



Expand Down
Loading