diff --git a/src/FGFDMExec.cpp b/src/FGFDMExec.cpp index 2a8357365b..755f348d2c 100644 --- a/src/FGFDMExec.cpp +++ b/src/FGFDMExec.cpp @@ -390,7 +390,6 @@ 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(); @@ -534,11 +533,12 @@ void FGFDMExec::LoadInputs(unsigned int idx) void FGFDMExec::LoadPlanetConstants(void) { - Propagate->in.vOmegaPlanet = Inertial->GetOmegaPlanet(); - Accelerations->in.vOmegaPlanet = Inertial->GetOmegaPlanet(); - Propagate->in.SemiMajor = Inertial->GetSemimajor(); - Propagate->in.SemiMinor = Inertial->GetSemiminor(); - Auxiliary->in.StandardGravity = Inertial->GetStandardGravity(); + Propagate->in.vOmegaPlanet = Inertial->GetOmegaPlanet(); + Accelerations->in.vOmegaPlanet = Inertial->GetOmegaPlanet(); + Propagate->in.SemiMajor = Inertial->GetSemimajor(); + Propagate->in.SemiMinor = Inertial->GetSemiminor(); + Auxiliary->in.StandardGravity = Inertial->GetStandardGravity(); + Auxiliary->in.StdDaySLsoundspeed = Atmosphere->StdDaySLsoundspeed; } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/src/initialization/FGInitialCondition.cpp b/src/initialization/FGInitialCondition.cpp index 23c21df9d4..949e0c8f2e 100644 --- a/src/initialization/FGInitialCondition.cpp +++ b/src/initialization/FGInitialCondition.cpp @@ -48,6 +48,7 @@ INCLUDES #include "models/FGInertial.h" #include "models/FGAtmosphere.h" #include "models/FGAccelerations.h" +#include "models/FGAuxiliary.h" #include "input_output/FGXMLFileRead.h" #include "FGTrim.h" @@ -64,6 +65,7 @@ FGInitialCondition::FGInitialCondition(FGFDMExec *FDMExec) : fdmex(FDMExec) if(FDMExec) { Atmosphere=fdmex->GetAtmosphere(); Aircraft=fdmex->GetAircraft(); + Auxiliary=fdmex->GetAuxiliary(); } else { cout << "FGInitialCondition: This class requires a pointer to a valid FGFDMExec object" << endl; } @@ -173,7 +175,8 @@ void FGInitialCondition::SetMachIC(double mach) void FGInitialCondition::SetVcalibratedKtsIC(double vcas) { double altitudeASL = GetAltitudeASLFtIC(); - double mach = Atmosphere->MachFromVcalibrated(fabs(vcas)*ktstofps, altitudeASL); + double pressure = Atmosphere->GetPressure(altitudeASL); + double mach = Auxiliary->MachFromVcalibrated(fabs(vcas)*ktstofps, pressure); double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL); SetVtrueFpsIC(mach * soundSpeed); @@ -675,12 +678,13 @@ 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 = Atmosphere->VcalibratedFromMach(mach0, altitudeASL); + double vc0 = Auxiliary->VcalibratedFromMach(mach0, pressure); double ve0 = vt * sqrt(rho/rhoSL); switch(lastLatitudeSet) { @@ -718,10 +722,11 @@ void FGInitialCondition::SetAltitudeAGLFtIC(double agl) altitudeASL = GetAltitudeASLFtIC(); soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL); rho = Atmosphere->GetDensity(altitudeASL); + pressure = Atmosphere->GetPressure(altitudeASL); switch(lastSpeedSet) { case setvc: - mach0 = Atmosphere->MachFromVcalibrated(vc0, altitudeASL); + mach0 = Auxiliary->MachFromVcalibrated(vc0, pressure); SetVtrueFpsIC(mach0 * soundSpeed); break; case setmach: @@ -745,12 +750,13 @@ 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 = Atmosphere->VcalibratedFromMach(mach0, altitudeASL); + double vc0 = Auxiliary->VcalibratedFromMach(mach0, pressure); double ve0 = vt * sqrt(rho/rhoSL); switch(lastLatitudeSet) { @@ -820,10 +826,11 @@ 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 = Atmosphere->MachFromVcalibrated(vc0, altitudeASL); + mach0 = Auxiliary->MachFromVcalibrated(vc0, pressure); SetVtrueFpsIC(mach0 * soundSpeed); break; case setmach: @@ -953,10 +960,11 @@ 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 * Atmosphere->VcalibratedFromMach(mach, altitudeASL); + return fpstokts * Auxiliary->VcalibratedFromMach(mach, pressure); } //****************************************************************************** diff --git a/src/initialization/FGInitialCondition.h b/src/initialization/FGInitialCondition.h index 6eb0c7fd14..860419548f 100644 --- a/src/initialization/FGInitialCondition.h +++ b/src/initialization/FGInitialCondition.h @@ -62,6 +62,7 @@ class FGMatrix33; class FGColumnVector3; class FGAtmosphere; class FGAircraft; +class FGAuxiliary; class FGPropertyManager; class Element; @@ -705,6 +706,7 @@ class JSBSIM_API FGInitialCondition : public FGJSBBase FGFDMExec *fdmex; FGAtmosphere* Atmosphere; FGAircraft* Aircraft; + FGAuxiliary* Auxiliary; bool Load_v1(Element* document); bool Load_v2(Element* document); diff --git a/src/models/FGAtmosphere.cpp b/src/models/FGAtmosphere.cpp index 2c7b29c699..b5aa81cd8d 100644 --- a/src/models/FGAtmosphere.cpp +++ b/src/models/FGAtmosphere.cpp @@ -309,79 +309,6 @@ 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 f5e545d7f7..95d54e486b 100644 --- a/src/models/FGAtmosphere.h +++ b/src/models/FGAtmosphere.h @@ -206,41 +206,6 @@ 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; @@ -248,6 +213,7 @@ class JSBSIM_API FGAtmosphere : public FGModel { static constexpr double StdDaySLtemperature = 518.67; static constexpr double StdDaySLpressure = 2116.228; const double StdDaySLsoundspeed; + static constexpr double SHRatio = 1.4; protected: // Sea level conditions @@ -323,8 +289,6 @@ class JSBSIM_API FGAtmosphere : public FGModel { static constexpr double Reng0 = Rstar / Mair; //@} - static constexpr double SHRatio = 1.4; - double Reng = Reng0; virtual void bind(void); diff --git a/src/models/FGAuxiliary.cpp b/src/models/FGAuxiliary.cpp index 9dbdc1bc1c..578341a4a6 100644 --- a/src/models/FGAuxiliary.cpp +++ b/src/models/FGAuxiliary.cpp @@ -61,9 +61,9 @@ CLASS IMPLEMENTATION FGAuxiliary::FGAuxiliary(FGFDMExec* fdmex) : FGModel(fdmex) { Name = "FGAuxiliary"; - pt = 2116.23; // ISA SL pressure - tatc = 15.0; // ISA SL temperature - tat = 518.67; + pt = FGAtmosphere::StdDaySLpressure; // ISA SL pressure + tat = FGAtmosphere::StdDaySLtemperature; // ISA SL temperature + tatc = RankineToCelsius(tat); vcas = veas = 0.0; qbar = qbarUW = qbarUV = 0.0; @@ -194,10 +194,10 @@ bool FGAuxiliary::Run(bool Holding) tat = in.Temperature*(1 + 0.2*Mach*Mach); // Total Temperature, isentropic flow tatc = RankineToCelsius(tat); - pt = FDMExec->GetAtmosphere()->PitotTotalPressure(Mach, in.Pressure); + pt = PitotTotalPressure(Mach, in.Pressure); if (abs(Mach) > 0.0) { - vcas = FDMExec->GetAtmosphere()->VcalibratedFromMach(Mach, in.AltitudeASL); + vcas = VcalibratedFromMach(Mach, in.Pressure); veas = sqrt(2 * qbar / in.DensitySL); } else @@ -229,6 +229,82 @@ bool FGAuxiliary::Run(bool Holding) return false; } +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +double FGAuxiliary::PitotTotalPressure(double mach, double pressure) const +{ + constexpr double SHRatio = FGAtmosphere::SHRatio; + 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 pressure; + if (mach < 1) //calculate total pressure assuming isentropic flow + return pressure*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 pressure*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 FGAuxiliary::MachFromImpactPressure(double qc, double pressure) const +{ + constexpr double SHRatio = FGAtmosphere::SHRatio; + 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 / pressure + 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 FGAuxiliary::VcalibratedFromMach(double mach, double pressure) const +{ + double qc = PitotTotalPressure(mach, pressure) - pressure; + return in.StdDaySLsoundspeed * MachFromImpactPressure(qc, FGAtmosphere::StdDaySLpressure); +} + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +double FGAuxiliary::MachFromVcalibrated(double vcas, double pressure) const +{ + constexpr double StdDaySLpressure = FGAtmosphere::StdDaySLpressure; + double qc = PitotTotalPressure(vcas / in.StdDaySLsoundspeed, StdDaySLpressure) - StdDaySLpressure; + return MachFromImpactPressure(qc, pressure); +} + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // // From Stevens and Lewis, "Aircraft Control and Simulation", 3rd Ed., the diff --git a/src/models/FGAuxiliary.h b/src/models/FGAuxiliary.h index 4a4f12c3cb..fc0bba4153 100644 --- a/src/models/FGAuxiliary.h +++ b/src/models/FGAuxiliary.h @@ -120,6 +120,41 @@ class JSBSIM_API FGAuxiliary : public FGModel { // GET functions + /** 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 pressure Pressure in psf + * @return The total pressure in front of the Pitot tube in psf */ + double PitotTotalPressure(double mach, double pressure) 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 pressure 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 pressure Pressure in psf + * @return The calibrated airspeed (CAS) in ft/s + * */ + double VcalibratedFromMach(double mach, double pressure) 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 pressure Pressure in psf + * @return The Mach number + * */ + double MachFromVcalibrated(double vcas, double pressure) const; + // Atmospheric parameters GET functions /** Returns Calibrated airspeed in feet/second.*/ double GetVcalibratedFPS(void) const { return vcas; } @@ -257,10 +292,10 @@ class JSBSIM_API FGAuxiliary : public FGModel { double DensitySL; double PressureSL; double Temperature; + double StdDaySLsoundspeed; double SoundSpeed; double KinematicViscosity; double DistanceAGL; - double AltitudeASL; double Wingspan; double Wingchord; double StandardGravity; diff --git a/tests/unit_tests/CMakeLists.txt b/tests/unit_tests/CMakeLists.txt index 2d9fbb9d97..626986b078 100644 --- a/tests/unit_tests/CMakeLists.txt +++ b/tests/unit_tests/CMakeLists.txt @@ -18,7 +18,8 @@ set(UNIT_TESTS FGColumnVector3Test FGParameterValueTest FGConditionTest FGPropertyManagerTest - FGAtmosphereTest) + FGAtmosphereTest + FGAuxiliaryTest) foreach(test ${UNIT_TESTS}) cxxtest_add_test(${test}1 ${test}.cpp ${CMAKE_CURRENT_SOURCE_DIR}/${test}.h) diff --git a/tests/unit_tests/FGAtmosphereTest.h b/tests/unit_tests/FGAtmosphereTest.h index b862d6dd9a..cd33973665 100644 --- a/tests/unit_tests/FGAtmosphereTest.h +++ b/tests/unit_tests/FGAtmosphereTest.h @@ -613,152 +613,4 @@ 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/FGAuxiliaryTest.h b/tests/unit_tests/FGAuxiliaryTest.h new file mode 100644 index 0000000000..1d4547ba70 --- /dev/null +++ b/tests/unit_tests/FGAuxiliaryTest.h @@ -0,0 +1,180 @@ +#include +#include + +#include +#include +#include + +const double epsilon = 100. * std::numeric_limits::epsilon(); + +using namespace JSBSim; + +class DummyAtmosphere : public FGAtmosphere +{ +public: + DummyAtmosphere(FGFDMExec* fdm) : FGAtmosphere(fdm) {} + + // Getters for the protected members + static constexpr double GetR(void) { return Reng0; } +}; + +constexpr double R = DummyAtmosphere::GetR(); + +class FGAtmosphereTest : public CxxTest::TestSuite +{ +public: + static constexpr double gama = FGAtmosphere::SHRatio; + + void testPitotTotalPressure() { + FGFDMExec fdmex; + auto atm = fdmex.GetAtmosphere(); + auto aux = FGAuxiliary(&fdmex); + 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 = aux.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 = fdmex.GetAtmosphere(); + auto aux = FGAuxiliary(&fdmex); + 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 = aux.MachFromImpactPressure(P2-p1, p1); + double a2 = sqrt(gama*R*t2); + double M2 = u2/a2; + double mach2 = aux.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 = fdmex.GetAtmosphere(); + auto aux = FGAuxiliary(&fdmex); + atm->InitModel(); + aux.in.StdDaySLsoundspeed = atm->StdDaySLsoundspeed; + + // Ambient conditions far upstream (i.e. upstream the normal schock + // in supersonic flight) + double t1 = atm->GetTemperatureSL(); + double p1 = atm->GetPressureSL(); + + TS_ASSERT_DELTA(aux.VcalibratedFromMach(0.0, p1), 0.0, epsilon); + TS_ASSERT_DELTA(aux.MachFromVcalibrated(0.0, p1), 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(aux.VcalibratedFromMach(M1, p1)/u1, 1.0, 1e-7); + TS_ASSERT_DELTA(aux.MachFromVcalibrated(u1, p1)/M1, 1.0, 1e-7); + } + + // Check the VCAS computation at an altitude of 1000 ft + double asl = atm->GetSoundSpeedSL(); + p1 = atm->GetPressure(1000.); + + TS_ASSERT_DELTA(aux.VcalibratedFromMach(0.0, p1), 0.0, epsilon); + TS_ASSERT_DELTA(aux.MachFromVcalibrated(0.0, p1), 0.0, epsilon); + + for(double M=0.1; M<3.0; M+=0.25) { + double vcas = M*asl; + double M1 = aux.MachFromVcalibrated(vcas, p1); + TS_ASSERT_DELTA(aux.VcalibratedFromMach(M1, p1)/vcas, 1.0, 1e-7); + } + + double psl = atm->GetPressureSL(); + 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 = aux.MachFromImpactPressure(P2-p1, psl); + + TS_ASSERT_DELTA(aux.VcalibratedFromMach(M1, p1)/(mach*asl), 1.0, 1e-8); + } + } +};