CableSegmentWaterForce.h

#pragma once
 
#include "cable/subroutines/CableHydrodynamics.h"
#include "cable/subroutines/CableSurface.h"
 
#include <algorithm>
#include <cmath>
#include <fhsim_coribo/eigen_matrix_defs.h>
#include <fhsim_environment/EnvironmentProvider.h>
 
namespace trawl_cable_water
{
using CoRiBoDynamics::vec3;
using CoRiBoDynamics::vec6;
 
struct Segment
{
   double radius;             
   double length;             
   double mass;               
   double displacementVolume; 
};
 
const double kRhoWater     = 1025.0;   
const double kKinViscWater = 1.005E-6; 
 
const double kAddedMassCoefficient = cable_hydrodynamics::kCylinderAddedMassCoefficient;
 
inline double WetFraction(environment::EnvironmentProvider& environment, double T, const vec3& P, const vec3& N, const Segment& segment)
{
   const vec3 endA = P - 0.5 * segment.length * N;
   const vec3 endB = P + 0.5 * segment.length * N;
   return cable_surface::CylinderWetFraction(environment, environment.MaxWaveElevation(), T, endA.data(), endB.data(), segment.radius);
}
 
inline vec6 InWaterForce(environment::EnvironmentProvider& environment, double T, const vec3& P, const vec3& V, const vec3& W, const vec3& N, const Segment& segment)
{
   vec6   Force = vec6::Zero();
   double currentVelocity[3];
   environment.GetParticleVelocity(T, P.data(), currentVelocity);
   vec3   V_relative = V - vec3(currentVelocity);
   double v_abs2     = V_relative.squaredNorm();
   const double f    = WetFraction(environment, T, P, N, segment);
   Force(2)          = 9.81 * (segment.mass - kRhoWater * segment.displacementVolume * f);
   if (v_abs2 > 0 && f > 0.0) {
      vec3   vz             = V_relative.dot(N) * N;
      vec3   vxy            = V_relative - vz;
      double vz_abs         = vz.norm();
      double vxy_abs        = vxy.norm();
      const double reynolds = cable_hydrodynamics::ReynoldsNumber(v_abs2, vxy_abs, 2.0 * segment.radius, kKinViscWater);
      const double Cd       = cable_hydrodynamics::NormalDragCoefficient(reynolds);
      const double Ct       = cable_hydrodynamics::TangentialDragCoefficient(reynolds);
      Force.segment<3>(0) += -kRhoWater * segment.radius * segment.length * (Cd * vxy_abs * vxy + Ct * vz_abs * vz) * f;
      Force.segment<3>(3) = -kRhoWater * segment.radius * segment.radius * segment.length * 0.1 * W.norm() * W.dot(N) * N * f;
   }
   const auto waves = environment.GetWaves();
   if (waves && f > 0.0) {
      const vec3 endA = P - 0.5 * segment.length * N;
      const vec3 endB = P + 0.5 * segment.length * N;
      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(kRhoWater, segment.displacementVolume, f, kAddedMassCoefficient), a.data(), N.data(), inertiaForce);
      Force.segment<3>(0) += vec3(inertiaForce);
   }
   return Force;
}
 
} // namespace trawl_cable_water