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 // ------------------------------------------------------------------------
50 #include <simgear/compiler.h>
60 // ------------------------------------------------------------------------
62 // ------------------------------------------------------------------------
65 // Calculate Engine RPM based on Propellor Lever Position
66 float FGEngine::Calc_Engine_RPM (float LeverPosition)
68 // Calculate RPM as set by Prop Lever Position. Assumes engine
69 // will run at 1000 RPM at full course
71 float RPM = LeverPosition * (Max_RPM - Min_RPM) /100 + Min_RPM ;
73 if ( RPM >= Max_RPM ) {
81 // Calculate Manifold Pressure based on Throttle lever Position
82 static float Calc_Manifold_Pressure ( float LeverPosn, float MaxMan)
88 Inches = LeverPosn * MaxMan / 100;
93 // set initial default values
94 void FGEngine::init() {
95 // Control and environment inputs
97 Throttle_Lever_Pos = 75;
98 Propeller_Lever_Pos = 75;
99 Mixture_Lever_Pos = 100;
101 // Engine Specific Variables used by this program that have limits.
102 // Will be set in a parameter file to be read in to create
103 // and instance for each engine.
104 Max_Manifold_Pressure = 29.50;
108 Mag_Derate_Percent = 5;
112 // Initialise Engine Variables used by this instance
113 Percentage_Power = 0;
114 Manifold_Pressure = 29.00; // Inches
116 Fuel_Flow = 0; // lbs/hour
120 Oil_Pressure = 0; // PSI
121 Oil_Temp = 85; // Deg C
124 Torque_Imbalance = 0;
127 // Initialise Propellor Variables used by this instance
128 FGProp1_Angular_V = 0;
129 FGProp1_Coef_Drag = 0.6;
133 FGProp1_Coef_Lift = 0.1;
135 FGProp1_Blade_Angle = 13.5;
136 FGProp_Fine_Pitch_Stop = 13.5;
137 FGProp_Course_Pitch_Stop = 55;
139 // Other internal values
144 // Calculate Oil Pressure
145 static float Oil_Press (float Oil_Temp, float Engine_RPM)
147 float Oil_Pressure = 0; //PSI
148 float Oil_Press_Relief_Valve = 60; //PSI
149 float Oil_Press_RPM_Max = 1800;
150 float Design_Oil_Temp = 85; //Celsius
151 float Oil_Viscosity_Index = 0.25; // PSI/Deg C
152 float Temp_Deviation = 0; // Deg C
154 Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
156 // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
157 if (Oil_Pressure >= Oil_Press_Relief_Valve)
159 Oil_Pressure = Oil_Press_Relief_Valve;
162 // Now adjust pressure according to Temp which affects the viscosity
164 Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
170 // Calculate Cylinder Head Temperature
171 static float Calc_CHT (float Fuel_Flow, float Mixture, float IAS)
179 // Calculate Density Ratio
180 static float Density_Ratio ( float x )
182 float y = ((3E-10 * x * x) - (3E-05 * x) + 0.9998);
187 // Calculate Air Density - Rho
188 static float Density ( float x )
190 float y = ((9E-08 * x * x) - (7E-08 * x) + 0.0024);
195 // Calculate Speed in FPS given Knots CAS
196 static float IAS_to_FPS (float ias)
198 return ias * 1.68888888;
202 // update the engine model based on current control positions
203 void FGEngine::update() {
204 // Declare local variables
205 int num = 0; // Not used. Counting variables
206 int num2 = 100; // Not used.
211 // Set up the new variables
212 float Blade_Station = 30;
213 float Rho = 0.002378;
214 float FGProp_Area = 1.405/3;
215 float PI = 3.1428571;
220 // 0 = Closed, 100 = Fully Open
221 // float Throttle_Lever_Pos = 75;
222 // 0 = Full Course 100 = Full Fine
223 // float Propeller_Lever_Pos = 75;
224 // 0 = Idle Cut Off 100 = Full Rich
225 // float Mixture_Lever_Pos = 100;
227 // Environmental Variables
229 // Temp Variation from ISA (Deg F)
230 float FG_ISA_VAR = 0;
231 // Pressure Altitude 1000's of Feet
232 float FG_Pressure_Ht = 0;
234 // Parameters that alter the operation of the engine.
235 // Yes = 1. Is there Fuel Available. Calculated elsewhere
236 int Fuel_Available = 1;
237 // Off = 0. Reduces power by 3 % for same throttle setting
238 int Alternate_Air_Pos =0;
239 // 1 = On. Reduces power by 5 % for same power lever settings
240 int Magneto_Left = 1;
241 // 1 = On. Ditto, Both of the above though do not alter fuel flow
242 int Magneto_Right = 1;
244 // There needs to be a section in here to trap silly values, like
245 // 0, otherwise they will crash the calculations
247 // cout << " Number of Iterations ";
251 // cout << " Throttle % ";
252 // cin >> Throttle_Lever_Pos;
255 // cout << " Prop % ";
256 // cin >> Propeller_Lever_Pos;
259 //==================================================================
260 // Engine & Environmental Inputs from elsewhere
262 // Calculate Air Density (Rho) - In FG this is calculated in
265 Rho = Density(FG_Pressure_Ht); // In FG FG_Pressure_Ht is "h"
266 // cout << "Rho = " << Rho << endl;
268 // Calculate Manifold Pressure (Engine 1) as set by throttle opening
271 Calc_Manifold_Pressure( Throttle_Lever_Pos, Max_Manifold_Pressure );
272 cout << "manifold pressure = " << Manifold_Pressure << endl;
275 RPM = Calc_Engine_RPM(Propeller_Lever_Pos);
276 // cout << "Engine RPM = " << RPM << endl;
279 cout << "Desired RPM = " << Desired_RPM << endl;
281 //==================================================================
282 // Engine Power & Torque Calculations
284 // Loop until stable - required for testing only
285 for (num = 0; num < num2; num++) {
286 // cout << endl << "====================" << endl;
287 // cout << "MP Inches = " << Manifold_Pressure << "\t";
288 // cout << " RPM = " << RPM << "\t";
290 // For a given Manifold Pressure and RPM calculate the % Power
291 // Multiply Manifold Pressure by RPM
292 ManXRPM = Manifold_Pressure * RPM;
293 // cout << ManXRPM << endl;
296 Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM)
297 + ( + 7E-04 * ManXRPM) - 0.1218;
298 // cout << "percent power = " << Percentage_Power << "%" << "\t";
300 // Adjust for Temperature - Temperature above Standard decrease
301 // power % by 7/120 per degree F increase, and incease power for
302 // temps below at the same ratio
303 Percentage_Power = Percentage_Power - (FG_ISA_VAR * 7 /120);
304 // cout << " adjusted T = " << Percentage_Power << "%" << "\t";
306 // Adjust for Altitude. In this version a linear variation is
307 // used. Decrease 1% for each 1000' increase in Altitde
308 Percentage_Power = Percentage_Power + (FG_Pressure_Ht * 12/10000);
309 // cout << " adjusted A = " << Percentage_Power << "%" << "\t";
311 // Now Calculate Fuel Flow based on % Power Best Power Mixture
312 Fuel_Flow = Percentage_Power * Max_Fuel_Flow / 100.0;
313 // cout << Fuel_Flow << " lbs/hr"<< endl;
315 // Now Derate engine for the effects of Bad/Switched off magnetos
316 if (Magneto_Left == 0 && Magneto_Right == 0) {
317 // cout << "Both OFF\n";
318 Percentage_Power = 0;
319 } else if (Magneto_Left && Magneto_Right) {
320 // cout << "Both On ";
321 } else if (Magneto_Left == 0 || Magneto_Right== 0) {
322 // cout << "1 Magneto Failed ";
324 Percentage_Power = Percentage_Power *
325 ((100.0 - Mag_Derate_Percent)/100.0);
327 // cout << "Final engine % power = " << Percentage_Power << "%" << endl;
329 // Calculate Engine Horsepower
331 HP = Percentage_Power * MaxHP / 100.0;
333 // Calculate Engine Torque
335 Torque = HP * 5252 / RPM;
336 // cout << Torque << "Ft/lbs" << "\t";
338 // Calculate Cylinder Head Temperature
339 CHT = Calc_CHT( Fuel_Flow, Mixture, IAS);
340 // cout << "Cylinder Head Temp (F) = " << CHT << endl;
342 // Calculate Oil Pressure
343 Oil_Pressure = Oil_Press( Oil_Temp, RPM );
344 // cout << "Oil Pressure (PSI) = " << Oil_Pressure << endl;
346 //==============================================================
348 // Now do the Propellor Calculations
351 FGProp1_RPS = RPM * Gear_Ratio / 60.0;
352 // cout << FGProp1_RPS << " RPS" << endl;
354 //Radial Flow Vector (V2) Ft/sec at Ref Blade Station (usually 30")
355 FGProp1_Angular_V = FGProp1_RPS * 2 * PI * (Blade_Station / 12);
356 // cout << "Angular Velocity " << FGProp1_Angular_V << endl;
358 // Axial Flow Vector (Vo) Ft/sec
359 // Some further work required here to allow for inflow at low speeds
360 // Vo = (IAS + 20) * 1.688888;
361 Vo = IAS_to_FPS(IAS + 20);
362 // cout << "Feet/sec = " << Vo << endl;
364 // cout << Vo << "Axial Velocity" << endl;
366 // Relative Velocity (V1)
367 V1 = sqrt((FGProp1_Angular_V * FGProp1_Angular_V) +
369 // cout << "Relative Velocity " << V1 << endl;
371 if ( FGProp1_Blade_Angle >= FGProp_Course_Pitch_Stop ) {
372 FGProp1_Blade_Angle = FGProp_Course_Pitch_Stop;
375 // cout << FGProp1_Blade_Angle << " Prop Blade Angle" << endl;
377 // Blade Angle of Attack (Alpha1)
379 Alpha1 = FGProp1_Blade_Angle -(atan(Vo / FGProp1_Angular_V) * (180/PI));
380 // cout << Alpha1 << " Alpha1" << endl;
382 // cout << " Alpha1 = " << Alpha1
383 // << " Blade angle = " << FGProp1_Blade_Angle
385 // << " FGProp1_Angular_V = " << FGProp1_Angular_V << endl;
387 // Calculate Coefficient of Drag at Alpha1
388 FGProp1_Coef_Drag = (0.0005 * (Alpha1 * Alpha1)) + (0.0003 * Alpha1)
390 // cout << FGProp1_Coef_Drag << " Coef Drag" << endl;
392 // Calculate Coefficient of Lift at Alpha1
393 FGProp1_Coef_Lift = -(0.0026 * (Alpha1 * Alpha1)) + (0.1027 * Alpha1)
395 // cout << FGProp1_Coef_Lift << " Coef Lift " << endl;
397 // Covert Alplha1 to Radians
398 // Alpha1 = Alpha1 * PI / 180;
400 // Calculate Prop Torque
401 FGProp1_Torque = (0.5 * Rho * (V1 * V1) * FGProp_Area
402 * ((FGProp1_Coef_Lift * sin(Alpha1 * PI / 180))
403 + (FGProp1_Coef_Drag * cos(Alpha1 * PI / 180))))
404 * (Blade_Station/12);
405 // cout << "Prop Torque = " << FGProp1_Torque << endl;
407 // Calculate Prop Thrust
408 // cout << " V1 = " << V1 << " Alpha1 = " << Alpha1 << endl;
409 FGProp1_Thrust = 0.5 * Rho * (V1 * V1) * FGProp_Area
410 * ((FGProp1_Coef_Lift * cos(Alpha1 * PI / 180))
411 - (FGProp1_Coef_Drag * sin(Alpha1 * PI / 180)));
412 // cout << " Prop Thrust = " << FGProp1_Thrust << endl;
414 // End of Propeller Calculations
415 //==============================================================
419 Torque_Imbalance = FGProp1_Torque - Torque;
420 // cout << Torque_Imbalance << endl;
422 if (Torque_Imbalance > 20) {
424 // FGProp1_RPM -= 25;
425 FGProp1_Blade_Angle -= 0.75;
428 if (FGProp1_Blade_Angle < FGProp_Fine_Pitch_Stop) {
429 FGProp1_Blade_Angle = FGProp_Fine_Pitch_Stop;
431 if (Torque_Imbalance < -20) {
433 // FGProp1_RPM += 25;
434 FGProp1_Blade_Angle += 0.75;
442 // cout << FGEng1_RPM << " Blade_Angle " << FGProp1_Blade_Angle << endl << endl;
452 // Calculate Oil Temperature
454 static float Oil_Temp (float Fuel_Flow, float Mixture, float IAS)