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 <math.h>
  92 #include <stdio.h>
  93
  94
  95 #define NCL 11
  96 #define DYN_ON_SPEED 33 /*20 knots*/
  97
  98
  99 #ifdef USENZ
 100         #define NZ generic_.n_cg_body_v[2]
 101 #else
 102         #define NZ 1
 103 #endif
 104
 105
 106 extern COCKPIT cockpit_;
 107 FILE *out;
 108
 109 SCALAR interp(SCALAR *y_table, SCALAR *x_table, int Ntable, SCALAR x)
 110 {
 111         SCALAR slope;
 112         int i=1;
 113         float y;
 114
 115
 116         /* if x is outside the table, return value at x[0] or x[Ntable-1]*/
 117         if(x <= x_table[0])
 118         {
 119                  y=y_table[0];
 120                  /* printf("x smaller than x_table[0]: %g %g\n",x,x_table[0]); */
 121         }
 122         else if(x >= x_table[Ntable-1])
 123         {
 124                  y=y_table[Ntable-1];
 125                  /* printf("x larger than x_table[N]: %g %g %d\n",x,x_table[NCL-1],Ntable-1); */
 126         }
 127         else /*x is within the table, interpolate linearly to find y value*/
 128         {
 129
 130             while(x_table[i] <= x) {i++;}
 131             slope=(y_table[i]-y_table[i-1])/(x_table[i]-x_table[i-1]);
 132                 /* 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); */
 133             y=slope*(x-x_table[i-1]) +y_table[i-1];
 134         }
 135         return y;
 136 }
 137
 138 void record()
 139 {
 140
 141         fprintf(out,"%g,%g,%g,%g,%g,%g,%g,%g,%g,",Long_control,Lat_control,Rudder_pedal,Aft_trim,Fwd_trim,V_rel_wind,Dynamic_pressure,P_body,R_body);
 142         fprintf(out,"%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,",Alpha,Cos_alpha,Sin_alpha,Alpha_dot,Q_body,Theta_dot,Sin_theta,Cos_theta,Beta,Cos_beta,Sin_beta);
 143         fprintf(out,"%g,%g,%g,%g,%g,%g,%g,%g\n",Sin_phi,Cos_phi,F_X_aero,F_Y_aero,F_Z_aero,M_l_aero,M_m_aero,M_n_aero);
 144     fflush(out);
 145 }
 146
 147 void aero( SCALAR dt, int Initialize ) {
 148   static int init = 0;
 149
 150
 151   static SCALAR trim_inc = 0.0002;
 152   SCALAR long_trim;
 153
 154
 155   SCALAR elevator, aileron, rudder;
 156
 157   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};
 158   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};
 159
 160
 161
 162
 163   /*Note that CLo,Cdo,Cmo will likely change with flap setting so
 164   they may not be declared static in the future */
 165
 166
 167   static SCALAR CLadot=1.7;
 168   static SCALAR CLq=3.9;
 169   static SCALAR CLde=0.43;
 170   static SCALAR CLo=0;
 171
 172
 173   static SCALAR Cdo=0.031;
 174   static SCALAR Cda=0.13;  /*Not used*/
 175   static SCALAR Cdde=0.06;
 176
 177   static SCALAR Cma=-0.89;
 178   static SCALAR Cmadot=-5.2;
 179   static SCALAR Cmq=-12.4;
 180   static SCALAR Cmo=-0.062;
 181   static SCALAR Cmde=-1.28;
 182
 183   static SCALAR Clbeta=-0.089;
 184   static SCALAR Clp=-0.47;
 185   static SCALAR Clr=0.096;
 186   static SCALAR Clda=0.178;
 187   static SCALAR Cldr=0.0147;
 188
 189   static SCALAR Cnbeta=0.065;
 190   static SCALAR Cnp=-0.03;
 191   static SCALAR Cnr=-0.099;
 192   static SCALAR Cnda=-0.053;
 193   static SCALAR Cndr=-0.0657;
 194
 195   static SCALAR Cybeta=-0.31;
 196   static SCALAR Cyp=-0.037;
 197   static SCALAR Cyr=0.21;
 198   static SCALAR Cyda=0.0;
 199   static SCALAR Cydr=0.187;
 200
 201   /*nondimensionalization quantities*/
 202   /*units here are ft and lbs */
 203   static SCALAR cbar=4.9; /*mean aero chord ft*/
 204   static SCALAR b=35.8; /*wing span ft */
 205   static SCALAR Sw=174; /*wing planform surface area ft^2*/
 206   static SCALAR rPiARe=0.054; /*reciprocal of Pi*AR*e*/
 207
 208   SCALAR W=Mass/INVG;
 209
 210   SCALAR CLwbh,CL,cm,cd,cn,cy,croll,cbar_2V,b_2V,qS,qScbar,qSb,ps,rs;
 211
 212   SCALAR F_X_wind,F_Y_wind,F_Z_wind,W_X,W_Y,W_Z;
 213
 214
 215
 216
 217
 218   if (Initialize != 0)
 219     {
 220
 221
 222       out=fopen("flight.csv","w");
 223           /* Initialize aero coefficients */
 224
 225
 226     }
 227
 228   record();
 229
 230   /*
 231   LaRCsim uses:
 232     Cm > 0 => ANU
 233         Cl > 0 => Right wing down
 234         Cn > 0 => ANL
 235   so:
 236     elevator > 0 => AND -- aircraft nose down
 237         aileron > 0  => right wing up
 238         rudder > 0   => ANL
 239   */
 240
 241   if(Aft_trim) long_trim = long_trim - trim_inc;
 242   if(Fwd_trim) long_trim = long_trim + trim_inc;
 243
 244   /*scale pct control to degrees deflection*/
 245   if ((Long_control+long_trim) <= 0)
 246         elevator=(Long_control+long_trim)*-28*DEG_TO_RAD;
 247   else
 248         elevator=(Long_control+long_trim)*23*DEG_TO_RAD;
 249
 250   aileron  = Lat_control*17.5*DEG_TO_RAD;
 251   rudder   = Rudder_pedal*16*DEG_TO_RAD;
 252
 253
 254
 255
 256
 257   /*check control surface travel limits*/
 258   /* if((elevator+long_trim) > 23)
 259      elevator=23;
 260   else if((elevator+long_trim) < -28)
 261          elevator=-23; */
 262
 263
 264   /*
 265     The aileron travel limits are 20 deg. TEU and 15 deg TED
 266     but since we don't distinguish between left and right we'll
 267         use the average here (17.5 deg)
 268   */
 269   /* if(fabs(aileron) > 17.5)
 270          aileron = 17.5;
 271   if(fabs(rudder) > 16)
 272          rudder = 16; */
 273
 274   /*calculate rate derivative nondimensionalization (is that a word?) factors */
 275   /*hack to avoid divide by zero*/
 276   /*the dynamic terms might be negligible at low ground speeds anyway*/
 277
 278   if(V_rel_wind > DYN_ON_SPEED)
 279   {
 280         cbar_2V=cbar/(2*V_rel_wind);
 281         b_2V=b/(2*V_rel_wind);
 282   }
 283   else
 284   {
 285         cbar_2V=0;
 286         b_2V=0;
 287   }
 288
 289   /*calcuate the qS nondimensionalization factors*/
 290
 291   qS=Dynamic_pressure*Sw;
 292   qScbar=qS*cbar;
 293   qSb=qS*b;
 294
 295   /*transform the aircraft rotation rates*/
 296   ps=-P_body*Cos_alpha + R_body*Sin_alpha;
 297   rs=-P_body*Sin_alpha + R_body*Cos_alpha;
 298
 299
 300   /* sum coefficients */
 301   CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
 302   CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
 303   cd = Cdo + rPiARe*CL*CL + Cdde*elevator;
 304   cy = Cybeta*Beta + (Cyp*ps + Cyr*rs)*b_2V + Cyda*aileron + Cydr*rudder;
 305
 306   cm = Cmo + Cma*Alpha + (Cmq*Theta_dot + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator+long_trim);
 307   cn = Cnbeta*Beta + (Cnp*ps + Cnr*rs)*b_2V + Cnda*aileron + Cndr*rudder;
 308   croll=Clbeta*Beta + (Clp*ps + Clr*rs)*b_2V + Clda*aileron + Cldr*rudder;
 309
 310   /*calculate wind axes forces*/
 311   F_X_wind=-1*cd*qS;
 312   F_Y_wind=cy*qS;
 313   F_Z_wind=-1*CL*qS;
 314
 315   /*calculate moments and body axis forces */
 316
 317   /*find body-axis components of weight*/
 318   /*with earth axis to body axis transform */
 319   W_X=-1*W*Sin_theta;
 320   W_Y=W*Sin_phi*Cos_theta;
 321   W_Z=W*Cos_phi*Cos_theta;
 322
 323   /* requires ugly wind-axes to body-axes transform */
 324   F_X_aero = W_X + F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
 325   F_Y_aero = W_Y + F_X_wind*Sin_beta + F_Z_wind*Cos_beta;
 326   F_Z_aero = W_Z*NZ + F_X_wind*Sin_alpha*Cos_beta - F_Y_wind*Sin_alpha*Sin_beta + F_Z_wind*Cos_alpha;
 327
 328   /*no axes transform here */
 329   M_l_aero = I_xx*croll*qSb;
 330   M_m_aero = I_yy*cm*qScbar;
 331   M_n_aero = I_zz*cn*qSb;
 332
 333 }
 334
 335