thermal_mohr_coulomb_yield_surface.h

Go to the documentation of this file.

// KRATOS  ___|  |                   |                   |
//       \___ \  __|  __| |   |  __| __| |   |  __| _` | |
//             | |   |    |   | (    |   |   | |   (   | |
//       _____/ \__|_|   \__,_|\___|\__|\__,_|_|  \__,_|_| MECHANICS
//
//  License:         BSD License
//                   license: structural_mechanics_application/license.txt
//
//  Main authors:    Alejandro Cornejo
//
 
#pragma once
 
// System includes
 
// Project includes
#include "custom_constitutive/auxiliary_files/yield_surfaces/mohr_coulomb_yield_surface.h"
#include "thermal_von_mises_yield_surface.h"
 
namespace Kratos
{
 
 
 
 
 
template <class TPlasticPotentialType>
class ThermalMohrCoulombYieldSurface
{
public:
 
    using PlasticPotentialType = TPlasticPotentialType;
 
    using BaseType = MohrCoulombYieldSurface<TPlasticPotentialType>;
 
    static constexpr SizeType Dimension = PlasticPotentialType::Dimension;
 
    static constexpr SizeType VoigtSize = PlasticPotentialType::VoigtSize;
 
    KRATOS_CLASS_POINTER_DEFINITION(ThermalMohrCoulombYieldSurface);
 
    static constexpr double tolerance = std::numeric_limits<double>::epsilon();
 
    using AdvCLutils = AdvancedConstitutiveLawUtilities<VoigtSize>;
 
 
    ThermalMohrCoulombYieldSurface()
    {
    }
 
    ThermalMohrCoulombYieldSurface(ThermalMohrCoulombYieldSurface const &rOther)
    {
    }
 
    ThermalMohrCoulombYieldSurface &operator=(ThermalMohrCoulombYieldSurface const &rOther)
    {
        return *this;
    }
 
    virtual ~ThermalMohrCoulombYieldSurface(){};
 
 
 
    static void CalculateEquivalentStress(
        const array_1d<double, VoigtSize>& rStressVector,
        const Vector& rStrainVector,
        double& rEquivalentStress,
        ConstitutiveLaw::Parameters& rValues
        )
    {
        double I1, J2, J3, lode_angle;
        array_1d<double, VoigtSize> deviator = ZeroVector(VoigtSize);
 
        AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateI1Invariant(rStressVector, I1);
        AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateJ2Invariant(rStressVector, I1, deviator, J2);
        AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateJ3Invariant(deviator, J3);
        AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateLodeAngle(J2, J3, lode_angle);
 
        const double friction_angle = AdvCLutils::GetMaterialPropertyThroughAccessor(FRICTION_ANGLE, rValues) * Globals::Pi / 180.0;
 
        rEquivalentStress = (std::cos(lode_angle) - std::sin(lode_angle) * std::sin(friction_angle) / std::sqrt(3.0)) * std::sqrt(J2) +
            I1 * std::sin(friction_angle) / 3.0;
    }
 
    static void GetInitialUniaxialThreshold(
        ConstitutiveLaw::Parameters& rValues,
        double& rThreshold
        )
    {
        const auto& r_material_properties = rValues.GetMaterialProperties();
        double friction_angle, cohesion;
        if (rValues.IsSetShapeFunctionsValues()) { // This is needed since at Initialize level the N are not set yet...
            friction_angle = AdvCLutils::GetMaterialPropertyThroughAccessor(FRICTION_ANGLE, rValues);
            cohesion = AdvCLutils::GetMaterialPropertyThroughAccessor(COHESION, rValues);
        } else {
            const double ref_temperature = r_material_properties.Has(REFERENCE_TEMPERATURE) ? r_material_properties[REFERENCE_TEMPERATURE] : rValues.GetElementGeometry().GetValue(REFERENCE_TEMPERATURE);
            friction_angle = AdvCLutils::GetPropertyFromTemperatureTable(FRICTION_ANGLE, rValues, ref_temperature);
            cohesion = AdvCLutils::GetPropertyFromTemperatureTable(COHESION, rValues, ref_temperature);
        }
        rThreshold = cohesion * std::cos(friction_angle);
    }
 
