feat(core): Implementation of the Syrope model for polyester ropes

* Syrope model implementation

* Optimization of the exponential working curve formula

* Fixed wrong initial mean tension and wrong working when Tmean < Tmax

* Documentation for Syrope extension

* Test for Syrope model with slow loading

* Three working curve formulas in one test

* Formatting

* Renamed the input file for the Syrope test

* Syrope test passed

* Added a reference for the Syrope model

* Removed bold texts in Syrope input description

* Keep the Syrope reference in Chicago style as the others

* Use (std::max) and (std::min) instead of std::max and std::min for Windows sake

* #define _USE_MATH_DEFINES before including cmath

* Disabled quadratic and exponential wc formulas, increased time step size and reduced IC duration

* formatted test syrope.cpp

* Put the references in chronological order

* Moved the nonlinear stiffness inputs to the Additional MoorDyn Files

* Read syrope working curve settings from a separate file

* Read initial Tmax and Tmean for Syrope lines in the LINES section

* Updated the test input files

* Updated the docs for Syrope model

* read Tmax and Tmean from a separate section instead from LINES section

* Throw a warning if a SYROPE IC is present but Syrope is not used

* Read Tmax and Tmean for Syrope lines from a separate section instead from LINES section

* Updated the input files for test syrope

* Updated documentation for Syrope IC inputs

* Excluded the syrope test from memcheck

* Removed Tmax0 and Tmean0 from LINES section

Author: zhilongwei

.github/workflows/memcheck.yml

@@ -36,4 +36,4 @@ jobs:
     - name: Test
       working-directory: ${{github.workspace}}/build
       # We are just testing in Linux
-      run: ctest -C ${{env.BUILD_TYPE}} -T memcheck --output-on-failure -j ${{env.PROCESSES}} -E "(seafloor|time_schemes|wavekin|wilson|lowe_and_langley_2006)"
+      run: ctest -C ${{env.BUILD_TYPE}} -T memcheck --output-on-failure -j ${{env.PROCESSES}} -E "(seafloor|time_schemes|wavekin|wilson|lowe_and_langley_2006|syrope)"

docs/inputs.rst

@@ -266,6 +266,7 @@ There is not a limit on the number of lines you can write here.
  --------------------- MoorDyn Input File ------------------------------------
  MoorDyn v2 sample input file
 
+.. _line-types-v2:
 Line Types
 ^^^^^^^^^^
 
@@ -285,7 +286,7 @@ The columns in order are as follows:
  - Diam –  the volume-equivalent diameter of the line – the diameter of a cylinder having the same 
    displacement per unit length (m)
  - MassDen –  the mass per unit length of the line (kg/m)
- - EA – the line stiffness, product of elasticity modulus and cross-sectional area (N)
+ - EA – the line stiffness, product of elasticity modulus and cross-sectional area (N), or nonlinear parameters (see descriptions below)
  - BA/-zeta –  the line internal damping (measured in N-s) or, if a negative value is entered, the 
    desired damping ratio (in fraction of critical) for the line type (and MoorDyn will set the BA 
    of each line accordingly)
@@ -306,24 +307,13 @@ The columns in order are as follows:
  - cF - OPTIONAL - the center of the range of non-dimensional frequencies for the CF VIV synchronization model. If it is not
    provided and VIV is enabled (Cl > 0) then it is default to 0.18 to align with the the theory found :ref:`here <version2>`
 
-Note: Non-linear values for the stiffness (EA) are an option in MoorDyn. For this, a file name can be provided instead of a number. This file 
-must be located in the same folder as the main MoorDyn input file for MoorDyn-C or for MoorDyn-F 
-in the same folder as the executable calling MoorDyn-F, unless a path is specified. Such file is a 
-tabulated file with 3 header lines and then a strain column and a tension column separated by a blank space:
-
-.. code-block:: none
-
-  ----Polyester----
-  Strain    Tension
-  (-)       (N)
-  0.0       0.0
-  ...       ...
-
 Note: MoorDyn has the ability to model the viscoelastic properties of synthetic lines in two ways. The first method, from work linked in the 
 :ref:`theory section <theory>`, allows a user to specify a bar-separated constant dynamic and static stiffness. The second method allows the user 
 to provide a constant static stiffness and two terms to determine the dynamic stiffness as a linear function of mean load. The equation is:
