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