CableRMWaterLoad.h

#pragma once
 
#include "Cable.h"
#include <fhsim_environment/hydrodynamics/BodySubmergence.h>
#include <fhsim_environment/hydrodynamics/WetCylinderLaw.h>
 
#include "cable/GaussPointDrag.h"
#include "cable/WetCylinderWater.h"
 
#include <sfh/constants.h>
 
#include <Eigen/Eigen>
 
#include <cmath>
#include <fhsim_environment/EnvironmentProvider.h>
#include <fhsim_environment/ElementFlowInteraction.h>
#include <marenv/marenv.h>
 
namespace RbCable
{
struct ElementProperties
{
   double length;     
   double radius;     
   bool   kelp;       
   double kelpWeight; 
};
 
const double kWaterDensity   = marenv::seawaterDensity; 
const double kAirDensity     = 1.025;    
 
inline environment::wet_cylinder::WetCylinder CylinderOf(const ElementProperties& element)
{
   environment::wet_cylinder::WetCylinder cylinder {};
   cylinder.length       = element.length;
   cylinder.radius       = element.radius;
   cylinder.volume       = sfh::pi * element.radius * element.radius * element.length;
   cylinder.waterDensity = kWaterDensity;
   cylinder.inertia      = environment::wet_cylinder::AddedMassInertia::NoCoriolis;
   return cylinder;
}
 
inline double AddAboveSeabedLoad(const environment::wake::ElementFlowInteraction& flow, double T, const Eigen::Vector3d& P, const Eigen::Vector3d& V,
      const Eigen::Vector3d& W, const Eigen::Vector3d& axis, const ElementProperties& element, int point, Eigen::Matrix<double, 6, 1>& load,
      const wet_cylinder_water::EndWater* gaussWater = nullptr)
{
   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 auto   water = gaussWater ? wet_cylinder_water::SampleWaterGivenVelocity(flow, T, P.data(), axis.data(), element.length, element.radius, gaussWater->Mean())
                                   : wet_cylinder_water::SampleWater(flow, T, P.data(), axis.data(), element.length, element.radius, point);
   const double f     = water.wetFraction;
   if (f > 0.0) {
      const double rho = kWaterDensity;
      load(2) -= element.length * element.radius * element.radius * sfh::pi * rho * marenv::gravity * f;
 
      vec3 v = V - vec3(water.waterVelocity.data());
 
      // The wet-cylinder law (fhsim_environment): the Coriolis force of the turning added mass, on the relative velocity
      // (MARE-0314, MARE-0323), the normal drag and the skin friction along the axis (R68) by the independence principle,
      // Cd on |v_n| D / ν, Ct on |v| D / ν (ENV-0104, ENV-0105), the skin-friction moment on the spin about the axis, relative
      // to the non-rotating water (MARE-0376), and the Morison inertia force (1 + Ca) rho f V a_n normal to the axis (owner
      // ruling R104, extending R67; MARE-0314).
      using environment::wet_cylinder::Vector;
      const environment::wet_cylinder::WetCylinderState state {f, Vector {V(0), V(1), V(2)}, water.waterVelocity, Vector {W(0), W(1), W(2)}, Vector {axis(0), axis(1), axis(2)},
            water.waveAcceleration};
      const auto cylinder = CylinderOf(element);
      const auto law      = gaussWater ? environment::wet_cylinder::WaterLoadWithoutDrag(cylinder, environment::wet_cylinder::ReynoldsDrag {}, state)
                                       : environment::wet_cylinder::WaterLoad(cylinder, environment::wet_cylinder::ReynoldsDrag {}, state);
      load.segment<3>(0) += vec3(law.coriolis.data());
 
      double v_abs2 = v.squaredNorm();
      if (v_abs2 > 0) {
         vec3 vz  = v.dot(axis) * axis;
         vec3 vxy = v - vz;
 
         double v_abs   = v.norm(); /*Kelp application v2*/
         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*/
 
         // 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*/
         }
         // Sugar kelp rope application - end
         if (!gaussWater)
            load.segment<3>(0) += vec3(law.drag.data());
      }
      if (gaussWater) {
         const auto gauss = wet_cylinder_water::GaussPointDragLoad(cylinder, f, V, W, axis, *gaussWater);
         load.segment<3>(0) += gauss.force;
         load.segment<3>(3) += gauss.moment;
      }
      load.segment<3>(3) += vec3(law.torsion.data());
      if (water.waveAcceleration)
         load.segment<3>(0) += vec3(law.morison.data());
   }
   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 * marenv::gravity * (1.0 - f);
      load.segment<3>(0) -= rho * element.radius * element.length * (Cd * vxy_abs * vxy + Ct * vz_abs * vz) * (1.0 - f);
   }
   return f;
}
 
} // namespace RbCable