12 #include <boost/math/special_functions/pow.hpp>
49 template <
int DisplacementDim>
52 return std::sqrt(2.0 / 3.0 * Invariants::FrobeniusNorm(eps_p.D.eval()));
61 template <
int DisplacementDim>
65 template <
int DisplacementDim>
97 :
value{1 + gamma_p * theta},
108 template <
int DisplacementDim>
111 double const sqrtPhi,
double const alpha_p,
double const beta_p,
112 double const delta_p,
double const epsilon_p)
116 4 * boost::math::pow<2>(delta_p) * boost::math::pow<3>(s.
I_1)) /
118 3 * beta_p + 6 * epsilon_p * s.
I_1;
121 template <
int DisplacementDim>
126 double const gamma_p,
double const m_p)
129 (s.
D + s.
J_2 * m_p * gamma_p * dtheta_dsigma / one_gt.
value)) /
133 template <
int DisplacementDim>
138 double const I_1_squared = boost::math::pow<2>(s.
I_1);
141 return std::sqrt(s.
J_2 * std::pow(1 + mp.
gamma * s.
J_3 /
142 (s.
J_2 * std::sqrt(s.
J_2)),
144 mp.
alpha / 2. * I_1_squared +
145 boost::math::pow<2>(mp.
delta) *
146 boost::math::pow<2>(I_1_squared)) +
150 template <
int DisplacementDim>
158 double const eps_p_V,
159 double const eps_p_V_dot,
160 double const eps_p_eff_dot,
165 static int const KelvinVectorSize =
171 auto const& P_dev = Invariants::deviatoric_projection;
172 auto const& identity2 = Invariants::identity2;
174 double const theta = s.
J_3 / (s.
J_2 * std::sqrt(s.
J_2));
178 residual.template segment<KelvinVectorSize>(0).noalias() =
179 s.
value / mp.
G - 2 * (eps_D - eps_p_D) -
180 mp.
K / mp.
G * (eps_V - eps_p_V) * identity2;
183 KelvinVector
const sigma_D_inverse_D =
185 KelvinVector
const dtheta_dsigma =
186 theta * sigma_D_inverse_D - 3. / 2. * theta / s.
J_2 * s.
D;
189 double const sqrtPhi = std::sqrt(
190 s.
J_2 * one_gt.pow_m_p + mp.
alpha_p / 2. * boost::math::pow<2>(s.
I_1) +
191 boost::math::pow<2>(mp.
delta_p) * boost::math::pow<4>(s.
I_1));
193 s, one_gt, sqrtPhi, dtheta_dsigma, mp.
gamma_p, mp.
m_p);
194 KelvinVector
const lambda_flow_D = lambda * flow_D;
196 residual.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
197 eps_p_D_dot - lambda_flow_D;
201 double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
203 residual(2 * KelvinVectorSize, 0) = eps_p_V_dot - lambda * flow_V;
207 residual(2 * KelvinVectorSize + 1) =
209 std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
212 residual(2 * KelvinVectorSize + 2) =
yieldFunction(mp, s, k) / mp.
G;
216 template <
int DisplacementDim>
223 static int const KelvinVectorSize =
231 auto const& P_dev = Invariants::deviatoric_projection;
232 auto const& identity2 = Invariants::identity2;
234 double const theta = s.
J_3 / (s.
J_2 * std::sqrt(s.
J_2));
238 if (Invariants::determinant(s.
D) == 0)
240 OGS_FATAL(
"Determinant is zero. Matrix is non-invertable.");
244 KelvinVector
const sigma_D_inverse_D = P_dev * sigma_D_inverse;
246 KelvinVector
const dtheta_dsigma =
247 theta * sigma_D_inverse_D - 3. / 2. * theta / s.
J_2 * s.
D;
250 double const sqrtPhi = std::sqrt(
251 s.
J_2 * one_gt.pow_m_p + mp.
alpha_p / 2. * boost::math::pow<2>(s.
I_1) +
252 boost::math::pow<2>(mp.
delta_p) * boost::math::pow<4>(s.
I_1));
254 s, one_gt, sqrtPhi, dtheta_dsigma, mp.
gamma_p, mp.
m_p);
255 KelvinVector
const lambda_flow_D = lambda * flow_D;
261 jacobian.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
262 .noalias() = KelvinMatrix::Identity();
266 .template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
267 .noalias() = 2 * KelvinMatrix::Identity();
270 jacobian.template block<KelvinVectorSize, 1>(0, 2 * KelvinVectorSize)
271 .noalias() = mp.
K / mp.
G * identity2;
279 KelvinVector
const M0 = s.
J_2 / one_gt.value * dtheta_dsigma;
281 KelvinVector
const dPhi_dsigma =
282 one_gt.pow_m_p * (s.
D + gm_p * M0) +
284 4 * boost::math::pow<2>(mp.
delta_p) * boost::math::pow<3>(s.
I_1)) *
288 KelvinMatrix
const M1 =
290 (s.
D * dPhi_dsigma.transpose() + gm_p * M0 * dPhi_dsigma.transpose());
292 KelvinMatrix
const M2 =
293 one_gt.pow_m_p * (P_dev + s.
D * gm_p * M0.transpose());
295 KelvinMatrix
const d2theta_dsigma2 =
296 theta * P_dev * sOdotS<DisplacementDim>(sigma_D_inverse) * P_dev +
297 sigma_D_inverse_D * dtheta_dsigma.transpose() -
298 3. / 2. * theta / s.
J_2 * P_dev -
299 3. / 2. * dtheta_dsigma / s.
J_2 * s.
D.transpose() +
300 3. / 2. * theta / boost::math::pow<2>(s.
J_2) * s.
D * s.
D.transpose();
303 KelvinMatrix
const M3 =
304 gm_p * one_gt.pow_m_p1 *
305 ((s.
D + (gm_p - mp.
gamma_p) * M0) * dtheta_dsigma.transpose() +
306 s.
J_2 * d2theta_dsigma2);
309 KelvinMatrix
const dflow_D_dsigma =
310 (-M1 / (4 * boost::math::pow<3>(sqrtPhi)) + (M2 + M3) / (2 * sqrtPhi)) *
313 .template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
314 .noalias() = -lambda * dflow_D_dsigma;
318 .template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
320 .noalias() = KelvinMatrix::Identity() / dt;
326 .template block<KelvinVectorSize, 1>(KelvinVectorSize,
327 2 * KelvinVectorSize + 2)
328 .noalias() = -flow_D;
333 KelvinVector
const dflow_V_dsigma =
336 boost::math::pow<3>(s.
I_1)) /
337 (4 * boost::math::pow<3>(sqrtPhi)) * dPhi_dsigma +
339 12 * boost::math::pow<2>(mp.
delta_p * s.
I_1) * identity2) /
343 jacobian.template block<1, KelvinVectorSize>(2 * KelvinVectorSize, 0)
344 .noalias() = -lambda * dflow_V_dsigma.transpose();
350 jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize) = 1. / dt;
356 double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
358 jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize + 2) = -flow_V;
362 double const eff_flow =
363 std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
368 KelvinVector
const eff_flow23_lambda_flow_D =
369 -2 / 3. / eff_flow * lambda_flow_D;
372 .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 1, 0)
373 .noalias() = lambda * dflow_D_dsigma * eff_flow23_lambda_flow_D;
375 jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 2) =
376 eff_flow23_lambda_flow_D.transpose() * flow_D;
382 jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 1) = 1. / dt;
386 double const one_gt_pow_m = std::pow(one_gt.value, mp.
m);
387 double const gm = mp.
gamma * mp.
m;
389 KelvinVector
const dF_dsigma =
391 (one_gt_pow_m * (s.
D + gm * M0) +
393 boost::math::pow<3>(s.
I_1)) *
399 .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 2, 0)
400 .noalias() = dF_dsigma.transpose() / mp.
G;
404 jacobian(2 * KelvinVectorSize + 2, 2 * KelvinVectorSize + 1) =
413 template <
int DisplacementDim>
415 double const K,
double const G)
417 static int const KelvinVectorSize =
421 auto const& P_dev = Invariants::deviatoric_projection;
422 auto const& P_sph = Invariants::spherical_projection;
426 return -2. * I * P_dev - 3. * K / G * I * P_sph;
430 double const hardening_coefficient,
431 double const eps_p_eff)
433 return kappa * (1. + eps_p_eff * hardening_coefficient);
436 template <
int DisplacementDim>
438 double const G,
double const K,
444 static int const KelvinVectorSize =
447 auto const& P_dev = Invariants::deviatoric_projection;
450 double const pressure_prev = Invariants::trace(sigma_prev) / (-3. * G);
452 double const e_prev = Invariants::trace(eps_prev);
454 double const pressure = pressure_prev - K / G * (eps_V - e_prev);
457 P_dev * sigma_prev / G;
460 sigma_D_prev + 2 * P_dev * (eps - eps_prev);
461 return sigma_D - pressure * Invariants::identity2;
466 template <
typename Res
idualVector,
typename KelvinVector>
467 std::tuple<KelvinVector, PlasticStrain<KelvinVector>,
double>
470 static auto const size = KelvinVector::SizeAtCompileTime;
471 return std::forward_as_tuple(
472 solution.template segment<size>(size * 0),
474 solution[size * 2], solution[size * 2 + 1]},
475 solution[size * 2 + 2]);
478 template <
int DisplacementDim>
482 std::unique_ptr<DamagePropertiesParameters>&& damage_properties,
484 : _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
485 _mp(std::move(material_properties)),
486 _damage_properties(std::move(damage_properties)),
487 _tangent_type(tangent_type)
491 template <
int DisplacementDim>
499 material_state_variables)
const
502 &material_state_variables) !=
nullptr);
505 material_state_variables)
508 auto const& identity2 = Invariants::identity2;
509 return (eps - eps_p.D - eps_p.V / 3 * identity2).dot(sigma) / 2;
512 template <
int DisplacementDim>
513 std::optional<std::tuple<typename SolidEhlers<DisplacementDim>::KelvinVector,
515 DisplacementDim>::MaterialStateVariables>,
522 material_state_variables)
const
524 auto const& eps_m = std::get<MPL::SymmetricTensor<DisplacementDim>>(
525 variable_array[
static_cast<int>(MPL::Variable::mechanical_strain)]);
526 auto const& eps_m_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
527 variable_array_prev[
static_cast<int>(
528 MPL::Variable::mechanical_strain)]);
529 auto const& sigma_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
530 variable_array_prev[
static_cast<int>(MPL::Variable::stress)]);
533 &material_state_variables) !=
nullptr);
537 material_state_variables);
543 double const eps_V = Invariants::trace(eps_m);
545 auto const& P_dev = Invariants::deviatoric_projection;
553 mp.
G, mp.
K, sigma_prev, eps_m, eps_m_prev, eps_V);
560 if ((sigma.squaredNorm() == 0 ||
564 state.
eps_p.eff)) < 0))
566 tangentStiffness = elasticTangentStiffness<DisplacementDim>(
567 mp.
K - 2. / 3 * mp.
G, mp.
G);
573 Eigen::FullPivLU<Eigen::Matrix<double, JacobianResidualSize,
574 JacobianResidualSize, Eigen::RowMajor>>
578 static int const KelvinVectorSize =
584 Eigen::Matrix<double, JacobianResidualSize, 1>;
586 Eigen::Matrix<double, JacobianResidualSize,
587 JacobianResidualSize, Eigen::RowMajor>;
598 auto const& eps_p_D =
599 solution.template segment<KelvinVectorSize>(
604 double const& eps_p_V = solution[KelvinVectorSize * 2];
605 double const eps_p_V_dot = (eps_p_V - state.
eps_p_prev.V) / dt;
607 double const& eps_p_eff = solution[KelvinVectorSize * 2 + 1];
608 double const eps_p_eff_dot =
613 solution[KelvinVectorSize * 2 + 1]);
614 residual = calculatePlasticResidual<DisplacementDim>(
616 solution.template segment<KelvinVectorSize>(
618 eps_p_D_dot, solution[KelvinVectorSize * 2], eps_p_V_dot,
619 eps_p_eff_dot, solution[KelvinVectorSize * 2 + 2],
625 jacobian = calculatePlasticJacobian<DisplacementDim>(
626 dt, s, solution[KelvinVectorSize * 2 + 2], mp);
629 auto const update_solution =
632 solution += increment;
634 mp.
G * solution.template segment<KelvinVectorSize>(0)};
640 decltype(update_residual), decltype(update_solution)>(
641 linear_solver, update_jacobian, update_residual,
642 update_solution, _nonlinear_solver_parameters);
644 auto const success_iterations = newton_solver.
solve(jacobian);
646 if (!success_iterations)
655 if (*success_iterations == 0)
657 linear_solver.compute(jacobian);
660 std::tie(sigma, state.
eps_p, std::ignore) =
661 splitSolutionVector<ResidualVectorType, KelvinVector>(solution);
665 Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
668 Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
669 Eigen::RowMajor>::Zero();
670 dresidual_deps.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
671 .noalias() = calculateDResidualDEps<DisplacementDim>(mp.
K, mp.
G);
676 elasticTangentStiffness<DisplacementDim>(mp.
K, mp.
G);
683 linear_solver.solve(-dresidual_deps)
684 .template block<KelvinVectorSize, KelvinVectorSize>(0, 0);
693 "Unimplemented tangent type behaviour for the tangent type "
701 return {std::make_tuple(
710 template <
int DisplacementDim>
711 std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
714 return {{
"damage.kappa_d", 1,
717 std::vector<double>& cache) -> std::vector<double>
const&
721 auto const& ehlers_state =
741 std::vector<double>& cache) -> std::vector<double>
const&
745 auto const& ehlers_state =
762 {
"eps_p.D", KelvinVector::RowsAtCompileTime,
765 std::vector<double>& cache) -> std::vector<double>
const&
769 auto const& ehlers_state =
772 cache.resize(KelvinVector::RowsAtCompileTime);
773 MathLib::toVector<KelvinVector>(
774 cache, KelvinVector::RowsAtCompileTime) =
776 ehlers_state.eps_p.D);
789 ehlers_state.
eps_p.D.data(),
790 static_cast<std::size_t
>(KelvinVector::RowsAtCompileTime)};
795 std::vector<double>& cache) -> std::vector<double>
const&
799 auto const& ehlers_state =
803 cache.front() = ehlers_state.
eps_p.V;
814 return {&ehlers_state.
eps_p.V, 1};
819 std::vector<double>& cache) -> std::vector<double>
const&
823 auto const& ehlers_state =
827 cache.front() = ehlers_state.
eps_p.eff;
838 return {&ehlers_state.
eps_p.eff, 1};
851 result(0, 0) = v(0) * v(0);
852 result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
853 result(0, 2) = result(2, 0) = v(5) * v(5) / 2.;
854 result(0, 3) = result(3, 0) = v(0) * v(3);
855 result(0, 4) = result(4, 0) = v(3) * v(5) / std::sqrt(2.);
856 result(0, 5) = result(5, 0) = v(0) * v(5);
858 result(1, 1) = v(1) * v(1);
859 result(1, 2) = result(2, 1) = v(4) * v(4) / 2.;
860 result(1, 3) = result(3, 1) = v(3) * v(1);
861 result(1, 4) = result(4, 1) = v(1) * v(4);
862 result(1, 5) = result(5, 1) = v(3) * v(4) / std::sqrt(2.);
864 result(2, 2) = v(2) * v(2);
865 result(2, 3) = result(3, 2) = v(5) * v(4) / std::sqrt(2.);
866 result(2, 4) = result(4, 2) = v(4) * v(2);
867 result(2, 5) = result(5, 2) = v(5) * v(2);
869 result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
870 result(3, 4) = result(4, 3) =
871 v(3) * v(4) / 2. + v(5) * v(1) / std::sqrt(2.);
872 result(3, 5) = result(5, 3) =
873 v(0) * v(4) / std::sqrt(2.) + v(3) * v(5) / 2.;
875 result(4, 4) = v(1) * v(2) + v(4) * v(4) / 2.;
876 result(4, 5) = result(5, 4) =
877 v(3) * v(2) / std::sqrt(2.) + v(5) * v(4) / 2.;
879 result(5, 5) = v(0) * v(2) + v(5) * v(5) / 2.;
889 result(0, 0) = v(0) * v(0);
890 result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
891 result(0, 2) = result(2, 0) = 0;
892 result(0, 3) = result(3, 0) = v(0) * v(3);
894 result(1, 1) = v(1) * v(1);
895 result(1, 2) = result(2, 1) = 0;
896 result(1, 3) = result(3, 1) = v(3) * v(1);
898 result(2, 2) = v(2) * v(2);
899 result(2, 3) = result(3, 2) = 0;
901 result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
std::vector< typename MechanicsBase< DisplacementDim >::InternalVariable > getInternalVariables() 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
Eigen::Matrix< double, JacobianResidualSize, JacobianResidualSize, Eigen::RowMajor > JacobianMatrix
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
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
std::optional< int > solve(JacobianMatrix &jacobian) const
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 > calculateDResidualDEps(double const K, double const G)
std::tuple< KelvinVector, PlasticStrain< KelvinVector >, double > splitSolutionVector(ResidualVector const &solution)
SolidEhlers< DisplacementDim >::JacobianMatrix calculatePlasticJacobian(double const dt, PhysicalStressWithInvariants< DisplacementDim > const &s, double const lambda, MaterialProperties const &mp)
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)
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< 3 > sOdotS< 3 >(MathLib::KelvinVector::KelvinVectorType< 3 > const &v)
double calculateIsotropicHardening(double const kappa, double const hardening_coefficient, double const eps_p_eff)
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > sOdotS(MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &v)
MathLib::KelvinVector::KelvinMatrixType< 2 > sOdotS< 2 >(MathLib::KelvinVector::KelvinVectorType< 2 > const &v)
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)
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
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)
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
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
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.