ScreenMF2022.h

#pragma once
 
#include "PanelFormulas.h"
#include "PanelLoadTypes.h"
#include "TwineFriction.h"
 
#include <cmath>
 
namespace hydrodynamics
{
 
struct ScreenMF2022
{
   static constexpr bool acceptsAtPanel = false;
 
   static constexpr double kMinValidSolidity = 0.18; 
   static constexpr double kMaxValidSolidity = 0.36; 
 
   static constexpr bool InValidSolidityRange(double solidity)
   {
      return kMinValidSolidity <= solidity && solidity <= kMaxValidSolidity;
   }
};
 
namespace screen_mf2022
{
 
constexpr double kMaxSolidity = 0.5;
 
constexpr double kMinSolidity = 0.10050014863866528;
 
constexpr double kSpeedEpsilon = 1e-9;
 
constexpr unsigned kFlagSolidityClamped      = 1u;
constexpr unsigned kFlagSpeedBelowEpsilon    = 4u;
constexpr unsigned kFlagVelocityRatioClamped = kPanelFlagVelocityRatioClamped; 
 
template<class T>
T DragCoefficientNormal(const T& solidity)
{
   return 1.782 * solidity * solidity + 1.057 * solidity - 0.053;
}
 
template<class T>
T DragCoefficient45(const T& solidity)
{
   return 1.165 * solidity - 0.0919;
}
 
template<class T>
T LiftCoefficient45(const T& solidity)
{
   return 1.693 * solidity * solidity - 0.217 * solidity + 0.022;
}
 
using panel_formulas::VelocityReductionFactor; 
 
template<class T>
T DragCoefficient(const T& cosTheta, const T& solidity)
{
   const T cd0          = DragCoefficientNormal(solidity);
   const T cd45         = DragCoefficient45(solidity);
   const T harmonic3    = 0.5 * (cd0 - std::sqrt(2.0) * cd45); // CD0 a3
   const T harmonic1    = cd0 - harmonic3;                     // CD0 (1 − a3)
   const T cosTripleArg = 4.0 * cosTheta * cosTheta * cosTheta - 3.0 * cosTheta;
   return harmonic1 * cosTheta + harmonic3 * cosTripleArg;
}
 
template<class T>
T ClampedSolidity(const T& solidity, unsigned& flags)
{
   if (solidity > kMaxSolidity) {
      flags |= kFlagSolidityClamped;
      return T(kMaxSolidity);
   }
   if (!(solidity >= kMinSolidity)) { // catches solidity < kMinSolidity and a NaN solidity (MARE-0145), as TwineCrossFlow's own floor does
      flags |= kFlagSolidityClamped;
      return T(kMinSolidity);
   }
   return solidity;
}
 
template<class T>
T SinDoubleAngle(const T& cosTheta)
{
   using std::sqrt;
   const T sinSquared = 1.0 - cosTheta * cosTheta;
   if (sinSquared > 0.0) {
      return 2.0 * cosTheta * sqrt(sinSquared);
   }
   return T(0.0);
}
 
template<class T>
T AngleFromCosine(const T& cosTheta)
{
   using std::acos;
   if (cosTheta < 1.0) {
      return acos(cosTheta);
   }
   return T(0.0);
}
 
} // namespace screen_mf2022
 
template<class T>
PanelLoad<T> Evaluate([[maybe_unused]] const ScreenMF2022& law, const PanelGeometry<T>& panel,
      const Netting<T>& netting, const Fluid& fluid, const FlowSample<T>& flow)
{
   using std::abs;
   using std::sqrt;
   namespace mf = screen_mf2022;
 
   PanelLoad<T> load {};
   load.flags = 0;
 
   const T* u        = flow.relativeVelocity;
   const T* n        = panel.normal;
   const T  speedSq  = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
   const T  solidity = mf::ClampedSolidity(netting.solidity, load.flags);
 
   T speed = T(mf::kSpeedEpsilon);
   if (speedSq > mf::kSpeedEpsilon * mf::kSpeedEpsilon) {
      speed = sqrt(speedSq);
   } else {
      load.flags |= mf::kFlagSpeedBelowEpsilon;
   }
 
   const T normalFlow = u[0] * n[0] + u[1] * n[1] + u[2] * n[2];
   const T cosSigned  = normalFlow / speed; // c = Û·n̂, either sign
   const T cosTheta   = abs(cosSigned);
 
   const T dragCoefficient = mf::DragCoefficient(cosTheta, solidity);
   const T lift45          = mf::LiftCoefficient45(solidity);
   const T pressureArea    = 0.5 * fluid.rho * speedSq * panel.area; // ½ ρ |U|² A
 
   // F = ½ρA|U|² (C_D Û + 2 CL45 (c n̂ − c² Û)).
   const T alongFlow   = pressureArea * (dragCoefficient - 2.0 * lift45 * cosSigned * cosSigned) / speed;
   const T alongNormal = pressureArea * 2.0 * lift45 * cosSigned;
   for (int i = 0; i < 3; ++i) {
      load.force[i] = alongFlow * u[i] + alongNormal * n[i];
   }
 
   const twine_friction::InPlaneFriction<T> friction = twine_friction::EvaluateTwineFriction(netting, fluid, solidity, u, n);
   for (int i = 0; i < 3; ++i) {
      load.force[i] = load.force[i] + friction.force[i];
   }
   const T frictionDragPerQ     = friction.drag / (0.5 * fluid.rho * speed * speed);
   const T totalDragCoefficient = dragCoefficient + frictionDragPerQ / panel.area;
 
   const T r         = mf::VelocityReductionFactor(solidity, load.flags);
   const T openRatio = 1.0 - solidity;
 
   load.drag             = pressureArea * dragCoefficient + friction.drag;
   load.lift             = pressureArea * lift45 * mf::SinDoubleAngle(cosTheta) + friction.lift;
   load.theta            = mf::AngleFromCosine(cosTheta);
   load.reynolds         = speed * netting.twineThickness / fluid.nu;
   load.momentumDeficit  = dragCoefficient * panel.area + frictionDragPerQ;
   load.throughFlowSpeed = speed * (1.0 + r) / (2.0 * openRatio);
   load.localCd          = totalDragCoefficient * 4.0 * openRatio * openRatio / (solidity * (1.0 + r) * (1.0 + r));
   return load;
}
 
} // namespace hydrodynamics