24 template <
int DisplacementDim>
25 Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
29 Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
32 Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
34 dGdE.template topLeftCorner<Lubby2<DisplacementDim>::KelvinVectorSize,
41 template <
int DisplacementDim,
typename LinearSolver>
43 double const GM0,
double const KM0, LinearSolver
const& linear_solver)
46 auto const dGdE = calculatedGdEBurgers<DisplacementDim>();
51 static int const KelvinVectorSize =
56 KelvinMatrix
const dzdE =
57 linear_solver.solve(-dGdE)
58 .template topLeftCorner<KelvinVectorSize, KelvinVectorSize>();
61 auto const& P_sph = Invariants::spherical_projection;
62 auto const& P_dev = Invariants::deviatoric_projection;
64 KelvinMatrix C = GM0 * dzdE * P_dev + 3. * KM0 * P_sph;
68 template <
int DisplacementDim>
69 std::optional<std::tuple<typename Lubby2<DisplacementDim>::KelvinVector,
71 DisplacementDim>::MaterialStateVariables>,
78 material_state_variables)
const
80 auto const& eps_m = std::get<MPL::SymmetricTensor<DisplacementDim>>(
81 variable_array[
static_cast<int>(MPL::Variable::mechanical_strain)]);
82 auto const& eps_m_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
83 variable_array_prev[
static_cast<int>(
84 MPL::Variable::mechanical_strain)]);
85 auto const& sigma_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
86 variable_array_prev[
static_cast<int>(MPL::Variable::stress)]);
91 &material_state_variables) !=
nullptr);
96 auto local_lubby2_properties =
100 auto const& P_dev = Invariants::deviatoric_projection;
107 KelvinVector sigd_t = P_dev * sigma_prev / local_lubby2_properties.GM0;
110 double sig_eff = Invariants::equivalentStress(sigd_j);
111 local_lubby2_properties.update(sig_eff);
113 using LocalJacobianMatrix =
114 Eigen::Matrix<double, KelvinVectorSize * 3, KelvinVectorSize * 3,
119 Eigen::FullPivLU<LocalJacobianMatrix> linear_solver;
132 LocalJacobianMatrix K_loc;
134 using LocalResidualVector =
135 Eigen::Matrix<double, KelvinVectorSize * 3, 1>;
137 auto const update_residual = [&](LocalResidualVector& residual)
139 calculateResidualBurgers(dt, eps_m_d_i, eps_m_d_t, sigd_j, sigd_t,
142 local_lubby2_properties);
145 auto const update_jacobian = [&](LocalJacobianMatrix& jacobian)
147 calculateJacobianBurgers(
148 t, x, dt, jacobian, sig_eff, sigd_j, state.
eps_K_j,
149 local_lubby2_properties);
152 auto const update_solution = [&](LocalResidualVector
const& increment)
155 sigd_j.noalias() += increment.template segment<KelvinVectorSize>(
156 KelvinVectorSize * 0);
158 increment.template segment<KelvinVectorSize>(KelvinVectorSize *
161 increment.template segment<KelvinVectorSize>(KelvinVectorSize *
166 KelvinVectorSize>::equivalentStress(sigd_j);
167 local_lubby2_properties.update(sig_eff);
171 decltype(linear_solver), LocalJacobianMatrix,
172 decltype(update_jacobian), LocalResidualVector,
173 decltype(update_residual), decltype(update_solution)>(
174 linear_solver, update_jacobian, update_residual, update_solution,
175 _nonlinear_solver_parameters);
177 auto const success_iterations = newton_solver.
solve(K_loc);
179 if (!success_iterations)
187 if (*success_iterations == 0)
189 linear_solver.compute(K_loc);
194 tangentStiffnessA<DisplacementDim>(local_lubby2_properties.GM0,
195 local_lubby2_properties.KM0,
199 double const delta_eps_m_trace = Invariants::trace(eps_m - eps_m_prev);
200 double const sigma_trace_prev = Invariants::trace(sigma_prev);
202 local_lubby2_properties.GM0 * sigd_j +
203 (local_lubby2_properties.KM0 * delta_eps_m_trace +
204 sigma_trace_prev / 3.) *
205 Invariants::identity2;
206 return {std::make_tuple(
214 template <
int DisplacementDim>
229 res.template segment<KelvinVectorSize>(0).noalias() =
230 (stress_curr - stress_t) -
231 2. * ((strain_curr - strain_t) - (strain_Kel_curr - strain_Kel_t) -
232 (strain_Max_curr - strain_Max_t));
235 res.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
236 (strain_Kel_curr - strain_Kel_t) -
237 dt / (2. * properties.
etaK) *
238 (properties.
GM0 * stress_curr -
239 2. * properties.
GK * strain_Kel_curr);
242 res.template segment<KelvinVectorSize>(2 * KelvinVectorSize).noalias() =
243 (strain_Max_curr - strain_Max_t) -
244 dt * 0.5 * properties.
GM0 / properties.
etaM * stress_curr;
247 template <
int DisplacementDim>
261 Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
266 Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
271 Jac.template block<KelvinVectorSize, KelvinVectorSize>(0,
272 2 * KelvinVectorSize)
277 Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
279 -0.5 * dt * properties.
GM0 / properties.
etaK * KelvinMatrix::Identity();
283 1. / (properties.
etaK * properties.
etaK) *
284 (properties.
GM0 * sig_i - 2. * properties.
GK * eps_K_i);
287 properties.
GM0 / s_eff * sig_i;
289 properties.
etaK / s_eff * sig_i;
290 Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
292 .noalias() += 0.5 * dt * eps_K_aid * dmu_vK.transpose() +
293 dt / properties.
etaK * eps_K_i * dG_K.transpose();
297 Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
300 .setConstant(1. + dt * properties.
GK / properties.
etaK);
305 Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
308 -0.5 * dt * properties.
GM0 / properties.
etaM * KelvinMatrix::Identity();
312 properties.
etaM / s_eff * sig_i;
313 Jac.template block<KelvinVectorSize, KelvinVectorSize>(
314 2 * KelvinVectorSize, 0)
315 .noalias() += 0.5 * dt * properties.
GM0 /
316 (properties.
etaM * properties.
etaM) * sig_i *
323 Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
324 2 * KelvinVectorSize)
Eigen::Matrix< double, JacobianResidualSize, 1 > ResidualVector
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
static int const KelvinVectorSize
Eigen::Matrix< double, JacobianResidualSize, JacobianResidualSize, Eigen::RowMajor > JacobianMatrix
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
std::optional< int > solve(JacobianMatrix &jacobian) const
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > tangentStiffnessA(double const GM0, double const KM0, LinearSolver const &linear_solver)
Eigen::Matrix< double, Lubby2< DisplacementDim >::JacobianResidualSize, Lubby2< DisplacementDim >::KelvinVectorSize > calculatedGdEBurgers()
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
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
void setInitialConditions()