From 0c519f82eef071e7b0ade403fdb909fc3e3673ab Mon Sep 17 00:00:00 2001 From: Bertrand Coconnier Date: Thu, 16 Feb 2023 09:22:42 +0100 Subject: [PATCH] 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 | 53 ++++++++++-- src/models/FGAtmosphere.h | 10 +++ .../atmosphere/FGStandardAtmosphere.cpp | 83 ++++++++++++++++--- tests/TestStdAtmosphere.py | 34 ++++++++ 4 files changed, 162 insertions(+), 18 deletions(-) diff --git a/src/models/FGAtmosphere.cpp b/src/models/FGAtmosphere.cpp index 219f9b68b5..baed3bbc0a 100644 --- a/src/models/FGAtmosphere.cpp +++ b/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"); } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/src/models/FGAtmosphere.h b/src/models/FGAtmosphere.h index 8b42bca077..e7f0320587 100644 --- a/src/models/FGAtmosphere.h +++ b/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 diff --git a/src/models/atmosphere/FGStandardAtmosphere.cpp b/src/models/atmosphere/FGStandardAtmosphere.cpp index 1077c52c65..45c9bf55c9 100644 --- a/src/models/atmosphere/FGStandardAtmosphere.cpp +++ b/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; + } } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/tests/TestStdAtmosphere.py b/tests/TestStdAtmosphere.py index 1e432a7424..e15f5df1ae 100644 --- a/tests/TestStdAtmosphere.py +++ b/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)