src/FDM/JSBSim/models/propulsion/FGPiston.cpp

   1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   2
   3  Module:       FGPiston.cpp
   4  Author:       Jon S. Berndt, JSBSim framework
   5                Dave Luff, Piston engine model
   6                Ronald Jensen, Piston engine model
   7  Date started: 09/12/2000
   8  Purpose:      This module models a Piston engine
   9
  10  ------------- Copyright (C) 2000  Jon S. Berndt (jsb@hal-pc.org) --------------
  11
  12  This program is free software; you can redistribute it and/or modify it under
  13  the terms of the GNU Lesser General Public License as published by the Free Software
  14  Foundation; either version 2 of the License, or (at your option) any later
  15  version.
  16
  17  This program is distributed in the hope that it will be useful, but WITHOUT
  18  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  19  FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  20  details.
  21
  22  You should have received a copy of the GNU Lesser General Public License along with
  23  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  24  Place - Suite 330, Boston, MA  02111-1307, USA.
  25
  26  Further information about the GNU Lesser General Public License can also be found on
  27  the world wide web at http://www.gnu.org.
  28
  29 FUNCTIONAL DESCRIPTION
  30 --------------------------------------------------------------------------------
  31
  32 This class descends from the FGEngine class and models a Piston engine based on
  33 parameters given in the engine config file for this class
  34
  35 HISTORY
  36 --------------------------------------------------------------------------------
  37 09/12/2000  JSB  Created
  38
  39 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  40 INCLUDES
  41 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  42
  43 #include <sstream>
  44
  45 #include "FGPiston.h"
  46 #include <models/FGPropulsion.h>
  47 #include "FGPropeller.h"
  48
  49 namespace JSBSim {
  50
  51 static const char *IdSrc = "$Id$";
  52 static const char *IdHdr = ID_PISTON;
  53
  54 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  55 CLASS IMPLEMENTATION
  56 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  57
  58 FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number)
  59   : FGEngine(exec, el, engine_number),
  60   R_air(287.3),                  // Gas constant for air J/Kg/K
  61   rho_fuel(800),                 // estimate
  62   calorific_value_fuel(47.3e6),
  63   Cp_air(1005),                  // Specific heat (constant pressure) J/Kg/K
  64   Cp_fuel(1700)
  65 {
  66   string token;
  67
  68   // Defaults and initializations
  69
  70   Type = etPiston;
  71   dt = State->Getdt();
  72
  73   // These items are read from the configuration file
  74
  75   Cycles = 2;
  76   IdleRPM = 600;
  77   MaxRPM = 2800;
  78   Displacement = 360;
  79   SparkFailDrop = 1.0;
  80   MaxHP = 200;
  81   MinManifoldPressure_inHg = 6.5;
  82   MaxManifoldPressure_inHg = 28.5;
  83   BSFC = 0.45;
  84
  85   // These are internal program variables
  86
  87   crank_counter = 0;
  88   Magnetos = 0;
  89   minMAP = 21950;
  90   maxMAP = 96250;
  91
  92   ResetToIC();
  93
  94   // Supercharging
  95   BoostSpeeds = 0;  // Default to no supercharging
  96   BoostSpeed = 0;
  97   Boosted = false;
  98   BoostOverride = 0;
  99   bBoostOverride = false;
 100   bTakeoffBoost = false;
 101   TakeoffBoost = 0.0;   // Default to no extra takeoff-boost
 102   int i;
 103   for (i=0; i<FG_MAX_BOOST_SPEEDS; i++) {
 104     RatedBoost[i] = 0.0;
 105     RatedPower[i] = 0.0;
 106     RatedAltitude[i] = 0.0;
 107     BoostMul[i] = 1.0;
 108     RatedMAP[i] = 100000;
 109     RatedRPM[i] = 2500;
 110     TakeoffMAP[i] = 100000;
 111   }
 112   for (i=0; i<FG_MAX_BOOST_SPEEDS-1; i++) {
 113     BoostSwitchAltitude[i] = 0.0;
 114     BoostSwitchPressure[i] = 0.0;
 115   }
 116   // Initialisation
 117   volumetric_efficiency = 0.8;  // Actually f(speed, load) but this will get us running
 118
 119   // First column is thi, second is neta (combustion efficiency)
 120   Lookup_Combustion_Efficiency = new FGTable(12);
 121   *Lookup_Combustion_Efficiency << 0.00 << 0.980;
 122   *Lookup_Combustion_Efficiency << 0.90 << 0.980;
 123   *Lookup_Combustion_Efficiency << 1.00 << 0.970;
 124   *Lookup_Combustion_Efficiency << 1.05 << 0.950;
 125   *Lookup_Combustion_Efficiency << 1.10 << 0.900;
 126   *Lookup_Combustion_Efficiency << 1.15 << 0.850;
 127   *Lookup_Combustion_Efficiency << 1.20 << 0.790;
 128   *Lookup_Combustion_Efficiency << 1.30 << 0.700;
 129   *Lookup_Combustion_Efficiency << 1.40 << 0.630;
 130   *Lookup_Combustion_Efficiency << 1.50 << 0.570;
 131   *Lookup_Combustion_Efficiency << 1.60 << 0.525;
 132   *Lookup_Combustion_Efficiency << 2.00 << 0.345;
 133
 134   Power_Mixture_Correlation = new FGTable(13);
 135   *Power_Mixture_Correlation << (14.7/1.6) << 0.780;
 136   *Power_Mixture_Correlation << 10 <<  0.860;
 137   *Power_Mixture_Correlation << 11 <<  0.935;
 138   *Power_Mixture_Correlation << 12 <<  0.980;
 139   *Power_Mixture_Correlation << 13 <<  1.000;
 140   *Power_Mixture_Correlation << 14 <<  0.990;
 141   *Power_Mixture_Correlation << 15 <<  0.964;
 142   *Power_Mixture_Correlation << 16 <<  0.925;
 143   *Power_Mixture_Correlation << 17 <<  0.880;
 144   *Power_Mixture_Correlation << 18 <<  0.830;
 145   *Power_Mixture_Correlation << 19 <<  0.785;
 146   *Power_Mixture_Correlation << 20 <<  0.740;
 147   *Power_Mixture_Correlation << (14.7/0.6) << 0.58;
 148
 149   Mixture_Efficiency_Correlation = new FGTable(15);
 150   *Mixture_Efficiency_Correlation << 0.05000 << 0.00000;
 151   *Mixture_Efficiency_Correlation << 0.05137 << 0.00862;
 152   *Mixture_Efficiency_Correlation << 0.05179 << 0.21552;
 153   *Mixture_Efficiency_Correlation << 0.05430 << 0.48276;
 154   *Mixture_Efficiency_Correlation << 0.05842 << 0.70690;
 155   *Mixture_Efficiency_Correlation << 0.06312 << 0.83621;
 156   *Mixture_Efficiency_Correlation << 0.06942 << 0.93103;
 157   *Mixture_Efficiency_Correlation << 0.07786 << 1.00000;
 158   *Mixture_Efficiency_Correlation << 0.08845 << 1.00000;
 159   *Mixture_Efficiency_Correlation << 0.09270 << 0.98276;
 160   *Mixture_Efficiency_Correlation << 0.10120 << 0.93103;
 161   *Mixture_Efficiency_Correlation << 0.11455 << 0.72414;
 162   *Mixture_Efficiency_Correlation << 0.12158 << 0.45690;
 163   *Mixture_Efficiency_Correlation << 0.12435 << 0.23276;
 164   *Mixture_Efficiency_Correlation << 0.12500 << 0.00000;
 165
 166
 167 /*
 168 Manifold_Pressure_Lookup = new
 169
 170         0       0.2      0.4      0.6      0.8     1
 171 0        1.0000    1.0000    1.0000    1.0000    1.0000    1.0000
 172 1000    0.7778    0.8212    0.8647    0.9081    0.9516    0.9950
 173 2000    0.5556    0.6424    0.7293    0.8162    0.9031    0.9900
 174 3000    0.3333    0.4637    0.5940    0.7243    0.8547    0.9850
 175 4000    0.2000    0.2849    0.4587    0.6324    0.8062    0.9800
 176 5000    0.2000    0.2000    0.3233    0.5406    0.7578    0.9750
 177 6000    0.2000    0.2000    0.2000    0.4487    0.7093    0.9700
 178 7000    0.2000    0.2000    0.2000    0.2000    0.4570    0.7611
 179 8000    0.2000    0.2000    0.2000    0.2000    0.2047    0.5522
 180 */
 181
 182   // Read inputs from engine data file where present.
 183
 184   if (el->FindElement("minmp")) // Should have ELSE statement telling default value used?
 185     MinManifoldPressure_inHg = el->FindElementValueAsNumberConvertTo("minmp","INHG");
 186   if (el->FindElement("maxmp"))
 187     MaxManifoldPressure_inHg = el->FindElementValueAsNumberConvertTo("maxmp","INHG");
 188   if (el->FindElement("displacement"))
 189     Displacement = el->FindElementValueAsNumberConvertTo("displacement","IN3");
 190   if (el->FindElement("maxhp"))
 191     MaxHP = el->FindElementValueAsNumberConvertTo("maxhp","HP");
 192   if (el->FindElement("sparkfaildrop"))
 193     SparkFailDrop = Constrain(0, 1 - el->FindElementValueAsNumber("sparkfaildrop"), 1);
 194   if (el->FindElement("cycles"))
 195     Cycles = el->FindElementValueAsNumber("cycles");
 196   if (el->FindElement("idlerpm"))
 197     IdleRPM = el->FindElementValueAsNumber("idlerpm");
 198   if (el->FindElement("maxrpm"))
 199     MaxRPM = el->FindElementValueAsNumber("maxrpm");
 200   if (el->FindElement("maxthrottle"))
 201     MaxThrottle = el->FindElementValueAsNumber("maxthrottle");
 202   if (el->FindElement("minthrottle"))
 203     MinThrottle = el->FindElementValueAsNumber("minthrottle");
 204   if (el->FindElement("bsfc"))
 205     BSFC = el->FindElementValueAsNumber("bsfc");
 206   if (el->FindElement("numboostspeeds")) { // Turbo- and super-charging parameters
 207     BoostSpeeds = (int)el->FindElementValueAsNumber("numboostspeeds");
 208     if (el->FindElement("boostoverride"))
 209       BoostOverride = (int)el->FindElementValueAsNumber("boostoverride");
 210     if (el->FindElement("takeoffboost"))
 211       TakeoffBoost = el->FindElementValueAsNumberConvertTo("takeoffboost", "PSI");
 212     if (el->FindElement("ratedboost1"))
 213       RatedBoost[0] = el->FindElementValueAsNumberConvertTo("ratedboost1", "PSI");
 214     if (el->FindElement("ratedboost2"))
 215       RatedBoost[1] = el->FindElementValueAsNumberConvertTo("ratedboost2", "PSI");
 216     if (el->FindElement("ratedboost3"))
 217       RatedBoost[2] = el->FindElementValueAsNumberConvertTo("ratedboost3", "PSI");
 218     if (el->FindElement("ratedpower1"))
 219       RatedPower[0] = el->FindElementValueAsNumberConvertTo("ratedpower1", "HP");
 220     if (el->FindElement("ratedpower2"))
 221       RatedPower[1] = el->FindElementValueAsNumberConvertTo("ratedpower2", "HP");
 222     if (el->FindElement("ratedpower3"))
 223       RatedPower[2] = el->FindElementValueAsNumberConvertTo("ratedpower3", "HP");
 224     if (el->FindElement("ratedrpm1"))
 225       RatedRPM[0] = el->FindElementValueAsNumber("ratedrpm1");
 226     if (el->FindElement("ratedrpm2"))
 227       RatedRPM[1] = el->FindElementValueAsNumber("ratedrpm2");
 228     if (el->FindElement("ratedrpm3"))
 229       RatedRPM[2] = el->FindElementValueAsNumber("ratedrpm3");
 230     if (el->FindElement("ratedaltitude1"))
 231       RatedAltitude[0] = el->FindElementValueAsNumberConvertTo("ratedaltitude1", "FT");
 232     if (el->FindElement("ratedaltitude2"))
 233       RatedAltitude[1] = el->FindElementValueAsNumberConvertTo("ratedaltitude2", "FT");
 234     if (el->FindElement("ratedaltitude3"))
 235       RatedAltitude[2] = el->FindElementValueAsNumberConvertTo("ratedaltitude3", "FT");
 236   }
 237   char property_name[80];
 238   snprintf(property_name, 80, "propulsion/engine[%d]/power_hp", EngineNumber);
 239   PropertyManager->Tie(property_name, &HP);
 240   snprintf(property_name, 80, "propulsion/engine[%d]/bsfc", EngineNumber);
 241   PropertyManager->Tie(property_name, &BSFC);
 242
 243   minMAP = MinManifoldPressure_inHg * inhgtopa;  // inHg to Pa
 244   maxMAP = MaxManifoldPressure_inHg * inhgtopa;
 245   StarterHP = sqrt(MaxHP) * 0.4;
 246
 247   // Set up and sanity-check the turbo/supercharging configuration based on the input values.
 248   if (TakeoffBoost > RatedBoost[0]) bTakeoffBoost = true;
 249   for (i=0; i<BoostSpeeds; ++i) {
 250     bool bad = false;
 251     if (RatedBoost[i] <= 0.0) bad = true;
 252     if (RatedPower[i] <= 0.0) bad = true;
 253     if (RatedAltitude[i] < 0.0) bad = true;  // 0.0 is deliberately allowed - this corresponds to unregulated supercharging.
 254     if (i > 0 && RatedAltitude[i] < RatedAltitude[i - 1]) bad = true;
 255     if (bad) {
 256       // We can't recover from the above - don't use this supercharger speed.
 257       BoostSpeeds--;
 258       // TODO - put out a massive error message!
 259       break;
 260     }
 261     // Now sanity-check stuff that is recoverable.
 262     if (i < BoostSpeeds - 1) {
 263       if (BoostSwitchAltitude[i] < RatedAltitude[i]) {
 264         // TODO - put out an error message
 265         // But we can also make a reasonable estimate, as below.
 266         BoostSwitchAltitude[i] = RatedAltitude[i] + 1000;
 267       }
 268       BoostSwitchPressure[i] = Atmosphere->GetPressure(BoostSwitchAltitude[i]) * psftopa;
 269       //cout << "BoostSwitchAlt = " << BoostSwitchAltitude[i] << ", pressure = " << BoostSwitchPressure[i] << '\n';
 270       // Assume there is some hysteresis on the supercharger gear switch, and guess the value for now
 271       BoostSwitchHysteresis = 1000;
 272     }
 273     // Now work out the supercharger pressure multiplier of this speed from the rated boost and altitude.
 274     RatedMAP[i] = Atmosphere->GetPressureSL() * psftopa + RatedBoost[i] * 6895;  // psi*6895 = Pa.
 275     // Sometimes a separate BCV setting for takeoff or extra power is fitted.
 276     if (TakeoffBoost > RatedBoost[0]) {
 277       // Assume that the effect on the BCV is the same whichever speed is in use.
 278       TakeoffMAP[i] = RatedMAP[i] + ((TakeoffBoost - RatedBoost[0]) * 6895);
 279       bTakeoffBoost = true;
 280     } else {
 281       TakeoffMAP[i] = RatedMAP[i];
 282       bTakeoffBoost = false;
 283     }
 284     BoostMul[i] = RatedMAP[i] / (Atmosphere->GetPressure(RatedAltitude[i]) * psftopa);
 285
 286   }
 287
 288   if (BoostSpeeds > 0) {
 289     Boosted = true;
 290     BoostSpeed = 0;
 291   }
 292   bBoostOverride = (BoostOverride == 1 ? true : false);
 293   if (MinThrottle < 0.001) MinThrottle = 0.001;  //MinThrottle is a denominator in a power equation so it can't be zero
 294   Debug(0); // Call Debug() routine from constructor if needed
 295 }
 296
 297 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 298
 299 FGPiston::~FGPiston()
 300 {
 301   char property_name[80];
 302
 303   delete Lookup_Combustion_Efficiency;
 304   delete Power_Mixture_Correlation;
 305   delete Mixture_Efficiency_Correlation;
 306   Debug(1); // Call Debug() routine from constructor if needed
 307 }
 308
 309 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 310
 311 void FGPiston::ResetToIC(void)
 312 {
 313   FGEngine::ResetToIC();
 314
 315   ManifoldPressure_inHg = Atmosphere->GetPressure() * psftoinhg; // psf to in Hg
 316   MAP = Atmosphere->GetPressure() * psftopa;
 317   double airTemperature_degK = RankineToKelvin(Atmosphere->GetTemperature());
 318   OilTemp_degK = airTemperature_degK;
 319   CylinderHeadTemp_degK = airTemperature_degK;
 320   ExhaustGasTemp_degK = airTemperature_degK;
 321   EGT_degC = ExhaustGasTemp_degK - 273;
 322   Thruster->SetRPM(0.0);
 323   RPM = 0.0;
 324 }
 325
 326 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 327
 328 double FGPiston::Calculate(void)
 329 {
 330   if (FuelFlow_gph > 0.0) ConsumeFuel();
 331
 332   Throttle = FCS->GetThrottlePos(EngineNumber);
 333   ThrottlePos = MinThrottle+((MaxThrottle-MinThrottle)*Throttle );
 334   Mixture = FCS->GetMixturePos(EngineNumber);
 335
 336   //
 337   // Input values.
 338   //
 339
 340   p_amb = Atmosphere->GetPressure() * psftopa;
 341   p_amb_sea_level = Atmosphere->GetPressureSL() * psftopa;
 342   T_amb = RankineToKelvin(Atmosphere->GetTemperature());
 343
 344   RPM = Thruster->GetRPM() * Thruster->GetGearRatio();
 345
 346   IAS = Auxiliary->GetVcalibratedKTS();
 347
 348   doEngineStartup();
 349   if (Boosted) doBoostControl();
 350   doMAP();
 351   doAirFlow();
 352   doFuelFlow();
 353
 354   //Now that the fuel flow is done check if the mixture is too lean to run the engine
 355   //Assume lean limit at 22 AFR for now - thats a thi of 0.668
 356   //This might be a bit generous, but since there's currently no audiable warning of impending
 357   //cutout in the form of misfiring and/or rough running its probably reasonable for now.
 358 //  if (equivalence_ratio < 0.668)
 359 //    Running = false;
 360
 361   doEnginePower();
 362 if(HP<0.1250)
 363   Running = false;
 364
 365   doEGT();
 366   doCHT();
 367   doOilTemperature();
 368   doOilPressure();
 369
 370   if (Thruster->GetType() == FGThruster::ttPropeller) {
 371     ((FGPropeller*)Thruster)->SetAdvance(FCS->GetPropAdvance(EngineNumber));
 372     ((FGPropeller*)Thruster)->SetFeather(FCS->GetPropFeather(EngineNumber));
 373   }
 374
 375   PowerAvailable = (HP * hptoftlbssec) - Thruster->GetPowerRequired();
 376
 377   return Thruster->Calculate(PowerAvailable);
 378 }
 379
 380 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 381
 382 double FGPiston::CalcFuelNeed(void)
 383 {
 384   return FuelFlow_gph / 3600 * 6 * State->Getdt() * Propulsion->GetRate();
 385 }
 386
 387 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 388
 389 int FGPiston::InitRunning(void) {
 390   Magnetos=3;
 391   //Thruster->SetRPM( 1.1*IdleRPM/Thruster->GetGearRatio() );
 392   Thruster->SetRPM( 1000 );
 393   Running=true;
 394   return 1;
 395 }
 396
 397 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 398 /**
 399  * Start or stop the engine.
 400  */
 401
 402 void FGPiston::doEngineStartup(void)
 403 {
 404   // Check parameters that may alter the operating state of the engine.
 405   // (spark, fuel, starter motor etc)
 406   bool spark;
 407   bool fuel;
 408
 409   // Check for spark
 410   Magneto_Left = false;
 411   Magneto_Right = false;
 412   // Magneto positions:
 413   // 0 -> off
 414   // 1 -> left only
 415   // 2 -> right only
 416   // 3 -> both
 417   if (Magnetos != 0) {
 418     spark = true;
 419   } else {
 420     spark = false;
 421   }  // neglects battery voltage, master on switch, etc for now.
 422
 423   if ((Magnetos == 1) || (Magnetos > 2)) Magneto_Left = true;
 424   if (Magnetos > 1)  Magneto_Right = true;
 425
 426   // Assume we have fuel for now
 427   fuel = !Starved;
 428
 429   // Check if we are turning the starter motor
 430   if (Cranking != Starter) {
 431     // This check saves .../cranking from getting updated every loop - they
 432     // only update when changed.
 433     Cranking = Starter;
 434     crank_counter = 0;
 435   }
 436
 437   if (Cranking) crank_counter++;  //Check mode of engine operation
 438
 439   if (!Running && spark && fuel) {  // start the engine if revs high enough
 440     if (Cranking) {
 441       if ((RPM > IdleRPM*0.8) && (crank_counter > 175)) // Add a little delay to startup
 442         Running = true;                         // on the starter
 443     } else {
 444       if (RPM > IdleRPM*0.8)                            // This allows us to in-air start
 445         Running = true;                         // when windmilling
 446     }
 447   }
 448
 449   // Cut the engine *power* - Note: the engine may continue to
 450   // spin if the prop is in a moving airstream
 451
 452   if ( Running && (!spark || !fuel) ) Running = false;
 453
 454   // Check for stalling (RPM = 0).
 455   if (Running) {
 456     if (RPM == 0) {
 457       Running = false;
 458     } else if ((RPM <= IdleRPM *0.8 ) && (Cranking)) {
 459       Running = false;
 460     }
 461   }
 462 }
 463
 464 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 465
 466 /**
 467  * Calculate the Current Boost Speed
 468  *
 469  * This function calculates the current turbo/supercharger boost speed
 470  * based on altitude and the (automatic) boost-speed control valve configuration.
 471  *
 472  * Inputs: p_amb, BoostSwitchPressure, BoostSwitchHysteresis
 473  *
 474  * Outputs: BoostSpeed
 475  */
 476
 477 void FGPiston::doBoostControl(void)
 478 {
 479   if(BoostSpeed < BoostSpeeds - 1) {
 480     // Check if we need to change to a higher boost speed
 481     if(p_amb < BoostSwitchPressure[BoostSpeed] - BoostSwitchHysteresis) {
 482       BoostSpeed++;
 483     }
 484   } else if(BoostSpeed > 0) {
 485     // Check if we need to change to a lower boost speed
 486     if(p_amb > BoostSwitchPressure[BoostSpeed - 1] + BoostSwitchHysteresis) {
 487       BoostSpeed--;
 488     }
 489   }
 490 }
 491
 492 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 493
 494 /**
 495  * Calculate the manifold absolute pressure (MAP) in inches hg
 496  *
 497  * This function calculates manifold absolute pressure (MAP)
 498  * from the throttle position, turbo/supercharger boost control
 499  * system, engine speed and local ambient air density.
 500  *
 501  * TODO: changes in MP should not be instantaneous -- introduce
 502  * a lag between throttle changes and MP changes, to allow pressure
 503  * to build up or disperse.
 504  *
 505  * Inputs: minMAP, maxMAP, p_amb, Throttle
 506  *
 507  * Outputs: MAP, ManifoldPressure_inHg
 508  */
 509
 510 void FGPiston::doMAP(void)
 511 {
 512     suction_loss = RPM > 0.0 ? ThrottlePos * MaxRPM / RPM : 1.0;
 513     if (suction_loss > 1.0) suction_loss = 1.0;
 514     MAP = p_amb * suction_loss;
 515
 516     if(Boosted) {
 517       // If takeoff boost is fitted, we currently assume the following throttle map:
 518       // (In throttle % - actual input is 0 -> 1)
 519       // 99 / 100 - Takeoff boost
 520       // 96 / 97 / 98 - Rated boost
 521       // 0 - 95 - Idle to Rated boost (MinManifoldPressure to MaxManifoldPressure)
 522       // In real life, most planes would be fitted with a mechanical 'gate' between
 523       // the rated boost and takeoff boost positions.
 524       double T = Throttle; // processed throttle value.
 525       bool bTakeoffPos = false;
 526       if(bTakeoffBoost) {
 527         if(Throttle > 0.98) {
 528           //cout << "Takeoff Boost!!!!\n";
 529           bTakeoffPos = true;
 530         } else if(Throttle <= 0.95) {
 531           bTakeoffPos = false;
 532           T *= 1.0 / 0.95;
 533         } else {
 534           bTakeoffPos = false;
 535           //cout << "Rated Boost!!\n";
 536           T = 1.0;
 537         }
 538       }
 539       // Boost the manifold pressure.
 540       MAP += MAP * BoostMul[BoostSpeed] * suction_loss * RPM/RatedRPM[BoostSpeed];
 541       // Now clip the manifold pressure to BCV or Wastegate setting.
 542       if(bTakeoffPos) {
 543         if(MAP > TakeoffMAP[BoostSpeed]) {
 544           MAP = TakeoffMAP[BoostSpeed];
 545         }
 546       } else {
 547         if(MAP > RatedMAP[BoostSpeed]) {
 548           MAP = RatedMAP[BoostSpeed];
 549         }
 550       }
 551     }
 552
 553   // And set the value in American units as well
 554   ManifoldPressure_inHg = MAP / inhgtopa;
 555 }
 556
 557 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 558 /**
 559  * Calculate the air flow through the engine.
 560  * Also calculates ambient air density
 561  * (used in CHT calculation for air-cooled engines).
 562  *
 563  * Inputs: p_amb, R_air, T_amb, MAP, Displacement,
 564  *   RPM, volumetric_efficiency, ThrottlePos
 565  *
 566  * TODO: Model inlet manifold air temperature.
 567  *
 568  * Outputs: rho_air, m_dot_air
 569  */
 570
 571 void FGPiston::doAirFlow(void)
 572 {
 573
 574 rho_air = p_amb / (R_air * T_amb);
 575   double displacement_SI = Displacement * in3tom3;
 576   double swept_volume = (displacement_SI * (RPM/60)) / 2;
 577   double v_dot_air = swept_volume * volumetric_efficiency * suction_loss;
 578
 579   double rho_air_manifold = MAP / (R_air * T_amb);
 580   m_dot_air = v_dot_air * rho_air_manifold;
 581 }
 582
 583 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 584 /**
 585  * Calculate the fuel flow into the engine.
 586  *
 587  * Inputs: Mixture, thi_sea_level, p_amb_sea_level, p_amb, m_dot_air
 588  *
 589  * Outputs: equivalence_ratio, m_dot_fuel
 590  */
 591
 592 void FGPiston::doFuelFlow(void)
 593 {
 594   double thi_sea_level = 1.3 * Mixture; // Allows an AFR of infinity:1 to 11.3075:1
 595   equivalence_ratio = thi_sea_level; // * p_amb_sea_level / p_amb;
 596   double AFR = 10+(12*(1-Mixture));// mixture 10:1 to 22:1
 597   m_dot_fuel = m_dot_air / AFR;
 598   FuelFlow_gph = m_dot_fuel
 599     * 3600            // seconds to hours
 600     * 2.2046            // kg to lb
 601     / 6.0;            // lb to gal_us of gasoline
 602 //    / 6.6;            // lb to gal_us of kerosene
 603 }
 604
 605 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 606 /**
 607  * Calculate the power produced by the engine.
 608  *
 609  * Currently, the JSBSim propellor model does not allow the
 610  * engine to produce enough RPMs to get up to a high horsepower.
 611  * When tested with sufficient RPM, it has no trouble reaching
 612  * 200HP.
 613  *
 614  * Inputs: ManifoldPressure_inHg, p_amb, RPM, T_amb,
 615  *   Mixture_Efficiency_Correlation, Cycles, MaxHP
 616  *
 617  * Outputs: Percentage_Power, HP
 618  */
 619
 620 void FGPiston::doEnginePower(void)
 621 {
 622   if (Running) {
 623     double T_amb_degF = KelvinToFahrenheit(T_amb);
 624     double T_amb_sea_lev_degF = KelvinToFahrenheit(288);
 625
 626     // FIXME: this needs to be generalized
 627     double ME, friction, percent_RPM;  // Convienience term for use in the calculations
 628     ME = Mixture_Efficiency_Correlation->GetValue(m_dot_fuel/m_dot_air);
 629
 630     percent_RPM = RPM/MaxRPM;
 631     friction = 1 - (percent_RPM * percent_RPM * percent_RPM * percent_RPM/10);
 632     if (friction < 0 ) friction = 0;
 633     Percentage_Power = friction;
 634
 635     if ( Magnetos != 3 ) Percentage_Power *= SparkFailDrop;
 636
 637
 638     HP = (FuelFlow_gph * 6.0 / BSFC )* ME * suction_loss * Percentage_Power;
 639
 640   } else {
 641
 642     // Power output when the engine is not running
 643     if (Cranking) {
 644       if (RPM < 10) {
 645         HP = StarterHP;
 646       } else if (RPM < IdleRPM*0.8) {
 647         HP = StarterHP + ((IdleRPM*0.8 - RPM) / 8.0);
 648         // This is a guess - would be nice to find a proper starter moter torque curve
 649       } else {
 650         HP = StarterHP;
 651       }
 652     } else {
 653       // Quick hack until we port the FMEP stuff
 654       if (RPM > 0.0)
 655         HP = -1.5;
 656       else
 657         HP = 0.0;
 658     }
 659   }
 660 //  cout << "Power = " << HP << "  RPM = " << RPM << "  Running = " << Running << "  Cranking = " << Cranking << endl;
 661 }
 662
 663 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 664 /**
 665  * Calculate the exhaust gas temperature.
 666  *
 667  * Inputs: equivalence_ratio, m_dot_fuel, calorific_value_fuel,
 668  *   Cp_air, m_dot_air, Cp_fuel, m_dot_fuel, T_amb, Percentage_Power
 669  *
 670  * Outputs: combustion_efficiency, ExhaustGasTemp_degK
 671  */
 672
 673 void FGPiston::doEGT(void)
 674 {
 675   double delta_T_exhaust;
 676   double enthalpy_exhaust;
 677   double heat_capacity_exhaust;
 678   double dEGTdt;
 679
 680   if ((Running) && (m_dot_air > 0.0)) {  // do the energy balance
 681     combustion_efficiency = Lookup_Combustion_Efficiency->GetValue(equivalence_ratio);
 682     enthalpy_exhaust = m_dot_fuel * calorific_value_fuel *
 683                               combustion_efficiency * 0.33;
 684     heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
 685     delta_T_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
 686     ExhaustGasTemp_degK = T_amb + delta_T_exhaust;
 687     ExhaustGasTemp_degK *= 0.444 + ((0.544 - 0.444) * Percentage_Power);
 688   } else {  // Drop towards ambient - guess an appropriate time constant for now
 689     combustion_efficiency = 0;
 690     dEGTdt = (RankineToKelvin(Atmosphere->GetTemperature()) - ExhaustGasTemp_degK) / 100.0;
 691     delta_T_exhaust = dEGTdt * dt;
 692     ExhaustGasTemp_degK += delta_T_exhaust;
 693   }
 694 }
 695
 696 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 697 /**
 698  * Calculate the cylinder head temperature.
 699  *
 700  * Inputs: T_amb, IAS, rho_air, m_dot_fuel, calorific_value_fuel,
 701  *   combustion_efficiency, RPM
 702  *
 703  * Outputs: CylinderHeadTemp_degK
 704  */
 705
 706 void FGPiston::doCHT(void)
 707 {
 708   double h1 = -95.0;
 709   double h2 = -3.95;
 710   double h3 = -0.05;
 711
 712   double arbitary_area = 1.0;
 713   double CpCylinderHead = 800.0;
 714   double MassCylinderHead = 8.0;
 715
 716   double temperature_difference = CylinderHeadTemp_degK - T_amb;
 717   double v_apparent = IAS * 0.5144444;
 718   double v_dot_cooling_air = arbitary_area * v_apparent;
 719   double m_dot_cooling_air = v_dot_cooling_air * rho_air;
 720   double dqdt_from_combustion =
 721     m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
 722   double dqdt_forced = (h2 * m_dot_cooling_air * temperature_difference) +
 723     (h3 * RPM * temperature_difference);
 724   double dqdt_free = h1 * temperature_difference;
 725   double dqdt_cylinder_head = dqdt_from_combustion + dqdt_forced + dqdt_free;
 726
 727   double HeatCapacityCylinderHead = CpCylinderHead * MassCylinderHead;
 728
 729   CylinderHeadTemp_degK +=
 730     (dqdt_cylinder_head / HeatCapacityCylinderHead) * dt;
 731 }
 732
 733 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 734 /**
 735  * Calculate the oil temperature.
 736  *
 737  * Inputs: Percentage_Power, running flag.
 738  *
 739  * Outputs: OilTemp_degK
 740  */
 741
 742 void FGPiston::doOilTemperature(void)
 743 {
 744   double idle_percentage_power = 0.023;        // approximately
 745   double target_oil_temp;        // Steady state oil temp at the current engine conditions
 746   double time_constant;          // The time constant for the differential equation
 747
 748   if (Running) {
 749     target_oil_temp = 363;
 750     time_constant = 500;        // Time constant for engine-on idling.
 751     if (Percentage_Power > idle_percentage_power) {
 752       time_constant /= ((Percentage_Power / idle_percentage_power) / 10.0); // adjust for power
 753     }
 754   } else {
 755     target_oil_temp = RankineToKelvin(Atmosphere->GetTemperature());
 756     time_constant = 1000;  // Time constant for engine-off; reflects the fact
 757                            // that oil is no longer getting circulated
 758   }
 759
 760   double dOilTempdt = (target_oil_temp - OilTemp_degK) / time_constant;
 761
 762   OilTemp_degK += (dOilTempdt * dt);
 763 }
 764
 765 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 766 /**
 767  * Calculate the oil pressure.
 768  *
 769  * Inputs: RPM
 770  *
 771  * Outputs: OilPressure_psi
 772  */
 773
 774 void FGPiston::doOilPressure(void)
 775 {
 776   double Oil_Press_Relief_Valve = 60; // FIXME: may vary by engine
 777   double Oil_Press_RPM_Max = MaxRPM * 0.75;    // 75% of max rpm FIXME: may vary by engine
 778   double Design_Oil_Temp = 358;          // degK; FIXME: may vary by engine
 779   double Oil_Viscosity_Index = 0.25;
 780
 781   OilPressure_psi = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * RPM;
 782
 783   if (OilPressure_psi >= Oil_Press_Relief_Valve) {
 784     OilPressure_psi = Oil_Press_Relief_Valve;
 785   }
 786
 787   OilPressure_psi += (Design_Oil_Temp - OilTemp_degK) * Oil_Viscosity_Index;
 788 }
 789
 790 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 791
 792 string FGPiston::GetEngineLabels(string delimeter)
 793 {
 794   std::ostringstream buf;
 795
 796   buf << Name << " Power Available (engine " << EngineNumber << " in HP)" << delimeter
 797       << Name << " HP (engine " << EngineNumber << ")" << delimeter
 798       << Name << " equivalent ratio (engine " << EngineNumber << ")" << delimeter
 799       << Name << " MAP (engine " << EngineNumber << ")" << delimeter
 800       << Thruster->GetThrusterLabels(EngineNumber, delimeter);
 801
 802   return buf.str();
 803 }
 804
 805 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 806
 807 string FGPiston::GetEngineValues(string delimeter)
 808 {
 809   std::ostringstream buf;
 810
 811   buf << PowerAvailable << delimeter << HP << delimeter
 812       << equivalence_ratio << delimeter << MAP << delimeter
 813       << Thruster->GetThrusterValues(EngineNumber, delimeter);
 814
 815   return buf.str();
 816 }
 817
 818 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 819 //
 820 //    The bitmasked value choices are as follows:
 821 //    unset: In this case (the default) JSBSim would only print
 822 //       out the normally expected messages, essentially echoing
 823 //       the config files as they are read. If the environment
 824 //       variable is not set, debug_lvl is set to 1 internally
 825 //    0: This requests JSBSim not to output any messages
 826 //       whatsoever.
 827 //    1: This value explicity requests the normal JSBSim
 828 //       startup messages
 829 //    2: This value asks for a message to be printed out when
 830 //       a class is instantiated
 831 //    4: When this value is set, a message is displayed when a
 832 //       FGModel object executes its Run() method
 833 //    8: When this value is set, various runtime state variables
 834 //       are printed out periodically
 835 //    16: When set various parameters are sanity checked and
 836 //       a message is printed out when they go out of bounds
 837
 838 void FGPiston::Debug(int from)
 839 {
 840   if (debug_lvl <= 0) return;
 841
 842   if (debug_lvl & 1) { // Standard console startup message output
 843     if (from == 0) { // Constructor
 844
 845       cout << "\n    Engine Name: "         << Name << endl;
 846       cout << "      MinManifoldPressure: " << MinManifoldPressure_inHg << endl;
 847       cout << "      MaxManifoldPressure: " << MaxManifoldPressure_inHg << endl;
 848       cout << "      MinMaP (Pa):         " << minMAP << endl;
 849       cout << "      MaxMaP (Pa): "         << maxMAP << endl;
 850       cout << "      Displacement: "        << Displacement             << endl;
 851       cout << "      MaxHP: "               << MaxHP                    << endl;
 852       cout << "      Cycles: "              << Cycles                   << endl;
 853       cout << "      IdleRPM: "             << IdleRPM                  << endl;
 854       cout << "      MaxThrottle: "         << MaxThrottle              << endl;
 855       cout << "      MinThrottle: "         << MinThrottle              << endl;
 856
 857       cout << endl;
 858       cout << "      Combustion Efficiency table:" << endl;
 859       Lookup_Combustion_Efficiency->Print();
 860       cout << endl;
 861
 862       cout << endl;
 863       cout << "      Power Mixture Correlation table:" << endl;
 864       Power_Mixture_Correlation->Print();
 865       cout << endl;
 866
 867       cout << endl;
 868       cout << "      Mixture Efficiency Correlation table:" << endl;
 869       Mixture_Efficiency_Correlation->Print();
 870       cout << endl;
 871
 872     }
 873   }
 874   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 875     if (from == 0) cout << "Instantiated: FGPiston" << endl;
 876     if (from == 1) cout << "Destroyed:    FGPiston" << endl;
 877   }
 878   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 879   }
 880   if (debug_lvl & 8 ) { // Runtime state variables
 881   }
 882   if (debug_lvl & 16) { // Sanity checking
 883   }
 884   if (debug_lvl & 64) {
 885     if (from == 0) { // Constructor
 886       cout << IdSrc << endl;
 887       cout << IdHdr << endl;
 888     }
 889   }
 890 }
 891 } // namespace JSBSim