CableSegmentWaterForce.h

#pragma once
 
#include "cable/GaussPointDrag.h"
#include "cable/WetCylinderWater.h"
 
#include <fhsim_environment/hydrodynamics/CylinderHydro.h>
#include <fhsim_environment/hydrodynamics/BodySubmergence.h>
#include <fhsim_environment/hydrodynamics/WetCylinderLaw.h>
 
#include <algorithm>
#include <cmath>
#include <numbers>
#include <limits>
#include <fhsim_coribo/eigen_matrix_defs.h>
#include <fhsim_environment/EnvironmentProvider.h>
#include <fhsim_environment/ElementFlowInteraction.h>
#include <marenv/marenv.h>
 
namespace trawl_cable_water
{
using CoRiBoDynamics::mat3;
using CoRiBoDynamics::vec3;
using CoRiBoDynamics::vec6;
 
struct Segment
{
   double radius;             
   double length;             
   double mass;               
   double displacementVolume; 
};
 
const double kRhoWater     = marenv::seawaterDensity; 
const double kGravity      = marenv::gravity;         
 
 
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 environment::body_submergence::CylinderWetFraction(environment, environment.MaxWaveElevation(), T, endA.data(), endB.data(), segment.radius);
}
 
inline environment::wet_cylinder::WetCylinder CylinderOf(const Segment& segment)
{
   environment::wet_cylinder::WetCylinder cylinder {};
   cylinder.length       = segment.length;
   cylinder.radius       = segment.radius;
   cylinder.volume       = segment.displacementVolume;
   cylinder.waterDensity = kRhoWater;
   cylinder.inertia      = environment::wet_cylinder::AddedMassInertia::CoriolisOnVelocity;
   return cylinder;
}
 
inline vec6 InWaterForce(const environment::wake::ElementFlowInteraction& flow, double T, const vec3& P, const vec3& V, const vec3& W, const vec3& N, const Segment& segment,
      int point, const wet_cylinder_water::EndWater* gaussWater = nullptr)
{
   using environment::wet_cylinder::Vector;
   const auto   water = gaussWater ? wet_cylinder_water::SampleWaterGivenVelocity(flow, T, P.data(), N.data(), segment.length, segment.radius, gaussWater->Mean())
                                   : wet_cylinder_water::SampleWater(flow, T, P.data(), N.data(), segment.length, segment.radius, point);
   const double f     = water.wetFraction;
   vec6         Force = vec6::Zero();
   const vec3   V_relative = V - vec3(water.waterVelocity.data());
   const double v_abs2     = V_relative.squaredNorm();
   Force(2)                = kGravity * (segment.mass - kRhoWater * segment.displacementVolume * f);
   if (f > 0.0) {
      // The wet-cylinder law: the Coriolis force of the added mass on the velocity relative to the water (MARE-0323), the
      // independence principle (Cd on the cross-flow Reynolds number, Ct on the total-speed one, ENV-0104, ENV-0105), the
      // skin-friction moment on the spin (MARE-0376) and the Morison inertia force (R104, MARE-0314).
      const environment::wet_cylinder::WetCylinderState state {f, Vector {V(0), V(1), V(2)}, water.waterVelocity, Vector {W(0), W(1), W(2)}, Vector {N(0), N(1), N(2)}, water.waveAcceleration};
      const auto cylinder = CylinderOf(segment);
      const auto load     = gaussWater ? environment::wet_cylinder::WaterLoadWithoutDrag(cylinder, environment::wet_cylinder::ReynoldsDrag {}, state)
                                       : environment::wet_cylinder::WaterLoad(cylinder, environment::wet_cylinder::ReynoldsDrag {}, state);
      Force.segment<3>(0) += vec3(load.coriolis.data());
      if (gaussWater) {
         const auto gauss = wet_cylinder_water::GaussPointDragLoad(cylinder, f, V, W, N, *gaussWater);
         Force.segment<3>(0) += gauss.force;
         Force.segment<3>(3) = gauss.moment;
      } else if (v_abs2 > 0)
         Force.segment<3>(0) += vec3(load.drag.data());
      Force.segment<3>(3) += vec3(load.torsion.data());
      if (water.waveAcceleration)
         Force.segment<3>(0) += vec3(load.morison.data());
   }
   return Force;
}
 
struct BodyState
{
   vec3 velocity;        
   vec3 angularVelocity; 
   mat3 rotation;        
};
 
