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.16 2011/09/17 16:39:19 bcoconni 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), MinimalRPM(0.0), MaximalRPM(0.0),
  81     ExternalRPM(0), 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   NominalRPM = Constrain(2.0, NominalRPM, 1e9);
 259
 260   MinimalRPM = ConfigValue(rotor_element, "minrpm", 1.0);
 261   MinimalRPM = Constrain(1.0, MinimalRPM, NominalRPM - 1.0);
 262
 263   MaximalRPM = ConfigValue(rotor_element, "maxrpm", 2.0*NominalRPM);
 264   MaximalRPM = Constrain(NominalRPM, MaximalRPM, 1e9);
 265
 266   estimate = Constrain(0.07, 2.0/Radius , 0.14); // guess solidity
 267   estimate = estimate * M_PI*Radius/BladeNum;
 268   BladeChord = ConfigValueConv(rotor_element, "chord", estimate, "FT", yell);
 269
 270   LiftCurveSlope = ConfigValue(rotor_element, "liftcurveslope", 6.0); // "1/RAD"
 271   BladeTwist = ConfigValueConv(rotor_element, "twist", -0.17, "RAD");
 272
 273   HingeOffset = ConfigValueConv(rotor_element, "hingeoffset", 0.05 * Radius, "FT" );
 274
 275   estimate = sqr(BladeChord) * sqr(Radius - HingeOffset) * 0.57;
 276   BladeFlappingMoment = ConfigValueConv(rotor_element, "flappingmoment", estimate, "SLUG*FT2");
 277   BladeFlappingMoment = Constrain(1.0e-6, BladeFlappingMoment, 1e9);
 278
 279   // guess mass from moment of a thin stick, and multiply by the blades cg distance
 280   estimate = ( 3.0 * BladeFlappingMoment / sqr(Radius) ) * (0.45 * Radius) ;
 281   BladeMassMoment = ConfigValue(rotor_element, "massmoment", estimate); // unit is slug-ft
 282   BladeMassMoment = Constrain(0.001, BladeMassMoment, 1e9);
 283
 284   estimate = 1.1 * BladeFlappingMoment * BladeNum;
 285   PolarMoment = ConfigValueConv(rotor_element, "polarmoment", estimate, "SLUG*FT2");
 286   PolarMoment = Constrain(1e-6, PolarMoment, 1e9);
 287
 288   // "inflowlag" is treated further down.
 289
 290   TipLossB = ConfigValue(rotor_element, "tiplossfactor", 1.0, silent);
 291
 292   estimate = 0.01 * PolarMoment ; // guesses for huey, bo105 20-30hp
 293   MaxBrakePower  = ConfigValueConv(rotor_element, "maxbrakepower", estimate, "HP");
 294   MaxBrakePower *= hptoftlbssec;
 295
 296   // ground effect
 297   if (rotor_element->FindElement("cgroundeffect")) {
 298     double cge,gee;
 299     cge = rotor_element->FindElementValueAsNumber("cgroundeffect");
 300     cge = Constrain(1e-9, cge, 1.0);
 301     gee = 1.0 / ( 2.0*Radius * cge );
 302     cerr << "# *** 'cgroundeffect' is defunct." << endl;
 303     cerr << "# *** use 'groundeffectexp' with: " << gee << endl;
 304   }
 305
 306   GroundEffectExp = ConfigValue(rotor_element, "groundeffectexp", 0.0);
 307   GroundEffectShift = ConfigValueConv(rotor_element, "groundeffectshift", 0.0, "FT");
 308
 309   // handle optional free-wheeling-unit (FWU)
 310   FreeWheelPresent = 0;
 311   FreeWheelTransmission = 1.0;
 312   if (rotor_element->FindElement("freewheelthresh")) {
 313     FreeWheelThresh = rotor_element->FindElementValueAsNumber("freewheelthresh");
 314     if (FreeWheelThresh > 1.0) {
 315       FreeWheelPresent = 1;
 316       FreeWheelTransmission = 0.0;
 317     }
 318   }
 319
 320   // precalc often used powers
 321   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];
 322   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];
 323
 324   // derived parameters
 325   LockNumberByRho = LiftCurveSlope * BladeChord * R[4] / BladeFlappingMoment;
 326   Solidity = BladeNum * BladeChord / (M_PI * Radius);
 327
 328   // estimate inflow lag, see /GE49/ eqn(1)
 329   double omega_tmp = (NominalRPM/60.0)*2.0*M_PI;
 330   estimate = 16.0/(LockNumberByRho*rho * omega_tmp ); // 16/(gamma*Omega)
 331   // printf("# Est. InflowLag: %f\n", estimate);
 332   InflowLag = ConfigValue(rotor_element, "inflowlag", estimate, yell);
 333   InflowLag = Constrain(1.0e-6, InflowLag, 2.0);
 334
 335   return;
 336 } // Configure
 337
 338 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 339
 340 // calculate control-axes components of total airspeed at the hub.
 341 // sets rotor orientation angle (beta) as side effect. /SH79/ eqn(19-22)
 342
 343 FGColumnVector3 FGRotor::hub_vel_body2ca( const FGColumnVector3 &uvw,
 344                                                  const FGColumnVector3 &pqr,
 345                                                  double a_ic, double b_ic)
 346 {
 347   FGColumnVector3  v_r, v_shaft, v_w;
 348   FGColumnVector3 pos;
 349
 350   pos = fdmex->GetMassBalance()->StructuralToBody(GetActingLocation());
 351
 352   v_r = uvw + pqr*pos;
 353   v_shaft = TboToHsr * InvTransform * v_r;
 354
 355   beta_orient = atan2(v_shaft(eV),v_shaft(eU));
 356
 357   v_w(eU) = v_shaft(eU)*cos(beta_orient) + v_shaft(eV)*sin(beta_orient);
 358   v_w(eV) = 0.0;
 359   v_w(eW) = v_shaft(eW) - b_ic*v_shaft(eU) - a_ic*v_shaft(eV);
 360
 361   return v_w;
 362 }
 363
 364 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 365
 366 // express fuselage angular velocity in control axes /SH79/ eqn(30,31)
 367
 368 FGColumnVector3 FGRotor::fus_angvel_body2ca( const FGColumnVector3 &pqr)
 369 {
 370   FGColumnVector3 av_s_fus, av_w_fus;
 371
 372   // for comparison:
 373   // av_s_fus = BodyToShaft * pqr; /SH79/
 374   // BodyToShaft = TboToHsr * InvTransform
 375   av_s_fus = TboToHsr * InvTransform * pqr;
 376
 377   av_w_fus(eP)=   av_s_fus(eP)*cos(beta_orient) + av_s_fus(eQ)*sin(beta_orient);
 378   av_w_fus(eQ)= - av_s_fus(eP)*sin(beta_orient) + av_s_fus(eQ)*cos(beta_orient);
 379   av_w_fus(eR)=   av_s_fus(eR);
 380
 381   return av_w_fus;
 382 }
 383
 384
 385 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 386
 387 // The calculation is a bit tricky because thrust depends on induced velocity,
 388 // and vice versa.
 389 //
 390 // The flow_scale parameter (ranging from 0.5-1.0) is used to approximate a
 391 // reduction of inflow if the helicopter is close to the ground, yielding to
 392 // higher thrust, see /TA77/ eqn(10a).
 393
 394 void FGRotor::calc_flow_and_thrust( double theta_0, double Uw, double Ww,
 395                                     double flow_scale)
 396 {
 397
 398   double ct_over_sigma = 0.0;
 399   double c0, ct_l, ct_t0, ct_t1;
 400   double mu2;
 401
 402   mu = Uw/(Omega*Radius); // /SH79/ eqn(24)
 403   if (mu > 0.7) mu = 0.7;
 404   mu2 = sqr(mu);
 405
 406   ct_t0 = (1.0/3.0*B[3] + 1.0/2.0 * TipLossB*mu2 - 4.0/(9.0*M_PI) * mu*mu2 ) * theta_0;
 407   ct_t1 = (1.0/4.0*B[4] + 1.0/4.0 * B[2]*mu2) * BladeTwist;
 408
 409   ct_l  = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda; // first time
 410
 411   c0 = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1) * Solidity;
 412   c0 = c0 / ( 2.0 * sqrt( sqr(mu) + sqr(lambda) ) + 1e-15);
 413
 414   // replacement for /SH79/ eqn(26).
 415   // ref: dnu/dt = 1/tau ( Ct / (2*sqrt(mu^2+lambda^2))  -  nu )
 416   // taking mu and lambda constant, this integrates to
 417
 418   nu  = flow_scale * ((nu - c0) * exp(-dt/InflowLag) + c0);
 419
 420   // now from nu to lambda, C_T, and Thrust
 421
 422   lambda = Ww/(Omega*Radius) - nu; // /SH79/ eqn(25)
 423
 424   ct_l  = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda;
 425
 426   ct_over_sigma = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1); // /SH79/ eqn(27)
 427
 428   Thrust = BladeNum*BladeChord*Radius*rho*sqr(Omega*Radius) * ct_over_sigma;
 429
 430   C_T = ct_over_sigma * Solidity;
 431   v_induced = nu * (Omega*Radius);
 432
 433 }
 434
 435
 436 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 437
 438 // The coning angle doesn't apply for teetering rotors, but calculating
 439 // doesn't hurt. /SH79/ eqn(29)
 440
 441 void FGRotor::calc_coning_angle(double theta_0)
 442 {
 443   double lock_gamma = LockNumberByRho * rho;
 444
 445   double a0_l  = (1.0/6.0  + 0.04 * mu*mu*mu) * lambda;
 446   double a0_t0 = (1.0/8.0  + 1.0/8.0  * mu*mu) * theta_0;
 447   double a0_t1 = (1.0/10.0 + 1.0/12.0 * mu*mu) * BladeTwist;
 448   a0 = lock_gamma * ( a0_l + a0_t0 + a0_t1);
 449   return;
 450 }
 451
 452 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 453
 454 // Flapping angles relative to control axes /SH79/ eqn(32)
 455
 456 void FGRotor::calc_flapping_angles(double theta_0, const FGColumnVector3 &pqr_fus_w)
 457 {
 458   double lock_gamma = LockNumberByRho * rho;
 459
 460
 461   double mu2_2 = sqr(mu)/2.0;
 462   double t075 = theta_0 + 0.75 * BladeTwist;  // common approximation for rectangular blades
 463
 464   a_1 = 1.0/(1.0 - mu2_2) * (
 465                                  (2.0*lambda + (8.0/3.0)*t075)*mu
 466                                + pqr_fus_w(eP)/Omega
 467                                - 16.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
 468                              );
 469
 470   b_1 = 1.0/(1.0 + mu2_2) * (
 471                                  (4.0/3.0)*mu*a0
 472                                - pqr_fus_w(eQ)/Omega
 473                                - 16.0 * pqr_fus_w(eP)/(lock_gamma*Omega)
 474                              );
 475
 476   // used in  force calc
 477   a_dw = 1.0/(1.0 - mu2_2) * (
 478                                  (2.0*lambda + (8.0/3.0)*t075)*mu
 479                                - 24.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
 480                                  * ( 1.0 - ( 0.29 * t075 / (C_T/Solidity) ) )
 481                              );
 482
 483   return;
 484 }
 485
 486 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 487
 488 // /SH79/ eqn(38,39)
 489
 490 void FGRotor::calc_drag_and_side_forces(double theta_0)
 491 {
 492   double cy_over_sigma;
 493   double t075 = theta_0 + 0.75 * BladeTwist;
 494
 495   H_drag = Thrust * a_dw;
 496
 497   cy_over_sigma = (
 498                       0.75*b_1*lambda - 1.5*a0*mu*lambda + 0.25*a_1*b_1*mu
 499                     - a0*a_1*sqr(mu) + (1.0/6.0)*a0*a_1
 500                     - (0.75*mu*a0 - (1.0/3.0)*b_1 - 0.5*sqr(mu)*b_1)*t075
 501                   );
 502   cy_over_sigma *= LiftCurveSlope/2.0;
 503
 504   J_side = BladeNum * BladeChord * Radius * rho * sqr(Omega*Radius) * cy_over_sigma;
 505
 506   return;
 507 }
 508
 509 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 510
 511 // Simplified version of /SH79/ eqn(36). Uses an estimate for blade drag
 512 // (a new config parameter to come...).
 513 // From "Bramwell's Helicopter Dynamics", second edition, eqn(3.43) and (3.44)
 514
 515 void FGRotor::calc_torque(double theta_0)
 516 {
 517   // estimate blade drag
 518   double delta_dr = 0.009 + 0.3*sqr(6.0*C_T/(LiftCurveSlope*Solidity));
 519
 520   Torque = rho * BladeNum * BladeChord * delta_dr * sqr(Omega*Radius) * R[2] *
 521            (1.0+4.5*sqr(mu))/8.0
 522                      - (Thrust*lambda + H_drag*mu)*Radius;
 523
 524   return;
 525 }
 526
 527 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 528
 529 // transform rotor forces from control axes to shaft axes, and express
 530 // in body axes /SH79/ eqn(40,41)
 531
 532 FGColumnVector3 FGRotor::body_forces(double a_ic, double b_ic)
 533 {
 534   FGColumnVector3 F_s(
 535         - H_drag*cos(beta_orient) - J_side*sin(beta_orient) + Thrust*b_ic,
 536         - H_drag*sin(beta_orient) + J_side*cos(beta_orient) + Thrust*a_ic,
 537         - Thrust);
 538
 539   return HsrToTbo * F_s;
 540 }
 541
 542 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 543
 544 // calculates the additional moments due to hinge offset and handles
 545 // torque and sense
 546
 547 FGColumnVector3 FGRotor::body_moments(double a_ic, double b_ic)
 548 {
 549   FGColumnVector3 M_s, M_hub, M_h;
 550   double mf;
 551
 552   // cyclic flapping relative to shaft axes /SH79/ eqn(43)
 553   a1s = a_1*cos(beta_orient) + b_1*sin(beta_orient) - b_ic;
 554   b1s = b_1*cos(beta_orient) - a_1*sin(beta_orient) + a_ic;
 555
 556   mf = 0.5 * HingeOffset * BladeNum * Omega*Omega * BladeMassMoment;
 557
 558   M_s(eL) = mf*b1s;
 559   M_s(eM) = mf*a1s;
 560   M_s(eN) = Torque * Sense ;
 561
 562   return HsrToTbo * M_s;
 563 }
 564
 565 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 566
 567 void FGRotor::CalcStatePart1(void)
 568 {
 569   double A_IC;       // lateral (roll) control in radians
 570   double B_IC;       // longitudinal (pitch) control in radians
 571   double theta_col;  // rotor collective pitch in radians
 572
 573   FGColumnVector3 vHub_ca, avFus_ca;
 574
 575   double filtered_hagl = 0.0;
 576   double ge_factor = 1.0;
 577
 578   // fetch needed values from environment
 579   rho = in.Density; // slugs/ft^3.
 580   double h_agl_ft = in.H_agl;
 581   // update InvTransform, the rotor orientation could have been altered
 582   InvTransform = Transform().Transposed();
 583
 584   // handle RPM requirements, calc omega.
 585   if (ExternalRPM && ExtRPMsource) {
 586     RPM = ExtRPMsource->getDoubleValue() / GearRatio;
 587   }
 588
 589   // MinimalRPM is always >= 1. MaximalRPM is always >= NominalRPM
 590   RPM = Constrain(MinimalRPM, RPM, MaximalRPM);
 591
 592   Omega = (RPM/60.0)*2.0*M_PI;
 593
 594   // set control inputs
 595   A_IC      = LateralCtrl;
 596   B_IC      = LongitudinalCtrl;
 597   theta_col = CollectiveCtrl;
 598
 599   // ground effect
 600   if (GroundEffectExp > 1e-5) {
 601     if (h_agl_ft<0.0) h_agl_ft = 0.0; // clamp
 602     filtered_hagl = damp_hagl.execute(h_agl_ft) + GroundEffectShift;
 603     // actual/nominal factor avoids absurd scales at startup
 604     ge_factor -= exp(-filtered_hagl*GroundEffectExp) * (RPM / NominalRPM);
 605     if (ge_factor<0.5) ge_factor=0.5; // clamp
 606   }
 607
 608   // all set, start calculations
 609
 610   vHub_ca  = hub_vel_body2ca(in.AeroUVW, in.AeroPQR, A_IC, B_IC);
 611
 612   avFus_ca = fus_angvel_body2ca(in.AeroPQR);
 613
 614   calc_flow_and_thrust(theta_col, vHub_ca(eU), vHub_ca(eW), ge_factor);
 615
 616   calc_coning_angle(theta_col);
 617
 618   calc_flapping_angles(theta_col, avFus_ca);
 619
 620   calc_drag_and_side_forces(theta_col);
 621
 622   calc_torque(theta_col);
 623
 624   // Fixme: only valid for a 'decent' rotor
 625   theta_downwash = atan2( -in.AeroUVW(eU), v_induced - in.AeroUVW(eW));
 626   phi_downwash   = atan2(  in.AeroUVW(eV), v_induced - in.AeroUVW(eW));
 627
 628   vFn = body_forces(A_IC, B_IC);
 629   vMn = Transform() * body_moments(A_IC, B_IC);
 630
 631 }
 632
 633 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 634
 635 void FGRotor::CalcStatePart2(double PowerAvailable)
 636 {
 637   if (! ExternalRPM) {
 638     // calculate new RPM
 639     double ExcessTorque = PowerAvailable / Omega;
 640     double deltaOmega   = ExcessTorque / PolarMoment * in.TotalDeltaT;
 641     RPM += deltaOmega/(2.0*M_PI) * 60.0;
 642   }
 643   RPM = Constrain(MinimalRPM, RPM, MaximalRPM); // trim again
 644 }
 645
 646 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 647
 648 // Simulation of a free-wheeling-unit (FWU). Might need improvements.
 649
 650 void FGRotor::calc_freewheel_state(double p_source, double p_load) {
 651
 652   // engine is off/detached, release.
 653   if (p_source<1e-3) {
 654     FreeWheelTransmission = 0.0;
 655     return;
 656   }
 657
 658   // engine is driving the rotor, engage.
 659   if (p_source >= p_load) {
 660     FreeWheelTransmission = 1.0;
 661     return;
 662   }
 663
 664   // releases if engine is detached, but stays calm if
 665   // the load changes due to rotor dynamics.
 666   if (p_source > 0.0 && p_load/(p_source+0.1) > FreeWheelThresh ) {
 667     FreeWheelTransmission = 0.0;
 668     return;
 669   }
 670
 671   return;
 672 }
 673
 674 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 675
 676 double FGRotor::Calculate(double EnginePower)
 677 {
 678   double FWmult = 1.0;
 679   double DeltaPower;
 680
 681   CalcStatePart1();
 682
 683   PowerRequired = Torque * Omega + BrakeCtrlNorm * MaxBrakePower;
 684
 685   if (FreeWheelPresent) {
 686     calc_freewheel_state(EnginePower * ClutchCtrlNorm, PowerRequired);
 687     FWmult = FreeWheelLag.execute(FreeWheelTransmission);
 688   }
 689
 690   DeltaPower = EnginePower * ClutchCtrlNorm * FWmult - PowerRequired;
 691
 692   CalcStatePart2(DeltaPower);
 693
 694   return Thrust;
 695 }
 696
 697
 698 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 699
 700
 701 bool FGRotor::BindModel(void)
 702 {
 703   string property_name, base_property_name;
 704   base_property_name = CreateIndexedPropertyName("propulsion/engine", EngineNum);
 705
 706   property_name = base_property_name + "/rotor-rpm";
 707   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetRPM );
 708
 709   property_name = base_property_name + "/x-engine-rpm"; // used for RPM eXchange
 710   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetEngineRPM );
 711
 712   property_name = base_property_name + "/rotor-thrust-lbs"; // might be redundant - check!
 713   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetThrust );
 714
 715   property_name = base_property_name + "/a0-rad";
 716   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA0 );
 717
 718   property_name = base_property_name + "/a1-rad";
 719   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA1 );
 720
 721   property_name = base_property_name + "/b1-rad";
 722   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetB1 );
 723
 724   property_name = base_property_name + "/inflow-ratio";
 725   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLambda );
 726
 727   property_name = base_property_name + "/advance-ratio";
 728   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetMu );
 729
 730   property_name = base_property_name + "/induced-inflow-ratio";
 731   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetNu );
 732
 733   property_name = base_property_name + "/vi-fps";
 734   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetVi );
 735
 736   property_name = base_property_name + "/thrust-coefficient";
 737   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCT );
 738
 739   property_name = base_property_name + "/torque-lbsft";
 740   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetTorque );
 741
 742   property_name = base_property_name + "/theta-downwash-rad";
 743   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetThetaDW );
 744
 745   property_name = base_property_name + "/phi-downwash-rad";
 746   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetPhiDW );
 747
 748   switch (ControlMap) {
 749     case eTailCtrl:
 750       property_name = base_property_name + "/antitorque-ctrl-rad";
 751       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
 752       break;
 753     case eTandemCtrl:
 754       property_name = base_property_name + "/tail-collective-ctrl-rad";
 755       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
 756       property_name = base_property_name + "/lateral-ctrl-rad";
 757       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
 758       property_name = base_property_name + "/longitudinal-ctrl-rad";
 759       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
 760       break;
 761     default: // eMainCtrl
 762       property_name = base_property_name + "/collective-ctrl-rad";
 763       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
 764       property_name = base_property_name + "/lateral-ctrl-rad";
 765       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
 766       property_name = base_property_name + "/longitudinal-ctrl-rad";
 767       PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
 768   }
 769
 770   property_name = base_property_name + "/brake-ctrl-norm";
 771   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetBrakeCtrl, &FGRotor::SetBrakeCtrl);
 772   property_name = base_property_name + "/free-wheel-transmission";
 773   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetFreeWheelTransmission);
 774
 775   if (ExternalRPM) {
 776     if (RPMdefinition == -1) {
 777       property_name = base_property_name + "/x-rpm-dict";
 778       ExtRPMsource = PropertyManager->GetNode(property_name, true);
 779     } else if (RPMdefinition >= 0 && RPMdefinition != EngineNum) {
 780       string ipn = CreateIndexedPropertyName("propulsion/engine", RPMdefinition);
 781       property_name = ipn + "/x-engine-rpm";
 782       ExtRPMsource = PropertyManager->GetNode(property_name, false);
 783       if (! ExtRPMsource) {
 784         cerr << "# Warning: Engine number " << EngineNum << "." << endl;
 785         cerr << "# No 'x-engine-rpm' property found for engine " << RPMdefinition << "." << endl;
 786         cerr << "# Please check order of engine definitons."  << endl;
 787       }
 788     } else {
 789       cerr << "# Engine number " << EngineNum;
 790       cerr << ", given ExternalRPM value '" << RPMdefinition << "' unhandled."  << endl;
 791     }
 792   }
 793
 794   return true;
 795 }
 796
 797 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 798
 799 string FGRotor::GetThrusterLabels(int id, string delimeter)
 800 {
 801
 802   ostringstream buf;
 803
 804   buf << Name << " RPM (engine " << id << ")";
 805
 806   return buf.str();
 807
 808 }
 809
 810 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 811
 812 string FGRotor::GetThrusterValues(int id, string delimeter)
 813 {
 814
 815   ostringstream buf;
 816
 817   buf << RPM;
 818
 819   return buf.str();
 820
 821 }
 822
 823 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 824 //    The bitmasked value choices are as follows:
 825 //    unset: In this case (the default) JSBSim would only print
 826 //       out the normally expected messages, essentially echoing
 827 //       the config files as they are read. If the environment
 828 //       variable is not set, debug_lvl is set to 1 internally
 829 //    0: This requests JSBSim not to output any messages
 830 //       whatsoever.
 831 //    1: This value explicity requests the normal JSBSim
 832 //       startup messages
 833 //    2: This value asks for a message to be printed out when
 834 //       a class is instantiated
 835 //    4: When this value is set, a message is displayed when a
 836 //       FGModel object executes its Run() method
 837 //    8: When this value is set, various runtime state variables
 838 //       are printed out periodically
 839 //    16: When set various parameters are sanity checked and
 840 //       a message is printed out when they go out of bounds
 841
 842 void FGRotor::Debug(int from)
 843 {
 844   string ControlMapName;
 845
 846   if (debug_lvl <= 0) return;
 847
 848   if (debug_lvl & 1) { // Standard console startup message output
 849     if (from == 0) { // Constructor
 850       cout << "\n    Rotor Name: " << Name << endl;
 851       cout << "      Diameter = " << 2.0 * Radius << " ft." << endl;
 852       cout << "      Number of Blades = " << BladeNum << endl;
 853       cout << "      Gear Ratio = " << GearRatio << endl;
 854       cout << "      Sense = " << Sense << endl;
 855       cout << "      Nominal RPM = " << NominalRPM << endl;
 856       cout << "      Minimal RPM = " << MinimalRPM << endl;
 857       cout << "      Maximal RPM = " << MaximalRPM << endl;
 858
 859       if (ExternalRPM) {
 860         if (RPMdefinition == -1) {
 861           cout << "      RPM is controlled externally" << endl;
 862         } else {
 863           cout << "      RPM source set to engine " << RPMdefinition << endl;
 864         }
 865       }
 866
 867       cout << "      Blade Chord = " << BladeChord << endl;
 868       cout << "      Lift Curve Slope = " << LiftCurveSlope << endl;
 869       cout << "      Blade Twist = " << BladeTwist << endl;
 870       cout << "      Hinge Offset = " << HingeOffset << endl;
 871       cout << "      Blade Flapping Moment = " << BladeFlappingMoment << endl;
 872       cout << "      Blade Mass Moment = " << BladeMassMoment << endl;
 873       cout << "      Polar Moment = " << PolarMoment << endl;
 874       cout << "      Inflow Lag = " << InflowLag << endl;
 875       cout << "      Tip Loss = " << TipLossB << endl;
 876       cout << "      Lock Number = " << LockNumberByRho * 0.002356 << " (SL)" << endl;
 877       cout << "      Solidity = " << Solidity << endl;
 878       cout << "      Max Brake Power = " << MaxBrakePower/hptoftlbssec << " HP" << endl;
 879
 880       switch (ControlMap) {
 881         case eTailCtrl:    ControlMapName = "Tail Rotor";   break;
 882         case eTandemCtrl:  ControlMapName = "Tandem Rotor"; break;
 883         default:           ControlMapName = "Main Rotor";
 884       }
 885       cout << "      Control Mapping = " << ControlMapName << endl;
 886
 887       if (FreeWheelPresent) {
 888         cout << "      Free Wheel Threshold = " << FreeWheelThresh << endl;
 889       } else {
 890         cout << "      No FWU present" << endl;
 891       }
 892
 893     }
 894   }
 895   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 896     if (from == 0) cout << "Instantiated: FGRotor" << endl;
 897     if (from == 1) cout << "Destroyed:    FGRotor" << endl;
 898   }
 899   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 900   }
 901   if (debug_lvl & 8 ) { // Runtime state variables
 902   }
 903   if (debug_lvl & 16) { // Sanity checking
 904   }
 905   if (debug_lvl & 64) {
 906     if (from == 0) { // Constructor
 907       cout << IdSrc << endl;
 908       cout << IdHdr << endl;
 909     }
 910   }
 911
 912 }
 913
 914
 915 } // namespace JSBSim
 916