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
 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
 139
 140 void aero( SCALAR dt, int Initialize ) {
 141
 142
 143   static int init = 0;
 144
 145
 146   static SCALAR trim_inc = 0.0002;
 147   SCALAR long_trim;
 148
 149
 150   SCALAR elevator, aileron, rudder;
 151
 152
 153
 154   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};
 155   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};
 156
 157
 158
 159
 160   /*Note that CLo,Cdo,Cmo will likely change with flap setting so
 161   they may not be declared static in the future */
 162
 163
 164   static SCALAR CLadot=1.7;
 165   static SCALAR CLq=3.9;
 166   static SCALAR CLde=0.43;
 167   static SCALAR CLo=0;
 168
 169
 170   static SCALAR Cdo=0.031;
 171   static SCALAR Cda=0.13;  /*Not used*/
 172   static SCALAR Cdde=0.06;
 173
 174   static SCALAR Cma=-0.89;
 175   static SCALAR Cmadot=-5.2;
 176   static SCALAR Cmq=-12.4;
 177   static SCALAR Cmo=-0.062;
 178   static SCALAR Cmde=-1.28;
 179
 180   static SCALAR Clbeta=-0.089;
 181   static SCALAR Clp=-0.47;
 182   static SCALAR Clr=0.096;
 183   static SCALAR Clda=-0.178;
 184   static SCALAR Cldr=0.0147;
 185
 186   static SCALAR Cnbeta=0.065;
 187   static SCALAR Cnp=-0.03;
 188   static SCALAR Cnr=-0.099;
 189   static SCALAR Cnda=-0.053;
 190   static SCALAR Cndr=-0.0657;
 191
 192   static SCALAR Cybeta=-0.31;
 193   static SCALAR Cyp=-0.037;
 194   static SCALAR Cyr=0.21;
 195   static SCALAR Cyda=0.0;
 196   static SCALAR Cydr=0.187;
 197
 198   /*nondimensionalization quantities*/
 199   /*units here are ft and lbs */
 200   static SCALAR cbar=4.9; /*mean aero chord ft*/
 201   static SCALAR b=35.8; /*wing span ft */
 202   static SCALAR Sw=174; /*wing planform surface area ft^2*/
 203   static SCALAR rPiARe=0.054; /*reciprocal of Pi*AR*e*/
 204
 205   SCALAR W=Mass/INVG;
 206
 207   SCALAR CLwbh,CL,cm,cd,cn,cy,croll,cbar_2V,b_2V,qS,qScbar,qSb,ps,rs;
 208
 209   SCALAR F_X_wind,F_Y_wind,F_Z_wind,W_X,W_Y,W_Z;
 210
 211
 212
 213
 214
 215   if (Initialize != 0)
 216     {
 217
 218
 219
 220
 221
 222     }
 223
 224
 225
 226   /*
 227   LaRCsim uses:
 228     Cm > 0 => ANU
 229         Cl > 0 => Right wing down
 230         Cn > 0 => ANL
 231   so:
 232     elevator > 0 => AND -- aircraft nose down
 233         aileron > 0  => right wing up
 234         rudder > 0   => ANL
 235   */
 236
 237   if(Aft_trim) long_trim = long_trim - trim_inc;
 238   if(Fwd_trim) long_trim = long_trim + trim_inc;
 239
 240   /*scale pct control to degrees deflection*/
 241   if ((Long_control+long_trim) <= 0)
 242         elevator=(Long_control+long_trim)*28*DEG_TO_RAD;
 243   else
 244         elevator=(Long_control+long_trim)*23*DEG_TO_RAD;
 245
 246   aileron  = Lat_control*17.5*DEG_TO_RAD;
 247   rudder   = Rudder_pedal*16*DEG_TO_RAD;
 248   /*
 249     The aileron travel limits are 20 deg. TEU and 15 deg TED
 250     but since we don't distinguish between left and right we'll
 251         use the average here (17.5 deg)
 252   */
 253
 254
 255   /*calculate rate derivative nondimensionalization (is that a word?) factors */
 256   /*hack to avoid divide by zero*/
 257   /*the dynamic terms might be negligible at low ground speeds anyway*/
 258
 259   if(V_rel_wind > DYN_ON_SPEED)
 260   {
 261         cbar_2V=cbar/(2*V_rel_wind);
 262         b_2V=b/(2*V_rel_wind);
 263   }
 264   else
 265   {
 266         cbar_2V=0;
 267         b_2V=0;
 268   }
 269
 270
 271   /*calcuate the qS nondimensionalization factors*/
 272
 273   qS=Dynamic_pressure*Sw;
 274   qScbar=qS*cbar;
 275   qSb=qS*b;
 276
 277   /*transform the aircraft rotation rates*/
 278   ps=-P_body*Cos_alpha + R_body*Sin_alpha;
 279   rs=-P_body*Sin_alpha + R_body*Cos_alpha;
 280
 281 /*   printf("Wb: %7.4f, Ub: %7.4f, Alpha: %7.4f, elev: %7.4f, ail: %7.4f, rud: %7.4f\n",W_body,U_body,Alpha*RAD_TO_DEG,elevator*RAD_TO_DEG,aileron*RAD_TO_DEG,rudder*RAD_TO_DEG);
 282   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);
 283  */
 284
 285   /* sum coefficients */
 286   CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
 287   CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
 288   cd = Cdo + rPiARe*CL*CL + Cdde*elevator;
 289   cy = Cybeta*Beta + (Cyp*ps + Cyr*rs)*b_2V + Cyda*aileron + Cydr*rudder;
 290
 291   cm = Cmo + Cma*Alpha + (Cmq*Theta_dot + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator+long_trim);
 292   cn = Cnbeta*Beta + (Cnp*ps + Cnr*rs)*b_2V + Cnda*aileron + Cndr*rudder;
 293   croll=Clbeta*Beta + (Clp*ps + Clr*rs)*b_2V + Clda*aileron + Cldr*rudder;
 294
 295 /*   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);
 296  */  /*calculate wind axes forces*/
 297   F_X_wind=-1*cd*qS;
 298   F_Y_wind=cy*qS;
 299   F_Z_wind=-1*CL*qS;
 300
 301 /*   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);
 302  */
 303   /*calculate moments and body axis forces */
 304
 305   /*find body-axis components of weight*/
 306   /*with earth axis to body axis transform */
 307   W_X=-1*W*Sin_theta;
 308   W_Y=W*Sin_phi*Cos_theta;
 309   W_Z=W*Cos_phi*Cos_theta;
 310
 311   /* requires ugly wind-axes to body-axes transform */
 312   F_X_aero = W_X + F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
 313   F_Y_aero = W_Y + F_X_wind*Sin_beta + F_Y_wind*Cos_beta;
 314   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;
 315
 316   /*no axes transform here */
 317   M_l_aero = croll*qSb;
 318   M_m_aero = cm*qScbar;
 319   M_n_aero = cn*qSb;
 320
 321 /*   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);
 322   printf("Fxbody: %7.4f Fybody: %7.4f Fzbody: %7.4f Weight: %7.4f\n",F_X_wind,F_Y_wind,F_Z_wind,W);
 323   printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);
 324  */
 325
 326 }
 327
 328