1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2014 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Authors: Justin Madsen
13 // =============================================================================
14 //
15 // Base class for a Pacejka type Magic formula 2002 tire model
16 //
17 // =============================================================================
18
19 #ifndef CH_PACEJKATIRE_H
20 #define CH_PACEJKATIRE_H
21
22 #include <fstream>
23 #include <sstream>
24 #include <string>
25 #include <vector>
26
27 #include "chrono/assets/ChCylinderShape.h"
28 #include "chrono/assets/ChTexture.h"
29 #include "chrono/physics/ChBody.h"
30
31 #include "chrono_vehicle/ChTerrain.h"
32 #include "chrono_vehicle/wheeled_vehicle/ChTire.h"
33 #include "chrono_vehicle/wheeled_vehicle/tire/ChPac2002_data.h"
34
35 namespace chrono {
36 namespace vehicle {
37
38 /// @addtogroup vehicle_wheeled_tire
39 /// @{
40
41 // Forward declarations for private structures
42 struct slips;
43 struct Pac2002_data;
44 struct pureLongCoefs;
45 struct pureLatCoefs;
46 struct pureTorqueCoefs;
47 struct combinedLongCoefs;
48 struct combinedLatCoefs;
49 struct combinedTorqueCoefs;
50 struct zetaCoefs;
51 struct relaxationL;
52 struct bessel;
53
54 /// Concrete tire class that implements the Pacejka tire model.
55 class CH_VEHICLE_API ChPacejkaTire : public ChTire {
56 public:
57 /// Default constructor for a Pacejka tire.
58 /// Construct a Pacejka tire for which the vertical load is calculated
59 /// internally. The model includes transient slip calculations.
60 /// chrono can suggest a time step for use with the ODE slips
61 ChPacejkaTire(const std::string& name, ///< [in] name of this tire
62 const std::string& pacTire_paramFile ///< [in] name of the parameter file
63 );
64
65 /// Construct a Pacejka tire with specified vertical load, for testing purposes
66 ChPacejkaTire(const std::string& name, ///< [in] name of this tire
67 const std::string& pacTire_paramFile, ///< [in] name of the parameter file
68 double Fz_override, ///< [in] prescribed vertical load
69 bool use_transient_slip = true ///< [in] indicate if using transient slip model
70 );
71
72 ~ChPacejkaTire();
73
74 /// Get the name of the vehicle subsystem template.
GetTemplateName()75 virtual std::string GetTemplateName() const override { return "PacejkaTire"; }
76
77 /// Specify whether or not the associated wheel is driven.
78 /// By default, the wheel is assumed not driven.
SetDrivenWheel(bool val)79 void SetDrivenWheel(bool val) { m_driven = val; }
80
81 /// Add visualization assets for the rigid tire subsystem.
82 virtual void AddVisualizationAssets(VisualizationType vis) override;
83
84 /// Remove visualization assets for the rigid tire subsystem.
85 virtual void RemoveVisualizationAssets() override;
86
87 /// Get the tire radius.
GetRadius()88 virtual double GetRadius() const override { return m_R_eff; }
89
90 /// Get tire width.
GetWidth()91 virtual double GetWidth() const override { return m_params->dimension.width; }
92
93 /// Get visualization tire width.
GetVisualizationWidth()94 virtual double GetVisualizationWidth() const { return 0.25; }
95
96 /// Report the tire force and moment.
97 virtual TerrainForce ReportTireForce(ChTerrain* terrain) const override;
98
99 /// Return the reactions for the pure slip EQs, in local or global coords
100 TerrainForce GetTireForce_pureSlip(const bool local = true) const;
101
102 /// Return the reactions for the combined slip EQs, in local or global coords
103 TerrainForce GetTireForce_combinedSlip(const bool local = true) const;
104
105 /// Get the tire slip angle computed internally by the Pacejka model (in radians).
106 /// The reported value will be the same as that reported by ChTire::GetSlipAngle.
GetSlipAngle_internal()107 double GetSlipAngle_internal() const { return m_slip->alpha; }
108
109 /// Get the tire longitudinal slip computed internally by the Pacejka model.
110 /// The reported value will be the same as that reported by ChTire::GetLongitudinalSlip.
GetLongitudinalSlip_internal()111 double GetLongitudinalSlip_internal() const { return m_slip->kappa; }
112
113 /// Get the camber angle for the Pacejka tire model (in radians).
114 /// The reported value will be the same as that reported by ChTire::GetCamberAngle.
GetCamberAngle_internal()115 double GetCamberAngle_internal() const { return m_slip->gamma; }
116
117 /// Write output data to a file.
118 void WriteOutData(double time, const std::string& outFilename);
119
120 /// Manually set the vertical wheel load as an input.
set_Fz_override(double Fz)121 void set_Fz_override(double Fz) { m_Fz_override = Fz; }
122
123 /// Return orientation, Vx (global) and omega/omega_y (global).
124 /// Assumes the tire is going straight forward (global x-dir), and the
125 /// returned state's orientation yields gamma and alpha, as x and z NASA angles
126 WheelState getState_from_KAG(double kappa, ///< [in] ...
127 double alpha, ///< [in] ...
128 double gamma, ///< [in] ...
129 double Vx ///< [in] tire forward velocity x-dir
130 );
131
132 /// Get the average simulation time per step spent in advance()
get_average_Advance_time()133 double get_average_Advance_time() { return m_sum_Advance_time / (double)m_num_Advance_calls; }
134
135 /// Get the average simulation time per step spent in calculating ODEs
get_average_ODE_time()136 double get_average_ODE_time() { return m_sum_ODE_time / (double)m_num_ODE_calls; }
137
138 /// Get current wheel longitudinal slip.
139 double get_kappa() const;
140
141 /// Get current wheel slip angle.
142 double get_alpha() const;
143
144 /// Get current wheel camber angle.
145 double get_gamma() const;
146
147 /// Get current long slip rate used in Magic Formula EQs.
148 double get_kappaPrime() const;
149
150 /// Get current slip angle using in Magic Formula.
151 double get_alphaPrime() const;
152
153 /// Get current camber angle used in Magic Formula.
154 double get_gammaPrime() const;
155
156 /// Get minimum longitudinal slip rate.
157 double get_min_long_slip() const;
158
159 /// Get maximum longitudinal slip rate.
160 double get_max_long_slip() const;
161
162 /// Get the minimum allowable lateral slip angle, alpha.
163 double get_min_lat_slip() const;
164
165 /// Get the maximum allowable lateral slip angle, alpha.
166 double get_max_lat_slip() const;
167
168 /// Get the longitudinal velocity.
169 double get_longvl() const;
170
171 /// Get the tire rolling radius, ideally updated each step.
get_tire_rolling_rad()172 double get_tire_rolling_rad() const { return m_R_eff; }
173
174 protected:
getPacTireParamFile()175 const std::string& getPacTireParamFile() const { return m_paramFile; }
176 std::string m_paramFile; // input parameter file
177
178 private:
Create(const rapidjson::Document & d)179 virtual void Create(const rapidjson::Document& d) override {}
180
181 /// Get the tire force and moment.
182 /// This represents the output from this tire system that is passed to the
183 /// vehicle system. Typically, the vehicle subsystem will pass the tire force
184 /// to the appropriate suspension subsystem which applies it as an external
185 /// force one the wheel body.
186 virtual TerrainForce GetTireForce() const override;
187
188 /// Initialize this tire by associating it to the specified wheel.
189 virtual void Initialize(std::shared_ptr<ChWheel> wheel) override;
190
191 /// Update the state of this tire system at the current time.
192 /// Set the PacTire spindle state data from the global wheel body state.
193 virtual void Synchronize(double time, ///< [in] current time
194 const ChTerrain& terrain ///< [in] reference to the terrain system
195 ) override;
196
197 /// Advance the state of this tire by the specified time step.
198 /// Use the new body state, calculate all the relevant quantities over the time increment.
199 virtual void Advance(double step) override;
200
201 // look for this data file
202 void loadPacTireParamFile();
203
204 // once Pac tire input text file has been successfully opened, read the input
205 // data, and populate the data struct
206 void readPacTireInput(std::ifstream& inFile);
207
208 // functions for reading each section in the parameter file
209 void readSection_UNITS(std::ifstream& inFile);
210 void readSection_MODEL(std::ifstream& inFile);
211 void readSection_DIMENSION(std::ifstream& inFile);
212 void readSection_SHAPE(std::ifstream& inFile);
213 void readSection_VERTICAL(std::ifstream& inFile);
214 void readSection_RANGES(std::ifstream& inFile);
215 void readSection_scaling(std::ifstream& inFile);
216 void readSection_longitudinal(std::ifstream& inFile);
217 void readSection_overturning(std::ifstream& inFile);
218 void readSection_lateral(std::ifstream& inFile);
219 void readSection_rolling(std::ifstream& inFile);
220 void readSection_aligning(std::ifstream& inFile);
221
222 /// update the tire contact coordinate system, TYDEX W-Axis
223 /// checks for contact, sets m_in_contact and m_depth
224 void update_W_frame(const ChTerrain& terrain);
225
226 // update the vertical load, tire deflection, and tire rolling radius
227 void update_verticalLoad(double step);
228
229 // find the vertical load, using a spring-damper model
230 double calc_Fz();
231
232 // calculate the various stiffness/relaxation lengths
233 void relaxationLengths();
234
235 // calculate the slips assume the steer/throttle/brake input to wheel has an
236 // instantaneous effect on contact patch slips
237 void slip_kinematic();
238
239 // set the tire m_slip vector to all zeros
240 void zero_slips();
241
242 // check if slips fall within a specified valid range
243 void evaluate_slips();
244
245 // after calculating all the reactions, evaluate output for any fishy business
246 // write_violations: any output threshold exceeded gets written to console
247 // enforce_threshold: enforce the limits when forces/moments exceed thresholds (except on Fz)
248 void evaluate_reactions(bool write_violations, bool enforce_threshold);
249
250 void advance_tire(double step);
251
252 // calculate transient slip properties, using first order ODEs to find slip
253 // displacements from velocities
254 // appends m_slips for the slip displacements, and integrated slip velocity terms
255 void advance_slip_transient(double step);
256
257 // calculate the increment delta_x using RK 45 integration for linear u, v_alpha
258 // Eq. 7.9 and 7.7, respectively
259 double ODE_RK_uv(double V_s, // slip velocity
260 double sigma, // either sigma_kappa or sigma_alpha
261 double V_cx, // wheel body center velocity
262 double step_size, // the simulation timestep size
263 double x_curr); // f(x_curr)
264
265 // calculate the linearized slope of Fy vs. alphaP at the current value
266 double get_dFy_dtan_alphaP(double x_curr);
267
268 // calculate v_alpha differently at low speeds
269 // relaxation length is non-linear
270 double ODE_RK_kappaAlpha(double V_sy, double V_cx, double C_Fy, double step_size, double x_curr);
271
272 // calculate the increment delta_gamma of the gamma ODE
273 // Note: gamma should be the kinematic camber, gamma
274 double ODE_RK_gamma(double C_Fgamma,
275 double C_Falpha,
276 double sigma_alpha,
277 double V_cx,
278 double step_size,
279 double gamma,
280 double v_gamma);
281
282 // calculate the ODE dphi/dt = f(phi), return the increment delta_phi
283 double ODE_RK_phi(double C_Fphi,
284 double C_Falpha,
285 double V_cx,
286 double psi_dot,
287 double omega,
288 double gamma,
289 double sigma_alpha,
290 double v_phi,
291 double eps_gamma,
292 double step_size);
293
294 // calculate and set the transient slip values (kappaP, alphaP, gammaP) from
295 // u, v deflections. optionally turn on/off Besselink low velocity damping
296 void slip_from_uv(bool use_besselink = true,
297 double bessel_Cx = 350.0,
298 double bessel_Cy = 200.0,
299 double V_low = 5.0);
300
301 /// calculate the reaction forces and moments, pure slip cases
302 /// assign longitudinal, lateral force, aligning moment:
303 /// Fx, Fy and Mz
304 void pureSlipReactions();
305
306 /// calculate combined slip reactions
307 /// assign Fx, Fy, Mz
308 void combinedSlipReactions();
309
310 /// longitudinal force, alpha ~= 0
311 /// assign to m_FM.force.x
312 /// assign m_pureLong, trigonometric function calculated constants
313 double Fx_pureLong(double gamma, double kappa);
314
315 /// lateral force, kappa ~= 0
316 /// assign to m_FM.force.y
317 /// assign m_pureLong, trigonometric function calculated constants
318 double Fy_pureLat(double alpha, double gamma);
319
320 /// aligning moment, kappa ~= 0
321 /// assign to m_FM.moment.z
322 /// assign m_pureLong, trigonometric function calculated constants
323 double Mz_pureLat(double alpha, double gamma, double Fy_pureSlip);
324
325 /// longitudinal force, combined slip (general case)
326 /// assign m_FM_combined.force.x
327 /// assign m_combinedLong
328 double Fx_combined(double alpha, double gamma, double kappa, double Fx_pureSlip);
329
330 /// calculate lateral force, combined slip (general case)
331 /// assign m_FM_combined.force.y
332 /// assign m_combinedLat
333 double Fy_combined(double alpha, double gamma, double kappa, double Fy_pureSlip);
334
335 // calculate aligning torque, combined slip (general case)
336 /// assign m_FM_combined.moment.z
337 /// assign m_combinedTorque
338 double Mz_combined(double alpha_r,
339 double alpha_t,
340 double gamma,
341 double kappa,
342 double Fx_combined,
343 double Fy_combined);
344
345 /// calculate the overturning couple moment
346 /// assign m_FM.moment.x and m_FM_combined.moment.x
347 double calc_Mx(double gamma, double Fy_combined);
348
349 /// calculate the rolling resistance moment,
350 /// assign m_FM.moment.y and m_FM_combined.moment.y
351 double calc_My(double Fx_combined);
352
353 // ----- Data members
354 bool m_use_transient_slip;
355 bool m_driven; // is this a driven tire?
356 int m_sameSide; // does parameter file side equal m_side? 1 = true, -1 opposite
357
358 WheelState m_tireState; // tire state, global coordinates
359 ChCoordsys<> m_W_frame; // tire contact coordinate system, TYDEX W-Axis
360 double m_simTime; // Chrono simulation time
361
362 bool m_in_contact; // indicates if there is tire-terrain contact
363 double m_depth; // if in contact, this is the contact depth
364
365 double m_R0; // unloaded radius
366 double m_R_eff; // current effect rolling radius
367 double m_R_l; // statically loaded radius
368
369 double m_Fz; // vertical load to use in Magic Formula
370 double m_dF_z; // (Fz - Fz,nom) / Fz,nom
371 bool m_use_Fz_override; // calculate Fz using collision, or user input
372 double m_Fz_override; // if manually inputting the vertical wheel load
373
374 double m_time_since_last_step; // init. to -1 in Initialize()
375 bool m_initial_step; // so Advance() gets called at time = 0
376 int m_num_ODE_calls;
377 double m_sum_ODE_time;
378 int m_num_Advance_calls;
379 double m_sum_Advance_time;
380
381 TerrainForce m_FM_pure; // output tire forces, based on pure slip
382 TerrainForce m_FM_combined; // output tire forces, based on combined slip
383 // previous steps calculated reaction
384 TerrainForce m_FM_pure_last;
385 TerrainForce m_FM_combined_last;
386
387 // TODO: could calculate these using sigma_kappa_adams and sigma_alpha_adams, in getRelaxationLengths()
388 // HARDCODED IN Initialize() for now
389 double m_C_Fx;
390 double m_C_Fy;
391
392 std::string m_outFilename; // output filename
393
394 int m_Num_WriteOutData; // number of times WriteOut was called
395
396 bool m_params_defined; // indicates if model params. have been defined/loaded
397
398 // MODEL PARAMETERS
399
400 // important slip quantities
401 slips* m_slip;
402
403 // model parameter factors stored here
404 Pac2002_data* m_params;
405
406 // for keeping track of intermediate factors in the PacTire model
407 pureLongCoefs* m_pureLong;
408 pureLatCoefs* m_pureLat;
409 pureTorqueCoefs* m_pureTorque;
410
411 combinedLongCoefs* m_combinedLong;
412 combinedLatCoefs* m_combinedLat;
413 combinedTorqueCoefs* m_combinedTorque;
414
415 zetaCoefs* m_zeta;
416
417 double m_mu; // current road friction coefficient
418 double m_mu0; // tire reference friction coeffient
419
420 ChFunction_Recorder m_areaDep; // lookup table for estimation of penetration depth from intersection area
421
422 // for transient contact point tire model
423 relaxationL* m_relaxation;
424 bessel* m_bessel;
425
426 std::shared_ptr<ChCylinderShape> m_cyl_shape; ///< visualization cylinder asset
427 std::shared_ptr<ChTexture> m_texture; ///< visualization texture asset
428 };
429
430 // -----------------------------------------------------------------------------
431 // String manipulation utility functions.
432 // -----------------------------------------------------------------------------
433
splitStr(const std::string & s,char delim,std::vector<std::string> & elems)434 inline std::vector<std::string>& splitStr(const std::string& s, char delim, std::vector<std::string>& elems) {
435 std::stringstream ss(s);
436 std::string item;
437 while (getline(ss, item, delim)) {
438 elems.push_back(item);
439 }
440 return elems;
441 }
442
splitStr(const std::string & s,char delim)443 inline std::vector<std::string> splitStr(const std::string& s, char delim) {
444 std::vector<std::string> elems;
445 return splitStr(s, delim, elems);
446 }
447
448 // convert a string to another type
449 template <class T>
fromString(const std::string & s)450 T fromString(const std::string& s) {
451 std::istringstream stream(s);
452 T t;
453 stream >> t;
454 return t;
455 }
456
457 // read an input line from tm.config, return the number that matters
458 template <class T>
fromTline(const std::string & tline)459 T fromTline(const std::string& tline) {
460 T t = fromString<T>(splitStr(splitStr(tline, '/')[0], '=')[1]);
461 return t;
462 }
463
464 /// @} vehicle_wheeled_tire
465
466 } // end namespace vehicle
467 } // end namespace chrono
468
469 #endif
470