2 // Author: Phil Schubert
3 // Date started: 12/03/99
4 // Purpose: Models a Continental IO-520-M Engine
5 // Called by: FGSimExec
7 // Copyright (C) 1999 Philip L. Schubert (philings@ozemail.com.au)
9 // This program is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU General Public License as
11 // published by the Free Software Foundation; either version 2 of the
12 // License, or (at your option) any later version.
14 // This program is distributed in the hope that it will be useful, but
15 // WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 // General Public License for more details.
19 // You should have received a copy of the GNU General Public License
20 // along with this program; if not, write to the Free Software
21 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24 // Further information about the GNU General Public License can also
25 // be found on the world wide web at http://www.gnu.org.
27 // FUNCTIONAL DESCRIPTION
28 // ------------------------------------------------------------------------
29 // Models a Continental IO-520-M engine. This engine is used in Cessna
30 // 210, 310, Beechcraft Bonaza and Baron C55. The equations used below
31 // were determined by a first and second order curve fits using Excel.
32 // The data is from the Cessna Aircraft Corporations Engine and Flight
33 // Computer for C310. Part Number D3500-13
36 // ------------------------------------------------------------------------
40 // ------------------------------------------------------------------------
41 // 12/03/99 PLS Created
42 // 07/03/99 PLS Added Calculation of Density, and Prop_Torque
43 // 07/03/99 PLS Restructered Variables to allow easier implementation
45 // 15/03/99 PLS Added Oil Pressure, Oil Temperature and CH Temp
46 // ------------------------------------------------------------------------
48 // ------------------------------------------------------------------------
51 /////////////////////////////////////////////////////////////////////
53 // Modified by Dave Luff (david.luff@nottingham.ac.uk) September 2000
55 // Altered manifold pressure range to add a minimum value at idle to simulate the throttle stop / idle bypass valve,
56 // and to reduce the maximum value whilst the engine is running to slightly below ambient to account for CdA losses across the throttle
58 // Altered it a bit to model an IO360 from C172 - 360 cubic inches, 180 HP max, fixed pitch prop
59 // Added a simple fixed pitch prop model by Nev Harbor - this is not intended as a final model but simply a hack to get it running for now
60 // I used Phil's ManXRPM correlation for power rather than do a new one for the C172 for now, but altered it a bit to reduce power at the low end
62 // Added EGT model based on combustion efficiency and an energy balance with the exhaust gases
64 // Added a mixture - power correlation based on a curve in the IO360 operating manual
66 // I've tried to match the prop and engine model to give roughly 600 RPM idle and 180 HP at 2700 RPM
67 // but it is by no means currently at a completed stage - DCL 15/9/00
69 // DCL 28/09/00 - Added estimate of engine and prop inertia and changed engine speed calculation to be calculated from Angular acceleration = Torque / Inertia.
70 // Requires a timestep to be passed to FGNewEngine::init and currently assumes this timestep does not change.
71 // Could easily be altered to pass a variable timestep to FGNewEngine::update every step instead if required.
73 // DCL 27/10/00 - Added first stab at cylinder head temperature model
74 // See the comment block in the code for details
76 // DCL 02/11/00 - Modified EGT code to reduce values to those more representative of a sensor downstream
78 // DCL 02/02/01 - Changed the prop model to one based on efficiency and co-efficient of power curves from McCormick instead of the
79 // blade element method we were using previously. This works much better, and is similar to how Jon is doing it in JSBSim.
81 // DCL 08/02/01 - Overhauled fuel consumption rate support.
83 // DCL 22/03/01 - Added input of actual air pressure and temperature (and hence density) to the model. Hence the power correlation
84 // with pressure height and temperature is no longer required since the power is based on the actual manifold pressure.
86 // DCL 22/03/01 - based on Riley's post on the list (25 rpm gain at 1000 rpm as lever is pulled out from full rich)
87 // I have reduced the sea level full rich mixture to thi = 1.3
89 // DCL 18/9/01 - Got the engine to start and stop in response to the magneto switch.
90 // Changed all PI to LS_PI (in ls_constants.h).
91 // Engine now checks for fuel and stops when not available.
93 //////////////////////////////////////////////////////////////////////
95 #include <simgear/compiler.h>
100 #include STL_IOSTREAM
102 #if !defined(SG_HAVE_NATIVE_SGI_COMPILERS)
107 #include "LaRCsim/ls_constants.h"
109 #include <Main/fg_props.hxx>
111 // Static utility functions
113 // Calculate Density Ratio
114 static float Density_Ratio ( float x )
117 y = ((3E-10 * x * x) - (3E-05 * x) + 0.9998);
122 // Calculate Air Density - Rho, using the ideal gas equation
123 // Takes and returns SI values
124 static float Density ( float temperature, float pressure )
130 float rho = pressure / (R * temperature);
135 // Calculate Speed in FPS given Knots CAS
136 static float IAS_to_FPS (float x)
143 // FGNewEngine member functions
145 float FGNewEngine::Lookup_Combustion_Efficiency(float thi_actual)
147 const int NUM_ELEMENTS = 11;
148 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
149 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
150 //combustion efficiency values from Heywood, "Internal Combustion Engine Fundamentals", ISBN 0-07-100499-8
151 float neta_comb_actual = 0.0f;
155 int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays
160 // Assume linear extrapolation of the slope between the last two points beyond the last point
161 float dydx = (neta_comb[i] - neta_comb[i-1]) / (thi[i] - thi[i-1]);
162 neta_comb_actual = neta_comb[i] + dydx * (thi_actual - thi[i]);
163 return neta_comb_actual;
165 if(thi_actual == thi[i]) {
166 neta_comb_actual = neta_comb[i];
167 return neta_comb_actual;
169 if((thi_actual > thi[i]) && (thi_actual < thi[i + 1])) {
170 //do linear interpolation between the two points
171 factor = (thi_actual - thi[i]) / (thi[i+1] - thi[i]);
172 neta_comb_actual = (factor * (neta_comb[i+1] - neta_comb[i])) + neta_comb[i];
173 return neta_comb_actual;
177 //if we get here something has gone badly wrong
178 cout << "ERROR: error in FGNewEngine::Lookup_Combustion_Efficiency\n";
179 return neta_comb_actual;
182 ////////////////////////////////////////////////////////////////////////////////////////////
183 // Return the percentage of best mixture power available at a given mixture strength
185 // Based on data from "Technical Considerations for Catalysts for the European Market"
186 // by H S Gandi, published 1988 by IMechE
188 // Note that currently no attempt is made to set a lean limit on stable combustion
189 ////////////////////////////////////////////////////////////////////////////////////////////
190 float FGNewEngine::Power_Mixture_Correlation(float thi_actual)
192 float AFR_actual = 14.7 / thi_actual;
193 // thi and thi_actual are equivalence ratio
194 const int NUM_ELEMENTS = 13;
195 // 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.
196 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
197 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
198 float mixPerPow_actual = 0.0f;
203 int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays
208 // Assume linear extrapolation of the slope between the last two points beyond the last point
209 dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
210 mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
211 return mixPerPow_actual;
213 if((i == 0) && (AFR_actual < AFR[i])) {
214 // Assume linear extrapolation of the slope between the first two points for points before the first point
215 dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
216 mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
217 return mixPerPow_actual;
219 if(AFR_actual == AFR[i]) {
220 mixPerPow_actual = mixPerPow[i];
221 return mixPerPow_actual;
223 if((AFR_actual > AFR[i]) && (AFR_actual < AFR[i + 1])) {
224 //do linear interpolation between the two points
225 factor = (AFR_actual - AFR[i]) / (AFR[i+1] - AFR[i]);
226 mixPerPow_actual = (factor * (mixPerPow[i+1] - mixPerPow[i])) + mixPerPow[i];
227 return mixPerPow_actual;
231 //if we get here something has gone badly wrong
232 cout << "ERROR: error in FGNewEngine::Power_Mixture_Correlation\n";
233 return mixPerPow_actual;
238 // Calculate Manifold Pressure based on Throttle lever Position
239 float FGNewEngine::Calc_Manifold_Pressure ( float LeverPosn, float MaxMan, float MinMan)
246 //Note that setting the manifold pressure as a function of lever position only is not strictly accurate
247 //MAP is also a function of engine speed. (and ambient pressure if we are going for an actual MAP model)
248 Inches = MinMan + (LeverPosn * (MaxMan - MinMan) / 100);
250 //allow for idle bypass valve or slightly open throttle stop
260 // Calculate Oil Temperature
261 float FGNewEngine::Calc_Oil_Temp (float Fuel_Flow, float Mixture, float IAS)
268 // Calculate Oil Pressure
269 float FGNewEngine::Calc_Oil_Press (float Oil_Temp, float Engine_RPM)
271 float Oil_Pressure = 0; //PSI
272 float Oil_Press_Relief_Valve = 60; //PSI
273 float Oil_Press_RPM_Max = 1800;
274 float Design_Oil_Temp = 85; //Celsius
275 float Oil_Viscosity_Index = 0.25; // PSI/Deg C
276 float Temp_Deviation = 0; // Deg C
278 Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
280 // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
281 if (Oil_Pressure >= Oil_Press_Relief_Valve)
283 Oil_Pressure = Oil_Press_Relief_Valve;
286 // Now adjust pressure according to Temp which affects the viscosity
288 Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
293 //*************************************************************************************
294 // Initialise the engine model
295 void FGNewEngine::init(double dt) {
297 // These constants should probably be moved eventually
298 CONVERT_CUBIC_INCHES_TO_METERS_CUBED = 1.638706e-5;
299 CONVERT_HP_TO_WATTS = 745.6999;
301 // Properties of working fluids
302 Cp_air = 1005; // J/KgK
303 Cp_fuel = 1700; // J/KgK
304 calorific_value_fuel = 47.3e6; // W/Kg Note that this is only an approximate value
305 rho_fuel = 800; // kg/m^3 - an estimate for now
307 p_amb_sea_level = 101325;
309 // Control and environment inputs
311 Throttle_Lever_Pos = 75;
312 Propeller_Lever_Pos = 75;
313 Mixture_Lever_Pos = 100;
317 // Engine Specific Variables.
318 // Will be set in a parameter file to be read in to create
319 // and instance for each engine.
320 Max_Manifold_Pressure = 28.50; //Inches Hg. An approximation - should be able to find it in the engine performance data
321 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
323 Min_RPM = 600; //Recommended idle from Continental data sheet
325 Mag_Derate_Percent = 5;
326 // MaxHP = 285; //Continental IO520-M
327 MaxHP = 200; //Lycoming IO360 -A-C-D series
328 // MaxHP = 180; //Current Lycoming IO360 ?
329 // displacement = 520; //Continental IO520-M
330 displacement = 360; //Lycoming IO360
331 displacement_SI = displacement * CONVERT_CUBIC_INCHES_TO_METERS_CUBED;
332 engine_inertia = 0.2; //kgm^2 - value taken from a popular family saloon car engine - need to find an aeroengine value !!!!!
333 prop_inertia = 0.05; //kgm^2 - this value is a total guess - dcl
335 n_R = 2; // Number of crank revolutions per power cycle - 2 for a 4 cylinder engine.
337 running = fgGetBool("/engines/engine[0]/running");
339 fgSetBool("/engines/engine[0]/cranking", false);
341 // Initialise Engine Variables used by this instance
346 Percentage_Power = 0;
347 Manifold_Pressure = 29.00; // Inches
348 Fuel_Flow_gals_hr = 0;
351 CHT = 298.0; //deg Kelvin
352 CHT_degF = (CHT * 1.8) - 459.67; //deg Fahrenheit
354 Oil_Pressure = 0; // PSI
355 Oil_Temp = 85; // Deg C
358 Torque_Imbalance = 0;
360 // Initialise Propellor Variables used by this instance
363 FGProp1_Blade_Angle = 13.5;
364 prop_diameter = 1.8; // meters
365 blade_angle = 23.0; // degrees
369 //*****************************************************************************
370 //*****************************************************************************
371 // update the engine model based on current control positions
372 void FGNewEngine::update() {
375 // Hack for testing - should output every 5 seconds
376 static int count1 = 0;
378 // cout << "P_atmos = " << p_amb << " T_atmos = " << T_amb << '\n';
379 // cout << "Manifold pressure = " << Manifold_Pressure << " True_Manifold_Pressure = " << True_Manifold_Pressure << '\n';
380 // cout << "p_amb_sea_level = " << p_amb_sea_level << '\n';
381 // cout << "equivalence_ratio = " << equivalence_ratio << '\n';
382 // cout << "combustion_efficiency = " << combustion_efficiency << '\n';
383 // cout << "AFR = " << 14.7 / equivalence_ratio << '\n';
384 // cout << "Mixture lever = " << Mixture_Lever_Pos << '\n';
385 // cout << "n = " << RPM << " rpm\n";
386 // cout << "T_amb = " << T_amb << '\n';
387 cout << "running = " << running << '\n';
388 cout << "fuel = " << fgGetFloat("/consumables/fuel/tank[0]/level-gal_us") << '\n';
399 // Parameters that alter the operation of the engine. (spark, fuel, starter motor etc)
401 int Magneto_Left = 0;
402 int Magneto_Right = 0;
403 int mag_pos = fgGetInt("/engines/engine[0]/magneto");
404 // Magneto positions:
413 } // neglects battery voltage, master on switch, etc for now.
414 if((mag_pos == 1) || (mag_pos > 2))
419 // crude check for fuel
420 if((fgGetFloat("/consumables/fuel/tank[0]/level-gal_us") > 0) || (fgGetFloat("/consumables/fuel/tank[1]/level-gal_us") > 0)) {
424 } // Need to make this better, eg position of fuel selector switch.
426 // Check if we are turning the starter motor
427 bool temp = fgGetBool("/engines/engine[0]/starter");
428 if(cranking != temp) {
429 // This check saves .../cranking from getting updated every loop - they only update when changed.
432 fgSetBool("/engines/engine[0]/cranking", true);
434 fgSetBool("/engines/engine[0]/cranking", false);
436 // Note that although /engines/engine[0]/starter and /engines/engine[0]/cranking might appear to be duplication it is
437 // not since the starter may be engaged with the battery voltage too low for cranking to occur (or perhaps the master
438 // switch just left off) and the sound manager will read .../cranking to determine wether to play a cranking sound.
439 // For now though none of that is implemented so cranking can be set equal to .../starter without further checks.
441 int Alternate_Air_Pos =0; // Off = 0. Reduces power by 3 % for same throttle setting
442 // DCL - don't know what this Alternate_Air_Pos is - this is a leftover from the Schubert code.
444 //Check mode of engine operation
451 // consider making a horrible noise if the starter is engaged with the engine running
454 if((!running) && (spark) && (fuel)) {
455 // start the engine if revs high enough
457 // For now just instantaneously start but later we should maybe crank for a bit
459 fgSetBool("/engines/engine[0]/running", true);
463 if( (running) && ((!spark)||(!fuel)) ) {
465 // note that we only cut the power - the engine may continue to spin if the prop is in a moving airstream
467 fgSetBool("/engines/engine[0]/running", false);
470 // Calculate Sea Level Manifold Pressure
471 Manifold_Pressure = Calc_Manifold_Pressure( Throttle_Lever_Pos, Max_Manifold_Pressure, Min_Manifold_Pressure );
472 // cout << "manifold pressure = " << Manifold_Pressure << endl;
474 //Then find the actual manifold pressure (the calculated one is the sea level pressure)
475 True_Manifold_Pressure = Manifold_Pressure * p_amb / p_amb_sea_level;
478 //DCL - next calculate m_dot_air and m_dot_fuel into engine
480 //DCL - calculate mass air flow into engine based on speed and load - separate this out into a function eventually
481 //t_amb is actual temperature calculated from altitude
482 //calculate density from ideal gas equation
483 rho_air = p_amb / ( R_air * T_amb );
484 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.
485 //calculate ideal engine volume inducted per second
486 swept_volume = (displacement_SI * (RPM / 60)) / 2; //This equation is only valid for a four stroke engine
487 //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
488 //Note that this is cylinder vol eff - the throttle drop is already accounted for in the MAP calculation
489 volumetric_efficiency = 0.8;
490 //Now use volumetric efficiency to calculate actual air volume inducted per second
491 v_dot_air = swept_volume * volumetric_efficiency;
492 //Now calculate mass flow rate of air into engine
493 m_dot_air = v_dot_air * rho_air_manifold;
497 //DCL - now calculate fuel flow into engine based on air flow and mixture lever position
498 //assume lever runs from no flow at fully out to thi = 1.3 at fully in at sea level
499 //also assume that the injector linkage is ideal - hence the set mixture is maintained at a given altitude throughout the speed and load range
500 thi_sea_level = 1.3 * ( Mixture_Lever_Pos / 100.0 );
501 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
502 m_dot_fuel = m_dot_air / 14.7 * equivalence_ratio;
503 Fuel_Flow_gals_hr = (m_dot_fuel / rho_fuel) * 264.172 * 3600.0; // Note this assumes US gallons
505 //***********************************************************************
506 //Engine power and torque calculations
508 // For a given Manifold Pressure and RPM calculate the % Power
509 // Multiply Manifold Pressure by RPM
510 ManXRPM = True_Manifold_Pressure * RPM;
511 // ManXRPM = Manifold_Pressure * RPM;
516 // Phil's %power correlation
518 Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM) + ( + 7E-04 * ManXRPM) - 0.1218;
519 // cout << Percentage_Power << "%" << "\t";
522 // DCL %power correlation - basically Phil's correlation modified to give slighty less power at the low end
523 // might need some adjustment as the prop model is adjusted
524 // 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
525 // Calculate % Power for Nev's prop model
526 //Percentage_Power = (+ 6E-09 * ManXRPM * ManXRPM) + ( + 8E-04 * ManXRPM) - 1.8524;
528 // Calculate %power for DCL prop model
529 Percentage_Power = (7e-9 * ManXRPM * ManXRPM) + (7e-4 * ManXRPM) - 1.0;
531 // Power de-rating for altitude has been removed now that we are basing the power
532 // on the actual manifold pressure, which takes air pressure into account. However - this fails to
533 // take the temperature into account - this is TODO.
535 // Adjust power for temperature - this is temporary until the power is done as a function of mass flow rate induced
536 // Adjust for Temperature - Temperature above Standard decrease
537 // power by 7/120 % per degree F increase, and incease power for
538 // temps below at the same ratio
539 float T_amb_degF = (T_amb * 1.8) - 459.67;
540 float T_amb_sea_lev_degF = (288 * 1.8) - 459.67;
541 Percentage_Power = Percentage_Power + ((T_amb_sea_lev_degF - T_amb_degF) * 7 /120);
543 //DCL - now adjust power to compensate for mixture
544 Percentage_of_best_power_mixture_power = Power_Mixture_Correlation(equivalence_ratio);
545 Percentage_Power = Percentage_Power * Percentage_of_best_power_mixture_power / 100.0;
547 // Now Derate engine for the effects of Bad/Switched off magnetos
548 //if (Magneto_Left == 0 && Magneto_Right == 0) {
550 // cout << "Both OFF\n";
551 Percentage_Power = 0;
552 } else if (Magneto_Left && Magneto_Right) {
553 // cout << "Both On ";
554 } else if (Magneto_Left == 0 || Magneto_Right== 0) {
555 // cout << "1 Magneto Failed ";
556 Percentage_Power = Percentage_Power * ((100.0 - Mag_Derate_Percent)/100.0);
557 // cout << FGEng1_Percentage_Power << "%" << "\t";
560 //DCL - stall the engine if RPM drops below 450 - this is possible if mixture lever is pulled right out
561 //This is a kludge that I should eliminate by adding a total fmep estimation.
563 Percentage_Power = 0;
565 if(Percentage_Power < 0)
566 Percentage_Power = 0;
568 // FMEP calculation. For now we will just use this during motored operation, ie when %Power == 0.
569 // Eventually we will calculate IMEP and use the FMEP all the time to give BMEP
571 if(Percentage_Power == 0) {
572 // This FMEP data is from the Patton paper, assumes fully warm conditions.
573 FMEP = 1e-12*pow(RPM,4) - 1e-8*pow(RPM,3) + 5e-5*pow(RPM,2) - 0.0722*RPM + 154.85;
574 // Gives FMEP in kPa - now convert to Pa
580 Torque_FMEP = (FMEP * displacement_SI) / (2.0 * LS_PI * n_R);
582 HP = Percentage_Power * MaxHP / 100.0;
584 Power_SI = HP * CONVERT_HP_TO_WATTS;
586 // Calculate Engine Torque. Check for div by zero since percentage power correlation does not guarantee zero power at zero rpm.
587 // However this is problematical since there is a resistance to movement even at rest
588 // Ie this is a dynamics equation not a statics one. This can be solved by going over to MEP based torque calculations.
590 Torque_SI = 0 - Torque_FMEP;
593 Torque_SI = ((Power_SI * 60.0) / (2.0 * LS_PI * RPM)) - Torque_FMEP; //Torque = power / angular velocity
594 // cout << Torque << " Nm\n";
597 //**********************************************************************
598 //Calculate Exhaust gas temperature
600 // cout << "Thi = " << equivalence_ratio << '\n';
602 combustion_efficiency = Lookup_Combustion_Efficiency(equivalence_ratio); //The combustion efficiency basically tells us what proportion of the fuels calorific value is released
604 // cout << "Combustion efficiency = " << combustion_efficiency << '\n';
606 //now calculate energy release to exhaust
607 //We will assume a three way split of fuel energy between useful work, the coolant system and the exhaust system
608 //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
609 //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.
610 enthalpy_exhaust = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
611 heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
612 delta_T_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
613 // delta_T_exhaust = Calculate_Delta_T_Exhaust();
615 // cout << "T_amb " << T_amb;
616 // cout << " dT exhaust = " << delta_T_exhaust;
618 EGT = T_amb + delta_T_exhaust;
620 //The above gives the exhaust temperature immediately prior to leaving the combustion chamber
621 //Now derate to give a more realistic figure as measured downstream
622 //For now we will aim for a peak of around 400 degC (750 degF)
624 EGT *= 0.444 + ((0.544 - 0.444) * Percentage_Power / 100.0);
626 EGT_degF = (EGT * 1.8) - 459.67;
628 //cout << " EGT = " << EGT << " degK " << EGT_degF << " degF";// << '\n';
630 //***************************************************************************************
631 // Calculate Cylinder Head Temperature
635 This is a somewhat rough first attempt at modelling cylinder head temperature. The cylinder head
636 is assumed to be at uniform temperature. Obviously this is incorrect, but it simplifies things a
637 lot, and we're just looking for the behaviour of CHT to be correct. Energy transfer to the cylinder
638 head is assumed to be one third of the energy released by combustion at all conditions. This is a
639 reasonable estimate, although obviously in real life it varies with different conditions and possibly
640 with CHT itself. I've split energy transfer from the cylinder head into 2 terms - free convection -
641 ie convection to stationary air, and forced convection, ie convection into flowing air. The basic
642 free convection equation is: dqdt = -hAdT Since we don't know A and are going to set h quite arbitarily
643 anyway I've knocked A out and just wrapped it up in h - the only real significance is that the units
644 of h will be different but that dosn't really matter to us anyway. In addition, we have the problem
645 that the prop model I'm currently using dosn't model the backwash from the prop which will add to the
646 velocity of the cooling air when the prop is turning, so I've added an extra term to try and cope
649 In real life, forced convection equations are genarally empirically derived, and are quite complicated
650 and generally contain such things as the Reynolds and Nusselt numbers to various powers. The best
651 course of action would probably to find an empirical correlation from the literature for a similar
652 situation and try and get it to fit well. However, for now I am using my own made up very simple
653 correlation for the energy transfer from the cylinder head:
655 dqdt = -(h1.dT) -(h2.m_dot.dT) -(h3.rpm.dT)
657 where dT is the temperature different between the cylinder head and the surrounding air, m_dot is the
658 mass flow rate of cooling air through an arbitary volume, rpm is the engine speed in rpm (this is the
659 backwash term), and h1, h2, h3 are co-efficients which we can play with to attempt to get the CHT
660 behaviour to match real life.
662 In order to change the values of CHT that the engine settles down at at various conditions,
663 have a play with h1, h2 and h3. In order to change the rate of heating/cooling without affecting
664 equilibrium values alter the cylinder head mass, which is really quite arbitary. Bear in mind that
665 altering h1, h2 and h3 will also alter the rate of heating or cooling as well as equilibrium values,
666 but altering the cylinder head mass will only alter the rate. It would I suppose be better to read
667 the values from file to avoid the necessity for re-compilation every time I change them.
670 //CHT = Calc_CHT( Fuel_Flow, Mixture, IAS);
671 // cout << "Cylinder Head Temp (F) = " << CHT << endl;
672 float h1 = -95.0; //co-efficient for free convection
673 float h2 = -3.95; //co-efficient for forced convection
674 float h3 = -0.05; //co-efficient for forced convection due to prop backwash
675 float v_apparent; //air velocity over cylinder head in m/s
676 float v_dot_cooling_air;
677 float m_dot_cooling_air;
678 float temperature_difference;
679 float arbitary_area = 1.0;
680 float dqdt_from_combustion;
681 float dqdt_forced; //Rate of energy transfer to/from cylinder head due to forced convection (Joules) (sign convention: to cylinder head is +ve)
682 float dqdt_free; //Rate of energy transfer to/from cylinder head due to free convection (Joules) (sign convention: to cylinder head is +ve)
683 float dqdt_cylinder_head; //Overall energy change in cylinder head
684 float CpCylinderHead = 800.0; //FIXME - this is a guess - I need to look up the correct value
685 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
686 float HeatCapacityCylinderHead;
689 temperature_difference = CHT - T_amb;
691 v_apparent = IAS * 0.5144444; //convert from knots to m/s
692 v_dot_cooling_air = arbitary_area * v_apparent;
693 m_dot_cooling_air = v_dot_cooling_air * rho_air;
695 //Calculate rate of energy transfer to cylinder head from combustion
696 dqdt_from_combustion = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
698 //Calculate rate of energy transfer from cylinder head due to cooling NOTE is calculated as rate to but negative
699 dqdt_forced = (h2 * m_dot_cooling_air * temperature_difference) + (h3 * RPM * temperature_difference);
700 dqdt_free = h1 * temperature_difference;
702 //Calculate net rate of energy transfer to or from cylinder head
703 dqdt_cylinder_head = dqdt_from_combustion + dqdt_forced + dqdt_free;
705 HeatCapacityCylinderHead = CpCylinderHead * MassCylinderHead;
707 dCHTdt = dqdt_cylinder_head / HeatCapacityCylinderHead;
709 CHT += (dCHTdt * time_step);
711 CHT_degF = (CHT * 1.8) - 459.67;
713 //cout << " CHT = " << CHT_degF << " degF\n";
716 // End calculate Cylinder Head Temperature
719 //***************************************************************************************
720 // Oil pressure calculation
722 // Calculate Oil Pressure
723 Oil_Pressure = Calc_Oil_Press( Oil_Temp, RPM );
724 // cout << "Oil Pressure (PSI) = " << Oil_Pressure << endl;
726 //**************************************************************************************
727 // Now do the Propeller Calculations
730 FGProp1_RPS = RPM * Gear_Ratio / 60.0; // Borrow this variable from Phils model for now !!
731 angular_velocity_SI = 2.0 * LS_PI * RPM / 60.0;
732 forward_velocity = IAS * 0.514444444444; // Convert to m/s
734 //cout << "Gear_Ratio = " << Gear_Ratio << '\n';
735 //cout << "IAS = " << IAS << '\n';
736 //cout << "forward_velocity = " << forward_velocity << '\n';
737 //cout << "FGProp1_RPS = " << FGProp1_RPS << '\n';
738 //cout << "prop_diameter = " << prop_diameter << '\n';
742 J = forward_velocity / (FGProp1_RPS * prop_diameter);
743 //cout << "advance_ratio = " << J << '\n';
745 //Cp correlations based on data from McCormick
746 Cp_20 = 0.0342*J*J*J*J - 0.1102*J*J*J + 0.0365*J*J - 0.0133*J + 0.064;
747 Cp_25 = 0.0119*J*J*J*J - 0.0652*J*J*J + 0.018*J*J - 0.0077*J + 0.0921;
749 //cout << "Cp_20 = " << Cp_20 << '\n';
750 //cout << "Cp_25 = " << Cp_25 << '\n';
752 //Assume that the blade angle is between 20 and 25 deg - reasonable for fixed pitch prop but won't hold for variable one !!!
753 Cp = Cp_20 + ( (Cp_25 - Cp_20) * ((blade_angle - 20)/(25 - 20)) );
754 //cout << "Cp = " << Cp << '\n';
755 //cout << "RPM = " << RPM << '\n';
756 //cout << "angular_velocity_SI = " << angular_velocity_SI << '\n';
758 prop_power_consumed_SI = Cp * rho_air * pow(FGProp1_RPS,3.0) * pow(prop_diameter,5.0);
759 //cout << "prop HP consumed = " << prop_power_consumed_SI / 745.699 << '\n';
760 if(angular_velocity_SI == 0)
762 // 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
763 // Unlikely to happen in practice - but I suppose someone could move the lever of a stopped large piston engine from feathered to windmilling.
764 // This *does* give problems - if the plane is put into a steep climb whilst windmilling the engine friction will eventually stop it spinning.
765 // 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!!!
767 prop_torque = prop_power_consumed_SI / angular_velocity_SI;
769 // calculate neta_prop here
770 neta_prop_20 = 0.1328*J*J*J*J - 1.3073*J*J*J + 0.3525*J*J + 1.5591*J + 0.0007;
771 neta_prop_25 = -0.3121*J*J*J*J + 0.4234*J*J*J - 0.7686*J*J + 1.5237*J - 0.0004;
772 neta_prop = neta_prop_20 + ( (neta_prop_25 - neta_prop_20) * ((blade_angle - 20)/(25 - 20)) );
774 // Check for zero forward velocity to avoid divide by zero
775 if(forward_velocity < 0.0001)
777 // I don't see how this works - how can the plane possibly start from rest ???
778 // Hmmmm - it works because the forward_velocity at present never drops below about 0.03 even at rest
779 // We can't really rely on this in the future - needs fixing !!!!
780 // 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.
781 // It would be far more natural to work out the power consumed from the thrust - FIXME!!!!!.
783 prop_thrust = neta_prop * prop_power_consumed_SI / forward_velocity; //TODO - rename forward_velocity to IAS_SI
784 //cout << "prop_thrust = " << prop_thrust << '\n';
786 //******************************************************************************
787 // Now do the engine - prop torque balance to calculate final RPM
789 //Calculate new RPM from torque balance and inertia.
790 Torque_Imbalance = Torque_SI - prop_torque; //This gives a +ve value when the engine torque exeeds the prop torque
791 // (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)
793 angular_acceleration = Torque_Imbalance / (engine_inertia + prop_inertia);
794 angular_velocity_SI += (angular_acceleration * time_step);
795 // Don't let the engine go into reverse
796 if(angular_velocity_SI < 0)
797 angular_velocity_SI = 0;
798 RPM = (angular_velocity_SI * 60) / (2.0 * LS_PI);
803 //DCL - stall the engine if RPM drops below 500 - this is possible if mixture lever is pulled right out