-//#include <ansi_c.h>\r
-//#include <math.h>\r
-//#include <stdio.h>\r
-//#include <stdlib.h>\r
-#include "uiuc_ice_rates.h"\r
-\r
-///////////////////////////////////////////////////////////////////////\r
-// Calculates shed rate depending on current aero loads, eta, temp, and freezing fraction\r
-// Code by Leia Blumenthal\r
-//\r
-// 13 Feb 02 - Created basic program with dummy variables and a constant shed rate (no dependency)\r
-//\r
-// Inputs:\r
-// aero_load - aerodynamic load\r
-// eta \r
-// T - Temperature in Farenheit\r
-// ff - freezing fraction\r
-//\r
-// Output:\r
-// rate - %eta shed/time\r
-//\r
-// Right now this is just a constant shed rate until we learn more...\r
-\r
-\r
-double shed(double aero_load, double eta, double T, double ff, double time_step)\r
-{\r
- double rate, eta_new;\r
-\r
- if (eta <= 0.0)\r
- rate = 0.0;\r
- else\r
- rate = 0.2;\r
-\r
- eta_new = eta-rate*eta*time_step;\r
- if (eta_new <= 0.0)\r
- eta_new = 0.0;\r
- \r
- return(eta_new);\r
-}\r
-\r
-\r
-///////////////////////////////////////////////////////////////////////////////////////////////////\r
-// Currently a simple linear approximation based on temperature and eta, but for next version, \r
-// should have so that it calculates sublimation rate depending on current temp,pressure, \r
-// dewpoint, radiation, and eta\r
-//\r
-// Code by Leia Blumenthal \r
-// 12 Feb 02 - Created basic program with linear rate for values when sublimation will occur\r
-// 16 May 02 - Modified so that outputs new eta as opposed to rate\r
-// Inputs:\r
-// T - temperature and must be input in Farenheit\r
-// P - pressure\r
-// Tdew - Dew point Temperature\r
-// rad - radiation\r
-// time_step- increment since last run\r
-//\r
-// Intermediate:\r
-// rate - sublimation rate (% eta change/time)\r
-//\r
-// Output:\r
-// eta_new- eta after sublimation has occurred\r
-//\r
-// This takes a simple approximation that the rate of sublimation will decrease\r
-// linearly with temperature increase.\r
-//\r
-// This code should be run every time step to every couple time steps \r
-//\r
-// If eta is less than zero, than there should be no sublimation\r
-\r
-double sublimation(double T, double eta, double time_step)\r
-{\r
- double rate, eta_new;\r
- \r
- if (eta <= 0.0) rate = 0;\r
- \r
- else{ \r
- // According to the Smithsonian Meteorological tables sublimation occurs\r
- // between -40 deg F < T < 32 deg F and between pressures of 0 atm < P < 0.00592 atm\r
- if (T < -40) rate = 0;\r
- else if (T >= -40 && T < 32)\r
- {\r
- // For a simple linear approximation, assume largest value is a rate of .2% per sec\r
- rate = 0.0028 * T + 0.0889;\r
- }\r
- else if (T >= 32) rate = 0;\r
- }\r
-\r
- eta_new = eta-rate*eta*time_step;\r
- if (eta_new <= 0.0)\r
- eta_new = 0.0;\r
- \r
- return(eta_new);\r
-} \r
+//#include <ansi_c.h>
+//#include <math.h>
+//#include <stdio.h>
+//#include <stdlib.h>
+#include "uiuc_ice_rates.h"
+
+///////////////////////////////////////////////////////////////////////
+// Calculates shed rate depending on current aero loads, eta, temp, and freezing fraction
+// Code by Leia Blumenthal
+//
+// 13 Feb 02 - Created basic program with dummy variables and a constant shed rate (no dependency)
+//
+// Inputs:
+// aero_load - aerodynamic load
+// eta
+// T - Temperature in Farenheit
+// ff - freezing fraction
+//
+// Output:
+// rate - %eta shed/time
+//
+// Right now this is just a constant shed rate until we learn more...
+
+
+double shed(double aero_load, double eta, double T, double ff, double time_step)
+{
+ double rate, eta_new;
+
+ if (eta <= 0.0)
+ rate = 0.0;
+ else
+ rate = 0.2;
+
+ eta_new = eta-rate*eta*time_step;
+ if (eta_new <= 0.0)
+ eta_new = 0.0;
+
+ return(eta_new);
+}
+
+
+///////////////////////////////////////////////////////////////////////////////////////////////////
+// Currently a simple linear approximation based on temperature and eta, but for next version,
+// should have so that it calculates sublimation rate depending on current temp,pressure,
+// dewpoint, radiation, and eta
+//
+// Code by Leia Blumenthal
+// 12 Feb 02 - Created basic program with linear rate for values when sublimation will occur
+// 16 May 02 - Modified so that outputs new eta as opposed to rate
+// Inputs:
+// T - temperature and must be input in Farenheit
+// P - pressure
+// Tdew - Dew point Temperature
+// rad - radiation
+// time_step- increment since last run
+//
+// Intermediate:
+// rate - sublimation rate (% eta change/time)
+//
+// Output:
+// eta_new- eta after sublimation has occurred
+//
+// This takes a simple approximation that the rate of sublimation will decrease
+// linearly with temperature increase.
+//
+// This code should be run every time step to every couple time steps
+//
+// If eta is less than zero, than there should be no sublimation
+
+double sublimation(double T, double eta, double time_step)
+{
+ double rate, eta_new;
+
+ if (eta <= 0.0) rate = 0;
+
+ else{
+ // According to the Smithsonian Meteorological tables sublimation occurs
+ // between -40 deg F < T < 32 deg F and between pressures of 0 atm < P < 0.00592 atm
+ if (T < -40) rate = 0;
+ else if (T >= -40 && T < 32)
+ {
+ // For a simple linear approximation, assume largest value is a rate of .2% per sec
+ rate = 0.0028 * T + 0.0889;
+ }
+ else if (T >= 32) rate = 0;
+ }
+
+ eta_new = eta-rate*eta*time_step;
+ if (eta_new <= 0.0)
+ eta_new = 0.0;
+
+ return(eta_new);
+}
-// SIS Twin Otter Iced aircraft Nonlinear model\r
-// Version 020409\r
-// read readme_020212.doc for information\r
-\r
-#include "uiuc_iced_nonlin.h"\r
-\r
-void Calc_Iced_Forces()\r
- {\r
- // alpha in deg\r
- double alpha;\r
- double de;\r
- double eta_ref_wing = 0.08; // eta of iced data used for curve fit\r
- double eta_ref_tail = 0.12;\r
- double eta_wing;\r
- //double delta_CL; // CL_clean - CL_iced;\r
- //double delta_CD; // CD_clean - CD_iced;\r
- //double delta_Cm; // CM_clean - CM_iced;\r
- double delta_Cm_a; // (Cm_clean - Cm_iced) as a function of AoA;\r
- double delta_Cm_de; // (Cm_clean - Cm_iced) as a function of de;\r
- double delta_Ch_a;\r
- double delta_Ch_e;\r
- double KCL;\r
- double KCD;\r
- double KCm_alpha;\r
- double KCm_de;\r
- double KCh;\r
- double CL_diff;\r
- \r
- \r
- \r
- alpha = Alpha*RAD_TO_DEG;\r
- de = elevator*RAD_TO_DEG;\r
- // lift fits\r
- if (alpha < 16)\r
- {\r
- delta_CL = (0.088449 + 0.004836*alpha - 0.0005459*alpha*alpha +\r
- 4.0859e-5*pow(alpha,3));\r
- }\r
- else\r
- {\r
- delta_CL = (-11.838 + 1.6861*alpha - 0.076707*alpha*alpha +\r
- 0.001142*pow(alpha,3));\r
- }\r
- KCL = -delta_CL/eta_ref_wing;\r
- eta_wing = 0.5*(eta_wing_left + eta_wing_right);\r
- delta_CL = eta_wing*KCL;\r
- \r
- \r
- // drag fit\r
- delta_CD = (-0.0089 + 0.001578*alpha - 0.00046253*pow(alpha,2) +\r
- -4.7511e-5*pow(alpha,3) + 2.3384e-6*pow(alpha,4));\r
- KCD = -delta_CD/eta_ref_wing;\r
- delta_CD = eta_wing*KCD;\r
- \r
- // pitching moment fit\r
- delta_Cm_a = (-0.01892 - 0.0056476*alpha + 1.0205e-5*pow(alpha,2)\r
- - 4.0692e-5*pow(alpha,3) + 1.7594e-6*pow(alpha,4));\r
- \r
- delta_Cm_de = (-0.014928 - 0.0037783*alpha + 0.00039086*pow(de,2)\r
- - 1.1304e-5*pow(de,3) - 1.8439e-6*pow(de,4));\r
- \r
- delta_Cm = delta_Cm_a + delta_Cm_de;\r
- KCm_alpha = delta_Cm_a/eta_ref_wing;\r
- KCm_de = delta_Cm_de/eta_ref_tail;\r
- delta_Cm = (0.75*eta_wing + 0.25*eta_tail)*KCm_alpha + (eta_tail)*KCm_de;\r
- \r
- // hinge moment\r
- if (alpha < 13)\r
- {\r
- delta_Ch_a = (-0.0012862 - 0.00022705*alpha + 1.5911e-5*pow(alpha,2)\r
- + 5.4536e-7*pow(alpha,3));\r
- }\r
- else\r
- {\r
- delta_Ch_a = 0;\r
- }\r
- delta_Ch_e = -0.0011851 - 0.00049924*de;\r
- delta_Ch = -(delta_Ch_a + delta_Ch_e);\r
- KCh = -delta_Ch/eta_ref_tail;\r
- delta_Ch = eta_tail*KCh;\r
- \r
- // rolling moment\r
- CL_diff = (eta_wing_left - eta_wing_right)*KCL;\r
- delta_Cl = CL_diff/4;\r
- \r
- }\r
-\r
-void add_ice_effects()\r
-{\r
- CL_clean = -1*CZ*cos(Alpha) + CX*sin(Alpha); //Check later\r
- CD_clean = -1*CZ*sin(Alpha) - CX*cos(Alpha);\r
- Cm_clean = Cm;\r
- Cl_clean = Cl;\r
- Ch_clean = Ch;\r
-\r
- CL_iced = CL_clean + delta_CL;\r
- CD_iced = CD_clean + delta_CD;\r
- Cm_iced = Cm_clean + delta_Cm;\r
- Cl_iced = Cl_clean + delta_Cl;\r
- //Ch_iced = Ch_clean + delta_Ch;\r
-\r
- CL = CL_iced;\r
- CD = CD_iced;\r
- Cm = Cm_iced;\r
- Cl = Cl_iced;\r
- //Ch = Ch_iced;\r
-\r
- CZ = -1*CL*cos(Alpha) - CD*sin(Alpha);\r
- CX = CL*sin(Alpha) - CD*cos(Alpha);\r
-\r
-}\r
+// SIS Twin Otter Iced aircraft Nonlinear model
+// Version 020409
+// read readme_020212.doc for information
+
+#include "uiuc_iced_nonlin.h"
+
+void Calc_Iced_Forces()
+ {
+ // alpha in deg
+ double alpha;
+ double de;
+ double eta_ref_wing = 0.08; // eta of iced data used for curve fit
+ double eta_ref_tail = 0.12;
+ double eta_wing;
+ //double delta_CL; // CL_clean - CL_iced;
+ //double delta_CD; // CD_clean - CD_iced;
+ //double delta_Cm; // CM_clean - CM_iced;
+ double delta_Cm_a; // (Cm_clean - Cm_iced) as a function of AoA;
+ double delta_Cm_de; // (Cm_clean - Cm_iced) as a function of de;
+ double delta_Ch_a;
+ double delta_Ch_e;
+ double KCL;
+ double KCD;
+ double KCm_alpha;
+ double KCm_de;
+ double KCh;
+ double CL_diff;
+
+
+
+ alpha = Alpha*RAD_TO_DEG;
+ de = elevator*RAD_TO_DEG;
+ // lift fits
+ if (alpha < 16)
+ {
+ delta_CL = (0.088449 + 0.004836*alpha - 0.0005459*alpha*alpha +
+ 4.0859e-5*pow(alpha,3));
+ }
+ else
+ {
+ delta_CL = (-11.838 + 1.6861*alpha - 0.076707*alpha*alpha +
+ 0.001142*pow(alpha,3));
+ }
+ KCL = -delta_CL/eta_ref_wing;
+ eta_wing = 0.5*(eta_wing_left + eta_wing_right);
+ delta_CL = eta_wing*KCL;
+
+
+ // drag fit
+ delta_CD = (-0.0089 + 0.001578*alpha - 0.00046253*pow(alpha,2) +
+ -4.7511e-5*pow(alpha,3) + 2.3384e-6*pow(alpha,4));
+ KCD = -delta_CD/eta_ref_wing;
+ delta_CD = eta_wing*KCD;
+
+ // pitching moment fit
+ delta_Cm_a = (-0.01892 - 0.0056476*alpha + 1.0205e-5*pow(alpha,2)
+ - 4.0692e-5*pow(alpha,3) + 1.7594e-6*pow(alpha,4));
+
+ delta_Cm_de = (-0.014928 - 0.0037783*alpha + 0.00039086*pow(de,2)
+ - 1.1304e-5*pow(de,3) - 1.8439e-6*pow(de,4));
+
+ delta_Cm = delta_Cm_a + delta_Cm_de;
+ KCm_alpha = delta_Cm_a/eta_ref_wing;
+ KCm_de = delta_Cm_de/eta_ref_tail;
+ delta_Cm = (0.75*eta_wing + 0.25*eta_tail)*KCm_alpha + (eta_tail)*KCm_de;
+
+ // hinge moment
+ if (alpha < 13)
+ {
+ delta_Ch_a = (-0.0012862 - 0.00022705*alpha + 1.5911e-5*pow(alpha,2)
+ + 5.4536e-7*pow(alpha,3));
+ }
+ else
+ {
+ delta_Ch_a = 0;
+ }
+ delta_Ch_e = -0.0011851 - 0.00049924*de;
+ delta_Ch = -(delta_Ch_a + delta_Ch_e);
+ KCh = -delta_Ch/eta_ref_tail;
+ delta_Ch = eta_tail*KCh;
+
+ // rolling moment
+ CL_diff = (eta_wing_left - eta_wing_right)*KCL;
+ delta_Cl = CL_diff/4;
+
+ }
+
+void add_ice_effects()
+{
+ CL_clean = -1*CZ*cos(Alpha) + CX*sin(Alpha); //Check later
+ CD_clean = -1*CZ*sin(Alpha) - CX*cos(Alpha);
+ Cm_clean = Cm;
+ Cl_clean = Cl;
+ Ch_clean = Ch;
+
+ CL_iced = CL_clean + delta_CL;
+ CD_iced = CD_clean + delta_CD;
+ Cm_iced = Cm_clean + delta_Cm;
+ Cl_iced = Cl_clean + delta_Cl;
+ //Ch_iced = Ch_clean + delta_Ch;
+
+ CL = CL_iced;
+ CD = CD_iced;
+ Cm = Cm_iced;
+ Cl = Cl_iced;
+ //Ch = Ch_iced;
+
+ CZ = -1*CL*cos(Alpha) - CD*sin(Alpha);
+ CX = CL*sin(Alpha) - CD*cos(Alpha);
+
+}
projects / flightgear.git / commitdiff