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 //////////////////////////////////////////////////////////////////////
75 #include <simgear/compiler.h>
85 // ------------------------------------------------------------------------
87 // ------------------------------------------------------------------------
90 // Calculate Engine RPM based on Propellor Lever Position
91 float FGNewEngine::Calc_Engine_RPM (float LeverPosition)
93 // Calculate RPM as set by Prop Lever Position. Assumes engine
94 // will run at 1000 RPM at full course
97 RPM = LeverPosition * Max_RPM / 100.0;
98 // * ((FGEng_Max_RPM + FGEng_Min_RPM) / 100);
100 if ( RPM >= Max_RPM ) {
108 float FGNewEngine::Lookup_Combustion_Efficiency(float thi_actual)
110 float thi[11]; //array of equivalence ratio values
111 float neta_comb[11]; //corresponding array of combustion efficiency values
112 float neta_comb_actual;
115 //thi = (0.0,0.9,1.0,1.05,1.1,1.15,1.2,1.3,1.4,1.5,1.6);
119 thi[3] = 1.05; //There must be an easier way of doing this !!!!!!!!
127 //neta_comb = (0.98,0.98,0.97,0.95,0.9,0.85,0.79,0.7,0.63,0.57,0.525);
138 neta_comb[10] = 0.525;
139 //combustion efficiency values from Heywood: ISBN 0-07-100499-8
143 j = 11; //This must be equal to the number of elements in the lookup table arrays
149 //this is just to avoid crashing the routine is we are bigger than the last element - for now just return the last element
150 //but at some point we will have to extrapolate further
151 neta_comb_actual = neta_comb[i];
152 return neta_comb_actual;
154 if(thi_actual == thi[i])
156 neta_comb_actual = neta_comb[i];
157 return neta_comb_actual;
159 if((thi_actual > thi[i]) && (thi_actual < thi[i + 1]))
161 //do linear interpolation between the two points
162 factor = (thi_actual - thi[i]) / (thi[i+1] - thi[i]);
163 neta_comb_actual = (factor * (neta_comb[i+1] - neta_comb[i])) + neta_comb[i];
164 return neta_comb_actual;
168 //if we get here something has gone badly wrong
169 cout << "ERROR: error in FGNewEngine::Lookup_Combustion_Efficiency\n";
171 return neta_comb_actual; //keeps the compiler happy
174 float FGNewEngine::Calculate_Delta_T_Exhaust(void)
177 heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
178 dT_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
184 // Calculate Manifold Pressure based on Throttle lever Position
185 static float Calc_Manifold_Pressure ( float LeverPosn, float MaxMan, float MinMan)
192 //Note that setting the manifold pressure as a function of lever position only is not strictly accurate
193 //MAP is also a function of engine speed.
194 Inches = MinMan + (LeverPosn * (MaxMan - MinMan) / 100);
196 //allow for idle bypass valve or slightly open throttle stop
204 // set initial default values
205 void FGNewEngine::init(double dt) {
207 CONVERT_CUBIC_INCHES_TO_METERS_CUBED = 1.638706e-5;
208 // Control and environment inputs
210 Throttle_Lever_Pos = 75;
211 Propeller_Lever_Pos = 75;
212 Mixture_Lever_Pos = 100;
213 Cp_air = 1005; // J/KgK
214 Cp_fuel = 1700; // J/KgK
215 calorific_value_fuel = 47.3e6; // W/Kg Note that this is only an approximate value
219 // Engine Specific Variables used by this program that have limits.
220 // Will be set in a parameter file to be read in to create
221 // and instance for each engine.
222 Max_Manifold_Pressure = 28.50; //Inches Hg. An approximation - should be able to find it in the engine performance data
223 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
225 Min_RPM = 600; //Recommended idle from Continental data sheet
227 Mag_Derate_Percent = 5;
228 // MaxHP = 285; //Continental IO520-M
229 MaxHP = 180; //Lycoming IO360
230 // displacement = 520; //Continental IO520-M
231 displacement = 360; //Lycoming IO360
232 engine_inertia = 0.2; //kgm^2 - value taken from a popular family saloon car engine - need to find an aeroengine value !!!!!
233 prop_inertia = 0.03; //kgm^2 - this value is a total guess - dcl
234 displacement_SI = displacement * CONVERT_CUBIC_INCHES_TO_METERS_CUBED;
240 CONVERT_HP_TO_WATTS = 745.6999;
242 // outfile.open(ios::out|ios::trunc);
244 // Initialise Engine Variables used by this instance
245 Percentage_Power = 0;
246 Manifold_Pressure = 29.00; // Inches
248 Fuel_Flow = 0; // lbs/hour
250 CHT = 298.0; //deg Kelvin
251 CHT_degF = (CHT * 1.8) - 459.67; //deg Fahrenheit
253 Oil_Pressure = 0; // PSI
254 Oil_Temp = 85; // Deg C
257 Torque_Imbalance = 0;
258 Desired_RPM = 2500; //Recommended cruise RPM from Continental datasheet
260 // Initialise Propellor Variables used by this instance
261 FGProp1_Angular_V = 0;
262 FGProp1_Coef_Drag = 0.6;
266 FGProp1_Coef_Lift = 0.1;
268 FGProp1_Blade_Angle = 13.5;
269 FGProp_Fine_Pitch_Stop = 13.5;
271 // Other internal values
276 // Calculate Oil Pressure
277 static float Oil_Press (float Oil_Temp, float Engine_RPM)
279 float Oil_Pressure = 0; //PSI
280 float Oil_Press_Relief_Valve = 60; //PSI
281 float Oil_Press_RPM_Max = 1800;
282 float Design_Oil_Temp = 85; //Celsius
283 float Oil_Viscosity_Index = 0.25; // PSI/Deg C
284 float Temp_Deviation = 0; // Deg C
286 Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
288 // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
289 if (Oil_Pressure >= Oil_Press_Relief_Valve)
291 Oil_Pressure = Oil_Press_Relief_Valve;
294 // Now adjust pressure according to Temp which affects the viscosity
296 Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
303 // Calculate Cylinder Head Temperature
304 static float Calc_CHT (float Fuel_Flow, float Mixture, float IAS, float rhoair, float tamb)
313 //Calculate Exhaust Gas Temperature
314 //For now we will simply adjust this as a function of mixture
315 //It may be necessary to consider fuel flow rates and CHT in the calculation in the future
316 static float Calc_EGT (float Mixture)
318 float EGT = 1000; //off the top of my head !!!!
319 //Now adjust for mixture strength
325 // Calculate Density Ratio
326 static float Density_Ratio ( float x )
329 y = ((3E-10 * x * x) - (3E-05 * x) + 0.9998);
334 // Calculate Air Density - Rho
335 static float Density ( float x )
338 y = ((9E-08 * x * x) - (7E-08 * x) + 0.0024);
343 // Calculate Speed in FPS given Knots CAS
344 static float IAS_to_FPS (float x)
352 //*****************************************************************************
353 //*****************************************************************************
354 // update the engine model based on current control positions
355 void FGNewEngine::update() {
356 // Declare local variables
358 // const int num2 = 500; // default is 100, number if iterations to run
359 // const int num2 = 5; // default is 100, number if iterations to run
365 // Set up the new variables
366 float Blade_Station = 30;
367 float FGProp_Area = 1.405/3;
368 float PI = 3.1428571;
372 // 0 = Closed, 100 = Fully Open
373 // float Throttle_Lever_Pos = 75;
374 // 0 = Full Course 100 = Full Fine
375 // float Propeller_Lever_Pos = 75;
376 // 0 = Idle Cut Off 100 = Full Rich
377 // float Mixture_Lever_Pos = 100;
379 // Environmental Variables
381 // Temp Variation from ISA (Deg F)
382 float FG_ISA_VAR = 0;
383 // Pressure Altitude 1000's of Feet
384 float FG_Pressure_Ht = 0;
386 // Parameters that alter the operation of the engine.
387 int Fuel_Available = 1; // Yes = 1. Is there Fuel Available. Calculated elsewhere
388 int Alternate_Air_Pos =0; // Off = 0. Reduces power by 3 % for same throttle setting
389 int Magneto_Left = 1; // 1 = On. Reduces power by 5 % for same power lever settings
390 int Magneto_Right = 1; // 1 = On. Ditto, Both of the above though do not alter fuel flow
393 //==================================================================
394 // Engine & Environmental Inputs from elsewhere
396 // Calculate Air Density (Rho) - In FG this is calculated in
399 Rho = Density(FG_Pressure_Ht); // In FG FG_Pressure_Ht is "h"
400 // cout << "Rho = " << Rho << endl;
402 // Calculate Manifold Pressure (Engine 1) as set by throttle opening
405 Calc_Manifold_Pressure( Throttle_Lever_Pos, Max_Manifold_Pressure, Min_Manifold_Pressure );
406 // cout << "manifold pressure = " << Manifold_Pressure << endl;
408 //**************************FIXME*******************************************
409 //DCL - hack for testing - fly at sea level
412 p_amb_sea_level = 101325;
414 //DCL - next calculate m_dot_air and m_dot_fuel into engine
416 //calculate actual ambient pressure and temperature from altitude
417 //Then find the actual manifold pressure (the calculated one is the sea level pressure)
418 True_Manifold_Pressure = Manifold_Pressure * p_amb / p_amb_sea_level;
420 // RPM = Calc_Engine_RPM(Propeller_Lever_Pos);
422 // cout << "Initial engine RPM = " << RPM << endl;
424 // Desired_RPM = RPM;
428 //DCL - calculate mass air flow into engine based on speed and load - separate this out into a function eventually
429 //t_amb is actual temperature calculated from altitude
430 //calculate density from ideal gas equation
431 rho_air = p_amb / ( R_air * T_amb );
432 rho_air_manifold = rho_air * Manifold_Pressure / 29.6;
433 //calculate ideal engine volume inducted per second
434 swept_volume = (displacement_SI * (RPM / 60)) / 2; //This equation is only valid for a four stroke engine
435 //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
436 volumetric_efficiency = 0.8;
437 //Now use volumetric efficiency to calculate actual air volume inducted per second
438 v_dot_air = swept_volume * volumetric_efficiency;
439 //Now calculate mass flow rate of air into engine
440 m_dot_air = v_dot_air * rho_air_manifold;
442 // cout << "rho air manifold " << rho_air_manifold << '\n';
443 // cout << "Swept volume " << swept_volume << '\n';
447 //DCL - now calculate fuel flow into engine based on air flow and mixture lever position
448 //assume lever runs from no flow at fully out to thi = 1.6 at fully in at sea level
449 //also assume that the injector linkage is ideal - hence the set mixture is maintained at a given altitude throughout the speed and load range
450 thi_sea_level = 1.6 * ( Mixture_Lever_Pos / 100.0 );
451 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
452 m_dot_fuel = m_dot_air / 14.7 * equivalence_ratio;
454 // cout << "fuel " << m_dot_fuel;
455 // cout << " air " << m_dot_air << '\n';
459 // cout << "Thi = " << equivalence_ratio << '\n';
461 combustion_efficiency = Lookup_Combustion_Efficiency(equivalence_ratio); //The combustion efficiency basically tells us what proportion of the fuels calorific value is released
463 // cout << "Combustion efficiency = " << combustion_efficiency << '\n';
465 //now calculate energy release to exhaust
466 //We will assume a three way split of fuel energy between useful work, the coolant system and the exhaust system
467 //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
468 //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.
469 enthalpy_exhaust = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
470 heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
471 delta_T_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
472 // delta_T_exhaust = Calculate_Delta_T_Exhaust();
474 // cout << "T_amb " << T_amb;
475 // cout << " dT exhaust = " << delta_T_exhaust;
477 EGT = T_amb + delta_T_exhaust;
479 // cout << " EGT = " << EGT << '\n';
481 //***************************************************************************************
482 // Calculate Cylinder Head Temperature
486 This is a somewhat rough first attempt at modelling cylinder head temperature. The cylinder head
487 is assumed to be at uniform temperature. Obviously this is incorrect, but it simplifies things a
488 lot, and we're just looking for the behaviour of CHT to be correct. Energy transfer to the cylinder
489 head is assumed to be one third of the energy released by combustion at all conditions. This is a
490 reasonable estimate, although obviously in real life it varies with different conditions and possibly
491 with CHT itself. I've split energy transfer from the cylinder head into 2 terms - free convection -
492 ie convection to stationary air, and forced convection, ie convection into flowing air. The basic
493 free convection equation is: dqdt = -hAdT Since we don't know A and are going to set h quite arbitarily
494 anyway I've knocked A out and just wrapped it up in h - the only real significance is that the units
495 of h will be different but that dosn't really matter to us anyway. In addition, we have the problem
496 that the prop model I'm currently using dosn't model the backwash from the prop which will add to the
497 velocity of the cooling air when the prop is turning, so I've added an extra term to try and cope
500 In real life, forced convection equations are genarally empirically derived, and are quite complicated
501 and generally contain such things as the Reynolds and Nusselt numbers to various powers. The best
502 course of action would probably to find an empirical correlation from the literature for a similar
503 situation and try and get it to fit well. However, for now I am using my own made up very simple
504 correlation for the energy transfer from the cylinder head:
506 dqdt = -(h1.dT) -(h2.m_dot.dT) -(h3.rpm.dT)
508 where dT is the temperature different between the cylinder head and the surrounding air, m_dot is the
509 mass flow rate of cooling air through an arbitary volume, rpm is the engine speed in rpm (this is the
510 backwash term), and h1, h2, h3 are co-efficients which we can play with to attempt to get the CHT
511 behaviour to match real life.
513 In order to change the values of CHT that the engine settles down at at various conditions,
514 have a play with h1, h2 and h3. In order to change the rate of heating/cooling without affecting
515 equilibrium values alter the cylinder head mass, which is really quite arbitary. Bear in mind that
516 altering h1, h2 and h3 will also alter the rate of heating or cooling as well as equilibrium values,
517 but altering the cylinder head mass will only alter the rate. It would I suppose be better to read
518 the values from file to avoid the necessity for re-compilation every time I change them.
521 //CHT = Calc_CHT( Fuel_Flow, Mixture, IAS);
522 // cout << "Cylinder Head Temp (F) = " << CHT << endl;
523 float h1 = -95.0; //co-efficient for free convection
524 float h2 = -3.95; //co-efficient for forced convection
525 float h3 = -0.05; //co-efficient for forced convection due to prop backwash
526 float v_apparent; //air velocity over cylinder head in m/s
527 float v_dot_cooling_air;
528 float m_dot_cooling_air;
529 float temperature_difference;
530 float arbitary_area = 1.0;
531 float dqdt_from_combustion;
532 float dqdt_forced; //Rate of energy transfer to/from cylinder head due to forced convection (Joules) (sign convention: to cylinder head is +ve)
533 float dqdt_free; //Rate of energy transfer to/from cylinder head due to free convection (Joules) (sign convention: to cylinder head is +ve)
534 float dqdt_cylinder_head; //Overall energy change in cylinder head
535 float CpCylinderHead = 800.0; //FIXME - this is a guess - I need to look up the correct value
536 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
537 float HeatCapacityCylinderHead;
540 temperature_difference = CHT - T_amb;
542 v_apparent = IAS * 0.5144444; //convert from knots to m/s
543 v_dot_cooling_air = arbitary_area * v_apparent;
544 m_dot_cooling_air = v_dot_cooling_air * rho_air;
546 //Calculate rate of energy transfer to cylinder head from combustion
547 dqdt_from_combustion = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
549 //Calculate rate of energy transfer from cylinder head due to cooling NOTE is calculated as rate to but negative
550 dqdt_forced = (h2 * m_dot_cooling_air * temperature_difference) + (h3 * RPM * temperature_difference);
551 dqdt_free = h1 * temperature_difference;
553 //Calculate net rate of energy transfer to or from cylinder head
554 dqdt_cylinder_head = dqdt_from_combustion + dqdt_forced + dqdt_free;
556 HeatCapacityCylinderHead = CpCylinderHead * MassCylinderHead;
558 dCHTdt = dqdt_cylinder_head / HeatCapacityCylinderHead;
560 CHT += (dCHTdt * time_step);
562 CHT_degF = (CHT * 1.8) - 459.67;
564 // cout << "CHT = " << CHT_degF << " degF\n";
567 // End calculate Cylinder Head Temperature
570 //***************************************************************************************
571 // Engine Power & Torque Calculations
574 // For a given Manifold Pressure and RPM calculate the % Power
575 // Multiply Manifold Pressure by RPM
576 ManXRPM = Manifold_Pressure * RPM;
581 // Phil's %power correlation
583 Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM) + ( + 7E-04 * ManXRPM) - 0.1218;
584 // cout << Percentage_Power << "%" << "\t";
587 // DCL %power correlation - basically Phil's correlation modified to give slighty less power at the low end
588 // might need some adjustment as the prop model is adjusted
589 // 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
591 Percentage_Power = (+ 6E-09 * ManXRPM * ManXRPM) + ( + 8E-04 * ManXRPM) - 1.8524;
592 // cout << Percentage_Power << "%" << "\t";
595 // Adjust for Temperature - Temperature above Standard decrease
596 // power % by 7/120 per degree F increase, and incease power for
597 // temps below at the same ratio
598 Percentage_Power = Percentage_Power - (FG_ISA_VAR * 7 /120);
599 // cout << Percentage_Power << "%" << "\t";
601 // Adjust for Altitude. In this version a linear variation is
602 // used. Decrease 1% for each 1000' increase in Altitde
603 Percentage_Power = Percentage_Power + (FG_Pressure_Ht * 12/10000);
604 // cout << Percentage_Power << "%" << "\t";
607 //DCL - now adjust power to compensate for mixture
608 //uses a curve fit to the data in the IO360 / O360 operating manual
609 //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,
610 //possibly by using separate fits for rich and lean of best power mixture
611 //first adjust actual mixture to abstract mixture - this is a temporary hack in order to account for the fact that the data I have
612 //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.
614 abstract_mixture = 10.0 * equivalence_ratio - 12.0;
615 float m = abstract_mixture; //to simplify writing the next equation
616 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);
617 Percentage_Power = Percentage_Power * Percentage_of_best_power_mixture_power / 100.0;
620 // Now Calculate Fuel Flow based on % Power Best Power Mixture
621 Fuel_Flow = Percentage_Power * Max_Fuel_Flow / 100.0;
622 // cout << Fuel_Flow << " lbs/hr"<< endl;
624 // Now Derate engine for the effects of Bad/Switched off magnetos
625 if (Magneto_Left == 0 && Magneto_Right == 0) {
626 // cout << "Both OFF\n";
627 Percentage_Power = 0;
628 } else if (Magneto_Left && Magneto_Right) {
629 // cout << "Both On ";
630 } else if (Magneto_Left == 0 || Magneto_Right== 0) {
631 // cout << "1 Magneto Failed ";
633 Percentage_Power = Percentage_Power *
634 ((100.0 - Mag_Derate_Percent)/100.0);
635 // cout << FGEng1_Percentage_Power << "%" << "\t";
638 // Calculate Engine Horsepower
640 HP = Percentage_Power * MaxHP / 100.0;
642 Power_SI = HP * CONVERT_HP_TO_WATTS;
644 // Calculate Engine Torque
646 Torque = HP * 5252 / RPM;
647 // cout << Torque << "Ft/lbs" << "\t";
649 Torque_SI = (Power_SI * 60.0) / (2.0 * PI * RPM); //Torque = power / angular velocity
650 // cout << Torque << " Nm\n";
653 // Calculate Oil Pressure
654 Oil_Pressure = Oil_Press( Oil_Temp, RPM );
655 // cout << "Oil Pressure (PSI) = " << Oil_Pressure << endl;
657 //==============================================================
659 // Now do the Propellor Calculations
661 #ifdef PHILS_PROP_MODEL
664 FGProp1_RPS = RPM * Gear_Ratio / 60.0;
665 // cout << FGProp1_RPS << " RPS" << endl;
667 //Radial Flow Vector (V2) Ft/sec at Ref Blade Station (usually 30")
668 FGProp1_Angular_V = FGProp1_RPS * 2 * PI * (Blade_Station / 12);
669 // cout << FGProp1_Angular_V << "Angular Velocity " << endl;
671 // Axial Flow Vector (Vo) Ft/sec
672 // Some further work required here to allow for inflow at low speeds
673 // Vo = (IAS + 20) * 1.688888;
674 Vo = IAS_to_FPS(IAS + 20);
675 // cout << "Feet/sec = " << Vo << endl;
677 // cout << Vo << "Axial Velocity" << endl;
679 // Relative Velocity (V1)
680 V1 = sqrt((FGProp1_Angular_V * FGProp1_Angular_V) +
682 // cout << V1 << "Relative Velocity " << endl;
684 // cout << FGProp1_Blade_Angle << " Prop Blade Angle" << endl;
686 // Blade Angle of Attack (Alpha1)
688 /* cout << " Alpha1 = " << Alpha1
689 << " Blade angle = " << FGProp1_Blade_Angle
691 << " FGProp1_Angular_V = " << FGProp1_Angular_V << endl;*/
692 Alpha1 = FGProp1_Blade_Angle -(atan(Vo / FGProp1_Angular_V) * (180/PI));
693 // cout << Alpha1 << " Alpha1" << endl;
695 // Calculate Coefficient of Drag at Alpha1
696 FGProp1_Coef_Drag = (0.0005 * (Alpha1 * Alpha1)) + (0.0003 * Alpha1)
698 // cout << FGProp1_Coef_Drag << " Coef Drag" << endl;
700 // Calculate Coefficient of Lift at Alpha1
701 FGProp1_Coef_Lift = -(0.0026 * (Alpha1 * Alpha1)) + (0.1027 * Alpha1)
703 // cout << FGProp1_Coef_Lift << " Coef Lift " << endl;
705 // Covert Alplha1 to Radians
706 // Alpha1 = Alpha1 * PI / 180;
708 // Calculate Prop Torque
709 FGProp1_Torque = (0.5 * Rho * (V1 * V1) * FGProp_Area
710 * ((FGProp1_Coef_Lift * sin(Alpha1 * PI / 180))
711 + (FGProp1_Coef_Drag * cos(Alpha1 * PI / 180))))
712 * (Blade_Station/12);
713 // cout << FGProp1_Torque << " Prop Torque" << endl;
715 // Calculate Prop Thrust
716 // cout << " V1 = " << V1 << " Alpha1 = " << Alpha1 << endl;
717 FGProp1_Thrust = 0.5 * Rho * (V1 * V1) * FGProp_Area
718 * ((FGProp1_Coef_Lift * cos(Alpha1 * PI / 180))
719 - (FGProp1_Coef_Drag * sin(Alpha1 * PI / 180)));
720 // cout << FGProp1_Thrust << " Prop Thrust " << endl;
722 // End of Propeller Calculations
723 //==============================================================
725 #endif //PHILS_PROP_MODEL
727 #ifdef NEVS_PROP_MODEL
732 number_of_blades = 2.0;
734 allowance_for_spinner = blade_length / 12.0;
735 prop_fudge_factor = 1.453401525;
736 forward_velocity = IAS;
745 angular_velocity_SI = 2.0 * PI * RPM / 60.0;
747 allowance_for_spinner = blade_length / 12.0;
748 //Calculate thrust and torque by summing the contributions from each of the blade elements
749 //Assumes equal length elements with numbered 1 inboard -> num_elements outboard
753 // outfile << "Rho = " << Rho << '\n\n';
754 // outfile << "Drag = ";
755 for(i=1;i<=num_elements;i++)
758 distance = (blade_length * (element / num_elements)) + allowance_for_spinner;
759 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))))))
760 * (0.1 * (blade_length / element)) * number_of_blades;
762 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)
763 * (0.1 * (blade_length / element)) * number_of_blades;
764 element_torque = element_drag * distance;
765 prop_torque += element_torque;
766 prop_thrust += element_lift;
767 // outfile << "Drag = " << element_drag << " n = " << element << '\n';
772 // outfile << "Angular velocity = " << angular_velocity_SI << " rad/s\n";
774 // cout << "Thrust = " << prop_thrust << '\n';
775 prop_thrust *= prop_fudge_factor;
776 prop_torque *= prop_fudge_factor;
777 prop_power_consumed_SI = prop_torque * angular_velocity_SI;
778 prop_power_consumed_HP = prop_power_consumed_SI / 745.699;
781 #endif //NEVS_PROP_MODEL
785 #ifdef PHILS_PROP_MODEL //Do Torque calculations in Ft/lbs - yuk :-(((
786 Torque_Imbalance = FGProp1_Torque - Torque;
788 if (Torque_Imbalance > 5) {
790 // FGProp1_RPM -= 25;
791 // FGProp1_Blade_Angle -= 0.75;
794 if (Torque_Imbalance < -5) {
796 // FGProp1_RPM += 25;
797 // FGProp1_Blade_Angle += 0.75;
802 #ifdef NEVS_PROP_MODEL //use proper units - Nm
803 Torque_Imbalance = Torque_SI - prop_torque; //This gives a +ve value when the engine torque exeeds the prop torque
805 angular_acceleration = Torque_Imbalance / (engine_inertia + prop_inertia);
806 angular_velocity_SI += (angular_acceleration * time_step);
807 RPM = (angular_velocity_SI * 60) / (2.0 * PI);
813 if( RPM > (Desired_RPM + 2)) {
814 FGProp1_Blade_Angle += 0.75; //This value could be altered depending on how far from the desired RPM we are
817 if( RPM < (Desired_RPM - 2)) {
818 FGProp1_Blade_Angle -= 0.75;
821 if (FGProp1_Blade_Angle < FGProp_Fine_Pitch_Stop) {
822 FGProp1_Blade_Angle = FGProp_Fine_Pitch_Stop;
829 //end constant speed prop
832 //DCL - stall the engine if RPM drops below 550 - this is possible if mixture lever is pulled right out
836 // outfile << "RPM = " << RPM << " Blade angle = " << FGProp1_Blade_Angle << " Engine torque = " << Torque << " Prop torque = " << FGProp1_Torque << '\n';
837 // outfile << "RPM = " << RPM << " Engine torque = " << Torque_SI << " Prop torque = " << prop_torque << '\n';
839 // cout << FGEng1_RPM << " Blade_Angle " << FGProp1_Blade_Angle << endl << endl;
843 // cout << "Final engine RPM = " << RPM << '\n';
851 // Calculate Oil Temperature
853 static float Oil_Temp (float Fuel_Flow, float Mixture, float IAS)