Avoid non physical values for the ambient pressure and temperature (#836

)

* Validate input values to avoid NaNs and other silly behaviors.

* Add tests for invalid input values

* Add a test for the sea level pressure

src/models/FGAtmosphere.cpp

@@ -102,20 +102,59 @@ bool FGAtmosphere::Run(bool Holding)
  return false;
 }
 
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+// Using pressure in Outer Space between stars in the Milky Way.
+
+double FGAtmosphere::ValidatePressure(double p, const string& msg, bool quiet) const
+{
+ const double MinPressure = ConvertToPSF(1E-15, ePascals);
+ if (p < MinPressure) {
+ if (!quiet) {
+ cerr << msg << " " << p << " is too low." << endl
+ << msg << " is capped to " << MinPressure << endl;
+ }
+ return MinPressure;
+ }
+ return p;
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ // Make sure that the ambient temperature never drops to zero.
+ // According to Wikipedia, 1K is the temperature at the coolest natural place
+ // currently (2023) known in the Universe: the Boomerang Nebula
+
+double FGAtmosphere::ValidateTemperature(double t, const string& msg, bool quiet) const
+{
+ // Using pressure in Outer Space between stars in the Milky Way.
+ const double MinTemperature = ConvertToRankine(1.0, eKelvin);
+ if (t < MinTemperature) {
+ if (!quiet) {
+ cerr << msg << " " << t << " is too low." << endl
+ << msg << " is capped to " << MinTemperature << endl;
+ }
+ return MinTemperature;
+ }
+ return t;
+}
+
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 void FGAtmosphere::Calculate(double altitude)
 {
  FGPropertyNode* node = PropertyManager->GetNode();
+ double t =0.0;
  if (!PropertyManager->HasNode("atmosphere/override/temperature"))
- Temperature = GetTemperature(altitude);
+ t = GetTemperature(altitude);
  else
- Temperature = node->GetDouble("atmosphere/override/temperature");
+ t = node->GetDouble("atmosphere/override/temperature");
+ Temperature = ValidateTemperature(t, "", true);
 
+ double p = 0.0;
  if (!PropertyManager->HasNode("atmosphere/override/pressure"))
- Pressure = GetPressure(altitude);
+ p = GetPressure(altitude);
  else
- Pressure = node->GetDouble("atmosphere/override/pressure");
+ p = node->GetDouble("atmosphere/override/pressure");
+ Pressure = ValidatePressure(p, "", true);
 
  if (!PropertyManager->HasNode("atmosphere/override/density"))
  Density = GetDensity(altitude);
@@ -136,7 +175,7 @@ void FGAtmosphere::SetPressureSL(ePressure unit, double pressure)
 {
  double press = ConvertToPSF(pressure, unit);
 
- SLpressure = press;
+ SLpressure = ValidatePressure(press, "Sea Level pressure");
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -161,7 +200,9 @@ double FGAtmosphere::GetSoundSpeed(double altitude) const
 
 void FGAtmosphere::SetTemperatureSL(double t, eTemperature unit)
 {
- SLtemperature = ConvertToRankine(t, unit);
+ double temp = ConvertToRankine(t, unit);
+
+ SLtemperature = ValidateTemperature(temp, "Sea Level temperature");
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

src/models/FGAtmosphere.h

@@ -254,6 +254,16 @@ class FGAtmosphere : public FGModel {
  /// Converts from PSF (pounds per square foot) to one of several unit systems.
  double ConvertFromPSF(double t, ePressure unit=ePSF) const;
 
+ /// Check that the pressure is within plausible boundaries.
+ /// @param msg Message to display if the pressure is out of boundaries
+ /// @param quiet Don't display the message if set to true
+ double ValidatePressure(double p, const string& msg, bool quiet=false) const;
+
+ /// Check that the temperature is within plausible boundaries.
+ /// @param msg Message to display if the pressure is out of boundaries
+ /// @param quiet Don't display the message if set to true
+ double ValidateTemperature(double t, const string& msg, bool quiet=false) const;
+
  /// @name ISA constants
  //@{
  /// Universal gas constant - ft*lbf/R/mol

src/models/atmosphere/FGStandardAtmosphere.cpp

@@ -224,7 +224,8 @@ double FGStandardAtmosphere::GetPressure(double altitude) const
 
 void FGStandardAtmosphere::SetPressureSL(ePressure unit, double pressure)
 {
- SLpressure = ConvertToPSF(pressure, unit);
+ double p = ConvertToPSF(pressure, unit);
+ SLpressure = ValidatePressure(p, "Sea Level pressure");
  CalculateSLDensity();
  CalculatePressureBreakpoints(SLpressure);
 }
@@ -328,12 +329,12 @@ void FGStandardAtmosphere::SetTemperature(double t, double h, eTemperature unit)
 {
  double targetTemp = ConvertToRankine(t, unit);
  double GeoPotAlt = GeopotentialAltitude(h);
-
- TemperatureBias = targetTemp - GetStdTemperature(h);
+ double bias = targetTemp - GetStdTemperature(h);
 
  if (GeoPotAlt <= GradientFadeoutAltitude)
- TemperatureBias -= TemperatureDeltaGradient * (GradientFadeoutAltitude - GeoPotAlt);
+ bias -= TemperatureDeltaGradient * (GradientFadeoutAltitude - GeoPotAlt);
 
+ SetTemperatureBias(eRankine, bias);
  CalculatePressureBreakpoints(SLpressure);
 
  SLtemperature = GetTemperature(0.0);
@@ -344,10 +345,20 @@ void FGStandardAtmosphere::SetTemperature(double t, double h, eTemperature unit)
 
 void FGStandardAtmosphere::SetTemperatureBias(eTemperature unit, double t)
 {
+ unsigned int numRows = StdAtmosTemperatureTable.GetNumRows();
+ const double minTemperature = StdAtmosTemperatureTable(numRows, 1) - 1.8;
+
  if (unit == eCelsius || unit == eKelvin)
  t *= 1.80; // If temp delta "t" is given in metric, scale up to English
 
  TemperatureBias = t;
+ if (TemperatureBias <= -minTemperature) {
+ cerr << "The temperature bias " << TemperatureBias << " R is too low. "
+ << "It could result in temperatures below the absolute zero." << endl
+ << "Temperature bias is therefore capped to " << -minTemperature << endl;
+ TemperatureBias = -minTemperature;
+ }
+
  CalculatePressureBreakpoints(SLpressure);
 
  SLtemperature = GetTemperature(0.0);
@@ -382,9 +393,20 @@ void FGStandardAtmosphere::SetSLTemperatureGradedDelta(eTemperature unit, double
 
 void FGStandardAtmosphere::SetTemperatureGradedDelta(double deltemp, double h, eTemperature unit)
 {
+ unsigned int numRows = StdAtmosTemperatureTable.GetNumRows();
+ const double minTemperature = StdAtmosTemperatureTable(numRows, 1);
+ const double minDeltaTemperature = minTemperature - StdSLtemperature;
+
  if (unit == eCelsius || unit == eKelvin)
  deltemp *= 1.80; // If temp delta "t" is given in metric, scale up to English
 
+ if (deltemp <= minDeltaTemperature) {
+ cerr << "The temperature delta " << deltemp << " R is too low. "
+ << "It could result in temperatures below the absolute zero." << endl
+ << "Temperature delta is therefore capped to " << minDeltaTemperature << endl;
+ deltemp = minDeltaTemperature;
+ }
+
  TemperatureDeltaGradient = deltemp/(GradientFadeoutAltitude - GeopotentialAltitude(h));
  CalculateLapseRates();
  CalculatePressureBreakpoints(SLpressure);
@@ -584,10 +606,24 @@ void FGStandardAtmosphere::ValidateVaporMassFraction(double h)
 
 void FGStandardAtmosphere::SetDewPoint(eTemperature unit, double dewpoint)
 {
- double altitude = CalculatePressureAltitude(Pressure, 0.0);
- double VaporPressure = CalculateVaporPressure(ConvertToRankine(dewpoint, unit));
- VaporMassFraction = Rdry * VaporPressure / (Rwater * (Pressure - VaporPressure));
- ValidateVaporMassFraction(altitude);
+ double dewPoint_R = ConvertToRankine(dewpoint, unit);
+ constexpr double minDewPoint = -CelsiusToRankine(c) + 1.0;
+
+ if (dewPoint_R <= minDewPoint) {
+ cerr << "The dew point temperature " << dewPoint_R << " is lower than "
+ << minDewPoint << " R." << endl
+ << "Dew point is therefore capped to " << minDewPoint << endl;
+ dewPoint_R = minDewPoint;
+ }
+
+ double VaporPressure = CalculateVaporPressure(dewPoint_R);
+ SetVaporPressure(ePSF, VaporPressure);
+
+ double finalizedDewPoint = GetDewPoint(eRankine);
+ if (finalizedDewPoint < dewPoint_R) {
+ cerr << "Dew point temperature has been capped to " << finalizedDewPoint
+ << endl;
+ }
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -604,7 +640,7 @@ double FGStandardAtmosphere::GetDewPoint(eTemperature to) const
  dewpoint_degC = c*x / (b - x);
  }
 
- return ConvertFromRankine(1.8*(dewpoint_degC + 273.15), to);
+ return ConvertFromRankine(CelsiusToRankine(dewpoint_degC), to);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -613,6 +649,16 @@ void FGStandardAtmosphere::SetVaporPressure(ePressure unit, double Pa)
 {
  double altitude = CalculatePressureAltitude(Pressure, 0.0);
  double VaporPressure = ConvertToPSF(Pa, unit);
+ if (VaporPressure < 0.0) {
+ cerr << "The vapor pressure cannot be negative." << endl
+ << "Vapor pressure is set to 0.0" << endl;
+ VaporPressure = 0.0;
+ } else if (VaporPressure >= Pressure) {
+ cerr << "The vapor pressure " << VaporPressure
+ << " PSF is higher than the ambient pressure." << endl
+ << "Vapor pressure is therefore capped to " << Pressure-1.0 << endl;
+ VaporPressure = Pressure - 1.0;
+ }
  VaporMassFraction = Rdry * VaporPressure / (Rwater * (Pressure - VaporPressure));
  ValidateVaporMassFraction(altitude);
 }
@@ -644,10 +690,18 @@ double FGStandardAtmosphere::GetRelativeHumidity(void) const
 
 void FGStandardAtmosphere::SetRelativeHumidity(double RH)
 {
- double altitude = CalculatePressureAltitude(Pressure, 0.0);
+ if (RH < 0.0) {
+ cerr << "The relative humidity cannot be negative." << endl
+ << "Relative humidity is set to 0%" << endl;
+ RH = 0.0;
+ } else if (RH > 100.0) {
+ cerr << "The relative humidity cannot be higher than 100%." << endl
+ << "Relative humidity is set to 100%" << endl;
+ RH = 100.0;
+ }
+
  double VaporPressure = 0.01*RH*SaturatedVaporPressure;
- VaporMassFraction = Rdry * VaporPressure / (Rwater * (Pressure - VaporPressure));
- ValidateVaporMassFraction(altitude);
+ SetVaporPressure(ePSF, VaporPressure);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -664,6 +718,11 @@ void FGStandardAtmosphere::SetVaporMassFractionPPM(double frac)
  double altitude = CalculatePressureAltitude(Pressure, 0.0);
  VaporMassFraction = frac*1E-6;
  ValidateVaporMassFraction(altitude);
+
+ if (fabs(VaporMassFraction*1E6-frac)>1E-2) {
+ cerr << "The vapor mass fraction " << frac << " has been capped to "
+ << VaporMassFraction*1E6 << "PPM." << endl;
+ }
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tests/TestStdAtmosphere.py

@@ -384,4 +384,38 @@ def test_max_HR(self):
  self.assertAlmostEqual(rhov/rhoa, ppm*1E-6)
  self.assertAlmostEqual(fdm['atmosphere/vapor-fraction-ppm']/ppm, 1.0)
 
+ def test_temperature_lower_limit(self):
+ fdm = self.create_fdm()
+ fdm.load_model('ball')
+ fdm['ic/h-sl-ft'] = 500000.
+ fdm.run_ic()
+
+ fdm['atmosphere/delta-T'] = -4000.
+ fdm.run()
+ self.assertAlmostEqual(fdm['atmosphere/T-R'], 1.8)
+
+ fdm['atmosphere/delta-T'] = 0.0
+ fdm['atmosphere/SL-graded-delta-T'] = -4000.
+ fdm.run()
+ self.assertAlmostEqual(fdm['atmosphere/T-sl-R'], fdm['atmosphere/T-R'])
+
+ fdm['atmosphere/delta-T'] = -4000.
+ fdm.run()
+ self.assertAlmostEqual(fdm['atmosphere/T-sl-R'], 1.8)
+
+ def test_pressure_lower_limit(self):
+ fdm = self.create_fdm()
+ fdm.load_model('ball')
+ fdm['ic/h-sl-ft'] = 1E7
+ fdm.run_ic()
+ self.assertAlmostEqual(fdm['atmosphere/P-psf']*1E17, 2.08854342)
+ self.assertGreater(fdm['atmosphere/pressure-altitude'], 900000.)
+ self.assertLess(fdm['atmosphere/pressure-altitude'], fdm['ic/h-sl-ft'])
+
+ fdm['atmosphere/P-sl-psf'] = 0.0
+ fdm.run()
+ self.assertAlmostEqual(fdm['atmosphere/P-sl-psf']*1E17, 2.08854342)
+ self.assertAlmostEqual(fdm['atmosphere/P-psf']*1E17, 2.08854342)
+
+
 RunTest(TestStdAtmosphere)