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 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 General Public License for more
  18  details.
  19
  20  You should have received a copy of the GNU 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 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   bind();
  80
  81   Debug(0);
  82 }
  83
  84 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  85
  86 FGAerodynamics::~FGAerodynamics()
  87 {
  88   unsigned int i,j;
  89
  90   unbind();
  91
  92   for (i=0; i<6; i++)
  93     for (j=0; j<Coeff[i].size(); j++)
  94       delete Coeff[i][j];
  95
  96   delete[] Coeff;
  97
  98   for (i=0; i<variables.size(); i++)
  99     delete variables[i];
 100
 101   Debug(1);
 102 }
 103
 104 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 105
 106 bool FGAerodynamics::Run(void)
 107 {
 108   unsigned int axis_ctr, ctr, i;
 109   double alpha, twovel;
 110
 111   if (FGModel::Run()) return true;
 112   if (FDMExec->Holding()) return false; // if paused don't execute
 113
 114   // calculate some oft-used quantities for speed
 115
 116   twovel = 2*Auxiliary->GetVt();
 117   if (twovel != 0) {
 118     bi2vel = Aircraft->GetWingSpan() / twovel;
 119     ci2vel = Aircraft->Getcbar() / twovel;
 120   }
 121   alphaw = Auxiliary->Getalpha() + Aircraft->GetWingIncidence();
 122   alpha = Auxiliary->Getalpha();
 123   qbar_area = Aircraft->GetWingArea() * Auxiliary->Getqbar();
 124
 125   if (alphaclmax != 0) {
 126     if (alpha > 0.85*alphaclmax) {
 127       impending_stall = 10*(alpha/alphaclmax - 0.85);
 128     } else {
 129       impending_stall = 0;
 130     }
 131   }
 132
 133   if (alphahystmax != 0.0 && alphahystmin != 0.0) {
 134     if (alpha > alphahystmax) {
 135        stall_hyst = 1;
 136     } else if (alpha < alphahystmin) {
 137        stall_hyst = 0;
 138     }
 139   }
 140
 141   vLastFs = vFs;
 142   vFs.InitMatrix();
 143
 144   // Tell the variable functions to cache their values, so while the aerodynamic
 145   // functions are being calculated for each axis, these functions do not get
 146   // calculated each time, but instead use the values that have already
 147   // been calculated for this frame.
 148   for (i=0; i<variables.size(); i++) variables[i]->cacheValue(true);
 149
 150   for (axis_ctr = 0; axis_ctr < 3; axis_ctr++) {
 151     for (ctr=0; ctr < Coeff[axis_ctr].size(); ctr++) {
 152       vFs(axis_ctr+1) += Coeff[axis_ctr][ctr]->GetValue();
 153     }
 154   }
 155
 156   // calculate lift coefficient squared
 157   if ( Auxiliary->Getqbar() > 0) {
 158     clsq = vFs(eLift) / (Aircraft->GetWingArea()*Auxiliary->Getqbar());
 159     clsq *= clsq;
 160   }
 161
 162   if ( vFs(eDrag)  > 0) {
 163     lod = vFs(eLift) / vFs(eDrag);
 164   }
 165
 166   //correct signs of drag and lift to wind axes convention
 167   //positive forward, right, down
 168   vFs(eDrag)*=-1; vFs(eLift)*=-1;
 169
 170   // transform stability axis forces into body axes
 171   vForces = State->GetTs2b()*vFs;
 172
 173   vDXYZcg = MassBalance->StructuralToBody(Aircraft->GetXYZrp());
 174
 175   vMoments = vDXYZcg*vForces; // M = r X F
 176
 177   for (axis_ctr = 0; axis_ctr < 3; axis_ctr++) {
 178     for (ctr = 0; ctr < Coeff[axis_ctr+3].size(); ctr++) {
 179       vMoments(axis_ctr+1) += Coeff[axis_ctr+3][ctr]->GetValue();
 180     }
 181   }
 182
 183   return false;
 184 }
 185
 186 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 187
 188 bool FGAerodynamics::Load(Element *element)
 189 {
 190   string parameter, axis, scratch;
 191   Element *temp_element, *axis_element, *function_element;
 192
 193   Debug(2);
 194
 195   if (temp_element = element->FindElement("alphalimits")) {
 196     alphaclmin = temp_element->FindElementValueAsNumberConvertTo("min", "DEG");
 197     alphaclmax = temp_element->FindElementValueAsNumberConvertTo("max", "DEG");
 198   }
 199
 200   if (temp_element = element->FindElement("hysteresis_limits")) {
 201     alphahystmin = temp_element->FindElementValueAsNumberConvertTo("min", "DEG");
 202     alphahystmax = temp_element->FindElementValueAsNumberConvertTo("max", "DEG");
 203   }
 204
 205   function_element = element->FindElement("function");
 206   while (function_element) {
 207     variables.push_back( new FGFunction(PropertyManager, function_element) );
 208     function_element = element->FindNextElement("function");
 209   }
 210
 211   axis_element = element->FindElement("axis");
 212   while (axis_element) {
 213     CoeffArray ca;
 214     axis = axis_element->GetAttributeValue("name");
 215     function_element = axis_element->FindElement("function");
 216     while (function_element) {
 217       ca.push_back( new FGFunction(PropertyManager, function_element) );
 218       function_element = axis_element->FindNextElement("function");
 219     }
 220     Coeff[AxisIdx[axis]] = ca;
 221     axis_element = element->FindNextElement("axis");
 222   }
 223
 224   return true;
 225 }
 226
 227 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 228
 229 string FGAerodynamics::GetCoefficientStrings(string delimeter)
 230 {
 231   string CoeffStrings = "";
 232   bool firstime = true;
 233   unsigned int axis, sd;
 234
 235   for (sd = 0; sd < variables.size(); sd++) {
 236     if (firstime) {
 237       firstime = false;
 238     } else {
 239       CoeffStrings += delimeter;
 240     }
 241     CoeffStrings += variables[sd]->GetName();
 242   }
 243
 244   for (axis = 0; axis < 6; axis++) {
 245     for (sd = 0; sd < Coeff[axis].size(); sd++) {
 246       CoeffStrings += delimeter;
 247       CoeffStrings += Coeff[axis][sd]->GetName();
 248     }
 249   }
 250   return CoeffStrings;
 251 }
 252
 253 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 254
 255 string FGAerodynamics::GetCoefficientValues(string delimeter)
 256 {
 257   string SDValues = "";
 258   bool firstime = true;
 259   unsigned int sd;
 260
 261   for (sd = 0; sd < variables.size(); sd++) {
 262     if (firstime) {
 263       firstime = false;
 264     } else {
 265       SDValues += delimeter;
 266     }
 267     SDValues += variables[sd]->GetValueAsString();
 268   }
 269
 270   for (unsigned int axis = 0; axis < 6; axis++) {
 271     for (unsigned int sd = 0; sd < Coeff[axis].size(); sd++) {
 272       SDValues += delimeter;
 273       SDValues += Coeff[axis][sd]->GetValueAsString();
 274     }
 275   }
 276
 277   return SDValues;
 278 }
 279
 280 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 281
 282 void FGAerodynamics::bind(void)
 283 {
 284   typedef double (FGAerodynamics::*PMF)(int) const;
 285
 286   PropertyManager->Tie("forces/fbx-aero-lbs", this,1,
 287                        (PMF)&FGAerodynamics::GetForces);
 288   PropertyManager->Tie("forces/fby-aero-lbs", this,2,
 289                        (PMF)&FGAerodynamics::GetForces);
 290   PropertyManager->Tie("forces/fbz-aero-lbs", this,3,
 291                        (PMF)&FGAerodynamics::GetForces);
 292   PropertyManager->Tie("moments/l-aero-lbsft", this,1,
 293                        (PMF)&FGAerodynamics::GetMoments);
 294   PropertyManager->Tie("moments/m-aero-lbsft", this,2,
 295                        (PMF)&FGAerodynamics::GetMoments);
 296   PropertyManager->Tie("moments/n-aero-lbsft", this,3,
 297                        (PMF)&FGAerodynamics::GetMoments);
 298   PropertyManager->Tie("forces/fwx-aero-lbs", this,1,
 299                        (PMF)&FGAerodynamics::GetvFs);
 300   PropertyManager->Tie("forces/fwy-aero-lbs", this,2,
 301                        (PMF)&FGAerodynamics::GetvFs);
 302   PropertyManager->Tie("forces/fwz-aero-lbs", this,3,
 303                        (PMF)&FGAerodynamics::GetvFs);
 304   PropertyManager->Tie("forces/lod-norm", this,
 305                        &FGAerodynamics::GetLoD);
 306   PropertyManager->Tie("aero/cl-squared", this,
 307                        &FGAerodynamics::GetClSquared);
 308   PropertyManager->Tie("aero/qbar-area", &qbar_area);
 309   PropertyManager->Tie("aero/alpha-max-deg", this,
 310                        &FGAerodynamics::GetAlphaCLMax,
 311                        &FGAerodynamics::SetAlphaCLMax,
 312                        true);
 313   PropertyManager->Tie("aero/alpha-min-deg", this,
 314                        &FGAerodynamics::GetAlphaCLMin,
 315                        &FGAerodynamics::SetAlphaCLMin,
 316                        true);
 317   PropertyManager->Tie("aero/bi2vel", this,
 318                        &FGAerodynamics::GetBI2Vel);
 319   PropertyManager->Tie("aero/ci2vel", this,
 320                        &FGAerodynamics::GetCI2Vel);
 321   PropertyManager->Tie("aero/alpha-wing-rad", this,
 322                        &FGAerodynamics::GetAlphaW);
 323   PropertyManager->Tie("systems/stall-warn-norm", this,
 324                         &FGAerodynamics::GetStallWarn);
 325   PropertyManager->Tie("aero/stall-hyst-norm", this,
 326                         &FGAerodynamics::GetHysteresisParm);
 327 }
 328
 329 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 330
 331 void FGAerodynamics::unbind(void)
 332 {
 333   unsigned i,j;
 334
 335   PropertyManager->Untie("forces/fbx-aero-lbs");
 336   PropertyManager->Untie("forces/fby-aero-lbs");
 337   PropertyManager->Untie("forces/fbz-aero-lbs");
 338   PropertyManager->Untie("moments/l-aero-lbsft");
 339   PropertyManager->Untie("moments/m-aero-lbsft");
 340   PropertyManager->Untie("moments/n-aero-lbsft");
 341   PropertyManager->Untie("forces/fwx-aero-lbs");
 342   PropertyManager->Untie("forces/fwy-aero-lbs");
 343   PropertyManager->Untie("forces/fwz-aero-lbs");
 344   PropertyManager->Untie("forces/lod-norm");
 345   PropertyManager->Untie("aero/cl-squared");
 346   PropertyManager->Untie("aero/qbar-area");
 347   PropertyManager->Untie("aero/alpha-max-deg");
 348   PropertyManager->Untie("aero/alpha-min-deg");
 349   PropertyManager->Untie("aero/bi2vel");
 350   PropertyManager->Untie("aero/ci2vel");
 351   PropertyManager->Untie("aero/alpha-wing-rad");
 352   PropertyManager->Untie("aero/stall-hyst-norm");
 353   PropertyManager->Untie("systems/stall-warn-norm");
 354
 355 }
 356
 357 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 358 //    The bitmasked value choices are as follows:
 359 //    unset: In this case (the default) JSBSim would only print
 360 //       out the normally expected messages, essentially echoing
 361 //       the config files as they are read. If the environment
 362 //       variable is not set, debug_lvl is set to 1 internally
 363 //    0: This requests JSBSim not to output any messages
 364 //       whatsoever.
 365 //    1: This value explicity requests the normal JSBSim
 366 //       startup messages
 367 //    2: This value asks for a message to be printed out when
 368 //       a class is instantiated
 369 //    4: When this value is set, a message is displayed when a
 370 //       FGModel object executes its Run() method
 371 //    8: When this value is set, various runtime state variables
 372 //       are printed out periodically
 373 //    16: When set various parameters are sanity checked and
 374 //       a message is printed out when they go out of bounds
 375
 376 void FGAerodynamics::Debug(int from)
 377 {
 378   if (debug_lvl <= 0) return;
 379
 380   if (debug_lvl & 1) { // Standard console startup message output
 381     if (from == 2) { // Loader
 382       cout << endl << "  Aerodynamics: " << endl;
 383     }
 384   }
 385   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 386     if (from == 0) cout << "Instantiated: FGAerodynamics" << endl;
 387     if (from == 1) cout << "Destroyed:    FGAerodynamics" << endl;
 388   }
 389   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 390   }
 391   if (debug_lvl & 8 ) { // Runtime state variables
 392   }
 393   if (debug_lvl & 16) { // Sanity checking
 394   }
 395   if (debug_lvl & 64) {
 396     if (from == 0) { // Constructor
 397       cout << IdSrc << endl;
 398       cout << IdHdr << endl;
 399     }
 400   }
 401 }
 402
 403 } // namespace JSBSim