11 #include <boost/math/special_functions/pow.hpp> 53 template <
int DisplacementDim>
57 template <
int DisplacementDim>
89 :
value{1 + gamma_p * theta},
90 pow_m_p{std::pow(
value, m_p)},
91 pow_m_p1{pow_m_p /
value}
100 template <
int DisplacementDim>
103 double const sqrtPhi,
double const alpha_p,
double const beta_p,
104 double const delta_p,
double const epsilon_p)
106 return 3 * (alpha_p * s.
I_1 +
107 4 * boost::math::pow<2>(delta_p) * boost::math::pow<3>(s.
I_1)) /
109 3 * beta_p + 6 * epsilon_p * s.
I_1;
112 template <
int DisplacementDim>
117 double const gamma_p,
double const m_p)
120 (s.
D + s.
J_2 * m_p * gamma_p * dtheta_dsigma / one_gt.
value)) /
124 template <
int DisplacementDim>
129 double const I_1_squared = boost::math::pow<2>(s.
I_1);
136 mp.
alpha / 2. * I_1_squared +
137 boost::math::pow<2>(mp.
delta) *
138 boost::math::pow<2>(I_1_squared)) +
142 template <
int DisplacementDim>
150 double const eps_p_V,
151 double const eps_p_V_dot,
152 double const eps_p_eff_dot,
166 double const theta = s.
J_3 / (s.
J_2 * std::sqrt(s.
J_2));
170 residual.template segment<KelvinVectorSize>(0).noalias() =
171 s.
value / mp.
G - 2 * (eps_D - eps_p_D) -
172 mp.
K / mp.
G * (eps_V - eps_p_V) * identity2;
178 theta * sigma_D_inverse_D - 3. / 2. * theta / s.
J_2 * s.
D;
181 double const sqrtPhi = std::sqrt(
182 s.
J_2 * one_gt.pow_m_p + mp.
alpha_p / 2. * boost::math::pow<2>(s.
I_1) +
183 boost::math::pow<2>(mp.
delta_p) * boost::math::pow<4>(s.
I_1));
185 s, one_gt, sqrtPhi, dtheta_dsigma, mp.
gamma_p, mp.
m_p);
189 eps_p_D_dot - lambda_flow_D;
193 double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
195 residual(2 * KelvinVectorSize, 0) = eps_p_V_dot - lambda * flow_V;
199 residual(2 * KelvinVectorSize + 1) =
201 std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
204 residual(2 * KelvinVectorSize + 2) =
yieldFunction(mp, s, k) / mp.
G;
208 template <
int DisplacementDim>
226 double const theta = s.
J_3 / (s.
J_2 * std::sqrt(s.
J_2));
232 OGS_FATAL(
"Determinant is zero. Matrix is non-invertable.");
236 KelvinVector const sigma_D_inverse_D = P_dev * sigma_D_inverse;
239 theta * sigma_D_inverse_D - 3. / 2. * theta / s.
J_2 * s.
D;
242 double const sqrtPhi = std::sqrt(
243 s.
J_2 * one_gt.pow_m_p + mp.
alpha_p / 2. * boost::math::pow<2>(s.
I_1) +
244 boost::math::pow<2>(mp.
delta_p) * boost::math::pow<4>(s.
I_1));
246 s, one_gt, sqrtPhi, dtheta_dsigma, mp.
gamma_p, mp.
m_p);
253 jacobian.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
254 .noalias() = KelvinMatrix::Identity();
259 .noalias() = 2 * KelvinMatrix::Identity();
263 .noalias() = mp.
K / mp.
G * identity2;
274 one_gt.pow_m_p * (s.
D + gm_p * M0) +
276 4 * boost::math::pow<2>(mp.
delta_p) * boost::math::pow<3>(s.
I_1)) *
280 KelvinMatrix
const M1 =
282 (s.
D * dPhi_dsigma.transpose() + gm_p * M0 * dPhi_dsigma.transpose());
284 KelvinMatrix
const M2 =
285 one_gt.pow_m_p * (P_dev + s.
D * gm_p * M0.transpose());
287 KelvinMatrix
const d2theta_dsigma2 =
288 theta * P_dev * sOdotS<DisplacementDim>(sigma_D_inverse) * P_dev +
289 sigma_D_inverse_D * dtheta_dsigma.transpose() -
290 3. / 2. * theta / s.
J_2 * P_dev -
291 3. / 2. * dtheta_dsigma / s.
J_2 * s.
D.transpose() +
292 3. / 2. * theta / boost::math::pow<2>(s.
J_2) * s.
D * s.
D.transpose();
295 KelvinMatrix
const M3 =
296 gm_p * one_gt.pow_m_p1 *
297 ((s.
D + (gm_p - mp.
gamma_p) * M0) * dtheta_dsigma.transpose() +
298 s.
J_2 * d2theta_dsigma2);
301 KelvinMatrix
const dflow_D_dsigma =
302 (-M1 / (4 * boost::math::pow<3>(sqrtPhi)) + (M2 + M3) / (2 * sqrtPhi)) *
305 .template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
306 .noalias() = -lambda * dflow_D_dsigma;
312 .noalias() = KelvinMatrix::Identity() / dt;
319 2 * KelvinVectorSize + 2)
320 .noalias() = -flow_D;
328 4 * boost::math::pow<2>(mp.
delta_p) *
329 boost::math::pow<3>(s.
I_1)) /
330 (4 * boost::math::pow<3>(sqrtPhi)) * dPhi_dsigma +
332 12 * boost::math::pow<2>(mp.
delta_p * s.
I_1) * identity2) /
337 .noalias() = -lambda * dflow_V_dsigma.transpose();
343 jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize) = 1. / dt;
349 double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
351 jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize + 2) = -flow_V;
355 double const eff_flow =
356 std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
362 -2 / 3. / eff_flow * lambda_flow_D;
365 .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 1, 0)
366 .noalias() = lambda * dflow_D_dsigma * eff_flow23_lambda_flow_D;
368 jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 2) =
369 eff_flow23_lambda_flow_D.transpose() * flow_D;
375 jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 1) = 1. / dt;
379 double const one_gt_pow_m = std::pow(one_gt.value, mp.
m);
380 double const gm = mp.
gamma * mp.
m;
383 mp.
G * (one_gt_pow_m * (s.
D + gm * M0) +
385 4 * boost::math::pow<2>(mp.
delta) *
386 boost::math::pow<3>(s.
I_1)) *
392 .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 2, 0)
393 .noalias() = dF_dsigma.transpose() / mp.
G;
397 jacobian(2 * KelvinVectorSize + 2, 2 * KelvinVectorSize + 1) =
406 template <
int DisplacementDim>
408 double const K,
double const G)
419 return -2. * I * P_dev - 3. * K / G * I * P_sph;
423 double const hardening_coefficient,
424 double const eps_p_eff)
426 return kappa * (1. + eps_p_eff * hardening_coefficient);
429 template <
int DisplacementDim>
431 double const G,
double const K,
447 double const pressure = pressure_prev - K / G * (eps_V - e_prev);
450 P_dev * sigma_prev / G;
453 sigma_D_prev + 2 * P_dev * (eps - eps_prev);
459 template <
typename Res
idualVector,
typename KelvinVector>
460 std::tuple<KelvinVector, PlasticStrain<KelvinVector>,
double>
463 static auto const size = KelvinVector::SizeAtCompileTime;
464 return std::forward_as_tuple(
465 solution.template segment<size>(size * 0),
467 solution[size * 2], solution[size * 2 + 1]},
468 solution[size * 2 + 2]);
471 template <
int DisplacementDim>
475 std::unique_ptr<DamagePropertiesParameters>&& damage_properties,
477 : _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
478 _mp(std::move(material_properties)),
479 _damage_properties(std::move(damage_properties)),
480 _tangent_type(tangent_type)
484 template <
int DisplacementDim>
492 material_state_variables)
const 495 &material_state_variables) !=
nullptr);
498 material_state_variables)
501 auto const& identity2 = Invariants::identity2;
502 return (eps - eps_p.D - eps_p.V / 3 * identity2).dot(sigma) / 2;
505 template <
int DisplacementDim>
506 boost::optional<std::tuple<typename SolidEhlers<DisplacementDim>::KelvinVector,
515 material_state_variables,
519 &material_state_variables) !=
nullptr);
523 material_state_variables);
529 double const eps_V = Invariants::trace(eps);
531 auto const& P_dev = Invariants::deviatoric_projection;
538 KelvinVector sigma = predict_sigma<DisplacementDim>(mp.
G, mp.
K, sigma_prev,
539 eps, eps_prev, eps_V);
546 if ((sigma.squaredNorm() == 0 ||
550 state.
eps_p.eff)) < 0))
552 tangentStiffness = elasticTangentStiffness<DisplacementDim>(
553 mp.
K - 2. / 3 * mp.
G, mp.
G);
566 DisplacementDim>::value;
570 Eigen::Matrix<double, JacobianResidualSize, 1>;
584 auto const& eps_p_D =
585 solution.template segment<KelvinVectorSize>(
590 double const& eps_p_V = solution[KelvinVectorSize * 2];
591 double const eps_p_V_dot =
594 double const& eps_p_eff = solution[KelvinVectorSize * 2 + 1];
595 double const eps_p_eff_dot =
600 solution[KelvinVectorSize * 2 + 1]);
601 residual = calculatePlasticResidual<DisplacementDim>(
603 solution.template segment<KelvinVectorSize>(
605 eps_p_D_dot, solution[KelvinVectorSize * 2], eps_p_V_dot,
606 eps_p_eff_dot, solution[KelvinVectorSize * 2 + 2],
611 jacobian = calculatePlasticJacobian<DisplacementDim>(
612 dt, s, solution[KelvinVectorSize * 2 + 2], mp);
615 auto const update_solution =
617 solution += increment;
619 mp.
G * solution.template segment<KelvinVectorSize>(0)};
625 decltype(update_residual), decltype(update_solution)>(
626 linear_solver, update_jacobian, update_residual,
629 auto const success_iterations = newton_solver.solve(jacobian);
631 if (!success_iterations)
638 if (*success_iterations == 0)
639 linear_solver.compute(jacobian);
641 std::tie(sigma, state.
eps_p, std::ignore) =
642 splitSolutionVector<ResidualVectorType, KelvinVector>(solution);
650 Eigen::RowMajor>::Zero();
651 dresidual_deps.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
652 .noalias() = calculateDResidualDEps<DisplacementDim>(mp.
K, mp.
G);
657 elasticTangentStiffness<DisplacementDim>(mp.
K, mp.
G);
664 linear_solver.solve(-dresidual_deps)
665 .template block<KelvinVectorSize, KelvinVectorSize>(0, 0);
668 tangentStiffness *= 1 - state.
damage.value();
674 "Unimplemented tangent type behaviour for the tangent type " 682 return {std::make_tuple(
691 template <
int DisplacementDim>
692 std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
696 {
"damage.kappa_d", 1,
698 DisplacementDim>::MaterialStateVariables
const& state,
699 std::vector<double>& cache) -> std::vector<double>
const& {
702 auto const& ehlers_state =
711 DisplacementDim>::MaterialStateVariables
const& state,
712 std::vector<double>& cache) -> std::vector<double>
const& {
715 auto const& ehlers_state =
722 {
"eps_p.D", KelvinVector::RowsAtCompileTime,
724 DisplacementDim>::MaterialStateVariables
const& state,
725 std::vector<double>& cache) -> std::vector<double>
const& {
728 auto const& ehlers_state =
731 cache.resize(KelvinVector::RowsAtCompileTime);
732 MathLib::toVector<KelvinVector>(cache,
733 KelvinVector::RowsAtCompileTime) =
735 ehlers_state.eps_p.D);
741 DisplacementDim>::MaterialStateVariables
const& state,
742 std::vector<double>& cache) -> std::vector<double>
const& {
745 auto const& ehlers_state =
749 cache.front() = ehlers_state.
eps_p.V;
754 DisplacementDim>::MaterialStateVariables
const& state,
755 std::vector<double>& cache) -> std::vector<double>
const& {
758 auto const& ehlers_state =
762 cache.front() = ehlers_state.
eps_p.eff;
776 result(0, 0) = v(0) * v(0);
777 result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
778 result(0, 2) = result(2, 0) = v(5) * v(5) / 2.;
779 result(0, 3) = result(3, 0) = v(0) * v(3);
780 result(0, 4) = result(4, 0) = v(3) * v(5) / std::sqrt(2.);
781 result(0, 5) = result(5, 0) = v(0) * v(5);
783 result(1, 1) = v(1) * v(1);
784 result(1, 2) = result(2, 1) = v(4) * v(4) / 2.;
785 result(1, 3) = result(3, 1) = v(3) * v(1);
786 result(1, 4) = result(4, 1) = v(1) * v(4);
787 result(1, 5) = result(5, 1) = v(3) * v(4) / std::sqrt(2.);
789 result(2, 2) = v(2) * v(2);
790 result(2, 3) = result(3, 2) = v(5) * v(4) / std::sqrt(2.);
791 result(2, 4) = result(4, 2) = v(4) * v(2);
792 result(2, 5) = result(5, 2) = v(5) * v(2);
794 result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
795 result(3, 4) = result(4, 3) =
796 v(3) * v(4) / 2. + v(5) * v(1) / std::sqrt(2.);
797 result(3, 5) = result(5, 3) =
798 v(0) * v(4) / std::sqrt(2.) + v(3) * v(5) / 2.;
800 result(4, 4) = v(1) * v(2) + v(4) * v(4) / 2.;
801 result(4, 5) = result(5, 4) =
802 v(3) * v(2) / std::sqrt(2.) + v(5) * v(4) / 2.;
804 result(5, 5) = v(0) * v(2) + v(5) * v(5) / 2.;
814 result(0, 0) = v(0) * v(0);
815 result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
816 result(0, 2) = result(2, 0) = 0;
817 result(0, 3) = result(3, 0) = v(0) * v(3);
819 result(1, 1) = v(1) * v(1);
820 result(1, 2) = result(2, 1) = 0;
821 result(1, 3) = result(3, 1) = v(3) * v(1);
823 result(2, 2) = v(2) * v(2);
824 result(2, 3) = result(3, 2) = 0;
826 result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
double calculateIsotropicHardening(double const kappa, double const hardening_coefficient, double const eps_p_eff)
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Kelvin vector dimensions for given displacement dimension.
MathLib::KelvinVector::KelvinMatrixType< 3 > sOdotS< 3 >(MathLib::KelvinVector::KelvinVectorType< 3 > const &v)
Damage damage
damage part of the state.
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > calculateDResidualDEps(double const K, double const G)
NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters
Eigen::Matrix< double, JacobianResidualSize, 1 > ResidualVectorType
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
static int const KelvinVectorSize
PlasticStrain< KelvinVector > eps_p
plastic part of the state.
Eigen::Matrix< double, JacobianResidualSize, JacobianResidualSize, Eigen::RowMajor > JacobianMatrix
MathLib::KelvinVector::KelvinMatrixType< 2 > sOdotS< 2 >(MathLib::KelvinVector::KelvinVectorType< 2 > const &v)
std::vector< typename MechanicsBase< DisplacementDim >::InternalVariable > getInternalVariables() const override
SolidEhlers(NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters, MaterialPropertiesParameters material_properties, std::unique_ptr< DamagePropertiesParameters > &&damage_properties, TangentType tangent_type)
SolidEhlers< DisplacementDim >::JacobianMatrix calculatePlasticJacobian(double const dt, PhysicalStressWithInvariants< DisplacementDim > const &s, double const lambda, MaterialProperties const &mp)
double const hardening_coefficient
static int const KelvinVectorSize
Holds powers of 1 + gamma_p*theta to base 0, m_p, and m_p-1.
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
static Eigen::Matrix< double, KelvinVectorSize, KelvinVectorSize > const deviatoric_projection
void setInitialConditions()
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
MathLib::KelvinVector::Invariants< KelvinVectorSize > Invariants
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.
MaterialPropertiesParameters _mp
PhysicalStressWithInvariants(KelvinVector const &stress)
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > sOdotS(MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &v)
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)
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
double yieldFunction(MaterialProperties const &mp, PhysicalStressWithInvariants< DisplacementDim > const &s, double const k)
double computeFreeEnergyDensity(double const, ProcessLib::SpatialPosition const &, double const, KelvinVector const &eps, KelvinVector const &sigma, typename MechanicsBase< DisplacementDim >::MaterialStateVariables const &material_state_variables) const override
std::tuple< KelvinVector, PlasticStrain< KelvinVector >, double > splitSolutionVector(ResidualVector const &solution)
static int const JacobianResidualSize
OnePlusGamma_pTheta(double const gamma_p, double const theta, double const m_p)
static double J3(Eigen::Matrix< double, KelvinVectorSize, 1 > const &deviatoric_v)
static double J2(Eigen::Matrix< double, KelvinVectorSize, 1 > const &deviatoric_v)
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)
#define OGS_FATAL(fmt,...)
static double determinant(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Determinant of a matrix in Kelvin vector representation.
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
PlasticStrain< KelvinVector > eps_p_prev
plastic part of the state.
boost::optional< std::tuple< KelvinVector, std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables >, KelvinMatrix > > integrateStress(double const t, ProcessLib::SpatialPosition const &x, double const dt, KelvinVector const &eps_prev, KelvinVector const &eps, KelvinVector const &sigma_prev, typename MechanicsBase< DisplacementDim >::MaterialStateVariables const &material_state_variables, double const T) const override
TangentType const _tangent_type
static Eigen::Matrix< double, KelvinVectorSize, KelvinVectorSize > const spherical_projection
Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > inverse(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
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)
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)
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, 1, Eigen::ColMajor > KelvinVectorType