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