thermal_simo_ju_yield_surface.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
//
 
#pragma once
 
// System includes
 
// Project includes
#include "custom_constitutive/auxiliary_files/yield_surfaces/simo_ju_yield_surface.h"
#include "thermal_von_mises_yield_surface.h"
 
namespace Kratos
{
 
 
 
 
 
template <class TPlasticPotentialType>
class ThermalSimoJuYieldSurface
{
  public:
 
    using PlasticPotentialType = TPlasticPotentialType;
 
    using BaseType = SimoJuYieldSurface<TPlasticPotentialType>;
    using ThermalVonMisesType = ThermalVonMisesYieldSurface<TPlasticPotentialType>;
 
    static constexpr SizeType Dimension = PlasticPotentialType::Dimension;
 
    static constexpr SizeType VoigtSize = PlasticPotentialType::VoigtSize;
 
    KRATOS_CLASS_POINTER_DEFINITION(ThermalSimoJuYieldSurface);
 
    static constexpr double tolerance = std::numeric_limits<double>::epsilon();
 
    using AdvCLutils = AdvancedConstitutiveLawUtilities<VoigtSize>;
 
 
    ThermalSimoJuYieldSurface()
    {
    }
 
    ThermalSimoJuYieldSurface(ThermalSimoJuYieldSurface const &rOther)
    {
    }
 
    ThermalSimoJuYieldSurface &operator=(ThermalSimoJuYieldSurface const &rOther)
    {
        return *this;
    }
 
    virtual ~ThermalSimoJuYieldSurface(){};
 
 
    static void CalculateEquivalentStress(
        const array_1d<double, VoigtSize>& rStressVector,
        const Vector& rStrainVector,
        double& rEquivalentStress,
        ConstitutiveLaw::Parameters& rValues
        )
    {
        BaseType::CalculateEquivalentStress(rStressVector, rStrainVector, rEquivalentStress, rValues);
    }
 
    static void GetInitialUniaxialThreshold(
        ConstitutiveLaw::Parameters& rValues,
        double& rThreshold
        )
    {
        const auto& r_props = rValues.GetMaterialProperties();
 
        double yield_compression, E;
        if (rValues.IsSetShapeFunctionsValues()) { // This is needed since at Initialize level the N are not set yet...
            E = AdvCLutils::GetMaterialPropertyThroughAccessor(YOUNG_MODULUS, rValues);
            yield_compression = r_props.Has(YIELD_STRESS) ? AdvCLutils::GetMaterialPropertyThroughAccessor(YIELD_STRESS, rValues) : AdvCLutils::GetMaterialPropertyThroughAccessor(YIELD_STRESS_COMPRESSION, rValues);
        } else {
            const double ref_temperature = r_props.Has(REFERENCE_TEMPERATURE) ? r_props[REFERENCE_TEMPERATURE] : rValues.GetElementGeometry().GetValue(REFERENCE_TEMPERATURE);
            E = AdvCLutils::GetPropertyFromTemperatureTable(YOUNG_MODULUS, rValues, ref_temperature);
            yield_compression = r_props.Has(YIELD_STRESS) ? AdvCLutils::GetPropertyFromTemperatureTable(YIELD_STRESS, rValues, ref_temperature) : AdvCLutils::GetPropertyFromTemperatureTable(YIELD_STRESS_COMPRESSION, rValues, ref_temperature);
        }
        rThreshold = yield_compression / std::sqrt(E);
    }
 
    static void CalculateDamageParameter(
        ConstitutiveLaw::Parameters& rValues,
        double& rAParameter,
        const double CharacteristicLength
        )
    {
        const auto& 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 yield_compression = r_material_properties.Has(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);
        const double yield_tension   = r_material_properties.Has(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);
        const double n = yield_compression / yield_tension;
 
        if (r_material_properties[SOFTENING_TYPE] == static_cast<int>(SofteningType::Exponential)) {
            rAParameter = 1.0 / (fracture_energy * n * n / (CharacteristicLength * std::pow(yield_compression, 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(yield_compression, 2) / (2.0 * fracture_energy * n * n / 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
        )
    {
        BaseType::CalculatePlasticPotentialDerivative(rStressVector, rDeviator, J2, rDerivativePlasticPotential, rValues);
    }
 
    static void CalculateYieldSurfaceDerivative(
        const array_1d<double, VoigtSize>& rPredictiveStressVector,
        const array_1d<double, VoigtSize>& Deviator,
        const double J2,
        array_1d<double, VoigtSize>& rFFlux,
        ConstitutiveLaw::Parameters& rValues
        )
    {
        KRATOS_ERROR << "Yield surface derivative not defined for SimoJu..." << std::endl;
    }
 
    static int Check(const Properties& rMaterialProperties)
    {
        return BaseType::Check(rMaterialProperties);
    }
 
    static bool IsWorkingWithTensionThreshold()
    {
        return BaseType::IsWorkingWithTensionThreshold();
    }
 
    static double GetScaleFactorTension(const Properties& rMaterialProperties)
    {
        // The reduction or increase of yield due to temperature is proportional in tension and compression so...
        return BaseType::GetScaleFactorTension(rMaterialProperties);
    }
 
 
 
 
 
 
  protected:
 
 
 
 
 
 
 
  private:
 
 
 
 
 
 
 
 
}; // Class SimoJuYieldSurface
 
 
 
 
 
} // namespace Kratos.