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 #ifndef CHCONTINUUMMATERIAL_H 16 #define CHCONTINUUMMATERIAL_H 17 18 #include "chrono/core/ChApiCE.h" 19 #include "chrono/core/ChMath.h" 20 #include "chrono/core/ChTensors.h" 21 22 namespace chrono { 23 namespace fea { 24 25 /// @addtogroup chrono_fea 26 /// @{ 27 28 // ----------------------------------------------------------------------------- 29 30 /// Base class for properties of materials in a continuum. 31 32 class ChApi ChContinuumMaterial { 33 protected: 34 double density; 35 36 public: density(mdensity)37 ChContinuumMaterial(double mdensity = 1000) : density(mdensity) {} 38 ChContinuumMaterial(const ChContinuumMaterial& other); ~ChContinuumMaterial()39 virtual ~ChContinuumMaterial() {} 40 41 /// Set the density of the material, in kg/m^2. Set_density(double m_density)42 void Set_density(double m_density) { density = m_density; } 43 /// Get the density of the material, in kg/m^2. Get_density()44 double Get_density() const { return density; } 45 46 virtual void ArchiveOUT(ChArchiveOut& marchive); 47 virtual void ArchiveIN(ChArchiveIn& marchive); 48 }; 49 50 /// Class for the basic properties of materials in an elastic continuum. 51 /// This is a base material with isotropic hookean elasticity. 52 53 class ChApi ChContinuumElastic : public ChContinuumMaterial { 54 private: 55 double E; ///< Young Modulus 56 double v; ///< Poisson ratio 57 double G; ///< shear modulus 58 double l; ///< Lame's modulus 59 ChMatrixDynamic<> StressStrainMatrix; ///< Elasticity (stiffness) matrix : = [E] ε 60 61 double damping_M; ///< Raleigh_damping, M proportional 62 double damping_K; ///< Raleigh_damping, K proportional 63 64 public: 65 /// Create a continuum isotropic hookean material. 66 /// Default value for Young elastic modulus is low (like a 67 /// rubber-type material), and same for density. 68 ChContinuumElastic(double myoung = 10000000, double mpoisson = 0.4, double mdensity = 1000); 69 ChContinuumElastic(const ChContinuumElastic& other); ~ChContinuumElastic()70 virtual ~ChContinuumElastic() {} 71 72 /// Set the Young E elastic modulus, in Pa (N/m^2), as the ratio of the uniaxial 73 /// stress over the uniaxial strain, for hookean materials. Intuitively, the 74 /// tensile pressure on a side of a parallelepiped in order to double its length. 75 /// Note that most metal materials require very high values, ex. steel 76 /// has E=210GPa (E=210e9), aluminium E=69e9, and this can cause numerical 77 /// problems if you do not set up the simulation integrator properly. 78 void Set_E(double m_E); 79 /// Get the Young E elastic modulus, in Pa (N/m^2). Get_E()80 double Get_E() const { return E; } 81 82 /// Set the Poisson v ratio, as v=-transverse_strain/axial_strain, so 83 /// takes into account the 'squeezing' effect of materials that are pulled (so, 84 /// if zero, when you push the two sizes of a cube, it won't inflate). Most 85 /// materials have some 0<v<0.5, for example steel has v=0.27..0.30, aluminium v=0.33, 86 /// rubber=0.49, etc. Note! v=0.5 means perfectly incompressible material, that 87 /// could give problems with some type of solvers. 88 /// Setting v also changes G. 89 void Set_v(double m_v); 90 /// Get the Young v ratio, as v=-transverse_strain/axial_strain. Get_v()91 double Get_v() const { return v; } 92 93 /// Set the shear modulus G, in Pa (N/m^2), as the ratio of shear stress to 94 /// the shear strain. Setting G also changes Poisson ratio v. 95 void Set_G(double m_G); 96 /// Get the shear modulus G, in Pa (N/m^2) Get_G()97 double Get_G() const { return G; } 98 99 /// Get Lamé first parameter (the second is shear modulus, so Get_G() ) Get_l()100 double Get_l() const { return l; } 101 102 /// Get bulk modulus (increase of pressure for decrease of volume), in Pa. Get_BulkModulus()103 double Get_BulkModulus() const { return E / (3 * (1 - 2 * v)); } 104 105 /// Get P-wave modulus (if V=speed of propagation of a P-wave, then (M/density)=V^2 ) Get_WaveModulus()106 double Get_WaveModulus() const { return E * ((1 - v) / (1 + v) * (1 - 2 * v)); } 107 108 /// Computes Elasticity matrix and stores the value in this->StressStrainMatrix 109 /// Note: is performed every time you change a material parameter 110 void ComputeStressStrainMatrix(); 111 /// Get the Elasticity matrix Get_StressStrainMatrix()112 ChMatrixDynamic<>& Get_StressStrainMatrix() { return StressStrainMatrix; } 113 114 /// Compute elastic stress from elastic strain 115 /// (using column tensors, in Voight notation) 116 void ComputeElasticStress(ChStressTensor<>& mstress, const ChStrainTensor<>& mstrain) const; 117 118 /// Compute elastic strain from elastic stress 119 /// (using column tensors, in Voight notation) 120 void ComputeElasticStrain(ChStrainTensor<>& mstrain, const ChStressTensor<>& mstress) const; 121 122 /// Set the Rayleigh mass-proportional damping factor alpha, to 123 /// build damping R as R=alpha*M + beta*K Set_RayleighDampingM(double m_d)124 void Set_RayleighDampingM(double m_d) { damping_M = m_d; } 125 /// Set the Rayleigh mass-proportional damping factor alpha, in R=alpha*M + beta*K Get_RayleighDampingM()126 double Get_RayleighDampingM() const { return damping_M; } 127 128 /// Set the Rayleigh stiffness-proportional damping factor beta, to 129 /// build damping R as R=alpha*M + beta*K Set_RayleighDampingK(double m_d)130 void Set_RayleighDampingK(double m_d) { damping_K = m_d; } 131 /// Set the Rayleigh stiffness-proportional damping factor beta, in R=alpha*M + beta*K Get_RayleighDampingK()132 double Get_RayleighDampingK() const { return damping_K; } 133 134 virtual void ArchiveOUT(ChArchiveOut& marchive) override; 135 virtual void ArchiveIN(ChArchiveIn& marchive) override; 136 }; 137 138 // ----------------------------------------------------------------------------- 139 140 /// Class for all elastic materials that can undergo plastic flow. 141 /// Defines simply some interface functions. 142 143 class ChApi ChContinuumElastoplastic : public ChContinuumElastic { 144 public: 145 ChContinuumElastoplastic(double myoung = 10000000, double mpoisson = 0.4, double mdensity = 1000) ChContinuumElastic(myoung,mpoisson,mdensity)146 : ChContinuumElastic(myoung, mpoisson, mdensity) {} 147 148 /// Return a scalar value that is 0 on the yield surface, <0 inside (elastic), >0 outside (incompatible->plastic 149 /// flow) 150 virtual double ComputeYeldFunction(const ChStressTensor<>& mstress) const = 0; 151 152 /// Compute plastic strain flow (flow derivative dE_plast/dt) from strain, 153 /// according to VonMises strain yield theory. 154 virtual void ComputePlasticStrainFlow(ChStrainTensor<>& mplasticstrainflow, 155 const ChStrainTensor<>& mtotstrain) const = 0; 156 157 /// Correct the strain-stress by enforcing that elastic stress must remain on the yield 158 /// surface, computing a plastic flow to be added to plastic strain while integrating. 159 virtual void ComputeReturnMapping(ChStrainTensor<>& mplasticstrainflow, 160 const ChStrainTensor<>& mincrementstrain, 161 const ChStrainTensor<>& mlastelasticstrain, 162 const ChStrainTensor<>& mlastplasticstrain) const = 0; 163 164 /// Set the plastic flow rate, i.e. the 'creep' speed. The lower the value, the slower 165 /// the plastic flow during dynamic simulations, with delayed plasticity. 166 virtual void Set_flow_rate(double mflow_rate) = 0; 167 /// Set the plastic flow rate. 168 virtual double Get_flow_rate() const = 0; 169 170 virtual void ArchiveOUT(ChArchiveOut& marchive) override; 171 virtual void ArchiveIN(ChArchiveIn& marchive) override; 172 }; 173 174 // ----------------------------------------------------------------------------- 175 176 /// Class for the basic properties of materials in an elastoplastic continuum, 177 /// with strain yield limit based on Von Mises yield. 178 179 class ChApi ChContinuumPlasticVonMises : public ChContinuumElastoplastic { 180 private: 181 double elastic_yeld; 182 double plastic_yeld; 183 184 double flow_rate; 185 186 public: 187 /// Create a continuum isotropic elastoplastic material, 188 /// where you can define also plastic and elastic max. stress (yield limits 189 /// for transition elastic->plastic and plastic->fracture). 190 ChContinuumPlasticVonMises(double myoung = 10000000, 191 double mpoisson = 0.4, 192 double mdensity = 1000, 193 double melastic_yeld = 0.1, 194 double mplastic_yeld = 0.2); 195 ChContinuumPlasticVonMises(const ChContinuumPlasticVonMises& other); ~ChContinuumPlasticVonMises()196 virtual ~ChContinuumPlasticVonMises() {} 197 198 /// Set the elastic yield modulus as the maximum VonMises 199 /// equivalent strain that can be withstood by material before 200 /// starting plastic flow. It defines the transition elastic->plastic. Set_elastic_yeld(double melastic_yeld)201 void Set_elastic_yeld(double melastic_yeld) { elastic_yeld = melastic_yeld; }; 202 /// Get the elastic yield modulus. Get_elastic_yeld()203 double Get_elastic_yeld() const { return elastic_yeld; } 204 205 /// Set the plastic yield modulus as the maximum VonMises 206 /// equivalent strain that can be withstood by material before 207 /// fracture. It defines the transition plastic->fracture. Set_plastic_yeld(double mplastic_yeld)208 void Set_plastic_yeld(double mplastic_yeld) { plastic_yeld = mplastic_yeld; }; 209 /// Get the plastic yield modulus. Get_plastic_yeld()210 double Get_plastic_yeld() const { return plastic_yeld; } 211 212 /// Set the plastic flow rate. The lower the value, the slower 213 /// the plastic flow during dynamic simulations. Set_flow_rate(double mflow_rate)214 virtual void Set_flow_rate(double mflow_rate) override { flow_rate = mflow_rate; }; 215 /// Set the plastic flow rate. Get_flow_rate()216 virtual double Get_flow_rate() const override { return flow_rate; } 217 218 virtual double ComputeYeldFunction(const ChStressTensor<>& mstress) const override; 219 220 virtual void ComputeReturnMapping(ChStrainTensor<>& mplasticstrainflow, 221 const ChStrainTensor<>& mincrementstrain, 222 const ChStrainTensor<>& mlastelasticstrain, 223 const ChStrainTensor<>& mlastplasticstrain) const override; 224 225 /// Compute plastic strain flow (flow derivative dE_plast/dt) from strain, 226 /// according to VonMises strain yield theory. 227 virtual void ComputePlasticStrainFlow(ChStrainTensor<>& mplasticstrainflow, 228 const ChStrainTensor<>& mestrain) const override; 229 230 virtual void ArchiveOUT(ChArchiveOut& marchive) override; 231 virtual void ArchiveIN(ChArchiveIn& marchive) override; 232 }; 233 234 // ----------------------------------------------------------------------------- 235 236 /// Class for the basic properties of elastoplastic materials of Drucker-Prager type, 237 /// that are useful for simulating soils. 238 239 class ChApi ChContinuumDruckerPrager : public ChContinuumElastoplastic { 240 private: 241 double elastic_yeld; 242 double alpha; 243 double dilatancy; 244 double hardening_speed; 245 double hardening_limit; 246 double flow_rate; 247 248 public: 249 /// Create a continuum isotropic Drucker-Prager material 250 ChContinuumDruckerPrager(double myoung = 10000000, 251 double mpoisson = 0.4, 252 double mdensity = 1000, 253 double melastic_yeld = 0.1, 254 double malpha = 0.5, 255 double mdilatancy = 0); 256 ChContinuumDruckerPrager(const ChContinuumDruckerPrager& other); ~ChContinuumDruckerPrager()257 virtual ~ChContinuumDruckerPrager() {} 258 259 /// Set the D-P yield modulus C, for Drucker-Prager 260 /// yield. It defines the transition elastic->plastic. Set_elastic_yeld(double melastic_yeld)261 void Set_elastic_yeld(double melastic_yeld) { elastic_yeld = melastic_yeld; } 262 /// Get the elastic yield modulus C Get_elastic_yeld()263 double Get_elastic_yeld() const { return elastic_yeld; } 264 265 /// Set the internal friction coefficient A Set_alpha(double malpha)266 void Set_alpha(double malpha) { alpha = malpha; } 267 /// Get the internal friction coefficient A Get_alpha()268 double Get_alpha() const { return alpha; } 269 270 /// Sets the C and A parameters of the Drucker-Prager model 271 /// starting from more 'practical' values of inner friction angle phi 272 /// and cohesion, as used in the faceted Mohr-Coulomb model. 273 /// Use the optional parameter inner_approx to set if the faceted 274 /// Mohr-Coulomb must be approximated with D-P inscribed (default) or circumscribed. 275 void Set_from_MohrCoulomb(double phi, double cohesion, bool inner_approx = true); 276 277 /// Set the plastic flow rate multiplier. The lower the value, the slower 278 /// the plastic flow during dynamic simulations. Set_flow_rate(double mflow_rate)279 virtual void Set_flow_rate(double mflow_rate) override { flow_rate = mflow_rate; } 280 /// Get the flow rate multiplier. Get_flow_rate()281 virtual double Get_flow_rate() const override { return flow_rate; } 282 283 /// Set the internal dilatation coefficient (usually 0.. < int.friction) Set_dilatancy(double mdilatancy)284 void Set_dilatancy(double mdilatancy) { dilatancy = mdilatancy; } 285 /// Get the internal dilatation coefficient Get_dilatancy()286 double Get_dilatancy() const { return dilatancy; } 287 288 /// Set the hardening limit (usually a bit larger than yield), or softening Set_hardening_limit(double mhl)289 void Set_hardening_limit(double mhl) { hardening_limit = mhl; } 290 /// Get the hardening limit Get_hardening_limit()291 double Get_hardening_limit() const { return hardening_limit; } 292 293 /// Set the hardening inverse speed coeff. for exponential hardening 294 /// (the larger, the slower the hardening or softening process that 295 /// will asymptotically make yield = hardening_limit ) Set_hardening_speed(double mhl)296 void Set_hardening_speed(double mhl) { hardening_speed = mhl; } 297 /// Get the hardening speed Get_hardening_speed()298 double Get_hardening_speed() const { return hardening_speed; } 299 300 virtual double ComputeYeldFunction(const ChStressTensor<>& mstress) const override; 301 302 virtual void ComputeReturnMapping(ChStrainTensor<>& mplasticstrainflow, 303 const ChStrainTensor<>& mincrementstrain, 304 const ChStrainTensor<>& mlastelasticstrain, 305 const ChStrainTensor<>& mlastplasticstrain) const override; 306 307 /// Compute plastic strain flow direction from strain 308 /// according to Drucker-Prager. 309 virtual void ComputePlasticStrainFlow(ChStrainTensor<>& mplasticstrainflow, 310 const ChStrainTensor<>& mestrain) const override; 311 312 virtual void ArchiveOUT(ChArchiveOut& marchive) override; 313 virtual void ArchiveIN(ChArchiveIn& marchive) override; 314 }; 315 316 /// @} chrono_fea 317 318 } // end namespace fea 319 320 CH_CLASS_VERSION(fea::ChContinuumMaterial, 0) 321 CH_CLASS_VERSION(fea::ChContinuumElastic, 0) 322 CH_CLASS_VERSION(fea::ChContinuumElastoplastic, 0) 323 CH_CLASS_VERSION(fea::ChContinuumPlasticVonMises, 0) 324 CH_CLASS_VERSION(fea::ChContinuumDruckerPrager, 0) 325 326 } // end namespace chrono 327 328 #endif 329