src/FDM/UIUCModel/uiuc_iced_nonlin.cpp

   1 //     SIS Twin Otter Iced aircraft Nonlinear model
   2 //     Version 020409
   3 //     read readme_020212.doc for information
   4
   5 #include "uiuc_iced_nonlin.h"
   6
   7 void Calc_Iced_Forces()
   8         {
   9         // alpha in deg
  10           double alpha;
  11           double de;
  12         double eta_ref_wing = 0.08;                      // eta of iced data used for curve fit
  13         double eta_ref_tail = 0.12;
  14         double eta_wing;
  15         //double delta_CL;                              // CL_clean - CL_iced;
  16         //double delta_CD;                              // CD_clean - CD_iced;
  17         //double delta_Cm;                              // CM_clean - CM_iced;
  18         double delta_Cm_a;                              // (Cm_clean - Cm_iced) as a function of AoA;
  19         double delta_Cm_de;                             // (Cm_clean - Cm_iced) as a function of de;
  20         double delta_Ch_a;
  21         double delta_Ch_e;
  22         double KCL;
  23         double KCD;
  24         double KCm_alpha;
  25         double KCm_de;
  26         double KCh;
  27         double CL_diff;
  28
  29
  30
  31         alpha = Alpha*RAD_TO_DEG;
  32         de = elevator*RAD_TO_DEG;
  33         // lift fits
  34         if (alpha < 16)
  35                 {
  36                 delta_CL = (0.088449 + 0.004836*alpha - 0.0005459*alpha*alpha +
  37                                         4.0859e-5*pow(alpha,3));
  38                 }
  39         else
  40                 {
  41                 delta_CL = (-11.838 + 1.6861*alpha - 0.076707*alpha*alpha +
  42                                         0.001142*pow(alpha,3));
  43                 }
  44         KCL = -delta_CL/eta_ref_wing;
  45         eta_wing = 0.5*(eta_wing_left + eta_wing_right);
  46         delta_CL = eta_wing*KCL;
  47
  48
  49         // drag fit
  50         delta_CD = (-0.0089 + 0.001578*alpha - 0.00046253*pow(alpha,2) +
  51                                         -4.7511e-5*pow(alpha,3) + 2.3384e-6*pow(alpha,4));
  52         KCD = -delta_CD/eta_ref_wing;
  53         delta_CD = eta_wing*KCD;
  54
  55         // pitching moment fit
  56         delta_Cm_a = (-0.01892 - 0.0056476*alpha + 1.0205e-5*pow(alpha,2)
  57                       - 4.0692e-5*pow(alpha,3) + 1.7594e-6*pow(alpha,4));
  58
  59         delta_Cm_de = (-0.014928 - 0.0037783*alpha + 0.00039086*pow(de,2)
  60                                         - 1.1304e-5*pow(de,3) - 1.8439e-6*pow(de,4));
  61
  62         delta_Cm = delta_Cm_a + delta_Cm_de;
  63         KCm_alpha = delta_Cm_a/eta_ref_wing;
  64         KCm_de = delta_Cm_de/eta_ref_tail;
  65         delta_Cm = (0.75*eta_wing + 0.25*eta_tail)*KCm_alpha + (eta_tail)*KCm_de;
  66
  67         // hinge moment
  68         if (alpha < 13)
  69                 {
  70                 delta_Ch_a = (-0.0012862 - 0.00022705*alpha + 1.5911e-5*pow(alpha,2)
  71                                                 + 5.4536e-7*pow(alpha,3));
  72                 }
  73         else
  74                 {
  75                 delta_Ch_a = 0;
  76                 }
  77         delta_Ch_e = -0.0011851 - 0.00049924*de;
  78         delta_Ch = -(delta_Ch_a + delta_Ch_e);
  79         KCh = -delta_Ch/eta_ref_tail;
  80         delta_Ch = eta_tail*KCh;
  81
  82         // rolling moment
  83         CL_diff = (eta_wing_left - eta_wing_right)*KCL;
  84         delta_Cl = CL_diff/4;
  85
  86         }
  87
  88 void add_ice_effects()
  89 {
  90   CL_clean = -1*CZ*cos(Alpha) + CX*sin(Alpha);  //Check later
  91   CD_clean = -1*CZ*sin(Alpha) - CX*cos(Alpha);
  92   Cm_clean = Cm;
  93   Cl_clean = Cl;
  94   Ch_clean = Ch;
  95
  96   CL_iced = CL_clean + delta_CL;
  97   CD_iced = CD_clean + delta_CD;
  98   Cm_iced = Cm_clean + delta_Cm;
  99   Cl_iced = Cl_clean + delta_Cl;
 100   //Ch_iced = Ch_clean + delta_Ch;
 101
 102   CL = CL_iced;
 103   CD = CD_iced;
 104   Cm = Cm_iced;
 105   Cl = Cl_iced;
 106   //Ch = Ch_iced;
 107
 108   CZ = -1*CL*cos(Alpha) - CD*sin(Alpha);
 109   CX = CL*sin(Alpha) - CD*cos(Alpha);
 110
 111 }