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 9
  98 #define Ndf 4
  99 #define DYN_ON_SPEED 33 /*20 knots*/
 100
 101
 102 #ifdef USENZ
 103         #define NZ generic_.n_cg_body_v[2]
 104 #else
 105         #define NZ 1
 106 #endif
 107
 108
 109 extern COCKPIT cockpit_;
 110
 111
 112 SCALAR interp(SCALAR *y_table, SCALAR *x_table, int Ntable, SCALAR x)
 113 {
 114         SCALAR slope;
 115         int i=1;
 116         float y;
 117
 118
 119         /* if x is outside the table, return value at x[0] or x[Ntable-1]*/
 120         if(x <= x_table[0])
 121         {
 122                  y=y_table[0];
 123                  /* printf("x smaller than x_table[0]: %g %g\n",x,x_table[0]); */
 124         }
 125         else if(x >= x_table[Ntable-1])
 126         {
 127                 slope=(y_table[Ntable-1]-y_table[Ntable-2])/(x_table[Ntable-1]-x_table[Ntable-2]);
 128             y=slope*(x-x_table[Ntable-1]) +y_table[Ntable-1];
 129
 130 /*               printf("x larger than x_table[N]: %g %g %d\n",x,x_table[NCL-1],Ntable-1);
 131  */     }
 132         else /*x is within the table, interpolate linearly to find y value*/
 133         {
 134
 135             while(x_table[i] <= x) {i++;}
 136             slope=(y_table[i]-y_table[i-1])/(x_table[i]-x_table[i-1]);
 137                 /* 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); */
 138             y=slope*(x-x_table[i-1]) +y_table[i-1];
 139         }
 140         return y;
 141 }
 142
 143
 144 void aero( SCALAR dt, int Initialize ) {
 145
 146
 147   static int init = 0;
 148   static int flap_dir=0;
 149   static SCALAR lastFlapHandle=0;
 150   static SCALAR Ai;
 151
 152   static SCALAR trim_inc = 0.0002;
 153
 154   static SCALAR alpha_ind[NCL]={-0.087,0,0.14,0.21,0.24,0.26,0.28,0.31,0.35};
 155   static SCALAR CLtable[NCL]={-0.22,0.25,1.02,1.252,1.354,1.44,1.466,1.298,0.97};
 156
 157   static SCALAR flap_ind[Ndf]={0,10,20,30};
 158   static SCALAR dCLf[Ndf]={0,0.20,0.30,0.35};
 159   static SCALAR dCdf[Ndf]={0,0.0021,0.0085,0.0191};
 160   static SCALAR dCmf[Ndf]={0,-0.0654,-0.0981,-0.114};
 161
 162   static SCALAR flap_transit_rate=2.5;
 163
 164
 165
 166
 167
 168    /* printf("Initialize= %d\n",Initialize); */
 169 /*         printf("Initializing aero model...Initialize= %d\n", Initialize);
 170  */
 171         /*nondimensionalization quantities*/
 172            /*units here are ft and lbs */
 173            cbar=4.9; /*mean aero chord ft*/
 174            b=35.8; /*wing span ft */
 175            Sw=174; /*wing planform surface area ft^2*/
 176            rPiARe=0.054; /*reciprocal of Pi*AR*e*/
 177            lbare=3.682;  /*elevator moment arm / MAC*/
 178
 179            CLadot=1.7;
 180            CLq=3.9;
 181
 182            CLob=0;
 183
 184        Ai=1.24;
 185            Cdob=0.036;
 186            Cda=0.13;  /*Not used*/
 187            Cdde=0.06;
 188
 189            Cma=-1.8;
 190            Cmadot=-5.2;
 191            Cmq=-12.4;
 192            Cmob=-0.02;
 193            Cmde=-1.28;
 194
 195            CLde=-Cmde / lbare; /* kinda backwards, huh? */
 196
 197            Clbeta=-0.089;
 198            Clp=-0.47;
 199            Clr=0.096;
 200            Clda=-0.09;
 201            Cldr=0.0147;
 202
 203            Cnbeta=0.065;
 204            Cnp=-0.03;
 205            Cnr=-0.099;
 206            Cnda=-0.0053;
 207            Cndr=-0.0657;
 208
 209            Cybeta=-0.31;
 210            Cyp=-0.037;
 211            Cyr=0.21;
 212            Cyda=0.0;
 213            Cydr=0.187;
 214
 215
 216
 217            MaxTakeoffWeight=2550;
 218            EmptyWeight=1500;
 219
 220            Zcg=0.51;
 221
 222   /*
 223   LaRCsim uses:
 224     Cm > 0 => ANU
 225         Cl > 0 => Right wing down
 226         Cn > 0 => ANL
 227   so:
 228     elevator > 0 => AND -- aircraft nose down
 229         aileron > 0  => right wing up
 230         rudder > 0   => ANL
 231   */
 232
 233   /*do weight & balance here since there is no better place*/
 234   Weight=Mass / INVG;
 235
 236   if(Weight > 2550)
 237   {  Weight=2550; }
 238   else if(Weight < 1500)
 239   {  Weight=1500; }
 240
 241
 242   if(Dx_cg > 0.5586)
 243   {  Dx_cg = 0.5586; }
 244   else if(Dx_cg < -0.4655)
 245   {  Dx_cg = -0.4655; }
 246
 247   Cg=Dx_cg/cbar +0.25;
 248
 249   Dz_cg=Zcg*cbar;
 250
 251
 252   if(Flap_handle < flap_ind[0])
 253   {
 254         Flap_handle=flap_ind[0];
 255         Flap_Position=flap_ind[0];
 256   }
 257   else if(Flap_handle > flap_ind[3])
 258   {
 259          Flap_handle=flap_ind[3];
 260          Flap_Position=flap_ind[3];
 261   }
 262   else
 263   {
 264
 265
 266          if((Flap_handle != lastFlapHandle) && (dt > 0))
 267          {
 268                 Flaps_In_Transit=1;
 269
 270          }
 271          else if(dt <= 0)
 272                 Flap_Position=Flap_handle;
 273
 274          lastFlapHandle=Flap_handle;
 275          if((Flaps_In_Transit) && (dt > 0))
 276          {
 277                 if(Flap_Position < 10)
 278                         flap_transit_rate = 2.5;
 279                 else
 280                         flap_transit_rate=5;
 281
 282                 if(Flaps_In_Transit)
 283                 {
 284                    if(Flap_Position < Flap_handle)
 285                         flap_dir=1;
 286                    else
 287                         flap_dir=-1;
 288
 289                    if(fabs(Flap_Position - Flap_handle) > dt*flap_transit_rate)
 290                            Flap_Position+=flap_dir*flap_transit_rate*dt;
 291
 292                    if(fabs(Flap_Position - Flap_handle) < dt*flap_transit_rate)
 293                    {
 294                            Flaps_In_Transit=0;
 295                            Flap_Position=Flap_handle;
 296                    }
 297         }
 298           }
 299   }
 300
 301   if(Aft_trim) long_trim = long_trim - trim_inc;
 302   if(Fwd_trim) long_trim = long_trim + trim_inc;
 303
 304 /*   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);
 305  */  /*scale pct control to degrees deflection*/
 306   if ((Long_control+Long_trim) <= 0)
 307         elevator=(Long_control+Long_trim)*28*DEG_TO_RAD;
 308   else
 309         elevator=(Long_control+Long_trim)*23*DEG_TO_RAD;
 310
 311   aileron  = -1*Lat_control*17.5*DEG_TO_RAD;
 312   rudder   = -1*Rudder_pedal*16*DEG_TO_RAD;
 313   /*
 314     The aileron travel limits are 20 deg. TEU and 15 deg TED
 315     but since we don't distinguish between left and right we'll
 316         use the average here (17.5 deg)
 317   */
 318
 319
 320   /*calculate rate derivative nondimensionalization (is that a word?) factors */
 321   /*hack to avoid divide by zero*/
 322   /*the dynamic terms are negligible at low ground speeds anyway*/
 323
 324 /*   printf("Vinf: %g, Span: %g\n",V_rel_wind,b);
 325  */
 326   if(V_rel_wind > DYN_ON_SPEED)
 327   {
 328         cbar_2V=cbar/(2*V_rel_wind);
 329         b_2V=b/(2*V_rel_wind);
 330   }
 331   else
 332   {
 333         cbar_2V=0;
 334         b_2V=0;
 335   }
 336
 337
 338   /*calcuate the qS nondimensionalization factors*/
 339
 340   qS=Dynamic_pressure*Sw;
 341   qScbar=qS*cbar;
 342   qSb=qS*b;
 343
 344
 345 /*   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);
 346   printf("aero: 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);
 347  */
 348
 349
 350
 351   /* sum coefficients */
 352   CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
 353 /*   printf("CLwbh: %g\n",CLwbh);
 354  */
 355   CLo = CLob + interp(dCLf,flap_ind,Ndf,Flap_Position);
 356   Cdo = Cdob + interp(dCdf,flap_ind,Ndf,Flap_Position);
 357   Cmo = Cmob + interp(dCmf,flap_ind,Ndf,Flap_Position);
 358
 359   /* printf("FP: %g\n",Flap_Position);
 360   printf("CLo: %g\n",CLo);
 361   printf("Cdo: %g\n",Cdo);
 362   printf("Cmo: %g\n",Cmo);        */
 363
 364
 365
 366
 367
 368   CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
 369   cd = Cdo + rPiARe*Ai*Ai*CL*CL + Cdde*elevator;
 370   cy = Cybeta*Beta + (Cyp*P_body + Cyr*R_body)*b_2V + Cyda*aileron + Cydr*rudder;
 371
 372   cm = Cmo + Cma*Alpha + (Cmq*Q_body + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator);
 373   cn = Cnbeta*Beta + (Cnp*P_body + Cnr*R_body)*b_2V + Cnda*aileron + Cndr*rudder;
 374   croll=Clbeta*Beta + (Clp*P_body + Clr*R_body)*b_2V + Clda*aileron + Cldr*rudder;
 375
 376 /*   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);
 377  */
 378
 379   /*calculate wind axes forces*/
 380   F_X_wind=-1*cd*qS;
 381   F_Y_wind=cy*qS;
 382   F_Z_wind=-1*CL*qS;
 383
 384 /*   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);
 385  */
 386
 387   /*calculate moments and body axis forces */
 388
 389
 390
 391   /* requires ugly wind-axes to body-axes transform */
 392   F_X_aero = F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
 393   F_Y_aero = F_X_wind*Sin_beta + F_Y_wind*Cos_beta;
 394   F_Z_aero = F_X_wind*Sin_alpha*Cos_beta - F_Y_wind*Sin_alpha*Sin_beta + F_Z_wind*Cos_alpha;
 395
 396   /*no axes transform here */
 397   M_l_aero = croll*qSb;
 398   M_m_aero = cm*qScbar;
 399   M_n_aero = cn*qSb;
 400
 401 /*   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);
 402
 403   printf("Fxaero: %7.4f Fyaero: %7.4f Fzaero: %7.4f Weight: %7.4f\n",F_X_aero,F_Y_aero,F_Z_aero,Weight);
 404
 405   printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);
 406   */
 407
 408 }
 409
 410