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: Antonio Recuero, Bryan Peterson
13 // =============================================================================
14 //
15 // FEA deformable terrain
16 //
17 // =============================================================================
18
19 #include <cstdio>
20 #include <cmath>
21
22 #include "chrono/physics/ChMaterialSurfaceNSC.h"
23 #include "chrono/physics/ChMaterialSurfaceSMC.h"
24
25 #include "chrono/fea/ChElementHexaANCF_3813_9.h"
26 #include "chrono/fea/ChContactSurfaceMesh.h"
27 #include "chrono/fea/ChVisualizationFEAmesh.h"
28
29 #include "chrono_vehicle/ChVehicleModelData.h"
30 #include "chrono_vehicle/ChWorldFrame.h"
31 #include "chrono_vehicle/terrain/FEADeformableTerrain.h"
32
33 using namespace chrono::fea;
34
35 namespace chrono {
36 namespace vehicle {
37
38 // -----------------------------------------------------------------------------
39 // Implementation of the FEADeformableTerrain wrapper class
40 // -----------------------------------------------------------------------------
FEADeformableTerrain(ChSystem * system)41 FEADeformableTerrain::FEADeformableTerrain(ChSystem* system)
42 : m_E(1.379e7),
43 m_nu(0.3),
44 m_rho(200.0),
45 m_yield_stress(10000.0),
46 m_hardening_slope(5000),
47 m_friction_angle(0.00001),
48 m_dilatancy_angle(0.00001) {
49 m_mesh = chrono_types::make_shared<fea::ChMesh>();
50 system->Add(m_mesh);
51 }
52
53 // Return the terrain height at the specified location
GetHeight(const ChVector<> & loc) const54 double FEADeformableTerrain::GetHeight(const ChVector<>& loc) const {
55 //// TODO
56 return 0;
57 }
58
59 // Return the terrain normal at the specified location
GetNormal(const ChVector<> & loc) const60 ChVector<> FEADeformableTerrain::GetNormal(const ChVector<>& loc) const {
61 //// TODO
62 return ChWorldFrame::Vertical();
63 }
64
65 // Return the terrain coefficient of friction at the specified location
GetCoefficientFriction(const ChVector<> & loc) const66 float FEADeformableTerrain::GetCoefficientFriction(const ChVector<>& loc) const {
67 return m_friction_fun ? (*m_friction_fun)(loc) : 0.8f;
68 }
69
70 // Set properties of the FEA soil model
SetSoilParametersFEA(double rho,double Emod,double nu,double yield_stress,double hardening_slope,double friction_angle,double dilatancy_angle)71 void FEADeformableTerrain::SetSoilParametersFEA(double rho, ///< Soil density
72 double Emod, ///< Soil modulus of elasticity
73 double nu, ///< Soil Poisson ratio
74 double yield_stress, ///< Soil yield stress, for plasticity
75 double hardening_slope, ///< Soil hardening slope, for plasticity
76 double friction_angle, ///< Soil internal friction angle
77 double dilatancy_angle ///< Soil dilatancy angle
78 ) {
79 m_rho = rho;
80 m_E = Emod;
81 m_nu = nu;
82 m_yield_stress = yield_stress;
83 m_hardening_slope = hardening_slope;
84 m_friction_angle = friction_angle;
85 m_dilatancy_angle = dilatancy_angle;
86 }
87
88 // Initialize the terrain as a box of 9-node brick elements of given dimensions.
Initialize(const ChVector<> & start_point,const ChVector<> & terrain_dimension,const ChVector<int> & terrain_discretization)89 void FEADeformableTerrain::Initialize(const ChVector<>& start_point,
90 const ChVector<>& terrain_dimension,
91 const ChVector<int>& terrain_discretization) {
92 // Specification of the mesh (40,20,6)
93 int numDiv_x = terrain_discretization.x();
94 int numDiv_y = terrain_discretization.y();
95 int numDiv_z = terrain_discretization.z();
96
97 int N_x = numDiv_x + 1;
98 int N_y = numDiv_y + 1;
99
100 // Calculate total number of elements based on user input
101 int TotalNumElements = numDiv_x * numDiv_y * numDiv_z;
102 int XYNumNodes = (numDiv_x + 1) * (numDiv_y + 1);
103
104 // For uniform mesh
105 double dx = terrain_dimension.x() / numDiv_x;
106 double dy = terrain_dimension.y() / numDiv_y;
107 double dz = terrain_dimension.z() / numDiv_z;
108
109 bool Plasticity = true;
110
111 // Define location of nodes and create/add the nodes
112 for (int j = 0; j <= numDiv_z; j++) {
113 for (int i = 0; i < XYNumNodes; i++) {
114 // Node location
115 double loc_x = start_point.x() + (i % (numDiv_x + 1)) * dx;
116 double loc_y = start_point.y() + (i / (numDiv_x + 1)) % (numDiv_y + 1) * dy;
117 double loc_z = start_point.z() + j * dz;
118
119 // Create the node
120 auto node = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(loc_x, loc_y, loc_z));
121 node->SetMass(0);
122 // Fix all nodes along the axis Z = 0
123 if (j == 0) {
124 node->SetFixed(true);
125 }
126 // Add node to mesh
127 m_mesh->AddNode(node);
128 }
129 }
130
131 // Fix nodes at the boundaries of the FEA 'box'
132 for (int iz = 0; iz < numDiv_z; iz++) {
133 for (int ix = 0; ix < N_x; ix++) {
134 auto sidenode = std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode(iz * N_x * N_y + ix));
135 sidenode->SetFixed(true);
136
137 auto farsidenode =
138 std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode((N_x * numDiv_y) + iz * N_x * N_y + ix));
139 farsidenode->SetFixed(true);
140 }
141 }
142 for (int iz = 0; iz < numDiv_z; iz++) {
143 for (int iy = 0; iy < N_y; iy++) {
144 auto sidenode = std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode((iy + iz * N_y) * N_y));
145 sidenode->SetFixed(true);
146
147 auto farsidenode =
148 std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode(numDiv_x + iz * N_x * N_y + iy * N_y));
149 farsidenode->SetFixed(true);
150 }
151 }
152 // Initialize coordinates for curvature (central) node
153 for (int i = 0; i < TotalNumElements; i++) {
154 auto node = chrono_types::make_shared<ChNodeFEAcurv>(ChVector<>(0.0, 0.0, 0.0), ChVector<>(0.0, 0.0, 0.0),
155 ChVector<>(0.0, 0.0, 0.0));
156 node->SetMass(0);
157 m_mesh->AddNode(node);
158 }
159
160 // Basic material properties for soil.
161 auto material = chrono_types::make_shared<ChContinuumElastic>();
162 material->Set_RayleighDampingK(0.0);
163 material->Set_RayleighDampingM(0.0);
164 material->Set_density(m_rho);
165 material->Set_E(m_E);
166 material->Set_v(m_nu);
167
168 // Initial plastic deformation tensor: Initially identity (elastic).
169 ChMatrixNM<double, 9, 8> CCPInitial;
170 for (int k = 0; k < 8; k++) {
171 CCPInitial(0, k) = 1;
172 CCPInitial(4, k) = 1;
173 CCPInitial(8, k) = 1;
174 }
175 int jj = -1;
176 int kk = 0;
177 // Create the elements
178 for (int i = 0; i < TotalNumElements; i++) {
179 if (i % (numDiv_x * numDiv_y) == 0) {
180 jj++;
181 kk = 0;
182 }
183 // Define node sequence for element node0 through node7 are corner nodes
184 // Node8 is the central curvature vector node.
185 int node0 = (kk / (numDiv_x)) * (N_x) + kk % numDiv_x + jj * (N_x * N_y);
186 int node1 = (kk / (numDiv_x)) * (N_x) + kk % numDiv_x + 1 + jj * (N_x * N_y);
187 int node2 = (kk / (numDiv_x)) * (N_x) + kk % numDiv_x + 1 + N_x + jj * (N_x * N_y);
188 int node3 = (kk / (numDiv_x)) * (N_x) + kk % numDiv_x + N_x + jj * (N_x * N_y);
189 int node4 = (kk / (numDiv_x)) * (N_x) + kk % numDiv_x + XYNumNodes + jj * (N_x * N_y);
190 int node5 = (kk / (numDiv_x)) * (N_x) + kk % numDiv_x + 1 + XYNumNodes + jj * (N_x * N_y);
191 int node6 = (kk / (numDiv_x)) * (N_x) + kk % numDiv_x + 1 + N_x + XYNumNodes + jj * (N_x * N_y);
192 int node7 = (kk / (numDiv_x)) * (N_x) + kk % numDiv_x + N_x + XYNumNodes + jj * (N_x * N_y);
193 int node8 = (numDiv_z + 1) * XYNumNodes + i;
194
195 // Create the element and set its nodes.
196 auto element = chrono_types::make_shared<ChElementHexaANCF_3813_9>();
197 element->SetNodes(std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode(node0)),
198 std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode(node1)),
199 std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode(node2)),
200 std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode(node3)),
201 std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode(node4)),
202 std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode(node5)),
203 std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode(node6)),
204 std::dynamic_pointer_cast<ChNodeFEAxyz>(m_mesh->GetNode(node7)),
205 std::dynamic_pointer_cast<ChNodeFEAcurv>(m_mesh->GetNode(node8)));
206
207 // Set element dimensions
208 element->SetDimensions(ChVector<>(dx, dy, dz));
209
210 // Add a single layers with a fiber angle of 0 degrees.
211 element->SetMaterial(material);
212
213 // Set other element properties
214 element->SetAlphaDamp(5e-4); // Structural damping for this element
215 element->SetDPIterationNo(50); // Set maximum number of iterations for Drucker-Prager Newton-Raphson
216 element->SetDPYieldTol(1e-8); // Set stop tolerance for Drucker-Prager Newton-Raphson
217 element->SetStrainFormulation(ChElementHexaANCF_3813_9::Hencky);
218 element->SetPlasticityFormulation(ChElementHexaANCF_3813_9::DruckerPrager);
219 if (element->GetStrainFormulation() == ChElementHexaANCF_3813_9::Hencky) {
220 element->SetPlasticity(Plasticity);
221 if (Plasticity) {
222 element->SetYieldStress(m_yield_stress);
223 element->SetHardeningSlope(m_hardening_slope);
224 element->SetCCPInitial(CCPInitial);
225 if (element->GetPlasticityFormulation() == ChElementHexaANCF_3813_9::DruckerPrager) {
226 element->SetFriction(m_friction_angle);
227 element->SetDilatancy(m_dilatancy_angle);
228 element->SetDPType(3);
229 }
230 }
231 }
232
233 // Add element to mesh
234 m_mesh->AddElement(element);
235 kk++;
236 }
237
238 // -------------------------------------
239 // Options for visualization in irrlicht
240 // -------------------------------------
241
242 auto mvisualizemesh = chrono_types::make_shared<ChVisualizationFEAmesh>(*(m_mesh.get()));
243 mvisualizemesh->SetFEMdataType(ChVisualizationFEAmesh::E_PLOT_NODE_SPEED_NORM);
244 mvisualizemesh->SetColorscaleMinMax(0.0, 5.50);
245 mvisualizemesh->SetShrinkElements(true, 0.995);
246 mvisualizemesh->SetSmoothFaces(false);
247 m_mesh->AddAsset(mvisualizemesh);
248 }
249 } // end namespace vehicle
250 } // end namespace chrono
251