CableRMWaterLoad.h

#pragma once
 
#include "Cable.h"
#include "cable/subroutines/CableSurface.h"
 
#include <sfh/constants.h>
 
#include <Eigen/Eigen>
 
#include <cmath>
#include <fhsim_environment/EnvironmentProvider.h>
 
namespace RbCable
{
struct ElementProperties
{
   double length;     
   double radius;     
   bool   kelp;       
   double kelpWeight; 
};
 
const double kWaterDensity   = 1025.0;   
const double kAirDensity     = 1.025;    
const double kWaterViscosity = 1.005E-6; 
 
inline void AddAboveSeabedLoad(environment::EnvironmentProvider& environment, double T, const Eigen::Vector3d& P, const Eigen::Vector3d& V,
      const Eigen::Vector3d& W, const Eigen::Vector3d& axis, const ElementProperties& element, Eigen::Matrix<double, 6, 1>& load)
{
   using vec3 = Eigen::Vector3d;
   // The wet fraction of the element's cylinder against the wave surface over its ends (owner
   // ruling R104, extending R66): the water part of the load takes f, the air part 1 - f. It was
   // all water or all air by the centre's side of the surface.
   const vec3   endA = P - 0.5 * element.length * axis;
   const vec3   endB = P + 0.5 * element.length * axis;
   const double f    = cable_surface::CylinderWetFraction(environment, environment.MaxWaveElevation(), T, endA.data(), endB.data(), element.radius);
   if (f > 0.0) {
      const double rho = kWaterDensity;
      load(2) -= element.length * element.radius * element.radius * sfh::pi * rho * 9.81 * f;
 
      double currentVelocity[3];
      environment.GetParticleVelocity(T, P.data(), currentVelocity);
      vec3 v = V - vec3(currentVelocity);
 
      double v_abs2 = v.squaredNorm();
      if (v_abs2 > 0) {
         vec3 vz  = v.dot(axis) * axis;
         vec3 vxy = v - vz;
 
         vec3 wz = W.dot(axis) * axis;
 
         double v_abs   = v.norm(); /*Kelp application v2*/
         double vz_abs  = vz.norm();
         double vxy_abs = vxy.norm();
 
         // v . v_xy = |v_xy|^2, so cos(theta) = |v_xy| / |v|, which tends to 0 as v_xy does:
         // 0 for flow along the axis, where the quotient below is 0/0 (MARE-0221).
         double cosTheta = vxy_abs > 0 ? v.dot(vxy) / (v_abs * vxy_abs) : 0.0; /*Kelp application v2*/
 
         const double           reynolds       = cable_hydrodynamics::ReynoldsNumber(v_abs2, vxy_abs, 2.0 * element.radius, kWaterViscosity);
         const DragCoefficients coefficients   = ElementDragCoefficients(reynolds, vxy_abs / sqrt(v_abs2));
         const double           Cd             = coefficients.normal;
         const double           CtTorsion      = coefficients.torsional;
         // Sugar kelp rope application - start
         if (element.kelp) {
            load.segment<3>(0) += -0.2 * (1 + 4 * cosTheta) * element.kelpWeight * element.length * (1.0731 * v_abs * v + 1.0509 * v) * f; /*Kelp application v2*/
         }
         const double Ct_xyz = coefficients.tangential; // skin friction along the axis (R68)
         // Sugar kelp rope application - end
         load.segment<3>(0) += -rho * element.radius * element.length * (Cd * vxy_abs * vxy + Ct_xyz * vz_abs * vz) * f;
         load.segment<3>(3) += -rho * 2 * sfh::pi * element.radius * element.length * element.radius * (CtTorsion * wz.norm() * wz) / 1000 * f;
      }
      // The Froude-Krylov force rho f V a_n normal to the axis (owner ruling R104, extending R67);
      // no added-mass term, as the element's inertia has none (MARE-0314).
      const auto waves = environment.GetWaves();
      if (waves) {
         double mean[3];
         cable_surface::MeanWaveAcceleration(waves.get(), T, endA.data(), endB.data(), mean);
         const vec3 a(mean);
         double     inertiaForce[3];
         cable_hydrodynamics::MorisonInertiaForce(
               cable_hydrodynamics::InertiaForceFactor(rho, sfh::pi * element.radius * element.radius * element.length, f, 0.0), a.data(), axis.data(), inertiaForce);
         load.segment<3>(0) += vec3(inertiaForce);
      }
   }
   if (f < 1.0) {
      const double rho     = kAirDensity;
      vec3         vz      = V.dot(axis) * axis;
      vec3         vxy     = V - vz;
      double       vz_abs  = vz.norm();
      double       vxy_abs = vxy.norm();
 
      double Cd = 1;
      double Ct = 1;
 
      load(2) -= element.length * element.radius * element.radius * sfh::pi * rho * 9.81 * (1.0 - f);
      load.segment<3>(0) -= rho * element.radius * element.length * (Cd * vxy_abs * vxy + Ct * vz_abs * vz) * (1.0 - f);
   }
}
 
} // namespace RbCable