MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim > Class Template Reference

Detailed Description

template<int DisplacementDim>
class MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >

Definition at line 104 of file Lubby2.h.

#include <Lubby2.h>

Inheritance diagram for MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >:

Collaboration diagram for MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >:

Classes
struct	MaterialStateVariables

Public Types
using	KelvinVector = MathLib::KelvinVector::KelvinVectorType< DisplacementDim >
using	KelvinMatrix = MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >
using	ResidualVector = Eigen::Matrix< double, JacobianResidualSize, 1 >
using	JacobianMatrix = Eigen::Matrix< double, JacobianResidualSize, JacobianResidualSize, Eigen::RowMajor >
Public Types inherited from MaterialLib::Solids::MechanicsBase< DisplacementDim >
using	KelvinVector = MathLib::KelvinVector::KelvinVectorType< DisplacementDim >
using	KelvinMatrix = MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >

Public Member Functions
std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables >	createMaterialStateVariables () const override
	Lubby2 (NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters, Lubby2MaterialProperties &material_properties)
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
double	getBulkModulus (double const t, ParameterLib::SpatialPosition const &x, KelvinMatrix const *const) 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
Public Member Functions inherited from MaterialLib::Solids::MechanicsBase< DisplacementDim >
virtual std::optional< std::tuple< KelvinVector, std::unique_ptr< MaterialStateVariables >, KelvinMatrix > >	integrateStress (MaterialPropertyLib::VariableArray const &variable_array_prev, MaterialPropertyLib::VariableArray const &variable_array, double const t, ParameterLib::SpatialPosition const &x, double const dt, MaterialStateVariables const &material_state_variables) const =0
virtual std::vector< InternalVariable >	getInternalVariables () const
virtual ConstitutiveModel	getConstitutiveModel () const
	Gets the type of constitutive model. More...
virtual double	getTemperatureRelatedCoefficient (double const, double const, ParameterLib::SpatialPosition const &, double const, double const) const
virtual double	computeFreeEnergyDensity (double const t, ParameterLib::SpatialPosition const &x, double const dt, KelvinVector const &eps, KelvinVector const &sigma, MaterialStateVariables const &material_state_variables) const =0
virtual	~MechanicsBase ()=default

Static Public Attributes
static int const	KelvinVectorSize
static int const	JacobianResidualSize

Private Member Functions
void	calculateResidualBurgers (double const dt, const KelvinVector &strain_curr, const KelvinVector &strain_t, const KelvinVector &stress_curr, const KelvinVector &stress_t, const KelvinVector &strain_Kel_curr, const KelvinVector &strain_Kel_t, const KelvinVector &strain_Max_curr, const KelvinVector &strain_Max_t, ResidualVector &res, detail::LocalLubby2Properties< DisplacementDim > const &properties) const
	Calculates the 18x1 residual vector. More...
void	calculateJacobianBurgers (double const t, ParameterLib::SpatialPosition const &x, double const dt, JacobianMatrix &Jac, double s_eff, const KelvinVector &sig_i, const KelvinVector &eps_K_i, detail::LocalLubby2Properties< DisplacementDim > const &properties) const
	Calculates the 18x18 Jacobian. More...

Private Attributes
NumLib::NewtonRaphsonSolverParameters const	_nonlinear_solver_parameters
Lubby2MaterialProperties	_mp

Member Typedef Documentation

JacobianMatrix

template<int DisplacementDim>

using MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::JacobianMatrix = Eigen::Matrix<double, JacobianResidualSize, JacobianResidualSize, Eigen::RowMajor>

Definition at line 170 of file Lubby2.h.

KelvinMatrix

template<int DisplacementDim>

using MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::KelvinMatrix = MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>

Definition at line 163 of file Lubby2.h.

KelvinVector

template<int DisplacementDim>

using MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::KelvinVector = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>

Definition at line 161 of file Lubby2.h.

ResidualVector

template<int DisplacementDim>

using MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::ResidualVector = Eigen::Matrix<double, JacobianResidualSize, 1>

Definition at line 169 of file Lubby2.h.

Constructor & Destructor Documentation

Lubby2()

template<int DisplacementDim>

MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::Lubby2 ( NumLib::NewtonRaphsonSolverParameters  nonlinear_solver_parameters, Lubby2MaterialProperties &  material_properties  )

inlineexplicit

Definition at line 176 of file Lubby2.h.

        : _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
          _mp(material_properties)
    {
    }

Member Function Documentation

calculateJacobianBurgers()

template<int DisplacementDim>

void MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::calculateJacobianBurgers ( double const  t, ParameterLib::SpatialPosition const &  x, double const  dt, JacobianMatrix &  Jac, double  s_eff, const KelvinVector &  sig_i, const KelvinVector &  eps_K_i, detail::LocalLubby2Properties< DisplacementDim > const &  properties  ) const

private

Calculates the 18x18 Jacobian.

Definition at line 248 of file Lubby2.cpp.

{
    Jac.setZero();
 
    // build G_11
    Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
        .diagonal()
        .setConstant(1.);
 
    // build G_12
    Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
        .diagonal()
        .setConstant(2.);
 
    // build G_13
    Jac.template block<KelvinVectorSize, KelvinVectorSize>(0,
                                                           2 * KelvinVectorSize)
        .diagonal()
        .setConstant(2.);
 
    // build G_21
    Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
        .noalias() =
        -0.5 * dt * properties.GM0 / properties.etaK * KelvinMatrix::Identity();
    if (s_eff > 0.)
    {
        KelvinVector const eps_K_aid =
            1. / (properties.etaK * properties.etaK) *
            (properties.GM0 * sig_i - 2. * properties.GK * eps_K_i);
 
        KelvinVector const dG_K = 1.5 * _mp.mK(t, x)[0] * properties.GK *
                                  properties.GM0 / s_eff * sig_i;
        KelvinVector const dmu_vK = 1.5 * _mp.mvK(t, x)[0] * properties.GM0 *
                                    properties.etaK / s_eff * sig_i;
        Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
                                                               0)
            .noalias() += 0.5 * dt * eps_K_aid * dmu_vK.transpose() +
                          dt / properties.etaK * eps_K_i * dG_K.transpose();
    }
 
    // build G_22
    Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
                                                           KelvinVectorSize)
        .diagonal()
        .setConstant(1. + dt * properties.GK / properties.etaK);
 
    // nothing to do for G_23
 
    // build G_31
    Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
                                                           0)
        .noalias() =
        -0.5 * dt * properties.GM0 / properties.etaM * KelvinMatrix::Identity();
    if (s_eff > 0.)
    {
        KelvinVector const dmu_vM = 1.5 * _mp.mvM(t, x)[0] * properties.GM0 *
                                    properties.etaM / s_eff * sig_i;
        Jac.template block<KelvinVectorSize, KelvinVectorSize>(
               2 * KelvinVectorSize, 0)
            .noalias() += 0.5 * dt * properties.GM0 /
                          (properties.etaM * properties.etaM) * sig_i *
                          dmu_vM.transpose();
    }
 
    // nothing to do for G_32
 
    // build G_33
    Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
                                                           2 * KelvinVectorSize)
        .diagonal()
        .setConstant(1.);
}

References MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::etaK, MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::etaM, MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::GK, and MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::GM0.

calculateResidualBurgers()

template<int DisplacementDim>

void MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::calculateResidualBurgers ( double const  dt, const KelvinVector &  strain_curr, const KelvinVector &  strain_t, const KelvinVector &  stress_curr, const KelvinVector &  stress_t, const KelvinVector &  strain_Kel_curr, const KelvinVector &  strain_Kel_t, const KelvinVector &  strain_Max_curr, const KelvinVector &  strain_Max_t, ResidualVector &  res, detail::LocalLubby2Properties< DisplacementDim > const &  properties  ) const

private

Calculates the 18x1 residual vector.

Definition at line 215 of file Lubby2.cpp.

{
    // calculate stress residual
    res.template segment<KelvinVectorSize>(0).noalias() =
        (stress_curr - stress_t) -
        2. * ((strain_curr - strain_t) - (strain_Kel_curr - strain_Kel_t) -
              (strain_Max_curr - strain_Max_t));
 
    // calculate Kelvin strain residual
    res.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
        (strain_Kel_curr - strain_Kel_t) -
        dt / (2. * properties.etaK) *
            (properties.GM0 * stress_curr -
             2. * properties.GK * strain_Kel_curr);
 
    // calculate Maxwell strain residual
    res.template segment<KelvinVectorSize>(2 * KelvinVectorSize).noalias() =
        (strain_Max_curr - strain_Max_t) -
        dt * 0.5 * properties.GM0 / properties.etaM * stress_curr;
}

References MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::etaK, MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::etaM, MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::GK, and MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::GM0.

computeFreeEnergyDensity()

template<int DisplacementDim>

double MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::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

inlineoverride