-`EA_d = EA_Dc + EA_D_Lm * mean_load` where `EA_D_Lm` is the slope of the load-stiffness curve. Both of these methods allow users to provide static 
-and dynamic damping coefficients as values separated by |. While the static damping can be described as a fraction of critical, the dynamic damping 
+`EA_d = EA_Dc + EA_D_Lm * mean_load` where `EA_D_Lm` is the slope of the load-stiffness curve. This method can be used standalone, 
+and is also part of the Syrope model described in :ref:`Syrope Model for Polyester Lines <syrope-model-polyester-lines>`. 
+Both of these methods allow users to provide static and dynamic damping coefficients as values separated by |. 
+While the static damping can be described as a fraction of critical, the dynamic damping 
 needs to be input as a value. Example inputs are below: 
 
 .. code-block:: none
@@ -497,9 +487,9 @@ The columns are as follows:
  - UnstrLen - the unstretched length of the line
  - NumSegs - how many segments the line is discretized into (it will have NumSegs+1 nodes total, 
    including its two end nodes)
- - LineOutputs - any data to be output in a dedicated output file for that line. 
+ - LineOutputs - any data to be output in a dedicated output file for that line.
 
-This last entry expects a string of one or more characters without spaces, each character 
+This `LineOutputs` entry expects a string of one or more characters without spaces, each character 
 activating a given output property. A placeholder character such as “-” should be used if no 
 outputs are wanted. Ten output properties are currently possible:
 
@@ -1008,6 +998,80 @@ data.
  5000            0.15         0.0
  --------------------- need this line ------------------
 
+Nonlinear Stiffness (EA) Inputs
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Non-linear values for the stiffness (EA) are an option in MoorDyn. For this, a file name can be provided instead of a number. This file 
+must be located in the same folder as the main MoorDyn input file for MoorDyn-C or for MoorDyn-F 
+in the same folder as the executable calling MoorDyn-F, unless a path is specified. Such file is a 
+tabulated file with 3 header lines and then a strain column and a tension column separated by a blank space:
+
+.. code-block:: none
+
+  ----Polyester----
+  Strain    Tension
+  (-)       (N)
+  0.0       0.0
+  ...       ...
+
+.. _syrope-model-polyester-lines:
+Syrope Model for Polyester Lines
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+MoorDyn supports the Syrope material model for polyester ropes. In this formulation, the
+tension–strain response is nonlinear and load-path dependent, meaning the unloading and
+reloading curves differ. See additional details in the :ref:`theory section <theory>`
+and the references therein.
+
+To enable the Syrope model for a given line type, specify the EA field using three
+bar-separated entries, analogous to the load-dependent dynamic stiffness method described in :ref:`Line Types <line-types-v2>`.
+To indicate Syrope, the first entry must start with ``SYROPE:``, followed by the name (or path)
+of a Syrope config file defining the working curves. The second and third entries are
+``EA_Dc`` and ``EA_D_Lm`` (as in the load-dependent dynamic stiffness method in :ref:`Line Types <line-types-v2>`). Static and dynamic
+damping may be provided in BA as ``BA_s|BA_d`` (same convention as above).
+
+Example:
+
+.. code-block:: none
+
+  TypeName   Diam    Mass/m     EA                                 BA
+  (name)     (m)     (kg/m)     (N)                               (N-s)
+  poly       ...     ...        SYROPE:syrope.dat|EA_Dc|EA_D_Lm   BA_s|BA_d
+
+The file ``syrope.dat`` defines 
+
+- OWC: a tabulated nonlinear original working curve
+- WCType: the working-curve type (options: ``LINEAR``, ``QUADRATIC``, ``EXP``)
+- k1: the first coefficient for the working curve formula (sets the strain offset where mean tension is zero)
+- k2: the second coefficient for the working curve formula (sets the curve shape)
+
+One example of such file is shown below.
+
+.. code-block:: none
+
+  ../owc.dat  OWC     Original working curve lookup table path, relative to the MoorDyn input file
+  LINEAR      WCType  Working curve formula: LINEAR, QUADRATIC, or EXP
+  0.25        k1      First parameter defining the working curve shape
+  1.00        k2      Second parameter defining the working curve shape
+
+The ``owc.dat`` file follows the same format as the tabulated non-linear stiffness files described
+above (three header lines, then columns for strain and tension).
+
+Initial conditions are specified for the Syrope lines in a separate ``SYROPE IC`` section placed after the ``LINES`` section.
+
+- Line(s) - the Syrope line number(s). This may be a single line number or a comma-separated list of line numbers. If non-Syrope lines are included, they are ignored and a warning is raised.
+- Tmax0 – the initial highest mean tension experienced by the line prior to the current time (N)
+- Tmean0 – the initial mean tension (N)
+
+One example of the ``SYROPE IC`` section is shown below.
+
+.. code-block:: none
+
+  ---------------------- SYROPE IC --------------------------------------
+  Line(s)   Tmax0    Tmean0
+  (-)        (N)       (N)
+  7,8,9    3.53e6    1.18e6
+
+In general, ``Tmax0 >= Tmean0``.
+
 MoorDyn-F with FAST.Farm - Inputs
 -------------------------------
 

