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: Radu Serban
13 // =============================================================================
14 //
15 // Base class for a tire.
16 // A tire subsystem is a force element. It is passed position and velocity
17 // information of the wheel body and it produces ground reaction forces and
18 // moments to be applied to the wheel body.
19 //
20 // =============================================================================
21
22 #include <cmath>
23
24 #include "chrono/physics/ChSystem.h"
25
26 #include "chrono_vehicle/ChVehicleModelData.h"
27 #include "chrono_vehicle/ChWorldFrame.h"
28 #include "chrono_vehicle/wheeled_vehicle/ChTire.h"
29
30 #include "chrono_thirdparty/filesystem/path.h"
31
32 namespace chrono {
33 namespace vehicle {
34
ChTire(const std::string & name)35 ChTire::ChTire(const std::string& name)
36 : ChPart(name),
37 m_collision_type(CollisionType::SINGLE_POINT),
38 m_stepsize(1e-3),
39 m_slip_angle(0),
40 m_longitudinal_slip(0),
41 m_camber_angle(0) {}
42
43 // -----------------------------------------------------------------------------
44 // Initialize this tire by associating it to the specified wheel.
45 // Increment the mass and inertia of the associated suspension spindle body to
46 // account for the tire mass and inertia.
47 // -----------------------------------------------------------------------------
Initialize(std::shared_ptr<ChWheel> wheel)48 void ChTire::Initialize(std::shared_ptr<ChWheel> wheel) {
49 m_wheel = wheel;
50
51 //// RADU
52 //// Todo: Properly account for offset in adjusting inertia.
53 //// This requires changing the spindle to a ChBodyAuxRef.
54 wheel->GetSpindle()->SetMass(wheel->GetSpindle()->GetMass() + GetMass());
55 wheel->GetSpindle()->SetInertiaXX(wheel->GetSpindle()->GetInertiaXX() + GetInertia());
56 }
57
58 // -----------------------------------------------------------------------------
59 // Default implementation of ReportMass simply returns the value from GetMass.
60 // However, concrete tire models which need to return 0 from GetMass (so that
61 // the mass of the tire is not double counted) will override this function.
62 // -----------------------------------------------------------------------------
ReportMass() const63 double ChTire::ReportMass() const {
64 return GetMass();
65 }
66
67 // -----------------------------------------------------------------------------
68 // Calculate kinematics quantities (slip angle, longitudinal slip, camber angle,
69 // and toe-in angle) using the given state of the associated wheel.
70 // -----------------------------------------------------------------------------
CalculateKinematics(double time,const WheelState & state,const ChTerrain & terrain)71 void ChTire::CalculateKinematics(double time, const WheelState& state, const ChTerrain& terrain) {
72 // Wheel normal (expressed in global frame)
73 ChVector<> wheel_normal = state.rot.GetYaxis();
74
75 // Terrain normal at wheel location (expressed in global frame)
76 ChVector<> Z_dir = terrain.GetNormal(state.pos);
77
78 // Longitudinal (heading) and lateral directions, in the terrain plane
79 ChVector<> X_dir = Vcross(wheel_normal, Z_dir);
80 X_dir.Normalize();
81 ChVector<> Y_dir = Vcross(Z_dir, X_dir);
82
83 // Tire reference coordinate system
84 ChMatrix33<> rot;
85 rot.Set_A_axis(X_dir, Y_dir, Z_dir);
86 ChCoordsys<> tire_csys(state.pos, rot.Get_A_quaternion());
87
88 // Express wheel linear velocity in tire frame
89 ChVector<> V = tire_csys.TransformDirectionParentToLocal(state.lin_vel);
90 // Express wheel normal in tire frame
91 ChVector<> n = tire_csys.TransformDirectionParentToLocal(wheel_normal);
92
93 // Slip angle (positive sign = left turn, negative sign = right turn)
94 double abs_Vx = std::abs(V.x());
95 double zero_Vx = 1e-4;
96 m_slip_angle = (abs_Vx > zero_Vx) ? std::atan(V.y() / abs_Vx) : 0;
97
98 // Longitudinal slip (positive sign = driving, negative sign = breaking)
99 m_longitudinal_slip = (abs_Vx > zero_Vx) ? -(V.x() - state.omega * GetRadius()) / abs_Vx : 0;
100
101 // Camber angle (positive sign = upper side tipping to the left, negative sign = upper side tipping to the right)
102 m_camber_angle = std::atan2(n.z(), n.y());
103 }
104
105 // -----------------------------------------------------------------------------
106 // Utility functions for adding/removing a mesh visualization asset to this tire
107 // -----------------------------------------------------------------------------
108
109 // Add visualization mesh: use one of the two provided OBJ files, depending on the side on which the tire is mounted.
AddVisualizationMesh(const std::string & mesh_file_left,const std::string & mesh_file_right)110 std::shared_ptr<ChTriangleMeshShape> ChTire::AddVisualizationMesh(const std::string& mesh_file_left,
111 const std::string& mesh_file_right) {
112 bool left = (m_wheel->GetSide() == VehicleSide::LEFT);
113 ChQuaternion<> rot = left ? Q_from_AngZ(0) : Q_from_AngZ(CH_C_PI);
114 m_vis_mesh_file = left ? mesh_file_left : mesh_file_right;
115
116 auto trimesh = chrono_types::make_shared<geometry::ChTriangleMeshConnected>();
117 trimesh->LoadWavefrontMesh(vehicle::GetDataFile(m_vis_mesh_file), false, false);
118 trimesh->Transform(ChVector<>(0, GetOffset(), 0), ChMatrix33<>(rot));
119
120 auto trimesh_shape = chrono_types::make_shared<ChTriangleMeshShape>();
121 trimesh_shape->Pos = ChVector<>(0, GetOffset(), 0);
122 trimesh_shape->Rot = ChMatrix33<>(rot);
123 trimesh_shape->SetMesh(trimesh);
124 trimesh_shape->SetName(filesystem::path(m_vis_mesh_file).stem());
125 trimesh_shape->SetStatic(true);
126
127 m_wheel->GetSpindle()->AddAsset(trimesh_shape);
128
129 return trimesh_shape;
130 }
131
132 // Remove visualization mesh shape. Make sure to only remove the specified asset. A associated body (a wheel spindle)
133 // may have additional assets.
RemoveVisualizationMesh(std::shared_ptr<ChTriangleMeshShape> trimesh_shape)134 void ChTire::RemoveVisualizationMesh(std::shared_ptr<ChTriangleMeshShape> trimesh_shape) {
135 auto& assets = m_wheel->GetSpindle()->GetAssets();
136 auto it = std::find(assets.begin(), assets.end(), trimesh_shape);
137 if (it != assets.end())
138 assets.erase(it);
139 }
140
141 // -----------------------------------------------------------------------------
142 // Utility functions for characterizing the geometric contact between a disc with
143 // specified center location, normal direction, and radius and the terrain,
144 // assumed to be specified as a height field (over the x-y domain).
145 // These functions returns false if no contact occurs.
146 // Otherwise, they set the contact points on the disc (ptD) and on the terrain (ptT),
147 // the normal contact direction, and the resulting penetration depth (a positive value).
148 //
149 // The first version uses a single point on the terrain (below the wheel center).
150 // The second version uses the average of four terrain heights below the wheel center.
151 // The third version uses the collision algorithm of Sui and Hirshey.
152 //
153 // NOTE: uses terrain normal at disc center for approximative calculation.
154 // Hence only valid for terrains with constant slope. A completely accurate
155 // solution would require an iterative calculation of the contact point.
156 // -----------------------------------------------------------------------------
DiscTerrainCollision(const ChTerrain & terrain,const ChVector<> & disc_center,const ChVector<> & disc_normal,double disc_radius,ChCoordsys<> & contact,double & depth)157 bool ChTire::DiscTerrainCollision(
158 const ChTerrain& terrain, // [in] reference to terrain system
159 const ChVector<>& disc_center, // [in] global location of the disc center
160 const ChVector<>& disc_normal, // [in] disc normal, expressed in the global frame
161 double disc_radius, // [in] disc radius
162 ChCoordsys<>& contact, // [out] contact coordinate system (relative to the global frame)
163 double& depth) // [out] penetration depth (positive if contact occurred)
164 {
165 // Find terrain height below disc center. There is no contact if the disc
166 // center is below the terrain or farther away by more than its radius.
167 double hc = terrain.GetHeight(disc_center);
168 double disc_height = ChWorldFrame::Height(disc_center);
169 if (disc_height <= hc || disc_height >= hc + disc_radius)
170 return false;
171
172 // Find the lowest point on the disc. There is no contact if the disc is (almost) horizontal.
173 ChVector<> nhelp = terrain.GetNormal(disc_center);
174 ChVector<> dir1 = Vcross(disc_normal, nhelp);
175 double sinTilt2 = dir1.Length2();
176
177 if (sinTilt2 < 1e-3)
178 return false;
179
180 // Contact point (lowest point on disc).
181 ChVector<> ptD = disc_center + disc_radius * Vcross(disc_normal, dir1 / sqrt(sinTilt2));
182
183 // Find terrain height at lowest point. No contact if lowest point is above the terrain.
184 double hp = terrain.GetHeight(ptD);
185 double ptD_height = ChWorldFrame::Height(ptD);
186 if (ptD_height > hp)
187 return false;
188
189 // Approximate the terrain with a plane. Define the projection of the lowest
190 // point onto this plane as the contact point on the terrain.
191 ChVector<> normal = terrain.GetNormal(ptD);
192 ChVector<> longitudinal = Vcross(disc_normal, normal);
193 longitudinal.Normalize();
194 ChVector<> lateral = Vcross(normal, longitudinal);
195 ChMatrix33<> rot;
196 rot.Set_A_axis(longitudinal, lateral, normal);
197
198 contact.pos = ptD;
199 contact.rot = rot.Get_A_quaternion();
200
201 depth = (hp - ptD_height) * (ChWorldFrame::Vertical() ^ normal);
202 assert(depth > 0);
203
204 return true;
205 }
206
DiscTerrainCollision4pt(const ChTerrain & terrain,const ChVector<> & disc_center,const ChVector<> & disc_normal,double disc_radius,double width,ChCoordsys<> & contact,double & depth,double & camber_angle)207 bool ChTire::DiscTerrainCollision4pt(
208 const ChTerrain& terrain, // [in] reference to terrain system
209 const ChVector<>& disc_center, // [in] global location of the disc center
210 const ChVector<>& disc_normal, // [in] disc normal, expressed in the global frame
211 double disc_radius, // [in] disc radius
212 double width, // [in] tire width
213 ChCoordsys<>& contact, // [out] contact coordinate system (relative to the global frame)
214 double& depth, // [out] penetration depth (positive if contact occurred),
215 double& camber_angle) // [out] camber angle
216 {
217 double dx = 0.1 * disc_radius;
218 double dy = 0.3 * width;
219
220 // Find terrain height below disc center. There is no contact if the disc
221 // center is below the terrain or farther away by more than its radius.
222 double hc = terrain.GetHeight(disc_center);
223 double disc_height = ChWorldFrame::Height(disc_center);
224 if (disc_height <= hc || disc_height >= hc + disc_radius)
225 return false;
226
227 // Find the lowest point on the disc. There is no contact if the disc is (almost) horizontal.
228 ChVector<> nhelp = terrain.GetNormal(disc_center);
229 ChVector<> dir1 = Vcross(disc_normal, nhelp);
230 double sinTilt2 = dir1.Length2();
231
232 if (sinTilt2 < 1e-3)
233 return false;
234
235 // Contact point (lowest point on disc).
236 ChVector<> ptD = disc_center + disc_radius * Vcross(disc_normal, dir1 / sqrt(sinTilt2));
237
238 // Approximate the terrain with a plane. Define the projection of the lowest
239 // point onto this plane as the contact point on the terrain.
240 ChVector<> normal = terrain.GetNormal(ptD);
241 ChVector<> longitudinal = Vcross(disc_normal, normal);
242 longitudinal.Normalize();
243 ChVector<> lateral = Vcross(normal, longitudinal);
244
245 // Calculate four contact points in the contact patch
246 ChVector<> ptQ1 = ptD + dx * longitudinal;
247 double hQ1 = terrain.GetHeight(ptQ1);
248 double ptQ1_height = ChWorldFrame::Height(ptQ1);
249 ptQ1 = ptQ1 - (ptQ1_height - hQ1) * ChWorldFrame::Vertical();
250
251 ChVector<> ptQ2 = ptD - dx * longitudinal;
252 double hQ2 = terrain.GetHeight(ptQ2);
253 double ptQ2_height = ChWorldFrame::Height(ptQ2);
254 ptQ2 = ptQ2 - (ptQ2_height - hQ2) * ChWorldFrame::Vertical();
255
256 ChVector<> ptQ3 = ptD + dy * lateral;
257 double hQ3 = terrain.GetHeight(ptQ3);
258 double ptQ3_height = ChWorldFrame::Height(ptQ3);
259 ptQ3 = ptQ3 - (ptQ3_height - hQ3) * ChWorldFrame::Vertical();
260
261 ChVector<> ptQ4 = ptD - dy * lateral;
262 double hQ4 = terrain.GetHeight(ptQ4);
263 double ptQ4_height = ChWorldFrame::Height(ptQ4);
264 ptQ4 = ptQ4 - (ptQ4_height - hQ4) * ChWorldFrame::Vertical();
265
266 // Calculate a smoothed road surface normal
267 ChVector<> rQ2Q1 = ptQ1 - ptQ2;
268 ChVector<> rQ4Q3 = ptQ3 - ptQ4;
269
270 ChVector<> terrain_normal = Vcross(rQ2Q1, rQ4Q3);
271 terrain_normal.Normalize();
272
273 // Find terrain height as average of four points. No contact if lowest point is above the terrain.
274 ptD = 0.25 * (ptQ1 + ptQ2 + ptQ3 + ptQ4);
275 ChVector<> d = ptD - disc_center;
276 double da = d.Length();
277
278 if (da >= disc_radius)
279 return false;
280
281 // Calculate an improved value for the camber angle
282 camber_angle = std::asin(Vdot(disc_normal, terrain_normal));
283
284 ChMatrix33<> rot;
285 rot.Set_A_axis(longitudinal, lateral, terrain_normal);
286
287 contact.pos = ptD;
288 contact.rot = rot.Get_A_quaternion();
289
290 depth = disc_radius - da;
291 assert(depth > 0);
292
293 return true;
294 }
295
ConstructAreaDepthTable(double disc_radius,ChFunction_Recorder & areaDep)296 void ChTire::ConstructAreaDepthTable(double disc_radius, ChFunction_Recorder& areaDep) {
297 const size_t n_lookup = 90;
298 double depMax = disc_radius; // should be high enough to avoid extrapolation
299 double depStep = depMax / double(n_lookup - 1);
300 for (size_t i = 0; i < n_lookup; i++) {
301 double dep = depStep * double(i);
302 double alpha = 2.0 * acos(1.0 - dep / disc_radius);
303 double area = 0.5 * disc_radius * disc_radius * (alpha - sin(alpha));
304 areaDep.AddPoint(area, dep);
305 }
306 }
307
DiscTerrainCollisionEnvelope(const ChTerrain & terrain,const ChVector<> & disc_center,const ChVector<> & disc_normal,double disc_radius,const ChFunction_Recorder & areaDep,ChCoordsys<> & contact,double & depth)308 bool ChTire::DiscTerrainCollisionEnvelope(
309 const ChTerrain& terrain, // [in] reference to terrain system
310 const ChVector<>& disc_center, // [in] global location of the disc center
311 const ChVector<>& disc_normal, // [in] disc normal, expressed in the global frame
312 double disc_radius, // [in] disc radius
313 const ChFunction_Recorder& areaDep, // [in] lookup table to calculate depth from intersection area
314 ChCoordsys<>& contact, // [out] contact coordinate system (relative to the global frame)
315 double& depth // [out] penetration depth (positive if contact occurred)
316 ) {
317 // Approximate the terrain with a plane. Define the projection of the lowest
318 // point onto this plane as the contact point on the terrain. We don't know
319 // where the equivalent contact point is exactly, so we use the intersection
320 // area to decide if there is contact or not.
321
322 ChVector<> normal = terrain.GetNormal(disc_center);
323 ChVector<> longitudinal = Vcross(disc_normal, normal);
324 longitudinal.Normalize();
325
326 const size_t n_div = 180;
327 double x_step = 2.0 * disc_radius / n_div;
328 double A = 0; // overlapping area of tire disc and road surface contour
329 for (size_t i = 1; i < n_div; i++) {
330 double x = -disc_radius + x_step * double(i);
331 ChVector<> pTest = disc_center + x * longitudinal;
332 double q = terrain.GetHeight(pTest);
333 double a = ChWorldFrame::Height(pTest) - sqrt(disc_radius * disc_radius - x * x);
334 if (q > a) {
335 A += q - a;
336 }
337 }
338 A *= x_step;
339
340 if (A == 0) {
341 return false;
342 }
343
344 // calculate equivalent depth from A
345 depth = areaDep.Get_y(A);
346
347 // Find the lowest point on the disc. There is no contact if the disc is (almost) horizontal.
348 ChVector<> nhelp = terrain.GetNormal(disc_center);
349 ChVector<> dir1 = Vcross(disc_normal, nhelp);
350 double sinTilt2 = dir1.Length2();
351
352 if (sinTilt2 < 1e-3)
353 return false;
354
355 // Contact point (lowest point on disc).
356 ChVector<> ptD = disc_center + (disc_radius - depth) * Vcross(disc_normal, dir1 / sqrt(sinTilt2));
357
358 // Find terrain height at lowest point. No contact if lowest point is above
359 // the terrain.
360
361 normal = terrain.GetNormal(ptD);
362 longitudinal = Vcross(disc_normal, normal);
363 longitudinal.Normalize();
364 ChVector<> lateral = Vcross(normal, longitudinal);
365 ChMatrix33<> rot;
366 rot.Set_A_axis(longitudinal, lateral, normal);
367
368 contact.pos = ptD;
369 contact.rot = rot.Get_A_quaternion();
370
371 return true;
372 }
373
374 // -----------------------------------------------------------------------------
375 // Estimate the tire moments of inertia, given the tire specification and mass.
376 // The tire is assumed to be composed of simple geometric shapes:
377 // - The sidewalls are treated as discs with an inner diameter equal to the
378 // wheel diameter, and an outer diameter equal to the static diameter of the
379 // tire.
380 // - The tread face is treated as a band with a diameter equal to the static
381 // radius of the tire.
382 // A rubber material is assumed, using a density of 1050 kg/m^3.
383 // -----------------------------------------------------------------------------
VolumeCyl(double r_outer,double r_inner,double cyl_height)384 double VolumeCyl(double r_outer, double r_inner, double cyl_height) {
385 return CH_C_PI * cyl_height * (r_outer * r_outer - r_inner * r_inner);
386 }
387
InertiaRotCyl(double mass,double r_outer,double r_inner)388 double InertiaRotCyl(double mass, double r_outer, double r_inner) {
389 return mass * (r_outer * r_outer + r_inner * r_inner) / 2.0;
390 }
391
InertiaTipCyl(double mass,double r_outer,double r_inner,double cyl_height)392 double InertiaTipCyl(double mass, double r_outer, double r_inner, double cyl_height) {
393 return mass * (3.0 * (r_outer * r_outer + r_inner * r_inner) + cyl_height * cyl_height) / 12.0;
394 }
395
EstimateInertia(double tire_width,double aspect_ratio,double rim_diameter,double tire_mass,double t_factor)396 ChVector<> ChTire::EstimateInertia(double tire_width, // tire width [mm]
397 double aspect_ratio, // aspect ratio: height to width [percentage]
398 double rim_diameter, // rim diameter [in]
399 double tire_mass, // mass of the tire [kg]
400 double t_factor // tread to sidewall thickness factor
401 ) {
402 double rho = 1050.0; // rubber density in kg/m^3
403
404 double width = tire_width / 1000; // tire width in meters
405 double secth = (aspect_ratio / 100) * width; // tire height in meters
406 double r_rim = (rim_diameter / 2) * 0.0254; // rim radius in meters
407 double r_tire = r_rim + secth; // tire radius in meters
408
409 // Estimate sidewall thickness.
410 double t = 0;
411 double m_test = 0;
412 while (tire_mass - m_test > 0.0) {
413 t += 1e-6;
414 m_test = rho * (VolumeCyl(r_tire, r_tire - t_factor * t, width - 2 * t) + 2 * VolumeCyl(r_tire - t, r_rim, t));
415 }
416
417 // Lateral sidewall offset.
418 double r_steiner = (width - t) / 2.0;
419
420 // Calculate mass and moments of inertia of tread section and sidewall, separately.
421 double m_tread = rho * VolumeCyl(r_tire, r_tire - t_factor * t, width - 2 * t);
422 double Irot_tread = InertiaRotCyl(m_tread, r_tire, r_tire - 2 * t);
423 double Itip_tread = InertiaTipCyl(m_tread, r_tire, r_tire - t_factor * t, width - 2 * t);
424
425 double m_sidewall = rho * VolumeCyl(r_tire - t, r_rim, t);
426 double Irot_sidewall = InertiaRotCyl(m_sidewall, r_tire - t, r_rim);
427 double Itip_sidewall = InertiaTipCyl(m_sidewall, r_tire - t, r_rim, t) + m_sidewall * pow(r_steiner, 2.0);
428
429 // Return composite tire inertia.
430 return ChVector<>(Itip_tread + 2 * Itip_sidewall, Irot_tread + 2 * Irot_sidewall, Itip_tread + 2 * Itip_sidewall);
431 }
432
433 } // end namespace vehicle
434 } // end namespace chrono
435