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