Discrepancy in Air Density and Atmospheric Pressure Calculation for Non-Standard Temperature

[viniciusrottolo]

I'm submitting a ...

• bug report
• feature request
• support request

Describe the issue
There is a discrepancy in the calculation of air density and atmospheric pressure when using the temperature deviation from standard atmosphere (atmosphere/delta-T) parameter, especially at altitudes above sea level. Notably, the atmospheric pressure value should remain constant for a fixed altitude regardless of temperature deviations, but JSBSim produces different pressure values.

What is the expected behavior?
The air density and atmospheric pressure values computed by JSBSim should closely align with standard reference values for the given temperature and altitude conditions, such as those provided by the Atmospheric Calculator.

What is the motivation / use case for changing the behavior?
Correct air density values are critical for accurate aerodynamic and flight dynamics simulations, especially for scenarios where temperature deviations from the standard atmosphere are involved.

code

import jsbsim

# Global variables that must be modified to match your particular need
# The aircraft name
# Note - It should match the exact spelling of the model file
AIRCRAFT_NAME = "737"
# Path to JSBSim files, location of the folders "aircraft", "engines" and "systems"
PATH_TO_JSBSIM_FILES = "../.."

# Avoid flooding the console with log messages
jsbsim.FGJSBBase().debug_lvl = 0

fdm = jsbsim.FGFDMExec(PATH_TO_JSBSIM_FILES)

# Load the aircraft model

# Load the aircraft model
fdm.load_model(AIRCRAFT_NAME)

fdm['atmosphere/delta-T'] =  0*1.8 # dISA °C to dISA R
fdm['ic/h-sl-ft'] = 0
fdm.run_ic()

print(" STANDARD ATMOSPHERE \n",
      f"Altitude: {fdm['position/h-sl-ft']:.0f} ft \n",
      f"Temperature_JSBSIM: {(fdm['atmosphere/T-R']-491.67)*5/9:.0f} °C \n",
      f"Density_JSBSIM: {fdm['atmosphere/rho-slugs_ft3']:.7f} slug/ft3 \n",
      f"Pressure_JSBSIM: {fdm['atmosphere/P-psf']:.2f} \n",
      "Temperature_ref:  15 °C \n",
      "Density_ref: 0.0023772 slug/ft3 \n",
      "Pressure_ref: 2116.20 psf \n")

fdm['atmosphere/delta-T'] =  25*1.8 # dISA °C to dISA R
fdm['ic/h-sl-ft'] = 0

fdm.run_ic() # Initialize the aircraft with initial conditions  
print(f" Altitude: {fdm['position/h-sl-ft']:.0f} ft \n",
      f"Temperature_JSBSIM: {(fdm['atmosphere/T-R']-491.67)*5/9:.0f} °C \n",
      f"Density_JSBSIM: {fdm['atmosphere/rho-slugs_ft3']:.7f} slug/ft3 \n",
      f"Pressure_JSBSIM: {fdm['atmosphere/P-psf']:.2f} \n",
      "Temperature_ref:  40 °C \n",
      "Density_ref: 0.0021874 slug/ft3 \n",
      "Pressure_ref: 2116.2 psf \n")

fdm['ic/h-sl-ft'] = 10000
fdm.run_ic() # Initialize the aircraft with initial conditions  

print(f" Altitude: {fdm['position/h-sl-ft']:.0f} ft \n",
      f"Temperature_JSBSIM: {(fdm['atmosphere/T-R']-491.67)*5/9:.4f} °C \n",
      f"Density_JSBSIM: {fdm['atmosphere/rho-slugs_ft3']:.7f} slug/ft3 \n",
      f"Pressure_JSBSIM: {fdm['atmosphere/P-psf']:.2f} \n",
      "Temperature_ref:  20.1975 °C \n",
      "Density_ref: 0.0016061 slug/ft3 \n",
      "Pressure_ref: 1455.60 psf \n")

fdm['ic/h-sl-ft'] = 20000
fdm.run_ic() # Initialize the aircraft with initial conditions  