docs/references.rst

@@ -82,16 +82,16 @@ Overview of MoorDyn v2 (bodies, rods, and line failures):
   `Hall, Matthew, “MoorDyn V2: New Capabilities in Mooring System Components and Load Cases.” In Proceedings of the ASME 2020 39th International 
   Conference on Ocean, Offshore and Arctic Engineering. virtual conference, 2020. <https://www.nrel.gov/docs/fy20osti/76555.pdf>`_
 
-Seabed friction and bathymetry approach used in v2:
-
-  `Housner, Stein, Ericka Lozon, Bruce Martin, Dorian Brefort, and Matthew Hall, “Seabed Bathymetry and Friction Modeling in MoorDyn.” Journal of 
-  Physics: Conference Series 2362, no. 1, Nov 2022: 012018. <https://doi.org/10.1088/1742-6596/2362/1/012018>`_
-
 Implementation of bending stiffness modeling for power cables:
 
   `Hall, Matthew, Senu Sirnivas, and Yi-Hsiang Yu, “Implementation and Verification of Cable Bending Stiffness in MoorDyn.” In ASME 2021 3rd International Offshore Wind 
   Technical Conference, V001T01A011. Virtual, Online: American Society of Mechanical Engineers, 2021. <https://doi.org/10.1115/IOWTC2021-3565>`_
 
+Seabed friction and bathymetry approach used in v2:
+
+  `Housner, Stein, Ericka Lozon, Bruce Martin, Dorian Brefort, and Matthew Hall, “Seabed Bathymetry and Friction Modeling in MoorDyn.” Journal of 
+  Physics: Conference Series 2362, no. 1, Nov 2022: 012018. <https://doi.org/10.1088/1742-6596/2362/1/012018>`_
+
 Non-linear line stiffness:
 
  `Lozon, Ericka, Matthew Hall, Paul McEvoy, Seojin Kim, and Bradley Ling, “Design and Analysis of a Floating-Wind Shallow-Water Mooring System 
@@ -122,6 +122,10 @@ Modeling of Bi-stable Nonlinear Energy Sinks in MoorDyn (most recent description
   `Anargyros Michaloliakos, Wei-Ying Wong, Ryan Davies, Malakonda Reddy Lekkala, Matthew Hall, Lei Zuo, Alexander F. Vakakis, "Stabilizing dynamic subsea power cables using 
   Bi-stable nonlinear energy sinks", Ocean Engineering, vol. 334, August 2025. <https://doi.org/10.1016/j.oceaneng.2025.121613>`_
 
+Syrope model for polyester ropes:
+
+  `Wei, Zhilong, Harry B. Bingham, and Yanlin Shao. 2026. “ESOMOOR D5.1: Extended Moordyn Solver and Validation Report”. Technical University of Denmark. <https://doi.org/10.11583/DTU.31408806>`_
+
 The Fortran version of MoorDyn is available as a module inside of OpenFAST:
 
   https://openfast.readthedocs.io/en/main/

source/Line.cpp

@@ -118,6 +118,50 @@ Line::getNonlinearBA(real ld_stretched, real l_unstretched) const
 	return Yi / Xi;
 }
 
