2 // Author: Phil Schubert
3 // Date started: 12/03/99
4 // Purpose: Models a Continental IO-520-M Engine
5 // Called by: FGSimExec
7 // Copyright (C) 1999 Philip L. Schubert (philings@ozemail.com.au)
9 // This program is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU General Public License as
11 // published by the Free Software Foundation; either version 2 of the
12 // License, or (at your option) any later version.
14 // This program is distributed in the hope that it will be useful, but
15 // WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 // General Public License for more details.
19 // You should have received a copy of the GNU General Public License
20 // along with this program; if not, write to the Free Software
21 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24 // Further information about the GNU General Public License can also
25 // be found on the world wide web at http://www.gnu.org.
27 // FUNCTIONAL DESCRIPTION
28 // ------------------------------------------------------------------------
29 // Models a Continental IO-520-M engine. This engine is used in Cessna
30 // 210, 310, Beechcraft Bonaza and Baron C55. The equations used below
31 // were determined by a first and second order curve fits using Excel.
32 // The data is from the Cessna Aircraft Corporations Engine and Flight
33 // Computer for C310. Part Number D3500-13
36 // ------------------------------------------------------------------------
40 // ------------------------------------------------------------------------
41 // 12/03/99 PLS Created
42 // 07/03/99 PLS Added Calculation of Density, and Prop_Torque
43 // 07/03/99 PLS Restructered Variables to allow easier implementation
45 // 15/03/99 PLS Added Oil Pressure, Oil Temperature and CH Temp
46 // ------------------------------------------------------------------------
48 // ------------------------------------------------------------------------
51 /////////////////////////////////////////////////////////////////////
53 // Modified by Dave Luff (david.luff@nottingham.ac.uk) September 2000
55 // Altered manifold pressure range to add a minimum value at idle to simulate the throttle stop / idle bypass valve,
56 // and to reduce the maximum value whilst the engine is running to slightly below ambient to account for CdA losses across the throttle
58 // Altered it a bit to model an IO360 from C172 - 360 cubic inches, 180 HP max, fixed pitch prop
59 // Added a simple fixed pitch prop model by Nev Harbor - this is not intended as a final model but simply a hack to get it running for now
60 // I used Phil's ManXRPM correlation for power rather than do a new one for the C172 for now, but altered it a bit to reduce power at the low end
62 // Added EGT model based on combustion efficiency and an energy balance with the exhaust gases
64 // Added a mixture - power correlation based on a curve in the IO360 operating manual
66 // I've tried to match the prop and engine model to give roughly 600 RPM idle and 180 HP at 2700 RPM
67 // but it is by no means currently at a completed stage - DCL 15/9/00
69 // DCL 28/9/00 - Added estimate of engine and prop inertia and changed engine speed calculation to be calculated from Angular acceleration = Torque / Inertia.
70 // Requires a timestep to be passed to FGNewEngine::init and currently assumes this timestep does not change.
71 // Could easily be altered to pass a variable timestep to FGNewEngine::update every step instead if required.
73 // DCL 27/10/00 - Added first stab at cylinder head temperature model
74 // See the comment block in the code for details
76 // DCL 02/11/00 - Modified EGT code to reduce values to those more representative of a sensor downstream
78 // DCL 02/02/01 - Changed the prop model to one based on efficiency and co-efficient of power curves from McCormick instead of the
79 // blade element method we were using previously. This works much better, and is similar to how Jon is doing it in JSBSim.
81 //////////////////////////////////////////////////////////////////////
83 #include <simgear/compiler.h>
93 float FGNewEngine::Lookup_Combustion_Efficiency(float thi_actual)
95 const int NUM_ELEMENTS = 11;
96 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
97 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
98 //combustion efficiency values from Heywood: ISBN 0-07-100499-8
99 float neta_comb_actual;
103 int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays
109 //this is just to avoid crashing the routine is we are bigger than the last element - for now just return the last element
110 //but at some point we will have to extrapolate further
111 neta_comb_actual = neta_comb[i];
112 return neta_comb_actual;
114 if(thi_actual == thi[i])
116 neta_comb_actual = neta_comb[i];
117 return neta_comb_actual;
119 if((thi_actual > thi[i]) && (thi_actual < thi[i + 1]))
121 //do linear interpolation between the two points
122 factor = (thi_actual - thi[i]) / (thi[i+1] - thi[i]);
123 neta_comb_actual = (factor * (neta_comb[i+1] - neta_comb[i])) + neta_comb[i];
124 return neta_comb_actual;
128 //if we get here something has gone badly wrong
129 cout << "ERROR: error in FGNewEngine::Lookup_Combustion_Efficiency\n";
130 return neta_comb_actual;
134 // Calculate Manifold Pressure based on Throttle lever Position
135 static float Calc_Manifold_Pressure ( float LeverPosn, float MaxMan, float MinMan)
142 //Note that setting the manifold pressure as a function of lever position only is not strictly accurate
143 //MAP is also a function of engine speed. (and ambient pressure if we are going for an actual MAP model)
144 Inches = MinMan + (LeverPosn * (MaxMan - MinMan) / 100);
146 //allow for idle bypass valve or slightly open throttle stop
156 // Calculate Oil Temperature
157 static float Oil_Temp (float Fuel_Flow, float Mixture, float IAS)
164 // Calculate Oil Pressure
165 static float Oil_Press (float Oil_Temp, float Engine_RPM)
167 float Oil_Pressure = 0; //PSI
168 float Oil_Press_Relief_Valve = 60; //PSI
169 float Oil_Press_RPM_Max = 1800;
170 float Design_Oil_Temp = 85; //Celsius
171 float Oil_Viscosity_Index = 0.25; // PSI/Deg C
172 float Temp_Deviation = 0; // Deg C
174 Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
176 // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
177 if (Oil_Pressure >= Oil_Press_Relief_Valve)
179 Oil_Pressure = Oil_Press_Relief_Valve;
182 // Now adjust pressure according to Temp which affects the viscosity
184 Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
190 // Calculate Density Ratio
191 static float Density_Ratio ( float x )
194 y = ((3E-10 * x * x) - (3E-05 * x) + 0.9998);
199 // Calculate Air Density - Rho
200 static float Density ( float x )
203 y = ((9E-08 * x * x) - (7E-08 * x) + 0.0024);
208 // Calculate Speed in FPS given Knots CAS
209 static float IAS_to_FPS (float x)
217 //*************************************************************************************
218 // Initialise the engine model
219 void FGNewEngine::init(double dt) {
221 // These constants should probably be moved eventually
222 CONVERT_CUBIC_INCHES_TO_METERS_CUBED = 1.638706e-5;
223 CONVERT_HP_TO_WATTS = 745.6999;
225 // Properties of working fluids
226 Cp_air = 1005; // J/KgK
227 Cp_fuel = 1700; // J/KgK
228 calorific_value_fuel = 47.3e6; // W/Kg Note that this is only an approximate value
229 rho_fuel = 800; // kg/m^3 - an estimate for now
232 // Control and environment inputs
234 Throttle_Lever_Pos = 75;
235 Propeller_Lever_Pos = 75;
236 Mixture_Lever_Pos = 100;
240 // Engine Specific Variables.
241 // Will be set in a parameter file to be read in to create
242 // and instance for each engine.
243 Max_Manifold_Pressure = 28.50; //Inches Hg. An approximation - should be able to find it in the engine performance data
244 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
246 Min_RPM = 600; //Recommended idle from Continental data sheet
248 Mag_Derate_Percent = 5;
249 // MaxHP = 285; //Continental IO520-M
250 MaxHP = 180; //Lycoming IO360
251 // displacement = 520; //Continental IO520-M
252 displacement = 360; //Lycoming IO360
253 displacement_SI = displacement * CONVERT_CUBIC_INCHES_TO_METERS_CUBED;
254 engine_inertia = 0.2; //kgm^2 - value taken from a popular family saloon car engine - need to find an aeroengine value !!!!!
255 prop_inertia = 0.03; //kgm^2 - this value is a total guess - dcl
262 // Initialise Engine Variables used by this instance
263 Percentage_Power = 0;
264 Manifold_Pressure = 29.00; // Inches
266 Fuel_Flow_gals_hr = 0;
269 CHT = 298.0; //deg Kelvin
270 CHT_degF = (CHT * 1.8) - 459.67; //deg Fahrenheit
272 Oil_Pressure = 0; // PSI
273 Oil_Temp = 85; // Deg C
276 Torque_Imbalance = 0;
278 // Initialise Propellor Variables used by this instance
279 FGProp1_Angular_V = 0;
280 FGProp1_Coef_Drag = 0.6;
284 FGProp1_Coef_Lift = 0.1;
286 FGProp1_Blade_Angle = 13.5;
287 FGProp_Fine_Pitch_Stop = 13.5;
289 // Other internal values
294 //*****************************************************************************
295 //*****************************************************************************
296 // update the engine model based on current control positions
297 void FGNewEngine::update() {
298 // Declare local variables
300 // const int num2 = 500; // default is 100, number if iterations to run
301 // const int num2 = 5; // default is 100, number if iterations to run
307 // Set up the new variables
308 float Blade_Station = 30;
309 float FGProp_Area = 1.405/3;
310 float PI = 3.1428571;
314 // 0 = Closed, 100 = Fully Open
315 // float Throttle_Lever_Pos = 75;
316 // 0 = Full Course 100 = Full Fine
317 // float Propeller_Lever_Pos = 75;
318 // 0 = Idle Cut Off 100 = Full Rich
319 // float Mixture_Lever_Pos = 100;
321 // Environmental Variables
323 // Temp Variation from ISA (Deg F)
324 float FG_ISA_VAR = 0;
325 // Pressure Altitude 1000's of Feet
326 float FG_Pressure_Ht = 0;
328 // Parameters that alter the operation of the engine.
329 int Fuel_Available = 1; // Yes = 1. Is there Fuel Available. Calculated elsewhere
330 int Alternate_Air_Pos =0; // Off = 0. Reduces power by 3 % for same throttle setting
331 int Magneto_Left = 1; // 1 = On. Reduces power by 5 % for same power lever settings
332 int Magneto_Right = 1; // 1 = On. Ditto, Both of the above though do not alter fuel flow
335 //==================================================================
336 // Engine & Environmental Inputs from elsewhere
338 // Calculate Air Density (Rho) - In FG this is calculated in
341 Rho = Density(FG_Pressure_Ht); // In FG FG_Pressure_Ht is "h"
342 // cout << "Rho = " << Rho << endl;
344 // Calculate Manifold Pressure (Engine 1) as set by throttle opening
347 Calc_Manifold_Pressure( Throttle_Lever_Pos, Max_Manifold_Pressure, Min_Manifold_Pressure );
348 // cout << "manifold pressure = " << Manifold_Pressure << endl;
350 //**************************FIXME*******************************************
351 //DCL - hack for testing - fly at sea level
354 p_amb_sea_level = 101325;
356 //DCL - next calculate m_dot_air and m_dot_fuel into engine
358 //calculate actual ambient pressure and temperature from altitude
359 //Then find the actual manifold pressure (the calculated one is the sea level pressure)
360 True_Manifold_Pressure = Manifold_Pressure * p_amb / p_amb_sea_level;
362 // RPM = Calc_Engine_RPM(Propeller_Lever_Pos);
364 // cout << "Initial engine RPM = " << RPM << endl;
366 // Desired_RPM = RPM;
370 //DCL - calculate mass air flow into engine based on speed and load - separate this out into a function eventually
371 //t_amb is actual temperature calculated from altitude
372 //calculate density from ideal gas equation
373 rho_air = p_amb / ( R_air * T_amb );
374 rho_air_manifold = rho_air * Manifold_Pressure / 29.6;
375 //calculate ideal engine volume inducted per second
376 swept_volume = (displacement_SI * (RPM / 60)) / 2; //This equation is only valid for a four stroke engine
377 //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
378 //Note that this is cylinder vol eff - the throttle drop is already accounted for in the MAP calculation
379 volumetric_efficiency = 0.8;
380 //Now use volumetric efficiency to calculate actual air volume inducted per second
381 v_dot_air = swept_volume * volumetric_efficiency;
382 //Now calculate mass flow rate of air into engine
383 m_dot_air = v_dot_air * rho_air_manifold;
385 //cout << "air = " << m_dot_air;
386 // cout << "rho air manifold " << rho_air_manifold << '\n';
387 // cout << "Swept volume " << swept_volume << '\n';
391 //DCL - now calculate fuel flow into engine based on air flow and mixture lever position
392 //assume lever runs from no flow at fully out to thi = 1.6 at fully in at sea level
393 //***TODO*** - MUST try and get some real idea of the actual full rich sea level mixture - this is important !!!
394 //also assume that the injector linkage is ideal - hence the set mixture is maintained at a given altitude throughout the speed and load range
395 thi_sea_level = 1.6 * ( Mixture_Lever_Pos / 100.0 );
396 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
397 m_dot_fuel = m_dot_air / 14.7 * equivalence_ratio;
398 Fuel_Flow_gals_hr = (m_dot_fuel / rho_fuel) * 264.172 * 3600.0; // Note this assumes US gallons
400 //cout << "fuel " << m_dot_fuel; << "kg/s " << Fuel_Flow_gals_hr << "gals/hr"
401 //cout << " air " << m_dot_air << '\n';
403 //***********************************************************************
404 //Engine power and torque calculations
406 // For a given Manifold Pressure and RPM calculate the % Power
407 // Multiply Manifold Pressure by RPM
408 ManXRPM = Manifold_Pressure * RPM;
413 // Phil's %power correlation
415 Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM) + ( + 7E-04 * ManXRPM) - 0.1218;
416 // cout << Percentage_Power << "%" << "\t";
419 // DCL %power correlation - basically Phil's correlation modified to give slighty less power at the low end
420 // might need some adjustment as the prop model is adjusted
421 // 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
422 // Calculate % Power for Nev's prop model
423 //Percentage_Power = (+ 6E-09 * ManXRPM * ManXRPM) + ( + 8E-04 * ManXRPM) - 1.8524;
425 // Calculate %power for DCL prop model
426 Percentage_Power = (7e-9 * ManXRPM * ManXRPM) + (7e-4 * ManXRPM) - 1.0;
428 // cout << Percentage_Power << "%" << "\t";
431 // Adjust for Temperature - Temperature above Standard decrease
432 // power % by 7/120 per degree F increase, and incease power for
433 // temps below at the same ratio
434 Percentage_Power = Percentage_Power - (FG_ISA_VAR * 7 /120);
435 // cout << Percentage_Power << "%" << "\t";
437 //******DCL - this bit will need altering or removing if I go to true altitude adjusted manifold pressure
438 // Adjust for Altitude. In this version a linear variation is
439 // used. Decrease 1% for each 1000' increase in Altitde
440 Percentage_Power = Percentage_Power + (FG_Pressure_Ht * 12/10000);
441 // cout << Percentage_Power << "%" << "\t";
444 //DCL - now adjust power to compensate for mixture
445 //uses a curve fit to the data in the IO360 / O360 operating manual
446 //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,
447 //possibly by using separate fits for rich and lean of best power mixture
448 //first adjust actual mixture to abstract mixture - this is a temporary hack in order to account for the fact that the data I have
449 //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.
451 abstract_mixture = 10.0 * equivalence_ratio - 12.0;
452 float m = abstract_mixture; //to simplify writing the next equation
453 Percentage_of_best_power_mixture_power = ((-0.0012*pow(m,6)) + (0.021*pow(m,5)) + (-0.1425*pow(m,4)) + (0.4395*pow(m,3)) + (-0.8909*m*m) + (-0.5155*m) + 100.03);
454 Percentage_Power = Percentage_Power * Percentage_of_best_power_mixture_power / 100.0;
456 //cout << " %POWER = " << Percentage_Power << '\n';
458 //***DCL - FIXME - this needs altering - for instance going richer than full power mixture decreases the %power but increases the fuel flow
459 // Now Calculate Fuel Flow based on % Power Best Power Mixture
460 // Fuel_Flow = Percentage_Power * Max_Fuel_Flow / 100.0;
461 // cout << Fuel_Flow << " lbs/hr"<< endl;
463 // Now Derate engine for the effects of Bad/Switched off magnetos
464 if (Magneto_Left == 0 && Magneto_Right == 0) {
465 // cout << "Both OFF\n";
466 Percentage_Power = 0;
467 } else if (Magneto_Left && Magneto_Right) {
468 // cout << "Both On ";
469 } else if (Magneto_Left == 0 || Magneto_Right== 0) {
470 // cout << "1 Magneto Failed ";
472 Percentage_Power = Percentage_Power *
473 ((100.0 - Mag_Derate_Percent)/100.0);
474 // cout << FGEng1_Percentage_Power << "%" << "\t";
477 HP = Percentage_Power * MaxHP / 100.0;
479 Power_SI = HP * CONVERT_HP_TO_WATTS;
481 // Calculate Engine Torque. Check for div by zero since percentage power correlation does not guarantee zero power at zero rpm.
486 Torque_SI = (Power_SI * 60.0) / (2.0 * PI * RPM); //Torque = power / angular velocity
487 // cout << Torque << " Nm\n";
490 //**********************************************************************
491 //Calculate Exhaust gas temperature
493 // cout << "Thi = " << equivalence_ratio << '\n';
495 combustion_efficiency = Lookup_Combustion_Efficiency(equivalence_ratio); //The combustion efficiency basically tells us what proportion of the fuels calorific value is released
497 // cout << "Combustion efficiency = " << combustion_efficiency << '\n';
499 //now calculate energy release to exhaust
500 //We will assume a three way split of fuel energy between useful work, the coolant system and the exhaust system
501 //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
502 //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.
503 enthalpy_exhaust = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
504 heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
505 delta_T_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
506 // delta_T_exhaust = Calculate_Delta_T_Exhaust();
508 // cout << "T_amb " << T_amb;
509 // cout << " dT exhaust = " << delta_T_exhaust;
511 EGT = T_amb + delta_T_exhaust;
513 //The above gives the exhaust temperature immediately prior to leaving the combustion chamber
514 //Now derate to give a more realistic figure as measured downstream
515 //For now we will aim for a peak of around 400 degC (750 degF)
517 EGT *= 0.444 + ((0.544 - 0.444) * Percentage_Power / 100.0);
519 EGT_degF = (EGT * 1.8) - 459.67;
521 //cout << " EGT = " << EGT << " degK " << EGT_degF << " degF";// << '\n';
523 //***************************************************************************************
524 // Calculate Cylinder Head Temperature
528 This is a somewhat rough first attempt at modelling cylinder head temperature. The cylinder head
529 is assumed to be at uniform temperature. Obviously this is incorrect, but it simplifies things a
530 lot, and we're just looking for the behaviour of CHT to be correct. Energy transfer to the cylinder
531 head is assumed to be one third of the energy released by combustion at all conditions. This is a
532 reasonable estimate, although obviously in real life it varies with different conditions and possibly
533 with CHT itself. I've split energy transfer from the cylinder head into 2 terms - free convection -
534 ie convection to stationary air, and forced convection, ie convection into flowing air. The basic
535 free convection equation is: dqdt = -hAdT Since we don't know A and are going to set h quite arbitarily
536 anyway I've knocked A out and just wrapped it up in h - the only real significance is that the units
537 of h will be different but that dosn't really matter to us anyway. In addition, we have the problem
538 that the prop model I'm currently using dosn't model the backwash from the prop which will add to the
539 velocity of the cooling air when the prop is turning, so I've added an extra term to try and cope
542 In real life, forced convection equations are genarally empirically derived, and are quite complicated
543 and generally contain such things as the Reynolds and Nusselt numbers to various powers. The best
544 course of action would probably to find an empirical correlation from the literature for a similar
545 situation and try and get it to fit well. However, for now I am using my own made up very simple
546 correlation for the energy transfer from the cylinder head:
548 dqdt = -(h1.dT) -(h2.m_dot.dT) -(h3.rpm.dT)
550 where dT is the temperature different between the cylinder head and the surrounding air, m_dot is the
551 mass flow rate of cooling air through an arbitary volume, rpm is the engine speed in rpm (this is the
552 backwash term), and h1, h2, h3 are co-efficients which we can play with to attempt to get the CHT
553 behaviour to match real life.
555 In order to change the values of CHT that the engine settles down at at various conditions,
556 have a play with h1, h2 and h3. In order to change the rate of heating/cooling without affecting
557 equilibrium values alter the cylinder head mass, which is really quite arbitary. Bear in mind that
558 altering h1, h2 and h3 will also alter the rate of heating or cooling as well as equilibrium values,
559 but altering the cylinder head mass will only alter the rate. It would I suppose be better to read
560 the values from file to avoid the necessity for re-compilation every time I change them.
563 //CHT = Calc_CHT( Fuel_Flow, Mixture, IAS);
564 // cout << "Cylinder Head Temp (F) = " << CHT << endl;
565 float h1 = -95.0; //co-efficient for free convection
566 float h2 = -3.95; //co-efficient for forced convection
567 float h3 = -0.05; //co-efficient for forced convection due to prop backwash
568 float v_apparent; //air velocity over cylinder head in m/s
569 float v_dot_cooling_air;
570 float m_dot_cooling_air;
571 float temperature_difference;
572 float arbitary_area = 1.0;
573 float dqdt_from_combustion;
574 float dqdt_forced; //Rate of energy transfer to/from cylinder head due to forced convection (Joules) (sign convention: to cylinder head is +ve)
575 float dqdt_free; //Rate of energy transfer to/from cylinder head due to free convection (Joules) (sign convention: to cylinder head is +ve)
576 float dqdt_cylinder_head; //Overall energy change in cylinder head
577 float CpCylinderHead = 800.0; //FIXME - this is a guess - I need to look up the correct value
578 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
579 float HeatCapacityCylinderHead;
582 temperature_difference = CHT - T_amb;
584 v_apparent = IAS * 0.5144444; //convert from knots to m/s
585 v_dot_cooling_air = arbitary_area * v_apparent;
586 m_dot_cooling_air = v_dot_cooling_air * rho_air;
588 //Calculate rate of energy transfer to cylinder head from combustion
589 dqdt_from_combustion = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
591 //Calculate rate of energy transfer from cylinder head due to cooling NOTE is calculated as rate to but negative
592 dqdt_forced = (h2 * m_dot_cooling_air * temperature_difference) + (h3 * RPM * temperature_difference);
593 dqdt_free = h1 * temperature_difference;
595 //Calculate net rate of energy transfer to or from cylinder head
596 dqdt_cylinder_head = dqdt_from_combustion + dqdt_forced + dqdt_free;
598 HeatCapacityCylinderHead = CpCylinderHead * MassCylinderHead;
600 dCHTdt = dqdt_cylinder_head / HeatCapacityCylinderHead;
602 CHT += (dCHTdt * time_step);
604 CHT_degF = (CHT * 1.8) - 459.67;
606 //cout << " CHT = " << CHT_degF << " degF\n";
609 // End calculate Cylinder Head Temperature
612 //***************************************************************************************
613 // Oil pressure calculation
615 // Calculate Oil Pressure
616 Oil_Pressure = Oil_Press( Oil_Temp, RPM );
617 // cout << "Oil Pressure (PSI) = " << Oil_Pressure << endl;
619 //**************************************************************************************
620 // Now do the Propeller Calculations
622 #ifdef NEVS_PROP_MODEL
627 number_of_blades = 2.0;
629 allowance_for_spinner = blade_length / 12.0;
630 prop_fudge_factor = 1.453401525;
631 forward_velocity = IAS;
640 angular_velocity_SI = 2.0 * PI * RPM / 60.0;
642 allowance_for_spinner = blade_length / 12.0;
643 //Calculate thrust and torque by summing the contributions from each of the blade elements
644 //Assumes equal length elements with numbered 1 inboard -> num_elements outboard
648 // outfile << "Rho = " << Rho << '\n\n';
649 // outfile << "Drag = ";
650 for(i=1;i<=num_elements;i++)
653 distance = (blade_length * (element / num_elements)) + allowance_for_spinner;
654 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))))))
655 * (0.1 * (blade_length / element)) * number_of_blades;
657 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)
658 * (0.1 * (blade_length / element)) * number_of_blades;
659 element_torque = element_drag * distance;
660 prop_torque += element_torque;
661 prop_thrust += element_lift;
662 // outfile << "Drag = " << element_drag << " n = " << element << '\n';
667 // outfile << "Angular velocity = " << angular_velocity_SI << " rad/s\n";
669 // cout << "Thrust = " << prop_thrust << '\n';
670 prop_thrust *= prop_fudge_factor;
671 prop_torque *= prop_fudge_factor;
672 prop_power_consumed_SI = prop_torque * angular_velocity_SI;
673 prop_power_consumed_HP = prop_power_consumed_SI / 745.699;
676 #endif //NEVS_PROP_MODEL
678 #ifdef DCL_PROP_MODEL
680 double prop_diameter = 1.8; // meters
681 double J; // advance ratio - dimensionless
682 double Cp_20; // coefficient of power for 20 degree blade angle
683 double Cp_25; // coefficient of power for 25 degree blade angle
684 double Cp; // our actual coefficient of power
685 double blade_angle = 23.0; // degrees
688 double neta_prop; // prop efficiency
691 FGProp1_RPS = RPM * Gear_Ratio / 60.0; // Borrow this variable from Phils model for now !!
692 angular_velocity_SI = 2.0 * PI * RPM / 60.0;
693 forward_velocity = IAS * 0.514444444444; // Convert to m/s
695 //cout << "Gear_Ratio = " << Gear_Ratio << '\n';
696 //cout << "IAS = " << IAS << '\n';
697 //cout << "forward_velocity = " << forward_velocity << '\n';
698 //cout << "FGProp1_RPS = " << FGProp1_RPS << '\n';
699 //cout << "prop_diameter = " << prop_diameter << '\n';
703 J = forward_velocity / (FGProp1_RPS * prop_diameter);
704 //cout << "advance_ratio = " << J << '\n';
706 //Cp correlations based on data from McCormick
707 Cp_20 = 0.0342*J*J*J*J - 0.1102*J*J*J + 0.0365*J*J - 0.0133*J + 0.064;
708 Cp_25 = 0.0119*J*J*J*J - 0.0652*J*J*J + 0.018*J*J - 0.0077*J + 0.0921;
710 //cout << "Cp_20 = " << Cp_20 << '\n';
711 //cout << "Cp_25 = " << Cp_25 << '\n';
713 //Assume that the blade angle is between 20 and 25 deg - reasonable for fixed pitch prop but won't hold for variable one !!!
714 Cp = Cp_20 + ( (Cp_25 - Cp_20) * ((blade_angle - 20)/(25 - 20)) );
715 //cout << "Cp = " << Cp << '\n';
716 //cout << "RPM = " << RPM << '\n';
717 //cout << "angular_velocity_SI = " << angular_velocity_SI << '\n';
719 prop_power_consumed_SI = Cp * rho_air * pow(FGProp1_RPS,3.0) * pow(prop_diameter,5.0);
720 //cout << "prop HP consumed = " << prop_power_consumed_SI / 745.699 << '\n';
721 if(angular_velocity_SI == 0)
724 prop_torque = prop_power_consumed_SI / angular_velocity_SI;
726 // calculate neta_prop here
727 neta_prop_20 = 0.1328*J*J*J*J - 1.3073*J*J*J + 0.3525*J*J + 1.5591*J + 0.0007;
728 neta_prop_25 = -0.3121*J*J*J*J + 0.4234*J*J*J - 0.7686*J*J + 1.5237*J - 0.0004;
729 neta_prop = neta_prop_20 + ( (neta_prop_25 - neta_prop_20) * ((blade_angle - 20)/(25 - 20)) );
731 //FIXME - need to check for zero forward velocity to avoid divide by zero
732 if(forward_velocity < 0.0001)
735 prop_thrust = neta_prop * prop_power_consumed_SI / forward_velocity; //TODO - rename forward_velocity to IAS_SI
736 //cout << "prop_thrust = " << prop_thrust << '\n';
738 #endif //DCL_PROP_MODEL
740 //Calculate new RPM from torque balance and inertia.
741 Torque_Imbalance = Torque_SI - prop_torque; //This gives a +ve value when the engine torque exeeds the prop torque
743 angular_acceleration = Torque_Imbalance / (engine_inertia + prop_inertia);
744 angular_velocity_SI += (angular_acceleration * time_step);
745 RPM = (angular_velocity_SI * 60) / (2.0 * PI);
747 //DCL - stall the engine if RPM drops below 550 - this is possible if mixture lever is pulled right out