Revision d8b826d344caa66abbcc3d86852ea99e6a14a329 authored by Jörg Buchwald on 08 March 2021, 20:38:20 UTC, committed by Dmitry Yu. Naumov on 14 March 2021, 16:08:33 UTC
1 parent 883db3d
Ehlers.cpp
/**
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
#include "Ehlers.h"
#include <boost/math/special_functions/pow.hpp>
#include "LinearElasticIsotropic.h"
#include "MaterialLib/MPL/Utils/GetSymmetricTensor.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
namespace MPL = MaterialPropertyLib;
/**
* Common convenitions for naming:
* x_D - deviatoric part of tensor x
* x_V - volumetric part of tensor x
* x_p - a variable related to plastic potential
* x_prev - value of x in previous time step
*
* Variables used in the code:
* eps_D - deviatoric strain
* eps_p_D_dot - deviatoric increment of plastic strain
* eps_p_eff_dot - increment of effective plastic strain
* eps_p_V_dot - volumetric increment of plastic strain
* sigma_D_inverse_D - deviatoric part of sigma_D_inverse
*
* derivation of the flow rule
* theta - J3 / J2^(3 / 2) from yield function
* dtheta_dsigma - derivative of theta
* sqrtPhi - square root of Phi from plastic potential
* flow_D - deviatoric part of flow
* flow_V - volumetric part of flow
* lambda_flow_D - deviatoric increment of plastic strain
*/
namespace MaterialLib
{
namespace Solids
{
namespace Ehlers
{
template <int DisplacementDim>
double StateVariables<DisplacementDim>::getEquivalentPlasticStrain() const
{
return std::sqrt(2.0 / 3.0 * Invariants::FrobeniusNorm(eps_p.D.eval()));
}
/// Special product of \c v with itself: \f$v \odot v\f$.
/// The tensor \c v is given in Kelvin mapping.
/// \note Implementation only for 2 and 3 dimensions.
/// \attention Pay attention to the sign of the result, which normally would be
/// negative, but the returned value is not negated. This has to do with \f$
/// d(A^{-1})/dA = -A^{-1} \odot A^{-1} \f$.
template <int DisplacementDim>
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> sOdotS(
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& v);
template <int DisplacementDim>
struct PhysicalStressWithInvariants final
{
static int const KelvinVectorSize =
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
explicit PhysicalStressWithInvariants(KelvinVector const& stress)
: value{stress},
D{Invariants::deviatoric_projection * stress},
I_1{Invariants::trace(stress)},
J_2{Invariants::J2(D)},
J_3{Invariants::J3(D)}
{
}
KelvinVector value;
KelvinVector D;
double I_1;
double J_2;
double J_3;
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
/// Holds powers of 1 + gamma_p*theta to base 0, m_p, and m_p-1.
struct OnePlusGamma_pTheta final
{
OnePlusGamma_pTheta(double const gamma_p, double const theta,
double const m_p)
: value{1 + gamma_p * theta},
pow_m_p{std::pow(value, m_p)},
pow_m_p1{pow_m_p / value}
{
}
double const value;
double const pow_m_p;
double const pow_m_p1;
};
template <int DisplacementDim>
double plasticFlowVolumetricPart(
PhysicalStressWithInvariants<DisplacementDim> const& s,
double const sqrtPhi, double const alpha_p, double const beta_p,
double const delta_p, double const epsilon_p)
{
return 3 * (alpha_p * s.I_1 +
4 * boost::math::pow<2>(delta_p) * boost::math::pow<3>(s.I_1)) /
(2 * sqrtPhi) +
3 * beta_p + 6 * epsilon_p * s.I_1;
}
template <int DisplacementDim>
typename SolidEhlers<DisplacementDim>::KelvinVector plasticFlowDeviatoricPart(
PhysicalStressWithInvariants<DisplacementDim> const& s,
OnePlusGamma_pTheta const& one_gt, double const sqrtPhi,
typename SolidEhlers<DisplacementDim>::KelvinVector const& dtheta_dsigma,
double const gamma_p, double const m_p)
{
return (one_gt.pow_m_p *
(s.D + s.J_2 * m_p * gamma_p * dtheta_dsigma / one_gt.value)) /
(2 * sqrtPhi);
}
template <int DisplacementDim>
double yieldFunction(MaterialProperties const& mp,
PhysicalStressWithInvariants<DisplacementDim> const& s,
double const k)
{
double const I_1_squared = boost::math::pow<2>(s.I_1);
assert(s.J_2 != 0);
return std::sqrt(
s.J_2 *
std::pow(1 + mp.gamma * s.J_3 / (s.J_2 * std::sqrt(s.J_2)),
mp.m) +
mp.alpha / 2. * I_1_squared +
boost::math::pow<2>(mp.delta) *
boost::math::pow<2>(I_1_squared)) +
mp.beta * s.I_1 + mp.epsilon * I_1_squared - k;
}
template <int DisplacementDim>
typename SolidEhlers<DisplacementDim>::ResidualVectorType
calculatePlasticResidual(
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps_D,
double const eps_V,
PhysicalStressWithInvariants<DisplacementDim> const& s,
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps_p_D,
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps_p_D_dot,
double const eps_p_V,
double const eps_p_V_dot,
double const eps_p_eff_dot,
double const lambda,
double const k,
MaterialProperties const& mp)
{
static int const KelvinVectorSize =
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
auto const& P_dev = Invariants::deviatoric_projection;
auto const& identity2 = Invariants::identity2;
double const theta = s.J_3 / (s.J_2 * std::sqrt(s.J_2));
typename SolidEhlers<DisplacementDim>::ResidualVectorType residual;
// calculate stress residual
residual.template segment<KelvinVectorSize>(0).noalias() =
s.value / mp.G - 2 * (eps_D - eps_p_D) -
mp.K / mp.G * (eps_V - eps_p_V) * identity2;
// deviatoric plastic strain
KelvinVector const sigma_D_inverse_D =
P_dev * MathLib::KelvinVector::inverse(s.D);
KelvinVector const dtheta_dsigma =
theta * sigma_D_inverse_D - 3. / 2. * theta / s.J_2 * s.D;
OnePlusGamma_pTheta const one_gt{mp.gamma_p, theta, mp.m_p};
double const sqrtPhi = std::sqrt(
s.J_2 * one_gt.pow_m_p + mp.alpha_p / 2. * boost::math::pow<2>(s.I_1) +
boost::math::pow<2>(mp.delta_p) * boost::math::pow<4>(s.I_1));
KelvinVector const flow_D = plasticFlowDeviatoricPart(
s, one_gt, sqrtPhi, dtheta_dsigma, mp.gamma_p, mp.m_p);
KelvinVector const lambda_flow_D = lambda * flow_D;
residual.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
eps_p_D_dot - lambda_flow_D;
// plastic volume strain
{
double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
s, sqrtPhi, mp.alpha_p, mp.beta_p, mp.delta_p, mp.epsilon_p);
residual(2 * KelvinVectorSize, 0) = eps_p_V_dot - lambda * flow_V;
}
// evolution of plastic equivalent strain
residual(2 * KelvinVectorSize + 1) =
eps_p_eff_dot -
std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
// yield function (for plastic multiplier)
residual(2 * KelvinVectorSize + 2) = yieldFunction(mp, s, k) / mp.G;
return residual;
}
template <int DisplacementDim>
typename SolidEhlers<DisplacementDim>::JacobianMatrix calculatePlasticJacobian(
double const dt,
PhysicalStressWithInvariants<DisplacementDim> const& s,
double const lambda,
MaterialProperties const& mp)
{
static int const KelvinVectorSize =
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
using KelvinMatrix =
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
auto const& P_dev = Invariants::deviatoric_projection;
auto const& identity2 = Invariants::identity2;
double const theta = s.J_3 / (s.J_2 * std::sqrt(s.J_2));
OnePlusGamma_pTheta const one_gt{mp.gamma_p, theta, mp.m_p};
// inverse of deviatoric stress tensor
if (Invariants::determinant(s.D) == 0)
{
OGS_FATAL("Determinant is zero. Matrix is non-invertable.");
}
// inverse of sigma_D
KelvinVector const sigma_D_inverse = MathLib::KelvinVector::inverse(s.D);
KelvinVector const sigma_D_inverse_D = P_dev * sigma_D_inverse;
KelvinVector const dtheta_dsigma =
theta * sigma_D_inverse_D - 3. / 2. * theta / s.J_2 * s.D;
// deviatoric flow
double const sqrtPhi = std::sqrt(
s.J_2 * one_gt.pow_m_p + mp.alpha_p / 2. * boost::math::pow<2>(s.I_1) +
boost::math::pow<2>(mp.delta_p) * boost::math::pow<4>(s.I_1));
KelvinVector const flow_D = plasticFlowDeviatoricPart(
s, one_gt, sqrtPhi, dtheta_dsigma, mp.gamma_p, mp.m_p);
KelvinVector const lambda_flow_D = lambda * flow_D;
typename SolidEhlers<DisplacementDim>::JacobianMatrix jacobian =
SolidEhlers<DisplacementDim>::JacobianMatrix::Zero();
// G_11
jacobian.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
.noalias() = KelvinMatrix::Identity();
// G_12
jacobian
.template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
.noalias() = 2 * KelvinMatrix::Identity();
// G_13
jacobian.template block<KelvinVectorSize, 1>(0, 2 * KelvinVectorSize)
.noalias() = mp.K / mp.G * identity2;
// G_14 and G_15 are zero
// G_21 -- derivative of deviatoric flow
double const gm_p = mp.gamma_p * mp.m_p;
// intermediate variable for derivative of deviatoric flow
KelvinVector const M0 = s.J_2 / one_gt.value * dtheta_dsigma;
// derivative of Phi w.r.t. sigma
KelvinVector const dPhi_dsigma =
one_gt.pow_m_p * (s.D + gm_p * M0) +
(mp.alpha_p * s.I_1 +
4 * boost::math::pow<2>(mp.delta_p) * boost::math::pow<3>(s.I_1)) *
identity2;
// intermediate variable for derivative of deviatoric flow
KelvinMatrix const M1 =
one_gt.pow_m_p *
(s.D * dPhi_dsigma.transpose() + gm_p * M0 * dPhi_dsigma.transpose());
// intermediate variable for derivative of deviatoric flow
KelvinMatrix const M2 =
one_gt.pow_m_p * (P_dev + s.D * gm_p * M0.transpose());
// second derivative of theta
KelvinMatrix const d2theta_dsigma2 =
theta * P_dev * sOdotS<DisplacementDim>(sigma_D_inverse) * P_dev +
sigma_D_inverse_D * dtheta_dsigma.transpose() -
3. / 2. * theta / s.J_2 * P_dev -
3. / 2. * dtheta_dsigma / s.J_2 * s.D.transpose() +
3. / 2. * theta / boost::math::pow<2>(s.J_2) * s.D * s.D.transpose();
// intermediate variable for derivative of deviatoric flow
KelvinMatrix const M3 =
gm_p * one_gt.pow_m_p1 *
((s.D + (gm_p - mp.gamma_p) * M0) * dtheta_dsigma.transpose() +
s.J_2 * d2theta_dsigma2);
// derivative of flow_D w.r.t. sigma
KelvinMatrix const dflow_D_dsigma =
(-M1 / (4 * boost::math::pow<3>(sqrtPhi)) + (M2 + M3) / (2 * sqrtPhi)) *
mp.G;
jacobian
.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
.noalias() = -lambda * dflow_D_dsigma;
// G_22
jacobian
.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
KelvinVectorSize)
.noalias() = KelvinMatrix::Identity() / dt;
// G_23 and G_24 are zero
// G_25
jacobian
.template block<KelvinVectorSize, 1>(KelvinVectorSize,
2 * KelvinVectorSize + 2)
.noalias() = -flow_D;
// G_31
{
// derivative of flow_V w.r.t. sigma
KelvinVector const dflow_V_dsigma =
3 * mp.G *
(-(mp.alpha_p * s.I_1 +
4 * boost::math::pow<2>(mp.delta_p) *
boost::math::pow<3>(s.I_1)) /
(4 * boost::math::pow<3>(sqrtPhi)) * dPhi_dsigma +
(mp.alpha_p * identity2 +
12 * boost::math::pow<2>(mp.delta_p * s.I_1) * identity2) /
(2 * sqrtPhi) +
2 * mp.epsilon_p * identity2);
jacobian.template block<1, KelvinVectorSize>(2 * KelvinVectorSize, 0)
.noalias() = -lambda * dflow_V_dsigma.transpose();
}
// G_32 is zero
// G_33
jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize) = 1. / dt;
// G_34 is zero
// G_35
{
double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
s, sqrtPhi, mp.alpha_p, mp.beta_p, mp.delta_p, mp.epsilon_p);
jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize + 2) = -flow_V;
}
// increment of effectiv plastic strain
double const eff_flow =
std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
if (eff_flow > 0)
{
// intermediate variable for derivative of plastic jacobian
KelvinVector const eff_flow23_lambda_flow_D =
-2 / 3. / eff_flow * lambda_flow_D;
// G_41
jacobian
.template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 1, 0)
.noalias() = lambda * dflow_D_dsigma * eff_flow23_lambda_flow_D;
// G_45
jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 2) =
eff_flow23_lambda_flow_D.transpose() * flow_D;
}
// G_42 and G_43 are zero
// G_44
jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 1) = 1. / dt;
// G_51
{
double const one_gt_pow_m = std::pow(one_gt.value, mp.m);
double const gm = mp.gamma * mp.m;
// derivative of yield function w.r.t. sigma
KelvinVector const dF_dsigma =
mp.G * (one_gt_pow_m * (s.D + gm * M0) +
(mp.alpha * s.I_1 +
4 * boost::math::pow<2>(mp.delta) *
boost::math::pow<3>(s.I_1)) *
identity2) /
(2. * sqrtPhi) +
mp.G * (mp.beta + 2 * mp.epsilon_p * s.I_1) * identity2;
jacobian
.template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 2, 0)
.noalias() = dF_dsigma.transpose() / mp.G;
}
// G_54
jacobian(2 * KelvinVectorSize + 2, 2 * KelvinVectorSize + 1) =
-mp.kappa * mp.hardening_coefficient / mp.G;
// G_52, G_53, G_55 are zero
return jacobian;
}
/// Calculates the derivative of the residuals with respect to total
/// strain. Implementation fully implicit only.
template <int DisplacementDim>
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> calculateDResidualDEps(
double const K, double const G)
{
static int const KelvinVectorSize =
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
auto const& P_dev = Invariants::deviatoric_projection;
auto const& P_sph = Invariants::spherical_projection;
auto const& I =
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>::Identity();
return -2. * I * P_dev - 3. * K / G * I * P_sph;
}
inline double calculateIsotropicHardening(double const kappa,
double const hardening_coefficient,
double const eps_p_eff)
{
return kappa * (1. + eps_p_eff * hardening_coefficient);
}
template <int DisplacementDim>
typename SolidEhlers<DisplacementDim>::KelvinVector predict_sigma(
double const G, double const K,
typename SolidEhlers<DisplacementDim>::KelvinVector const& sigma_prev,
typename SolidEhlers<DisplacementDim>::KelvinVector const& eps,
typename SolidEhlers<DisplacementDim>::KelvinVector const& eps_prev,
double const eps_V)
{
static int const KelvinVectorSize =
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
auto const& P_dev = Invariants::deviatoric_projection;
// dimensionless initial hydrostatic pressure
double const pressure_prev = Invariants::trace(sigma_prev) / (-3. * G);
// initial strain invariant
double const e_prev = Invariants::trace(eps_prev);
// dimensioness hydrostatic stress increment
double const pressure = pressure_prev - K / G * (eps_V - e_prev);
// dimensionless deviatoric initial stress
typename SolidEhlers<DisplacementDim>::KelvinVector const sigma_D_prev =
P_dev * sigma_prev / G;
// dimensionless deviatoric stress
typename SolidEhlers<DisplacementDim>::KelvinVector const sigma_D =
sigma_D_prev + 2 * P_dev * (eps - eps_prev);
return sigma_D - pressure * Invariants::identity2;
}
/// Split the agglomerated solution vector in separate items. The arrangement
/// must be the same as in the newton() function.
template <typename ResidualVector, typename KelvinVector>
std::tuple<KelvinVector, PlasticStrain<KelvinVector>, double>
splitSolutionVector(ResidualVector const& solution)
{
static auto const size = KelvinVector::SizeAtCompileTime;
return std::forward_as_tuple(
solution.template segment<size>(size * 0),
PlasticStrain<KelvinVector>{solution.template segment<size>(size * 1),
solution[size * 2], solution[size * 2 + 1]},
solution[size * 2 + 2]);
}
template <int DisplacementDim>
SolidEhlers<DisplacementDim>::SolidEhlers(
NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters,
MaterialPropertiesParameters material_properties,
std::unique_ptr<DamagePropertiesParameters>&& damage_properties,
TangentType tangent_type)
: _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
_mp(std::move(material_properties)),
_damage_properties(std::move(damage_properties)),
_tangent_type(tangent_type)
{
}
template <int DisplacementDim>
double SolidEhlers<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
{
assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
&material_state_variables) != nullptr);
auto const& eps_p = static_cast<StateVariables<DisplacementDim> const&>(
material_state_variables)
.eps_p;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
auto const& identity2 = Invariants::identity2;
return (eps - eps_p.D - eps_p.V / 3 * identity2).dot(sigma) / 2;
}
template <int DisplacementDim>
std::optional<std::tuple<typename SolidEhlers<DisplacementDim>::KelvinVector,
std::unique_ptr<typename MechanicsBase<
DisplacementDim>::MaterialStateVariables>,
typename SolidEhlers<DisplacementDim>::KelvinMatrix>>
SolidEhlers<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
{
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)]);
assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
&material_state_variables) != nullptr);
StateVariables<DisplacementDim> state =
static_cast<StateVariables<DisplacementDim> const&>(
material_state_variables);
state.setInitialConditions();
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
// volumetric strain
double const eps_V = Invariants::trace(eps_m);
auto const& P_dev = Invariants::deviatoric_projection;
// deviatoric strain
KelvinVector const eps_D = P_dev * eps_m;
// do the evaluation once per function call.
MaterialProperties const mp(t, x, _mp);
KelvinVector sigma = predict_sigma<DisplacementDim>(
mp.G, mp.K, sigma_prev, eps_m, eps_m_prev, eps_V);
KelvinMatrix tangentStiffness;
PhysicalStressWithInvariants<DisplacementDim> s{mp.G * sigma};
// Quit early if sigma is zero (nothing to do) or if we are still in elastic
// zone.
if ((sigma.squaredNorm() == 0 ||
yieldFunction(
mp, s,
calculateIsotropicHardening(mp.kappa, mp.hardening_coefficient,
state.eps_p.eff)) < 0))
{
tangentStiffness = elasticTangentStiffness<DisplacementDim>(
mp.K - 2. / 3 * mp.G, mp.G);
}
else
{
// Linear solver for the newton loop is required after the loop with the
// same matrix. This saves one decomposition.
Eigen::FullPivLU<Eigen::Matrix<double, JacobianResidualSize,
JacobianResidualSize, Eigen::RowMajor>>
linear_solver;
{
static int const KelvinVectorSize =
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value;
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
using ResidualVectorType =
Eigen::Matrix<double, JacobianResidualSize, 1>;
using JacobianMatrix =
Eigen::Matrix<double, JacobianResidualSize,
JacobianResidualSize, Eigen::RowMajor>;
JacobianMatrix jacobian;
// Agglomerated solution vector construction. It is later split
// into individual parts by splitSolutionVector().
ResidualVectorType solution;
solution << sigma, state.eps_p.D, state.eps_p.V, state.eps_p.eff, 0;
auto const update_residual = [&](ResidualVectorType& residual) {
auto const& eps_p_D =
solution.template segment<KelvinVectorSize>(
KelvinVectorSize);
KelvinVector const eps_p_D_dot =
(eps_p_D - state.eps_p_prev.D) / dt;
double const& eps_p_V = solution[KelvinVectorSize * 2];
double const eps_p_V_dot =
(eps_p_V - state.eps_p_prev.V) / dt;
double const& eps_p_eff = solution[KelvinVectorSize * 2 + 1];
double const eps_p_eff_dot =
(eps_p_eff - state.eps_p_prev.eff) / dt;
double const k_hardening = calculateIsotropicHardening(
mp.kappa, mp.hardening_coefficient,
solution[KelvinVectorSize * 2 + 1]);
residual = calculatePlasticResidual<DisplacementDim>(
eps_D, eps_V, s,
solution.template segment<KelvinVectorSize>(
KelvinVectorSize),
eps_p_D_dot, solution[KelvinVectorSize * 2], eps_p_V_dot,
eps_p_eff_dot, solution[KelvinVectorSize * 2 + 2],
k_hardening, mp);
};
auto const update_jacobian = [&](JacobianMatrix& jacobian) {
jacobian = calculatePlasticJacobian<DisplacementDim>(
dt, s, solution[KelvinVectorSize * 2 + 2], mp);
};
auto const update_solution =
[&](ResidualVectorType const& increment) {
solution += increment;
s = PhysicalStressWithInvariants<DisplacementDim>{
mp.G * solution.template segment<KelvinVectorSize>(0)};
};
auto newton_solver = NumLib::NewtonRaphson<
decltype(linear_solver), JacobianMatrix,
decltype(update_jacobian), ResidualVectorType,
decltype(update_residual), decltype(update_solution)>(
linear_solver, update_jacobian, update_residual,
update_solution, _nonlinear_solver_parameters);
auto const success_iterations = newton_solver.solve(jacobian);
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(jacobian);
}
std::tie(sigma, state.eps_p, std::ignore) =
splitSolutionVector<ResidualVectorType, KelvinVector>(solution);
}
// Calculate residual derivative w.r.t. strain
Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
Eigen::RowMajor>
dresidual_deps =
Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
Eigen::RowMajor>::Zero();
dresidual_deps.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
.noalias() = calculateDResidualDEps<DisplacementDim>(mp.K, mp.G);
if (_tangent_type == TangentType::Elastic)
{
tangentStiffness =
elasticTangentStiffness<DisplacementDim>(mp.K, mp.G);
}
else if (_tangent_type == TangentType::Plastic ||
_tangent_type == TangentType::PlasticDamageSecant)
{
tangentStiffness =
mp.G *
linear_solver.solve(-dresidual_deps)
.template block<KelvinVectorSize, KelvinVectorSize>(0, 0);
if (_tangent_type == TangentType::PlasticDamageSecant)
{
tangentStiffness *= 1 - state.damage.value();
}
}
else
{
OGS_FATAL(
"Unimplemented tangent type behaviour for the tangent type "
"'{:d}'.",
_tangent_type);
}
}
KelvinVector sigma_final = mp.G * sigma;
return {std::make_tuple(
sigma_final,
std::unique_ptr<
typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
new StateVariables<DisplacementDim>{
static_cast<StateVariables<DisplacementDim> const&>(state)}},
tangentStiffness)};
}
template <int DisplacementDim>
std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
SolidEhlers<DisplacementDim>::getInternalVariables() const
{
return {
{"damage.kappa_d", 1,
[](typename MechanicsBase<
DisplacementDim>::MaterialStateVariables const& state,
std::vector<double>& cache) -> std::vector<double> const& {
assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
&state) != nullptr);
auto const& ehlers_state =
static_cast<StateVariables<DisplacementDim> const&>(state);
cache.resize(1);
cache.front() = ehlers_state.damage.kappa_d();
return cache;
}},
{"damage.value", 1,
[](typename MechanicsBase<
DisplacementDim>::MaterialStateVariables const& state,
std::vector<double>& cache) -> std::vector<double> const& {
assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
&state) != nullptr);
auto const& ehlers_state =
static_cast<StateVariables<DisplacementDim> const&>(state);
cache.resize(1);
cache.front() = ehlers_state.damage.value();
return cache;
}},
{"eps_p.D", KelvinVector::RowsAtCompileTime,
[](typename MechanicsBase<
DisplacementDim>::MaterialStateVariables const& state,
std::vector<double>& cache) -> std::vector<double> const& {
assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
&state) != nullptr);
auto const& ehlers_state =
static_cast<StateVariables<DisplacementDim> const&>(state);
cache.resize(KelvinVector::RowsAtCompileTime);
MathLib::toVector<KelvinVector>(cache,
KelvinVector::RowsAtCompileTime) =
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
ehlers_state.eps_p.D);
return cache;
}},
{"eps_p.V", 1,
[](typename MechanicsBase<
DisplacementDim>::MaterialStateVariables const& state,
std::vector<double>& cache) -> std::vector<double> const& {
assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
&state) != nullptr);
auto const& ehlers_state =
static_cast<StateVariables<DisplacementDim> const&>(state);
cache.resize(1);
cache.front() = ehlers_state.eps_p.V;
return cache;
}},
{"eps_p.eff", 1,
[](typename MechanicsBase<
DisplacementDim>::MaterialStateVariables const& state,
std::vector<double>& cache) -> std::vector<double> const& {
assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
&state) != nullptr);
auto const& ehlers_state =
static_cast<StateVariables<DisplacementDim> const&>(state);
cache.resize(1);
cache.front() = ehlers_state.eps_p.eff;
return cache;
}}};
}
template class SolidEhlers<2>;
template class SolidEhlers<3>;
template <>
MathLib::KelvinVector::KelvinMatrixType<3> sOdotS<3>(
MathLib::KelvinVector::KelvinVectorType<3> const& v)
{
MathLib::KelvinVector::KelvinMatrixType<3> result;
result(0, 0) = v(0) * v(0);
result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
result(0, 2) = result(2, 0) = v(5) * v(5) / 2.;
result(0, 3) = result(3, 0) = v(0) * v(3);
result(0, 4) = result(4, 0) = v(3) * v(5) / std::sqrt(2.);
result(0, 5) = result(5, 0) = v(0) * v(5);
result(1, 1) = v(1) * v(1);
result(1, 2) = result(2, 1) = v(4) * v(4) / 2.;
result(1, 3) = result(3, 1) = v(3) * v(1);
result(1, 4) = result(4, 1) = v(1) * v(4);
result(1, 5) = result(5, 1) = v(3) * v(4) / std::sqrt(2.);
result(2, 2) = v(2) * v(2);
result(2, 3) = result(3, 2) = v(5) * v(4) / std::sqrt(2.);
result(2, 4) = result(4, 2) = v(4) * v(2);
result(2, 5) = result(5, 2) = v(5) * v(2);
result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
result(3, 4) = result(4, 3) =
v(3) * v(4) / 2. + v(5) * v(1) / std::sqrt(2.);
result(3, 5) = result(5, 3) =
v(0) * v(4) / std::sqrt(2.) + v(3) * v(5) / 2.;
result(4, 4) = v(1) * v(2) + v(4) * v(4) / 2.;
result(4, 5) = result(5, 4) =
v(3) * v(2) / std::sqrt(2.) + v(5) * v(4) / 2.;
result(5, 5) = v(0) * v(2) + v(5) * v(5) / 2.;
return result;
}
template <>
MathLib::KelvinVector::KelvinMatrixType<2> sOdotS<2>(
MathLib::KelvinVector::KelvinVectorType<2> const& v)
{
MathLib::KelvinVector::KelvinMatrixType<2> result;
result(0, 0) = v(0) * v(0);
result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
result(0, 2) = result(2, 0) = 0;
result(0, 3) = result(3, 0) = v(0) * v(3);
result(1, 1) = v(1) * v(1);
result(1, 2) = result(2, 1) = 0;
result(1, 3) = result(3, 1) = v(3) * v(1);
result(2, 2) = v(2) * v(2);
result(2, 3) = result(3, 2) = 0;
result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
return result;
}
} // namespace Ehlers
} // namespace Solids
} // namespace MaterialLib
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...