JSBSim Flight Dynamics Model  1.0 (02 March 2017)
An Open Source Flight Dynamics and Control Software Library in C++
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 <string>
47 #include <sstream>
48 
49 #include "FGRotor.h"
50 #include "models/FGMassBalance.h"
51 #include "models/FGPropulsion.h" // to get the GearRatio from a linked rotor
52 #include "input_output/FGXMLElement.h"
53 
54 using std::cerr;
55 using std::cout;
56 using std::endl;
57 using std::string;
58 using std::ostringstream;
59 
60 namespace JSBSim {
61 
62 IDENT(IdSrc,"$Id: FGRotor.cpp,v 1.25 2015/09/28 08:57:32 ehofman Exp $");
63 IDENT(IdHdr,ID_ROTOR);
64 
65 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
66 MISC
67 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
68 
69 static inline double sqr(double x) { return x*x; }
70 
71 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
72 CLASS IMPLEMENTATION
73 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
74 
75 // Constructor
76 
77 FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
78  : FGThruster(exec, rotor_element, num),
79  rho(0.002356), // environment
80  Radius(0.0), BladeNum(0), // configuration parameters
81  Sense(1.0), NominalRPM(0.0), MinimalRPM(0.0), MaximalRPM(0.0),
82  ExternalRPM(0), RPMdefinition(0), ExtRPMsource(NULL), SourceGearRatio(1.0),
83  BladeChord(0.0), LiftCurveSlope(0.0), BladeTwist(0.0), HingeOffset(0.0),
84  BladeFlappingMoment(0.0), BladeMassMoment(0.0), PolarMoment(0.0),
85  InflowLag(0.0), TipLossB(0.0),
86  GroundEffectExp(0.0), GroundEffectShift(0.0), GroundEffectScaleNorm(1.0),
87  LockNumberByRho(0.0), Solidity(0.0), // derived parameters
88  RPM(0.0), Omega(0.0), // dynamic values
89  beta_orient(0.0),
90  a0(0.0), a_1(0.0), b_1(0.0), a_dw(0.0),
91  a1s(0.0), b1s(0.0),
92  H_drag(0.0), J_side(0.0), Torque(0.0), C_T(0.0),
93  lambda(-0.001), mu(0.0), nu(0.001), v_induced(0.0),
94  theta_downwash(0.0), phi_downwash(0.0),
95  ControlMap(eMainCtrl), // control
96  CollectiveCtrl(0.0), LateralCtrl(0.0), LongitudinalCtrl(0.0),
97  Transmission(NULL), // interaction with engine
98  EngineRPM(0.0), MaxBrakePower(0.0), GearLoss(0.0), GearMoment(0.0)
99 {
100  FGColumnVector3 location(0.0, 0.0, 0.0), orientation(0.0, 0.0, 0.0);
101  Element *thruster_element;
102  double engine_power_est = 0.0;
103 
104  // initialise/set remaining variables
105  SetTransformType(FGForce::tCustom);
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(); // body to custom/native
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  SourceGearRatio = 1.0;
162  RPMdefinition = (int) rotor_element->FindElementValueAsNumber("ExternalRPM");
163  int rdef = RPMdefinition;
164  if (RPMdefinition>=0) {
165  // avoid ourself and (still) unknown engines.
166  if (!exec->GetPropulsion()->GetEngine(RPMdefinition) || RPMdefinition==num) {
167  RPMdefinition = -1;
168  } else {
169  FGThruster *tr = exec->GetPropulsion()->GetEngine(RPMdefinition)->GetThruster();
170  SourceGearRatio = tr->GetGearRatio();
171  // cout << "# got sources' GearRatio: " << SourceGearRatio << endl;
172  }
173  }
174  if (RPMdefinition != rdef) {
175  cerr << "# discarded given RPM source (" << rdef << ") and switched to external control (-1)." << endl;
176  }
177  }
178 
179  // process rotor parameters
180  engine_power_est = Configure(rotor_element);
181 
182  // setup transmission if needed
183  if (!ExternalRPM) {
184 
185  Transmission = new FGTransmission(exec, num, dt);
186 
187  Transmission->SetThrusterMoment(PolarMoment);
188 
189  // The MOI sensed behind the gear ( MOI_engine*sqr(GearRatio) ).
190  GearMoment = ConfigValueConv(rotor_element, "gearmoment", 0.1*PolarMoment, "SLUG*FT2");
191  GearMoment = Constrain(1e-6, GearMoment, 1e9);
192  Transmission->SetEngineMoment(GearMoment);
193 
194  Transmission->SetMaxBrakePower(MaxBrakePower);
195 
196  GearLoss = ConfigValueConv(rotor_element, "gearloss", 0.0025 * engine_power_est, "HP");
197  GearLoss = Constrain(0.0, GearLoss, 1e9);
198  GearLoss *= hptoftlbssec;
199  Transmission->SetEngineFriction(GearLoss);
200 
201  }
202 
203  // shaft representation - a rather simple transform,
204  // but using a matrix is safer.
205  TboToHsr.InitMatrix( 0.0, 0.0, 1.0,
206  0.0, 1.0, 0.0,
207  -1.0, 0.0, 0.0 );
208  HsrToTbo = TboToHsr.Transposed();
209 
210  // smooth out jumps in hagl reported, otherwise the ground effect
211  // calculation would cause jumps too. 1Hz seems sufficient.
212  damp_hagl = Filter(1.0, dt);
213 
214  // enable import-export
215  bindmodel(exec->GetPropertyManager());
216 
217  Debug(0);
218 
219 } // Constructor
220 
221 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
222 
224  if (Transmission) delete Transmission;
225  Debug(1);
226 }
227 
228 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
229 
230 // 5in1: value-fetch-convert-default-return function
231 
232 double FGRotor::ConfigValueConv( Element* el, const string& ename, double default_val,
233  const string& unit, bool tell)
234 {
235 
236  Element *e = NULL;
237  double val = default_val;
238 
239  string pname = "*No parent element*";
240 
241  if (el) {
242  e = el->FindElement(ename);
243  pname = el->GetName();
244  }
245 
246  if (e) {
247  if (unit.empty()) {
248  val = e->GetDataAsNumber();
249  } else {
250  val = el->FindElementValueAsNumberConvertTo(ename,unit);
251  }
252  } else {
253  if (tell) {
254  cerr << pname << ": missing element '" << ename <<
255  "' using estimated value: " << default_val << endl;
256  }
257  }
258 
259  return val;
260 }
261 
262 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
263 
264 double FGRotor::ConfigValue(Element* el, const string& ename, double default_val, bool tell)
265 {
266  return ConfigValueConv(el, ename, default_val, "", tell);
267 }
268 
269 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
270 
271 // 1. read configuration and try to fill holes, ymmv
272 // 2. calculate derived parameters
273 double FGRotor::Configure(Element* rotor_element)
274 {
275 
276  double estimate, engine_power_est=0.0;
277  const bool yell = true;
278  const bool silent = false;
279 
280 
281  Radius = 0.5 * ConfigValueConv(rotor_element, "diameter", 42.0, "FT", yell);
282  Radius = Constrain(1e-3, Radius, 1e9);
283 
284  BladeNum = (int) ConfigValue(rotor_element, "numblades", 3 , yell);
285 
286  GearRatio = ConfigValue(rotor_element, "gearratio", 1.0, yell);
287  GearRatio = Constrain(1e-9, GearRatio, 1e9);
288 
289  // make sure that v_tip (omega*r) is below 0.7mach ~ 750ft/s
290  estimate = (750.0/Radius)/(2.0*M_PI) * 60.0; // 7160/Radius
291  NominalRPM = ConfigValue(rotor_element, "nominalrpm", estimate, yell);
292  NominalRPM = Constrain(2.0, NominalRPM, 1e9);
293 
294  MinimalRPM = ConfigValue(rotor_element, "minrpm", 1.0);
295  MinimalRPM = Constrain(1.0, MinimalRPM, NominalRPM - 1.0);
296 
297  MaximalRPM = ConfigValue(rotor_element, "maxrpm", 2.0*NominalRPM);
298  MaximalRPM = Constrain(NominalRPM, MaximalRPM, 1e9);
299 
300  estimate = Constrain(0.07, 2.0/Radius , 0.14); // guess solidity
301  estimate = estimate * M_PI*Radius/BladeNum;
302  BladeChord = ConfigValueConv(rotor_element, "chord", estimate, "FT", yell);
303 
304  LiftCurveSlope = ConfigValue(rotor_element, "liftcurveslope", 6.0); // "1/RAD"
305  BladeTwist = ConfigValueConv(rotor_element, "twist", -0.17, "RAD");
306 
307  HingeOffset = ConfigValueConv(rotor_element, "hingeoffset", 0.05 * Radius, "FT" );
308 
309  estimate = sqr(BladeChord) * sqr(Radius - HingeOffset) * 0.57;
310  BladeFlappingMoment = ConfigValueConv(rotor_element, "flappingmoment", estimate, "SLUG*FT2");
311  BladeFlappingMoment = Constrain(1e-9, BladeFlappingMoment, 1e9);
312 
313  // guess mass from moment of a thin stick, and multiply by the blades cg distance
314  estimate = ( 3.0 * BladeFlappingMoment / sqr(Radius) ) * (0.45 * Radius) ;
315  BladeMassMoment = ConfigValue(rotor_element, "massmoment", estimate); // unit is slug-ft
316  BladeMassMoment = Constrain(1e-9, BladeMassMoment, 1e9);
317 
318  estimate = 1.1 * BladeFlappingMoment * BladeNum;
319  PolarMoment = ConfigValueConv(rotor_element, "polarmoment", estimate, "SLUG*FT2");
320  PolarMoment = Constrain(1e-9, PolarMoment, 1e9);
321 
322  // "inflowlag" is treated further down.
323 
324  TipLossB = ConfigValue(rotor_element, "tiplossfactor", 1.0, silent);
325 
326  // estimate engine power (bit of a pity, cause our caller already knows)
327  engine_power_est = 0.5 * BladeNum*BladeChord*Radius*Radius;
328 
329  estimate = engine_power_est / 30.0;
330  MaxBrakePower = ConfigValueConv(rotor_element, "maxbrakepower", estimate, "HP");
331  MaxBrakePower *= hptoftlbssec;
332 
333  GroundEffectExp = ConfigValue(rotor_element, "groundeffectexp", 0.0);
334  GroundEffectShift = ConfigValueConv(rotor_element, "groundeffectshift", 0.0, "FT");
335 
336  // precalc often used powers
337  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];
338  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];
339 
340  // derived parameters
341  LockNumberByRho = LiftCurveSlope * BladeChord * R[4] / BladeFlappingMoment;
342  Solidity = BladeNum * BladeChord / (M_PI * Radius);
343 
344  // estimate inflow lag, see /GE49/ eqn(1)
345  double omega_tmp = (NominalRPM/60.0)*2.0*M_PI;
346  estimate = 16.0/(LockNumberByRho*rho * omega_tmp ); // 16/(gamma*Omega)
347  // printf("# Est. InflowLag: %f\n", estimate);
348  InflowLag = ConfigValue(rotor_element, "inflowlag", estimate, yell);
349  InflowLag = Constrain(1e-6, InflowLag, 2.0);
350 
351  return engine_power_est;
352 } // Configure
353 
354 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
355 
356 // calculate control-axes components of total airspeed at the hub.
357 // sets rotor orientation angle (beta) as side effect. /SH79/ eqn(19-22)
358 
359 FGColumnVector3 FGRotor::hub_vel_body2ca( const FGColumnVector3 &uvw,
360  const FGColumnVector3 &pqr,
361  double a_ic, double b_ic)
362 {
363  FGColumnVector3 v_r, v_shaft, v_w;
364  FGColumnVector3 pos;
365 
366  pos = fdmex->GetMassBalance()->StructuralToBody(GetActingLocation());
367 
368  v_r = uvw + pqr*pos;
369  v_shaft = TboToHsr * InvTransform * v_r;
370 
371  beta_orient = atan2(v_shaft(eV),v_shaft(eU));
372 
373  v_w(eU) = v_shaft(eU)*cos(beta_orient) + v_shaft(eV)*sin(beta_orient);
374  v_w(eV) = 0.0;
375  v_w(eW) = v_shaft(eW) - b_ic*v_shaft(eU) - a_ic*v_shaft(eV);
376 
377  return v_w;
378 }
379 
380 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
381 
382 // express fuselage angular velocity in control axes /SH79/ eqn(30,31)
383 
384 FGColumnVector3 FGRotor::fus_angvel_body2ca( const FGColumnVector3 &pqr)
385 {
386  FGColumnVector3 av_s_fus, av_w_fus;
387 
388  // for comparison:
389  // av_s_fus = BodyToShaft * pqr; /SH79/
390  // BodyToShaft = TboToHsr * InvTransform
391  av_s_fus = TboToHsr * InvTransform * pqr;
392 
393  av_w_fus(eP)= av_s_fus(eP)*cos(beta_orient) + av_s_fus(eQ)*sin(beta_orient);
394  av_w_fus(eQ)= - av_s_fus(eP)*sin(beta_orient) + av_s_fus(eQ)*cos(beta_orient);
395  av_w_fus(eR)= av_s_fus(eR);
396 
397  return av_w_fus;
398 }
399 
400 
401 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
402 
403 // The calculation is a bit tricky because thrust depends on induced velocity,
404 // and vice versa.
405 //
406 // The flow_scale parameter (ranging from 0.5-1.0) is used to approximate a
407 // reduction of inflow if the helicopter is close to the ground, yielding to
408 // higher thrust, see /TA77/ eqn(10a).
409 
410 void FGRotor::calc_flow_and_thrust( double theta_0, double Uw, double Ww,
411  double flow_scale)
412 {
413 
414  double ct_over_sigma = 0.0;
415  double c0, ct_l, ct_t0, ct_t1;
416  double mu2;
417 
418  mu = Uw/(Omega*Radius); // /SH79/ eqn(24)
419  if (mu > 0.7) mu = 0.7;
420  mu2 = sqr(mu);
421 
422  ct_t0 = (1.0/3.0*B[3] + 1.0/2.0 * TipLossB*mu2 - 4.0/(9.0*M_PI) * mu*mu2 ) * theta_0;
423  ct_t1 = (1.0/4.0*B[4] + 1.0/4.0 * B[2]*mu2) * BladeTwist;
424 
425  ct_l = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda; // first time
426 
427  c0 = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1) * Solidity;
428  c0 = c0 / ( 2.0 * sqrt( sqr(mu) + sqr(lambda) ) + 1e-15);
429 
430  // replacement for /SH79/ eqn(26).
431  // ref: dnu/dt = 1/tau ( Ct / (2*sqrt(mu^2+lambda^2)) - nu )
432  // taking mu and lambda constant, this integrates to
433 
434  nu = flow_scale * ((nu - c0) * exp(-dt/InflowLag) + c0);
435 
436  // now from nu to lambda, C_T, and Thrust
437 
438  lambda = Ww/(Omega*Radius) - nu; // /SH79/ eqn(25)
439 
440  ct_l = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda;
441 
442  ct_over_sigma = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1); // /SH79/ eqn(27)
443 
444  Thrust = BladeNum*BladeChord*Radius*rho*sqr(Omega*Radius) * ct_over_sigma;
445 
446  C_T = ct_over_sigma * Solidity;
447  v_induced = nu * (Omega*Radius);
448 
449 }
450 
451 
452 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
453 
454 // Two blade teetering rotors are often 'preconed' to a fixed angle, but the
455 // calculated value is pretty close to the real one. /SH79/ eqn(29)
456 
457 void FGRotor::calc_coning_angle(double theta_0)
458 {
459  double lock_gamma = LockNumberByRho * rho;
460 
461  double a0_l = (1.0/6.0 + 0.04 * mu*mu*mu) * lambda;
462  double a0_t0 = (1.0/8.0 + 1.0/8.0 * mu*mu) * theta_0;
463  double a0_t1 = (1.0/10.0 + 1.0/12.0 * mu*mu) * BladeTwist;
464  a0 = lock_gamma * ( a0_l + a0_t0 + a0_t1);
465  return;
466 }
467 
468 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
469 
470 // Flapping angles relative to control axes /SH79/ eqn(32)
471 
472 void FGRotor::calc_flapping_angles(double theta_0, const FGColumnVector3 &pqr_fus_w)
473 {
474  double lock_gamma = LockNumberByRho * rho;
475 
476 
477  double mu2_2 = sqr(mu)/2.0;
478  double t075 = theta_0 + 0.75 * BladeTwist; // common approximation for rectangular blades
479 
480  a_1 = 1.0/(1.0 - mu2_2) * (
481  (2.0*lambda + (8.0/3.0)*t075)*mu
482  + pqr_fus_w(eP)/Omega
483  - 16.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
484  );
485 
486  b_1 = 1.0/(1.0 + mu2_2) * (
487  (4.0/3.0)*mu*a0
488  - pqr_fus_w(eQ)/Omega
489  - 16.0 * pqr_fus_w(eP)/(lock_gamma*Omega)
490  );
491 
492  // used in force calc
493  a_dw = 1.0/(1.0 - mu2_2) * (
494  (2.0*lambda + (8.0/3.0)*t075)*mu
495  - 24.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
496  * ( 1.0 - ( 0.29 * t075 / (C_T/Solidity) ) )
497  );
498 
499  return;
500 }
501 
502 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
503 
504 // /SH79/ eqn(38,39)
505 
506 void FGRotor::calc_drag_and_side_forces(double theta_0)
507 {
508  double cy_over_sigma;
509  double t075 = theta_0 + 0.75 * BladeTwist;
510 
511  H_drag = Thrust * a_dw;
512 
513  cy_over_sigma = (
514  0.75*b_1*lambda - 1.5*a0*mu*lambda + 0.25*a_1*b_1*mu
515  - a0*a_1*sqr(mu) + (1.0/6.0)*a0*a_1
516  - (0.75*mu*a0 - (1.0/3.0)*b_1 - 0.5*sqr(mu)*b_1)*t075
517  );
518  cy_over_sigma *= LiftCurveSlope/2.0;
519 
520  J_side = BladeNum * BladeChord * Radius * rho * sqr(Omega*Radius) * cy_over_sigma;
521 
522  return;
523 }
524 
525 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
526 
527 // Simplified version of /SH79/ eqn(36). Uses an estimate for blade drag
528 // (a new config parameter to come...).
529 // From "Bramwell's Helicopter Dynamics", second edition, eqn(3.43) and (3.44)
530 
531 void FGRotor::calc_torque(double theta_0)
532 {
533  // estimate blade drag
534  double delta_dr = 0.009 + 0.3*sqr(6.0*C_T/(LiftCurveSlope*Solidity));
535 
536  Torque = rho * BladeNum * BladeChord * delta_dr * sqr(Omega*Radius) * R[2] *
537  (1.0+4.5*sqr(mu))/8.0
538  - (Thrust*lambda + H_drag*mu)*Radius;
539 
540  return;
541 }
542 
543 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
544 
545 // Get the downwash angles with respect to the shaft axis.
546 // Given a 'regular' main rotor, the angles are zero when the downwash points
547 // down, positive theta values mean that the downwash turns towards the nose,
548 // and positive phi values mean it turns to the left side. (Note: only airspeed
549 // is transformed, the rotational speed contribution is ignored.)
550 
551 void FGRotor::calc_downwash_angles()
552 {
553  FGColumnVector3 v_shaft;
554  v_shaft = TboToHsr * InvTransform * in.AeroUVW;
555 
556  theta_downwash = atan2( -v_shaft(eU), v_induced - v_shaft(eW)) + a1s;
557  phi_downwash = atan2( v_shaft(eV), v_induced - v_shaft(eW)) + b1s;
558 
559  return;
560 }
561 
562 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
563 
564 // transform rotor forces from control axes to shaft axes, and express
565 // in body axes /SH79/ eqn(40,41)
566 
567 FGColumnVector3 FGRotor::body_forces(double a_ic, double b_ic)
568 {
569  FGColumnVector3 F_s(
570  - H_drag*cos(beta_orient) - J_side*sin(beta_orient) + Thrust*b_ic,
571  - H_drag*sin(beta_orient) + J_side*cos(beta_orient) + Thrust*a_ic,
572  - Thrust);
573 
574  return HsrToTbo * F_s;
575 }
576 
577 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
578 
579 // calculates the additional moments due to hinge offset and handles
580 // torque and sense
581 
582 FGColumnVector3 FGRotor::body_moments(double a_ic, double b_ic)
583 {
584  FGColumnVector3 M_s, M_hub, M_h;
585  double mf;
586 
587  // cyclic flapping relative to shaft axes /SH79/ eqn(43)
588  a1s = a_1*cos(beta_orient) + b_1*sin(beta_orient) - b_ic;
589  b1s = b_1*cos(beta_orient) - a_1*sin(beta_orient) + a_ic;
590 
591  mf = 0.5 * HingeOffset * BladeNum * Omega*Omega * BladeMassMoment;
592 
593  M_s(eL) = mf*b1s;
594  M_s(eM) = mf*a1s;
595  M_s(eN) = Torque * Sense ;
596 
597  return HsrToTbo * M_s;
598 }
599 
600 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
601 
602 void FGRotor::CalcRotorState(void)
603 {
604  double A_IC; // lateral (roll) control in radians
605  double B_IC; // longitudinal (pitch) control in radians
606  double theta_col; // rotor collective pitch in radians
607 
608  FGColumnVector3 vHub_ca, avFus_ca;
609 
610  double filtered_hagl = 0.0;
611  double ge_factor = 1.0;
612 
613  // fetch needed values from environment
614  rho = in.Density; // slugs/ft^3.
615  double h_agl_ft = in.H_agl;
616 
617  // update InvTransform, the rotor orientation could have been altered
618  InvTransform = Transform().Transposed();
619 
620  // handle RPM requirements, calc omega.
621  if (ExternalRPM && ExtRPMsource) {
622  RPM = ExtRPMsource->getDoubleValue() * ( SourceGearRatio / GearRatio );
623  }
624 
625  // MinimalRPM is always >= 1. MaximalRPM is always >= NominalRPM
626  RPM = Constrain(MinimalRPM, RPM, MaximalRPM);
627 
628  Omega = (RPM/60.0)*2.0*M_PI;
629 
630  // set control inputs
631  A_IC = LateralCtrl;
632  B_IC = LongitudinalCtrl;
633  theta_col = CollectiveCtrl;
634 
635  // optional ground effect, a ge_factor of 1.0 gives no effect
636  // and 0.5 yields to maximal influence.
637  if (GroundEffectExp > 1e-5) {
638  if (h_agl_ft<0.0) h_agl_ft = 0.0; // clamp
639  filtered_hagl = damp_hagl.execute(h_agl_ft) + GroundEffectShift;
640  // actual/nominal factor avoids absurd scales at startup
641  ge_factor -= GroundEffectScaleNorm *
642  ( exp(-filtered_hagl*GroundEffectExp) * (RPM / NominalRPM) );
643  ge_factor = Constrain(0.5, ge_factor, 1.0);
644  }
645 
646  // all set, start calculations ...
647 
648  vHub_ca = hub_vel_body2ca(in.AeroUVW, in.AeroPQR, A_IC, B_IC);
649 
650  avFus_ca = fus_angvel_body2ca(in.AeroPQR);
651 
652  calc_flow_and_thrust(theta_col, vHub_ca(eU), vHub_ca(eW), ge_factor);
653 
654  calc_coning_angle(theta_col);
655 
656  calc_flapping_angles(theta_col, avFus_ca);
657 
658  calc_drag_and_side_forces(theta_col);
659 
660  calc_torque(theta_col);
661 
662  calc_downwash_angles();
663 
664  // ... and assign to inherited vFn and vMn members
665  // (for processing see FGForce::GetBodyForces).
666  vFn = body_forces(A_IC, B_IC);
667  vMn = Transform() * body_moments(A_IC, B_IC);
668 
669 }
670 
671 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
672 
673 double FGRotor::Calculate(double EnginePower)
674 {
675 
676  CalcRotorState();
677 
678  if (! ExternalRPM) {
679  // the RPM values are handled inside Transmission
680  Transmission->Calculate(EnginePower, Torque, in.TotalDeltaT);
681 
682  EngineRPM = Transmission->GetEngineRPM() * GearRatio;
683  RPM = Transmission->GetThrusterRPM();
684  } else {
685  EngineRPM = RPM * GearRatio;
686  }
687 
688  RPM = Constrain(MinimalRPM, RPM, MaximalRPM); // trim again
689 
690  return Thrust;
691 }
692 
693 
694 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
695 
696 
697 bool FGRotor::bindmodel(FGPropertyManager* PropertyManager)
698 {
699  string property_name, base_property_name;
700  base_property_name = CreateIndexedPropertyName("propulsion/engine", EngineNum);
701 
702  property_name = base_property_name + "/rotor-rpm";
703  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetRPM );
704 
705  property_name = base_property_name + "/engine-rpm";
706  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetEngineRPM );
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  property_name = base_property_name + "/groundeffect-scale-norm";
742  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetGroundEffectScaleNorm,
744 
745  switch (ControlMap) {
746  case eTailCtrl:
747  property_name = base_property_name + "/antitorque-ctrl-rad";
748  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
749  break;
750  case eTandemCtrl:
751  property_name = base_property_name + "/tail-collective-ctrl-rad";
752  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
753  property_name = base_property_name + "/lateral-ctrl-rad";
754  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
755  property_name = base_property_name + "/longitudinal-ctrl-rad";
756  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
757  break;
758  default: // eMainCtrl
759  property_name = base_property_name + "/collective-ctrl-rad";
760  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
761  property_name = base_property_name + "/lateral-ctrl-rad";
762  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
763  property_name = base_property_name + "/longitudinal-ctrl-rad";
764  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
765  }
766 
767  if (ExternalRPM) {
768  if (RPMdefinition == -1) {
769  property_name = base_property_name + "/x-rpm-dict";
770  ExtRPMsource = PropertyManager->GetNode(property_name, true);
771  } else if (RPMdefinition >= 0 && RPMdefinition != EngineNum) {
772  string ipn = CreateIndexedPropertyName("propulsion/engine", RPMdefinition);
773  property_name = ipn + "/rotor-rpm";
774  ExtRPMsource = PropertyManager->GetNode(property_name, false);
775  if (! ExtRPMsource) {
776  cerr << "# Warning: Engine number " << EngineNum << "." << endl;
777  cerr << "# No 'rotor-rpm' property found for engine " << RPMdefinition << "." << endl;
778  cerr << "# Please check order of engine definitons." << endl;
779  }
780  } else {
781  cerr << "# Engine number " << EngineNum;
782  cerr << ", given ExternalRPM value '" << RPMdefinition << "' unhandled." << endl;
783  }
784  }
785 
786  return true;
787 }
788 
789 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
790 
791 string FGRotor::GetThrusterLabels(int id, const string& delimeter)
792 {
793 
794  ostringstream buf;
795 
796  buf << Name << " RPM (engine " << id << ")";
797 
798  return buf.str();
799 
800 }
801 
802 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
803 
804 string FGRotor::GetThrusterValues(int id, const string& delimeter)
805 {
806 
807  ostringstream buf;
808 
809  buf << RPM;
810 
811  return buf.str();
812 
813 }
814 
815 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
816 // The bitmasked value choices are as follows:
817 // unset: In this case (the default) JSBSim would only print
818 // out the normally expected messages, essentially echoing
819 // the config files as they are read. If the environment
820 // variable is not set, debug_lvl is set to 1 internally
821 // 0: This requests JSBSim not to output any messages
822 // whatsoever.
823 // 1: This value explicity requests the normal JSBSim
824 // startup messages
825 // 2: This value asks for a message to be printed out when
826 // a class is instantiated
827 // 4: When this value is set, a message is displayed when a
828 // FGModel object executes its Run() method
829 // 8: When this value is set, various runtime state variables
830 // are printed out periodically
831 // 16: When set various parameters are sanity checked and
832 // a message is printed out when they go out of bounds
833 
834 void FGRotor::Debug(int from)
835 {
836  string ControlMapName;
837 
838  if (debug_lvl <= 0) return;
839 
840  if (debug_lvl & 1) { // Standard console startup message output
841  if (from == 0) { // Constructor
842  cout << "\n Rotor Name: " << Name << endl;
843  cout << " Diameter = " << 2.0 * Radius << " ft." << endl;
844  cout << " Number of Blades = " << BladeNum << endl;
845  cout << " Gear Ratio = " << GearRatio << endl;
846  cout << " Sense = " << Sense << endl;
847  cout << " Nominal RPM = " << NominalRPM << endl;
848  cout << " Minimal RPM = " << MinimalRPM << endl;
849  cout << " Maximal RPM = " << MaximalRPM << endl;
850 
851  if (ExternalRPM) {
852  if (RPMdefinition == -1) {
853  cout << " RPM is controlled externally" << endl;
854  } else {
855  cout << " RPM source set to thruster " << RPMdefinition << endl;
856  }
857  }
858 
859  cout << " Blade Chord = " << BladeChord << endl;
860  cout << " Lift Curve Slope = " << LiftCurveSlope << endl;
861  cout << " Blade Twist = " << BladeTwist << endl;
862  cout << " Hinge Offset = " << HingeOffset << endl;
863  cout << " Blade Flapping Moment = " << BladeFlappingMoment << endl;
864  cout << " Blade Mass Moment = " << BladeMassMoment << endl;
865  cout << " Polar Moment = " << PolarMoment << endl;
866  cout << " Inflow Lag = " << InflowLag << endl;
867  cout << " Tip Loss = " << TipLossB << endl;
868  cout << " Lock Number = " << LockNumberByRho * 0.002356 << " (SL)" << endl;
869  cout << " Solidity = " << Solidity << endl;
870  cout << " Max Brake Power = " << MaxBrakePower/hptoftlbssec << " HP" << endl;
871  cout << " Gear Loss = " << GearLoss/hptoftlbssec << " HP" << endl;
872  cout << " Gear Moment = " << GearMoment << endl;
873 
874  switch (ControlMap) {
875  case eTailCtrl: ControlMapName = "Tail Rotor"; break;
876  case eTandemCtrl: ControlMapName = "Tandem Rotor"; break;
877  default: ControlMapName = "Main Rotor";
878  }
879  cout << " Control Mapping = " << ControlMapName << endl;
880 
881  }
882  }
883  if (debug_lvl & 2 ) { // Instantiation/Destruction notification
884  if (from == 0) cout << "Instantiated: FGRotor" << endl;
885  if (from == 1) cout << "Destroyed: FGRotor" << endl;
886  }
887  if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
888  }
889  if (debug_lvl & 8 ) { // Runtime state variables
890  }
891  if (debug_lvl & 16) { // Sanity checking
892  }
893  if (debug_lvl & 64) {
894  if (from == 0) { // Constructor
895  cout << IdSrc << endl;
896  cout << IdHdr << endl;
897  }
898  }
899 
900 }
901 
902 
903 } // namespace JSBSim
904 
Element * GetParent(void)
Returns a pointer to the parent of an element.
Definition: FGXMLElement.h:231
double GetGroundEffectScaleNorm(void) const
Retrieves the ground effect scaling factor.
Definition: FGRotor.h:302
double Calculate(double EnginePower)
Returns the scalar thrust of the rotor, and adjusts the RPM value.
Definition: FGRotor.cpp:673
double GetThetaDW(void) const
Downwash angle - positive values point forward (given a horizontal spinning rotor) ...
Definition: FGRotor.h:297
~FGRotor()
Destructor for FGRotor.
Definition: FGRotor.cpp:223
double GetVi(void) const
Retrieves the induced velocity.
Definition: FGRotor.h:290
static double Constrain(double min, double value, double max)
Constrain a value between a minimum and a maximum value.
Definition: FGJSBBase.h:332
void SetLateralCtrl(double c)
Sets the lateral control input in radians.
Definition: FGRotor.h:316
double FindElementValueAsNumberConvertTo(const std::string &el, const std::string &target_units)
Searches for the named element and converts and returns the data belonging to it. ...
void SetLongitudinalCtrl(double c)
Sets the longitudinal control input in radians.
Definition: FGRotor.h:318
FGEngine * GetEngine(unsigned int index) const
Retrieves an engine object pointer from the list of engines.
Definition: FGPropulsion.h:140
double GetCT(void) const
Retrieves the thrust coefficient.
Definition: FGRotor.h:292
Element * FindElement(const std::string &el="")
Searches for a specified element.
double FindElementValueAsNumber(const std::string &el="")
Searches for the named element and returns the data belonging to it as a number.
Utility class that handles power transmission in conjunction with FGRotor.
Base class for specific thrusting devices such as propellers, nozzles, etc.
Definition: FGThruster.h:84
double GetRPM(void) const
Retrieves the RPMs of the rotor.
Definition: FGRotor.h:265
double GetA1(void) const
Retrieves the longitudinal flapping angle with respect to the rotor shaft.
Definition: FGRotor.h:279
FGPropertyManager * GetPropertyManager(void)
Returns a pointer to the property manager object.
Definition: FGFDMExec.cpp:1099
double GetEngineRPM(void) const
Retrieves the RPMs of the Engine, as seen from this rotor.
Definition: FGRotor.h:269
void SetCollectiveCtrl(double c)
Sets the collective control input in radians.
Definition: FGRotor.h:314
First order, (low pass / lag) filter.
Definition: FGJSBBase.h:102
double GetLongitudinalCtrl(void) const
Retrieves the longitudinal control input in radians.
Definition: FGRotor.h:311
void Tie(const std::string &name, bool *pointer, bool useDefault=true)
Tie a property to an external bool variable.
double GetPhiDW(void) const
Downwash angle - positive values point leftward (given a horizontal spinning rotor) ...
Definition: FGRotor.h:299
void SetGroundEffectScaleNorm(double g)
Sets the ground effect scaling factor.
Definition: FGRotor.h:304
double GetDataAsNumber(void)
Converts the element data to a number.
const std::string & GetName(void) const
Retrieves the element name.
Definition: FGXMLElement.h:186
std::string FindElementValue(const std::string &el="")
Searches for the named element and returns the string data belonging to it.
double GetCollectiveCtrl(void) const
Retrieves the collective control input in radians.
Definition: FGRotor.h:307
double GetNu(void) const
Retrieves the induced inflow ratio.
Definition: FGRotor.h:288
double GetA0(void) const
Retrieves the rotor&#39;s coning angle.
Definition: FGRotor.h:277
This class implements a 3 element column vector.
double GetLambda(void) const
Retrieves the inflow ratio.
Definition: FGRotor.h:284
double GetLateralCtrl(void) const
Retrieves the lateral control input in radians.
Definition: FGRotor.h:309
FGRotor(FGFDMExec *exec, Element *rotor_element, int num)
Constructor for FGRotor.
Definition: FGRotor.cpp:77
double GetB1(void) const
Retrieves the lateral flapping angle with respect to the rotor shaft.
Definition: FGRotor.h:281
double GetTorque(void) const
Retrieves the torque.
Definition: FGRotor.h:294
FGMatrix33 Transposed(void) const
Transposed matrix.
Definition: FGMatrix33.h:243
void InitMatrix(void)
Initialize the matrix.
Definition: FGMatrix33.cpp:257
double GetDeltaT(void) const
Returns the simulation delta T.
Definition: FGFDMExec.h:536
FGPropulsion * GetPropulsion(void)
Returns the FGPropulsion pointer.
Definition: FGFDMExec.h:353
Encapsulates the JSBSim simulation executive.
Definition: FGFDMExec.h:189
double GetMu(void) const
Retrieves the tip-speed (aka advance) ratio.
Definition: FGRotor.h:286
FGColumnVector3 StructuralToBody(const FGColumnVector3 &r) const
Conversion from the structural frame to the body frame.
FGColumnVector3 FindElementTripletConvertTo(const std::string &target_units)
Composes a 3-element column vector for the supplied location or orientation.
FGMassBalance * GetMassBalance(void)
Returns the FGAircraft pointer.
Definition: FGFDMExec.h:355