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