print(f" Altitude: {fdm['position/h-sl-ft']:.0f} ft \n",
      f"Temperature_JSBSIM: {(fdm['atmosphere/T-R']-491.67)*5/9:.4f} °C \n",
      f"Density_JSBSIM: {fdm['atmosphere/rho-slugs_ft3']:.7f} slug/ft3 \n",
      f"Pressure_JSBSIM: {fdm['atmosphere/P-psf']:.2f} \n",
      "Temperature_ref:  0.4140 °C \n",
      "Density_ref: 0.0011516 slug/ft3 \n",
      "Pressure_ref: 973.27 psf \n")

fdm['ic/h-sl-ft'] = 30000
fdm.run_ic() # Initialize the aircraft with initial conditions  

print(f" Altitude: {fdm['position/h-sl-ft']:.0f} ft \n",
      f"Temperature_JSBSIM: {(fdm['atmosphere/T-R']-491.67)*5/9:.4f} °C \n",
      f"Density_JSBSIM: {fdm['atmosphere/rho-slugs_ft3']:.7f} slug/ft3 \n",
      f"Pressure_JSBSIM: {fdm['atmosphere/P-psf']:.2f} \n",
      "Temperature_ref: -19.3506 °C \n",
      "Density_ref: 0.00080305 slug/ft3 \n",
      "Pressure_ref: 629.67 psf \n")

fdm['ic/h-sl-ft'] = 40000
fdm.run_ic() # Initialize the aircraft with initial conditions  

print(f" Altitude: {fdm['position/h-sl-ft']:.0f} ft \n",
      f"Temperature_JSBSIM: {(fdm['atmosphere/T-R']-491.67)*5/9:.4f} °C \n",
      f"Density_JSBSIM: {fdm['atmosphere/rho-slugs_ft3']:.7f} slug/ft3 \n",
      f"Pressure_JSBSIM: {fdm['atmosphere/P-psf']:.2f} \n",
      "Temperature_ref:  -31.5 °C \n",
      "Density_ref: 0.00052659 slug/ft3 \n",
      "Pressure_ref: 393.13 psf \n")

fdm['atmosphere/delta-T'] =  35*1.8 # dISA °C to dISA R
fdm.run_ic() # Initialize the aircraft with initial conditions  

print(f" Altitude: {fdm['position/h-sl-ft']:.0f} ft \n",
      f"Temperature_JSBSIM: {(fdm['atmosphere/T-R']-491.67)*5/9:.4f} °C \n",
      f"Density_JSBSIM: {fdm['atmosphere/rho-slugs_ft3']:.7f} slug/ft3 \n",
      f"Pressure_JSBSIM: {fdm['atmosphere/P-psf']:.2f} \n",
      "Temperature_ref:  -21.5 °C \n",
      "Density_ref: 0.00050566 slug/ft3 \n",
      "Pressure_ref: 393.13 psf \n")

output

 STANDARD ATMOSPHERE 
 Altitude: 0 ft 
 Temperature_JSBSIM: 15 °C 
 Density_JSBSIM: 0.0023769 slug/ft3 
 Pressure_JSBIM: 2116.23 
 Temperature_ref:  15 °C 
 Density_ref: 0.0023772 slug/ft3 
 Pressure_ref: 2116.20 psf 

 Altitude: 0 ft 
 Temperature_JSBSIM: 40 °C 
 Density_JSBSIM: 0.0021872 slug/ft3 
 Pressure_JSBSIM: 2116.23 
 Temperature_ref:  40 °C 
 Density_ref: 0.0021874 slug/ft3 
 Pressure_ref: 2116.2 psf 

 Altitude: 10000 ft 
 Temperature_JSBSIM: 20.1975 °C 
 Density_JSBSIM: 0.0016563 slug/ft3 
 Pressure_JSBSIM: 1501.25 
 Temperature_ref:  20.1975 °C 
 Density_ref: 0.0016061 slug/ft3 
 Pressure_ref: 1455.60 psf 

 Altitude: 20000 ft 
 Temperature_JSBSIM: 0.4140 °C 
 Density_JSBSIM: 0.0012305 slug/ft3 
 Pressure_JSBSIM: 1040.10 
 Temperature_ref:  0.4140 °C 
 Density_ref: 0.0011516 slug/ft3 
 Pressure_ref: 973.27 psf 

 Altitude: 30000 ft 
 Temperature_JSBSIM: -19.3506 °C 
 Density_JSBSIM: 0.0008943 slug/ft3 
 Pressure_JSBSIM: 701.29 
 Temperature_ref: -19.3506 °C 
 Density_ref: 0.00080305 slug/ft3 
 Pressure_ref: 629.67 psf 

 Altitude: 40000 ft 
 Temperature_JSBSIM: -31.5000 °C 
 Density_JSBSIM: 0.0006153 slug/ft3 
 Pressure_JSBSIM: 459.38 
 Temperature_ref:  -31.5 °C 
 Density_ref: 0.00052659 slug/ft3 
 Pressure_ref: 393.13 psf 

 Altitude: 40000 ft 
 Temperature_JSBSIM: -21.5000 °C 
 Density_JSBSIM: 0.0006239 slug/ft3 
 Pressure_JSBSIM: 485.10 
 Temperature_ref:  -21.5 °C 
 Density_ref: 0.00050566 slug/ft3 
 Pressure_ref: 393.13 psf 

