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 02/05/12  T.Kreitler brake, clutch, and FWU now in FGTransmission,
  40                      downwash angles relate to shaft orientation
  41
  42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  43 INCLUDES
  44 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  45
  46 #include <sstream>
  47
  48 #include "FGRotor.h"
  49 #include "models/FGMassBalance.h"
  50 #include "models/FGPropulsion.h" // to get the GearRatio from a linked rotor
  51
  52 using std::cerr;
  53 using std::cout;
  54 using std::endl;
  55 using std::ostringstream;
  56
  57 namespace JSBSim {
  58
  59 static const char *IdSrc = "$Id: FGRotor.cpp,v 1.20 2012/03/18 15:48:36 jentron 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 // Constructor
  73
  74 FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
  75   : FGThruster(exec, rotor_element, num),
  76     rho(0.002356),                                  // environment
  77     Radius(0.0), BladeNum(0),                       // configuration parameters
  78     Sense(1.0), NominalRPM(0.0), MinimalRPM(0.0), MaximalRPM(0.0),
  79     ExternalRPM(0), RPMdefinition(0), ExtRPMsource(NULL), SourceGearRatio(1.0),
  80     BladeChord(0.0), LiftCurveSlope(0.0), BladeTwist(0.0), HingeOffset(0.0),
  81     BladeFlappingMoment(0.0), BladeMassMoment(0.0), PolarMoment(0.0),
  82     InflowLag(0.0), TipLossB(0.0),
  83     GroundEffectExp(0.0), GroundEffectShift(0.0), GroundEffectScaleNorm(1.0),
  84     LockNumberByRho(0.0), Solidity(0.0),            // derived parameters
  85     RPM(0.0), Omega(0.0),                           // dynamic values
  86     beta_orient(0.0),
  87     a0(0.0), a_1(0.0), b_1(0.0), a_dw(0.0),
  88     a1s(0.0), b1s(0.0),
  89     H_drag(0.0), J_side(0.0), Torque(0.0), C_T(0.0),
  90     lambda(-0.001), mu(0.0), nu(0.001), v_induced(0.0),
  91     theta_downwash(0.0), phi_downwash(0.0),
  92     ControlMap(eMainCtrl),                          // control
  93     CollectiveCtrl(0.0), LateralCtrl(0.0), LongitudinalCtrl(0.0),
  94     Transmission(NULL),                             // interaction with engine
  95     EngineRPM(0.0), MaxBrakePower(0.0), GearLoss(0.0), GearMoment(0.0)
  96 {
  97   FGColumnVector3 location(0.0, 0.0, 0.0), orientation(0.0, 0.0, 0.0);
  98   Element *thruster_element;
  99   double engine_power_est = 0.0;
 100
 101   // initialise/set remaining variables
 102   SetTransformType(FGForce::tCustom);
 103   PropertyManager = exec->GetPropertyManager();
 104   Type = ttRotor;
 105   GearRatio = 1.0;
 106
 107   dt = exec->GetDeltaT();
 108   for (int i=0; i<5; i++) R[i] = 0.0;
 109   for (int i=0; i<5; i++) B[i] = 0.0;
 110
 111   // get positions
 112   thruster_element = rotor_element->GetParent()->FindElement("sense");
 113   if (thruster_element) {
 114     double s = thruster_element->GetDataAsNumber();
 115     if (s < -0.1) {
 116       Sense = -1.0; // 'CW' as seen from above
 117     } else if (s < 0.1) {
 118       Sense = 0.0;  // 'coaxial'
 119     } else {
 120       Sense = 1.0; // 'CCW' as seen from above
 121     }
 122   }
 123
 124   thruster_element = rotor_element->GetParent()->FindElement("location");
 125   if (thruster_element) {
 126     location = thruster_element->FindElementTripletConvertTo("IN");
 127   } else {
 128     cerr << "No thruster location found." << endl;
 129   }
 130
 131   thruster_element = rotor_element->GetParent()->FindElement("orient");
 132   if (thruster_element) {
 133     orientation = thruster_element->FindElementTripletConvertTo("RAD");
 134   } else {
 135     cerr << "No thruster orientation found." << endl;
 136   }
 137
 138   SetLocation(location);
 139   SetAnglesToBody(orientation);
 140   InvTransform = Transform().Transposed(); // body to custom/native
 141
 142   // wire controls
 143   ControlMap = eMainCtrl;
 144   if (rotor_element->FindElement("controlmap")) {
 145     string cm = rotor_element->FindElementValue("controlmap");
 146     cm = to_upper(cm);
 147     if (cm == "TAIL") {
 148       ControlMap = eTailCtrl;
 149     } else if (cm == "TANDEM") {
 150       ControlMap = eTandemCtrl;
 151     } else {
 152       cerr << "# found unknown controlmap: '" << cm << "' using main rotor config."  << endl;
 153     }
 154   }
 155
 156   // ExternalRPM -- is the RPM dictated ?
 157   if (rotor_element->FindElement("ExternalRPM")) {
 158     ExternalRPM = 1;
 159     SourceGearRatio = 1.0;
 160     RPMdefinition = (int) rotor_element->FindElementValueAsNumber("ExternalRPM");
 161     int rdef = RPMdefinition;
 162     if (RPMdefinition>=0) {
 163       // avoid ourself and (still) unknown engines.
 164       if (!exec->GetPropulsion()->GetEngine(RPMdefinition) || RPMdefinition==num) {
 165         RPMdefinition = -1;
 166       } else {
 167         FGThruster *tr = exec->GetPropulsion()->GetEngine(RPMdefinition)->GetThruster();
 168         SourceGearRatio = tr->GetGearRatio();
 169         // cout << "# got sources' GearRatio: " << SourceGearRatio << endl;
 170       }
 171     }
 172     if (RPMdefinition != rdef) {
 173       cerr << "# discarded given RPM source (" << rdef << ") and switched to external control (-1)." << endl;
 174     }
 175   }
 176
 177   // process rotor parameters
 178   engine_power_est = Configure(rotor_element);
 179
 180   // setup transmission if needed
 181   if (!ExternalRPM) {
 182
 183     Transmission = new FGTransmission(exec, num, dt);
 184
 185     Transmission->SetThrusterMoment(PolarMoment);
 186
 187     // The MOI sensed behind the gear ( MOI_engine*sqr(GearRatio) ).
 188     GearMoment = ConfigValueConv(rotor_element, "gearmoment", 0.1*PolarMoment, "SLUG*FT2");
 189     GearMoment = Constrain(1e-6, GearMoment, 1e9);
 190     Transmission->SetEngineMoment(GearMoment);
 191
 192     Transmission->SetMaxBrakePower(MaxBrakePower);
 193
 194     GearLoss = ConfigValueConv(rotor_element, "gearloss", 0.0025 * engine_power_est, "HP");
 195     GearLoss = Constrain(0.0, GearLoss, 1e9);
 196     GearLoss *= hptoftlbssec;
 197     Transmission->SetEngineFriction(GearLoss);
 198
 199   }
 200
 201   // shaft representation - a rather simple transform,
 202   // but using a matrix is safer.
 203   TboToHsr.InitMatrix(   0.0, 0.0, 1.0,
 204                          0.0, 1.0, 0.0,
 205                         -1.0, 0.0, 0.0  );
 206   HsrToTbo  =  TboToHsr.Transposed();
 207
 208   // smooth out jumps in hagl reported, otherwise the ground effect
 209   // calculation would cause jumps too. 1Hz seems sufficient.
 210   damp_hagl = Filter(1.0, dt);
 211
 212   // enable import-export
 213   BindModel();
 214
 215   Debug(0);
 216
 217 }  // Constructor
 218
 219 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 220
 221 FGRotor::~FGRotor(){
 222   if (Transmission) delete Transmission;
 223   Debug(1);
 224 }
 225
 226 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 227
 228 // 5in1: value-fetch-convert-default-return function
 229
 230 double FGRotor::ConfigValueConv( Element* el, const string& ename, double default_val,
 231                                   const string& unit, bool tell)
 232 {
 233
 234   Element *e = NULL;
 235   double val = default_val;
 236
 237   string pname = "*No parent element*";
 238
 239   if (el) {
 240     e = el->FindElement(ename);
 241     pname = el->GetName();
 242   }
 243
 244   if (e) {
 245     if (unit.empty()) {
 246       val = e->GetDataAsNumber();
 247     } else {
 248       val = el->FindElementValueAsNumberConvertTo(ename,unit);
 249     }
 250   } else {
 251     if (tell) {
 252       cerr << pname << ": missing element '" << ename <<
 253                        "' using estimated value: " << default_val << endl;
 254     }
 255   }
 256
 257   return val;
 258 }
 259
 260 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 261
 262 double FGRotor::ConfigValue(Element* el, const string& ename, double default_val, bool tell)
 263 {
 264   return ConfigValueConv(el, ename, default_val, "", tell);
 265 }
 266
 267 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 268
 269 // 1. read configuration and try to fill holes, ymmv
 270 // 2. calculate derived parameters
 271 double FGRotor::Configure(Element* rotor_element)
 272 {
 273
 274   double estimate, engine_power_est=0.0;
 275   const bool yell   = true;
 276   const bool silent = false;
 277
 278
 279   Radius = 0.5 * ConfigValueConv(rotor_element, "diameter", 42.0, "FT", yell);
 280   Radius = Constrain(1e-3, Radius, 1e9);
 281
 282   BladeNum = (int) ConfigValue(rotor_element, "numblades", 3 , yell);
 283
 284   GearRatio = ConfigValue(rotor_element, "gearratio", 1.0, yell);
 285   GearRatio = Constrain(1e-9, GearRatio, 1e9);
 286
 287   // make sure that v_tip (omega*r) is below 0.7mach ~ 750ft/s
 288   estimate = (750.0/Radius)/(2.0*M_PI) * 60.0;  // 7160/Radius
 289   NominalRPM = ConfigValue(rotor_element, "nominalrpm", estimate, yell);
 290   NominalRPM = Constrain(2.0, NominalRPM, 1e9);
 291
 292   MinimalRPM = ConfigValue(rotor_element, "minrpm", 1.0);
 293   MinimalRPM = Constrain(1.0, MinimalRPM, NominalRPM - 1.0);
 294
 295   MaximalRPM = ConfigValue(rotor_element, "maxrpm", 2.0*NominalRPM);
 296   MaximalRPM = Constrain(NominalRPM, MaximalRPM, 1e9);
 297
 298   estimate = Constrain(0.07, 2.0/Radius , 0.14); // guess solidity
 299   estimate = estimate * M_PI*Radius/BladeNum;
 300   BladeChord = ConfigValueConv(rotor_element, "chord", estimate, "FT", yell);
 301
 302   LiftCurveSlope = ConfigValue(rotor_element, "liftcurveslope", 6.0); // "1/RAD"
 303   BladeTwist = ConfigValueConv(rotor_element, "twist", -0.17, "RAD");
 304
 305   HingeOffset = ConfigValueConv(rotor_element, "hingeoffset", 0.05 * Radius, "FT" );
 306
 307   estimate = sqr(BladeChord) * sqr(Radius - HingeOffset) * 0.57;
 308   BladeFlappingMoment = ConfigValueConv(rotor_element, "flappingmoment", estimate, "SLUG*FT2");
 309   BladeFlappingMoment = Constrain(1e-9, BladeFlappingMoment, 1e9);
 310
 311   // guess mass from moment of a thin stick, and multiply by the blades cg distance
 312   estimate = ( 3.0 * BladeFlappingMoment / sqr(Radius) ) * (0.45 * Radius) ;
 313   BladeMassMoment = ConfigValue(rotor_element, "massmoment", estimate); // unit is slug-ft
 314   BladeMassMoment = Constrain(1e-9, BladeMassMoment, 1e9);
 315
 316   estimate = 1.1 * BladeFlappingMoment * BladeNum;
 317   PolarMoment = ConfigValueConv(rotor_element, "polarmoment", estimate, "SLUG*FT2");
 318   PolarMoment = Constrain(1e-9, PolarMoment, 1e9);
 319
 320   // "inflowlag" is treated further down.
 321
 322   TipLossB = ConfigValue(rotor_element, "tiplossfactor", 1.0, silent);
 323
 324   // estimate engine power (bit of a pity, cause our caller already knows)
 325   engine_power_est = 0.5 * BladeNum*BladeChord*Radius*Radius;
 326
 327   estimate = engine_power_est / 30.0;
 328   MaxBrakePower  = ConfigValueConv(rotor_element, "maxbrakepower", estimate, "HP");
 329   MaxBrakePower *= hptoftlbssec;
 330
 331   GroundEffectExp = ConfigValue(rotor_element, "groundeffectexp", 0.0);
 332   GroundEffectShift = ConfigValueConv(rotor_element, "groundeffectshift", 0.0, "FT");
 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(1e-6, InflowLag, 2.0);
 348
 349   return engine_power_est;
 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   if (mu > 0.7) mu = 0.7;
 418   mu2 = sqr(mu);
 419
 420   ct_t0 = (1.0/3.0*B[3] + 1.0/2.0 * TipLossB*mu2 - 4.0/(9.0*M_PI) * mu*mu2 ) * theta_0;
 421   ct_t1 = (1.0/4.0*B[4] + 1.0/4.0 * B[2]*mu2) * BladeTwist;
 422
 423   ct_l  = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda; // first time
 424
 425   c0 = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1) * Solidity;
 426   c0 = c0 / ( 2.0 * sqrt( sqr(mu) + sqr(lambda) ) + 1e-15);
 427
 428   // replacement for /SH79/ eqn(26).
 429   // ref: dnu/dt = 1/tau ( Ct / (2*sqrt(mu^2+lambda^2))  -  nu )
 430   // taking mu and lambda constant, this integrates to
 431
 432   nu  = flow_scale * ((nu - c0) * exp(-dt/InflowLag) + c0);
 433
 434   // now from nu to lambda, C_T, and Thrust
 435
 436   lambda = Ww/(Omega*Radius) - nu; // /SH79/ eqn(25)
 437
 438   ct_l  = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda;
 439
 440   ct_over_sigma = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1); // /SH79/ eqn(27)
 441
 442   Thrust = BladeNum*BladeChord*Radius*rho*sqr(Omega*Radius) * ct_over_sigma;
 443
 444   C_T = ct_over_sigma * Solidity;
 445   v_induced = nu * (Omega*Radius);
 446
 447 }
 448
 449
 450 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 451
 452 // Two blade teetering rotors are often 'preconed' to a fixed angle, but the
 453 // calculated value is pretty close to the real one. /SH79/ eqn(29)
 454
 455 void FGRotor::calc_coning_angle(double theta_0)
 456 {
 457   double lock_gamma = LockNumberByRho * rho;
 458
 459   double a0_l  = (1.0/6.0  + 0.04 * mu*mu*mu) * lambda;
 460   double a0_t0 = (1.0/8.0  + 1.0/8.0  * mu*mu) * theta_0;
 461   double a0_t1 = (1.0/10.0 + 1.0/12.0 * mu*mu) * BladeTwist;
 462   a0 = lock_gamma * ( a0_l + a0_t0 + a0_t1);
 463   return;
 464 }
 465
 466 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 467
 468 // Flapping angles relative to control axes /SH79/ eqn(32)
 469
 470 void FGRotor::calc_flapping_angles(double theta_0, const FGColumnVector3 &pqr_fus_w)
 471 {
 472   double lock_gamma = LockNumberByRho * rho;
 473
 474
 475   double mu2_2 = sqr(mu)/2.0;
 476   double t075 = theta_0 + 0.75 * BladeTwist;  // common approximation for rectangular blades
 477
 478   a_1 = 1.0/(1.0 - mu2_2) * (
 479                                  (2.0*lambda + (8.0/3.0)*t075)*mu
 480                                + pqr_fus_w(eP)/Omega
 481                                - 16.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
 482                              );
 483
 484   b_1 = 1.0/(1.0 + mu2_2) * (
 485                                  (4.0/3.0)*mu*a0
 486                                - pqr_fus_w(eQ)/Omega
 487                                - 16.0 * pqr_fus_w(eP)/(lock_gamma*Omega)
 488                              );
 489
 490   // used in  force calc
 491   a_dw = 1.0/(1.0 - mu2_2) * (
 492                                  (2.0*lambda + (8.0/3.0)*t075)*mu
 493                                - 24.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
 494                                  * ( 1.0 - ( 0.29 * t075 / (C_T/Solidity) ) )
 495                              );
 496
 497   return;
 498 }
 499
 500 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 501
 502 // /SH79/ eqn(38,39)
 503
 504 void FGRotor::calc_drag_and_side_forces(double theta_0)
 505 {
 506   double cy_over_sigma;
 507   double t075 = theta_0 + 0.75 * BladeTwist;
 508
 509   H_drag = Thrust * a_dw;
 510
 511   cy_over_sigma = (
 512                       0.75*b_1*lambda - 1.5*a0*mu*lambda + 0.25*a_1*b_1*mu
 513                     - a0*a_1*sqr(mu) + (1.0/6.0)*a0*a_1
 514                     - (0.75*mu*a0 - (1.0/3.0)*b_1 - 0.5*sqr(mu)*b_1)*t075
 515                   );
 516   cy_over_sigma *= LiftCurveSlope/2.0;
 517
 518   J_side = BladeNum * BladeChord * Radius * rho * sqr(Omega*Radius) * cy_over_sigma;
 519
 520   return;
 521 }
 522
 523 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 524
 525 // Simplified version of /SH79/ eqn(36). Uses an estimate for blade drag
 526 // (a new config parameter to come...).
 527 // From "Bramwell's Helicopter Dynamics", second edition, eqn(3.43) and (3.44)
 528
 529 void FGRotor::calc_torque(double theta_0)
 530 {
 531   // estimate blade drag
 532   double delta_dr = 0.009 + 0.3*sqr(6.0*C_T/(LiftCurveSlope*Solidity));
 533
 534   Torque = rho * BladeNum * BladeChord * delta_dr * sqr(Omega*Radius) * R[2] *
 535            (1.0+4.5*sqr(mu))/8.0
 536                      - (Thrust*lambda + H_drag*mu)*Radius;
 537
 538   return;
 539 }
 540
 541 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 542
 543 // Get the downwash angles with respect to the shaft axis.
 544 // Given a 'regular'  main rotor, the angles are zero when the downwash points
 545 // down, positive theta values mean that the downwash turns towards the nose,
 546 // and positive phi values mean it turns to the left side. (Note: only airspeed
 547 // is transformed, the rotational speed contribution is ignored.)
 548
 549 void FGRotor::calc_downwash_angles()
 550 {
 551   FGColumnVector3 v_shaft;
 552   v_shaft = TboToHsr * InvTransform * in.AeroUVW;
 553
 554   theta_downwash = atan2( -v_shaft(eU), v_induced - v_shaft(eW)) + a1s;
 555   phi_downwash   = atan2(  v_shaft(eV), v_induced - v_shaft(eW)) + b1s;
 556
 557   return;
 558 }
 559
 560 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 561
 562 // transform rotor forces from control axes to shaft axes, and express
 563 // in body axes /SH79/ eqn(40,41)
 564
 565 FGColumnVector3 FGRotor::body_forces(double a_ic, double b_ic)
 566 {
 567   FGColumnVector3 F_s(
 568         - H_drag*cos(beta_orient) - J_side*sin(beta_orient) + Thrust*b_ic,
 569         - H_drag*sin(beta_orient) + J_side*cos(beta_orient) + Thrust*a_ic,
 570         - Thrust);
 571
 572   return HsrToTbo * F_s;
 573 }
 574
 575 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 576
 577 // calculates the additional moments due to hinge offset and handles
 578 // torque and sense
 579
 580 FGColumnVector3 FGRotor::body_moments(double a_ic, double b_ic)
 581 {
 582   FGColumnVector3 M_s, M_hub, M_h;
 583   double mf;
 584
 585   // cyclic flapping relative to shaft axes /SH79/ eqn(43)
 586   a1s = a_1*cos(beta_orient) + b_1*sin(beta_orient) - b_ic;
 587   b1s = b_1*cos(beta_orient) - a_1*sin(beta_orient) + a_ic;
 588
 589   mf = 0.5 * HingeOffset * BladeNum * Omega*Omega * BladeMassMoment;
 590
 591   M_s(eL) = mf*b1s;
 592   M_s(eM) = mf*a1s;
 593   M_s(eN) = Torque * Sense ;
 594
 595   return HsrToTbo * M_s;
 596 }
 597
 598 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 599
 600 void FGRotor::CalcRotorState(void)
 601 {
 602   double A_IC;       // lateral (roll) control in radians
 603   double B_IC;       // longitudinal (pitch) control in radians
 604   double theta_col;  // rotor collective pitch in radians
 605
 606   FGColumnVector3 vHub_ca, avFus_ca;
 607
 608   double filtered_hagl = 0.0;
 609   double ge_factor = 1.0;
 610
 611   // fetch needed values from environment
 612   rho = in.Density; // slugs/ft^3.
 613   double h_agl_ft = in.H_agl;
 614
 615   // update InvTransform, the rotor orientation could have been altered
 616   InvTransform = Transform().Transposed();
 617
 618   // handle RPM requirements, calc omega.
 619   if (ExternalRPM && ExtRPMsource) {
 620     RPM = ExtRPMsource->getDoubleValue() * ( SourceGearRatio / GearRatio );
 621   }
 622
 623   // MinimalRPM is always >= 1. MaximalRPM is always >= NominalRPM
 624   RPM = Constrain(MinimalRPM, RPM, MaximalRPM);
 625
 626   Omega = (RPM/60.0)*2.0*M_PI;
 627
 628   // set control inputs
 629   A_IC      = LateralCtrl;
 630   B_IC      = LongitudinalCtrl;
 631   theta_col = CollectiveCtrl;
 632
 633   // optional ground effect, a ge_factor of 1.0 gives no effect
 634   // and 0.5 yields to maximal influence.
 635   if (GroundEffectExp > 1e-5) {
 636     if (h_agl_ft<0.0) h_agl_ft = 0.0; // clamp
 637     filtered_hagl = damp_hagl.execute(h_agl_ft) + GroundEffectShift;
 638     // actual/nominal factor avoids absurd scales at startup
 639     ge_factor -= GroundEffectScaleNorm *
 640                  ( exp(-filtered_hagl*GroundEffectExp) * (RPM / NominalRPM) );
 641     ge_factor = Constrain(0.5, ge_factor, 1.0);
 642   }
 643
 644   // all set, start calculations ...
 645
 646   vHub_ca  = hub_vel_body2ca(in.AeroUVW, in.AeroPQR, A_IC, B_IC);
 647
 648   avFus_ca = fus_angvel_body2ca(in.AeroPQR);
 649
 650   calc_flow_and_thrust(theta_col, vHub_ca(eU), vHub_ca(eW), ge_factor);
 651
 652   calc_coning_angle(theta_col);
 653
 654   calc_flapping_angles(theta_col, avFus_ca);
 655
 656   calc_drag_and_side_forces(theta_col);
 657
 658   calc_torque(theta_col);
 659
 660   calc_downwash_angles();
 661
 662   // ... and assign to inherited vFn and vMn members
 663   //     (for processing see FGForce::GetBodyForces).
 664   vFn = body_forces(A_IC, B_IC);
 665   vMn = Transform() * body_moments(A_IC, B_IC);
 666
 667 }
 668
 669 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 670
 671 double FGRotor::Calculate(double EnginePower)
 672 {
 673
 674   CalcRotorState();
 675
 676   if (! ExternalRPM) {
 677     // the RPM values are handled inside Transmission
 678     Transmission->Calculate(EnginePower, Torque, in.TotalDeltaT);
 679
 680     EngineRPM = Transmission->GetEngineRPM() * GearRatio;
 681     RPM = Transmission->GetThrusterRPM();
 682   } else {
 683     EngineRPM = RPM * GearRatio;
 684   }
 685
 686   RPM = Constrain(MinimalRPM, RPM, MaximalRPM); // trim again
 687
 688   return Thrust;
 689 }
 690
 691
 692 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 693
 694
 695 bool FGRotor::BindModel(void)
 696 {
 697   string property_name, base_property_name;
 698   base_property_name = CreateIndexedPropertyName("propulsion/engine", EngineNum);
 699
 700   property_name = base_property_name + "/rotor-rpm";
 701   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetRPM );
 702
 703   property_name = base_property_name + "/engine-rpm";
 704   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetEngineRPM );
 705
 706   property_name = base_property_name + "/a0-rad";
 707   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA0 );
 708
 709   property_name = base_property_name + "/a1-rad";
 710   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA1 );
 711
 712   property_name = base_property_name + "/b1-rad";
 713   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetB1 );
 714
 715   property_name = base_property_name + "/inflow-ratio";
 716   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLambda );
 717
 718   property_name = base_property_name + "/advance-ratio";
 719   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetMu );
 720
 721   property_name = base_property_name + "/induced-inflow-ratio";
 722   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetNu );
 723
 724   property_name = base_property_name + "/vi-fps";
 725   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetVi );
 726
 727   property_name = base_property_name + "/thrust-coefficient";
 728   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCT );
 729
 730   property_name = base_property_name + "/torque-lbsft";
 731   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetTorque );
 732
 733   property_name = base_property_name + "/theta-downwash-rad";
 734   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetThetaDW );
 735
 736   property_name = base_property_name + "/phi-downwash-rad";
 737   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetPhiDW );
 738
 739   property_name = base_property_name + "/groundeffect-scale-norm";
 740   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetGroundEffectScaleNorm,
 741                                                      &FGRotor::SetGroundEffectScaleNorm );
 742
 743   switch (ControlMap) {
 744     case eTailCtrl:
 745       property_name = base_property_name + "/antitorque-ctrl-rad";
 746       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
 747       break;
 748     case eTandemCtrl:
 749       property_name = base_property_name + "/tail-collective-ctrl-rad";
 750       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
 751       property_name = base_property_name + "/lateral-ctrl-rad";
 752       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
 753       property_name = base_property_name + "/longitudinal-ctrl-rad";
 754       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
 755       break;
 756     default: // eMainCtrl
 757       property_name = base_property_name + "/collective-ctrl-rad";
 758       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
 759       property_name = base_property_name + "/lateral-ctrl-rad";
 760       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
 761       property_name = base_property_name + "/longitudinal-ctrl-rad";
 762       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
 763   }
 764
 765   if (ExternalRPM) {
 766     if (RPMdefinition == -1) {
 767       property_name = base_property_name + "/x-rpm-dict";
 768       ExtRPMsource = PropertyManager->GetNode(property_name, true);
 769     } else if (RPMdefinition >= 0 && RPMdefinition != EngineNum) {
 770       string ipn = CreateIndexedPropertyName("propulsion/engine", RPMdefinition);
 771       property_name = ipn + "/rotor-rpm";
 772       ExtRPMsource = PropertyManager->GetNode(property_name, false);
 773       if (! ExtRPMsource) {
 774         cerr << "# Warning: Engine number " << EngineNum << "." << endl;
 775         cerr << "# No 'rotor-rpm' property found for engine " << RPMdefinition << "." << endl;
 776         cerr << "# Please check order of engine definitons."  << endl;
 777       }
 778     } else {
 779       cerr << "# Engine number " << EngineNum;
 780       cerr << ", given ExternalRPM value '" << RPMdefinition << "' unhandled."  << endl;
 781     }
 782   }
 783
 784   return true;
 785 }
 786
 787 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 788
 789 string FGRotor::GetThrusterLabels(int id, const string& delimeter)
 790 {
 791
 792   ostringstream buf;
 793
 794   buf << Name << " RPM (engine " << id << ")";
 795
 796   return buf.str();
 797
 798 }
 799
 800 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 801
 802 string FGRotor::GetThrusterValues(int id, const string& delimeter)
 803 {
 804
 805   ostringstream buf;
 806
 807   buf << RPM;
 808
 809   return buf.str();
 810
 811 }
 812
 813 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 814 //    The bitmasked value choices are as follows:
 815 //    unset: In this case (the default) JSBSim would only print
 816 //       out the normally expected messages, essentially echoing
 817 //       the config files as they are read. If the environment
 818 //       variable is not set, debug_lvl is set to 1 internally
 819 //    0: This requests JSBSim not to output any messages
 820 //       whatsoever.
 821 //    1: This value explicity requests the normal JSBSim
 822 //       startup messages
 823 //    2: This value asks for a message to be printed out when
 824 //       a class is instantiated
 825 //    4: When this value is set, a message is displayed when a
 826 //       FGModel object executes its Run() method
 827 //    8: When this value is set, various runtime state variables
 828 //       are printed out periodically
 829 //    16: When set various parameters are sanity checked and
 830 //       a message is printed out when they go out of bounds
 831
 832 void FGRotor::Debug(int from)
 833 {
 834   string ControlMapName;
 835
 836   if (debug_lvl <= 0) return;
 837
 838   if (debug_lvl & 1) { // Standard console startup message output
 839     if (from == 0) { // Constructor
 840       cout << "\n    Rotor Name: " << Name << endl;
 841       cout << "      Diameter = " << 2.0 * Radius << " ft." << endl;
 842       cout << "      Number of Blades = " << BladeNum << endl;
 843       cout << "      Gear Ratio = " << GearRatio << endl;
 844       cout << "      Sense = " << Sense << endl;
 845       cout << "      Nominal RPM = " << NominalRPM << endl;
 846       cout << "      Minimal RPM = " << MinimalRPM << endl;
 847       cout << "      Maximal RPM = " << MaximalRPM << endl;
 848
 849       if (ExternalRPM) {
 850         if (RPMdefinition == -1) {
 851           cout << "      RPM is controlled externally" << endl;
 852         } else {
 853           cout << "      RPM source set to thruster " << RPMdefinition << endl;
 854         }
 855       }
 856
 857       cout << "      Blade Chord = " << BladeChord << endl;
 858       cout << "      Lift Curve Slope = " << LiftCurveSlope << endl;
 859       cout << "      Blade Twist = " << BladeTwist << endl;
 860       cout << "      Hinge Offset = " << HingeOffset << endl;
 861       cout << "      Blade Flapping Moment = " << BladeFlappingMoment << endl;
 862       cout << "      Blade Mass Moment = " << BladeMassMoment << endl;
 863       cout << "      Polar Moment = " << PolarMoment << endl;
 864       cout << "      Inflow Lag = " << InflowLag << endl;
 865       cout << "      Tip Loss = " << TipLossB << endl;
 866       cout << "      Lock Number = " << LockNumberByRho * 0.002356 << " (SL)" << endl;
 867       cout << "      Solidity = " << Solidity << endl;
 868       cout << "      Max Brake Power = " << MaxBrakePower/hptoftlbssec << " HP" << endl;
 869       cout << "      Gear Loss = " << GearLoss/hptoftlbssec << " HP" << endl;
 870       cout << "      Gear Moment = " << GearMoment << endl;
 871
 872       switch (ControlMap) {
 873         case eTailCtrl:    ControlMapName = "Tail Rotor";   break;
 874         case eTandemCtrl:  ControlMapName = "Tandem Rotor"; break;
 875         default:           ControlMapName = "Main Rotor";
 876       }
 877       cout << "      Control Mapping = " << ControlMapName << endl;
 878
 879     }
 880   }
 881   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 882     if (from == 0) cout << "Instantiated: FGRotor" << endl;
 883     if (from == 1) cout << "Destroyed:    FGRotor" << endl;
 884   }
 885   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 886   }
 887   if (debug_lvl & 8 ) { // Runtime state variables
 888   }
 889   if (debug_lvl & 16) { // Sanity checking
 890   }
 891   if (debug_lvl & 64) {
 892     if (from == 0) { // Constructor
 893       cout << IdSrc << endl;
 894       cout << IdHdr << endl;
 895     }
 896   }
 897
 898 }
 899
 900
 901 } // namespace JSBSim
 902