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 (philip@zedley.com)
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
46 // 19/8/2000 PLS Updated E-mail address - This version compiles
47 // 19/8/2000 PLS Set Max Prop blade angle to prevent prop exeeding 90
48 // ------------------------------------------------------------------------
50 // ------------------------------------------------------------------------
52 #include <simgear/compiler.h>
58 #if !defined(SG_HAVE_NATIVE_SGI_COMPILERS)
63 // ------------------------------------------------------------------------
65 // ------------------------------------------------------------------------
67 // prototype definitions
68 // These should be in a header file 10520c.h
70 float Density (float x);
71 void ShowRho (float x);
73 float IAS_to_FPS (float x);
74 void ShowFPS(float x);
76 float Get_Throttle (float x);
77 void Show_Throttle (float x);
79 float Manifold_Pressure (float x, float z);
80 void Show_Manifold_Pressure (float x);
82 float CHT (float Fuel_Flow, float Mixture, float IAS);
83 void Show_CHT (float x);
85 float Oil_Temp (float x, float y);
86 void Show_Oil_Temp (float x);
88 float Oil_Press (float Oil_Temp, float Engine_RPM);
89 void Show_Oil_Press (float x);
94 // Declare local variables
95 int num = 0; // Not used. Counting variables
96 int num2 = 100; // Not used.
102 // Set up the new variables
103 float Blade_Station = 30;
104 float Rho = 0.002378;
105 float FGProp_Area = 1.405/3;
106 float PI = 3.1428571;
110 cout << "Enter IAS ";
116 // 0 = Closed, 100 = Fully Open
117 float FGEng1_Throttle_Lever_Pos = 75;
118 // 0 = Full Course 100 = Full Fine
119 float FGEng1_Propeller_Lever_Pos = 75;
120 // 0 = Idle Cut Off 100 = Full Rich
121 float FGEng1_Mixture_Lever_Pos = 100;
123 // Environmental Variables
125 // Temp Variation from ISA (Deg F)
126 float FG_ISA_VAR = 0;
127 // Pressure Altitude 1000's of Feet
128 float FG_Pressure_Ht = 0;
130 // Parameters that alter the operation of the engine.
131 // Yes = 1. Is there Fuel Available. Calculated elsewhere
132 int FGEng1_Fuel_Available = 1;
133 // Off = 0. Reduces power by 3 % for same throttle setting
134 int FGEng1_Alternate_Air_Pos =0;
135 // 1 = On. Reduces power by 5 % for same power lever settings
136 int FGEng1_Magneto_Left = 1;
137 // 1 = On. Ditto, Both of the above though do not alter fuel flow
138 int FGEng1_Magneto_Right = 1;
140 // There needs to be a section in here to trap silly values, like
141 // 0, otherwise they will crash the calculations
143 // Engine Specific Variables used by this program that have limits.
144 // Will be set in a parameter file to be read in to create
145 // and instance for each engine.
146 float FGEng_Max_Manifold_Pressure = 29.50;
147 float FGEng_Max_RPM = 2700;
148 float FGEng_Min_RPM = 1000;
149 float FGEng_Max_Fuel_Flow = 130;
150 float FGEng_Mag_Derate_Percent = 5;
151 float FGEng_MaxHP = 285;
152 float FGEng_Gear_Ratio = 1;
154 // Initialise Engine Variables used by this instance
155 float FGEng1_Percentage_Power = 0;
156 float FGEng1_Manifold_Pressure = 29.00; // Inches
157 float FGEng1_RPM = 500;
158 float FGEng1_Fuel_Flow = 0; // lbs/hour
159 float FGEng1_Torque = 0;
160 float FGEng1_CHT = 370;
161 float FGEng1_Mixture = 14;
162 float FGEng1_Oil_Pressure = 0; // PSI
163 float FGEng1_Oil_Temp = 85; // Deg C
165 float FGEng1_RPS = 0;
166 float Torque_Imbalance = 0;
167 float FGEng1_Desired_RPM = 0;
169 // Initialise Propellor Variables used by this instance
170 float FGProp1_Angular_V = 0;
171 float FGProp1_Coef_Drag = 0.6;
172 float FGProp1_Torque = 0;
173 float FGProp1_Thrust = 0;
174 float FGProp1_RPS = 0;
175 float FGProp1_Coef_Lift = 0.1;
177 float FGProp1_Blade_Angle = 13.5;
178 float FGProp_Fine_Pitch_Stop = 13.5;
179 float FGProp_Course_Pitch_Stop = 55;
181 // cout << "Enter Blade Angle ";
182 // cin >> FGProp1_Blade_Angle;
185 cout << " Number of Iterations ";
190 cout << " Throttle % ";
191 // cin >> FGEng1_Throttle_Lever_Pos;
192 FGEng1_Throttle_Lever_Pos = 50;
196 // cin >> FGEng1_Propeller_Lever_Pos;
197 FGEng1_Propeller_Lever_Pos = 100;
200 //==================================================================
201 // Engine & Environmental Inputs from elsewhere
203 // Calculate Air Density (Rho) - In FG this is calculated in
206 Rho = Density(FG_Pressure_Ht); // In FG FG_Pressure_Ht is "h"
210 // Calculate Manifold Pressure (Engine 1) as set by throttle opening
212 FGEng1_Manifold_Pressure = Manifold_Pressure(FGEng1_Throttle_Lever_Pos,
213 FGEng1_Manifold_Pressure );
214 Show_Manifold_Pressure(FGEng1_Manifold_Pressure);
216 // Calculate Desired RPM as set by Prop Lever Position.
217 // Actual engine RPM may be different
218 // The governed max RPM at 100% Prop Lever Position = FGEng_MaxRPM
219 // The governed minimum RPM at 0% Prop Lever Position = FGEng_Min_RPM
220 // The actual minimum RPM of the engine can be < FGEng_Min_RPM if there is insufficient
221 // engine torque to counter act the propeller torque at FGProp_Fine_Pitch_Stop
223 FGEng1_RPM = (FGEng1_Propeller_Lever_Pos * (FGEng_Max_RPM - FGEng_Min_RPM) /100)
226 // * ((FGEng_Max_RPM + FGEng_Min_RPM) / 100);
228 if (FGEng1_RPM >= 2700) {
231 FGEng1_Desired_RPM = FGEng1_RPM;
233 cout << "Desired RPM = " << FGEng1_Desired_RPM << endl;
235 //==================================================================
236 // Engine Power & Torque Calculations
238 // Loop until stable - required for testing only
239 for (num = 1; num < num2; num++) {
240 cout << endl << "====================" << endl;
241 cout << "MP Inches = " << FGEng1_Manifold_Pressure << "\t";
242 cout << FGEng1_RPM << " RPM" << "\t";
244 // For a givem Manifold Pressure and RPM calculate the % Power
245 // Multiply Manifold Pressure by RPM
246 ManXRPM = FGEng1_Manifold_Pressure * FGEng1_RPM;
247 cout << ManXRPM << endl;
250 FGEng1_Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM)
251 + ( + 7E-04 * ManXRPM) - 0.1218;
252 cout << "percent power = " << FGEng1_Percentage_Power << "%" << "\t";
254 // Adjust for Temperature - Temperature above Standard decrease
255 // power % by 7/120 per degree F increase, and incease power for
256 // temps below at the same ratio
257 FGEng1_Percentage_Power = FGEng1_Percentage_Power - (FG_ISA_VAR * 7 /120);
258 cout << " adjusted T = " << FGEng1_Percentage_Power << "%" << "\t";
260 // Adjust for Altitude. In this version a linear variation is
261 // used. Decrease 1% for each 1000' increase in Altitde
262 FGEng1_Percentage_Power = FGEng1_Percentage_Power
263 + (FG_Pressure_Ht * 12/10000);
264 cout << " adjusted A = " << FGEng1_Percentage_Power << "%" << "\t";
266 // Now Calculate Fuel Flow based on % Power Best Power Mixture
267 FGEng1_Fuel_Flow = FGEng1_Percentage_Power
268 * FGEng_Max_Fuel_Flow / 100;
269 // cout << FGEng1_Fuel_Flow << " lbs/hr"<< endl;
271 // Now Derate engine for the effects of Bad/Switched off magnetos
272 if (FGEng1_Magneto_Left == 0 && FGEng1_Magneto_Right == 0) {
273 // cout << "Both OFF\n";
274 FGEng1_Percentage_Power = 0;
275 } else if (FGEng1_Magneto_Left && FGEng1_Magneto_Right) {
276 // cout << "Both On ";
277 } else if (FGEng1_Magneto_Left == 0 || FGEng1_Magneto_Right== 0) {
278 // cout << "1 Magneto Failed ";
280 FGEng1_Percentage_Power = FGEng1_Percentage_Power *
281 ((100 - FGEng_Mag_Derate_Percent)/100);
282 // cout << FGEng1_Percentage_Power << "%" << "\t";
285 // Calculate Engine Horsepower
287 FGEng1_HP = FGEng1_Percentage_Power * FGEng_MaxHP/100;
289 // Calculate Engine Torque
291 FGEng1_Torque = FGEng1_HP * 5252 / FGEng1_RPM;
292 cout << FGEng1_Torque << "Ft/lbs" << "\t";
294 // Calculate Cylinder Head Temperature
295 FGEng1_CHT = CHT (FGEng1_Fuel_Flow, FGEng1_Mixture, IAS);
296 // Show_CHT (FGEng1_CHT);
298 // Calculate Oil Pressure
299 FGEng1_Oil_Pressure = Oil_Press (FGEng1_Oil_Temp, FGEng1_RPM);
300 // Show_Oil_Press(FGEng1_Oil_Pressure);
303 //==============================================================
305 // Now do the Propellor Calculations
308 FGProp1_RPS = FGEng1_RPM * FGEng_Gear_Ratio/60;
309 // cout << FGProp1_RPS << " RPS" << endl;
311 //Radial Flow Vector (V2) Ft/sec at Ref Blade Station (usually 30")
312 FGProp1_Angular_V = FGProp1_RPS * 2 * PI * (Blade_Station / 12);
313 cout << "Angular Velocity " << FGProp1_Angular_V << endl;
315 // Axial Flow Vector (Vo) Ft/sec
316 // Some further work required here to allow for inflow at low speeds
317 // Vo = (IAS + 20) * 1.688888;
318 Vo = IAS_to_FPS(IAS + 20);
321 // cout << Vo << "Axial Velocity" << endl;
323 // Relative Velocity (V1)
324 V1 = sqrt((FGProp1_Angular_V * FGProp1_Angular_V) +
326 cout << "Relative Velocity " << V1 << endl;
328 if ( FGProp1_Blade_Angle >= FGProp_Course_Pitch_Stop ) {
329 FGProp1_Blade_Angle = FGProp_Course_Pitch_Stop;
332 cout << FGProp1_Blade_Angle << " Prop Blade Angle" << endl;
334 // Blade Angle of Attack (Alpha1)
336 Alpha1 = FGProp1_Blade_Angle -(atan(Vo / FGProp1_Angular_V) * (180/PI));
337 // cout << Alpha1 << " Alpha1" << endl;
339 cout << " Alpha1 = " << Alpha1
340 << " Blade angle = " << FGProp1_Blade_Angle
342 << " FGProp1_Angular_V = " << FGProp1_Angular_V << endl;
344 // Calculate Coefficient of Drag at Alpha1
345 FGProp1_Coef_Drag = (0.0005 * (Alpha1 * Alpha1)) + (0.0003 * Alpha1)
347 // cout << FGProp1_Coef_Drag << " Coef Drag" << endl;
349 // Calculate Coefficient of Lift at Alpha1
350 FGProp1_Coef_Lift = -(0.0026 * (Alpha1 * Alpha1)) + (0.1027 * Alpha1)
352 // cout << FGProp1_Coef_Lift << " Coef Lift " << endl;
354 // Covert Alplha1 to Radians
355 // Alpha1 = Alpha1 * PI / 180;
357 // Calculate Prop Torque
358 FGProp1_Torque = (0.5 * Rho * (V1 * V1) * FGProp_Area
359 * ((FGProp1_Coef_Lift * sin(Alpha1 * PI / 180))
360 + (FGProp1_Coef_Drag * cos(Alpha1 * PI / 180))))
361 * (Blade_Station/12);
362 cout << "Prop Torque = " << FGProp1_Torque << endl;
364 // Calculate Prop Thrust
365 FGProp1_Thrust = 0.5 * Rho * (V1 * V1) * FGProp_Area
366 * ((FGProp1_Coef_Lift * cos(Alpha1 * PI / 180))
367 - (FGProp1_Coef_Drag * sin(Alpha1 * PI / 180)));
368 cout << " Prop Thrust = " << FGProp1_Thrust << endl;
370 // End of Propeller Calculations
371 //==============================================================
375 Torque_Imbalance = FGProp1_Torque - FGEng1_Torque;
376 // cout << Torque_Imbalance << endl;
378 if (Torque_Imbalance > 20) {
380 // FGProp1_RPM -= 25;
381 FGProp1_Blade_Angle -= 0.75;
384 if (FGProp1_Blade_Angle < FGProp_Fine_Pitch_Stop) {
385 FGProp1_Blade_Angle = FGProp_Fine_Pitch_Stop;
387 if (Torque_Imbalance < -20) {
389 // FGProp1_RPM += 25;
390 FGProp1_Blade_Angle += 0.75;
393 if (FGEng1_RPM >= 2700) {
398 // cout << FGEng1_RPM << " Blade_Angle " << FGProp1_Blade_Angle << endl << endl;
411 // Calculate Air Density - Rho
412 float Density ( float x )
415 y = ((9E-08 * x * x) - (7E-08 * x) + 0.0024);
419 // Show Air Density Calculations
420 void ShowRho (float x)
430 // Calculate Speed in FPS given Knots CAS
431 float IAS_to_FPS (float x)
438 // Show Feet per Second
439 void ShowFPS (float x)
441 cout << "Feet/sec = ";
447 // Calculate Manifold Pressure based on Throttle lever Position
449 float Manifold_Pressure ( float x, float z)
460 // Show Manifold Pressure
461 void Show_Manifold_Pressure (float x)
463 cout << "Manifold Pressure = ";
467 // Calculate Oil Temperature
469 float Oil_Temp (float Fuel_Flow, float Mixture, float IAS)
476 // Show Oil Temperature
478 void Show_Oil_Temp (float x)
480 cout << "Oil Temperature (F) = ";
485 // Calculate Oil Pressure
487 float Oil_Press (float Oil_Temp, float Engine_RPM)
489 float Oil_Pressure = 0; //PSI
490 float Oil_Press_Relief_Valve = 60; //PSI
491 float Oil_Press_RPM_Max = 1800;
492 float Design_Oil_Temp = 85; //Celsius
493 float Oil_Viscosity_Index = 0.25; // PSI/Deg C
494 float Temp_Deviation = 0; // Deg C
496 Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
498 // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
499 if (Oil_Pressure >= Oil_Press_Relief_Valve)
501 Oil_Pressure = Oil_Press_Relief_Valve;
504 // Now adjust pressure according to Temp which affects the viscosity
506 Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
508 return (Oil_Pressure);
512 void Show_Oil_Press (float x)
514 cout << "Oil Pressure (PSI) = ";
520 // Calculate Cylinder Head Temperature
522 float CHT (float Fuel_Flow, float Mixture, float IAS)
529 // Show Cyl Head Temperature
531 void Show_CHT (float x)
533 cout << "CHT (F) = ";