[Warning39]

Non-standard (hot or cold) days are modeled by adding a specified temperature delta to the standard temperature at altitude, but pressure is taken as the standard day value. Density and viscosity are recalculated at the resultant temperature and pressure using the ideal gas equation of state.

From - https://en.wikipedia.org/wiki/International_Standard_Atmosphere

So it looks like TemperatureBias shouldn't be used here in CalculatePressureBreakpoints().

jsbsim/src/models/atmosphere/FGStandardAtmosphere.cpp

Lines 486 to 488 in 318c804

	double Tmb = BaseTemp
	+ TemperatureBias
	+ (GradientFadeoutAltitude - BaseAlt)*TemperatureDeltaGradient;

And instead the TemperatureBias needs to be used in:

jsbsim/src/models/atmosphere/FGStandardAtmosphere.cpp

Lines 523 to 528 in 318c804

	void FGStandardAtmosphere::CalculateStdDensityBreakpoints()
	{
	StdDensityBreakpoints.clear();
	for (unsigned int i = 0; i < StdPressureBreakpoints.size(); i++)
	StdDensityBreakpoints.push_back(StdPressureBreakpoints[i] / (Rdry * StdAtmosTemperatureTable(i + 1, 1)));
	}

And also in:

jsbsim/src/models/atmosphere/FGStandardAtmosphere.cpp

Lines 238 to 242 in 318c804

	void FGStandardAtmosphere::CalculateSLSoundSpeedAndDensity(void)
	{
	SLsoundspeed = sqrt(SHRatio*Reng*SLtemperature);
	CalculateSLDensity();
	}

And:

jsbsim/src/models/atmosphere/FGStandardAtmosphere.h

Line 352 in 318c804

void CalculateSLDensity(void) { SLdensity = SLpressure / (Reng * SLtemperature); }

[bcoconni]

The air density and atmospheric pressure values computed by JSBSim should closely align with standard reference values for the given temperature and altitude condition

Well, according to ISA 1976 chapter 1.2.2 and its equations (4) and (5), we have:
$dP=-g\rho dZ$
and
$\dfrac{dP}{P}=-\dfrac{-g}{RT}dZ$
So if the temperature is changed and the pressure is not, then the equations above are no longer met.

Unfortunately, the ISA document does not explicitly tell how a temperature deviation from standard should be handled. However I would say that our current implementation is making sense because it complies with ISA 1976 chapters 1.2.2 and 1.3.1. Unless proven otherwise.

[bcoconni]

Non-standard (hot or cold) days are modeled by adding a specified temperature delta to the standard temperature at altitude, but pressure is taken as the standard day value. Density and viscosity are recalculated at the resultant temperature and pressure using the ideal gas equation of state.

From - https://en.wikipedia.org/wiki/International_Standard_Atmosphere

So it looks like TemperatureBias shouldn't be used here in CalculatePressureBreakpoints().

jsbsim/src/models/atmosphere/FGStandardAtmosphere.cpp

Lines 486 to 488 in 318c804

	double Tmb = BaseTemp
	+ TemperatureBias
	+ (GradientFadeoutAltitude - BaseAlt)*TemperatureDeltaGradient;

