src/FDM/LaRCsim/basic_aero.c

   1 /***************************************************************************
   2
   3   TITLE:        basic_aero
   4
   5 ----------------------------------------------------------------------------
   6
   7   FUNCTION:      aerodynamics model based on stability derivatives
   8                  + tweaks for nonlinear aero
   9
  10 ----------------------------------------------------------------------------
  11
  12   MODULE STATUS: developmental
  13
  14 ----------------------------------------------------------------------------
  15
  16   GENEALOGY:     based on aero model from crrcsim code (Drela's aero model)
  17
  18 ----------------------------------------------------------------------------
  19
  20   DESIGNED BY:
  21
  22   CODED BY:             Michael Selig
  23
  24   MAINTAINED BY:        Michael Selig
  25
  26 ----------------------------------------------------------------------------
  27
  28   MODIFICATION HISTORY:
  29
  30   DATE          PURPOSE                                                                                         BY
  31   7/25/03       LaRCsim debugging purposes
  32
  33 ----------------------------------------------------------------------------
  34
  35   REFERENCES:
  36
  37 ----------------------------------------------------------------------------
  38
  39   CALLED BY:
  40
  41 ----------------------------------------------------------------------------
  42
  43   CALLS TO:
  44
  45 ----------------------------------------------------------------------------
  46
  47   INPUTS:
  48
  49 ----------------------------------------------------------------------------
  50
  51   OUTPUTS:
  52
  53 --------------------------------------------------------------------------*/
  54
  55
  56
  57 #include "ls_generic.h"
  58 #include "ls_cockpit.h"
  59 #include "ls_constants.h"
  60 #include "ls_types.h"
  61 #include "basic_aero.h"
  62
  63 #include <math.h>
  64 #include <stdio.h>
  65
  66
  67 #ifdef USENZ
  68 #define NZ generic_.n_cg_body_v[2]
  69 #else
  70 #define NZ 1
  71 #endif
  72
  73
  74 extern COCKPIT cockpit_;
  75
  76
  77 void basic_aero(SCALAR dt, int Initialize)
  78 // Calculate forces and moments for the current time step.  If Initialize is
  79 // zero, then re-initialize coefficients by reading in the coefficient file.
  80 {
  81  static int init = 0;
  82  //static SCALAR elevator_drela, aileron_drela, rudder_drela;
  83
  84  SCALAR C_ref;
  85  SCALAR B_ref;
  86  SCALAR S_ref;
  87  SCALAR U_ref;
  88
  89  /*  SCALAR Mass; */
  90  /*  SCALAR I_xx; */
  91  /*  SCALAR I_yy; */
  92  /*  SCALAR I_zz; */
  93  /*  SCALAR I_xz; */
  94
  95  SCALAR Alpha_0;
  96  SCALAR Cm_0;
  97  SCALAR CL_0;
  98  SCALAR CL_max;
  99  SCALAR CL_min;
 100
 101  SCALAR CD_prof;
 102  SCALAR Uexp_CD;
 103
 104  SCALAR CL_a;
 105  SCALAR Cm_a;
 106
 107  SCALAR CY_b;
 108  SCALAR Cl_b;
 109  SCALAR Cn_b;
 110
 111  SCALAR CL_q;
 112  SCALAR Cm_q;
 113
 114  SCALAR CY_p;
 115  SCALAR Cl_p;
 116  SCALAR Cn_p;
 117
 118  SCALAR CY_r;
 119  SCALAR Cl_r;
 120  SCALAR Cn_r;
 121
 122  SCALAR CL_de;
 123  SCALAR Cm_de;
 124
 125  SCALAR CY_dr;
 126  SCALAR Cl_dr;
 127  SCALAR Cn_dr;
 128
 129  //SCALAR CY_da;
 130  SCALAR Cl_da;
 131  SCALAR Cn_da;
 132
 133  SCALAR eta_loc;
 134  SCALAR CG_arm;
 135  SCALAR CL_drop;
 136  SCALAR CD_stall = 0.05;
 137  //int stalling;
 138
 139  SCALAR span_eff;
 140  SCALAR CL_CD0;
 141  SCALAR CD_CLsq;
 142  SCALAR CD_AIsq;
 143  SCALAR CD_ELsq;
 144
 145  SCALAR Phat, Qhat, Rhat;  // dimensionless rotation rates
 146  SCALAR CL_left, CL_cent, CL_right; // CL values across the span
 147  SCALAR dCL_left,dCL_cent,dCL_right; // stall-induced CL changes
 148  SCALAR dCD_left,dCD_cent,dCD_right; // stall-induced CD changes
 149  SCALAR dCl,dCn,dCl_stall,dCn_stall;  // stall-induced roll,yaw moments
 150  SCALAR dCm_stall;  // Stall induced pitching moment.
 151  SCALAR CL_wing, CD_all, CD_scaled, Cl_w;
 152  SCALAR Cl_r_mod,Cn_p_mod;
 153  SCALAR CL_drela,CD_drela,Cx_drela,Cy_drela,Cz_drela,Cl_drela,Cm_drela,Cn_drela;
 154  SCALAR QS;
 155  /*
 156  SCALAR G_11,G_12,G_13;
 157  SCALAR G_21,G_22,G_23;
 158  SCALAR G_31,G_32,G_33;
 159  SCALAR U_body_X,U_body_Y,U_body_Z;
 160  SCALAR V_body_X,V_body_Y,V_body_Z;
 161  SCALAR W_body_X,W_body_Y,W_body_Z;
 162  */
 163  SCALAR P_atmo,Q_atmo,R_atmo;
 164
 165  // set the parameters
 166
 167  C_ref = 0.551667;
 168  B_ref = 6.55000;
 169  S_ref = 3.61111;
 170  U_ref = 19.6850;
 171  Alpha_0 = 0.349066E-01;
 172  Cm_0 = -0.112663E-01;
 173  CL_0 = 0.563172;
 174  CL_max = 1.10000;
 175  CL_min = -0.600000;
 176  CD_prof = 0.200000E-01;
 177  Uexp_CD = -0.500000;
 178  CL_a = 5.50360;
 179  Cm_a = -0.575335;
 180  CY_b = -0.415610;
 181  Cl_b = -0.250926;
 182  Cn_b = 0.567069E-01;
 183  CL_q = 7.50999;
 184  Cm_q = -11.4975;
 185  CY_p = -0.423820;
 186  Cl_p = -0.611798;
 187  Cn_p = -0.740898E-01;
 188  CY_r = 0.297540;
 189  Cl_r = 0.139581;
 190  Cn_r = -0.687755E-01;
 191  CL_de = 0.162000;
 192  Cm_de = -0.597537;
 193  CY_dr = 0.000000E+00;
 194  Cl_dr = 0.000000E+00;
 195  Cn_dr = 0.000000E+00;
 196  //CY_da = -0.135890;
 197  Cl_da = -0.307921E-02;
 198  Cn_da = 0.527143E-01;
 199  span_eff = 0.95;
 200  CL_CD0 = 0.0;
 201  CD_CLsq = 0.01;
 202  CD_AIsq = 0.0;
 203  CD_ELsq = 0.0;
 204  eta_loc = 0.3;
 205  CG_arm = 0.25;
 206  CL_drop = 0.5;
 207
 208  if (!init)
 209   {
 210    init = -1;
 211   }
 212
 213  // jan's data goes -.5 to .5 while
 214  // fgfs data goes +- 1.
 215  // so I need to divide by 2 below.
 216   elevator = Long_control + Long_trim;
 217   aileron  = Lat_control;
 218   rudder   = Rudder_pedal;
 219
 220   elevator = elevator * 0.5;
 221   aileron  = aileron * 0.5;
 222   rudder   = rudder * 0.5;
 223
 224
 225
 226 //  printf("%f\n",V_rel_wind);
 227
 228   /* compute gradients of Local velocities w.r.t. Body coordinates
 229       G_11  =  dU_local/dx_body
 230       G_12  =  dU_local/dy_body   etc.  */
 231
 232   /*
 233   G_11 = U_atmo_X*T_local_to_body_11
 234        + U_atmo_Y*T_local_to_body_12
 235        + U_atmo_Z*T_local_to_body_13;
 236   G_12 = U_atmo_X*T_local_to_body_21
 237        + U_atmo_Y*T_local_to_body_22
 238        + U_atmo_Z*T_local_to_body_23;
 239   G_13 = U_atmo_X*T_local_to_body_31
 240        + U_atmo_Y*T_local_to_body_32
 241        + U_atmo_Z*T_local_to_body_33;
 242
 243   G_21 = V_atmo_X*T_local_to_body_11
 244        + V_atmo_Y*T_local_to_body_12
 245        + V_atmo_Z*T_local_to_body_13;
 246   G_22 = V_atmo_X*T_local_to_body_21
 247        + V_atmo_Y*T_local_to_body_22
 248        + V_atmo_Z*T_local_to_body_23;
 249   G_23 = V_atmo_X*T_local_to_body_31
 250        + V_atmo_Y*T_local_to_body_32
 251        + V_atmo_Z*T_local_to_body_33;
 252
 253   G_31 = W_atmo_X*T_local_to_body_11
 254        + W_atmo_Y*T_local_to_body_12
 255        + W_atmo_Z*T_local_to_body_13;
 256   G_32 = W_atmo_X*T_local_to_body_21
 257        + W_atmo_Y*T_local_to_body_22
 258        + W_atmo_Z*T_local_to_body_23;
 259   G_33 = W_atmo_X*T_local_to_body_31
 260        + W_atmo_Y*T_local_to_body_32
 261        + W_atmo_Z*T_local_to_body_33;
 262   */
 263
 264   //printf("%f %f %f %f\n",W_atmo_X,W_atmo_Y,G_31,G_32);
 265
 266   /* now compute gradients of Body velocities w.r.t. Body coordinates */
 267   /*  U_body_x  =  dU_body/dx_body   etc.  */
 268
 269   /*
 270   U_body_X = T_local_to_body_11*G_11
 271            + T_local_to_body_12*G_21
 272            + T_local_to_body_13*G_31;
 273   U_body_Y = T_local_to_body_11*G_12
 274            + T_local_to_body_12*G_22
 275            + T_local_to_body_13*G_32;
 276   U_body_Z = T_local_to_body_11*G_13
 277            + T_local_to_body_12*G_23
 278            + T_local_to_body_13*G_33;
 279
 280   V_body_X = T_local_to_body_21*G_11
 281            + T_local_to_body_22*G_21
 282            + T_local_to_body_23*G_31;
 283   V_body_Y = T_local_to_body_21*G_12
 284            + T_local_to_body_22*G_22
 285            + T_local_to_body_23*G_32;
 286   V_body_Z = T_local_to_body_21*G_13
 287            + T_local_to_body_22*G_23
 288            + T_local_to_body_23*G_33;
 289
 290   W_body_X = T_local_to_body_31*G_11
 291            + T_local_to_body_32*G_21
 292            + T_local_to_body_33*G_31;
 293   W_body_Y = T_local_to_body_31*G_12
 294            + T_local_to_body_32*G_22
 295            + T_local_to_body_33*G_32;
 296   W_body_Z = T_local_to_body_31*G_13
 297            + T_local_to_body_32*G_23
 298            + T_local_to_body_33*G_33;
 299   */
 300
 301   /* set rotation rates of airmass motion */
 302   /*   BUG
 303         P_atmo =  W_body_X;
 304         Q_atmo = -W_body_Y;
 305         R_atmo =  V_body_Z;
 306   */
 307
 308   /*
 309         P_atmo =  W_body_Y;
 310         Q_atmo = -W_body_X;
 311         R_atmo =  V_body_X;
 312   */
 313
 314         P_atmo = 0;
 315         Q_atmo = 0;
 316         R_atmo = 0;
 317
 318   //  printf("%f %f %f\n",P_atmo,Q_atmo,R_atmo);
 319   if (V_rel_wind != 0)
 320    {
 321   /* set net effective dimensionless rotation rates */
 322     Phat = (P_body - P_atmo) * B_ref / (2.0*V_rel_wind);
 323     Qhat = (Q_body - Q_atmo) * C_ref / (2.0*V_rel_wind);
 324     Rhat = (R_body - R_atmo) * B_ref / (2.0*V_rel_wind);
 325    }
 326   else
 327    {
 328     Phat=0;
 329     Qhat=0;
 330     Rhat=0;
 331    }
 332
 333 //  printf("Phat: %f Qhat: %f Rhat: %f\n",Phat,Qhat,Rhat);
 334
 335     /* compute local CL at three spanwise locations */
 336   CL_left  = CL_0 + CL_a*(Std_Alpha - Alpha_0 - Phat*eta_loc);
 337   CL_cent  = CL_0 + CL_a*(Std_Alpha - Alpha_0               );
 338   CL_right = CL_0 + CL_a*(Std_Alpha - Alpha_0 + Phat*eta_loc);
 339
 340 // printf("CL_left: %f CL_cent: %f CL_right: %f\n",CL_left,CL_cent,CL_right);
 341
 342     /* set CL-limit changes */
 343   dCL_left  = 0.;
 344   dCL_cent  = 0.;
 345   dCL_right = 0.;
 346
 347   //stalling=0;
 348   if (CL_left  > CL_max)
 349     {
 350       dCL_left  = CL_max-CL_left -CL_drop;
 351       //stalling=1;
 352     }
 353
 354   if (CL_cent  > CL_max)
 355     {
 356       dCL_cent  = CL_max-CL_cent -CL_drop;
 357       //stalling=1;
 358     }
 359
 360   if (CL_right > CL_max)
 361     {
 362       dCL_right = CL_max-CL_right -CL_drop;
 363       //stalling=1;
 364     }
 365
 366   if (CL_left  < CL_min)
 367     {
 368       dCL_left  = CL_min-CL_left -CL_drop;
 369       //stalling=1;
 370     }
 371
 372   if (CL_cent  < CL_min)
 373     {
 374       dCL_cent  = CL_min-CL_cent -CL_drop;
 375       //stalling=1;
 376     }
 377
 378   if (CL_right < CL_min)
 379     {
 380       dCL_right = CL_min-CL_right -CL_drop;
 381       //stalling=1;
 382     }
 383
 384   /* set average wing CL */
 385   CL_wing = CL_0 + CL_a*(Std_Alpha-Alpha_0)
 386           + 0.25*dCL_left + 0.5*dCL_cent + 0.25*dCL_right;
 387
 388 //  printf("CL_wing: %f\n",CL_wing);
 389
 390
 391   /* correct profile CD for CL dependence and aileron dependence */
 392   CD_all = CD_prof
 393          + CD_CLsq*(CL_wing-CL_CD0)*(CL_wing-CL_CD0)
 394          + CD_AIsq*aileron*aileron
 395          + CD_ELsq*elevator*elevator;
 396
 397 //  printf("CD_all:%lf\n",CD_all);
 398
 399     /* scale profile CD with Reynolds number via simple power law */
 400  if (V_rel_wind > 0.1)
 401  {
 402   CD_scaled = CD_all*pow(((double)V_rel_wind/(double)U_ref),Uexp_CD);
 403  }
 404  else
 405  {
 406         CD_scaled=CD_all;
 407  }
 408
 409 //  printf("CD_scaled:%lf\n",CD_scaled);
 410
 411
 412
 413   /* Scale lateral cross-coupling derivatives with wing CL */
 414   Cl_r_mod = Cl_r*CL_wing/CL_0;
 415   Cn_p_mod = Cn_p*CL_wing/CL_0;
 416
 417   //  printf("Cl_r_mod: %f Cn_p_mod: %f\n",Cl_r_mod,Cn_p_mod);
 418
 419     /* total average CD with induced and stall contributions */
 420   dCD_left  = CD_stall*dCL_left *dCL_left ;
 421   dCD_cent  = CD_stall*dCL_cent *dCL_cent ;
 422   dCD_right = CD_stall*dCL_right*dCL_right;
 423
 424   /* total CL, with pitch rate and elevator contributions */
 425   CL_drela = (CL_wing + CL_q*Qhat + CL_de*elevator)*Cos_alpha;
 426
 427 //  printf("CL:%f\n",CL);
 428
 429     /* assymetric stall will cause roll and yaw moments */
 430   dCl =  0.45*-1*(dCL_right-dCL_left)*0.5*eta_loc;
 431   dCn =  0.45*(dCD_right-dCD_left)*0.5*eta_loc;
 432   dCm_stall = (0.25*dCL_left + 0.5*dCL_cent + 0.25*dCL_right)*CG_arm;
 433
 434 //  printf("dCl: %f dCn:%f\n",dCl,dCn);
 435
 436     /* stall-caused moments in body axes */
 437   dCl_stall = dCl*Cos_alpha - dCn*Sin_alpha;
 438   dCn_stall = dCl*Sin_alpha + dCn*Cos_alpha;
 439
 440     /* total CD, with induced and stall contributions */
 441
 442   Cl_w = Cl_b*Std_Beta  + Cl_p*Phat + Cl_r_mod*Rhat
 443        + dCl_stall  + Cl_da*aileron;
 444   CD_drela = CD_scaled
 445      + (CL*CL + 32.0*Cl_w*Cl_w)*S_ref/(B_ref*B_ref*3.14159*span_eff)
 446      + 0.25*dCD_left + 0.5*dCD_cent + 0.25*dCD_right;
 447
 448   //printf("CL: %f CD:%f L/D:%f\n",CL,CD,CL/CD);
 449
 450     /* total forces in body axes */
 451   Cx_drela = -CD_drela*Cos_alpha + CL_drela*Sin_alpha*Cos_beta*Cos_beta;
 452   Cz_drela = -CD_drela*Sin_alpha - CL_drela*Cos_alpha*Cos_beta*Cos_beta;
 453   Cy_drela = CY_b*Std_Beta  + CY_p*Phat + CY_r*Rhat + CY_dr*rudder;
 454
 455     /* total moments in body axes */
 456   Cl_drela =        Cl_b*Std_Beta  + Cl_p*Phat + Cl_r_mod*Rhat + Cl_dr*rudder
 457            + dCl_stall                               + Cl_da*aileron;
 458   Cn_drela =        Cn_b*Std_Beta  + Cn_p_mod*Phat + Cn_r*Rhat + Cn_dr*rudder
 459            + dCn_stall                               + Cn_da*aileron;
 460   Cm_drela = Cm_0 + Cm_a*(Std_Alpha-Alpha_0) +dCm_stall
 461                          + Cm_q*Qhat                 + Cm_de*elevator;
 462
 463     /* set dimensional forces and moments */
 464   QS = 0.5*Density*V_rel_wind*V_rel_wind * S_ref;
 465
 466   F_X_aero = Cx_drela * QS;
 467   F_Y_aero = Cy_drela * QS;
 468   F_Z_aero = Cz_drela * QS;
 469
 470   M_l_aero = Cl_drela * QS * B_ref;
 471   M_m_aero = Cm_drela * QS * C_ref;
 472   M_n_aero = Cn_drela * QS * B_ref;
 473 //  printf("%f %f %f %f %f %f\n",F_X_aero,F_Y_aero,F_Z_aero,M_l_aero,M_m_aero,M_n_aero);
 474 }
 475