Ehlers.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause

 
#pragma once
 
#ifndef NDEBUG
#include <iosfwd>
#endif
 
#include "BaseLib/Error.h"
#include "MathLib/KelvinVector.h"
#include "MechanicsBase.h"
#include "NumLib/NewtonRaphson.h"
#include "ParameterLib/Parameter.h"
 
namespace MaterialLib
{
namespace Solids
{
namespace Ehlers
{
enum class TangentType
{
    Elastic,
    PlasticDamageSecant,
    Plastic
};
inline TangentType makeTangentType(std::string const& s)
{
    if (s == "Elastic")
    {
        return TangentType::Elastic;
    }
    if (s == "PlasticDamageSecant")
    {
        return TangentType::PlasticDamageSecant;
    }
    if (s == "Plastic")
    {
        return TangentType::Plastic;
    }
    OGS_FATAL("Not valid string '{:s}' to create a tangent type from.", s);
}

struct MaterialPropertiesParameters
{
    using P = ParameterLib::Parameter<double>;
 
    MaterialPropertiesParameters(P const& G_, P const& K_, P const& alpha_,
                                 P const& beta_, P const& gamma_,
                                 P const& delta_, P const& epsilon_,
                                 P const& m_, P const& alpha_p_,
                                 P const& beta_p_, P const& gamma_p_,
                                 P const& delta_p_, P const& epsilon_p_,
                                 P const& m_p_, P const& kappa_,
                                 P const& hardening_coefficient_)
        : G(G_),
          K(K_),
          alpha(alpha_),
          beta(beta_),
          gamma(gamma_),
          delta(delta_),
          epsilon(epsilon_),
          m(m_),
          alpha_p(alpha_p_),
          beta_p(beta_p_),
          gamma_p(gamma_p_),
          delta_p(delta_p_),
          epsilon_p(epsilon_p_),
          m_p(m_p_),
          kappa(kappa_),
          hardening_coefficient(hardening_coefficient_)
    {
    }
    // basic material parameters
    P const& G;  
    P const& K;  
 
    P const& alpha;    
    P const& beta;     
    P const& gamma;    
    P const& delta;    
    P const& epsilon;  
    P const& m;        
 
    P const& alpha_p;    
    P const& beta_p;     
    P const& gamma_p;    
    P const& delta_p;    
    P const& epsilon_p;  
    P const& m_p;        
 
    P const& kappa;  
    P const& hardening_coefficient;
};
 
struct DamagePropertiesParameters
{
    using P = ParameterLib::Parameter<double>;
    P const& alpha_d;
    P const& beta_d;
    P const& h_d;
};

struct MaterialProperties final
{
    MaterialProperties(double const t, ParameterLib::SpatialPosition const& x,
                       MaterialPropertiesParameters const& mp)
        : G(mp.G(t, x)[0]),
          K(mp.K(t, x)[0]),
          alpha(mp.alpha(t, x)[0]),
          beta(mp.beta(t, x)[0]),
          gamma(mp.gamma(t, x)[0]),
          delta(mp.delta(t, x)[0]),
          epsilon(mp.epsilon(t, x)[0]),
          m(mp.m(t, x)[0]),
          alpha_p(mp.alpha_p(t, x)[0]),
          beta_p(mp.beta_p(t, x)[0]),
          gamma_p(mp.gamma_p(t, x)[0]),
          delta_p(mp.delta_p(t, x)[0]),
          epsilon_p(mp.epsilon_p(t, x)[0]),
          m_p(mp.m_p(t, x)[0]),
          kappa(mp.kappa(t, x)[0]),
          hardening_coefficient(mp.hardening_coefficient(t, x)[0])
    {
    }
    double const G;
    double const K;
 
    double const alpha;
    double const beta;
    double const gamma;
    double const delta;
    double const epsilon;
    double const m;
 
    double const alpha_p;
    double const beta_p;
    double const gamma_p;
    double const delta_p;
    double const epsilon_p;
    double const m_p;
 
    double const kappa;
    double const hardening_coefficient;
};

struct DamageProperties
{
    DamageProperties(double const t,
                     ParameterLib::SpatialPosition const& x,
                     DamagePropertiesParameters const& dp)
        : alpha_d(dp.alpha_d(t, x)[0]),
          beta_d(dp.beta_d(t, x)[0]),
          h_d(dp.h_d(t, x)[0])
    {
    }
    double const alpha_d;
    double const beta_d;
    double const h_d;
};
 
template <typename KelvinVector>
struct PlasticStrain final
{
    PlasticStrain() : D(KelvinVector::Zero()) {}
    PlasticStrain(KelvinVector eps_p_D_, double const eps_p_V_,
                  double const eps_p_eff_)
        : D(std::move(eps_p_D_)), V(eps_p_V_), eff(eps_p_eff_){};
 
