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