src/FDM/YASim/Jet.cpp

   1 #include "Atmosphere.hpp"
   2 #include "Math.hpp"
   3 #include "Jet.hpp"
   4 namespace yasim {
   5
   6 Jet::Jet()
   7 {
   8     _maxThrust = 0;
   9     _abFactor = 1;
  10     _reheat = 0;
  11     _rotControl = 0;
  12     _maxRot = 0;
  13
  14     // Initialize parameters for an early-ish subsonic turbojet.  More
  15     // recent turbofans will typically have a lower vMax, epr0, and
  16     // tsfc.
  17     _vMax = 800;
  18     _epr0 = 3.0;
  19     _tsfc = 0.8f;
  20     _egt0 = 1050;
  21     _n1Min = 55;
  22     _n1Max = 102;
  23     _n2Min = 73;
  24     _n2Max = 103;
  25     setSpooling(4); // 4 second spool time? s'bout right.
  26
  27     // And initialize to an engine that is idling
  28     _n1 = _n1Min;
  29     _n2 = _n2Min;
  30
  31     // And sanify the remaining junk, just in case.
  32     _running = true;
  33     _cranking = false;
  34     _epr = 1;
  35     _fuelFlow = 0;
  36     _egt = 273;
  37     _tempCorrect = 1;
  38     _pressureCorrect = 1;
  39 }
  40
  41 void Jet::stabilize()
  42 {
  43     // Just run it for an hour, there's no need to iterate given the
  44     // algorithms used.
  45     integrate(3600);
  46 }
  47
  48 void Jet::setMaxThrust(float thrust, float afterburner)
  49 {
  50     _maxThrust = thrust;
  51     if(afterburner == 0) _abFactor = 1;
  52     else                 _abFactor = afterburner/thrust;
  53 }
  54
  55 void Jet::setVMax(float spd)
  56 {
  57     _vMax = spd;
  58 }
  59
  60 void Jet::setTSFC(float tsfc)
  61 {
  62     _tsfc = tsfc;
  63 }
  64
  65 void Jet::setRPMs(float idleN1, float maxN1, float idleN2, float maxN2)
  66 {
  67     _n1Min = idleN1;
  68     _n1Max = maxN1;
  69     _n2Min = idleN2;
  70     _n2Max = maxN2;
  71 }
  72
  73 void Jet::setEGT(float takeoffEGT)
  74 {
  75     _egt0 = takeoffEGT;
  76 }
  77
  78 void Jet::setEPR(float takeoffEPR)
  79 {
  80     _epr0 = takeoffEPR;
  81 }
  82
  83 void Jet::setSpooling(float time)
  84 {
  85     // 2.3 = -ln(0.1), i.e. x=2.3 is the 90% point we're defining
  86     // The extra fudge factor is there because the N1 speed (which
  87     // determines thrust) lags the N2 speed.
  88     _decay = 1.5f * 2.3f / time;
  89 }
  90
  91 void Jet::setVectorAngle(float angle)
  92 {
  93     _maxRot = angle;
  94 }
  95
  96 void Jet::setReheat(float reheat)
  97 {
  98     _reheat = Math::clamp(reheat, 0, 1);
  99 }
 100
 101 void Jet::setRotation(float rot)
 102 {
 103     if(rot < 0) rot = 0;
 104     if(rot > 1) rot = 1;
 105     _rotControl = rot;
 106 }
 107
 108 float Jet::getN1()
 109 {
 110     return _n1 * _tempCorrect;
 111 }
 112
 113 float Jet::getN2()
 114 {
 115     return _n2 * _tempCorrect;
 116 }
 117
 118 float Jet::getEPR()
 119 {
 120     return _epr;
 121 }
 122
 123 float Jet::getEGT()
 124 {
 125     // Exactly zero means "off" -- return the ambient temperature
 126     if(_egt == 0) return _temp;
 127
 128     return _egt * _tempCorrect * _tempCorrect;
 129 }
 130
 131 float Jet::getFuelFlow()
 132 {
 133     return _fuelFlow * _pressureCorrect;
 134 }
 135
 136 void Jet::integrate(float dt)
 137 {
 138     // Sea-level values
 139     const static float P0 = Atmosphere::getStdPressure(0);
 140     const static float T0 = Atmosphere::getStdTemperature(0);
 141     const static float D0 = Atmosphere::getStdDensity(0);
 142
 143     float speed = -Math::dot3(_wind, _dir);
 144
 145     float statT, statP, statD;
 146     Atmosphere::calcStaticAir(_pressure, _temp, _rho, speed,
 147                               &statP, &statT, &statD);
 148     _pressureCorrect = statP/P0;
 149     _tempCorrect = Math::sqrt(statT/T0);
 150
 151     // Linearly taper maxThrust to zero at vMax
 152     float vCorr = 1 - (speed/_vMax);
 153
 154     float maxThrust = _maxThrust * vCorr * (statD/D0);
 155     _thrust = maxThrust * _throttle;
 156
 157     // Now get a "beta" (i.e. EPR - 1) value.  The output values are
 158     // expressed as functions of beta.
 159     float ibeta0 = 1/(_epr0 - 1);
 160     float betaTarget = (_epr0 - 1) * (_thrust/_maxThrust) * (P0/_pressure)
 161         * (_temp/statT);
 162     float n2Target = _n2Min + (betaTarget*ibeta0) * (_n2Max - _n2Min);
 163
 164     // Note that this "first" beta value is used to compute a target
 165     // for N2 only Integrate the N2 speed and back-calculate a beta1
 166     // target.  The N1 speed will seek to this.
 167     _n2 = (_n2 + dt*_decay * n2Target) / (1 + dt*_decay);
 168
 169     float betaN2 = (_epr0-1) * (_n2 - _n2Min) / (_n2Max - _n2Min);
 170     float n1Target = _n1Min + betaN2*ibeta0 * (_n1Max - _n1Min);
 171     _n1 = (_n1 + dt*_decay * n1Target) / (1 + dt*_decay);
 172
 173     // The actual thrust produced is keyed to the N1 speed.  Add the
 174     // afterburners in at the end.
 175     float betaN1 =  (_epr0-1) * (_n1 - _n1Min) / (_n1Max - _n1Min);
 176     _thrust *= betaN1/(betaTarget+0.00001f); // blowup protection
 177     _thrust *= 1 + _reheat*(_abFactor-1);
 178
 179     // Finally, calculate the output variables.   Use a 80/20 mix of
 180     // the N2/N1 speeds as the key.
 181     float beta = 0.8f*betaN2 + 0.2f*betaN1;
 182     _epr = beta + 1;
 183     float ff0 = _maxThrust*_tsfc*(1/(3600.0f*9.8f)); // takeoff fuel flow, kg/s
 184     _fuelFlow = ff0 * beta*ibeta0;
 185     _fuelFlow *= 1 + (3.5f * _reheat * _abFactor); // Afterburners take
 186                                                   // 3.5 times as much
 187                                                   // fuel per thrust unit
 188     _egt = T0 + beta*ibeta0 * (_egt0 - T0);
 189 }
 190
 191 bool Jet::isRunning()
 192 {
 193     return _running;
 194 }
 195
 196 bool Jet::isCranking()
 197 {
 198     return _cranking;
 199 }
 200
 201 void Jet::getThrust(float* out)
 202 {
 203     Math::mul3(_thrust, _dir, out);
 204
 205     // Rotate about the Y axis for thrust vectoring
 206     float angle = _rotControl * _maxRot;
 207     float s = Math::sin(angle);
 208     float c = Math::cos(angle);
 209     float o0 = out[0];
 210     out[0] =  c * o0 + s * out[2];
 211     out[2] = -s * o0 + c * out[2];
 212 }
 213
 214 void Jet::getTorque(float* out)
 215 {
 216     out[0] = out[1] = out[2] = 0;
 217     return;
 218 }
 219
 220 void Jet::getGyro(float* out)
 221 {
 222     out[0] = out[1] = out[2] = 0;
 223     return;
 224 }
 225
 226 }; // namespace yasim