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