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.
tresca_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 
99  virtual ~TrescaYieldSurface(){};
100 
104 
108 
117  const array_1d<double, VoigtSize>& rPredictiveStressVector,
118  const Vector& rStrainVector,
119  double& rEquivalentStress,
121  )
122  {
123  double I1, J2, J3, lode_angle;
125 
127  AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateJ2Invariant(rPredictiveStressVector, I1, deviator, J2);
130 
131  rEquivalentStress = 2.0 * std::cos(lode_angle) * std::sqrt(J2);
132  }
133 
141  double& rThreshold
142  )
143  {
144  const Properties& r_material_properties = rValues.GetMaterialProperties();
145 
146  const double yield_tension = r_material_properties.Has(YIELD_STRESS) ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_TENSION];
147  rThreshold = std::abs(yield_tension);
148  }
149 
158  double& rAParameter,
159  const double CharacteristicLength
160  )
161  {
162  const Properties& r_material_properties = rValues.GetMaterialProperties();
163 
164  const double fracture_energy = r_material_properties[FRACTURE_ENERGY];
165  const double young_modulus = r_material_properties[YOUNG_MODULUS];
166  const bool has_symmetric_yield_stress = r_material_properties.Has(YIELD_STRESS);
167  const double yield_compression = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_COMPRESSION];
168  const double yield_tension = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_TENSION];
169  const double n = yield_compression / yield_tension;
170 
171  if (r_material_properties[SOFTENING_TYPE] == static_cast<int>(SofteningType::Exponential)) {
172  rAParameter = 1.00 / (fracture_energy * n * n * young_modulus / (CharacteristicLength * std::pow(yield_compression, 2)) - 0.5);
173  KRATOS_ERROR_IF(rAParameter < 0.0) << "Fracture energy is too low, increase FRACTURE_ENERGY..." << std::endl;
174  } else { // linear
175  rAParameter = -std::pow(yield_compression, 2) / (2.0 * young_modulus * fracture_energy * n * n / CharacteristicLength);
176  }
177  }
178 
188  const array_1d<double, VoigtSize>& rPredictiveStressVector,
189  const array_1d<double, VoigtSize>& rDeviator,
190  const double J2,
191  array_1d<double, VoigtSize>& rDerivativePlasticPotential,
193  )
194  {
195  TPlasticPotentialType::CalculatePlasticPotentialDerivative(rPredictiveStressVector, rDeviator, J2, rDerivativePlasticPotential, rValues);
196  }
197 
210  const array_1d<double, VoigtSize>& rPredictiveStressVector,
211  const array_1d<double, VoigtSize>& rDeviator,
212  const double J2,
215  )
216  {
217  array_1d<double, VoigtSize> second_vector, third_vector;
218 
221 
222  double J3, lode_angle;
225 
226  const double checker = std::abs(lode_angle * 180.0 / Globals::Pi);
227 
228  double c2, c3;
229  if (std::abs(checker) < 29.0) { // the lode_angle cannot be greater than pi/6
230  c2 = 2.0 * (std::cos(lode_angle) + std::sin(lode_angle) * std::tan(3.0 * lode_angle));
231  c3 = std::sqrt(3.0) * std::sin(lode_angle) / (J2 * std::cos(3.0 * lode_angle));
232  } else {
233  c2 = std::sqrt(3.0);
234  c3 = 0.0;
235  }
236 
237  noalias(rFFlux) = c2 * second_vector + c3 * third_vector;
238  }
239 
244  static int Check(const Properties& rMaterialProperties)
245  {
246  if (!rMaterialProperties.Has(YIELD_STRESS)) {
247  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(YIELD_STRESS_TENSION)) << "YIELD_STRESS_TENSION is not a defined value" << std::endl;
248  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(YIELD_STRESS_COMPRESSION)) << "YIELD_STRESS_COMPRESSION is not a defined value" << std::endl;
249 
250  const double yield_compression = rMaterialProperties[YIELD_STRESS_COMPRESSION];
251  const double yield_tension = rMaterialProperties[YIELD_STRESS_TENSION];
252 
253  KRATOS_ERROR_IF(yield_compression < tolerance) << "Yield stress in compression almost zero or negative, include YIELD_STRESS_COMPRESSION in definition";
254  KRATOS_ERROR_IF(yield_tension < tolerance) << "Yield stress in tension almost zero or negative, include YIELD_STRESS_TENSION in definition";
255  } else {
256  const double yield_stress = rMaterialProperties[YIELD_STRESS];
257 
258  KRATOS_ERROR_IF(yield_stress < tolerance) << "Yield stress almost zero or negative, include YIELD_STRESS in definition";
259  }
260  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(FRACTURE_ENERGY)) << "FRACTURE_ENERGY is not a defined value" << std::endl;
261  KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(YOUNG_MODULUS)) << "YOUNG_MODULUS is not a defined value" << std::endl;
262 
263  return TPlasticPotentialType::Check(rMaterialProperties);
264  }
265 
270  {
271  return true;
272  }
273 
277  static double GetScaleFactorTension(const Properties& rMaterialProperties)
278  {
279  return 1.0;
280  }
281 
285 
289 
293 
297 
299 
300 protected:
303 
307 
311 
315 
319 
323 
327 
329 private:
332 
336 
340 
344 
348 
352 
356 
358 
359 }; // Class TrescaYieldSurface
360 
362 
365 
369 
371 
372 } // 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 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
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
This class defines a yield surface according to Tresca theory.
Definition: tresca_yield_surface.h:58
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: tresca_yield_surface.h:209
TrescaYieldSurface()
Initialization constructor.
Definition: tresca_yield_surface.h:83
virtual ~TrescaYieldSurface()
Destructor.
Definition: tresca_yield_surface.h:99
TrescaYieldSurface(TrescaYieldSurface const &rOther)
Copy constructor.
Definition: tresca_yield_surface.h:88
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: tresca_yield_surface.h:156
KRATOS_CLASS_POINTER_DEFINITION(TrescaYieldSurface)
Counted pointer of TrescaYieldSurface.
static bool IsWorkingWithTensionThreshold()
This method returns true if the yield surfacecompares with the tension tield stress.
Definition: tresca_yield_surface.h:269
static constexpr double tolerance
The machine precision zero tolerance.
Definition: tresca_yield_surface.h:76
TPlasticPotentialType PlasticPotentialType
The type of potential plasticity.
Definition: tresca_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 > &rDerivativePlasticPotential, ConstitutiveLaw::Parameters &rValues)
This method calculates the derivative of the plastic potential DG/DS.
Definition: tresca_yield_surface.h:187
static void CalculateEquivalentStress(const array_1d< double, VoigtSize > &rPredictiveStressVector, const Vector &rStrainVector, double &rEquivalentStress, ConstitutiveLaw::Parameters &rValues)
This method the uniaxial equivalent stress.
Definition: tresca_yield_surface.h:116
static constexpr SizeType Dimension
The Plastic potential already defines the working simension size.
Definition: tresca_yield_surface.h:67
static constexpr SizeType VoigtSize
The Plastic potential already defines the Voigt size.
Definition: tresca_yield_surface.h:70
static int Check(const Properties &rMaterialProperties)
This method defines the check to be performed in the yield surface.
Definition: tresca_yield_surface.h:244
static void GetInitialUniaxialThreshold(ConstitutiveLaw::Parameters &rValues, double &rThreshold)
This method returns the initial uniaxial stress threshold.
Definition: tresca_yield_surface.h:139
TrescaYieldSurface & operator=(TrescaYieldSurface const &rOther)
Assignment operator.
Definition: tresca_yield_surface.h:93
static double GetScaleFactorTension(const Properties &rMaterialProperties)
This method returns the scaling factor of the yield surface surfacecompares with the tension tield st...
Definition: tresca_yield_surface.h:277
#define KRATOS_ERROR_IF_NOT(conditional)
Definition: exception.h:163
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
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
float J2
Definition: isotropic_damage_automatic_differentiation.py:133
I1
Definition: isotropic_damage_automatic_differentiation.py:230
def J3
Definition: isotropic_damage_automatic_differentiation.py:176
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