Repository: FloatingArrayDesign/MoorDyn (language: python)

Issue:
Body negative local drag coefficients
Strictly speaking, this is a bug that I encountered on our internal MoorDyn fork, which is somewhat out of date. I don't see anything in the latest public code that would address this, but I haven't directly tested this behavior with the latest branch of this repo.

When bodies are in certain orientations, the computed global drag coefficients can end up having negative values.

The code at fault is
https://github.com/FloatingArrayDesign/MoorDyn/blob/99cd2b6e6afa7e5d58e2ac7ac3a4b05e127f1f72/source/Body.cpp#L496-L501

An example of a case where this fails to work is when the body's pitch is -90 degrees. The rotation matrix in this case will be
```
[  0.0,  0.0, -1.0;
   0.0,  1.0,  0.0;
   1.0,  0.0,  0.0 ]
```

The `cda` variable will look like `bodyCda * [1, 1, -1, 1, 1, -1]`. Which leads to a negative drag coefficient in the z axis, which is a destabilizing force that is both not physical and will eventually cause the simulation to error out.

I'm making this issue because I'm not sure what the correct solution to this issue is. I think it's valid to just take the absolute value of the rotated drag coefficients, but I'll defer to others for how they want to address this.


Discussion / hints:

Fix patch:
diff --git a/source/Body.cpp b/source/Body.cpp
index 836c9b69..5d7203f3 100644
--- a/source/Body.cpp
+++ b/source/Body.cpp
@@ -581,7 +581,7 @@ Body::doRHS()
 	cda(Eigen::seqN(0, 3)) = OrMat.transpose() * bodyCdA.head<3>();
 	cda(Eigen::seqN(3, 3)) = OrMat.transpose() * bodyCdA.tail<3>();
 	F6net +=
