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.
drucker_prager_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 
57 template<class TPlasticPotentialType>
59 {
60 public:
63 
65  typedef TPlasticPotentialType PlasticPotentialType;
66 
68  static constexpr SizeType Dimension = PlasticPotentialType::Dimension;
69 
71  static constexpr SizeType VoigtSize = PlasticPotentialType::VoigtSize;
72 
75 
77  static constexpr double tolerance = std::numeric_limits<double>::epsilon();
78 
82 
85  {
86  }
87 
90  {
91  }
92 
95  {
96  return *this;
97  }
98 
101 
108 
117  array_1d<double, VoigtSize>& rPredictiveStressVector,
118  const Vector& rStrainVector,
119  double& rEquivalentStress,
121  )
122  {
123  const Properties& r_material_properties = rValues.GetMaterialProperties();
124 
125  double friction_angle = r_material_properties[FRICTION_ANGLE] * Globals::Pi / 180.0; // In radians!
126  const double sin_phi = std::sin(friction_angle);
127  const double root_3 = std::sqrt(3.0);
128 
129  // Check input variables
130  if (friction_angle < tolerance) {
131  friction_angle = 32.0 * Globals::Pi / 180.0;
132  KRATOS_WARNING("DruckerPragerYieldSurface") << "Friction Angle not defined, assumed equal to 32 " << std::endl;
133  }
134 
135  double I1, J2;
138  AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateJ2Invariant(rPredictiveStressVector, I1, deviator, J2);
139 
140 
141  const double CFL = -root_3 * (3.0 - sin_phi) / (3.0 * sin_phi - 3.0);
142  const double TEN0 = 2.0 * I1 * sin_phi / (root_3 * (3.0 - sin_phi)) + std::sqrt(J2);
143  rEquivalentStress = (CFL * TEN0);
144  }
145 
153  double& rThreshold
154  )
155  {
156  const Properties& r_material_properties = rValues.GetMaterialProperties();
157 
158  const double yield_tension = r_material_properties.Has(YIELD_STRESS) ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_TENSION];
159  const double friction_angle = r_material_properties[FRICTION_ANGLE] * Globals::Pi / 180.0; // In radians!
160  const double sin_phi = std::sin(friction_angle);
161  rThreshold = std::abs(yield_tension * (3.0 + sin_phi) / (3.0 * sin_phi - 3.0));
162  }
163 
172  double& rAParameter,
173  const double CharacteristicLength
174  )
175  {
176  const Properties& r_material_properties = rValues.GetMaterialProperties();
177 
178  const double Gf = r_material_properties[FRACTURE_ENERGY];
179  const double E = r_material_properties[YOUNG_MODULUS];
180  const bool has_symmetric_yield_stress = r_material_properties.Has(YIELD_STRESS);
181  const double yield_compression = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_COMPRESSION];
182  const double yield_tension = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_TENSION];
183  const double n = yield_compression / yield_tension;
184 
185  if (r_material_properties[SOFTENING_TYPE] == static_cast<int>(SofteningType::Exponential)) {
186  rAParameter = 1.00 / (Gf * n * n * E / (CharacteristicLength * std::pow(yield_compression, 2)) - 0.5);
187  KRATOS_ERROR_IF(rAParameter < 0.0) << "Fracture energy is too low, increase FRACTURE_ENERGY..." << std::endl;
188  } else { // linear
189  rAParameter = -std::pow(yield_compression, 2) / (2.0 * E * Gf * n * n / CharacteristicLength);
190  }
191  }
192 
202  const array_1d<double, VoigtSize>& rPredictiveStressVector,
203  const array_1d<double, VoigtSize>& rDeviator,
204  const double J2,
207  )
208  {
209  TPlasticPotentialType::CalculatePlasticPotentialDerivative(rPredictiveStressVector, rDeviator, J2, rGFlux, rValues);
210  }
211 
224  const array_1d<double, VoigtSize>& rPredictiveStressVector,
225  const array_1d<double, VoigtSize>& rDeviator,
226  const double J2,
229  )
230  {
231  const Properties& r_material_properties = rValues.GetMaterialProperties();
232 
233  array_1d<double, VoigtSize> first_vector, second_vector, third_vector;
236 
237  const double friction_angle = r_material_properties[FRICTION_ANGLE] * Globals::Pi / 180.0;;
238  const double sin_phi = std::sin(friction_angle);
239  const double Root3 = std::sqrt(3.0);
240 
241  const double CFL = -Root3 * (3.0 - sin_phi) / (3.0 * sin_phi - 3.0);
242  const double c1 = CFL * 2.0 * sin_phi / (Root3 * (3.0 - sin_phi));
243  const double c2 = CFL;
244 
245  noalias(rFFlux) = c1 * first_vector + c2 * second_vector;
246  }
247 
252  static int Check(const Properties& rMaterialProperties)
253  {
254  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(FRICTION_ANGLE)) << "FRICTION_ANGLE is not a defined value" << std::endl;
255  if (!rMaterialProperties.Has(YIELD_STRESS)) {
256  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(YIELD_STRESS_TENSION)) << "YIELD_STRESS_TENSION is not a defined value" << std::endl;
257  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(YIELD_STRESS_COMPRESSION)) << "YIELD_STRESS_COMPRESSION is not a defined value" << std::endl;
258 
259  const double yield_compression = rMaterialProperties[YIELD_STRESS_COMPRESSION];
260  const double yield_tension = rMaterialProperties[YIELD_STRESS_TENSION];
261 
262  KRATOS_ERROR_IF(yield_compression < tolerance) << "Yield stress in compression almost zero or negative, include YIELD_STRESS_COMPRESSION in definition";
263  KRATOS_ERROR_IF(yield_tension < tolerance) << "Yield stress in tension almost zero or negative, include YIELD_STRESS_TENSION in definition";
264  } else {
265  const double yield_stress = rMaterialProperties[YIELD_STRESS];
266 
267  KRATOS_ERROR_IF(yield_stress < tolerance) << "Yield stress almost zero or negative, include YIELD_STRESS in definition";
268  }
269  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(FRACTURE_ENERGY)) << "FRACTURE_ENERGY is not a defined value" << std::endl;
270  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(YOUNG_MODULUS)) << "YOUNG_MODULUS is not a defined value" << std::endl;
271 
272  return TPlasticPotentialType::Check(rMaterialProperties);
273  }
274 
280  {
281  return true;
282  }
283 
289  static double GetScaleFactorTension(const Properties& rMaterialProperties)
290  {
291  const double friction_angle = rMaterialProperties[FRICTION_ANGLE] * Globals::Pi / 180.0; // In radians!
292  const double sin_phi = std::sin(friction_angle);
293  return 1.0 / std::abs((3.0 + sin_phi) / (3.0 * sin_phi - 3.0));
294  }
295 
299 
303 
307 
311 
313 
314 protected:
317 
321 
325 
329 
333 
337 
341 
343 private:
346 
350 
354 
358 
362 
366 
370 
372 
373 }; // Class DruckerPragerYieldSurface
374 
376 
379 
383 
385 
386 } // 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 CalculateI1Invariant(const TVector &rStressVector, double &rI1)
This method computes the first invariant from a given stress vector.
Definition: advanced_constitutive_law_utilities.h:116
This class defines a yield surface according to Drucker-Prager theory.
Definition: drucker_prager_yield_surface.h:59
static void GetInitialUniaxialThreshold(ConstitutiveLaw::Parameters &rValues, double &rThreshold)
This method returns the initial uniaxial stress threshold.
Definition: drucker_prager_yield_surface.h:151
static bool IsWorkingWithTensionThreshold()
This method returns true if the yield surfacecompares with the tension tield stress.
Definition: drucker_prager_yield_surface.h:279
DruckerPragerYieldSurface & operator=(DruckerPragerYieldSurface const &rOther)
Assignment operator.
Definition: drucker_prager_yield_surface.h:94
static void CalculateEquivalentStress(array_1d< double, VoigtSize > &rPredictiveStressVector, const Vector &rStrainVector, double &rEquivalentStress, ConstitutiveLaw::Parameters &rValues)
This method the uniaxial equivalent stress.
Definition: drucker_prager_yield_surface.h:116
DruckerPragerYieldSurface()
Initialization constructor.
Definition: drucker_prager_yield_surface.h:84
static constexpr SizeType VoigtSize
The Plastic potential already defines the Voigt size.
Definition: drucker_prager_yield_surface.h:71
static int Check(const Properties &rMaterialProperties)
This method defines the check to be performed in the yield surface.
Definition: drucker_prager_yield_surface.h:252
static constexpr double tolerance
The machine precision zero tolerance.
Definition: drucker_prager_yield_surface.h:77
TPlasticPotentialType PlasticPotentialType
The type of potential plasticity.
Definition: drucker_prager_yield_surface.h:65
static constexpr SizeType Dimension
The Plastic potential already defines the working simension size.
Definition: drucker_prager_yield_surface.h:68
KRATOS_CLASS_POINTER_DEFINITION(DruckerPragerYieldSurface)
Counted pointer of DruckerPragerYieldSurface.
virtual ~DruckerPragerYieldSurface()
Destructor.
Definition: drucker_prager_yield_surface.h:100
static double GetScaleFactorTension(const Properties &rMaterialProperties)
This method returns the scaling factor of the yield surface surfacecompares with the tension tield st...
Definition: drucker_prager_yield_surface.h:289
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: drucker_prager_yield_surface.h:223
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: drucker_prager_yield_surface.h:170
static void CalculatePlasticPotentialDerivative(const array_1d< double, VoigtSize > &rPredictiveStressVector, const array_1d< double, VoigtSize > &rDeviator, const double J2, array_1d< double, VoigtSize > &rGFlux, ConstitutiveLaw::Parameters &rValues)
This method calculates the derivative of the plastic potential DG/DS.
Definition: drucker_prager_yield_surface.h:201
DruckerPragerYieldSurface(DruckerPragerYieldSurface const &rOther)
Copy constructor.
Definition: drucker_prager_yield_surface.h:89
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_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
std::size_t SizeType
The definition of the size type.
Definition: mortar_classes.h:43
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
E
Definition: generate_hyper_elastic_simo_taylor_neo_hookean.py:26
float J2
Definition: isotropic_damage_automatic_differentiation.py:133
I1
Definition: isotropic_damage_automatic_differentiation.py:230
CFL
Definition: isotropic_damage_automatic_differentiation.py:156
root_3
Definition: isotropic_damage_automatic_differentiation.py:155
Gf
Definition: isotropic_damage_automatic_differentiation.py:135
float TEN0
Definition: isotropic_damage_automatic_differentiation.py:157
sin_phi
Definition: isotropic_damage_automatic_differentiation.py:153
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