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 // Demo code illustrating the SCM semi-empirical model for deformable soil
13 // =============================================================================
14
15 #include "chrono/geometry/ChTriangleMeshConnected.h"
16 #include "chrono/physics/ChLinkMotorRotationAngle.h"
17 #include "chrono/physics/ChLoadContainer.h"
18 #include "chrono/physics/ChSystemSMC.h"
19 #include "chrono/utils/ChUtilsInputOutput.h"
20 #ifdef CHRONO_COLLISION
21 #include "chrono/collision/ChCollisionSystemChrono.h"
22 #endif
23
24 #include "chrono_irrlicht/ChIrrApp.h"
25
26 #include "chrono_vehicle/ChVehicleModelData.h"
27 #include "chrono_vehicle/terrain/SCMDeformableTerrain.h"
28
29 #include "chrono_thirdparty/filesystem/path.h"
30
31 using namespace chrono;
32 using namespace chrono::irrlicht;
33
34 using namespace irr;
35
36 bool output = false;
37 const std::string out_dir = GetChronoOutputPath() + "SCM_DEF_SOIL";
38
39 // Type of tire (controls both contact and visualization)
40 enum class TireType { CYLINDRICAL, LUGGED };
41 TireType tire_type = TireType::LUGGED;
42
43 // SCM grid spacing
44 double mesh_resolution = 0.04;
45
46 // Enable/disable bulldozing effects
47 bool enable_bulldozing = true;
48
49 // Enable/disable moving patch feature
50 bool enable_moving_patch = true;
51
52 // If true, use provided callback to change soil properties based on location
53 bool var_params = true;
54
55 // Custom callback for setting location-dependent soil properties.
56 // Note that the location is given in the SCM reference frame.
57 // Here, the vehicle moves in the terrain's negative y direction!
58 class MySoilParams : public vehicle::SCMDeformableTerrain::SoilParametersCallback {
59 public:
Set(const ChVector<> & loc,double & Bekker_Kphi,double & Bekker_Kc,double & Bekker_n,double & Mohr_cohesion,double & Mohr_friction,double & Janosi_shear,double & elastic_K,double & damping_R)60 virtual void Set(const ChVector<>& loc,
61 double& Bekker_Kphi,
62 double& Bekker_Kc,
63 double& Bekker_n,
64 double& Mohr_cohesion,
65 double& Mohr_friction,
66 double& Janosi_shear,
67 double& elastic_K,
68 double& damping_R) override {
69 if (loc.y() > 0) {
70 Bekker_Kphi = 0.2e6;
71 Bekker_Kc = 0;
72 Bekker_n = 1.1;
73 Mohr_cohesion = 0;
74 Mohr_friction = 30;
75 Janosi_shear = 0.01;
76 elastic_K = 4e7;
77 damping_R = 3e4;
78 } else {
79 Bekker_Kphi = 5301e3;
80 Bekker_Kc = 102e3;
81 Bekker_n = 0.793;
82 Mohr_cohesion = 1.3e3;
83 Mohr_friction = 31.1;
84 Janosi_shear = 1.2e-2;
85 elastic_K = 4e8;
86 damping_R = 3e4;
87 }
88 }
89 };
90
main(int argc,char * argv[])91 int main(int argc, char* argv[]) {
92 GetLog() << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << "\n\n";
93
94 // Set world frame with Y up
95 vehicle::ChWorldFrame::SetYUP();
96
97 // Global parameter for tire:
98 double tire_rad = 0.8;
99 ChVector<> tire_center(0, 0.02 + tire_rad, -1.5);
100
101 // Create a Chrono::Engine physical system
102 auto collsys_type = collision::ChCollisionSystemType::BULLET;
103 ChSystemSMC my_system;
104 my_system.SetNumThreads(4, 8, 1);
105 if (collsys_type == collision::ChCollisionSystemType::CHRONO) {
106 #ifdef CHRONO_COLLISION
107 auto collsys = chrono_types::make_shared<collision::ChCollisionSystemChrono>();
108 collsys->SetBroadphaseGridResolution(ChVector<int>(20, 20, 10));
109 my_system.SetCollisionSystem(collsys);
110 #endif
111 }
112
113 // Create the Irrlicht visualization (open the Irrlicht device,
114 // bind a simple user interface, etc. etc.)
115 ChIrrApp application(&my_system, L"Deformable soil", core::dimension2d<u32>(1280, 720), VerticalDir::Y, false, true);
116
117 // Easy shortcuts to add camera, lights, logo and sky in Irrlicht scene:
118 application.AddTypicalLogo();
119 application.AddTypicalSky();
120 application.AddTypicalLights();
121 application.AddTypicalCamera(core::vector3df(2.0f, 1.4f, 0.0f), core::vector3df(0, (f32)tire_rad, 0));
122 application.AddLightWithShadow(core::vector3df(1.5f, 5.5f, -2.5f), core::vector3df(0, 0, 0), 3, 2.2, 7.2, 40, 512,
123 video::SColorf(0.8f, 0.8f, 1.0f));
124
125 auto mtruss = chrono_types::make_shared<ChBody>(collsys_type);
126 mtruss->SetBodyFixed(true);
127 my_system.Add(mtruss);
128
129 // Initialize output
130 if (output) {
131 if (!filesystem::create_directory(filesystem::path(out_dir))) {
132 std::cout << "Error creating directory " << out_dir << std::endl;
133 return 1;
134 }
135 }
136 utils::CSV_writer csv(" ");
137
138 //
139 // Create a rigid body with a mesh or a cylinder collision shape
140 //
141
142 auto mrigidbody = chrono_types::make_shared<ChBody>(collsys_type);
143 my_system.Add(mrigidbody);
144 mrigidbody->SetMass(500);
145 mrigidbody->SetInertiaXX(ChVector<>(20, 20, 20));
146 mrigidbody->SetPos(tire_center + ChVector<>(0, 0.3, 0));
147
148 auto material = chrono_types::make_shared<ChMaterialSurfaceSMC>();
149 mrigidbody->GetCollisionModel()->ClearModel();
150 switch (tire_type) {
151 case TireType::LUGGED: {
152 auto trimesh = chrono_types::make_shared<geometry::ChTriangleMeshConnected>();
153 trimesh->LoadWavefrontMesh(GetChronoDataFile("models/tractor_wheel/tractor_wheel.obj"));
154
155 std::shared_ptr<ChTriangleMeshShape> mrigidmesh(new ChTriangleMeshShape);
156 mrigidmesh->SetMesh(trimesh);
157 mrigidbody->AddAsset(mrigidmesh);
158
159 mrigidbody->GetCollisionModel()->AddTriangleMesh(material, trimesh, false, false, VNULL, ChMatrix33<>(1),
160 0.01);
161 break;
162 }
163 case TireType::CYLINDRICAL: {
164 double radius = 0.5;
165 double width = 0.4;
166 mrigidbody->GetCollisionModel()->AddCylinder(material, radius, radius, width / 2, ChVector<>(0), Q_from_AngZ(CH_C_PI_2));
167
168 auto cyl_shape = chrono_types::make_shared<ChCylinderShape>();
169 cyl_shape->GetCylinderGeometry().rad = radius;
170 cyl_shape->GetCylinderGeometry().p1 = ChVector<>(+width / 2, 0, 0);
171 cyl_shape->GetCylinderGeometry().p2 = ChVector<>(-width / 2, 0, 0);
172 mrigidbody->AddAsset(cyl_shape);
173
174 break;
175 }
176 }
177 mrigidbody->GetCollisionModel()->BuildModel();
178 mrigidbody->SetCollide(true);
179
180 std::shared_ptr<ChColorAsset> mcol(new ChColorAsset);
181 mcol->SetColor(ChColor(0.3f, 0.3f, 0.3f));
182 mrigidbody->AddAsset(mcol);
183
184 auto motor = chrono_types::make_shared<ChLinkMotorRotationAngle>();
185 motor->SetSpindleConstraint(ChLinkMotorRotation::SpindleConstraint::OLDHAM);
186 motor->SetAngleFunction(chrono_types::make_shared<ChFunction_Ramp>(0, CH_C_PI / 4.0));
187 motor->Initialize(mrigidbody, mtruss, ChFrame<>(tire_center, Q_from_AngY(CH_C_PI_2)));
188 my_system.Add(motor);
189
190 //
191 // THE DEFORMABLE TERRAIN
192 //
193
194 // Create the 'deformable terrain' object
195 vehicle::SCMDeformableTerrain mterrain(&my_system);
196
197 // Displace/rotate the terrain reference plane.
198 // Note that SCMDeformableTerrain uses a default ISO reference frame (Z up). Since the mechanism is modeled here in
199 // a Y-up global frame, we rotate the terrain plane by -90 degrees about the X axis.
200 mterrain.SetPlane(ChCoordsys<>(ChVector<>(0, 0.2, 0), Q_from_AngX(-CH_C_PI_2)));
201
202 // Initialize the geometry of the soil
203
204 // Use either a regular grid:
205 double length = 6;
206 double width = 2;
207 mterrain.Initialize(width, length, mesh_resolution);
208
209 // Or use a height map:
210 ////mterrain.Initialize(vehicle::GetDataFile("terrain/height_maps/test64.bmp"), "test64", 1.6, 1.6, 0, 0.3);
211
212 // Set the soil terramechanical parameters
213 if (var_params) {
214 // Location-dependent soil properties
215 auto my_params = chrono_types::make_shared<MySoilParams>();
216 mterrain.RegisterSoilParametersCallback(my_params);
217 } else {
218 // Constant soil properties
219 mterrain.SetSoilParameters(0.2e6, // Bekker Kphi
220 0, // Bekker Kc
221 1.1, // Bekker n exponent
222 0, // Mohr cohesive limit (Pa)
223 30, // Mohr friction limit (degrees)
224 0.01, // Janosi shear coefficient (m)
225 4e7, // Elastic stiffness (Pa/m), before plastic yield, must be > Kphi
226 3e4 // Damping (Pa s/m), proportional to negative vertical speed (optional)
227 );
228
229 // LETE sand parameters
230 ////mterrain.SetSoilParameters(5301e3, // Bekker Kphi
231 //// 102e3, // Bekker Kc
232 //// 0.793, // Bekker n exponent
233 //// 1.3e3, // Mohr cohesive limit (Pa)
234 //// 31.1, // Mohr friction limit (degrees)
235 //// 1.2e-2, // Janosi shear coefficient (m)
236 //// 4e8, // Elastic stiffness (Pa/m), before plastic yield, must be > Kphi
237 //// 3e4 // Damping (Pa s/m), proportional to negative vertical speed (optional)
238 ////);
239 }
240
241 if (enable_bulldozing) {
242 mterrain.EnableBulldozing(true); // inflate soil at the border of the rut
243 mterrain.SetBulldozingParameters(
244 55, // angle of friction for erosion of displaced material at the border of the rut
245 1, // displaced material vs downward pressed material.
246 5, // number of erosion refinements per timestep
247 6); // number of concentric vertex selections subject to erosion
248 }
249
250 // Optionally, enable moving patch feature (reduces number of ray casts)
251 if (enable_moving_patch) {
252 mterrain.AddMovingPatch(mrigidbody, ChVector<>(0, 0, 0), ChVector<>(0.5, 2 * tire_rad, 2 * tire_rad));
253 }
254
255 // Set some visualization parameters: either with a texture, or with falsecolor plot, etc.
256 ////mterrain.SetTexture(vehicle::GetDataFile("terrain/textures/grass.jpg"), 16, 16);
257 mterrain.SetPlotType(vehicle::SCMDeformableTerrain::PLOT_PRESSURE, 0, 30000.2);
258 ////mterrain.SetPlotType(vehicle::SCMDeformableTerrain::PLOT_PRESSURE_YELD, 0, 30000.2);
259 ////mterrain.SetPlotType(vehicle::SCMDeformableTerrain::PLOT_SINKAGE, 0, 0.15);
260 ////mterrain.SetPlotType(vehicle::SCMDeformableTerrain::PLOT_SINKAGE_PLASTIC, 0, 0.15);
261 ////mterrain.SetPlotType(vehicle::SCMDeformableTerrain::PLOT_SINKAGE_ELASTIC, 0, 0.05);
262 ////mterrain.SetPlotType(vehicle::SCMDeformableTerrain::PLOT_STEP_PLASTIC_FLOW, 0, 0.0001);
263 ////mterrain.SetPlotType(vehicle::SCMDeformableTerrain::PLOT_ISLAND_ID, 0, 8);
264 ////mterrain.SetPlotType(vehicle::SCMDeformableTerrain::PLOT_IS_TOUCHED, 0, 8);
265 mterrain.SetMeshWireframe(true);
266
267 // ==IMPORTANT!== Use this function for adding a ChIrrNodeAsset to all items
268 application.AssetBindAll();
269
270 // ==IMPORTANT!== Use this function for 'converting' into Irrlicht meshes the assets
271 application.AssetUpdateAll();
272
273 // Use shadows in realtime view
274 application.AddShadowAll();
275
276 //
277 // THE SOFT-REAL-TIME CYCLE
278 //
279 /*
280 // Change the timestepper to HHT:
281 my_system.SetTimestepperType(ChTimestepper::Type::HHT);
282 auto integrator = std::static_pointer_cast<ChTimestepperHHT>(my_system.GetTimestepper());
283 integrator->SetAlpha(-0.2);
284 integrator->SetMaxiters(8);
285 integrator->SetAbsTolerances(1e-05, 1.8e00);
286 integrator->SetMode(ChTimestepperHHT::POSITION);
287 integrator->SetModifiedNewton(true);
288 integrator->SetScaling(true);
289 integrator->SetVerbose(true);
290 */
291 /*
292 my_system.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT);
293 */
294
295 application.SetTimestep(0.002);
296
297 while (application.GetDevice()->run()) {
298 double time = my_system.GetChTime();
299 if (output) {
300 vehicle::TerrainForce frc = mterrain.GetContactForce(mrigidbody);
301 csv << time << frc.force << frc.moment << frc.point << std::endl;
302 }
303
304 ////std::cout << "\nTime: " << time << std::endl;
305 ////std::cout << "Wheel pos: " << mrigidbody->GetPos() << std::endl;
306 ////std::cout << "Wheel rot: " << mrigidbody->GetRot() << std::endl;
307
308 application.BeginScene();
309 application.GetActiveCamera()->setTarget(core::vector3dfCH(mrigidbody->GetPos()));
310 application.DrawAll();
311 application.DoStep();
312 tools::drawColorbar(0, 30000, "Pressure yield [Pa]", application.GetDevice(), 1180);
313 application.EndScene();
314
315 ////mterrain.PrintStepStatistics(std::cout);
316 }
317
318 if (output) {
319 csv.write_to_file(out_dir + "/output.dat");
320 }
321
322 return 0;
323 }
324
