Weather Balloon script crash

[bcoconni]

Discussed in #643

^{Originally posted by @dpculp June 5, 2022}
When executing the script weather-balloon.xml as supplied by the repo I get a core dump "Floating point exception" after about 4980 seconds of virtual time. This is at the point in the flight around 92K feet where the balloon is deflating and beginning descent.

I changed the termination event to stop execution at 80K feet of altitude and the script ran perfectly.

It seems the transition from ascent to descent is not working presently in the weather balloon (gas cell) model. I'll investigate further.

[bcoconni]

The program crashes at src/models/FGGasCell.cpp line 323 due a division by zero because GasVolume is 0.0.

jsbsim/src/models/FGGasCell.cpp

Lines 317 to 330 in 33735db

	// First compute the difference in pressure between the gas in the
	// cell and the air above it.
	// FixMe: CellHeight should depend on current volume.
	const double CellHeight = 2 * Zradius + Zwidth; // [ft]
	const double GasMass = Contents * M_gas(); // [slug]
	const double GasVolume = Contents * R * Temperature / Pressure; // [ft^3]
	const double GasDensity = GasMass / GasVolume;
	const double DeltaPressure =
	Pressure + CellHeight * g * (AirDensity - GasDensity) - AirPressure;
	const double VolumeValved =
	ValveOpen * ValveCoefficient * DeltaPressure * dt;
	Contents =
	max(0.0, Contents - Pressure * VolumeValved / (R * Temperature));
	}

By setting a minimum value of 1E-8 for Contents the crash is fixed:

    Contents =
-      max(0.0, Contents - Pressure * VolumeValved / (R * Temperature));
+      max(1E-8, Contents - Pressure * VolumeValved / (R * Temperature));

[bcoconni]

Setting the minimum content to 10^-8 is a workaround to avoid the crash. The mere fact that this formula can empty the balloon until its interior volume is zero makes it questionable.

@andgi I know you wrote a disclaimer in the code about the formula you used to compute the valved flow of a balloon, but I am surprised that the weight of the column of air CellHeight * g * (AirDensity - GasDensity) is involved in such a formula ?

In addition, I have the feeling that the sign of AirDensity - GasDensity is incorrect: the air column weight should push the ambient air inside the balloon because the ambient air is heavier:

$$P_{air}=P_0+\rho_{air}gh$$

$$P_{gas}=P_0+\rho_{gas}gh$$

since $\rho_{air}>\rho_{gas}$ then $P_{air}>P_{gas}$ whatever $P_0$ is. As the air moves from high pressures to low pressures, this formula should result in air moving from the outside to the inside of the balloon which is exactly the opposite of what the valve is expected to do.

[bcoconni]

Actually I am now certain the formula is wrong because AirDensity should have the same sign than AirPressure in the code below:

jsbsim/src/models/FGGasCell.cpp

Lines 324 to 325 in 33735db

	const double DeltaPressure =
	Pressure + CellHeight * g * (AirDensity - GasDensity) - AirPressure;

Indeed the computation of $\Delta P$ is subtracting two $P+\rho gh$:

$$P_{air}=P^0_{air}+\rho_{air}gh$$

$$P_{gas}=P^0_{gas}+\rho_{gas}gh$$

which gives:

$$\Delta P=P_{gas}-P_{air}=P^0_{gas}+\rho_{gas}gh-P^0_{air}-\rho_{air}gh$$

$$\Delta P=(P^0_{gas}-P^0_{air})+(\rho_{gas}-\rho_{air})gh$$

So the air quantities $P^0_{air}$ and $\rho_{air}$ have the same sign which they don't in FGGasCell.cpp as shown above. And as I wrote previously the air column weight should be pushing ambient air inside the balloon.

[andgi]

Hi Bertrand, It is some years since a last looked at this code (seems to be in 2008), but I'm not sure the computation of the delta pressure at the manual control release-valve at the top of the gas cell is wrong. The "standard" pressure in the gas cell is computed for the bottom of the cell (and at least as large as the air pressure). The pressure difference at the top is the difference in gas and air pressure at the bottom plus the difference from the height of the gas cell. Cheers, Anders

…

On Sun, 12 Jun 2022 at 13:30, Bertrand Coconnier ***@***.***> wrote: Actually I am now certain the formula is wrong because AirDensity should have the same sign than AirPressure in the code below: https://github.com/JSBSim-Team/jsbsim/blob/33735db0fb1e4fdf47f253a4215dbe3bb63a8e92/src/models/FGGasCell.cpp#L324-L325 Indeed the computation of $\Delta P$ is subtracting two $P+\rho gh$: $$P_{air}=P^0_{air}+\rho_{air}gh$$ $$P_{gas}=P^0_{gas}+\rho_{gas}gh$$ which gives: $$\Delta P=P_{gas}-P_{air}=P^0_{gas}+\rho_{gas}gh-P^0_{air}-\rho_{air}gh$$ $$\Delta P=(P^0_{gas}-P^0_{air})+(\rho_{gas}-\rho_{air})gh$$ So the air quantities $P^0_{air}$ and $\rho_{air}$ have the same sign which they don't in FGGasCell.cpp as shown above. And as I wrote previously the air column weight is pushing ambient air inside the balloon. — Reply to this email directly, view it on GitHub <#654 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACDXMKRKIPLMKJFFEWBKTODVOXC2BANCNFSM5YQRCHVA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[bcoconni]

Hi Anders,

The "standard" pressure in the gas cell is computed for the bottom of the cell (and at least as large as the air pressure). The pressure difference at the top is the difference in gas and air pressure at the bottom plus the difference from the height of the gas cell.

Thanks for the explanation. I think I got the idea: CellHeight is the absolute value of the column height i.e. the code uses the formula $P=P^0-\rho\times g\times(-h)$ where $-h$ is stored in CellHeight which is positive since $h&lt;0$ because the valve is assumed to be at the top.

This makes me realize that I had assumed that the valve was located at a position below the reference pressure $P^0$. This must be a reminiscence of the draining water tank exercise that is taught in fluid mechanics 😄.

It remains that the physics is somehow wrongly simulated in FGGasCell since balloon deflation cannot reach a state where not a single molecule of gas is remaining inside the ballon envelope (which is what Contents == 0.0means).

Since the mass flow is generated by the pressure difference $\Delta P$ between the ambient air and the gas inside the balloon envelope, at some point $\Delta P$ should vanish before Contents reaches zero. So I still have this intuition that there is something wrong with the computation of DeltaPressure.