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 - 2011  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/FGMassBalance.h"
  40 #include "FGGasCell.h"
  41 #include "input_output/FGXMLElement.h"
  42 #include <iostream>
  43 #include <cstdlib>
  44
  45 using std::cerr;
  46 using std::endl;
  47 using std::cout;
  48 using std::string;
  49 using std::max;
  50
  51 namespace JSBSim {
  52
  53 static const char *IdSrc = "$Id: FGGasCell.cpp,v 1.15 2011/08/06 13:47:59 jberndt Exp $";
  54 static const char *IdHdr = ID_GASCELL;
  55
  56 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  57 CLASS IMPLEMENTATION
  58 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  59 /* Constants. */
  60 const double FGGasCell::R = 3.4071;              // [lbs ft/(mol Rankine)]
  61 const double FGGasCell::M_air = 0.0019186;       // [slug/mol]
  62 const double FGGasCell::M_hydrogen = 0.00013841; // [slug/mol]
  63 const double FGGasCell::M_helium = 0.00027409;   // [slug/mol]
  64
  65 FGGasCell::FGGasCell(FGFDMExec* exec, Element* el, int num, const struct Inputs& input)
  66   : FGForce(exec), in(input)
  67 {
  68   string token;
  69   Element* element;
  70
  71   PropertyManager = exec->GetPropertyManager();
  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 = in.Temperature;
 174   }
 175   if (Pressure == 0.0) {
 176     Pressure = in.Pressure;
 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, false );
 205   PropertyManager->SetWritable( property_name, false );
 206   property_name = base_property_name + "/temp-R";
 207   PropertyManager->Tie( property_name.c_str(), &Temperature, false );
 208   property_name = base_property_name + "/pressure-psf";
 209   PropertyManager->Tie( property_name.c_str(), &Pressure, false );
 210   property_name = base_property_name + "/volume-ft3";
 211   PropertyManager->Tie( property_name.c_str(), &Volume, false );
 212   property_name = base_property_name + "/buoyancy-lbs";
 213   PropertyManager->Tie( property_name.c_str(), &Buoyancy, false );
 214   property_name = base_property_name + "/contents-mol";
 215   PropertyManager->Tie( property_name.c_str(), &Contents, false );
 216   property_name = base_property_name + "/valve_open";
 217   PropertyManager->Tie( property_name.c_str(), &ValveOpen, false );
 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, in));
 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 = in.Temperature;  // [Rankine]
 264   const double AirPressure    = in.Pressure;     // [lbs/ft^2]
 265   const double AirDensity     = in.Density;      // [slug/ft^3]
 266   const double g = in.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   // Transform the moments of inertia to the body frame.
 404   gasCellJ += MassBalance->GetPointmassInertia(Mass, GetXYZ());
 405
 406   gasCellM.InitMatrix();
 407   gasCellM(eX) +=
 408     GetXYZ(eX) * Mass*slugtolb;
 409   gasCellM(eY) +=
 410     GetXYZ(eY) * Mass*slugtolb;
 411   gasCellM(eZ) +=
 412     GetXYZ(eZ) * Mass*slugtolb;
 413
 414   if (no_ballonets > 0) {
 415     // Add the mass, moment and inertia of any ballonets.
 416     for (i = 0; i < no_ballonets; i++) {
 417       Mass += Ballonet[i]->GetMass();
 418
 419       // Add ballonet moments due to mass (in the structural frame).
 420       gasCellM(eX) +=
 421         Ballonet[i]->GetXYZ(eX) * Ballonet[i]->GetMass()*slugtolb;
 422       gasCellM(eY) +=
 423         Ballonet[i]->GetXYZ(eY) * Ballonet[i]->GetMass()*slugtolb;
 424       gasCellM(eZ) +=
 425         Ballonet[i]->GetXYZ(eZ) * Ballonet[i]->GetMass()*slugtolb;
 426
 427       gasCellJ += Ballonet[i]->GetInertia();
 428     }
 429   }
 430 }
 431
 432 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 433 //    The bitmasked value choices are as follows:
 434 //    unset: In this case (the default) JSBSim would only print
 435 //       out the normally expected messages, essentially echoing
 436 //       the config files as they are read. If the environment
 437 //       variable is not set, debug_lvl is set to 1 internally
 438 //    0: This requests JSBSim not to output any messages
 439 //       whatsoever.
 440 //    1: This value explicity requests the normal JSBSim
 441 //       startup messages
 442 //    2: This value asks for a message to be printed out when
 443 //       a class is instantiated
 444 //    4: When this value is set, a message is displayed when a
 445 //       FGModel object executes its Run() method
 446 //    8: When this value is set, various runtime state variables
 447 //       are printed out periodically
 448 //    16: When set various parameters are sanity checked and
 449 //       a message is printed out when they go out of bounds
 450
 451 void FGGasCell::Debug(int from)
 452 {
 453   if (debug_lvl <= 0) return;
 454
 455   if (debug_lvl & 1) { // Standard console startup message output
 456     if (from == 0) { // Constructor
 457       cout << "    Gas cell holds " << Contents << " mol " <<
 458         type << endl;
 459       cout << "      Cell location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
 460         vXYZ(eY) << ", " << vXYZ(eZ) << endl;
 461       cout << "      Maximum volume: " << MaxVolume << " ft3" << endl;
 462       cout << "      Relief valve release pressure: " << MaxOverpressure <<
 463         " lbs/ft2" << endl;
 464       cout << "      Manual valve coefficient: " << ValveCoefficient <<
 465         " ft4*sec/slug" << endl;
 466       cout << "      Initial temperature: " << Temperature << " Rankine" <<
 467         endl;
 468       cout << "      Initial pressure: " << Pressure << " lbs/ft2" << endl;
 469       cout << "      Initial volume: " << Volume << " ft3" << endl;
 470       cout << "      Initial mass: " << GetMass() << " slug mass" << endl;
 471       cout << "      Initial weight: " << GetMass()*lbtoslug << " lbs force" <<
 472         endl;
 473       cout << "      Heat transfer: " << endl;
 474     }
 475   }
 476   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 477     if (from == 0) cout << "Instantiated: FGGasCell" << endl;
 478     if (from == 1) cout << "Destroyed:    FGGasCell" << endl;
 479   }
 480   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 481   }
 482   if (debug_lvl & 8 ) { // Runtime state variables
 483       cout << "      " << type << " cell holds " << Contents << " mol " << endl;
 484       cout << "      Temperature: " << Temperature << " Rankine" << endl;
 485       cout << "      Pressure: " << Pressure << " lbs/ft2" << endl;
 486       cout << "      Volume: " << Volume << " ft3" << endl;
 487       cout << "      Mass: " << GetMass() << " slug mass" << endl;
 488       cout << "      Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
 489   }
 490   if (debug_lvl & 16) { // Sanity checking
 491   }
 492   if (debug_lvl & 64) {
 493     if (from == 0) { // Constructor
 494       cout << IdSrc << endl;
 495       cout << IdHdr << endl;
 496     }
 497   }
 498 }
 499
 500 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 501 CLASS IMPLEMENTATION
 502 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 503 const double FGBallonet::R = 3.4071;              // [lbs ft/(mol Rankine)]
 504 const double FGBallonet::M_air = 0.0019186;       // [slug/mol]
 505 const double FGBallonet::Cv_air = 5.0/2.0;        // [??]
 506
 507 FGBallonet::FGBallonet(FGFDMExec* exec, Element* el, int num, FGGasCell* parent, const struct FGGasCell::Inputs& input)
 508   : in(input)
 509 {
 510   string token;
 511   Element* element;
 512
 513   PropertyManager = exec->GetPropertyManager();
 514   MassBalance = exec->GetMassBalance();
 515
 516   ballonetJ = FGMatrix33();
 517
 518   MaxVolume = MaxOverpressure = Temperature = Pressure =
 519     Contents = Volume = dVolumeIdeal = dU = 0.0;
 520   Xradius = Yradius = Zradius = Xwidth = Ywidth = Zwidth = 0.0;
 521   ValveCoefficient = ValveOpen = 0.0;
 522   BlowerInput = NULL;
 523   CellNum = num;
 524   Parent = parent;
 525
 526   // NOTE: In the local system X points north, Y points east and Z points down.
 527   element = el->FindElement("location");
 528   if (element) {
 529     vXYZ = element->FindElementTripletConvertTo("IN");
 530   } else {
 531     cerr << "Fatal Error: No location found for this ballonet." << endl;
 532     exit(-1);
 533   }
 534   if ((el->FindElement("x_radius") || el->FindElement("x_width")) &&
 535       (el->FindElement("y_radius") || el->FindElement("y_width")) &&
 536       (el->FindElement("z_radius") || el->FindElement("z_width"))) {
 537
 538     if (el->FindElement("x_radius")) {
 539       Xradius = el->FindElementValueAsNumberConvertTo("x_radius", "FT");
 540     }
 541     if (el->FindElement("y_radius")) {
 542       Yradius = el->FindElementValueAsNumberConvertTo("y_radius", "FT");
 543     }
 544     if (el->FindElement("z_radius")) {
 545       Zradius = el->FindElementValueAsNumberConvertTo("z_radius", "FT");
 546     }
 547
 548     if (el->FindElement("x_width")) {
 549       Xwidth = el->FindElementValueAsNumberConvertTo("x_width", "FT");
 550     }
 551     if (el->FindElement("y_width")) {
 552       Ywidth = el->FindElementValueAsNumberConvertTo("y_width", "FT");
 553     }
 554     if (el->FindElement("z_width")) {
 555       Zwidth = el->FindElementValueAsNumberConvertTo("z_width", "FT");
 556     }
 557
 558     // The volume is a (potentially) extruded ellipsoid.
 559     // FIXME: However, currently only a few combinations of radius and
 560     //        width are fully supported.
 561     if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 562         (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 563       // Ellipsoid volume.
 564       MaxVolume = 4.0  * M_PI * Xradius * Yradius * Zradius / 3.0;
 565     } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 566                 (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 567       // Cylindrical volume.
 568       MaxVolume = M_PI * Yradius * Zradius * Xwidth;
 569     } else {
 570       cerr << "Warning: Unsupported ballonet shape." << endl;
 571       MaxVolume =
 572         (4.0  * M_PI * Xradius * Yradius * Zradius / 3.0 +
 573          M_PI * Yradius * Zradius * Xwidth +
 574          M_PI * Xradius * Zradius * Ywidth +
 575          M_PI * Xradius * Yradius * Zwidth +
 576          2.0  * Xradius * Ywidth * Zwidth +
 577          2.0  * Yradius * Xwidth * Zwidth +
 578          2.0  * Zradius * Xwidth * Ywidth +
 579          Xwidth * Ywidth * Zwidth);
 580     }
 581   } else {
 582     cerr << "Fatal Error: Ballonet shape must be given." << endl;
 583     exit(-1);
 584   }
 585   if (el->FindElement("max_overpressure")) {
 586     MaxOverpressure = el->FindElementValueAsNumberConvertTo("max_overpressure",
 587                                                             "LBS/FT2");
 588   }
 589   if (el->FindElement("fullness")) {
 590     const double Fullness = el->FindElementValueAsNumber("fullness");
 591     if (0 <= Fullness) {
 592       Volume = Fullness * MaxVolume;
 593     } else {
 594       cerr << "Warning: Invalid initial ballonet fullness value." << endl;
 595     }
 596   }
 597   if (el->FindElement("valve_coefficient")) {
 598     ValveCoefficient =
 599       el->FindElementValueAsNumberConvertTo("valve_coefficient",
 600                                             "FT4*SEC/SLUG");
 601     ValveCoefficient = max(ValveCoefficient, 0.0);
 602   }
 603
 604   // Initialize state
 605   if (Temperature == 0.0) {
 606     Temperature = Parent->GetTemperature();
 607   }
 608   if (Pressure == 0.0) {
 609     Pressure = Parent->GetPressure();
 610   }
 611   if (Volume != 0.0) {
 612     // Calculate initial air content.
 613     Contents = Pressure * Volume / (R * Temperature);
 614
 615     // Clip to max allowed value.
 616     const double IdealPressure = Contents * R * Temperature / MaxVolume;
 617     if (IdealPressure > Pressure + MaxOverpressure) {
 618       Contents = (Pressure + MaxOverpressure) * MaxVolume / (R * Temperature);
 619       Pressure = Pressure + MaxOverpressure;
 620     } else {
 621       Pressure = max(IdealPressure, Pressure);
 622     }
 623   } else {
 624     // Calculate initial air content.
 625     Contents = Pressure * MaxVolume / (R * Temperature);
 626   }
 627
 628   Volume = Contents * R * Temperature / Pressure;
 629
 630   // Bind relevant properties
 631   string property_name, base_property_name;
 632   base_property_name = CreateIndexedPropertyName("buoyant_forces/gas-cell", Parent->GetIndex());
 633   base_property_name = CreateIndexedPropertyName(base_property_name + "/ballonet", CellNum);
 634
 635   property_name = base_property_name + "/max_volume-ft3";
 636   PropertyManager->Tie( property_name, &MaxVolume, false );
 637   PropertyManager->SetWritable( property_name, false );
 638
 639   property_name = base_property_name + "/temp-R";
 640   PropertyManager->Tie( property_name, &Temperature, false );
 641
 642   property_name = base_property_name + "/pressure-psf";
 643   PropertyManager->Tie( property_name, &Pressure, false );
 644
 645   property_name = base_property_name + "/volume-ft3";
 646   PropertyManager->Tie( property_name, &Volume, false );
 647
 648   property_name = base_property_name + "/contents-mol";
 649   PropertyManager->Tie( property_name, &Contents, false );
 650
 651   property_name = base_property_name + "/valve_open";
 652   PropertyManager->Tie( property_name, &ValveOpen, false );
 653
 654   Debug(0);
 655
 656   // Read heat transfer coefficients
 657   if (Element* heat = el->FindElement("heat")) {
 658     Element* function_element = heat->FindElement("function");
 659     while (function_element) {
 660       HeatTransferCoeff.push_back(new FGFunction(PropertyManager,
 661                                                  function_element));
 662       function_element = heat->FindNextElement("function");
 663     }
 664   }
 665   // Read blower input function
 666   if (Element* blower = el->FindElement("blower_input")) {
 667     Element* function_element = blower->FindElement("function");
 668     BlowerInput = new FGFunction(PropertyManager,
 669                                  function_element);
 670   }
 671 }
 672
 673 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 674
 675 FGBallonet::~FGBallonet()
 676 {
 677   unsigned int i;
 678
 679   for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
 680   HeatTransferCoeff.clear();
 681
 682   delete BlowerInput;
 683   BlowerInput = NULL;
 684
 685   Debug(1);
 686 }
 687
 688 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 689
 690 void FGBallonet::Calculate(double dt)
 691 {
 692   const double ParentPressure = Parent->GetPressure(); // [lbs/ft^2]
 693   const double AirPressure    = in.Pressure;           // [lbs/ft^2]
 694
 695   const double OldTemperature = Temperature;
 696   const double OldPressure    = Pressure;
 697   unsigned int i;
 698
 699   //-- Gas temperature --
 700
 701   // The model is based on the ideal gas law.
 702   // However, it does look a bit fishy. Please verify.
 703   //   dT/dt = dU / (Cv n R)
 704   dU = 0.0;
 705   for (i = 0; i < HeatTransferCoeff.size(); i++) {
 706     dU += HeatTransferCoeff[i]->GetValue();
 707   }
 708   // dt is already accounted for in dVolumeIdeal.
 709   if (Contents > 0) {
 710     Temperature +=
 711       (dU * dt - Pressure * dVolumeIdeal) / (Cv_air * Contents * R);
 712   } else {
 713     Temperature = Parent->GetTemperature();
 714   }
 715
 716   //-- Pressure --
 717   const double IdealPressure = Contents * R * Temperature / MaxVolume;
 718   // The pressure is at least that of the parent gas cell.
 719   Pressure = max(IdealPressure, ParentPressure);
 720
 721   //-- Blower input --
 722   if (BlowerInput) {
 723     const double AddedVolume = BlowerInput->GetValue() * dt;
 724     if (AddedVolume > 0.0) {
 725       Contents += Pressure * AddedVolume / (R * Temperature);
 726     }
 727   }
 728
 729   //-- Pressure relief and manual valving --
 730   // FIXME: Presently the effect of valving is computed using
 731   //        an ad hoc formula which might not be a good representation
 732   //        of reality.
 733   if ((ValveCoefficient > 0.0) &&
 734       ((ValveOpen > 0.0) || (Pressure > AirPressure + MaxOverpressure))) {
 735     const double DeltaPressure = Pressure - AirPressure;
 736     const double VolumeValved =
 737       ((Pressure > AirPressure + MaxOverpressure) ? 1.0 : ValveOpen) *
 738       ValveCoefficient * DeltaPressure * dt;
 739     // FIXME: Too small values of Contents sometimes leads to NaN.
 740     //        Currently the minimum is restricted to a safe value.
 741     Contents =
 742       max(1.0, Contents - Pressure * VolumeValved / (R * Temperature));
 743   }
 744
 745   //-- Volume --
 746   Volume = Contents * R * Temperature / Pressure;
 747   dVolumeIdeal =
 748     Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);
 749
 750   // Compute the inertia of the ballonet.
 751   // Consider the ballonet as a shape of uniform density.
 752   // FIXME: If the ballonet isn't ellipsoid or cylindrical the inertia will
 753   //        be wrong.
 754   ballonetJ = FGMatrix33();
 755   const double mass = Contents * M_air;
 756   double Ixx, Iyy, Izz;
 757   if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 758       (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 759     // Ellipsoid volume.
 760     Ixx = (1.0 / 5.0) * mass * (Yradius*Yradius + Zradius*Zradius);
 761     Iyy = (1.0 / 5.0) * mass * (Xradius*Xradius + Zradius*Zradius);
 762     Izz = (1.0 / 5.0) * mass * (Xradius*Xradius + Yradius*Yradius);
 763   } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
 764               (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
 765     // Cylindrical volume (might not be valid with an elliptical cross-section).
 766     Ixx = (1.0 / 2.0) * mass * Yradius * Zradius;
 767     Iyy =
 768       (1.0 / 4.0) * mass * Yradius * Zradius +
 769       (1.0 / 12.0) * mass * Xwidth * Xwidth;
 770     Izz =
 771       (1.0 / 4.0) * mass * Yradius * Zradius +
 772       (1.0 / 12.0) * mass * Xwidth * Xwidth;
 773   } else {
 774     // Not supported. Revert to pointmass model.
 775     Ixx = Iyy = Izz = 0.0;
 776   }
 777   // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
 778   ballonetJ(1,1) = Ixx;
 779   ballonetJ(2,2) = Iyy;
 780   ballonetJ(3,3) = Izz;
 781   // Transform the moments of inertia to the body frame.
 782   ballonetJ += MassBalance->GetPointmassInertia(GetMass(), GetXYZ());
 783 }
 784
 785 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 786 //    The bitmasked value choices are as follows:
 787 //    unset: In this case (the default) JSBSim would only print
 788 //       out the normally expected messages, essentially echoing
 789 //       the config files as they are read. If the environment
 790 //       variable is not set, debug_lvl is set to 1 internally
 791 //    0: This requests JSBSim not to output any messages
 792 //       whatsoever.
 793 //    1: This value explicity requests the normal JSBSim
 794 //       startup messages
 795 //    2: This value asks for a message to be printed out when
 796 //       a class is instantiated
 797 //    4: When this value is set, a message is displayed when a
 798 //       FGModel object executes its Run() method
 799 //    8: When this value is set, various runtime state variables
 800 //       are printed out periodically
 801 //    16: When set various parameters are sanity checked and
 802 //       a message is printed out when they go out of bounds
 803
 804 void FGBallonet::Debug(int from)
 805 {
 806   if (debug_lvl <= 0) return;
 807
 808   if (debug_lvl & 1) { // Standard console startup message output
 809     if (from == 0) { // Constructor
 810       cout << "      Ballonet holds " << Contents << " mol air" << endl;
 811       cout << "        Location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
 812         vXYZ(eY) << ", " << vXYZ(eZ) << endl;
 813       cout << "        Maximum volume: " << MaxVolume << " ft3" << endl;
 814       cout << "        Relief valve release pressure: " << MaxOverpressure <<
 815         " lbs/ft2" << endl;
 816       cout << "        Relief valve coefficient: " << ValveCoefficient <<
 817         " ft4*sec/slug" << endl;
 818       cout << "        Initial temperature: " << Temperature << " Rankine" <<
 819         endl;
 820       cout << "        Initial pressure: " << Pressure << " lbs/ft2" << endl;
 821       cout << "        Initial volume: " << Volume << " ft3" << endl;
 822       cout << "        Initial mass: " << GetMass() << " slug mass" << endl;
 823       cout << "        Initial weight: " << GetMass()*lbtoslug <<
 824         " lbs force" << endl;
 825       cout << "        Heat transfer: " << endl;
 826     }
 827   }
 828   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 829     if (from == 0) cout << "Instantiated: FGBallonet" << endl;
 830     if (from == 1) cout << "Destroyed:    FGBallonet" << endl;
 831   }
 832   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 833   }
 834   if (debug_lvl & 8 ) { // Runtime state variables
 835       cout << "        Ballonet holds " << Contents <<
 836         " mol air" << endl;
 837       cout << "        Temperature: " << Temperature << " Rankine" << endl;
 838       cout << "        Pressure: " << Pressure << " lbs/ft2" << endl;
 839       cout << "        Volume: " << Volume << " ft3" << endl;
 840       cout << "        Mass: " << GetMass() << " slug mass" << endl;
 841       cout << "        Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
 842   }
 843   if (debug_lvl & 16) { // Sanity checking
 844   }
 845   if (debug_lvl & 64) {
 846     if (from == 0) { // Constructor
 847       cout << IdSrc << endl;
 848       cout << IdHdr << endl;
 849     }
 850   }
 851 }
 852 }