X-Git-Url: https://git.mxchange.org/?a=blobdiff_plain;f=src%2FFDM%2FJSBSim%2Fmodels%2FFGAtmosphere.cpp;h=12fb084a70f879c384bcf25f3962d8d0f4f37a9f;hb=0917a5e062b531963f9f3d16bb0f95f769d34f61;hp=8c47b022676420c76c3ae9a6f3ebf23aefed24ed;hpb=27a730573637faa26c38c82d4b4b2a09f1079684;p=flightgear.git diff --git a/src/FDM/JSBSim/models/FGAtmosphere.cpp b/src/FDM/JSBSim/models/FGAtmosphere.cpp index 8c47b0226..12fb084a7 100644 --- a/src/FDM/JSBSim/models/FGAtmosphere.cpp +++ b/src/FDM/JSBSim/models/FGAtmosphere.cpp @@ -7,7 +7,7 @@ Purpose: Models the atmosphere Called by: FGSimExec - ------------- Copyright (C) 1999 Jon S. Berndt (jsb@hal-pc.org) ------------- + ------------- Copyright (C) 1999 Jon S. Berndt (jon@jsbsim.org) ------------- This program is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software @@ -48,16 +48,20 @@ INCLUDES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ #include "FGAtmosphere.h" -#include -#include #include "FGAircraft.h" #include "FGPropagate.h" #include "FGInertial.h" -#include +#include "FGAuxiliary.h" +#include "FGFDMExec.h" +#include "input_output/FGPropertyManager.h" +#include +#include + +using namespace std; namespace JSBSim { -static const char *IdSrc = "$Id$"; +static const char *IdSrc = "$Id: FGAtmosphere.cpp,v 1.38 2010/09/16 11:01:24 jberndt Exp $"; static const char *IdHdr = ID_ATMOSPHERE; /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -82,8 +86,8 @@ FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex) MagnitudedAccelDt = MagnitudeAccel = Magnitude = 0.0; // SetTurbType( ttCulp ); SetTurbType( ttNone ); - TurbGain = 0.0; - TurbRate = 1.7; + TurbGain = 1.0; + TurbRate = 10.0; Rhythmicity = 0.1; spike = target_time = strength = 0.0; wind_from_clockwise = 0.0; @@ -96,6 +100,22 @@ FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex) vGustNED.InitMatrix(); vTurbulenceNED.InitMatrix(); + // Milspec turbulence model + windspeed_at_20ft = 0.; + probability_of_exceedence_index = 0; + POE_Table = new FGTable(7,12); + // this is Figure 7 from p. 49 of MIL-F-8785C + // rows: probability of exceedance curve index, cols: altitude in ft + *POE_Table + << 500.0 << 1750.0 << 3750.0 << 7500.0 << 15000.0 << 25000.0 << 35000.0 << 45000.0 << 55000.0 << 65000.0 << 75000.0 << 80000.0 + << 1 << 3.2 << 2.2 << 1.5 << 0.0 << 0.0 << 0.0 << 0.0 << 0.0 << 0.0 << 0.0 << 0.0 << 0.0 + << 2 << 4.2 << 3.6 << 3.3 << 1.6 << 0.0 << 0.0 << 0.0 << 0.0 << 0.0 << 0.0 << 0.0 << 0.0 + << 3 << 6.6 << 6.9 << 7.4 << 6.7 << 4.6 << 2.7 << 0.4 << 0.0 << 0.0 << 0.0 << 0.0 << 0.0 + << 4 << 8.6 << 9.6 << 10.6 << 10.1 << 8.0 << 6.6 << 5.0 << 4.2 << 2.7 << 0.0 << 0.0 << 0.0 + << 5 << 11.8 << 13.0 << 16.0 << 15.1 << 11.6 << 9.7 << 8.1 << 8.2 << 7.9 << 4.9 << 3.2 << 2.1 + << 6 << 15.6 << 17.6 << 23.0 << 23.6 << 22.1 << 20.0 << 16.0 << 15.1 << 12.1 << 7.9 << 6.2 << 5.1 + << 7 << 18.7 << 21.5 << 28.4 << 30.2 << 30.7 << 31.0 << 25.2 << 23.1 << 17.5 << 10.7 << 8.4 << 7.2; + bind(); Debug(0); } @@ -135,8 +155,10 @@ bool FGAtmosphere::Run(void) if (FGModel::Run()) return true; if (FDMExec->Holding()) return false; + RunPreFunctions(); + T_dev = 0.0; - h = Propagate->Geth(); + h = Propagate->GetAltitudeASL(); if (!useExternal) { Calculate(h); @@ -145,6 +167,8 @@ bool FGAtmosphere::Run(void) CalculateDerived(); } + RunPostFunctions(); + Debug(2); return false; } @@ -269,7 +293,9 @@ void FGAtmosphere::Calculate(double altitude) void FGAtmosphere::CalculateDerived(void) { T_dev = (*temperature) - GetTemperature(h); - density_altitude = h + T_dev * 66.7; + + if (T_dev == 0.0) density_altitude = h; + else density_altitude = 518.67/0.00356616 * (1.0 - pow(GetDensityRatio(),0.235)); if (turbType != ttNone) Turbulence(); @@ -337,6 +363,9 @@ static inline double square_signed (double value) return value * value; } +/// simply square a value +static inline double sqr(double x) { return x*x; } + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // // psi is the angle that the wind is blowing *towards* @@ -375,9 +404,11 @@ void FGAtmosphere::SetWindPsi(double dir) void FGAtmosphere::Turbulence(void) { + double DeltaT = rate*FDMExec->GetDeltaT(); + switch (turbType) { case ttStandard: { - TurbGain = TurbGain * TurbGain * 100.0; + // TurbGain = TurbGain * TurbGain * 100.0; // what is this!? vDirectiondAccelDt(eX) = 1 - 2.0*(double(rand())/double(RAND_MAX)); vDirectiondAccelDt(eY) = 1 - 2.0*(double(rand())/double(RAND_MAX)); @@ -388,8 +419,8 @@ void FGAtmosphere::Turbulence(void) // away from the peaks MagnitudedAccelDt = ((MagnitudedAccelDt - Magnitude) / (1 + fabs(Magnitude))); - MagnitudeAccel += MagnitudedAccelDt*rate*TurbRate*State->Getdt(); - Magnitude += MagnitudeAccel*rate*State->Getdt(); + MagnitudeAccel += MagnitudedAccelDt*TurbRate*DeltaT; + Magnitude += MagnitudeAccel*DeltaT; Magnitude = fabs(Magnitude); vDirectiondAccelDt.Normalize(); @@ -398,9 +429,9 @@ void FGAtmosphere::Turbulence(void) vDirectiondAccelDt(eX) = square_signed(vDirectiondAccelDt(eX)); vDirectiondAccelDt(eY) = square_signed(vDirectiondAccelDt(eY)); - vDirectionAccel += vDirectiondAccelDt*rate*TurbRate*State->Getdt(); + vDirectionAccel += vDirectiondAccelDt*TurbRate*DeltaT; vDirectionAccel.Normalize(); - vDirection += vDirectionAccel*rate*State->Getdt(); + vDirection += vDirectionAccel*DeltaT; vDirection.Normalize(); @@ -449,42 +480,39 @@ void FGAtmosphere::Turbulence(void) } case ttBerndt: { // This is very experimental and incomplete at the moment. - TurbGain = TurbGain * TurbGain * 100.0; - - vDirectiondAccelDt(eX) = 1 - 2.0*(double(rand())/double(RAND_MAX)); - vDirectiondAccelDt(eY) = 1 - 2.0*(double(rand())/double(RAND_MAX)); - vDirectiondAccelDt(eZ) = 1 - 2.0*(double(rand())/double(RAND_MAX)); - - - MagnitudedAccelDt = 1 - 2.0*(double(rand())/double(RAND_MAX)) - Magnitude; - MagnitudeAccel += MagnitudedAccelDt*rate*State->Getdt(); - Magnitude += MagnitudeAccel*rate*State->Getdt(); + vDirectiondAccelDt(eX) = GaussianRandomNumber(); + vDirectiondAccelDt(eY) = GaussianRandomNumber(); + vDirectiondAccelDt(eZ) = GaussianRandomNumber(); +/* + MagnitudedAccelDt = GaussianRandomNumber(); + MagnitudeAccel += MagnitudedAccelDt * DeltaT; + Magnitude += MagnitudeAccel * DeltaT; +*/ + Magnitude += GaussianRandomNumber() * DeltaT; vDirectiondAccelDt.Normalize(); - vDirectionAccel += vDirectiondAccelDt*rate*State->Getdt(); + vDirectionAccel += TurbRate * vDirectiondAccelDt * DeltaT; vDirectionAccel.Normalize(); - vDirection += vDirectionAccel*rate*State->Getdt(); + vDirection += vDirectionAccel*DeltaT; - // Diminish z-vector within two wingspans - // of the ground + // Diminish z-vector within two wingspans of the ground double HOverBMAC = Auxiliary->GetHOverBMAC(); - if (HOverBMAC < 2.0) - vDirection(eZ) *= HOverBMAC / 2.0; + if (HOverBMAC < 2.0) vDirection(eZ) *= HOverBMAC / 2.0; vDirection.Normalize(); vTurbulenceNED = TurbGain*Magnitude * vDirection; vTurbulenceGrad = TurbGain*MagnitudeAccel * vDirection; - vBodyTurbGrad = Propagate->GetTl2b()*vTurbulenceGrad; - vTurbPQR(eP) = vBodyTurbGrad(eY)/Aircraft->GetWingSpan(); + vBodyTurbGrad = Propagate->GetTl2b() * vTurbulenceGrad; + vTurbPQR(eP) = vBodyTurbGrad(eY) / Aircraft->GetWingSpan(); if (Aircraft->GetHTailArm() > 0) - vTurbPQR(eQ) = vBodyTurbGrad(eZ)/Aircraft->GetHTailArm(); + vTurbPQR(eQ) = vBodyTurbGrad(eZ) / Aircraft->GetHTailArm(); else - vTurbPQR(eQ) = vBodyTurbGrad(eZ)/10.0; + vTurbPQR(eQ) = vBodyTurbGrad(eZ) / 10.0; if (Aircraft->GetVTailArm() > 0) - vTurbPQR(eR) = vBodyTurbGrad(eX)/Aircraft->GetVTailArm(); + vTurbPQR(eR) = vBodyTurbGrad(eX) / Aircraft->GetVTailArm(); else vTurbPQR(eR) = vBodyTurbGrad(eX)/10.0; @@ -540,6 +568,134 @@ void FGAtmosphere::Turbulence(void) spike = spike * 0.9; break; } + case ttMilspec: + case ttTustin: { + // an index of zero means turbulence is disabled + if (probability_of_exceedence_index == 0) { + vTurbulenceNED(1) = vTurbulenceNED(2) = vTurbulenceNED(3) = 0.0; + vTurbPQR(1) = vTurbPQR(2) = vTurbPQR(3) = 0.0; + return; + } + + // Turbulence model according to MIL-F-8785C (Flying Qualities of Piloted Aircraft) + double + h = Propagate->GetDistanceAGL(), + V = Auxiliary->GetVt(), // true airspeed in ft/s + b_w = Aircraft->GetWingSpan(), + L_u, L_w, sig_u, sig_w; + + // clip height functions at 10 ft + if (h <= 10.) h = 10; + + // Scale lengths L and amplitudes sigma as function of height + if (h <= 1000) { + L_u = h/pow(0.177 + 0.000823*h, 1.2); // MIL-F-8785c, Fig. 10, p. 55 + L_w = h; + sig_w = 0.1*windspeed_at_20ft; + sig_u = sig_w/pow(0.177 + 0.000823*h, 0.4); // MIL-F-8785c, Fig. 11, p. 56 + } else if (h <= 2000) { + // linear interpolation between low altitude and high altitude models + L_u = L_w = 1000 + (h-1000.)/1000.*750.; + sig_u = sig_w = 0.1*windspeed_at_20ft + + (h-1000.)/1000.*(POE_Table->GetValue(probability_of_exceedence_index, h) - 0.1*windspeed_at_20ft); + } else { + L_u = L_w = 1750.; // MIL-F-8785c, Sec. 3.7.2.1, p. 48 + sig_u = sig_w = POE_Table->GetValue(probability_of_exceedence_index, h); + } + + // keep values from last timesteps + // TODO maybe use deque? + static double + xi_u_km1 = 0, nu_u_km1 = 0, + xi_v_km1 = 0, xi_v_km2 = 0, nu_v_km1 = 0, nu_v_km2 = 0, + xi_w_km1 = 0, xi_w_km2 = 0, nu_w_km1 = 0, nu_w_km2 = 0, + xi_p_km1 = 0, nu_p_km1 = 0, + xi_q_km1 = 0, xi_r_km1 = 0; + + + double + T_V = DeltaT, // for compatibility of nomenclature + sig_p = 1.9/sqrt(L_w*b_w)*sig_w, // Yeager1998, eq. (8) + sig_q = sqrt(M_PI/2/L_w/b_w), // eq. (14) + sig_r = sqrt(2*M_PI/3/L_w/b_w), // eq. (17) + L_p = sqrt(L_w*b_w)/2.6, // eq. (10) + tau_u = L_u/V, // eq. (6) + tau_w = L_w/V, // eq. (3) + tau_p = L_p/V, // eq. (9) + tau_q = 4*b_w/M_PI/V, // eq. (13) + tau_r =3*b_w/M_PI/V, // eq. (17) + nu_u = GaussianRandomNumber(), + nu_v = GaussianRandomNumber(), + nu_w = GaussianRandomNumber(), + nu_p = GaussianRandomNumber(), + xi_u, xi_v, xi_w, xi_p, xi_q, xi_r; + + // values of turbulence NED velocities + + if (turbType == ttTustin) { + // the following is the Tustin formulation of Yeager's report + double + omega_w = V/L_w, // hidden in nomenclature p. 3 + omega_v = V/L_u, // this is defined nowhere + C_BL = 1/tau_u/tan(T_V/2/tau_u), // eq. (19) + C_BLp = 1/tau_p/tan(T_V/2/tau_p), // eq. (22) + C_BLq = 1/tau_q/tan(T_V/2/tau_q), // eq. (24) + C_BLr = 1/tau_r/tan(T_V/2/tau_r); // eq. (26) + + xi_u = -(1 - C_BL*tau_u)/(1 + C_BL*tau_u)*xi_u_km1 + + sig_u*sqrt(2*tau_u/T_V)/(1 + C_BL*tau_u)*(nu_u + nu_u_km1); // eq. (18) + xi_v = -2*(sqr(omega_v) - sqr(C_BL))/sqr(omega_v + C_BL)*xi_v_km1 + - sqr(omega_v - C_BL)/sqr(omega_v + C_BL) * xi_v_km2 + + sig_u*sqrt(3*omega_v/T_V)/sqr(omega_v + C_BL)*( + (C_BL + omega_v/sqrt(3.))*nu_v + + 2/sqrt(3.)*omega_v*nu_v_km1 + + (omega_v/sqrt(3.) - C_BL)*nu_v_km2); // eq. (20) for v + xi_w = -2*(sqr(omega_w) - sqr(C_BL))/sqr(omega_w + C_BL)*xi_w_km1 + - sqr(omega_w - C_BL)/sqr(omega_w + C_BL) * xi_w_km2 + + sig_w*sqrt(3*omega_w/T_V)/sqr(omega_w + C_BL)*( + (C_BL + omega_w/sqrt(3.))*nu_w + + 2/sqrt(3.)*omega_w*nu_w_km1 + + (omega_w/sqrt(3.) - C_BL)*nu_w_km2); // eq. (20) for w + xi_p = -(1 - C_BLp*tau_p)/(1 + C_BLp*tau_p)*xi_p_km1 + + sig_p*sqrt(2*tau_p/T_V)/(1 + C_BLp*tau_p) * (nu_p + nu_p_km1); // eq. (21) + xi_q = -(1 - 4*b_w*C_BLq/M_PI/V)/(1 + 4*b_w*C_BLq/M_PI/V) * xi_q_km1 + + C_BLq/V/(1 + 4*b_w*C_BLq/M_PI/V) * (xi_w - xi_w_km1); // eq. (23) + xi_r = - (1 - 3*b_w*C_BLr/M_PI/V)/(1 + 3*b_w*C_BLr/M_PI/V) * xi_r_km1 + + C_BLr/V/(1 + 3*b_w*C_BLr/M_PI/V) * (xi_v - xi_v_km1); // eq. (25) + + } else if (turbType == ttMilspec) { + // the following is the MIL-STD-1797A formulation + // as cited in Yeager's report + xi_u = (1 - T_V/tau_u) *xi_u_km1 + sig_u*sqrt(2*T_V/tau_u)*nu_u; // eq. (30) + xi_v = (1 - 2*T_V/tau_u)*xi_v_km1 + sig_u*sqrt(4*T_V/tau_u)*nu_v; // eq. (31) + xi_w = (1 - 2*T_V/tau_w)*xi_w_km1 + sig_w*sqrt(4*T_V/tau_w)*nu_w; // eq. (32) + xi_p = (1 - T_V/tau_p) *xi_p_km1 + sig_p*sqrt(2*T_V/tau_p)*nu_p; // eq. (33) + xi_q = (1 - T_V/tau_q) *xi_q_km1 + M_PI/4/b_w*(xi_w - xi_w_km1); // eq. (34) + xi_r = (1 - T_V/tau_r) *xi_r_km1 + M_PI/3/b_w*(xi_v - xi_v_km1); // eq. (35) + } + + // rotate by wind azimuth and assign the velocities + double cospsi = cos(psiw), sinpsi = sin(psiw); + vTurbulenceNED(1) = cospsi*xi_u + sinpsi*xi_v; + vTurbulenceNED(2) = -sinpsi*xi_u + cospsi*xi_v; + vTurbulenceNED(3) = xi_w; + + vTurbPQR(1) = cospsi*xi_p + sinpsi*xi_q; + vTurbPQR(2) = -sinpsi*xi_p + cospsi*xi_q; + vTurbPQR(3) = xi_r; + + // vTurbPQR is in the body fixed frame, not NED + vTurbPQR = Propagate->GetTl2b()*vTurbPQR; + + // hand on the values for the next timestep + xi_u_km1 = xi_u; nu_u_km1 = nu_u; + xi_v_km2 = xi_v_km1; xi_v_km1 = xi_v; nu_v_km2 = nu_v_km1; nu_v_km1 = nu_v; + xi_w_km2 = xi_w_km1; xi_w_km1 = xi_w; nu_w_km2 = nu_w_km1; nu_w_km1 = nu_w; + xi_p_km1 = xi_p; nu_p_km1 = nu_p; + xi_q_km1 = xi_q; + xi_r_km1 = xi_r; + + } default: break; } @@ -571,7 +727,9 @@ void FGAtmosphere::bind(void) { typedef double (FGAtmosphere::*PMF)(int) const; typedef double (FGAtmosphere::*PMFv)(void) const; + typedef int (FGAtmosphere::*PMFt)(void) const; typedef void (FGAtmosphere::*PMFd)(int,double); + typedef void (FGAtmosphere::*PMFi)(int); PropertyManager->Tie("atmosphere/T-R", this, (PMFv)&FGAtmosphere::GetTemperature); PropertyManager->Tie("atmosphere/rho-slugs_ft3", this, (PMFv)&FGAtmosphere::GetDensity); PropertyManager->Tie("atmosphere/P-psf", this, (PMFv)&FGAtmosphere::GetPressure); @@ -608,13 +766,29 @@ void FGAtmosphere::bind(void) PropertyManager->Tie("atmosphere/gust-down-fps", this, eDown, (PMF)&FGAtmosphere::GetGustNED, (PMFd)&FGAtmosphere::SetGustNED); + PropertyManager->Tie("atmosphere/turb-north-fps", this, eNorth, (PMF)&FGAtmosphere::GetTurbNED, + (PMFd)&FGAtmosphere::SetTurbNED); + PropertyManager->Tie("atmosphere/turb-east-fps", this, eEast, (PMF)&FGAtmosphere::GetTurbNED, + (PMFd)&FGAtmosphere::SetTurbNED); + PropertyManager->Tie("atmosphere/turb-down-fps", this, eDown, (PMF)&FGAtmosphere::GetTurbNED, + (PMFd)&FGAtmosphere::SetTurbNED); + PropertyManager->Tie("atmosphere/p-turb-rad_sec", this,1, (PMF)&FGAtmosphere::GetTurbPQR); PropertyManager->Tie("atmosphere/q-turb-rad_sec", this,2, (PMF)&FGAtmosphere::GetTurbPQR); PropertyManager->Tie("atmosphere/r-turb-rad_sec", this,3, (PMF)&FGAtmosphere::GetTurbPQR); + PropertyManager->Tie("atmosphere/turb-type", this, (PMFt)&FGAtmosphere::GetTurbType, (PMFi)&FGAtmosphere::SetTurbType); PropertyManager->Tie("atmosphere/turb-rate", this, &FGAtmosphere::GetTurbRate, &FGAtmosphere::SetTurbRate); PropertyManager->Tie("atmosphere/turb-gain", this, &FGAtmosphere::GetTurbGain, &FGAtmosphere::SetTurbGain); PropertyManager->Tie("atmosphere/turb-rhythmicity", this, &FGAtmosphere::GetRhythmicity, &FGAtmosphere::SetRhythmicity); + + PropertyManager->Tie("atmosphere/turbulence/milspec/windspeed_at_20ft_AGL-fps", + this, &FGAtmosphere::GetWindspeed20ft, + &FGAtmosphere::SetWindspeed20ft); + PropertyManager->Tie("atmosphere/turbulence/milspec/severity", + this, &FGAtmosphere::GetProbabilityOfExceedence, + &FGAtmosphere::SetProbabilityOfExceedence); + } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%