KF2012Common.h

#pragma once
 
#include "PanelFormulas.h"
#include "PanelLoadTypes.h"
#include "TwineFriction.h"
 
#include <cmath>
 
namespace hydrodynamics
{
namespace kf2012
{
 
const double kMinReynolds = 31.622776601683793;
const double kMaxReynolds = 1.0e4;
 
const double kMaxSolidity = 0.5;
const double kMinSolidity = 1.0e-4;
 
const double kMinSpeed = 1.0e-10;
 
const unsigned kFlagSolidityClamped = 1u;
const unsigned kFlagReynoldsClamped = 2u;
const unsigned kFlagSpeedBelowEpsilon = 4u;
const unsigned kFlagVelocityRatioClamped = kPanelFlagVelocityRatioClamped;
 
const double kB2 = 1.0;
 
const double kMaxA3 = 0.25;
 
const double kMinA3 = -0.125;
 
const double kMaxAbsB4 = 0.5;
 
 
using panel_formulas::CylinderDragCoefficient; 
 
 
template<class T>
T ScreenDragCoefficient(const T& cylinderCd, const T& solidity)
{
   const T open = 1.0 - solidity;
   return cylinderCd * solidity * (2.0 - solidity) / (2.0 * open * open);
}
 
 
template<class T>
T SchubauerLiftCoefficient(const T& cd)
{
   using std::sqrt;
 
   const double quarterPi  = 0.78539816339744830962;
   const T      normal     = 0.5 * cd;
   const T      tangential = quarterPi * 4.0 * normal / (8.0 + normal);
   return (normal - tangential) / sqrt(2.0);
}
 
 
using panel_formulas::VelocityReductionFactor; 
 
 
template<class T>
T InductionFactor(const T& solidity)
{
   return 0.5 * (1.0 + VelocityReductionFactor(solidity));
}
 
 
template<class T>
T LocalDragCoefficient(const T& panelCd, const T& solidity)
{
   const T open     = 1.0 - solidity;
   const T onePlusR = 1.0 + VelocityReductionFactor(solidity);
   return panelCd * 4.0 * open * open / (solidity * onePlusR * onePlusR);
}
 
 
template<class T>
T ClampAndFlag(const T& value, double lower, double upper, unsigned flag, unsigned& flags)
{
   if (!(value >= lower)) { // catches value < lower and a NaN value (MARE-0145), as TwineCrossFlow's own floor does
      flags |= flag;
      return T(lower);
   }
   if (value > upper) {
      flags |= flag;
      return T(upper);
   }
   return value;
}
 
 
template<class T>
struct Inflow
{
   T speed;        
   T direction[3]; 
   T cosTheta;     
   T theta;        
 
