src/FDM/JSBSim/models/FGAerodynamics.cpp

   1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   2
   3  Module:       FGAerodynamics.cpp
   4  Author:       Jon S. Berndt
   5  Date started: 09/13/00
   6  Purpose:      Encapsulates the aerodynamic forces
   7
   8  ------------- Copyright (C) 2000  Jon S. Berndt (jsb@hal-pc.org) -------------
   9
  10  This program is free software; you can redistribute it and/or modify it under
  11  the terms of the GNU Lesser General Public License as published by the Free Software
  12  Foundation; either version 2 of the License, or (at your option) any later
  13  version.
  14
  15  This program is distributed in the hope that it will be useful, but WITHOUT
  16  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  17  FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  18  details.
  19
  20  You should have received a copy of the GNU Lesser General Public License along with
  21  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  22  Place - Suite 330, Boston, MA  02111-1307, USA.
  23
  24  Further information about the GNU Lesser General Public License can also be found on
  25  the world wide web at http://www.gnu.org.
  26
  27 FUNCTIONAL DESCRIPTION
  28 --------------------------------------------------------------------------------
  29
  30 HISTORY
  31 --------------------------------------------------------------------------------
  32 09/13/00   JSB   Created
  33 04/22/01   JSB   Moved code into here from FGAircraft
  34
  35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  36 INCLUDES
  37 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  38
  39 #include "FGAerodynamics.h"
  40 #include "FGPropagate.h"
  41 #include "FGAircraft.h"
  42 #include "FGState.h"
  43 #include "FGMassBalance.h"
  44 #include <input_output/FGPropertyManager.h>
  45
  46 namespace JSBSim {
  47
  48 static const char *IdSrc = "$Id$";
  49 static const char *IdHdr = ID_AERODYNAMICS;
  50
  51 const unsigned NAxes=6;
  52 const char* AxisNames[] = { "drag", "side-force", "lift", "rolling-moment",
  53                             "pitching-moment","yawing-moment" };
  54
  55 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  56 CLASS IMPLEMENTATION
  57 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  58
  59
  60 FGAerodynamics::FGAerodynamics(FGFDMExec* FDMExec) : FGModel(FDMExec)
  61 {
  62   Name = "FGAerodynamics";
  63
  64   AxisIdx["DRAG"]  = 0;
  65   AxisIdx["SIDE"]  = 1;
  66   AxisIdx["LIFT"]  = 2;
  67   AxisIdx["ROLL"]  = 3;
  68   AxisIdx["PITCH"] = 4;
  69   AxisIdx["YAW"]   = 5;
  70
  71   Coeff = new CoeffArray[6];
  72
  73   impending_stall = stall_hyst = 0.0;
  74   alphaclmin = alphaclmax = 0.0;
  75   alphahystmin = alphahystmax = 0.0;
  76   clsq = lod = 0.0;
  77   alphaw = 0.0;
  78   bi2vel = ci2vel = 0.0;
  79   AeroRPShift = 0;
  80   vDeltaRP.InitMatrix();
  81
  82   bind();
  83
  84   Debug(0);
  85 }
  86
  87 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  88
  89 FGAerodynamics::~FGAerodynamics()
  90 {
  91   unsigned int i,j;
  92
  93   for (i=0; i<6; i++)
  94     for (j=0; j<Coeff[i].size(); j++)
  95       delete Coeff[i][j];
  96
  97   delete[] Coeff;
  98
  99   for (i=0; i<variables.size(); i++)
 100     delete variables[i];
 101
 102   if (AeroRPShift) delete AeroRPShift;
 103
 104   unbind();
 105
 106   Debug(1);
 107 }
 108
 109 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 110
 111 bool FGAerodynamics::Run(void)
 112 {
 113   unsigned int axis_ctr, ctr, i;
 114   double alpha, twovel;
 115
 116   if (FGModel::Run()) return true;
 117   if (FDMExec->Holding()) return false; // if paused don't execute
 118
 119   // calculate some oft-used quantities for speed
 120
 121   twovel = 2*Auxiliary->GetVt();
 122   if (twovel != 0) {
 123     bi2vel = Aircraft->GetWingSpan() / twovel;
 124     ci2vel = Aircraft->Getcbar() / twovel;
 125   }
 126   alphaw = Auxiliary->Getalpha() + Aircraft->GetWingIncidence();
 127   alpha = Auxiliary->Getalpha();
 128   qbar_area = Aircraft->GetWingArea() * Auxiliary->Getqbar();
 129
 130   if (alphaclmax != 0) {
 131     if (alpha > 0.85*alphaclmax) {
 132       impending_stall = 10*(alpha/alphaclmax - 0.85);
 133     } else {
 134       impending_stall = 0;
 135     }
 136   }
 137
 138   if (alphahystmax != 0.0 && alphahystmin != 0.0) {
 139     if (alpha > alphahystmax) {
 140        stall_hyst = 1;
 141     } else if (alpha < alphahystmin) {
 142        stall_hyst = 0;
 143     }
 144   }
 145
 146   vLastFs = vFs;
 147   vFs.InitMatrix();
 148
 149   // Tell the variable functions to cache their values, so while the aerodynamic
 150   // functions are being calculated for each axis, these functions do not get
 151   // calculated each time, but instead use the values that have already
 152   // been calculated for this frame.
 153   for (i=0; i<variables.size(); i++) variables[i]->cacheValue(true);
 154
 155   for (axis_ctr = 0; axis_ctr < 3; axis_ctr++) {
 156     for (ctr=0; ctr < Coeff[axis_ctr].size(); ctr++) {
 157       vFs(axis_ctr+1) += Coeff[axis_ctr][ctr]->GetValue();
 158     }
 159   }
 160
 161   // Calculate aerodynamic reference point shift, if any
 162   if (AeroRPShift) {
 163     vDeltaRP(eX) = AeroRPShift->GetValue()*Aircraft->Getcbar()*12.0;
 164   }
 165
 166   // calculate lift coefficient squared
 167   if ( Auxiliary->Getqbar() > 0) {
 168     clsq = vFs(eLift) / (Aircraft->GetWingArea()*Auxiliary->Getqbar());
 169     clsq *= clsq;
 170   }
 171
 172   if ( vFs(eDrag)  > 0) {
 173     lod = vFs(eLift) / vFs(eDrag);
 174   }
 175
 176   //correct signs of drag and lift to wind axes convention
 177   //positive forward, right, down
 178   vFs(eDrag)*=-1; vFs(eLift)*=-1;
 179
 180   // transform stability axis forces into body axes
 181   vForces = State->GetTs2b()*vFs;
 182
 183   vDXYZcg = MassBalance->StructuralToBody(Aircraft->GetXYZrp() + vDeltaRP);
 184
 185   vMoments = vDXYZcg*vForces; // M = r X F
 186
 187   for (axis_ctr = 0; axis_ctr < 3; axis_ctr++) {
 188     for (ctr = 0; ctr < Coeff[axis_ctr+3].size(); ctr++) {
 189       vMoments(axis_ctr+1) += Coeff[axis_ctr+3][ctr]->GetValue();
 190     }
 191   }
 192
 193   return false;
 194 }
 195
 196 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 197
 198 bool FGAerodynamics::Load(Element *element)
 199 {
 200   string parameter, axis, scratch;
 201   Element *temp_element, *axis_element, *function_element;
 202
 203   Debug(2);
 204
 205   if (temp_element = element->FindElement("alphalimits")) {
 206     alphaclmin = temp_element->FindElementValueAsNumberConvertTo("min", "DEG");
 207     alphaclmax = temp_element->FindElementValueAsNumberConvertTo("max", "DEG");
 208   }
 209
 210   if (temp_element = element->FindElement("hysteresis_limits")) {
 211     alphahystmin = temp_element->FindElementValueAsNumberConvertTo("min", "DEG");
 212     alphahystmax = temp_element->FindElementValueAsNumberConvertTo("max", "DEG");
 213   }
 214
 215   if (temp_element = element->FindElement("aero_ref_pt_shift_x")) {
 216     function_element = temp_element->FindElement("function");
 217     AeroRPShift = new FGFunction(PropertyManager, function_element);
 218   }
 219
 220   function_element = element->FindElement("function");
 221   while (function_element) {
 222     variables.push_back( new FGFunction(PropertyManager, function_element) );
 223     function_element = element->FindNextElement("function");
 224   }
 225
 226   axis_element = element->FindElement("axis");
 227   while (axis_element) {
 228     CoeffArray ca;
 229     axis = axis_element->GetAttributeValue("name");
 230     function_element = axis_element->FindElement("function");
 231     while (function_element) {
 232       ca.push_back( new FGFunction(PropertyManager, function_element) );
 233       function_element = axis_element->FindNextElement("function");
 234     }
 235     Coeff[AxisIdx[axis]] = ca;
 236     axis_element = element->FindNextElement("axis");
 237   }
 238
 239   return true;
 240 }
 241
 242 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 243
 244 string FGAerodynamics::GetCoefficientStrings(string delimeter)
 245 {
 246   string CoeffStrings = "";
 247   bool firstime = true;
 248   unsigned int axis, sd;
 249
 250   for (sd = 0; sd < variables.size(); sd++) {
 251     if (firstime) {
 252       firstime = false;
 253     } else {
 254       CoeffStrings += delimeter;
 255     }
 256     CoeffStrings += variables[sd]->GetName();
 257   }
 258
 259   for (axis = 0; axis < 6; axis++) {
 260     for (sd = 0; sd < Coeff[axis].size(); sd++) {
 261       if (firstime) {
 262         firstime = false;
 263       } else {
 264         CoeffStrings += delimeter;
 265       }
 266       CoeffStrings += Coeff[axis][sd]->GetName();
 267     }
 268   }
 269   return CoeffStrings;
 270 }
 271
 272 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 273
 274 string FGAerodynamics::GetCoefficientValues(string delimeter)
 275 {
 276   string SDValues = "";
 277   bool firstime = true;
 278   unsigned int sd;
 279
 280   for (sd = 0; sd < variables.size(); sd++) {
 281     if (firstime) {
 282       firstime = false;
 283     } else {
 284       SDValues += delimeter;
 285     }
 286     SDValues += variables[sd]->GetValueAsString();
 287   }
 288
 289   for (unsigned int axis = 0; axis < 6; axis++) {
 290     for (unsigned int sd = 0; sd < Coeff[axis].size(); sd++) {
 291       if (firstime) {
 292         firstime = false;
 293       } else {
 294         SDValues += delimeter;
 295       }
 296       SDValues += Coeff[axis][sd]->GetValueAsString();
 297     }
 298   }
 299
 300   return SDValues;
 301 }
 302
 303 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 304
 305 void FGAerodynamics::bind(void)
 306 {
 307   typedef double (FGAerodynamics::*PMF)(int) const;
 308
 309   PropertyManager->Tie("forces/fbx-aero-lbs", this,1,
 310                        (PMF)&FGAerodynamics::GetForces);
 311   PropertyManager->Tie("forces/fby-aero-lbs", this,2,
 312                        (PMF)&FGAerodynamics::GetForces);
 313   PropertyManager->Tie("forces/fbz-aero-lbs", this,3,
 314                        (PMF)&FGAerodynamics::GetForces);
 315   PropertyManager->Tie("moments/l-aero-lbsft", this,1,
 316                        (PMF)&FGAerodynamics::GetMoments);
 317   PropertyManager->Tie("moments/m-aero-lbsft", this,2,
 318                        (PMF)&FGAerodynamics::GetMoments);
 319   PropertyManager->Tie("moments/n-aero-lbsft", this,3,
 320                        (PMF)&FGAerodynamics::GetMoments);
 321   PropertyManager->Tie("forces/fwx-aero-lbs", this,1,
 322                        (PMF)&FGAerodynamics::GetvFs);
 323   PropertyManager->Tie("forces/fwy-aero-lbs", this,2,
 324                        (PMF)&FGAerodynamics::GetvFs);
 325   PropertyManager->Tie("forces/fwz-aero-lbs", this,3,
 326                        (PMF)&FGAerodynamics::GetvFs);
 327   PropertyManager->Tie("forces/lod-norm", this,
 328                        &FGAerodynamics::GetLoD);
 329   PropertyManager->Tie("aero/cl-squared", this,
 330                        &FGAerodynamics::GetClSquared);
 331   PropertyManager->Tie("aero/qbar-area", &qbar_area);
 332   PropertyManager->Tie("aero/alpha-max-deg", this,
 333                        &FGAerodynamics::GetAlphaCLMax,
 334                        &FGAerodynamics::SetAlphaCLMax,
 335                        true);
 336   PropertyManager->Tie("aero/alpha-min-deg", this,
 337                        &FGAerodynamics::GetAlphaCLMin,
 338                        &FGAerodynamics::SetAlphaCLMin,
 339                        true);
 340   PropertyManager->Tie("aero/bi2vel", this,
 341                        &FGAerodynamics::GetBI2Vel);
 342   PropertyManager->Tie("aero/ci2vel", this,
 343                        &FGAerodynamics::GetCI2Vel);
 344   PropertyManager->Tie("aero/alpha-wing-rad", this,
 345                        &FGAerodynamics::GetAlphaW);
 346   PropertyManager->Tie("systems/stall-warn-norm", this,
 347                         &FGAerodynamics::GetStallWarn);
 348   PropertyManager->Tie("aero/stall-hyst-norm", this,
 349                         &FGAerodynamics::GetHysteresisParm);
 350 }
 351
 352 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 353
 354 void FGAerodynamics::unbind(void)
 355 {
 356   PropertyManager->Untie("forces/fbx-aero-lbs");
 357   PropertyManager->Untie("forces/fby-aero-lbs");
 358   PropertyManager->Untie("forces/fbz-aero-lbs");
 359   PropertyManager->Untie("moments/l-aero-lbsft");
 360   PropertyManager->Untie("moments/m-aero-lbsft");
 361   PropertyManager->Untie("moments/n-aero-lbsft");
 362   PropertyManager->Untie("forces/fwx-aero-lbs");
 363   PropertyManager->Untie("forces/fwy-aero-lbs");
 364   PropertyManager->Untie("forces/fwz-aero-lbs");
 365   PropertyManager->Untie("forces/lod-norm");
 366   PropertyManager->Untie("aero/cl-squared");
 367   PropertyManager->Untie("aero/qbar-area");
 368   PropertyManager->Untie("aero/alpha-max-deg");
 369   PropertyManager->Untie("aero/alpha-min-deg");
 370   PropertyManager->Untie("aero/bi2vel");
 371   PropertyManager->Untie("aero/ci2vel");
 372   PropertyManager->Untie("aero/alpha-wing-rad");
 373   PropertyManager->Untie("aero/stall-hyst-norm");
 374   PropertyManager->Untie("systems/stall-warn-norm");
 375 }
 376
 377 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 378 //    The bitmasked value choices are as follows:
 379 //    unset: In this case (the default) JSBSim would only print
 380 //       out the normally expected messages, essentially echoing
 381 //       the config files as they are read. If the environment
 382 //       variable is not set, debug_lvl is set to 1 internally
 383 //    0: This requests JSBSim not to output any messages
 384 //       whatsoever.
 385 //    1: This value explicity requests the normal JSBSim
 386 //       startup messages
 387 //    2: This value asks for a message to be printed out when
 388 //       a class is instantiated
 389 //    4: When this value is set, a message is displayed when a
 390 //       FGModel object executes its Run() method
 391 //    8: When this value is set, various runtime state variables
 392 //       are printed out periodically
 393 //    16: When set various parameters are sanity checked and
 394 //       a message is printed out when they go out of bounds
 395
 396 void FGAerodynamics::Debug(int from)
 397 {
 398   if (debug_lvl <= 0) return;
 399
 400   if (debug_lvl & 1) { // Standard console startup message output
 401     if (from == 2) { // Loader
 402       cout << endl << "  Aerodynamics: " << endl;
 403     }
 404   }
 405   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 406     if (from == 0) cout << "Instantiated: FGAerodynamics" << endl;
 407     if (from == 1) cout << "Destroyed:    FGAerodynamics" << endl;
 408   }
 409   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 410   }
 411   if (debug_lvl & 8 ) { // Runtime state variables
 412   }
 413   if (debug_lvl & 16) { // Sanity checking
 414   }
 415   if (debug_lvl & 64) {
 416     if (from == 0) { // Constructor
 417       cout << IdSrc << endl;
 418       cout << IdHdr << endl;
 419     }
 420   }
 421 }
 422
 423 } // namespace JSBSim