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