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
 151   static SCALAR trim_inc = 0.0002;
 152
 153   static SCALAR alpha_ind[NCL]={-0.087,0,0.14,0.21,0.24,0.26,0.28,0.31,0.35};
 154   static SCALAR CLtable[NCL]={-0.22,0.25,1.02,1.252,1.354,1.44,1.466,1.298,0.97};
 155
 156   static SCALAR flap_ind[Ndf]={0,10,20,30};
 157   static SCALAR dCLf[Ndf]={0,0.20,0.30,0.35};
 158   static SCALAR dCdf[Ndf]={0,0.0021,0.0085,0.0191};
 159   static SCALAR dCmf[Ndf]={0,-0.0654,-0.0981,-0.114};
 160
 161   static SCALAR flap_transit_rate=2.5;
 162
 163
 164
 165
 166
 167    /* printf("Initialize= %d\n",Initialize); */
 168 /*         printf("Initializing aero model...Initialize= %d\n", Initialize);
 169  */
 170         /*nondimensionalization quantities*/
 171            /*units here are ft and lbs */
 172            cbar=4.9; /*mean aero chord ft*/
 173            b=35.8; /*wing span ft */
 174            Sw=174; /*wing planform surface area ft^2*/
 175            rPiARe=0.054; /*reciprocal of Pi*AR*e*/
 176            lbare=3.682;  /*elevator moment arm / MAC*/
 177
 178            CLadot=1.7;
 179            CLq=3.9;
 180
 181            CLob=0;
 182
 183
 184            /* original */ /* Cdob=0.031; */
 185            Cdob=0.046;
 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                 Flaps_In_Transit=1;
 268          } else if(dt <= 0) {
 269                 Flap_Position=Flap_handle;
 270          }
 271
 272          lastFlapHandle=Flap_handle;
 273          if((Flaps_In_Transit) && (dt > 0))
 274          {
 275                 if(Flap_Position < 10)
 276                         flap_transit_rate = 2.5;
 277                 else
 278                         flap_transit_rate=5;
 279
 280                 if(Flaps_In_Transit)
 281                 {
 282                    if(Flap_Position < Flap_handle)
 283                         flap_dir=1;
 284                    else
 285                         flap_dir=-1;
 286
 287                    if(fabs(Flap_Position - Flap_handle) > dt*flap_transit_rate)
 288                            Flap_Position+=flap_dir*flap_transit_rate*dt;
 289
 290                    if(fabs(Flap_Position - Flap_handle) < dt*flap_transit_rate)
 291                    {
 292                            Flaps_In_Transit=0;
 293                            Flap_Position=Flap_handle;
 294                    }
 295         }
 296           }
 297   }
 298
 299   if(Aft_trim) long_trim = long_trim - trim_inc;
 300   if(Fwd_trim) long_trim = long_trim + trim_inc;
 301
 302 /*   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);
 303  */  /*scale pct control to degrees deflection*/
 304   if ((Long_control+Long_trim) <= 0)
 305         elevator=(Long_control+Long_trim)*28*DEG_TO_RAD;
 306   else
 307         elevator=(Long_control+Long_trim)*23*DEG_TO_RAD;
 308
 309   aileron  = -1*Lat_control*17.5*DEG_TO_RAD;
 310   rudder   = -1*Rudder_pedal*16*DEG_TO_RAD;
 311   /*
 312     The aileron travel limits are 20 deg. TEU and 15 deg TED
 313     but since we don't distinguish between left and right we'll
 314         use the average here (17.5 deg)
 315   */
 316
 317
 318   /*calculate rate derivative nondimensionalization (is that a word?) factors */
 319   /*hack to avoid divide by zero*/
 320   /*the dynamic terms might be negligible at low ground speeds anyway*/
 321
 322 /*   printf("Vinf: %g, Span: %g\n",V_rel_wind,b);
 323  */
 324   if(V_rel_wind > DYN_ON_SPEED)
 325   {
 326         cbar_2V=cbar/(2*V_rel_wind);
 327         b_2V=b/(2*V_rel_wind);
 328   }
 329   else
 330   {
 331         cbar_2V=0;
 332         b_2V=0;
 333   }
 334
 335
 336   /*calcuate the qS nondimensionalization factors*/
 337
 338   qS=Dynamic_pressure*Sw;
 339   qScbar=qS*cbar;
 340   qSb=qS*b;
 341
 342
 343 /*   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);
 344   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);
 345  */
 346
 347
 348
 349   /* sum coefficients */
 350   CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
 351 /*   printf("CLwbh: %g\n",CLwbh);
 352  */
 353   CLo = CLob + interp(dCLf,flap_ind,Ndf,Flap_Position);
 354   Cdo = Cdob + interp(dCdf,flap_ind,Ndf,Flap_Position);
 355   Cmo = Cmob + interp(dCmf,flap_ind,Ndf,Flap_Position);
 356
 357   /* printf("FP: %g\n",Flap_Position);
 358   printf("CLo: %g\n",CLo);
 359   printf("Cdo: %g\n",Cdo);
 360   printf("Cmo: %g\n",Cmo); */
 361
 362
 363
 364   CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
 365   cd = Cdo + rPiARe*CL*CL + Cdde*elevator;
 366   cy = Cybeta*Beta + (Cyp*P_body + Cyr*R_body)*b_2V + Cyda*aileron + Cydr*rudder;
 367
 368   cm = Cmo + Cma*Alpha + (Cmq*Q_body + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator);
 369   cn = Cnbeta*Beta + (Cnp*P_body + Cnr*R_body)*b_2V + Cnda*aileron + Cndr*rudder;
 370   croll=Clbeta*Beta + (Clp*P_body + Clr*R_body)*b_2V + Clda*aileron + Cldr*rudder;
 371
 372 /*   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);
 373  */
 374
 375   /*calculate wind axes forces*/
 376   F_X_wind=-1*cd*qS;
 377   F_Y_wind=cy*qS;
 378   F_Z_wind=-1*CL*qS;
 379
 380 /*   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);
 381  */
 382
 383   /*calculate moments and body axis forces */
 384
 385
 386
 387   /* requires ugly wind-axes to body-axes transform */
 388   F_X_aero = F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
 389   F_Y_aero = F_X_wind*Sin_beta + F_Y_wind*Cos_beta;
 390   F_Z_aero = F_X_wind*Sin_alpha*Cos_beta - F_Y_wind*Sin_alpha*Sin_beta + F_Z_wind*Cos_alpha;
 391
 392   /*no axes transform here */
 393   M_l_aero = croll*qSb;
 394   M_m_aero = cm*qScbar;
 395   M_n_aero = cn*qSb;
 396
 397 /*   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);
 398
 399   printf("Fxaero: %7.4f Fyaero: %7.4f Fzaero: %7.4f Weight: %7.4f\n",F_X_aero,F_Y_aero,F_Z_aero,Weight);
 400
 401   printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);
 402   */
 403
 404 }
 405
 406