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 ) {
107 float FGNewEngine::Lookup_Combustion_Efficiency(float thi_actual)
109 float thi[11]; //array of equivalence ratio values
110 float neta_comb[11]; //corresponding array of combustion efficiency values
111 float neta_comb_actual;
114 //thi = (0.0,0.9,1.0,1.05,1.1,1.15,1.2,1.3,1.4,1.5,1.6);
118 thi[3] = 1.05; //There must be an easier way of doing this !!!!!!!!
126 //neta_comb = (0.98,0.98,0.97,0.95,0.9,0.85,0.79,0.7,0.63,0.57,0.525);
137 neta_comb[10] = 0.525;
138 //combustion efficiency values from Heywood: ISBN 0-07-100499-8
142 j = 11; //This must be equal to the number of elements in the lookup table arrays
148 //this is just to avoid crashing the routine is we are bigger than the last element - for now just return the last element
149 //but at some point we will have to extrapolate further
150 neta_comb_actual = neta_comb[i];
151 return neta_comb_actual;
153 if(thi_actual == thi[i])
155 neta_comb_actual = neta_comb[i];
156 return neta_comb_actual;
158 if((thi_actual > thi[i]) && (thi_actual < thi[i + 1]))
160 //do linear interpolation between the two points
161 factor = (thi_actual - thi[i]) / (thi[i+1] - thi[i]);
162 neta_comb_actual = (factor * (neta_comb[i+1] - neta_comb[i])) + neta_comb[i];
163 return neta_comb_actual;
167 //if we get here something has gone badly wrong
168 cout << "ERROR: error in FGNewEngine::Lookup_Combustion_Efficiency\n";
170 return neta_comb_actual; //keeps the compiler happy
173 float FGNewEngine::Calculate_Delta_T_Exhaust(void)
176 heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
177 dT_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
183 // Calculate Manifold Pressure based on Throttle lever Position
184 static float Calc_Manifold_Pressure ( float LeverPosn, float MaxMan, float MinMan)
191 //Note that setting the manifold pressure as a function of lever position only is not strictly accurate
192 //MAP is also a function of engine speed.
193 Inches = MinMan + (LeverPosn * (MaxMan - MinMan) / 100);
195 //allow for idle bypass valve or slightly open throttle stop
203 // set initial default values
204 void FGNewEngine::init(double dt) {
206 CONVERT_CUBIC_INCHES_TO_METERS_CUBED = 1.638706e-5;
207 // Control and environment inputs
209 Throttle_Lever_Pos = 75;
210 Propeller_Lever_Pos = 75;
211 Mixture_Lever_Pos = 100;
212 Cp_air = 1005; // J/KgK
213 Cp_fuel = 1700; // J/KgK
214 calorific_value_fuel = 47.3e6; // W/Kg Note that this is only an approximate value
218 // Engine Specific Variables used by this program that have limits.
219 // Will be set in a parameter file to be read in to create
220 // and instance for each engine.
221 Max_Manifold_Pressure = 28.50; //Inches Hg. An approximation - should be able to find it in the engine performance data
222 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
224 Min_RPM = 600; //Recommended idle from Continental data sheet
226 Mag_Derate_Percent = 5;
227 // MaxHP = 285; //Continental IO520-M
228 MaxHP = 180; //Lycoming IO360
229 // displacement = 520; //Continental IO520-M
230 displacement = 360; //Lycoming IO360
231 engine_inertia = 0.2; //kgm^2 - value taken from a popular family saloon car engine - need to find an aeroengine value !!!!!
232 prop_inertia = 0.03; //kgm^2 - this value is a total guess - dcl
233 displacement_SI = displacement * CONVERT_CUBIC_INCHES_TO_METERS_CUBED;
239 CONVERT_HP_TO_WATTS = 745.6999;
241 // outfile.open(ios::out|ios::trunc);
243 // Initialise Engine Variables used by this instance
244 Percentage_Power = 0;
245 Manifold_Pressure = 29.00; // Inches
247 Fuel_Flow = 0; // lbs/hour
251 Oil_Pressure = 0; // PSI
252 Oil_Temp = 85; // Deg C
255 Torque_Imbalance = 0;
256 Desired_RPM = 2500; //Recommended cruise RPM from Continental datasheet
258 // Initialise Propellor Variables used by this instance
259 FGProp1_Angular_V = 0;
260 FGProp1_Coef_Drag = 0.6;
264 FGProp1_Coef_Lift = 0.1;
266 FGProp1_Blade_Angle = 13.5;
267 FGProp_Fine_Pitch_Stop = 13.5;
269 // Other internal values
274 // Calculate Oil Pressure
275 static float Oil_Press (float Oil_Temp, float Engine_RPM)
277 float Oil_Pressure = 0; //PSI
278 float Oil_Press_Relief_Valve = 60; //PSI
279 float Oil_Press_RPM_Max = 1800;
280 float Design_Oil_Temp = 85; //Celsius
281 float Oil_Viscosity_Index = 0.25; // PSI/Deg C
282 float Temp_Deviation = 0; // Deg C
284 Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
286 // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
287 if (Oil_Pressure >= Oil_Press_Relief_Valve)
289 Oil_Pressure = Oil_Press_Relief_Valve;
292 // Now adjust pressure according to Temp which affects the viscosity
294 Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
300 // Calculate Cylinder Head Temperature
301 static float Calc_CHT (float Fuel_Flow, float Mixture, float IAS)
309 //Calculate Exhaust Gas Temperature
310 //For now we will simply adjust this as a function of mixture
311 //It may be necessary to consider fuel flow rates and CHT in the calculation in the future
312 static float Calc_EGT (float Mixture)
314 float EGT = 1000; //off the top of my head !!!!
315 //Now adjust for mixture strength
321 // Calculate Density Ratio
322 static float Density_Ratio ( float x )
325 y = ((3E-10 * x * x) - (3E-05 * x) + 0.9998);
330 // Calculate Air Density - Rho
331 static float Density ( float x )
334 y = ((9E-08 * x * x) - (7E-08 * x) + 0.0024);
339 // Calculate Speed in FPS given Knots CAS
340 static float IAS_to_FPS (float x)
348 // update the engine model based on current control positions
349 void FGNewEngine::update() {
350 // Declare local variables
352 // const int num2 = 500; // default is 100, number if iterations to run
353 const int num2 = 5; // default is 100, number if iterations to run
359 // Set up the new variables
360 float Blade_Station = 30;
361 float FGProp_Area = 1.405/3;
362 float PI = 3.1428571;
366 // 0 = Closed, 100 = Fully Open
367 // float Throttle_Lever_Pos = 75;
368 // 0 = Full Course 100 = Full Fine
369 // float Propeller_Lever_Pos = 75;
370 // 0 = Idle Cut Off 100 = Full Rich
371 // float Mixture_Lever_Pos = 100;
373 // Environmental Variables
375 // Temp Variation from ISA (Deg F)
376 float FG_ISA_VAR = 0;
377 // Pressure Altitude 1000's of Feet
378 float FG_Pressure_Ht = 0;
380 // Parameters that alter the operation of the engine.
381 // Yes = 1. Is there Fuel Available. Calculated elsewhere
382 int Fuel_Available = 1;
383 // Off = 0. Reduces power by 3 % for same throttle setting
384 int Alternate_Air_Pos =0;
385 // 1 = On. Reduces power by 5 % for same power lever settings
386 int Magneto_Left = 1;
387 // 1 = On. Ditto, Both of the above though do not alter fuel flow
388 int Magneto_Right = 1;
390 // There needs to be a section in here to trap silly values, like
391 // 0, otherwise they will crash the calculations
393 // cout << " Number of Iterations ";
397 // cout << " Throttle % ";
398 // cin >> Throttle_Lever_Pos;
401 // cout << " Prop % ";
402 // cin >> Propeller_Lever_Pos;
405 //==================================================================
406 // Engine & Environmental Inputs from elsewhere
408 // Calculate Air Density (Rho) - In FG this is calculated in
411 Rho = Density(FG_Pressure_Ht); // In FG FG_Pressure_Ht is "h"
412 // cout << "Rho = " << Rho << endl;
414 // Calculate Manifold Pressure (Engine 1) as set by throttle opening
417 Calc_Manifold_Pressure( Throttle_Lever_Pos, Max_Manifold_Pressure, Min_Manifold_Pressure );
418 // cout << "manifold pressure = " << Manifold_Pressure << endl;
420 //DCL - hack for testing - fly at sea level
423 p_amb_sea_level = 101325;
425 //DCL - next calculate m_dot_air and m_dot_fuel into engine
427 //calculate actual ambient pressure and temperature from altitude
428 //Then find the actual manifold pressure (the calculated one is the sea level pressure)
429 True_Manifold_Pressure = Manifold_Pressure * p_amb / p_amb_sea_level;
431 // RPM = Calc_Engine_RPM(Propeller_Lever_Pos);
433 // cout << "Initial engine RPM = " << RPM << endl;
435 // Desired_RPM = RPM;
439 //DCL - calculate mass air flow into engine based on speed and load - separate this out into a function eventually
440 //t_amb is actual temperature calculated from altitude
441 //calculate density from ideal gas equation
442 rho_air = p_amb / ( R_air * T_amb );
443 rho_air_manifold = rho_air * Manifold_Pressure / 29.6;
444 //calculate ideal engine volume inducted per second
445 swept_volume = (displacement_SI * (RPM / 60)) / 2; //This equation is only valid for a four stroke engine
446 //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
447 volumetric_efficiency = 0.8;
448 //Now use volumetric efficiency to calculate actual air volume inducted per second
449 v_dot_air = swept_volume * volumetric_efficiency;
450 //Now calculate mass flow rate of air into engine
451 m_dot_air = v_dot_air * rho_air_manifold;
453 // cout << "rho air manifold " << rho_air_manifold << '\n';
454 // cout << "Swept volume " << swept_volume << '\n';
458 //DCL - now calculate fuel flow into engine based on air flow and mixture lever position
459 //assume lever runs from no flow at fully out to thi = 1.6 at fully in at sea level
460 //also assume that the injector linkage is ideal - hence the set mixture is maintained at a given altitude throughout the speed and load range
461 thi_sea_level = 1.6 * ( Mixture_Lever_Pos / 100.0 );
462 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
463 m_dot_fuel = m_dot_air / 14.7 * equivalence_ratio;
465 // cout << "fuel " << m_dot_fuel;
466 // cout << " air " << m_dot_air << '\n';
470 // cout << "Thi = " << equivalence_ratio << '\n';
472 combustion_efficiency = Lookup_Combustion_Efficiency(equivalence_ratio); //The combustion efficiency basically tells us what proportion of the fuels calorific value is released
474 // cout << "Combustion efficiency = " << combustion_efficiency << '\n';
476 //now calculate energy release to exhaust
477 //We will assume a three way split of fuel energy between useful work, the coolant system and the exhaust system
478 //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
479 //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.
480 enthalpy_exhaust = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
481 heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
482 delta_T_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
483 // delta_T_exhaust = Calculate_Delta_T_Exhaust();
485 // cout << "T_amb " << T_amb;
486 // cout << " dT exhaust = " << delta_T_exhaust;
488 EGT = T_amb + delta_T_exhaust;
490 // cout << " EGT = " << EGT << '\n';
493 // Calculate Manifold Pressure (Engine 2) as set by throttle opening
495 // FGEng2_Manifold_Pressure = Manifold_Pressure(FGEng2_Throttle_Lever_Pos, FGEng2_Manifold_Pressure);
496 // Show_Manifold_Pressure(FGEng2_Manifold_Pressure);
500 //==================================================================
501 // Engine Power & Torque Calculations
503 // Loop until stable - required for testing only
504 // for (num = 0; num < num2; num++) {
505 // cout << Manifold_Pressure << " Inches" << "\t";
506 // cout << RPM << " RPM" << "\t";
508 // For a given Manifold Pressure and RPM calculate the % Power
509 // Multiply Manifold Pressure by RPM
510 ManXRPM = Manifold_Pressure * RPM;
514 // Phil's %power correlation
515 /* // Calculate % Power
516 Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM)
517 + ( + 7E-04 * ManXRPM) - 0.1218;
518 // cout << Percentage_Power << "%" << "\t"; */
520 // DCL %power correlation - basically Phil's correlation modified to give slighty less power at the low end
521 // might need some adjustment as the prop model is adjusted
522 // 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
524 Percentage_Power = (+ 6E-09 * ManXRPM * ManXRPM)
525 + ( + 8E-04 * ManXRPM) - 1.8524;
526 // cout << Percentage_Power << "%" << "\t";
528 // Adjust for Temperature - Temperature above Standard decrease
529 // power % by 7/120 per degree F increase, and incease power for
530 // temps below at the same ratio
531 Percentage_Power = Percentage_Power - (FG_ISA_VAR * 7 /120);
532 // cout << Percentage_Power << "%" << "\t";
534 // Adjust for Altitude. In this version a linear variation is
535 // used. Decrease 1% for each 1000' increase in Altitde
536 Percentage_Power = Percentage_Power + (FG_Pressure_Ht * 12/10000);
537 // cout << Percentage_Power << "%" << "\t";
539 //DCL - now adjust power to compensate for mixture
540 //uses a curve fit to the data in the IO360 / O360 operating manual
541 //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,
542 //possibly by using separate fits for rich and lean of best power mixture
543 //first adjust actual mixture to abstract mixture - this is a temporary hack in order to account for the fact that the data I have
544 //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.
546 abstract_mixture = 10.0 * equivalence_ratio - 12.0;
547 float m = abstract_mixture; //to simplify writing the next equation
548 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);
549 Percentage_Power = Percentage_Power * Percentage_of_best_power_mixture_power / 100.0;
552 // Now Calculate Fuel Flow based on % Power Best Power Mixture
553 Fuel_Flow = Percentage_Power * Max_Fuel_Flow / 100.0;
554 // cout << Fuel_Flow << " lbs/hr"<< endl;
556 // Now Derate engine for the effects of Bad/Switched off magnetos
557 if (Magneto_Left == 0 && Magneto_Right == 0) {
558 // cout << "Both OFF\n";
559 Percentage_Power = 0;
560 } else if (Magneto_Left && Magneto_Right) {
561 // cout << "Both On ";
562 } else if (Magneto_Left == 0 || Magneto_Right== 0) {
563 // cout << "1 Magneto Failed ";
565 Percentage_Power = Percentage_Power *
566 ((100.0 - Mag_Derate_Percent)/100.0);
567 // cout << FGEng1_Percentage_Power << "%" << "\t";
570 // Calculate Engine Horsepower
572 HP = Percentage_Power * MaxHP / 100.0;
574 Power_SI = HP * CONVERT_HP_TO_WATTS;
576 // Calculate Engine Torque
578 Torque = HP * 5252 / RPM;
579 // cout << Torque << "Ft/lbs" << "\t";
581 Torque_SI = (Power_SI * 60.0) / (2.0 * PI * RPM); //Torque = power / angular velocity
582 // cout << Torque << " Nm\n";
584 // Calculate Cylinder Head Temperature
585 CHT = Calc_CHT( Fuel_Flow, Mixture, IAS);
586 // cout << "Cylinder Head Temp (F) = " << CHT << endl;
588 // EGT = Calc_EGT( Mixture );
590 // Calculate Oil Pressure
591 Oil_Pressure = Oil_Press( Oil_Temp, RPM );
592 // cout << "Oil Pressure (PSI) = " << Oil_Pressure << endl;
594 //==============================================================
596 // Now do the Propellor Calculations
598 #ifdef PHILS_PROP_MODEL
601 FGProp1_RPS = RPM * Gear_Ratio / 60.0;
602 // cout << FGProp1_RPS << " RPS" << endl;
604 //Radial Flow Vector (V2) Ft/sec at Ref Blade Station (usually 30")
605 FGProp1_Angular_V = FGProp1_RPS * 2 * PI * (Blade_Station / 12);
606 // cout << FGProp1_Angular_V << "Angular Velocity " << endl;
608 // Axial Flow Vector (Vo) Ft/sec
609 // Some further work required here to allow for inflow at low speeds
610 // Vo = (IAS + 20) * 1.688888;
611 Vo = IAS_to_FPS(IAS + 20);
612 // cout << "Feet/sec = " << Vo << endl;
614 // cout << Vo << "Axial Velocity" << endl;
616 // Relative Velocity (V1)
617 V1 = sqrt((FGProp1_Angular_V * FGProp1_Angular_V) +
619 // cout << V1 << "Relative Velocity " << endl;
621 // cout << FGProp1_Blade_Angle << " Prop Blade Angle" << endl;
623 // Blade Angle of Attack (Alpha1)
625 /* cout << " Alpha1 = " << Alpha1
626 << " Blade angle = " << FGProp1_Blade_Angle
628 << " FGProp1_Angular_V = " << FGProp1_Angular_V << endl;*/
629 Alpha1 = FGProp1_Blade_Angle -(atan(Vo / FGProp1_Angular_V) * (180/PI));
630 // cout << Alpha1 << " Alpha1" << endl;
632 // Calculate Coefficient of Drag at Alpha1
633 FGProp1_Coef_Drag = (0.0005 * (Alpha1 * Alpha1)) + (0.0003 * Alpha1)
635 // cout << FGProp1_Coef_Drag << " Coef Drag" << endl;
637 // Calculate Coefficient of Lift at Alpha1
638 FGProp1_Coef_Lift = -(0.0026 * (Alpha1 * Alpha1)) + (0.1027 * Alpha1)
640 // cout << FGProp1_Coef_Lift << " Coef Lift " << endl;
642 // Covert Alplha1 to Radians
643 // Alpha1 = Alpha1 * PI / 180;
645 // Calculate Prop Torque
646 FGProp1_Torque = (0.5 * Rho * (V1 * V1) * FGProp_Area
647 * ((FGProp1_Coef_Lift * sin(Alpha1 * PI / 180))
648 + (FGProp1_Coef_Drag * cos(Alpha1 * PI / 180))))
649 * (Blade_Station/12);
650 // cout << FGProp1_Torque << " Prop Torque" << endl;
652 // Calculate Prop Thrust
653 // cout << " V1 = " << V1 << " Alpha1 = " << Alpha1 << endl;
654 FGProp1_Thrust = 0.5 * Rho * (V1 * V1) * FGProp_Area
655 * ((FGProp1_Coef_Lift * cos(Alpha1 * PI / 180))
656 - (FGProp1_Coef_Drag * sin(Alpha1 * PI / 180)));
657 // cout << FGProp1_Thrust << " Prop Thrust " << endl;
659 // End of Propeller Calculations
660 //==============================================================
662 #endif //PHILS_PROP_MODEL
664 #ifdef NEVS_PROP_MODEL
669 number_of_blades = 2.0;
671 allowance_for_spinner = blade_length / 12.0;
672 prop_fudge_factor = 1.453401525;
673 forward_velocity = IAS;
682 angular_velocity_SI = 2.0 * PI * RPM / 60.0;
684 allowance_for_spinner = blade_length / 12.0;
685 //Calculate thrust and torque by summing the contributions from each of the blade elements
686 //Assumes equal length elements with numbered 1 inboard -> num_elements outboard
690 // outfile << "Rho = " << Rho << '\n\n';
691 // outfile << "Drag = ";
692 for(i=1;i<=num_elements;i++)
695 distance = (blade_length * (element / num_elements)) + allowance_for_spinner;
696 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))))))
697 * (0.1 * (blade_length / element)) * number_of_blades;
699 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)
700 * (0.1 * (blade_length / element)) * number_of_blades;
701 element_torque = element_drag * distance;
702 prop_torque += element_torque;
703 prop_thrust += element_lift;
704 // outfile << "Drag = " << element_drag << " n = " << element << '\n';
709 // outfile << "Angular velocity = " << angular_velocity_SI << " rad/s\n";
711 // cout << "Thrust = " << prop_thrust << '\n';
712 prop_thrust *= prop_fudge_factor;
713 prop_torque *= prop_fudge_factor;
714 prop_power_consumed_SI = prop_torque * angular_velocity_SI;
715 prop_power_consumed_HP = prop_power_consumed_SI / 745.699;
718 #endif //NEVS_PROP_MODEL
722 #ifdef PHILS_PROP_MODEL //Do Torque calculations in Ft/lbs - yuk :-(((
723 Torque_Imbalance = FGProp1_Torque - Torque;
725 if (Torque_Imbalance > 5) {
727 // FGProp1_RPM -= 25;
728 // FGProp1_Blade_Angle -= 0.75;
731 if (Torque_Imbalance < -5) {
733 // FGProp1_RPM += 25;
734 // FGProp1_Blade_Angle += 0.75;
739 #ifdef NEVS_PROP_MODEL //use proper units - Nm
740 Torque_Imbalance = Torque_SI - prop_torque; //This gives a +ve value when the engine torque exeeds the prop torque
742 angular_acceleration = Torque_Imbalance / (engine_inertia + prop_inertia);
743 angular_velocity_SI += (angular_acceleration * time_step);
744 RPM = (angular_velocity_SI * 60) / (2.0 * PI);
750 if( RPM > (Desired_RPM + 2)) {
751 FGProp1_Blade_Angle += 0.75; //This value could be altered depending on how far from the desired RPM we are
754 if( RPM < (Desired_RPM - 2)) {
755 FGProp1_Blade_Angle -= 0.75;
758 if (FGProp1_Blade_Angle < FGProp_Fine_Pitch_Stop) {
759 FGProp1_Blade_Angle = FGProp_Fine_Pitch_Stop;
766 //end constant speed prop
769 //DCL - stall the engine if RPM drops below 550 - this is possible if mixture lever is pulled right out
773 // outfile << "RPM = " << RPM << " Blade angle = " << FGProp1_Blade_Angle << " Engine torque = " << Torque << " Prop torque = " << FGProp1_Torque << '\n';
774 // outfile << "RPM = " << RPM << " Engine torque = " << Torque_SI << " Prop torque = " << prop_torque << '\n';
776 // cout << FGEng1_RPM << " Blade_Angle " << FGProp1_Blade_Angle << endl << endl;
780 // cout << "Final engine RPM = " << RPM << '\n';
788 // Calculate Oil Temperature
790 static float Oil_Temp (float Fuel_Flow, float Mixture, float IAS)