   T liftBasis[3];
};
 
 
template<class T>
Inflow<T> ComputeInflow(const T velocity[3], const T normal[3], unsigned& flags)
{
   using std::abs;
   using std::acos;
   using std::sqrt;
 
   Inflow<T> inflow;
   inflow.speed =
         sqrt(velocity[0] * velocity[0] + velocity[1] * velocity[1] + velocity[2] * velocity[2]);
 
   T guardedSpeed = inflow.speed;
   if (inflow.speed < kMinSpeed) {
      flags |= kFlagSpeedBelowEpsilon;
      guardedSpeed = T(kMinSpeed);
   }
 
   T normalComponent = T(0.0);
   for (int i = 0; i < 3; i++) {
      inflow.direction[i] = velocity[i] / guardedSpeed;
      normalComponent     = normalComponent + inflow.direction[i] * normal[i];
   }
 
   inflow.cosTheta = abs(normalComponent);
   if (inflow.cosTheta > 1.0) {
      inflow.cosTheta = T(1.0);
   }
   inflow.theta = acos(inflow.cosTheta);
 
   for (int i = 0; i < 3; i++) {
      inflow.liftBasis[i] =
            normalComponent * normal[i] - normalComponent * normalComponent * inflow.direction[i];
   }
   return inflow;
}
 
 
template<class T>
T DragAngleFactor(const T& cosTheta, double a3)
{
   const T cos3Theta = (4.0 * cosTheta * cosTheta - 3.0) * cosTheta;
   return (1.0 - a3) * cosTheta + a3 * cos3Theta;
}
 
 
template<class T>
T LiftAngleFactorPerSinCos(const T& cosTheta, double b4)
{
   const T cos2Theta = 2.0 * cosTheta * cosTheta - 1.0;
   return 2.0 * kB2 + 4.0 * b4 * cos2Theta;
}
 
 
template<class T>
PanelLoad<T> EvaluateScreen(double a3, double b4, bool withInduction,
      const PanelGeometry<T>& panel, const Netting<T>& netting, const Fluid& fluid,
      const FlowSample<T>& flow)
{
   using std::sqrt;
 
   PanelLoad<T> load {};
   load.flags = 0u;
 
   const T solidity =
         ClampAndFlag(netting.solidity, kMinSolidity, kMaxSolidity, kFlagSolidityClamped, load.flags);
   const Inflow<T> inflow = ComputeInflow(flow.relativeVelocity, panel.normal, load.flags);
 
   const T velocityRatio = VelocityReductionFactor(solidity, load.flags);
   const T induction     = withInduction ? 0.5 * (1.0 + velocityRatio) : T(1.0);
 
   // U_t = f |U| / (1 - Sn): MF2022 Eq. 8 (Eq. 6 when f = 1), KF2012 Eq. 5 at theta = 0.
   load.throughFlowSpeed = induction * inflow.speed / (1.0 - solidity);
   load.reynolds         = load.throughFlowSpeed * netting.twineThickness / fluid.nu;
 
   const T reynolds =
         ClampAndFlag(load.reynolds, kMinReynolds, kMaxReynolds, kFlagReynoldsClamped, load.flags);
   // f^2 scales c_d and c_l alike (lead rulings R1 and R59): MF2022 Eq. 9 normalises the lift with
   // the same factor as the drag (p. 041301-7), so c_l = f^2 Schubauer(c_d0), not Schubauer(f^2 c_d0).
   const T cdWithoutInduction = ScreenDragCoefficient(CylinderDragCoefficient(reynolds), solidity);
   const T inductionSquared   = induction * induction;
   const T cd                 = cdWithoutInduction * inductionSquared;
   const T cl                 = inductionSquared * SchubauerLiftCoefficient(cdWithoutInduction);
 
   const T dragCoefficient = cd * DragAngleFactor(inflow.cosTheta, a3);
   const T liftPerSinCos   = cl * LiftAngleFactorPerSinCos(inflow.cosTheta, b4);
   const T sinTheta        = sqrt(1.0 - inflow.cosTheta * inflow.cosTheta);
   const T liftCoefficient = liftPerSinCos * sinTheta * inflow.cosTheta;
   const T dynamicForce    = 0.5 * fluid.rho * panel.area * inflow.speed * inflow.speed;
 
   const twine_friction::InPlaneFriction<T> friction = twine_friction::EvaluateTwineFriction(
         netting, fluid, solidity, flow.relativeVelocity, panel.normal);
   const T guardedSpeed       = inflow.speed < kMinSpeed ? T(kMinSpeed) : inflow.speed;
   const T frictionDragPerQ   = friction.drag / (0.5 * fluid.rho * guardedSpeed * guardedSpeed);
   const T totalDragPerQ      = panel.area * dragCoefficient + frictionDragPerQ;
 
   for (int i = 0; i < 3; i++) {
      load.force[i] = dynamicForce * (dragCoefficient * inflow.direction[i] + liftPerSinCos * inflow.liftBasis[i])
            + friction.force[i];
   }
   load.drag            = dynamicForce * dragCoefficient + friction.drag;
   load.lift            = dynamicForce * liftCoefficient + friction.lift;
   load.theta           = inflow.theta;
   load.momentumDeficit = totalDragPerQ; // drag / (1/2 rho |U|^2); the screen part without dividing by |U|
   load.localCd         = LocalDragCoefficient(dragCoefficient + frictionDragPerQ / panel.area, solidity);
   return load;
}
 
} // namespace kf2012
} // namespace hydrodynamics