src/FDM/JSBSim/models/FGGasCell.cpp

   1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   2
   3  Header:       FGGasCell.h
   4  Author:       Anders Gidenstam
   5  Date started: 01/21/2006
   6
   7  ----- Copyright (C) 2006 - 2008  Anders Gidenstam (anders(at)gidenstam.org) --
   8
   9  This program is free software; you can redistribute it and/or modify it under
  10  the terms of the GNU Lesser General Public License as published by the Free Software
  11  Foundation; either version 2 of the License, or (at your option) any later
  12  version.
  13
  14  This program is distributed in the hope that it will be useful, but WITHOUT
  15  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  16  FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  17  details.
  18
  19  You should have received a copy of the GNU Lesser General Public License along with
  20  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  21  Place - Suite 330, Boston, MA  02111-1307, USA.
  22
  23  Further information about the GNU Lesser General Public License can also be found on
  24  the world wide web at http://www.gnu.org.
  25
  26 FUNCTIONAL DESCRIPTION
  27 --------------------------------------------------------------------------------
  28 See header file.
  29
  30 HISTORY
  31 --------------------------------------------------------------------------------
  32 01/21/2006  AG   Created
  33
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 INCLUDES
  36 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  37
  38 #include <FGFDMExec.h>
  39 #include <models/FGAuxiliary.h>
  40 #include <models/FGAtmosphere.h>
  41 #include <models/FGInertial.h>
  42 #include <models/FGMassBalance.h>
  43 #include "FGGasCell.h"
  44
  45 #if !defined ( sgi ) || defined( __GNUC__ ) && (_COMPILER_VERSION < 740)
  46 using std::cerr;
  47 using std::endl;
  48 using std::cout;
  49 #endif
  50
  51 namespace JSBSim {
  52
  53 static const char *IdSrc = "$Id$";
  54 static const char *IdHdr = ID_GASCELL;
  55
  56 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  57 CLASS IMPLEMENTATION
  58 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  59 /* Constants. */
  60 const double FGGasCell::R = 3.4071;              // [lbs ft/(mol Rankine)]
  61 const double FGGasCell::M_air = 0.0019186;       // [slug/mol]
  62 const double FGGasCell::M_hydrogen = 0.00013841; // [slug/mol]
  63 const double FGGasCell::M_helium = 0.00027409;   // [slug/mol]
  64
  65 FGGasCell::FGGasCell(FGFDMExec* exec, Element* el, int num) : FGForce(exec)
  66 {
  67   string token;
  68   Element* element;
  69
  70   Auxiliary = exec->GetAuxiliary();
  71   Atmosphere = exec->GetAtmosphere();
  72   PropertyManager = exec->GetPropertyManager();
  73   Inertial = exec->GetInertial();
  74   MassBalance = exec->GetMassBalance();
  75
  76   gasCellJ = FGMatrix33();
  77   gasCellM = FGColumnVector3();
  78
  79   Buoyancy = MaxVolume = MaxOverpressure = Temperature = Pressure =
  80     Contents = Volume = dVolumeIdeal = 0.0;
  81   Xradius = Yradius = Zradius = Xwidth = Ywidth = Zwidth = 0.0;
  82   ValveCoefficient = ValveOpen = 0.0;
  83   CellNum = num;
  84
  85   // NOTE: In the local system X points north, Y points east and Z points down.
  86   SetTransformType(FGForce::tLocalBody);
  87
  88   type = el->GetAttributeValue("type");
  89   if      (type == "HYDROGEN") Type = ttHYDROGEN;
  90   else if (type == "HELIUM")   Type = ttHELIUM;
  91   else if (type == "AIR")      Type = ttAIR;
  92   else                         Type = ttUNKNOWN;
  93
  94   element = el->FindElement("location");
  95   if (element) {
  96     vXYZ = element->FindElementTripletConvertTo("IN");
  97   } else {
  98     cerr << "Fatal Error: No location found for this gas cell." << endl;
  99     exit(-1);
 100   }
 101   if ((el->FindElement("x_radius") || el->FindElement("x_width")) &&
 102       (el->FindElement("y_radius") || el->FindElement("y_width")) &&
 103       (el->FindElement("z_radius") || el->FindElement("z_width"))) {
 104
 105     if (el->FindElement("x_radius")) {
 106       Xradius = el->FindElementValueAsNumberConvertTo("x_radius", "FT");
 107     }
 108     if (el->FindElement("y_radius")) {
 109       Yradius = el->FindElementValueAsNumberConvertTo("y_radius", "FT");
 110     }
 111     if (el->FindElement("z_radius")) {
 112       Zradius = el->FindElementValueAsNumberConvertTo("z_radius", "FT");
 113     }
 114
 115     if (el->FindElement("x_width")) {
 116       Xwidth = el->FindElementValueAsNumberConvertTo("x_width", "FT");
 117     }
 118     if (el->FindElement("y_width")) {
 119       Ywidth = el->FindElementValueAsNumberConvertTo("y_width", "FT");
 120     }
 121     if (el->FindElement("z_width")) {
 122       Zwidth = el->FindElementValueAsNumberConvertTo("z_width", "FT");
 123     }
 124
 125     // The volume is a (potentially) extruded ellipsoid.
 126     // However, currently only a few combinations of radius and width are
 127     // fully supported.
 128     if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 129         (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 130       // Ellipsoid volume.
 131       MaxVolume = 4.0  * M_PI * Xradius * Yradius * Zradius / 3.0;
 132     } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 133                 (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 134       // Cylindrical volume.
 135       MaxVolume = M_PI * Yradius * Zradius * Xwidth;
 136     } else {
 137       cerr << "Warning: Unsupported gas cell shape." << endl;
 138       MaxVolume =
 139         (4.0  * M_PI * Xradius * Yradius * Zradius / 3.0 +
 140          M_PI * Yradius * Zradius * Xwidth +
 141          M_PI * Xradius * Zradius * Ywidth +
 142          M_PI * Xradius * Yradius * Zwidth +
 143          2.0  * Xradius * Ywidth * Zwidth +
 144          2.0  * Yradius * Xwidth * Zwidth +
 145          2.0  * Zradius * Xwidth * Ywidth +
 146          Xwidth * Ywidth * Zwidth);
 147     }
 148   } else {
 149     cerr << "Fatal Error: Gas cell shape must be given." << endl;
 150     exit(-1);
 151   }
 152   if (el->FindElement("max_overpressure")) {
 153     MaxOverpressure = el->FindElementValueAsNumberConvertTo("max_overpressure",
 154                                                             "LBS/FT2");
 155   }
 156   if (el->FindElement("fullness")) {
 157     const double Fullness = el->FindElementValueAsNumber("fullness");
 158     if (0 <= Fullness) {
 159       Volume = Fullness * MaxVolume;
 160     } else {
 161       cerr << "Warning: Invalid initial gas cell fullness value." << endl;
 162     }
 163   }
 164   if (el->FindElement("valve_coefficient")) {
 165     ValveCoefficient =
 166       el->FindElementValueAsNumberConvertTo("valve_coefficient",
 167                                             "FT4*SEC/SLUG");
 168     ValveCoefficient = max(ValveCoefficient, 0.0);
 169   }
 170
 171   // Initialize state
 172   SetLocation(vXYZ);
 173
 174   if (Temperature == 0.0) {
 175     Temperature = Atmosphere->GetTemperature();
 176   }
 177   if (Pressure == 0.0) {
 178     Pressure = Atmosphere->GetPressure();
 179   }
 180   if (Volume != 0.0) {
 181     // Calculate initial gas content.
 182     Contents = Pressure * Volume / (R * Temperature);
 183
 184     // Clip to max allowed value.
 185     const double IdealPressure = Contents * R * Temperature / MaxVolume;
 186     if (IdealPressure > Pressure + MaxOverpressure) {
 187       Contents = (Pressure + MaxOverpressure) * MaxVolume / (R * Temperature);
 188       Pressure = Pressure + MaxOverpressure;
 189     } else {
 190       Pressure = max(IdealPressure, Pressure);
 191     }
 192   } else {
 193     // Calculate initial gas content.
 194     Contents = Pressure * MaxVolume / (R * Temperature);
 195   }
 196
 197   Volume = Contents * R * Temperature / Pressure;
 198   Mass = Contents * M_gas();
 199
 200   // Bind relevant properties
 201   char property_name[80];
 202   snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/max_volume-ft3",
 203            CellNum);
 204   PropertyManager->Tie( property_name, &MaxVolume );
 205   PropertyManager->SetWritable( property_name, false );
 206   snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/temp-R",
 207            CellNum);
 208   PropertyManager->Tie( property_name, &Temperature );
 209   snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/pressure-psf",
 210            CellNum);
 211   PropertyManager->Tie( property_name, &Pressure );
 212   snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/volume-ft3",
 213            CellNum);
 214   PropertyManager->Tie( property_name, &Volume );
 215   snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/buoyancy-lbs",
 216            CellNum);
 217   PropertyManager->Tie( property_name, &Buoyancy );
 218   snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/contents-mol",
 219            CellNum);
 220   PropertyManager->Tie( property_name, &Contents );
 221   snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/valve_open",
 222            CellNum);
 223   PropertyManager->Tie( property_name, &ValveOpen );
 224
 225   Debug(0);
 226
 227   // Read heat transfer coefficients
 228   if (Element* heat = el->FindElement("heat")) {
 229     Element* function_element = heat->FindElement("function");
 230     while (function_element) {
 231       HeatTransferCoeff.push_back(new FGFunction(PropertyManager,
 232                                                  function_element));
 233       function_element = heat->FindNextElement("function");
 234     }
 235   }
 236
 237   // Load ballonets if there are any
 238   if (Element* ballonet_element = el->FindElement("ballonet")) {
 239     while (ballonet_element) {
 240       Ballonet.push_back(new FGBallonet(exec,
 241                                         ballonet_element,
 242                                         Ballonet.size(),
 243                                         this));
 244       ballonet_element = el->FindNextElement("ballonet");
 245     }
 246   }
 247
 248 }
 249
 250 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 251
 252 FGGasCell::~FGGasCell()
 253 {
 254   unsigned int i;
 255
 256   for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
 257   HeatTransferCoeff.clear();
 258
 259   for (i = 0; i < Ballonet.size(); i++) delete Ballonet[i];
 260   Ballonet.clear();
 261
 262   Debug(1);
 263 }
 264
 265 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 266
 267 void FGGasCell::Calculate(double dt)
 268 {
 269   const double AirTemperature = Atmosphere->GetTemperature();  // [Rankine]
 270   const double AirPressure    = Atmosphere->GetPressure();     // [lbs/ft²]
 271   const double AirDensity     = Atmosphere->GetDensity();      // [slug/ft³]
 272   const double g = Inertial->gravity();                        // [lbs/slug]
 273
 274   const double OldTemperature = Temperature;
 275   const double OldPressure    = Pressure;
 276   unsigned int i;
 277   const unsigned int no_ballonets = Ballonet.size();
 278
 279   //-- Read ballonet state --
 280   // NOTE: This model might need a more proper integration technique.
 281   double BallonetsVolume = 0.0;
 282   double BallonetsHeatFlow = 0.0;
 283   for (i = 0; i < no_ballonets; i++) {
 284     BallonetsVolume   += Ballonet[i]->GetVolume();
 285     BallonetsHeatFlow += Ballonet[i]->GetHeatFlow();
 286   }
 287
 288   //-- Gas temperature --
 289
 290   if (HeatTransferCoeff.size() > 0) {
 291     // The model is based on the ideal gas law.
 292     // However, it does look a bit fishy. Please verify.
 293     //   dT/dt = dU / (Cv n R)
 294     double dU = 0.0;
 295     for (i = 0; i < HeatTransferCoeff.size(); i++) {
 296       dU += HeatTransferCoeff[i]->GetValue();
 297     }
 298     // Don't include dt when accounting for adiabatic expansion/contraction.
 299     // The rate of adiabatic cooling looks about right: ~5.4 Rankine/1000ft.
 300     if (Contents > 0) {
 301       Temperature +=
 302         (dU * dt - Pressure * dVolumeIdeal - BallonetsHeatFlow) /
 303         (Cv_gas() * Contents * R);
 304     } else {
 305       Temperature = AirTemperature;
 306     }
 307   } else {
 308     // No simulation of complex temperature changes.
 309     // Note: Making the gas cell behave adiabatically might be a better
 310     // option.
 311     Temperature = AirTemperature;
 312   }
 313
 314   //-- Pressure --
 315   const double IdealPressure =
 316     Contents * R * Temperature / (MaxVolume - BallonetsVolume);
 317   if (IdealPressure > AirPressure + MaxOverpressure) {
 318     Pressure = AirPressure + MaxOverpressure;
 319   } else {
 320     Pressure = max(IdealPressure, AirPressure);
 321   }
 322
 323   //-- Manual valving --
 324
 325   // FIXME: Presently the effect of manual valving is computed using
 326   //        an ad hoc formula which might not be a good representation
 327   //        of reality.
 328   if ((ValveCoefficient > 0.0) && (ValveOpen > 0.0)) {
 329     // First compute the difference in pressure between the gas in the
 330     // cell and the air above it.
 331     // FixMe: CellHeight should depend on current volume.
 332     const double CellHeight = 2 * Zradius + Zwidth;                   // [ft]
 333     const double GasMass    = Contents * M_gas();                     // [slug]
 334     const double GasVolume  = Contents * R * Temperature / Pressure;  // [ft³]
 335     const double GasDensity = GasMass / GasVolume;
 336     const double DeltaPressure =
 337       Pressure + CellHeight * g * (AirDensity - GasDensity) - AirPressure;
 338     const double VolumeValved =
 339       ValveOpen * ValveCoefficient * DeltaPressure * dt;
 340     Contents =
 341       max(0.0, Contents - Pressure * VolumeValved / (R * Temperature));
 342   }
 343
 344   //-- Update ballonets. --
 345   // Doing that here should give them the opportunity to react to the
 346   // new pressure.
 347   BallonetsVolume = 0.0;
 348   for (i = 0; i < no_ballonets; i++) {
 349     Ballonet[i]->Calculate(dt);
 350     BallonetsVolume += Ballonet[i]->GetVolume();
 351   }
 352
 353   //-- Automatic safety valving. --
 354   if (Contents * R * Temperature / (MaxVolume - BallonetsVolume) >
 355       AirPressure + MaxOverpressure) {
 356     // Gas is automatically valved. Valving capacity is assumed to be infinite.
 357     // FIXME: This could/should be replaced by damage to the gas cell envelope.
 358     Contents =
 359       (AirPressure + MaxOverpressure) *
 360       (MaxVolume - BallonetsVolume) / (R * Temperature);
 361   }
 362
 363   //-- Volume --
 364   Volume = Contents * R * Temperature / Pressure + BallonetsVolume;
 365   dVolumeIdeal =
 366     Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);
 367
 368   //-- Current buoyancy --
 369   // The buoyancy is computed using the atmospheres local density.
 370   Buoyancy = Volume * AirDensity * g;
 371
 372   // Note: This is gross buoyancy. The weight of the gas itself and
 373   // any ballonets is not deducted here as the effects of the gas mass
 374   // is handled by FGMassBalance.
 375   vFn.InitMatrix(0.0, 0.0, - Buoyancy);
 376
 377   // Compute the inertia of the gas cell.
 378   // Consider the gas cell as a shape of uniform density.
 379   // FIXME: If the cell isn't ellipsoid or cylindrical the inertia will
 380   //        be wrong.
 381   gasCellJ = FGMatrix33();
 382   const double mass = Contents * M_gas();
 383   double Ixx, Iyy, Izz;
 384   if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 385       (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 386     // Ellipsoid volume.
 387     Ixx = (1.0 / 5.0) * mass * (Yradius*Yradius + Zradius*Zradius);
 388     Iyy = (1.0 / 5.0) * mass * (Xradius*Xradius + Zradius*Zradius);
 389     Izz = (1.0 / 5.0) * mass * (Xradius*Xradius + Yradius*Yradius);
 390   } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 391               (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 392     // Cylindrical volume (might not be valid with an elliptical cross-section).
 393     Ixx = (1.0 / 2.0) * mass * Yradius * Zradius;
 394     Iyy =
 395       (1.0 / 4.0) * mass * Yradius * Zradius +
 396       (1.0 / 12.0) * mass * Xwidth * Xwidth;
 397     Izz =
 398       (1.0 / 4.0) * mass * Yradius * Zradius +
 399       (1.0 / 12.0) * mass * Xwidth * Xwidth;
 400   } else {
 401     // Not supported. Revert to pointmass model.
 402     Ixx = Iyy = Izz = 0.0;
 403   }
 404   // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
 405   gasCellJ(1,1) = Ixx;
 406   gasCellJ(2,2) = Iyy;
 407   gasCellJ(3,3) = Izz;
 408   Mass = mass;
 409   gasCellM.InitMatrix();
 410   gasCellM(eX) +=
 411     GetXYZ(eX) * Mass*slugtolb;
 412   gasCellM(eY) +=
 413     GetXYZ(eY) * Mass*slugtolb;
 414   gasCellM(eZ) +=
 415     GetXYZ(eZ) * Mass*slugtolb;
 416
 417   if (no_ballonets > 0) {
 418     // Add the mass, moment and inertia of any ballonets.
 419     const FGColumnVector3 p = MassBalance->StructuralToBody( GetXYZ() );
 420
 421     for (i = 0; i < no_ballonets; i++) {
 422       Mass += Ballonet[i]->GetMass();
 423
 424       // Add ballonet moments.
 425       gasCellM(eX) +=
 426         Ballonet[i]->GetXYZ(eX) * Ballonet[i]->GetMass()*slugtolb;
 427       gasCellM(eY) +=
 428         Ballonet[i]->GetXYZ(eY) * Ballonet[i]->GetMass()*slugtolb;
 429       gasCellM(eZ) +=
 430         Ballonet[i]->GetXYZ(eZ) * Ballonet[i]->GetMass()*slugtolb;
 431
 432       // Moments of inertia must be converted to the gas cell frame here.
 433       FGColumnVector3 v =
 434         MassBalance->StructuralToBody( Ballonet[i]->GetXYZ() ) - p;
 435       // Body basis is in FT.
 436       const double mass = Ballonet[i]->GetMass();
 437       gasCellJ += Ballonet[i]->GetInertia() +
 438         FGMatrix33( 0,                - mass*v(1)*v(2), - mass*v(1)*v(3),
 439                     - mass*v(2)*v(1), 0,                - mass*v(2)*v(3),
 440                     - mass*v(3)*v(1), - mass*v(3)*v(2), 0 );
 441     }
 442   }
 443 }
 444
 445 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 446 //    The bitmasked value choices are as follows:
 447 //    unset: In this case (the default) JSBSim would only print
 448 //       out the normally expected messages, essentially echoing
 449 //       the config files as they are read. If the environment
 450 //       variable is not set, debug_lvl is set to 1 internally
 451 //    0: This requests JSBSim not to output any messages
 452 //       whatsoever.
 453 //    1: This value explicity requests the normal JSBSim
 454 //       startup messages
 455 //    2: This value asks for a message to be printed out when
 456 //       a class is instantiated
 457 //    4: When this value is set, a message is displayed when a
 458 //       FGModel object executes its Run() method
 459 //    8: When this value is set, various runtime state variables
 460 //       are printed out periodically
 461 //    16: When set various parameters are sanity checked and
 462 //       a message is printed out when they go out of bounds
 463
 464 void FGGasCell::Debug(int from)
 465 {
 466   if (debug_lvl <= 0) return;
 467
 468   if (debug_lvl & 1) { // Standard console startup message output
 469     if (from == 0) { // Constructor
 470       cout << "    Gas cell holds " << Contents << " mol " <<
 471         type << endl;
 472       cout << "      Cell location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
 473         vXYZ(eY) << ", " << vXYZ(eZ) << endl;
 474       cout << "      Maximum volume: " << MaxVolume << " ft3" << endl;
 475       cout << "      Relief valve release pressure: " << MaxOverpressure <<
 476         " lbs/ft2" << endl;
 477       cout << "      Manual valve coefficient: " << ValveCoefficient <<
 478         " ft4*sec/slug" << endl;
 479       cout << "      Initial temperature: " << Temperature << " Rankine" <<
 480         endl;
 481       cout << "      Initial pressure: " << Pressure << " lbs/ft2" << endl;
 482       cout << "      Initial volume: " << Volume << " ft3" << endl;
 483       cout << "      Initial mass: " << GetMass() << " slug mass" << endl;
 484       cout << "      Initial weight: " << GetMass()*lbtoslug << " lbs force" <<
 485         endl;
 486       cout << "      Heat transfer: " << endl;
 487     }
 488   }
 489   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 490     if (from == 0) cout << "Instantiated: FGGasCell" << endl;
 491     if (from == 1) cout << "Destroyed:    FGGasCell" << endl;
 492   }
 493   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 494   }
 495   if (debug_lvl & 8 ) { // Runtime state variables
 496       cout << "      " << type << " cell holds " << Contents << " mol " << endl;
 497       cout << "      Temperature: " << Temperature << " Rankine" << endl;
 498       cout << "      Pressure: " << Pressure << " lbs/ft2" << endl;
 499       cout << "      Volume: " << Volume << " ft3" << endl;
 500       cout << "      Mass: " << GetMass() << " slug mass" << endl;
 501       cout << "      Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
 502   }
 503   if (debug_lvl & 16) { // Sanity checking
 504   }
 505   if (debug_lvl & 64) {
 506     if (from == 0) { // Constructor
 507       cout << IdSrc << endl;
 508       cout << IdHdr << endl;
 509     }
 510   }
 511 }
 512
 513 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 514 CLASS IMPLEMENTATION
 515 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 516 const double FGBallonet::R = 3.4071;              // [lbs ft/(mol Rankine)]
 517 const double FGBallonet::M_air = 0.0019186;       // [slug/mol]
 518 const double FGBallonet::Cv_air = 5.0/2.0;        // [??]
 519
 520 FGBallonet::FGBallonet(FGFDMExec* exec, Element* el, int num, FGGasCell* parent)
 521 {
 522   string token;
 523   Element* element;
 524
 525   Auxiliary = exec->GetAuxiliary();
 526   Atmosphere = exec->GetAtmosphere();
 527   PropertyManager = exec->GetPropertyManager();
 528   Inertial = exec->GetInertial();
 529
 530   ballonetJ = FGMatrix33();
 531
 532   MaxVolume = MaxOverpressure = Temperature = Pressure =
 533     Contents = Volume = dVolumeIdeal = dU = 0.0;
 534   Xradius = Yradius = Zradius = Xwidth = Ywidth = Zwidth = 0.0;
 535   ValveCoefficient = ValveOpen = 0.0;
 536   BlowerInput = NULL;
 537   CellNum = num;
 538   Parent = parent;
 539
 540   // NOTE: In the local system X points north, Y points east and Z points down.
 541   element = el->FindElement("location");
 542   if (element) {
 543     vXYZ = element->FindElementTripletConvertTo("IN");
 544   } else {
 545     cerr << "Fatal Error: No location found for this ballonet." << endl;
 546     exit(-1);
 547   }
 548   if ((el->FindElement("x_radius") || el->FindElement("x_width")) &&
 549       (el->FindElement("y_radius") || el->FindElement("y_width")) &&
 550       (el->FindElement("z_radius") || el->FindElement("z_width"))) {
 551
 552     if (el->FindElement("x_radius")) {
 553       Xradius = el->FindElementValueAsNumberConvertTo("x_radius", "FT");
 554     }
 555     if (el->FindElement("y_radius")) {
 556       Yradius = el->FindElementValueAsNumberConvertTo("y_radius", "FT");
 557     }
 558     if (el->FindElement("z_radius")) {
 559       Zradius = el->FindElementValueAsNumberConvertTo("z_radius", "FT");
 560     }
 561
 562     if (el->FindElement("x_width")) {
 563       Xwidth = el->FindElementValueAsNumberConvertTo("x_width", "FT");
 564     }
 565     if (el->FindElement("y_width")) {
 566       Ywidth = el->FindElementValueAsNumberConvertTo("y_width", "FT");
 567     }
 568     if (el->FindElement("z_width")) {
 569       Zwidth = el->FindElementValueAsNumberConvertTo("z_width", "FT");
 570     }
 571
 572     // The volume is a (potentially) extruded ellipsoid.
 573     // FIXME: However, currently only a few combinations of radius and
 574     //        width are fully supported.
 575     if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 576         (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 577       // Ellipsoid volume.
 578       MaxVolume = 4.0  * M_PI * Xradius * Yradius * Zradius / 3.0;
 579     } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 580                 (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 581       // Cylindrical volume.
 582       MaxVolume = M_PI * Yradius * Zradius * Xwidth;
 583     } else {
 584       cerr << "Warning: Unsupported ballonet shape." << endl;
 585       MaxVolume =
 586         (4.0  * M_PI * Xradius * Yradius * Zradius / 3.0 +
 587          M_PI * Yradius * Zradius * Xwidth +
 588          M_PI * Xradius * Zradius * Ywidth +
 589          M_PI * Xradius * Yradius * Zwidth +
 590          2.0  * Xradius * Ywidth * Zwidth +
 591          2.0  * Yradius * Xwidth * Zwidth +
 592          2.0  * Zradius * Xwidth * Ywidth +
 593          Xwidth * Ywidth * Zwidth);
 594     }
 595   } else {
 596     cerr << "Fatal Error: Ballonet shape must be given." << endl;
 597     exit(-1);
 598   }
 599   if (el->FindElement("max_overpressure")) {
 600     MaxOverpressure = el->FindElementValueAsNumberConvertTo("max_overpressure",
 601                                                             "LBS/FT2");
 602   }
 603   if (el->FindElement("fullness")) {
 604     const double Fullness = el->FindElementValueAsNumber("fullness");
 605     if (0 <= Fullness) {
 606       Volume = Fullness * MaxVolume;
 607     } else {
 608       cerr << "Warning: Invalid initial ballonet fullness value." << endl;
 609     }
 610   }
 611   if (el->FindElement("valve_coefficient")) {
 612     ValveCoefficient =
 613       el->FindElementValueAsNumberConvertTo("valve_coefficient",
 614                                             "FT4*SEC/SLUG");
 615     ValveCoefficient = max(ValveCoefficient, 0.0);
 616   }
 617
 618   // Initialize state
 619   if (Temperature == 0.0) {
 620     Temperature = Parent->GetTemperature();
 621   }
 622   if (Pressure == 0.0) {
 623     Pressure = Parent->GetPressure();
 624   }
 625   if (Volume != 0.0) {
 626     // Calculate initial air content.
 627     Contents = Pressure * Volume / (R * Temperature);
 628
 629     // Clip to max allowed value.
 630     const double IdealPressure = Contents * R * Temperature / MaxVolume;
 631     if (IdealPressure > Pressure + MaxOverpressure) {
 632       Contents = (Pressure + MaxOverpressure) * MaxVolume / (R * Temperature);
 633       Pressure = Pressure + MaxOverpressure;
 634     } else {
 635       Pressure = max(IdealPressure, Pressure);
 636     }
 637   } else {
 638     // Calculate initial air content.
 639     Contents = Pressure * MaxVolume / (R * Temperature);
 640   }
 641
 642   Volume = Contents * R * Temperature / Pressure;
 643
 644   // Bind relevant properties
 645   char property_name[80];
 646   snprintf(property_name, 80,
 647            "buoyant_forces/gas-cell[%d]/ballonet[%d]/max_volume-ft3",
 648            Parent->GetIndex(),
 649            CellNum);
 650   PropertyManager->Tie( property_name, &MaxVolume );
 651   PropertyManager->SetWritable( property_name, false );
 652   snprintf(property_name, 80,
 653            "buoyant_forces/gas-cell[%d]/ballonet[%d]/temp-R",
 654            Parent->GetIndex(),
 655            CellNum);
 656   PropertyManager->Tie( property_name, &Temperature );
 657   snprintf(property_name, 80,
 658            "buoyant_forces/gas-cell[%d]/ballonet[%d]/pressure-psf",
 659            Parent->GetIndex(),
 660            CellNum);
 661   PropertyManager->Tie( property_name, &Pressure );
 662   snprintf(property_name, 80,
 663            "buoyant_forces/gas-cell[%d]/ballonet[%d]/volume-ft3",
 664            Parent->GetIndex(),
 665            CellNum);
 666   PropertyManager->Tie( property_name, &Volume );
 667   snprintf(property_name, 80,
 668            "buoyant_forces/gas-cell[%d]/ballonet[%d]/contents-mol",
 669            Parent->GetIndex(),
 670            CellNum);
 671   PropertyManager->Tie( property_name, &Contents );
 672   snprintf(property_name, 80,
 673            "buoyant_forces/gas-cell[%d]/ballonet[%d]/valve_open",
 674            Parent->GetIndex(),
 675            CellNum);
 676   PropertyManager->Tie( property_name, &ValveOpen );
 677
 678   Debug(0);
 679
 680   // Read heat transfer coefficients
 681   if (Element* heat = el->FindElement("heat")) {
 682     Element* function_element = heat->FindElement("function");
 683     while (function_element) {
 684       HeatTransferCoeff.push_back(new FGFunction(PropertyManager,
 685                                                  function_element));
 686       function_element = heat->FindNextElement("function");
 687     }
 688   }
 689   // Read blower input function
 690   if (Element* blower = el->FindElement("blower_input")) {
 691     Element* function_element = blower->FindElement("function");
 692     BlowerInput = new FGFunction(PropertyManager,
 693                                  function_element);
 694   }
 695 }
 696
 697 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 698
 699 FGBallonet::~FGBallonet()
 700 {
 701   unsigned int i;
 702
 703   for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
 704   HeatTransferCoeff.clear();
 705
 706   delete BlowerInput;
 707   BlowerInput = NULL;
 708
 709   Debug(1);
 710 }
 711
 712 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 713
 714 void FGBallonet::Calculate(double dt)
 715 {
 716   const double ParentPressure = Parent->GetPressure();         // [lbs/ft²]
 717   const double AirPressure    = Atmosphere->GetPressure();     // [lbs/ft²]
 718
 719   const double OldTemperature = Temperature;
 720   const double OldPressure    = Pressure;
 721   unsigned int i;
 722
 723   //-- Gas temperature --
 724
 725   // The model is based on the ideal gas law.
 726   // However, it does look a bit fishy. Please verify.
 727   //   dT/dt = dU / (Cv n R)
 728   dU = 0.0;
 729   for (i = 0; i < HeatTransferCoeff.size(); i++) {
 730     dU += HeatTransferCoeff[i]->GetValue();
 731   }
 732   // dt is already accounted for in dVolumeIdeal.
 733   if (Contents > 0) {
 734     Temperature +=
 735       (dU * dt - Pressure * dVolumeIdeal) / (Cv_air * Contents * R);
 736   } else {
 737     Temperature = Parent->GetTemperature();
 738   }
 739
 740   //-- Pressure --
 741   const double IdealPressure = Contents * R * Temperature / MaxVolume;
 742   // The pressure is at least that of the parent gas cell.
 743   Pressure = max(IdealPressure, ParentPressure);
 744
 745   //-- Blower input --
 746   if (BlowerInput) {
 747     const double AddedVolume = BlowerInput->GetValue() * dt;
 748     if (AddedVolume > 0.0) {
 749       Contents += Pressure * AddedVolume / (R * Temperature);
 750     }
 751   }
 752
 753   //-- Pressure relief and manual valving --
 754   // FIXME: Presently the effect of valving is computed using
 755   //        an ad hoc formula which might not be a good representation
 756   //        of reality.
 757   if ((ValveCoefficient > 0.0) &&
 758       ((ValveOpen > 0.0) || (Pressure > AirPressure + MaxOverpressure))) {
 759     const double DeltaPressure = Pressure - AirPressure;
 760     const double VolumeValved =
 761       ((Pressure > AirPressure + MaxOverpressure) ? 1.0 : ValveOpen) *
 762       ValveCoefficient * DeltaPressure * dt;
 763     // FIXME: Too small values of Contents sometimes leads to NaN.
 764     //        Currently the minimum is restricted to a safe value.
 765     Contents =
 766       max(1.0, Contents - Pressure * VolumeValved / (R * Temperature));
 767   }
 768
 769   //-- Volume --
 770   Volume = Contents * R * Temperature / Pressure;
 771   dVolumeIdeal =
 772     Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);
 773
 774   // Compute the inertia of the ballonet.
 775   // Consider the ballonet as a shape of uniform density.
 776   // FIXME: If the ballonet isn't ellipsoid or cylindrical the inertia will
 777   //        be wrong.
 778   ballonetJ = FGMatrix33();
 779   const double mass = Contents * M_air;
 780   double Ixx, Iyy, Izz;
 781   if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 782       (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 783     // Ellipsoid volume.
 784     Ixx = (1.0 / 5.0) * mass * (Yradius*Yradius + Zradius*Zradius);
 785     Iyy = (1.0 / 5.0) * mass * (Xradius*Xradius + Zradius*Zradius);
 786     Izz = (1.0 / 5.0) * mass * (Xradius*Xradius + Yradius*Yradius);
 787   } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 788               (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 789     // Cylindrical volume (might not be valid with an elliptical cross-section).
 790     Ixx = (1.0 / 2.0) * mass * Yradius * Zradius;
 791     Iyy =
 792       (1.0 / 4.0) * mass * Yradius * Zradius +
 793       (1.0 / 12.0) * mass * Xwidth * Xwidth;
 794     Izz =
 795       (1.0 / 4.0) * mass * Yradius * Zradius +
 796       (1.0 / 12.0) * mass * Xwidth * Xwidth;
 797   } else {
 798     // Not supported. Revert to pointmass model.
 799     Ixx = Iyy = Izz = 0.0;
 800   }
 801   // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
 802   ballonetJ(1,1) = Ixx;
 803   ballonetJ(2,2) = Iyy;
 804   ballonetJ(3,3) = Izz;
 805 }
 806
 807 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 808 //    The bitmasked value choices are as follows:
 809 //    unset: In this case (the default) JSBSim would only print
 810 //       out the normally expected messages, essentially echoing
 811 //       the config files as they are read. If the environment
 812 //       variable is not set, debug_lvl is set to 1 internally
 813 //    0: This requests JSBSim not to output any messages
 814 //       whatsoever.
 815 //    1: This value explicity requests the normal JSBSim
 816 //       startup messages
 817 //    2: This value asks for a message to be printed out when
 818 //       a class is instantiated
 819 //    4: When this value is set, a message is displayed when a
 820 //       FGModel object executes its Run() method
 821 //    8: When this value is set, various runtime state variables
 822 //       are printed out periodically
 823 //    16: When set various parameters are sanity checked and
 824 //       a message is printed out when they go out of bounds
 825
 826 void FGBallonet::Debug(int from)
 827 {
 828   if (debug_lvl <= 0) return;
 829
 830   if (debug_lvl & 1) { // Standard console startup message output
 831     if (from == 0) { // Constructor
 832       cout << "      Ballonet holds " << Contents << " mol air" << endl;
 833       cout << "        Location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
 834         vXYZ(eY) << ", " << vXYZ(eZ) << endl;
 835       cout << "        Maximum volume: " << MaxVolume << " ft3" << endl;
 836       cout << "        Relief valve release pressure: " << MaxOverpressure <<
 837         " lbs/ft2" << endl;
 838       cout << "        Relief valve coefficient: " << ValveCoefficient <<
 839         " ft4*sec/slug" << endl;
 840       cout << "        Initial temperature: " << Temperature << " Rankine" <<
 841         endl;
 842       cout << "        Initial pressure: " << Pressure << " lbs/ft2" << endl;
 843       cout << "        Initial volume: " << Volume << " ft3" << endl;
 844       cout << "        Initial mass: " << GetMass() << " slug mass" << endl;
 845       cout << "        Initial weight: " << GetMass()*lbtoslug <<
 846         " lbs force" << endl;
 847       cout << "        Heat transfer: " << endl;
 848     }
 849   }
 850   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 851     if (from == 0) cout << "Instantiated: FGBallonet" << endl;
 852     if (from == 1) cout << "Destroyed:    FGBallonet" << endl;
 853   }
 854   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 855   }
 856   if (debug_lvl & 8 ) { // Runtime state variables
 857       cout << "        Ballonet holds " << Contents <<
 858         " mol air" << endl;
 859       cout << "        Temperature: " << Temperature << " Rankine" << endl;
 860       cout << "        Pressure: " << Pressure << " lbs/ft2" << endl;
 861       cout << "        Volume: " << Volume << " ft3" << endl;
 862       cout << "        Mass: " << GetMass() << " slug mass" << endl;
 863       cout << "        Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
 864   }
 865   if (debug_lvl & 16) { // Sanity checking
 866   }
 867   if (debug_lvl & 64) {
 868     if (from == 0) { // Constructor
 869       cout << IdSrc << endl;
 870       cout << IdHdr << endl;
 871     }
 872   }
 873 }
 874 }