diff --git a/src/FGFDMExec.cpp b/src/FGFDMExec.cpp index 9d08dbc319..02cd5185d7 100644 --- a/src/FGFDMExec.cpp +++ b/src/FGFDMExec.cpp @@ -470,6 +470,7 @@ void FGFDMExec::LoadInputs(unsigned int idx) Auxiliary->in.SoundSpeed = Atmosphere->GetSoundSpeed(); Auxiliary->in.KinematicViscosity = Atmosphere->GetKinematicViscosity(); Auxiliary->in.DistanceAGL = Propagate->GetDistanceAGL(); + Auxiliary->in.AltitudeASL = Propagate->GetAltitudeASL(); Auxiliary->in.Mass = MassBalance->GetMass(); Auxiliary->in.Tl2b = Propagate->GetTl2b(); Auxiliary->in.Tb2l = Propagate->GetTb2l(); diff --git a/src/FGJSBBase.cpp b/src/FGJSBBase.cpp index 30b39688be..43f98f0d9f 100644 --- a/src/FGJSBBase.cpp +++ b/src/FGJSBBase.cpp @@ -102,72 +102,6 @@ string FGJSBBase::CreateIndexedPropertyName(const string& Property, int index) buf << Property << '[' << index << ']'; return buf.str(); } - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -double FGJSBBase::PitotTotalPressure(double mach, double p) -{ - if (mach < 0) return p; - if (mach < 1) //calculate total pressure assuming isentropic flow - return p*pow((1 + 0.2*mach*mach),3.5); - else { - // shock in front of pitot tube, we'll assume its normal and use - // the Rayleigh Pitot Tube Formula, i.e. the ratio of total - // pressure behind the shock to the static pressure in front of - // the normal shock assumption should not be a bad one -- most supersonic - // aircraft place the pitot probe out front so that it is the forward - // most point on the aircraft. The real shock would, of course, take - // on something like the shape of a rounded-off cone but, here again, - // the assumption should be good since the opening of the pitot probe - // is very small and, therefore, the effects of the shock curvature - // should be small as well. AFAIK, this approach is fairly well accepted - // within the aerospace community - - // The denominator below is zero for Mach ~ 0.38, for which - // we'll never be here, so we're safe - - return p*166.92158009316827*pow(mach,7.0)/pow(7*mach*mach-1,2.5); - } -} - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -// Based on the formulas in the US Air Force Aircraft Performance Flight Testing -// Manual (AFFTC-TIH-99-01). In particular sections 4.6 to 4.8. - -double FGJSBBase::MachFromImpactPressure(double qc, double p) -{ - double A = qc / p + 1; - double M = sqrt(5.0*(pow(A, 1. / 3.5) - 1)); // Equation (4.12) - - if (M > 1.0) - for (unsigned int i = 0; i<10; i++) - M = 0.8812848543473311*sqrt(A*pow(1 - 1.0 / (7.0*M*M), 2.5)); // Equation (4.17) - - return M; -} - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -double FGJSBBase::VcalibratedFromMach(double mach, double p) -{ - double asl = FGAtmosphere::StdDaySLsoundspeed; - double psl = FGAtmosphere::StdDaySLpressure; - double qc = PitotTotalPressure(mach, p) - p; - - return asl * MachFromImpactPressure(qc, psl); -} - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -double FGJSBBase::MachFromVcalibrated(double vcas, double p) -{ - double asl = FGAtmosphere::StdDaySLsoundspeed; - double psl = FGAtmosphere::StdDaySLpressure; - double qc = PitotTotalPressure(vcas / asl, psl) - psl; - - return MachFromImpactPressure(qc, p); -} //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% } // namespace JSBSim diff --git a/src/FGJSBBase.h b/src/FGJSBBase.h index 22327a6bfb..f203b211c0 100644 --- a/src/FGJSBBase.h +++ b/src/FGJSBBase.h @@ -250,41 +250,6 @@ class JSBSIM_API FGJSBBase { return measure*0.3048; } - /** Compute the total pressure in front of the Pitot tube. It uses the - * Rayleigh formula for supersonic speeds (See "Introduction to Aerodynamics - * of a Compressible Fluid - H.W. Liepmann, A.E. Puckett - Wiley & sons - * (1947)" §5.4 pp 75-80) - * @param mach The Mach number - * @param p Pressure in psf - * @return The total pressure in front of the Pitot tube in psf */ - static double PitotTotalPressure(double mach, double p); - - /** Compute the Mach number from the differential pressure (qc) and the - * static pressure. Based on the formulas in the US Air Force Aircraft - * Performance Flight Testing Manual (AFFTC-TIH-99-01). - * @param qc The differential/impact pressure - * @param p Pressure in psf - * @return The Mach number */ - static double MachFromImpactPressure(double qc, double p); - - /** Calculate the calibrated airspeed from the Mach number. Based on the - * formulas in the US Air Force Aircraft Performance Flight Testing - * Manual (AFFTC-TIH-99-01). - * @param mach The Mach number - * @param p Pressure in psf - * @return The calibrated airspeed (CAS) in ft/s - * */ - static double VcalibratedFromMach(double mach, double p); - - /** Calculate the Mach number from the calibrated airspeed.Based on the - * formulas in the US Air Force Aircraft Performance Flight Testing - * Manual (AFFTC-TIH-99-01). - * @param vcas The calibrated airspeed (CAS) in ft/s - * @param p Pressure in psf - * @return The Mach number - * */ - static double MachFromVcalibrated(double vcas, double p); - /** Finite precision comparison. @param a first value to compare @param b second value to compare diff --git a/src/initialization/FGInitialCondition.cpp b/src/initialization/FGInitialCondition.cpp index 2778eed19b..f463825a57 100644 --- a/src/initialization/FGInitialCondition.cpp +++ b/src/initialization/FGInitialCondition.cpp @@ -174,8 +174,7 @@ void FGInitialCondition::SetMachIC(double mach) void FGInitialCondition::SetVcalibratedKtsIC(double vcas) { double altitudeASL = GetAltitudeASLFtIC(); - double pressure = Atmosphere->GetPressure(altitudeASL); - double mach = MachFromVcalibrated(fabs(vcas)*ktstofps, pressure); + double mach = Atmosphere->MachFromVcalibrated(fabs(vcas)*ktstofps, altitudeASL); double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL); SetVtrueFpsIC(mach * soundSpeed); @@ -677,13 +676,12 @@ double FGInitialCondition::GetTerrainElevationFtIC(void) const void FGInitialCondition::SetAltitudeAGLFtIC(double agl) { double altitudeASL = GetAltitudeASLFtIC(); - double pressure = Atmosphere->GetPressure(altitudeASL); double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL); double rho = Atmosphere->GetDensity(altitudeASL); double rhoSL = Atmosphere->GetDensitySL(); double mach0 = vt / soundSpeed; - double vc0 = VcalibratedFromMach(mach0, pressure); + double vc0 = Atmosphere->VcalibratedFromMach(mach0, altitudeASL); double ve0 = vt * sqrt(rho/rhoSL); switch(lastLatitudeSet) { @@ -721,11 +719,10 @@ void FGInitialCondition::SetAltitudeAGLFtIC(double agl) altitudeASL = GetAltitudeASLFtIC(); soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL); rho = Atmosphere->GetDensity(altitudeASL); - pressure = Atmosphere->GetPressure(altitudeASL); switch(lastSpeedSet) { case setvc: - mach0 = MachFromVcalibrated(vc0, pressure); + mach0 = Atmosphere->MachFromVcalibrated(vc0, altitudeASL); SetVtrueFpsIC(mach0 * soundSpeed); break; case setmach: @@ -749,13 +746,12 @@ void FGInitialCondition::SetAltitudeAGLFtIC(double agl) void FGInitialCondition::SetAltitudeASLFtIC(double alt) { double altitudeASL = GetAltitudeASLFtIC(); - double pressure = Atmosphere->GetPressure(altitudeASL); double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL); double rho = Atmosphere->GetDensity(altitudeASL); double rhoSL = Atmosphere->GetDensitySL(); double mach0 = vt / soundSpeed; - double vc0 = VcalibratedFromMach(mach0, pressure); + double vc0 = Atmosphere->VcalibratedFromMach(mach0, altitudeASL); double ve0 = vt * sqrt(rho/rhoSL); switch(lastLatitudeSet) { @@ -825,11 +821,10 @@ void FGInitialCondition::SetAltitudeASLFtIC(double alt) altitudeASL = position.GetGeodAltitude(); soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL); rho = Atmosphere->GetDensity(altitudeASL); - pressure = Atmosphere->GetPressure(altitudeASL); switch(lastSpeedSet) { case setvc: - mach0 = MachFromVcalibrated(vc0, pressure); + mach0 = Atmosphere->MachFromVcalibrated(vc0, altitudeASL); SetVtrueFpsIC(mach0 * soundSpeed); break; case setmach: @@ -959,11 +954,10 @@ double FGInitialCondition::GetBodyWindFpsIC(int idx) const double FGInitialCondition::GetVcalibratedKtsIC(void) const { double altitudeASL = GetAltitudeASLFtIC(); - double pressure = Atmosphere->GetPressure(altitudeASL); double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL); double mach = vt / soundSpeed; - return fpstokts * VcalibratedFromMach(mach, pressure); + return fpstokts * Atmosphere->VcalibratedFromMach(mach, altitudeASL); } //****************************************************************************** diff --git a/src/models/FGAtmosphere.cpp b/src/models/FGAtmosphere.cpp index 4816402e74..2c7b29c699 100644 --- a/src/models/FGAtmosphere.cpp +++ b/src/models/FGAtmosphere.cpp @@ -51,11 +51,9 @@ namespace JSBSim { CLASS IMPLEMENTATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ -// Atmosphere constants in British units converted from the SI values specified in the -// ISA document - https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19770009539.pdf -const double FGAtmosphere::StdDaySLsoundspeed = sqrt(SHRatio*Reng0*StdDaySLtemperature); - -FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex) +FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) + : FGModel(fdmex), + StdDaySLsoundspeed(sqrt(SHRatio*Reng0*StdDaySLtemperature)) { Name = "FGAtmosphere"; @@ -283,6 +281,8 @@ double FGAtmosphere::ConvertToPSF(double p, ePressure unit) const return targetPressure; } +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + double FGAtmosphere::ConvertFromPSF(double p, ePressure unit) const { double targetPressure=0; // Pressure @@ -309,6 +309,79 @@ double FGAtmosphere::ConvertFromPSF(double p, ePressure unit) const //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +double FGAtmosphere::PitotTotalPressure(double mach, double p) const +{ + constexpr double a = (SHRatio-1.0) / 2.0; + constexpr double b = SHRatio / (SHRatio-1.0); + constexpr double c = 2.0*b; + constexpr double d = 1.0 / (SHRatio-1.0); + const double coeff = pow(0.5*(SHRatio+1.0), b)*pow((SHRatio+1.0)/(SHRatio-1.0), d); + + if (mach < 0) return p; + if (mach < 1) //calculate total pressure assuming isentropic flow + return p*pow((1.0 + a*mach*mach), b); + else { + // Shock in front of pitot tube, we'll assume its normal and use the + // Rayleigh Pitot Tube Formula, i.e. the ratio of total pressure behind the + // shock to the static pressure in front of the normal shock assumption + // should not be a bad one -- most supersonic aircraft place the pitot probe + // out front so that it is the forward most point on the aircraft. + // The real shock would, of course, take on something like the shape of a + // rounded-off cone but, here again, the assumption should be good since the + // opening of the pitot probe is very small and, therefore, the effects of + // the shock curvature should be small as well. AFAIK, this approach is + // fairly well accepted within the aerospace community + + // The denominator below is zero for Mach ~ 0.38, for which + // we'll never be here, so we're safe + + return p*coeff*pow(mach, c)/pow(c*mach*mach-1.0, d); + } +} + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +// Based on the formulas in the US Air Force Aircraft Performance Flight Testing +// Manual (AFFTC-TIH-99-01). In particular sections 4.6 to 4.8. + +double FGAtmosphere::MachFromImpactPressure(double qc, double p) const +{ + constexpr double a = 2.0/(SHRatio-1.0); + constexpr double b = (SHRatio-1.0)/SHRatio; + constexpr double c = 2.0/b; + constexpr double d = 0.5*a; + const double coeff = pow(0.5*(SHRatio+1.0), -0.25*c)*pow(0.5*(SHRatio+1.0)/SHRatio, -0.5*d); + + double A = qc / p + 1; + double M = sqrt(a*(pow(A, b) - 1.0)); // Equation (4.12) + + if (M > 1.0) + for (unsigned int i = 0; i<10; i++) + M = coeff*sqrt(A*pow(1 - 1.0 / (c*M*M), d)); // Equation (4.17) + + return M; +} + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +double FGAtmosphere::VcalibratedFromMach(double mach, double altitude) const +{ + double p = GetPressure(altitude); + double qc = PitotTotalPressure(mach, p) - p; + return StdDaySLsoundspeed * MachFromImpactPressure(qc, StdDaySLpressure); +} + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +double FGAtmosphere::MachFromVcalibrated(double vcas, double altitude) const +{ + double p = GetPressure(altitude); + double qc = PitotTotalPressure(vcas / StdDaySLsoundspeed, StdDaySLpressure) - StdDaySLpressure; + return MachFromImpactPressure(qc, p); +} + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + void FGAtmosphere::bind(void) { PropertyManager->Tie("atmosphere/T-R", this, &FGAtmosphere::GetTemperature); diff --git a/src/models/FGAtmosphere.h b/src/models/FGAtmosphere.h index 6a4f70feca..f5e545d7f7 100644 --- a/src/models/FGAtmosphere.h +++ b/src/models/FGAtmosphere.h @@ -206,13 +206,48 @@ class JSBSIM_API FGAtmosphere : public FGModel { virtual double GetPressureAltitude() const {return PressureAltitude;} + /** Compute the total pressure in front of the Pitot tube. It uses the + * Rayleigh formula for supersonic speeds (See "Introduction to Aerodynamics + * of a Compressible Fluid - H.W. Liepmann, A.E. Puckett - Wiley & sons + * (1947)" §5.4 pp 75-80) + * @param mach The Mach number + * @param p Pressure in psf + * @return The total pressure in front of the Pitot tube in psf */ + double PitotTotalPressure(double mach, double p) const; + + /** Compute the Mach number from the differential pressure (qc) and the + * static pressure. Based on the formulas in the US Air Force Aircraft + * Performance Flight Testing Manual (AFFTC-TIH-99-01). + * @param qc The differential/impact pressure + * @param p Pressure in psf + * @return The Mach number */ + double MachFromImpactPressure(double qc, double p) const; + + /** Calculate the calibrated airspeed from the Mach number. Based on the + * formulas in the US Air Force Aircraft Performance Flight Testing + * Manual (AFFTC-TIH-99-01). + * @param mach The Mach number + * @param altitude The altitude above sea level (ASL) in feet. + * @return The calibrated airspeed (CAS) in ft/s + * */ + double VcalibratedFromMach(double mach, double altitude) const; + + /** Calculate the Mach number from the calibrated airspeed.Based on the + * formulas in the US Air Force Aircraft Performance Flight Testing + * Manual (AFFTC-TIH-99-01). + * @param vcas The calibrated airspeed (CAS) in ft/s + * @param altitude The altitude above sea level (ASL) in feet. + * @return The Mach number + * */ + double MachFromVcalibrated(double vcas, double altitude) const; + struct Inputs { double altitudeASL; } in; static constexpr double StdDaySLtemperature = 518.67; static constexpr double StdDaySLpressure = 2116.228; - static const double StdDaySLsoundspeed; + const double StdDaySLsoundspeed; protected: // Sea level conditions diff --git a/src/models/FGAuxiliary.cpp b/src/models/FGAuxiliary.cpp index c0057788c8..70ddca6b64 100644 --- a/src/models/FGAuxiliary.cpp +++ b/src/models/FGAuxiliary.cpp @@ -47,6 +47,7 @@ INCLUDES #include "FGFDMExec.h" #include "input_output/FGPropertyManager.h" #include "FGInertial.h" +#include "FGAtmosphere.h" using namespace std; @@ -193,10 +194,10 @@ bool FGAuxiliary::Run(bool Holding) tat = in.Temperature*(1 + 0.2*Mach*Mach); // Total Temperature, isentropic flow tatc = RankineToCelsius(tat); - pt = PitotTotalPressure(Mach, in.Pressure); + pt = FDMExec->GetAtmosphere()->PitotTotalPressure(Mach, in.Pressure); if (abs(Mach) > 0.0) { - vcas = VcalibratedFromMach(Mach, in.Pressure); + vcas = FDMExec->GetAtmosphere()->VcalibratedFromMach(Mach, in.AltitudeASL); veas = sqrt(2 * qbar / in.DensitySL); } else diff --git a/src/models/FGAuxiliary.h b/src/models/FGAuxiliary.h index f1e3fa1b34..4a4f12c3cb 100644 --- a/src/models/FGAuxiliary.h +++ b/src/models/FGAuxiliary.h @@ -260,6 +260,7 @@ class JSBSIM_API FGAuxiliary : public FGModel { double SoundSpeed; double KinematicViscosity; double DistanceAGL; + double AltitudeASL; double Wingspan; double Wingchord; double StandardGravity; diff --git a/tests/unit_tests/FGAtmosphereTest.h b/tests/unit_tests/FGAtmosphereTest.h index 2bf5fb78b2..b862d6dd9a 100644 --- a/tests/unit_tests/FGAtmosphereTest.h +++ b/tests/unit_tests/FGAtmosphereTest.h @@ -86,6 +86,8 @@ class FGAtmosphereTest : public CxxTest::TestSuite TS_ASSERT_EQUALS(atm.GetAbsoluteViscosity(), 0.0); TS_ASSERT_EQUALS(atm.GetKinematicViscosity(), 0.0); + + FGJSBBase::debug_lvl = 0; } void testDefaultValuesAfterInit() @@ -611,4 +613,152 @@ class FGAtmosphereTest : public CxxTest::TestSuite TS_ASSERT_EQUALS(atm.GetTemperature(1000.), 1.8); TS_ASSERT_DELTA(atm.GetPressure(1000.)*psftopa*1e15, 1.0, 1e-5); } + + void testPitotTotalPressure() { + FGFDMExec fdmex; + auto atm = DummyAtmosphere(&fdmex, -1.0, -100.0); + atm.InitModel(); + + // Ambient conditions far upstream (i.e. upstream the normal schock + // in supersonic flight) + double p1 = atm.GetPressureSL(); + double t1 = atm.GetTemperatureSL(); + double rho1 = atm.GetDensitySL(); + constexpr double Cp = gama*R/(gama-1.0); + + // Based on formulas from Modern Compressible Flow (3rd edition) + // - John D. Anderson + for(double M1=0; M1<3.0; M1+=0.25) { + double a1 = sqrt(gama*R*t1); + double u1 = M1*a1; + // Total temperature + double T0 = t1+u1*u1/(2.0*Cp); + // Compute conditions downstream (at the pitot tube) + double u2 = u1; + if (M1 >= 1.0){ + // Assess the normal shock effect knowing that a_star=u1*u2 + double a_star = sqrt((a1*a1/(gama-1.0)+0.5*u1*u1)*2*(gama-1.0)/(gama+1.0)); // equation (3.32) p.81 + u2 = a_star*a_star/u1;// equation (3.47) p.89 + } + double t2 = T0-u2*u2/(2*Cp); + double P2 = atm.PitotTotalPressure(M1, p1); + double p2 = P2*pow(t2/T0, gama/(gama-1.0)); + double rho2 = p2/(R*t2); + + // mass conservation + TS_ASSERT_DELTA(rho1*u1, rho2*u2, epsilon); + // momentum conservation + TS_ASSERT_DELTA(p1+rho1*u1*u1, p2+rho2*u2*u2, 1000.*epsilon); + // energy conservation + TS_ASSERT_DELTA(Cp*t1+0.5*u1*u1, Cp*t2+0.5*u2*u2, epsilon); + } + } + + void testMachFromImpactPressure() { + FGFDMExec fdmex; + auto atm = DummyAtmosphere(&fdmex, -1.0, -100.0); + atm.InitModel(); + + // Ambient conditions far upstream (i.e. upstream the normal schock + // in supersonic flight) + double p1 = atm.GetPressureSL(); + double t1 = atm.GetTemperatureSL(); + double rho1 = atm.GetDensitySL(); + constexpr double Cp = gama*R/(gama-1.0); + + // Based on formulas from Modern Compressible Flow (3rd edition) + // - John D. Anderson + for(double M1=0; M1<3.0; M1+=0.25) { + double a1 = sqrt(gama*R*t1); + double u1 = M1*a1; + // Total temperature + double T0 = t1+u1*u1/(2.0*Cp); + // Compute conditions downstream (at the pitot tube) + double u2 = u1; + if (M1 >= 1.0) { + // Assess the normal shock effect knowing that a_star=u1*u2 + double a_star = sqrt((a1*a1/(gama-1.0)+0.5*u1*u1)*2*(gama-1.0)/(gama+1.0)); // equation (3.32) p.81 + u2 = a_star*a_star/u1;// equation (3.47) p.89 + } + double t2 = T0-u2*u2/(2*Cp); + double rho2 = M1 == 0.0 ? rho1 : rho1*u1/u2; + double p2 = rho2*R*t2; + double P2 = p2*pow(T0/t2, gama/(gama-1.0)); + double mach1 = atm.MachFromImpactPressure(P2-p1, p1); + double a2 = sqrt(gama*R*t2); + double M2 = u2/a2; + double mach2 = atm.MachFromImpactPressure(P2-p2, p2); + + // mass conservation + TS_ASSERT_DELTA(rho1*u1, rho2*u2, epsilon); + // momentum conservation + TS_ASSERT_DELTA(p1+rho1*u1*u1, p2+rho2*u2*u2, 1000.*epsilon); + // energy conservation + TS_ASSERT_DELTA(Cp*t1+0.5*u1*u1, Cp*t2+0.5*u2*u2, epsilon); + // Check the Mach computations + TS_ASSERT_DELTA(mach1, M1, 1e-7); + TS_ASSERT_DELTA(mach2, M2, 1e-7); + } + } + + void testCASConversion() { + FGFDMExec fdmex; + auto atm = DummyAtmosphere(&fdmex, -0.1, -1.0); + atm.InitModel(); + + // Ambient conditions far upstream (i.e. upstream the normal schock + // in supersonic flight) + double t1 = atm.GetTemperatureSL(); + + TS_ASSERT_DELTA(atm.VcalibratedFromMach(0.0, 0.0), 0.0, epsilon); + TS_ASSERT_DELTA(atm.MachFromVcalibrated(0.0, 0.0), 0.0, epsilon); + + // Check that VCAS match the true airspeed at sea level + for(double M1=0.1; M1<3.0; M1+=0.25) { + double u1 = M1*sqrt(gama*R*t1); + TS_ASSERT_DELTA(atm.VcalibratedFromMach(M1, 0.0)/u1, 1.0, 1e-7); + TS_ASSERT_DELTA(atm.MachFromVcalibrated(u1, 0.0)/M1, 1.0, 1e-7); + } + + // Check the VCAS computation at an altitude of 1000 ft + double asl = atm.GetSoundSpeedSL(); + + TS_ASSERT_DELTA(atm.VcalibratedFromMach(0.0, 1000.), 0.0, epsilon); + TS_ASSERT_DELTA(atm.MachFromVcalibrated(0.0, 1000.), 0.0, epsilon); + + for(double M=0.1; M<3.0; M+=0.25) { + double vcas = M*asl; + double M1 = atm.MachFromVcalibrated(vcas, 1000.); + TS_ASSERT_DELTA(atm.VcalibratedFromMach(M1, 1000.)/vcas, 1.0, 1e-7); + } + + double psl = atm.GetPressureSL(); + double p1 = atm.GetPressure(1000.); + t1 = atm.GetTemperature(1000.); + double rho1 = atm.GetDensity(1000.); + constexpr double Cp = gama*R/(gama-1.0); + + // Based on formulas from Modern Compressible Flow (3rd edition) + // - John D. Anderson + for(double M1=0.1; M1<3.0; M1+=0.25) { + double a1 = sqrt(gama*R*t1); + double u1 = M1*a1; + // Total temperature + double T0 = t1+u1*u1/(2.0*Cp); + // Compute conditions downstream (at the pitot tube) + double u2 = u1; + if (M1 >= 1.0) { + // Assess the normal shock effect knowing that a_star=u1*u2 + double a_star = sqrt((a1*a1/(gama-1.0)+0.5*u1*u1)*2*(gama-1.0)/(gama+1.0)); // equation (3.32) p.81 + u2 = a_star*a_star/u1;// equation (3.47) p.89 + } + double t2 = T0-u2*u2/(2*Cp); + double rho2 = M1 == 0.0 ? rho1 : rho1*u1/u2; + double p2 = rho2*R*t2; + double P2 = p2*pow(T0/t2, gama/(gama-1.0)); + double mach = atm.MachFromImpactPressure(P2-p1, psl); + + TS_ASSERT_DELTA(atm.VcalibratedFromMach(M1, 1000.)/(mach*asl), 1.0, 1e-8); + } + } }; diff --git a/tests/unit_tests/FGJSBBaseTest.h b/tests/unit_tests/FGJSBBaseTest.h index 732fe72cbf..dc456479dc 100644 --- a/tests/unit_tests/FGJSBBaseTest.h +++ b/tests/unit_tests/FGJSBBaseTest.h @@ -6,19 +6,6 @@ class FGJSBBaseTest : public CxxTest::TestSuite, public JSBSim::FGJSBBase { public: - void testCASConversion() { - double p = 2116.228; - TS_ASSERT_EQUALS(VcalibratedFromMach(-0.1, p), 0.0); - TS_ASSERT_EQUALS(VcalibratedFromMach(0, p), 0.0); - TS_ASSERT_DELTA(VcalibratedFromMach(0.5, p), 558.2243, 1E-4); - TS_ASSERT_DELTA(VcalibratedFromMach(1.0, p), 1116.4486, 1E-4); - TS_ASSERT_DELTA(VcalibratedFromMach(1.5, p), 1674.6728, 1E-4); - TS_ASSERT_EQUALS(MachFromVcalibrated(0.0, p), 0.0); - TS_ASSERT_DELTA(MachFromVcalibrated(558.2243, p), 0.5, 1E-4); - TS_ASSERT_DELTA(MachFromVcalibrated(1116.4486, p), 1.0, 1E-4); - TS_ASSERT_DELTA(MachFromVcalibrated(1674.6728, p), 1.5, 1E-4); - } - void testNumericRoutines() { double dx = 1.0; float fx = 1.0;