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.43)
 243   {  Dx_cg = 0.43; }
 244   else if(Dx_cg < -0.6)
 245   {  Dx_cg = -0.6; }
 246
 247   Cg=0.25 - Dx_cg/cbar;
 248
 249   Dz_cg=Zcg*cbar;
 250   Dy_cg=0;
 251
 252
 253   if(Flap_handle < flap_ind[0])
 254   {
 255         Flap_handle=flap_ind[0];
 256         Flap_Position=flap_ind[0];
 257   }
 258   else if(Flap_handle > flap_ind[3])
 259   {
 260          Flap_handle=flap_ind[3];
 261          Flap_Position=flap_ind[3];
 262   }
 263   else
 264   {
 265
 266
 267          if((Flap_handle != lastFlapHandle) && (dt > 0))
 268          {
 269                 Flaps_In_Transit=1;
 270
 271          }
 272          else if(dt <= 0)
 273                 Flap_Position=Flap_handle;
 274
 275          lastFlapHandle=Flap_handle;
 276          if((Flaps_In_Transit) && (dt > 0))
 277          {
 278                 if(Flap_Position < 10)
 279                         flap_transit_rate = 2.5;
 280                 else
 281                         flap_transit_rate=5;
 282
 283                 if(Flaps_In_Transit)
 284                 {
 285                    if(Flap_Position < Flap_handle)
 286                         flap_dir=1;
 287                    else
 288                         flap_dir=-1;
 289
 290                    if(fabs(Flap_Position - Flap_handle) > dt*flap_transit_rate)
 291                            Flap_Position+=flap_dir*flap_transit_rate*dt;
 292
 293                    if(fabs(Flap_Position - Flap_handle) < dt*flap_transit_rate)
 294                    {
 295                            Flaps_In_Transit=0;
 296                            Flap_Position=Flap_handle;
 297                    }
 298         }
 299           }
 300   }
 301
 302   if(Aft_trim) long_trim = long_trim - trim_inc;
 303   if(Fwd_trim) long_trim = long_trim + trim_inc;
 304
 305 /*   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);
 306  */  /*scale pct control to degrees deflection*/
 307   if ((Long_control+Long_trim) <= 0)
 308         elevator=(Long_control+Long_trim)*28*DEG_TO_RAD;
 309   else
 310         elevator=(Long_control+Long_trim)*23*DEG_TO_RAD;
 311
 312   aileron  = -1*Lat_control*17.5*DEG_TO_RAD;
 313   rudder   = -1*Rudder_pedal*16*DEG_TO_RAD;
 314   /*
 315     The aileron travel limits are 20 deg. TEU and 15 deg TED
 316     but since we don't distinguish between left and right we'll
 317         use the average here (17.5 deg)
 318   */
 319
 320
 321   /*calculate rate derivative nondimensionalization (is that a word?) factors */
 322   /*hack to avoid divide by zero*/
 323   /*the dynamic terms are negligible at low ground speeds anyway*/
 324
 325 /*   printf("Vinf: %g, Span: %g\n",V_rel_wind,b);
 326  */
 327   if(V_rel_wind > DYN_ON_SPEED)
 328   {
 329         cbar_2V=cbar/(2*V_rel_wind);
 330         b_2V=b/(2*V_rel_wind);
 331   }
 332   else
 333   {
 334         cbar_2V=0;
 335         b_2V=0;
 336   }
 337
 338
 339   /*calcuate the qS nondimensionalization factors*/
 340
 341   qS=Dynamic_pressure*Sw;
 342   qScbar=qS*cbar;
 343   qSb=qS*b;
 344
 345
 346 /*   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);
 347   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);
 348  */
 349
 350
 351
 352   /* sum coefficients */
 353   CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
 354 /*   printf("CLwbh: %g\n",CLwbh);
 355  */
 356   CLo = CLob + interp(dCLf,flap_ind,Ndf,Flap_Position);
 357   Cdo = Cdob + interp(dCdf,flap_ind,Ndf,Flap_Position);
 358   Cmo = Cmob + interp(dCmf,flap_ind,Ndf,Flap_Position);
 359
 360   /* printf("FP: %g\n",Flap_Position);
 361   printf("CLo: %g\n",CLo);
 362   printf("Cdo: %g\n",Cdo);
 363   printf("Cmo: %g\n",Cmo);        */
 364
 365
 366
 367
 368
 369   CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
 370   cd = Cdo + rPiARe*Ai*Ai*CL*CL + Cdde*elevator;
 371   cy = Cybeta*Beta + (Cyp*P_body + Cyr*R_body)*b_2V + Cyda*aileron + Cydr*rudder;
 372
 373   cm = Cmo + Cma*Alpha + (Cmq*Q_body + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator);
 374   cn = Cnbeta*Beta + (Cnp*P_body + Cnr*R_body)*b_2V + Cnda*aileron + Cndr*rudder;
 375   croll=Clbeta*Beta + (Clp*P_body + Clr*R_body)*b_2V + Clda*aileron + Cldr*rudder;
 376
 377 /*   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);
 378  */
 379
 380   /*calculate wind axes forces*/
 381   F_X_wind=-1*cd*qS;
 382   F_Y_wind=cy*qS;
 383   F_Z_wind=-1*CL*qS;
 384
 385 /*   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);
 386  */
 387
 388   /*calculate moments and body axis forces */
 389
 390
 391
 392   /* requires ugly wind-axes to body-axes transform */
 393   F_X_aero = F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
 394   F_Y_aero = F_X_wind*Sin_beta + F_Y_wind*Cos_beta;
 395   F_Z_aero = F_X_wind*Sin_alpha*Cos_beta - F_Y_wind*Sin_alpha*Sin_beta + F_Z_wind*Cos_alpha;
 396
 397   /*no axes transform here */
 398   M_l_aero = croll*qSb;
 399   M_m_aero = cm*qScbar;
 400   M_n_aero = cn*qSb;
 401
 402 /*   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);
 403
 404   printf("Fxaero: %7.4f Fyaero: %7.4f Fzaero: %7.4f Weight: %7.4f\n",F_X_aero,F_Y_aero,F_Z_aero,Weight);
 405
 406   printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);
 407   */
 408
 409 }
 410
 411