    KelvinVector D;  
    double V = 0;    
    double eff = 0;  
 
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
 
class Damage final
{
public:
    Damage() = default;
    Damage(double const kappa_d, double const value)
        : _kappa_d(kappa_d), _value(value)
    {
    }
 
    double kappa_d() const { return _kappa_d; }
    double& kappa_d() { return _kappa_d; }
    double value() const { return _value; }
    double& value() { return _value; }
 
private:
    double _kappa_d = 0;  
    double _value = 0;    
};
 
template <int DisplacementDim>
struct StateVariables
    : public MechanicsBase<DisplacementDim>::MaterialStateVariables
{
    void setInitialConditions()
    {
        eps_p = eps_p_prev;
        damage = damage_prev;
    }
 
    void pushBackState() override
    {
        eps_p_prev = eps_p;
        damage_prev = damage;
    }
 
    double getEquivalentPlasticStrain() const override;
 
    static int const KelvinVectorSize =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
    using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
    using KelvinVector =
        MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
 
    PlasticStrain<KelvinVector> eps_p;  
    Damage damage;                      
 
    // Initial values from previous timestep
    PlasticStrain<KelvinVector> eps_p_prev;  
    Damage damage_prev;                      
 
#ifndef NDEBUG
    friend std::ostream& operator<<(std::ostream& os,
                                    StateVariables<DisplacementDim> const& m)
    {
        os << "State:\n"
           << "eps_p_D: " << m.eps_p.D << "\n"
           << "eps_p_eff: " << m.eps_p.eff << "\n"
           << "kappa_d: " << m.damage.kappa_d() << "\n"
           << "damage: " << m.damage.value() << "\n"
           << "eps_p_D_prev: " << m.eps_p_prev.D << "\n"
           << "eps_p_eff_prev: " << m.eps_p_prev.eff << "\n"
           << "kappa_d_prev: " << m.damage_prev.kappa_d() << "\n"
           << "damage_prev: " << m.damage_prev.value() << "\n";
        return os;
    }
#endif  // NDEBUG
 
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
 
template <int DisplacementDim>
class SolidEhlers final : public MechanicsBase<DisplacementDim>
{
public:
    static int const KelvinVectorSize =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
    static int const JacobianResidualSize =
        2 * KelvinVectorSize + 3;  // 2 is the number of components in the
                                   // jacobian/residual, not the space
                                   // dimension. And 3 is for additional
                                   // variables.
    using ResidualVectorType = Eigen::Matrix<double, JacobianResidualSize, 1>;
    using JacobianMatrix = Eigen::Matrix<double, JacobianResidualSize,
                                         JacobianResidualSize, Eigen::RowMajor>;
 
public:
    std::unique_ptr<
        typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
    createMaterialStateVariables() const override
    {
        return std::make_unique<StateVariables<DisplacementDim>>();
    }
 
public:
    using KelvinVector =
        MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
    using KelvinMatrix =
        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
 
public:
    SolidEhlers(
        NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters,
        MaterialPropertiesParameters material_properties,
        std::unique_ptr<DamagePropertiesParameters>&& damage_properties,
        TangentType tangent_type);
 
    double computeFreeEnergyDensity(
        double const /*t*/,
        ParameterLib::SpatialPosition const& /*x*/,
        double const /*dt*/,
        KelvinVector const& eps,
        KelvinVector const& sigma,
        typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
            material_state_variables) const override;
 
    std::optional<std::tuple<KelvinVector,
                             std::unique_ptr<typename MechanicsBase<
                                 DisplacementDim>::MaterialStateVariables>,
                             KelvinMatrix>>
    integrateStress(
        MaterialPropertyLib::VariableArray const& variable_array_prev,
        MaterialPropertyLib::VariableArray const& variable_array,
        double const t, ParameterLib::SpatialPosition const& x, double const dt,
        typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
            material_state_variables) const override;
 
    std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
    getInternalVariables() const override;
 
    MaterialProperties evaluatedMaterialProperties(
        double const t, ParameterLib::SpatialPosition const& x) const
    {
        return MaterialProperties(t, x, _mp);
    }
 
    double getBulkModulus(double const t,
                          ParameterLib::SpatialPosition const& x,
                          KelvinMatrix const* const /*C*/) const override
    {
        return _mp.K(t, x)[0];
    }
 
    DamageProperties evaluatedDamageProperties(
        double const t, ParameterLib::SpatialPosition const& x) const
    {
        return DamageProperties(t, x, *_damage_properties);
    }
 
private:
    NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters;
 
    MaterialPropertiesParameters _mp;
    std::unique_ptr<DamagePropertiesParameters> _damage_properties;
    TangentType const _tangent_type;
};
 
extern template class SolidEhlers<2>;
extern template class SolidEhlers<3>;
 
extern template struct StateVariables<2>;
extern template struct StateVariables<3>;
}  // namespace Ehlers
}  // namespace Solids
}  // namespace MaterialLib