JSBSim Flight Dynamics Model  1.0 (02 March 2017)
An Open Source Flight Dynamics and Control Software Library in C++
FGLGear.cpp
1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 
3  Module: FGLGear.cpp
4  Author: Jon S. Berndt
5  Norman H. Princen
6  Bertrand Coconnier
7  Date started: 11/18/99
8  Purpose: Encapsulates the landing gear elements
9  Called by: FGAircraft
10 
11  ------------- Copyright (C) 1999 Jon S. Berndt (jon@jsbsim.org) -------------
12 
13  This program is free software; you can redistribute it and/or modify it under
14  the terms of the GNU Lesser General Public License as published by the Free Software
15  Foundation; either version 2 of the License, or (at your option) any later
16  version.
17 
18  This program is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20  FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
21  details.
22 
23  You should have received a copy of the GNU Lesser General Public License along with
24  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
25  Place - Suite 330, Boston, MA 02111-1307, USA.
26 
27  Further information about the GNU Lesser General Public License can also be found on
28  the world wide web at http://www.gnu.org.
29 
30 FUNCTIONAL DESCRIPTION
31 --------------------------------------------------------------------------------
32 
33 HISTORY
34 --------------------------------------------------------------------------------
35 11/18/99 JSB Created
36 01/30/01 NHP Extended gear model to properly simulate steering and braking
37 07/08/09 BC Modified gear model to support large angles between aircraft and ground
38 
39 /%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40 INCLUDES
41 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
42 
43 #include <cstdlib>
44 #include <cstring>
45 
46 #include "math/FGFunction.h"
47 #include "FGLGear.h"
48 #include "models/FGGroundReactions.h"
49 #include "math/FGTable.h"
50 #include "input_output/FGXMLElement.h"
51 
52 using namespace std;
53 
54 namespace JSBSim {
55 
56 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57 DEFINITIONS
58 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
59 
60 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61 GLOBAL DATA
62 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
63 
64 IDENT(IdSrc,"$Id: FGLGear.cpp,v 1.125 2017/02/21 21:14:13 bcoconni Exp $");
65 IDENT(IdHdr,ID_LGEAR);
66 
67 // Body To Structural (body frame is rotated 180 deg about Y and lengths are given in
68 // ft instead of inches)
69 const FGMatrix33 FGLGear::Tb2s(-1./inchtoft, 0., 0., 0., 1./inchtoft, 0., 0., 0., -1./inchtoft);
70 const FGMatrix33 FGLGear::Ts2b(-inchtoft, 0., 0., 0., inchtoft, 0., 0., 0., -inchtoft);
71 
72 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
73 CLASS IMPLEMENTATION
74 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
75 
76 FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number, const struct Inputs& inputs) :
77  FGSurface(fdmex, number),
78  FGForce(fdmex),
79  in(inputs),
80  GearNumber(number),
81  SteerAngle(0.0),
82  Castered(false),
83  StaticFriction(false),
84  eSteerType(stSteer)
85 {
86  kSpring = bDamp = bDampRebound = dynamicFCoeff = staticFCoeff = rollingFCoeff = maxSteerAngle = 0;
87  isRetractable = false;
88  eDampType = dtLinear;
89  eDampTypeRebound = dtLinear;
90 
91  name = el->GetAttributeValue("name");
92  string sContactType = el->GetAttributeValue("type");
93  if (sContactType == "BOGEY") {
94  eContactType = ctBOGEY;
95  } else if (sContactType == "STRUCTURE") {
96  eContactType = ctSTRUCTURE;
97  } else {
98  // Unknown contact point types will be treated as STRUCTURE.
99  eContactType = ctSTRUCTURE;
100  }
101 
102  // Default values for structural contact points
103  if (eContactType == ctSTRUCTURE) {
104  kSpring = in.EmptyWeight;
105  bDamp = kSpring;
106  bDampRebound = kSpring * 10;
107  staticFCoeff = 1.0;
108  dynamicFCoeff = 1.0;
109  }
110 
111  PropertyManager = fdmex->GetPropertyManager();
112 
113  fStrutForce = 0;
114  Element* strutForce = el->FindElement("strut_force");
115  if (strutForce) {
116  Element* springFunc = strutForce->FindElement("function");
117  fStrutForce = new FGFunction(PropertyManager, springFunc);
118  }
119  else {
120  if (el->FindElement("spring_coeff"))
121  kSpring = el->FindElementValueAsNumberConvertTo("spring_coeff", "LBS/FT");
122  if (el->FindElement("damping_coeff")) {
123  Element* dampCoeff = el->FindElement("damping_coeff");
124  if (dampCoeff->GetAttributeValue("type") == "SQUARE") {
125  eDampType = dtSquare;
126  bDamp = el->FindElementValueAsNumberConvertTo("damping_coeff", "LBS/FT2/SEC2");
127  } else {
128  bDamp = el->FindElementValueAsNumberConvertTo("damping_coeff", "LBS/FT/SEC");
129  }
130  }
131 
132  if (el->FindElement("damping_coeff_rebound")) {
133  Element* dampCoeffRebound = el->FindElement("damping_coeff_rebound");
134  if (dampCoeffRebound->GetAttributeValue("type") == "SQUARE") {
135  eDampTypeRebound = dtSquare;
136  bDampRebound = el->FindElementValueAsNumberConvertTo("damping_coeff_rebound", "LBS/FT2/SEC2");
137  } else {
138  bDampRebound = el->FindElementValueAsNumberConvertTo("damping_coeff_rebound", "LBS/FT/SEC");
139  }
140  } else {
141  bDampRebound = bDamp;
142  eDampTypeRebound = eDampType;
143  }
144  }
145 
146  if (el->FindElement("dynamic_friction"))
147  dynamicFCoeff = el->FindElementValueAsNumber("dynamic_friction");
148  if (el->FindElement("static_friction"))
149  staticFCoeff = el->FindElementValueAsNumber("static_friction");
150  if (el->FindElement("rolling_friction"))
151  rollingFCoeff = el->FindElementValueAsNumber("rolling_friction");
152  if (el->FindElement("retractable"))
153  isRetractable = ((unsigned int)el->FindElementValueAsNumber("retractable"))>0.0?true:false;
154 
155  if (el->FindElement("max_steer"))
156  maxSteerAngle = el->FindElementValueAsNumberConvertTo("max_steer", "DEG");
157 
158  Element* castered_el = el->FindElement("castered");
159 
160  if ((maxSteerAngle == 360 && !castered_el)
161  || (castered_el && castered_el->GetDataAsNumber() != 0.0)) {
162  eSteerType = stCaster;
163  Castered = true;
164  }
165  else if (maxSteerAngle == 0.0) {
166  eSteerType = stFixed;
167  }
168  else
169  eSteerType = stSteer;
170 
171  GroundReactions = fdmex->GetGroundReactions();
172 
173  ForceY_Table = 0;
174  Element* force_table = el->FindElement("table");
175  while (force_table) {
176  string force_type = force_table->GetAttributeValue("type");
177  if (force_type == "CORNERING_COEFF") {
178  ForceY_Table = new FGTable(PropertyManager, force_table);
179  break;
180  } else {
181  cerr << "Undefined force table for " << name << " contact point" << endl;
182  }
183  force_table = el->FindNextElement("table");
184  }
185 
186  Element* element = el->FindElement("location");
187  if (element) vXYZn = element->FindElementTripletConvertTo("IN");
188  else {cerr << "No location given for contact " << name << endl; exit(-1);}
189  SetTransformType(FGForce::tCustom);
190 
191  element = el->FindElement("orientation");
192  if (element && (eContactType == ctBOGEY)) {
193  FGQuaternion quatFromEuler(element->FindElementTripletConvertTo("RAD"));
194 
195  mTGear = quatFromEuler.GetT();
196  }
197  else {
198  mTGear(1,1) = 1.;
199  mTGear(2,2) = 1.;
200  mTGear(3,3) = 1.;
201  }
202 
203  string sBrakeGroup = el->FindElementValue("brake_group");
204 
205  if (sBrakeGroup == "LEFT" ) eBrakeGrp = bgLeft;
206  else if (sBrakeGroup == "RIGHT" ) eBrakeGrp = bgRight;
207  else if (sBrakeGroup == "CENTER") eBrakeGrp = bgCenter;
208  else if (sBrakeGroup == "NOSE" ) eBrakeGrp = bgCenter; // Nose brake is not supported by FGFCS
209  else if (sBrakeGroup == "TAIL" ) eBrakeGrp = bgCenter; // Tail brake is not supported by FGFCS
210  else if (sBrakeGroup == "NONE" ) eBrakeGrp = bgNone;
211  else if (sBrakeGroup.empty() ) eBrakeGrp = bgNone;
212  else {
213  cerr << "Improper braking group specification in config file: "
214  << sBrakeGroup << " is undefined." << endl;
215  }
216 
217 // Add some AI here to determine if gear is located properly according to its
218 // brake group type ??
219 
220  useFCSGearPos = false;
221  ReportEnable = true;
222  TakeoffReported = LandingReported = false;
223 
224  // Set Pacejka terms
225 
226  Stiffness = 0.06;
227  Shape = 2.8;
228  Peak = staticFCoeff;
229  Curvature = 1.03;
230 
231  ResetToIC();
232 
233  Debug(0);
234 }
235 
236 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
237 
239 {
240  delete ForceY_Table;
241  delete fStrutForce;
242  Debug(1);
243 }
244 
245 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
246 
247 void FGLGear::ResetToIC(void)
248 {
249  GearPos = 1.0;
250 
251  WOW = lastWOW = false;
252  FirstContact = false;
253  StartedGroundRun = false;
254  LandingDistanceTraveled = TakeoffDistanceTraveled = TakeoffDistanceTraveled50ft = 0.0;
255  MaximumStrutForce = MaximumStrutTravel = 0.0;
256  SinkRate = GroundSpeed = 0.0;
257  SteerAngle = 0.0;
258 
259  vWhlVelVec.InitMatrix();
260 
261  compressLength = 0.0;
262  compressSpeed = 0.0;
263  maxCompLen = 0.0;
264 
265  WheelSlip = 0.0;
266 
267  // Initialize Lagrange multipliers
268  for (int i=0; i < 3; i++) {
269  LMultiplier[i].ForceJacobian.InitMatrix();
270  LMultiplier[i].MomentJacobian.InitMatrix();
271  LMultiplier[i].Min = 0.0;
272  LMultiplier[i].Max = 0.0;
273  LMultiplier[i].value = 0.0;
274  }
275 }
276 
277 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
278 
279 const FGColumnVector3& FGLGear::GetBodyForces(FGSurface *surface)
280 {
281  double gearPos = 1.0;
282 
283  vFn.InitMatrix();
284 
285  if (isRetractable) gearPos = GetGearUnitPos();
286 
287  if (gearPos > 0.99) { // Gear DOWN
288  FGColumnVector3 normal, terrainVel, dummy;
289  FGLocation gearLoc, contact;
290  FGColumnVector3 vWhlBodyVec = Ts2b * (vXYZn - in.vXYZcg);
291 
292  vLocalGear = in.Tb2l * vWhlBodyVec; // Get local frame wheel location
293  gearLoc = in.Location.LocalToLocation(vLocalGear);
294 
295  // Compute the height of the theoretical location of the wheel (if strut is
296  // not compressed) with respect to the ground level
297  double height = gearLoc.GetContactPoint(contact, normal, terrainVel, dummy);
298 
299  // Does this surface contact point interact with another surface?
300  if (surface) {
301  height -= (*surface).GetBumpHeight();
302  staticFFactor = (*surface).GetStaticFFactor();
303  rollingFFactor = (*surface).GetRollingFFactor();
304  maximumForce = (*surface).GetMaximumForce();
305  isSolid = (*surface).GetSolid();
306  }
307 
308  if (height < 0.0) {
309  WOW = isSolid;
310  vGroundNormal = in.Tec2b * normal;
311 
312  // The height returned by GetGroundCallback() is the AGL and is expressed
313  // in the Z direction of the local coordinate frame. We now need to transform
314  // this height in actual compression of the strut (BOGEY) or in the normal
315  // direction to the ground (STRUCTURE)
316  double normalZ = (in.Tec2l*normal)(eZ);
317  double LGearProj = -(mTGear.Transposed() * vGroundNormal)(eZ);
318  FGColumnVector3 vWhlDisplVec;
319 
320  // The following equations use the vector to the tire contact patch
321  // including the strut compression.
322  switch(eContactType) {
323  case ctBOGEY:
324  if (isSolid) {
325  compressLength = LGearProj > 0.0 ? height * normalZ / LGearProj : 0.0;
326  vWhlDisplVec = mTGear * FGColumnVector3(0., 0., -compressLength);
327  } else {
328  // Gears don't (or hardly) compress in liquids
329  compressLength = 0.0;
330  vWhlDisplVec = 0.0 * vGroundNormal;
331  }
332  break;
333  case ctSTRUCTURE:
334  compressLength = height * normalZ / DotProduct(normal, normal);
335  vWhlDisplVec = compressLength * vGroundNormal;
336  break;
337  }
338 
339  FGColumnVector3 vWhlContactVec = vWhlBodyVec + vWhlDisplVec;
340  vActingXYZn = vXYZn + Tb2s * vWhlDisplVec;
341  FGColumnVector3 vBodyWhlVel = in.PQR * vWhlContactVec;
342  vBodyWhlVel += in.UVW - in.Tec2b * terrainVel;
343 
344  if (isSolid) {
345  vWhlVelVec = mTGear.Transposed() * vBodyWhlVel;
346  } else {
347  // wheels don't spin up in liquids: let wheel spin down slowly
348  vWhlVelVec(eX) -= 13.0 * in.TotalDeltaT;
349  if (vWhlVelVec(eX) < 0.0) vWhlVelVec(eX) = 0.0;
350  }
351 
352  InitializeReporting();
353  ComputeSteeringAngle();
354  ComputeGroundFrame();
355 
356  vGroundWhlVel = mT.Transposed() * vBodyWhlVel;
357 
358  if (fdmex->GetTrimStatus())
359  compressSpeed = 0.0; // Steady state is sought during trimming
360  else {
361  compressSpeed = -vGroundWhlVel(eZ);
362  if (eContactType == ctBOGEY)
363  compressSpeed /= LGearProj;
364  }
365 
366  ComputeVerticalStrutForce();
367 
368  // Compute the friction coefficients in the wheel ground plane.
369  if (eContactType == ctBOGEY) {
370  ComputeSlipAngle();
371  ComputeBrakeForceCoefficient();
372  ComputeSideForceCoefficient();
373  }
374 
375  // Prepare the Jacobians and the Lagrange multipliers for later friction
376  // forces calculations.
377  ComputeJacobian(vWhlContactVec);
378 
379  } else { // Gear is NOT compressed
380 
381  WOW = false;
382  compressLength = 0.0;
383  compressSpeed = 0.0;
384  WheelSlip = 0.0;
385  StrutForce = 0.0;
386 
387  LMultiplier[ftRoll].value = 0.0;
388  LMultiplier[ftSide].value = 0.0;
389  LMultiplier[ftDynamic].value = 0.0;
390 
391  // Let wheel spin down slowly
392  vWhlVelVec(eX) -= 13.0 * in.TotalDeltaT;
393  if (vWhlVelVec(eX) < 0.0) vWhlVelVec(eX) = 0.0;
394 
395  // Return to neutral position between 1.0 and 0.8 gear pos.
396  SteerAngle *= max(gearPos-0.8, 0.0)/0.2;
397 
398  ResetReporting();
399  }
400 
401  } else if (gearPos < 0.01) { // Gear UP
402 
403  WOW = false;
404  vWhlVelVec.InitMatrix();
405  }
406 
407  if (!fdmex->GetTrimStatus()) {
408  ReportTakeoffOrLanding();
409 
410  // Require both WOW and LastWOW to be true before checking crash conditions
411  // to allow the WOW flag to be used in terminating a scripted run.
412  if (WOW && lastWOW) CrashDetect();
413 
414  lastWOW = WOW;
415  }
416 
417  return FGForce::GetBodyForces();
418 }
419 
420 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
421 // Build a local "ground" coordinate system defined by
422 // eX : projection of the rolling direction on the ground
423 // eY : projection of the sliping direction on the ground
424 // eZ : normal to the ground
425 
426 void FGLGear::ComputeGroundFrame(void)
427 {
428  FGColumnVector3 roll = mTGear * FGColumnVector3(cos(SteerAngle), sin(SteerAngle), 0.);
429  FGColumnVector3 side = vGroundNormal * roll;
430 
431  roll -= DotProduct(roll, vGroundNormal) * vGroundNormal;
432  roll.Normalize();
433  side.Normalize();
434 
435  mT(eX,eX) = roll(eX);
436  mT(eY,eX) = roll(eY);
437  mT(eZ,eX) = roll(eZ);
438  mT(eX,eY) = side(eX);
439  mT(eY,eY) = side(eY);
440  mT(eZ,eY) = side(eZ);
441  mT(eX,eZ) = vGroundNormal(eX);
442  mT(eY,eZ) = vGroundNormal(eY);
443  mT(eZ,eZ) = vGroundNormal(eZ);
444 }
445 
446 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
447 // Calculate tire slip angle.
448 
449 void FGLGear::ComputeSlipAngle(void)
450 {
451 // Check that the speed is non-null otherwise keep the current angle
452  if (vGroundWhlVel.Magnitude(eX,eY) > 1E-3)
453  WheelSlip = -atan2(vGroundWhlVel(eY), fabs(vGroundWhlVel(eX)))*radtodeg;
454 }
455 
456 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
457 // Compute the steering angle in any case.
458 // This will also make sure that animations will look right.
459 
460 void FGLGear::ComputeSteeringAngle(void)
461 {
462  if (Castered) {
463  // Check that the speed is non-null otherwise keep the current angle
464  if (vWhlVelVec.Magnitude(eX,eY) > 0.1)
465  SteerAngle = atan2(vWhlVelVec(eY), fabs(vWhlVelVec(eX)));
466  }
467 }
468 
469 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
470 // Reset reporting functionality after takeoff
471 
472 void FGLGear::ResetReporting(void)
473 {
474  if (in.DistanceAGL > 200.0) {
475  FirstContact = false;
476  StartedGroundRun = false;
477  LandingReported = false;
478  TakeoffReported = true;
479  LandingDistanceTraveled = 0.0;
480  MaximumStrutForce = MaximumStrutTravel = 0.0;
481  }
482 }
483 
484 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
485 
486 void FGLGear::InitializeReporting(void)
487 {
488  // If this is the first time the wheel has made contact, remember some values
489  // for later printout.
490 
491  if (!FirstContact) {
492  FirstContact = true;
493  SinkRate = compressSpeed;
494  GroundSpeed = in.Vground;
495  TakeoffReported = false;
496  }
497 
498  // If the takeoff run is starting, initialize.
499 
500  if ((in.Vground > 0.1) &&
501  (in.BrakePos[bgLeft] == 0) &&
502  (in.BrakePos[bgRight] == 0) &&
503  (in.TakeoffThrottle && !StartedGroundRun))
504  {
505  TakeoffDistanceTraveled = 0;
506  TakeoffDistanceTraveled50ft = 0;
507  StartedGroundRun = true;
508  }
509 }
510 
511 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
512 // Takeoff and landing reporting functionality
513 
514 void FGLGear::ReportTakeoffOrLanding(void)
515 {
516  if (FirstContact)
517  LandingDistanceTraveled += in.Vground * in.TotalDeltaT;
518 
519  if (StartedGroundRun) {
520  TakeoffDistanceTraveled50ft += in.Vground * in.TotalDeltaT;
521  if (WOW) TakeoffDistanceTraveled += in.Vground * in.TotalDeltaT;
522  }
523 
524  if ( ReportEnable
525  && in.Vground <= 0.05
526  && !LandingReported
527  && in.WOW)
528  {
529  if (debug_lvl > 0) Report(erLand);
530  }
531 
532  if ( ReportEnable
533  && !TakeoffReported
534  && (in.DistanceAGL - vLocalGear(eZ)) > 50.0
535  && !in.WOW)
536  {
537  if (debug_lvl > 0) Report(erTakeoff);
538  }
539 
540  if (lastWOW != WOW)
541  {
542  ostringstream buf;
543  buf << "GEAR_CONTACT: " << fdmex->GetSimTime() << " seconds: " << name;
544  PutMessage(buf.str(), WOW);
545  }
546 }
547 
548 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
549 // Crash detection logic (really out-of-bounds detection)
550 
551 void FGLGear::CrashDetect(void)
552 {
553  if ( (compressLength > 500.0 ||
554  vFn.Magnitude() > 100000000.0 ||
555  GetMoments().Magnitude() > 5000000000.0 ||
556  SinkRate > 1.4666*30 ) && !fdmex->IntegrationSuspended())
557  {
558  ostringstream buf;
559  buf << "*CRASH DETECTED* " << fdmex->GetSimTime() << " seconds: " << name;
560  PutMessage(buf.str());
561  // fdmex->SuspendIntegration();
562  }
563 }
564 
565 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
566 // The following needs work regarding friction coefficients and braking and
567 // steering The BrakeFCoeff formula assumes that an anti-skid system is used.
568 // It also assumes that we won't be turning and braking at the same time.
569 // Will fix this later.
570 // [JSB] The braking force coefficients include normal rolling coefficient +
571 // a percentage of the static friction coefficient based on braking applied.
572 
573 void FGLGear::ComputeBrakeForceCoefficient(void)
574 {
575  BrakeFCoeff = rollingFFactor * rollingFCoeff;
576 
577  if (eBrakeGrp != bgNone)
578  BrakeFCoeff += in.BrakePos[eBrakeGrp] * staticFFactor * (staticFCoeff - rollingFCoeff);
579 }
580 
581 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
582 // Compute the sideforce coefficients using Pacejka's Magic Formula.
583 //
584 // y(x) = D sin {C arctan [Bx - E(Bx - arctan Bx)]}
585 //
586 // Where: B = Stiffness Factor (0.06, here)
587 // C = Shape Factor (2.8, here)
588 // D = Peak Factor (0.8, here)
589 // E = Curvature Factor (1.03, here)
590 
591 void FGLGear::ComputeSideForceCoefficient(void)
592 {
593  if (ForceY_Table) {
594  FCoeff = ForceY_Table->GetValue(WheelSlip);
595  } else {
596  double StiffSlip = Stiffness*WheelSlip;
597  FCoeff = Peak * sin(Shape*atan(StiffSlip - Curvature*(StiffSlip - atan(StiffSlip))));
598  }
599  FCoeff *= staticFFactor;
600 }
601 
602 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
603 // Compute the vertical force on the wheel using square-law damping (per comment
604 // in paper AIAA-2000-4303 - see header prologue comments). We might consider
605 // allowing for both square and linear damping force calculation. Also need to
606 // possibly give a "rebound damping factor" that differs from the compression
607 // case.
608 
609 void FGLGear::ComputeVerticalStrutForce()
610 {
611  double springForce = 0;
612  double dampForce = 0;
613 
614  if (fStrutForce)
615  StrutForce = min(fStrutForce->GetValue(), (double)0.0);
616  else {
617  springForce = -compressLength * kSpring;
618 
619  if (compressSpeed >= 0.0) {
620 
621  if (eDampType == dtLinear)
622  dampForce = -compressSpeed * bDamp;
623  else
624  dampForce = -compressSpeed * compressSpeed * bDamp;
625 
626  } else {
627 
628  if (eDampTypeRebound == dtLinear)
629  dampForce = -compressSpeed * bDampRebound;
630  else
631  dampForce = compressSpeed * compressSpeed * bDampRebound;
632 
633  }
634 
635  StrutForce = min(springForce + dampForce, (double)0.0);
636  if (StrutForce > maximumForce) {
637  StrutForce = maximumForce;
638  compressLength = -StrutForce / kSpring;
639  }
640  }
641 
642  // The reaction force of the wheel is always normal to the ground
643  switch (eContactType) {
644  case ctBOGEY:
645  // Project back the strut force in the local coordinate frame of the ground
646  vFn(eZ) = StrutForce / (mTGear.Transposed()*vGroundNormal)(eZ);
647  break;
648  case ctSTRUCTURE:
649  vFn(eZ) = -StrutForce;
650  break;
651  }
652 
653  // Remember these values for reporting
654  MaximumStrutForce = max(MaximumStrutForce, fabs(StrutForce));
655  MaximumStrutTravel = max(MaximumStrutTravel, fabs(compressLength));
656 }
657 
658 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
659 
660 double FGLGear::GetGearUnitPos(void) const
661 {
662  // hack to provide backward compatibility to gear/gear-pos-norm property
663  if( useFCSGearPos || in.FCSGearPos != 1.0 ) {
664  useFCSGearPos = true;
665  return in.FCSGearPos;
666  }
667  return GearPos;
668 }
669 
670 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
671 // Compute the jacobian entries for the friction forces resolution later
672 // in FGPropagate
673 
674 void FGLGear::ComputeJacobian(const FGColumnVector3& vWhlContactVec)
675 {
676  // When the point of contact is moving, dynamic friction is used
677  // This type of friction is limited to ctSTRUCTURE elements because their
678  // friction coefficient is the same in every directions
679  if ((eContactType == ctSTRUCTURE) && (vGroundWhlVel.Magnitude(eX,eY) > 1E-3)) {
680 
681  FGColumnVector3 velocityDirection = vGroundWhlVel;
682 
683  StaticFriction = false;
684 
685  velocityDirection(eZ) = 0.;
686  velocityDirection.Normalize();
687 
688  LMultiplier[ftDynamic].ForceJacobian = mT * velocityDirection;
689  LMultiplier[ftDynamic].MomentJacobian = vWhlContactVec * LMultiplier[ftDynamic].ForceJacobian;
690  LMultiplier[ftDynamic].Max = 0.;
691  LMultiplier[ftDynamic].Min = -fabs(staticFFactor * dynamicFCoeff * vFn(eZ));
692 
693  // The Lagrange multiplier value obtained from the previous iteration is kept
694  // This is supposed to accelerate the convergence of the projected Gauss-Seidel
695  // algorithm. The code just below is to make sure that the initial value
696  // is consistent with the current friction coefficient and normal reaction.
697  LMultiplier[ftDynamic].value = Constrain(LMultiplier[ftDynamic].Min, LMultiplier[ftDynamic].value, LMultiplier[ftDynamic].Max);
698 
699  GroundReactions->RegisterLagrangeMultiplier(&LMultiplier[ftDynamic]);
700  }
701  else {
702  // Static friction is used for ctSTRUCTURE when the contact point is not moving.
703  // It is always used for ctBOGEY elements because the friction coefficients
704  // of a tyre depend on the direction of the movement (roll & side directions).
705  // This cannot be handled properly by the so-called "dynamic friction".
706  StaticFriction = true;
707 
708  LMultiplier[ftRoll].ForceJacobian = mT * FGColumnVector3(1.,0.,0.);
709  LMultiplier[ftSide].ForceJacobian = mT * FGColumnVector3(0.,1.,0.);
710  LMultiplier[ftRoll].MomentJacobian = vWhlContactVec * LMultiplier[ftRoll].ForceJacobian;
711  LMultiplier[ftSide].MomentJacobian = vWhlContactVec * LMultiplier[ftSide].ForceJacobian;
712 
713  switch(eContactType) {
714  case ctBOGEY:
715  LMultiplier[ftRoll].Max = fabs(BrakeFCoeff * vFn(eZ));
716  LMultiplier[ftSide].Max = fabs(FCoeff * vFn(eZ));
717  break;
718  case ctSTRUCTURE:
719  LMultiplier[ftRoll].Max = fabs(staticFFactor * staticFCoeff * vFn(eZ));
720  LMultiplier[ftSide].Max = LMultiplier[ftRoll].Max;
721  break;
722  }
723 
724  LMultiplier[ftRoll].Min = -LMultiplier[ftRoll].Max;
725  LMultiplier[ftSide].Min = -LMultiplier[ftSide].Max;
726 
727  // The Lagrange multiplier value obtained from the previous iteration is kept
728  // This is supposed to accelerate the convergence of the projected Gauss-Seidel
729  // algorithm. The code just below is to make sure that the initial value
730  // is consistent with the current friction coefficient and normal reaction.
731  LMultiplier[ftRoll].value = Constrain(LMultiplier[ftRoll].Min, LMultiplier[ftRoll].value, LMultiplier[ftRoll].Max);
732  LMultiplier[ftSide].value = Constrain(LMultiplier[ftSide].Min, LMultiplier[ftSide].value, LMultiplier[ftSide].Max);
733 
734  GroundReactions->RegisterLagrangeMultiplier(&LMultiplier[ftRoll]);
735  GroundReactions->RegisterLagrangeMultiplier(&LMultiplier[ftSide]);
736  }
737 }
738 
739 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
740 // This routine is called after the Lagrange multiplier has been computed in
741 // the FGAccelerations class. The friction forces of the landing gear are then
742 // updated accordingly.
743 void FGLGear::UpdateForces(void)
744 {
745  if (StaticFriction) {
746  vFn(eX) = LMultiplier[ftRoll].value;
747  vFn(eY) = LMultiplier[ftSide].value;
748  }
749  else {
750  FGColumnVector3 forceDir = mT.Transposed() * LMultiplier[ftDynamic].ForceJacobian;
751  vFn(eX) = LMultiplier[ftDynamic].value * forceDir(eX);
752  vFn(eY) = LMultiplier[ftDynamic].value * forceDir(eY);
753  }
754 }
755 
756 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
757 
758 void FGLGear::SetstaticFCoeff(double coeff)
759 {
760  staticFCoeff = coeff;
761  Peak = coeff;
762 }
763 
764 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
765 
766 void FGLGear::bind(void)
767 {
768  string property_name;
769  string base_property_name;
770 
771  switch(eContactType) {
772  case ctBOGEY:
773  eSurfaceType = FGSurface::ctBOGEY;
774  base_property_name = CreateIndexedPropertyName("gear/unit", GearNumber);
775  break;
776  case ctSTRUCTURE:
777  eSurfaceType = FGSurface::ctSTRUCTURE;
778  base_property_name = CreateIndexedPropertyName("contact/unit", GearNumber);
779  break;
780  default:
781  return;
782  }
783  FGSurface::bind();
784 
785  property_name = base_property_name + "/WOW";
786  PropertyManager->Tie( property_name.c_str(), &WOW );
787  property_name = base_property_name + "/x-position";
788  PropertyManager->Tie( property_name.c_str(), (FGForce*)this,
789  &FGForce::GetLocationX, &FGForce::SetLocationX);
790  property_name = base_property_name + "/y-position";
791  PropertyManager->Tie( property_name.c_str(), (FGForce*)this,
792  &FGForce::GetLocationY, &FGForce::SetLocationY);
793  property_name = base_property_name + "/z-position";
794  PropertyManager->Tie( property_name.c_str(), (FGForce*)this,
795  &FGForce::GetLocationZ, &FGForce::SetLocationZ);
796  property_name = base_property_name + "/compression-ft";
797  PropertyManager->Tie( property_name.c_str(), &compressLength );
798  property_name = base_property_name + "/compression-velocity-fps";
799  PropertyManager->Tie( property_name.c_str(), &compressSpeed );
800  property_name = base_property_name + "/static_friction_coeff";
801  PropertyManager->Tie( property_name.c_str(), (FGLGear*)this,
802  &FGLGear::GetstaticFCoeff, &FGLGear::SetstaticFCoeff);
803  property_name = base_property_name + "/dynamic_friction_coeff";
804  PropertyManager->Tie( property_name.c_str(), &dynamicFCoeff );
805 
806  if (eContactType == ctBOGEY) {
807  property_name = base_property_name + "/slip-angle-deg";
808  PropertyManager->Tie( property_name.c_str(), &WheelSlip );
809  property_name = base_property_name + "/wheel-speed-fps";
810  PropertyManager->Tie( property_name.c_str(), (FGLGear*)this,
811  &FGLGear::GetWheelRollVel);
812  property_name = base_property_name + "/side_friction_coeff";
813  PropertyManager->Tie( property_name.c_str(), &FCoeff );
814  property_name = base_property_name + "/rolling_friction_coeff";
815  PropertyManager->Tie( property_name.c_str(), &rollingFCoeff );
816 
817  if (eSteerType == stCaster) {
818  property_name = base_property_name + "/steering-angle-deg";
819  PropertyManager->Tie( property_name.c_str(), this, &FGLGear::GetSteerAngleDeg );
820  property_name = base_property_name + "/castered";
821  PropertyManager->Tie( property_name.c_str(), &Castered);
822  }
823  }
824 
825  if( isRetractable ) {
826  property_name = base_property_name + "/pos-norm";
827  PropertyManager->Tie( property_name.c_str(), &GearPos );
828  }
829 
830  if (eSteerType != stFixed) {
831  // This property allows the FCS to override the steering position angle that
832  // is set by the property fcs/steer-cmd-norm. The prefix fcs/ has been kept
833  // for backward compatibility.
834  string tmp = CreateIndexedPropertyName("fcs/steer-pos-deg", GearNumber);
835  PropertyManager->Tie(tmp.c_str(), this, &FGLGear::GetSteerAngleDeg, &FGLGear::SetSteerAngleDeg);
836  }
837 }
838 
839 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
840 
841 void FGLGear::Report(ReportType repType)
842 {
843  if (fabs(TakeoffDistanceTraveled) < 0.001) return; // Don't print superfluous reports
844 
845  switch(repType) {
846  case erLand:
847  cout << endl << "Touchdown report for " << name << " (WOW at time: "
848  << fdmex->GetSimTime() << " seconds)" << endl;
849  cout << " Sink rate at contact: " << SinkRate << " fps, "
850  << SinkRate*0.3048 << " mps" << endl;
851  cout << " Contact ground speed: " << GroundSpeed*.5925 << " knots, "
852  << GroundSpeed*0.3048 << " mps" << endl;
853  cout << " Maximum contact force: " << MaximumStrutForce << " lbs, "
854  << MaximumStrutForce*4.448 << " Newtons" << endl;
855  cout << " Maximum strut travel: " << MaximumStrutTravel*12.0 << " inches, "
856  << MaximumStrutTravel*30.48 << " cm" << endl;
857  cout << " Distance traveled: " << LandingDistanceTraveled << " ft, "
858  << LandingDistanceTraveled*0.3048 << " meters" << endl;
859  LandingReported = true;
860  break;
861  case erTakeoff:
862  cout << endl << "Takeoff report for " << name << " (Liftoff at time: "
863  << fdmex->GetSimTime() << " seconds)" << endl;
864  cout << " Distance traveled: " << TakeoffDistanceTraveled
865  << " ft, " << TakeoffDistanceTraveled*0.3048 << " meters" << endl;
866  cout << " Distance traveled (over 50'): " << TakeoffDistanceTraveled50ft
867  << " ft, " << TakeoffDistanceTraveled50ft*0.3048 << " meters" << endl;
868  cout << " [Altitude (ASL): " << in.DistanceASL << " ft. / "
869  << in.DistanceASL*FGJSBBase::fttom << " m | Temperature: "
870  << in.Temperature - 459.67 << " F / "
871  << RankineToCelsius(in.Temperature) << " C]" << endl;
872  cout << " [Velocity (KCAS): " << in.VcalibratedKts << "]" << endl;
873  TakeoffReported = true;
874  break;
875  case erNone:
876  break;
877  }
878 }
879 
880 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
881 // The bitmasked value choices are as follows:
882 // unset: In this case (the default) JSBSim would only print
883 // out the normally expected messages, essentially echoing
884 // the config files as they are read. If the environment
885 // variable is not set, debug_lvl is set to 1 internally
886 // 0: This requests JSBSim not to output any messages
887 // whatsoever.
888 // 1: This value explicity requests the normal JSBSim
889 // startup messages
890 // 2: This value asks for a message to be printed out when
891 // a class is instantiated
892 // 4: When this value is set, a message is displayed when a
893 // FGModel object executes its Run() method
894 // 8: When this value is set, various runtime state variables
895 // are printed out periodically
896 // 16: When set various parameters are sanity checked and
897 // a message is printed out when they go out of bounds
898 
899 void FGLGear::Debug(int from)
900 {
901  static const char* sSteerType[] = {"STEERABLE", "FIXED", "CASTERED" };
902  static const char* sBrakeGroup[] = {"NONE", "LEFT", "RIGHT", "CENTER", "NOSE", "TAIL"};
903  static const char* sContactType[] = {"BOGEY", "STRUCTURE" };
904 
905  if (debug_lvl <= 0) return;
906 
907  if (debug_lvl & 1) { // Standard console startup message output
908  if (from == 0) { // Constructor - loading and initialization
909  cout << " " << sContactType[eContactType] << " " << name << endl;
910  cout << " Location: " << vXYZn << endl;
911  cout << " Spring Constant: " << kSpring << endl;
912 
913  if (eDampType == dtLinear)
914  cout << " Damping Constant: " << bDamp << " (linear)" << endl;
915  else
916  cout << " Damping Constant: " << bDamp << " (square law)" << endl;
917 
918  if (eDampTypeRebound == dtLinear)
919  cout << " Rebound Damping Constant: " << bDampRebound << " (linear)" << endl;
920  else
921  cout << " Rebound Damping Constant: " << bDampRebound << " (square law)" << endl;
922 
923  cout << " Dynamic Friction: " << dynamicFCoeff << endl;
924  cout << " Static Friction: " << staticFCoeff << endl;
925  if (eContactType == ctBOGEY) {
926  cout << " Rolling Friction: " << rollingFCoeff << endl;
927  cout << " Steering Type: " << sSteerType[eSteerType] << endl;
928  cout << " Grouping: " << sBrakeGroup[eBrakeGrp] << endl;
929  cout << " Max Steer Angle: " << maxSteerAngle << endl;
930  cout << " Retractable: " << isRetractable << endl;
931  }
932  }
933  }
934  if (debug_lvl & 2 ) { // Instantiation/Destruction notification
935  if (from == 0) cout << "Instantiated: FGLGear" << endl;
936  if (from == 1) cout << "Destroyed: FGLGear" << endl;
937  }
938  if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
939  }
940  if (debug_lvl & 8 ) { // Runtime state variables
941  }
942  if (debug_lvl & 16) { // Sanity checking
943  }
944  if (debug_lvl & 64) {
945  if (from == 0) { // Constructor
946  cout << IdSrc << endl;
947  cout << IdHdr << endl;
948  }
949  }
950 }
951 
952 } // namespace JSBSim
Models the Quaternion representation of rotations.
Definition: FGQuaternion.h:92
ReportType
Report type enumerators.
Definition: FGLGear.h:225
static double Constrain(double min, double value, double max)
Constrain a value between a minimum and a maximum value.
Definition: FGJSBBase.h:332
std::string GetAttributeValue(const std::string &key)
Retrieves an attribute.
static double RankineToCelsius(double rankine)
Converts from degrees Rankine to degrees Celsius.
Definition: FGJSBBase.h:211
double FindElementValueAsNumberConvertTo(const std::string &el, const std::string &target_units)
Searches for the named element and converts and returns the data belonging to it. ...
FGLocation holds an arbitrary location in the Earth centered Earth fixed reference frame (ECEF)...
Definition: FGLocation.h:160
STL namespace.
FGColumnVector3 & Normalize(void)
Normalize.
Element * FindElement(const std::string &el="")
Searches for a specified element.
double FindElementValueAsNumber(const std::string &el="")
Searches for the named element and returns the data belonging to it as a number.
double GetContactPoint(FGLocation &contact, FGColumnVector3 &normal, FGColumnVector3 &v, FGColumnVector3 &w) const
Get terrain contact point information below the current location.
Definition: FGLocation.h:378
void PutMessage(const Message &msg)
Places a Message structure on the Message queue.
Definition: FGJSBBase.cpp:123
double GetValue(void) const
Retrieves the value of the function object.
Definition: FGFunction.cpp:364
FGPropertyManager * GetPropertyManager(void)
Returns a pointer to the property manager object.
Definition: FGFDMExec.cpp:1099
void Tie(const std::string &name, bool *pointer, bool useDefault=true)
Tie a property to an external bool variable.
double GetDataAsNumber(void)
Converts the element data to a number.
FGGroundReactions * GetGroundReactions(void)
Returns the FGGroundReactions pointer.
Definition: FGFDMExec.h:361
FGLocation LocalToLocation(const FGColumnVector3 &lvec) const
Conversion from Local frame coordinates to a location in the earth centered and fixed frame...
Definition: FGLocation.h:467
std::string FindElementValue(const std::string &el="")
Searches for the named element and returns the string data belonging to it.
Represents a mathematical function.
Definition: FGFunction.h:699
This class implements a 3 element column vector.
bool IntegrationSuspended(void) const
Returns the simulation suspension state.
Definition: FGFDMExec.h:546
double GetSimTime(void) const
Returns the cumulative simulation time in seconds.
Definition: FGFDMExec.h:533
Element * FindNextElement(const std::string &el="")
Searches for the next element as specified.
FGMatrix33 Transposed(void) const
Transposed matrix.
Definition: FGMatrix33.h:243
double Magnitude(void) const
Length of the vector.
Base class for all surface properties.
Definition: FGSurface.h:67
Encapsulates the JSBSim simulation executive.
Definition: FGFDMExec.h:189
Utility class that aids in the conversion of forces between coordinate systems and calculation of mom...
Definition: FGForce.h:225
Lookup table class.
Definition: FGTable.h:243
FGColumnVector3 FindElementTripletConvertTo(const std::string &target_units)
Composes a 3-element column vector for the supplied location or orientation.
~FGLGear()
Destructor.
Definition: FGLGear.cpp:238
Landing gear model.
Definition: FGLGear.h:194