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 = -1;
  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
 238   // Create a BSFC to match the engine if not provided
 239   // The 0.8 in the equation below is volumetric efficiency
 240   if (BSFC < 0) {
 241       BSFC = ( Displacement * MaxRPM * 0.8 ) / (9411 * MaxHP);
 242   }
 243   char property_name[80];
 244   snprintf(property_name, 80, "propulsion/engine[%d]/power_hp", EngineNumber);
 245   PropertyManager->Tie(property_name, &HP);
 246   snprintf(property_name, 80, "propulsion/engine[%d]/bsfc", EngineNumber);
 247   PropertyManager->Tie(property_name, &BSFC);
 248
 249   minMAP = MinManifoldPressure_inHg * inhgtopa;  // inHg to Pa
 250   maxMAP = MaxManifoldPressure_inHg * inhgtopa;
 251   StarterHP = sqrt(MaxHP) * 0.4;
 252
 253   // Set up and sanity-check the turbo/supercharging configuration based on the input values.
 254   if (TakeoffBoost > RatedBoost[0]) bTakeoffBoost = true;
 255   for (i=0; i<BoostSpeeds; ++i) {
 256     bool bad = false;
 257     if (RatedBoost[i] <= 0.0) bad = true;
 258     if (RatedPower[i] <= 0.0) bad = true;
 259     if (RatedAltitude[i] < 0.0) bad = true;  // 0.0 is deliberately allowed - this corresponds to unregulated supercharging.
 260     if (i > 0 && RatedAltitude[i] < RatedAltitude[i - 1]) bad = true;
 261     if (bad) {
 262       // We can't recover from the above - don't use this supercharger speed.
 263       BoostSpeeds--;
 264       // TODO - put out a massive error message!
 265       break;
 266     }
 267     // Now sanity-check stuff that is recoverable.
 268     if (i < BoostSpeeds - 1) {
 269       if (BoostSwitchAltitude[i] < RatedAltitude[i]) {
 270         // TODO - put out an error message
 271         // But we can also make a reasonable estimate, as below.
 272         BoostSwitchAltitude[i] = RatedAltitude[i] + 1000;
 273       }
 274       BoostSwitchPressure[i] = Atmosphere->GetPressure(BoostSwitchAltitude[i]) * psftopa;
 275       //cout << "BoostSwitchAlt = " << BoostSwitchAltitude[i] << ", pressure = " << BoostSwitchPressure[i] << '\n';
 276       // Assume there is some hysteresis on the supercharger gear switch, and guess the value for now
 277       BoostSwitchHysteresis = 1000;
 278     }
 279     // Now work out the supercharger pressure multiplier of this speed from the rated boost and altitude.
 280     RatedMAP[i] = Atmosphere->GetPressureSL() * psftopa + RatedBoost[i] * 6895;  // psi*6895 = Pa.
 281     // Sometimes a separate BCV setting for takeoff or extra power is fitted.
 282     if (TakeoffBoost > RatedBoost[0]) {
 283       // Assume that the effect on the BCV is the same whichever speed is in use.
 284       TakeoffMAP[i] = RatedMAP[i] + ((TakeoffBoost - RatedBoost[0]) * 6895);
 285       bTakeoffBoost = true;
 286     } else {
 287       TakeoffMAP[i] = RatedMAP[i];
 288       bTakeoffBoost = false;
 289     }
 290     BoostMul[i] = RatedMAP[i] / (Atmosphere->GetPressure(RatedAltitude[i]) * psftopa);
 291
 292   }
 293
 294   if (BoostSpeeds > 0) {
 295     Boosted = true;
 296     BoostSpeed = 0;
 297   }
 298   bBoostOverride = (BoostOverride == 1 ? true : false);
 299   if (MinThrottle < 0.001) MinThrottle = 0.001;  //MinThrottle is a denominator in a power equation so it can't be zero
 300   Debug(0); // Call Debug() routine from constructor if needed
 301 }
 302
 303 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 304
 305 FGPiston::~FGPiston()
 306 {
 307   delete Lookup_Combustion_Efficiency;
 308   delete Power_Mixture_Correlation;
 309   delete Mixture_Efficiency_Correlation;
 310   Debug(1); // Call Debug() routine from constructor if needed
 311 }
 312
 313 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 314
 315 void FGPiston::ResetToIC(void)
 316 {
 317   FGEngine::ResetToIC();
 318
 319   ManifoldPressure_inHg = Atmosphere->GetPressure() * psftoinhg; // psf to in Hg
 320   MAP = Atmosphere->GetPressure() * psftopa;
 321   double airTemperature_degK = RankineToKelvin(Atmosphere->GetTemperature());
 322   OilTemp_degK = airTemperature_degK;
 323   CylinderHeadTemp_degK = airTemperature_degK;
 324   ExhaustGasTemp_degK = airTemperature_degK;
 325   EGT_degC = ExhaustGasTemp_degK - 273;
 326   Thruster->SetRPM(0.0);
 327   RPM = 0.0;
 328 }
 329
 330 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 331
 332 double FGPiston::Calculate(void)
 333 {
 334   if (FuelFlow_gph > 0.0) ConsumeFuel();
 335
 336   Throttle = FCS->GetThrottlePos(EngineNumber);
 337   ThrottlePos = MinThrottle+((MaxThrottle-MinThrottle)*Throttle );
 338   Mixture = FCS->GetMixturePos(EngineNumber);
 339
 340   //
 341   // Input values.
 342   //
 343
 344   p_amb = Atmosphere->GetPressure() * psftopa;
 345   p_amb_sea_level = Atmosphere->GetPressureSL() * psftopa;
 346   T_amb = RankineToKelvin(Atmosphere->GetTemperature());
 347
 348   RPM = Thruster->GetRPM() * Thruster->GetGearRatio();
 349
 350   IAS = Auxiliary->GetVcalibratedKTS();
 351
 352   doEngineStartup();
 353   if (Boosted) doBoostControl();
 354   doMAP();
 355   doAirFlow();
 356   doFuelFlow();
 357
 358   //Now that the fuel flow is done check if the mixture is too lean to run the engine
 359   //Assume lean limit at 22 AFR for now - thats a thi of 0.668
 360   //This might be a bit generous, but since there's currently no audiable warning of impending
 361   //cutout in the form of misfiring and/or rough running its probably reasonable for now.
 362 //  if (equivalence_ratio < 0.668)
 363 //    Running = false;
 364
 365   doEnginePower();
 366 if(HP<0.1250)
 367   Running = false;
 368
 369   doEGT();
 370   doCHT();
 371   doOilTemperature();
 372   doOilPressure();
 373
 374   if (Thruster->GetType() == FGThruster::ttPropeller) {
 375     ((FGPropeller*)Thruster)->SetAdvance(FCS->GetPropAdvance(EngineNumber));
 376     ((FGPropeller*)Thruster)->SetFeather(FCS->GetPropFeather(EngineNumber));
 377   }
 378
 379   PowerAvailable = (HP * hptoftlbssec) - Thruster->GetPowerRequired();
 380
 381   return Thruster->Calculate(PowerAvailable);
 382 }
 383
 384 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 385
 386 double FGPiston::CalcFuelNeed(void)
 387 {
 388   return FuelFlow_gph / 3600 * 6 * State->Getdt() * Propulsion->GetRate();
 389 }
 390
 391 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 392
 393 int FGPiston::InitRunning(void) {
 394   Magnetos=3;
 395   //Thruster->SetRPM( 1.1*IdleRPM/Thruster->GetGearRatio() );
 396   Thruster->SetRPM( 1000 );
 397   Running=true;
 398   return 1;
 399 }
 400
 401 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 402 /**
 403  * Start or stop the engine.
 404  */
 405
 406 void FGPiston::doEngineStartup(void)
 407 {
 408   // Check parameters that may alter the operating state of the engine.
 409   // (spark, fuel, starter motor etc)
 410   bool spark;
 411   bool fuel;
 412
 413   // Check for spark
 414   Magneto_Left = false;
 415   Magneto_Right = false;
 416   // Magneto positions:
 417   // 0 -> off
 418   // 1 -> left only
 419   // 2 -> right only
 420   // 3 -> both
 421   if (Magnetos != 0) {
 422     spark = true;
 423   } else {
 424     spark = false;
 425   }  // neglects battery voltage, master on switch, etc for now.
 426
 427   if ((Magnetos == 1) || (Magnetos > 2)) Magneto_Left = true;
 428   if (Magnetos > 1)  Magneto_Right = true;
 429
 430   // Assume we have fuel for now
 431   fuel = !Starved;
 432
 433   // Check if we are turning the starter motor
 434   if (Cranking != Starter) {
 435     // This check saves .../cranking from getting updated every loop - they
 436     // only update when changed.
 437     Cranking = Starter;
 438     crank_counter = 0;
 439   }
 440
 441   if (Cranking) crank_counter++;  //Check mode of engine operation
 442
 443   if (!Running && spark && fuel) {  // start the engine if revs high enough
 444     if (Cranking) {
 445       if ((RPM > IdleRPM*0.8) && (crank_counter > 175)) // Add a little delay to startup
 446         Running = true;                         // on the starter
 447     } else {
 448       if (RPM > IdleRPM*0.8)                            // This allows us to in-air start
 449         Running = true;                         // when windmilling
 450     }
 451   }
 452
 453   // Cut the engine *power* - Note: the engine may continue to
 454   // spin if the prop is in a moving airstream
 455
 456   if ( Running && (!spark || !fuel) ) Running = false;
 457
 458   // Check for stalling (RPM = 0).
 459   if (Running) {
 460     if (RPM == 0) {
 461       Running = false;
 462     } else if ((RPM <= IdleRPM *0.8 ) && (Cranking)) {
 463       Running = false;
 464     }
 465   }
 466 }
 467
 468 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 469
 470 /**
 471  * Calculate the Current Boost Speed
 472  *
 473  * This function calculates the current turbo/supercharger boost speed
 474  * based on altitude and the (automatic) boost-speed control valve configuration.
 475  *
 476  * Inputs: p_amb, BoostSwitchPressure, BoostSwitchHysteresis
 477  *
 478  * Outputs: BoostSpeed
 479  */
 480
 481 void FGPiston::doBoostControl(void)
 482 {
 483   if(BoostSpeed < BoostSpeeds - 1) {
 484     // Check if we need to change to a higher boost speed
 485     if(p_amb < BoostSwitchPressure[BoostSpeed] - BoostSwitchHysteresis) {
 486       BoostSpeed++;
 487     }
 488   } else if(BoostSpeed > 0) {
 489     // Check if we need to change to a lower boost speed
 490     if(p_amb > BoostSwitchPressure[BoostSpeed - 1] + BoostSwitchHysteresis) {
 491       BoostSpeed--;
 492     }
 493   }
 494 }
 495
 496 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 497
 498 /**
 499  * Calculate the manifold absolute pressure (MAP) in inches hg
 500  *
 501  * This function calculates manifold absolute pressure (MAP)
 502  * from the throttle position, turbo/supercharger boost control
 503  * system, engine speed and local ambient air density.
 504  *
 505  * TODO: changes in MP should not be instantaneous -- introduce
 506  * a lag between throttle changes and MP changes, to allow pressure
 507  * to build up or disperse.
 508  *
 509  * Inputs: minMAP, maxMAP, p_amb, Throttle
 510  *
 511  * Outputs: MAP, ManifoldPressure_inHg
 512  */
 513
 514 void FGPiston::doMAP(void)
 515 {
 516     suction_loss = RPM > 0.0 ? ThrottlePos * MaxRPM / RPM : 1.0;
 517     if (suction_loss > 1.0) suction_loss = 1.0;
 518     MAP = p_amb * suction_loss;
 519
 520     if(Boosted) {
 521       // If takeoff boost is fitted, we currently assume the following throttle map:
 522       // (In throttle % - actual input is 0 -> 1)
 523       // 99 / 100 - Takeoff boost
 524       // 96 / 97 / 98 - Rated boost
 525       // 0 - 95 - Idle to Rated boost (MinManifoldPressure to MaxManifoldPressure)
 526       // In real life, most planes would be fitted with a mechanical 'gate' between
 527       // the rated boost and takeoff boost positions.
 528       double T = Throttle; // processed throttle value.
 529       bool bTakeoffPos = false;
 530       if(bTakeoffBoost) {
 531         if(Throttle > 0.98) {
 532           //cout << "Takeoff Boost!!!!\n";
 533           bTakeoffPos = true;
 534         } else if(Throttle <= 0.95) {
 535           bTakeoffPos = false;
 536           T *= 1.0 / 0.95;
 537         } else {
 538           bTakeoffPos = false;
 539           //cout << "Rated Boost!!\n";
 540           T = 1.0;
 541         }
 542       }
 543       // Boost the manifold pressure.
 544       MAP += MAP * BoostMul[BoostSpeed] * suction_loss * RPM/RatedRPM[BoostSpeed];
 545       // Now clip the manifold pressure to BCV or Wastegate setting.
 546       if(bTakeoffPos) {
 547         if(MAP > TakeoffMAP[BoostSpeed]) {
 548           MAP = TakeoffMAP[BoostSpeed];
 549         }
 550       } else {
 551         if(MAP > RatedMAP[BoostSpeed]) {
 552           MAP = RatedMAP[BoostSpeed];
 553         }
 554       }
 555     }
 556
 557   // And set the value in American units as well
 558   ManifoldPressure_inHg = MAP / inhgtopa;
 559 }
 560
 561 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 562 /**
 563  * Calculate the air flow through the engine.
 564  * Also calculates ambient air density
 565  * (used in CHT calculation for air-cooled engines).
 566  *
 567  * Inputs: p_amb, R_air, T_amb, MAP, Displacement,
 568  *   RPM, volumetric_efficiency, ThrottlePos
 569  *
 570  * TODO: Model inlet manifold air temperature.
 571  *
 572  * Outputs: rho_air, m_dot_air
 573  */
 574
 575 void FGPiston::doAirFlow(void)
 576 {
 577
 578 rho_air = p_amb / (R_air * T_amb);
 579   double displacement_SI = Displacement * in3tom3;
 580   double swept_volume = (displacement_SI * (RPM/60)) / 2;
 581   double v_dot_air = swept_volume * volumetric_efficiency * suction_loss;
 582
 583   double rho_air_manifold = MAP / (R_air * T_amb);
 584   m_dot_air = v_dot_air * rho_air_manifold;
 585 }
 586
 587 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 588 /**
 589  * Calculate the fuel flow into the engine.
 590  *
 591  * Inputs: Mixture, thi_sea_level, p_amb_sea_level, p_amb, m_dot_air
 592  *
 593  * Outputs: equivalence_ratio, m_dot_fuel
 594  */
 595
 596 void FGPiston::doFuelFlow(void)
 597 {
 598   double thi_sea_level = 1.3 * Mixture; // Allows an AFR of infinity:1 to 11.3075:1
 599   equivalence_ratio = thi_sea_level; // * p_amb_sea_level / p_amb;
 600   double AFR = 10+(12*(1-Mixture));// mixture 10:1 to 22:1
 601   m_dot_fuel = m_dot_air / AFR;
 602   FuelFlow_gph = m_dot_fuel
 603     * 3600            // seconds to hours
 604     * 2.2046            // kg to lb
 605     / 6.0;            // lb to gal_us of gasoline
 606 //    / 6.6;            // lb to gal_us of kerosene
 607 }
 608
 609 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 610 /**
 611  * Calculate the power produced by the engine.
 612  *
 613  * Currently, the JSBSim propellor model does not allow the
 614  * engine to produce enough RPMs to get up to a high horsepower.
 615  * When tested with sufficient RPM, it has no trouble reaching
 616  * 200HP.
 617  *
 618  * Inputs: ManifoldPressure_inHg, p_amb, RPM, T_amb,
 619  *   Mixture_Efficiency_Correlation, Cycles, MaxHP
 620  *
 621  * Outputs: Percentage_Power, HP
 622  */
 623
 624 void FGPiston::doEnginePower(void)
 625 {
 626   if (Running) {
 627     double T_amb_degF = KelvinToFahrenheit(T_amb);
 628     double T_amb_sea_lev_degF = KelvinToFahrenheit(288);
 629
 630     // FIXME: this needs to be generalized
 631     double ME, friction, percent_RPM, power;  // Convienience term for use in the calculations
 632     ME = Mixture_Efficiency_Correlation->GetValue(m_dot_fuel/m_dot_air);
 633
 634     percent_RPM = RPM/MaxRPM;
 635     friction = 1 - (percent_RPM * percent_RPM * percent_RPM * percent_RPM/10);
 636     if (friction < 0 ) friction = 0;
 637     power = friction;
 638
 639     if ( Magnetos != 3 ) power *= SparkFailDrop;
 640
 641
 642     HP = (FuelFlow_gph * 6.0 / BSFC )* ME * suction_loss * power;
 643
 644   } else {
 645
 646     // Power output when the engine is not running
 647     if (Cranking) {
 648       if (RPM < 10) {
 649         HP = StarterHP;
 650       } else if (RPM < IdleRPM*0.8) {
 651         HP = StarterHP + ((IdleRPM*0.8 - RPM) / 8.0);
 652         // This is a guess - would be nice to find a proper starter moter torque curve
 653       } else {
 654         HP = StarterHP;
 655       }
 656     } else {
 657       // Quick hack until we port the FMEP stuff
 658       if (RPM > 0.0)
 659         HP = -1.5;
 660       else
 661         HP = 0.0;
 662     }
 663   }
 664   Percentage_Power = HP / MaxHP ;
 665 //  cout << "Power = " << HP << "  RPM = " << RPM << "  Running = " << Running << "  Cranking = " << Cranking << endl;
 666 }
 667
 668 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 669 /**
 670  * Calculate the exhaust gas temperature.
 671  *
 672  * Inputs: equivalence_ratio, m_dot_fuel, calorific_value_fuel,
 673  *   Cp_air, m_dot_air, Cp_fuel, m_dot_fuel, T_amb, Percentage_Power
 674  *
 675  * Outputs: combustion_efficiency, ExhaustGasTemp_degK
 676  */
 677
 678 void FGPiston::doEGT(void)
 679 {
 680   double delta_T_exhaust;
 681   double enthalpy_exhaust;
 682   double heat_capacity_exhaust;
 683   double dEGTdt;
 684
 685   if ((Running) && (m_dot_air > 0.0)) {  // do the energy balance
 686     combustion_efficiency = Lookup_Combustion_Efficiency->GetValue(equivalence_ratio);
 687     enthalpy_exhaust = m_dot_fuel * calorific_value_fuel *
 688                               combustion_efficiency * 0.33;
 689     heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
 690     delta_T_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
 691     ExhaustGasTemp_degK = T_amb + delta_T_exhaust;
 692     ExhaustGasTemp_degK *= 0.444 + ((0.544 - 0.444) * Percentage_Power);
 693   } else {  // Drop towards ambient - guess an appropriate time constant for now
 694     combustion_efficiency = 0;
 695     dEGTdt = (RankineToKelvin(Atmosphere->GetTemperature()) - ExhaustGasTemp_degK) / 100.0;
 696     delta_T_exhaust = dEGTdt * dt;
 697     ExhaustGasTemp_degK += delta_T_exhaust;
 698   }
 699 }
 700
 701 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 702 /**
 703  * Calculate the cylinder head temperature.
 704  *
 705  * Inputs: T_amb, IAS, rho_air, m_dot_fuel, calorific_value_fuel,
 706  *   combustion_efficiency, RPM
 707  *
 708  * Outputs: CylinderHeadTemp_degK
 709  */
 710
 711 void FGPiston::doCHT(void)
 712 {
 713   double h1 = -95.0;
 714   double h2 = -3.95;
 715   double h3 = -0.05;
 716
 717   double arbitary_area = 1.0;
 718   double CpCylinderHead = 800.0;
 719   double MassCylinderHead = 8.0;
 720
 721   double temperature_difference = CylinderHeadTemp_degK - T_amb;
 722   double v_apparent = IAS * 0.5144444;
 723   double v_dot_cooling_air = arbitary_area * v_apparent;
 724   double m_dot_cooling_air = v_dot_cooling_air * rho_air;
 725   double dqdt_from_combustion =
 726     m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
 727   double dqdt_forced = (h2 * m_dot_cooling_air * temperature_difference) +
 728     (h3 * RPM * temperature_difference);
 729   double dqdt_free = h1 * temperature_difference;
 730   double dqdt_cylinder_head = dqdt_from_combustion + dqdt_forced + dqdt_free;
 731
 732   double HeatCapacityCylinderHead = CpCylinderHead * MassCylinderHead;
 733
 734   CylinderHeadTemp_degK +=
 735     (dqdt_cylinder_head / HeatCapacityCylinderHead) * dt;
 736 }
 737
 738 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 739 /**
 740  * Calculate the oil temperature.
 741  *
 742  * Inputs: Percentage_Power, running flag.
 743  *
 744  * Outputs: OilTemp_degK
 745  */
 746
 747 void FGPiston::doOilTemperature(void)
 748 {
 749   double idle_percentage_power = 0.023;        // approximately
 750   double target_oil_temp;        // Steady state oil temp at the current engine conditions
 751   double time_constant;          // The time constant for the differential equation
 752
 753   if (Running) {
 754     target_oil_temp = 363;
 755     time_constant = 500;        // Time constant for engine-on idling.
 756     if (Percentage_Power > idle_percentage_power) {
 757       time_constant /= ((Percentage_Power / idle_percentage_power) / 10.0); // adjust for power
 758     }
 759   } else {
 760     target_oil_temp = RankineToKelvin(Atmosphere->GetTemperature());
 761     time_constant = 1000;  // Time constant for engine-off; reflects the fact
 762                            // that oil is no longer getting circulated
 763   }
 764
 765   double dOilTempdt = (target_oil_temp - OilTemp_degK) / time_constant;
 766
 767   OilTemp_degK += (dOilTempdt * dt);
 768 }
 769
 770 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 771 /**
 772  * Calculate the oil pressure.
 773  *
 774  * Inputs: RPM
 775  *
 776  * Outputs: OilPressure_psi
 777  */
 778
 779 void FGPiston::doOilPressure(void)
 780 {
 781   double Oil_Press_Relief_Valve = 60; // FIXME: may vary by engine
 782   double Oil_Press_RPM_Max = MaxRPM * 0.75;    // 75% of max rpm FIXME: may vary by engine
 783   double Design_Oil_Temp = 358;          // degK; FIXME: may vary by engine
 784   double Oil_Viscosity_Index = 0.25;
 785
 786   OilPressure_psi = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * RPM;
 787
 788   if (OilPressure_psi >= Oil_Press_Relief_Valve) {
 789     OilPressure_psi = Oil_Press_Relief_Valve;
 790   }
 791
 792   OilPressure_psi += (Design_Oil_Temp - OilTemp_degK) * Oil_Viscosity_Index * OilPressure_psi / Oil_Press_Relief_Valve;
 793 }
 794
 795 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 796
 797 string FGPiston::GetEngineLabels(string delimeter)
 798 {
 799   std::ostringstream buf;
 800
 801   buf << Name << " Power Available (engine " << EngineNumber << " in HP)" << delimeter
 802       << Name << " HP (engine " << EngineNumber << ")" << delimeter
 803       << Name << " equivalent ratio (engine " << EngineNumber << ")" << delimeter
 804       << Name << " MAP (engine " << EngineNumber << ")" << delimeter
 805       << Thruster->GetThrusterLabels(EngineNumber, delimeter);
 806
 807   return buf.str();
 808 }
 809
 810 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 811
 812 string FGPiston::GetEngineValues(string delimeter)
 813 {
 814   std::ostringstream buf;
 815
 816   buf << PowerAvailable << delimeter << HP << delimeter
 817       << equivalence_ratio << delimeter << MAP << delimeter
 818       << Thruster->GetThrusterValues(EngineNumber, delimeter);
 819
 820   return buf.str();
 821 }
 822
 823 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 824 //
 825 //    The bitmasked value choices are as follows:
 826 //    unset: In this case (the default) JSBSim would only print
 827 //       out the normally expected messages, essentially echoing
 828 //       the config files as they are read. If the environment
 829 //       variable is not set, debug_lvl is set to 1 internally
 830 //    0: This requests JSBSim not to output any messages
 831 //       whatsoever.
 832 //    1: This value explicity requests the normal JSBSim
 833 //       startup messages
 834 //    2: This value asks for a message to be printed out when
 835 //       a class is instantiated
 836 //    4: When this value is set, a message is displayed when a
 837 //       FGModel object executes its Run() method
 838 //    8: When this value is set, various runtime state variables
 839 //       are printed out periodically
 840 //    16: When set various parameters are sanity checked and
 841 //       a message is printed out when they go out of bounds
 842
 843 void FGPiston::Debug(int from)
 844 {
 845   if (debug_lvl <= 0) return;
 846
 847   if (debug_lvl & 1) { // Standard console startup message output
 848     if (from == 0) { // Constructor
 849
 850       cout << "\n    Engine Name: "         << Name << endl;
 851       cout << "      MinManifoldPressure: " << MinManifoldPressure_inHg << endl;
 852       cout << "      MaxManifoldPressure: " << MaxManifoldPressure_inHg << endl;
 853       cout << "      MinMaP (Pa):         " << minMAP << endl;
 854       cout << "      MaxMaP (Pa): "         << maxMAP << endl;
 855       cout << "      Displacement: "        << Displacement             << endl;
 856       cout << "      MaxHP: "               << MaxHP                    << endl;
 857       cout << "      Cycles: "              << Cycles                   << endl;
 858       cout << "      IdleRPM: "             << IdleRPM                  << endl;
 859       cout << "      MaxThrottle: "         << MaxThrottle              << endl;
 860       cout << "      MinThrottle: "         << MinThrottle              << endl;
 861
 862       cout << endl;
 863       cout << "      Combustion Efficiency table:" << endl;
 864       Lookup_Combustion_Efficiency->Print();
 865       cout << endl;
 866
 867       cout << endl;
 868       cout << "      Power Mixture Correlation table:" << endl;
 869       Power_Mixture_Correlation->Print();
 870       cout << endl;
 871
 872       cout << endl;
 873       cout << "      Mixture Efficiency Correlation table:" << endl;
 874       Mixture_Efficiency_Correlation->Print();
 875       cout << endl;
 876
 877     }
 878   }
 879   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 880     if (from == 0) cout << "Instantiated: FGPiston" << endl;
 881     if (from == 1) cout << "Destroyed:    FGPiston" << endl;
 882   }
 883   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 884   }
 885   if (debug_lvl & 8 ) { // Runtime state variables
 886   }
 887   if (debug_lvl & 16) { // Sanity checking
 888   }
 889   if (debug_lvl & 64) {
 890     if (from == 0) { // Constructor
 891       cout << IdSrc << endl;
 892       cout << IdHdr << endl;
 893     }
 894   }
 895 }
 896 } // namespace JSBSim