src/FDM/JSBSim/models/propulsion/FGRotor.cpp

   1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   2
   3  Module:       FGRotor.cpp
   4  Author:       Jon S. Berndt
   5  Date started: 08/24/00
   6  Purpose:      Encapsulates the rotor object
   7
   8  ------------- Copyright (C) 2000  Jon S. Berndt (jon@jsbsim.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 08/24/00  JSB  Created
  33 01/01/10  T.Kreitler test implementation
  34 11/15/10  T.Kreitler treated flow solver bug, flow and torque calculations
  35                      simplified, tiploss influence removed from flapping angles
  36 01/10/11  T.Kreitler changed to single rotor model
  37
  38 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  39 INCLUDES
  40 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  41
  42 #include <iostream>
  43 #include <sstream>
  44
  45 #include "FGRotor.h"
  46
  47 #include "models/FGPropagate.h"
  48 #include "models/FGAtmosphere.h"
  49 #include "models/FGAuxiliary.h"
  50 #include "models/FGMassBalance.h"
  51
  52 #include "input_output/FGXMLElement.h"
  53
  54
  55 using namespace std;
  56
  57 namespace JSBSim {
  58
  59 static const char *IdSrc = "$Id: FGRotor.cpp,v 1.11 2011/01/17 22:09:59 jberndt Exp $";
  60 static const char *IdHdr = ID_ROTOR;
  61
  62 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  63 MISC
  64 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  65
  66 static inline double sqr(double x) { return x*x; }
  67
  68 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  69 CLASS IMPLEMENTATION
  70 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  71
  72
  73 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  74
  75 // Constructor
  76
  77 FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
  78                     : FGThruster(exec, rotor_element, num),
  79
  80
  81   // environment
  82   dt(0.0), rho(0.002356),
  83
  84   // configuration parameters
  85   Radius(0.0), BladeNum(0),
  86
  87   Sense(1.0), NominalRPM(0.0), ExternalRPM(0), RPMdefinition(0), ExtRPMsource(NULL),
  88
  89   BladeChord(0.0), LiftCurveSlope(0.0), BladeTwist(0.0), HingeOffset(0.0),
  90   BladeFlappingMoment(0.0), BladeMassMoment(0.0), PolarMoment(0.0),
  91   InflowLag(0.0),
  92   TipLossB(0.0),
  93
  94   GroundEffectExp(0.0), GroundEffectShift(0.0),
  95
  96   // derived parameters
  97   LockNumberByRho(0.0), Solidity(0.0),
  98
  99   // dynamic values
 100   RPM(0.0), Omega(0.0),
 101
 102   beta_orient(0.0),
 103   a0(0.0), a_1(0.0), b_1(0.0), a_dw(0.0), a1s(0.0), b1s(0.0),
 104
 105   H_drag(0.0), J_side(0.0), Torque(0.0), C_T(0.0),
 106
 107   lambda(-0.001), mu(0.0), nu(0.001), v_induced(0.0),
 108   theta_downwash(0.0), phi_downwash(0.0),
 109
 110   // control
 111   ControlMap(eMainCtrl),
 112   CollectiveCtrl(0.0), LateralCtrl(0.0), LongitudinalCtrl(0.0)
 113
 114 {
 115   FGColumnVector3 location(0.0, 0.0, 0.0), orientation(0.0, 0.0, 0.0);
 116   Element *thruster_element;
 117
 118   // initialise/set remaining variables
 119   SetTransformType(FGForce::tCustom);
 120   PropertyManager = exec->GetPropertyManager();
 121   Type = ttRotor;
 122   GearRatio = 1.0;
 123
 124   dt = exec->GetDeltaT();
 125   for (int i=0; i<5; i++) R[i] = 0.0;
 126   for (int i=0; i<5; i++) B[i] = 0.0;
 127
 128   // get positions
 129   thruster_element = rotor_element->GetParent()->FindElement("sense");
 130   if (thruster_element) {
 131     double s = thruster_element->GetDataAsNumber();
 132     if (s < -0.1) {
 133       Sense = -1.0; // 'CW' as seen from above
 134     } else if (s < 0.1) {
 135       Sense = 0.0;  // 'coaxial'
 136     } else {
 137       Sense = 1.0; // 'CCW' as seen from above
 138     }
 139   }
 140
 141   thruster_element = rotor_element->GetParent()->FindElement("location");
 142   if (thruster_element) {
 143     location = thruster_element->FindElementTripletConvertTo("IN");
 144   } else {
 145     cerr << "No thruster location found." << endl;
 146   }
 147
 148   thruster_element = rotor_element->GetParent()->FindElement("orient");
 149   if (thruster_element) {
 150     orientation = thruster_element->FindElementTripletConvertTo("RAD");
 151   } else {
 152     cerr << "No thruster orientation found." << endl;
 153   }
 154
 155   SetLocation(location);
 156   SetAnglesToBody(orientation);
 157   InvTransform = Transform().Transposed();
 158
 159   // wire controls
 160   ControlMap = eMainCtrl;
 161   if (rotor_element->FindElement("controlmap")) {
 162     string cm = rotor_element->FindElementValue("controlmap");
 163     cm = to_upper(cm);
 164     if (cm == "TAIL") {
 165       ControlMap = eTailCtrl;
 166     } else if (cm == "TANDEM") {
 167       ControlMap = eTandemCtrl;
 168     } else {
 169       cerr << "# found unknown controlmap: '" << cm << "' using main rotor config."  << endl;
 170     }
 171   }
 172
 173   // ExternalRPM -- is the RPM dictated ?
 174   if (rotor_element->FindElement("ExternalRPM")) {
 175     ExternalRPM = 1;
 176     RPMdefinition = (int) rotor_element->FindElementValueAsNumber("ExternalRPM");
 177   }
 178
 179   // configure the rotor parameters
 180   Configure(rotor_element);
 181
 182   // shaft representation - a rather simple transform,
 183   // but using a matrix is safer.
 184   TboToHsr.InitMatrix(   0.0, 0.0, 1.0,
 185                          0.0, 1.0, 0.0,
 186                         -1.0, 0.0, 0.0  );
 187   HsrToTbo  =  TboToHsr.Transposed();
 188
 189   // smooth out jumps in hagl reported, otherwise the ground effect
 190   // calculation would cause jumps too. 1Hz seems sufficient.
 191   damp_hagl = Filter(1.0,dt);
 192
 193   // enable import-export
 194   BindModel();
 195
 196   Debug(0);
 197
 198 }  // Constructor
 199
 200 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 201
 202 FGRotor::~FGRotor(){
 203   Debug(1);
 204 }
 205
 206
 207 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 208
 209 // 5in1: value-fetch-convert-default-return function
 210
 211 double FGRotor::ConfigValueConv( Element* el, const string& ename, double default_val,
 212                                   const string& unit, bool tell)
 213 {
 214
 215   Element *e = NULL;
 216   double val = default_val;
 217
 218   string pname = "*No parent element*";
 219
 220   if (el) {
 221     e = el->FindElement(ename);
 222     pname = el->GetName();
 223   }
 224
 225   if (e) {
 226     if (unit.empty()) {
 227       val = e->GetDataAsNumber();
 228     } else {
 229       val = el->FindElementValueAsNumberConvertTo(ename,unit);
 230     }
 231   } else {
 232     if (tell) {
 233       cerr << pname << ": missing element '" << ename <<
 234                        "' using estimated value: " << default_val << endl;
 235     }
 236   }
 237
 238   return val;
 239 }
 240
 241 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 242
 243 double FGRotor::ConfigValue(Element* el, const string& ename, double default_val, bool tell)
 244 {
 245   return ConfigValueConv(el, ename, default_val, "", tell);
 246 }
 247
 248 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 249
 250 // 1. read configuration and try to fill holes, ymmv
 251 // 2. calculate derived parameters and transforms
 252 void FGRotor::Configure(Element* rotor_element)
 253 {
 254
 255   double estimate;
 256   const bool yell   = true;
 257   const bool silent = false;
 258
 259
 260   Radius = 0.5 * ConfigValueConv(rotor_element, "diameter", 42.0, "FT", yell);
 261   Radius = Constrain(1e-3, Radius, 1e9);
 262
 263   BladeNum = (int) ConfigValue(rotor_element, "numblades", 3 , yell);
 264
 265   GearRatio = ConfigValue(rotor_element, "gearratio", 1.0, yell);
 266
 267   // make sure that v_tip (omega*r) is below 0.7mach ~ 750ft/s
 268   estimate = (750.0/Radius)/(2.0*M_PI) * 60.0;  // 7160/Radius
 269   NominalRPM = ConfigValue(rotor_element, "nominalrpm", estimate, yell);
 270
 271   estimate = Constrain(0.07, 2.0/Radius , 0.14); // guess solidity
 272   estimate = estimate * M_PI*Radius/BladeNum;
 273   BladeChord = ConfigValueConv(rotor_element, "chord", estimate, "FT", yell);
 274
 275   LiftCurveSlope = ConfigValue(rotor_element, "liftcurveslope", 6.0); // "1/RAD"
 276   BladeTwist = ConfigValueConv(rotor_element, "twist", -0.17, "RAD");
 277
 278   HingeOffset = ConfigValueConv(rotor_element, "hingeoffset", 0.05 * Radius, "FT" );
 279
 280   estimate = sqr(BladeChord) * sqr(Radius - HingeOffset) * 0.57;
 281   BladeFlappingMoment = ConfigValueConv(rotor_element, "flappingmoment", estimate, "SLUG*FT2");
 282   BladeFlappingMoment = Constrain(0.001, BladeFlappingMoment, 1e9);
 283
 284   // guess mass from moment of a thin stick, and multiply by the blades cg distance
 285   estimate = ( 3.0 * BladeFlappingMoment / sqr(Radius) ) * (0.45 * Radius) ;
 286   BladeMassMoment = ConfigValue(rotor_element, "massmoment", estimate); // unit is slug-ft
 287   BladeMassMoment = Constrain(0.001, BladeMassMoment, 1e9);
 288
 289   TipLossB = ConfigValue(rotor_element, "tiplossfactor", 1.0, silent);
 290
 291   estimate = 1.1 * BladeFlappingMoment * BladeNum;
 292   PolarMoment = ConfigValueConv(rotor_element, "polarmoment", estimate, "SLUG*FT2");
 293   PolarMoment = Constrain(0.001, PolarMoment, 1e9);
 294
 295   InflowLag = ConfigValue(rotor_element, "inflowlag", 0.2, yell); // fixme, depends on size
 296   InflowLag = Constrain(1e-6, InflowLag, 2.0);
 297
 298
 299   // ground effect
 300   if (rotor_element->FindElement("cgroundeffect")) {
 301     double cge,gee;
 302     cge = rotor_element->FindElementValueAsNumber("cgroundeffect");
 303     cge = Constrain(1e-9, cge, 1.0);
 304     gee = 1.0 / ( 2.0*Radius * cge );
 305     cerr << "# *** 'cgroundeffect' is defunct." << endl;
 306     cerr << "# *** use 'groundeffectexp' with: " << gee << endl;
 307   }
 308
 309   GroundEffectExp = ConfigValue(rotor_element, "groundeffectexp", 0.0);
 310   GroundEffectShift = ConfigValueConv(rotor_element, "groundeffectshift", 0.0, "FT");
 311
 312   // precalc often used powers
 313   R[0]=1.0; R[1]=Radius;   R[2]=R[1]*R[1]; R[3]=R[2]*R[1]; R[4]=R[3]*R[1];
 314   B[0]=1.0; B[1]=TipLossB; B[2]=B[1]*B[1]; B[3]=B[2]*B[1]; B[4]=B[3]*B[1];
 315
 316   // derived parameters
 317   LockNumberByRho = LiftCurveSlope * BladeChord * R[4] / BladeFlappingMoment;
 318   Solidity = BladeNum * BladeChord / (M_PI * Radius);
 319
 320   return;
 321 } // Configure
 322
 323 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 324
 325 // calculate control-axes components of total airspeed at the hub.
 326 // sets rotor orientation angle (beta) as side effect. /SH79/ eqn(19-22)
 327
 328 FGColumnVector3 FGRotor::hub_vel_body2ca( const FGColumnVector3 &uvw,
 329                                                  const FGColumnVector3 &pqr,
 330                                                  double a_ic, double b_ic)
 331 {
 332   FGColumnVector3  v_r, v_shaft, v_w;
 333   FGColumnVector3 pos;
 334
 335   pos = fdmex->GetMassBalance()->StructuralToBody(GetActingLocation());
 336
 337   v_r = uvw + pqr*pos;
 338   v_shaft = TboToHsr * InvTransform * v_r;
 339
 340   beta_orient = atan2(v_shaft(eV),v_shaft(eU));
 341
 342   v_w(eU) = v_shaft(eU)*cos(beta_orient) + v_shaft(eV)*sin(beta_orient);
 343   v_w(eV) = 0.0;
 344   v_w(eW) = v_shaft(eW) - b_ic*v_shaft(eU) - a_ic*v_shaft(eV);
 345
 346   return v_w;
 347 }
 348
 349 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 350
 351 // express fuselage angular velocity in control axes /SH79/ eqn(30,31)
 352
 353 FGColumnVector3 FGRotor::fus_angvel_body2ca( const FGColumnVector3 &pqr)
 354 {
 355   FGColumnVector3 av_s_fus, av_w_fus;
 356
 357   // for comparison:
 358   // av_s_fus = BodyToShaft * pqr; /SH79/
 359   // BodyToShaft = TboToHsr * InvTransform
 360   av_s_fus = TboToHsr * InvTransform * pqr;
 361
 362   av_w_fus(eP)=   av_s_fus(eP)*cos(beta_orient) + av_s_fus(eQ)*sin(beta_orient);
 363   av_w_fus(eQ)= - av_s_fus(eP)*sin(beta_orient) + av_s_fus(eQ)*cos(beta_orient);
 364   av_w_fus(eR)=   av_s_fus(eR);
 365
 366   return av_w_fus;
 367 }
 368
 369
 370 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 371
 372 // The calculation is a bit tricky because thrust depends on induced velocity,
 373 // and vice versa.
 374 //
 375 // The flow_scale parameter (ranging from 0.5-1.0) is used to approximate a
 376 // reduction of inflow if the helicopter is close to the ground, yielding to
 377 // higher thrust, see /TA77/ eqn(10a).
 378
 379 void FGRotor::calc_flow_and_thrust( double theta_0, double Uw, double Ww,
 380                                     double flow_scale)
 381 {
 382
 383   double ct_over_sigma = 0.0;
 384   double c0, ct_l, ct_t0, ct_t1;
 385   double mu2;
 386
 387   mu = Uw/(Omega*Radius); // /SH79/ eqn(24)
 388   mu2 = sqr(mu);
 389
 390   ct_t0 = (1.0/3.0*B[3] + 1.0/2.0 * TipLossB*mu2 - 4.0/(9.0*M_PI) * mu*mu2 ) * theta_0;
 391   ct_t1 = (1.0/4.0*B[4] + 1.0/4.0 * B[2]*mu2) * BladeTwist;
 392
 393   ct_l  = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda; // first time
 394
 395   c0 = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1) * Solidity;
 396   c0 = c0 / ( 2.0 * sqrt( sqr(mu) + sqr(lambda) ) + 1e-15);
 397
 398   // replacement for /SH79/ eqn(26).
 399   // ref: dnu/dt = 1/tau ( Ct / (2*sqrt(mu^2+lambda^2))  -  nu )
 400   // taking mu and lambda constant, this integrates to
 401
 402   nu  = flow_scale * ((nu - c0) * exp(-dt/InflowLag) + c0);
 403
 404   // now from nu to lambda, C_T, and Thrust
 405
 406   lambda = Ww/(Omega*Radius) - nu; // /SH79/ eqn(25)
 407
 408   ct_l  = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda;
 409
 410   ct_over_sigma = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1); // /SH79/ eqn(27)
 411
 412   Thrust = BladeNum*BladeChord*Radius*rho*sqr(Omega*Radius) * ct_over_sigma;
 413
 414   C_T = ct_over_sigma * Solidity;
 415   v_induced = nu * (Omega*Radius);
 416
 417 }
 418
 419
 420 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 421
 422 // The coning angle doesn't apply for teetering rotors, but calculating
 423 // doesn't hurt. /SH79/ eqn(29)
 424
 425 void FGRotor::calc_coning_angle(double theta_0)
 426 {
 427   double lock_gamma = LockNumberByRho * rho;
 428
 429   double a0_l  = (1.0/6.0  + 0.04 * mu*mu*mu) * lambda;
 430   double a0_t0 = (1.0/8.0  + 1.0/8.0  * mu*mu) * theta_0;
 431   double a0_t1 = (1.0/10.0 + 1.0/12.0 * mu*mu) * BladeTwist;
 432   a0 = lock_gamma * ( a0_l + a0_t0 + a0_t1);
 433   return;
 434 }
 435
 436 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 437
 438 // Flapping angles relative to control axes /SH79/ eqn(32)
 439
 440 void FGRotor::calc_flapping_angles(double theta_0, const FGColumnVector3 &pqr_fus_w)
 441 {
 442   double lock_gamma = LockNumberByRho * rho;
 443
 444
 445   double mu2_2 = sqr(mu)/2.0;
 446   double t075 = theta_0 + 0.75 * BladeTwist;  // common approximation for rectangular blades
 447
 448   a_1 = 1.0/(1.0 - mu2_2) * (
 449                                  (2.0*lambda + (8.0/3.0)*t075)*mu
 450                                + pqr_fus_w(eP)/Omega
 451                                - 16.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
 452                              );
 453
 454   b_1 = 1.0/(1.0 + mu2_2) * (
 455                                  (4.0/3.0)*mu*a0
 456                                - pqr_fus_w(eQ)/Omega
 457                                - 16.0 * pqr_fus_w(eP)/(lock_gamma*Omega)
 458                              );
 459
 460   // used in  force calc
 461   a_dw = 1.0/(1.0 - mu2_2) * (
 462                                  (2.0*lambda + (8.0/3.0)*t075)*mu
 463                                - 24.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
 464                                  * ( 1.0 - ( 0.29 * t075 / (C_T/Solidity) ) )
 465                              );
 466
 467   return;
 468 }
 469
 470 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 471
 472 // /SH79/ eqn(38,39)
 473
 474 void FGRotor::calc_drag_and_side_forces(double theta_0)
 475 {
 476   double cy_over_sigma  ;
 477   double t075 = theta_0 + 0.75 * BladeTwist;
 478
 479   H_drag = Thrust * a_dw;
 480
 481   cy_over_sigma = (
 482                       0.75*b_1*lambda - 1.5*a0*mu*lambda + 0.25*a_1*b_1*mu
 483                     - a0*a_1*sqr(mu) + (1.0/6.0)*a0*a_1
 484                     - (0.75*mu*a0 - (1.0/3.0)*b_1 - 0.5*sqr(mu)*b_1)*t075
 485                   );
 486   cy_over_sigma *= LiftCurveSlope/2.0;
 487
 488   J_side = BladeNum * BladeChord * Radius * rho * sqr(Omega*Radius) * cy_over_sigma;
 489
 490   return;
 491 }
 492
 493 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 494
 495 // Simplified version of /SH79/ eqn(36). Uses an estimate for blade drag
 496 // (a new config parameter to come...).
 497 // From "Bramwell's Helicopter Dynamics"  second edition, eqn(3.43) and (3.44)
 498
 499 void FGRotor::calc_torque(double theta_0)
 500 {
 501   // estimate blade drag
 502   double delta_dr = 0.009 + 0.3*sqr(6.0*C_T/(LiftCurveSlope*Solidity));
 503
 504   Torque = rho * BladeNum * BladeChord * delta_dr * sqr(Omega*Radius) * R[2] *
 505            (1.0+4.5*sqr(mu))/8.0
 506                      - (Thrust*lambda + H_drag*mu)*Radius;
 507
 508   return;
 509 }
 510
 511 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 512
 513 // transform rotor forces from control axes to shaft axes, and express
 514 // in body axes /SH79/ eqn(40,41)
 515
 516 FGColumnVector3 FGRotor::body_forces(double a_ic, double b_ic)
 517 {
 518   FGColumnVector3 F_s(
 519         - H_drag*cos(beta_orient) - J_side*sin(beta_orient) + Thrust*b_ic,
 520         - H_drag*sin(beta_orient) + J_side*cos(beta_orient) + Thrust*a_ic,
 521         - Thrust);
 522
 523   return HsrToTbo * F_s;
 524 }
 525
 526 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 527
 528 // calculates the additional moments due to hinge offset and handles
 529 // torque and sense
 530
 531 FGColumnVector3 FGRotor::body_moments(double a_ic, double b_ic)
 532 {
 533   FGColumnVector3 M_s, M_hub, M_h;
 534   double mf;
 535
 536   // cyclic flapping relative to shaft axes /SH79/ eqn(43)
 537   a1s = a_1*cos(beta_orient) + b_1*sin(beta_orient) - b_ic;
 538   b1s = b_1*cos(beta_orient) - a_1*sin(beta_orient) + a_ic;
 539
 540   mf = 0.5 * HingeOffset * BladeNum * Omega*Omega * BladeMassMoment;
 541
 542   M_s(eL) = mf*b1s;
 543   M_s(eM) = mf*a1s;
 544   M_s(eN) = Torque * Sense ;
 545
 546   return HsrToTbo * M_s;
 547 }
 548
 549 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 550
 551 void FGRotor::CalcStatePart1(void)
 552 {
 553   double A_IC;       // lateral (roll) control in radians
 554   double B_IC;       // longitudinal (pitch) control in radians
 555   double theta_col;  // rotor collective pitch in radians
 556
 557   double Vt ;
 558
 559   FGColumnVector3 UVW_h, PQR_h;
 560   FGColumnVector3 vHub_ca, avFus_ca;
 561
 562   double h_agl_ft, filtered_hagl = 0.0;
 563   double ge_factor = 1.0;
 564
 565   // fetch needed values from environment
 566   Vt = fdmex->GetAuxiliary()->GetVt(); // total vehicle velocity including wind
 567   dt = fdmex->GetDeltaT();
 568   rho = fdmex->GetAtmosphere()->GetDensity(); // slugs/ft^3.
 569   UVW_h = fdmex->GetAuxiliary()->GetAeroUVW();
 570   PQR_h = fdmex->GetAuxiliary()->GetAeroPQR();
 571   h_agl_ft = fdmex->GetPropagate()->GetDistanceAGL();
 572   // update InvTransform, the rotor orientation could have been altered
 573   InvTransform = Transform().Transposed();
 574
 575   // handle RPM requirements, calc omega.
 576   if (ExternalRPM && ExtRPMsource) {
 577     RPM = ExtRPMsource->getDoubleValue() / GearRatio;
 578   }
 579
 580   if (RPM < 1.0) { // kludge, otherwise calculations go bananas
 581     RPM = 1.0;
 582   }
 583
 584   Omega = (RPM/60.0)*2.0*M_PI;
 585
 586   // set control inputs
 587   A_IC      = LateralCtrl;
 588   B_IC      = LongitudinalCtrl;
 589   theta_col = CollectiveCtrl;
 590
 591   // ground effect
 592   if (GroundEffectExp > 1e-5) {
 593     if (h_agl_ft<0.0) h_agl_ft = 0.0; // clamp
 594     filtered_hagl = damp_hagl.execute(h_agl_ft) + GroundEffectShift;
 595     // actual/nominal factor avoids absurd scales at startup
 596     ge_factor -= exp(-filtered_hagl*GroundEffectExp) * (RPM / NominalRPM);
 597     if (ge_factor<0.5) ge_factor=0.5; // clamp
 598   }
 599
 600   // all set, start calculations
 601
 602   vHub_ca  = hub_vel_body2ca(UVW_h, PQR_h, A_IC, B_IC);
 603
 604   avFus_ca = fus_angvel_body2ca(PQR_h);
 605
 606   calc_flow_and_thrust(theta_col, vHub_ca(eU), vHub_ca(eW), ge_factor);
 607
 608   calc_coning_angle(theta_col);
 609
 610   calc_flapping_angles(theta_col, avFus_ca);
 611
 612   calc_drag_and_side_forces(theta_col);
 613
 614   calc_torque(theta_col);
 615
 616   // Fixme: only valid for a 'decent' rotor
 617   theta_downwash = atan2( - UVW_h(eU), v_induced - UVW_h(eW));
 618   phi_downwash   = atan2(   UVW_h(eV), v_induced - UVW_h(eW));
 619
 620   vFn = body_forces(A_IC, B_IC);
 621   vMn = Transform() * body_moments(A_IC, B_IC);
 622
 623 }
 624
 625 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 626
 627 void FGRotor::CalcStatePart2(double PowerAvailable)
 628 {
 629   if (! ExternalRPM) {
 630     // calculate new RPM
 631     double ExcessTorque = PowerAvailable / Omega;
 632     double deltaOmega   = ExcessTorque / PolarMoment * dt;
 633     RPM += deltaOmega/(2.0*M_PI) * 60.0;
 634     if (RPM < 0.0) RPM = 0.0; // Engine won't turn backwards
 635   }
 636 }
 637
 638 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 639
 640 double FGRotor::GetPowerRequired(void)
 641 {
 642   CalcStatePart1();
 643   PowerRequired = Torque * Omega;
 644   return PowerRequired;
 645 }
 646
 647 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 648
 649 double FGRotor::Calculate(double PowerAvailable)
 650 {
 651   CalcStatePart2(PowerAvailable);
 652   return Thrust;
 653 }
 654
 655
 656 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 657
 658
 659 bool FGRotor::BindModel(void)
 660 {
 661   string property_name, base_property_name;
 662   base_property_name = CreateIndexedPropertyName("propulsion/engine", EngineNum);
 663
 664   property_name = base_property_name + "/rotor-rpm";
 665   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetRPM );
 666
 667   property_name = base_property_name + "/x-engine-rpm"; // used for RPM eXchange
 668   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetEngineRPM );
 669
 670   property_name = base_property_name + "/rotor-thrust-lbs"; // might be redundant - check!
 671   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetThrust );
 672
 673   property_name = base_property_name + "/a0-rad";
 674   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA0 );
 675
 676   property_name = base_property_name + "/a1-rad";
 677   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA1 );
 678
 679   property_name = base_property_name + "/b1-rad";
 680   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetB1 );
 681
 682   property_name = base_property_name + "/inflow-ratio";
 683   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLambda );
 684
 685   property_name = base_property_name + "/advance-ratio";
 686   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetMu );
 687
 688   property_name = base_property_name + "/induced-inflow-ratio";
 689   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetNu );
 690
 691   property_name = base_property_name + "/vi-fps";
 692   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetVi );
 693
 694   property_name = base_property_name + "/thrust-coefficient";
 695   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCT );
 696
 697   property_name = base_property_name + "/torque-lbsft";
 698   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetTorque );
 699
 700   property_name = base_property_name + "/theta-downwash-rad";
 701   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetThetaDW );
 702
 703   property_name = base_property_name + "/phi-downwash-rad";
 704   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetPhiDW );
 705
 706   switch (ControlMap) {
 707     case eTailCtrl:
 708       property_name = base_property_name + "/antitorque-ctrl-rad";
 709       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
 710       break;
 711     case eTandemCtrl:
 712       property_name = base_property_name + "/tail-collective-ctrl-rad";
 713       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
 714       property_name = base_property_name + "/lateral-ctrl-rad";
 715       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
 716       property_name = base_property_name + "/longitudinal-ctrl-rad";
 717       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
 718       break;
 719     default: // eMainCtrl
 720       property_name = base_property_name + "/collective-ctrl-rad";
 721       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
 722       property_name = base_property_name + "/lateral-ctrl-rad";
 723       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
 724       property_name = base_property_name + "/longitudinal-ctrl-rad";
 725       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
 726   }
 727
 728   if (ExternalRPM) {
 729     if (RPMdefinition == -1) {
 730       property_name = base_property_name + "/x-rpm-dict";
 731       ExtRPMsource = PropertyManager->GetNode(property_name, true);
 732     } else if (RPMdefinition >= 0 && RPMdefinition != EngineNum) {
 733       string ipn = CreateIndexedPropertyName("propulsion/engine", RPMdefinition);
 734       property_name = ipn + "/x-engine-rpm";
 735       ExtRPMsource = PropertyManager->GetNode(property_name, false);
 736       if (! ExtRPMsource) {
 737         cerr << "# Warning: Engine number " << EngineNum << "." << endl;
 738         cerr << "# No 'x-engine-rpm' property found for engine " << RPMdefinition << "." << endl;
 739         cerr << "# Please check order of engine definitons."  << endl;
 740       }
 741     } else {
 742       cerr << "# Engine number " << EngineNum;
 743       cerr << ", given ExternalRPM value '" << RPMdefinition << "' unhandled."  << endl;
 744     }
 745   }
 746
 747   return true;
 748 }
 749
 750 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 751
 752 string FGRotor::GetThrusterLabels(int id, string delimeter)
 753 {
 754
 755   ostringstream buf;
 756
 757   buf << Name << " RPM (engine " << id << ")";
 758
 759   return buf.str();
 760
 761 }
 762
 763 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 764
 765 string FGRotor::GetThrusterValues(int id, string delimeter)
 766 {
 767
 768   ostringstream buf;
 769
 770   buf << RPM;
 771
 772   return buf.str();
 773
 774 }
 775
 776 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 777 //    The bitmasked value choices are as follows:
 778 //    unset: In this case (the default) JSBSim would only print
 779 //       out the normally expected messages, essentially echoing
 780 //       the config files as they are read. If the environment
 781 //       variable is not set, debug_lvl is set to 1 internally
 782 //    0: This requests JSBSim not to output any messages
 783 //       whatsoever.
 784 //    1: This value explicity requests the normal JSBSim
 785 //       startup messages
 786 //    2: This value asks for a message to be printed out when
 787 //       a class is instantiated
 788 //    4: When this value is set, a message is displayed when a
 789 //       FGModel object executes its Run() method
 790 //    8: When this value is set, various runtime state variables
 791 //       are printed out periodically
 792 //    16: When set various parameters are sanity checked and
 793 //       a message is printed out when they go out of bounds
 794
 795 void FGRotor::Debug(int from)
 796 {
 797   string ControlMapName;
 798
 799   if (debug_lvl <= 0) return;
 800
 801   if (debug_lvl & 1) { // Standard console startup message output
 802     if (from == 0) { // Constructor
 803       cout << "\n    Rotor Name: " << Name << endl;
 804       cout << "      Diameter = " << 2.0 * Radius << " ft." << endl;
 805       cout << "      Number of Blades = " << BladeNum << endl;
 806       cout << "      Gear Ratio = " << GearRatio << endl;
 807       cout << "      Sense = " << Sense << endl;
 808       cout << "      Nominal RPM = " << NominalRPM << endl;
 809
 810       if (ExternalRPM) {
 811         if (RPMdefinition == -1) {
 812           cout << "      RPM is controlled externally" << endl;
 813         } else {
 814           cout << "      RPM source set to engine " << RPMdefinition << endl;
 815         }
 816       }
 817
 818       cout << "      Blade Chord = " << BladeChord << endl;
 819       cout << "      Lift Curve Slope = " << LiftCurveSlope << endl;
 820       cout << "      Blade Twist = " << BladeTwist << endl;
 821       cout << "      Hinge Offset = " << HingeOffset << endl;
 822       cout << "      Blade Flapping Moment = " << BladeFlappingMoment << endl;
 823       cout << "      Blade Mass Moment = " << BladeMassMoment << endl;
 824       cout << "      Polar Moment = " << PolarMoment << endl;
 825       cout << "      Inflow Lag = " << InflowLag << endl;
 826       cout << "      Tip Loss = " << TipLossB << endl;
 827       cout << "      Lock Number = " << LockNumberByRho * 0.002356 << " (SL)" << endl;
 828       cout << "      Solidity = " << Solidity << endl;
 829
 830       switch (ControlMap) {
 831         case eTailCtrl:    ControlMapName = "Tail Rotor";   break;
 832         case eTandemCtrl:  ControlMapName = "Tandem Rotor"; break;
 833         default:           ControlMapName = "Main Rotor";
 834       }
 835       cout << "      Control Mapping = " << ControlMapName << endl;
 836
 837     }
 838   }
 839   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 840     if (from == 0) cout << "Instantiated: FGRotor" << endl;
 841     if (from == 1) cout << "Destroyed:    FGRotor" << endl;
 842   }
 843   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 844   }
 845   if (debug_lvl & 8 ) { // Runtime state variables
 846   }
 847   if (debug_lvl & 16) { // Sanity checking
 848   }
 849   if (debug_lvl & 64) {
 850     if (from == 0) { // Constructor
 851       cout << IdSrc << endl;
 852       cout << IdHdr << endl;
 853     }
 854   }
 855
 856 }
 857
 858
 859 } // namespace JSBSim
 860