src/FDM/YASim/Propeller.cpp

   1 #include <stdio.h>
   2
   3 #include "Atmosphere.hpp"
   4 #include "Math.hpp"
   5 #include "Propeller.hpp"
   6 namespace yasim {
   7
   8 Propeller::Propeller(float radius, float v, float omega,
   9                      float rho, float power)
  10 {
  11     // Initialize numeric constants:
  12     _lambdaPeak = Math::pow(5.0, -1.0/4.0);
  13     _beta = 1.0f/(Math::pow(5.0f, -1.0f/4.0f) - Math::pow(5.0f, -5.0f/4.0f));
  14
  15     _r = radius;
  16     _etaC = 0.85f; // make this settable?
  17
  18     _j0 = v/(omega*_lambdaPeak);
  19     _baseJ0 = _j0;
  20
  21     float V2 = v*v + (_r*omega)*(_r*omega);
  22     _f0 = 2*_etaC*power/(rho*v*V2);
  23
  24     _matchTakeoff = false;
  25 }
  26
  27 void Propeller::setTakeoff(float omega0, float power0)
  28 {
  29     // Takeoff thrust coefficient at lambda==0
  30     _matchTakeoff = true;
  31     float V2 = _r*omega0 * _r*omega0;
  32     float gamma = _etaC * _beta / _j0;
  33     float torque = power0 / omega0;
  34     float density = Atmosphere::getStdDensity(0);
  35     _tc0 = (torque * gamma) / (0.5f * density * V2 * _f0);
  36 }
  37
  38 void Propeller::modPitch(float mod)
  39 {
  40     _j0 *= mod;
  41     if(_j0 < 0.25f*_baseJ0) _j0 = 0.25f*_baseJ0;
  42     if(_j0 > 4*_baseJ0)     _j0 = 4*_baseJ0;
  43 }
  44
  45
  46 void Propeller::calc(float density, float v, float omega,
  47                      float* thrustOut, float* torqueOut)
  48 {
  49     float tipspd = _r*omega;
  50     float V2 = v*v + tipspd*tipspd;
  51
  52     // Clamp v (forward velocity) to zero, now that we've used it to
  53     // calculate V (propeller "speed")
  54     if(v < 0) v = 0;
  55
  56     float J = v/omega;
  57     float lambda = J/_j0;
  58
  59     float torque = 0;
  60     if(lambda > 1) {
  61         lambda = 1.0f/lambda;
  62         torque = (density*V2*_f0*_j0)/(4*_etaC*_beta*(1-_lambdaPeak));
  63     }
  64
  65     // There's an undefined point at 1.  Just offset by a tiny bit to
  66     // fix (note: the discontinuity is at EXACTLY one, this is about
  67     // the only time in history you'll see me use == on a floating
  68     // point number!)
  69     if(lambda == 1.0) lambda = 0.9999f;
  70
  71     // Calculate lambda^4
  72     float l4 = lambda*lambda; l4 = l4*l4;
  73
  74     // thrust/torque ratio
  75     float gamma = (_etaC*_beta/_j0)*(1-l4);
  76
  77     // Compute a thrust, clamp to takeoff thrust to prevend huge
  78     // numbers at slow speeds.
  79     float tc = (1 - lambda) / (1 - _lambdaPeak);
  80     if(_matchTakeoff && tc > _tc0) tc = _tc0;
  81
  82     float thrust = 0.5f * density * V2 * _f0 * tc;
  83
  84     if(torque > 0) {
  85         torque -= thrust/gamma;
  86         thrust = -thrust;
  87     } else {
  88         torque = thrust/gamma;
  89     }
  90
  91     *thrustOut = thrust;
  92     *torqueOut = torque;
  93 }
  94
  95 }; // namespace yasim