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: Alessandro Tasora, Radu Serban
13 // =============================================================================
14
15 #include "chrono/fea/ChContinuumMaterial.h"
16
17 namespace chrono {
18 namespace fea {
19
20 // -----------------------------------------------------------------------------
21
ChContinuumMaterial(const ChContinuumMaterial & other)22 ChContinuumMaterial::ChContinuumMaterial(const ChContinuumMaterial& other) {
23 density = other.density;
24 }
25
26 // Register into the object factory, to enable run-time dynamic creation and persistence
CH_FACTORY_REGISTER(ChContinuumMaterial)27 CH_FACTORY_REGISTER(ChContinuumMaterial)
28
29 void ChContinuumMaterial::ArchiveOUT(ChArchiveOut& marchive) {
30 // version number
31 marchive.VersionWrite<ChContinuumMaterial>();
32 // serialize parent class
33 // serialize all member data:
34 marchive << CHNVP(density, "density");
35 }
36
ArchiveIN(ChArchiveIn & marchive)37 void ChContinuumMaterial::ArchiveIN(ChArchiveIn& marchive) {
38 // version number
39 /*int version =*/ marchive.VersionRead<ChContinuumMaterial>();
40 // deserialize parent class
41 // stream in all member data:
42 marchive >> CHNVP(density, "density");
43 }
44
45 // -----------------------------------------------------------------------------
46
47 // Register into the object factory, to enable run-time dynamic creation and persistence
CH_FACTORY_REGISTER(ChContinuumElastic)48 CH_FACTORY_REGISTER(ChContinuumElastic)
49
50 ChContinuumElastic::ChContinuumElastic(double myoung, double mpoisson, double mdensity)
51 : ChContinuumMaterial(mdensity) {
52 E = myoung;
53 Set_v(mpoisson); // sets also G and l
54 ComputeStressStrainMatrix(); // sets Elasticity matrix
55 this->damping_M = 0;
56 this->damping_K = 0;
57 }
58
ChContinuumElastic(const ChContinuumElastic & other)59 ChContinuumElastic::ChContinuumElastic(const ChContinuumElastic& other) : ChContinuumMaterial(other) {
60 E = other.E;
61 v = other.v;
62 G = other.G;
63 l = other.l;
64 damping_M = other.damping_M;
65 damping_K = other.damping_K;
66
67 StressStrainMatrix = other.StressStrainMatrix;
68 }
69
Set_E(double m_E)70 void ChContinuumElastic::Set_E(double m_E) {
71 E = m_E;
72 G = E / (2 * (1 + v)); // fixed v, E, get G
73 l = (v * E) / ((1 + v) * (1 - 2 * v)); // Lame's constant l
74 ComputeStressStrainMatrix(); // updates Elasticity matrix
75 }
76
Set_v(double m_v)77 void ChContinuumElastic::Set_v(double m_v) {
78 v = m_v;
79 G = E / (2 * (1 + v)); // fixed v, E, get G
80 l = (v * E) / ((1 + v) * (1 - 2 * v)); // Lame's constant l
81 ComputeStressStrainMatrix(); // updates Elasticity matrix
82 }
83
Set_G(double m_G)84 void ChContinuumElastic::Set_G(double m_G) {
85 G = m_G;
86 v = (E / (2 * G)) - 1; // fixed G, E, get v
87 l = (v * E) / ((1 + v) * (1 - 2 * v)); // Lame's constant l
88 ComputeStressStrainMatrix(); // updates Elasticity matrix
89 }
90
ComputeStressStrainMatrix()91 void ChContinuumElastic::ComputeStressStrainMatrix() {
92 StressStrainMatrix.setZero(6, 6);
93 StressStrainMatrix(0, 0) = (E * (1 - v)) / (1 + v) / (1 - 2 * v);
94 // StressStrainMatrix(1,1)=StressStrainMatrix(0,0); //
95 // StressStrainMatrix(2,2)=StressStrainMatrix(0,0); //per non ricalcolare; qual'� meglio?
96 StressStrainMatrix(1, 1) = (E * (1 - v)) / (1 + v) / (1 - 2 * v);
97 StressStrainMatrix(2, 2) = (E * (1 - v)) / (1 + v) / (1 - 2 * v);
98 StressStrainMatrix(0, 1) = (E * (v)) / (1 + v) / (1 - 2 * v);
99 StressStrainMatrix(0, 2) = (E * (v)) / (1 + v) / (1 - 2 * v);
100 StressStrainMatrix(1, 0) = (E * (v)) / (1 + v) / (1 - 2 * v);
101 StressStrainMatrix(1, 2) = (E * (v)) / (1 + v) / (1 - 2 * v);
102 StressStrainMatrix(2, 0) = (E * (v)) / (1 + v) / (1 - 2 * v);
103 StressStrainMatrix(2, 1) = (E * (v)) / (1 + v) / (1 - 2 * v);
104 StressStrainMatrix(3, 3) = (E * (1 - 2 * v)) / (1 + v) / (1 - 2 * v) / 2;
105 StressStrainMatrix(4, 4) = (E * (1 - 2 * v)) / (1 + v) / (1 - 2 * v) / 2;
106 StressStrainMatrix(5, 5) = (E * (1 - 2 * v)) / (1 + v) / (1 - 2 * v) / 2;
107 }
108
ComputeElasticStress(ChStressTensor<> & mstress,const ChStrainTensor<> & mstrain) const109 void ChContinuumElastic::ComputeElasticStress(ChStressTensor<>& mstress, const ChStrainTensor<>& mstrain) const {
110 mstress.XX() = mstrain.XX() * (l + 2 * G) + mstrain.YY() * l + mstrain.ZZ() * l;
111 mstress.YY() = mstrain.XX() * l + mstrain.YY() * (l + 2 * G) + mstrain.ZZ() * l;
112 mstress.ZZ() = mstrain.XX() * l + mstrain.YY() * l + mstrain.ZZ() * (l + 2 * G);
113 mstress.XY() = mstrain.XY() * 2 * G;
114 mstress.XZ() = mstrain.XZ() * 2 * G;
115 mstress.YZ() = mstrain.YZ() * 2 * G;
116 }
117
ComputeElasticStrain(ChStrainTensor<> & mstrain,const ChStressTensor<> & mstress) const118 void ChContinuumElastic::ComputeElasticStrain(ChStrainTensor<>& mstrain, const ChStressTensor<>& mstress) const {
119 double invE = 1. / E;
120 double invhG = 0.5 / G;
121 mstrain.XX() = invE * (mstress.XX() - mstress.YY() * v - mstress.ZZ() * v);
122 mstrain.YY() = invE * (-mstress.XX() * v + mstress.YY() - mstress.ZZ() * v);
123 mstrain.ZZ() = invE * (-mstress.XX() * v - mstress.YY() * v + mstress.ZZ());
124 mstrain.XY() = mstress.XY() * invhG;
125 mstrain.XZ() = mstress.XZ() * invhG;
126 mstrain.YZ() = mstress.YZ() * invhG;
127 }
128
ArchiveOUT(ChArchiveOut & marchive)129 void ChContinuumElastic::ArchiveOUT(ChArchiveOut& marchive) {
130 // version number
131 marchive.VersionWrite<ChContinuumElastic>();
132 // serialize parent class
133 ChContinuumMaterial::ArchiveOUT(marchive);
134 // serialize all member data:
135 marchive << CHNVP(this->E);
136 marchive << CHNVP(this->v);
137 marchive << CHNVP(this->damping_M);
138 marchive << CHNVP(this->damping_K);
139 }
140
ArchiveIN(ChArchiveIn & marchive)141 void ChContinuumElastic::ArchiveIN(ChArchiveIn& marchive) {
142 // version number
143 /*int version =*/ marchive.VersionRead<ChContinuumElastic>();
144 // deserialize parent class
145 ChContinuumMaterial::ArchiveIN(marchive);
146 // stream in all member data:
147 marchive >> CHNVP(this->E);
148 marchive >> CHNVP(this->v);
149 this->Set_v(this->v); // G and l from v
150 marchive >> CHNVP(this->damping_M);
151 marchive >> CHNVP(this->damping_K);
152 }
153
154 // -----------------------------------------------------------------------------
155
ArchiveOUT(ChArchiveOut & marchive)156 void ChContinuumElastoplastic::ArchiveOUT(ChArchiveOut& marchive) {
157 // version number
158 marchive.VersionWrite<ChContinuumElastoplastic>();
159 // serialize parent class
160 ChContinuumElastic::ArchiveOUT(marchive);
161 // serialize all member data:
162 }
163
ArchiveIN(ChArchiveIn & marchive)164 void ChContinuumElastoplastic::ArchiveIN(ChArchiveIn& marchive) {
165 // version number
166 /*int version =*/ marchive.VersionRead<ChContinuumElastoplastic>();
167 // deserialize parent class
168 ChContinuumElastic::ArchiveIN(marchive);
169 // stream in all member data:
170 }
171
172 // -----------------------------------------------------------------------------
173
174 // Register into the object factory, to enable run-time dynamic creation and persistence
CH_FACTORY_REGISTER(ChContinuumPlasticVonMises)175 CH_FACTORY_REGISTER(ChContinuumPlasticVonMises)
176
177 ChContinuumPlasticVonMises::ChContinuumPlasticVonMises(double myoung,
178 double mpoisson,
179 double mdensity,
180 double melastic_yeld,
181 double mplastic_yeld)
182 : ChContinuumElastoplastic(myoung, mpoisson, mdensity),
183 elastic_yeld(melastic_yeld),
184 plastic_yeld(mplastic_yeld),
185 flow_rate(1) {}
186
ChContinuumPlasticVonMises(const ChContinuumPlasticVonMises & other)187 ChContinuumPlasticVonMises::ChContinuumPlasticVonMises(const ChContinuumPlasticVonMises& other)
188 : ChContinuumElastoplastic(other) {
189 elastic_yeld = other.elastic_yeld;
190 plastic_yeld = other.plastic_yeld;
191 flow_rate = other.flow_rate;
192 }
193
ComputeYeldFunction(const ChStressTensor<> & mstress) const194 double ChContinuumPlasticVonMises::ComputeYeldFunction(const ChStressTensor<>& mstress) const {
195 return (mstress.GetEquivalentVonMises() - this->elastic_yeld);
196 }
197
ComputeReturnMapping(ChStrainTensor<> & mplasticstrainflow,const ChStrainTensor<> & mincrementstrain,const ChStrainTensor<> & mlastelasticstrain,const ChStrainTensor<> & mlastplasticstrain) const198 void ChContinuumPlasticVonMises::ComputeReturnMapping(ChStrainTensor<>& mplasticstrainflow,
199 const ChStrainTensor<>& mincrementstrain,
200 const ChStrainTensor<>& mlastelasticstrain,
201 const ChStrainTensor<>& mlastplasticstrain) const {
202 ChStrainTensor<> guesselstrain(mlastelasticstrain);
203 guesselstrain += mincrementstrain; // assume increment is all elastic
204
205 double vonm = guesselstrain.GetEquivalentVonMises();
206 if (vonm > this->elastic_yeld) {
207 ChVoightTensor<> mdev;
208 guesselstrain.GetDeviatoricPart(mdev);
209 mplasticstrainflow = mdev * ((vonm - this->elastic_yeld) / (vonm));
210 } else {
211 mplasticstrainflow.setZero();
212 }
213 }
214
ComputePlasticStrainFlow(ChStrainTensor<> & mplasticstrainflow,const ChStrainTensor<> & mtotstrain) const215 void ChContinuumPlasticVonMises::ComputePlasticStrainFlow(ChStrainTensor<>& mplasticstrainflow,
216 const ChStrainTensor<>& mtotstrain) const {
217 double vonm = mtotstrain.GetEquivalentVonMises();
218 if (vonm > this->elastic_yeld) {
219 ChVoightTensor<> mdev;
220 mtotstrain.GetDeviatoricPart(mdev);
221 mplasticstrainflow = mdev * ((vonm - this->elastic_yeld) / (vonm));
222 } else {
223 mplasticstrainflow.setZero();
224 }
225 }
226
ArchiveOUT(ChArchiveOut & marchive)227 void ChContinuumPlasticVonMises::ArchiveOUT(ChArchiveOut& marchive) {
228 // version number
229 marchive.VersionWrite<ChContinuumPlasticVonMises>();
230 // serialize parent class
231 ChContinuumElastoplastic::ArchiveOUT(marchive);
232 // serialize all member data:
233 marchive << CHNVP(this->elastic_yeld);
234 marchive << CHNVP(this->plastic_yeld);
235 marchive << CHNVP(this->flow_rate);
236 }
237
ArchiveIN(ChArchiveIn & marchive)238 void ChContinuumPlasticVonMises::ArchiveIN(ChArchiveIn& marchive) {
239 // version number
240 /*int version =*/ marchive.VersionRead<ChContinuumPlasticVonMises>();
241 // deserialize parent class
242 ChContinuumElastoplastic::ArchiveIN(marchive);
243 // stream in all member data:
244 marchive >> CHNVP(this->elastic_yeld);
245 marchive >> CHNVP(this->plastic_yeld);
246 marchive >> CHNVP(this->flow_rate);
247 }
248
249 // -----------------------------------------------------------------------------
250
251 // Register into the object factory, to enable run-time dynamic creation and persistence
CH_FACTORY_REGISTER(ChContinuumDruckerPrager)252 CH_FACTORY_REGISTER(ChContinuumDruckerPrager)
253
254 ChContinuumDruckerPrager::ChContinuumDruckerPrager(double myoung,
255 double mpoisson,
256 double mdensity,
257 double melastic_yeld,
258 double malpha,
259 double mdilatancy)
260 : ChContinuumElastoplastic(myoung, mpoisson, mdensity),
261 elastic_yeld(melastic_yeld),
262 alpha(malpha),
263 dilatancy(mdilatancy),
264 hardening_limit(elastic_yeld),
265 hardening_speed(0),
266 flow_rate(1) {}
267
ChContinuumDruckerPrager(const ChContinuumDruckerPrager & other)268 ChContinuumDruckerPrager::ChContinuumDruckerPrager(const ChContinuumDruckerPrager& other)
269 : ChContinuumElastoplastic(other) {
270 elastic_yeld = other.elastic_yeld;
271 alpha = other.alpha;
272 dilatancy = other.dilatancy;
273 hardening_speed = other.hardening_speed;
274 hardening_limit = other.hardening_limit;
275 flow_rate = other.flow_rate;
276 }
277
Set_from_MohrCoulomb(double phi,double cohesion,bool inner_approx)278 void ChContinuumDruckerPrager::Set_from_MohrCoulomb(double phi, double cohesion, bool inner_approx) {
279 if (inner_approx) {
280 alpha = (2 * sin(phi)) / (sqrt(3.0) * (3.0 - sin(phi)));
281 elastic_yeld = (6 * cohesion * cos(phi)) / (sqrt(3.0) * (3.0 - sin(phi)));
282 } else {
283 alpha = (2 * sin(phi)) / (sqrt(3.0) * (3.0 + sin(phi)));
284 elastic_yeld = (6 * cohesion * cos(phi)) / (sqrt(3.0) * (3.0 + sin(phi)));
285 }
286 }
287
ComputeYeldFunction(const ChStressTensor<> & mstress) const288 double ChContinuumDruckerPrager::ComputeYeldFunction(const ChStressTensor<>& mstress) const {
289 return (mstress.GetInvariant_I1() * this->alpha + sqrt(mstress.GetInvariant_J2()) - this->elastic_yeld);
290 }
291
ComputeReturnMapping(ChStrainTensor<> & mplasticstrainflow,const ChStrainTensor<> & mincrementstrain,const ChStrainTensor<> & mlastelasticstrain,const ChStrainTensor<> & mlastplasticstrain) const292 void ChContinuumDruckerPrager::ComputeReturnMapping(ChStrainTensor<>& mplasticstrainflow,
293 const ChStrainTensor<>& mincrementstrain,
294 const ChStrainTensor<>& mlastelasticstrain,
295 const ChStrainTensor<>& mlastplasticstrain) const {
296 ChStrainTensor<> guesselstrain(mlastelasticstrain);
297 guesselstrain += mincrementstrain; // assume increment is all elastic
298
299 ChStressTensor<> mstress;
300 this->ComputeElasticStress(mstress, guesselstrain);
301 double fprager = this->ComputeYeldFunction(mstress);
302
303 if (fprager > 0) {
304 if (mstress.GetInvariant_I1() * this->alpha - sqrt(mstress.GetInvariant_J2()) * this->alpha * this->alpha -
305 this->elastic_yeld >
306 0) {
307 // Case: tentative stress is in polar cone; a singular region where the gradient of
308 // the yield function (or flow potential) is not defined. Just project to vertex.
309 ChStressTensor<> vertexstress;
310 double vertcoord = this->elastic_yeld / (3 * this->alpha);
311 vertexstress.XX() = vertcoord;
312 vertexstress.YY() = vertcoord;
313 vertexstress.ZZ() = vertcoord;
314 ChStrainTensor<> vertexstrain;
315 this->ComputeElasticStrain(vertexstrain, vertexstress);
316 mplasticstrainflow = guesselstrain - vertexstrain;
317 } else {
318 // Case: tentative stress is out of the yield cone.
319 // Just project using the yield (or flow potential) gradient.
320 ChStrainTensor<> dFdS;
321 ChStrainTensor<> dGdS;
322 double devsq = sqrt(mstress.GetInvariant_J2());
323 if (devsq > 10e-16) {
324 double sixdevsq = 6 * devsq;
325
326 dFdS.XX() = this->alpha + (2 * mstress.XX() - mstress.YY() - mstress.ZZ()) / sixdevsq;
327 dFdS.YY() = this->alpha + (-mstress.XX() + 2 * mstress.YY() - mstress.ZZ()) / sixdevsq;
328 dFdS.ZZ() = this->alpha + (-mstress.XX() - mstress.YY() + 2 * mstress.ZZ()) / sixdevsq;
329 dFdS.XY() = mstress.XY() / devsq;
330 dFdS.YZ() = mstress.YZ() / devsq;
331 dFdS.XZ() = mstress.XZ() / devsq;
332
333 dGdS.XX() = this->dilatancy + (2 * mstress.XX() - mstress.YY() - mstress.ZZ()) / sixdevsq;
334 dGdS.YY() = this->dilatancy + (-mstress.XX() + 2 * mstress.YY() - mstress.ZZ()) / sixdevsq;
335 dGdS.ZZ() = this->dilatancy + (-mstress.XX() - mstress.YY() + 2 * mstress.ZZ()) / sixdevsq;
336 dGdS.XY() = mstress.XY() / devsq;
337 dGdS.YZ() = mstress.YZ() / devsq;
338 dGdS.XZ() = mstress.XZ() / devsq;
339 } else {
340 GetLog() << " ... axial singularity - SHOULD NEVER OCCUR - handled by polar cone\n";
341 dFdS.setZero();
342 dFdS.XX() = 1;
343 dFdS.YY() = 1;
344 dFdS.ZZ() = 1;
345 dGdS.setZero();
346 dGdS.XX() = 1;
347 dGdS.YY() = 1;
348 dGdS.ZZ() = 1;
349 }
350 ChStressTensor<> aux_dFdS_C;
351 this->ComputeElasticStress(aux_dFdS_C, dFdS);
352
353 ChMatrixNM<double, 1, 1> inner_up;
354 inner_up = aux_dFdS_C.transpose() * mincrementstrain;
355 ChMatrixNM<double, 1, 1> inner_dw;
356 inner_dw = aux_dFdS_C.transpose() * dGdS;
357
358 mplasticstrainflow = dGdS;
359 mplasticstrainflow *= inner_up(0) / inner_dw(0);
360 }
361 } else {
362 mplasticstrainflow.setZero();
363 }
364 }
365
366 //***OBSOLETE***
ComputePlasticStrainFlow(ChStrainTensor<> & mplasticstrainflow,const ChStrainTensor<> & mestrain) const367 void ChContinuumDruckerPrager::ComputePlasticStrainFlow(ChStrainTensor<>& mplasticstrainflow,
368 const ChStrainTensor<>& mestrain) const {
369 ChStressTensor<> mstress;
370 this->ComputeElasticStress(mstress, mestrain);
371 double prager = mstress.GetInvariant_I1() * this->alpha + sqrt(mstress.GetInvariant_J2());
372 if (prager > this->elastic_yeld) {
373 ChVoightTensor<> mdev;
374 mstress.GetDeviatoricPart(mdev);
375 double divisor = 2. * sqrt(mstress.GetInvariant_J2());
376 if (divisor > 10e-20)
377 mdev *= 1. / divisor;
378 mdev.XX() += this->dilatancy;
379 mdev.YY() += this->dilatancy;
380 mdev.ZZ() += this->dilatancy;
381 mplasticstrainflow = mdev;
382 } else {
383 mplasticstrainflow.setZero();
384 }
385 }
386
ArchiveOUT(ChArchiveOut & marchive)387 void ChContinuumDruckerPrager::ArchiveOUT(ChArchiveOut& marchive) {
388 // version number
389 marchive.VersionWrite<ChContinuumDruckerPrager>();
390 // serialize parent class
391 ChContinuumElastoplastic::ArchiveOUT(marchive);
392 // serialize all member data:
393 marchive << CHNVP(this->elastic_yeld);
394 marchive << CHNVP(this->alpha);
395 marchive << CHNVP(this->dilatancy);
396 marchive << CHNVP(this->hardening_speed);
397 marchive << CHNVP(this->hardening_limit);
398 marchive << CHNVP(this->flow_rate);
399 }
400
ArchiveIN(ChArchiveIn & marchive)401 void ChContinuumDruckerPrager::ArchiveIN(ChArchiveIn& marchive) {
402 // version number
403 /*int version =*/ marchive.VersionRead<ChContinuumDruckerPrager>();
404 // deserialize parent class
405 ChContinuumElastoplastic::ArchiveIN(marchive);
406 // stream in all member data:
407 marchive >> CHNVP(this->elastic_yeld);
408 marchive >> CHNVP(this->alpha);
409 marchive >> CHNVP(this->dilatancy);
410 marchive >> CHNVP(this->hardening_speed);
411 marchive >> CHNVP(this->hardening_limit);
412 marchive >> CHNVP(this->flow_rate);
413 }
414
415 } // end namespace fea
416 } // end namespace chrono
417