src/FDM/YASim/PistonEngine.cpp

   1 #include "Atmosphere.hpp"
   2 #include "Math.hpp"
   3 #include "PistonEngine.hpp"
   4 namespace yasim {
   5
   6 const static float HP2W = 745.7f;
   7 const static float CIN2CM = 1.6387064e-5f;
   8 const static float RPM2RADPS = 0.1047198f;
   9
  10 PistonEngine::PistonEngine(float power, float speed)
  11 {
  12     _boost = 1;
  13     _running = false;
  14     _fuel = true;
  15     _boostPressure = 0;
  16     _hasSuper = false;
  17
  18     _oilTemp = Atmosphere::getStdTemperature(0);
  19     _oilTempTarget = _oilTemp;
  20     _dOilTempdt = 0;
  21
  22     // Presume a BSFC (in lb/hour per HP) of 0.45.  In SI that becomes
  23     // (2.2 lb/kg, 745.7 W/hp, 3600 sec/hour) 7.62e-08 kg/Ws.
  24     _f0 = power * 7.62e-08f;
  25
  26     _power0 = power;
  27     _omega0 = speed;
  28
  29     // We must be at sea level under standard conditions
  30     _rho0 = Atmosphere::getStdDensity(0);
  31
  32     // Further presume that takeoff is (duh) full throttle and
  33     // peak-power, that means that by our efficiency function, we are
  34     // at 11/8 of "ideal" fuel flow.
  35     float realFlow = _f0 * (11.0f/8.0f);
  36     _mixCoeff = realFlow * 1.1f / _omega0;
  37
  38     _turbo = 1;
  39     _maxMP = 1e6; // No waste gate on non-turbo engines.
  40     _wastegate = 1;
  41     _charge = 1;
  42     _chargeTarget = 1;
  43     _turboLag = 2;
  44
  45     // Guess at reasonable values for these guys.  Displacements run
  46     // at about 2 cubic inches per horsepower or so, at least for
  47     // non-turbocharged engines.
  48     _compression = 8;
  49     _displacement = power * (2*CIN2CM/HP2W);
  50 }
  51
  52 void PistonEngine::setTurboParams(float turbo, float maxMP)
  53 {
  54     _turbo = turbo;
  55     _maxMP = maxMP;
  56
  57     // This changes the "sea level" manifold air density
  58     float P0 = Atmosphere::getStdPressure(0);
  59     float P = P0 * (1 + _boost * (_turbo - 1));
  60     if(P > _maxMP) P = _maxMP;
  61     float T = Atmosphere::getStdTemperature(0) * Math::pow(P/P0, 2./7.);
  62     _rho0 = P / (287.1f * T);
  63 }
  64
  65 void PistonEngine::setDisplacement(float d)
  66 {
  67     _displacement = d;
  68 }
  69
  70 void PistonEngine::setCompression(float c)
  71 {
  72     _compression = c;
  73 }
  74
  75 float PistonEngine::getMaxPower()
  76 {
  77     return _power0;
  78 }
  79
  80 bool PistonEngine::isCranking()
  81 {
  82     return _starter;
  83 }
  84
  85 float PistonEngine::getTorque()
  86 {
  87     return _torque;
  88 }
  89
  90 float PistonEngine::getFuelFlow()
  91 {
  92     return _fuelFlow;
  93 }
  94
  95 float PistonEngine::getMP()
  96 {
  97     return _mp;
  98 }
  99
 100 float PistonEngine::getEGT()
 101 {
 102     return _egt;
 103 }
 104
 105 void PistonEngine::stabilize()
 106 {
 107     _oilTemp = _oilTempTarget;
 108     _charge = _chargeTarget;
 109 }
 110
 111 void PistonEngine::integrate(float dt)
 112 {
 113     _oilTemp += (_dOilTempdt * dt);
 114
 115     // See comments in Jet.cpp for how this decay constant works
 116     float decay = 2.3f / _turboLag;
 117     _charge = (_charge + dt*decay * _chargeTarget) / (1 + dt*decay);
 118 }
 119
 120 void PistonEngine::calc(float pressure, float temp, float speed)
 121 {
 122     _running = _magnetos && _fuel && (speed > 60*RPM2RADPS);
 123
 124     // Calculate the factor required to modify supercharger output for
 125     // rpm. Assume that the normalized supercharger output ~= 1 when
 126     // the engine is at the nominal peak-power rpm.  A power equation
 127     // of the form (A * B^x * x^C) has been derived empirically from
 128     // some representative supercharger data.  This provides
 129     // near-linear output over the normal operating range, with
 130     // fall-off in the over-speed situation.
 131     float rpm_norm = (speed / _omega0);
 132     float A = 1.795206541;
 133     float B = 0.55620178;
 134     float C = 1.246708471;
 135     float rpm_factor = A * Math::pow(B, rpm_norm) * Math::pow(rpm_norm, C);
 136     _chargeTarget = 1 + (_boost * (_turbo-1) * rpm_factor);
 137
 138     if(_hasSuper) {
 139         // Superchargers have no lag
 140         _charge = _chargeTarget;
 141     } else if(!_running) {
 142         // Turbochargers only work well when the engine is actually
 143         // running.  The 25% number is a guesstimate from Vivian.
 144         _chargeTarget = 1 + (_chargeTarget - 1) * 0.25;
 145     }
 146
 147     // We need to adjust the minimum manifold pressure to get a
 148     // reasonable idle speed (a "closed" throttle doesn't suck a total
 149     // vacuum in real manifolds).  This is a hack.
 150     float _minMP = (-0.008 * _turbo ) + 0.1;
 151
 152     _mp = pressure * _charge;
 153
 154     // Scale to throttle setting, clamp to wastegate
 155     if(_running)
 156         _mp *= _minMP + (1 -_minMP) * _throttle;
 157
 158     // Scale the max MP according to the WASTEGATE control input.  Use
 159     // the un-supercharged MP as the bottom limit.
 160     float max = _wastegate * _maxMP;
 161     if(max < _mp/_charge) max = _mp/_charge;
 162     if(_mp > max) _mp = max;
 163
 164
 165     // The "boost" is the delta above ambient
 166     _boostPressure = _mp - pressure;
 167
 168     // Air entering the manifold does so rapidly, and thus the
 169     // pressure change can be assumed to be adiabatic.  Calculate a
 170     // temperature change, and use that to get the density.
 171     // Note: need to model intercoolers here...
 172     float T = temp * Math::pow(_mp/pressure, 2.0/7.0);
 173     float rho = _mp / (287.1f * T);
 174
 175     // The actual fuel flow is determined only by engine RPM and the
 176     // mixture setting.  Not all of this will burn with the same
 177     // efficiency.
 178     _fuelFlow = _mixture * speed * _mixCoeff;
 179     if(_fuel == false) _fuelFlow = 0;
 180
 181     // How much fuel could be burned with ideal (i.e. uncorrected!)
 182     // combustion.
 183     float burnable = _f0 * (rho/_rho0) * (speed/_omega0);
 184
 185     // Calculate the fuel that actually burns to produce work.  The
 186     // idea is that less than 5/8 of ideal, we get complete
 187     // combustion.  We use up all the oxygen at 1 3/8 of ideal (that
 188     // is, you need to waste fuel to use all your O2).  In between,
 189     // interpolate.  This vaguely matches a curve I copied out of a
 190     // book for a single engine.  Shrug.
 191     float burned;
 192     float r = _fuelFlow/burnable;
 193     if     (burnable == 0) burned = 0;
 194     else if(r < .625)      burned = _fuelFlow;
 195     else if(r > 1.375)     burned = burnable;
 196     else
 197         burned = _fuelFlow + (burnable-_fuelFlow)*(r-0.625f)*(4.0f/3.0f);
 198
 199     // Correct for engine control state
 200     if(!_running)
 201         burned = 0;
 202     if(_magnetos < 3)
 203         burned *= 0.9f;
 204
 205     // And finally the power is just the reference power scaled by the
 206     // amount of fuel burned, and torque is that divided by RPM.
 207     float power = _power0 * burned/_f0;
 208     _torque = power/speed;
 209
 210     // Figure that the starter motor produces 15% of the engine's
 211     // cruise torque.  Assuming 60RPM starter speed vs. 1800RPM cruise
 212     // speed on a 160HP engine, that comes out to about 160*.15/30 ==
 213     // 0.8 HP starter motor.  Which sounds about right to me.  I think
 214     // I've finally got this tuned. :)
 215     if(_starter && !_running)
 216         _torque += 0.15f * _power0/_omega0;
 217
 218     // Also, add a negative torque of 8% of cruise, to represent
 219     // internal friction.  Propeller aerodynamic friction is too low
 220     // at low RPMs to provide a good deceleration.  Interpolate it
 221     // away as we approach cruise RPMs (full at 50%, zero at 100%),
 222     // though, to prevent interaction with the power computations.
 223     // Ugly.
 224     if(speed > 0 && speed < _omega0) {
 225         float interp = 2 - 2*speed/_omega0;
 226         interp = (interp > 1) ? 1 : interp;
 227         _torque -= 0.08f * (_power0/_omega0) * interp;
 228     }
 229
 230     // Now EGT.  This one gets a little goofy.  We can calculate the
 231     // work done by an isentropically expanding exhaust gas as the
 232     // mass of the gas times the specific heat times the change in
 233     // temperature.  The mass is just the engine displacement times
 234     // the manifold density, plus the mass of the fuel, which we know.
 235     // The change in temperature can be calculated adiabatically as a
 236     // function of the exhaust gas temperature and the compression
 237     // ratio (which we know).  So just rearrange the equation to get
 238     // EGT as a function of engine power.  Cool.  I'm using a value of
 239     // 1300 J/(kg*K) for the exhaust gas specific heat.  I found this
 240     // on a web page somewhere; no idea if it's accurate.  Also,
 241     // remember that four stroke engines do one combustion cycle every
 242     // TWO revolutions, so the displacement per revolution is half of
 243     // what we'd expect.  And diddle the work done by the gas a bit to
 244     // account for non-thermodynamic losses like internal friction;
 245     // 10% should do it.
 246     float massFlow = _fuelFlow + (rho * 0.5f * _displacement * speed);
 247     float specHeat = 1300;
 248     float corr = 1.0f/(Math::pow(_compression, 0.4f) - 1.0f);
 249     _egt = corr * (power * 1.1f) / (massFlow * specHeat);
 250     if(_egt < temp) _egt = temp;
 251
 252
 253     // Oil temperature.
 254     // Assume a linear variation between ~90degC at idle and ~120degC
 255     // at full power.  No attempt to correct for airflow over the
 256     // engine is made.  Make the time constant to attain target steady-
 257     // state oil temp greater at engine off than on to reflect no
 258     // circulation.  Nothing fancy, but populates the guage with a
 259     // plausible value.
 260     float tau;  // secs
 261     if(_running) {
 262         _oilTempTarget = 363.0f + (30.0f * (power/_power0));
 263         tau = 600;
 264         // Reduce tau linearly to 300 at max power
 265         tau -= (power/_power0) * 300.0f;
 266     } else {
 267         _oilTempTarget = temp;
 268         tau = 1500;
 269     }
 270     _dOilTempdt = (_oilTempTarget - _oilTemp) / tau;
 271 }
 272
 273 }; // namespace yasim