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
40 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
49 #include "models/FGPropagate.h"
50 #include "models/FGAtmosphere.h"
51 #include "models/FGAuxiliary.h"
52 #include "models/FGMassBalance.h"
54 #include "input_output/FGXMLElement.h"
61 static const char *IdSrc = "$Id: FGRotor.cpp,v 1.12 2011/03/10 01:35:25 dpculp Exp $";
62 static const char *IdHdr = ID_ROTOR;
64 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
66 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
68 static inline double sqr(double x) { return x*x; }
70 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
72 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
75 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
79 FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
80 : FGThruster(exec, rotor_element, num),
84 dt(0.0), rho(0.002356),
86 // configuration parameters
87 Radius(0.0), BladeNum(0),
89 Sense(1.0), NominalRPM(0.0), ExternalRPM(0), RPMdefinition(0), ExtRPMsource(NULL),
91 BladeChord(0.0), LiftCurveSlope(0.0), BladeTwist(0.0), HingeOffset(0.0),
92 BladeFlappingMoment(0.0), BladeMassMoment(0.0), PolarMoment(0.0),
96 GroundEffectExp(0.0), GroundEffectShift(0.0),
99 LockNumberByRho(0.0), Solidity(0.0),
102 RPM(0.0), Omega(0.0),
105 a0(0.0), a_1(0.0), b_1(0.0), a_dw(0.0), a1s(0.0), b1s(0.0),
107 H_drag(0.0), J_side(0.0), Torque(0.0), C_T(0.0),
109 lambda(-0.001), mu(0.0), nu(0.001), v_induced(0.0),
110 theta_downwash(0.0), phi_downwash(0.0),
113 ControlMap(eMainCtrl),
114 CollectiveCtrl(0.0), LateralCtrl(0.0), LongitudinalCtrl(0.0),
115 BrakeCtrlNorm(0.0), MaxBrakePower(0.0),
117 // free-wheeling-unit (FWU)
118 FreeWheelPresent(0), FreeWheelThresh(0.0), FreeWheelTransmission(0.0)
121 FGColumnVector3 location(0.0, 0.0, 0.0), orientation(0.0, 0.0, 0.0);
122 Element *thruster_element;
124 // initialise/set remaining variables
125 SetTransformType(FGForce::tCustom);
126 PropertyManager = exec->GetPropertyManager();
130 dt = exec->GetDeltaT();
131 for (int i=0; i<5; i++) R[i] = 0.0;
132 for (int i=0; i<5; i++) B[i] = 0.0;
135 thruster_element = rotor_element->GetParent()->FindElement("sense");
136 if (thruster_element) {
137 double s = thruster_element->GetDataAsNumber();
139 Sense = -1.0; // 'CW' as seen from above
140 } else if (s < 0.1) {
141 Sense = 0.0; // 'coaxial'
143 Sense = 1.0; // 'CCW' as seen from above
147 thruster_element = rotor_element->GetParent()->FindElement("location");
148 if (thruster_element) {
149 location = thruster_element->FindElementTripletConvertTo("IN");
151 cerr << "No thruster location found." << endl;
154 thruster_element = rotor_element->GetParent()->FindElement("orient");
155 if (thruster_element) {
156 orientation = thruster_element->FindElementTripletConvertTo("RAD");
158 cerr << "No thruster orientation found." << endl;
161 SetLocation(location);
162 SetAnglesToBody(orientation);
163 InvTransform = Transform().Transposed();
166 ControlMap = eMainCtrl;
167 if (rotor_element->FindElement("controlmap")) {
168 string cm = rotor_element->FindElementValue("controlmap");
171 ControlMap = eTailCtrl;
172 } else if (cm == "TANDEM") {
173 ControlMap = eTandemCtrl;
175 cerr << "# found unknown controlmap: '" << cm << "' using main rotor config." << endl;
179 // ExternalRPM -- is the RPM dictated ?
180 if (rotor_element->FindElement("ExternalRPM")) {
182 RPMdefinition = (int) rotor_element->FindElementValueAsNumber("ExternalRPM");
185 // configure the rotor parameters
186 Configure(rotor_element);
188 // shaft representation - a rather simple transform,
189 // but using a matrix is safer.
190 TboToHsr.InitMatrix( 0.0, 0.0, 1.0,
193 HsrToTbo = TboToHsr.Transposed();
195 // smooth out jumps in hagl reported, otherwise the ground effect
196 // calculation would cause jumps too. 1Hz seems sufficient.
197 damp_hagl = Filter(1.0,dt);
199 // avoid too abrupt changes in power transmission
200 FreeWheelLag = Filter(200.0,dt);
202 // enable import-export
209 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
216 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
218 // 5in1: value-fetch-convert-default-return function
220 double FGRotor::ConfigValueConv( Element* el, const string& ename, double default_val,
221 const string& unit, bool tell)
225 double val = default_val;
227 string pname = "*No parent element*";
230 e = el->FindElement(ename);
231 pname = el->GetName();
236 val = e->GetDataAsNumber();
238 val = el->FindElementValueAsNumberConvertTo(ename,unit);
242 cerr << pname << ": missing element '" << ename <<
243 "' using estimated value: " << default_val << endl;
250 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
252 double FGRotor::ConfigValue(Element* el, const string& ename, double default_val, bool tell)
254 return ConfigValueConv(el, ename, default_val, "", tell);
257 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
259 // 1. read configuration and try to fill holes, ymmv
260 // 2. calculate derived parameters
261 void FGRotor::Configure(Element* rotor_element)
265 const bool yell = true;
266 const bool silent = false;
269 Radius = 0.5 * ConfigValueConv(rotor_element, "diameter", 42.0, "FT", yell);
270 Radius = Constrain(1e-3, Radius, 1e9);
272 BladeNum = (int) ConfigValue(rotor_element, "numblades", 3 , yell);
274 GearRatio = ConfigValue(rotor_element, "gearratio", 1.0, yell);
276 // make sure that v_tip (omega*r) is below 0.7mach ~ 750ft/s
277 estimate = (750.0/Radius)/(2.0*M_PI) * 60.0; // 7160/Radius
278 NominalRPM = ConfigValue(rotor_element, "nominalrpm", estimate, yell);
280 estimate = Constrain(0.07, 2.0/Radius , 0.14); // guess solidity
281 estimate = estimate * M_PI*Radius/BladeNum;
282 BladeChord = ConfigValueConv(rotor_element, "chord", estimate, "FT", yell);
284 LiftCurveSlope = ConfigValue(rotor_element, "liftcurveslope", 6.0); // "1/RAD"
285 BladeTwist = ConfigValueConv(rotor_element, "twist", -0.17, "RAD");
287 HingeOffset = ConfigValueConv(rotor_element, "hingeoffset", 0.05 * Radius, "FT" );
289 estimate = sqr(BladeChord) * sqr(Radius - HingeOffset) * 0.57;
290 BladeFlappingMoment = ConfigValueConv(rotor_element, "flappingmoment", estimate, "SLUG*FT2");
291 BladeFlappingMoment = Constrain(1.0e-6, BladeFlappingMoment, 1e9);
293 // guess mass from moment of a thin stick, and multiply by the blades cg distance
294 estimate = ( 3.0 * BladeFlappingMoment / sqr(Radius) ) * (0.45 * Radius) ;
295 BladeMassMoment = ConfigValue(rotor_element, "massmoment", estimate); // unit is slug-ft
296 BladeMassMoment = Constrain(0.001, BladeMassMoment, 1e9);
298 estimate = 1.1 * BladeFlappingMoment * BladeNum;
299 PolarMoment = ConfigValueConv(rotor_element, "polarmoment", estimate, "SLUG*FT2");
300 PolarMoment = Constrain(1e-6, PolarMoment, 1e9);
302 // "inflowlag" is treated further down.
304 TipLossB = ConfigValue(rotor_element, "tiplossfactor", 1.0, silent);
306 estimate = 0.01 * PolarMoment ; // guesses for huey, bo105 20-30hp
307 MaxBrakePower = ConfigValueConv(rotor_element, "maxbrakepower", estimate, "HP");
308 MaxBrakePower *= hptoftlbssec;
311 if (rotor_element->FindElement("cgroundeffect")) {
313 cge = rotor_element->FindElementValueAsNumber("cgroundeffect");
314 cge = Constrain(1e-9, cge, 1.0);
315 gee = 1.0 / ( 2.0*Radius * cge );
316 cerr << "# *** 'cgroundeffect' is defunct." << endl;
317 cerr << "# *** use 'groundeffectexp' with: " << gee << endl;
320 GroundEffectExp = ConfigValue(rotor_element, "groundeffectexp", 0.0);
321 GroundEffectShift = ConfigValueConv(rotor_element, "groundeffectshift", 0.0, "FT");
323 // handle optional free-wheeling-unit (FWU)
324 FreeWheelPresent = 0;
325 FreeWheelTransmission = 1.0;
326 if (rotor_element->FindElement("freewheelthresh")) {
327 FreeWheelThresh = rotor_element->FindElementValueAsNumber("freewheelthresh");
328 if (FreeWheelThresh > 1.0) {
329 FreeWheelPresent = 1;
330 FreeWheelTransmission = 0.0;
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(1.0e-6, InflowLag, 2.0);
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)
419 ct_t0 = (1.0/3.0*B[3] + 1.0/2.0 * TipLossB*mu2 - 4.0/(9.0*M_PI) * mu*mu2 ) * theta_0;
420 ct_t1 = (1.0/4.0*B[4] + 1.0/4.0 * B[2]*mu2) * BladeTwist;
422 ct_l = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda; // first time
424 c0 = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1) * Solidity;
425 c0 = c0 / ( 2.0 * sqrt( sqr(mu) + sqr(lambda) ) + 1e-15);
427 // replacement for /SH79/ eqn(26).
428 // ref: dnu/dt = 1/tau ( Ct / (2*sqrt(mu^2+lambda^2)) - nu )
429 // taking mu and lambda constant, this integrates to
431 nu = flow_scale * ((nu - c0) * exp(-dt/InflowLag) + c0);
433 // now from nu to lambda, C_T, and Thrust
435 lambda = Ww/(Omega*Radius) - nu; // /SH79/ eqn(25)
437 ct_l = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda;
439 ct_over_sigma = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1); // /SH79/ eqn(27)
441 Thrust = BladeNum*BladeChord*Radius*rho*sqr(Omega*Radius) * ct_over_sigma;
443 C_T = ct_over_sigma * Solidity;
444 v_induced = nu * (Omega*Radius);
449 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
451 // The coning angle doesn't apply for teetering rotors, but calculating
452 // doesn't hurt. /SH79/ eqn(29)
454 void FGRotor::calc_coning_angle(double theta_0)
456 double lock_gamma = LockNumberByRho * rho;
458 double a0_l = (1.0/6.0 + 0.04 * mu*mu*mu) * lambda;
459 double a0_t0 = (1.0/8.0 + 1.0/8.0 * mu*mu) * theta_0;
460 double a0_t1 = (1.0/10.0 + 1.0/12.0 * mu*mu) * BladeTwist;
461 a0 = lock_gamma * ( a0_l + a0_t0 + a0_t1);
465 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
467 // Flapping angles relative to control axes /SH79/ eqn(32)
469 void FGRotor::calc_flapping_angles(double theta_0, const FGColumnVector3 &pqr_fus_w)
471 double lock_gamma = LockNumberByRho * rho;
474 double mu2_2 = sqr(mu)/2.0;
475 double t075 = theta_0 + 0.75 * BladeTwist; // common approximation for rectangular blades
477 a_1 = 1.0/(1.0 - mu2_2) * (
478 (2.0*lambda + (8.0/3.0)*t075)*mu
479 + pqr_fus_w(eP)/Omega
480 - 16.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
483 b_1 = 1.0/(1.0 + mu2_2) * (
485 - pqr_fus_w(eQ)/Omega
486 - 16.0 * pqr_fus_w(eP)/(lock_gamma*Omega)
489 // used in force calc
490 a_dw = 1.0/(1.0 - mu2_2) * (
491 (2.0*lambda + (8.0/3.0)*t075)*mu
492 - 24.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
493 * ( 1.0 - ( 0.29 * t075 / (C_T/Solidity) ) )
499 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
503 void FGRotor::calc_drag_and_side_forces(double theta_0)
505 double cy_over_sigma;
506 double t075 = theta_0 + 0.75 * BladeTwist;
508 H_drag = Thrust * a_dw;
511 0.75*b_1*lambda - 1.5*a0*mu*lambda + 0.25*a_1*b_1*mu
512 - a0*a_1*sqr(mu) + (1.0/6.0)*a0*a_1
513 - (0.75*mu*a0 - (1.0/3.0)*b_1 - 0.5*sqr(mu)*b_1)*t075
515 cy_over_sigma *= LiftCurveSlope/2.0;
517 J_side = BladeNum * BladeChord * Radius * rho * sqr(Omega*Radius) * cy_over_sigma;
522 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
524 // Simplified version of /SH79/ eqn(36). Uses an estimate for blade drag
525 // (a new config parameter to come...).
526 // From "Bramwell's Helicopter Dynamics", second edition, eqn(3.43) and (3.44)
528 void FGRotor::calc_torque(double theta_0)
530 // estimate blade drag
531 double delta_dr = 0.009 + 0.3*sqr(6.0*C_T/(LiftCurveSlope*Solidity));
533 Torque = rho * BladeNum * BladeChord * delta_dr * sqr(Omega*Radius) * R[2] *
534 (1.0+4.5*sqr(mu))/8.0
535 - (Thrust*lambda + H_drag*mu)*Radius;
540 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
542 // transform rotor forces from control axes to shaft axes, and express
543 // in body axes /SH79/ eqn(40,41)
545 FGColumnVector3 FGRotor::body_forces(double a_ic, double b_ic)
548 - H_drag*cos(beta_orient) - J_side*sin(beta_orient) + Thrust*b_ic,
549 - H_drag*sin(beta_orient) + J_side*cos(beta_orient) + Thrust*a_ic,
552 return HsrToTbo * F_s;
555 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
557 // calculates the additional moments due to hinge offset and handles
560 FGColumnVector3 FGRotor::body_moments(double a_ic, double b_ic)
562 FGColumnVector3 M_s, M_hub, M_h;
565 // cyclic flapping relative to shaft axes /SH79/ eqn(43)
566 a1s = a_1*cos(beta_orient) + b_1*sin(beta_orient) - b_ic;
567 b1s = b_1*cos(beta_orient) - a_1*sin(beta_orient) + a_ic;
569 mf = 0.5 * HingeOffset * BladeNum * Omega*Omega * BladeMassMoment;
573 M_s(eN) = Torque * Sense ;
575 return HsrToTbo * M_s;
578 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
580 void FGRotor::CalcStatePart1(void)
582 double A_IC; // lateral (roll) control in radians
583 double B_IC; // longitudinal (pitch) control in radians
584 double theta_col; // rotor collective pitch in radians
588 FGColumnVector3 UVW_h, PQR_h;
589 FGColumnVector3 vHub_ca, avFus_ca;
591 double h_agl_ft, filtered_hagl = 0.0;
592 double ge_factor = 1.0;
594 // fetch needed values from environment
595 Vt = fdmex->GetAuxiliary()->GetVt(); // total vehicle velocity including wind
596 dt = fdmex->GetDeltaT();
597 rho = fdmex->GetAtmosphere()->GetDensity(); // slugs/ft^3.
598 UVW_h = fdmex->GetAuxiliary()->GetAeroUVW();
599 PQR_h = fdmex->GetAuxiliary()->GetAeroPQR();
600 h_agl_ft = fdmex->GetPropagate()->GetDistanceAGL();
601 // update InvTransform, the rotor orientation could have been altered
602 InvTransform = Transform().Transposed();
604 // handle RPM requirements, calc omega.
605 if (ExternalRPM && ExtRPMsource) {
606 RPM = ExtRPMsource->getDoubleValue() / GearRatio;
609 if (RPM < 1.0) { // kludge, otherwise calculations go bananas
613 Omega = (RPM/60.0)*2.0*M_PI;
615 // set control inputs
617 B_IC = LongitudinalCtrl;
618 theta_col = CollectiveCtrl;
621 if (GroundEffectExp > 1e-5) {
622 if (h_agl_ft<0.0) h_agl_ft = 0.0; // clamp
623 filtered_hagl = damp_hagl.execute(h_agl_ft) + GroundEffectShift;
624 // actual/nominal factor avoids absurd scales at startup
625 ge_factor -= exp(-filtered_hagl*GroundEffectExp) * (RPM / NominalRPM);
626 if (ge_factor<0.5) ge_factor=0.5; // clamp
629 // all set, start calculations
631 vHub_ca = hub_vel_body2ca(UVW_h, PQR_h, A_IC, B_IC);
633 avFus_ca = fus_angvel_body2ca(PQR_h);
635 calc_flow_and_thrust(theta_col, vHub_ca(eU), vHub_ca(eW), ge_factor);
637 calc_coning_angle(theta_col);
639 calc_flapping_angles(theta_col, avFus_ca);
641 calc_drag_and_side_forces(theta_col);
643 calc_torque(theta_col);
645 // Fixme: only valid for a 'decent' rotor
646 theta_downwash = atan2( - UVW_h(eU), v_induced - UVW_h(eW));
647 phi_downwash = atan2( UVW_h(eV), v_induced - UVW_h(eW));
649 vFn = body_forces(A_IC, B_IC);
650 vMn = Transform() * body_moments(A_IC, B_IC);
654 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
656 void FGRotor::CalcStatePart2(double PowerAvailable)
660 double ExcessTorque = PowerAvailable / Omega;
661 double deltaOmega = ExcessTorque / PolarMoment * dt;
662 RPM += deltaOmega/(2.0*M_PI) * 60.0;
663 if (RPM < 0.0) RPM = 0.0; // Engine won't turn backwards
667 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
669 // Simulation of a free-wheeling-unit (FWU). Might need improvements.
671 void FGRotor::calc_freewheel_state(double p_source, double p_load) {
673 // engine is off/detached, release.
675 FreeWheelTransmission = 0.0;
679 // engine is driving the rotor, engage.
680 if (p_source >= p_load) {
681 FreeWheelTransmission = 1.0;
685 // releases if engine is detached, but stays calm if
686 // the load changes due to rotor dynamics.
687 if (p_source > 0.0 && p_load/(p_source+0.1) > FreeWheelThresh ) {
688 FreeWheelTransmission = 0.0;
695 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
697 double FGRotor::Calculate(double EnginePower)
704 PowerRequired = Torque * Omega + BrakeCtrlNorm * MaxBrakePower;
706 if (FreeWheelPresent) {
707 calc_freewheel_state(EnginePower * ClutchCtrlNorm, PowerRequired);
708 FWmult = FreeWheelLag.execute(FreeWheelTransmission);
711 DeltaPower = EnginePower * ClutchCtrlNorm * FWmult - PowerRequired;
713 CalcStatePart2(DeltaPower);
719 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
722 bool FGRotor::BindModel(void)
724 string property_name, base_property_name;
725 base_property_name = CreateIndexedPropertyName("propulsion/engine", EngineNum);
727 property_name = base_property_name + "/rotor-rpm";
728 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetRPM );
730 property_name = base_property_name + "/x-engine-rpm"; // used for RPM eXchange
731 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetEngineRPM );
733 property_name = base_property_name + "/rotor-thrust-lbs"; // might be redundant - check!
734 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetThrust );
736 property_name = base_property_name + "/a0-rad";
737 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA0 );
739 property_name = base_property_name + "/a1-rad";
740 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA1 );
742 property_name = base_property_name + "/b1-rad";
743 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetB1 );
745 property_name = base_property_name + "/inflow-ratio";
746 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLambda );
748 property_name = base_property_name + "/advance-ratio";
749 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetMu );
751 property_name = base_property_name + "/induced-inflow-ratio";
752 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetNu );
754 property_name = base_property_name + "/vi-fps";
755 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetVi );
757 property_name = base_property_name + "/thrust-coefficient";
758 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCT );
760 property_name = base_property_name + "/torque-lbsft";
761 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetTorque );
763 property_name = base_property_name + "/theta-downwash-rad";
764 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetThetaDW );
766 property_name = base_property_name + "/phi-downwash-rad";
767 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetPhiDW );
769 switch (ControlMap) {
771 property_name = base_property_name + "/antitorque-ctrl-rad";
772 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
775 property_name = base_property_name + "/tail-collective-ctrl-rad";
776 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
777 property_name = base_property_name + "/lateral-ctrl-rad";
778 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
779 property_name = base_property_name + "/longitudinal-ctrl-rad";
780 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
782 default: // eMainCtrl
783 property_name = base_property_name + "/collective-ctrl-rad";
784 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
785 property_name = base_property_name + "/lateral-ctrl-rad";
786 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
787 property_name = base_property_name + "/longitudinal-ctrl-rad";
788 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
791 property_name = base_property_name + "/brake-ctrl-norm";
792 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetBrakeCtrl, &FGRotor::SetBrakeCtrl);
793 property_name = base_property_name + "/free-wheel-transmission";
794 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetFreeWheelTransmission);
797 if (RPMdefinition == -1) {
798 property_name = base_property_name + "/x-rpm-dict";
799 ExtRPMsource = PropertyManager->GetNode(property_name, true);
800 } else if (RPMdefinition >= 0 && RPMdefinition != EngineNum) {
801 string ipn = CreateIndexedPropertyName("propulsion/engine", RPMdefinition);
802 property_name = ipn + "/x-engine-rpm";
803 ExtRPMsource = PropertyManager->GetNode(property_name, false);
804 if (! ExtRPMsource) {
805 cerr << "# Warning: Engine number " << EngineNum << "." << endl;
806 cerr << "# No 'x-engine-rpm' property found for engine " << RPMdefinition << "." << endl;
807 cerr << "# Please check order of engine definitons." << endl;
810 cerr << "# Engine number " << EngineNum;
811 cerr << ", given ExternalRPM value '" << RPMdefinition << "' unhandled." << endl;
818 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
820 string FGRotor::GetThrusterLabels(int id, string delimeter)
825 buf << Name << " RPM (engine " << id << ")";
831 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
833 string FGRotor::GetThrusterValues(int id, string delimeter)
844 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
845 // The bitmasked value choices are as follows:
846 // unset: In this case (the default) JSBSim would only print
847 // out the normally expected messages, essentially echoing
848 // the config files as they are read. If the environment
849 // variable is not set, debug_lvl is set to 1 internally
850 // 0: This requests JSBSim not to output any messages
852 // 1: This value explicity requests the normal JSBSim
854 // 2: This value asks for a message to be printed out when
855 // a class is instantiated
856 // 4: When this value is set, a message is displayed when a
857 // FGModel object executes its Run() method
858 // 8: When this value is set, various runtime state variables
859 // are printed out periodically
860 // 16: When set various parameters are sanity checked and
861 // a message is printed out when they go out of bounds
863 void FGRotor::Debug(int from)
865 string ControlMapName;
867 if (debug_lvl <= 0) return;
869 if (debug_lvl & 1) { // Standard console startup message output
870 if (from == 0) { // Constructor
871 cout << "\n Rotor Name: " << Name << endl;
872 cout << " Diameter = " << 2.0 * Radius << " ft." << endl;
873 cout << " Number of Blades = " << BladeNum << endl;
874 cout << " Gear Ratio = " << GearRatio << endl;
875 cout << " Sense = " << Sense << endl;
876 cout << " Nominal RPM = " << NominalRPM << endl;
879 if (RPMdefinition == -1) {
880 cout << " RPM is controlled externally" << endl;
882 cout << " RPM source set to engine " << RPMdefinition << endl;
886 cout << " Blade Chord = " << BladeChord << endl;
887 cout << " Lift Curve Slope = " << LiftCurveSlope << endl;
888 cout << " Blade Twist = " << BladeTwist << endl;
889 cout << " Hinge Offset = " << HingeOffset << endl;
890 cout << " Blade Flapping Moment = " << BladeFlappingMoment << endl;
891 cout << " Blade Mass Moment = " << BladeMassMoment << endl;
892 cout << " Polar Moment = " << PolarMoment << endl;
893 cout << " Inflow Lag = " << InflowLag << endl;
894 cout << " Tip Loss = " << TipLossB << endl;
895 cout << " Lock Number = " << LockNumberByRho * 0.002356 << " (SL)" << endl;
896 cout << " Solidity = " << Solidity << endl;
897 cout << " Max Brake Power = " << MaxBrakePower/hptoftlbssec << " HP" << endl;
899 switch (ControlMap) {
900 case eTailCtrl: ControlMapName = "Tail Rotor"; break;
901 case eTandemCtrl: ControlMapName = "Tandem Rotor"; break;
902 default: ControlMapName = "Main Rotor";
904 cout << " Control Mapping = " << ControlMapName << endl;
906 if (FreeWheelPresent) {
907 cout << " Free Wheel Threshold = " << FreeWheelThresh << endl;
909 cout << " No FWU present" << endl;
914 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
915 if (from == 0) cout << "Instantiated: FGRotor" << endl;
916 if (from == 1) cout << "Destroyed: FGRotor" << endl;
918 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
920 if (debug_lvl & 8 ) { // Runtime state variables
922 if (debug_lvl & 16) { // Sanity checking
924 if (debug_lvl & 64) {
925 if (from == 0) { // Constructor
926 cout << IdSrc << endl;
927 cout << IdHdr << endl;
934 } // namespace JSBSim