1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 Author: Anders Gidenstam
5 Date started: 01/21/2006
7 ----- Copyright (C) 2006 - 2008 Anders Gidenstam (anders(at)gidenstam.org) --
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
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
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.
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.
26 FUNCTIONAL DESCRIPTION
27 --------------------------------------------------------------------------------
31 --------------------------------------------------------------------------------
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
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"
51 static const char *IdSrc = "$Id$";
52 static const char *IdHdr = ID_GASCELL;
54 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
58 const double FGGasCell::R = 3.4071; // [lbs ft/(mol Rankine)]
59 const double FGGasCell::M_air = 0.0019186; // [slug/mol]
60 const double FGGasCell::M_hydrogen = 0.00013841; // [slug/mol]
61 const double FGGasCell::M_helium = 0.00027409; // [slug/mol]
63 FGGasCell::FGGasCell(FGFDMExec* exec, Element* el, int num) : FGForce(exec)
68 Auxiliary = exec->GetAuxiliary();
69 Atmosphere = exec->GetAtmosphere();
70 PropertyManager = exec->GetPropertyManager();
71 Inertial = exec->GetInertial();
72 MassBalance = exec->GetMassBalance();
74 gasCellJ = FGMatrix33();
75 gasCellM = FGColumnVector3();
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;
83 // NOTE: In the local system X points north, Y points east and Z points down.
84 SetTransformType(FGForce::tLocalBody);
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;
92 element = el->FindElement("location");
94 vXYZ = element->FindElementTripletConvertTo("IN");
96 cerr << "Fatal Error: No location found for this gas cell." << endl;
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"))) {
103 if (el->FindElement("x_radius")) {
104 Xradius = el->FindElementValueAsNumberConvertTo("x_radius", "FT");
106 if (el->FindElement("y_radius")) {
107 Yradius = el->FindElementValueAsNumberConvertTo("y_radius", "FT");
109 if (el->FindElement("z_radius")) {
110 Zradius = el->FindElementValueAsNumberConvertTo("z_radius", "FT");
113 if (el->FindElement("x_width")) {
114 Xwidth = el->FindElementValueAsNumberConvertTo("x_width", "FT");
116 if (el->FindElement("y_width")) {
117 Ywidth = el->FindElementValueAsNumberConvertTo("y_width", "FT");
119 if (el->FindElement("z_width")) {
120 Zwidth = el->FindElementValueAsNumberConvertTo("z_width", "FT");
123 // The volume is a (potentially) extruded ellipsoid.
124 // However, currently only a few combinations of radius and width are
126 if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
127 (Xwidth == 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
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;
135 cerr << "Warning: Unsupported gas cell shape." << endl;
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);
147 cerr << "Fatal Error: Gas cell shape must be given." << endl;
150 if (el->FindElement("max_overpressure")) {
151 MaxOverpressure = el->FindElementValueAsNumberConvertTo("max_overpressure",
154 if (el->FindElement("fullness")) {
155 const double Fullness = el->FindElementValueAsNumber("fullness");
157 Volume = Fullness * MaxVolume;
159 cerr << "Warning: Invalid initial gas cell fullness value." << endl;
162 if (el->FindElement("valve_coefficient")) {
164 el->FindElementValueAsNumberConvertTo("valve_coefficient",
166 ValveCoefficient = max(ValveCoefficient, 0.0);
172 if (Temperature == 0.0) {
173 Temperature = Atmosphere->GetTemperature();
175 if (Pressure == 0.0) {
176 Pressure = Atmosphere->GetPressure();
179 // Calculate initial gas content.
180 Contents = Pressure * Volume / (R * Temperature);
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;
188 Pressure = max(IdealPressure, Pressure);
191 // Calculate initial gas content.
192 Contents = Pressure * MaxVolume / (R * Temperature);
195 Volume = Contents * R * Temperature / Pressure;
196 Mass = Contents * M_gas();
198 // Bind relevant properties
199 char property_name[80];
200 snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/max_volume-ft3",
202 PropertyManager->Tie( property_name, &MaxVolume );
203 PropertyManager->SetWritable( property_name, false );
204 snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/temp-R",
206 PropertyManager->Tie( property_name, &Temperature );
207 snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/pressure-psf",
209 PropertyManager->Tie( property_name, &Pressure );
210 snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/volume-ft3",
212 PropertyManager->Tie( property_name, &Volume );
213 snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/buoyancy-lbs",
215 PropertyManager->Tie( property_name, &Buoyancy );
216 snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/contents-mol",
218 PropertyManager->Tie( property_name, &Contents );
219 snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/valve_open",
221 PropertyManager->Tie( property_name, &ValveOpen );
225 // Read heat transfer coefficients
226 if (Element* heat = el->FindElement("heat")) {
227 Element* function_element = heat->FindElement("function");
228 while (function_element) {
229 HeatTransferCoeff.push_back(new FGFunction(PropertyManager,
231 function_element = heat->FindNextElement("function");
235 // Load ballonets if there are any
236 if (Element* ballonet_element = el->FindElement("ballonet")) {
237 while (ballonet_element) {
238 Ballonet.push_back(new FGBallonet(exec,
242 ballonet_element = el->FindNextElement("ballonet");
248 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
250 FGGasCell::~FGGasCell()
254 for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
255 HeatTransferCoeff.clear();
257 for (i = 0; i < Ballonet.size(); i++) delete Ballonet[i];
263 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
265 void FGGasCell::Calculate(double dt)
267 const double AirTemperature = Atmosphere->GetTemperature(); // [Rankine]
268 const double AirPressure = Atmosphere->GetPressure(); // [lbs/ft²]
269 const double AirDensity = Atmosphere->GetDensity(); // [slug/ft³]
270 const double g = Inertial->gravity(); // [lbs/slug]
272 const double OldTemperature = Temperature;
273 const double OldPressure = Pressure;
275 const unsigned int no_ballonets = Ballonet.size();
277 //-- Read ballonet state --
278 // NOTE: This model might need a more proper integration technique.
279 double BallonetsVolume = 0.0;
280 double BallonetsHeatFlow = 0.0;
281 for (i = 0; i < no_ballonets; i++) {
282 BallonetsVolume += Ballonet[i]->GetVolume();
283 BallonetsHeatFlow += Ballonet[i]->GetHeatFlow();
286 //-- Gas temperature --
288 if (HeatTransferCoeff.size() > 0) {
289 // The model is based on the ideal gas law.
290 // However, it does look a bit fishy. Please verify.
291 // dT/dt = dU / (Cv n R)
293 for (i = 0; i < HeatTransferCoeff.size(); i++) {
294 dU += HeatTransferCoeff[i]->GetValue();
296 // Don't include dt when accounting for adiabatic expansion/contraction.
297 // The rate of adiabatic cooling looks about right: ~5.4 Rankine/1000ft.
300 (dU * dt - Pressure * dVolumeIdeal - BallonetsHeatFlow) /
301 (Cv_gas() * Contents * R);
303 Temperature = AirTemperature;
306 // No simulation of complex temperature changes.
307 // Note: Making the gas cell behave adiabatically might be a better
309 Temperature = AirTemperature;
313 const double IdealPressure =
314 Contents * R * Temperature / (MaxVolume - BallonetsVolume);
315 if (IdealPressure > AirPressure + MaxOverpressure) {
316 Pressure = AirPressure + MaxOverpressure;
318 Pressure = max(IdealPressure, AirPressure);
321 //-- Manual valving --
323 // FIXME: Presently the effect of manual valving is computed using
324 // an ad hoc formula which might not be a good representation
326 if ((ValveCoefficient > 0.0) && (ValveOpen > 0.0)) {
327 // First compute the difference in pressure between the gas in the
328 // cell and the air above it.
329 // FixMe: CellHeight should depend on current volume.
330 const double CellHeight = 2 * Zradius + Zwidth; // [ft]
331 const double GasMass = Contents * M_gas(); // [slug]
332 const double GasVolume = Contents * R * Temperature / Pressure; // [ft³]
333 const double GasDensity = GasMass / GasVolume;
334 const double DeltaPressure =
335 Pressure + CellHeight * g * (AirDensity - GasDensity) - AirPressure;
336 const double VolumeValved =
337 ValveOpen * ValveCoefficient * DeltaPressure * dt;
339 max(0.0, Contents - Pressure * VolumeValved / (R * Temperature));
342 //-- Update ballonets. --
343 // Doing that here should give them the opportunity to react to the
345 BallonetsVolume = 0.0;
346 for (i = 0; i < no_ballonets; i++) {
347 Ballonet[i]->Calculate(dt);
348 BallonetsVolume += Ballonet[i]->GetVolume();
351 //-- Automatic safety valving. --
352 if (Contents * R * Temperature / (MaxVolume - BallonetsVolume) >
353 AirPressure + MaxOverpressure) {
354 // Gas is automatically valved. Valving capacity is assumed to be infinite.
355 // FIXME: This could/should be replaced by damage to the gas cell envelope.
357 (AirPressure + MaxOverpressure) *
358 (MaxVolume - BallonetsVolume) / (R * Temperature);
362 Volume = Contents * R * Temperature / Pressure + BallonetsVolume;
364 Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);
366 //-- Current buoyancy --
367 // The buoyancy is computed using the atmospheres local density.
368 Buoyancy = Volume * AirDensity * g;
370 // Note: This is gross buoyancy. The weight of the gas itself and
371 // any ballonets is not deducted here as the effects of the gas mass
372 // is handled by FGMassBalance.
373 vFn.InitMatrix(0.0, 0.0, - Buoyancy);
375 // Compute the inertia of the gas cell.
376 // Consider the gas cell as a shape of uniform density.
377 // FIXME: If the cell isn't ellipsoid or cylindrical the inertia will
379 gasCellJ = FGMatrix33();
380 const double mass = Contents * M_gas();
381 double Ixx, Iyy, Izz;
382 if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
383 (Xwidth == 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
385 Ixx = (1.0 / 5.0) * mass * (Yradius*Yradius + Zradius*Zradius);
386 Iyy = (1.0 / 5.0) * mass * (Xradius*Xradius + Zradius*Zradius);
387 Izz = (1.0 / 5.0) * mass * (Xradius*Xradius + Yradius*Yradius);
388 } else if ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
389 (Xwidth != 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
390 // Cylindrical volume (might not be valid with an elliptical cross-section).
391 Ixx = (1.0 / 2.0) * mass * Yradius * Zradius;
393 (1.0 / 4.0) * mass * Yradius * Zradius +
394 (1.0 / 12.0) * mass * Xwidth * Xwidth;
396 (1.0 / 4.0) * mass * Yradius * Zradius +
397 (1.0 / 12.0) * mass * Xwidth * Xwidth;
399 // Not supported. Revert to pointmass model.
400 Ixx = Iyy = Izz = 0.0;
402 // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
407 gasCellM.InitMatrix();
409 GetXYZ(eX) * Mass*slugtolb;
411 GetXYZ(eY) * Mass*slugtolb;
413 GetXYZ(eZ) * Mass*slugtolb;
415 if (no_ballonets > 0) {
416 // Add the mass, moment and inertia of any ballonets.
417 const FGColumnVector3 p = MassBalance->StructuralToBody( GetXYZ() );
419 for (i = 0; i < no_ballonets; i++) {
420 Mass += Ballonet[i]->GetMass();
422 // Add ballonet moments.
424 Ballonet[i]->GetXYZ(eX) * Ballonet[i]->GetMass()*slugtolb;
426 Ballonet[i]->GetXYZ(eY) * Ballonet[i]->GetMass()*slugtolb;
428 Ballonet[i]->GetXYZ(eZ) * Ballonet[i]->GetMass()*slugtolb;
430 // Moments of inertia must be converted to the gas cell frame here.
432 MassBalance->StructuralToBody( Ballonet[i]->GetXYZ() ) - p;
433 // Body basis is in FT.
434 const double mass = Ballonet[i]->GetMass();
435 gasCellJ += Ballonet[i]->GetInertia() +
436 FGMatrix33( 0, - mass*v(1)*v(2), - mass*v(1)*v(3),
437 - mass*v(2)*v(1), 0, - mass*v(2)*v(3),
438 - mass*v(3)*v(1), - mass*v(3)*v(2), 0 );
443 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
444 // The bitmasked value choices are as follows:
445 // unset: In this case (the default) JSBSim would only print
446 // out the normally expected messages, essentially echoing
447 // the config files as they are read. If the environment
448 // variable is not set, debug_lvl is set to 1 internally
449 // 0: This requests JSBSim not to output any messages
451 // 1: This value explicity requests the normal JSBSim
453 // 2: This value asks for a message to be printed out when
454 // a class is instantiated
455 // 4: When this value is set, a message is displayed when a
456 // FGModel object executes its Run() method
457 // 8: When this value is set, various runtime state variables
458 // are printed out periodically
459 // 16: When set various parameters are sanity checked and
460 // a message is printed out when they go out of bounds
462 void FGGasCell::Debug(int from)
464 if (debug_lvl <= 0) return;
466 if (debug_lvl & 1) { // Standard console startup message output
467 if (from == 0) { // Constructor
468 cout << " Gas cell holds " << Contents << " mol " <<
470 cout << " Cell location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
471 vXYZ(eY) << ", " << vXYZ(eZ) << endl;
472 cout << " Maximum volume: " << MaxVolume << " ft3" << endl;
473 cout << " Relief valve release pressure: " << MaxOverpressure <<
475 cout << " Manual valve coefficient: " << ValveCoefficient <<
476 " ft4*sec/slug" << endl;
477 cout << " Initial temperature: " << Temperature << " Rankine" <<
479 cout << " Initial pressure: " << Pressure << " lbs/ft2" << endl;
480 cout << " Initial volume: " << Volume << " ft3" << endl;
481 cout << " Initial mass: " << GetMass() << " slug mass" << endl;
482 cout << " Initial weight: " << GetMass()*lbtoslug << " lbs force" <<
484 cout << " Heat transfer: " << endl;
487 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
488 if (from == 0) cout << "Instantiated: FGGasCell" << endl;
489 if (from == 1) cout << "Destroyed: FGGasCell" << endl;
491 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
493 if (debug_lvl & 8 ) { // Runtime state variables
494 cout << " " << type << " cell holds " << Contents << " mol " << endl;
495 cout << " Temperature: " << Temperature << " Rankine" << endl;
496 cout << " Pressure: " << Pressure << " lbs/ft2" << endl;
497 cout << " Volume: " << Volume << " ft3" << endl;
498 cout << " Mass: " << GetMass() << " slug mass" << endl;
499 cout << " Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
501 if (debug_lvl & 16) { // Sanity checking
503 if (debug_lvl & 64) {
504 if (from == 0) { // Constructor
505 cout << IdSrc << endl;
506 cout << IdHdr << endl;
511 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
513 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
514 const double FGBallonet::R = 3.4071; // [lbs ft/(mol Rankine)]
515 const double FGBallonet::M_air = 0.0019186; // [slug/mol]
516 const double FGBallonet::Cv_air = 5.0/2.0; // [??]
518 FGBallonet::FGBallonet(FGFDMExec* exec, Element* el, int num, FGGasCell* parent)
523 Auxiliary = exec->GetAuxiliary();
524 Atmosphere = exec->GetAtmosphere();
525 PropertyManager = exec->GetPropertyManager();
526 Inertial = exec->GetInertial();
528 ballonetJ = FGMatrix33();
530 MaxVolume = MaxOverpressure = Temperature = Pressure =
531 Contents = Volume = dVolumeIdeal = dU = 0.0;
532 Xradius = Yradius = Zradius = Xwidth = Ywidth = Zwidth = 0.0;
533 ValveCoefficient = ValveOpen = 0.0;
538 // NOTE: In the local system X points north, Y points east and Z points down.
539 element = el->FindElement("location");
541 vXYZ = element->FindElementTripletConvertTo("IN");
543 cerr << "Fatal Error: No location found for this ballonet." << endl;
546 if ((el->FindElement("x_radius") || el->FindElement("x_width")) &&
547 (el->FindElement("y_radius") || el->FindElement("y_width")) &&
548 (el->FindElement("z_radius") || el->FindElement("z_width"))) {
550 if (el->FindElement("x_radius")) {
551 Xradius = el->FindElementValueAsNumberConvertTo("x_radius", "FT");
553 if (el->FindElement("y_radius")) {
554 Yradius = el->FindElementValueAsNumberConvertTo("y_radius", "FT");
556 if (el->FindElement("z_radius")) {
557 Zradius = el->FindElementValueAsNumberConvertTo("z_radius", "FT");
560 if (el->FindElement("x_width")) {
561 Xwidth = el->FindElementValueAsNumberConvertTo("x_width", "FT");
563 if (el->FindElement("y_width")) {
564 Ywidth = el->FindElementValueAsNumberConvertTo("y_width", "FT");
566 if (el->FindElement("z_width")) {
567 Zwidth = el->FindElementValueAsNumberConvertTo("z_width", "FT");
570 // The volume is a (potentially) extruded ellipsoid.
571 // FIXME: However, currently only a few combinations of radius and
572 // width are fully supported.
573 if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
574 (Xwidth == 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
576 MaxVolume = 4.0 * M_PI * Xradius * Yradius * Zradius / 3.0;
577 } else if ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
578 (Xwidth != 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
579 // Cylindrical volume.
580 MaxVolume = M_PI * Yradius * Zradius * Xwidth;
582 cerr << "Warning: Unsupported ballonet shape." << endl;
584 (4.0 * M_PI * Xradius * Yradius * Zradius / 3.0 +
585 M_PI * Yradius * Zradius * Xwidth +
586 M_PI * Xradius * Zradius * Ywidth +
587 M_PI * Xradius * Yradius * Zwidth +
588 2.0 * Xradius * Ywidth * Zwidth +
589 2.0 * Yradius * Xwidth * Zwidth +
590 2.0 * Zradius * Xwidth * Ywidth +
591 Xwidth * Ywidth * Zwidth);
594 cerr << "Fatal Error: Ballonet shape must be given." << endl;
597 if (el->FindElement("max_overpressure")) {
598 MaxOverpressure = el->FindElementValueAsNumberConvertTo("max_overpressure",
601 if (el->FindElement("fullness")) {
602 const double Fullness = el->FindElementValueAsNumber("fullness");
604 Volume = Fullness * MaxVolume;
606 cerr << "Warning: Invalid initial ballonet fullness value." << endl;
609 if (el->FindElement("valve_coefficient")) {
611 el->FindElementValueAsNumberConvertTo("valve_coefficient",
613 ValveCoefficient = max(ValveCoefficient, 0.0);
617 if (Temperature == 0.0) {
618 Temperature = Parent->GetTemperature();
620 if (Pressure == 0.0) {
621 Pressure = Parent->GetPressure();
624 // Calculate initial air content.
625 Contents = Pressure * Volume / (R * Temperature);
627 // Clip to max allowed value.
628 const double IdealPressure = Contents * R * Temperature / MaxVolume;
629 if (IdealPressure > Pressure + MaxOverpressure) {
630 Contents = (Pressure + MaxOverpressure) * MaxVolume / (R * Temperature);
631 Pressure = Pressure + MaxOverpressure;
633 Pressure = max(IdealPressure, Pressure);
636 // Calculate initial air content.
637 Contents = Pressure * MaxVolume / (R * Temperature);
640 Volume = Contents * R * Temperature / Pressure;
642 // Bind relevant properties
643 char property_name[80];
644 snprintf(property_name, 80,
645 "buoyant_forces/gas-cell[%d]/ballonet[%d]/max_volume-ft3",
648 PropertyManager->Tie( property_name, &MaxVolume );
649 PropertyManager->SetWritable( property_name, false );
650 snprintf(property_name, 80,
651 "buoyant_forces/gas-cell[%d]/ballonet[%d]/temp-R",
654 PropertyManager->Tie( property_name, &Temperature );
655 snprintf(property_name, 80,
656 "buoyant_forces/gas-cell[%d]/ballonet[%d]/pressure-psf",
659 PropertyManager->Tie( property_name, &Pressure );
660 snprintf(property_name, 80,
661 "buoyant_forces/gas-cell[%d]/ballonet[%d]/volume-ft3",
664 PropertyManager->Tie( property_name, &Volume );
665 snprintf(property_name, 80,
666 "buoyant_forces/gas-cell[%d]/ballonet[%d]/contents-mol",
669 PropertyManager->Tie( property_name, &Contents );
670 snprintf(property_name, 80,
671 "buoyant_forces/gas-cell[%d]/ballonet[%d]/valve_open",
674 PropertyManager->Tie( property_name, &ValveOpen );
678 // Read heat transfer coefficients
679 if (Element* heat = el->FindElement("heat")) {
680 Element* function_element = heat->FindElement("function");
681 while (function_element) {
682 HeatTransferCoeff.push_back(new FGFunction(PropertyManager,
684 function_element = heat->FindNextElement("function");
687 // Read blower input function
688 if (Element* blower = el->FindElement("blower_input")) {
689 Element* function_element = blower->FindElement("function");
690 BlowerInput = new FGFunction(PropertyManager,
695 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
697 FGBallonet::~FGBallonet()
701 for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
702 HeatTransferCoeff.clear();
710 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
712 void FGBallonet::Calculate(double dt)
714 const double ParentPressure = Parent->GetPressure(); // [lbs/ft²]
715 const double AirPressure = Atmosphere->GetPressure(); // [lbs/ft²]
717 const double OldTemperature = Temperature;
718 const double OldPressure = Pressure;
721 //-- Gas temperature --
723 // The model is based on the ideal gas law.
724 // However, it does look a bit fishy. Please verify.
725 // dT/dt = dU / (Cv n R)
727 for (i = 0; i < HeatTransferCoeff.size(); i++) {
728 dU += HeatTransferCoeff[i]->GetValue();
730 // dt is already accounted for in dVolumeIdeal.
733 (dU * dt - Pressure * dVolumeIdeal) / (Cv_air * Contents * R);
735 Temperature = Parent->GetTemperature();
739 const double IdealPressure = Contents * R * Temperature / MaxVolume;
740 // The pressure is at least that of the parent gas cell.
741 Pressure = max(IdealPressure, ParentPressure);
745 const double AddedVolume = BlowerInput->GetValue() * dt;
746 if (AddedVolume > 0.0) {
747 Contents += Pressure * AddedVolume / (R * Temperature);
751 //-- Pressure relief and manual valving --
752 // FIXME: Presently the effect of valving is computed using
753 // an ad hoc formula which might not be a good representation
755 if ((ValveCoefficient > 0.0) &&
756 ((ValveOpen > 0.0) || (Pressure > AirPressure + MaxOverpressure))) {
757 const double DeltaPressure = Pressure - AirPressure;
758 const double VolumeValved =
759 ((Pressure > AirPressure + MaxOverpressure) ? 1.0 : ValveOpen) *
760 ValveCoefficient * DeltaPressure * dt;
761 // FIXME: Too small values of Contents sometimes leads to NaN.
762 // Currently the minimum is restricted to a safe value.
764 max(1.0, Contents - Pressure * VolumeValved / (R * Temperature));
768 Volume = Contents * R * Temperature / Pressure;
770 Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);
772 // Compute the inertia of the ballonet.
773 // Consider the ballonet as a shape of uniform density.
774 // FIXME: If the ballonet isn't ellipsoid or cylindrical the inertia will
776 ballonetJ = FGMatrix33();
777 const double mass = Contents * M_air;
778 double Ixx, Iyy, Izz;
779 if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
780 (Xwidth == 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
782 Ixx = (1.0 / 5.0) * mass * (Yradius*Yradius + Zradius*Zradius);
783 Iyy = (1.0 / 5.0) * mass * (Xradius*Xradius + Zradius*Zradius);
784 Izz = (1.0 / 5.0) * mass * (Xradius*Xradius + Yradius*Yradius);
785 } else if ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
786 (Xwidth != 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
787 // Cylindrical volume (might not be valid with an elliptical cross-section).
788 Ixx = (1.0 / 2.0) * mass * Yradius * Zradius;
790 (1.0 / 4.0) * mass * Yradius * Zradius +
791 (1.0 / 12.0) * mass * Xwidth * Xwidth;
793 (1.0 / 4.0) * mass * Yradius * Zradius +
794 (1.0 / 12.0) * mass * Xwidth * Xwidth;
796 // Not supported. Revert to pointmass model.
797 Ixx = Iyy = Izz = 0.0;
799 // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
800 ballonetJ(1,1) = Ixx;
801 ballonetJ(2,2) = Iyy;
802 ballonetJ(3,3) = Izz;
805 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
806 // The bitmasked value choices are as follows:
807 // unset: In this case (the default) JSBSim would only print
808 // out the normally expected messages, essentially echoing
809 // the config files as they are read. If the environment
810 // variable is not set, debug_lvl is set to 1 internally
811 // 0: This requests JSBSim not to output any messages
813 // 1: This value explicity requests the normal JSBSim
815 // 2: This value asks for a message to be printed out when
816 // a class is instantiated
817 // 4: When this value is set, a message is displayed when a
818 // FGModel object executes its Run() method
819 // 8: When this value is set, various runtime state variables
820 // are printed out periodically
821 // 16: When set various parameters are sanity checked and
822 // a message is printed out when they go out of bounds
824 void FGBallonet::Debug(int from)
826 if (debug_lvl <= 0) return;
828 if (debug_lvl & 1) { // Standard console startup message output
829 if (from == 0) { // Constructor
830 cout << " Ballonet holds " << Contents << " mol air" << endl;
831 cout << " Location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
832 vXYZ(eY) << ", " << vXYZ(eZ) << endl;
833 cout << " Maximum volume: " << MaxVolume << " ft3" << endl;
834 cout << " Relief valve release pressure: " << MaxOverpressure <<
836 cout << " Relief valve coefficient: " << ValveCoefficient <<
837 " ft4*sec/slug" << endl;
838 cout << " Initial temperature: " << Temperature << " Rankine" <<
840 cout << " Initial pressure: " << Pressure << " lbs/ft2" << endl;
841 cout << " Initial volume: " << Volume << " ft3" << endl;
842 cout << " Initial mass: " << GetMass() << " slug mass" << endl;
843 cout << " Initial weight: " << GetMass()*lbtoslug <<
844 " lbs force" << endl;
845 cout << " Heat transfer: " << endl;
848 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
849 if (from == 0) cout << "Instantiated: FGBallonet" << endl;
850 if (from == 1) cout << "Destroyed: FGBallonet" << endl;
852 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
854 if (debug_lvl & 8 ) { // Runtime state variables
855 cout << " Ballonet holds " << Contents <<
857 cout << " Temperature: " << Temperature << " Rankine" << endl;
858 cout << " Pressure: " << Pressure << " lbs/ft2" << endl;
859 cout << " Volume: " << Volume << " ft3" << endl;
860 cout << " Mass: " << GetMass() << " slug mass" << endl;
861 cout << " Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
863 if (debug_lvl & 16) { // Sanity checking
865 if (debug_lvl & 64) {
866 if (from == 0) { // Constructor
867 cout << IdSrc << endl;
868 cout << IdHdr << endl;