- // Adjust for Altitude. In this version a linear variation is
- // used. Decrease 1% for each 1000' increase in Altitde
- Percentage_Power = Percentage_Power + (FG_Pressure_Ht * 12/10000);
- // cout << Percentage_Power << "%" << "\t";
-
- //DCL - now adjust power to compensate for mixture
- //uses a curve fit to the data in the IO360 / O360 operating manual
- //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,
- //possibly by using separate fits for rich and lean of best power mixture
- //first adjust actual mixture to abstract mixture - this is a temporary hack in order to account for the fact that the data I have
- //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.
- //y=10x-12 for now
- abstract_mixture = 10.0 * equivalence_ratio - 12.0;
- float m = abstract_mixture; //to simplify writing the next equation
- 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);
- Percentage_Power = Percentage_Power * Percentage_of_best_power_mixture_power / 100.0;
-
-
- // Now Calculate Fuel Flow based on % Power Best Power Mixture
- Fuel_Flow = Percentage_Power * Max_Fuel_Flow / 100.0;
- // cout << Fuel_Flow << " lbs/hr"<< endl;
-
- // Now Derate engine for the effects of Bad/Switched off magnetos
- if (Magneto_Left == 0 && Magneto_Right == 0) {
- // cout << "Both OFF\n";
- Percentage_Power = 0;
- } else if (Magneto_Left && Magneto_Right) {
- // cout << "Both On ";
- } else if (Magneto_Left == 0 || Magneto_Right== 0) {
- // cout << "1 Magneto Failed ";
-
- Percentage_Power = Percentage_Power *
- ((100.0 - Mag_Derate_Percent)/100.0);
- // cout << FGEng1_Percentage_Power << "%" << "\t";
- }
-
- // Calculate Engine Horsepower
-
- HP = Percentage_Power * MaxHP / 100.0;
-
- Power_SI = HP * CONVERT_HP_TO_WATTS;
-
- // Calculate Engine Torque
-
- Torque = HP * 5252 / RPM;
- // cout << Torque << "Ft/lbs" << "\t";
-
- Torque_SI = (Power_SI * 60.0) / (2.0 * PI * RPM); //Torque = power / angular velocity
- // cout << Torque << " Nm\n";
-
- // Calculate Cylinder Head Temperature
- CHT = Calc_CHT( Fuel_Flow, Mixture, IAS);
- // cout << "Cylinder Head Temp (F) = " << CHT << endl;
-
-// EGT = Calc_EGT( Mixture );
-
- // Calculate Oil Pressure
- Oil_Pressure = Oil_Press( Oil_Temp, RPM );
- // cout << "Oil Pressure (PSI) = " << Oil_Pressure << endl;
-
- //==============================================================
-
- // Now do the Propellor Calculations
-
-#ifdef PHILS_PROP_MODEL
-
- // Revs per second
- FGProp1_RPS = RPM * Gear_Ratio / 60.0;
- // cout << FGProp1_RPS << " RPS" << endl;
-
- //Radial Flow Vector (V2) Ft/sec at Ref Blade Station (usually 30")
- FGProp1_Angular_V = FGProp1_RPS * 2 * PI * (Blade_Station / 12);
- // cout << FGProp1_Angular_V << "Angular Velocity " << endl;
-
- // Axial Flow Vector (Vo) Ft/sec
- // Some further work required here to allow for inflow at low speeds
- // Vo = (IAS + 20) * 1.688888;
- Vo = IAS_to_FPS(IAS + 20);
- // cout << "Feet/sec = " << Vo << endl;
-
- // cout << Vo << "Axial Velocity" << endl;
-
- // Relative Velocity (V1)
- V1 = sqrt((FGProp1_Angular_V * FGProp1_Angular_V) +
- (Vo * Vo));
- // cout << V1 << "Relative Velocity " << endl;
-
- // cout << FGProp1_Blade_Angle << " Prop Blade Angle" << endl;
-
- // Blade Angle of Attack (Alpha1)
-
-/* cout << " Alpha1 = " << Alpha1
- << " Blade angle = " << FGProp1_Blade_Angle
- << " Vo = " << Vo
- << " FGProp1_Angular_V = " << FGProp1_Angular_V << endl;*/
- Alpha1 = FGProp1_Blade_Angle -(atan(Vo / FGProp1_Angular_V) * (180/PI));
- // cout << Alpha1 << " Alpha1" << endl;
-
- // Calculate Coefficient of Drag at Alpha1
- FGProp1_Coef_Drag = (0.0005 * (Alpha1 * Alpha1)) + (0.0003 * Alpha1)
- + 0.0094;
- // cout << FGProp1_Coef_Drag << " Coef Drag" << endl;
-
- // Calculate Coefficient of Lift at Alpha1
- FGProp1_Coef_Lift = -(0.0026 * (Alpha1 * Alpha1)) + (0.1027 * Alpha1)
- + 0.2295;
- // cout << FGProp1_Coef_Lift << " Coef Lift " << endl;
-
- // Covert Alplha1 to Radians
- // Alpha1 = Alpha1 * PI / 180;
-
- // Calculate Prop Torque
- FGProp1_Torque = (0.5 * Rho * (V1 * V1) * FGProp_Area
- * ((FGProp1_Coef_Lift * sin(Alpha1 * PI / 180))
- + (FGProp1_Coef_Drag * cos(Alpha1 * PI / 180))))
- * (Blade_Station/12);
- // cout << FGProp1_Torque << " Prop Torque" << endl;
-
- // Calculate Prop Thrust
- // cout << " V1 = " << V1 << " Alpha1 = " << Alpha1 << endl;
- FGProp1_Thrust = 0.5 * Rho * (V1 * V1) * FGProp_Area
- * ((FGProp1_Coef_Lift * cos(Alpha1 * PI / 180))
- - (FGProp1_Coef_Drag * sin(Alpha1 * PI / 180)));
- // cout << FGProp1_Thrust << " Prop Thrust " << endl;
-
- // End of Propeller Calculations
- //==============================================================
-
-#endif //PHILS_PROP_MODEL
-
-#ifdef NEVS_PROP_MODEL
-
- // Nev's prop model
-
- num_elements = 6.0;
- number_of_blades = 2.0;
- blade_length = 0.95;
- allowance_for_spinner = blade_length / 12.0;
- prop_fudge_factor = 1.453401525;
- forward_velocity = IAS;
-
- theta[0] = 25.0;
- theta[1] = 20.0;
- theta[2] = 15.0;
- theta[3] = 10.0;
- theta[4] = 5.0;
- theta[5] = 0.0;
-
- angular_velocity_SI = 2.0 * PI * RPM / 60.0;
-
- allowance_for_spinner = blade_length / 12.0;
- //Calculate thrust and torque by summing the contributions from each of the blade elements
- //Assumes equal length elements with numbered 1 inboard -> num_elements outboard
- prop_torque = 0.0;
- prop_thrust = 0.0;
- int i;
-// outfile << "Rho = " << Rho << '\n\n';
-// outfile << "Drag = ";
- for(i=1;i<=num_elements;i++)
- {
- element = float(i);
- distance = (blade_length * (element / num_elements)) + allowance_for_spinner;
- 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))))))
- * (0.1 * (blade_length / element)) * number_of_blades;
-
- 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)
- * (0.1 * (blade_length / element)) * number_of_blades;
- element_torque = element_drag * distance;
- prop_torque += element_torque;
- prop_thrust += element_lift;
-// outfile << "Drag = " << element_drag << " n = " << element << '\n';
- }