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