TwineFriction.h

#pragma once
 
#include "PanelLoadTypes.h"
#include "TwineCrossFlow.h"
 
#include <sfh/constants.h>
#include <sfh/math/math.h>
 
#include <cmath>
 
namespace hydrodynamics
{
namespace twine_friction
{
 
template<class T>
struct InPlaneFriction
{
   T force[3]; 
   T drag;     
   T lift;     
};
 
 
template<class T>
void AddBarFamilyFriction(double dragFactor, double frictionCoefficient, const T barDirection[3],
      const T effVel[3], const T& effSpeed, T force[3])
{
   const T alongBar = twine_cross_flow_detail::Dot(effVel, barDirection);
   const T factor   = dragFactor * frictionCoefficient * effSpeed * alongBar;
   for (int i = 0; i < 3; i++) {
      force[i] = force[i] + factor * barDirection[i];
   }
}
 
 
template<class T>
InPlaneFriction<T> EvaluateTwineFriction(const Netting<T>& netting, const Fluid& fluid, const T& solidity,
      const T velocity[3], const T normal[3])
{
   using namespace twine_cross_flow_detail;
   using std::sqrt;
 
   InPlaneFriction<T> friction;
   for (int i = 0; i < 3; i++) {
      friction.force[i] = T(0);
   }
   friction.drag = T(0);
   friction.lift = T(0);
   if (!netting.hasBarDirections) {
      return friction;
   }
 
   // The speed-up of TwineCrossFlow, on the normal component only (R10).
   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));
 
   const double twineLength         = sfh::math::Max(netting.barLength - netting.knotDiameter, 0.0);
   const double barDragFactor       = 0.5 * fluid.rho * netting.twineThickness * twineLength;
   const double frictionCoefficient = TwineCrossFlow {}.ct * sfh::pi;
   AddBarFamilyFriction(barDragFactor * netting.numBarsU, frictionCoefficient, netting.barU, effVel, effSpeed, friction.force);
   AddBarFamilyFriction(barDragFactor * netting.numBarsV, frictionCoefficient, netting.barV, effVel, effSpeed, friction.force);
 
   const T boundedSpeed = sfh::math::Max(sfh::math::Norm(velocity, 3), 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;
   }
   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));
 
   friction.drag = Dot(friction.force, flowDirection);
   friction.lift = Dot(friction.force, liftDirection) / liftDirectionNorm;
   return friction;
}
 
} // namespace twine_friction
} // namespace hydrodynamics