]> git.mxchange.org Git - flightgear.git/blob - src/FDM/LaRCsim/IO360.cxx
Adds a basic FDM model for LaRCsim debugging purposes.
[flightgear.git] / src / FDM / LaRCsim / IO360.cxx
1 // IO360.cxx - a piston engine model currently for the IO360 engine fitted to the C172
2 //             but with the potential to model other naturally aspirated piston engines
3 //             given appropriate config input.
4 //
5 // Written by David Luff, started 2000.
6 // Based on code by Phil Schubert, started 1999.
7 //
8 // This program is free software; you can redistribute it and/or
9 // modify it under the terms of the GNU General Public License as
10 // published by the Free Software Foundation; either version 2 of the
11 // License, or (at your option) any later version.
12 //
13 // This program is distributed in the hope that it will be useful, but
14 // WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16 // General Public License for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with this program; if not, write to the Free Software
20 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
21
22 #include <simgear/compiler.h>
23
24 #include <math.h>
25
26 #include STL_FSTREAM
27 #include STL_IOSTREAM
28
29 #include <Main/fg_props.hxx>
30
31 #include "IO360.hxx"
32 #include "ls_constants.h"
33
34
35 //*************************************************************************************
36 // Initialise the engine model
37 void FGNewEngine::init(double dt) {
38
39     // These constants should probably be moved eventually
40     CONVERT_CUBIC_INCHES_TO_METERS_CUBED = 1.638706e-5;
41     CONVERT_HP_TO_WATTS = 745.6999;
42
43     // Properties of working fluids
44     Cp_air = 1005;      // J/KgK
45     Cp_fuel = 1700;     // J/KgK
46     calorific_value_fuel = 47.3e6;  // W/Kg  Note that this is only an approximate value
47     rho_fuel = 800;     // kg/m^3 - an estimate for now
48     R_air = 287.3;
49
50     // environment inputs
51     p_amb_sea_level = 101325;   // Pascals              
52
53     // Control inputs  - ARE THESE NEEDED HERE???
54     Throttle_Lever_Pos = 75;
55     Propeller_Lever_Pos = 75;
56     Mixture_Lever_Pos = 100;
57
58     //misc
59     IAS = 0;
60     time_step = dt;
61
62     // Engine Specific Variables that should be read in from a config file
63     MaxHP = 200;    //Lycoming IO360 -A-C-D series
64 //  MaxHP = 180;    //Current Lycoming IO360 ?
65 //  displacement = 520;  //Continental IO520-M
66     displacement = 360;  //Lycoming IO360
67     displacement_SI = displacement * CONVERT_CUBIC_INCHES_TO_METERS_CUBED;
68     engine_inertia = 0.2;  //kgm^2 - value taken from a popular family saloon car engine - need to find an aeroengine value !!!!!
69     prop_inertia = 0.05;  //kgm^2 - this value is a total guess - dcl
70     Max_Fuel_Flow = 130;  // Units??? Do we need this variable any more??
71
72     // Engine specific variables that maybe should be read in from config but are pretty generic and won't vary much for a naturally aspirated piston engine.
73     Max_Manifold_Pressure = 28.50;  //Inches Hg. An approximation - should be able to find it in the engine performance data
74     Min_Manifold_Pressure = 6.5;    //Inches Hg. This is a guess corresponding to approx 0.24 bar MAP (7 in Hg) - need to find some proper data for this
75     Max_RPM = 2700;
76     Min_RPM = 600;                  //Recommended idle from Continental data sheet
77     Mag_Derate_Percent = 5;
78     Gear_Ratio = 1;
79     n_R = 2;         // Number of crank revolutions per power cycle - 2 for a 4 stroke engine.
80
81     // Various bits of housekeeping describing the engines initial state.
82     running = false;
83     cranking = false;
84     crank_counter = false;
85     starter = false;
86
87     // Initialise Engine Variables used by this instance
88     if(running)
89         RPM = 600;
90     else
91         RPM = 0;
92     Percentage_Power = 0;
93     Manifold_Pressure = 29.96; // Inches
94     Fuel_Flow_gals_hr = 0;
95 //    Torque = 0;
96     Torque_SI = 0;
97     CHT = 298.0;                        //deg Kelvin
98     CHT_degF = (CHT * 1.8) - 459.67;    //deg Fahrenheit
99     Mixture = 14;
100     Oil_Pressure = 0;   // PSI
101     Oil_Temp = 85;      // Deg C
102     current_oil_temp = 298.0;   //deg Kelvin
103     /**** one of these is superfluous !!!!***/
104     HP = 0;
105     RPS = 0;
106     Torque_Imbalance = 0;
107
108     // Initialise Propellor Variables used by this instance
109     FGProp1_RPS = 0;
110     // Hardcode propellor for now - the following two should be read from config eventually
111     prop_diameter = 1.8;         // meters
112     blade_angle = 23.0;          // degrees
113 }
114
115 //*****************************************************************************
116 // update the engine model based on current control positions
117 void FGNewEngine::update() {
118
119 /*
120     // Hack for testing - should output every 5 seconds
121     static int count1 = 0;
122     if(count1 == 0) {
123 //      cout << "P_atmos = " << p_amb << "  T_atmos = " << T_amb << '\n';
124 //      cout << "Manifold pressure = " << Manifold_Pressure << "  True_Manifold_Pressure = " << True_Manifold_Pressure << '\n';
125 //      cout << "p_amb_sea_level = " << p_amb_sea_level << '\n';
126 //      cout << "equivalence_ratio = " << equivalence_ratio << '\n';
127 //      cout << "combustion_efficiency = " << combustion_efficiency << '\n';
128 //      cout << "AFR = " << 14.7 / equivalence_ratio << '\n';
129 //      cout << "Mixture lever = " << Mixture_Lever_Pos << '\n';
130 //      cout << "n = " << RPM << " rpm\n";
131 //      cout << "T_amb = " << T_amb << '\n';
132 //      cout << "running = " << running << '\n';
133 //      cout << "fuel = " << fgGetFloat("/consumables/fuel/tank[0]/level-gal_us") << '\n';
134 //      cout << "Percentage_Power = " << Percentage_Power << '\n';
135 //      cout << "current_oil_temp = " << current_oil_temp << '\n';
136 //      cout << "EGT = " << EGT << '\n';
137     }
138     count1++;
139     if(count1 == 100)
140         count1 = 0;
141 */
142
143     // Check parameters that may alter the operating state of the engine. 
144     // (spark, fuel, starter motor etc)
145
146     // Check for spark
147     bool Magneto_Left = false;
148     bool Magneto_Right = false;
149     // Magneto positions:
150     // 0 -> off
151     // 1 -> left only
152     // 2 -> right only
153     // 3 -> both
154     if(mag_pos != 0) {
155         spark = true;
156     } else {
157         spark = false;
158     }  // neglects battery voltage, master on switch, etc for now.
159     if((mag_pos == 1) || (mag_pos > 2)) 
160         Magneto_Left = true;
161     if(mag_pos > 1)
162         Magneto_Right = true;
163  
164     // crude check for fuel
165     if((fgGetFloat("/consumables/fuel/tank[0]/level-gal_us") > 0) || (fgGetFloat("/consumables/fuel/tank[1]/level-gal_us") > 0)) {
166         fuel = true;
167     } else {
168         fuel = false;
169     }  // Need to make this better, eg position of fuel selector switch.
170
171     // Check if we are turning the starter motor
172     if(cranking != starter) {
173         // This check saves .../cranking from getting updated every loop - they only update when changed.
174         cranking = starter;
175         crank_counter = 0;
176     }
177     // Note that although /engines/engine[0]/starter and /engines/engine[0]/cranking might appear to be duplication it is
178     // not since the starter may be engaged with the battery voltage too low for cranking to occur (or perhaps the master 
179     // switch just left off) and the sound manager will read .../cranking to determine wether to play a cranking sound.
180     // For now though none of that is implemented so cranking can be set equal to .../starter without further checks.
181
182 //    int Alternate_Air_Pos =0; // Off = 0. Reduces power by 3 % for same throttle setting
183     // DCL - don't know what this Alternate_Air_Pos is - this is a leftover from the Schubert code.
184
185     //Check mode of engine operation
186     if(cranking) {
187         crank_counter++;
188         if(RPM <= 480) {
189             RPM += 100;
190             if(RPM > 480)
191                 RPM = 480;
192         } else {
193             // consider making a horrible noise if the starter is engaged with the engine running
194         }
195     }
196     if((!running) && (spark) && (fuel) && (crank_counter > 120)) {
197         // start the engine if revs high enough
198         if(RPM > 450) {
199             // For now just instantaneously start but later we should maybe crank for a bit
200             running = true;
201 //          RPM = 600;
202         }
203     }
204     if( (running) && ((!spark)||(!fuel)) ) {
205         // Cut the engine
206         // note that we only cut the power - the engine may continue to spin if the prop is in a moving airstream
207         running = false;
208     }
209
210     // Now we've ascertained whether the engine is running or not we can start to do the engine calculations 'proper'
211
212     // Calculate Sea Level Manifold Pressure
213     Manifold_Pressure = Calc_Manifold_Pressure( Throttle_Lever_Pos, Max_Manifold_Pressure, Min_Manifold_Pressure );
214     // cout << "manifold pressure = " << Manifold_Pressure << endl;
215
216     //Then find the actual manifold pressure (the calculated one is the sea level pressure)
217     True_Manifold_Pressure = Manifold_Pressure * p_amb / p_amb_sea_level;
218
219     //Do the fuel flow calculations
220     Calc_Fuel_Flow_Gals_Hr();
221
222     //Calculate engine power
223     Calc_Percentage_Power(Magneto_Left, Magneto_Right);
224     HP = Percentage_Power * MaxHP / 100.0;
225     Power_SI = HP * CONVERT_HP_TO_WATTS;
226
227     // FMEP calculation.  For now we will just use this during motored operation.
228     // Eventually we will calculate IMEP and use the FMEP all the time to give BMEP (maybe!)
229     if(!running) {
230         // This FMEP data is from the Patton paper, assumes fully warm conditions.
231         FMEP = 1e-12*pow(RPM,4) - 1e-8*pow(RPM,3) + 5e-5*pow(RPM,2) - 0.0722*RPM + 154.85;
232         // Gives FMEP in kPa - now convert to Pa
233         FMEP *= 1000;
234     } else {
235         FMEP = 0.0;
236     }
237     // Is this total FMEP or friction FMEP ???
238
239     Torque_FMEP = (FMEP * displacement_SI) / (2.0 * LS_PI * n_R);
240
241     // Calculate Engine Torque. Check for div by zero since percentage power correlation does not guarantee zero power at zero rpm.
242     // However this is problematical since there is a resistance to movement even at rest
243     // Ie this is a dynamics equation not a statics one.  This can be solved by going over to MEP based torque calculations.
244     if(RPM == 0) {
245         Torque_SI = 0 - Torque_FMEP;
246     }
247     else {
248         Torque_SI = ((Power_SI * 60.0) / (2.0 * LS_PI * RPM)) - Torque_FMEP;  //Torque = power / angular velocity
249         // cout << Torque << " Nm\n";
250     }
251
252     //Calculate Exhaust gas temperature
253     if(running)
254         Calc_EGT();
255     else
256         EGT = 298.0;
257
258     // Calculate Cylinder Head Temperature
259     Calc_CHT();
260     
261     // Calculate oil temperature
262     current_oil_temp = Calc_Oil_Temp(current_oil_temp);
263     
264     // Calculate Oil Pressure
265     Oil_Pressure = Calc_Oil_Press( Oil_Temp, RPM );
266     
267     // Now do the Propeller Calculations
268     Do_Prop_Calcs();
269     
270 // Now do the engine - prop torque balance to calculate final RPM
271     
272     //Calculate new RPM from torque balance and inertia.
273     Torque_Imbalance = Torque_SI - prop_torque;  //This gives a +ve value when the engine torque exeeds the prop torque
274     // (Engine torque is +ve when it acts in the direction of engine revolution, prop torque is +ve when it opposes the direction of engine revolution)
275     
276     angular_acceleration = Torque_Imbalance / (engine_inertia + prop_inertia);
277     angular_velocity_SI += (angular_acceleration * time_step);
278     // Don't let the engine go into reverse
279     if(angular_velocity_SI < 0)
280         angular_velocity_SI = 0;
281     RPM = (angular_velocity_SI * 60) / (2.0 * LS_PI);
282
283     // And finally a last check on the engine state after doing the torque balance with the prop - have we stalled?
284     if(running) { 
285         //Check if we have stalled the engine
286         if (RPM == 0) {
287             running = false;
288         } else if((RPM <= 480) && (cranking)) {
289             //Make sure the engine noise dosn't play if the engine won't start due to eg mixture lever pulled out.
290             running = false;
291             EGT = 298.0;
292         }
293     }
294
295     // And finally, do any unit conversions from internal units to output units
296     EGT_degF = (EGT * 1.8) - 459.67;
297     CHT_degF = (CHT * 1.8) - 459.67;
298 }
299
300 //*****************************************************************************************************
301
302 // FGNewEngine member functions
303
304 ////////////////////////////////////////////////////////////////////////////////////////////
305 // Return the combustion efficiency as a function of equivalence ratio
306 //
307 // Combustion efficiency values from Heywood, 
308 // "Internal Combustion Engine Fundamentals", ISBN 0-07-100499-8
309 ////////////////////////////////////////////////////////////////////////////////////////////
310 float FGNewEngine::Lookup_Combustion_Efficiency(float thi_actual)
311 {
312     const int NUM_ELEMENTS = 11;
313     float thi[NUM_ELEMENTS] = {0.0, 0.9, 1.0, 1.05, 1.1, 1.15, 1.2, 1.3, 1.4, 1.5, 1.6};  //array of equivalence ratio values
314     float neta_comb[NUM_ELEMENTS] = {0.98, 0.98, 0.97, 0.95, 0.9, 0.85, 0.79, 0.7, 0.63, 0.57, 0.525};  //corresponding array of combustion efficiency values
315     float neta_comb_actual = 0.0f;
316     float factor;
317
318     int i;
319     int j = NUM_ELEMENTS;  //This must be equal to the number of elements in the lookup table arrays
320
321     for(i=0;i<j;i++)
322     {
323         if(i == (j-1)) {
324             // Assume linear extrapolation of the slope between the last two points beyond the last point
325             float dydx = (neta_comb[i] - neta_comb[i-1]) / (thi[i] - thi[i-1]);
326             neta_comb_actual = neta_comb[i] + dydx * (thi_actual - thi[i]);
327             return neta_comb_actual;
328         }
329         if(thi_actual == thi[i]) {
330             neta_comb_actual = neta_comb[i];
331             return neta_comb_actual;
332         }
333         if((thi_actual > thi[i]) && (thi_actual < thi[i + 1])) {
334             //do linear interpolation between the two points
335             factor = (thi_actual - thi[i]) / (thi[i+1] - thi[i]);
336             neta_comb_actual = (factor * (neta_comb[i+1] - neta_comb[i])) + neta_comb[i];
337             return neta_comb_actual;
338         }
339     }
340
341     //if we get here something has gone badly wrong
342     SG_LOG(SG_AIRCRAFT, SG_ALERT, "error in FGNewEngine::Lookup_Combustion_Efficiency");
343     return neta_comb_actual;
344 }
345
346 ////////////////////////////////////////////////////////////////////////////////////////////
347 // Return the percentage of best mixture power available at a given mixture strength
348 //
349 // Based on data from "Technical Considerations for Catalysts for the European Market"
350 // by H S Gandi, published 1988 by IMechE
351 //
352 // Note that currently no attempt is made to set a lean limit on stable combustion
353 ////////////////////////////////////////////////////////////////////////////////////////////
354 float FGNewEngine::Power_Mixture_Correlation(float thi_actual)
355 {
356     float AFR_actual = 14.7 / thi_actual;
357     // thi and thi_actual are equivalence ratio
358     const int NUM_ELEMENTS = 13;
359     // The lookup table is in AFR because the source data was.  I added the two end elements to make sure we are almost always in it.
360     float AFR[NUM_ELEMENTS] = {(14.7/1.6), 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, (14.7/0.6)};  //array of equivalence ratio values
361     float mixPerPow[NUM_ELEMENTS] = {78, 86, 93.5, 98, 100, 99, 96.4, 92.5, 88, 83, 78.5, 74, 58};  //corresponding array of combustion efficiency values
362     float mixPerPow_actual = 0.0f;
363     float factor;
364     float dydx;
365
366     int i;
367     int j = NUM_ELEMENTS;  //This must be equal to the number of elements in the lookup table arrays
368
369     for(i=0;i<j;i++)
370     {
371         if(i == (j-1)) {
372             // Assume linear extrapolation of the slope between the last two points beyond the last point
373             dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
374             mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
375             return mixPerPow_actual;
376         }
377         if((i == 0) && (AFR_actual < AFR[i])) {
378             // Assume linear extrapolation of the slope between the first two points for points before the first point
379             dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
380             mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
381             return mixPerPow_actual;
382         }
383         if(AFR_actual == AFR[i]) {
384             mixPerPow_actual = mixPerPow[i];
385             return mixPerPow_actual;
386         }
387         if((AFR_actual > AFR[i]) && (AFR_actual < AFR[i + 1])) {
388             //do linear interpolation between the two points
389             factor = (AFR_actual - AFR[i]) / (AFR[i+1] - AFR[i]);
390             mixPerPow_actual = (factor * (mixPerPow[i+1] - mixPerPow[i])) + mixPerPow[i];
391             return mixPerPow_actual;
392         }
393     }
394
395     //if we get here something has gone badly wrong
396     SG_LOG(SG_AIRCRAFT, SG_ALERT, "error in FGNewEngine::Power_Mixture_Correlation");
397     return mixPerPow_actual;
398 }
399
400 // Calculate Cylinder Head Temperature
401 // Crudely models the cylinder head as an arbitary lump of arbitary size and area with one third of combustion energy
402 // as heat input and heat output as a function of airspeed and temperature.  Could be improved!!!
403 void FGNewEngine::Calc_CHT()
404 {
405     float h1 = -95.0;   //co-efficient for free convection
406     float h2 = -3.95;   //co-efficient for forced convection
407     float h3 = -0.05;   //co-efficient for forced convection due to prop backwash
408     float v_apparent;   //air velocity over cylinder head in m/s
409     float v_dot_cooling_air;
410     float m_dot_cooling_air;
411     float temperature_difference;
412     float arbitary_area = 1.0;
413     float dqdt_from_combustion;
414     float dqdt_forced;      //Rate of energy transfer to/from cylinder head due to forced convection (Joules) (sign convention: to cylinder head is +ve)
415     float dqdt_free;        //Rate of energy transfer to/from cylinder head due to free convection (Joules) (sign convention: to cylinder head is +ve)
416     float dqdt_cylinder_head;       //Overall energy change in cylinder head
417     float CpCylinderHead = 800.0;   //FIXME - this is a guess - I need to look up the correct value
418     float MassCylinderHead = 8.0;   //Kg - this is a guess - it dosn't have to be absolutely accurate but can be varied to alter the warm-up rate
419     float HeatCapacityCylinderHead;
420     float dCHTdt;
421
422     // The above values are hardwired to give reasonable results for an IO360 (C172 engine)
423     // Now adjust to try to compensate for arbitary sized engines - this is currently untested
424     arbitary_area *= (MaxHP / 180.0);
425     MassCylinderHead *= (MaxHP / 180.0);
426
427     temperature_difference = CHT - T_amb;
428
429     v_apparent = IAS * 0.5144444;  //convert from knots to m/s
430     v_dot_cooling_air = arbitary_area * v_apparent;
431     m_dot_cooling_air = v_dot_cooling_air * rho_air;
432
433     //Calculate rate of energy transfer to cylinder head from combustion
434     dqdt_from_combustion = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
435     
436     //Calculate rate of energy transfer from cylinder head due to cooling  NOTE is calculated as rate to but negative
437     dqdt_forced = (h2 * m_dot_cooling_air * temperature_difference) + (h3 * RPM * temperature_difference);
438     dqdt_free = h1 * temperature_difference;
439     
440     //Calculate net rate of energy transfer to or from cylinder head
441     dqdt_cylinder_head = dqdt_from_combustion + dqdt_forced + dqdt_free;
442     
443     HeatCapacityCylinderHead = CpCylinderHead * MassCylinderHead;
444     
445     dCHTdt = dqdt_cylinder_head / HeatCapacityCylinderHead;
446     
447     CHT += (dCHTdt * time_step);
448 }
449
450 // Calculate exhaust gas temperature
451 void FGNewEngine::Calc_EGT()
452 {
453     combustion_efficiency = Lookup_Combustion_Efficiency(equivalence_ratio);  //The combustion efficiency basically tells us what proportion of the fuels calorific value is released
454
455     //now calculate energy release to exhaust
456     //We will assume a three way split of fuel energy between useful work, the coolant system and the exhaust system
457     //This is a reasonable first suck of the thumb estimate for a water cooled automotive engine - whether it holds for an air cooled aero engine is probably open to question
458     //Regardless - it won't affect the variation of EGT with mixture, and we can always put a multiplier on EGT to get a reasonable peak value.
459     enthalpy_exhaust = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
460     heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
461     delta_T_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
462
463     EGT = T_amb + delta_T_exhaust;
464
465     //The above gives the exhaust temperature immediately prior to leaving the combustion chamber
466     //Now derate to give a more realistic figure as measured downstream
467     //For now we will aim for a peak of around 400 degC (750 degF)
468
469     EGT *= 0.444 + ((0.544 - 0.444) * Percentage_Power / 100.0);
470 }
471
472 // Calculate Manifold Pressure based on Throttle lever Position
473 float FGNewEngine::Calc_Manifold_Pressure ( float LeverPosn, float MaxMan, float MinMan)
474 {
475     float Inches;
476
477     //Note that setting the manifold pressure as a function of lever position only is not strictly accurate
478     //MAP is also a function of engine speed. (and ambient pressure if we are going for an actual MAP model)
479     Inches = MinMan + (LeverPosn * (MaxMan - MinMan) / 100);
480
481     //allow for idle bypass valve or slightly open throttle stop
482     if(Inches < MinMan)
483         Inches = MinMan;
484
485     //Check for stopped engine (crudest way of compensating for engine speed)
486     if(RPM == 0)
487         Inches = 29.92;
488
489     return Inches;
490 }
491
492 // Calculate fuel flow in gals/hr
493 void FGNewEngine::Calc_Fuel_Flow_Gals_Hr()
494 {
495     //DCL - calculate mass air flow into engine based on speed and load - separate this out into a function eventually
496     //t_amb is actual temperature calculated from altitude
497     //calculate density from ideal gas equation
498     rho_air = p_amb / ( R_air * T_amb );
499     rho_air_manifold = rho_air * Manifold_Pressure / 29.6;  //This is a bit of a roundabout way of calculating this but it works !!  If we put manifold pressure into SI units we could do it simpler.
500     //calculate ideal engine volume inducted per second
501     swept_volume = (displacement_SI * (RPM / 60)) / 2;  //This equation is only valid for a four stroke engine
502     //calculate volumetric efficiency - for now we will just use 0.8, but actually it is a function of engine speed and the exhaust to manifold pressure ratio
503     //Note that this is cylinder vol eff - the throttle drop is already accounted for in the MAP calculation
504     volumetric_efficiency = 0.8;
505     //Now use volumetric efficiency to calculate actual air volume inducted per second
506     v_dot_air = swept_volume * volumetric_efficiency;
507     //Now calculate mass flow rate of air into engine
508     m_dot_air = v_dot_air * rho_air_manifold;
509
510 //**************
511
512     //DCL - now calculate fuel flow into engine based on air flow and mixture lever position
513     //assume lever runs from no flow at fully out to thi = 1.3 at fully in at sea level
514     //also assume that the injector linkage is ideal - hence the set mixture is maintained at a given altitude throughout the speed and load range
515     thi_sea_level = 1.3 * ( Mixture_Lever_Pos / 100.0 );
516     equivalence_ratio = thi_sea_level * p_amb_sea_level / p_amb; //ie as we go higher the mixture gets richer for a given lever position
517     m_dot_fuel = m_dot_air / 14.7 * equivalence_ratio;
518     Fuel_Flow_gals_hr = (m_dot_fuel / rho_fuel) * 264.172 * 3600.0;  // Note this assumes US gallons
519 }
520
521 // Calculate the percentage of maximum rated power delivered as a function of Manifold pressure multiplied by engine speed (rpm)
522 // This is not necessarilly the best approach but seems to work for now.
523 // May well need tweaking at the bottom end if the prop model is changed.
524 void FGNewEngine::Calc_Percentage_Power(bool mag_left, bool mag_right)
525 {
526     // For a given Manifold Pressure and RPM calculate the % Power
527     // Multiply Manifold Pressure by RPM
528     float ManXRPM = True_Manifold_Pressure * RPM;
529
530 /*
531 //  Phil's %power correlation
532     //  Calculate % Power
533     Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM) + ( + 7E-04 * ManXRPM) - 0.1218;
534     // cout << Percentage_Power <<  "%" << "\t";
535 */
536
537 // DCL %power correlation - basically Phil's correlation modified to give slighty less power at the low end
538 // might need some adjustment as the prop model is adjusted
539 // My aim is to match the prop model and engine model at the low end to give the manufacturer's recommended idle speed with the throttle closed - 600rpm for the Continental IO520
540     //  Calculate % Power for Nev's prop model
541     //Percentage_Power = (+ 6E-09 * ManXRPM * ManXRPM) + ( + 8E-04 * ManXRPM) - 1.8524;
542
543     // Calculate %power for DCL prop model
544     Percentage_Power = (7e-9 * ManXRPM * ManXRPM) + (7e-4 * ManXRPM) - 1.0;
545
546     // Power de-rating for altitude has been removed now that we are basing the power
547     // on the actual manifold pressure, which takes air pressure into account.  However - this fails to
548     // take the temperature into account - this is TODO.
549
550     // Adjust power for temperature - this is temporary until the power is done as a function of mass flow rate induced
551     // Adjust for Temperature - Temperature above Standard decrease
552     // power by 7/120 % per degree F increase, and incease power for
553     // temps below at the same ratio
554     float T_amb_degF = (T_amb * 1.8) - 459.67;
555     float T_amb_sea_lev_degF = (288 * 1.8) - 459.67; 
556     Percentage_Power = Percentage_Power + ((T_amb_sea_lev_degF - T_amb_degF) * 7 /120);
557
558     //DCL - now adjust power to compensate for mixture
559     Percentage_of_best_power_mixture_power = Power_Mixture_Correlation(equivalence_ratio);
560     Percentage_Power = Percentage_Power * Percentage_of_best_power_mixture_power / 100.0;
561
562     // Now Derate engine for the effects of Bad/Switched off magnetos
563     //if (Magneto_Left == 0 && Magneto_Right == 0) {
564     if(!running) {
565         // cout << "Both OFF\n";
566         Percentage_Power = 0;
567     } else if (mag_left && mag_right) {
568         // cout << "Both On    ";
569     } else if (mag_left == 0 || mag_right== 0) {
570         // cout << "1 Magneto Failed   ";
571         Percentage_Power = Percentage_Power * ((100.0 - Mag_Derate_Percent)/100.0);
572         //  cout << FGEng1_Percentage_Power <<  "%" << "\t";
573     }
574 /*
575     //DCL - stall the engine if RPM drops below 450 - this is possible if mixture lever is pulled right out
576     //This is a kludge that I should eliminate by adding a total fmep estimation.
577     if(RPM < 450)
578         Percentage_Power = 0;
579 */
580     if(Percentage_Power < 0)
581         Percentage_Power = 0;
582 }
583
584 // Calculate Oil Temperature in degrees Kelvin
585 float FGNewEngine::Calc_Oil_Temp (float oil_temp)
586 {
587     float idle_percentage_power = 2.3;  // approximately
588     float target_oil_temp;          // Steady state oil temp at the current engine conditions
589     float time_constant;            // The time constant for the differential equation
590     if(running) {
591         target_oil_temp = 363;
592         time_constant = 500;        // Time constant for engine-on idling.
593         if(Percentage_Power > idle_percentage_power) {
594             time_constant /= ((Percentage_Power / idle_percentage_power) / 10.0);       // adjust for power 
595         }
596     } else {
597         target_oil_temp = 298;
598         time_constant = 1000;  // Time constant for engine-off; reflects the fact that oil is no longer getting circulated
599     }
600
601     float dOilTempdt = (target_oil_temp - oil_temp) / time_constant;
602
603     oil_temp += (dOilTempdt * time_step);
604
605     return (oil_temp);
606 }
607
608 // Calculate Oil Pressure
609 float FGNewEngine::Calc_Oil_Press (float Oil_Temp, float Engine_RPM)
610 {
611     float Oil_Pressure = 0;                     //PSI
612     float Oil_Press_Relief_Valve = 60;  //PSI
613     float Oil_Press_RPM_Max = 1800;
614     float Design_Oil_Temp = 85;         //Celsius
615     float Oil_Viscosity_Index = 0.25;   // PSI/Deg C
616 //    float Temp_Deviation = 0;         // Deg C
617
618     Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
619
620     // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
621     if (Oil_Pressure >= Oil_Press_Relief_Valve) {
622         Oil_Pressure = Oil_Press_Relief_Valve;
623     }
624
625     // Now adjust pressure according to Temp which affects the viscosity
626
627     Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
628
629     return Oil_Pressure;
630 }
631
632
633 // Propeller calculations.
634 void FGNewEngine::Do_Prop_Calcs()
635 {
636     float Gear_Ratio = 1.0;
637     float forward_velocity;             // m/s
638     float prop_power_consumed_SI;       // Watts
639     double J;                           // advance ratio - dimensionless
640     double Cp_20;                   // coefficient of power for 20 degree blade angle
641     double Cp_25;                   // coefficient of power for 25 degree blade angle
642     double Cp;                      // Our actual coefficient of power
643     double neta_prop_20;
644     double neta_prop_25;
645     double neta_prop;               // prop efficiency
646
647     FGProp1_RPS = RPM * Gear_Ratio / 60.0; 
648     angular_velocity_SI = 2.0 * LS_PI * RPM / 60.0;
649     forward_velocity = IAS * 0.514444444444;        // Convert to m/s
650     
651     if(FGProp1_RPS == 0)
652         J = 0;
653     else
654         J = forward_velocity / (FGProp1_RPS * prop_diameter);
655     //cout << "advance_ratio = " << J << '\n';
656     
657     //Cp correlations based on data from McCormick
658     Cp_20 = 0.0342*J*J*J*J - 0.1102*J*J*J + 0.0365*J*J - 0.0133*J + 0.064;
659     Cp_25 = 0.0119*J*J*J*J - 0.0652*J*J*J + 0.018*J*J - 0.0077*J + 0.0921;
660     
661     //cout << "Cp_20 = " << Cp_20 << '\n';
662     //cout << "Cp_25 = " << Cp_25 << '\n';
663     
664     //Assume that the blade angle is between 20 and 25 deg - reasonable for fixed pitch prop but won't hold for variable one !!!
665     Cp = Cp_20 + ( (Cp_25 - Cp_20) * ((blade_angle - 20)/(25 - 20)) );
666     //cout << "Cp = " << Cp << '\n';
667     //cout << "RPM = " << RPM << '\n';
668     //cout << "angular_velocity_SI = " << angular_velocity_SI << '\n';
669     
670     prop_power_consumed_SI = Cp * rho_air * pow(FGProp1_RPS,3.0f) * pow(float(prop_diameter),5.0f);
671     //cout << "prop HP consumed = " << prop_power_consumed_SI / 745.699 << '\n';
672     if(angular_velocity_SI == 0)
673         prop_torque = 0;
674         // However this can give problems - if rpm == 0 but forward velocity increases the prop should be able to generate a torque to start the engine spinning
675         // Unlikely to happen in practice - but I suppose someone could move the lever of a stopped large piston engine from feathered to windmilling.
676         // This *does* give problems - if the plane is put into a steep climb whilst windmilling the engine friction will eventually stop it spinning.
677         // When put back into a dive it never starts re-spinning again.  Although it is unlikely that anyone would do this in real life, they might well do it in a sim!!!
678     else
679         prop_torque = prop_power_consumed_SI / angular_velocity_SI;
680     
681     // calculate neta_prop here
682     neta_prop_20 = 0.1328*J*J*J*J - 1.3073*J*J*J + 0.3525*J*J + 1.5591*J + 0.0007;
683     neta_prop_25 = -0.3121*J*J*J*J + 0.4234*J*J*J - 0.7686*J*J + 1.5237*J - 0.0004;
684     neta_prop = neta_prop_20 + ( (neta_prop_25 - neta_prop_20) * ((blade_angle - 20)/(25 - 20)) );
685     
686     // Check for zero forward velocity to avoid divide by zero
687     if(forward_velocity < 0.0001)
688         prop_thrust = 0.0;
689         // I don't see how this works - how can the plane possibly start from rest ???
690         // Hmmmm - it works because the forward_velocity at present never drops below about 0.03 even at rest
691         // We can't really rely on this in the future - needs fixing !!!!
692         // The problem is that we're doing this calculation backwards - we're working out the thrust from the power consumed and the velocity, which becomes invalid as velocity goes to zero.
693         // It would be far more natural to work out the power consumed from the thrust - FIXME!!!!!.
694     else
695         prop_thrust = neta_prop * prop_power_consumed_SI / forward_velocity;       //TODO - rename forward_velocity to IAS_SI
696     //cout << "prop_thrust = " << prop_thrust << '\n';
697 }