src/FDM/10520d.cxx

   1 // Module:        10520c.c
   2 //  Author:       Phil Schubert
   3 //  Date started: 12/03/99
   4 //  Purpose:      Models a Continental IO-520-M Engine
   5 //  Called by:    FGSimExec
   6 //
   7 //  Copyright (C) 1999  Philip L. Schubert (philings@ozemail.com.au)
   8 //
   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.
  13 //
  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.
  18 //
  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
  22 // 02111-1307, USA.
  23 //
  24 // Further information about the GNU General Public License can also
  25 // be found on the world wide web at http://www.gnu.org.
  26 //
  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
  34 //
  35 // ARGUMENTS
  36 // ------------------------------------------------------------------------
  37 //
  38 //
  39 // HISTORY
  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
  44 //                      of Classes
  45 // 15/03/99     PLS     Added Oil Pressure, Oil Temperature and CH Temp
  46 // ------------------------------------------------------------------------
  47 // INCLUDES
  48 // ------------------------------------------------------------------------
  49
  50 #include <simgear/compiler.h>
  51
  52 #include <iostream>
  53 #include <math.h>
  54
  55 #include "10520d.hxx"
  56
  57 SG_USING_STD(cout);
  58 SG_USING_STD(endl);
  59
  60 // ------------------------------------------------------------------------
  61 // CODE
  62 // ------------------------------------------------------------------------
  63
  64
  65 // Calculate Engine RPM based on Propellor Lever Position
  66 float FGEngine::Calc_Engine_RPM (float LeverPosition)
  67 {
  68     // Calculate RPM as set by Prop Lever Position. Assumes engine
  69     // will run at 1000 RPM at full course
  70
  71     float RPM = LeverPosition * (Max_RPM - Min_RPM) /100 + Min_RPM ;
  72
  73     if ( RPM >= Max_RPM ) {
  74         RPM = Max_RPM;
  75     }
  76
  77     return RPM;
  78 }
  79
  80
  81 // Calculate Manifold Pressure based on Throttle lever Position
  82 static float Calc_Manifold_Pressure ( float LeverPosn, float MaxMan)
  83 {
  84     float Inches;
  85     // if ( x < = 0 ) {
  86     //   x = 0.00001;
  87     // }
  88     Inches = LeverPosn * MaxMan / 100;
  89     return Inches;
  90 }
  91
  92
  93 // set initial default values
  94 void FGEngine::init() {
  95     // Control and environment inputs
  96     IAS = 0;
  97     Throttle_Lever_Pos = 75;
  98     Propeller_Lever_Pos = 75;
  99     Mixture_Lever_Pos = 100;
 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;
 105     Max_RPM = 2700;
 106     Min_RPM = 1000;
 107     Max_Fuel_Flow = 130;
 108     Mag_Derate_Percent = 5;
 109     MaxHP = 285;
 110     Gear_Ratio = 1;
 111
 112     // Initialise Engine Variables used by this instance
 113     Percentage_Power = 0;
 114     Manifold_Pressure = 29.00; // Inches
 115     RPM = 500;
 116     Fuel_Flow = 0;      // lbs/hour
 117     Torque = 0;
 118     CHT = 370;
 119     Mixture = 14;
 120     Oil_Pressure = 0;   // PSI
 121     Oil_Temp = 85;      // Deg C
 122     HP = 0;
 123     RPS = 0;
 124     Torque_Imbalance = 0;
 125     Desired_RPM = 0;
 126
 127     // Initialise Propellor Variables used by this instance
 128     FGProp1_Angular_V = 0;
 129     FGProp1_Coef_Drag =  0.6;
 130     FGProp1_Torque = 0;
 131     FGProp1_Thrust = 0;
 132     FGProp1_RPS = 0;
 133     FGProp1_Coef_Lift = 0.1;
 134     Alpha1 = 13.5;
 135     FGProp1_Blade_Angle = 13.5;
 136     FGProp_Fine_Pitch_Stop = 13.5;
 137     FGProp_Course_Pitch_Stop = 55;
 138
 139     // Other internal values
 140     Rho = 0.002378;
 141 }
 142
 143
 144 // Calculate Oil Pressure
 145 static float Oil_Press (float Oil_Temp, float Engine_RPM)
 146 {
 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
 153
 154     Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
 155
 156     // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
 157     if (Oil_Pressure >= Oil_Press_Relief_Valve)
 158         {
 159             Oil_Pressure = Oil_Press_Relief_Valve;
 160         }
 161
 162     // Now adjust pressure according to Temp which affects the viscosity
 163
 164     Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
 165
 166     return Oil_Pressure;
 167 }
 168
 169
 170 // Calculate Cylinder Head Temperature
 171 static float Calc_CHT (float Fuel_Flow, float Mixture, float IAS)
 172 {
 173     float CHT = 350;
 174
 175     return CHT;
 176 }
 177
 178
 179 // Calculate Density Ratio
 180 static float Density_Ratio ( float x )
 181 {
 182     float y = ((3E-10 * x * x) - (3E-05 * x) + 0.9998);
 183     return y;
 184 }
 185
 186
 187 // Calculate Air Density - Rho
 188 static float Density ( float x )
 189 {
 190     float y = ((9E-08 * x * x) - (7E-08 * x) + 0.0024);
 191     return y;
 192 }
 193
 194
 195 // Calculate Speed in FPS given Knots CAS
 196 static float IAS_to_FPS (float ias)
 197 {
 198     return ias * 1.68888888;
 199 }
 200
 201
 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.
 207     float ManXRPM = 0;
 208     float Vo = 0;
 209     float V1 = 0;
 210
 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;
 216
 217     // Input Variables
 218     // float IAS = 0;
 219
 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;
 226
 227     // Environmental Variables
 228
 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;
 233
 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;
 243
 244     // There needs to be a section in here to trap silly values, like
 245     // 0, otherwise they will crash the calculations
 246
 247     // cout << " Number of Iterations ";
 248     // cin >> num2;
 249     // cout << endl;
 250
 251     // cout << " Throttle % ";
 252     // cin >> Throttle_Lever_Pos;
 253     // cout << endl;
 254
 255     // cout << " Prop % ";
 256     // cin >> Propeller_Lever_Pos;
 257     // cout << endl;
 258
 259     //==================================================================
 260     // Engine & Environmental Inputs from elsewhere
 261
 262     // Calculate Air Density (Rho) - In FG this is calculated in
 263     // FG_Atomoshere.cxx
 264
 265     Rho = Density(FG_Pressure_Ht); // In FG FG_Pressure_Ht is "h"
 266     // cout << "Rho = " << Rho << endl;
 267
 268     // Calculate Manifold Pressure (Engine 1) as set by throttle opening
 269
 270     Manifold_Pressure =
 271         Calc_Manifold_Pressure( Throttle_Lever_Pos, Max_Manifold_Pressure );
 272     cout << "manifold pressure = " << Manifold_Pressure << endl;
 273
 274
 275     RPM = Calc_Engine_RPM(Propeller_Lever_Pos);
 276     // cout << "Engine RPM = " << RPM << endl;
 277
 278     Desired_RPM = RPM;
 279     cout << "Desired RPM = " << Desired_RPM << endl;
 280
 281     //==================================================================
 282     // Engine Power & Torque Calculations
 283
 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";
 289
 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;
 294
 295         //  Calculate % Power
 296         Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM)
 297             + ( + 7E-04 * ManXRPM) - 0.1218;
 298         // cout << "percent power = " << Percentage_Power <<  "%" << "\t";
 299
 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";
 305
 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";
 310
 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;
 314
 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   ";
 323
 324             Percentage_Power = Percentage_Power *
 325                 ((100.0 - Mag_Derate_Percent)/100.0);
 326         }
 327         // cout << "Final engine % power = " <<  Percentage_Power << "%" << endl;
 328
 329         // Calculate Engine Horsepower
 330
 331         HP = Percentage_Power * MaxHP / 100.0;
 332
 333         // Calculate Engine Torque
 334
 335         Torque = HP * 5252 / RPM;
 336         // cout << Torque << "Ft/lbs" << "\t";
 337
 338         // Calculate Cylinder Head Temperature
 339         CHT = Calc_CHT( Fuel_Flow, Mixture, IAS);
 340         // cout << "Cylinder Head Temp (F) = " << CHT << endl;
 341
 342         // Calculate Oil Pressure
 343         Oil_Pressure = Oil_Press( Oil_Temp, RPM );
 344         // cout << "Oil Pressure (PSI) = " << Oil_Pressure << endl;
 345
 346         //==============================================================
 347
 348         // Now do the Propellor Calculations
 349
 350         // Revs per second
 351         FGProp1_RPS = RPM * Gear_Ratio / 60.0;
 352         // cout << FGProp1_RPS << " RPS" <<  endl;
 353
 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;
 357
 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;
 363
 364         // cout << Vo << "Axial Velocity" << endl;
 365
 366         // Relative Velocity (V1)
 367         V1 = sqrt((FGProp1_Angular_V * FGProp1_Angular_V) +
 368                   (Vo * Vo));
 369         // cout << "Relative Velocity " << V1 << endl;
 370
 371         if ( FGProp1_Blade_Angle >= FGProp_Course_Pitch_Stop )  {
 372             FGProp1_Blade_Angle = FGProp_Course_Pitch_Stop;
 373         }
 374
 375         // cout << FGProp1_Blade_Angle << " Prop Blade Angle" << endl;
 376
 377         // Blade Angle of Attack (Alpha1)
 378
 379         Alpha1 = FGProp1_Blade_Angle -(atan(Vo / FGProp1_Angular_V) * (180/PI));
 380         // cout << Alpha1 << " Alpha1" << endl;
 381
 382         // cout << "  Alpha1 = " << Alpha1
 383         //      << "  Blade angle = " << FGProp1_Blade_Angle
 384         //      << "  Vo = " << Vo
 385         //      << "  FGProp1_Angular_V = " << FGProp1_Angular_V << endl;
 386
 387         // Calculate Coefficient of Drag at Alpha1
 388         FGProp1_Coef_Drag = (0.0005 * (Alpha1 * Alpha1)) + (0.0003 * Alpha1)
 389             + 0.0094;
 390         //      cout << FGProp1_Coef_Drag << " Coef Drag" << endl;
 391
 392         // Calculate Coefficient of Lift at Alpha1
 393         FGProp1_Coef_Lift = -(0.0026 * (Alpha1 * Alpha1)) + (0.1027 * Alpha1)
 394             + 0.2295;
 395         // cout << FGProp1_Coef_Lift << " Coef Lift " << endl;
 396
 397         // Covert Alplha1 to Radians
 398         // Alpha1 = Alpha1 * PI / 180;
 399
 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;
 406
 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;
 413
 414         // End of Propeller Calculations
 415         //==============================================================
 416
 417
 418
 419         Torque_Imbalance = FGProp1_Torque - Torque;
 420         //  cout <<  Torque_Imbalance << endl;
 421
 422         if (Torque_Imbalance > 20) {
 423             RPM -= 14.5;
 424             // FGProp1_RPM -= 25;
 425             FGProp1_Blade_Angle -= 0.75;
 426         }
 427
 428         if (FGProp1_Blade_Angle < FGProp_Fine_Pitch_Stop) {
 429             FGProp1_Blade_Angle = FGProp_Fine_Pitch_Stop;
 430         }
 431         if (Torque_Imbalance < -20) {
 432             RPM += 14.5;
 433             // FGProp1_RPM += 25;
 434             FGProp1_Blade_Angle += 0.75;
 435         }
 436
 437         if (RPM >= 2700) {
 438             RPM = 2700;
 439         }
 440
 441
 442         // cout << FGEng1_RPM << " Blade_Angle  " << FGProp1_Blade_Angle << endl << endl;
 443
 444     }
 445 }
 446
 447
 448
 449
 450 // Functions
 451
 452 // Calculate Oil Temperature
 453
 454 static float Oil_Temp (float Fuel_Flow, float Mixture, float IAS)
 455 {
 456     float Oil_Temp = 85;
 457
 458     return (Oil_Temp);
 459 }