modified_mohr_coulomb_plastic_potential.h

Go to the documentation of this file.

// KRATOS ___                _   _ _         _   _             __                       _
//       / __\___  _ __  ___| |_(_) |_ _   _| |_(_)_   _____  / /  __ ___      _____   /_\  _ __  _ __
//      / /  / _ \| '_ \/ __| __| | __| | | | __| \ \ / / _ \/ /  / _` \ \ /\ / / __| //_\\| '_ \| '_  |
//     / /__| (_) | | | \__ \ |_| | |_| |_| | |_| |\ V /  __/ /__| (_| |\ V  V /\__ \/  _  \ |_) | |_) |
//     \____/\___/|_| |_|___/\__|_|\__|\__,_|\__|_| \_/ \___\____/\__,_| \_/\_/ |___/\_/ \_/ .__/| .__/
//                                                                                         |_|   |_|
//
//  License:         BSD License
//                   license: structural_mechanics_application/license.txt
//
//  Main authors:    Alejandro Cornejo & Lucia Barbu
//
 
#pragma once
 
// System includes
 
// Project includes
#include "includes/checks.h"
#include "generic_plastic_potential.h"
 
namespace Kratos
{
 
 
    // The size type definition
    typedef std::size_t SizeType;
 
 
 
 
template <SizeType TVoigtSize = 6>
class ModifiedMohrCoulombPlasticPotential
{
  public:
 
    static constexpr SizeType Dimension = TVoigtSize == 6 ? 3 : 2;
 
    static constexpr SizeType VoigtSize = TVoigtSize;
 
    KRATOS_CLASS_POINTER_DEFINITION(ModifiedMohrCoulombPlasticPotential);
 
    static constexpr double tolerance = std::numeric_limits<double>::epsilon();
 
 
    ModifiedMohrCoulombPlasticPotential()
    {
    }
 
    ModifiedMohrCoulombPlasticPotential(ModifiedMohrCoulombPlasticPotential const &rOther)
    {
    }
 
    ModifiedMohrCoulombPlasticPotential &operator=(ModifiedMohrCoulombPlasticPotential const &rOther)
    {
        return *this;
    }
 
    virtual ~ModifiedMohrCoulombPlasticPotential(){};
 
 
 
    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
        )
    {
        const Properties& r_material_properties = rValues.GetMaterialProperties();
 
        array_1d<double, VoigtSize> first_vector, second_vector, third_vector;
 
        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);
 
        const double checker = std::abs(lode_angle * 180.0 / Globals::Pi);
 
        const double dilatancy = r_material_properties[DILATANCY_ANGLE] * Globals::Pi / 180.0;
        const double sin_dil = std::sin(dilatancy);
        const double cos_dil = std::cos(dilatancy);
        const double sin_theta = std::sin(lode_angle);
        const double cos_theta = std::cos(lode_angle);
        const double cos_3theta = std::cos(3.0 * lode_angle);
        const double tan_theta = std::tan(lode_angle);
        const double tan_3theta = std::tan(3.0 * lode_angle);
        const double Root3 = std::sqrt(3.0);
 
        const bool has_symmetric_yield_stress = r_material_properties.Has(YIELD_STRESS);
        const double compr_yield = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_COMPRESSION];
        const double tensi_yield = has_symmetric_yield_stress ? r_material_properties[YIELD_STRESS] : r_material_properties[YIELD_STRESS_TENSION];
        const double n = compr_yield / tensi_yield;
 
        const double angle_phi = (Globals::Pi * 0.25) + dilatancy * 0.5;
        const double alpha = n / (std::tan(angle_phi) * std::tan(angle_phi));
 
        const double CFL = 2.0 * std::tan(angle_phi) / cos_dil;
 
        const double K1 = 0.5 * (1 + alpha) - 0.5 * (1 - alpha) * sin_dil;
        const double K2 = 0.5 * (1 + alpha) - 0.5 * (1 - alpha) / sin_dil;
        const double K3 = 0.5 * (1 + alpha) * sin_dil - 0.5 * (1 - alpha);
 
        double c1, c2, c3;
        if (std::abs(sin_dil) > tolerance)
            c1 = CFL * K3 / 3.0;
        else
            c1 = 0.0; // check
 
        if (checker < 29.0) {
            c2 = cos_theta * CFL * (K1 * (1 + tan_theta * tan_3theta) + K2 * sin_dil * (tan_3theta - tan_theta) / Root3);
            c3 = CFL * (K1 * Root3 * sin_theta + K2 * sin_dil * cos_theta) / (2.0 * J2 * cos_3theta);
        } else {
            c3 = 0.0;
            double Aux = 1.0;
            if (std::abs(lode_angle) > tolerance)
                Aux = -1.0;
            c2 = 0.5 * CFL * (K1 * Root3 + Aux * K2 * sin_dil / Root3);
        }
 
        noalias(rGFlux) = c1 * first_vector + c2 * second_vector + c3 * third_vector;
    }
 
    static int Check(const Properties& rMaterialProperties)
    {
        KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(DILATANCY_ANGLE)) << "DILATANCY_ANGLE is not a defined value" << std::endl;
        if (!rMaterialProperties.Has(YIELD_STRESS)) {
            KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(YIELD_STRESS_TENSION)) << "YIELD_STRESS_TENSION is not a defined value" << std::endl;
            KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(YIELD_STRESS_COMPRESSION)) << "YIELD_STRESS_COMPRESSION is not a defined value" << std::endl;
 
            const double yield_compression = rMaterialProperties[YIELD_STRESS_COMPRESSION];
            const double yield_tension = rMaterialProperties[YIELD_STRESS_TENSION];
 
            KRATOS_ERROR_IF(yield_compression < tolerance) << "Yield stress in compression almost zero or negative, include YIELD_STRESS_COMPRESSION in definition";
            KRATOS_ERROR_IF(yield_tension < tolerance) << "Yield stress in tension almost zero or negative, include YIELD_STRESS_TENSION in definition";
        } else {
            const double yield_stress = rMaterialProperties[YIELD_STRESS];
 
            KRATOS_ERROR_IF(yield_stress < tolerance) << "Yield stress almost zero or negative, include YIELD_STRESS in definition";
        }
 
        return 0;
    }
 
 
 
 
 
 
protected:
 
 
 
 
 
 
 
 
private:
 
 
 
 
 
 
 
 
}; // Class GenericYieldSurface
 
 
 
 
 
} // namespace Kratos.