LinearElasticIsotropicPhaseField.h

Go to the documentation of this file.

 
#pragma once
 
#include <Eigen/Eigenvalues>
#include <boost/math/special_functions/pow.hpp>
 
#include "MathLib/KelvinVector.h"
 
namespace MaterialLib
{
namespace Solids
{
namespace Phasefield
{
inline double heaviside(double const v)
{
    return (v < 0) ? 0.0 : 1.0;
}
 
inline double macaulayTensile(double const v)
{
    return v * heaviside(v);
}
inline double macaulayCompressive(double v)
{
    return v * (1 - heaviside(v));
}
 
template <int DisplacementDim>
std::tuple<MathLib::KelvinVector::KelvinVectorType<
               DisplacementDim> /* sigma_real */,
           MathLib::KelvinVector::KelvinVectorType<
               DisplacementDim> /* sigma_tensile */,
           MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> /* D */,
           double /* strain_energy_tensile */, double /* elastic_energy */,
           MathLib::KelvinVector::KelvinMatrixType<
               DisplacementDim> /* C_tensile */,
           MathLib::KelvinVector::KelvinMatrixType<
               DisplacementDim> /* C_compressive
                                 */
           >
calculateVolDevDegradedStress(
    double const degradation, double const bulk_modulus, double const mu,
    MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps)
{
    static int const KelvinVectorSize =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
    using KelvinVector =
        MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
    using KelvinMatrix =
        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
    using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
    // calculation of deviatoric parts
    auto const& P_dev = Invariants::deviatoric_projection;
    KelvinVector const epsd_curr = P_dev * eps;
 
    // Hydrostatic part for the stress and the tangent.
    double const eps_curr_trace = Invariants::trace(eps);
 
    KelvinMatrix C_tensile = KelvinMatrix::Zero();
    KelvinMatrix C_compressive = KelvinMatrix::Zero();
 
    auto strain_energy_computation_vol = [&](auto&& macaulay)
    {
        auto macaulay_squared = [&macaulay](double x)
        { return boost::math::pow<2>(macaulay(x)); };
        return bulk_modulus / 2 * macaulay_squared(eps_curr_trace);
    };
 
    auto stress_computation_vol = [&](auto&& macaulay)
    { return bulk_modulus * macaulay(eps_curr_trace) * Invariants::identity2; };
 
    auto hs = [&](double const v) { return heaviside(v); };
 
    auto mt = [&](double const v) { return macaulayTensile(v); };
 
    auto mc = [&](double const v) { return macaulayCompressive(v); };
 
    double const strain_energy_tensile = strain_energy_computation_vol(mt) +
                                         mu * epsd_curr.transpose() * epsd_curr;
 
    KelvinVector const sigma_tensile =
        stress_computation_vol(mt) + 2 * mu * epsd_curr;
 
    KelvinVector const sigma_compressive = stress_computation_vol(mc);
 
    C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus *
                                                         hs(eps_curr_trace));
    C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
 
    C_compressive.template topLeftCorner<3, 3>().setConstant(
        bulk_modulus * (1 - hs(eps_curr_trace)));
 
    double const elastic_energy =
        bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
        mu * epsd_curr.transpose() * epsd_curr;
    KelvinVector const sigma_real =
        degradation * sigma_tensile + sigma_compressive;
 
    KelvinMatrix const D = degradation * C_tensile + C_compressive;
 
