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   long_trim=0;
 298   if(Aft_trim) long_trim = long_trim - trim_inc;
 299   if(Fwd_trim) long_trim = long_trim + trim_inc;
 300
 301 /*   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);
 302  */  /*scale pct control to degrees deflection*/
 303   if ((Long_control+long_trim) <= 0)
 304         elevator=(Long_control+long_trim)*28*DEG_TO_RAD;
 305   else
 306         elevator=(Long_control+long_trim)*23*DEG_TO_RAD;
 307
 308   aileron  = -1*Lat_control*17.5*DEG_TO_RAD;
 309   rudder   = -1*Rudder_pedal*16*DEG_TO_RAD;
 310   /*
 311     The aileron travel limits are 20 deg. TEU and 15 deg TED
 312     but since we don't distinguish between left and right we'll
 313         use the average here (17.5 deg)
 314   */
 315
 316
 317   /*calculate rate derivative nondimensionalization (is that a word?) factors */
 318   /*hack to avoid divide by zero*/
 319   /*the dynamic terms might be negligible at low ground speeds anyway*/
 320
 321 /*   printf("Vinf: %g, Span: %g\n",V_rel_wind,b);
 322  */
 323   if(V_rel_wind > DYN_ON_SPEED)
 324   {
 325         cbar_2V=cbar/(2*V_rel_wind);
 326         b_2V=b/(2*V_rel_wind);
 327   }
 328   else
 329   {
 330         cbar_2V=0;
 331         b_2V=0;
 332   }
 333
 334
 335   /*calcuate the qS nondimensionalization factors*/
 336
 337   qS=Dynamic_pressure*Sw;
 338   qScbar=qS*cbar;
 339   qSb=qS*b;
 340
 341
 342 /*   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);
 343   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);
 344  */
 345
 346
 347
 348   /* sum coefficients */
 349   CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
 350 /*   printf("CLwbh: %g\n",CLwbh);
 351  */
 352   CLo = CLob + interp(dCLf,flap_ind,Ndf,Flap_Position);
 353   Cdo = Cdob + interp(dCdf,flap_ind,Ndf,Flap_Position);
 354   Cmo = Cmob + interp(dCmf,flap_ind,Ndf,Flap_Position);
 355
 356   /* printf("FP: %g\n",Flap_Position);
 357   printf("CLo: %g\n",CLo);
 358   printf("Cdo: %g\n",Cdo);
 359   printf("Cmo: %g\n",Cmo); */
 360
 361
 362
 363   CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
 364   cd = Cdo + rPiARe*CL*CL + Cdde*elevator;
 365   cy = Cybeta*Beta + (Cyp*P_body + Cyr*R_body)*b_2V + Cyda*aileron + Cydr*rudder;
 366
 367   cm = Cmo + Cma*Alpha + (Cmq*Q_body + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator);
 368   cn = Cnbeta*Beta + (Cnp*P_body + Cnr*R_body)*b_2V + Cnda*aileron + Cndr*rudder;
 369   croll=Clbeta*Beta + (Clp*P_body + Clr*R_body)*b_2V + Clda*aileron + Cldr*rudder;
 370
 371 /*   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);
 372  */
 373
 374   /*calculate wind axes forces*/
 375   F_X_wind=-1*cd*qS;
 376   F_Y_wind=cy*qS;
 377   F_Z_wind=-1*CL*qS;
 378
 379 /*   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);
 380  */
 381
 382   /*calculate moments and body axis forces */
 383
 384
 385
 386   /* requires ugly wind-axes to body-axes transform */
 387   F_X_aero = F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
 388   F_Y_aero = F_X_wind*Sin_beta + F_Y_wind*Cos_beta;
 389   F_Z_aero = F_X_wind*Sin_alpha*Cos_beta - F_Y_wind*Sin_alpha*Sin_beta + F_Z_wind*Cos_alpha;
 390
 391   /*no axes transform here */
 392   M_l_aero = croll*qSb;
 393   M_m_aero = cm*qScbar;
 394   M_n_aero = cn*qSb;
 395
 396 /*   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);
 397
 398   printf("Fxaero: %7.4f Fyaero: %7.4f Fzaero: %7.4f Weight: %7.4f\n",F_X_aero,F_Y_aero,F_Z_aero,Weight);
 399
 400   printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);
 401   */
 402
 403 }
 404
 405