Definition at line 184 of file Lubby2.h.

    {
        assert(dynamic_cast<MaterialStateVariables const*>(
                   &material_state_variables) != nullptr);
        MaterialStateVariables state(static_cast<MaterialStateVariables const&>(
            material_state_variables));
 
        auto const& eps_K = state.eps_K_j;
        auto const& eps_K_prev = state.eps_K_t;
 
        auto const& eps_M = state.eps_M_j;
 
        auto local_lubby2_properties =
            detail::LocalLubby2Properties<DisplacementDim>{t, x, _mp};
 
        // calculation of deviatoric parts
        using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
        auto const& P_dev = Invariants::deviatoric_projection;
        KelvinVector const epsd_i = P_dev * eps;
 
        // initial guess as elastic predictor
        KelvinVector sigd_j = 2.0 * (epsd_i - state.eps_M_t - state.eps_K_t);
 
        // Calculate effective stress and update material properties
        double sig_eff = Invariants::equivalentStress(sigd_j);
        local_lubby2_properties.update(sig_eff);
 
        auto const& eta_K = local_lubby2_properties.etaK;
 
        // This is specific to the backward Euler time scheme and needs to be
        // updated if other time schemes are used.
        return (eps - eps_K - eps_M).dot(sigma) / 2 +
               eps_K.dot(sigma - eta_K * (eps_K - eps_K_prev) / dt) / 2;
    }

References MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::_mp, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_K_j, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_K_t, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_M_j, and MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_M_t.

createMaterialStateVariables()

template<int DisplacementDim>

std::unique_ptr< typename MechanicsBase<DisplacementDim>::MaterialStateVariables> MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::createMaterialStateVariables ( ) const

inlineoverridevirtual

Polymorphic creator for MaterialStateVariables objects specific for a material model.

Reimplemented from MaterialLib::Solids::MechanicsBase< DisplacementDim >.

Definition at line 151 of file Lubby2.h.

    {
        return std::unique_ptr<
            typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
            new MaterialStateVariables};
    }

getBulkModulus()

template<int DisplacementDim>

double MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::getBulkModulus ( double const  t, ParameterLib::SpatialPosition const &  x, KelvinMatrix const * const    ) const

inlineoverridevirtual

Reimplemented from MaterialLib::Solids::MechanicsBase< DisplacementDim >.

Definition at line 226 of file Lubby2.h.

    {
        return _mp.KM0(t, x)[0];
    }

References MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::_mp, and MaterialLib::Solids::Lubby2::Lubby2MaterialProperties::KM0.

integrateStress()

template<int DisplacementDim>

std::optional< std::tuple< typename Lubby2< DisplacementDim >::KelvinVector, std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables >, typename Lubby2< DisplacementDim >::KelvinMatrix > > MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::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

Definition at line 73 of file Lubby2.cpp.

