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.
thermal_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
10 //
11 
12 #pragma once
13 
14 // System includes
15 
16 // Project includes
18 
19 namespace Kratos
20 {
23 
27 
31 
35 
39 
53 template<class TPlasticPotentialType>
55 {
56 public:
59 
61  using PlasticPotentialType = TPlasticPotentialType;
62 
64 
66  static constexpr SizeType Dimension = PlasticPotentialType::Dimension;
67 
69  static constexpr SizeType VoigtSize = PlasticPotentialType::VoigtSize;
70 
73 
75  static constexpr double tolerance = std::numeric_limits<double>::epsilon();
76 
79 
83 
86  {
87  }
88 
91  {
92  }
93 
96  {
97  return *this;
98  }
99 
102 
109 
118  array_1d<double, VoigtSize>& rStressVector,
119  const Vector& rStrainVector,
120  double& rEquivalentStress,
122  )
123  {
124 
125  double friction_angle = AdvCLutils::GetMaterialPropertyThroughAccessor(FRICTION_ANGLE, rValues) * 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;
136  AdvCLutils::CalculateI1Invariant(rStressVector, I1);
138  AdvCLutils::CalculateJ2Invariant(rStressVector, 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 auto& r_material_properties = rValues.GetMaterialProperties();
157 
158  double yield_tension, friction_angle;
159  if (rValues.IsSetShapeFunctionsValues()) { // This is needed since at Initialize level the N are not set yet...
160  yield_tension = r_material_properties.Has(YIELD_STRESS) ? AdvCLutils::GetMaterialPropertyThroughAccessor(YIELD_STRESS, rValues) : AdvCLutils::GetMaterialPropertyThroughAccessor(YIELD_STRESS_TENSION, rValues);
161  friction_angle = AdvCLutils::GetMaterialPropertyThroughAccessor(FRICTION_ANGLE, rValues) * Globals::Pi / 180.0;
162  } else {
163  const double ref_temperature = r_material_properties.Has(REFERENCE_TEMPERATURE) ? r_material_properties[REFERENCE_TEMPERATURE] : rValues.GetElementGeometry().GetValue(REFERENCE_TEMPERATURE);
164  yield_tension = r_material_properties.Has(YIELD_STRESS) ? AdvCLutils::GetPropertyFromTemperatureTable(YIELD_STRESS, rValues, ref_temperature) : AdvCLutils::GetPropertyFromTemperatureTable(YIELD_STRESS_TENSION, rValues, ref_temperature);
165  friction_angle = AdvCLutils::GetPropertyFromTemperatureTable(FRICTION_ANGLE, rValues, ref_temperature) * Globals::Pi / 180.0;
166  }
167  const double sin_phi = std::sin(friction_angle);
168  rThreshold = std::abs(yield_tension * (3.0 + sin_phi) / (3.0 * sin_phi - 3.0));
169  }
170 
179  double& rAParameter,
180  const double CharacteristicLength
181  )
182  {
183  const auto& r_material_properties = rValues.GetMaterialProperties();
184 
185  const auto &r_geom = rValues.GetElementGeometry();
186  const auto &r_N = rValues.GetShapeFunctionsValues();
187  const auto &r_process_info = rValues.GetProcessInfo();
188 
189  const double fracture_energy = r_material_properties.GetValue(FRACTURE_ENERGY, r_geom, r_N, r_process_info);
190  const double young_modulus = r_material_properties.GetValue(YOUNG_MODULUS, r_geom, r_N, r_process_info);
191  const bool has_symmetric_yield_stress = r_material_properties.Has(YIELD_STRESS);
192 
193  const double yield_compression = has_symmetric_yield_stress ? r_material_properties.GetValue(YIELD_STRESS, r_geom, r_N, r_process_info) : r_material_properties.GetValue(YIELD_STRESS_COMPRESSION, r_geom, r_N, r_process_info);
194  const double yield_tension = has_symmetric_yield_stress ? r_material_properties.GetValue(YIELD_STRESS, r_geom, r_N, r_process_info) : r_material_properties.GetValue(YIELD_STRESS_TENSION, r_geom, r_N, r_process_info);
195  const double n = yield_compression / yield_tension;
196 
197  if (r_material_properties[SOFTENING_TYPE] == static_cast<int>(SofteningType::Exponential)) {
198  rAParameter = 1.00 / (fracture_energy * n * n * young_modulus / (CharacteristicLength * std::pow(yield_compression, 2)) - 0.5);
199  KRATOS_ERROR_IF(rAParameter < 0.0) << "Fracture energy is too low, increase FRACTURE_ENERGY..." << std::endl;
200  } else if (r_material_properties[SOFTENING_TYPE] == static_cast<int>(SofteningType::Linear)) { // linear
201  rAParameter = -std::pow(yield_compression, 2) / (2.0 * young_modulus * fracture_energy * n * n / CharacteristicLength);
202  }
203  }
204 
214  const array_1d<double, VoigtSize>& rPredictiveStressVector,
215  const array_1d<double, VoigtSize>& rDeviator,
216  const double J2,
219  )
220  {
221  TPlasticPotentialType::CalculatePlasticPotentialDerivative(rPredictiveStressVector, rDeviator, J2, rGFlux, rValues);
222  }
223 
236  const array_1d<double, VoigtSize>& rPredictiveStressVector,
237  const array_1d<double, VoigtSize>& rDeviator,
238  const double J2,
241  )
242  {
243  array_1d<double, VoigtSize> first_vector, second_vector, third_vector;
246 
247  const double friction_angle = AdvCLutils::GetMaterialPropertyThroughAccessor(FRICTION_ANGLE, rValues) * Globals::Pi / 180.0;;
248  const double sin_phi = std::sin(friction_angle);
249  const double Root3 = std::sqrt(3.0);
250 
251  const double CFL = -Root3 * (3.0 - sin_phi) / (3.0 * sin_phi - 3.0);
252  const double c1 = CFL * 2.0 * sin_phi / (Root3 * (3.0 - sin_phi));
253  const double c2 = CFL;
254 
255  noalias(rFFlux) = c1 * first_vector + c2 * second_vector;
256  }
257 
262  static int Check(const Properties& rMaterialProperties)
263  {
264  return BaseType::Check(rMaterialProperties);
265  }
266 
272  {
274  }
275 
280  static double GetScaleFactorTension(const Properties& rMaterialProperties)
281  {
283  }
284 
288 
292 
296 
300 
302 
303 protected:
306 
310 
314 
318 
322 
326 
330 
332 private:
335 
339 
343 
347 
351 
355 
359 
361 
362 }; // Class DruckerPragerYieldSurface
363 
365 
368 
372 
374 
375 } // namespace Kratos.
This class includes several utilities necessaries for the computation of the constitutive law.
Definition: advanced_constitutive_law_utilities.h:59
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 double GetPropertyFromTemperatureTable(const Variable< double > &rVariable, ConstitutiveLaw::Parameters &rValues, const double Temperature)
This retrieves a double type variable from a table if exists, assumes TEMPERATURE to be the independe...
Definition: advanced_constitutive_law_utilities.cpp:833
static double GetMaterialPropertyThroughAccessor(const Variable< double > &rVariable, ConstitutiveLaw::Parameters &rValues)
This retrieves a double type variable checking the accessor.
Definition: advanced_constitutive_law_utilities.cpp:818
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
TDataType & GetValue(const Variable< TDataType > &rThisVariable)
Gets the value associated with a given variable.
Definition: data_value_container.h:268
This class defines a yield surface according to Drucker-Prager theory.
Definition: drucker_prager_yield_surface.h:59
static bool IsWorkingWithTensionThreshold()
This method returns true if the yield surfacecompares with the tension tield stress.
Definition: drucker_prager_yield_surface.h:279
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 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
TVariableType::Type & GetValue(const TVariableType &rThisVariable)
Definition: geometry.h:627
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 Drucker-Prager theory.
Definition: thermal_drucker_prager_yield_surface.h:55
TPlasticPotentialType PlasticPotentialType
The type of potential plasticity.
Definition: thermal_drucker_prager_yield_surface.h:61
static double GetScaleFactorTension(const Properties &rMaterialProperties)
This method returns the scaling factor of the yield surface compares with the tension yield stress.
Definition: thermal_drucker_prager_yield_surface.h:280
virtual ~ThermalDruckerPragerYieldSurface()
Destructor.
Definition: thermal_drucker_prager_yield_surface.h:101
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: thermal_drucker_prager_yield_surface.h:213
static constexpr SizeType Dimension
The Plastic potential already defines the working simension size.
Definition: thermal_drucker_prager_yield_surface.h:66
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: thermal_drucker_prager_yield_surface.h:177
ThermalDruckerPragerYieldSurface & operator=(ThermalDruckerPragerYieldSurface const &rOther)
Assignment operator.
Definition: thermal_drucker_prager_yield_surface.h:95
ThermalDruckerPragerYieldSurface()
Initialization constructor.
Definition: thermal_drucker_prager_yield_surface.h:85
static bool IsWorkingWithTensionThreshold()
This method returns true if the yield surfacecompares with the tension tield stress.
Definition: thermal_drucker_prager_yield_surface.h:271
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: thermal_drucker_prager_yield_surface.h:235
KRATOS_CLASS_POINTER_DEFINITION(ThermalDruckerPragerYieldSurface)
Counted pointer of ThermalDruckerPragerYieldSurface.
static int Check(const Properties &rMaterialProperties)
This method defines the check to be performed in the yield surface.
Definition: thermal_drucker_prager_yield_surface.h:262
static constexpr double tolerance
The machine precision zero tolerance.
Definition: thermal_drucker_prager_yield_surface.h:75
static void GetInitialUniaxialThreshold(ConstitutiveLaw::Parameters &rValues, double &rThreshold)
This method returns the initial uniaxial stress threshold.
Definition: thermal_drucker_prager_yield_surface.h:151
ThermalDruckerPragerYieldSurface(ThermalDruckerPragerYieldSurface const &rOther)
Copy constructor.
Definition: thermal_drucker_prager_yield_surface.h:90
static constexpr SizeType VoigtSize
The Plastic potential already defines the Voigt size.
Definition: thermal_drucker_prager_yield_surface.h:69
static void CalculateEquivalentStress(array_1d< double, VoigtSize > &rStressVector, const Vector &rStrainVector, double &rEquivalentStress, ConstitutiveLaw::Parameters &rValues)
This method the uniaxial equivalent stress.
Definition: thermal_drucker_prager_yield_surface.h:117
#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
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
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 GeometryType & GetElementGeometry()
Definition: constitutive_law.h:462
const Vector & GetShapeFunctionsValues()
Definition: constitutive_law.h:419
bool IsSetShapeFunctionsValues()
Definition: constitutive_law.h:483
const ProcessInfo & GetProcessInfo()
Definition: constitutive_law.h:452
const Properties & GetMaterialProperties()
Definition: constitutive_law.h:457