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"
44 #include "input_output/FGXMLElement.h"
56 static const char *IdSrc = "$Id: FGGasCell.cpp,v 1.13 2010/12/29 22:39:25 andgi Exp $";
57 static const char *IdHdr = ID_GASCELL;
59 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
63 const double FGGasCell::R = 3.4071; // [lbs ft/(mol Rankine)]
64 const double FGGasCell::M_air = 0.0019186; // [slug/mol]
65 const double FGGasCell::M_hydrogen = 0.00013841; // [slug/mol]
66 const double FGGasCell::M_helium = 0.00027409; // [slug/mol]
68 FGGasCell::FGGasCell(FGFDMExec* exec, Element* el, int num) : FGForce(exec)
73 Auxiliary = exec->GetAuxiliary();
74 Atmosphere = exec->GetAtmosphere();
75 PropertyManager = exec->GetPropertyManager();
76 Inertial = exec->GetInertial();
77 MassBalance = exec->GetMassBalance();
79 gasCellJ = FGMatrix33();
80 gasCellM = FGColumnVector3();
82 Buoyancy = MaxVolume = MaxOverpressure = Temperature = Pressure =
83 Contents = Volume = dVolumeIdeal = 0.0;
84 Xradius = Yradius = Zradius = Xwidth = Ywidth = Zwidth = 0.0;
85 ValveCoefficient = ValveOpen = 0.0;
88 // NOTE: In the local system X points north, Y points east and Z points down.
89 SetTransformType(FGForce::tLocalBody);
91 type = el->GetAttributeValue("type");
92 if (type == "HYDROGEN") Type = ttHYDROGEN;
93 else if (type == "HELIUM") Type = ttHELIUM;
94 else if (type == "AIR") Type = ttAIR;
95 else Type = ttUNKNOWN;
97 element = el->FindElement("location");
99 vXYZ = element->FindElementTripletConvertTo("IN");
101 cerr << "Fatal Error: No location found for this gas cell." << endl;
104 if ((el->FindElement("x_radius") || el->FindElement("x_width")) &&
105 (el->FindElement("y_radius") || el->FindElement("y_width")) &&
106 (el->FindElement("z_radius") || el->FindElement("z_width"))) {
108 if (el->FindElement("x_radius")) {
109 Xradius = el->FindElementValueAsNumberConvertTo("x_radius", "FT");
111 if (el->FindElement("y_radius")) {
112 Yradius = el->FindElementValueAsNumberConvertTo("y_radius", "FT");
114 if (el->FindElement("z_radius")) {
115 Zradius = el->FindElementValueAsNumberConvertTo("z_radius", "FT");
118 if (el->FindElement("x_width")) {
119 Xwidth = el->FindElementValueAsNumberConvertTo("x_width", "FT");
121 if (el->FindElement("y_width")) {
122 Ywidth = el->FindElementValueAsNumberConvertTo("y_width", "FT");
124 if (el->FindElement("z_width")) {
125 Zwidth = el->FindElementValueAsNumberConvertTo("z_width", "FT");
128 // The volume is a (potentially) extruded ellipsoid.
129 // However, currently only a few combinations of radius and width are
131 if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
132 (Xwidth == 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
134 MaxVolume = 4.0 * M_PI * Xradius * Yradius * Zradius / 3.0;
135 } else if ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
136 (Xwidth != 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
137 // Cylindrical volume.
138 MaxVolume = M_PI * Yradius * Zradius * Xwidth;
140 cerr << "Warning: Unsupported gas cell shape." << endl;
142 (4.0 * M_PI * Xradius * Yradius * Zradius / 3.0 +
143 M_PI * Yradius * Zradius * Xwidth +
144 M_PI * Xradius * Zradius * Ywidth +
145 M_PI * Xradius * Yradius * Zwidth +
146 2.0 * Xradius * Ywidth * Zwidth +
147 2.0 * Yradius * Xwidth * Zwidth +
148 2.0 * Zradius * Xwidth * Ywidth +
149 Xwidth * Ywidth * Zwidth);
152 cerr << "Fatal Error: Gas cell shape must be given." << endl;
155 if (el->FindElement("max_overpressure")) {
156 MaxOverpressure = el->FindElementValueAsNumberConvertTo("max_overpressure",
159 if (el->FindElement("fullness")) {
160 const double Fullness = el->FindElementValueAsNumber("fullness");
162 Volume = Fullness * MaxVolume;
164 cerr << "Warning: Invalid initial gas cell fullness value." << endl;
167 if (el->FindElement("valve_coefficient")) {
169 el->FindElementValueAsNumberConvertTo("valve_coefficient",
171 ValveCoefficient = max(ValveCoefficient, 0.0);
177 if (Temperature == 0.0) {
178 Temperature = Atmosphere->GetTemperature();
180 if (Pressure == 0.0) {
181 Pressure = Atmosphere->GetPressure();
184 // Calculate initial gas content.
185 Contents = Pressure * Volume / (R * Temperature);
187 // Clip to max allowed value.
188 const double IdealPressure = Contents * R * Temperature / MaxVolume;
189 if (IdealPressure > Pressure + MaxOverpressure) {
190 Contents = (Pressure + MaxOverpressure) * MaxVolume / (R * Temperature);
191 Pressure = Pressure + MaxOverpressure;
193 Pressure = max(IdealPressure, Pressure);
196 // Calculate initial gas content.
197 Contents = Pressure * MaxVolume / (R * Temperature);
200 Volume = Contents * R * Temperature / Pressure;
201 Mass = Contents * M_gas();
203 // Bind relevant properties
204 string property_name, base_property_name;
206 base_property_name = CreateIndexedPropertyName("buoyant_forces/gas-cell", CellNum);
208 property_name = base_property_name + "/max_volume-ft3";
209 PropertyManager->Tie( property_name.c_str(), &MaxVolume, false );
210 PropertyManager->SetWritable( property_name, false );
211 property_name = base_property_name + "/temp-R";
212 PropertyManager->Tie( property_name.c_str(), &Temperature, false );
213 property_name = base_property_name + "/pressure-psf";
214 PropertyManager->Tie( property_name.c_str(), &Pressure, false );
215 property_name = base_property_name + "/volume-ft3";
216 PropertyManager->Tie( property_name.c_str(), &Volume, false );
217 property_name = base_property_name + "/buoyancy-lbs";
218 PropertyManager->Tie( property_name.c_str(), &Buoyancy, false );
219 property_name = base_property_name + "/contents-mol";
220 PropertyManager->Tie( property_name.c_str(), &Contents, false );
221 property_name = base_property_name + "/valve_open";
222 PropertyManager->Tie( property_name.c_str(), &ValveOpen, false );
226 // Read heat transfer coefficients
227 if (Element* heat = el->FindElement("heat")) {
228 Element* function_element = heat->FindElement("function");
229 while (function_element) {
230 HeatTransferCoeff.push_back(new FGFunction(PropertyManager,
232 function_element = heat->FindNextElement("function");
236 // Load ballonets if there are any
237 if (Element* ballonet_element = el->FindElement("ballonet")) {
238 while (ballonet_element) {
239 Ballonet.push_back(new FGBallonet(exec,
243 ballonet_element = el->FindNextElement("ballonet");
249 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
251 FGGasCell::~FGGasCell()
255 for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
256 HeatTransferCoeff.clear();
258 for (i = 0; i < Ballonet.size(); i++) delete Ballonet[i];
264 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
266 void FGGasCell::Calculate(double dt)
268 const double AirTemperature = Atmosphere->GetTemperature(); // [Rankine]
269 const double AirPressure = Atmosphere->GetPressure(); // [lbs/ft²]
270 const double AirDensity = Atmosphere->GetDensity(); // [slug/ft³]
271 const double g = Inertial->gravity(); // [lbs/slug]
273 const double OldTemperature = Temperature;
274 const double OldPressure = Pressure;
276 const unsigned int no_ballonets = Ballonet.size();
278 //-- Read ballonet state --
279 // NOTE: This model might need a more proper integration technique.
280 double BallonetsVolume = 0.0;
281 double BallonetsHeatFlow = 0.0;
282 for (i = 0; i < no_ballonets; i++) {
283 BallonetsVolume += Ballonet[i]->GetVolume();
284 BallonetsHeatFlow += Ballonet[i]->GetHeatFlow();
287 //-- Gas temperature --
289 if (HeatTransferCoeff.size() > 0) {
290 // The model is based on the ideal gas law.
291 // However, it does look a bit fishy. Please verify.
292 // dT/dt = dU / (Cv n R)
294 for (i = 0; i < HeatTransferCoeff.size(); i++) {
295 dU += HeatTransferCoeff[i]->GetValue();
297 // Don't include dt when accounting for adiabatic expansion/contraction.
298 // The rate of adiabatic cooling looks about right: ~5.4 Rankine/1000ft.
301 (dU * dt - Pressure * dVolumeIdeal - BallonetsHeatFlow) /
302 (Cv_gas() * Contents * R);
304 Temperature = AirTemperature;
307 // No simulation of complex temperature changes.
308 // Note: Making the gas cell behave adiabatically might be a better
310 Temperature = AirTemperature;
314 const double IdealPressure =
315 Contents * R * Temperature / (MaxVolume - BallonetsVolume);
316 if (IdealPressure > AirPressure + MaxOverpressure) {
317 Pressure = AirPressure + MaxOverpressure;
319 Pressure = max(IdealPressure, AirPressure);
322 //-- Manual valving --
324 // FIXME: Presently the effect of manual valving is computed using
325 // an ad hoc formula which might not be a good representation
327 if ((ValveCoefficient > 0.0) && (ValveOpen > 0.0)) {
328 // First compute the difference in pressure between the gas in the
329 // cell and the air above it.
330 // FixMe: CellHeight should depend on current volume.
331 const double CellHeight = 2 * Zradius + Zwidth; // [ft]
332 const double GasMass = Contents * M_gas(); // [slug]
333 const double GasVolume = Contents * R * Temperature / Pressure; // [ft³]
334 const double GasDensity = GasMass / GasVolume;
335 const double DeltaPressure =
336 Pressure + CellHeight * g * (AirDensity - GasDensity) - AirPressure;
337 const double VolumeValved =
338 ValveOpen * ValveCoefficient * DeltaPressure * dt;
340 max(0.0, Contents - Pressure * VolumeValved / (R * Temperature));
343 //-- Update ballonets. --
344 // Doing that here should give them the opportunity to react to the
346 BallonetsVolume = 0.0;
347 for (i = 0; i < no_ballonets; i++) {
348 Ballonet[i]->Calculate(dt);
349 BallonetsVolume += Ballonet[i]->GetVolume();
352 //-- Automatic safety valving. --
353 if (Contents * R * Temperature / (MaxVolume - BallonetsVolume) >
354 AirPressure + MaxOverpressure) {
355 // Gas is automatically valved. Valving capacity is assumed to be infinite.
356 // FIXME: This could/should be replaced by damage to the gas cell envelope.
358 (AirPressure + MaxOverpressure) *
359 (MaxVolume - BallonetsVolume) / (R * Temperature);
363 Volume = Contents * R * Temperature / Pressure + BallonetsVolume;
365 Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);
367 //-- Current buoyancy --
368 // The buoyancy is computed using the atmospheres local density.
369 Buoyancy = Volume * AirDensity * g;
371 // Note: This is gross buoyancy. The weight of the gas itself and
372 // any ballonets is not deducted here as the effects of the gas mass
373 // is handled by FGMassBalance.
374 vFn.InitMatrix(0.0, 0.0, - Buoyancy);
376 // Compute the inertia of the gas cell.
377 // Consider the gas cell as a shape of uniform density.
378 // FIXME: If the cell isn't ellipsoid or cylindrical the inertia will
380 gasCellJ = FGMatrix33();
381 const double mass = Contents * M_gas();
382 double Ixx, Iyy, Izz;
383 if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
384 (Xwidth == 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
386 Ixx = (1.0 / 5.0) * mass * (Yradius*Yradius + Zradius*Zradius);
387 Iyy = (1.0 / 5.0) * mass * (Xradius*Xradius + Zradius*Zradius);
388 Izz = (1.0 / 5.0) * mass * (Xradius*Xradius + Yradius*Yradius);
389 } else if ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
390 (Xwidth != 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
391 // Cylindrical volume (might not be valid with an elliptical cross-section).
392 Ixx = (1.0 / 2.0) * mass * Yradius * Zradius;
394 (1.0 / 4.0) * mass * Yradius * Zradius +
395 (1.0 / 12.0) * mass * Xwidth * Xwidth;
397 (1.0 / 4.0) * mass * Yradius * Zradius +
398 (1.0 / 12.0) * mass * Xwidth * Xwidth;
400 // Not supported. Revert to pointmass model.
401 Ixx = Iyy = Izz = 0.0;
403 // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
408 gasCellM.InitMatrix();
410 GetXYZ(eX) * Mass*slugtolb;
412 GetXYZ(eY) * Mass*slugtolb;
414 GetXYZ(eZ) * Mass*slugtolb;
416 if (no_ballonets > 0) {
417 // Add the mass, moment and inertia of any ballonets.
418 const FGColumnVector3 p = MassBalance->StructuralToBody( GetXYZ() );
420 for (i = 0; i < no_ballonets; i++) {
421 Mass += Ballonet[i]->GetMass();
423 // Add ballonet moments.
425 Ballonet[i]->GetXYZ(eX) * Ballonet[i]->GetMass()*slugtolb;
427 Ballonet[i]->GetXYZ(eY) * Ballonet[i]->GetMass()*slugtolb;
429 Ballonet[i]->GetXYZ(eZ) * Ballonet[i]->GetMass()*slugtolb;
431 // Moments of inertia must be converted to the gas cell frame here.
433 MassBalance->StructuralToBody( Ballonet[i]->GetXYZ() ) - p;
434 // Body basis is in FT.
435 const double mass = Ballonet[i]->GetMass();
436 gasCellJ += Ballonet[i]->GetInertia() +
437 FGMatrix33( 0, - mass*v(1)*v(2), - mass*v(1)*v(3),
438 - mass*v(2)*v(1), 0, - mass*v(2)*v(3),
439 - mass*v(3)*v(1), - mass*v(3)*v(2), 0 );
444 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
445 // The bitmasked value choices are as follows:
446 // unset: In this case (the default) JSBSim would only print
447 // out the normally expected messages, essentially echoing
448 // the config files as they are read. If the environment
449 // variable is not set, debug_lvl is set to 1 internally
450 // 0: This requests JSBSim not to output any messages
452 // 1: This value explicity requests the normal JSBSim
454 // 2: This value asks for a message to be printed out when
455 // a class is instantiated
456 // 4: When this value is set, a message is displayed when a
457 // FGModel object executes its Run() method
458 // 8: When this value is set, various runtime state variables
459 // are printed out periodically
460 // 16: When set various parameters are sanity checked and
461 // a message is printed out when they go out of bounds
463 void FGGasCell::Debug(int from)
465 if (debug_lvl <= 0) return;
467 if (debug_lvl & 1) { // Standard console startup message output
468 if (from == 0) { // Constructor
469 cout << " Gas cell holds " << Contents << " mol " <<
471 cout << " Cell location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
472 vXYZ(eY) << ", " << vXYZ(eZ) << endl;
473 cout << " Maximum volume: " << MaxVolume << " ft3" << endl;
474 cout << " Relief valve release pressure: " << MaxOverpressure <<
476 cout << " Manual valve coefficient: " << ValveCoefficient <<
477 " ft4*sec/slug" << endl;
478 cout << " Initial temperature: " << Temperature << " Rankine" <<
480 cout << " Initial pressure: " << Pressure << " lbs/ft2" << endl;
481 cout << " Initial volume: " << Volume << " ft3" << endl;
482 cout << " Initial mass: " << GetMass() << " slug mass" << endl;
483 cout << " Initial weight: " << GetMass()*lbtoslug << " lbs force" <<
485 cout << " Heat transfer: " << endl;
488 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
489 if (from == 0) cout << "Instantiated: FGGasCell" << endl;
490 if (from == 1) cout << "Destroyed: FGGasCell" << endl;
492 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
494 if (debug_lvl & 8 ) { // Runtime state variables
495 cout << " " << type << " cell holds " << Contents << " mol " << endl;
496 cout << " Temperature: " << Temperature << " Rankine" << endl;
497 cout << " Pressure: " << Pressure << " lbs/ft2" << endl;
498 cout << " Volume: " << Volume << " ft3" << endl;
499 cout << " Mass: " << GetMass() << " slug mass" << endl;
500 cout << " Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
502 if (debug_lvl & 16) { // Sanity checking
504 if (debug_lvl & 64) {
505 if (from == 0) { // Constructor
506 cout << IdSrc << endl;
507 cout << IdHdr << endl;
512 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
514 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
515 const double FGBallonet::R = 3.4071; // [lbs ft/(mol Rankine)]
516 const double FGBallonet::M_air = 0.0019186; // [slug/mol]
517 const double FGBallonet::Cv_air = 5.0/2.0; // [??]
519 FGBallonet::FGBallonet(FGFDMExec* exec, Element* el, int num, FGGasCell* parent)
524 Auxiliary = exec->GetAuxiliary();
525 Atmosphere = exec->GetAtmosphere();
526 PropertyManager = exec->GetPropertyManager();
527 Inertial = exec->GetInertial();
529 ballonetJ = FGMatrix33();
531 MaxVolume = MaxOverpressure = Temperature = Pressure =
532 Contents = Volume = dVolumeIdeal = dU = 0.0;
533 Xradius = Yradius = Zradius = Xwidth = Ywidth = Zwidth = 0.0;
534 ValveCoefficient = ValveOpen = 0.0;
539 // NOTE: In the local system X points north, Y points east and Z points down.
540 element = el->FindElement("location");
542 vXYZ = element->FindElementTripletConvertTo("IN");
544 cerr << "Fatal Error: No location found for this ballonet." << endl;
547 if ((el->FindElement("x_radius") || el->FindElement("x_width")) &&
548 (el->FindElement("y_radius") || el->FindElement("y_width")) &&
549 (el->FindElement("z_radius") || el->FindElement("z_width"))) {
551 if (el->FindElement("x_radius")) {
552 Xradius = el->FindElementValueAsNumberConvertTo("x_radius", "FT");
554 if (el->FindElement("y_radius")) {
555 Yradius = el->FindElementValueAsNumberConvertTo("y_radius", "FT");
557 if (el->FindElement("z_radius")) {
558 Zradius = el->FindElementValueAsNumberConvertTo("z_radius", "FT");
561 if (el->FindElement("x_width")) {
562 Xwidth = el->FindElementValueAsNumberConvertTo("x_width", "FT");
564 if (el->FindElement("y_width")) {
565 Ywidth = el->FindElementValueAsNumberConvertTo("y_width", "FT");
567 if (el->FindElement("z_width")) {
568 Zwidth = el->FindElementValueAsNumberConvertTo("z_width", "FT");
571 // The volume is a (potentially) extruded ellipsoid.
572 // FIXME: However, currently only a few combinations of radius and
573 // width are fully supported.
574 if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
575 (Xwidth == 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
577 MaxVolume = 4.0 * M_PI * Xradius * Yradius * Zradius / 3.0;
578 } else if ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
579 (Xwidth != 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
580 // Cylindrical volume.
581 MaxVolume = M_PI * Yradius * Zradius * Xwidth;
583 cerr << "Warning: Unsupported ballonet shape." << endl;
585 (4.0 * M_PI * Xradius * Yradius * Zradius / 3.0 +
586 M_PI * Yradius * Zradius * Xwidth +
587 M_PI * Xradius * Zradius * Ywidth +
588 M_PI * Xradius * Yradius * Zwidth +
589 2.0 * Xradius * Ywidth * Zwidth +
590 2.0 * Yradius * Xwidth * Zwidth +
591 2.0 * Zradius * Xwidth * Ywidth +
592 Xwidth * Ywidth * Zwidth);
595 cerr << "Fatal Error: Ballonet shape must be given." << endl;
598 if (el->FindElement("max_overpressure")) {
599 MaxOverpressure = el->FindElementValueAsNumberConvertTo("max_overpressure",
602 if (el->FindElement("fullness")) {
603 const double Fullness = el->FindElementValueAsNumber("fullness");
605 Volume = Fullness * MaxVolume;
607 cerr << "Warning: Invalid initial ballonet fullness value." << endl;
610 if (el->FindElement("valve_coefficient")) {
612 el->FindElementValueAsNumberConvertTo("valve_coefficient",
614 ValveCoefficient = max(ValveCoefficient, 0.0);
618 if (Temperature == 0.0) {
619 Temperature = Parent->GetTemperature();
621 if (Pressure == 0.0) {
622 Pressure = Parent->GetPressure();
625 // Calculate initial air content.
626 Contents = Pressure * Volume / (R * Temperature);
628 // Clip to max allowed value.
629 const double IdealPressure = Contents * R * Temperature / MaxVolume;
630 if (IdealPressure > Pressure + MaxOverpressure) {
631 Contents = (Pressure + MaxOverpressure) * MaxVolume / (R * Temperature);
632 Pressure = Pressure + MaxOverpressure;
634 Pressure = max(IdealPressure, Pressure);
637 // Calculate initial air content.
638 Contents = Pressure * MaxVolume / (R * Temperature);
641 Volume = Contents * R * Temperature / Pressure;
643 // Bind relevant properties
644 string property_name, base_property_name;
645 base_property_name = CreateIndexedPropertyName("buoyant_forces/gas-cell", Parent->GetIndex());
646 base_property_name = CreateIndexedPropertyName(base_property_name + "/ballonet", CellNum);
648 property_name = base_property_name + "/max_volume-ft3";
649 PropertyManager->Tie( property_name, &MaxVolume, false );
650 PropertyManager->SetWritable( property_name, false );
652 property_name = base_property_name + "/temp-R";
653 PropertyManager->Tie( property_name, &Temperature, false );
655 property_name = base_property_name + "/pressure-psf";
656 PropertyManager->Tie( property_name, &Pressure, false );
658 property_name = base_property_name + "/volume-ft3";
659 PropertyManager->Tie( property_name, &Volume, false );
661 property_name = base_property_name + "/contents-mol";
662 PropertyManager->Tie( property_name, &Contents, false );
664 property_name = base_property_name + "/valve_open";
665 PropertyManager->Tie( property_name, &ValveOpen, false );
669 // Read heat transfer coefficients
670 if (Element* heat = el->FindElement("heat")) {
671 Element* function_element = heat->FindElement("function");
672 while (function_element) {
673 HeatTransferCoeff.push_back(new FGFunction(PropertyManager,
675 function_element = heat->FindNextElement("function");
678 // Read blower input function
679 if (Element* blower = el->FindElement("blower_input")) {
680 Element* function_element = blower->FindElement("function");
681 BlowerInput = new FGFunction(PropertyManager,
686 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
688 FGBallonet::~FGBallonet()
692 for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
693 HeatTransferCoeff.clear();
701 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
703 void FGBallonet::Calculate(double dt)
705 const double ParentPressure = Parent->GetPressure(); // [lbs/ft²]
706 const double AirPressure = Atmosphere->GetPressure(); // [lbs/ft²]
708 const double OldTemperature = Temperature;
709 const double OldPressure = Pressure;
712 //-- Gas temperature --
714 // The model is based on the ideal gas law.
715 // However, it does look a bit fishy. Please verify.
716 // dT/dt = dU / (Cv n R)
718 for (i = 0; i < HeatTransferCoeff.size(); i++) {
719 dU += HeatTransferCoeff[i]->GetValue();
721 // dt is already accounted for in dVolumeIdeal.
724 (dU * dt - Pressure * dVolumeIdeal) / (Cv_air * Contents * R);
726 Temperature = Parent->GetTemperature();
730 const double IdealPressure = Contents * R * Temperature / MaxVolume;
731 // The pressure is at least that of the parent gas cell.
732 Pressure = max(IdealPressure, ParentPressure);
736 const double AddedVolume = BlowerInput->GetValue() * dt;
737 if (AddedVolume > 0.0) {
738 Contents += Pressure * AddedVolume / (R * Temperature);
742 //-- Pressure relief and manual valving --
743 // FIXME: Presently the effect of valving is computed using
744 // an ad hoc formula which might not be a good representation
746 if ((ValveCoefficient > 0.0) &&
747 ((ValveOpen > 0.0) || (Pressure > AirPressure + MaxOverpressure))) {
748 const double DeltaPressure = Pressure - AirPressure;
749 const double VolumeValved =
750 ((Pressure > AirPressure + MaxOverpressure) ? 1.0 : ValveOpen) *
751 ValveCoefficient * DeltaPressure * dt;
752 // FIXME: Too small values of Contents sometimes leads to NaN.
753 // Currently the minimum is restricted to a safe value.
755 max(1.0, Contents - Pressure * VolumeValved / (R * Temperature));
759 Volume = Contents * R * Temperature / Pressure;
761 Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);
763 // Compute the inertia of the ballonet.
764 // Consider the ballonet as a shape of uniform density.
765 // FIXME: If the ballonet isn't ellipsoid or cylindrical the inertia will
767 ballonetJ = FGMatrix33();
768 const double mass = Contents * M_air;
769 double Ixx, Iyy, Izz;
770 if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
771 (Xwidth == 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
773 Ixx = (1.0 / 5.0) * mass * (Yradius*Yradius + Zradius*Zradius);
774 Iyy = (1.0 / 5.0) * mass * (Xradius*Xradius + Zradius*Zradius);
775 Izz = (1.0 / 5.0) * mass * (Xradius*Xradius + Yradius*Yradius);
776 } else if ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
777 (Xwidth != 0.0) && (Ywidth == 0.0) && (Zwidth == 0.0)) {
778 // Cylindrical volume (might not be valid with an elliptical cross-section).
779 Ixx = (1.0 / 2.0) * mass * Yradius * Zradius;
781 (1.0 / 4.0) * mass * Yradius * Zradius +
782 (1.0 / 12.0) * mass * Xwidth * Xwidth;
784 (1.0 / 4.0) * mass * Yradius * Zradius +
785 (1.0 / 12.0) * mass * Xwidth * Xwidth;
787 // Not supported. Revert to pointmass model.
788 Ixx = Iyy = Izz = 0.0;
790 // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
791 ballonetJ(1,1) = Ixx;
792 ballonetJ(2,2) = Iyy;
793 ballonetJ(3,3) = Izz;
796 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
797 // The bitmasked value choices are as follows:
798 // unset: In this case (the default) JSBSim would only print
799 // out the normally expected messages, essentially echoing
800 // the config files as they are read. If the environment
801 // variable is not set, debug_lvl is set to 1 internally
802 // 0: This requests JSBSim not to output any messages
804 // 1: This value explicity requests the normal JSBSim
806 // 2: This value asks for a message to be printed out when
807 // a class is instantiated
808 // 4: When this value is set, a message is displayed when a
809 // FGModel object executes its Run() method
810 // 8: When this value is set, various runtime state variables
811 // are printed out periodically
812 // 16: When set various parameters are sanity checked and
813 // a message is printed out when they go out of bounds
815 void FGBallonet::Debug(int from)
817 if (debug_lvl <= 0) return;
819 if (debug_lvl & 1) { // Standard console startup message output
820 if (from == 0) { // Constructor
821 cout << " Ballonet holds " << Contents << " mol air" << endl;
822 cout << " Location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
823 vXYZ(eY) << ", " << vXYZ(eZ) << endl;
824 cout << " Maximum volume: " << MaxVolume << " ft3" << endl;
825 cout << " Relief valve release pressure: " << MaxOverpressure <<
827 cout << " Relief valve coefficient: " << ValveCoefficient <<
828 " ft4*sec/slug" << endl;
829 cout << " Initial temperature: " << Temperature << " Rankine" <<
831 cout << " Initial pressure: " << Pressure << " lbs/ft2" << endl;
832 cout << " Initial volume: " << Volume << " ft3" << endl;
833 cout << " Initial mass: " << GetMass() << " slug mass" << endl;
834 cout << " Initial weight: " << GetMass()*lbtoslug <<
835 " lbs force" << endl;
836 cout << " Heat transfer: " << endl;
839 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
840 if (from == 0) cout << "Instantiated: FGBallonet" << endl;
841 if (from == 1) cout << "Destroyed: FGBallonet" << endl;
843 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
845 if (debug_lvl & 8 ) { // Runtime state variables
846 cout << " Ballonet holds " << Contents <<
848 cout << " Temperature: " << Temperature << " Rankine" << endl;
849 cout << " Pressure: " << Pressure << " lbs/ft2" << endl;
850 cout << " Volume: " << Volume << " ft3" << endl;
851 cout << " Mass: " << GetMass() << " slug mass" << endl;
852 cout << " Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
854 if (debug_lvl & 16) { // Sanity checking
856 if (debug_lvl & 64) {
857 if (from == 0) { // Constructor
858 cout << IdSrc << endl;
859 cout << IdHdr << endl;