+void
+Line::setWorkingCurve(real Tmax)
+{
+	// Check if Tmax is within the valid range
+	if (Tmax < stiffYs.front() || Tmax > stiffYs.back()) {
+		LOGERR << "Line " << lineId << ": Tmax value " << Tmax
+		       << " is outside the working curve range [" << stiffYs.front()
+		       << ", " << stiffYs.back() << "]." << endl;
+		throw moordyn::invalid_value_error("Invalid Tmax for working curve");
+	}
+	// maximum and minimum strains for the working curve
+	real eps_max = interp(stiffYs, stiffXs, Tmax);
+	real eps_min = stiffXs.front() + k1 * (eps_max - stiffXs.front());
+
+	// For linear working curve, if k1 >= 1.0 (usually much larger than unity),
+	// it is taken as the slope (or dimensional stiffness)
+	if (syropeWCForm == SYROPE_LINEAR_WC && k1 >= 1.0)
+		eps_min = eps_max - Tmax / k1;
+
+	for (unsigned int I = 0; I < nEApointsWC; I++) {
+		// Normalize the strain points between eps_min and eps_max
+		stiffxs_[I] = eps_min + (eps_max - eps_min) * I /
+		                            (nEApointsWC - 1); // total strain
+		real xi = (stiffxs_[I] - eps_min) / (eps_max - eps_min);
+
+		if (syropeWCForm == SYROPE_LINEAR_WC)
+			stiffys_[I] = Tmax * xi;
+		else if (syropeWCForm == SYROPE_QUADRATIC_WC)
+			stiffys_[I] = Tmax * xi * (k2 * xi + (1.0 - k2));
+		else if (syropeWCForm == SYROPE_EXP_WC)
+			stiffys_[I] = Tmax * (1.0 - exp(k2 * xi)) / (1.0 - exp(k2));
+		else {
+			LOGERR << "Line " << lineId
+			       << ": Invalid Syrope working curve formula " << syropeWCForm
+			       << "." << endl;
+			throw moordyn::invalid_value_error(
+			    "Invalid Syrope working curve formula");
+		}
+
+		stiffzs_[I] = stiffxs_[I] -
+		              1.0 / vbeta * log(1.0 + vbeta / alphaMBL * stiffys_[I]);
+	}
+}
+
 void
 Line::setup(int number_in,
             LineProps* props,
@@ -158,16 +202,50 @@ Line::setup(int number_in,
 	Cl = props->Cl;
 	dF = props->dF;
 	cF = props->cF;
+	syropeWCForm = (syrope_wc_formula)props->SyropeWCForm;
+	k1 = props->k1;
+	k2 = props->k2;
 
 	// copy in nonlinear stress-strain data if applicable
 	stiffXs.clear();
 	stiffYs.clear();
+	stiffZs.clear();
 	nEApoints = props->nEApoints;
 	for (unsigned int I = 0; I < nEApoints; I++) {
 		stiffXs.push_back(props->stiffXs[I]);
 		stiffYs.push_back(props->stiffYs[I]);
 	}
 
+	// Compute the slow strain
+	if (ElasticMod == ELASTIC_SYROPE) {
+		for (unsigned int I = 0; I < nEApoints; I++) {
+			stiffZs.push_back(
+			    props->stiffXs[I] -
+			    1.0 / vbeta * log(1.0 + vbeta / alphaMBL * props->stiffYs[I]));
+		}
+
+		// Check if stiffZs is strictly increasing
+		for (unsigned int I = 1; I < nEApoints; I++) {
+			if (stiffZs[I] <= stiffZs[I - 1]) {
+				LOGERR << "Line " << number
+				       << ": The slow strain values for the Syrope model are "
+				          "not strictly increasing. Please check the "
+				          "stress-strain data provided."
+				       << endl;
+				throw moordyn::invalid_value_error(
+				    "Invalid slow strain data for Syrope model");
+			}
+		}
+
+		// Initialize working curve arrays
+		stiffxs_.clear();
+		stiffys_.clear();
+		stiffzs_.clear();
+		stiffxs_.assign(nEApointsWC, 0.0);
+		stiffys_.assign(nEApointsWC, 0.0);
+		stiffzs_.assign(nEApointsWC, 0.0);
+	}
+
 	// Use the last entry on the lookup table. see Line::initialize()
 	const real EA = nEApoints ? stiffYs.back() / stiffXs.back() : props->EA;
 	NatFreqCFL::length(UnstrLen / N);
@@ -606,13 +684,24 @@ Line::initialize()
 	// 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 (ElasticMod != ELASTIC_CONSTANT) {
+	if (ElasticMod != ELASTIC_CONSTANT && ElasticMod != ELASTIC_SYROPE) {
 		for (unsigned int i = 0; i < N; i++) {
 			lstr[i] = unitvector(qs[i], r[i], r[i + 1]);
 			dl_1[i] = lstr[i] - l[i]; // delta l of the segment
 		}
 	}
 
+	// If using the Syrope model, initalize the deltaL_1 to slow stretch based
+	// on the mean tension and the active working curve
+	if (ElasticMod == ELASTIC_SYROPE) {
+		for (unsigned int i = 0; i < N; i++) {
+			lstr[i] = unitvector(qs[i], r[i], r[i + 1]);
+			setWorkingCurve(Tmax[i]);
+			dl_1[i] = interp(stiffys_, stiffzs_, Tmean[i]) *
+			          l[i]; // stretch instead of strain
+		}
+	}
+
 	// If using the VIV model, initalize as a distribution between 0 and 2pi for
 	// inital phase of lift force to avoid initial transient
 	if (Cl > 0)
@@ -1019,7 +1108,8 @@ Line::getStateDeriv(InstanceStateVarView drdt)
 				BA = getNonlinearBA(ldstr[i], l[i]);
 			Td[i] = BA * ldstr[i] / l[i] * qs[i];
 
-		} else {
+		} else if (ElasticMod == ELASTIC_VISCO_CTE ||
+		           ElasticMod == ELASTIC_VISCO_MEAN) {
 			// viscoelastic model from
 			// https://asmedigitalcollection.asme.org/OMAE/proceedings/IOWTC2023/87578/V001T01A029/1195018
 			// note that dl_1[i] is the same as Line%dl_1 in MD-F. This is the
@@ -1094,6 +1184,42 @@ Line::getStateDeriv(InstanceStateVarView drdt)
 			// update state derivative for static stiffness stretch. The visco
 			// state is always the last element in the row.
 			drdt.row(i).tail<1>()[0] = ld_1[i];
+		} else if (ElasticMod == ELASTIC_SYROPE) {
+			const real dl = lstr[i] - l[i];
+
+			if (moordyn::EqualRealNos(Tmean[i], Tmax[i])) {
+				// on the original working curve
+				Tmean[i] = interp(stiffZs, stiffYs, dl_1[i] / l[i]);
+				if (dl >= 0.0)
+					T[i] = Tmean[i] * qs[i];
+				else
+					T[i] = vec::Zero();
+
+				ld_1[i] = ((dl - interp(stiffYs, stiffXs, Tmean[i]) * l[i]) *
+				               (alphaMBL + vbeta * Tmean[i]) +
+				           BA_D * ldstr[i]) /
+				          (BA + BA_D);
+				drdt.row(i).tail<1>()[0] = ld_1[i];
+				Td[i] = ld_1[i] / l[i] * BA * qs[i];
+			} else { // on the working curve
+				Tmean[i] = interp(stiffzs_, stiffys_, dl_1[i] / l[i]);
+				if (dl >= 0.0)
+					T[i] = Tmean[i] * qs[i];
+				else
+					T[i] = vec::Zero();
+
+				ld_1[i] = ((dl - interp(stiffys_, stiffxs_, Tmean[i]) * l[i]) *
+				               (alphaMBL + vbeta * Tmean[i]) +
+				           BA_D * ldstr[i]) /
+				          (BA + BA_D);
+				drdt.row(i).tail<1>()[0] = ld_1[i];
+				Td[i] = ld_1[i] / l[i] * BA * qs[i];
+			}
+
+			if (Tmean[i] > Tmax[i]) {
+				Tmax[i] = Tmean[i];
+				setWorkingCurve(Tmax[i]);
+			}
 		}
 	}