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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
48 #include "input_output/FGXMLElement.h"
49 #include "models/FGMassBalance.h"
53 using std::ostringstream;
58 static const char *IdSrc = "$Id: FGRotor.cpp,v 1.17 2011/09/24 14:26:46 jentron Exp $";
59 static const char *IdHdr = ID_ROTOR;
61 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
65 static inline double sqr(double x) { return x*x; }
67 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
69 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
72 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
76 FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
77 : FGThruster(exec, rotor_element, num),
78 rho(0.002356), // environment
79 Radius(0.0), BladeNum(0), // configuration parameters
80 Sense(1.0), NominalRPM(0.0), MinimalRPM(0.0), MaximalRPM(0.0),
81 ExternalRPM(0), RPMdefinition(0), ExtRPMsource(NULL),
82 BladeChord(0.0), LiftCurveSlope(0.0), BladeTwist(0.0), HingeOffset(0.0),
83 BladeFlappingMoment(0.0), BladeMassMoment(0.0), PolarMoment(0.0),
84 InflowLag(0.0), TipLossB(0.0),
85 GroundEffectExp(0.0), GroundEffectShift(0.0),
86 LockNumberByRho(0.0), Solidity(0.0), // derived parameters
87 RPM(0.0), Omega(0.0), // dynamic values
89 a0(0.0), a_1(0.0), b_1(0.0), a_dw(0.0),
91 H_drag(0.0), J_side(0.0), Torque(0.0), C_T(0.0),
92 lambda(-0.001), mu(0.0), nu(0.001), v_induced(0.0),
93 theta_downwash(0.0), phi_downwash(0.0),
94 ControlMap(eMainCtrl), // control
95 CollectiveCtrl(0.0), LateralCtrl(0.0), LongitudinalCtrl(0.0),
96 BrakeCtrlNorm(0.0), MaxBrakePower(0.0),
97 FreeWheelPresent(0), FreeWheelThresh(0.0), // free-wheeling-unit (FWU)
98 FreeWheelTransmission(0.0)
100 FGColumnVector3 location(0.0, 0.0, 0.0), orientation(0.0, 0.0, 0.0);
101 Element *thruster_element;
103 // initialise/set remaining variables
104 SetTransformType(FGForce::tCustom);
105 PropertyManager = exec->GetPropertyManager();
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;
114 thruster_element = rotor_element->GetParent()->FindElement("sense");
115 if (thruster_element) {
116 double s = thruster_element->GetDataAsNumber();
118 Sense = -1.0; // 'CW' as seen from above
119 } else if (s < 0.1) {
120 Sense = 0.0; // 'coaxial'
122 Sense = 1.0; // 'CCW' as seen from above
126 thruster_element = rotor_element->GetParent()->FindElement("location");
127 if (thruster_element) {
128 location = thruster_element->FindElementTripletConvertTo("IN");
130 cerr << "No thruster location found." << endl;
133 thruster_element = rotor_element->GetParent()->FindElement("orient");
134 if (thruster_element) {
135 orientation = thruster_element->FindElementTripletConvertTo("RAD");
137 cerr << "No thruster orientation found." << endl;
140 SetLocation(location);
141 SetAnglesToBody(orientation);
142 InvTransform = Transform().Transposed();
145 ControlMap = eMainCtrl;
146 if (rotor_element->FindElement("controlmap")) {
147 string cm = rotor_element->FindElementValue("controlmap");
150 ControlMap = eTailCtrl;
151 } else if (cm == "TANDEM") {
152 ControlMap = eTandemCtrl;
154 cerr << "# found unknown controlmap: '" << cm << "' using main rotor config." << endl;
158 // ExternalRPM -- is the RPM dictated ?
159 if (rotor_element->FindElement("ExternalRPM")) {
161 RPMdefinition = (int) rotor_element->FindElementValueAsNumber("ExternalRPM");
164 // configure the rotor parameters
165 Configure(rotor_element);
167 // shaft representation - a rather simple transform,
168 // but using a matrix is safer.
169 TboToHsr.InitMatrix( 0.0, 0.0, 1.0,
172 HsrToTbo = TboToHsr.Transposed();
174 // smooth out jumps in hagl reported, otherwise the ground effect
175 // calculation would cause jumps too. 1Hz seems sufficient.
176 damp_hagl = Filter(1.0, dt);
178 // avoid too abrupt changes in power transmission
179 FreeWheelLag = Filter(200.0,dt);
181 // enable import-export
188 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
195 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
197 // 5in1: value-fetch-convert-default-return function
199 double FGRotor::ConfigValueConv( Element* el, const string& ename, double default_val,
200 const string& unit, bool tell)
204 double val = default_val;
206 string pname = "*No parent element*";
209 e = el->FindElement(ename);
210 pname = el->GetName();
215 val = e->GetDataAsNumber();
217 val = el->FindElementValueAsNumberConvertTo(ename,unit);
221 cerr << pname << ": missing element '" << ename <<
222 "' using estimated value: " << default_val << endl;
229 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
231 double FGRotor::ConfigValue(Element* el, const string& ename, double default_val, bool tell)
233 return ConfigValueConv(el, ename, default_val, "", tell);
236 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
238 // 1. read configuration and try to fill holes, ymmv
239 // 2. calculate derived parameters
240 void FGRotor::Configure(Element* rotor_element)
244 const bool yell = true;
245 const bool silent = false;
248 Radius = 0.5 * ConfigValueConv(rotor_element, "diameter", 42.0, "FT", yell);
249 Radius = Constrain(1e-3, Radius, 1e9);
251 BladeNum = (int) ConfigValue(rotor_element, "numblades", 3 , yell);
253 GearRatio = ConfigValue(rotor_element, "gearratio", 1.0, yell);
255 // make sure that v_tip (omega*r) is below 0.7mach ~ 750ft/s
256 estimate = (750.0/Radius)/(2.0*M_PI) * 60.0; // 7160/Radius
257 NominalRPM = ConfigValue(rotor_element, "nominalrpm", estimate, yell);
258 NominalRPM = Constrain(2.0, NominalRPM, 1e9);
260 MinimalRPM = ConfigValue(rotor_element, "minrpm", 1.0);
261 MinimalRPM = Constrain(1.0, MinimalRPM, NominalRPM - 1.0);
263 MaximalRPM = ConfigValue(rotor_element, "maxrpm", 2.0*NominalRPM);
264 MaximalRPM = Constrain(NominalRPM, MaximalRPM, 1e9);
266 estimate = Constrain(0.07, 2.0/Radius , 0.14); // guess solidity
267 estimate = estimate * M_PI*Radius/BladeNum;
268 BladeChord = ConfigValueConv(rotor_element, "chord", estimate, "FT", yell);
270 LiftCurveSlope = ConfigValue(rotor_element, "liftcurveslope", 6.0); // "1/RAD"
271 BladeTwist = ConfigValueConv(rotor_element, "twist", -0.17, "RAD");
273 HingeOffset = ConfigValueConv(rotor_element, "hingeoffset", 0.05 * Radius, "FT" );
275 estimate = sqr(BladeChord) * sqr(Radius - HingeOffset) * 0.57;
276 BladeFlappingMoment = ConfigValueConv(rotor_element, "flappingmoment", estimate, "SLUG*FT2");
277 BladeFlappingMoment = Constrain(1.0e-6, BladeFlappingMoment, 1e9);
279 // guess mass from moment of a thin stick, and multiply by the blades cg distance
280 estimate = ( 3.0 * BladeFlappingMoment / sqr(Radius) ) * (0.45 * Radius) ;
281 BladeMassMoment = ConfigValue(rotor_element, "massmoment", estimate); // unit is slug-ft
282 BladeMassMoment = Constrain(0.001, BladeMassMoment, 1e9);
284 estimate = 1.1 * BladeFlappingMoment * BladeNum;
285 PolarMoment = ConfigValueConv(rotor_element, "polarmoment", estimate, "SLUG*FT2");
286 PolarMoment = Constrain(1e-6, PolarMoment, 1e9);
288 // "inflowlag" is treated further down.
290 TipLossB = ConfigValue(rotor_element, "tiplossfactor", 1.0, silent);
292 estimate = 0.01 * PolarMoment ; // guesses for huey, bo105 20-30hp
293 MaxBrakePower = ConfigValueConv(rotor_element, "maxbrakepower", estimate, "HP");
294 MaxBrakePower *= hptoftlbssec;
297 if (rotor_element->FindElement("cgroundeffect")) {
299 cge = rotor_element->FindElementValueAsNumber("cgroundeffect");
300 cge = Constrain(1e-9, cge, 1.0);
301 gee = 1.0 / ( 2.0*Radius * cge );
302 cerr << "# *** 'cgroundeffect' is defunct." << endl;
303 cerr << "# *** use 'groundeffectexp' with: " << gee << endl;
306 GroundEffectExp = ConfigValue(rotor_element, "groundeffectexp", 0.0);
307 GroundEffectShift = ConfigValueConv(rotor_element, "groundeffectshift", 0.0, "FT");
309 // handle optional free-wheeling-unit (FWU)
310 FreeWheelPresent = 0;
311 FreeWheelTransmission = 1.0;
312 if (rotor_element->FindElement("freewheelthresh")) {
313 FreeWheelThresh = rotor_element->FindElementValueAsNumber("freewheelthresh");
314 if (FreeWheelThresh > 1.0) {
315 FreeWheelPresent = 1;
316 FreeWheelTransmission = 0.0;
320 // precalc often used powers
321 R[0]=1.0; R[1]=Radius; R[2]=R[1]*R[1]; R[3]=R[2]*R[1]; R[4]=R[3]*R[1];
322 B[0]=1.0; B[1]=TipLossB; B[2]=B[1]*B[1]; B[3]=B[2]*B[1]; B[4]=B[3]*B[1];
324 // derived parameters
325 LockNumberByRho = LiftCurveSlope * BladeChord * R[4] / BladeFlappingMoment;
326 Solidity = BladeNum * BladeChord / (M_PI * Radius);
328 // estimate inflow lag, see /GE49/ eqn(1)
329 double omega_tmp = (NominalRPM/60.0)*2.0*M_PI;
330 estimate = 16.0/(LockNumberByRho*rho * omega_tmp ); // 16/(gamma*Omega)
331 // printf("# Est. InflowLag: %f\n", estimate);
332 InflowLag = ConfigValue(rotor_element, "inflowlag", estimate, yell);
333 InflowLag = Constrain(1.0e-6, InflowLag, 2.0);
338 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
340 // calculate control-axes components of total airspeed at the hub.
341 // sets rotor orientation angle (beta) as side effect. /SH79/ eqn(19-22)
343 FGColumnVector3 FGRotor::hub_vel_body2ca( const FGColumnVector3 &uvw,
344 const FGColumnVector3 &pqr,
345 double a_ic, double b_ic)
347 FGColumnVector3 v_r, v_shaft, v_w;
350 pos = fdmex->GetMassBalance()->StructuralToBody(GetActingLocation());
353 v_shaft = TboToHsr * InvTransform * v_r;
355 beta_orient = atan2(v_shaft(eV),v_shaft(eU));
357 v_w(eU) = v_shaft(eU)*cos(beta_orient) + v_shaft(eV)*sin(beta_orient);
359 v_w(eW) = v_shaft(eW) - b_ic*v_shaft(eU) - a_ic*v_shaft(eV);
364 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
366 // express fuselage angular velocity in control axes /SH79/ eqn(30,31)
368 FGColumnVector3 FGRotor::fus_angvel_body2ca( const FGColumnVector3 &pqr)
370 FGColumnVector3 av_s_fus, av_w_fus;
373 // av_s_fus = BodyToShaft * pqr; /SH79/
374 // BodyToShaft = TboToHsr * InvTransform
375 av_s_fus = TboToHsr * InvTransform * pqr;
377 av_w_fus(eP)= av_s_fus(eP)*cos(beta_orient) + av_s_fus(eQ)*sin(beta_orient);
378 av_w_fus(eQ)= - av_s_fus(eP)*sin(beta_orient) + av_s_fus(eQ)*cos(beta_orient);
379 av_w_fus(eR)= av_s_fus(eR);
385 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
387 // The calculation is a bit tricky because thrust depends on induced velocity,
390 // The flow_scale parameter (ranging from 0.5-1.0) is used to approximate a
391 // reduction of inflow if the helicopter is close to the ground, yielding to
392 // higher thrust, see /TA77/ eqn(10a).
394 void FGRotor::calc_flow_and_thrust( double theta_0, double Uw, double Ww,
398 double ct_over_sigma = 0.0;
399 double c0, ct_l, ct_t0, ct_t1;
402 mu = Uw/(Omega*Radius); // /SH79/ eqn(24)
403 if (mu > 0.7) mu = 0.7;
406 ct_t0 = (1.0/3.0*B[3] + 1.0/2.0 * TipLossB*mu2 - 4.0/(9.0*M_PI) * mu*mu2 ) * theta_0;
407 ct_t1 = (1.0/4.0*B[4] + 1.0/4.0 * B[2]*mu2) * BladeTwist;
409 ct_l = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda; // first time
411 c0 = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1) * Solidity;
412 c0 = c0 / ( 2.0 * sqrt( sqr(mu) + sqr(lambda) ) + 1e-15);
414 // replacement for /SH79/ eqn(26).
415 // ref: dnu/dt = 1/tau ( Ct / (2*sqrt(mu^2+lambda^2)) - nu )
416 // taking mu and lambda constant, this integrates to
418 nu = flow_scale * ((nu - c0) * exp(-dt/InflowLag) + c0);
420 // now from nu to lambda, C_T, and Thrust
422 lambda = Ww/(Omega*Radius) - nu; // /SH79/ eqn(25)
424 ct_l = (1.0/2.0*B[2] + 1.0/4.0 * mu2) * lambda;
426 ct_over_sigma = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1); // /SH79/ eqn(27)
428 Thrust = BladeNum*BladeChord*Radius*rho*sqr(Omega*Radius) * ct_over_sigma;
430 C_T = ct_over_sigma * Solidity;
431 v_induced = nu * (Omega*Radius);
436 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
438 // The coning angle doesn't apply for teetering rotors, but calculating
439 // doesn't hurt. /SH79/ eqn(29)
441 void FGRotor::calc_coning_angle(double theta_0)
443 double lock_gamma = LockNumberByRho * rho;
445 double a0_l = (1.0/6.0 + 0.04 * mu*mu*mu) * lambda;
446 double a0_t0 = (1.0/8.0 + 1.0/8.0 * mu*mu) * theta_0;
447 double a0_t1 = (1.0/10.0 + 1.0/12.0 * mu*mu) * BladeTwist;
448 a0 = lock_gamma * ( a0_l + a0_t0 + a0_t1);
452 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
454 // Flapping angles relative to control axes /SH79/ eqn(32)
456 void FGRotor::calc_flapping_angles(double theta_0, const FGColumnVector3 &pqr_fus_w)
458 double lock_gamma = LockNumberByRho * rho;
461 double mu2_2 = sqr(mu)/2.0;
462 double t075 = theta_0 + 0.75 * BladeTwist; // common approximation for rectangular blades
464 a_1 = 1.0/(1.0 - mu2_2) * (
465 (2.0*lambda + (8.0/3.0)*t075)*mu
466 + pqr_fus_w(eP)/Omega
467 - 16.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
470 b_1 = 1.0/(1.0 + mu2_2) * (
472 - pqr_fus_w(eQ)/Omega
473 - 16.0 * pqr_fus_w(eP)/(lock_gamma*Omega)
476 // used in force calc
477 a_dw = 1.0/(1.0 - mu2_2) * (
478 (2.0*lambda + (8.0/3.0)*t075)*mu
479 - 24.0 * pqr_fus_w(eQ)/(lock_gamma*Omega)
480 * ( 1.0 - ( 0.29 * t075 / (C_T/Solidity) ) )
486 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
490 void FGRotor::calc_drag_and_side_forces(double theta_0)
492 double cy_over_sigma;
493 double t075 = theta_0 + 0.75 * BladeTwist;
495 H_drag = Thrust * a_dw;
498 0.75*b_1*lambda - 1.5*a0*mu*lambda + 0.25*a_1*b_1*mu
499 - a0*a_1*sqr(mu) + (1.0/6.0)*a0*a_1
500 - (0.75*mu*a0 - (1.0/3.0)*b_1 - 0.5*sqr(mu)*b_1)*t075
502 cy_over_sigma *= LiftCurveSlope/2.0;
504 J_side = BladeNum * BladeChord * Radius * rho * sqr(Omega*Radius) * cy_over_sigma;
509 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
511 // Simplified version of /SH79/ eqn(36). Uses an estimate for blade drag
512 // (a new config parameter to come...).
513 // From "Bramwell's Helicopter Dynamics", second edition, eqn(3.43) and (3.44)
515 void FGRotor::calc_torque(double theta_0)
517 // estimate blade drag
518 double delta_dr = 0.009 + 0.3*sqr(6.0*C_T/(LiftCurveSlope*Solidity));
520 Torque = rho * BladeNum * BladeChord * delta_dr * sqr(Omega*Radius) * R[2] *
521 (1.0+4.5*sqr(mu))/8.0
522 - (Thrust*lambda + H_drag*mu)*Radius;
527 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
529 // transform rotor forces from control axes to shaft axes, and express
530 // in body axes /SH79/ eqn(40,41)
532 FGColumnVector3 FGRotor::body_forces(double a_ic, double b_ic)
535 - H_drag*cos(beta_orient) - J_side*sin(beta_orient) + Thrust*b_ic,
536 - H_drag*sin(beta_orient) + J_side*cos(beta_orient) + Thrust*a_ic,
539 return HsrToTbo * F_s;
542 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
544 // calculates the additional moments due to hinge offset and handles
547 FGColumnVector3 FGRotor::body_moments(double a_ic, double b_ic)
549 FGColumnVector3 M_s, M_hub, M_h;
552 // cyclic flapping relative to shaft axes /SH79/ eqn(43)
553 a1s = a_1*cos(beta_orient) + b_1*sin(beta_orient) - b_ic;
554 b1s = b_1*cos(beta_orient) - a_1*sin(beta_orient) + a_ic;
556 mf = 0.5 * HingeOffset * BladeNum * Omega*Omega * BladeMassMoment;
560 M_s(eN) = Torque * Sense ;
562 return HsrToTbo * M_s;
565 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
567 void FGRotor::CalcStatePart1(void)
569 double A_IC; // lateral (roll) control in radians
570 double B_IC; // longitudinal (pitch) control in radians
571 double theta_col; // rotor collective pitch in radians
573 FGColumnVector3 vHub_ca, avFus_ca;
575 double filtered_hagl = 0.0;
576 double ge_factor = 1.0;
578 // fetch needed values from environment
579 rho = in.Density; // slugs/ft^3.
580 double h_agl_ft = in.H_agl;
581 // update InvTransform, the rotor orientation could have been altered
582 InvTransform = Transform().Transposed();
584 // handle RPM requirements, calc omega.
585 if (ExternalRPM && ExtRPMsource) {
586 RPM = ExtRPMsource->getDoubleValue() / GearRatio;
589 // MinimalRPM is always >= 1. MaximalRPM is always >= NominalRPM
590 RPM = Constrain(MinimalRPM, RPM, MaximalRPM);
592 Omega = (RPM/60.0)*2.0*M_PI;
594 // set control inputs
596 B_IC = LongitudinalCtrl;
597 theta_col = CollectiveCtrl;
600 if (GroundEffectExp > 1e-5) {
601 if (h_agl_ft<0.0) h_agl_ft = 0.0; // clamp
602 filtered_hagl = damp_hagl.execute(h_agl_ft) + GroundEffectShift;
603 // actual/nominal factor avoids absurd scales at startup
604 ge_factor -= exp(-filtered_hagl*GroundEffectExp) * (RPM / NominalRPM);
605 if (ge_factor<0.5) ge_factor=0.5; // clamp
608 // all set, start calculations
610 vHub_ca = hub_vel_body2ca(in.AeroUVW, in.AeroPQR, A_IC, B_IC);
612 avFus_ca = fus_angvel_body2ca(in.AeroPQR);
614 calc_flow_and_thrust(theta_col, vHub_ca(eU), vHub_ca(eW), ge_factor);
616 calc_coning_angle(theta_col);
618 calc_flapping_angles(theta_col, avFus_ca);
620 calc_drag_and_side_forces(theta_col);
622 calc_torque(theta_col);
624 // Fixme: only valid for a 'decent' rotor
625 theta_downwash = atan2( -in.AeroUVW(eU), v_induced - in.AeroUVW(eW));
626 phi_downwash = atan2( in.AeroUVW(eV), v_induced - in.AeroUVW(eW));
628 vFn = body_forces(A_IC, B_IC);
629 vMn = Transform() * body_moments(A_IC, B_IC);
633 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
635 void FGRotor::CalcStatePart2(double PowerAvailable)
639 double ExcessTorque = PowerAvailable / Omega;
640 double deltaOmega = ExcessTorque / PolarMoment * in.TotalDeltaT;
641 RPM += deltaOmega/(2.0*M_PI) * 60.0;
643 RPM = Constrain(MinimalRPM, RPM, MaximalRPM); // trim again
646 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
648 // Simulation of a free-wheeling-unit (FWU). Might need improvements.
650 void FGRotor::calc_freewheel_state(double p_source, double p_load) {
652 // engine is off/detached, release.
654 FreeWheelTransmission = 0.0;
658 // engine is driving the rotor, engage.
659 if (p_source >= p_load) {
660 FreeWheelTransmission = 1.0;
664 // releases if engine is detached, but stays calm if
665 // the load changes due to rotor dynamics.
666 if (p_source > 0.0 && p_load/(p_source+0.1) > FreeWheelThresh ) {
667 FreeWheelTransmission = 0.0;
674 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
676 double FGRotor::Calculate(double EnginePower)
683 PowerRequired = Torque * Omega + BrakeCtrlNorm * MaxBrakePower;
685 if (FreeWheelPresent) {
686 calc_freewheel_state(EnginePower * ClutchCtrlNorm, PowerRequired);
687 FWmult = FreeWheelLag.execute(FreeWheelTransmission);
690 DeltaPower = EnginePower * ClutchCtrlNorm * FWmult - PowerRequired;
692 CalcStatePart2(DeltaPower);
698 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
701 bool FGRotor::BindModel(void)
703 string property_name, base_property_name;
704 base_property_name = CreateIndexedPropertyName("propulsion/engine", EngineNum);
706 property_name = base_property_name + "/rotor-rpm";
707 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetRPM );
709 property_name = base_property_name + "/engine-rpm";
710 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetEngineRPM );
712 property_name = base_property_name + "/rotor-thrust-lbs"; // might be redundant - check!
713 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetThrust );
715 property_name = base_property_name + "/a0-rad";
716 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA0 );
718 property_name = base_property_name + "/a1-rad";
719 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA1 );
721 property_name = base_property_name + "/b1-rad";
722 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetB1 );
724 property_name = base_property_name + "/inflow-ratio";
725 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLambda );
727 property_name = base_property_name + "/advance-ratio";
728 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetMu );
730 property_name = base_property_name + "/induced-inflow-ratio";
731 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetNu );
733 property_name = base_property_name + "/vi-fps";
734 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetVi );
736 property_name = base_property_name + "/thrust-coefficient";
737 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCT );
739 property_name = base_property_name + "/torque-lbsft";
740 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetTorque );
742 property_name = base_property_name + "/theta-downwash-rad";
743 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetThetaDW );
745 property_name = base_property_name + "/phi-downwash-rad";
746 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetPhiDW );
748 switch (ControlMap) {
750 property_name = base_property_name + "/antitorque-ctrl-rad";
751 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
754 property_name = base_property_name + "/tail-collective-ctrl-rad";
755 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
756 property_name = base_property_name + "/lateral-ctrl-rad";
757 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
758 property_name = base_property_name + "/longitudinal-ctrl-rad";
759 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
761 default: // eMainCtrl
762 property_name = base_property_name + "/collective-ctrl-rad";
763 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetCollectiveCtrl, &FGRotor::SetCollectiveCtrl);
764 property_name = base_property_name + "/lateral-ctrl-rad";
765 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLateralCtrl, &FGRotor::SetLateralCtrl);
766 property_name = base_property_name + "/longitudinal-ctrl-rad";
767 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetLongitudinalCtrl, &FGRotor::SetLongitudinalCtrl);
770 property_name = base_property_name + "/brake-ctrl-norm";
771 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetBrakeCtrl, &FGRotor::SetBrakeCtrl);
772 property_name = base_property_name + "/free-wheel-transmission";
773 PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetFreeWheelTransmission);
776 if (RPMdefinition == -1) {
777 property_name = base_property_name + "/x-rpm-dict";
778 ExtRPMsource = PropertyManager->GetNode(property_name, true);
779 } else if (RPMdefinition >= 0 && RPMdefinition != EngineNum) {
780 string ipn = CreateIndexedPropertyName("propulsion/engine", RPMdefinition);
781 property_name = ipn + "/engine-rpm";
782 ExtRPMsource = PropertyManager->GetNode(property_name, false);
783 if (! ExtRPMsource) {
784 cerr << "# Warning: Engine number " << EngineNum << "." << endl;
785 cerr << "# No 'engine-rpm' property found for engine " << RPMdefinition << "." << endl;
786 cerr << "# Please check order of engine definitons." << endl;
789 cerr << "# Engine number " << EngineNum;
790 cerr << ", given ExternalRPM value '" << RPMdefinition << "' unhandled." << endl;
797 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
799 string FGRotor::GetThrusterLabels(int id, string delimeter)
804 buf << Name << " RPM (engine " << id << ")";
810 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
812 string FGRotor::GetThrusterValues(int id, string delimeter)
823 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
824 // The bitmasked value choices are as follows:
825 // unset: In this case (the default) JSBSim would only print
826 // out the normally expected messages, essentially echoing
827 // the config files as they are read. If the environment
828 // variable is not set, debug_lvl is set to 1 internally
829 // 0: This requests JSBSim not to output any messages
831 // 1: This value explicity requests the normal JSBSim
833 // 2: This value asks for a message to be printed out when
834 // a class is instantiated
835 // 4: When this value is set, a message is displayed when a
836 // FGModel object executes its Run() method
837 // 8: When this value is set, various runtime state variables
838 // are printed out periodically
839 // 16: When set various parameters are sanity checked and
840 // a message is printed out when they go out of bounds
842 void FGRotor::Debug(int from)
844 string ControlMapName;
846 if (debug_lvl <= 0) return;
848 if (debug_lvl & 1) { // Standard console startup message output
849 if (from == 0) { // Constructor
850 cout << "\n Rotor Name: " << Name << endl;
851 cout << " Diameter = " << 2.0 * Radius << " ft." << endl;
852 cout << " Number of Blades = " << BladeNum << endl;
853 cout << " Gear Ratio = " << GearRatio << endl;
854 cout << " Sense = " << Sense << endl;
855 cout << " Nominal RPM = " << NominalRPM << endl;
856 cout << " Minimal RPM = " << MinimalRPM << endl;
857 cout << " Maximal RPM = " << MaximalRPM << endl;
860 if (RPMdefinition == -1) {
861 cout << " RPM is controlled externally" << endl;
863 cout << " RPM source set to engine " << RPMdefinition << endl;
867 cout << " Blade Chord = " << BladeChord << endl;
868 cout << " Lift Curve Slope = " << LiftCurveSlope << endl;
869 cout << " Blade Twist = " << BladeTwist << endl;
870 cout << " Hinge Offset = " << HingeOffset << endl;
871 cout << " Blade Flapping Moment = " << BladeFlappingMoment << endl;
872 cout << " Blade Mass Moment = " << BladeMassMoment << endl;
873 cout << " Polar Moment = " << PolarMoment << endl;
874 cout << " Inflow Lag = " << InflowLag << endl;
875 cout << " Tip Loss = " << TipLossB << endl;
876 cout << " Lock Number = " << LockNumberByRho * 0.002356 << " (SL)" << endl;
877 cout << " Solidity = " << Solidity << endl;
878 cout << " Max Brake Power = " << MaxBrakePower/hptoftlbssec << " HP" << endl;
880 switch (ControlMap) {
881 case eTailCtrl: ControlMapName = "Tail Rotor"; break;
882 case eTandemCtrl: ControlMapName = "Tandem Rotor"; break;
883 default: ControlMapName = "Main Rotor";
885 cout << " Control Mapping = " << ControlMapName << endl;
887 if (FreeWheelPresent) {
888 cout << " Free Wheel Threshold = " << FreeWheelThresh << endl;
890 cout << " No FWU present" << endl;
895 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
896 if (from == 0) cout << "Instantiated: FGRotor" << endl;
897 if (from == 1) cout << "Destroyed: FGRotor" << endl;
899 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
901 if (debug_lvl & 8 ) { // Runtime state variables
903 if (debug_lvl & 16) { // Sanity checking
905 if (debug_lvl & 64) {
906 if (from == 0) { // Constructor
907 cout << IdSrc << endl;
908 cout << IdHdr << endl;
915 } // namespace JSBSim