X-Git-Url: https://git.mxchange.org/?a=blobdiff_plain;f=src%2FFDM%2FIO360.cxx;h=c75b9b3f2fea84f37ca9d4ac62a46a1444b9474d;hb=ee98995d30e75cda88c9866f3cb6a761fda7d078;hp=22e7f97f9d3ccef7cb6559adeb325117968b2f29;hpb=6d4e03361a8b0f8f10909fd7542694a83f478256;p=flightgear.git diff --git a/src/FDM/IO360.cxx b/src/FDM/IO360.cxx index 22e7f97f9..c75b9b3f2 100644 --- a/src/FDM/IO360.cxx +++ b/src/FDM/IO360.cxx @@ -1,10 +1,9 @@ -// 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 @@ -18,775 +17,684 @@ // // 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/9/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 FGEngine::init and currently assumes this timestep does not change. -// Could easily be altered to pass a variable timestep to FGEngine::update every step instead if required. -// -////////////////////////////////////////////////////////////////////// - -#include -#include -#include - -#include "IO360.hxx" - - -// ------------------------------------------------------------------------ -// CODE -// ------------------------------------------------------------------------ - - -// Calculate Engine RPM based on Propellor Lever Position -float FGEngine::Calc_Engine_RPM (float LeverPosition) -{ - // Calculate RPM as set by Prop Lever Position. Assumes engine - // will run at 1000 RPM at full course - - float RPM; - RPM = LeverPosition * Max_RPM / 100.0; - // * ((FGEng_Max_RPM + FGEng_Min_RPM) / 100); - - if ( RPM >= Max_RPM ) { - RPM = Max_RPM; - } - - return RPM; -} - -float FGEngine::Lookup_Combustion_Efficiency(float thi_actual) -{ - float thi[11]; //array of equivalence ratio values - float neta_comb[11]; //corresponding array of combustion efficiency values - float neta_comb_actual; - float factor; - - //thi = (0.0,0.9,1.0,1.05,1.1,1.15,1.2,1.3,1.4,1.5,1.6); - thi[0] = 0.0; - thi[1] = 0.9; - thi[2] = 1.0; - thi[3] = 1.05; //There must be an easier way of doing this !!!!!!!! - thi[4] = 1.1; - thi[5] = 1.15; - thi[6] = 1.2; - thi[7] = 1.3; - thi[8] = 1.4; - thi[9] = 1.5; - thi[10] = 1.6; - //neta_comb = (0.98,0.98,0.97,0.95,0.9,0.85,0.79,0.7,0.63,0.57,0.525); - neta_comb[0] = 0.98; - neta_comb[1] = 0.98; - neta_comb[2] = 0.97; - neta_comb[3] = 0.95; - neta_comb[4] = 0.9; - neta_comb[5] = 0.85; - neta_comb[6] = 0.79; - neta_comb[7] = 0.7; - neta_comb[8] = 0.63; - neta_comb[9] = 0.57; - neta_comb[10] = 0.525; - //combustion efficiency values from Heywood: ISBN 0-07-100499-8 - - int i; - int j; - j = 11; //This must be equal to the number of elements in the lookup table arrays +// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. - for(i=0;i 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; - } - } +#include - //if we get here something has gone badly wrong - cout << "ERROR: error in FGEngine::Lookup_Combustion_Efficiency\n"; - //exit(-1); - return neta_comb_actual; //keeps the compiler happy -} -/* -float FGEngine::Calculate_Delta_T_Exhaust(void) -{ - float dT_exhaust; - heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel); - dT_exhaust = enthalpy_exhaust / heat_capacity_exhaust; - - return(dT_exhaust); -} -*/ - -// Calculate Manifold Pressure based on Throttle lever Position -static float Calc_Manifold_Pressure ( float LeverPosn, float MaxMan, float MinMan) -{ - float Inches; - // if ( x < = 0 ) { - // x = 0.00001; - // } +#include - //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. - Inches = MinMan + (LeverPosn * (MaxMan - MinMan) / 100); +#include STL_FSTREAM +#include STL_IOSTREAM - //allow for idle bypass valve or slightly open throttle stop - if(Inches < MinMan) - Inches = MinMan; +#if !defined(SG_HAVE_NATIVE_SGI_COMPILERS) +SG_USING_STD(cout); +#endif - return Inches; -} +#include "IO360.hxx" +#include "LaRCsim/ls_constants.h" +#include
-// set initial default values -void FGEngine::init(double dt) { +//************************************************************************************* +// Initialise the engine model +void FGNewEngine::init(double dt) { + // These constants should probably be moved eventually CONVERT_CUBIC_INCHES_TO_METERS_CUBED = 1.638706e-5; - // Control and environment inputs - IAS = 0; - Throttle_Lever_Pos = 75; - Propeller_Lever_Pos = 75; - Mixture_Lever_Pos = 100; + CONVERT_HP_TO_WATTS = 745.6999; + + // Properties of working fluids Cp_air = 1005; // J/KgK Cp_fuel = 1700; // J/KgK 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; + + // 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 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. + // 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 = 360; //Lycoming IO360 + 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 - Max_Fuel_Flow = 130; Mag_Derate_Percent = 5; -// MaxHP = 285; //Continental IO520-M - MaxHP = 180; //Lycoming IO360 -// displacement = 520; //Continental IO520-M - displacement = 360; //Lycoming IO360 - engine_inertia = 0.2; //kgm^2 - value taken from a popular family saloon car engine - need to find an aeroengine value !!!!! - prop_inertia = 0.03; //kgm^2 - this value is a total guess - dcl - displacement_SI = displacement * CONVERT_CUBIC_INCHES_TO_METERS_CUBED; - Gear_Ratio = 1; - started = true; - cranking = false; + n_R = 2; // Number of crank revolutions per power cycle - 2 for a 4 stroke engine. - CONVERT_HP_TO_WATTS = 745.6999; -// ofstream outfile; - // outfile.open(ios::out|ios::trunc); + // Various bits of housekeeping describing the engines initial state. + running = false; + cranking = false; + crank_counter = false; + starter = false; // Initialise Engine Variables used by this instance + if(running) + RPM = 600; + else + RPM = 0; Percentage_Power = 0; - Manifold_Pressure = 29.00; // Inches - RPM = 600; - Fuel_Flow = 0; // lbs/hour - Torque = 0; - CHT = 370; + Manifold_Pressure = 29.96; // Inches + Fuel_Flow_gals_hr = 0; +// Torque = 0; + Torque_SI = 0; + CHT = 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; - Desired_RPM = 2500; //Recommended cruise RPM from Continental datasheet // Initialise Propellor Variables used by this instance - FGProp1_Angular_V = 0; - FGProp1_Coef_Drag = 0.6; - FGProp1_Torque = 0; - FGProp1_Thrust = 0; FGProp1_RPS = 0; - FGProp1_Coef_Lift = 0.1; - Alpha1 = 13.5; - FGProp1_Blade_Angle = 13.5; - FGProp_Fine_Pitch_Stop = 13.5; - - // Other internal values - Rho = 0.002378; -} - - -// Calculate Oil Pressure -static float 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; + // Hardcode propellor for now - the following two should be read from config eventually + prop_diameter = 1.8; // meters + blade_angle = 23.0; // degrees } - -// Calculate Cylinder Head Temperature -static float Calc_CHT (float Fuel_Flow, float Mixture, float IAS) -{ - float CHT = 350; - - return CHT; -} +//***************************************************************************** +// update the engine model based on current control positions +void FGNewEngine::update() { /* -//Calculate Exhaust Gas Temperature -//For now we will simply adjust this as a function of mixture -//It may be necessary to consider fuel flow rates and CHT in the calculation in the future -static float Calc_EGT (float Mixture) -{ - float EGT = 1000; //off the top of my head !!!! - //Now adjust for mixture strength - - return EGT; -}*/ - - -// 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 -static float Density ( float x ) -{ - float y ; - y = ((9E-08 * x * x) - (7E-08 * x) + 0.0024); - return(y); -} - + // Hack for testing - should output every 5 seconds + static int count1 = 0; + if(count1 == 0) { +// cout << "P_atmos = " << p_amb << " T_atmos = " << T_amb << '\n'; +// cout << "Manifold pressure = " << Manifold_Pressure << " True_Manifold_Pressure = " << True_Manifold_Pressure << '\n'; +// cout << "p_amb_sea_level = " << p_amb_sea_level << '\n'; +// cout << "equivalence_ratio = " << equivalence_ratio << '\n'; +// cout << "combustion_efficiency = " << combustion_efficiency << '\n'; +// cout << "AFR = " << 14.7 / equivalence_ratio << '\n'; +// 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 << "Percentage_Power = " << Percentage_Power << '\n'; +// cout << "current_oil_temp = " << current_oil_temp << '\n'; + cout << "EGT = " << EGT << '\n'; + } + count1++; + if(count1 == 100) + count1 = 0; +*/ -// Calculate Speed in FPS given Knots CAS -static float IAS_to_FPS (float x) -{ - float y; - y = x * 1.68888888; - return y; -} + // Check parameters that may alter the operating state of the engine. + // (spark, fuel, starter motor etc) + + // Check for spark + bool Magneto_Left = false; + bool Magneto_Right = false; + // Magneto positions: + // 0 -> off + // 1 -> left only + // 2 -> right only + // 3 -> both + if(mag_pos != 0) { + spark = true; + } else { + spark = false; + } // neglects battery voltage, master on switch, etc for now. + if((mag_pos == 1) || (mag_pos > 2)) + Magneto_Left = true; + if(mag_pos > 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)) { + fuel = true; + } else { + fuel = false; + } // Need to make this better, eg position of fuel selector switch. + + // Check if we are turning the starter motor + if(cranking != starter) { + // This check saves .../cranking from getting updated every loop - they only update when changed. + cranking = starter; + 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 + // 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) + RPM = 480; + } else { + // consider making a horrible noise if the starter is engaged with the engine running + } + } + 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; +// RPM = 600; + } + } + if( (running) && ((!spark)||(!fuel)) ) { + // Cut the engine + // note that we only cut the power - the engine may continue to spin if the prop is in a moving airstream + running = false; + } + // Now we've ascertained whether the engine is running or not we can start to do the engine calculations 'proper' -// update the engine model based on current control positions -void FGEngine::update() { - // Declare local variables - int num = 0; - // const int num2 = 500; // default is 100, number if iterations to run - const int num2 = 5; // default is 100, number if iterations to run - float ManXRPM = 0; - float Vo = 0; - float V1 = 0; - - - // Set up the new variables - float Blade_Station = 30; - float FGProp_Area = 1.405/3; - float PI = 3.1428571; - - // Input Variables - - // 0 = Closed, 100 = Fully Open - // float Throttle_Lever_Pos = 75; - // 0 = Full Course 100 = Full Fine - // float Propeller_Lever_Pos = 75; - // 0 = Idle Cut Off 100 = Full Rich - // float Mixture_Lever_Pos = 100; - - // Environmental Variables - - // Temp Variation from ISA (Deg F) - float FG_ISA_VAR = 0; - // Pressure Altitude 1000's of Feet - float FG_Pressure_Ht = 0; - - // Parameters that alter the operation of the engine. - // Yes = 1. Is there Fuel Available. Calculated elsewhere - int Fuel_Available = 1; - // Off = 0. Reduces power by 3 % for same throttle setting - int Alternate_Air_Pos =0; - // 1 = On. Reduces power by 5 % for same power lever settings - int Magneto_Left = 1; - // 1 = On. Ditto, Both of the above though do not alter fuel flow - int Magneto_Right = 1; - - // There needs to be a section in here to trap silly values, like - // 0, otherwise they will crash the calculations - - // cout << " Number of Iterations "; - // cin >> num2; - // cout << endl; - - // cout << " Throttle % "; - // cin >> Throttle_Lever_Pos; - // cout << endl; - - // cout << " Prop % "; - // cin >> Propeller_Lever_Pos; - // cout << endl; - - //================================================================== - // Engine & Environmental Inputs from elsewhere - - // Calculate Air Density (Rho) - In FG this is calculated in - // FG_Atomoshere.cxx - - Rho = Density(FG_Pressure_Ht); // In FG FG_Pressure_Ht is "h" - // cout << "Rho = " << Rho << endl; - - // Calculate Manifold Pressure (Engine 1) as set by throttle opening - - Manifold_Pressure = - Calc_Manifold_Pressure( Throttle_Lever_Pos, Max_Manifold_Pressure, Min_Manifold_Pressure ); + // 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; - //DCL - hack for testing - fly at sea level - T_amb = 298.0; - p_amb = 101325; - p_amb_sea_level = 101325; - - //DCL - next calculate m_dot_air and m_dot_fuel into engine - - //calculate actual ambient pressure and temperature from altitude //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; - // RPM = Calc_Engine_RPM(Propeller_Lever_Pos); - // RPM = 600; - // cout << "Initial engine RPM = " << RPM << endl; - -// Desired_RPM = RPM; - -//************** - - //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; - //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 - 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; - - // cout << "rho air manifold " << rho_air_manifold << '\n'; - // cout << "Swept volume " << swept_volume << '\n'; - -//************** - - //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.6 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.6 * ( 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; - - // cout << "fuel " << m_dot_fuel; - // cout << " air " << m_dot_air << '\n'; - -//************** + //Do the fuel flow calculations + Calc_Fuel_Flow_Gals_Hr(); + + //Calculate engine power + Calc_Percentage_Power(Magneto_Left, Magneto_Right); + HP = Percentage_Power * MaxHP / 100.0; + Power_SI = HP * CONVERT_HP_TO_WATTS; + + // 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 ??? - // cout << "Thi = " << equivalence_ratio << '\n'; + Torque_FMEP = (FMEP * displacement_SI) / (2.0 * LS_PI * n_R); - combustion_efficiency = Lookup_Combustion_Efficiency(equivalence_ratio); //The combustion efficiency basically tells us what proportion of the fuels calorific value is released + // 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"; + } - // cout << "Combustion efficiency = " << combustion_efficiency << '\n'; + //Calculate Exhaust gas temperature + if(running) + Calc_EGT(); + else + EGT = 298.0; - //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 con 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(); + // Calculate Cylinder Head Temperature + Calc_CHT(); + + // 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); + + // 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; + } 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; + EGT = 298.0; + } + } - // cout << "T_amb " << T_amb; - // cout << " dT exhaust = " << delta_T_exhaust; + // And finally, do any unit conversions from internal units to output units + EGT_degF = (EGT * 1.8) - 459.67; + CHT_degF = (CHT * 1.8) - 459.67; +} - EGT = T_amb + delta_T_exhaust; +//***************************************************************************************************** - // cout << " EGT = " << EGT << '\n'; +// FGNewEngine member functions - - // Calculate Manifold Pressure (Engine 2) as set by throttle opening +//////////////////////////////////////////////////////////////////////////////////////////// +// 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; - // FGEng2_Manifold_Pressure = Manifold_Pressure(FGEng2_Throttle_Lever_Pos, FGEng2_Manifold_Pressure); - // Show_Manifold_Pressure(FGEng2_Manifold_Pressure); + int i; + int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays + for(i=0;i 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; +} - //================================================================== - // Engine Power & Torque Calculations +//////////////////////////////////////////////////////////////////////////////////////////// +// 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; - // Loop until stable - required for testing only -// for (num = 0; num < num2; num++) { - // cout << Manifold_Pressure << " Inches" << "\t"; - // cout << RPM << " RPM" << "\t"; + int i; + int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays - // For a given Manifold Pressure and RPM calculate the % Power - // Multiply Manifold Pressure by RPM - ManXRPM = Manifold_Pressure * RPM; - // cout << ManXRPM; - // cout << endl; + for(i=0;i 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; + } + } -// Phil's %power correlation -/* // Calculate % Power - Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM) - + ( + 7E-04 * ManXRPM) - 0.1218; - // cout << Percentage_Power << "%" << "\t"; */ + //if we get here something has gone badly wrong + cout << "ERROR: error in FGNewEngine::Power_Mixture_Correlation\n"; + return mixPerPow_actual; +} -// 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 - Percentage_Power = (+ 6E-09 * ManXRPM * ManXRPM) - + ( + 8E-04 * ManXRPM) - 1.8524; - // cout << Percentage_Power << "%" << "\t"; - - // 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 - Percentage_Power = Percentage_Power - (FG_ISA_VAR * 7 /120); - // cout << Percentage_Power << "%" << "\t"; +// Calculate Cylinder Head Temperature +// 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!!! +void FGNewEngine::Calc_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 v_apparent; //air velocity over cylinder head in m/s + float v_dot_cooling_air; + float m_dot_cooling_air; + float temperature_difference; + float arbitary_area = 1.0; + float dqdt_from_combustion; + float dqdt_forced; //Rate of energy transfer to/from cylinder head due to forced convection (Joules) (sign convention: to cylinder head is +ve) + float dqdt_free; //Rate of energy transfer to/from cylinder head due to free convection (Joules) (sign convention: to cylinder head is +ve) + float dqdt_cylinder_head; //Overall energy change in cylinder head + float CpCylinderHead = 800.0; //FIXME - this is a guess - I need to look up the correct value + 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 + 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 + v_dot_cooling_air = arbitary_area * v_apparent; + m_dot_cooling_air = v_dot_cooling_air * rho_air; + + //Calculate rate of energy transfer to cylinder head from combustion + dqdt_from_combustion = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33; - // Adjust for Altitude. In this version a linear variation is - // used. Decrease 1% for each 1000' increase in Altitde - Percentage_Power = Percentage_Power + (FG_Pressure_Ht * 12/10000); - // cout << Percentage_Power << "%" << "\t"; - - //DCL - now adjust power to compensate for mixture - //uses a curve fit to the data in the IO360 / O360 operating manual - //due to the shape of the curve I had to use a 6th order fit - I am sure it must be possible to reduce this in future, - //possibly by using separate fits for rich and lean of best power mixture - //first adjust actual mixture to abstract mixture - this is a temporary hack in order to account for the fact that the data I have - //dosn't specify actual mixtures and I want to be able to change what I think they are without redoing the curve fit each time. - //y=10x-12 for now - abstract_mixture = 10.0 * equivalence_ratio - 12.0; - float m = abstract_mixture; //to simplify writing the next equation - Percentage_of_best_power_mixture_power = ((-0.0012*m*m*m*m*m*m) + (0.021*m*m*m*m*m) + (-0.1425*m*m*m*m) + (0.4395*m*m*m) + (-0.8909*m*m) + (-0.5155*m) + 100.03); - Percentage_Power = Percentage_Power * Percentage_of_best_power_mixture_power / 100.0; - - - // Now Calculate Fuel Flow based on % Power Best Power Mixture - Fuel_Flow = Percentage_Power * Max_Fuel_Flow / 100.0; - // cout << Fuel_Flow << " lbs/hr"<< endl; - - // Now Derate engine for the effects of Bad/Switched off magnetos - if (Magneto_Left == 0 && Magneto_Right == 0) { - // 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"; - } - - // Calculate Engine Horsepower - - HP = Percentage_Power * MaxHP / 100.0; - - Power_SI = HP * CONVERT_HP_TO_WATTS; - - // Calculate Engine Torque - - Torque = HP * 5252 / RPM; - // cout << Torque << "Ft/lbs" << "\t"; - - Torque_SI = (Power_SI * 60.0) / (2.0 * PI * RPM); //Torque = power / angular velocity - // cout << Torque << " Nm\n"; - - // Calculate Cylinder Head Temperature - CHT = Calc_CHT( Fuel_Flow, Mixture, IAS); - // cout << "Cylinder Head Temp (F) = " << CHT << endl; - -// EGT = Calc_EGT( Mixture ); - - // Calculate Oil Pressure - Oil_Pressure = Oil_Press( Oil_Temp, RPM ); - // cout << "Oil Pressure (PSI) = " << Oil_Pressure << endl; - - //============================================================== - - // Now do the Propellor Calculations - -#ifdef PHILS_PROP_MODEL - - // Revs per second - FGProp1_RPS = RPM * Gear_Ratio / 60.0; - // cout << FGProp1_RPS << " RPS" << endl; - - //Radial Flow Vector (V2) Ft/sec at Ref Blade Station (usually 30") - FGProp1_Angular_V = FGProp1_RPS * 2 * PI * (Blade_Station / 12); - // cout << FGProp1_Angular_V << "Angular Velocity " << endl; - - // Axial Flow Vector (Vo) Ft/sec - // Some further work required here to allow for inflow at low speeds - // Vo = (IAS + 20) * 1.688888; - Vo = IAS_to_FPS(IAS + 20); - // cout << "Feet/sec = " << Vo << endl; - - // cout << Vo << "Axial Velocity" << endl; - - // Relative Velocity (V1) - V1 = sqrt((FGProp1_Angular_V * FGProp1_Angular_V) + - (Vo * Vo)); - // cout << V1 << "Relative Velocity " << endl; - - // cout << FGProp1_Blade_Angle << " Prop Blade Angle" << endl; - - // Blade Angle of Attack (Alpha1) - -/* cout << " Alpha1 = " << Alpha1 - << " Blade angle = " << FGProp1_Blade_Angle - << " Vo = " << Vo - << " FGProp1_Angular_V = " << FGProp1_Angular_V << endl;*/ - Alpha1 = FGProp1_Blade_Angle -(atan(Vo / FGProp1_Angular_V) * (180/PI)); - // cout << Alpha1 << " Alpha1" << endl; - - // Calculate Coefficient of Drag at Alpha1 - FGProp1_Coef_Drag = (0.0005 * (Alpha1 * Alpha1)) + (0.0003 * Alpha1) - + 0.0094; - // cout << FGProp1_Coef_Drag << " Coef Drag" << endl; - - // Calculate Coefficient of Lift at Alpha1 - FGProp1_Coef_Lift = -(0.0026 * (Alpha1 * Alpha1)) + (0.1027 * Alpha1) - + 0.2295; - // cout << FGProp1_Coef_Lift << " Coef Lift " << endl; - - // Covert Alplha1 to Radians - // Alpha1 = Alpha1 * PI / 180; - - // Calculate Prop Torque - FGProp1_Torque = (0.5 * Rho * (V1 * V1) * FGProp_Area - * ((FGProp1_Coef_Lift * sin(Alpha1 * PI / 180)) - + (FGProp1_Coef_Drag * cos(Alpha1 * PI / 180)))) - * (Blade_Station/12); - // cout << FGProp1_Torque << " Prop Torque" << endl; - - // Calculate Prop Thrust - // cout << " V1 = " << V1 << " Alpha1 = " << Alpha1 << endl; - FGProp1_Thrust = 0.5 * Rho * (V1 * V1) * FGProp_Area - * ((FGProp1_Coef_Lift * cos(Alpha1 * PI / 180)) - - (FGProp1_Coef_Drag * sin(Alpha1 * PI / 180))); - // cout << FGProp1_Thrust << " Prop Thrust " << endl; - - // End of Propeller Calculations - //============================================================== - -#endif //PHILS_PROP_MODEL - -#ifdef NEVS_PROP_MODEL - - // Nev's prop model - - num_elements = 6.0; - number_of_blades = 2.0; - blade_length = 0.95; - allowance_for_spinner = blade_length / 12.0; - prop_fudge_factor = 1.453401525; - forward_velocity = IAS; - - theta[0] = 25.0; - theta[1] = 20.0; - theta[2] = 15.0; - theta[3] = 10.0; - theta[4] = 5.0; - theta[5] = 0.0; - - angular_velocity_SI = 2.0 * PI * RPM / 60.0; - - allowance_for_spinner = blade_length / 12.0; - //Calculate thrust and torque by summing the contributions from each of the blade elements - //Assumes equal length elements with numbered 1 inboard -> num_elements outboard - prop_torque = 0.0; - prop_thrust = 0.0; - int i; -// outfile << "Rho = " << Rho << '\n\n'; -// outfile << "Drag = "; - for(i=1;i<=num_elements;i++) - { - element = float(i); - distance = (blade_length * (element / num_elements)) + allowance_for_spinner; - element_drag = 0.5 * rho_air * ((distance * angular_velocity_SI)*(distance * angular_velocity_SI)) * (0.000833 * ((theta[int(element-1)] - (atan(forward_velocity/(distance * angular_velocity_SI))))*(theta[int(element-1)] - (atan(forward_velocity/(distance * angular_velocity_SI)))))) - * (0.1 * (blade_length / element)) * number_of_blades; - - element_lift = 0.5 * rho_air * ((distance * angular_velocity_SI)*(distance * angular_velocity_SI)) * (0.036 * (theta[int(element-1)] - (atan(forward_velocity/(distance * angular_velocity_SI))))+0.4) - * (0.1 * (blade_length / element)) * number_of_blades; - element_torque = element_drag * distance; - prop_torque += element_torque; - prop_thrust += element_lift; -// outfile << "Drag = " << element_drag << " n = " << element << '\n'; - } + //Calculate rate of energy transfer from cylinder head due to cooling NOTE is calculated as rate to but negative + dqdt_forced = (h2 * m_dot_cooling_air * temperature_difference) + (h3 * RPM * temperature_difference); + dqdt_free = h1 * temperature_difference; + + //Calculate net rate of energy transfer to or from cylinder head + dqdt_cylinder_head = dqdt_from_combustion + dqdt_forced + dqdt_free; + + HeatCapacityCylinderHead = CpCylinderHead * MassCylinderHead; + + dCHTdt = dqdt_cylinder_head / HeatCapacityCylinderHead; + + CHT += (dCHTdt * time_step); +} -// outfile << '\n'; +// 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 -// outfile << "Angular velocity = " << angular_velocity_SI << " rad/s\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; - // cout << "Thrust = " << prop_thrust << '\n'; - prop_thrust *= prop_fudge_factor; - prop_torque *= prop_fudge_factor; - prop_power_consumed_SI = prop_torque * angular_velocity_SI; - prop_power_consumed_HP = prop_power_consumed_SI / 745.699; + 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) -#endif //NEVS_PROP_MODEL + EGT *= 0.444 + ((0.544 - 0.444) * Percentage_Power / 100.0); +} +// Calculate Manifold Pressure based on Throttle lever Position +float FGNewEngine::Calc_Manifold_Pressure ( float LeverPosn, float MaxMan, float MinMan) +{ + float Inches; -//#if 0 -#ifdef PHILS_PROP_MODEL //Do Torque calculations in Ft/lbs - yuk :-((( - Torque_Imbalance = FGProp1_Torque - Torque; + //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); - if (Torque_Imbalance > 5) { - RPM -= 14.5; - // FGProp1_RPM -= 25; -// FGProp1_Blade_Angle -= 0.75; - } + //allow for idle bypass valve or slightly open throttle stop + if(Inches < MinMan) + Inches = MinMan; - if (Torque_Imbalance < -5) { - RPM += 14.5; - // FGProp1_RPM += 25; -// FGProp1_Blade_Angle += 0.75; - } -#endif + //Check for stopped engine (crudest way of compensating for engine speed) + if(RPM == 0) + Inches = 29.92; + return Inches; +} -#ifdef NEVS_PROP_MODEL //use proper units - Nm - Torque_Imbalance = Torque_SI - prop_torque; //This gives a +ve value when the engine torque exeeds the prop torque +// 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; - angular_acceleration = Torque_Imbalance / (engine_inertia + prop_inertia); - angular_velocity_SI += (angular_acceleration * time_step); - RPM = (angular_velocity_SI * 60) / (2.0 * PI); -#endif +//************** + //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; /* - if( RPM > (Desired_RPM + 2)) { - FGProp1_Blade_Angle += 0.75; //This value could be altered depending on how far from the desired RPM we are - } - - if( RPM < (Desired_RPM - 2)) { - FGProp1_Blade_Angle -= 0.75; - } - - if (FGProp1_Blade_Angle < FGProp_Fine_Pitch_Stop) { - FGProp1_Blade_Angle = FGProp_Fine_Pitch_Stop; - } +// Phil's %power correlation + // Calculate % Power + Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM) + ( + 7E-04 * ManXRPM) - 0.1218; + // cout << Percentage_Power << "%" << "\t"; +*/ - if (RPM >= 2700) { - RPM = 2700; - } +// 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; */ - //end constant speed prop -//#endif + if(Percentage_Power < 0) + Percentage_Power = 0; +} - //DCL - stall the engine if RPM drops below 550 - this is possible if mixture lever is pulled right out - if(RPM < 550) - RPM = 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 + } -// outfile << "RPM = " << RPM << " Blade angle = " << FGProp1_Blade_Angle << " Engine torque = " << Torque << " Prop torque = " << FGProp1_Torque << '\n'; -// outfile << "RPM = " << RPM << " Engine torque = " << Torque_SI << " Prop torque = " << prop_torque << '\n'; + float dOilTempdt = (target_oil_temp - oil_temp) / time_constant; - // cout << FGEng1_RPM << " Blade_Angle " << FGProp1_Blade_Angle << endl << endl; + 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 - // cout << "Final engine RPM = " << RPM << '\n'; -} + 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; -// Functions + return Oil_Pressure; +} -// Calculate Oil Temperature -static float Oil_Temp (float Fuel_Flow, float Mixture, float IAS) +// Propeller calculations. +void FGNewEngine::Do_Prop_Calcs() { - float Oil_Temp = 85; - - return (Oil_Temp); + float Gear_Ratio = 1.0; + float forward_velocity; // m/s + float prop_power_consumed_SI; // Watts + 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 + + if(FGProp1_RPS == 0) + J = 0; + else + J = forward_velocity / (FGProp1_RPS * prop_diameter); + //cout << "advance_ratio = " << J << '\n'; + + //Cp correlations based on data from McCormick + Cp_20 = 0.0342*J*J*J*J - 0.1102*J*J*J + 0.0365*J*J - 0.0133*J + 0.064; + Cp_25 = 0.0119*J*J*J*J - 0.0652*J*J*J + 0.018*J*J - 0.0077*J + 0.0921; + + //cout << "Cp_20 = " << Cp_20 << '\n'; + //cout << "Cp_25 = " << Cp_25 << '\n'; + + //Assume that the blade angle is between 20 and 25 deg - reasonable for fixed pitch prop but won't hold for variable one !!! + Cp = Cp_20 + ( (Cp_25 - Cp_20) * ((blade_angle - 20)/(25 - 20)) ); + //cout << "Cp = " << Cp << '\n'; + //cout << "RPM = " << RPM << '\n'; + //cout << "angular_velocity_SI = " << angular_velocity_SI << '\n'; + + prop_power_consumed_SI = Cp * rho_air * pow(FGProp1_RPS,3.0f) * pow(float(prop_diameter),5.0f); + //cout << "prop HP consumed = " << prop_power_consumed_SI / 745.699 << '\n'; + if(angular_velocity_SI == 0) + prop_torque = 0; + // 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 + // Unlikely to happen in practice - but I suppose someone could move the lever of a stopped large piston engine from feathered to windmilling. + // This *does* give problems - if the plane is put into a steep climb whilst windmilling the engine friction will eventually stop it spinning. + // 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!!! + else + prop_torque = prop_power_consumed_SI / angular_velocity_SI; + + // calculate neta_prop here + neta_prop_20 = 0.1328*J*J*J*J - 1.3073*J*J*J + 0.3525*J*J + 1.5591*J + 0.0007; + neta_prop_25 = -0.3121*J*J*J*J + 0.4234*J*J*J - 0.7686*J*J + 1.5237*J - 0.0004; + neta_prop = neta_prop_20 + ( (neta_prop_25 - neta_prop_20) * ((blade_angle - 20)/(25 - 20)) ); + + // Check for zero forward velocity to avoid divide by zero + if(forward_velocity < 0.0001) + prop_thrust = 0.0; + // I don't see how this works - how can the plane possibly start from rest ??? + // Hmmmm - it works because the forward_velocity at present never drops below about 0.03 even at rest + // We can't really rely on this in the future - needs fixing !!!! + // 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. + // It would be far more natural to work out the power consumed from the thrust - FIXME!!!!!. + else + prop_thrust = neta_prop * prop_power_consumed_SI / forward_velocity; //TODO - rename forward_velocity to IAS_SI + //cout << "prop_thrust = " << prop_thrust << '\n'; }