@Warning39, if the pressure is independent of the temperature I don't see why the delta gradient should be kept ?

And instead the TemperatureBias needs to be used in:

jsbsim/src/models/atmosphere/FGStandardAtmosphere.cpp

Lines 523 to 528 in 318c804

	void FGStandardAtmosphere::CalculateStdDensityBreakpoints()
	{
	StdDensityBreakpoints.clear();
	for (unsigned int i = 0; i < StdPressureBreakpoints.size(); i++)
	StdDensityBreakpoints.push_back(StdPressureBreakpoints[i] / (Rdry * StdAtmosTemperatureTable(i + 1, 1)));
	}

I don't see why standard density breakpoints should depend on TemperatureBias ?

And also in:

jsbsim/src/models/atmosphere/FGStandardAtmosphere.cpp

Lines 238 to 242 in 318c804

	void FGStandardAtmosphere::CalculateSLSoundSpeedAndDensity(void)
	{
	SLsoundspeed = sqrt(SHRatio*Reng*SLtemperature);
	CalculateSLDensity();
	}

And:

jsbsim/src/models/atmosphere/FGStandardAtmosphere.h

Line 352 in 318c804

void CalculateSLDensity(void) { SLdensity = SLpressure / (Reng * SLtemperature); }

SLtemperature already accounts for TemperatureBias because it is computed from FGStandardAtmosphere::GetTemperature which adds TemperatureBias to its result.

jsbsim/src/models/atmosphere/FGStandardAtmosphere.cpp

Line 343 in d1bedfb

SLtemperature = GetTemperature(0.0);

jsbsim/src/models/atmosphere/FGStandardAtmosphere.cpp

Line 376 in d1bedfb

SLtemperature = GetTemperature(0.0);

jsbsim/src/models/atmosphere/FGStandardAtmosphere.cpp

Line 429 in d1bedfb

SLtemperature = GetTemperature(0.0);

jsbsim/src/models/atmosphere/FGStandardAtmosphere.cpp

Lines 248 to 272 in d1bedfb

	double FGStandardAtmosphere::GetTemperature(double altitude) const
	{
	double GeoPotAlt = GeopotentialAltitude(altitude);
	double T;
	if (GeoPotAlt >= 0.0) {
	T = StdAtmosTemperatureTable.GetValue(GeoPotAlt);
	if (GeoPotAlt <= GradientFadeoutAltitude)
	T -= TemperatureDeltaGradient * GeoPotAlt;
	}
	else {
	// We don't need to add TemperatureDeltaGradient*GeoPotAlt here because
	// the lapse rate vector already accounts for the temperature gradient.
	T = StdAtmosTemperatureTable.GetValue(0.0) + GeoPotAlt*LapseRates[0];
	}
	T += TemperatureBias;
	if (GeoPotAlt <= GradientFadeoutAltitude)
	T += TemperatureDeltaGradient * GradientFadeoutAltitude;
	return T;
	}

[Warning39]

Unfortunately, the ISA document does not explicitly tell how a temperature deviation from standard should be handled.

Yep, just glanced through it now. Previously when I did a quick internet search I came across the Wikipedia entry:

Non-standard (hot or cold) days are modeled by adding a specified temperature delta to the standard temperature at altitude, but pressure is taken as the standard day value. Density and viscosity are recalculated at the resultant temperature and pressure using the ideal gas equation of state.

From - https://en.wikipedia.org/wiki/International_Standard_Atmosphere and I assumed the description was based on a description in the ISA document.

The online ISA calculator - https://aerospaceweb.org/design/scripts/atmosphere/ also appears to implement the description from the Wikipedia entry.

Doing a bit more searching now I came across - https://www.pdas.com/atmosnonstd.html

The committee that defined the standard atmosphere did not define any such non-standard atmospheres. As a result, there may be more than one such definition in practice today.

Regardless of the definition, all of the non-standard atmosphere definitions that I have examined are simply a different table of temperature versus altitude. At a given altitude, the computational algorithm is

• Compute the pressure at this altitude from the standard atmosphere. Assume this does not change.
• Compute the temperature from the revised schedule.
• Compute density from the pressure and temperature using the perfect gas law.

So the initial question I guess is whether JSBSim should follow this same approach. Which seems somewhat common?