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