KratosMultiphysics
KRATOS Multiphysics (Kratos) is a framework for building parallel, multi-disciplinary simulation software, aiming at modularity, extensibility, and high performance. Kratos is written in C++, and counts with an extensive Python interface.
modified_mohr_coulomb_yield_surface.h
Go to the documentation of this file.
1 // KRATOS ___| | | |
2 // \___ \ __| __| | | __| __| | | __| _` | |
3 // | | | | | ( | | | | ( | |
4 // _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS
5 //
6 // License: BSD License
7 // license: structural_mechanics_application/license.txt
8 //
9 // Main authors: Alejandro Cornejo & Lucia Barbu
10 //
11 
12 #pragma once
13 
14 // System includes
15 
16 // Project includes
17 #include "includes/checks.h"
18 #include "generic_yield_surface.h"
19 
20 namespace Kratos
21 {
24 
28 
29  // The size type definition
30  typedef std::size_t SizeType;
31 
35 
39 
43 
56 template <class TPlasticPotentialType>
58 {
59  public:
62 
64  typedef TPlasticPotentialType PlasticPotentialType;
65 
67  static constexpr SizeType Dimension = PlasticPotentialType::Dimension;
68 
70  static constexpr SizeType VoigtSize = PlasticPotentialType::VoigtSize;
71 
74 
76  static constexpr double tolerance = std::numeric_limits<double>::epsilon();
77 
81 
84  {
85  }
86 
89  {
90  }
91 
94  {
95  return *this;
96  }
97 
100 
107 
116  const array_1d<double, VoigtSize>& rPredictiveStressVector,
117  const Vector& rStrainVector,
118  double& rEquivalentStress,
120  )
121  {
122  const Properties& r_material_properties = rValues.GetMaterialProperties();
123 
124  const bool has_symmetric_yield_stress = r_material_properties.Has(YIELD_STRESS);
125  const double yield_compression = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_COMPRESSION];
126  const double yield_tension = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_TENSION];
127  double friction_angle = r_material_properties[FRICTION_ANGLE] * Globals::Pi / 180.0; // In radians!
128 
129  // Check input variables
130  if (friction_angle < tolerance) {
131  friction_angle = 32.0 * Globals::Pi / 180.0;
132  KRATOS_WARNING("ModifiedMohrCoulombYieldSurface") << "Friction Angle not defined, assumed equal to 32 deg " << std::endl;
133  }
134 
135  double theta;
136  const double R = std::abs(yield_compression / yield_tension);
137  const double Rmorh = std::pow(std::tan((Globals::Pi / 4.0) + friction_angle / 2.0), 2);
138  const double alpha_r = R / Rmorh;
139  const double sin_phi = std::sin(friction_angle);
140 
141  double I1, J2, J3;
144  AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateJ2Invariant(rPredictiveStressVector, I1, deviator, J2);
146 
147  const double K1 = 0.5 * (1.0 + alpha_r) - 0.5 * (1.0 - alpha_r) * sin_phi;
148  const double K2 = 0.5 * (1.0 + alpha_r) - 0.5 * (1.0 - alpha_r) / sin_phi;
149  const double K3 = 0.5 * (1.0 + alpha_r) * sin_phi - 0.5 * (1.0 - alpha_r);
150 
151  // Check Modified Mohr-Coulomb criterion
152  if (std::abs(I1) < tolerance) {
153  rEquivalentStress = 0.0;
154  } else {
156  rEquivalentStress = (2.0 * std::tan(Globals::Pi * 0.25 + friction_angle * 0.5) / std::cos(friction_angle)) * ((I1 * K3 / 3.0) +
157  std::sqrt(J2) * (K1 * std::cos(theta) - K2 * std::sin(theta) * sin_phi / std::sqrt(3.0)));
158  }
159  }
160 
166  static void GetInitialUniaxialThreshold(ConstitutiveLaw::Parameters& rValues, double& rThreshold)
167  {
168  const Properties& r_material_properties = rValues.GetMaterialProperties();
169 
170  const double yield_compression = r_material_properties.Has(YIELD_STRESS) ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_COMPRESSION];
171  rThreshold = std::abs(yield_compression);
172  }
173 
183  const array_1d<double, VoigtSize>& rPredictiveStressVector,
184  const array_1d<double, VoigtSize>& rDeviator,
185  const double J2,
188  )
189  {
190  TPlasticPotentialType::CalculatePlasticPotentialDerivative(rPredictiveStressVector, rDeviator, J2, GFlux, rValues);
191  }
192 
201  double& rAParameter,
202  const double CharacteristicLength
203  )
204  {
205  const Properties& r_material_properties = rValues.GetMaterialProperties();
206 
207  const double fracture_energy = r_material_properties[FRACTURE_ENERGY];
208  const double young_modulus = r_material_properties[YOUNG_MODULUS];
209  const bool has_symmetric_yield_stress = r_material_properties.Has(YIELD_STRESS);
210  const double yield_compression = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_COMPRESSION];
211  const double yield_tension = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_TENSION];
212  const double n = yield_compression / yield_tension;
213 
214  if (r_material_properties[SOFTENING_TYPE] == static_cast<int>(SofteningType::Exponential)) {
215  rAParameter = 1.00 / (fracture_energy * n * n * young_modulus / (CharacteristicLength * std::pow(yield_compression, 2)) - 0.5);
216  KRATOS_DEBUG_ERROR_IF(rAParameter < 0.0) << "Fracture energy is too low, increase FRACTURE_ENERGY..." << std::endl;
217  } else { // linear
218  rAParameter = -std::pow(yield_compression, 2) / (2.0 * young_modulus * fracture_energy * n * n / CharacteristicLength);
219  }
220  }
221 
234  const array_1d<double, VoigtSize>& rPredictiveStressVector,
235  const array_1d<double, VoigtSize>& rDeviator,
236  const double J2,
239  {
240  const Properties& r_material_properties = rValues.GetMaterialProperties();
241 
242  array_1d<double, VoigtSize> first_vector, second_vector, third_vector;
243 
247 
248  double J3, lode_angle;
251 
252  const double checker = std::abs(lode_angle * 180.0 / Globals::Pi);
253 
254  double c1, c2, c3;
255  double friction_angle = r_material_properties[FRICTION_ANGLE] * Globals::Pi / 180.0;
256 
257  // Check input variables
258  if (friction_angle < tolerance) {
259  friction_angle = 32.0 * Globals::Pi / 180.0;
260  KRATOS_WARNING("ModifiedMohrCoulombYieldSurface") << "Friction Angle not defined, assumed equal to 32 deg " << std::endl;
261  }
262 
263  const double sin_phi = std::sin(friction_angle);
264  const double cons_phi = std::cos(friction_angle);
265  const double sin_theta = std::sin(lode_angle);
266  const double cos_theta = std::cos(lode_angle);
267  const double cos_3theta = std::cos(3.0 * lode_angle);
268  const double tan_theta = std::tan(lode_angle);
269  const double tan_3theta = std::tan(3.0 * lode_angle);
270  const double Root3 = std::sqrt(3.0);
271 
272  const bool has_symmetric_yield_stress = r_material_properties.Has(YIELD_STRESS);
273  const double compr_yield = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_COMPRESSION];
274  const double tens_yield = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_TENSION];
275  const double n = compr_yield / tens_yield;
276 
277  const double angle_phi = (Globals::Pi * 0.25) + friction_angle * 0.5;
278  const double alpha = n / (std::tan(angle_phi) * std::tan(angle_phi));
279 
280  const double CFL = 2.0 * std::tan(angle_phi) / cons_phi;
281 
282  const double K1 = 0.5 * (1.0 + alpha) - 0.5 * (1.0 - alpha) * sin_phi;
283  const double K2 = 0.5 * (1.0 + alpha) - 0.5 * (1.0 - alpha) / sin_phi;
284  const double K3 = 0.5 * (1.0 + alpha) * sin_phi - 0.5 * (1.0 - alpha);
285 
286  if (std::abs(sin_phi) > tolerance)
287  c1 = CFL * K3 / 3.0;
288  else
289  c1 = 0.0; // check
290 
291  if (std::abs(checker) < 29.0) { // Lode angle not greater than pi/6
292  c2 = cos_theta * CFL * (K1 * (1.0 + tan_theta * tan_3theta) + K2 * sin_phi * (tan_3theta - tan_theta) / Root3);
293  c3 = CFL * (K1 * Root3 * sin_theta + K2 * sin_phi * cos_theta) / (2.0 * J2 * cos_3theta);
294  } else {
295  c3 = 0.0;
296  double aux = 1.0;
297  if (lode_angle > tolerance)
298  aux = -1.0;
299  c2 = 0.5 * CFL * (K1 * Root3 + aux * K2 * sin_phi / Root3);
300  }
301  noalias(rFFlux) = c1 * first_vector + c2 * second_vector + c3 * third_vector;
302  }
303 
308  static int Check(const Properties& rMaterialProperties)
309  {
310  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(FRICTION_ANGLE)) << "FRICTION_ANGLE is not a defined value" << std::endl;
311  if (!rMaterialProperties.Has(YIELD_STRESS)) {
312  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(YIELD_STRESS_TENSION)) << "YIELD_STRESS_TENSION is not a defined value" << std::endl;
313  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(YIELD_STRESS_COMPRESSION)) << "YIELD_STRESS_COMPRESSION is not a defined value" << std::endl;
314 
315  const double yield_compression = rMaterialProperties[YIELD_STRESS_COMPRESSION];
316  const double yield_tension = rMaterialProperties[YIELD_STRESS_TENSION];
317 
318  KRATOS_ERROR_IF(yield_compression < tolerance) << "Yield stress in compression almost zero or negative, include YIELD_STRESS_COMPRESSION in definition";
319  KRATOS_ERROR_IF(yield_tension < tolerance) << "Yield stress in tension almost zero or negative, include YIELD_STRESS_TENSION in definition";
320  } else {
321  const double yield_stress = rMaterialProperties[YIELD_STRESS];
322 
323  KRATOS_ERROR_IF(yield_stress < tolerance) << "Yield stress almost zero or negative, include YIELD_STRESS in definition";
324  }
325  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(FRACTURE_ENERGY)) << "FRACTURE_ENERGY is not a defined value" << std::endl;
326  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(YOUNG_MODULUS)) << "YOUNG_MODULUS is not a defined value" << std::endl;
327 
328  return TPlasticPotentialType::Check(rMaterialProperties);
329  }
330 
336  {
337  return false;
338  }
339 
345  static double GetScaleFactorTension(const Properties& rMaterialProperties)
346  {
347  const double yield_compression = rMaterialProperties[YIELD_STRESS_COMPRESSION];
348  const double yield_tension = rMaterialProperties[YIELD_STRESS_TENSION];
349  return yield_compression / yield_tension;
350  }
351 
352 
356 
360 
364 
368 
370 
371  protected:
374 
378 
382 
386 
390 
394 
398 
400  private:
403 
407 
411 
415 
419 
423 
427 
429 
430 }; // Class ModifiedMohrCoulombYieldSurface
431 
433 
436 
440 
442 
443 } // namespace Kratos.
static void CalculateJ2Invariant(const TVector &rStressVector, const double I1, BoundedVectorType &rDeviator, double &rJ2)
This method computes the second invariant of J.
Definition: advanced_constitutive_law_utilities.h:157
static void CalculateSecondVector(const BoundedVectorType &rDeviator, const double J2, BoundedVectorType &rSecondVector)
This method computes the first vector to be used in the derivative of the yield surface.
Definition: advanced_constitutive_law_utilities.cpp:100
static void CalculateFirstVector(BoundedVectorType &rFirstVector)
This method computes the first vector to be used in the derivative of the yield surface.
Definition: advanced_constitutive_law_utilities.cpp:80
static void CalculateThirdVector(const BoundedVectorType &rDeviator, const double J2, BoundedVectorType &rThirdVector)
This method computes the third vector to be used in the derivative of the yield surface.
Definition: advanced_constitutive_law_utilities.cpp:131
static void CalculateLodeAngle(const double J2, const double J3, double &rLodeAngle)
This method computes the lode angle.
Definition: advanced_constitutive_law_utilities.cpp:158
static void CalculateI1Invariant(const TVector &rStressVector, double &rI1)
This method computes the first invariant from a given stress vector.
Definition: advanced_constitutive_law_utilities.h:116
static void CalculateJ3Invariant(const BoundedVectorType &rDeviator, double &rJ3)
This method computes the third invariant of J.
Definition: advanced_constitutive_law_utilities.cpp:62
This class defines a yield surface according to Modified Mohr-Coulumb theory.
Definition: modified_mohr_coulomb_yield_surface.h:58
KRATOS_CLASS_POINTER_DEFINITION(ModifiedMohrCoulombYieldSurface)
Counted pointer of ModifiedMohrCoulombYieldSurface.
static bool IsWorkingWithTensionThreshold()
This method returns true if the yield surfacecompares with the tension tield stress.
Definition: modified_mohr_coulomb_yield_surface.h:335
ModifiedMohrCoulombYieldSurface()
Initialization constructor.
Definition: modified_mohr_coulomb_yield_surface.h:83
virtual ~ModifiedMohrCoulombYieldSurface()
Destructor.
Definition: modified_mohr_coulomb_yield_surface.h:99
static constexpr double tolerance
The machine precision zero tolerance.
Definition: modified_mohr_coulomb_yield_surface.h:76
static constexpr SizeType VoigtSize
The Plastic potential already defines the Voigt size.
Definition: modified_mohr_coulomb_yield_surface.h:70
static void GetInitialUniaxialThreshold(ConstitutiveLaw::Parameters &rValues, double &rThreshold)
This method returns the initial uniaxial stress threshold.
Definition: modified_mohr_coulomb_yield_surface.h:166
static constexpr SizeType Dimension
The Plastic potential already defines the working simension size.
Definition: modified_mohr_coulomb_yield_surface.h:67
static void CalculateYieldSurfaceDerivative(const array_1d< double, VoigtSize > &rPredictiveStressVector, const array_1d< double, VoigtSize > &rDeviator, const double J2, array_1d< double, VoigtSize > &rFFlux, ConstitutiveLaw::Parameters &rValues)
This script calculates the derivatives of the Yield Surf according to NAYAK-ZIENKIEWICZ paper Interna...
Definition: modified_mohr_coulomb_yield_surface.h:233
TPlasticPotentialType PlasticPotentialType
The type of potential plasticity.
Definition: modified_mohr_coulomb_yield_surface.h:64
static void CalculatePlasticPotentialDerivative(const array_1d< double, VoigtSize > &rPredictiveStressVector, const array_1d< double, VoigtSize > &rDeviator, const double J2, array_1d< double, VoigtSize > &GFlux, ConstitutiveLaw::Parameters &rValues)
This method calculates the derivative of the plastic potential DG/DS.
Definition: modified_mohr_coulomb_yield_surface.h:182
ModifiedMohrCoulombYieldSurface(ModifiedMohrCoulombYieldSurface const &rOther)
Copy constructor.
Definition: modified_mohr_coulomb_yield_surface.h:88
static double GetScaleFactorTension(const Properties &rMaterialProperties)
This method returns the scaling factor of the yield surface surfacecompares with the tension tield st...
Definition: modified_mohr_coulomb_yield_surface.h:345
static void CalculateEquivalentStress(const array_1d< double, VoigtSize > &rPredictiveStressVector, const Vector &rStrainVector, double &rEquivalentStress, ConstitutiveLaw::Parameters &rValues)
This method the uniaxial equivalent stress.
Definition: modified_mohr_coulomb_yield_surface.h:115
static int Check(const Properties &rMaterialProperties)
This method defines the check to be performed in the yield surface.
Definition: modified_mohr_coulomb_yield_surface.h:308
static void CalculateDamageParameter(ConstitutiveLaw::Parameters &rValues, double &rAParameter, const double CharacteristicLength)
This method returns the damage parameter needed in the exp/linear expressions of damage.
Definition: modified_mohr_coulomb_yield_surface.h:199
ModifiedMohrCoulombYieldSurface & operator=(ModifiedMohrCoulombYieldSurface const &rOther)
Assignment operator.
Definition: modified_mohr_coulomb_yield_surface.h:93
Properties encapsulates data shared by different Elements or Conditions. It can store any type of dat...
Definition: properties.h:69
bool Has(TVariableType const &rThisVariable) const
Definition: properties.h:578
#define KRATOS_ERROR_IF_NOT(conditional)
Definition: exception.h:163
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
#define KRATOS_DEBUG_ERROR_IF(conditional)
Definition: exception.h:171
#define KRATOS_WARNING(label)
Definition: logger.h:265
constexpr double Pi
Definition: global_variables.h:25
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
std::size_t SizeType
The definition of the size type.
Definition: mortar_classes.h:43
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
alpha
Definition: generate_convection_diffusion_explicit_element.py:113
float J2
Definition: isotropic_damage_automatic_differentiation.py:133
float K2
Definition: isotropic_damage_automatic_differentiation.py:178
float K3
Definition: isotropic_damage_automatic_differentiation.py:179
I1
Definition: isotropic_damage_automatic_differentiation.py:230
tuple alpha_r
Definition: isotropic_damage_automatic_differentiation.py:174
R
Definition: isotropic_damage_automatic_differentiation.py:172
CFL
Definition: isotropic_damage_automatic_differentiation.py:156
def J3
Definition: isotropic_damage_automatic_differentiation.py:176
sin_phi
Definition: isotropic_damage_automatic_differentiation.py:153
float K1
Definition: isotropic_damage_automatic_differentiation.py:177
int n
manufactured solution and derivatives (u=0 at z=0 dudz=0 at z=domain_height)
Definition: ode_solve.py:402
Definition: constitutive_law.h:189
const Properties & GetMaterialProperties()
Definition: constitutive_law.h:457