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    /* printf("Initialize= %d\n",Initialize); */
 156    /* if (Initialize != 0)
 157     { */
 158 /*         printf("Initializing aero model...Initialize= %d\n", Initialize);
 159  */        CLadot=1.7;
 160            CLq=3.9;
 161            CLde=0.43;
 162            CLo=0;
 163
 164
 165            Cdo=0.031;
 166            Cda=0.13;  /*Not used*/
 167            Cdde=0.06;
 168
 169            Cma=-0.89;
 170            Cmadot=-5.2;
 171            Cmq=-12.4;
 172            Cmo=-0.062;
 173            Cmde=-1.28;
 174
 175            Clbeta=-0.089;
 176            Clp=-0.47;
 177            Clr=0.096;
 178            Clda=-0.178;
 179            Cldr=0.0147;
 180
 181            Cnbeta=0.065;
 182            Cnp=-0.03;
 183            Cnr=-0.099;
 184            Cnda=-0.053;
 185            Cndr=-0.0657;
 186
 187            Cybeta=-0.31;
 188            Cyp=-0.037;
 189            Cyr=0.21;
 190            Cyda=0.0;
 191            Cydr=0.187;
 192
 193            /*nondimensionalization quantities*/
 194            /*units here are ft and lbs */
 195            cbar=4.9; /*mean aero chord ft*/
 196            b=35.8; /*wing span ft */
 197            Sw=174; /*wing planform surface area ft^2*/
 198            rPiARe=0.054; /*reciprocal of Pi*AR*e*/
 199     /* } */
 200
 201   /*
 202   LaRCsim uses:
 203     Cm > 0 => ANU
 204         Cl > 0 => Right wing down
 205         Cn > 0 => ANL
 206   so:
 207     elevator > 0 => AND -- aircraft nose down
 208         aileron > 0  => right wing up
 209         rudder > 0   => ANL
 210   */
 211   long_trim=0;
 212   if(Aft_trim) long_trim = long_trim - trim_inc;
 213   if(Fwd_trim) long_trim = long_trim + trim_inc;
 214
 215 /*   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);
 216  */  /*scale pct control to degrees deflection*/
 217   if ((Long_control+long_trim) <= 0)
 218         elevator=(Long_control+long_trim)*28*DEG_TO_RAD;
 219   else
 220         elevator=(Long_control+long_trim)*23*DEG_TO_RAD;
 221
 222   aileron  = Lat_control*17.5*DEG_TO_RAD;
 223   rudder   = Rudder_pedal*16*DEG_TO_RAD;
 224   /*
 225     The aileron travel limits are 20 deg. TEU and 15 deg TED
 226     but since we don't distinguish between left and right we'll
 227         use the average here (17.5 deg)
 228   */
 229
 230
 231   /*calculate rate derivative nondimensionalization (is that a word?) factors */
 232   /*hack to avoid divide by zero*/
 233   /*the dynamic terms might be negligible at low ground speeds anyway*/
 234
 235   if(V_rel_wind > DYN_ON_SPEED)
 236   {
 237         cbar_2V=cbar/(2*V_rel_wind);
 238         b_2V=b/(2*V_rel_wind);
 239   }
 240   else
 241   {
 242         cbar_2V=0;
 243         b_2V=0;
 244   }
 245
 246
 247   /*calcuate the qS nondimensionalization factors*/
 248
 249   qS=Dynamic_pressure*Sw;
 250   qScbar=qS*cbar;
 251   qSb=qS*b;
 252
 253
 254 /*   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);
 255  */  //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);
 256
 257
 258   /* sum coefficients */
 259   CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
 260   CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
 261   cd = Cdo + rPiARe*CL*CL + Cdde*elevator;
 262   cy = Cybeta*Beta + (Cyp*P_body + Cyr*R_body)*b_2V + Cyda*aileron + Cydr*rudder;
 263
 264   cm = Cmo + Cma*Alpha + (Cmq*Q_body + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator+long_trim);
 265   cn = Cnbeta*Beta + (Cnp*P_body + Cnr*R_body)*b_2V + Cnda*aileron + Cndr*rudder;
 266   croll=Clbeta*Beta + (Clp*P_body + Clr*R_body)*b_2V + Clda*aileron + Cldr*rudder;
 267
 268 /*   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);
 269  */  /*calculate wind axes forces*/
 270   F_X_wind=-1*cd*qS;
 271   F_Y_wind=cy*qS;
 272   F_Z_wind=-1*CL*qS;
 273
 274 /*   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);
 275  */
 276   /*calculate moments and body axis forces */
 277
 278
 279
 280   /* requires ugly wind-axes to body-axes transform */
 281   F_X_aero = F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
 282   F_Y_aero = F_X_wind*Sin_beta + F_Y_wind*Cos_beta;
 283   F_Z_aero = F_X_wind*Sin_alpha*Cos_beta - F_Y_wind*Sin_alpha*Sin_beta + F_Z_wind*Cos_alpha;
 284
 285   /*no axes transform here */
 286   M_l_aero = croll*qSb;
 287   M_m_aero = cm*qScbar;
 288   M_n_aero = cn*qSb;
 289
 290 /*   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);
 291  */
 292 /*  printf("Fxaero: %7.4f Fyaero: %7.4f Fzaero: %7.4f Weight: %7.4f\n",F_X_aero,F_Y_aero,F_Z_aero,W);
 293  *//*  printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);
 294  */
 295 }