    return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
                           elastic_energy, C_tensile, C_compressive);
}
 
template <int DisplacementDim>
std::tuple<MathLib::KelvinVector::KelvinVectorType<
               DisplacementDim> /* sigma_real */,
           MathLib::KelvinVector::KelvinVectorType<
               DisplacementDim> /* sigma_tensile */,
           MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> /* D */,
           double /* strain_energy_tensile */, double /* elastic_energy */,
           MathLib::KelvinVector::KelvinMatrixType<
               DisplacementDim> /* C_tensile */,
           MathLib::KelvinVector::KelvinMatrixType<
               DisplacementDim> /* C_compressive
                                 */
           >
calculateIsotropicDegradedStress(
    double const degradation,
    double const bulk_modulus,
    double const mu,
    MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps)
{
    static int const KelvinVectorSize =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
    using KelvinVector =
        MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
    using KelvinMatrix =
        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
    using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
    // calculation of deviatoric parts
    auto const& P_dev = Invariants::deviatoric_projection;
    KelvinVector const epsd_curr = P_dev * eps;
 
    // Hydrostatic part for the stress and the tangent.
    double const eps_curr_trace = Invariants::trace(eps);
 
    KelvinMatrix C_tensile = KelvinMatrix::Zero();
    KelvinMatrix C_compressive = KelvinMatrix::Zero();
 
    double const strain_energy_tensile =
        bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
        mu * epsd_curr.transpose() * epsd_curr;
    double const elastic_energy = degradation * strain_energy_tensile;
    KelvinVector const sigma_tensile =
        bulk_modulus * eps_curr_trace * Invariants::identity2 +
        2 * mu * epsd_curr;
    C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
    C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
    KelvinVector const sigma_real = degradation * sigma_tensile;
 
    KelvinMatrix const D = degradation * C_tensile + C_compressive;
 
    return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
                           elastic_energy, C_tensile, C_compressive);
}
 
template <int DisplacementDim>
std::tuple<MathLib::KelvinVector::KelvinVectorType<
               DisplacementDim> /* sigma_real */,
           MathLib::KelvinVector::KelvinVectorType<
               DisplacementDim> /* sigma_tensile */,
           MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> /* D */,
           double /* strain_energy_tensile */, double /* elastic_energy */,
           MathLib::KelvinVector::KelvinMatrixType<
               DisplacementDim> /* C_tensile */,
           MathLib::KelvinVector::KelvinMatrixType<
               DisplacementDim> /* C_compressive
                                 */
           >
calculateIsotropicDegradedStressWithRankineEnergy(
    double const degradation,
    double const bulk_modulus,
    double const mu,
    MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps)
{
    static int const KelvinVectorSize =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
    using KelvinVector =
        MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
    using KelvinMatrix =
        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
    using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
    // calculation of deviatoric parts
    auto const& P_dev = Invariants::deviatoric_projection;
    KelvinVector const epsd_curr = P_dev * eps;
 
    // Hydrostatic part for the stress and the tangent.
    double const eps_curr_trace = Invariants::trace(eps);
 
    KelvinMatrix C_tensile = KelvinMatrix::Zero();
    KelvinMatrix C_compressive = KelvinMatrix::Zero();
 
    KelvinVector const sigma_tensile =
        bulk_modulus * eps_curr_trace * Invariants::identity2 +
        2 * mu * epsd_curr;
    // The maximum principal stress
    auto sigma_tensor =
        MathLib::KelvinVector::kelvinVectorToTensor(sigma_tensile);
    Eigen::SelfAdjointEigenSolver<decltype(sigma_tensor)>
        selfAdjoint_eigen_solver(sigma_tensor);
    double const max_principle_stress =
        selfAdjoint_eigen_solver
            .eigenvalues()[2];  // selfAdjoint_eigen_solver function gives an
                                // ascending sequence.
    //
    double const E0 = 9.0 * mu * bulk_modulus / (3.0 * bulk_modulus + mu);
    auto hs = [&](double const v) { return heaviside(v); };
    //
    double const strain_energy_tensile = 0.5 * max_principle_stress *
                                         max_principle_stress *
                                         hs(max_principle_stress) / E0;
    double const elastic_energy =
        degradation * (bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
                       mu * epsd_curr.transpose() * epsd_curr);
    //
    C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
    C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
    KelvinVector const sigma_real = degradation * sigma_tensile;
 
    KelvinMatrix const D = degradation * C_tensile + C_compressive;
    return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
                           elastic_energy, C_tensile, C_compressive);
}
 
}  // namespace Phasefield
}  // namespace Solids
}  // namespace MaterialLib