-// Module: 10520c.c
-// Author: Phil Schubert
-// Date started: 12/03/99
-// Purpose: Models a Continental IO-520-M Engine
-// Called by: FGSimExec
-//
-// Copyright (C) 1999 Philip L. Schubert (philings@ozemail.com.au)
+// IO360.cxx - a piston engine model currently for the IO360 engine fitted to the C172
+// but with the potential to model other naturally aspirated piston engines
+// given appropriate config input.
+//
+// Written by David Luff, started 2000.
+// Based on code by Phil Schubert, started 1999.
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
-// 02111-1307, USA.
-//
-// Further information about the GNU General Public License can also
-// be found on the world wide web at http://www.gnu.org.
-//
-// FUNCTIONAL DESCRIPTION
-// ------------------------------------------------------------------------
-// Models a Continental IO-520-M engine. This engine is used in Cessna
-// 210, 310, Beechcraft Bonaza and Baron C55. The equations used below
-// were determined by a first and second order curve fits using Excel.
-// The data is from the Cessna Aircraft Corporations Engine and Flight
-// Computer for C310. Part Number D3500-13
-//
-// ARGUMENTS
-// ------------------------------------------------------------------------
-//
-//
-// HISTORY
-// ------------------------------------------------------------------------
-// 12/03/99 PLS Created
-// 07/03/99 PLS Added Calculation of Density, and Prop_Torque
-// 07/03/99 PLS Restructered Variables to allow easier implementation
-// of Classes
-// 15/03/99 PLS Added Oil Pressure, Oil Temperature and CH Temp
-// ------------------------------------------------------------------------
-// INCLUDES
-// ------------------------------------------------------------------------
-//
-//
-/////////////////////////////////////////////////////////////////////
-//
-// Modified by Dave Luff (david.luff@nottingham.ac.uk) September 2000
-//
-// Altered manifold pressure range to add a minimum value at idle to simulate the throttle stop / idle bypass valve,
-// and to reduce the maximum value whilst the engine is running to slightly below ambient to account for CdA losses across the throttle
-//
-// Altered it a bit to model an IO360 from C172 - 360 cubic inches, 180 HP max, fixed pitch prop
-// 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
-// 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
-//
-// Added EGT model based on combustion efficiency and an energy balance with the exhaust gases
-//
-// Added a mixture - power correlation based on a curve in the IO360 operating manual
-//
-// I've tried to match the prop and engine model to give roughly 600 RPM idle and 180 HP at 2700 RPM
-// but it is by no means currently at a completed stage - DCL 15/9/00
-//
-// DCL 28/09/00 - Added estimate of engine and prop inertia and changed engine speed calculation to be calculated from Angular acceleration = Torque / Inertia.
-// Requires a timestep to be passed to FGNewEngine::init and currently assumes this timestep does not change.
-// Could easily be altered to pass a variable timestep to FGNewEngine::update every step instead if required.
-//
-// DCL 27/10/00 - Added first stab at cylinder head temperature model
-// See the comment block in the code for details
-//
-// DCL 02/11/00 - Modified EGT code to reduce values to those more representative of a sensor downstream
-//
-// DCL 02/02/01 - Changed the prop model to one based on efficiency and co-efficient of power curves from McCormick instead of the
-// blade element method we were using previously. This works much better, and is similar to how Jon is doing it in JSBSim.
-//
-// DCL 08/02/01 - Overhauled fuel consumption rate support.
-//
-// DCL 22/03/01 - Added input of actual air pressure and temperature (and hence density) to the model. Hence the power correlation
-// with pressure height and temperature is no longer required since the power is based on the actual manifold pressure.
-//
-// 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)
-// I have reduced the sea level full rich mixture to thi = 1.3
-//
-// DCL 18/9/01 - Got the engine to start and stop in response to the magneto switch.
-// Changed all PI to LS_PI (in ls_constants.h).
-// Engine now checks for fuel and stops when not available.
-//
-//////////////////////////////////////////////////////////////////////
+// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#include <simgear/compiler.h>
#include <Main/fg_props.hxx>
-// Static utility functions
-
-// Calculate Density Ratio
-static float Density_Ratio ( float x )
-{
- float y ;
- y = ((3E-10 * x * x) - (3E-05 * x) + 0.9998);
- return(y);
-}
-
-
-// Calculate Air Density - Rho, using the ideal gas equation
-// Takes and returns SI values
-static float Density ( float temperature, float pressure )
-{
- // rho = P / RT
- // R = 287.3 for air
-
- float R = 287.3;
- float rho = pressure / (R * temperature);
- return(rho);
-}
-
-
-// Calculate Speed in FPS given Knots CAS
-static float IAS_to_FPS (float x)
-{
- float y;
- y = x * 1.68888888;
- return y;
-}
-
-// FGNewEngine member functions
-
-float FGNewEngine::Lookup_Combustion_Efficiency(float thi_actual)
-{
- const int NUM_ELEMENTS = 11;
- 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
- 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
- //combustion efficiency values from Heywood, "Internal Combustion Engine Fundamentals", ISBN 0-07-100499-8
- float neta_comb_actual = 0.0f;
- float factor;
-
- int i;
- int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays
-
- for(i=0;i<j;i++)
- {
- if(i == (j-1)) {
- // Assume linear extrapolation of the slope between the last two points beyond the last point
- float dydx = (neta_comb[i] - neta_comb[i-1]) / (thi[i] - thi[i-1]);
- neta_comb_actual = neta_comb[i] + dydx * (thi_actual - thi[i]);
- return neta_comb_actual;
- }
- if(thi_actual == thi[i]) {
- neta_comb_actual = neta_comb[i];
- return neta_comb_actual;
- }
- if((thi_actual > thi[i]) && (thi_actual < thi[i + 1])) {
- //do linear interpolation between the two points
- factor = (thi_actual - thi[i]) / (thi[i+1] - thi[i]);
- neta_comb_actual = (factor * (neta_comb[i+1] - neta_comb[i])) + neta_comb[i];
- return neta_comb_actual;
- }
- }
-
- //if we get here something has gone badly wrong
- cout << "ERROR: error in FGNewEngine::Lookup_Combustion_Efficiency\n";
- return neta_comb_actual;
-}
-
-////////////////////////////////////////////////////////////////////////////////////////////
-// Return the percentage of best mixture power available at a given mixture strength
-//
-// Based on data from "Technical Considerations for Catalysts for the European Market"
-// by H S Gandi, published 1988 by IMechE
-//
-// Note that currently no attempt is made to set a lean limit on stable combustion
-////////////////////////////////////////////////////////////////////////////////////////////
-float FGNewEngine::Power_Mixture_Correlation(float thi_actual)
-{
- float AFR_actual = 14.7 / thi_actual;
- // thi and thi_actual are equivalence ratio
- const int NUM_ELEMENTS = 13;
- // 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.
- 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
- 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
- float mixPerPow_actual = 0.0f;
- float factor;
- float dydx;
-
- int i;
- int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays
-
- for(i=0;i<j;i++)
- {
- if(i == (j-1)) {
- // Assume linear extrapolation of the slope between the last two points beyond the last point
- dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
- mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
- return mixPerPow_actual;
- }
- if((i == 0) && (AFR_actual < AFR[i])) {
- // Assume linear extrapolation of the slope between the first two points for points before the first point
- dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
- mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
- return mixPerPow_actual;
- }
- if(AFR_actual == AFR[i]) {
- mixPerPow_actual = mixPerPow[i];
- return mixPerPow_actual;
- }
- if((AFR_actual > AFR[i]) && (AFR_actual < AFR[i + 1])) {
- //do linear interpolation between the two points
- factor = (AFR_actual - AFR[i]) / (AFR[i+1] - AFR[i]);
- mixPerPow_actual = (factor * (mixPerPow[i+1] - mixPerPow[i])) + mixPerPow[i];
- return mixPerPow_actual;
- }
- }
-
- //if we get here something has gone badly wrong
- cout << "ERROR: error in FGNewEngine::Power_Mixture_Correlation\n";
- return mixPerPow_actual;
-}
-
-
-
-// Calculate Manifold Pressure based on Throttle lever Position
-float FGNewEngine::Calc_Manifold_Pressure ( float LeverPosn, float MaxMan, float MinMan)
-{
- float Inches;
- // if ( x < = 0 ) {
- // x = 0.00001;
- // }
-
- //Note that setting the manifold pressure as a function of lever position only is not strictly accurate
- //MAP is also a function of engine speed. (and ambient pressure if we are going for an actual MAP model)
- Inches = MinMan + (LeverPosn * (MaxMan - MinMan) / 100);
-
- //allow for idle bypass valve or slightly open throttle stop
- if(Inches < MinMan)
- Inches = MinMan;
-
- return Inches;
-}
-
-
-
-
-// Calculate Oil Temperature
-float FGNewEngine::Calc_Oil_Temp (float Fuel_Flow, float Mixture, float IAS)
-{
- float Oil_Temp = 85;
-
- return (Oil_Temp);
-}
-
-// Calculate Oil Pressure
-float FGNewEngine::Calc_Oil_Press (float Oil_Temp, float Engine_RPM)
-{
- float Oil_Pressure = 0; //PSI
- float Oil_Press_Relief_Valve = 60; //PSI
- float Oil_Press_RPM_Max = 1800;
- float Design_Oil_Temp = 85; //Celsius
- float Oil_Viscosity_Index = 0.25; // PSI/Deg C
- float Temp_Deviation = 0; // Deg C
-
- Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
-
- // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
- if (Oil_Pressure >= Oil_Press_Relief_Valve)
- {
- Oil_Pressure = Oil_Press_Relief_Valve;
- }
-
- // Now adjust pressure according to Temp which affects the viscosity
-
- Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
-
- return Oil_Pressure;
-}
-
//*************************************************************************************
// Initialise the engine model
void FGNewEngine::init(double dt) {
calorific_value_fuel = 47.3e6; // W/Kg Note that this is only an approximate value
rho_fuel = 800; // kg/m^3 - an estimate for now
R_air = 287.3;
- p_amb_sea_level = 101325;
- // Control and environment inputs
- IAS = 0;
+ // environment inputs
+ p_amb_sea_level = 101325; // Pascals
+
+ // Control inputs - ARE THESE NEEDED HERE???
Throttle_Lever_Pos = 75;
Propeller_Lever_Pos = 75;
Mixture_Lever_Pos = 100;
+ //misc
+ IAS = 0;
time_step = dt;
- // Engine Specific Variables.
- // Will be set in a parameter file to be read in to create
- // and instance for each engine.
- Max_Manifold_Pressure = 28.50; //Inches Hg. An approximation - should be able to find it in the engine performance data
- 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
- Max_RPM = 2700;
- Min_RPM = 600; //Recommended idle from Continental data sheet
- Max_Fuel_Flow = 130;
- Mag_Derate_Percent = 5;
-// MaxHP = 285; //Continental IO520-M
+ // Engine Specific Variables that should be read in from a config file
MaxHP = 200; //Lycoming IO360 -A-C-D series
// MaxHP = 180; //Current Lycoming IO360 ?
// displacement = 520; //Continental IO520-M
displacement_SI = displacement * CONVERT_CUBIC_INCHES_TO_METERS_CUBED;
engine_inertia = 0.2; //kgm^2 - value taken from a popular family saloon car engine - need to find an aeroengine value !!!!!
prop_inertia = 0.05; //kgm^2 - this value is a total guess - dcl
+ Max_Fuel_Flow = 130; // Units??? Do we need this variable any more??
+
+ // 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.
+ Max_Manifold_Pressure = 28.50; //Inches Hg. An approximation - should be able to find it in the engine performance data
+ 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
+ Max_RPM = 2700;
+ Min_RPM = 600; //Recommended idle from Continental data sheet
+ Mag_Derate_Percent = 5;
Gear_Ratio = 1;
- n_R = 2; // Number of crank revolutions per power cycle - 2 for a 4 cylinder engine.
+ n_R = 2; // Number of crank revolutions per power cycle - 2 for a 4 stroke engine.
+ // Various bits of housekeeping describing the engines initial state.
running = fgGetBool("/engines/engine[0]/running");
cranking = false;
+ crank_counter = false;
fgSetBool("/engines/engine[0]/cranking", false);
// Initialise Engine Variables used by this instance
else
RPM = 0;
Percentage_Power = 0;
- Manifold_Pressure = 29.00; // Inches
+ Manifold_Pressure = 29.96; // Inches
Fuel_Flow_gals_hr = 0;
- Torque = 0;
+// Torque = 0;
Torque_SI = 0;
- CHT = 298.0; //deg Kelvin
- CHT_degF = (CHT * 1.8) - 459.67; //deg Fahrenheit
+ CHT_degK = 298.0; //deg Kelvin
+ CHT_degF = (CHT_degF * 1.8) - 459.67; //deg Fahrenheit
Mixture = 14;
Oil_Pressure = 0; // PSI
Oil_Temp = 85; // Deg C
+ current_oil_temp = 298.0; //deg Kelvin
+ /**** one of these is superfluous !!!!***/
HP = 0;
RPS = 0;
Torque_Imbalance = 0;
// Initialise Propellor Variables used by this instance
- FGProp1_Thrust = 0;
FGProp1_RPS = 0;
- FGProp1_Blade_Angle = 13.5;
+ // Hardcode propellor for now - the following two should be read from config eventually
prop_diameter = 1.8; // meters
blade_angle = 23.0; // degrees
}
-
-//*****************************************************************************
//*****************************************************************************
// update the engine model based on current control positions
void FGNewEngine::update() {
// cout << "Mixture lever = " << Mixture_Lever_Pos << '\n';
// cout << "n = " << RPM << " rpm\n";
// cout << "T_amb = " << T_amb << '\n';
- cout << "running = " << running << '\n';
- cout << "fuel = " << fgGetFloat("/consumables/fuel/tank[0]/level-gal_us") << '\n';
+// cout << "running = " << running << '\n';
+// cout << "fuel = " << fgGetFloat("/consumables/fuel/tank[0]/level-gal_us") << '\n';
+ cout << "Percentage_Power = " << Percentage_Power << '\n';
+ cout << "current_oil_temp = " << current_oil_temp << '\n';
}
count1++;
if(count1 == 600)
count1 = 0;
*/
- float ManXRPM = 0;
- float Vo = 0;
- float V1 = 0;
+ // Check parameters that may alter the operating state of the engine.
+ // (spark, fuel, starter motor etc)
- // Parameters that alter the operation of the engine. (spark, fuel, starter motor etc)
// Check for spark
- int Magneto_Left = 0;
- int Magneto_Right = 0;
+ bool Magneto_Left = false;
+ bool Magneto_Right = false;
int mag_pos = fgGetInt("/engines/engine[0]/magneto");
// Magneto positions:
// 0 -> off
spark = false;
} // neglects battery voltage, master on switch, etc for now.
if((mag_pos == 1) || (mag_pos > 2))
- Magneto_Left = 1;
+ Magneto_Left = true;
if(mag_pos > 1)
- Magneto_Right = 1;
+ Magneto_Right = true;
// crude check for fuel
if((fgGetFloat("/consumables/fuel/tank[0]/level-gal_us") > 0) || (fgGetFloat("/consumables/fuel/tank[1]/level-gal_us") > 0)) {
fgSetBool("/engines/engine[0]/cranking", true);
else
fgSetBool("/engines/engine[0]/cranking", false);
+ crank_counter = 0;
}
// Note that although /engines/engine[0]/starter and /engines/engine[0]/cranking might appear to be duplication it is
// not since the starter may be engaged with the battery voltage too low for cranking to occur (or perhaps the master
// switch just left off) and the sound manager will read .../cranking to determine wether to play a cranking sound.
// For now though none of that is implemented so cranking can be set equal to .../starter without further checks.
- int Alternate_Air_Pos =0; // Off = 0. Reduces power by 3 % for same throttle setting
+// int Alternate_Air_Pos =0; // Off = 0. Reduces power by 3 % for same throttle setting
// DCL - don't know what this Alternate_Air_Pos is - this is a leftover from the Schubert code.
//Check mode of engine operation
if(cranking) {
+ crank_counter++;
if(RPM <= 480) {
RPM += 100;
if(RPM > 480)
// consider making a horrible noise if the starter is engaged with the engine running
}
}
- if((!running) && (spark) && (fuel)) {
+ if((!running) && (spark) && (fuel) && (crank_counter > 120)) {
// start the engine if revs high enough
if(RPM > 450) {
// For now just instantaneously start but later we should maybe crank for a bit
running = true;
fgSetBool("/engines/engine[0]/running", true);
- RPM = 600;
+// RPM = 600;
}
}
if( (running) && ((!spark)||(!fuel)) ) {
fgSetBool("/engines/engine[0]/running", false);
}
+ // Now we've ascertained whether the engine is running or not we can start to do the engine calculations 'proper'
+
// Calculate Sea Level Manifold Pressure
Manifold_Pressure = Calc_Manifold_Pressure( Throttle_Lever_Pos, Max_Manifold_Pressure, Min_Manifold_Pressure );
// cout << "manifold pressure = " << Manifold_Pressure << endl;
//Then find the actual manifold pressure (the calculated one is the sea level pressure)
True_Manifold_Pressure = Manifold_Pressure * p_amb / p_amb_sea_level;
-//*************
-//DCL - next calculate m_dot_air and m_dot_fuel into engine
-
- //DCL - calculate mass air flow into engine based on speed and load - separate this out into a function eventually
- //t_amb is actual temperature calculated from altitude
- //calculate density from ideal gas equation
- rho_air = p_amb / ( R_air * T_amb );
- 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.
- //calculate ideal engine volume inducted per second
- swept_volume = (displacement_SI * (RPM / 60)) / 2; //This equation is only valid for a four stroke engine
- //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
- //Note that this is cylinder vol eff - the throttle drop is already accounted for in the MAP calculation
- volumetric_efficiency = 0.8;
- //Now use volumetric efficiency to calculate actual air volume inducted per second
- v_dot_air = swept_volume * volumetric_efficiency;
- //Now calculate mass flow rate of air into engine
- m_dot_air = v_dot_air * rho_air_manifold;
-
-//**************
-
- //DCL - now calculate fuel flow into engine based on air flow and mixture lever position
- //assume lever runs from no flow at fully out to thi = 1.3 at fully in at sea level
- //also assume that the injector linkage is ideal - hence the set mixture is maintained at a given altitude throughout the speed and load range
- thi_sea_level = 1.3 * ( Mixture_Lever_Pos / 100.0 );
- 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
- m_dot_fuel = m_dot_air / 14.7 * equivalence_ratio;
- Fuel_Flow_gals_hr = (m_dot_fuel / rho_fuel) * 264.172 * 3600.0; // Note this assumes US gallons
-
-//***********************************************************************
-//Engine power and torque calculations
+ //Do the fuel flow calculations
+ Calc_Fuel_Flow_Gals_Hr();
- // For a given Manifold Pressure and RPM calculate the % Power
- // Multiply Manifold Pressure by RPM
- ManXRPM = True_Manifold_Pressure * RPM;
-// ManXRPM = Manifold_Pressure * RPM;
- // cout << ManXRPM;
- // cout << endl;
+ //Calculate engine power
+ Calc_Percentage_Power(Magneto_Left, Magneto_Right);
+ HP = Percentage_Power * MaxHP / 100.0;
+ Power_SI = HP * CONVERT_HP_TO_WATTS;
-/*
-// Phil's %power correlation
- // Calculate % Power
- Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM) + ( + 7E-04 * ManXRPM) - 0.1218;
- // cout << Percentage_Power << "%" << "\t";
-*/
+ // FMEP calculation. For now we will just use this during motored operation.
+ // Eventually we will calculate IMEP and use the FMEP all the time to give BMEP (maybe!)
+ if(!running) {
+ // This FMEP data is from the Patton paper, assumes fully warm conditions.
+ FMEP = 1e-12*pow(RPM,4) - 1e-8*pow(RPM,3) + 5e-5*pow(RPM,2) - 0.0722*RPM + 154.85;
+ // Gives FMEP in kPa - now convert to Pa
+ FMEP *= 1000;
+ } else {
+ FMEP = 0.0;
+ }
+ // Is this total FMEP or friction FMEP ???
-// DCL %power correlation - basically Phil's correlation modified to give slighty less power at the low end
-// might need some adjustment as the prop model is adjusted
-// 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
- // Calculate % Power for Nev's prop model
- //Percentage_Power = (+ 6E-09 * ManXRPM * ManXRPM) + ( + 8E-04 * ManXRPM) - 1.8524;
+ Torque_FMEP = (FMEP * displacement_SI) / (2.0 * LS_PI * n_R);
- // Calculate %power for DCL prop model
- Percentage_Power = (7e-9 * ManXRPM * ManXRPM) + (7e-4 * ManXRPM) - 1.0;
+ // Calculate Engine Torque. Check for div by zero since percentage power correlation does not guarantee zero power at zero rpm.
+ // However this is problematical since there is a resistance to movement even at rest
+ // Ie this is a dynamics equation not a statics one. This can be solved by going over to MEP based torque calculations.
+ if(RPM == 0) {
+ Torque_SI = 0 - Torque_FMEP;
+ }
+ else {
+ Torque_SI = ((Power_SI * 60.0) / (2.0 * LS_PI * RPM)) - Torque_FMEP; //Torque = power / angular velocity
+ // cout << Torque << " Nm\n";
+ }
- // Power de-rating for altitude has been removed now that we are basing the power
- // on the actual manifold pressure, which takes air pressure into account. However - this fails to
- // take the temperature into account - this is TODO.
+ //Calculate Exhaust gas temperature
+ Calc_EGT();
- // Adjust power for temperature - this is temporary until the power is done as a function of mass flow rate induced
- // Adjust for Temperature - Temperature above Standard decrease
- // power by 7/120 % per degree F increase, and incease power for
- // temps below at the same ratio
- float T_amb_degF = (T_amb * 1.8) - 459.67;
- float T_amb_sea_lev_degF = (288 * 1.8) - 459.67;
- Percentage_Power = Percentage_Power + ((T_amb_sea_lev_degF - T_amb_degF) * 7 /120);
-
- //DCL - now adjust power to compensate for mixture
- Percentage_of_best_power_mixture_power = Power_Mixture_Correlation(equivalence_ratio);
- Percentage_Power = Percentage_Power * Percentage_of_best_power_mixture_power / 100.0;
+ // Calculate Cylinder Head Temperature
+ CHT_degK = Calc_CHT(CHT_degK);
+ CHT_degF = (CHT_degK * 1.8) - 459.67;
+
+ // Calculate oil temperature
+ current_oil_temp = Calc_Oil_Temp(current_oil_temp);
+
+ // Calculate Oil Pressure
+ Oil_Pressure = Calc_Oil_Press( Oil_Temp, RPM );
+
+ // Now do the Propeller Calculations
+ Do_Prop_Calcs();
+
+// Now do the engine - prop torque balance to calculate final RPM
+
+ //Calculate new RPM from torque balance and inertia.
+ Torque_Imbalance = Torque_SI - prop_torque; //This gives a +ve value when the engine torque exeeds the prop torque
+ // (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)
+
+ angular_acceleration = Torque_Imbalance / (engine_inertia + prop_inertia);
+ angular_velocity_SI += (angular_acceleration * time_step);
+ // Don't let the engine go into reverse
+ if(angular_velocity_SI < 0)
+ angular_velocity_SI = 0;
+ RPM = (angular_velocity_SI * 60) / (2.0 * LS_PI);
- // Now Derate engine for the effects of Bad/Switched off magnetos
- //if (Magneto_Left == 0 && Magneto_Right == 0) {
- if(!running) {
- // cout << "Both OFF\n";
- Percentage_Power = 0;
- } else if (Magneto_Left && Magneto_Right) {
- // cout << "Both On ";
- } else if (Magneto_Left == 0 || Magneto_Right== 0) {
- // cout << "1 Magneto Failed ";
- Percentage_Power = Percentage_Power * ((100.0 - Mag_Derate_Percent)/100.0);
- // cout << FGEng1_Percentage_Power << "%" << "\t";
+ // And finally a last check on the engine state after doing the torque balance with the prop - have we stalled?
+ if(running) {
+ //Check if we have stalled the engine
+ if (RPM == 0) {
+ running = false;
+ fgSetBool("/engines/engine[0]/running", false);
+ } else if((RPM <= 480) && (cranking)) {
+ //Make sure the engine noise dosn't play if the engine won't start due to eg mixture lever pulled out.
+ running = false;
+ fgSetBool("/engines/engine[0]/running", false);
+ }
}
+}
- //DCL - stall the engine if RPM drops below 450 - this is possible if mixture lever is pulled right out
- //This is a kludge that I should eliminate by adding a total fmep estimation.
- if(RPM < 450)
- Percentage_Power = 0;
-
- if(Percentage_Power < 0)
- Percentage_Power = 0;
-
- // FMEP calculation. For now we will just use this during motored operation, ie when %Power == 0.
- // Eventually we will calculate IMEP and use the FMEP all the time to give BMEP
- //
- if(Percentage_Power == 0) {
- // This FMEP data is from the Patton paper, assumes fully warm conditions.
- FMEP = 1e-12*pow(RPM,4) - 1e-8*pow(RPM,3) + 5e-5*pow(RPM,2) - 0.0722*RPM + 154.85;
- // Gives FMEP in kPa - now convert to Pa
- FMEP *= 1000;
- } else {
- FMEP = 0.0;
- }
+//*****************************************************************************************************
- Torque_FMEP = (FMEP * displacement_SI) / (2.0 * LS_PI * n_R);
+// FGNewEngine member functions
- HP = Percentage_Power * MaxHP / 100.0;
+////////////////////////////////////////////////////////////////////////////////////////////
+// Return the combustion efficiency as a function of equivalence ratio
+//
+// Combustion efficiency values from Heywood,
+// "Internal Combustion Engine Fundamentals", ISBN 0-07-100499-8
+////////////////////////////////////////////////////////////////////////////////////////////
+float FGNewEngine::Lookup_Combustion_Efficiency(float thi_actual)
+{
+ const int NUM_ELEMENTS = 11;
+ 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
+ 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
+ float neta_comb_actual = 0.0f;
+ float factor;
- Power_SI = HP * CONVERT_HP_TO_WATTS;
+ int i;
+ int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays
- // Calculate Engine Torque. Check for div by zero since percentage power correlation does not guarantee zero power at zero rpm.
- // However this is problematical since there is a resistance to movement even at rest
- // Ie this is a dynamics equation not a statics one. This can be solved by going over to MEP based torque calculations.
- if(RPM == 0) {
- Torque_SI = 0 - Torque_FMEP;
- }
- else {
- Torque_SI = ((Power_SI * 60.0) / (2.0 * LS_PI * RPM)) - Torque_FMEP; //Torque = power / angular velocity
- // cout << Torque << " Nm\n";
+ for(i=0;i<j;i++)
+ {
+ if(i == (j-1)) {
+ // Assume linear extrapolation of the slope between the last two points beyond the last point
+ float dydx = (neta_comb[i] - neta_comb[i-1]) / (thi[i] - thi[i-1]);
+ neta_comb_actual = neta_comb[i] + dydx * (thi_actual - thi[i]);
+ return neta_comb_actual;
+ }
+ if(thi_actual == thi[i]) {
+ neta_comb_actual = neta_comb[i];
+ return neta_comb_actual;
+ }
+ if((thi_actual > thi[i]) && (thi_actual < thi[i + 1])) {
+ //do linear interpolation between the two points
+ factor = (thi_actual - thi[i]) / (thi[i+1] - thi[i]);
+ neta_comb_actual = (factor * (neta_comb[i+1] - neta_comb[i])) + neta_comb[i];
+ return neta_comb_actual;
+ }
}
-//**********************************************************************
-//Calculate Exhaust gas temperature
-
- // cout << "Thi = " << equivalence_ratio << '\n';
-
- combustion_efficiency = Lookup_Combustion_Efficiency(equivalence_ratio); //The combustion efficiency basically tells us what proportion of the fuels calorific value is released
-
- // cout << "Combustion efficiency = " << combustion_efficiency << '\n';
-
- //now calculate energy release to exhaust
- //We will assume a three way split of fuel energy between useful work, the coolant system and the exhaust system
- //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
- //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.
- enthalpy_exhaust = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
- heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
- delta_T_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
-// delta_T_exhaust = Calculate_Delta_T_Exhaust();
-
- // cout << "T_amb " << T_amb;
- // cout << " dT exhaust = " << delta_T_exhaust;
-
- EGT = T_amb + delta_T_exhaust;
+ //if we get here something has gone badly wrong
+ cout << "ERROR: error in FGNewEngine::Lookup_Combustion_Efficiency\n";
+ return neta_comb_actual;
+}
- //The above gives the exhaust temperature immediately prior to leaving the combustion chamber
- //Now derate to give a more realistic figure as measured downstream
- //For now we will aim for a peak of around 400 degC (750 degF)
+////////////////////////////////////////////////////////////////////////////////////////////
+// Return the percentage of best mixture power available at a given mixture strength
+//
+// Based on data from "Technical Considerations for Catalysts for the European Market"
+// by H S Gandi, published 1988 by IMechE
+//
+// Note that currently no attempt is made to set a lean limit on stable combustion
+////////////////////////////////////////////////////////////////////////////////////////////
+float FGNewEngine::Power_Mixture_Correlation(float thi_actual)
+{
+ float AFR_actual = 14.7 / thi_actual;
+ // thi and thi_actual are equivalence ratio
+ const int NUM_ELEMENTS = 13;
+ // 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.
+ 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
+ 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
+ float mixPerPow_actual = 0.0f;
+ float factor;
+ float dydx;
- EGT *= 0.444 + ((0.544 - 0.444) * Percentage_Power / 100.0);
+ int i;
+ int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays
- EGT_degF = (EGT * 1.8) - 459.67;
+ for(i=0;i<j;i++)
+ {
+ if(i == (j-1)) {
+ // Assume linear extrapolation of the slope between the last two points beyond the last point
+ dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
+ mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
+ return mixPerPow_actual;
+ }
+ if((i == 0) && (AFR_actual < AFR[i])) {
+ // Assume linear extrapolation of the slope between the first two points for points before the first point
+ dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
+ mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
+ return mixPerPow_actual;
+ }
+ if(AFR_actual == AFR[i]) {
+ mixPerPow_actual = mixPerPow[i];
+ return mixPerPow_actual;
+ }
+ if((AFR_actual > AFR[i]) && (AFR_actual < AFR[i + 1])) {
+ //do linear interpolation between the two points
+ factor = (AFR_actual - AFR[i]) / (AFR[i+1] - AFR[i]);
+ mixPerPow_actual = (factor * (mixPerPow[i+1] - mixPerPow[i])) + mixPerPow[i];
+ return mixPerPow_actual;
+ }
+ }
- //cout << " EGT = " << EGT << " degK " << EGT_degF << " degF";// << '\n';
+ //if we get here something has gone badly wrong
+ cout << "ERROR: error in FGNewEngine::Power_Mixture_Correlation\n";
+ return mixPerPow_actual;
+}
-//***************************************************************************************
// Calculate Cylinder Head Temperature
-
-/* DCL 27/10/00
-
-This is a somewhat rough first attempt at modelling cylinder head temperature. The cylinder head
-is assumed to be at uniform temperature. Obviously this is incorrect, but it simplifies things a
-lot, and we're just looking for the behaviour of CHT to be correct. Energy transfer to the cylinder
-head is assumed to be one third of the energy released by combustion at all conditions. This is a
-reasonable estimate, although obviously in real life it varies with different conditions and possibly
-with CHT itself. I've split energy transfer from the cylinder head into 2 terms - free convection -
-ie convection to stationary air, and forced convection, ie convection into flowing air. The basic
-free convection equation is: dqdt = -hAdT Since we don't know A and are going to set h quite arbitarily
-anyway I've knocked A out and just wrapped it up in h - the only real significance is that the units
-of h will be different but that dosn't really matter to us anyway. In addition, we have the problem
-that the prop model I'm currently using dosn't model the backwash from the prop which will add to the
-velocity of the cooling air when the prop is turning, so I've added an extra term to try and cope
-with this.
-
-In real life, forced convection equations are genarally empirically derived, and are quite complicated
-and generally contain such things as the Reynolds and Nusselt numbers to various powers. The best
-course of action would probably to find an empirical correlation from the literature for a similar
-situation and try and get it to fit well. However, for now I am using my own made up very simple
-correlation for the energy transfer from the cylinder head:
-
-dqdt = -(h1.dT) -(h2.m_dot.dT) -(h3.rpm.dT)
-
-where dT is the temperature different between the cylinder head and the surrounding air, m_dot is the
-mass flow rate of cooling air through an arbitary volume, rpm is the engine speed in rpm (this is the
-backwash term), and h1, h2, h3 are co-efficients which we can play with to attempt to get the CHT
-behaviour to match real life.
-
-In order to change the values of CHT that the engine settles down at at various conditions,
-have a play with h1, h2 and h3. In order to change the rate of heating/cooling without affecting
-equilibrium values alter the cylinder head mass, which is really quite arbitary. Bear in mind that
-altering h1, h2 and h3 will also alter the rate of heating or cooling as well as equilibrium values,
-but altering the cylinder head mass will only alter the rate. It would I suppose be better to read
-the values from file to avoid the necessity for re-compilation every time I change them.
-
-*/
- //CHT = Calc_CHT( Fuel_Flow, Mixture, IAS);
- // cout << "Cylinder Head Temp (F) = " << CHT << endl;
+// Crudely models the cylinder head as an arbitary lump of arbitary size and area with one third of combustion energy
+// as heat input and heat output as a function of airspeed and temperature. Could be improved!!!
+float FGNewEngine::Calc_CHT(float CHT)
+{
float h1 = -95.0; //co-efficient for free convection
float h2 = -3.95; //co-efficient for forced convection
float h3 = -0.05; //co-efficient for forced convection due to prop backwash
float HeatCapacityCylinderHead;
float dCHTdt;
+ // The above values are hardwired to give reasonable results for an IO360 (C172 engine)
+ // Now adjust to try to compensate for arbitary sized engines - this is currently untested
+ arbitary_area *= (MaxHP / 180.0);
+ MassCylinderHead *= (MaxHP / 180.0);
+
temperature_difference = CHT - T_amb;
v_apparent = IAS * 0.5144444; //convert from knots to m/s
dCHTdt = dqdt_cylinder_head / HeatCapacityCylinderHead;
CHT += (dCHTdt * time_step);
-
- CHT_degF = (CHT * 1.8) - 459.67;
-
- //cout << " CHT = " << CHT_degF << " degF\n";
-
-
- // End calculate Cylinder Head Temperature
-
-
-//***************************************************************************************
-// Oil pressure calculation
-
- // Calculate Oil Pressure
- Oil_Pressure = Calc_Oil_Press( Oil_Temp, RPM );
- // cout << "Oil Pressure (PSI) = " << Oil_Pressure << endl;
-
-//**************************************************************************************
-// Now do the Propeller Calculations
-
- Gear_Ratio = 1.0;
- FGProp1_RPS = RPM * Gear_Ratio / 60.0; // Borrow this variable from Phils model for now !!
+
+ return(CHT);
+}
+
+// Calculate exhaust gas temperature
+void FGNewEngine::Calc_EGT()
+{
+ combustion_efficiency = Lookup_Combustion_Efficiency(equivalence_ratio); //The combustion efficiency basically tells us what proportion of the fuels calorific value is released
+
+ //now calculate energy release to exhaust
+ //We will assume a three way split of fuel energy between useful work, the coolant system and the exhaust system
+ //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
+ //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.
+ enthalpy_exhaust = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
+ heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
+ delta_T_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
+
+ EGT = T_amb + delta_T_exhaust;
+
+ //The above gives the exhaust temperature immediately prior to leaving the combustion chamber
+ //Now derate to give a more realistic figure as measured downstream
+ //For now we will aim for a peak of around 400 degC (750 degF)
+
+ EGT *= 0.444 + ((0.544 - 0.444) * Percentage_Power / 100.0);
+
+ EGT_degF = (EGT * 1.8) - 459.67;
+}
+
+// Calculate Manifold Pressure based on Throttle lever Position
+float FGNewEngine::Calc_Manifold_Pressure ( float LeverPosn, float MaxMan, float MinMan)
+{
+ float Inches;
+
+ //Note that setting the manifold pressure as a function of lever position only is not strictly accurate
+ //MAP is also a function of engine speed. (and ambient pressure if we are going for an actual MAP model)
+ Inches = MinMan + (LeverPosn * (MaxMan - MinMan) / 100);
+
+ //allow for idle bypass valve or slightly open throttle stop
+ if(Inches < MinMan)
+ Inches = MinMan;
+
+ //Check for stopped engine (crudest way of compensating for engine speed)
+ if(RPM == 0)
+ Inches = 29.92;
+
+ return Inches;
+}
+
+// Calculate fuel flow in gals/hr
+void FGNewEngine::Calc_Fuel_Flow_Gals_Hr()
+{
+ //DCL - calculate mass air flow into engine based on speed and load - separate this out into a function eventually
+ //t_amb is actual temperature calculated from altitude
+ //calculate density from ideal gas equation
+ rho_air = p_amb / ( R_air * T_amb );
+ 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.
+ //calculate ideal engine volume inducted per second
+ swept_volume = (displacement_SI * (RPM / 60)) / 2; //This equation is only valid for a four stroke engine
+ //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
+ //Note that this is cylinder vol eff - the throttle drop is already accounted for in the MAP calculation
+ volumetric_efficiency = 0.8;
+ //Now use volumetric efficiency to calculate actual air volume inducted per second
+ v_dot_air = swept_volume * volumetric_efficiency;
+ //Now calculate mass flow rate of air into engine
+ m_dot_air = v_dot_air * rho_air_manifold;
+
+//**************
+
+ //DCL - now calculate fuel flow into engine based on air flow and mixture lever position
+ //assume lever runs from no flow at fully out to thi = 1.3 at fully in at sea level
+ //also assume that the injector linkage is ideal - hence the set mixture is maintained at a given altitude throughout the speed and load range
+ thi_sea_level = 1.3 * ( Mixture_Lever_Pos / 100.0 );
+ 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
+ m_dot_fuel = m_dot_air / 14.7 * equivalence_ratio;
+ Fuel_Flow_gals_hr = (m_dot_fuel / rho_fuel) * 264.172 * 3600.0; // Note this assumes US gallons
+}
+
+// Calculate the percentage of maximum rated power delivered as a function of Manifold pressure multiplied by engine speed (rpm)
+// This is not necessarilly the best approach but seems to work for now.
+// May well need tweaking at the bottom end if the prop model is changed.
+void FGNewEngine::Calc_Percentage_Power(bool mag_left, bool mag_right)
+{
+ // For a given Manifold Pressure and RPM calculate the % Power
+ // Multiply Manifold Pressure by RPM
+ float ManXRPM = True_Manifold_Pressure * RPM;
+
+/*
+// Phil's %power correlation
+ // Calculate % Power
+ Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM) + ( + 7E-04 * ManXRPM) - 0.1218;
+ // cout << Percentage_Power << "%" << "\t";
+*/
+
+// DCL %power correlation - basically Phil's correlation modified to give slighty less power at the low end
+// might need some adjustment as the prop model is adjusted
+// 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
+ // Calculate % Power for Nev's prop model
+ //Percentage_Power = (+ 6E-09 * ManXRPM * ManXRPM) + ( + 8E-04 * ManXRPM) - 1.8524;
+
+ // Calculate %power for DCL prop model
+ Percentage_Power = (7e-9 * ManXRPM * ManXRPM) + (7e-4 * ManXRPM) - 1.0;
+
+ // Power de-rating for altitude has been removed now that we are basing the power
+ // on the actual manifold pressure, which takes air pressure into account. However - this fails to
+ // take the temperature into account - this is TODO.
+
+ // Adjust power for temperature - this is temporary until the power is done as a function of mass flow rate induced
+ // Adjust for Temperature - Temperature above Standard decrease
+ // power by 7/120 % per degree F increase, and incease power for
+ // temps below at the same ratio
+ float T_amb_degF = (T_amb * 1.8) - 459.67;
+ float T_amb_sea_lev_degF = (288 * 1.8) - 459.67;
+ Percentage_Power = Percentage_Power + ((T_amb_sea_lev_degF - T_amb_degF) * 7 /120);
+
+ //DCL - now adjust power to compensate for mixture
+ Percentage_of_best_power_mixture_power = Power_Mixture_Correlation(equivalence_ratio);
+ Percentage_Power = Percentage_Power * Percentage_of_best_power_mixture_power / 100.0;
+
+ // Now Derate engine for the effects of Bad/Switched off magnetos
+ //if (Magneto_Left == 0 && Magneto_Right == 0) {
+ if(!running) {
+ // cout << "Both OFF\n";
+ Percentage_Power = 0;
+ } else if (mag_left && mag_right) {
+ // cout << "Both On ";
+ } else if (mag_left == 0 || mag_right== 0) {
+ // cout << "1 Magneto Failed ";
+ Percentage_Power = Percentage_Power * ((100.0 - Mag_Derate_Percent)/100.0);
+ // cout << FGEng1_Percentage_Power << "%" << "\t";
+ }
+/*
+ //DCL - stall the engine if RPM drops below 450 - this is possible if mixture lever is pulled right out
+ //This is a kludge that I should eliminate by adding a total fmep estimation.
+ if(RPM < 450)
+ Percentage_Power = 0;
+*/
+ if(Percentage_Power < 0)
+ Percentage_Power = 0;
+}
+
+// Calculate Oil Temperature in degrees Kelvin
+float FGNewEngine::Calc_Oil_Temp (float oil_temp)
+{
+ float idle_percentage_power = 2.3; // approximately
+ float target_oil_temp; // Steady state oil temp at the current engine conditions
+ float time_constant; // The time constant for the differential equation
+ if(running) {
+ target_oil_temp = 363;
+ time_constant = 500; // Time constant for engine-on idling.
+ if(Percentage_Power > idle_percentage_power) {
+ time_constant /= ((Percentage_Power / idle_percentage_power) / 10.0); // adjust for power
+ }
+ } else {
+ target_oil_temp = 298;
+ time_constant = 1000; // Time constant for engine-off; reflects the fact that oil is no longer getting circulated
+ }
+
+ float dOilTempdt = (target_oil_temp - oil_temp) / time_constant;
+
+ oil_temp += (dOilTempdt * time_step);
+
+ return (oil_temp);
+}
+
+// Calculate Oil Pressure
+float FGNewEngine::Calc_Oil_Press (float Oil_Temp, float Engine_RPM)
+{
+ float Oil_Pressure = 0; //PSI
+ float Oil_Press_Relief_Valve = 60; //PSI
+ float Oil_Press_RPM_Max = 1800;
+ float Design_Oil_Temp = 85; //Celsius
+ float Oil_Viscosity_Index = 0.25; // PSI/Deg C
+// float Temp_Deviation = 0; // Deg C
+
+ Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
+
+ // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
+ if (Oil_Pressure >= Oil_Press_Relief_Valve) {
+ Oil_Pressure = Oil_Press_Relief_Valve;
+ }
+
+ // Now adjust pressure according to Temp which affects the viscosity
+
+ Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
+
+ return Oil_Pressure;
+}
+
+
+// Propeller calculations.
+void FGNewEngine::Do_Prop_Calcs()
+{
+ float Gear_Ratio = 1.0;
+ float blade_length; // meters
+ float forward_velocity; // m/s
+ float prop_power_consumed_SI; // Watts
+ float prop_power_consumed_HP; // HP
+ double J; // advance ratio - dimensionless
+ double Cp_20; // coefficient of power for 20 degree blade angle
+ double Cp_25; // coefficient of power for 25 degree blade angle
+ double Cp; // Our actual coefficient of power
+ double neta_prop_20;
+ double neta_prop_25;
+ double neta_prop; // prop efficiency
+
+ FGProp1_RPS = RPM * Gear_Ratio / 60.0;
angular_velocity_SI = 2.0 * LS_PI * RPM / 60.0;
forward_velocity = IAS * 0.514444444444; // Convert to m/s
- //cout << "Gear_Ratio = " << Gear_Ratio << '\n';
- //cout << "IAS = " << IAS << '\n';
- //cout << "forward_velocity = " << forward_velocity << '\n';
- //cout << "FGProp1_RPS = " << FGProp1_RPS << '\n';
- //cout << "prop_diameter = " << prop_diameter << '\n';
if(FGProp1_RPS == 0)
J = 0;
else
else
prop_thrust = neta_prop * prop_power_consumed_SI / forward_velocity; //TODO - rename forward_velocity to IAS_SI
//cout << "prop_thrust = " << prop_thrust << '\n';
-
-//******************************************************************************
-// Now do the engine - prop torque balance to calculate final RPM
-
- //Calculate new RPM from torque balance and inertia.
- Torque_Imbalance = Torque_SI - prop_torque; //This gives a +ve value when the engine torque exeeds the prop torque
- // (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)
-
- angular_acceleration = Torque_Imbalance / (engine_inertia + prop_inertia);
- angular_velocity_SI += (angular_acceleration * time_step);
- // Don't let the engine go into reverse
- if(angular_velocity_SI < 0)
- angular_velocity_SI = 0;
- RPM = (angular_velocity_SI * 60) / (2.0 * LS_PI);
-
-// if(RPM < 0)
-// RPM = 0;
-
- //DCL - stall the engine if RPM drops below 500 - this is possible if mixture lever is pulled right out
-// if(RPM < 500)
-// RPM = 0;
-
}
-// Module: 10520c.c
-// Author: Phil Schubert
-// Date started: 12/03/99
-// Purpose: Models a Continental IO-520-M Engine
-// Called by: FGSimExec
-//
-// Copyright (C) 1999 Philip L. Schubert (philings@ozemail.com.au)
+// IO360.hxx - a piston engine model currently for the IO360 engine fitted to the C172
+// but with the potential to model other naturally aspirated piston engines
+// given appropriate config input.
+//
+// Written by David Luff, started 2000.
+// Based on code by Phil Schubert, started 1999.
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
-// 02111-1307, USA.
-//
-// Further information about the GNU General Public License can also
-// be found on the world wide web at http://www.gnu.org.
-//
-// FUNCTIONAL DESCRIPTION
-// ------------------------------------------------------------------------
-// Models a Continental IO-520-M engine. This engine is used in Cessna
-// 210, 310, Beechcraft Bonaza and Baron C55. The equations used below
-// were determined by a first and second order curve fits using Excel.
-// The data is from the Cessna Aircraft Corporations Engine and Flight
-// Computer for C310. Part Number D3500-13
-//
-// ARGUMENTS
-// ------------------------------------------------------------------------
-//
-//
-// HISTORY
-// ------------------------------------------------------------------------
-// 12/03/99 PLS Created
-// 07/03/99 PLS Added Calculation of Density, and Prop_Torque
-// 07/03/99 PLS Restructered Variables to allow easier implementation
-// of Classes
-// 15/03/99 PLS Added Oil Pressure, Oil Temperature and CH Temp
-// ------------------------------------------------------------------------
-// INCLUDES
-// ------------------------------------------------------------------------
+// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
#ifndef _IO360_HXX_
#define _IO360_HXX_
private:
- float CONVERT_HP_TO_WATTS;
+ // These constants should probably be moved eventually
float CONVERT_CUBIC_INCHES_TO_METERS_CUBED;
+ float CONVERT_HP_TO_WATTS;
+
+ // Properties of working fluids
+ float Cp_air; // J/KgK
+ float Cp_fuel; // J/KgK
+ float calorific_value_fuel; // W/Kg
+ float rho_fuel; // kg/m^3
+ float rho_air; // kg/m^3
- // Control and environment inputs
+ // environment inputs
+ float p_amb_sea_level; // Sea level ambient pressure in Pascals
+ float p_amb; // Ambient pressure at current altitude in Pascals
+ float T_amb; // ditto deg Kelvin
+
+ // Control inputs
+ float Throttle_Lever_Pos; // 0 = Closed, 100 = Fully Open
+ float Propeller_Lever_Pos; // 0 = Full Course 100 = Full Fine
+ float Mixture_Lever_Pos; // 0 = Idle Cut Off 100 = Full Rich
+
+ //misc
float IAS;
- // 0 = Closed, 100 = Fully Open
- float Throttle_Lever_Pos;
- // 0 = Full Course 100 = Full Fine
- float Propeller_Lever_Pos;
- // 0 = Idle Cut Off 100 = Full Rich
- float Mixture_Lever_Pos;
-
- // Engine Specific Variables used by this program that have limits.
- // Will be set in a parameter file to be read in to create
- // and instance for each engine.
- float Max_Manifold_Pressure; //will be lower than ambient pressure for a non turbo/super charged engine due to losses through the throttle. This is the sea level full throttle value.
- float Min_Manifold_Pressure; //Closed throttle valueat idle - governed by the idle bypass valve
- float Max_RPM;
- float Min_RPM;
- float Max_Fuel_Flow;
- float Mag_Derate_Percent;
- float MaxHP;
- float Gear_Ratio;
-
- // Initialise Engine Variables used by this instance
+ double time_step;
+
+ // Engine Specific Variables that should be read in from a config file
+ float MaxHP; // Horsepower
+ float displacement; // Cubic inches
+ float displacement_SI; //m^3 (derived from above rather than read in)
+ float engine_inertia; //kg.m^2
+ float prop_inertia; //kg.m^2
+ float Max_Fuel_Flow; // Units??? Do we need this variable any more??
+
+ // 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.
+ float Max_Manifold_Pressure; // inches Hg - typical manifold pressure at full throttle and typical max rpm
+ float Min_Manifold_Pressure; // inches Hg - typical manifold pressure at idle (closed throttle)
+ float Max_RPM; // rpm - this is really a bit of a hack and could be make redundant if the prop model works properly and takes tips at speed of sound into account.
+ float Min_RPM; // rpm - possibly redundant ???
+ float Mag_Derate_Percent; // Percentage reduction in power when mags are switched from 'both' to either 'L' or 'R'
+ float Gear_Ratio; // Gearing between engine and propellor
+ float n_R; // Number of cycles per power stroke - 2 for a 4 stroke engine.
+
+ // Engine Variables not read in from config file
+ float RPM; // rpm
float Percentage_Power; // Power output as percentage of maximum power output
- float Manifold_Pressure; // Inches
- float RPM;
- float Fuel_Flow_gals_hr; // gals/hour
- float Torque;
- float CHT; // Cylinder head temperature deg K
+ float Manifold_Pressure; // Inches Hg
+ float Fuel_Flow_gals_hr; // USgals/hour
+ float Torque_lbft; // lb-ft
+ float Torque_SI; // Nm
+ float CHT_degK; // Cylinder head temperature deg K
float CHT_degF; // Ditto in deg Fahrenheit
- float EGT; // Exhaust gas temperature deg K
- float EGT_degF; // Exhaust gas temperature deg Fahrenheit
float Mixture;
float Oil_Pressure; // PSI
float Oil_Temp; // Deg C
+ float current_oil_temp; // deg K
+ /**** one of these is superfluous !!!!***/
float HP; // Current power output in HP
float Power_SI; // Current power output in Watts
- float Torque_SI; // Torque in Nm
- float Torque_FMEP; // The component of Engine torque due to FMEP (Nm)
float RPS;
- float Torque_Imbalance;
- bool running; //flag to indicate the engine is running self sustaining
- bool cranking; //flag to indicate the engine is being cranked
- bool spark; //flag to indicate a spark is available
- bool fuel; //flag to indicate fuel is available
-
- //DCL
+ float angular_velocity_SI; // rad/s
+ float Torque_FMEP; // The component of Engine torque due to FMEP (Nm)
+ float Torque_Imbalance; // difference between engine and prop torque
+ float EGT; // Exhaust gas temperature deg K
+ float EGT_degF; // Exhaust gas temperature deg Fahrenheit
float volumetric_efficiency;
float combustion_efficiency;
- float equivalence_ratio;
- float v_dot_air;
- float m_dot_air;
- float m_dot_fuel;
- float swept_volume;
- float True_Manifold_Pressure; //in Hg
- float rho_air_manifold;
- float R_air;
- float p_amb_sea_level; // Pascals
- float p_amb; // Pascals
- float T_amb; // deg Kelvin
- float calorific_value_fuel;
- float rho_air;
- float rho_fuel; // kg/m^3
- float thi_sea_level;
- float delta_T_exhaust;
- float displacement; // Engine displacement in cubic inches - to be read in from config file for each engine
- float displacement_SI; // ditto in meters cubed
- float Cp_air; // J/KgK
- float Cp_fuel; // J/KgK
- float heat_capacity_exhaust;
- float enthalpy_exhaust;
- float Percentage_of_best_power_mixture_power;
+ float equivalence_ratio; // ratio of stoichiometric AFR over actual AFR
+ float thi_sea_level; // the equivalence ratio we would be running at assuming sea level air denisity in the manifold
+ float v_dot_air; // volume flow rate of air into engine - m^3/s
+ float m_dot_air; // mass flow rate of air into engine - kg/s
+ float m_dot_fuel; // mass flow rate of fuel into engine - kg/s
+ float swept_volume; // total engine swept volume - m^3
+ /********* swept volume or the geometry used to calculate it should be in the config read section surely ??? ******/
+ float True_Manifold_Pressure; //in Hg - actual manifold gauge pressure
+ float rho_air_manifold; // denisty of air in the manifold - kg/m^3
+ float R_air; // Gas constant of air (287) UNITS???
+ float delta_T_exhaust; // Temperature change of exhaust this time step - degK
+ float heat_capacity_exhaust; // exhaust gas specific heat capacity - J/kgK
+ float enthalpy_exhaust; // Enthalpy at current exhaust gas conditions - UNITS???
+ float Percentage_of_best_power_mixture_power; // Current power as a percentage of what power we would have at the same conditions but at best power mixture
float abstract_mixture; //temporary hack
- float engine_inertia; //kg.m^2
- float prop_inertia; //kg.m^2
float angular_acceleration; //rad/s^2
- float n_R; //Number of cycles per power stroke
float FMEP; //Friction Mean Effective Pressure (Pa)
- double time_step;
+
+ // Various bits of housekeeping describing the engines state.
+ bool running; // flag to indicate the engine is running self sustaining
+ bool cranking; // flag to indicate the engine is being cranked
+ int crank_counter; // Number of iterations that the engine has been cranked non-stop
+ bool spark; // flag to indicate a spark is available
+ bool fuel; // flag to indicate fuel is available
// Propellor Variables
- float FGProp1_Thrust;
- float FGProp1_RPS;
- float FGProp1_Blade_Angle;
+ float FGProp1_RPS; // rps
float prop_torque; // Nm
float prop_thrust; // Newtons
- float blade_length; // meters
- float forward_velocity; // m/s
- float angular_velocity_SI; // rad/s
- float prop_power_consumed_SI; // Watts
- float prop_power_consumed_HP; // HP
- double prop_diameter; // meters
- double J; // advance ratio - dimensionless
- double Cp_20; // coefficient of power for 20 degree blade angle
- double Cp_25; // coefficient of power for 25 degree blade angle
- double Cp; // Our actual coefficient of power
- double blade_angle; // degrees
- double neta_prop_20;
- double neta_prop_25;
- double neta_prop; // prop efficiency
+ double prop_diameter; // meters
+ double blade_angle; // degrees
+
+// MEMBER FUNCTIONS
+
// Calculate Engine RPM based on Propellor Lever Position
float Calc_Engine_RPM(float Position);
// Calculate percentage of best power mixture power based on equivalence ratio
float Power_Mixture_Correlation(float thi_actual);
- // Calculate exhaust gas temperature rise
+ // Calculate exhaust gas temperature change
float Calculate_Delta_T_Exhaust(void);
+ // Calculate cylinder head temperature
+ float FGNewEngine::Calc_CHT(float CHT);
+
+ void FGNewEngine::Calc_EGT(void);
+
+ // Calculate fuel flow in gals/hr
+ void FGNewEngine::Calc_Fuel_Flow_Gals_Hr(void);
+
+ // Calculate current percentage power
+ void FGNewEngine::Calc_Percentage_Power(bool mag_left, bool mag_right);
+
// Calculate Oil Temperature
- float Calc_Oil_Temp (float Fuel_Flow, float Mixture, float IAS);
+ float Calc_Oil_Temp (float oil_temp);
// Calculate Oil Pressure
float Calc_Oil_Press (float Oil_Temp, float Engine_RPM);
+ // Propeller calculations.
+ void FGNewEngine::Do_Prop_Calcs(void);
+
public:
- ofstream outfile;
+// ofstream outfile;
//constructor
FGNewEngine() {
// accessors
inline float get_RPM() const { return RPM; }
inline float get_Manifold_Pressure() const { return True_Manifold_Pressure; }
- inline float get_FGProp1_Thrust() const { return FGProp1_Thrust; }
- inline float get_FGProp1_Blade_Angle() const { return FGProp1_Blade_Angle; }
-
// inline float get_Rho() const { return Rho; }
inline float get_MaxHP() const { return MaxHP; }
inline float get_Percentage_Power() const { return Percentage_Power; }
inline float get_prop_thrust_SI() const { return prop_thrust; }
inline float get_prop_thrust_lbs() const { return (prop_thrust * 0.2248); }
inline float get_fuel_flow_gals_hr() const { return (Fuel_Flow_gals_hr); }
+ inline float get_oil_temp() const { return ((current_oil_temp * 1.8) - 459.67); }
};
projects / flightgear.git / commitdiff