struct Door
{
   double length;             
   double width;              
   double mass;               
   double displacementVolume; 
};
 
struct ClumpWeight
{
   double length; 
   double radius; 
   double mass;   
};
 
const double kDoorStallAngle = 0.66; 
 
inline double DoorDragCoefficient(double alpha)
{
   if (std::abs(alpha) > kDoorStallAngle) {
      return 1.28;
   }
   return 0.3 + std::abs(alpha) / kDoorStallAngle * 0.98;
}
 
inline double DoorLiftCoefficient(double alpha)
{
   if (alpha > kDoorStallAngle) {
      return std::tanh((0.0523 + kDoorStallAngle) * 3) * 2.5 * std::exp((kDoorStallAngle - alpha) * 3.0);
   }
   if (alpha > -0.1) {
      return std::tanh((0.0523 + alpha) * 3) * 2.5;
   }
   return std::tanh((0.0523 - 0.1) * 3) * 2.5 * std::exp((alpha + 0.1) * 3.0);
}
 
inline vec6 DoorHydrodynamicForce(const BodyState& state, const Door& door, const vec3& waterVelocity)
{
   const vec3 Nx = state.rotation.col(0);
   const vec3 Ny = state.rotation.col(1);
   const vec3 Nz = state.rotation.col(2);
   const vec3 W  = state.angularVelocity;
 
   const vec3   V     = state.velocity - waterVelocity;
   const double Vabs  = V.norm();
   const double VxAbs = V.dot(Nx);
   const vec3   Vy    = V.dot(Ny) * Ny;
   const double VzAbs = V.dot(Nz);
 
   const vec3 Vxz           = V - Vy;
   vec3       liftDirection = (Vxz.cross(Nz)).cross(Vxz);
   liftDirection /= liftDirection.norm() + std::numeric_limits<double>::epsilon();
 
   const double alpha = std::atan(-VzAbs / (VxAbs + std::numeric_limits<double>::epsilon()));
   const double Cd    = DoorDragCoefficient(alpha);
   const double Cl    = DoorLiftCoefficient(alpha);
 
   vec6         Force = vec6::Zero();
   const double Area  = door.length * door.width;
   const vec3   lift  = 0.5 * kRhoWater * Area * Cl * Vabs * Vabs * liftDirection;
   Force.segment<3>(0) += -0.5 * kRhoWater * Area * Vabs * (Cd * V) + lift;
   const vec3 Wz  = W.dot(Nz) * Nz;
   const vec3 Wxy = W - Wz;
   Force.segment<3>(3) += -0.5 * kRhoWater * Area * 0.7 * Area * W.norm() * (1.28 * Wxy + 0.1 * Wz);
   return Force;
}
 
inline vec6 DoorForce(const BodyState& state, const Door& door, const vec3& waterVelocity)
{
   vec6 Force = DoorHydrodynamicForce(state, door, waterVelocity);
   Force(2) += kGravity * (door.mass - door.displacementVolume * kRhoWater);
   return Force;
}
 
inline vec6 ClumpWeightHydrodynamicForce(const BodyState& state, const ClumpWeight& clump, const vec3& waterVelocity)
{
   const vec3 Nz = state.rotation.col(2);
   const vec3 W  = state.angularVelocity;
 
   const vec3   V    = state.velocity - waterVelocity;
   const double Vabs = V.norm();
   const vec3   Vz   = V.dot(Nz) * Nz;
   const vec3   Vxy  = V - Vz;
 
   vec6         Force = vec6::Zero();
   const double Area  = clump.length * 2 * clump.radius;
   Force.segment<3>(0) += -0.5 * kRhoWater * Area * Vabs * (1.1 * Vxy + 0.1 * Vz);
   const vec3 Wz  = W.dot(Nz) * Nz;
   const vec3 Wxy = W - Wz;
   Force.segment<3>(3) += -0.125 * kRhoWater * Area * 0.7 * Area * W.norm() * (1.28 * Wxy + 0.1 * Wz);
   return Force;
}
 
inline vec6 ClumpWeightForce(const BodyState& state, const ClumpWeight& clump, const vec3& waterVelocity)
{
   vec6 Force = ClumpWeightHydrodynamicForce(state, clump, waterVelocity);
   Force(2) += kGravity * (clump.mass - clump.length * clump.radius * clump.radius * std::numbers::pi * kRhoWater);
   return Force;
}
 
} // namespace trawl_cable_water