fix/feat: Fixed viscoelastic bug by initalizing sub-segment stretch as non-zero, added viscoelastic test based on polyester test

Author: RyanDavies19

source/Line.cpp

@@ -247,9 +247,6 @@ Line::setup(int number_in,
 		phi_dot.assign(N + 1, 0.0);     // synch model frequency 
 		yd_rms_old.assign(N + 1, 0.0);  // node old yd_rms
 		ydd_rms_old.assign(N + 1, 0.0); // node old ydd_rms
-
-		// give a random distribution between 0 and 2pi for inital phase of lift force to avoid initial transient
-		for (unsigned int i = 1; i < N; i++) phi[i] = (i/moordyn::real(N))*2*pi; // internal nodes only. Change if end node viv included
 	}
 
 	// ensure end moments start at zero
@@ -580,6 +577,20 @@ Line::initialize()
 		LOGDBG << "Vertical linear initial profile for Line " << number << endl;
 	}
 
+	// If using the viscoelastic model, 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 (ElasticMod > 1) {
+		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 VIV model, initalize as a distribution between 0 and 2pi for inital phase of lift force to avoid initial transient
+	if (Cl > 0)
+		for (unsigned int i = 1; i < N; i++) phi[i] = (i/moordyn::real(N))*2*pi; // internal nodes only. Change if end node viv included
+
 	LOGMSG << "Initialized Line " << number << endl;
 
 	// NOTE: becasue Line.hpp is the only user of this function, no need to return any state information
