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 SampleWaterGivenVelocity(const environment::wake::ElementFlowInteraction& flow, double T, const double centre[3], const double axis[3], double length, double radius,
      const environment::wet_cylinder::Vector& waterVelocity)
{
   environment::EnvironmentProvider& environment = *flow.Environment();
   WaterAtElement                    water;
   water.waterVelocity = waterVelocity;
 
   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;
}
 
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::wet_cylinder::Vector waterVelocity = {0.0, 0.0, 0.0};
   flow.WaterVelocity(T, point, centre, waterVelocity.data());
   return SampleWaterGivenVelocity(flow, T, centre, axis, length, radius, waterVelocity);
}
 
struct FixedDragElement
{
   double length;                                 
   double radius;                                 
   double tangential;                             
   double normal;                                 
   double waterDensity = marenv::seawaterDensity; 
   double addedMassCoefficient = 0.0;             
};
 
struct FixedDragLoad
{
   Eigen::Vector3d force         = Eigen::Vector3d::Zero(); 
   Eigen::Vector3d moment        = Eigen::Vector3d::Zero(); 
   double          displacedMass = 0.0;                     
   double          addedMass     = 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, const Eigen::Vector3d& angularVelocity = Eigen::Vector3d::Zero(),
      const environment::wet_cylinder::EndWater* gaussWater = nullptr)
{
   using environment::wet_cylinder::Vector;
   FixedDragLoad load;
   const auto    water = gaussWater ? SampleWaterGivenVelocity(flow, T, P.data(), axis.data(), element.length, element.radius, gaussWater->Mean())
                                    : 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 = element.addedMassCoefficient;
   load.displacedMass            = element.waterDensity * cylinder.volume * f;
   load.addedMass                = environment::wet_cylinder::AddedMass(cylinder, 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]};
   if (element.addedMassCoefficient != 0.0) {
      const environment::wet_cylinder::WetCylinderState state {f, {V(0), V(1), V(2)}, water.waterVelocity, {angularVelocity(0), angularVelocity(1), angularVelocity(2)}, unitAxis, water.waveAcceleration};
      load.force += Eigen::Vector3d(environment::wet_cylinder::CoriolisForce(cylinder, state).data());
   }
   const environment::wet_cylinder::DragCoefficients model {element.normal, element.tangential};
   if (gaussWater) {
      const auto gauss = environment::wet_cylinder::GaussPointDragLoad(cylinder, model, f, {V(0), V(1), V(2)}, {angularVelocity(0), angularVelocity(1), angularVelocity(2)}, unitAxis, *gaussWater);
      load.force += Eigen::Vector3d(gauss.force.data());
      load.moment += Eigen::Vector3d(gauss.moment.data());
   } else {
      const Vector drag = environment::wet_cylinder::Drag(cylinder, model, 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