{
    auto const& eps_m = std::get<MPL::SymmetricTensor<DisplacementDim>>(
        variable_array[static_cast<int>(MPL::Variable::mechanical_strain)]);
    auto const& eps_m_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
        variable_array_prev[static_cast<int>(
            MPL::Variable::mechanical_strain)]);
    auto const& sigma_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
        variable_array_prev[static_cast<int>(MPL::Variable::stress)]);
 
    using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
 
    assert(dynamic_cast<MaterialStateVariables const*>(
               &material_state_variables) != nullptr);
    MaterialStateVariables state(
        static_cast<MaterialStateVariables const&>(material_state_variables));
    state.setInitialConditions();
 
    auto local_lubby2_properties =
        detail::LocalLubby2Properties<DisplacementDim>{t, x, _mp};
 
    // calculation of deviatoric parts
    auto const& P_dev = Invariants::deviatoric_projection;
    KelvinVector const eps_m_d_i = P_dev * eps_m;
    KelvinVector const eps_m_d_t = P_dev * eps_m_prev;
 
    // initial guess as elastic predictor.
    KelvinVector sigd_j = 2.0 * (eps_m_d_i - state.eps_M_t - state.eps_K_t);
    // Note: sigd_t contains dimensionless stresses!
    KelvinVector sigd_t = P_dev * sigma_prev / local_lubby2_properties.GM0;
 
    // Calculate effective stress and update material properties
    double sig_eff = Invariants::equivalentStress(sigd_j);
    local_lubby2_properties.update(sig_eff);
 
    using LocalJacobianMatrix =
        Eigen::Matrix<double, KelvinVectorSize * 3, KelvinVectorSize * 3,
                      Eigen::RowMajor>;
 
    // Linear solver for the newton loop is required after the loop with the
    // same matrix. This saves one decomposition.
    Eigen::FullPivLU<LocalJacobianMatrix> linear_solver;
 
    // Different solvers are available for the solution of the local system.
    // TODO Make the following choice of linear solvers available from the
    // input file configuration:
    //      K_loc.partialPivLu().solve(-res_loc);
    //      K_loc.fullPivLu().solve(-res_loc);
    //      K_loc.householderQr().solve(-res_loc);
    //      K_loc.colPivHouseholderQr().solve(res_loc);
    //      K_loc.fullPivHouseholderQr().solve(-res_loc);
    //      K_loc.llt().solve(-res_loc);
    //      K_loc.ldlt().solve(-res_loc);
 
    LocalJacobianMatrix K_loc;
    {  // Local Newton solver
        using LocalResidualVector =
            Eigen::Matrix<double, KelvinVectorSize * 3, 1>;
 
        auto const update_residual = [&](LocalResidualVector& residual)
        {
            calculateResidualBurgers(dt, eps_m_d_i, eps_m_d_t, sigd_j, sigd_t,
                                     state.eps_K_j, state.eps_K_t,
                                     state.eps_M_j, state.eps_M_t, residual,
                                     local_lubby2_properties);
        };
 
        auto const update_jacobian = [&](LocalJacobianMatrix& jacobian)
        {
            calculateJacobianBurgers(
                t, x, dt, jacobian, sig_eff, sigd_j, state.eps_K_j,
                local_lubby2_properties);  // for solution dependent Jacobians
        };
 
        auto const update_solution = [&](LocalResidualVector const& increment)
        {
            // increment solution vectors
            sigd_j.noalias() += increment.template segment<KelvinVectorSize>(
                KelvinVectorSize * 0);
            state.eps_K_j.noalias() +=
                increment.template segment<KelvinVectorSize>(KelvinVectorSize *
                                                             1);
            state.eps_M_j.noalias() +=
                increment.template segment<KelvinVectorSize>(KelvinVectorSize *
                                                             2);
 
            // Calculate effective stress and update material properties
            sig_eff = MathLib::KelvinVector::Invariants<
                KelvinVectorSize>::equivalentStress(sigd_j);
            local_lubby2_properties.update(sig_eff);
        };
 
        auto newton_solver = NumLib::NewtonRaphson<
            decltype(linear_solver), LocalJacobianMatrix,
            decltype(update_jacobian), LocalResidualVector,
            decltype(update_residual), decltype(update_solution)>(
            linear_solver, update_jacobian, update_residual, update_solution,
            _nonlinear_solver_parameters);
 
        auto const success_iterations = newton_solver.solve(K_loc);
 
        if (!success_iterations)
        {
            return {};
        }
 
        // If the Newton loop didn't run, the linear solver will not be
        // initialized.
        // This happens usually for the first iteration of the first timestep.
        if (*success_iterations == 0)
        {
            linear_solver.compute(K_loc);
        }
    }
 
    KelvinMatrix C =
        tangentStiffnessA<DisplacementDim>(local_lubby2_properties.GM0,
                                           local_lubby2_properties.KM0,
                                           linear_solver);
 
    // Hydrostatic part for the stress and the tangent.
    double const delta_eps_m_trace = Invariants::trace(eps_m - eps_m_prev);
    double const sigma_trace_prev = Invariants::trace(sigma_prev);
    KelvinVector const sigma =
        local_lubby2_properties.GM0 * sigd_j +
        (local_lubby2_properties.KM0 * delta_eps_m_trace +
         sigma_trace_prev / 3.) *
            Invariants::identity2;
    return {std::make_tuple(
        sigma,
        std::unique_ptr<
            typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
            new MaterialStateVariables{state}},
        C)};
}

References MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_K_j, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_K_t, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_M_j, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_M_t, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::setInitialConditions(), and NumLib::NewtonRaphson< LinearSolver, JacobianMatrix, JacobianMatrixUpdate, ResidualVector, ResidualUpdate, SolutionUpdate >::solve().

Member Data Documentation

_mp

template<int DisplacementDim>

Lubby2MaterialProperties MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::_mp

private

Definition at line 272 of file Lubby2.h.

Referenced by MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::computeFreeEnergyDensity(), and MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::getBulkModulus().

_nonlinear_solver_parameters

template<int DisplacementDim>

NumLib::NewtonRaphsonSolverParameters const MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::_nonlinear_solver_parameters

private

Definition at line 271 of file Lubby2.h.

JacobianResidualSize

template<int DisplacementDim>

int const MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::JacobianResidualSize

static

Initial value:

=
        3 * KelvinVectorSize

Definition at line 166 of file Lubby2.h.

KelvinVectorSize

template<int DisplacementDim>

int const MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::KelvinVectorSize

static

Initial value:

=
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)

Definition at line 159 of file Lubby2.h.

Referenced by MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::MaterialStateVariables().

---

The documentation for this class was generated from the following files:

• MaterialLib/SolidModels/Lubby2.h
• MaterialLib/SolidModels/Lubby2.cpp