-	    0.5 * env->rho_w * vi.cwiseProduct(vi.cwiseAbs()).cwiseProduct(cda);
+	    0.5 * env->rho_w * vi.cwiseProduct((vi.cwiseProduct(cda)).cwiseAbs());
 
 	// Get contributions from any points attached to the body
 	for (auto attached : attachedP) {


Test patch:
diff --git a/tests/Mooring/body_tests/bodyDrag.txt b/tests/Mooring/body_tests/bodyDrag.txt
new file mode 100644
index 00000000..60b24cda
--- /dev/null
+++ b/tests/Mooring/body_tests/bodyDrag.txt
@@ -0,0 +1,30 @@
+A simple test to ensure that body drag coefficients don't go negative in certain orientations
+------------------------- LINE TYPES --------------------------------------------------
+LineType    Diam      MassDenInAir    EA        BA/-zeta    EI      tbd   Can   Cat  Cdn   Cdt  
+(-)         (m)       (kg/m)          (N)       (Pa-s/-)    (n-m^2) (-)   (-)   (-)  (-)   (-)  
+spring      0.1       10.0            500.0     -0.4        0.0     0.0   0.0   0.0  0.0   0.00
+---------------------------- BODIES -----------------------------------------------------
+ID   Attachment  X0     Y0    Z0     r0      p0     y0     Mass  CG*   I*      Volume   CdA*   Ca
+(#)     (-)      (m)    (m)   (m)   (deg)   (deg)  (deg)   (kg)  (m)  (kg-m^2)  (m^3)   (m^2)  (-)
+1       free      0     0      0     0       -90      0      1e1     0    5        0     0.1    0
+----------------------- POINTS ----------------------------------------------
+Node      Type       X         Y          Z        M        V          CdA    CA
+(-)       (-)        (m)       (m)        (m)      (kg)     (m^3)      (m^2)  (-)
+1         Body1      0.0       0          0        0        0          0      0
+2         Fixed      0.0       0         -5        0        0          0      0
+-------------------------- LINES -------------------------------------------------
+Line     LineType NodeA     NodeB  UnstrLen  NumSegs     Flags/Outputs
+(-)      (-)       (-)       (-)   (m)         (-)          (-)
+1        spring      1         2    4.5        1           -
+-------------------------- SOLVER OPTIONS---------------------------------------------------
+2        writeLog     - Write a log file
+9.81        g            - No gravity
+0.002   dtM          - time step to use in mooring integration
+3.0e6    kb           - bottom stiffness
+3.0e5    cb           - bottom damping
+70       WtrDpth      - water depth
+3.0      ICDfac       - factor by which to scale drag coefficients during dynamic relaxation IC gen
+0.015    threshIC     - threshold for IC convergence
+0.0    TmaxIC       - threshold for IC convergence
+0.01     dtIC         - Time lapse between convergence tests (s)
+--------------------------- need this line -------------------------------------------------
diff --git a/tests/bodies.cpp b/tests/bodies.cpp
index a02e1fa6..83928458 100644
--- a/tests/bodies.cpp
+++ b/tests/bodies.cpp
@@ -386,7 +386,9 @@ rotatingBody(SeriesWriter* series_writer)
 	double angle = moordyn::rad2deg * 2 * acos(q.w());
 	double denom = (sqrt(1 - q.w() * q.w()));
 	moordyn::vec3 axis{ q.x() / denom, q.y() / denom, q.z() / denom };
-	if (!(abs(axis.x()) > 0.85 && abs(axis.y()) < 0.2 && abs(axis.z()) < 0.2)) { // if we are just checking axis direction then +/- does not matter
+	if (!(abs(axis.x()) > 0.85 && abs(axis.y()) < 0.2 &&
+	      abs(axis.z()) < 0.2)) { // if we are just checking axis direction then
+		                          // +/- does not matter
 		cerr << "The final rotation of the body in angle axis form should "
 		        "have an axis in the x direction, but axis is "
 		     << axis.transpose() << endl;
@@ -421,19 +423,19 @@ rotatingBody(SeriesWriter* series_writer)
 /**
  * @brief Compares inertial deflection of a pinned body to analytical solution
  *
- * 
- * The coupled pinned body that is massless and volumeless with a rod fixed to it
- * is moved with constant acceleration in a vaccum (0 water density). 
- * The resulting avg inertial deflection should match an analytical solution of 
+ *
+ * The coupled pinned body that is massless and volumeless with a rod fixed to
+ * it is moved with constant acceleration in a vaccum (0 water density). The
+ * resulting avg inertial deflection should match an analytical solution of
  * theta = arctan(-accel/gravity)
- * 
- * This solution was derived both with Newtonian and LaGrangian mechanics. 
- * It is the same as the pendulum on an accelerating cart problem. 
- * 
- * This only tests the inertial properties of pinned bodies, other tests deal 
+ *
+ * This solution was derived both with Newtonian and LaGrangian mechanics.
+ * It is the same as the pendulum on an accelerating cart problem.
+ *
+ * This only tests the inertial properties of pinned bodies, other tests deal
  * with hydrodynamics and general body properties
  *
- * 
+ *
  * @param series_writer
  * @return true
  * @return false
@@ -456,13 +458,14 @@ pinnedBody(SeriesWriter* series_writer)
 		cerr << "Failure getting NCoupledDOF: " << err << endl;
 		return false;
 	}
-	if (n_dof != 6) { // rotational DOF are ignored by MDC, same as a coupled pinned rods
+	if (n_dof !=
+	    6) { // rotational DOF are ignored by MDC, same as a coupled pinned rods
 		cerr << "Expected 6 DOFs but got " << n_dof << endl;
 		return false;
 	}
 
-	moordyn::vec6 x{ 0, 0, -5, 0, 0, 0};
-	moordyn::vec6 xd{ 0, 0, 0 , 0, 0, 0};
+	moordyn::vec6 x{ 0, 0, -5, 0, 0, 0 };
+	moordyn::vec6 xd{ 0, 0, 0, 0, 0, 0 };
 	double f[6];
 
 	err = MoorDyn_Init(system, x.data(), xd.data());
@@ -488,8 +491,8 @@ pinnedBody(SeriesWriter* series_writer)
 
 	while (t < 50.0) {
 
-		x[1] = 0.5*accel*pow(t, 2);
-		xd[1] = accel*t;
+		x[1] = 0.5 * accel * pow(t, 2);
+		xd[1] = accel * t;
 		err = MoorDyn_Step(system, x.data(), xd.data(), f, &t, &dt);
 		if (err != MOORDYN_SUCCESS) {
 			cerr << "Failure during the mooring initialization: " << err
@@ -504,20 +507,24 @@ pinnedBody(SeriesWriter* series_writer)
 		roll.push_back(r[3]);
 
 		if (i >= 30) { // after the simulation has run for a few time steps
-			// When local min or max of oscillation, indicates half of an oscialltion has occured
-			local_min_max = (((roll[i]-roll[i-1])/dt) * ((roll[i-1]-roll[i-2])/dt)) < 0; 
-			if (local_min_max) j++;
+			// When local min or max of oscillation, indicates half of an
+			// oscialltion has occured
+			local_min_max = (((roll[i] - roll[i - 1]) / dt) *
+			                 ((roll[i - 1] - roll[i - 2]) / dt)) < 0;
+			if (local_min_max)
+				j++;
 		}
-		if (j > 3) break; // after 2 full oscillations
+		if (j > 3)
+			break; // after 2 full oscillations
 
 		t = t + dt;
 		i++;
 	}
-	double theta = atan(-accel/9.80665);
+	double theta = atan(-accel / 9.80665);
 	double average = reduce(roll.begin(), roll.end()) / roll.size();
 	if (abs(average - theta) > 0.001) {
-		cerr << "Pinned body inertial deflection should be "
-		<< theta << " but it is " << average << endl;
+		cerr << "Pinned body inertial deflection should be " << theta
+		     << " but it is " << average << endl;
 		return false;
 	}
 
@@ -533,6 +540,65 @@ pinnedBody(SeriesWriter* series_writer)
 	return true;
 }
 
+bool
+bodyDrag(SeriesWriter* series_writer)
+{
+	int err;
+	cout << endl << " => " << __PRETTY_FUNC_NAME__ << "..." << endl;
+
+	MoorDyn system = MoorDyn_Create("Mooring/body_tests/bodyDrag.txt");
+	if (!system) {
+		cerr << "Failure Creating the Mooring system" << endl;
+		return false;
+	}
+
+	unsigned int n_dof;
+	err = MoorDyn_NCoupledDOF(system, &n_dof);
+	if (err != MOORDYN_SUCCESS) {
+		cerr << "Failure getting NCoupledDOF: " << err << endl;
+		return false;
+	}
+	if (n_dof != 0) {
+		cerr << "Expected 0 DOFs but got " << n_dof << endl;
+		return false;
+	}
+
+	double f[3];
+
+	err = MoorDyn_Init(system, nullptr, nullptr);
+	if (err != MOORDYN_SUCCESS) {
+		cerr << "Failure during the mooring initialization: " << err << endl;
+		return false;
+	}
+
+	if (!write_system_vtk(system, 0, series_writer)) {
+		return false;
+	}
+
+	double t = 0, dt = 0.1;
+	double max_t = 5;
+	while (t < max_t) {
+		// do one outer time step just to make sure everything is settled
+		err = MoorDyn_Step(system, nullptr, nullptr, f, &t, &dt);
+		if (err != MOORDYN_SUCCESS) {
+			cerr << "Failure during the mooring dynamics: " << err << endl;
+			return false;
+		}
+
+		if (!write_system_vtk(system, t, series_writer)) {
+			return false;
+		}
+	}
+
+	err = MoorDyn_Close(system);
+	if (err != MOORDYN_SUCCESS) {
+		cerr << "Failure closing Moordyn: " << err << endl;
+		return false;
+	}
+
+	return true;
+}
+
 /** @brief Runs all the test
  * @return 0 if the tests have ran just fine. The index of the failing test
  * otherwise
@@ -566,6 +632,10 @@ main(int, char**)
 		return 3;
 	}
 
+	if (!bodyDrag(NULL)) {
+		return 2;
+	}
+
 	cout << "bodies.cpp passed successfully" << endl;
 	return 0;
 }