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