TwineCrossFlow.h

#pragma once
 
#include "PanelFormulas.h"
#include "PanelLoadTypes.h"
 
#include <sfh/constants.h>
#include <sfh/math/math.h>
 
#include <cmath>
 
namespace hydrodynamics
{
 
struct TwineCrossFlow
{
   double ct      = 0.01; 
   double cnKnots = 0.4;  
};
 
namespace twine_cross_flow_detail
{
 
const double kMinSpeed = 1.0e-10;
const double kMaxSpeed = 1.0e10;
const double kMaxDrag = 1.0e100;
const double kMinSolidity = 1.0e-4;
const double kMaxSolidity = 0.5;
const double kMinLiftDirectionNorm = 1.0e-10;
 
const double kMinReynolds = 32.0;
const double kMaxReynolds = 1.0e4;
 
const unsigned kFlagSolidityClamped      = 1u << 0;
const unsigned kFlagReynoldsClamped      = 1u << 1;
const unsigned kFlagSpeedBelowEpsilon    = 1u << 2;
const unsigned kFlagVelocityRatioClamped = kPanelFlagVelocityRatioClamped; 
 
template<class T>
T Dot(const T a[3], const T b[3])
{
   return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
 
using panel_formulas::CylinderDragCoefficient; 
 
template<class T>
T PolynomialReynoldsNumber(const T& reynoldsNumber, unsigned& flags)
{
   if (reynoldsNumber < kMinReynolds) {
      flags |= kFlagReynoldsClamped;
      return T(kMinReynolds);
   }
   if (reynoldsNumber > kMaxReynolds) {
      flags |= kFlagReynoldsClamped;
      return T(kMaxReynolds);
   }
   return reynoldsNumber;
}
 
template<class T>
void BarFamilyDrag(const TwineCrossFlow& law, double dragFactor, const T barDirection[3],
      const T effVel[3], const T& effSpeed, const T& cylinderCd, T force[3])
{
   const T alongBar = Dot(effVel, barDirection);
   T       tangentialVel[3];
   T       normalVel[3];
   for (int i = 0; i < 3; i++) {
      tangentialVel[i] = alongBar * barDirection[i];
      normalVel[i]     = effVel[i] - tangentialVel[i];
   }
   const T normalSpeed = sfh::math::Max(T(kMinSpeed), sfh::math::Norm(normalVel, 3));
 
   // The angle of attack between the effective velocity and the bars.
   T sinAngle = normalSpeed / effSpeed;
   if (sinAngle > 1) {
      sinAngle = T(1);
   }
 
   const T cn = cylinderCd * sinAngle;
 
   const T normalDrag = sfh::math::Bound(T(dragFactor * effSpeed * cn), T(-kMaxDrag), T(kMaxDrag));
   const T tangentialDrag =
         sfh::math::Bound(T(dragFactor * effSpeed * (law.ct * sfh::pi)), T(-kMaxDrag), T(kMaxDrag));
 
   for (int i = 0; i < 3; i++) {
      force[i] = normalDrag * normalVel[i] + tangentialDrag * tangentialVel[i];
   }
}
 
template<class T>
PanelLoad<T> ZeroLoad(unsigned flags)
{
   PanelLoad<T> load {};
   for (int i = 0; i < 3; i++) {
      load.force[i] = T(0);
   }
   load.drag             = T(0);
   load.lift             = T(0);
   load.theta            = T(0);
   load.reynolds         = T(0);
   load.momentumDeficit  = T(0);
   load.throughFlowSpeed = T(0);
   load.localCd          = T(0);
   load.flags            = flags;
   return load;
}
 
} // namespace twine_cross_flow_detail
 
template<class T>
PanelLoad<T> Evaluate(const TwineCrossFlow& law, const PanelGeometry<T>& geometry,
      const Netting<T>& netting, const Fluid& fluid, const FlowSample<T>& flow)
{
   using namespace twine_cross_flow_detail;
   using std::acos;
   using std::sqrt;
 
   const T* const velocity = flow.relativeVelocity;
   const T* const normal   = geometry.normal;
 
   unsigned flags = 0;
   const T  speed = sfh::math::Norm(velocity, 3);
   if (speed < kMinSpeed) {
      flags |= kFlagSpeedBelowEpsilon;
   }
 
   T solidity = netting.solidity;
   if (!(solidity > kMinSolidity)) { // Also catches a NaN.
      solidity = T(kMinSolidity);
      flags |= kFlagSolidityClamped;
   } else if (solidity > kMaxSolidity) {
      solidity = T(kMaxSolidity);
      flags |= kFlagSolidityClamped;
   }
 
   if (!netting.hasBarDirections) {
      return ZeroLoad<T>(flags);
   }
 
   // The speed-up of Kristiansen and Faltinsen (2012), on the normal component only.
   const T speedUp     = sqrt(2 - solidity) / (sqrt(2.0) * (1 - solidity));
   const T normalSpeed = Dot(velocity, normal);
   T       effVel[3];
   for (int i = 0; i < 3; i++) {
      effVel[i] = velocity[i] + (speedUp - 1) * normalSpeed * normal[i];
   }
   const T effSpeed = sfh::math::Bound(sfh::math::Norm(effVel, 3), T(kMinSpeed), T(kMaxSpeed));
 
   // Reynolds number from the flow through the meshes and the twine diameter (MARE-0085).
   const T meshFlowSpeed  = effSpeed * sqrt(2.0) / sqrt(2 - solidity);
   const T reynoldsNumber = meshFlowSpeed * netting.twineThickness / fluid.nu;
 
   const double twineLength   = sfh::math::Max(netting.barLength - netting.knotDiameter, 0.0);
   const double barDragFactor = 0.5 * fluid.rho * netting.twineThickness * twineLength;
   const double dragFactorU   = barDragFactor * netting.numBarsU;
   const double dragFactorV   = barDragFactor * netting.numBarsV;
 
   T forceU[3];
   T forceV[3];
   const T cylinderCd = CylinderDragCoefficient(PolynomialReynoldsNumber(reynoldsNumber, flags));
   BarFamilyDrag(law, dragFactorU, netting.barU, effVel, effSpeed, cylinderCd, forceU);
   BarFamilyDrag(law, dragFactorV, netting.barV, effVel, effSpeed, cylinderCd, forceV);
 
   // Drag on the knots acts along the effective velocity, with no tangential component.
   const double knotArea       = sfh::pi * netting.knotDiameter * netting.knotDiameter / 4;
   const double knotDragFactor = 0.5 * fluid.rho * netting.numKnots * law.cnKnots * knotArea;
   const T      knotDrag       = knotDragFactor * effSpeed * effSpeed;
 
   PanelLoad<T> load {};
   for (int i = 0; i < 3; i++) {
      load.force[i] = forceU[i] + forceV[i] + knotDrag * effVel[i] / effSpeed;
   }
 
   // Drag and lift directions, and the angle between U and the normal folded into [0, pi/2].
   const T boundedSpeed = sfh::math::Max(speed, T(kMinSpeed));
   T       flowDirection[3];
   for (int i = 0; i < 3; i++) {
      flowDirection[i] = velocity[i] / boundedSpeed;
   }
   T cosTheta = Dot(flowDirection, normal);
   T normalDownstream[3];
   for (int i = 0; i < 3; i++) {
      normalDownstream[i] = (cosTheta < 0) ? -normal[i] : normal[i];
   }
   if (cosTheta < 0) {
      cosTheta = -cosTheta;
   }
   // (U-hat x n) x U-hat = n - (n . U-hat) U-hat, for a unit U-hat.
   T liftDirection[3];
   for (int i = 0; i < 3; i++) {
      liftDirection[i] = normalDownstream[i] - cosTheta * flowDirection[i];
   }
   const T liftDirectionNorm =
         sfh::math::Max(sfh::math::Norm(liftDirection, 3), T(kMinLiftDirectionNorm));
 
   load.drag  = Dot(load.force, flowDirection);
   load.lift  = Dot(load.force, liftDirection) / liftDirectionNorm;
   load.theta = acos(sfh::math::Min(cosTheta, T(1)));
 
   const T dynamicPressure = 0.5 * fluid.rho * boundedSpeed * boundedSpeed;
   load.reynolds           = reynoldsNumber;
   load.throughFlowSpeed   = meshFlowSpeed;
   load.momentumDeficit    = load.drag / dynamicPressure;
 
   // MF2022 Eq. 9, with the velocity ratio r = min(1.08 - 0.97 Sn, 1) of MF2022 Eq. 11 (MARE-0109).
   const T dragCoefficient = load.momentumDeficit / geometry.area;
   const T reductionFactor = panel_formulas::VelocityReductionFactor(solidity, flags);
   load.localCd            = dragCoefficient * 4 * (1 - solidity) * (1 - solidity) / (solidity * (1 + reductionFactor) * (1 + reductionFactor));
   load.flags              = flags;
   return load;
}
 
} // namespace hydrodynamics