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            Cdob=0.031;
 185            Cda=0.13;  /*Not used*/
 186            Cdde=0.06;
 187
 188            Cma=-1.8;
 189            Cmadot=-5.2;
 190            Cmq=-12.4;
 191            Cmob=-0.02;
 192            Cmde=-1.28;
 193
 194            CLde=-Cmde / lbare; /* kinda backwards, huh? */
 195
 196            Clbeta=-0.089;
 197            Clp=-0.47;
 198            Clr=0.096;
 199            Clda=-0.09;
 200            Cldr=0.0147;
 201
 202            Cnbeta=0.065;
 203            Cnp=-0.03;
 204            Cnr=-0.099;
 205            Cnda=-0.0053;
 206            Cndr=-0.0657;
 207
 208            Cybeta=-0.31;
 209            Cyp=-0.037;
 210            Cyr=0.21;
 211            Cyda=0.0;
 212            Cydr=0.187;
 213
 214
 215
 216            MaxTakeoffWeight=2550;
 217            EmptyWeight=1500;
 218
 219            Zcg=0.51;
 220
 221   /*
 222   LaRCsim uses:
 223     Cm > 0 => ANU
 224         Cl > 0 => Right wing down
 225         Cn > 0 => ANL
 226   so:
 227     elevator > 0 => AND -- aircraft nose down
 228         aileron > 0  => right wing up
 229         rudder > 0   => ANL
 230   */
 231
 232   /*do weight & balance here since there is no better place*/
 233   Weight=Mass / INVG;
 234
 235   if(Weight > 2550)
 236   {  Weight=2550; }
 237   else if(Weight < 1500)
 238   {  Weight=1500; }
 239
 240
 241   if(Dx_cg > 0.5586)
 242   {  Dx_cg = 0.5586; }
 243   else if(Dx_cg < -0.4655)
 244   {  Dx_cg = -0.4655; }
 245
 246   Cg=Dx_cg/cbar +0.25;
 247
 248   Dz_cg=Zcg*cbar;
 249
 250
 251   if(Flap_handle < flap_ind[0])
 252   {
 253         Flap_handle=flap_ind[0];
 254         Flap_Position=flap_ind[0];
 255   }
 256   else if(Flap_handle > flap_ind[3])
 257   {
 258          Flap_handle=flap_ind[3];
 259          Flap_Position=flap_ind[3];
 260   }
 261   else
 262   {
 263
 264
 265          if((Flap_handle != lastFlapHandle) && (dt > 0))
 266                 Flaps_In_Transit=1;
 267          else if(dt <= 0)
 268                 Flap_Position=Flap_handle;
 269
 270          lastFlapHandle=Flap_handle;
 271          if((Flaps_In_Transit) && (dt > 0))
 272          {
 273                 if(Flap_Position < 10)
 274                         flap_transit_rate = 2.5;
 275                 else
 276                         flap_transit_rate=5;
 277
 278                 if(Flaps_In_Transit)
 279                 {
 280                    if(Flap_Position < Flap_handle)
 281                         flap_dir=1;
 282                    else
 283                         flap_dir=-1;
 284
 285                    if(fabs(Flap_Position - Flap_handle) > dt*flap_transit_rate)
 286                            Flap_Position+=flap_dir*flap_transit_rate*dt;
 287
 288                    if(fabs(Flap_Position - Flap_handle) < dt*flap_transit_rate)
 289                    {
 290                            Flaps_In_Transit=0;
 291                            Flap_Position=Flap_handle;
 292                    }
 293         }
 294           }
 295   }
 296
 297   if(Aft_trim) long_trim = long_trim - trim_inc;
 298   if(Fwd_trim) long_trim = long_trim + trim_inc;
 299
 300 /*   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);
 301  */  /*scale pct control to degrees deflection*/
 302   if ((Long_control+Long_trim) <= 0)
 303         elevator=(Long_control+Long_trim)*28*DEG_TO_RAD;
 304   else
 305         elevator=(Long_control+Long_trim)*23*DEG_TO_RAD;
 306
 307   aileron  = -1*Lat_control*17.5*DEG_TO_RAD;
 308   rudder   = -1*Rudder_pedal*16*DEG_TO_RAD;
 309   /*
 310     The aileron travel limits are 20 deg. TEU and 15 deg TED
 311     but since we don't distinguish between left and right we'll
 312         use the average here (17.5 deg)
 313   */
 314
 315
 316   /*calculate rate derivative nondimensionalization (is that a word?) factors */
 317   /*hack to avoid divide by zero*/
 318   /*the dynamic terms might be negligible at low ground speeds anyway*/
 319
 320 /*   printf("Vinf: %g, Span: %g\n",V_rel_wind,b);
 321  */
 322   if(V_rel_wind > DYN_ON_SPEED)
 323   {
 324         cbar_2V=cbar/(2*V_rel_wind);
 325         b_2V=b/(2*V_rel_wind);
 326   }
 327   else
 328   {
 329         cbar_2V=0;
 330         b_2V=0;
 331   }
 332
 333
 334   /*calcuate the qS nondimensionalization factors*/
 335
 336   qS=Dynamic_pressure*Sw;
 337   qScbar=qS*cbar;
 338   qSb=qS*b;
 339
 340
 341 /*   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);
 342   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);
 343  */
 344
 345
 346
 347   /* sum coefficients */
 348   CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
 349 /*   printf("CLwbh: %g\n",CLwbh);
 350  */
 351   CLo = CLob + interp(dCLf,flap_ind,Ndf,Flap_Position);
 352   Cdo = Cdob + interp(dCdf,flap_ind,Ndf,Flap_Position);
 353   Cmo = Cmob + interp(dCmf,flap_ind,Ndf,Flap_Position);
 354
 355   /* printf("FP: %g\n",Flap_Position);
 356   printf("CLo: %g\n",CLo);
 357   printf("Cdo: %g\n",Cdo);
 358   printf("Cmo: %g\n",Cmo); */
 359
 360
 361
 362   CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
 363   cd = Cdo + rPiARe*CL*CL + Cdde*elevator;
 364   cy = Cybeta*Beta + (Cyp*P_body + Cyr*R_body)*b_2V + Cyda*aileron + Cydr*rudder;
 365
 366   cm = Cmo + Cma*Alpha + (Cmq*Q_body + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator);
 367   cn = Cnbeta*Beta + (Cnp*P_body + Cnr*R_body)*b_2V + Cnda*aileron + Cndr*rudder;
 368   croll=Clbeta*Beta + (Clp*P_body + Clr*R_body)*b_2V + Clda*aileron + Cldr*rudder;
 369
 370 /*   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);
 371  */
 372
 373   /*calculate wind axes forces*/
 374   F_X_wind=-1*cd*qS;
 375   F_Y_wind=cy*qS;
 376   F_Z_wind=-1*CL*qS;
 377
 378 /*   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);
 379  */
 380
 381   /*calculate moments and body axis forces */
 382
 383
 384
 385   /* requires ugly wind-axes to body-axes transform */
 386   F_X_aero = F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
 387   F_Y_aero = F_X_wind*Sin_beta + F_Y_wind*Cos_beta;
 388   F_Z_aero = F_X_wind*Sin_alpha*Cos_beta - F_Y_wind*Sin_alpha*Sin_beta + F_Z_wind*Cos_alpha;
 389
 390   /*no axes transform here */
 391   M_l_aero = croll*qSb;
 392   M_m_aero = cm*qScbar;
 393   M_n_aero = cn*qSb;
 394
 395 /*   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);
 396
 397   printf("Fxaero: %7.4f Fyaero: %7.4f Fzaero: %7.4f Weight: %7.4f\n",F_X_aero,F_Y_aero,F_Z_aero,Weight);
 398
 399   printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);
 400   */
 401
 402 }
 403
 404