    static void CalculateDamageParameter(
        ConstitutiveLaw::Parameters& rValues,
        double& rAParameter,
        const double CharacteristicLength
        )
    {
        const Properties& r_material_properties = rValues.GetMaterialProperties();
 
        const auto &r_geom = rValues.GetElementGeometry();
        const auto &r_N = rValues.GetShapeFunctionsValues();
        const auto &r_process_info = rValues.GetProcessInfo();
 
        const double fracture_energy = r_material_properties.GetValue(FRACTURE_ENERGY, r_geom, r_N, r_process_info);
        const double young_modulus   = r_material_properties.GetValue(YOUNG_MODULUS, r_geom, r_N, r_process_info);
 
        double equivalent_yield;
        GetInitialUniaxialThreshold(rValues, equivalent_yield);
        if (r_material_properties[SOFTENING_TYPE] == static_cast<int>(SofteningType::Exponential)) {
            rAParameter = 1.00 / (fracture_energy * young_modulus / (CharacteristicLength * std::pow(equivalent_yield, 2)) - 0.5);
            KRATOS_ERROR_IF(rAParameter < 0.0) << "Fracture Energy is too low, increase FRACTURE_ENERGY..." << std::endl;
        } else if (r_material_properties[SOFTENING_TYPE] == static_cast<int>(SofteningType::Linear)) { // linear
            rAParameter = -std::pow(equivalent_yield, 2) / (2.0 * young_modulus * fracture_energy / CharacteristicLength);
        } else {
            rAParameter = 0.0;
        }
    }
 
    static void CalculatePlasticPotentialDerivative(
        const array_1d<double, VoigtSize>& rStressVector,
        const array_1d<double, VoigtSize>& rDeviator,
        const double J2,
        array_1d<double, VoigtSize>& rDerivativePlasticPotential,
        ConstitutiveLaw::Parameters& rValues
        )
    {
        TPlasticPotentialType::CalculatePlasticPotentialDerivative(rStressVector, rDeviator, J2, rDerivativePlasticPotential, rValues);
    }
 
    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
        )
    {
        array_1d<double, VoigtSize> first_vector, second_vector, third_vector;
        const double friction_angle = AdvCLutils::GetMaterialPropertyThroughAccessor(FRICTION_ANGLE, rValues) * Globals::Pi / 180.0;
 
        AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateFirstVector(first_vector);
        AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateSecondVector(rDeviator, J2, second_vector);
        AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateThirdVector(rDeviator, J2, third_vector);
 
        double J3, lode_angle;
        AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateJ3Invariant(rDeviator, J3);
        AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateLodeAngle(J2, J3, lode_angle);
 
        double c1, c3, c2;
        double checker = std::abs(lode_angle * 180.0 / Globals::Pi);
 
        if (std::abs(checker) < 29.0) { // If it is not the edge
            c1 = std::sin(friction_angle) / 3.0;
            c3 = (std::sqrt(3.0) * std::sin(lode_angle) + std::sin(friction_angle) * std::cos(lode_angle)) /
                (2.0 * J2 * std::cos(3.0 * lode_angle));
            c2 = 0.5 * std::cos(lode_angle)*(1.0 + std::tan(lode_angle) * std::sin(3.0 * lode_angle) +
                std::sin(friction_angle) * (std::tan(3.0 * lode_angle) - std::tan(lode_angle)) / std::sqrt(3.0));
        } else { // smoothing with drucker-prager
            c1 = 3.0 * (2.0 * std::sin(friction_angle) / (std::sqrt(3.0) * (3.0 - std::sin(friction_angle))));
            c2 = 1.0;
            c3 = 0.0;
        }
        noalias(rFFlux) = c1 * first_vector + c2 * second_vector + c3 * third_vector;
    }
 
    static int Check(const Properties& rMaterialProperties)
    {
        return BaseType::Check(rMaterialProperties);
    }
 
    static bool IsWorkingWithTensionThreshold()
    {
        return BaseType::IsWorkingWithTensionThreshold();
    }
 
    static double GetScaleFactorTension(const Properties& rMaterialProperties)
    {
        return BaseType::GetScaleFactorTension();
    }
 
 
 
 
 
 
protected:
 
 
 
 
 
 
 
private:
 
 
 
 
 
 
 
 
}; // Class MohrCoulombYieldSurface
 
 
 
 
 
} // namespace Kratos.