diff --git a/docs/inputs.rst b/docs/inputs.rst index 4ac84b86..b31d4899 100644 --- a/docs/inputs.rst +++ b/docs/inputs.rst @@ -333,6 +333,65 @@ needs to be input as a value. Example inputs are below: Polyester ... ... EA_s|EA_d BA_s|BA_d <-- Constant dynamic stiffness method with static and dynamic damping Polyester ... ... EA_s|EA_Dc|EA_D_Lm BA_s|BA_d <-- Load dependent dynamic stiffness method with static and dynamic damping +MoorDyn also 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 (hysteresis). See additional details in the :ref:`theory section ` +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 above. +To indicate Syrope, the first entry must start with ``SYROPE:``, followed by the name (or path) +of a tabulated file defining the original working curve. The second and third entries are +``EA_Dc`` and ``EA_D_Lm`` (as in the load-dependent dynamic stiffness method). 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:owc.txt|EA_Dc|EA_D_Lm BA_s|BA_d + +The file ``owc.txt`` follows the same format as the tabulated non-linear stiffness files described +above (three header lines, then columns for strain and tension). + +In addition to the line-type inputs above, users must also provide (i) the working-curve +parameterization and (ii) Syrope initial conditions for each Syrope line. + +The working-curve parameters are provided in a dedicated table. The columns are: + +- LineType – a string matching a Line Dictionary entry that uses the Syrope model +- WCType – the working-curve type (options: ``LINEAR``, ``QUADRATIC``, ``EXP``) +- k1 – the primary coefficient of the working curve (sets the strain offset where mean tension is zero) +- k2 – the secondary coefficient of the working curve (sets the curve shape) + +Example: + +.. code-block:: none + + -------------------------- SYROPE WORKING CURVES ------------------------------ + LineType WCType k1 k2 + (-) (-) (-) (-) + poly LINEAR k1 k2 + +Initial conditions are specified per line. The Line ID must match a line defined in the +*Lines* section and must reference a line type that uses the Syrope model. + +- Tmax – the highest mean tension experienced by the line prior to the current time (N) +- Tmean – the current mean tension at initialization (N) + +In general, ``Tmax0 >= Tmean0``. + +Example: + +.. code-block:: none + + -------------------------- SYROPE IC ------------------------------------------ + Line Tmax Tmean + (-) (N) (N) + 1 Tmax0 Tmean0 + Rod Types ^^^^^^^^^ diff --git a/docs/references.rst b/docs/references.rst index a6e24059..b531c35d 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -107,6 +107,10 @@ Viscoelastic approach for non-linear rope behavior: `Hall, Matthew, Brian Duong, and Ericka Lozon, “Streamlined Loads Analysis of Floating Wind Turbines With Fiber Rope Mooring Lines.” In ASME 2023 5th International Offshore Wind Technical Conference, V001T01A029. Exeter, UK: American Society of Mechanical Engineers, 2023. `_ +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. `_ + Updated MoorDyn-OpenFOAM Coupling: `Haifei Chen, Tanausú Almeida Medina, and Jose Luis Cercos-Pita, "CFD simulation of multiple moored floating structures using OpenFOAM: An open-access mooring restraints diff --git a/source/Line.cpp b/source/Line.cpp index f14b743c..27b424b4 100644 --- a/source/Line.cpp +++ b/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]); + } } } diff --git a/source/Line.hpp b/source/Line.hpp index 9ee1b5f5..4691374c 100644 --- a/source/Line.hpp +++ b/source/Line.hpp @@ -90,8 +90,20 @@ class DECLDIR Line final ELASTIC_VISCO_CTE = 2, /// mean load dependent dynamic stiffness ELASTIC_VISCO_MEAN = 3, + /// Syrope model + ELASTIC_SYROPE = 4, } elastic_model; + /// @brief Syrope working curve formulas + typedef enum + { + /// linear working curve + SYROPE_LINEAR_WC = 1, + /// quadratic working curve + SYROPE_QUADRATIC_WC = 2, + /// exponential working curve + SYROPE_EXP_WC = 3, + } syrope_wc_formula; /** @brief Get the non-linear stiffness. This is interpolated from a * curve provided in the input file. @@ -130,6 +142,12 @@ class DECLDIR Line final { return seafloor ? seafloor->getAverageDepth() : -env->WtrDpth; } + + /** @brief Set up the active working curve with given Tmax + * @param Tmax preceding highest mean tension + */ + inline void setWorkingCurve(real Tmax); + // ENVIRONMENTAL STUFF /// Global struct that holds environmental settings EnvCondRef env; @@ -218,6 +236,9 @@ class DECLDIR Line final std::vector stiffXs; /// y array for stress-strain lookup table std::vector stiffYs; + /// z array for stress-strain original working curves (used only in Syrope + /// model) + std::vector stiffZs; /// number of values in bent stiffness lookup table (0 for constant EI) unsigned int nEIpoints; /// x array for bent stiffness lookup table @@ -231,6 +252,13 @@ class DECLDIR Line final /// y array for stress-strainrate lookup table std::vector dampYs; + /// auxiliary x array for Syrope working curve + std::vector stiffxs_; + /// auxiliary y array for Syrope working curve + std::vector stiffys_; + /// auxiliary z array for Syrope working curve + std::vector stiffzs_; + // Externally provided data /** true if pressure bending forces shall be considered, false otherwise * @see @@ -265,6 +293,20 @@ class DECLDIR Line final // line segment volumes std::vector V; + // Syrope variables + /// Working curve formula. See ::syrope_wc_formula + syrope_wc_formula syropeWCForm; + /// first coefficient for the Syrope working curve formula + moordyn::real k1; + /// second coefficient for the Syrope working curve formula + moordyn::real k2; + /// preceding highest mean tensions + std::vector Tmax; + /// mean tensions + std::vector Tmean; + /// number of interpolation points for Syrope working curve + const unsigned int nEApointsWC = 30; + // forces /// segment tensions std::vector T; @@ -487,6 +529,24 @@ class DECLDIR Line final * @} */ + /** @brief Set the initial Tmax values + * @param Tmax0 Initial Tmax value + */ + inline void setInitialTmax(moordyn::real Tmax0) { Tmax.assign(N, Tmax0); } + + /** @brief Set the initial mean tension values + * @param Tmean0 Initial Tmean value + */ + inline void setInitialTmean(moordyn::real Tmean0) + { + Tmean.assign(N, Tmean0); + } + + /** @brief Get the elastic model used by the line + * @return The elastic model. See ::elastic_model + */ + inline elastic_model getElasticModel() const { return ElasticMod; } + /** @brief Number of segments * * The number of nodes can be computed as moordyn::Line::getN() + 1 diff --git a/source/Misc.hpp b/source/Misc.hpp index 45c75e7f..d93013fd 100644 --- a/source/Misc.hpp +++ b/source/Misc.hpp @@ -1178,6 +1178,10 @@ typedef struct _LineProps // (matching Line Dictionary inputs) double bstiffXs[nCoef]; // x array for stress-strain lookup table (up to nCoef) double bstiffYs[nCoef]; // y array for stress-strain lookup table + int SyropeWCForm; // Syrope working curve formula (1=linear, 2=quadratic, + // 3=exponential) + double k1; // first coefficient for the Syrope working curve formula + double k2; // second coefficient for the Syrope working curve formula } LineProps; typedef struct _RodProps // (matching Rod Dictionary inputs) diff --git a/source/MoorDyn2.cpp b/source/MoorDyn2.cpp index 5bd39417..c930cdab 100644 --- a/source/MoorDyn2.cpp +++ b/source/MoorDyn2.cpp @@ -996,6 +996,63 @@ moordyn::MoorDyn::ReadInFile() } } + // If Syrope working curve is defined, read it + if ((i = findStartOfSection(in_txt, { "SYROPE WORKING CURVES" })) != -1) { + LOGDBG << " Reading Syrope working curve formulas:" << endl; + + // parse until the next header or the end of the file + while ((in_txt[i].find("---") == string::npos) && + (i < (int)in_txt.size())) { + vector entries = moordyn::str::split(in_txt[i], ' '); + // Expecting 4 entries: LineType (string), WCType (string), k1 + // (float), and k2 (float) + if (entries.size() != 4) { + LOGERR << "Error in " << _filepath << " at line " << i + 1 + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "4 fields are required, but just " << entries.size() + << " are provided" << endl; + return MOORDYN_INVALID_INPUT; + } + + string line_prop_type = entries[0]; + string wc_formula = entries[1]; + double k1 = atof(entries[2].c_str()); + double k2 = atof(entries[3].c_str()); + + // Loop all line properties to find the matching one + bool line_type_found = false; + for (auto& line_prop : LinePropList) { + if (line_prop->type == line_prop_type) { + line_type_found = true; + + if (wc_formula == "LINEAR") + line_prop->SyropeWCForm = 1; + else if (wc_formula == "QUADRATIC") + line_prop->SyropeWCForm = 2; + else if (wc_formula == "EXP") + line_prop->SyropeWCForm = 3; + else { + LOGERR + << "Error in " << _filepath << " at line " << i + 1 + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "Invalid working curve formula: " << wc_formula + << endl; + return MOORDYN_INVALID_INPUT; + } + line_prop->k1 = k1; + line_prop->k2 = k2; + } + } + + LOGDBG << " Line Type: " << line_prop_type + << ", WC Type: " << wc_formula << ", k1: " << k1 + << ", k2: " << k2 << endl; + i++; + } + } + if ((i = findStartOfSection(in_txt, { "ROD DICTIONARY", "ROD TYPES" })) != -1) { LOGDBG << " Reading rod types:" << endl; @@ -1332,6 +1389,44 @@ moordyn::MoorDyn::ReadInFile() } } + if ((i = findStartOfSection(in_txt, { "SYROPE IC" })) != -1) { + LOGDBG << " Reading Syrope line initial conditions:" << endl; + + // parse until the next header or the end of the file + while ((in_txt[i].find("---") == string::npos) && + (i < (int)in_txt.size())) { + vector entries = moordyn::str::split(in_txt[i], ' '); + if (entries.size() < 3) { + LOGERR << "Error in " << _filepath << " at line " << i + 1 + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "3 fields are required, but just " << entries.size() + << " are provided" << endl; + return MOORDYN_INVALID_INPUT; + } + + int line_id = atoi(entries[0].c_str()); + real Tmax0 = atof(entries[1].c_str()); + real Tmean0 = atof(entries[2].c_str()); + + if (line_id < 1 || line_id > (int)LineList.size()) { + LOGERR << "Error in " << _filepath << " at line " << i + 1 + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "There is no line with ID " << line_id << "!" << endl; + return MOORDYN_INVALID_INPUT; + } + + auto* line = LineList[line_id - 1]; + if (line->getElasticModel() == 4) { + line->setInitialTmax(Tmax0); + line->setInitialTmean(Tmean0); + } + + i++; + } + } + if ((i = findStartOfSection(in_txt, { "FAILURE" })) != -1) { LOGDBG << " Reading failure conditions:" << endl; // parse until the next header or the end of the file @@ -1922,29 +2017,62 @@ moordyn::MoorDyn::readLineProps(string inputText, int lineNum) moordyn::error_id err; vector EA_stuff = moordyn::str::split(entries[3], '|'); const int EA_N = EA_stuff.size(); - if (EA_N == 1) { - obj->ElasticMod = 1; // normal case - } else if (EA_N == 2) { - obj->ElasticMod = 2; // viscoelastic model, constant dynamic stiffness - obj->EA_D = atof(EA_stuff[1].c_str()); - } else if (EA_N == 3) { - obj->ElasticMod = - 3; // viscoelastic model load dependent dynamic stiffness + + // Check if Syrope model is being used + const std::string& ea0 = EA_stuff.empty() ? std::string() : EA_stuff[0]; + const bool is_syrope = (!ea0.empty() && ea0.rfind("SYROPE:", 0) == 0); + + if (is_syrope) { + obj->ElasticMod = 4; // Syrope model + LOGDBG + << "Line type '" << obj->type + << "' is using the Syrope model for axial tension-strain behavior." + << endl; + + if (EA_N != 3) { + LOGERR << "A line type EA entry with Syrope model must have 3 " + "(bar-separated) values." + << endl; + return nullptr; + } + + const std::string owc = ea0.substr(7); obj->alphaMBL = atof(EA_stuff[1].c_str()); obj->vbeta = atof(EA_stuff[2].c_str()); - } else { - LOGERR - << "A line type EA entry can have at most 3 (bar-separated) values." - << endl; - return nullptr; + + err = read_curve(owc.c_str(), + &(obj->EA), + &(obj->nEApoints), + obj->stiffXs, + obj->stiffYs); + if (err) + return nullptr; + } else { // Not Syrope, other models + if (EA_N == 1) { + obj->ElasticMod = 1; // normal case + } else if (EA_N == 2) { + obj->ElasticMod = + 2; // viscoelastic model, constant dynamic stiffness + obj->EA_D = atof(EA_stuff[1].c_str()); + } else if (EA_N == 3) { + obj->ElasticMod = + 3; // viscoelastic model load dependent dynamic stiffness + obj->alphaMBL = atof(EA_stuff[1].c_str()); + obj->vbeta = atof(EA_stuff[2].c_str()); + } else { + LOGERR << "A line type EA entry can have at most 3 (bar-separated) " + "values." + << endl; + return nullptr; + } + err = read_curve(EA_stuff[0].c_str(), + &(obj->EA), + &(obj->nEApoints), + obj->stiffXs, + obj->stiffYs); + if (err) + return nullptr; } - err = read_curve(EA_stuff[0].c_str(), - &(obj->EA), - &(obj->nEApoints), - obj->stiffXs, - obj->stiffYs); - if (err) - return nullptr; vector BA_stuff = moordyn::str::split(entries[4], '|'); unsigned int BA_N = BA_stuff.size(); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 108f526e..c85a02dd 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -43,6 +43,7 @@ set(CATCH2_TESTS math_tests state_tests testDecomposeString + syrope bodies beam conveying_fluid diff --git a/tests/Mooring/syrope/exp_wc/line.txt b/tests/Mooring/syrope/exp_wc/line.txt new file mode 100644 index 00000000..d841c79d --- /dev/null +++ b/tests/Mooring/syrope/exp_wc/line.txt @@ -0,0 +1,35 @@ +One single segment. +------------------------- LINE TYPES -------------------------------------------------- +LineType Diam MassDenInAir EA BA/-zeta EI Can Cat Cdn Cdt +(-) (m) (kg/m) (-) (Pa-s/-) (n-m^2) (-) (-) (-) (-) +poly 0.203 40 SYROPE:../owc.dat|1.53e8|23.12 5.0e10|1.0e5 0.0 0.0 0.0 0.0 0.0 +-------------------------- SYROPE WORKING CURVES ------------------------------------------------- +LineType WCType k1 k2 +(-) (-) (-) (-) +poly EXP 0.2 1.5 +----------------------- POINTS ---------------------------------------------- +Node Type X Y Z M V CdA CA +(-) (-) (m) (m) (m) (kg) (m^3) (m^2) (-) +1 Fixed 0.0 0 -5.0 0 0 0 0 +2 Coupled 1.0124577 0 -5.0 0 0 0 0 +-------------------------- LINES ------------------------------------------------- +Line LineType NodeA NodeB UnstrLen NumSegs Flags/Outputs +(-) (-) (-) (-) (m) (-) (-) +1 poly 2 1 1.0 1 tsc +-------------------------- SYROPE IC ------------------------------------------------- +Line Tmax Tmean +(-) (N) (N) +1 1.3886e6 4.339e5 +-------------------------- SOLVER OPTIONS--------------------------------------------------- +2 writeLog - Write a log file +0 g - gravitational acceleration +1025 rho - water density +1.0e-2 dtM - time step to use in mooring integration +3.0e6 kb - bottom stiffness +3.0e5 cb - bottom damping +70 WtrDpth - water depth +5.0 ICDfac - factor by which to scale drag coefficients during dynamic relaxation IC gen +1.0e-6 threshIC - threshold for IC convergence +100.0 TmaxIC - threshold for IC convergence +0.01 dtIC - Time lapse between convergence tests (s) +--------------------------- need this line ------------------------------------------------- \ No newline at end of file diff --git a/tests/Mooring/syrope/linear_wc/line.txt b/tests/Mooring/syrope/linear_wc/line.txt new file mode 100644 index 00000000..2ecc4e8e --- /dev/null +++ b/tests/Mooring/syrope/linear_wc/line.txt @@ -0,0 +1,35 @@ +One single segment. +------------------------- LINE TYPES -------------------------------------------------- +LineType Diam MassDenInAir EA BA/-zeta EI Can Cat Cdn Cdt +(-) (m) (kg/m) (-) (Pa-s/-) (n-m^2) (-) (-) (-) (-) +poly 0.203 40 SYROPE:../owc.dat|1.53e8|23.12 5.0e10|1.0e5 0.0 0.0 0.0 0.0 0.0 +-------------------------- SYROPE WORKING CURVES ------------------------------------------------- +LineType WCType k1 k2 +(-) (-) (-) (-) +poly LINEAR 1.25e8 0.0 +----------------------- POINTS ---------------------------------------------- +Node Type X Y Z M V CdA CA +(-) (-) (m) (m) (m) (kg) (m^3) (m^2) (-) +1 Fixed 0.0 0 -5.0 0 0 0 0 +2 Coupled 1.0148657 0 -5.0 0 0 0 0 +-------------------------- LINES ------------------------------------------------- +Line LineType NodeA NodeB UnstrLen NumSegs Flags/Outputs +(-) (-) (-) (-) (m) (-) (-) +1 poly 2 1 1.0 1 tsc +-------------------------- SYROPE IC ------------------------------------------------- +Line Tmax Tmean +(-) (N) (N) +1 1.3886e6 4.3394e5 +-------------------------- SOLVER OPTIONS--------------------------------------------------- +2 writeLog - Write a log file +0 g - gravitational acceleration +1025 rho - water density +1.0e-2 dtM - time step to use in mooring integration +3.0e6 kb - bottom stiffness +3.0e5 cb - bottom damping +70 WtrDpth - water depth +5.0 ICDfac - factor by which to scale drag coefficients during dynamic relaxation IC gen +1.0e-6 threshIC - threshold for IC convergence +100.0 TmaxIC - threshold for IC convergence +0.01 dtIC - Time lapse between convergence tests (s) +--------------------------- need this line ------------------------------------------------- \ No newline at end of file diff --git a/tests/Mooring/syrope/owc.dat b/tests/Mooring/syrope/owc.dat new file mode 100644 index 00000000..090a733d --- /dev/null +++ b/tests/Mooring/syrope/owc.dat @@ -0,0 +1,32 @@ +Strain Tension +(-) (N) +0.00000e+00 0.00000e+00 +2.06897e-03 1.71768e+05 +4.13793e-03 3.30952e+05 +6.20690e-03 4.78788e+05 +8.27586e-03 6.16510e+05 +1.03448e-02 7.45355e+05 +1.24138e-02 8.66556e+05 +1.44828e-02 9.81351e+05 +1.65517e-02 1.09097e+06 +1.86207e-02 1.19666e+06 +2.06897e-02 1.29964e+06 +2.27586e-02 1.40116e+06 +2.48276e-02 1.50245e+06 +2.68966e-02 1.60474e+06 +2.89655e-02 1.70927e+06 +3.10345e-02 1.81728e+06 +3.31034e-02 1.93000e+06 +3.51724e-02 2.04866e+06 +3.72414e-02 2.17450e+06 +3.93103e-02 2.30876e+06 +4.13793e-02 2.45267e+06 +4.34483e-02 2.60747e+06 +4.55172e-02 2.77439e+06 +4.75862e-02 2.95467e+06 +4.96552e-02 3.14954e+06 +5.17241e-02 3.36024e+06 +5.37931e-02 3.58800e+06 +5.58621e-02 3.83406e+06 +5.79310e-02 4.09965e+06 +6.00000e-02 4.38601e+06 \ No newline at end of file diff --git a/tests/Mooring/syrope/quadratic_wc/line.txt b/tests/Mooring/syrope/quadratic_wc/line.txt new file mode 100644 index 00000000..482cdeba --- /dev/null +++ b/tests/Mooring/syrope/quadratic_wc/line.txt @@ -0,0 +1,35 @@ +One single segment. +------------------------- LINE TYPES -------------------------------------------------- +LineType Diam MassDenInAir EA BA/-zeta EI Can Cat Cdn Cdt +(-) (m) (kg/m) (-) (Pa-s/-) (n-m^2) (-) (-) (-) (-) +poly 0.203 40 SYROPE:../owc.dat|1.53e8|23.12 5.0e10|1.0e5 0.0 0.0 0.0 0.0 0.0 +-------------------------- SYROPE WORKING CURVES ------------------------------------------------- +LineType WCType k1 k2 +(-) (-) (-) (-) +poly QUADRATIC 0.25 1.0 +----------------------- POINTS ---------------------------------------------- +Node Type X Y Z M V CdA CA +(-) (-) (m) (m) (m) (kg) (m^3) (m^2) (-) +1 Fixed 0.0 0 -5.0 0 0 0 0 +2 Coupled 1.0150575 0 -5.0 0 0 0 0 +-------------------------- LINES ------------------------------------------------- +Line LineType NodeA NodeB UnstrLen NumSegs Flags/Outputs +(-) (-) (-) (-) (m) (-) (-) +1 poly 2 1 1.0 1 tsc +-------------------------- SYROPE IC ------------------------------------------------- +Line Tmax Tmean +(-) (N) (N) +1 1.3886e6 4.339e5 +-------------------------- SOLVER OPTIONS--------------------------------------------------- +2 writeLog - Write a log file +0 g - gravitational acceleration +1025 rho - water density +1.0e-2 dtM - time step to use in mooring integration +3.0e6 kb - bottom stiffness +3.0e5 cb - bottom damping +70 WtrDpth - water depth +5.0 ICDfac - factor by which to scale drag coefficients during dynamic relaxation IC gen +1.0e-6 threshIC - threshold for IC convergence +100.0 TmaxIC - threshold for IC convergence +0.01 dtIC - Time lapse between convergence tests (s) +--------------------------- need this line ------------------------------------------------- \ No newline at end of file diff --git a/tests/syrope.cpp b/tests/syrope.cpp new file mode 100644 index 00000000..f110efdc --- /dev/null +++ b/tests/syrope.cpp @@ -0,0 +1,629 @@ +/* + * This file is based on the C++ driver to verify the Syrope implementation in MoorDyn-C, + * which is available at https://github.com/zhilongwei/moordyn-syrope-tests.git. + * A Python script is used there to check the L2-error and plot the results, + * whereas here we use Catch2 to check the L2-error only. + */ + +#define _USE_MATH_DEFINES +#include "MoorDyn2.h" + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +constexpr double kG = 9.80665; +constexpr double kMBL = 885.0e3 * kG; +constexpr double kTmax0 = 0.16 * kMBL; +constexpr double kL2Tol = 0.05; + +enum class WorkingCurveForm +{ + Linear, + Quadratic, + Exponential +}; + +constexpr double kLinearK1 = 1.25e8; +constexpr double kLinearK2 = 0.0; +constexpr double kQuadraticK1 = 0.25; +constexpr double kQuadraticK2 = 1.00; +constexpr double kExpK1 = 0.20; +constexpr double kExpK2 = 1.50; + +struct WcCase +{ + std::string name; + std::string input_file; + WorkingCurveForm form; + double eps_0; +}; + +struct SimulationResult +{ + Eigen::VectorXd times; + Eigen::VectorXd strain; + Eigen::VectorXd tension_output; +}; + +inline void +check_md(const int err, const std::string& msg) +{ + REQUIRE(err == MOORDYN_SUCCESS); + if (err != MOORDYN_SUCCESS) { + throw std::runtime_error(msg + " (err=" + std::to_string(err) + ")"); + } +} + +struct MoorDynRAII +{ + explicit MoorDynRAII(const std::string& input_file) + : sys(MoorDyn_Create(input_file.c_str())) + { + REQUIRE(sys != nullptr); + if (!sys) { + throw std::runtime_error("MoorDyn_Create failed for input: " + + input_file); + } + } + + ~MoorDynRAII() + { + if (sys) { + (void)MoorDyn_Close(sys); + } + } + + MoorDynRAII(const MoorDynRAII&) = delete; + MoorDynRAII& operator=(const MoorDynRAII&) = delete; + + MoorDyn sys = nullptr; +}; + +double +interpolate_clamped(double x, + const Eigen::VectorXd& xdata, + const Eigen::VectorXd& ydata) +{ + if (xdata.size() != ydata.size()) { + throw std::invalid_argument("interpolate_clamped: size mismatch"); + } + if (xdata.size() < 2) { + throw std::invalid_argument("interpolate_clamped: need >= 2 points"); + } + + for (Eigen::Index i = 1; i < xdata.size(); ++i) { + if (xdata[i] < xdata[i - 1]) { + throw std::invalid_argument( + "interpolate_clamped: xdata must be monotone increasing"); + } + } + + if (x <= xdata[0]) + return ydata[0]; + if (x >= xdata[xdata.size() - 1]) + return ydata[ydata.size() - 1]; + + Eigen::Index lo = 0; + Eigen::Index hi = xdata.size() - 1; + + while (hi - lo > 1) { + const Eigen::Index mid = lo + (hi - lo) / 2; + if (xdata[mid] <= x) + lo = mid; + else + hi = mid; + } + + const double x0 = xdata[lo], x1 = xdata[hi]; + const double y0 = ydata[lo], y1 = ydata[hi]; + const double dx = x1 - x0; + + if (std::abs(dx) <= std::numeric_limits::epsilon()) { + return 0.5 * (y0 + y1); + } + + const double t = (x - x0) / dx; + return y0 + t * (y1 - y0); +} + +double +find_mean_tension(double strain, + double Tmean, + double Tmax, + const Eigen::VectorXd& owc_strains, + const Eigen::VectorXd& owc_tensions, + WorkingCurveForm wc_form) +{ + (void)Tmean; + + double k1 = 0.0; + double k2 = 0.0; + switch (wc_form) { + case WorkingCurveForm::Linear: + k1 = kLinearK1; + k2 = kLinearK2; + break; + case WorkingCurveForm::Quadratic: + k1 = kQuadraticK1; + k2 = kQuadraticK2; + break; + case WorkingCurveForm::Exponential: + k1 = kExpK1; + k2 = kExpK2; + break; + default: + throw std::invalid_argument( + "find_mean_tension: unknown WorkingCurveForm"); + } + + const double eps_max = interpolate_clamped(Tmax, owc_tensions, owc_strains); + const double eps0 = owc_strains[0]; + + double eps_min = eps0 + k1 * (eps_max - eps0); + if (wc_form == WorkingCurveForm::Linear && k1 >= 1.0) { + eps_min = eps_max - Tmax / k1; + } + + const double denom = eps_max - eps_min; + if (std::abs(denom) <= std::numeric_limits::epsilon()) { + return interpolate_clamped(strain, owc_strains, owc_tensions); + } + + double xi = (strain - eps_min) / denom; + xi = std::clamp(xi, 0.0, 1.0); + + double tension_wc = 0.0; + + switch (wc_form) { + case WorkingCurveForm::Linear: + tension_wc = Tmax * xi; + break; + case WorkingCurveForm::Quadratic: + tension_wc = Tmax * xi * (k2 * xi + (1.0 - k2)); + break; + case WorkingCurveForm::Exponential: { + const double den = std::exp(k2) - 1.0; + if (std::abs(den) < 1e-14) + tension_wc = Tmax * xi; + else + tension_wc = Tmax * (std::exp(k2 * xi) - 1.0) / den; + break; + } + default: + throw std::invalid_argument( + "find_mean_tension: unknown WorkingCurveForm"); + } + + if (std::abs(tension_wc - Tmax) <= 1e-12 * (std::max)(1.0, Tmax)) { + return interpolate_clamped(strain, owc_strains, owc_tensions); + } + + return tension_wc; +} + +std::pair +read_strain_tension_table(const std::string& path) +{ + std::ifstream in(path); + if (!in) { + throw std::runtime_error("Cannot open file: " + path); + } + + std::string line; + if (!std::getline(in, line)) + throw std::runtime_error("Missing header line 1 in: " + path); + if (!std::getline(in, line)) + throw std::runtime_error("Missing header line 2 in: " + path); + + std::vector strain, tension; + strain.reserve(1024); + tension.reserve(1024); + + std::size_t lineno = 2; + while (std::getline(in, line)) { + ++lineno; + if (line.empty()) + continue; + + const auto first = line.find_first_not_of(" \t\r"); + if (first == std::string::npos) + continue; + const char c0 = line[first]; + if (c0 == '#' || c0 == '!') + continue; + + std::istringstream iss(line); + double eps = 0.0, ten = 0.0; + if (!(iss >> eps >> ten)) { + throw std::runtime_error("Failed to parse line " + + std::to_string(lineno) + " in " + path + + ": '" + line + "'"); + } + strain.push_back(eps); + tension.push_back(ten); + } + + if (strain.empty()) { + throw std::runtime_error("No data rows found in: " + path); + } + + Eigen::VectorXd eps(static_cast(strain.size())); + Eigen::VectorXd ten(static_cast(tension.size())); + for (Eigen::Index i = 0; i < eps.size(); ++i) { + eps[i] = strain[static_cast(i)]; + ten[i] = tension[static_cast(i)]; + } + + return { eps, ten }; +} + +Eigen::VectorXd +load_seg_te_column(const std::string& path) +{ + std::ifstream in(path); + if (!in) { + throw std::runtime_error("Cannot open " + path); + } + + std::string line; + if (!std::getline(in, line)) { + throw std::runtime_error("Missing header line in " + path); + } + + std::istringstream hdr(line); + std::vector headers; + for (std::string tok; hdr >> tok;) + headers.push_back(tok); + if (headers.empty()) { + throw std::runtime_error("Empty header line in " + path); + } + + const std::regex re(R"(^Seg\d+Te$)"); + int col = -1; + for (size_t i = 0; i < headers.size(); ++i) { + if (std::regex_match(headers[i], re)) { + col = static_cast(i); + break; + } + } + if (col < 0) { + throw std::runtime_error("No SegTe column found in " + path); + } + + if (!std::getline(in, line)) { + throw std::runtime_error("Missing units line in " + path); + } + + std::istringstream units_iss(line); + std::vector units; + for (std::string tok; units_iss >> tok;) + units.push_back(tok); + + if (static_cast(units.size()) <= col) { + throw std::runtime_error( + "Units line has fewer columns than header in " + path); + } + if (units[static_cast(col)] != "(N)") { + throw std::runtime_error("Unexpected unit for " + + headers[static_cast(col)] + ": got '" + + units[static_cast(col)] + + "', expected '(N)' in " + path); + } + + std::vector vals; + vals.reserve(2048); + + std::size_t lineno = 2; + while (std::getline(in, line)) { + ++lineno; + if (line.empty()) + continue; + + std::istringstream row(line); + double v = 0.0; + for (int i = 0; i <= col; ++i) { + if (!(row >> v)) { + throw std::runtime_error( + "Row has fewer numeric columns than expected at line " + + std::to_string(lineno) + " in " + path + ": '" + line + "'"); + } + } + vals.push_back(v); + } + + Eigen::VectorXd out(static_cast(vals.size())); + for (Eigen::Index i = 0; i < out.size(); ++i) { + out[i] = vals[static_cast(i)]; + } + return out; +} + +double +jonswap_psd(double Hs, + double Tp, + double gamma, + double omega, + double sigma_a = 0.07, + double sigma_b = 0.09) +{ + if (!(omega > 0.0) || !(Tp > 0.0) || !(gamma > 0.0) || !(Hs >= 0.0)) { + return 0.0; + } + + const double wp = 2.0 * M_PI / Tp; + const double A = 1.0 - 0.287 * std::log(gamma); + const double sigma = (omega > wp) ? sigma_b : sigma_a; + const double rw = (omega - wp) / (sigma * wp); + const double a = std::pow(wp / omega, 4.0); + + return A * (5.0 / 16.0) * Hs * Hs * (a / omega) * std::exp(-1.25 * a) * + std::pow(gamma, std::exp(-0.5 * rw * rw)); +} + +double +slow_strain(double t, double seg_dur, double eps0) +{ + if (seg_dur <= 0.0) + return eps0; + + if (t < 1.0 * seg_dur) + return eps0; + if (t < 2.0 * seg_dur) + return eps0 + 2.0 * eps0 * (t - 1.0 * seg_dur) / seg_dur; + if (t < 3.0 * seg_dur) + return eps0 + 2.0 * eps0; + if (t < 4.0 * seg_dur) + return eps0 + 2.0 * eps0 - 1.0 * eps0 * (t - 3.0 * seg_dur) / seg_dur; + if (t < 5.0 * seg_dur) + return eps0 + 1.0 * eps0; + return eps0 + 1.0 * eps0 + 1.5 * eps0 * (t - 5.0 * seg_dur) / seg_dur; +} + +Eigen::VectorXd +surface_elevation_from_jonswap(double Hs, + double Tp, + double gamma, + const Eigen::VectorXd& times, + unsigned int n_omega_edges, + double omega_min, + double omega_max, + std::mt19937& rng) +{ + if (n_omega_edges < 2) { + throw std::runtime_error("n_omega_edges must be >= 2"); + } + if (!(omega_min > 0.0 && omega_max > omega_min)) { + throw std::runtime_error("Invalid omega_min/omega_max"); + } + + const Eigen::VectorXd omega_edges = + Eigen::VectorXd::LinSpaced(n_omega_edges, omega_min, omega_max); + const double domega = omega_edges[1] - omega_edges[0]; + + const Eigen::VectorXd omega_mid = + 0.5 * (omega_edges.head(n_omega_edges - 1) + + omega_edges.tail(n_omega_edges - 1)); + + const Eigen::VectorXd S = omega_mid.unaryExpr( + [&](double om) { return jonswap_psd(Hs, Tp, gamma, om); }); + + const Eigen::VectorXd a = (2.0 * S * domega).cwiseSqrt(); + + std::uniform_real_distribution unif(0.0, 2.0 * M_PI); + Eigen::VectorXd phi(n_omega_edges - 1); + for (Eigen::Index i = 0; i < phi.size(); ++i) { + phi[i] = unif(rng); + } + + Eigen::VectorXd eta = Eigen::VectorXd::Zero(times.size()); + for (Eigen::Index i = 0; i < times.size(); ++i) { + const double t = times[i]; + eta[i] = ((omega_mid * t + phi).array().cos() * a.array()).sum(); + } + return eta; +} + +std::string +line1_output_from_input(const std::string& input_path) +{ + const std::size_t sep = input_path.find_last_of("/\\"); + const std::size_t base = (sep == std::string::npos) ? 0 : (sep + 1); + + const std::size_t dot = input_path.find_last_of('.'); + const bool has_ext = (dot != std::string::npos) && (dot > base); + + const std::string stem = has_ext ? input_path.substr(0, dot) : input_path; + return stem + "_Line1.out"; +} + +SimulationResult +run_case(const WcCase& c, + bool superimpose_fast) +{ + MoorDynRAII md(c.input_file); + + unsigned int ndof = 0; + check_md(MoorDyn_NCoupledDOF(md.sys, &ndof), + "MoorDyn_NCoupledDOF failed for: " + c.input_file); + + auto fairlead = MoorDyn_GetPoint(md.sys, 2); + REQUIRE(fairlead != nullptr); + + double r[3] = { 0.0, 0.0, 0.0 }; + double dr[3] = { 0.0, 0.0, 0.0 }; + check_md(MoorDyn_GetPointPos(fairlead, r), + "MoorDyn_GetPointPos failed for: " + c.input_file); + + const double seg_dur = 3.0 * 3600.0; + const double total_dur = 6.0 * seg_dur; + const double dt0 = 0.1; + const unsigned int nsteps = static_cast(total_dur / dt0); + + const Eigen::VectorXd times = Eigen::VectorXd::LinSpaced( + static_cast(nsteps) + 1, 0.0, total_dur); + + constexpr unsigned int n_omega_edges = 300; + std::mt19937 rng_wf(12345); + std::mt19937 rng_lf(54321); + + const double Hs_WF = 5.0; + const double Tp_WF = 12.0; + const double Tz_WF = Tp_WF / 1.402; + const double omega_min_WF = 1.0 / Tz_WF; + const double omega_max_WF = 20.0 / Tz_WF; + + const Eigen::VectorXd eta_WF = surface_elevation_from_jonswap(Hs_WF, + Tp_WF, + 3.3, + times, + n_omega_edges, + omega_min_WF, + omega_max_WF, + rng_wf); + + const double Hs_LF = 20.0; + const double Tp_LF = 120.0; + const double Tz_LF = Tp_LF / 1.402; + const double omega_min_LF = 1.0 / Tz_LF; + const double omega_max_LF = 20.0 / Tz_LF; + + const Eigen::VectorXd eta_LF = surface_elevation_from_jonswap(Hs_LF, + Tp_LF, + 3.3, + times, + n_omega_edges, + omega_min_LF, + omega_max_LF, + rng_lf); + + const double x0 = 1.0; + const double scale_WF = 4.00e-4; + const double scale_LF = 0.80e-4; + + Eigen::VectorXd strain_slow(times.size()); + for (Eigen::Index i = 0; i < times.size(); ++i) { + strain_slow[i] = slow_strain(times[i], seg_dur, c.eps_0); + } + + const Eigen::VectorXd strain = superimpose_fast + ? (strain_slow + scale_WF * eta_WF + + scale_LF * eta_LF) + : strain_slow; + const Eigen::VectorXd x = + x0 * (Eigen::VectorXd::Ones(strain.size()) + strain); + + r[0] = x[0]; + dr[0] = (x[1] - x[0]) / dt0; + dr[1] = 0.0; + dr[2] = 0.0; + + check_md(MoorDyn_Init(md.sys, r, dr), "MoorDyn_Init failed for: " + c.input_file); + + double t = 0.0; + double f[3] = { 0.0, 0.0, 0.0 }; + + for (unsigned int i = 0; i < nsteps; ++i) { + r[0] = x[static_cast(i)]; + dr[0] = (x[static_cast(i) + 1] - + x[static_cast(i)]) / + dt0; + + double t_in = t; + double dt_in = dt0; + + check_md(MoorDyn_Step(md.sys, r, dr, f, &t_in, &dt_in), + "MoorDyn_Step failed for: " + c.input_file); + + REQUIRE(std::abs(dt_in - dt0) <= 1e-12); + t = t_in; + } + + check_md(MoorDyn_Close(md.sys), "MoorDyn_Close failed for: " + c.input_file); + md.sys = nullptr; + + SimulationResult out; + out.times = times; + out.strain = strain; + out.tension_output = load_seg_te_column(line1_output_from_input(c.input_file)); + return out; +} + +double +relative_l2(const Eigen::VectorXd& ref, const Eigen::VectorXd& val) +{ + const Eigen::Index n = (std::min)(ref.size(), val.size()); + REQUIRE(n >= 2); + const double denom = ref.head(n).squaredNorm(); + REQUIRE(denom > 0.0); + return std::sqrt((ref.head(n) - val.head(n)).squaredNorm() / denom); +} + +} // namespace + +TEST_CASE("Syrope tests", "[syrope][working-curve]") +{ + const auto [owc_strains, owc_tensions] = + read_strain_tension_table("Mooring/syrope/owc.dat"); + + const std::vector cases = { + { "Linear", + "Mooring/syrope/linear_wc/line.txt", + WorkingCurveForm::Linear, + 1.48657e-02 }, + { "Quadratic", + "Mooring/syrope/quadratic_wc/line.txt", + WorkingCurveForm::Quadratic, + 1.50575e-02 }, + { "Exponential", + "Mooring/syrope/exp_wc/line.txt", + WorkingCurveForm::Exponential, + 1.24577e-02 }, + }; + + const bool superimpose_fast = true; + + for (const auto& c : cases) { + DYNAMIC_SECTION("Case: " << c.name) + { + const auto sim = run_case(c, superimpose_fast); + const Eigen::Index n = + (std::min)(sim.tension_output.size(), sim.strain.size()); + REQUIRE(n >= 2); + + Eigen::VectorXd tmax_mean(n); + tmax_mean[0] = kTmax0; + for (Eigen::Index i = 1; i < n; ++i) { + tmax_mean[i] = (std::max)(tmax_mean[i - 1], sim.tension_output[i]); + } + + Eigen::VectorXd tension_analytical(n); + for (Eigen::Index i = 0; i < n; ++i) { + tension_analytical[i] = find_mean_tension(sim.strain[i], + sim.tension_output[i], + tmax_mean[i], + owc_strains, + owc_tensions, + c.form); + } + + const double l2 = relative_l2(tension_analytical, sim.tension_output); + INFO("Relative L2 error = " << l2); + REQUIRE(l2 < kL2Tol); + } + } +}