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  SCALAR G_11,G_12,G_13;
 156  SCALAR G_21,G_22,G_23;
 157  SCALAR G_31,G_32,G_33;
 158  SCALAR U_body_X,U_body_Y,U_body_Z;
 159  SCALAR V_body_X,V_body_Y,V_body_Z;
 160  SCALAR W_body_X,W_body_Y,W_body_Z;
 161  SCALAR P_atmo,Q_atmo,R_atmo;
 162
 163  // set the parameters
 164
 165  C_ref = 0.551667;
 166  B_ref = 6.55000;
 167  S_ref = 3.61111;
 168  U_ref = 19.6850;
 169  Alpha_0 = 0.349066E-01;
 170  Cm_0 = -0.112663E-01;
 171  CL_0 = 0.563172;
 172  CL_max = 1.10000;
 173  CL_min = -0.600000;
 174  CD_prof = 0.200000E-01;
 175  Uexp_CD = -0.500000;
 176  CL_a = 5.50360;
 177  Cm_a = -0.575335;
 178  CY_b = -0.415610;
 179  Cl_b = -0.250926;
 180  Cn_b = 0.567069E-01;
 181  CL_q = 7.50999;
 182  Cm_q = -11.4975;
 183  CY_p = -0.423820;
 184  Cl_p = -0.611798;
 185  Cn_p = -0.740898E-01;
 186  CY_r = 0.297540;
 187  Cl_r = 0.139581;
 188  Cn_r = -0.687755E-01;
 189  CL_de = 0.162000;
 190  Cm_de = -0.597537;
 191  CY_dr = 0.000000E+00;
 192  Cl_dr = 0.000000E+00;
 193  Cn_dr = 0.000000E+00;
 194  CY_da = -0.135890;
 195  Cl_da = -0.307921E-02;
 196  Cn_da = 0.527143E-01;
 197  span_eff = 0.95;
 198  CL_CD0 = 0.0;
 199  CD_CLsq = 0.01;
 200  CD_AIsq = 0.0;
 201  CD_ELsq = 0.0;
 202  eta_loc = 0.3;
 203  CG_arm = 0.25;
 204  CL_drop = 0.5;
 205
 206  if (!init)
 207   {
 208    init = -1;
 209   }
 210
 211  // jan's data goes -.5 to .5 while
 212  // fgfs data goes +- 1.
 213  // so I need to divide by 2 below.
 214   elevator = Long_control + Long_trim;
 215   aileron  = Lat_control;
 216   rudder   = Rudder_pedal;
 217
 218   elevator = elevator * 0.5;
 219   aileron  = aileron * 0.5;
 220   rudder   = rudder * 0.5;
 221
 222
 223
 224 //  printf("%f\n",V_rel_wind);
 225
 226   /* compute gradients of Local velocities w.r.t. Body coordinates
 227       G_11  =  dU_local/dx_body
 228       G_12  =  dU_local/dy_body   etc.  */
 229
 230   /*
 231   G_11 = U_atmo_X*T_local_to_body_11
 232        + U_atmo_Y*T_local_to_body_12
 233        + U_atmo_Z*T_local_to_body_13;
 234   G_12 = U_atmo_X*T_local_to_body_21
 235        + U_atmo_Y*T_local_to_body_22
 236        + U_atmo_Z*T_local_to_body_23;
 237   G_13 = U_atmo_X*T_local_to_body_31
 238        + U_atmo_Y*T_local_to_body_32
 239        + U_atmo_Z*T_local_to_body_33;
 240
 241   G_21 = V_atmo_X*T_local_to_body_11
 242        + V_atmo_Y*T_local_to_body_12
 243        + V_atmo_Z*T_local_to_body_13;
 244   G_22 = V_atmo_X*T_local_to_body_21
 245        + V_atmo_Y*T_local_to_body_22
 246        + V_atmo_Z*T_local_to_body_23;
 247   G_23 = V_atmo_X*T_local_to_body_31
 248        + V_atmo_Y*T_local_to_body_32
 249        + V_atmo_Z*T_local_to_body_33;
 250
 251   G_31 = W_atmo_X*T_local_to_body_11
 252        + W_atmo_Y*T_local_to_body_12
 253        + W_atmo_Z*T_local_to_body_13;
 254   G_32 = W_atmo_X*T_local_to_body_21
 255        + W_atmo_Y*T_local_to_body_22
 256        + W_atmo_Z*T_local_to_body_23;
 257   G_33 = W_atmo_X*T_local_to_body_31
 258        + W_atmo_Y*T_local_to_body_32
 259        + W_atmo_Z*T_local_to_body_33;
 260   */
 261
 262   //printf("%f %f %f %f\n",W_atmo_X,W_atmo_Y,G_31,G_32);
 263
 264   /* now compute gradients of Body velocities w.r.t. Body coordinates */
 265   /*  U_body_x  =  dU_body/dx_body   etc.  */
 266
 267   /*
 268   U_body_X = T_local_to_body_11*G_11
 269            + T_local_to_body_12*G_21
 270            + T_local_to_body_13*G_31;
 271   U_body_Y = T_local_to_body_11*G_12
 272            + T_local_to_body_12*G_22
 273            + T_local_to_body_13*G_32;
 274   U_body_Z = T_local_to_body_11*G_13
 275            + T_local_to_body_12*G_23
 276            + T_local_to_body_13*G_33;
 277
 278   V_body_X = T_local_to_body_21*G_11
 279            + T_local_to_body_22*G_21
 280            + T_local_to_body_23*G_31;
 281   V_body_Y = T_local_to_body_21*G_12
 282            + T_local_to_body_22*G_22
 283            + T_local_to_body_23*G_32;
 284   V_body_Z = T_local_to_body_21*G_13
 285            + T_local_to_body_22*G_23
 286            + T_local_to_body_23*G_33;
 287
 288   W_body_X = T_local_to_body_31*G_11
 289            + T_local_to_body_32*G_21
 290            + T_local_to_body_33*G_31;
 291   W_body_Y = T_local_to_body_31*G_12
 292            + T_local_to_body_32*G_22
 293            + T_local_to_body_33*G_32;
 294   W_body_Z = T_local_to_body_31*G_13
 295            + T_local_to_body_32*G_23
 296            + T_local_to_body_33*G_33;
 297   */
 298
 299   /* set rotation rates of airmass motion */
 300   /*   BUG
 301         P_atmo =  W_body_X;
 302         Q_atmo = -W_body_Y;
 303         R_atmo =  V_body_Z;
 304   */
 305
 306   /*
 307         P_atmo =  W_body_Y;
 308         Q_atmo = -W_body_X;
 309         R_atmo =  V_body_X;
 310   */
 311
 312         P_atmo = 0;
 313         Q_atmo = 0;
 314         R_atmo = 0;
 315
 316   //  printf("%f %f %f\n",P_atmo,Q_atmo,R_atmo);
 317   if (V_rel_wind != 0)
 318    {
 319   /* set net effective dimensionless rotation rates */
 320     Phat = (P_body - P_atmo) * B_ref / (2.0*V_rel_wind);
 321     Qhat = (Q_body - Q_atmo) * C_ref / (2.0*V_rel_wind);
 322     Rhat = (R_body - R_atmo) * B_ref / (2.0*V_rel_wind);
 323    }
 324   else
 325    {
 326     Phat=0;
 327     Qhat=0;
 328     Rhat=0;
 329    }
 330
 331 //  printf("Phat: %f Qhat: %f Rhat: %f\n",Phat,Qhat,Rhat);
 332
 333     /* compute local CL at three spanwise locations */
 334   CL_left  = CL_0 + CL_a*(Std_Alpha - Alpha_0 - Phat*eta_loc);
 335   CL_cent  = CL_0 + CL_a*(Std_Alpha - Alpha_0               );
 336   CL_right = CL_0 + CL_a*(Std_Alpha - Alpha_0 + Phat*eta_loc);
 337
 338 // printf("CL_left: %f CL_cent: %f CL_right: %f\n",CL_left,CL_cent,CL_right);
 339
 340     /* set CL-limit changes */
 341   dCL_left  = 0.;
 342   dCL_cent  = 0.;
 343   dCL_right = 0.;
 344
 345   stalling=0;
 346   if (CL_left  > CL_max)
 347     {
 348       dCL_left  = CL_max-CL_left -CL_drop;
 349       stalling=1;
 350     }
 351
 352   if (CL_cent  > CL_max)
 353     {
 354       dCL_cent  = CL_max-CL_cent -CL_drop;
 355       stalling=1;
 356     }
 357
 358   if (CL_right > CL_max)
 359     {
 360       dCL_right = CL_max-CL_right -CL_drop;
 361       stalling=1;
 362     }
 363
 364   if (CL_left  < CL_min)
 365     {
 366       dCL_left  = CL_min-CL_left -CL_drop;
 367       stalling=1;
 368     }
 369
 370   if (CL_cent  < CL_min)
 371     {
 372       dCL_cent  = CL_min-CL_cent -CL_drop;
 373       stalling=1;
 374     }
 375
 376   if (CL_right < CL_min)
 377     {
 378       dCL_right = CL_min-CL_right -CL_drop;
 379       stalling=1;
 380     }
 381
 382   /* set average wing CL */
 383   CL_wing = CL_0 + CL_a*(Std_Alpha-Alpha_0)
 384           + 0.25*dCL_left + 0.5*dCL_cent + 0.25*dCL_right;
 385
 386 //  printf("CL_wing: %f\n",CL_wing);
 387
 388
 389   /* correct profile CD for CL dependence and aileron dependence */
 390   CD_all = CD_prof
 391          + CD_CLsq*(CL_wing-CL_CD0)*(CL_wing-CL_CD0)
 392          + CD_AIsq*aileron*aileron
 393          + CD_ELsq*elevator*elevator;
 394
 395 //  printf("CD_all:%lf\n",CD_all);
 396
 397     /* scale profile CD with Reynolds number via simple power law */
 398  if (V_rel_wind > 0.1)
 399  {
 400   CD_scaled = CD_all*pow(((double)V_rel_wind/(double)U_ref),Uexp_CD);
 401  }
 402  else
 403  {
 404         CD_scaled=CD_all;
 405  }
 406
 407 //  printf("CD_scaled:%lf\n",CD_scaled);
 408
 409
 410
 411   /* Scale lateral cross-coupling derivatives with wing CL */
 412   Cl_r_mod = Cl_r*CL_wing/CL_0;
 413   Cn_p_mod = Cn_p*CL_wing/CL_0;
 414
 415   //  printf("Cl_r_mod: %f Cn_p_mod: %f\n",Cl_r_mod,Cn_p_mod);
 416
 417     /* total average CD with induced and stall contributions */
 418   dCD_left  = CD_stall*dCL_left *dCL_left ;
 419   dCD_cent  = CD_stall*dCL_cent *dCL_cent ;
 420   dCD_right = CD_stall*dCL_right*dCL_right;
 421
 422   /* total CL, with pitch rate and elevator contributions */
 423   CL_drela = (CL_wing + CL_q*Qhat + CL_de*elevator)*Cos_alpha;
 424
 425 //  printf("CL:%f\n",CL);
 426
 427     /* assymetric stall will cause roll and yaw moments */
 428   dCl =  0.45*-1*(dCL_right-dCL_left)*0.5*eta_loc;
 429   dCn =  0.45*(dCD_right-dCD_left)*0.5*eta_loc;
 430   dCm_stall = (0.25*dCL_left + 0.5*dCL_cent + 0.25*dCL_right)*CG_arm;
 431
 432 //  printf("dCl: %f dCn:%f\n",dCl,dCn);
 433
 434     /* stall-caused moments in body axes */
 435   dCl_stall = dCl*Cos_alpha - dCn*Sin_alpha;
 436   dCn_stall = dCl*Sin_alpha + dCn*Cos_alpha;
 437
 438     /* total CD, with induced and stall contributions */
 439
 440   Cl_w = Cl_b*Std_Beta  + Cl_p*Phat + Cl_r_mod*Rhat
 441        + dCl_stall  + Cl_da*aileron;
 442   CD_drela = CD_scaled
 443      + (CL*CL + 32.0*Cl_w*Cl_w)*S_ref/(B_ref*B_ref*3.14159*span_eff)
 444      + 0.25*dCD_left + 0.5*dCD_cent + 0.25*dCD_right;
 445
 446   //printf("CL: %f CD:%f L/D:%f\n",CL,CD,CL/CD);
 447
 448     /* total forces in body axes */
 449   Cx_drela = -CD_drela*Cos_alpha + CL_drela*Sin_alpha*Cos_beta*Cos_beta;
 450   Cz_drela = -CD_drela*Sin_alpha - CL_drela*Cos_alpha*Cos_beta*Cos_beta;
 451   Cy_drela = CY_b*Std_Beta  + CY_p*Phat + CY_r*Rhat + CY_dr*rudder;
 452
 453     /* total moments in body axes */
 454   Cl_drela =        Cl_b*Std_Beta  + Cl_p*Phat + Cl_r_mod*Rhat + Cl_dr*rudder
 455            + dCl_stall                               + Cl_da*aileron;
 456   Cn_drela =        Cn_b*Std_Beta  + Cn_p_mod*Phat + Cn_r*Rhat + Cn_dr*rudder
 457            + dCn_stall                               + Cn_da*aileron;
 458   Cm_drela = Cm_0 + Cm_a*(Std_Alpha-Alpha_0) +dCm_stall
 459                          + Cm_q*Qhat                 + Cm_de*elevator;
 460
 461     /* set dimensional forces and moments */
 462   QS = 0.5*Density*V_rel_wind*V_rel_wind * S_ref;
 463
 464   F_X_aero = Cx_drela * QS;
 465   F_Y_aero = Cy_drela * QS;
 466   F_Z_aero = Cz_drela * QS;
 467
 468   M_l_aero = Cl_drela * QS * B_ref;
 469   M_m_aero = Cm_drela * QS * C_ref;
 470   M_n_aero = Cn_drela * QS * B_ref;
 471 //  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);
 472 }
 473