@@ -755,8 +766,10 @@ Line::setState(const InstanceStateVarView state)
 	for (unsigned int i = 0; i < (ElasticMod > 1 ? N : N-1); i++) {
 		// node number is i+1
 		// segment number is i
-		r[i+1] = state.row(i).head<3>();
-		rd[i+1] = state.row(i).segment<3>(3);
+		if (i < N-1) { // only assign the internal nodes
+			r[i+1] = state.row(i).head<3>();
+			rd[i+1] = state.row(i).segment<3>(3);
+		}
 
 		if (ElasticMod > 1) dl_1[i] = state.row(i).tail<1>()[0]; // [0] needed becasue tail<1> returns a one element vector. Viscoelastic state is always the last element in the row
 
@@ -771,6 +784,12 @@ Line::setState(const InstanceStateVarView state)
 void
 Line::setEndKinematics(vec pos, vec vel, EndPoints end_point)
 {
+	if (t < 5 && int(t / dtm0) % 20 == 0) {
+		LOGDBG << "------" << endl;
+    	LOGDBG << "Setting end kinematics for end point " << end_point << endl;
+    	LOGDBG << "Position: " << pos.transpose() << ", Velocity: " << vel.transpose() << endl;
+	}
+
 	switch (end_point) {
 		case ENDPOINT_TOP:
 			endTypeB = PINNED; // indicate pinned
@@ -925,7 +944,8 @@ Line::getStateDeriv(InstanceStateVarView drdt)
 	for (unsigned int i = 0; i < N; i++) {
 		// calculate current (Stretched) segment lengths and unit tangent
 		// vectors (qs) for each segment (this is used for bending and stiffness calculations)
-		lstr[i] = unitvector(qs[i], r[i], r[i + 1]);
+
+		lstr[i] = unitvector(qs[i], r[i], r[i + 1]); // if using the viscoelastic model this is redundant for the first time step, as it is also called by Line::initalize
 
 		ldstr[i] = qs[i].dot(rd[i + 1] - rd[i]); // strain rate of segment
 

tests/Mooring/polyester/visco.txt

@@ -11,7 +11,7 @@ vbeta = KRD2 * MBL = 0.35 * 25.3e6 = 8.855e6
 ----------------------- LINE TYPES ------------------------------------------
 TypeName         Diam    Mass/m     EA                      BA/-zeta   EI         Cd     Ca     CdAx    CaAx
 (name)           (m)     (kg/m)     (N)                     (N-s/-)    (N-m^2)    (-)    (-)    (-)     (-)
-Polyester_290mm  0.217   52.0       3.39e8|4.048e8|8.855e6  -1.0|1e11  0.0        1.6    1.79   0.008   0.0    
+Polyester_290mm  0.217   52.0       3.39e8|4.048e8|8.855e6  -1|1e10  0.0        1.6    1.79   0.008   0.0    
 ---------------------- POINT PROPERTIES --------------------------------
 ID    Type      X       Y       Z         Mass   Volume  CdA    Ca
 (#)   (-)       (m)     (m)     (m)       (kg)   (mˆ3)   (m^2)  (-)
@@ -20,23 +20,19 @@ ID    Type      X       Y       Z         Mass   Volume  CdA    Ca
 ---------------------- LINES ----------------------------------------
 ID   LineType         AttachA  AttachB   UnstrLen  NumSegs  LineOutputs
 (#)  (name)           (#)      (#)       (m)       (-)      (-)
-1    Polyester_290mm  1        2         1400.0    70       p
+1    Polyester_290mm  1        2         1400.0    70       -
 ---------------------- OPTIONS -----------------------------------------
 2             writeLog             Write a log file
-0.0005         dtM                  time step to use in mooring integration (s)
+0.00005         dtM                  time step to use in mooring integration (s)
 1.0e5         kBot                 bottom stiffness (Pa/m)
 1.0e4         cBot                 bottom damping (Pa-s/m)
 9.80665       g                    gravitational constant (m/s^2)
 1025.0        WtrDnsty             water density (kg/m^3)
 372.3         WtrDpth              water depth (m)
 4.0           CdScaleIC            factor by which to scale drag coefficients during dynamic relaxation (-)
-1             ICgenDynamic
 0.5           FrictionCoefficient  general bottom friction coefficient, as a start (-)
-5.0           dtIC                 time interval for analyzing convergence during IC gen (s)
-00           TmaxIC               max time for ic gen (s)
-0.0000001          threshIC             threshold for IC convergence (-)
-0.001          dtOut
-------------------------- Outputs -------------------------
-Fairten1
-POINT1PX
+0.1           dtIC                 time interval for analyzing convergence during IC gen (s)
+100           TmaxIC               max time for ic gen (s)
+0.01          threshIC             threshold for IC convergence (-)
+0.01          dtOut
 ------------------------- need this line -------------------------------------- 

tests/polyester.cpp

@@ -170,12 +170,85 @@ TEST_CASE("Ramp up, stabilization and cycles")
 }
 
 
-// TEST_CASE("Visco-elastic testing")
-// {
-// 	// TODO: static dynamic stiffness
+TEST_CASE("Visco-elastic testing")
+{
+	// this is the same as the above test, but the static and dynamic stiffness coefficients have been added to MD instead of calculated externally. 
+	MoorDyn system = MoorDyn_Create("Mooring/polyester/visco.txt");
+	REQUIRE(system);
+	unsigned int n_dof;
+	REQUIRE(MoorDyn_NCoupledDOF(system, &n_dof) == MOORDYN_SUCCESS);
+	REQUIRE(n_dof == 3);
+
+	double r[3], dr[3];
+	auto anchor = MoorDyn_GetPoint(system, 1);
+	REQUIRE(anchor);
+	auto fairlead = MoorDyn_GetPoint(system, 2);
+	REQUIRE(fairlead);
+	auto line = MoorDyn_GetLine(system, 1);
+	REQUIRE(line);
+
+	double l0;
+	REQUIRE(MoorDyn_GetLineUnstretchedLength(line, &l0) == MOORDYN_SUCCESS);
+
+	REQUIRE(MoorDyn_GetPointPos(fairlead, r) == MOORDYN_SUCCESS);
+	std::fill(dr, dr + 3, 0.0);
+	REQUIRE(MoorDyn_Init(system, r, dr) == MOORDYN_SUCCESS);
 
-// 	// TODO: load dependent dynamic stiffness
+	// Read the csv with the motions
+	std::ifstream csv("Mooring/polyester/motion.csv");
+	aria::csv::CsvParser csv_parser(csv);
 
-// 	REQUIRE(MoorDyn_Close(system) == MOORDYN_SUCCESS);
-// }
+	std::vector<double> tdata, xdata;
+	unsigned int skip = 2;
+	for (auto& row : csv_parser) {
+		if (skip > 0) {
+			skip--;
+			continue;
+		}
+		if (row.size() < 2)
+			continue;
+		tdata.push_back(std::stod(row.at(0)));
+		xdata.push_back(std::stod(row.at(1)));
+	}
+
+	// Time to move
+	double t = 0.0, dt;
+	// For checking the performance
+	std::deque<double> times;  
+	std::deque<double> tensions;  
+	std::deque<double> ttimes(TTIMES);
+	std::deque<double> tmeans(TMEANS); // tension means
+	for (unsigned int i = 0; i < tdata.size(); i++) {
+		REQUIRE(MoorDyn_GetPointPos(fairlead, r) == MOORDYN_SUCCESS);
+		double t_dst = tdata.at(i);
+		double dt = t_dst - t;
+		double x_dst = xdata.at(i);
+		dr[0] = (x_dst - r[0]) / dt;
+		double f[3];
+		REQUIRE(MoorDyn_Step(system, r, dr, f, &t, &dt) == MOORDYN_SUCCESS);
+		times.push_back(t);
+		tensions.push_back(get_average_tension(line, anchor, fairlead));
+
+		// Let's check and tweak the line if there is info enough
+		if (times.back() - times.front() < TC)
+			continue;
+
+		double tension = 0.0;
+		for (auto f : tensions)
+			tension += f;
+		tension /= tensions.size();
+		times.pop_front();
+		tensions.pop_front();
+
+		if (t >= ttimes.front()) {
+			const double tmean = tmeans.front();
+			ttimes.pop_front();
+			tmeans.pop_front();
+			REQUIRE(
+				fabs(tension / MBL - tmean) < 0.025);
+		}
+	}
+
+	REQUIRE(MoorDyn_Close(system) == MOORDYN_SUCCESS);
+}