Branch data Line data Source code
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 (jsb@hal-pc.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 : :
35 : : %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36 : : INCLUDES
37 : : %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
38 : :
39 : : #include <iostream>
40 : : #include <sstream>
41 : :
42 : : #include "FGRotor.h"
43 : :
44 : : #include "models/FGPropagate.h"
45 : : #include "models/FGAtmosphere.h"
46 : : #include "models/FGAuxiliary.h"
47 : : #include "models/FGMassBalance.h"
48 : :
49 : : #include "input_output/FGXMLElement.h"
50 : :
51 : : #include "math/FGRungeKutta.h"
52 : :
53 : : using namespace std;
54 : :
55 : : namespace JSBSim {
56 : :
57 : : static const char *IdSrc = "$Id: FGRotor.cpp,v 1.9 2010/06/05 12:12:34 jberndt Exp $";
58 : : static const char *IdHdr = ID_ROTOR;
59 : :
60 : : /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61 : : MISC
62 : : %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
63 : :
64 : : static int dump_req; // debug schwafel flag
65 : :
66 : 0 : static inline double sqr(double x) { return x*x; }
67 : :
68 : : /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
69 : : CLASS IMPLEMENTATION
70 : : %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
71 : :
72 : : // starting with 'inner' rotor, FGRotor constructor is further down
73 : :
74 [ # # ][ # # ]: 0 : FGRotor::rotor::~rotor() { }
[ # # ]
75 : :
76 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77 : :
78 : : // hmm, not a real alternative to a pretty long initializer list
79 : :
80 : 0 : void FGRotor::rotor::zero() {
81 : : FGColumnVector3 zero_vec(0.0, 0.0, 0.0);
82 : :
83 : 0 : flags = 0;
84 : 0 : parent = NULL ;
85 : :
86 : 0 : reports = 0;
87 : :
88 : : // configuration
89 : 0 : Radius = 0.0 ;
90 : 0 : BladeNum = 0 ;
91 : 0 : RelDistance_xhub = 0.0 ;
92 : 0 : RelShift_yhub = 0.0 ;
93 : 0 : RelHeight_zhub = 0.0 ;
94 : 0 : NominalRPM = 0.0 ;
95 : 0 : MinRPM = 0.0 ;
96 : 0 : BladeChord = 0.0 ;
97 : 0 : LiftCurveSlope = 0.0 ;
98 : 0 : BladeFlappingMoment = 0.0 ;
99 : 0 : BladeTwist = 0.0 ;
100 : 0 : BladeMassMoment = 0.0 ;
101 : 0 : TipLossB = 0.0 ;
102 : 0 : PolarMoment = 0.0 ;
103 : 0 : InflowLag = 0.0 ;
104 : 0 : ShaftTilt = 0.0 ;
105 : 0 : HingeOffset = 0.0 ;
106 : 0 : HingeOffset_hover = 0.0 ;
107 : 0 : CantAngleD3 = 0.0 ;
108 : :
109 : 0 : theta_shaft = 0.0 ;
110 : 0 : phi_shaft = 0.0 ;
111 : :
112 : : // derived parameters
113 : 0 : LockNumberByRho = 0.0 ;
114 : 0 : solidity = 0.0 ;
115 : 0 : RpmRatio = 0.0 ;
116 : :
117 [ # # ]: 0 : for (int i=0; i<5; i++) R[i] = 0.0;
118 [ # # ]: 0 : for (int i=0; i<6; i++) B[i] = 0.0;
119 : :
120 : 0 : BodyToShaft.InitMatrix();
121 : 0 : ShaftToBody.InitMatrix();
122 : :
123 : : // dynamic values
124 : 0 : ActualRPM = 0.0 ;
125 : 0 : Omega = 0.0 ;
126 : 0 : beta_orient = 0.0 ;
127 : 0 : a0 = 0.0 ;
128 : 0 : a_1 = b_1 = a_dw = 0.0 ;
129 : 0 : a1s = b1s = 0.0 ;
130 : 0 : H_drag = J_side = 0.0 ;
131 : :
132 : 0 : Torque = 0.0 ;
133 : 0 : Thrust = 0.0 ;
134 : 0 : Ct = 0.0 ;
135 : 0 : lambda = 0.0 ;
136 : 0 : mu = 0.0 ;
137 : 0 : nu = 0.0 ;
138 : 0 : v_induced = 0.0 ;
139 : :
140 : 0 : force = zero_vec ;
141 : 0 : moment = zero_vec ;
142 : :
143 : 0 : }
144 : :
145 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
146 : :
147 : : // 5in1: value-fetch-convert-default-return function
148 : :
149 : : double FGRotor::rotor::cnf_elem( const string& ename, double default_val,
150 : 0 : const string& unit, bool tell)
151 : : {
152 : :
153 : 0 : Element *e = NULL;
154 : 0 : double val=default_val;
155 : :
156 : 0 : std::string pname = "*No parent element*";
157 : :
158 [ # # ]: 0 : if (parent) {
159 : 0 : e = parent->FindElement(ename);
160 : 0 : pname = parent->GetName();
161 : : }
162 : :
163 [ # # ]: 0 : if (e) {
164 [ # # ]: 0 : if (unit.empty()) {
165 : : // val = e->FindElementValueAsNumber(ename);
166 : : // yields to: Attempting to get single data value from multiple lines
167 : 0 : val = parent->FindElementValueAsNumber(ename);
168 : : } else {
169 : : // val = e->FindElementValueAsNumberConvertTo(ename,unit);
170 : : // yields to: Attempting to get non-existent element diameter + crash, why ?
171 : 0 : val = parent->FindElementValueAsNumberConvertTo(ename,unit);
172 : : }
173 : : } else {
174 [ # # ]: 0 : if (tell) {
175 : 0 : cerr << pname << ": missing element '" << ename <<"' using estimated value: " << default_val << endl;
176 : : }
177 : : }
178 : :
179 : 0 : return val;
180 : : }
181 : :
182 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
183 : :
184 : 0 : double FGRotor::rotor::cnf_elem(const string& ename, double default_val, bool tell)
185 : : {
186 : 0 : return cnf_elem(ename, default_val, "", tell);
187 : : }
188 : :
189 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
190 : :
191 : : // 1. read configuration and try to fill holes, ymmv
192 : : // 2. calculate derived parameters and transforms
193 : 0 : void FGRotor::rotor::configure(int f, const rotor *xmain)
194 : : {
195 : :
196 : : double estimate;
197 : 0 : const bool yell = true;
198 : 0 : const bool silent = false;
199 : :
200 : 0 : flags = f;
201 : :
202 [ # # ]: 0 : estimate = (xmain) ? 2.0 * xmain->Radius * 0.2 : 42.0;
203 : 0 : Radius = 0.5 * cnf_elem("diameter", estimate, "FT", yell);
204 : :
205 [ # # ]: 0 : estimate = (xmain) ? xmain->BladeNum : 2.0;
206 : 0 : estimate = Constrain(1.0,estimate,4.0);
207 : 0 : BladeNum = (int) cnf_elem("numblades", estimate, yell);
208 : :
209 [ # # ]: 0 : estimate = (xmain) ? - xmain->Radius * 1.05 - Radius : - 0.025 * Radius ;
210 : 0 : RelDistance_xhub = cnf_elem("xhub", estimate, "FT", yell);
211 : :
212 : 0 : RelShift_yhub = cnf_elem("yhub", 0.0, "FT", silent);
213 : :
214 : 0 : estimate = - 0.1 * Radius - 4.0;
215 : 0 : RelHeight_zhub = cnf_elem("zhub", estimate, "FT", yell);
216 : :
217 : : // make sure that v_tip (omega*r) is below 0.7mach ~ 750ft/s
218 : 0 : estimate = (750.0/Radius)/(2.0*M_PI) * 60.0; // 7160/Radius
219 : 0 : NominalRPM = cnf_elem("nominalrpm", estimate, yell);
220 : :
221 : 0 : MinRPM = cnf_elem("minrpm", 1.0, silent);
222 : 0 : MinRPM = Constrain(1.0, MinRPM, NominalRPM-1.0);
223 : :
224 [ # # ]: 0 : estimate = (xmain) ? 0.12 : 0.07; // guess solidity
225 : 0 : estimate = estimate * M_PI*Radius/BladeNum;
226 : 0 : BladeChord = cnf_elem("chord", estimate, "FT", yell);
227 : :
228 : 0 : LiftCurveSlope = cnf_elem("liftcurveslope", 6.0, yell); // "1/RAD"
229 : :
230 : 0 : estimate = sqr(BladeChord) * sqr(Radius) * 0.57;
231 : 0 : BladeFlappingMoment = cnf_elem("flappingmoment", estimate, "SLUG*FT2", yell);
232 : 0 : BladeFlappingMoment = Constrain(0.1, BladeFlappingMoment, 1e9);
233 : :
234 : 0 : BladeTwist = cnf_elem("twist", -0.17, "RAD", yell);
235 : :
236 : 0 : estimate = sqr(BladeChord) * BladeChord * 15.66; // might be really wrong!
237 : 0 : BladeMassMoment = cnf_elem("massmoment", estimate, yell); // slug-ft
238 : 0 : BladeMassMoment = Constrain(0.1, BladeMassMoment, 1e9);
239 : :
240 : 0 : TipLossB = cnf_elem("tiplossfactor", 0.98, silent);
241 : :
242 : 0 : estimate = 1.1 * BladeFlappingMoment * BladeNum;
243 : 0 : PolarMoment = cnf_elem("polarmoment", estimate, "SLUG*FT2", silent);
244 : 0 : PolarMoment = Constrain(0.1, PolarMoment, 1e9);
245 : :
246 : 0 : InflowLag = cnf_elem("inflowlag", 0.2, silent); // fixme, depends on size
247 : :
248 [ # # ]: 0 : estimate = (xmain) ? 0.0 : -0.06;
249 : 0 : ShaftTilt = cnf_elem("shafttilt", estimate, "RAD", silent);
250 : :
251 : : // ignore differences between teeter/hingeless/fully-articulated constructions
252 : 0 : estimate = 0.05 * Radius;
253 [ # # ]: 0 : HingeOffset = cnf_elem("hingeoffset", estimate, "FT", (xmain) ? silent : yell);
254 : :
255 : 0 : CantAngleD3 = cnf_elem("cantangle", 0.0, "RAD", silent);
256 : :
257 : : // derived parameters
258 : :
259 : : // precalc often used powers
260 : 0 : 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];
261 : 0 : 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]; B[5]=B[4]*B[1];
262 : :
263 : 0 : LockNumberByRho = LiftCurveSlope * BladeChord * R[4] / BladeFlappingMoment;
264 : 0 : solidity = BladeNum * BladeChord / (M_PI * Radius);
265 : :
266 : : // use simple orientations at the moment
267 [ # # ]: 0 : if (flags & eTail) { // axis parallel to Y_body
268 : 0 : theta_shaft = 0.0; // no tilt
269 : 0 : phi_shaft = 0.5*M_PI;
270 : :
271 : : // opposite direction if main rotor is spinning CW
272 [ # # ][ # # ]: 0 : if (xmain && (xmain->flags & eRotCW) ) {
273 : 0 : phi_shaft = -phi_shaft;
274 : : }
275 : : } else { // more or less upright
276 : 0 : theta_shaft = ShaftTilt;
277 : 0 : phi_shaft = 0.0;
278 : : }
279 : :
280 : : // setup Shaft-Body transforms, see /SH79/ eqn(17,18)
281 : 0 : double st = sin(theta_shaft);
282 : 0 : double ct = cos(theta_shaft);
283 : 0 : double sp = sin(phi_shaft);
284 : 0 : double cp = cos(phi_shaft);
285 : :
286 : : ShaftToBody.InitMatrix( ct, st*sp, st*cp,
287 : : 0.0, cp, -sp,
288 : 0 : -st, ct*sp, ct*cp );
289 : :
290 : 0 : BodyToShaft = ShaftToBody.Inverse();
291 : :
292 : : // misc defaults
293 : 0 : nu = 0.001; // help the flow solver by providing some moving molecules
294 : :
295 : : return;
296 : : }
297 : :
298 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
299 : :
300 : : // calculate control-axes components of total airspeed at the hub.
301 : : // sets rotor orientation angle (beta) as side effect. /SH79/ eqn(19-22)
302 : :
303 : : FGColumnVector3 FGRotor::rotor::hub_vel_body2ca( const FGColumnVector3 &uvw,
304 : : const FGColumnVector3 &pqr,
305 : 0 : double a_ic, double b_ic)
306 : : {
307 : :
308 : 0 : FGColumnVector3 v_r, v_shaft, v_w;
309 : 0 : FGColumnVector3 pos(RelDistance_xhub,0.0,RelHeight_zhub);
310 : :
311 : : v_r = uvw + pqr*pos;
312 : 0 : v_shaft = BodyToShaft * v_r;
313 : :
314 : 0 : beta_orient = atan2(v_shaft(eV),v_shaft(eU));
315 : :
316 : 0 : v_w(eU) = v_shaft(eU)*cos(beta_orient) + v_shaft(eV)*sin(beta_orient);
317 : 0 : v_w(eV) = 0.0;
318 : 0 : v_w(eW) = v_shaft(eW) - b_ic*v_shaft(eU) - a_ic*v_shaft(eV);
319 : :
320 : 0 : return v_w;
321 : : }
322 : :
323 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
324 : :
325 : : // express fuselage angular velocity in control axes /SH79/ eqn(30,31)
326 : :
327 : 0 : FGColumnVector3 FGRotor::rotor::fus_angvel_body2ca( const FGColumnVector3 &pqr)
328 : : {
329 : 0 : FGColumnVector3 av_s_fus, av_w_fus;
330 : :
331 : 0 : av_s_fus = BodyToShaft * pqr;
332 : :
333 : 0 : av_w_fus(eP)= av_s_fus(eP)*cos(beta_orient) + av_s_fus(eQ)*sin(beta_orient);
334 : 0 : av_w_fus(eQ)= - av_s_fus(eP)*sin(beta_orient) + av_s_fus(eQ)*cos(beta_orient);
335 : 0 : av_w_fus(eR)= av_s_fus(eR);
336 : :
337 : 0 : return av_w_fus;
338 : : }
339 : :
340 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
341 : :
342 : : // problem function passed to rk solver
343 : :
344 : 0 : double FGRotor::rotor::dnuFunction::pFunc(double x, double nu) {
345 : : double d_nu;
346 : : d_nu = k_sat * (ct_lambda * (k_wor - nu) + k_theta) /
347 : 0 : (2.0 * sqrt( mu2 + sqr(k_wor - nu)) );
348 : 0 : d_nu = d_nu * k_flowscale - nu;
349 : 0 : return d_nu;
350 : : };
351 : :
352 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
353 : :
354 : : // merge params to keep the equation short
355 : 0 : void FGRotor::rotor::dnuFunction::update_params(rotor *r, double ct_t01, double fs, double w) {
356 : 0 : k_sat = 0.5* r->solidity * r->LiftCurveSlope;
357 : 0 : ct_lambda = 1.0/2.0*r->B[2] + 1.0/4.0 * r->mu*r->mu;
358 : 0 : k_wor = w/(r->Omega*r->Radius);
359 : 0 : k_theta = ct_t01;
360 : 0 : mu2 = r->mu * r->mu;
361 : 0 : k_flowscale = fs;
362 : 0 : };
363 : :
364 : :
365 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
366 : :
367 : : // Calculate rotor thrust and inflow-ratio (lambda), this is achieved by
368 : : // approximating a solution for the differential equation:
369 : : //
370 : : // dnu/dt = 1/tau ( Ct / (2*sqrt(mu^2+lambda^2)) - nu ) , /SH79/ eqn(26).
371 : : //
372 : : // Propper calculation of the inflow-ratio (lambda) is vital for the
373 : : // following calculations. Simple implementations (i.e. Newton-Raphson w/o
374 : : // checking) tend to oscillate or overshoot in the low speed region,
375 : : // therefore a more expensive solver is used.
376 : : //
377 : : // The flow_scale parameter is used to approximate a reduction of inflow
378 : : // if the helicopter is close to the ground, yielding to higher thrust,
379 : : // see /TA77/ eqn(10a). Doing the ground effect calculation here seems
380 : : // more favorable then to code it in the fdm_config.
381 : :
382 : : void FGRotor::rotor::calc_flow_and_thrust(
383 : : double dt, double rho, double theta_0,
384 : 0 : double Uw, double Ww, double flow_scale)
385 : : {
386 : :
387 : 0 : double ct_over_sigma = 0.0;
388 : : double ct_l, ct_t0, ct_t1;
389 : 0 : double nu_ret = 0.0;
390 : :
391 : 0 : mu = Uw/(Omega*Radius); // /SH79/ eqn(24)
392 : :
393 : 0 : ct_t0 = (1.0/3.0*B[3] + 1.0/2.0 * TipLossB*mu*mu - 4.0/(9.0*M_PI) * mu*mu*mu )*theta_0;
394 : 0 : ct_t1 = (1.0/4.0*B[4] + 1.0/4.0 * B[2]*mu*mu)*BladeTwist;
395 : :
396 : : // merge params together
397 : 0 : flowEquation.update_params(this, ct_t0+ct_t1, flow_scale, Ww);
398 : :
399 : 0 : nu_ret = rk.evolve(nu, &flowEquation);
400 : :
401 [ # # ]: 0 : if (rk.getStatus() != FGRungeKutta::eNoError) { // never observed so far
402 : 0 : cerr << "# IEHHHH [" << flags << "]: Solver Error - resetting!" << endl;
403 : 0 : rk.clearStatus();
404 : 0 : nu_ret = nu; // use old value and keep fingers crossed.
405 : : }
406 : :
407 : : // keep an eye on the solver, but be quiet after a hundred messages
408 [ # # ][ # # ]: 0 : if (reports < 100 && rk.getIterations()>6) {
[ # # ]
409 : : cerr << "# LOOK [" << flags << "]: Solver took "
410 : 0 : << rk.getIterations() << " rounds." << endl;
411 : 0 : reports++;
412 [ # # ]: 0 : if (reports==100) {
413 : 0 : cerr << "# stopped babbling after 100 notifications." << endl;
414 : : }
415 : : }
416 : :
417 : : // now from nu to lambda, Ct, and Thrust
418 : :
419 : 0 : nu = nu_ret;
420 : 0 : lambda = Ww/(Omega*Radius) - nu; // /SH79/ eqn(25)
421 : :
422 : 0 : ct_l = (1.0/2.0*B[2] + 1.0/4.0 * mu*mu)*lambda;
423 : 0 : ct_over_sigma = (LiftCurveSlope/2.0)*(ct_l + ct_t0 + ct_t1); // /SH79/ eqn(27)
424 : :
425 : 0 : Thrust = BladeNum*BladeChord*Radius*rho*sqr(Omega*Radius) * ct_over_sigma;
426 : :
427 : 0 : Ct = ct_over_sigma * solidity;
428 : 0 : v_induced = nu * (Omega*Radius);
429 : :
430 [ # # ][ # # ]: 0 : if (dump_req && (flags & eMain) ) {
431 : 0 : printf("# mu %f : nu %f lambda %f vi %f\n", mu, nu, lambda, v_induced);
432 : : }
433 : :
434 : 0 : }
435 : :
436 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
437 : :
438 : : // this is the most arcane part in the calculation chain the
439 : : // constants used should be reverted to more general parameters.
440 : : // otoh: it works also for smaller rotors, sigh!
441 : : // See /SH79/ eqn(36), and /BA41/ for a richer set of equations.
442 : :
443 : 0 : void FGRotor::rotor::calc_torque(double rho, double theta_0)
444 : : {
445 : : double Qa0;
446 : : double cq_s_m[5], cq_over_sigma;
447 : : double l,m,t075; // shortcuts
448 : :
449 : 0 : t075 = theta_0 + 0.75 * BladeTwist;
450 : :
451 : 0 : m = mu;
452 : 0 : l = lambda;
453 : :
454 : 0 : cq_s_m[0] = 0.00109 - 0.0036*l - 0.0027*t075 - 1.10*sqr(l) - 0.545*l*t075 + 0.122*sqr(t075);
455 : 0 : cq_s_m[2] = ( 0.00109 - 0.0027*t075 - 3.13*sqr(l) - 6.35*l*t075 - 1.93*sqr(t075) ) * sqr(m);
456 : 0 : cq_s_m[3] = - 0.133*l*t075 * sqr(m)*m;
457 : 0 : cq_s_m[4] = ( - 0.976*sqr(l) - 6.38*l*t075 - 5.26*sqr(t075) ) * sqr(m)*sqr(m);
458 : :
459 : 0 : cq_over_sigma = cq_s_m[0] + cq_s_m[2] + cq_s_m[3] + cq_s_m[4];
460 : : // guess an a (LiftCurveSlope) is included in eqn above, so check if there is a large
461 : : // influence when a_'other-model'/ a_'ch54' diverts from 1.0.
462 : :
463 : 0 : Qa0 = BladeNum * BladeChord * R[2] * rho * sqr(Omega*Radius);
464 : :
465 : : // TODO: figure out how to handle negative cq_over_sigma/torque
466 : :
467 : 0 : Torque = Qa0 * cq_over_sigma;
468 : :
469 : : return;
470 : : }
471 : :
472 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
473 : :
474 : : // The coning angle doesn't apply for teetering rotors, but calculating
475 : : // doesn't hurt. /SH79/ eqn(29)
476 : :
477 : 0 : void FGRotor::rotor::calc_coning_angle(double rho, double theta_0)
478 : : {
479 : 0 : double lock_gamma = LockNumberByRho * rho;
480 : :
481 : 0 : double a0_l = (1.0/6.0 * B[3] + 0.04 * mu*mu*mu) * lambda;
482 : 0 : double a0_t0 = (1.0/8.0 * B[4] + 1.0/8.0 * B[2]*mu*mu) * theta_0;
483 : 0 : double a0_t1 = (1.0/10.0 * B[5] + 1.0/12.0 * B[3]*mu*mu) * BladeTwist;
484 : 0 : a0 = lock_gamma * ( a0_l + a0_t0 + a0_t1);
485 : : return;
486 : : }
487 : :
488 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
489 : :
490 : : // Flapping angles relative to control axes /SH79/ eqn(32)
491 : :
492 : : void FGRotor::rotor::calc_flapping_angles( double rho, double theta_0,
493 : 0 : const FGColumnVector3 &pqr_fus_w)
494 : : {
495 : 0 : double lock_gamma = LockNumberByRho * rho;
496 : :
497 : 0 : double mu2_2B2 = sqr(mu)/(2.0*B[2]);
498 : 0 : double t075 = theta_0 + 0.75 * BladeTwist; // common approximation for rectangular blades
499 : :
500 : : a_1 = 1.0/(1.0 - mu2_2B2) * (
501 : : (2.0*lambda + (8.0/3.0)*t075)*mu
502 : : + pqr_fus_w(eP)/Omega
503 : : - 16.0 * pqr_fus_w(eQ)/(B[4]*lock_gamma*Omega)
504 : 0 : );
505 : :
506 : : b_1 = 1.0/(1.0 + mu2_2B2) * (
507 : : (4.0/3.0)*mu*a0
508 : : - pqr_fus_w(eQ)/Omega
509 : : - 16.0 * pqr_fus_w(eP)/(B[4]*lock_gamma*Omega)
510 : 0 : );
511 : :
512 : : // used in force calc
513 : : a_dw = 1.0/(1.0 - mu2_2B2) * (
514 : : (2.0*lambda + (8.0/3.0)*t075)*mu
515 : : - 24.0 * pqr_fus_w(eQ)/(B[4]*lock_gamma*Omega)
516 : : * ( 1.0 - ( 0.29 * t075 / (Ct/solidity) ) )
517 : 0 : );
518 : :
519 : : return;
520 : : }
521 : :
522 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
523 : :
524 : : // /SH79/ eqn(38,39)
525 : :
526 : 0 : void FGRotor::rotor::calc_drag_and_side_forces(double rho, double theta_0)
527 : : {
528 : : double cy_over_sigma ;
529 : 0 : double t075 = theta_0 + 0.75 * BladeTwist;
530 : :
531 : 0 : H_drag = Thrust * a_dw;
532 : :
533 : : cy_over_sigma = (
534 : : 0.75*b_1*lambda - 1.5*a0*mu*lambda + 0.25*a_1*b_1*mu
535 : : - a0*a_1*sqr(mu) + (1.0/6.0)*a0*a_1
536 : : - (0.75*mu*a0 - (1.0/3.0)*b_1 - 0.5*sqr(mu)*b_1)*t075
537 : 0 : );
538 : 0 : cy_over_sigma *= LiftCurveSlope/2.0;
539 : :
540 : 0 : J_side = BladeNum * BladeChord * Radius * rho * sqr(Omega*Radius) * cy_over_sigma;
541 : :
542 : : return;
543 : : }
544 : :
545 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
546 : :
547 : : // transform rotor forces from control axes to shaft axes, and express
548 : : // in body axes /SH79/ eqn(40,41)
549 : :
550 : 0 : FGColumnVector3 FGRotor::rotor::body_forces(double a_ic, double b_ic)
551 : : {
552 : : FGColumnVector3 F_s(
553 : : - H_drag*cos(beta_orient) - J_side*sin(beta_orient) + Thrust*b_ic,
554 : : - H_drag*sin(beta_orient) + J_side*cos(beta_orient) + Thrust*a_ic,
555 : 0 : - Thrust);
556 : :
557 [ # # ][ # # ]: 0 : if (dump_req && (flags & eMain) ) {
558 : 0 : printf("# abß: % f % f % f\n", a_ic, b_ic, beta_orient );
559 : 0 : printf("# HJT: % .2f % .2f % .2f\n", H_drag, J_side, Thrust );
560 : 0 : printf("# F_s: % .2f % .2f % .2f\n", F_s(1), F_s(2), F_s(3) );
561 : 0 : FGColumnVector3 F_h;
562 : 0 : F_h = ShaftToBody * F_s;
563 : 0 : printf("# F_h: % .2f % .2f % .2f\n", F_h(1), F_h(2), F_h(3) );
564 : : }
565 : :
566 : 0 : return ShaftToBody * F_s;
567 : : }
568 : :
569 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
570 : :
571 : : // rotational sense is handled here
572 : : // still a todo: how to get propper values for 'BladeMassMoment'
573 : : // here might be a good place to tweak hovering stability, check /AM50/.
574 : :
575 : 0 : FGColumnVector3 FGRotor::rotor::body_moments(FGColumnVector3 F_h, double a_ic, double b_ic)
576 : : {
577 : 0 : FGColumnVector3 M_s, M_hub, M_h;
578 : :
579 : 0 : FGColumnVector3 h_pos(RelDistance_xhub, 0.0, RelHeight_zhub);
580 : :
581 : : // vermutlich ein biege moment, bzw.widerstands moment ~ d^3
582 : 0 : double M_w_tilde = 0.0 ;
583 : 0 : double mf = 0.0 ;
584 : :
585 : 0 : M_w_tilde = BladeMassMoment;
586 : :
587 : : // cyclic flapping relative to shaft axes /SH79/ eqn(43)
588 : 0 : a1s = a_1*cos(beta_orient) + b_1*sin(beta_orient) - b_ic;
589 : 0 : b1s = b_1*cos(beta_orient) - a_1*sin(beta_orient) + a_ic;
590 : :
591 : : // mind this: no HingeOffset, no additional pitch/roll moments
592 : 0 : mf = 0.5 * (HingeOffset+HingeOffset_hover) * BladeNum * Omega*Omega * M_w_tilde;
593 : 0 : M_s(eL) = mf*b1s;
594 : 0 : M_s(eM) = mf*a1s;
595 : 0 : M_s(eN) = Torque;
596 : :
597 [ # # ]: 0 : if (flags & eRotCW) {
598 : 0 : M_s(eN) = -M_s(eN);
599 : : }
600 : :
601 [ # # ]: 0 : if (flags & eCoaxial) {
602 : 0 : M_s(eN) = 0.0;
603 : : }
604 : :
605 : 0 : M_hub = ShaftToBody * M_s;
606 : :
607 : 0 : M_h = M_hub + (h_pos * F_h);
608 : :
609 : 0 : return M_h;
610 : : }
611 : :
612 : :
613 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
614 : :
615 : : // Constructor
616 : :
617 : 0 : FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
618 : 0 : : FGThruster(exec, rotor_element, num)
619 : : {
620 : :
621 : 0 : FGColumnVector3 location, orientation;
622 : : Element *thruster_element;
623 : :
624 : 0 : PropertyManager = fdmex->GetPropertyManager();
625 : 0 : dt = fdmex->GetDeltaT();
626 : :
627 : : /* apply defaults */
628 : :
629 : 0 : rho = 0.002356; // just a sane value
630 : :
631 : 0 : RPM = 0.0;
632 : 0 : Sense = 1.0;
633 : 0 : tailRotorPresent = false;
634 : :
635 : 0 : effective_tail_col = 0.001; // just a sane value
636 : :
637 : 0 : prop_inflow_ratio_lambda = 0.0;
638 : 0 : prop_advance_ratio_mu = 0.0;
639 : 0 : prop_inflow_ratio_induced_nu = 0.0;
640 : 0 : prop_mr_torque = 0.0;
641 : 0 : prop_thrust_coefficient = 0.0;
642 : 0 : prop_coning_angle = 0.0;
643 : :
644 : 0 : prop_theta_downwash = prop_phi_downwash = 0.0;
645 : :
646 : 0 : hover_threshold = 0.0;
647 : 0 : hover_scale = 0.0;
648 : :
649 : 0 : mr.zero();
650 : 0 : tr.zero();
651 : :
652 : : // debug stuff
653 : 0 : prop_DumpFlag = 0;
654 : :
655 : : /* configure */
656 : :
657 : 0 : Type = ttRotor;
658 : 0 : SetTransformType(FGForce::tCustom);
659 : :
660 : : // get data from parent and 'mount' the rotor system
661 : :
662 : 0 : thruster_element = rotor_element->GetParent()->FindElement("sense");
663 [ # # # # ]: 0 : if (thruster_element) {
664 [ # # ][ # # ]: 0 : Sense = thruster_element->GetDataAsNumber() >= 0.0 ? 1.0: -1.0;
665 : : }
666 : :
667 : 0 : thruster_element = rotor_element->GetParent()->FindElement("location");
668 [ # # ][ # # ]: 0 : if (thruster_element) location = thruster_element->FindElementTripletConvertTo("IN");
669 : 0 : else cerr << "No thruster location found." << endl;
670 : :
671 : 0 : thruster_element = rotor_element->GetParent()->FindElement("orient");
672 [ # # ][ # # ]: 0 : if (thruster_element) orientation = thruster_element->FindElementTripletConvertTo("RAD");
673 : 0 : else cerr << "No thruster orientation found." << endl;
674 : :
675 : 0 : SetLocation(location);
676 : 0 : SetAnglesToBody(orientation);
677 : :
678 : : // get main rotor parameters
679 : 0 : mr.parent = rotor_element;
680 : :
681 : 0 : int flags = eMain;
682 : :
683 : 0 : string a_val="";
684 : 0 : a_val = rotor_element->GetAttributeValue("variant");
685 [ # # # # ]: 0 : if ( a_val == "coaxial" ) {
686 : 0 : flags += eCoaxial;
687 : 0 : cerr << "# found 'coaxial' variant" << endl;
688 : : }
689 : :
690 [ # # ][ # # ]: 0 : if (Sense<0.0) {
691 : 0 : flags += eRotCW;
692 : : }
693 : :
694 : 0 : mr.configure(flags);
695 : :
696 : 0 : mr.rk.init(0,dt,6);
697 : :
698 : : // get tail rotor parameters
699 : 0 : tr.parent=rotor_element->FindElement("tailrotor");
700 [ # # # # ]: 0 : if (tr.parent) {
701 : 0 : tailRotorPresent = true;
702 : : } else {
703 : 0 : tailRotorPresent = false;
704 : 0 : cerr << "# No tailrotor found, assuming a single rotor." << endl;
705 : : }
706 : :
707 [ # # ][ # # ]: 0 : if (tailRotorPresent) {
708 : 0 : int flags = eTail;
709 [ # # ][ # # ]: 0 : if (Sense<0.0) {
710 : 0 : flags += eRotCW;
711 : : }
712 : 0 : tr.configure(flags, &mr);
713 : 0 : tr.rk.init(0,dt,6);
714 : 0 : tr.RpmRatio = tr.NominalRPM/mr.NominalRPM; // 'connect'
715 : : }
716 : :
717 : : /* remaining parameters */
718 : :
719 : : // ground effect
720 : 0 : double c_ground_effect = 0.0; // uh1 ~ 0.28 , larger values increase the effect
721 : 0 : ground_effect_exp = 0.0;
722 : 0 : ground_effect_shift = 0.0;
723 : :
724 [ # # ][ # # ]: 0 : if (rotor_element->FindElement("cgroundeffect"))
725 : 0 : c_ground_effect = rotor_element->FindElementValueAsNumber("cgroundeffect");
726 : :
727 [ # # ][ # # ]: 0 : if (rotor_element->FindElement("groundeffectshift"))
728 : 0 : ground_effect_shift = rotor_element->FindElementValueAsNumberConvertTo("groundeffectshift","FT");
729 : :
730 : : // prepare calculations, see /TA77/
731 [ # # ][ # # ]: 0 : if (c_ground_effect > 1e-5) {
732 : 0 : ground_effect_exp = 1.0 / ( 2.0*mr.Radius * c_ground_effect );
733 : : } else {
734 : 0 : ground_effect_exp = 0.0; // disable
735 : : }
736 : :
737 : : // smooth out jumps in hagl reported, otherwise the ground effect
738 : : // calculation would cause jumps too. 1Hz seems sufficient.
739 : 0 : damp_hagl = Filter(1.0,dt);
740 : :
741 : :
742 : : // misc, experimental
743 [ # # ][ # # ]: 0 : if (rotor_element->FindElement("hoverthreshold"))
744 : 0 : hover_threshold = rotor_element->FindElementValueAsNumberConvertTo("hoverthreshold", "FT/SEC");
745 : :
746 [ # # ][ # # ]: 0 : if (rotor_element->FindElement("hoverscale"))
747 : 0 : hover_scale = rotor_element->FindElementValueAsNumber("hoverscale");
748 : :
749 : : // enable import-export
750 : 0 : bind();
751 : :
752 : : // unused right now
753 : 0 : prop_rotorbrake->setDoubleValue(0.0);
754 : 0 : prop_freewheel_factor->setDoubleValue(1.0);
755 : :
756 : 0 : Debug(0);
757 : :
758 : 0 : } // Constructor
759 : :
760 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
761 : :
762 : 0 : FGRotor::~FGRotor()
763 : : {
764 : 0 : Debug(1);
765 [ # # ][ # # ]: 0 : }
[ # # ]
766 : :
767 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
768 : :
769 : : // mea-culpa - the connection to the engine might be wrong, but the calling
770 : : // convention appears to be 'variable' too.
771 : : // piston call:
772 : : // return Thruster->Calculate((Eng_HP * hptoftlbssec)-Thruster->GetPowerRequired());
773 : : // turbine call:
774 : : // Thrust = Thruster->Calculate(Thrust); // allow thruster to modify thrust (i.e. reversing)
775 : : //
776 : : // Here 'Calculate' takes thrust and estimates the power provided.
777 : :
778 : 0 : double FGRotor::Calculate(double PowerAvailable)
779 : : {
780 : : // controls
781 : : double A_IC; // lateral (roll) control in radians
782 : : double B_IC; // longitudinal (pitch) control in radians
783 : : double theta_col; // main rotor collective pitch in radians
784 : : double tail_col; // tail rotor collective in radians
785 : :
786 : : // state
787 : 0 : double h_agl_ft = 0.0;
788 : : double Vt ;
789 : :
790 : 0 : FGColumnVector3 UVW_h;
791 : 0 : FGColumnVector3 PQR_h;
792 : :
793 : : /* total vehicle velocity including wind effects in feet per second. */
794 : 0 : Vt = fdmex->GetAuxiliary()->GetVt();
795 : :
796 : 0 : dt = fdmex->GetDeltaT(); // might be variable ?
797 : :
798 : 0 : dump_req = prop_DumpFlag;
799 : 0 : prop_DumpFlag = 0;
800 : :
801 : : // fetch often used values
802 : 0 : rho = fdmex->GetAtmosphere()->GetDensity(); // slugs/ft^3.
803 : :
804 : 0 : UVW_h = fdmex->GetAuxiliary()->GetAeroUVW();
805 : 0 : PQR_h = fdmex->GetAuxiliary()->GetAeroPQR();
806 : :
807 : : // handle present RPM now, calc omega values.
808 : :
809 [ # # ]: 0 : if (RPM < mr.MinRPM) { // kludge, otherwise calculations go bananas
810 : 0 : RPM = mr.MinRPM;
811 : : }
812 : :
813 : 0 : mr.ActualRPM = RPM;
814 : 0 : mr.Omega = (RPM/60.0)*2.0*M_PI;
815 : :
816 [ # # ]: 0 : if (tailRotorPresent) {
817 : 0 : tr.ActualRPM = RPM*tr.RpmRatio;
818 : 0 : tr.Omega = (RPM*tr.RpmRatio/60.0)*2.0*M_PI;
819 : : }
820 : :
821 : : // read control inputs
822 : :
823 : 0 : A_IC = prop_lateral_ctrl->getDoubleValue();
824 : 0 : B_IC = prop_longitudinal_ctrl->getDoubleValue();
825 : 0 : theta_col = prop_collective_ctrl->getDoubleValue();
826 : 0 : tail_col = 0.0;
827 [ # # ]: 0 : if (tailRotorPresent) {
828 : 0 : tail_col = prop_antitorque_ctrl->getDoubleValue();
829 : : }
830 : :
831 : :
832 : 0 : FGColumnVector3 vHub_ca = mr.hub_vel_body2ca(UVW_h,PQR_h,A_IC,B_IC);
833 : 0 : FGColumnVector3 avFus_ca = mr.fus_angvel_body2ca(PQR_h);
834 : :
835 : :
836 : 0 : h_agl_ft = fdmex->GetPropagate()->GetDistanceAGL();
837 : :
838 : : double filtered_hagl;
839 : 0 : filtered_hagl = damp_hagl.execute( h_agl_ft + ground_effect_shift );
840 : :
841 : : // gnuplot> plot [-1:50] 1 - exp((-x/44)/.28)
842 : 0 : double ge_factor = 1.0;
843 [ # # ]: 0 : if (ground_effect_exp > 1e-5) {
844 : 0 : ge_factor -= exp(-filtered_hagl*ground_effect_exp);
845 : : }
846 : : // clamp
847 [ # # ]: 0 : if (ge_factor<0.5) ge_factor=0.5;
848 : :
849 [ # # ]: 0 : if (dump_req) {
850 : 0 : printf("# GE h: %.3f (%.3f) f: %f\n", filtered_hagl, h_agl_ft + ground_effect_shift, ge_factor);
851 : : }
852 : :
853 : :
854 : : // EXPERIMENTAL: modify rotor for hover, see rotor::body_moments for the consequences
855 [ # # ][ # # ]: 0 : if (hover_threshold > 1e-5 && Vt < hover_threshold) {
856 : 0 : double scale = 1.0 - Vt/hover_threshold;
857 : 0 : mr.HingeOffset_hover = scale*hover_scale*mr.Radius;
858 : : } else {
859 : 0 : mr.HingeOffset_hover = 0.0;
860 : : }
861 : :
862 : : // all set, start calculations
863 : :
864 : : /* MAIN ROTOR */
865 : :
866 : 0 : mr.calc_flow_and_thrust(dt, rho, theta_col, vHub_ca(eU), vHub_ca(eW), ge_factor);
867 : :
868 : 0 : prop_inflow_ratio_lambda = mr.lambda;
869 : 0 : prop_advance_ratio_mu = mr.mu;
870 : 0 : prop_inflow_ratio_induced_nu = mr.nu;
871 : 0 : prop_thrust_coefficient = mr.Ct;
872 : :
873 : 0 : mr.calc_coning_angle(rho, theta_col);
874 : 0 : prop_coning_angle = mr.a0;
875 : :
876 : 0 : mr.calc_torque(rho, theta_col);
877 : 0 : prop_mr_torque = mr.Torque;
878 : :
879 : 0 : mr.calc_flapping_angles(rho, theta_col, avFus_ca);
880 : 0 : mr.calc_drag_and_side_forces(rho, theta_col);
881 : :
882 : 0 : FGColumnVector3 F_h_mr = mr.body_forces(A_IC,B_IC);
883 : 0 : FGColumnVector3 M_h_mr = mr.body_moments(F_h_mr, A_IC, B_IC);
884 : :
885 : : // export downwash angles
886 : : // theta: positive for downwash moving from +z_h towards +x_h
887 : : // phi: positive for downwash moving from +z_h towards -y_h
888 : :
889 : 0 : prop_theta_downwash = atan2( - UVW_h(eU), mr.v_induced - UVW_h(eW));
890 : 0 : prop_phi_downwash = atan2( UVW_h(eV), mr.v_induced - UVW_h(eW));
891 : :
892 : 0 : mr.force(eX) = F_h_mr(1);
893 : 0 : mr.force(eY) = F_h_mr(2);
894 : 0 : mr.force(eZ) = F_h_mr(3);
895 : :
896 : 0 : mr.moment(eL) = M_h_mr(1);
897 : 0 : mr.moment(eM) = M_h_mr(2);
898 : 0 : mr.moment(eN) = M_h_mr(3);
899 : :
900 : : /* TAIL ROTOR */
901 : :
902 : : FGColumnVector3 F_h_tr(0.0, 0.0, 0.0);
903 : : FGColumnVector3 M_h_tr(0.0, 0.0, 0.0);
904 : :
905 [ # # ]: 0 : if (tailRotorPresent) {
906 : 0 : FGColumnVector3 vHub_ca_tr = tr.hub_vel_body2ca(UVW_h,PQR_h);
907 : 0 : FGColumnVector3 avFus_ca_tr = tr.fus_angvel_body2ca(PQR_h);
908 : :
909 : 0 : tr.calc_flow_and_thrust(dt, rho, tail_col, vHub_ca_tr(eU), vHub_ca_tr(eW));
910 : 0 : tr.calc_coning_angle(rho, tail_col);
911 : :
912 : : // test code for cantered tail rotor, see if it has a notable effect. /SH79/ eqn(47)
913 [ # # ]: 0 : if (fabs(tr.CantAngleD3)>1e-5) {
914 : 0 : double tan_d3 = tan(tr.CantAngleD3);
915 : : double d_t0t;
916 : 0 : double ca_dt = dt/12.0;
917 [ # # ]: 0 : for (int i = 0; i<12; i++) {
918 : 0 : d_t0t = 1/0.1*(tail_col - tr.a0 * tan_d3 - effective_tail_col);
919 : 0 : effective_tail_col += d_t0t*ca_dt;
920 : : }
921 : : } else {
922 : 0 : effective_tail_col = tail_col;
923 : : }
924 : :
925 : 0 : tr.calc_torque(rho, effective_tail_col);
926 : 0 : tr.calc_flapping_angles(rho, effective_tail_col, avFus_ca_tr);
927 : 0 : tr.calc_drag_and_side_forces(rho, effective_tail_col);
928 : :
929 : 0 : F_h_tr = tr.body_forces();
930 : 0 : M_h_tr = tr.body_moments(F_h_tr);
931 : : }
932 : :
933 : 0 : tr.force(eX) = F_h_tr(1) ;
934 : 0 : tr.force(eY) = F_h_tr(2) ;
935 : 0 : tr.force(eZ) = F_h_tr(3) ;
936 : 0 : tr.moment(eL) = M_h_tr(1) ;
937 : 0 : tr.moment(eM) = M_h_tr(2) ;
938 : 0 : tr.moment(eN) = M_h_tr(3) ;
939 : :
940 : : /*
941 : : TODO:
942 : : check negative mr.Torque conditions
943 : : freewheel factor: assure [0..1] just multiply with available power
944 : : rotorbrake: just steal from available power
945 : :
946 : : */
947 : :
948 : : // calculate new RPM, assuming a stiff connection between engine and rotor.
949 : :
950 : 0 : double engine_hp = PowerAvailable/2.24; // 'undo' force via the estimation factor used in aeromatic
951 : 0 : double engine_torque = 550.0*engine_hp/mr.Omega;
952 : 0 : double Omega_dot = (engine_torque - mr.Torque) / mr.PolarMoment;
953 : :
954 : 0 : RPM += ( Omega_dot * dt )/(2.0*M_PI) * 60.0;
955 : :
956 : : if (0 && dump_req) {
957 : : printf("# SENSE : % d % d\n", mr.flags & eRotCW ? -1 : 1, tr.flags & eRotCW ? -1 : 1);
958 : : printf("# vi : % f % f\n", mr.v_induced, tr.v_induced);
959 : : printf("# a0 a1 b1 : % f % f % f\n", mr.a0, mr.a_1, mr.b_1 );
960 : : printf("# m forces : % f % f % f\n", mr.force(eX), mr.force(eY), mr.force(eZ) );
961 : : printf("# m moments : % f % f % f\n", mr.moment(eL), mr.moment(eM), mr.moment(eN) );
962 : : printf("# t forces : % f % f % f\n", tr.force(eX), tr.force(eY), tr.force(eZ) );
963 : : printf("# t moments : % f % f % f\n", tr.moment(eL), tr.moment(eM), tr.moment(eN) );
964 : : }
965 : :
966 : : // finally set vFn & vMn
967 : 0 : vFn = mr.force + tr.force;
968 : 0 : vMn = mr.moment + tr.moment;
969 : :
970 : : // and just lie here
971 : 0 : Thrust = 0.0;
972 : :
973 : : // return unmodified thrust to the turbine.
974 : : // :TK: As far as I can see the return value is unused.
975 : 0 : return PowerAvailable;
976 : :
977 : : } // Calculate
978 : :
979 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
980 : :
981 : : // FGThruster does return 0.0 (the implicit direct thruster)
982 : : // piston CALL: return Thruster->Calculate((Eng_HP * hptoftlbssec)-Thruster->GetPowerRequired());
983 : :
984 : 0 : double FGRotor::GetPowerRequired(void)
985 : : {
986 : 0 : PowerRequired = 0.0;
987 : 0 : return PowerRequired;
988 : : }
989 : :
990 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
991 : :
992 : 0 : bool FGRotor::bind(void) {
993 : :
994 : 0 : string property_name, base_property_name;
995 : 0 : base_property_name = CreateIndexedPropertyName("propulsion/engine", EngineNum);
996 : :
997 : 0 : PropertyManager->Tie( base_property_name + "/rotor-rpm", this, &FGRotor::GetRPM );
998 : 0 : PropertyManager->Tie( base_property_name + "/thrust-mr-lbs", &mr.Thrust );
999 : 0 : PropertyManager->Tie( base_property_name + "/vi-mr-fps", &mr.v_induced );
1000 : 0 : PropertyManager->Tie( base_property_name + "/a0-mr-rad", &mr.a0 );
1001 : 0 : PropertyManager->Tie( base_property_name + "/a1-mr-rad", &mr.a1s ); // s means shaft axes
1002 : 0 : PropertyManager->Tie( base_property_name + "/b1-mr-rad", &mr.b1s );
1003 : 0 : PropertyManager->Tie( base_property_name + "/thrust-tr-lbs", &tr.Thrust );
1004 : 0 : PropertyManager->Tie( base_property_name + "/vi-tr-fps", &tr.v_induced );
1005 : :
1006 : : // lambda
1007 : 0 : PropertyManager->Tie( base_property_name + "/inflow-ratio", &prop_inflow_ratio_lambda );
1008 : : // mu
1009 : 0 : PropertyManager->Tie( base_property_name + "/advance-ratio", &prop_advance_ratio_mu );
1010 : : // nu
1011 : 0 : PropertyManager->Tie( base_property_name + "/induced-inflow-ratio", &prop_inflow_ratio_induced_nu );
1012 : :
1013 : 0 : PropertyManager->Tie( base_property_name + "/torque-mr-lbsft", &prop_mr_torque );
1014 : 0 : PropertyManager->Tie( base_property_name + "/thrust-coefficient", &prop_thrust_coefficient );
1015 : 0 : PropertyManager->Tie( base_property_name + "/main-rotor-rpm", &mr.ActualRPM );
1016 : 0 : PropertyManager->Tie( base_property_name + "/tail-rotor-rpm", &tr.ActualRPM );
1017 : :
1018 : : // position of the downwash
1019 : 0 : PropertyManager->Tie( base_property_name + "/theta-downwash-rad", &prop_theta_downwash );
1020 : 0 : PropertyManager->Tie( base_property_name + "/phi-downwash-rad", &prop_phi_downwash );
1021 : :
1022 : : // nodes to use via get<xyz>Value
1023 : 0 : prop_collective_ctrl = PropertyManager->GetNode(base_property_name + "/collective-ctrl-rad",true);
1024 : 0 : prop_lateral_ctrl = PropertyManager->GetNode(base_property_name + "/lateral-ctrl-rad",true);
1025 : 0 : prop_longitudinal_ctrl = PropertyManager->GetNode(base_property_name + "/longitudinal-ctrl-rad",true);
1026 : 0 : prop_antitorque_ctrl = PropertyManager->GetNode(base_property_name + "/antitorque-ctrl-rad",true);
1027 : :
1028 : 0 : prop_rotorbrake = PropertyManager->GetNode(base_property_name + "/rotorbrake-hp", true);
1029 : 0 : prop_freewheel_factor = PropertyManager->GetNode(base_property_name + "/freewheel-factor", true);
1030 : :
1031 : 0 : PropertyManager->Tie( base_property_name + "/dump-flag", &prop_DumpFlag );
1032 : :
1033 : 0 : return true;
1034 : : }
1035 : :
1036 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1037 : :
1038 : 0 : string FGRotor::GetThrusterLabels(int id, string delimeter)
1039 : : {
1040 : :
1041 : 0 : std::ostringstream buf;
1042 : :
1043 : 0 : buf << Name << " RPM (engine " << id << ")";
1044 : :
1045 : 0 : return buf.str();
1046 : :
1047 : : }
1048 : :
1049 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1050 : :
1051 : 0 : string FGRotor::GetThrusterValues(int id, string delimeter)
1052 : : {
1053 : 0 : std::ostringstream buf;
1054 : :
1055 : 0 : buf << RPM;
1056 : :
1057 : 0 : return buf.str();
1058 : : }
1059 : :
1060 : : //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1061 : : // The bitmasked value choices are as follows:
1062 : : // unset: In this case (the default) JSBSim would only print
1063 : : // out the normally expected messages, essentially echoing
1064 : : // the config files as they are read. If the environment
1065 : : // variable is not set, debug_lvl is set to 1 internally
1066 : : // 0: This requests JSBSim not to output any messages
1067 : : // whatsoever.
1068 : : // 1: This value explicity requests the normal JSBSim
1069 : : // startup messages
1070 : : // 2: This value asks for a message to be printed out when
1071 : : // a class is instantiated
1072 : : // 4: When this value is set, a message is displayed when a
1073 : : // FGModel object executes its Run() method
1074 : : // 8: When this value is set, various runtime state variables
1075 : : // are printed out periodically
1076 : : // 16: When set various parameters are sanity checked and
1077 : : // a message is printed out when they go out of bounds
1078 : :
1079 : 0 : void FGRotor::Debug(int from)
1080 : : {
1081 [ # # ]: 0 : if (debug_lvl <= 0) return;
1082 : :
1083 [ # # ]: 0 : if (debug_lvl & 1) { // Standard console startup message output
1084 [ # # ]: 0 : if (from == 0) { // Constructor
1085 : 0 : cout << "\n Rotor Name: " << Name << endl;
1086 : : }
1087 : : }
1088 [ # # ]: 0 : if (debug_lvl & 2 ) { // Instantiation/Destruction notification
1089 [ # # ]: 0 : if (from == 0) cout << "Instantiated: FGRotor" << endl;
1090 [ # # ]: 0 : if (from == 1) cout << "Destroyed: FGRotor" << endl;
1091 : : }
1092 : 0 : if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
1093 : : }
1094 : 0 : if (debug_lvl & 8 ) { // Runtime state variables
1095 : : }
1096 : 0 : if (debug_lvl & 16) { // Sanity checking
1097 : : }
1098 [ # # ]: 0 : if (debug_lvl & 64) {
1099 [ # # ]: 0 : if (from == 0) { // Constructor
1100 : 0 : cout << IdSrc << endl;
1101 : 0 : cout << IdHdr << endl;
1102 : : }
1103 : : }
1104 : : }
1105 : :
1106 [ + + ][ + - ]: 12 : } // namespace JSBSim
1107 : :
|
