1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 Purpose: Encapsulates the rotor object
8 ------------- Copyright (C) 2000 Jon S. Berndt (jon@jsbsim.org) -------------
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
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
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.
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.
27 FUNCTIONAL DESCRIPTION
28 --------------------------------------------------------------------------------
31 --------------------------------------------------------------------------------
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
42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
49 #include "models/FGMassBalance.h"
50 #include "models/FGPropulsion.h" // to get the GearRatio from a linked rotor
55 using std::ostringstream;
59 static const char *IdSrc = "$Id: FGRotor.cpp,v 1.20 2012/03/18 15:48:36 jentron Exp $";
60 static const char *IdHdr = ID_ROTOR;
62 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
66 static inline double sqr(double x) { return x*x; }
68 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
70 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
74 FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
75 : FGThruster(exec, rotor_element, num),
76 rho(0.002356), // environment
77 Radius(0.0), BladeNum(0), // configuration parameters
78 Sense(1.0), NominalRPM(0.0), MinimalRPM(0.0), MaximalRPM(0.0),
79 ExternalRPM(0), RPMdefinition(0), ExtRPMsource(NULL), SourceGearRatio(1.0),
80 BladeChord(0.0), LiftCurveSlope(0.0), BladeTwist(0.0), HingeOffset(0.0),
81 BladeFlappingMoment(0.0), BladeMassMoment(0.0), PolarMoment(0.0),
82 InflowLag(0.0), TipLossB(0.0),
83 GroundEffectExp(0.0), GroundEffectShift(0.0), GroundEffectScaleNorm(1.0),
84 LockNumberByRho(0.0), Solidity(0.0), // derived parameters
85 RPM(0.0), Omega(0.0), // dynamic values
87 a0(0.0), a_1(0.0), b_1(0.0), a_dw(0.0),
89 H_drag(0.0), J_side(0.0), Torque(0.0), C_T(0.0),
90 lambda(-0.001), mu(0.0), nu(0.001), v_induced(0.0),
91 theta_downwash(0.0), phi_downwash(0.0),
92 ControlMap(eMainCtrl), // control
93 CollectiveCtrl(0.0), LateralCtrl(0.0), LongitudinalCtrl(0.0),
94 Transmission(NULL), // interaction with engine
95 EngineRPM(0.0), MaxBrakePower(0.0), GearLoss(0.0), GearMoment(0.0)
97 FGColumnVector3 location(0.0, 0.0, 0.0), orientation(0.0, 0.0, 0.0);
98 Element *thruster_element;
99 double engine_power_est = 0.0;
101 // initialise/set remaining variables
102 SetTransformType(FGForce::tCustom);
103 PropertyManager = exec->GetPropertyManager();
107 dt = exec->GetDeltaT();
108 for (int i=0; i<5; i++) R[i] = 0.0;
109 for (int i=0; i<5; i++) B[i] = 0.0;
112 thruster_element = rotor_element->GetParent()->FindElement("sense");
113 if (thruster_element) {
114 double s = thruster_element->GetDataAsNumber();
116 Sense = -1.0; // 'CW' as seen from above
117 } else if (s < 0.1) {
118 Sense = 0.0; // 'coaxial'
120 Sense = 1.0; // 'CCW' as seen from above
124 thruster_element = rotor_element->GetParent()->FindElement("location");
125 if (thruster_element) {
126 location = thruster_element->FindElementTripletConvertTo("IN");
128 cerr << "No thruster location found." << endl;
131 thruster_element = rotor_element->GetParent()->FindElement("orient");
132 if (thruster_element) {
133 orientation = thruster_element->FindElementTripletConvertTo("RAD");
135 cerr << "No thruster orientation found." << endl;
138 SetLocation(location);
139 SetAnglesToBody(orientation);
140 InvTransform = Transform().Transposed(); // body to custom/native
143 ControlMap = eMainCtrl;
144 if (rotor_element->FindElement("controlmap")) {
145 string cm = rotor_element->FindElementValue("controlmap");
148 ControlMap = eTailCtrl;
149 } else if (cm == "TANDEM") {
150 ControlMap = eTandemCtrl;
152 cerr << "# found unknown controlmap: '" << cm << "' using main rotor config." << endl;
156 // ExternalRPM -- is the RPM dictated ?
157 if (rotor_element->FindElement("ExternalRPM")) {
159 SourceGearRatio = 1.0;
160 RPMdefinition = (int) rotor_element->FindElementValueAsNumber("ExternalRPM");
161 int rdef = RPMdefinition;
162 if (RPMdefinition>=0) {
163 // avoid ourself and (still) unknown engines.
164 if (!exec->GetPropulsion()->GetEngine(RPMdefinition) || RPMdefinition==num) {
167 FGThruster *tr = exec->GetPropulsion()->GetEngine(RPMdefinition)->GetThruster();
168 SourceGearRatio = tr->GetGearRatio();
169 // cout << "# got sources' GearRatio: " << SourceGearRatio << endl;
172 if (RPMdefinition != rdef) {
173 cerr << "# discarded given RPM source (" << rdef << ") and switched to external control (-1)." << endl;
177 // process rotor parameters
178 engine_power_est = Configure(rotor_element);
180 // setup transmission if needed
183 Transmission = new FGTransmission(exec, num, dt);
185 Transmission->SetThrusterMoment(PolarMoment);
187 // The MOI sensed behind the gear ( MOI_engine*sqr(GearRatio) ).
188 GearMoment = ConfigValueConv(rotor_element, "gearmoment", 0.1*PolarMoment, "SLUG*FT2");
189 GearMoment = Constrain(1e-6, GearMoment, 1e9);
190 Transmission->SetEngineMoment(GearMoment);
192 Transmission->SetMaxBrakePower(MaxBrakePower);
194 GearLoss = ConfigValueConv(rotor_element, "gearloss", 0.0025 * engine_power_est, "HP");
195 GearLoss = Constrain(0.0, GearLoss, 1e9);
196 GearLoss *= hptoftlbssec;
197 Transmission->SetEngineFriction(GearLoss);
201 // shaft representation - a rather simple transform,
202 // but using a matrix is safer.
203 TboToHsr.InitMatrix( 0.0, 0.0, 1.0,
206 HsrToTbo = TboToHsr.Transposed();
208 // smooth out jumps in hagl reported, otherwise the ground effect
209 // calculation would cause jumps too. 1Hz seems sufficient.
210 damp_hagl = Filter(1.0, dt);
212 // enable import-export
219 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
222 if (Transmission) delete Transmission;
226 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
228 // 5in1: value-fetch-convert-default-return function
230 double FGRotor::ConfigValueConv( Element* el, const string& ename, double default_val,
231 const string& unit, bool tell)
235 double val = default_val;
237 string pname = "*No parent element*";
240 e = el->FindElement(ename);
241 pname = el->GetName();
246 val = e->GetDataAsNumber();
248 val = el->FindElementValueAsNumberConvertTo(ename,unit);
252 cerr << pname << ": missing element '" << ename <<
253 "' using estimated value: " << default_val << endl;
260 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
262 double FGRotor::ConfigValue(Element* el, const string& ename, double default_val, bool tell)
264 return ConfigValueConv(el, ename, default_val, "", tell);
267 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
269 // 1. read configuration and try to fill holes, ymmv
270 // 2. calculate derived parameters
271 double FGRotor::Configure(Element* rotor_element)
274 double estimate, engine_power_est=0.0;
275 const bool yell = true;
276 const bool silent = false;
279 Radius = 0.5 * ConfigValueConv(rotor_element, "diameter", 42.0, "FT", yell);
280 Radius = Constrain(1e-3, Radius, 1e9);
282 BladeNum = (int) ConfigValue(rotor_element, "numblades", 3 , yell);
284 GearRatio = ConfigValue(rotor_element, "gearratio", 1.0, yell);
285 GearRatio = Constrain(1e-9, GearRatio, 1e9);
287 // make sure that v_tip (omega*r) is below 0.7mach ~ 750ft/s
288 estimate = (750.0/Radius)/(2.0*M_PI) * 60.0; // 7160/Radius
289 NominalRPM = ConfigValue(rotor_element, "nominalrpm", estimate, yell);
290 NominalRPM = Constrain(2.0, NominalRPM, 1e9);
292 MinimalRPM = ConfigValue(rotor_element, "minrpm", 1.0);
293 MinimalRPM = Constrain(1.0, MinimalRPM, NominalRPM - 1.0);
295 MaximalRPM = ConfigValue(rotor_element, "maxrpm", 2.0*NominalRPM);
296 MaximalRPM = Constrain(NominalRPM, MaximalRPM, 1e9);
298 estimate = Constrain(0.07, 2.0/Radius , 0.14); // guess solidity
299 estimate = estimate * M_PI*Radius/BladeNum;
300 BladeChord = ConfigValueConv(rotor_element, "chord", estimate, "FT", yell);
302 LiftCurveSlope = ConfigValue(rotor_element, "liftcurveslope", 6.0); // "1/RAD"
303 BladeTwist = ConfigValueConv(rotor_element, "twist", -0.17, "RAD");
305 HingeOffset = ConfigValueConv(rotor_element, "hingeoffset", 0.05 * Radius, "FT" );
307 estimate = sqr(BladeChord) * sqr(Radius - HingeOffset) * 0.57;
308 BladeFlappingMoment = ConfigValueConv(rotor_element, "flappingmoment", estimate, "SLUG*FT2");
309 BladeFlappingMoment = Constrain(1e-9, BladeFlappingMoment, 1e9);
311 // guess mass from moment of a thin stick, and multiply by the blades cg distance
312 estimate = ( 3.0 * BladeFlappingMoment / sqr(Radius) ) * (0.45 * Radius) ;
313 BladeMassMoment = ConfigValue(rotor_element, "massmoment", estimate); // unit is slug-ft
314 BladeMassMoment = Constrain(1e-9, BladeMassMoment, 1e9);
316 estimate = 1.1 * BladeFlappingMoment * BladeNum;
317 PolarMoment = ConfigValueConv(rotor_element, "polarmoment", estimate, "SLUG*FT2");
318 PolarMoment = Constrain(1e-9, PolarMoment, 1e9);
320 // "inflowlag" is treated further down.
322 TipLossB = ConfigValue(rotor_element, "tiplossfactor", 1.0, silent);
324 // estimate engine power (bit of a pity, cause our caller already knows)
325 engine_power_est = 0.5 * BladeNum*BladeChord*Radius*Radius;
327 estimate = engine_power_est / 30.0;
328 MaxBrakePower = ConfigValueConv(rotor_element, "maxbrakepower", estimate, "HP");
329 MaxBrakePower *= hptoftlbssec;
331 GroundEffectExp = ConfigValue(rotor_element, "groundeffectexp", 0.0);
332 GroundEffectShift = ConfigValueConv(rotor_element, "groundeffectshift", 0.0, "FT");
334 // precalc often used powers
335 R[0]=1.0; R[1]=Radius; R[2]=R[1]*R[1]; R[3]=R[2]*R[1]; R[4]=R[3]*R[1];
336 B[0]=1.0; B[1]=TipLossB; B[2]=B[1]*B[1]; B[3]=B[2]*B[1]; B[4]=B[3]*B[1];
338 // derived parameters
339 LockNumberByRho = LiftCurveSlope * BladeChord * R[4] / BladeFlappingMoment;
340 Solidity = BladeNum * BladeChord / (M_PI * Radius);
342 // estimate inflow lag, see /GE49/ eqn(1)
343 double omega_tmp = (NominalRPM/60.0)*2.0*M_PI;
344 estimate = 16.0/(LockNumberByRho*rho * omega_tmp ); // 16/(gamma*Omega)
345 // printf("# Est. InflowLag: %f\n", estimate);
346 InflowLag = ConfigValue(rotor_element, "inflowlag", estimate, yell);
347 InflowLag = Constrain(1e-6, InflowLag, 2.0);
349 return engine_power_est;
352 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
354 // calculate control-axes components of total airspeed at the hub.
355 // sets rotor orientation angle (beta) as side effect. /SH79/ eqn(19-22)
357 FGColumnVector3 FGRotor::hub_vel_body2ca( const FGColumnVector3 &uvw,
358 const FGColumnVector3 &pqr,
359 double a_ic, double b_ic)
361 FGColumnVector3 v_r, v_shaft, v_w;
364 pos = fdmex->GetMassBalance()->StructuralToBody(GetActingLocation());
367 v_shaft = TboToHsr * InvTransform * v_r;
369 beta_orient = atan2(v_shaft(eV),v_shaft(eU));
371 v_w(eU) = v_shaft(eU)*cos(beta_orient) + v_shaft(eV)*sin(beta_orient);
373 v_w(eW) = v_shaft(eW) - b_ic*v_shaft(eU) - a_ic*v_shaft(eV);
378 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
380 // express fuselage angular velocity in control axes /SH79/ eqn(30,31)
382 FGColumnVector3 FGRotor::fus_angvel_body2ca( const FGColumnVector3 &pqr)
384 FGColumnVector3 av_s_fus, av_w_fus;
387 // av_s_fus = BodyToShaft * pqr; /SH79/
388 // BodyToShaft = TboToHsr * InvTransform
389 av_s_fus = TboToHsr * InvTransform * pqr;
391 av_w_fus(eP)= av_s_fus(eP)*cos(beta_orient) + av_s_fus(eQ)*sin(beta_orient);
392 av_w_fus(eQ)= - av_s_fus(eP)*sin(beta_orient) + av_s_fus(eQ)*cos(beta_orient);
393 av_w_fus(eR)= av_s_fus(eR);
399 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
401 // The calculation is a bit tricky because thrust depends on induced velocity,
404 // The flow_scale parameter (ranging from 0.5-1.0) is used to approximate a
405 // reduction of inflow if the helicopter is close to the ground, yielding to
406 // higher thrust, see /TA77/ eqn(10a).
408 void FGRotor::calc_flow_and_thrust( double theta_0, double Uw, double Ww,
412 double ct_over_sigma = 0.0;
413 double c0, ct_l, ct_t0, ct_t1;
416 mu = Uw/(Omega*Radius); // /SH79/ eqn(24)
417 if (mu > 0.7) mu = 0.7;
420 ct_t0 = (1.0/3.0*B[3] + 1.0/2.0 * TipLossB*mu2 - 4.0/(9.0*M_PI) * mu*mu2 ) * theta_0;
421 ct_t1 = (1.0/4.0*B[4] + 1.0/4.0 * B[2]*mu2) * BladeTwist;
423 ct_l = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda; // first time
425 c0 = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1) * Solidity;
426 c0 = c0 / ( 2.0 * sqrt( sqr(mu) + sqr(lambda) ) + 1e-15);
428 // replacement for /SH79/ eqn(26).
429 // ref: dnu/dt = 1/tau ( Ct / (2*sqrt(mu^2+lambda^2)) - nu )
430 // taking mu and lambda constant, this integrates to
432 nu = flow_scale * ((nu - c0) * exp(-dt/InflowLag) + c0);
434 // now from nu to lambda, C_T, and Thrust
436 lambda = Ww/(Omega*Radius) - nu; // /SH79/ eqn(25)
438 ct_l = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda;
440 ct_over_sigma = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1); // /SH79/ eqn(27)
442 Thrust = BladeNum*BladeChord*Radius*rho*sqr(Omega*Radius) * ct_over_sigma;
444 C_T = ct_over_sigma * Solidity;
445 v_induced = nu * (Omega*Radius);
450 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
452 // Two blade teetering rotors are often 'preconed' to a fixed angle, but the
453 // calculated value is pretty close to the real one. /SH79/ eqn(29)
455 void FGRotor::calc_coning_angle(double theta_0)
457 double lock_gamma = LockNumberByRho * rho;
459 double a0_l = (1.0/6.0 + 0.04 * mu*mu*mu) * lambda;
460 double a0_t0 = (1.0/8.0 + 1.0/8.0 * mu*mu) * theta_0;
461 double a0_t1 = (1.0/10.0 + 1.0/12.0 * mu*mu) * BladeTwist;
462 a0 = lock_gamma * ( a0_l + a0_t0 + a0_t1);
466 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
468 // Flapping angles relative to control axes /SH79/ eqn(32)
470 void FGRotor::calc_flapping_angles(double theta_0, const FGColumnVector3 &pqr_fus_w)
472 double lock_gamma = LockNumberByRho * rho;
475 double mu2_2 = sqr(mu)/2.0;
476 double t075 = theta_0 + 0.75 * BladeTwist; // common approximation for rectangular blades
478 a_1 = 1.0/(1.0 - mu2_2) * (
479 (2.0*lambda + (8.0/3.0)*t075)*mu
480 + pqr_fus_w(eP)/Omega
481 - 16.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
484 b_1 = 1.0/(1.0 + mu2_2) * (
486 - pqr_fus_w(eQ)/Omega
487 - 16.0 * pqr_fus_w(eP)/(lock_gamma*Omega)
490 // used in force calc
491 a_dw = 1.0/(1.0 - mu2_2) * (
492 (2.0*lambda + (8.0/3.0)*t075)*mu
493 - 24.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
494 * ( 1.0 - ( 0.29 * t075 / (C_T/Solidity) ) )
500 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
504 void FGRotor::calc_drag_and_side_forces(double theta_0)
506 double cy_over_sigma;
507 double t075 = theta_0 + 0.75 * BladeTwist;
509 H_drag = Thrust * a_dw;
512 0.75*b_1*lambda - 1.5*a0*mu*lambda + 0.25*a_1*b_1*mu
513 - a0*a_1*sqr(mu) + (1.0/6.0)*a0*a_1
514 - (0.75*mu*a0 - (1.0/3.0)*b_1 - 0.5*sqr(mu)*b_1)*t075
516 cy_over_sigma *= LiftCurveSlope/2.0;
518 J_side = BladeNum * BladeChord * Radius * rho * sqr(Omega*Radius) * cy_over_sigma;
523 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
525 // Simplified version of /SH79/ eqn(36). Uses an estimate for blade drag
526 // (a new config parameter to come...).
527 // From "Bramwell's Helicopter Dynamics", second edition, eqn(3.43) and (3.44)
529 void FGRotor::calc_torque(double theta_0)
531 // estimate blade drag
532 double delta_dr = 0.009 + 0.3*sqr(6.0*C_T/(LiftCurveSlope*Solidity));
534 Torque = rho * BladeNum * BladeChord * delta_dr * sqr(Omega*Radius) * R[2] *
535 (1.0+4.5*sqr(mu))/8.0
536 - (Thrust*lambda + H_drag*mu)*Radius;
541 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
543 // Get the downwash angles with respect to the shaft axis.
544 // Given a 'regular' main rotor, the angles are zero when the downwash points
545 // down, positive theta values mean that the downwash turns towards the nose,
546 // and positive phi values mean it turns to the left side. (Note: only airspeed
547 // is transformed, the rotational speed contribution is ignored.)
549 void FGRotor::calc_downwash_angles()
551 FGColumnVector3 v_shaft;
552 v_shaft = TboToHsr * InvTransform * in.AeroUVW;
554 theta_downwash = atan2( -v_shaft(eU), v_induced - v_shaft(eW)) + a1s;
555 phi_downwash = atan2( v_shaft(eV), v_induced - v_shaft(eW)) + b1s;
560 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
562 // transform rotor forces from control axes to shaft axes, and express
563 // in body axes /SH79/ eqn(40,41)
565 FGColumnVector3 FGRotor::body_forces(double a_ic, double b_ic)
568 - H_drag*cos(beta_orient) - J_side*sin(beta_orient) + Thrust*b_ic,
569 - H_drag*sin(beta_orient) + J_side*cos(beta_orient) + Thrust*a_ic,
572 return HsrToTbo * F_s;
575 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
577 // calculates the additional moments due to hinge offset and handles
580 FGColumnVector3 FGRotor::body_moments(double a_ic, double b_ic)
582 FGColumnVector3 M_s, M_hub, M_h;
585 // cyclic flapping relative to shaft axes /SH79/ eqn(43)
586 a1s = a_1*cos(beta_orient) + b_1*sin(beta_orient) - b_ic;
587 b1s = b_1*cos(beta_orient) - a_1*sin(beta_orient) + a_ic;
589 mf = 0.5 * HingeOffset * BladeNum * Omega*Omega * BladeMassMoment;
593 M_s(eN) = Torque * Sense ;
595 return HsrToTbo * M_s;
598 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
600 void FGRotor::CalcRotorState(void)
602 double A_IC; // lateral (roll) control in radians
603 double B_IC; // longitudinal (pitch) control in radians
604 double theta_col; // rotor collective pitch in radians
606 FGColumnVector3 vHub_ca, avFus_ca;
608 double filtered_hagl = 0.0;
609 double ge_factor = 1.0;
611 // fetch needed values from environment
612 rho = in.Density; // slugs/ft^3.
613 double h_agl_ft = in.H_agl;
615 // update InvTransform, the rotor orientation could have been altered
616 InvTransform = Transform().Transposed();
618 // handle RPM requirements, calc omega.
619 if (ExternalRPM && ExtRPMsource) {
620 RPM = ExtRPMsource->getDoubleValue() * ( SourceGearRatio / GearRatio );
623 // MinimalRPM is always >= 1. MaximalRPM is always >= NominalRPM
624 RPM = Constrain(MinimalRPM, RPM, MaximalRPM);
626 Omega = (RPM/60.0)*2.0*M_PI;
628 // set control inputs
630 B_IC = LongitudinalCtrl;
631 theta_col = CollectiveCtrl;
633 // optional ground effect, a ge_factor of 1.0 gives no effect
634 // and 0.5 yields to maximal influence.
635 if (GroundEffectExp > 1e-5) {
636 if (h_agl_ft<0.0) h_agl_ft = 0.0; // clamp
637 filtered_hagl = damp_hagl.execute(h_agl_ft) + GroundEffectShift;
638 // actual/nominal factor avoids absurd scales at startup
639 ge_factor -= GroundEffectScaleNorm *
640 ( exp(-filtered_hagl*GroundEffectExp) * (RPM / NominalRPM) );
641 ge_factor = Constrain(0.5, ge_factor, 1.0);
644 // all set, start calculations ...
646 vHub_ca = hub_vel_body2ca(in.AeroUVW, in.AeroPQR, A_IC, B_IC);
648 avFus_ca = fus_angvel_body2ca(in.AeroPQR);
650 calc_flow_and_thrust(theta_col, vHub_ca(eU), vHub_ca(eW), ge_factor);
652 calc_coning_angle(theta_col);
654 calc_flapping_angles(theta_col, avFus_ca);
656 calc_drag_and_side_forces(theta_col);
658 calc_torque(theta_col);
660 calc_downwash_angles();
662 // ... and assign to inherited vFn and vMn members
663 // (for processing see FGForce::GetBodyForces).
664 vFn = body_forces(A_IC, B_IC);
665 vMn = Transform() * body_moments(A_IC, B_IC);
669 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
671 double FGRotor::Calculate(double EnginePower)
677 // the RPM values are handled inside Transmission
678 Transmission->Calculate(EnginePower, Torque, in.TotalDeltaT);
680 EngineRPM = Transmission->GetEngineRPM() * GearRatio;
681 RPM = Transmission->GetThrusterRPM();
683 EngineRPM = RPM * GearRatio;
686 RPM = Constrain(MinimalRPM, RPM, MaximalRPM); // trim again
692 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
695 bool FGRotor::BindModel(void)
697 string property_name, base_property_name;
698 base_property_name = CreateIndexedPropertyName("propulsion/engine", EngineNum);
700 property_name = base_property_name + "/rotor-rpm";
701 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetRPM );
703 property_name = base_property_name + "/engine-rpm";
704 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetEngineRPM );
706 property_name = base_property_name + "/a0-rad";
707 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA0 );
709 property_name = base_property_name + "/a1-rad";
710 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA1 );
712 property_name = base_property_name + "/b1-rad";
713 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetB1 );
715 property_name = base_property_name + "/inflow-ratio";
716 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLambda );
718 property_name = base_property_name + "/advance-ratio";
719 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetMu );
721 property_name = base_property_name + "/induced-inflow-ratio";
722 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetNu );
724 property_name = base_property_name + "/vi-fps";
725 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetVi );
727 property_name = base_property_name + "/thrust-coefficient";
728 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCT );
730 property_name = base_property_name + "/torque-lbsft";
731 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetTorque );
733 property_name = base_property_name + "/theta-downwash-rad";
734 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetThetaDW );
736 property_name = base_property_name + "/phi-downwash-rad";
737 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetPhiDW );
739 property_name = base_property_name + "/groundeffect-scale-norm";
740 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetGroundEffectScaleNorm,
741 &FGRotor::SetGroundEffectScaleNorm );
743 switch (ControlMap) {
745 property_name = base_property_name + "/antitorque-ctrl-rad";
746 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
749 property_name = base_property_name + "/tail-collective-ctrl-rad";
750 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
751 property_name = base_property_name + "/lateral-ctrl-rad";
752 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
753 property_name = base_property_name + "/longitudinal-ctrl-rad";
754 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
756 default: // eMainCtrl
757 property_name = base_property_name + "/collective-ctrl-rad";
758 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
759 property_name = base_property_name + "/lateral-ctrl-rad";
760 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
761 property_name = base_property_name + "/longitudinal-ctrl-rad";
762 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
766 if (RPMdefinition == -1) {
767 property_name = base_property_name + "/x-rpm-dict";
768 ExtRPMsource = PropertyManager->GetNode(property_name, true);
769 } else if (RPMdefinition >= 0 && RPMdefinition != EngineNum) {
770 string ipn = CreateIndexedPropertyName("propulsion/engine", RPMdefinition);
771 property_name = ipn + "/rotor-rpm";
772 ExtRPMsource = PropertyManager->GetNode(property_name, false);
773 if (! ExtRPMsource) {
774 cerr << "# Warning: Engine number " << EngineNum << "." << endl;
775 cerr << "# No 'rotor-rpm' property found for engine " << RPMdefinition << "." << endl;
776 cerr << "# Please check order of engine definitons." << endl;
779 cerr << "# Engine number " << EngineNum;
780 cerr << ", given ExternalRPM value '" << RPMdefinition << "' unhandled." << endl;
787 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
789 string FGRotor::GetThrusterLabels(int id, const string& delimeter)
794 buf << Name << " RPM (engine " << id << ")";
800 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
802 string FGRotor::GetThrusterValues(int id, const string& delimeter)
813 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
814 // The bitmasked value choices are as follows:
815 // unset: In this case (the default) JSBSim would only print
816 // out the normally expected messages, essentially echoing
817 // the config files as they are read. If the environment
818 // variable is not set, debug_lvl is set to 1 internally
819 // 0: This requests JSBSim not to output any messages
821 // 1: This value explicity requests the normal JSBSim
823 // 2: This value asks for a message to be printed out when
824 // a class is instantiated
825 // 4: When this value is set, a message is displayed when a
826 // FGModel object executes its Run() method
827 // 8: When this value is set, various runtime state variables
828 // are printed out periodically
829 // 16: When set various parameters are sanity checked and
830 // a message is printed out when they go out of bounds
832 void FGRotor::Debug(int from)
834 string ControlMapName;
836 if (debug_lvl <= 0) return;
838 if (debug_lvl & 1) { // Standard console startup message output
839 if (from == 0) { // Constructor
840 cout << "\n Rotor Name: " << Name << endl;
841 cout << " Diameter = " << 2.0 * Radius << " ft." << endl;
842 cout << " Number of Blades = " << BladeNum << endl;
843 cout << " Gear Ratio = " << GearRatio << endl;
844 cout << " Sense = " << Sense << endl;
845 cout << " Nominal RPM = " << NominalRPM << endl;
846 cout << " Minimal RPM = " << MinimalRPM << endl;
847 cout << " Maximal RPM = " << MaximalRPM << endl;
850 if (RPMdefinition == -1) {
851 cout << " RPM is controlled externally" << endl;
853 cout << " RPM source set to thruster " << RPMdefinition << endl;
857 cout << " Blade Chord = " << BladeChord << endl;
858 cout << " Lift Curve Slope = " << LiftCurveSlope << endl;
859 cout << " Blade Twist = " << BladeTwist << endl;
860 cout << " Hinge Offset = " << HingeOffset << endl;
861 cout << " Blade Flapping Moment = " << BladeFlappingMoment << endl;
862 cout << " Blade Mass Moment = " << BladeMassMoment << endl;
863 cout << " Polar Moment = " << PolarMoment << endl;
864 cout << " Inflow Lag = " << InflowLag << endl;
865 cout << " Tip Loss = " << TipLossB << endl;
866 cout << " Lock Number = " << LockNumberByRho * 0.002356 << " (SL)" << endl;
867 cout << " Solidity = " << Solidity << endl;
868 cout << " Max Brake Power = " << MaxBrakePower/hptoftlbssec << " HP" << endl;
869 cout << " Gear Loss = " << GearLoss/hptoftlbssec << " HP" << endl;
870 cout << " Gear Moment = " << GearMoment << endl;
872 switch (ControlMap) {
873 case eTailCtrl: ControlMapName = "Tail Rotor"; break;
874 case eTandemCtrl: ControlMapName = "Tandem Rotor"; break;
875 default: ControlMapName = "Main Rotor";
877 cout << " Control Mapping = " << ControlMapName << endl;
881 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
882 if (from == 0) cout << "Instantiated: FGRotor" << endl;
883 if (from == 1) cout << "Destroyed: FGRotor" << endl;
885 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
887 if (debug_lvl & 8 ) { // Runtime state variables
889 if (debug_lvl & 16) { // Sanity checking
891 if (debug_lvl & 64) {
892 if (from == 0) { // Constructor
893 cout << IdSrc << endl;
894 cout << IdHdr << endl;
901 } // namespace JSBSim