src/FDM/LaRCsim/c172_aero.c

   1 /***************************************************************************
   2
   3   TITLE:        c172_aero
   4
   5 ----------------------------------------------------------------------------
   6
   7   FUNCTION:     aerodynamics model based on constant stability derivatives
   8
   9 ----------------------------------------------------------------------------
  10
  11   MODULE STATUS:        developmental
  12
  13 ----------------------------------------------------------------------------
  14
  15   GENEALOGY:    Based on data from:
  16                                 Part 1 of Roskam's S&C text
  17                                 The FAA type certificate data sheet for the 172
  18                                 Various sources on the net
  19                                 John D. Anderson's Intro to Flight text (NACA 2412 data)
  20                                 UIUC's airfoil data web site
  21
  22 ----------------------------------------------------------------------------
  23
  24   DESIGNED BY:  Tony Peden
  25
  26   CODED BY:             Tony Peden
  27
  28   MAINTAINED BY:        Tony Peden
  29
  30 ----------------------------------------------------------------------------
  31
  32   MODIFICATION HISTORY:
  33
  34   DATE          PURPOSE                                                                                         BY
  35   6/10/99   Initial test release
  36
  37
  38 ----------------------------------------------------------------------------
  39
  40   REFERENCES:
  41
  42   Aero Coeffs:
  43         CL                      lift
  44         Cd                      drag
  45         Cm                      pitching moment
  46         Cy                      sideforce
  47         Cn                      yawing moment
  48         Croll,Cl        rolling moment (yeah, I know.  Shoot me.)
  49
  50   Subscripts
  51         o               constant i.e. not a function of alpha or beta
  52         a               alpha
  53         adot    d(alpha)/dt
  54         q       pitch rate
  55         qdot    d(q)/dt
  56         beta    sideslip angle
  57         p               roll rate
  58         r               yaw rate
  59         da      aileron deflection
  60         de      elevator deflection
  61         dr      rudder deflection
  62
  63         s               stability axes
  64
  65
  66
  67 ----------------------------------------------------------------------------
  68
  69   CALLED BY:
  70
  71 ----------------------------------------------------------------------------
  72
  73   CALLS TO:
  74
  75 ----------------------------------------------------------------------------
  76
  77   INPUTS:
  78
  79 ----------------------------------------------------------------------------
  80
  81   OUTPUTS:
  82
  83 --------------------------------------------------------------------------*/
  84
  85
  86
  87 #include "ls_generic.h"
  88 #include "ls_cockpit.h"
  89 #include "ls_constants.h"
  90 #include "ls_types.h"
  91 #include "c172_aero.h"
  92
  93 #include <math.h>
  94 #include <stdio.h>
  95
  96
  97 #define NCL 11
  98 #define DYN_ON_SPEED 33 /*20 knots*/
  99
 100
 101 #ifdef USENZ
 102         #define NZ generic_.n_cg_body_v[2]
 103 #else
 104         #define NZ 1
 105 #endif
 106
 107
 108 extern COCKPIT cockpit_;
 109
 110
 111 SCALAR interp(SCALAR *y_table, SCALAR *x_table, int Ntable, SCALAR x)
 112 {
 113         SCALAR slope;
 114         int i=1;
 115         float y;
 116
 117
 118         /* if x is outside the table, return value at x[0] or x[Ntable-1]*/
 119         if(x <= x_table[0])
 120         {
 121                  y=y_table[0];
 122                  /* printf("x smaller than x_table[0]: %g %g\n",x,x_table[0]); */
 123         }
 124         else if(x >= x_table[Ntable-1])
 125         {
 126                  y=y_table[Ntable-1];
 127                  /* printf("x larger than x_table[N]: %g %g %d\n",x,x_table[NCL-1],Ntable-1); */
 128         }
 129         else /*x is within the table, interpolate linearly to find y value*/
 130         {
 131
 132             while(x_table[i] <= x) {i++;}
 133             slope=(y_table[i]-y_table[i-1])/(x_table[i]-x_table[i-1]);
 134                 /* printf("x: %g, i: %d, cl[i]: %g, cl[i-1]: %g, slope: %g\n",x,i,y_table[i],y_table[i-1],slope); */
 135             y=slope*(x-x_table[i-1]) +y_table[i-1];
 136         }
 137         return y;
 138 }
 139
 140
 141
 142 void aero( SCALAR dt, int Initialize ) {
 143
 144
 145   static int init = 0;
 146
 147
 148   static SCALAR trim_inc = 0.0002;
 149
 150   static SCALAR alpha_ind[NCL]={-0.087,0,0.175,0.209,0.24,0.262,0.278,0.303,0.314,0.332,0.367};
 151   static SCALAR CLtable[NCL]={-0.14,0.31,1.21,1.376,1.51249,1.591,1.63,1.60878,1.53712,1.376,1.142};
 152
 153  /*Note that CLo,Cdo,Cmo will likely change with flap setting so
 154   they may not be declared static in the future */
 155    if (Initialize != 0)
 156     {
 157            CLadot=1.7;
 158            CLq=3.9;
 159            CLde=0.43;
 160            CLo=0;
 161
 162
 163            Cdo=0.031;
 164            Cda=0.13;  /*Not used*/
 165            Cdde=0.06;
 166
 167            Cma=-0.89;
 168            Cmadot=-5.2;
 169            Cmq=-12.4;
 170            Cmo=-0.062;
 171            Cmde=-1.28;
 172
 173            Clbeta=-0.089;
 174            Clp=-0.47;
 175            Clr=0.096;
 176            Clda=-0.178;
 177            Cldr=0.0147;
 178
 179            Cnbeta=0.065;
 180            Cnp=-0.03;
 181            Cnr=-0.099;
 182            Cnda=-0.053;
 183            Cndr=-0.0657;
 184
 185            Cybeta=-0.31;
 186            Cyp=-0.037;
 187            Cyr=0.21;
 188            Cyda=0.0;
 189            Cydr=0.187;
 190
 191            /*nondimensionalization quantities*/
 192            /*units here are ft and lbs */
 193            cbar=4.9; /*mean aero chord ft*/
 194            b=35.8; /*wing span ft */
 195            Sw=174; /*wing planform surface area ft^2*/
 196            rPiARe=0.054; /*reciprocal of Pi*AR*e*/
 197     }
 198
 199   /*
 200   LaRCsim uses:
 201     Cm > 0 => ANU
 202         Cl > 0 => Right wing down
 203         Cn > 0 => ANL
 204   so:
 205     elevator > 0 => AND -- aircraft nose down
 206         aileron > 0  => right wing up
 207         rudder > 0   => ANL
 208   */
 209   long_trim=0;
 210   if(Aft_trim) long_trim = long_trim - trim_inc;
 211   if(Fwd_trim) long_trim = long_trim + trim_inc;
 212
 213 /*   printf("Long_control: %7.4f, long_trim: %7.4f,DEG_TO_RAD: %7.4f, RAD_TO_DEG: %7.4f\n",Long_control,long_trim,DEG_TO_RAD,RAD_TO_DEG);
 214  */  /*scale pct control to degrees deflection*/
 215   if ((Long_control+long_trim) <= 0)
 216         elevator=(Long_control+long_trim)*28*DEG_TO_RAD;
 217   else
 218         elevator=(Long_control+long_trim)*23*DEG_TO_RAD;
 219
 220   aileron  = Lat_control*17.5*DEG_TO_RAD;
 221   rudder   = Rudder_pedal*16*DEG_TO_RAD;
 222   /*
 223     The aileron travel limits are 20 deg. TEU and 15 deg TED
 224     but since we don't distinguish between left and right we'll
 225         use the average here (17.5 deg)
 226   */
 227
 228
 229   /*calculate rate derivative nondimensionalization (is that a word?) factors */
 230   /*hack to avoid divide by zero*/
 231   /*the dynamic terms might be negligible at low ground speeds anyway*/
 232
 233   if(V_rel_wind > DYN_ON_SPEED)
 234   {
 235         cbar_2V=cbar/(2*V_rel_wind);
 236         b_2V=b/(2*V_rel_wind);
 237   }
 238   else
 239   {
 240         cbar_2V=0;
 241         b_2V=0;
 242   }
 243
 244
 245   /*calcuate the qS nondimensionalization factors*/
 246
 247   qS=Dynamic_pressure*Sw;
 248   qScbar=qS*cbar;
 249   qSb=qS*b;
 250
 251
 252 /*   printf("Wb: %7.4f, Ub: %7.4f, Alpha: %7.4f, elev: %7.4f, ail: %7.4f, rud: %7.4f, long_trim: %7.4f\n",W_body,U_body,Alpha*RAD_TO_DEG,elevator*RAD_TO_DEG,aileron*RAD_TO_DEG,rudder*RAD_TO_DEG,long_trim*RAD_TO_DEG);
 253   printf("Theta: %7.4f, Gamma: %7.4f, Beta: %7.4f, Phi: %7.4f, Psi: %7.4f\n",Theta*RAD_TO_DEG,Gamma_vert_rad*RAD_TO_DEG,Beta*RAD_TO_DEG,Phi*RAD_TO_DEG,Psi*RAD_TO_DEG);
 254  */
 255
 256   /* sum coefficients */
 257   CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
 258   CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
 259   cd = Cdo + rPiARe*CL*CL + Cdde*elevator;
 260   cy = Cybeta*Beta + (Cyp*P_body + Cyr*R_body)*b_2V + Cyda*aileron + Cydr*rudder;
 261
 262   cm = Cmo + Cma*Alpha + (Cmq*Q_body + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator+long_trim);
 263   cn = Cnbeta*Beta + (Cnp*P_body + Cnr*R_body)*b_2V + Cnda*aileron + Cndr*rudder;
 264   croll=Clbeta*Beta + (Clp*P_body + Clr*R_body)*b_2V + Clda*aileron + Cldr*rudder;
 265
 266 /*   printf("CL: %7.4f, Cd: %7.4f, Cm: %7.4f, Cy: %7.4f, Cn: %7.4f, Cl: %7.4f\n",CL,cd,cm,cy,cn,croll);
 267  */  /*calculate wind axes forces*/
 268   F_X_wind=-1*cd*qS;
 269   F_Y_wind=cy*qS;
 270   F_Z_wind=-1*CL*qS;
 271
 272 /*   printf("V_rel_wind: %7.4f, Fxwind: %7.4f Fywind: %7.4f Fzwind: %7.4f\n",V_rel_wind,F_X_wind,F_Y_wind,F_Z_wind);
 273  */
 274   /*calculate moments and body axis forces */
 275
 276
 277
 278   /* requires ugly wind-axes to body-axes transform */
 279   F_X_aero = F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
 280   F_Y_aero = F_X_wind*Sin_beta + F_Y_wind*Cos_beta;
 281   F_Z_aero = F_X_wind*Sin_alpha*Cos_beta - F_Y_wind*Sin_alpha*Sin_beta + F_Z_wind*Cos_alpha;
 282
 283   /*no axes transform here */
 284   M_l_aero = croll*qSb;
 285   M_m_aero = cm*qScbar;
 286   M_n_aero = cn*qSb;
 287
 288 /*   printf("I_yy: %7.4f, qScbar: %7.4f, qbar: %7.4f, Sw: %7.4f, cbar: %7.4f, 0.5*rho*V^2: %7.4f\n",I_yy,qScbar,Dynamic_pressure,Sw,cbar,0.5*0.0023081*V_rel_wind*V_rel_wind);
 289  */
 290 /*  printf("Fxaero: %7.4f Fyaero: %7.4f Fzaero: %7.4f Weight: %7.4f\n",F_X_aero,F_Y_aero,F_Z_aero,W);
 291  *//*  printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);
 292  */
 293 }
 294