ElementWaterLoad.h

#pragma once
 
#include "cable/subroutines/CableHydrodynamics.h"
#include "cable/subroutines/CableSurface.h"
 
#include <sfh/constants.h>
 
#include <Eigen/Eigen>
 
#include <cmath>
#include <fhsim_environment/EnvironmentProvider.h>
 
namespace RbCable
{
struct CylinderElement
{
   double length;     
   double radius;     
   double tangential; 
   double normal;     
   double waterDensity = 1025.0; 
};
 
struct ElementWaterLoad
{
   Eigen::Vector3d force = Eigen::Vector3d::Zero(); 
   double          displacedMass = 0.0;             
};
 
inline ElementWaterLoad WaterLoadOnElement(environment::EnvironmentProvider& environment, double T, const Eigen::Vector3d& P,
      const Eigen::Vector3d& V, const Eigen::Vector3d& axis, const CylinderElement& element)
{
   using vec3 = Eigen::Vector3d;
   ElementWaterLoad load;
 
   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))
      return load;
 
   const double volume  = sfh::pi * element.radius * element.radius * element.length;
   load.displacedMass   = element.waterDensity * volume * f;
 
   double particleVelocity[3];
   environment.GetParticleVelocity(T, P.data(), particleVelocity);
   const vec3 relative = V - vec3(particleVelocity);
   const vec3 tangent  = relative.dot(axis) * axis;
   const vec3 normal   = relative - tangent;
   load.force -= element.waterDensity * element.radius * element.length * f
         * (sfh::pi * element.tangential * tangent.norm() * tangent + element.normal * normal.norm() * normal);
 
   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(element.waterDensity, volume, f, 0.0), a.data(), axis.data(), inertiaForce);
      load.force += vec3(inertiaForce);
   }
   return load;
}
 
} // namespace RbCable