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 string property_name, base_property_name;
201 base_property_name = CreateIndexedPropertyName("buoyant_forces/gas-cell", CellNum);
203 property_name = base_property_name + "/max_volume-ft3";
204 PropertyManager->Tie( property_name.c_str(), &MaxVolume );
205 PropertyManager->SetWritable( property_name, false );
206 property_name = base_property_name + "/temp-R";
207 PropertyManager->Tie( property_name.c_str(), &Temperature );
208 property_name = base_property_name + "/pressure-psf";
209 PropertyManager->Tie( property_name.c_str(), &Pressure );
210 property_name = base_property_name + "/volume-ft3";
211 PropertyManager->Tie( property_name.c_str(), &Volume );
212 property_name = base_property_name + "/buoyancy-lbs";
213 PropertyManager->Tie( property_name.c_str(), &Buoyancy );
214 property_name = base_property_name + "/contents-mol";
215 PropertyManager->Tie( property_name.c_str(), &Contents );
216 property_name = base_property_name + "/valve_open";
217 PropertyManager->Tie( property_name.c_str(), &ValveOpen );
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,
227 function_element = heat->FindNextElement("function");
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,
238 ballonet_element = el->FindNextElement("ballonet");
244 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
246 FGGasCell::~FGGasCell()
250 for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
251 HeatTransferCoeff.clear();
253 for (i = 0; i < Ballonet.size(); i++) delete Ballonet[i];
259 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
261 void FGGasCell::Calculate(double dt)
263 const double AirTemperature = Atmosphere->GetTemperature(); // [Rankine]
264 const double AirPressure = Atmosphere->GetPressure(); // [lbs/ft²]
265 const double AirDensity = Atmosphere->GetDensity(); // [slug/ft³]
266 const double g = Inertial->gravity(); // [lbs/slug]
268 const double OldTemperature = Temperature;
269 const double OldPressure = Pressure;
271 const unsigned int no_ballonets = Ballonet.size();
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();
282 //-- Gas temperature --
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)
289 for (i = 0; i < HeatTransferCoeff.size(); i++) {
290 dU += HeatTransferCoeff[i]->GetValue();
292 // Don't include dt when accounting for adiabatic expansion/contraction.
293 // The rate of adiabatic cooling looks about right: ~5.4 Rankine/1000ft.
296 (dU * dt - Pressure * dVolumeIdeal - BallonetsHeatFlow) /
297 (Cv_gas() * Contents * R);
299 Temperature = AirTemperature;
302 // No simulation of complex temperature changes.
303 // Note: Making the gas cell behave adiabatically might be a better
305 Temperature = AirTemperature;
309 const double IdealPressure =
310 Contents * R * Temperature / (MaxVolume - BallonetsVolume);
311 if (IdealPressure > AirPressure + MaxOverpressure) {
312 Pressure = AirPressure + MaxOverpressure;
314 Pressure = max(IdealPressure, AirPressure);
317 //-- Manual valving --
319 // FIXME: Presently the effect of manual valving is computed using
320 // an ad hoc formula which might not be a good representation
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;
335 max(0.0, Contents - Pressure * VolumeValved / (R * Temperature));
338 //-- Update ballonets. --
339 // Doing that here should give them the opportunity to react to the
341 BallonetsVolume = 0.0;
342 for (i = 0; i < no_ballonets; i++) {
343 Ballonet[i]->Calculate(dt);
344 BallonetsVolume += Ballonet[i]->GetVolume();
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.
353 (AirPressure + MaxOverpressure) *
354 (MaxVolume - BallonetsVolume) / (R * Temperature);
358 Volume = Contents * R * Temperature / Pressure + BallonetsVolume;
360 Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);
362 //-- Current buoyancy --
363 // The buoyancy is computed using the atmospheres local density.
364 Buoyancy = Volume * AirDensity * g;
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);
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
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)) {
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;
389 (1.0 / 4.0) * mass * Yradius * Zradius +
390 (1.0 / 12.0) * mass * Xwidth * Xwidth;
392 (1.0 / 4.0) * mass * Yradius * Zradius +
393 (1.0 / 12.0) * mass * Xwidth * Xwidth;
395 // Not supported. Revert to pointmass model.
396 Ixx = Iyy = Izz = 0.0;
398 // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
403 gasCellM.InitMatrix();
405 GetXYZ(eX) * Mass*slugtolb;
407 GetXYZ(eY) * Mass*slugtolb;
409 GetXYZ(eZ) * Mass*slugtolb;
411 if (no_ballonets > 0) {
412 // Add the mass, moment and inertia of any ballonets.
413 const FGColumnVector3 p = MassBalance->StructuralToBody( GetXYZ() );
415 for (i = 0; i < no_ballonets; i++) {
416 Mass += Ballonet[i]->GetMass();
418 // Add ballonet moments.
420 Ballonet[i]->GetXYZ(eX) * Ballonet[i]->GetMass()*slugtolb;
422 Ballonet[i]->GetXYZ(eY) * Ballonet[i]->GetMass()*slugtolb;
424 Ballonet[i]->GetXYZ(eZ) * Ballonet[i]->GetMass()*slugtolb;
426 // Moments of inertia must be converted to the gas cell frame here.
428 MassBalance->StructuralToBody( Ballonet[i]->GetXYZ() ) - p;
429 // Body basis is in FT.
430 const double mass = Ballonet[i]->GetMass();
431 gasCellJ += Ballonet[i]->GetInertia() +
432 FGMatrix33( 0, - mass*v(1)*v(2), - mass*v(1)*v(3),
433 - mass*v(2)*v(1), 0, - mass*v(2)*v(3),
434 - mass*v(3)*v(1), - mass*v(3)*v(2), 0 );
439 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
440 // The bitmasked value choices are as follows:
441 // unset: In this case (the default) JSBSim would only print
442 // out the normally expected messages, essentially echoing
443 // the config files as they are read. If the environment
444 // variable is not set, debug_lvl is set to 1 internally
445 // 0: This requests JSBSim not to output any messages
447 // 1: This value explicity requests the normal JSBSim
449 // 2: This value asks for a message to be printed out when
450 // a class is instantiated
451 // 4: When this value is set, a message is displayed when a
452 // FGModel object executes its Run() method
453 // 8: When this value is set, various runtime state variables
454 // are printed out periodically
455 // 16: When set various parameters are sanity checked and
456 // a message is printed out when they go out of bounds
458 void FGGasCell::Debug(int from)
460 if (debug_lvl <= 0) return;
462 if (debug_lvl & 1) { // Standard console startup message output
463 if (from == 0) { // Constructor
464 cout << " Gas cell holds " << Contents << " mol " <<
466 cout << " Cell location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
467 vXYZ(eY) << ", " << vXYZ(eZ) << endl;
468 cout << " Maximum volume: " << MaxVolume << " ft3" << endl;
469 cout << " Relief valve release pressure: " << MaxOverpressure <<
471 cout << " Manual valve coefficient: " << ValveCoefficient <<
472 " ft4*sec/slug" << endl;
473 cout << " Initial temperature: " << Temperature << " Rankine" <<
475 cout << " Initial pressure: " << Pressure << " lbs/ft2" << endl;
476 cout << " Initial volume: " << Volume << " ft3" << endl;
477 cout << " Initial mass: " << GetMass() << " slug mass" << endl;
478 cout << " Initial weight: " << GetMass()*lbtoslug << " lbs force" <<
480 cout << " Heat transfer: " << endl;
483 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
484 if (from == 0) cout << "Instantiated: FGGasCell" << endl;
485 if (from == 1) cout << "Destroyed: FGGasCell" << endl;
487 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
489 if (debug_lvl & 8 ) { // Runtime state variables
490 cout << " " << type << " cell holds " << Contents << " mol " << endl;
491 cout << " Temperature: " << Temperature << " Rankine" << endl;
492 cout << " Pressure: " << Pressure << " lbs/ft2" << endl;
493 cout << " Volume: " << Volume << " ft3" << endl;
494 cout << " Mass: " << GetMass() << " slug mass" << endl;
495 cout << " Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
497 if (debug_lvl & 16) { // Sanity checking
499 if (debug_lvl & 64) {
500 if (from == 0) { // Constructor
501 cout << IdSrc << endl;
502 cout << IdHdr << endl;
507 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
509 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
510 const double FGBallonet::R = 3.4071; // [lbs ft/(mol Rankine)]
511 const double FGBallonet::M_air = 0.0019186; // [slug/mol]
512 const double FGBallonet::Cv_air = 5.0/2.0; // [??]
514 FGBallonet::FGBallonet(FGFDMExec* exec, Element* el, int num, FGGasCell* parent)
519 Auxiliary = exec->GetAuxiliary();
520 Atmosphere = exec->GetAtmosphere();
521 PropertyManager = exec->GetPropertyManager();
522 Inertial = exec->GetInertial();
524 ballonetJ = FGMatrix33();
526 MaxVolume = MaxOverpressure = Temperature = Pressure =
527 Contents = Volume = dVolumeIdeal = dU = 0.0;
528 Xradius = Yradius = Zradius = Xwidth = Ywidth = Zwidth = 0.0;
529 ValveCoefficient = ValveOpen = 0.0;
534 // NOTE: In the local system X points north, Y points east and Z points down.
535 element = el->FindElement("location");
537 vXYZ = element->FindElementTripletConvertTo("IN");
539 cerr << "Fatal Error: No location found for this ballonet." << endl;
542 if ((el->FindElement("x_radius") || el->FindElement("x_width")) &&
543 (el->FindElement("y_radius") || el->FindElement("y_width")) &&
544 (el->FindElement("z_radius") || el->FindElement("z_width"))) {
546 if (el->FindElement("x_radius")) {
547 Xradius = el->FindElementValueAsNumberConvertTo("x_radius", "FT");
549 if (el->FindElement("y_radius")) {
550 Yradius = el->FindElementValueAsNumberConvertTo("y_radius", "FT");
552 if (el->FindElement("z_radius")) {
553 Zradius = el->FindElementValueAsNumberConvertTo("z_radius", "FT");
556 if (el->FindElement("x_width")) {
557 Xwidth = el->FindElementValueAsNumberConvertTo("x_width", "FT");
559 if (el->FindElement("y_width")) {
560 Ywidth = el->FindElementValueAsNumberConvertTo("y_width", "FT");
562 if (el->FindElement("z_width")) {
563 Zwidth = el->FindElementValueAsNumberConvertTo("z_width", "FT");
566 // The volume is a (potentially) extruded ellipsoid.
567 // FIXME: However, currently only a few combinations of radius and
568 // width are fully supported.
569 if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
570 (Xwidth == 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
572 MaxVolume = 4.0 * M_PI * Xradius * Yradius * Zradius / 3.0;
573 } else if ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
574 (Xwidth != 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
575 // Cylindrical volume.
576 MaxVolume = M_PI * Yradius * Zradius * Xwidth;
578 cerr << "Warning: Unsupported ballonet shape." << endl;
580 (4.0 * M_PI * Xradius * Yradius * Zradius / 3.0 +
581 M_PI * Yradius * Zradius * Xwidth +
582 M_PI * Xradius * Zradius * Ywidth +
583 M_PI * Xradius * Yradius * Zwidth +
584 2.0 * Xradius * Ywidth * Zwidth +
585 2.0 * Yradius * Xwidth * Zwidth +
586 2.0 * Zradius * Xwidth * Ywidth +
587 Xwidth * Ywidth * Zwidth);
590 cerr << "Fatal Error: Ballonet shape must be given." << endl;
593 if (el->FindElement("max_overpressure")) {
594 MaxOverpressure = el->FindElementValueAsNumberConvertTo("max_overpressure",
597 if (el->FindElement("fullness")) {
598 const double Fullness = el->FindElementValueAsNumber("fullness");
600 Volume = Fullness * MaxVolume;
602 cerr << "Warning: Invalid initial ballonet fullness value." << endl;
605 if (el->FindElement("valve_coefficient")) {
607 el->FindElementValueAsNumberConvertTo("valve_coefficient",
609 ValveCoefficient = max(ValveCoefficient, 0.0);
613 if (Temperature == 0.0) {
614 Temperature = Parent->GetTemperature();
616 if (Pressure == 0.0) {
617 Pressure = Parent->GetPressure();
620 // Calculate initial air content.
621 Contents = Pressure * Volume / (R * Temperature);
623 // Clip to max allowed value.
624 const double IdealPressure = Contents * R * Temperature / MaxVolume;
625 if (IdealPressure > Pressure + MaxOverpressure) {
626 Contents = (Pressure + MaxOverpressure) * MaxVolume / (R * Temperature);
627 Pressure = Pressure + MaxOverpressure;
629 Pressure = max(IdealPressure, Pressure);
632 // Calculate initial air content.
633 Contents = Pressure * MaxVolume / (R * Temperature);
636 Volume = Contents * R * Temperature / Pressure;
638 // Bind relevant properties
639 string property_name, base_property_name;
640 base_property_name = CreateIndexedPropertyName("buoyant_forces/gas-cell", Parent->GetIndex());
641 base_property_name = CreateIndexedPropertyName(base_property_name + "/ballonet", CellNum);
643 property_name = base_property_name + "/max_volume-ft3";
644 PropertyManager->Tie( property_name, &MaxVolume );
645 PropertyManager->SetWritable( property_name, false );
647 property_name = base_property_name + "/temp-R";
648 PropertyManager->Tie( property_name, &Temperature );
650 property_name = base_property_name + "/pressure-psf";
651 PropertyManager->Tie( property_name, &Pressure );
653 property_name = base_property_name + "/volume-ft3";
654 PropertyManager->Tie( property_name, &Volume );
656 property_name = base_property_name + "/contents-mol";
657 PropertyManager->Tie( property_name, &Contents );
659 property_name = base_property_name + "/valve_open";
660 PropertyManager->Tie( property_name, &ValveOpen );
664 // Read heat transfer coefficients
665 if (Element* heat = el->FindElement("heat")) {
666 Element* function_element = heat->FindElement("function");
667 while (function_element) {
668 HeatTransferCoeff.push_back(new FGFunction(PropertyManager,
670 function_element = heat->FindNextElement("function");
673 // Read blower input function
674 if (Element* blower = el->FindElement("blower_input")) {
675 Element* function_element = blower->FindElement("function");
676 BlowerInput = new FGFunction(PropertyManager,
681 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
683 FGBallonet::~FGBallonet()
687 for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
688 HeatTransferCoeff.clear();
696 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
698 void FGBallonet::Calculate(double dt)
700 const double ParentPressure = Parent->GetPressure(); // [lbs/ft²]
701 const double AirPressure = Atmosphere->GetPressure(); // [lbs/ft²]
703 const double OldTemperature = Temperature;
704 const double OldPressure = Pressure;
707 //-- Gas temperature --
709 // The model is based on the ideal gas law.
710 // However, it does look a bit fishy. Please verify.
711 // dT/dt = dU / (Cv n R)
713 for (i = 0; i < HeatTransferCoeff.size(); i++) {
714 dU += HeatTransferCoeff[i]->GetValue();
716 // dt is already accounted for in dVolumeIdeal.
719 (dU * dt - Pressure * dVolumeIdeal) / (Cv_air * Contents * R);
721 Temperature = Parent->GetTemperature();
725 const double IdealPressure = Contents * R * Temperature / MaxVolume;
726 // The pressure is at least that of the parent gas cell.
727 Pressure = max(IdealPressure, ParentPressure);
731 const double AddedVolume = BlowerInput->GetValue() * dt;
732 if (AddedVolume > 0.0) {
733 Contents += Pressure * AddedVolume / (R * Temperature);
737 //-- Pressure relief and manual valving --
738 // FIXME: Presently the effect of valving is computed using
739 // an ad hoc formula which might not be a good representation
741 if ((ValveCoefficient > 0.0) &&
742 ((ValveOpen > 0.0) || (Pressure > AirPressure + MaxOverpressure))) {
743 const double DeltaPressure = Pressure - AirPressure;
744 const double VolumeValved =
745 ((Pressure > AirPressure + MaxOverpressure) ? 1.0 : ValveOpen) *
746 ValveCoefficient * DeltaPressure * dt;
747 // FIXME: Too small values of Contents sometimes leads to NaN.
748 // Currently the minimum is restricted to a safe value.
750 max(1.0, Contents - Pressure * VolumeValved / (R * Temperature));
754 Volume = Contents * R * Temperature / Pressure;
756 Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);
758 // Compute the inertia of the ballonet.
759 // Consider the ballonet as a shape of uniform density.
760 // FIXME: If the ballonet isn't ellipsoid or cylindrical the inertia will
762 ballonetJ = FGMatrix33();
763 const double mass = Contents * M_air;
764 double Ixx, Iyy, Izz;
765 if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
766 (Xwidth == 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
768 Ixx = (1.0 / 5.0) * mass * (Yradius*Yradius + Zradius*Zradius);
769 Iyy = (1.0 / 5.0) * mass * (Xradius*Xradius + Zradius*Zradius);
770 Izz = (1.0 / 5.0) * mass * (Xradius*Xradius + Yradius*Yradius);
771 } else if ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
772 (Xwidth != 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
773 // Cylindrical volume (might not be valid with an elliptical cross-section).
774 Ixx = (1.0 / 2.0) * mass * Yradius * Zradius;
776 (1.0 / 4.0) * mass * Yradius * Zradius +
777 (1.0 / 12.0) * mass * Xwidth * Xwidth;
779 (1.0 / 4.0) * mass * Yradius * Zradius +
780 (1.0 / 12.0) * mass * Xwidth * Xwidth;
782 // Not supported. Revert to pointmass model.
783 Ixx = Iyy = Izz = 0.0;
785 // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
786 ballonetJ(1,1) = Ixx;
787 ballonetJ(2,2) = Iyy;
788 ballonetJ(3,3) = Izz;
791 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
792 // The bitmasked value choices are as follows:
793 // unset: In this case (the default) JSBSim would only print
794 // out the normally expected messages, essentially echoing
795 // the config files as they are read. If the environment
796 // variable is not set, debug_lvl is set to 1 internally
797 // 0: This requests JSBSim not to output any messages
799 // 1: This value explicity requests the normal JSBSim
801 // 2: This value asks for a message to be printed out when
802 // a class is instantiated
803 // 4: When this value is set, a message is displayed when a
804 // FGModel object executes its Run() method
805 // 8: When this value is set, various runtime state variables
806 // are printed out periodically
807 // 16: When set various parameters are sanity checked and
808 // a message is printed out when they go out of bounds
810 void FGBallonet::Debug(int from)
812 if (debug_lvl <= 0) return;
814 if (debug_lvl & 1) { // Standard console startup message output
815 if (from == 0) { // Constructor
816 cout << " Ballonet holds " << Contents << " mol air" << endl;
817 cout << " Location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
818 vXYZ(eY) << ", " << vXYZ(eZ) << endl;
819 cout << " Maximum volume: " << MaxVolume << " ft3" << endl;
820 cout << " Relief valve release pressure: " << MaxOverpressure <<
822 cout << " Relief valve coefficient: " << ValveCoefficient <<
823 " ft4*sec/slug" << endl;
824 cout << " Initial temperature: " << Temperature << " Rankine" <<
826 cout << " Initial pressure: " << Pressure << " lbs/ft2" << endl;
827 cout << " Initial volume: " << Volume << " ft3" << endl;
828 cout << " Initial mass: " << GetMass() << " slug mass" << endl;
829 cout << " Initial weight: " << GetMass()*lbtoslug <<
830 " lbs force" << endl;
831 cout << " Heat transfer: " << endl;
834 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
835 if (from == 0) cout << "Instantiated: FGBallonet" << endl;
836 if (from == 1) cout << "Destroyed: FGBallonet" << endl;
838 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
840 if (debug_lvl & 8 ) { // Runtime state variables
841 cout << " Ballonet holds " << Contents <<
843 cout << " Temperature: " << Temperature << " Rankine" << endl;
844 cout << " Pressure: " << Pressure << " lbs/ft2" << endl;
845 cout << " Volume: " << Volume << " ft3" << endl;
846 cout << " Mass: " << GetMass() << " slug mass" << endl;
847 cout << " Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
849 if (debug_lvl & 16) { // Sanity checking
851 if (debug_lvl & 64) {
852 if (from == 0) { // Constructor
853 cout << IdSrc << endl;
854 cout << IdHdr << endl;