WetCylinderWater.h

#pragma once
 
#include <fhsim_environment/ElementFlowInteraction.h>
#include <fhsim_environment/EnvironmentProvider.h>
#include <fhsim_environment/hydrodynamics/BodySubmergence.h>
#include <fhsim_environment/hydrodynamics/WetCylinderLaw.h>
 
#include <sfh/constants.h>
 
#include <marenv/marenv.h>
 
#include <Eigen/Eigen>
 
#include <optional>
 
namespace wet_cylinder_water
{
struct WaterAtElement
{
   double                                         wetFraction;      
   environment::wet_cylinder::Vector              waterVelocity;    
   std::optional<environment::wet_cylinder::Vector> waveAcceleration; 
};
 
inline WaterAtElement SampleWater(const environment::wake::ElementFlowInteraction& flow, double T, const double centre[3], const double axis[3], double length, double radius, int point)
{
   environment::EnvironmentProvider& environment = *flow.Environment();
   WaterAtElement                    water;
   water.waterVelocity = {0.0, 0.0, 0.0};
   flow.WaterVelocity(T, point, centre, water.waterVelocity.data());
 
   const double half     = 0.5 * length;
   const double endA[3]  = {centre[0] - half * axis[0], centre[1] - half * axis[1], centre[2] - half * axis[2]};
   const double endB[3]  = {centre[0] + half * axis[0], centre[1] + half * axis[1], centre[2] + half * axis[2]};
   water.wetFraction     = environment::body_submergence::CylinderWetFraction(environment, environment.MaxWaveElevation(), T, endA, endB, radius);
   const auto waves = environment.GetWaves();
   if (waves && water.wetFraction > 0.0) {
      environment::wet_cylinder::Vector mean;
      environment::body_submergence::MeanWaveAcceleration(waves.get(), T, endA, endB, mean.data());
      water.waveAcceleration = mean;
   }
   return water;
}
 
struct FixedDragElement
{
   double length;                                 
   double radius;                                 
   double tangential;                             
   double normal;                                 
   double waterDensity = marenv::seawaterDensity; 
};
 
struct FixedDragLoad
{
   Eigen::Vector3d force         = Eigen::Vector3d::Zero(); 
   double          displacedMass = 0.0;                     
};
 
inline FixedDragLoad FixedDragLoadOnElement(const environment::wake::ElementFlowInteraction& flow, double T, const Eigen::Vector3d& P, const Eigen::Vector3d& V,
      const Eigen::Vector3d& axis, const FixedDragElement& element, int point)
{
   using environment::wet_cylinder::Vector;
   FixedDragLoad load;
   const auto    water = SampleWater(flow, T, P.data(), axis.data(), element.length, element.radius, point);
   const double  f     = water.wetFraction;
   if (!(f > 0.0))
      return load;
 
   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         = element.waterDensity;
   cylinder.addedMassCoefficient = 0.0;
   load.displacedMass            = element.waterDensity * cylinder.volume * f;
 
   const Vector unitAxis {axis(0), axis(1), axis(2)};
   const Vector relative {V(0) - water.waterVelocity[0], V(1) - water.waterVelocity[1], V(2) - water.waterVelocity[2]};
   const Vector drag = environment::wet_cylinder::Drag(cylinder, environment::wet_cylinder::DragCoefficients {element.normal, element.tangential}, f, relative, unitAxis);
   load.force += Eigen::Vector3d(drag.data());
   if (water.waveAcceleration)
      load.force += Eigen::Vector3d(environment::wet_cylinder::MorisonForce(cylinder, f, *water.waveAcceleration, unitAxis).data());
   return load;
}
} // namespace wet_cylinder_water