13#include <boost/math/special_functions/pow.hpp>
51template <
int DisplacementDim>
54 return std::sqrt(2.0 / 3.0 * Invariants::FrobeniusNorm(eps_p.D.eval()));
63template <
int DisplacementDim>
67template <
int DisplacementDim>
99 :
value{1 + gamma_p * theta},
110template <
int DisplacementDim>
113 double const sqrtPhi,
double const alpha_p,
double const beta_p,
114 double const delta_p,
double const epsilon_p)
118 4 * boost::math::pow<2>(delta_p) * boost::math::pow<3>(s.
I_1)) /
120 3 * beta_p + 6 * epsilon_p * s.
I_1;
123template <
int DisplacementDim>
128 double const gamma_p,
double const m_p)
131 (s.
D + s.
J_2 * m_p * gamma_p * dtheta_dsigma / one_gt.
value)) /
135template <
int DisplacementDim>
140 double const I_1_squared = boost::math::pow<2>(s.
I_1);
143 return std::sqrt(s.
J_2 * std::pow(1 + mp.
gamma * s.
J_3 /
144 (s.
J_2 * std::sqrt(s.
J_2)),
146 mp.
alpha / 2. * I_1_squared +
147 boost::math::pow<2>(mp.
delta) *
148 boost::math::pow<2>(I_1_squared)) +
152template <
int DisplacementDim>
153std::pair<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>,
158 constexpr int KelvinVectorSize =
168 return {KelvinVector::Zero(), KelvinMatrix::Zero()};
171 auto const& P_dev = Invariants::deviatoric_projection;
174 if (Invariants::determinant(s.
D) == 0)
176 OGS_FATAL(
"Determinant is zero. Matrix is non-invertable.");
180 KelvinVector
const sigma_D_inverse_D = P_dev * sigma_D_inverse;
182 KelvinVector
const dtheta_dsigma =
183 theta * sigma_D_inverse_D - 3. / 2. * theta / s.
J_2 * s.
D;
184 KelvinMatrix
const d2theta_dsigma2 =
186 sigma_D_inverse_D * dtheta_dsigma.transpose() -
187 3. / 2. * theta / s.
J_2 * P_dev -
188 3. / 2. * dtheta_dsigma / s.
J_2 * s.
D.transpose() +
189 3. / 2. * theta / boost::math::pow<2>(s.
J_2) * s.
D * s.
D.transpose();
191 return {dtheta_dsigma, d2theta_dsigma2};
194template <
int DisplacementDim>
202 double const eps_p_V,
203 double const eps_p_V_dot,
204 double const eps_p_eff_dot,
209 static int const KelvinVectorSize =
215 auto const& identity2 = Invariants::identity2;
217 double const theta = s.
J_3 / (s.
J_2 * std::sqrt(s.
J_2));
221 residual.template segment<KelvinVectorSize>(0).noalias() =
222 s.
value / mp.
G - 2 * (eps_D - eps_p_D) -
223 mp.
K / mp.
G * (eps_V - eps_p_V) * identity2;
225 auto const [dtheta_dsigma, d2theta_dsigma2] =
229 double const sqrtPhi = std::sqrt(
230 s.
J_2 * one_gt.pow_m_p + mp.
alpha_p / 2. * boost::math::pow<2>(s.
I_1) +
231 boost::math::pow<2>(mp.
delta_p) * boost::math::pow<4>(s.
I_1));
233 s, one_gt, sqrtPhi, dtheta_dsigma, mp.
gamma_p, mp.
m_p);
234 KelvinVector
const lambda_flow_D = lambda * flow_D;
236 residual.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
237 eps_p_D_dot - lambda_flow_D;
243 residual(2 * KelvinVectorSize, 0) = eps_p_V_dot - lambda * flow_V;
247 residual(2 * KelvinVectorSize + 1) =
249 std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
252 residual(2 * KelvinVectorSize + 2) =
yieldFunction(mp, s, k) / mp.
G;
256template <
int DisplacementDim>
263 static int const KelvinVectorSize =
271 auto const& P_dev = Invariants::deviatoric_projection;
272 auto const& identity2 = Invariants::identity2;
274 double const theta = s.
J_3 / (s.
J_2 * std::sqrt(s.
J_2));
277 auto const [dtheta_dsigma, d2theta_dsigma2] =
281 double const sqrtPhi = std::sqrt(
282 s.
J_2 * one_gt.pow_m_p + mp.
alpha_p / 2. * boost::math::pow<2>(s.
I_1) +
283 boost::math::pow<2>(mp.
delta_p) * boost::math::pow<4>(s.
I_1));
285 s, one_gt, sqrtPhi, dtheta_dsigma, mp.
gamma_p, mp.
m_p);
286 KelvinVector
const lambda_flow_D = lambda * flow_D;
292 jacobian.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
293 .noalias() = KelvinMatrix::Identity();
297 .template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
298 .noalias() = 2 * KelvinMatrix::Identity();
301 jacobian.template block<KelvinVectorSize, 1>(0, 2 * KelvinVectorSize)
302 .noalias() = mp.
K / mp.
G * identity2;
310 KelvinVector
const M0 = s.
J_2 / one_gt.value * dtheta_dsigma;
312 KelvinVector
const dPhi_dsigma =
313 one_gt.pow_m_p * (s.
D + gm_p * M0) +
315 4 * boost::math::pow<2>(mp.
delta_p) * boost::math::pow<3>(s.
I_1)) *
319 KelvinMatrix
const M1 =
321 (s.
D * dPhi_dsigma.transpose() + gm_p * M0 * dPhi_dsigma.transpose());
323 KelvinMatrix
const M2 =
324 one_gt.pow_m_p * (P_dev + s.
D * gm_p * M0.transpose());
327 KelvinMatrix
const M3 =
328 gm_p * one_gt.pow_m_p1 *
329 ((s.
D + (gm_p - mp.
gamma_p) * M0) * dtheta_dsigma.transpose() +
330 s.
J_2 * d2theta_dsigma2);
333 KelvinMatrix
const dflow_D_dsigma =
334 (-M1 / (4 * boost::math::pow<3>(sqrtPhi)) + (M2 + M3) / (2 * sqrtPhi)) *
337 .template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
338 .noalias() = -lambda * dflow_D_dsigma;
342 .template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
344 .noalias() = KelvinMatrix::Identity() / dt;
350 .template block<KelvinVectorSize, 1>(KelvinVectorSize,
351 2 * KelvinVectorSize + 2)
352 .noalias() = -flow_D;
357 KelvinVector
const dflow_V_dsigma =
360 boost::math::pow<3>(s.
I_1)) /
361 (4 * boost::math::pow<3>(sqrtPhi)) * dPhi_dsigma +
363 12 * boost::math::pow<2>(mp.
delta_p * s.
I_1) * identity2) /
367 jacobian.template block<1, KelvinVectorSize>(2 * KelvinVectorSize, 0)
368 .noalias() = -lambda * dflow_V_dsigma.transpose();
374 jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize) = 1. / dt;
382 jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize + 2) = -flow_V;
386 double const eff_flow =
387 std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
392 KelvinVector
const eff_flow23_lambda_flow_D =
393 -2 / 3. / eff_flow * lambda_flow_D;
396 .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 1, 0)
397 .noalias() = lambda * dflow_D_dsigma * eff_flow23_lambda_flow_D;
399 jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 2) =
400 eff_flow23_lambda_flow_D.transpose() * flow_D;
406 jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 1) = 1. / dt;
410 double const one_gt_pow_m = std::pow(one_gt.value, mp.
m);
411 double const gm = mp.
gamma * mp.
m;
413 KelvinVector
const dF_dsigma =
415 (one_gt_pow_m * (s.
D + gm * M0) +
417 boost::math::pow<3>(s.
I_1)) *
423 .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 2, 0)
424 .noalias() = dF_dsigma.transpose() / mp.
G;
428 jacobian(2 * KelvinVectorSize + 2, 2 * KelvinVectorSize + 1) =
437template <
int DisplacementDim>
439 double const K,
double const G)
441 static int const KelvinVectorSize =
445 auto const& P_dev = Invariants::deviatoric_projection;
446 auto const& P_sph = Invariants::spherical_projection;
450 return -2. * I * P_dev - 3. * K / G * I * P_sph;
454 double const hardening_coefficient,
455 double const eps_p_eff)
457 return kappa * (1. + eps_p_eff * hardening_coefficient);
460template <
int DisplacementDim>
462 double const G,
double const K,
468 static int const KelvinVectorSize =
471 auto const& P_dev = Invariants::deviatoric_projection;
474 double const pressure_prev = Invariants::trace(sigma_prev) / (-3. * G);
476 double const e_prev = Invariants::trace(eps_prev);
478 double const pressure = pressure_prev - K / G * (eps_V - e_prev);
481 P_dev * sigma_prev / G;
484 sigma_D_prev + 2 * P_dev * (eps - eps_prev);
485 return sigma_D - pressure * Invariants::identity2;
490template <
typename Res
idualVector,
typename KelvinVector>
491std::tuple<KelvinVector, PlasticStrain<KelvinVector>,
double>
494 static auto const size = KelvinVector::SizeAtCompileTime;
495 return std::forward_as_tuple(
496 solution.template segment<size>(size * 0),
498 solution[size * 2], solution[size * 2 + 1]},
499 solution[size * 2 + 2]);
502template <
int DisplacementDim>
506 std::unique_ptr<DamagePropertiesParameters>&& damage_properties,
508 : _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
509 _mp(std::move(material_properties)),
510 _damage_properties(std::move(damage_properties)),
511 _tangent_type(tangent_type)
515template <
int DisplacementDim>
523 material_state_variables)
const
526 &material_state_variables) !=
nullptr);
529 material_state_variables)
532 auto const& identity2 = Invariants::identity2;
533 return (eps - eps_p.D - eps_p.V / 3 * identity2).dot(sigma) / 2;
536template <
int DisplacementDim>
537std::optional<std::tuple<typename SolidEhlers<DisplacementDim>::KelvinVector,
539 DisplacementDim>::MaterialStateVariables>,
546 material_state_variables)
const
548 auto const& eps_m = std::get<MPL::SymmetricTensor<DisplacementDim>>(
550 auto const& eps_m_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
552 auto const& sigma_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
553 variable_array_prev.
stress);
556 &material_state_variables) !=
nullptr);
560 material_state_variables);
566 double const eps_V = Invariants::trace(eps_m);
568 auto const& P_dev = Invariants::deviatoric_projection;
576 mp.
G, mp.
K, sigma_prev, eps_m, eps_m_prev, eps_V);
586 if ((sigma.squaredNorm() == 0. || dt == 0. ||
590 state.
eps_p.eff)) < 0.))
593 mp.
K - 2. / 3 * mp.
G, mp.
G);
599 Eigen::FullPivLU<Eigen::Matrix<double, JacobianResidualSize,
600 JacobianResidualSize, Eigen::RowMajor>>
607 Eigen::Matrix<double, JacobianResidualSize, 1>;
609 Eigen::Matrix<double, JacobianResidualSize,
610 JacobianResidualSize, Eigen::RowMajor>;
621 auto const& eps_p_D =
622 solution.template segment<KelvinVectorSize>(
627 double const& eps_p_V = solution[KelvinVectorSize * 2];
628 double const eps_p_V_dot = (eps_p_V - state.
eps_p_prev.V) / dt;
630 double const& eps_p_eff = solution[KelvinVectorSize * 2 + 1];
631 double const eps_p_eff_dot =
636 solution[KelvinVectorSize * 2 + 1]);
639 solution.template segment<KelvinVectorSize>(
641 eps_p_D_dot, solution[KelvinVectorSize * 2], eps_p_V_dot,
642 eps_p_eff_dot, solution[KelvinVectorSize * 2 + 2],
649 dt, s, solution[KelvinVectorSize * 2 + 2], mp);
652 auto const update_solution =
655 solution += increment;
657 mp.
G * solution.template segment<KelvinVectorSize>(0)};
661 linear_solver, update_jacobian, update_residual,
662 update_solution, _nonlinear_solver_parameters);
664 auto const success_iterations = newton_solver.solve(jacobian);
666 if (!success_iterations)
675 if (*success_iterations == 0)
677 linear_solver.compute(jacobian);
680 std::tie(sigma, state.
eps_p, std::ignore) =
685 Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
688 Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
689 Eigen::RowMajor>::Zero();
690 dresidual_deps.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
703 linear_solver.solve(-dresidual_deps)
704 .template block<KelvinVectorSize, KelvinVectorSize>(0, 0);
713 "Unimplemented tangent type behaviour for the tangent type "
721 return {std::make_tuple(
730template <
int DisplacementDim>
731std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
734 return {{
"damage.kappa_d", 1,
736 DisplacementDim>::MaterialStateVariables
const& state,
737 std::vector<double>& cache) -> std::vector<double>
const&
741 auto const& ehlers_state =
749 state) -> std::span<double>
760 DisplacementDim>::MaterialStateVariables
const& state,
761 std::vector<double>& cache) -> std::vector<double>
const&
765 auto const& ehlers_state =
773 state) -> std::span<double>
782 {
"eps_p.D", KelvinVector::RowsAtCompileTime,
784 DisplacementDim>::MaterialStateVariables
const& state,
785 std::vector<double>& cache) -> std::vector<double>
const&
789 auto const& ehlers_state =
792 cache.resize(KelvinVector::RowsAtCompileTime);
794 cache, KelvinVector::RowsAtCompileTime) =
796 ehlers_state.eps_p.D);
801 state) -> std::span<double>
809 ehlers_state.
eps_p.D.data(),
810 static_cast<std::size_t
>(KelvinVector::RowsAtCompileTime)};
814 DisplacementDim>::MaterialStateVariables
const& state,
815 std::vector<double>& cache) -> std::vector<double>
const&
819 auto const& ehlers_state =
823 cache.front() = ehlers_state.
eps_p.V;
827 state) -> std::span<double>
834 return {&ehlers_state.
eps_p.V, 1};
838 DisplacementDim>::MaterialStateVariables
const& state,
839 std::vector<double>& cache) -> std::vector<double>
const&
843 auto const& ehlers_state =
847 cache.front() = ehlers_state.
eps_p.eff;
851 state) -> std::span<double>
858 return {&ehlers_state.
eps_p.eff, 1};
874 result(0, 0) = v(0) * v(0);
875 result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
876 result(0, 2) = result(2, 0) = v(5) * v(5) / 2.;
877 result(0, 3) = result(3, 0) = v(0) * v(3);
878 result(0, 4) = result(4, 0) = v(3) * v(5) / std::sqrt(2.);
879 result(0, 5) = result(5, 0) = v(0) * v(5);
881 result(1, 1) = v(1) * v(1);
882 result(1, 2) = result(2, 1) = v(4) * v(4) / 2.;
883 result(1, 3) = result(3, 1) = v(3) * v(1);
884 result(1, 4) = result(4, 1) = v(1) * v(4);
885 result(1, 5) = result(5, 1) = v(3) * v(4) / std::sqrt(2.);
887 result(2, 2) = v(2) * v(2);
888 result(2, 3) = result(3, 2) = v(5) * v(4) / std::sqrt(2.);
889 result(2, 4) = result(4, 2) = v(4) * v(2);
890 result(2, 5) = result(5, 2) = v(5) * v(2);
892 result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
893 result(3, 4) = result(4, 3) =
894 v(3) * v(4) / 2. + v(5) * v(1) / std::sqrt(2.);
895 result(3, 5) = result(5, 3) =
896 v(0) * v(4) / std::sqrt(2.) + v(3) * v(5) / 2.;
898 result(4, 4) = v(1) * v(2) + v(4) * v(4) / 2.;
899 result(4, 5) = result(5, 4) =
900 v(3) * v(2) / std::sqrt(2.) + v(5) * v(4) / 2.;
902 result(5, 5) = v(0) * v(2) + v(5) * v(5) / 2.;
912 result(0, 0) = v(0) * v(0);
913 result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
914 result(0, 2) = result(2, 0) = 0;
915 result(0, 3) = result(3, 0) = v(0) * v(3);
917 result(1, 1) = v(1) * v(1);
918 result(1, 2) = result(2, 1) = 0;
919 result(1, 3) = result(3, 1) = v(3) * v(1);
921 result(2, 2) = v(2) * v(2);
922 result(2, 3) = result(3, 2) = 0;
924 result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
std::vector< typename MechanicsBase< DisplacementDim >::InternalVariable > getInternalVariables() const override
Eigen::Matrix< double, JacobianResidualSize, JacobianResidualSize, Eigen::RowMajor > JacobianMatrix
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
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
Eigen::Matrix< double, JacobianResidualSize, 1 > ResidualVectorType
SolidEhlers(NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters, MaterialPropertiesParameters material_properties, std::unique_ptr< DamagePropertiesParameters > &&damage_properties, TangentType tangent_type)
double computeFreeEnergyDensity(double const, ParameterLib::SpatialPosition const &, double const, KelvinVector const &eps, KelvinVector const &sigma, typename MechanicsBase< DisplacementDim >::MaterialStateVariables const &material_state_variables) const override
KelvinVector mechanical_strain
constexpr auto to_underlying(E e) noexcept
Converts an enumeration to its underlying type.
MathLib::KelvinVector::KelvinMatrixType< 3 > sOdotS< 3 >(MathLib::KelvinVector::KelvinVectorType< 3 > const &v)
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > calculateDResidualDEps(double const K, double const G)
SolidEhlers< DisplacementDim >::JacobianMatrix calculatePlasticJacobian(double const dt, PhysicalStressWithInvariants< DisplacementDim > const &s, double const lambda, MaterialProperties const &mp)
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)
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > sOdotS(MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &v)
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)
double calculateIsotropicHardening(double const kappa, double const hardening_coefficient, double const eps_p_eff)
MathLib::KelvinVector::KelvinMatrixType< 2 > sOdotS< 2 >(MathLib::KelvinVector::KelvinVectorType< 2 > const &v)
std::tuple< KelvinVector, PlasticStrain< KelvinVector >, double > splitSolutionVector(ResidualVector const &solution)
std::pair< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > thetaSigmaDerivatives(double theta, PhysicalStressWithInvariants< DisplacementDim > const &s)
double yieldFunction(MaterialProperties const &mp, PhysicalStressWithInvariants< DisplacementDim > const &s, double const k)
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)
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)
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > elasticTangentStiffness(double const first_lame_parameter, double const shear_modulus)
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > inverse(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
double const hardening_coefficient
Holds powers of 1 + gamma_p*theta to base 0, m_p, and m_p-1.
OnePlusGamma_pTheta(double const gamma_p, double const theta, double const m_p)
PhysicalStressWithInvariants(KelvinVector const &stress)
static int const KelvinVectorSize
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
void setInitialConditions()
Damage damage
damage part of the state.
PlasticStrain< KelvinVector > eps_p
plastic part of the state.
double getEquivalentPlasticStrain() const override
PlasticStrain< KelvinVector > eps_p_prev
plastic part of the state.