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