OGS
MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim > Class Template Referencefinal

Detailed Description

template<int DisplacementDim>
class MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >

Definition at line 104 of file Lubby2.h.

#include <Lubby2.h>

Inheritance diagram for MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >:
[legend]
Collaboration diagram for MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >:
[legend]

Classes

struct  MaterialStateVariables
 

Public Types

using KelvinVector
 
using KelvinMatrix
 
using ResidualVector = Eigen::Matrix<double, JacobianResidualSize, 1>
 
using JacobianMatrix
 
- Public Types inherited from MaterialLib::Solids::MechanicsBase< DisplacementDim >
using KelvinVector
 
using KelvinMatrix
 

Public Member Functions

std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables > createMaterialStateVariables () const override
 
 Lubby2 (NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters, Lubby2MaterialProperties &material_properties)
 
double 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 override
 
double getBulkModulus (double const t, ParameterLib::SpatialPosition const &x, KelvinMatrix const *const) 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
 
- Public Member Functions inherited from MaterialLib::Solids::MechanicsBase< DisplacementDim >
virtual void initializeInternalStateVariables (double const, ParameterLib::SpatialPosition const &, typename MechanicsBase< DisplacementDim >::MaterialStateVariables &) const
 
virtual std::optional< std::tuple< KelvinVector, std::unique_ptr< MaterialStateVariables >, KelvinMatrix > > integrateStress (MaterialPropertyLib::VariableArray const &variable_array_prev, MaterialPropertyLib::VariableArray const &variable_array, double const t, ParameterLib::SpatialPosition const &x, double const dt, MaterialStateVariables const &material_state_variables) const =0
 
virtual std::vector< InternalVariablegetInternalVariables () const
 
virtual ConstitutiveModel getConstitutiveModel () const
 Gets the type of constitutive model.
 
virtual double getTemperatureRelatedCoefficient (double const, double const, ParameterLib::SpatialPosition const &, double const, double const) const
 
virtual double computeFreeEnergyDensity (double const t, ParameterLib::SpatialPosition const &x, double const dt, KelvinVector const &eps, KelvinVector const &sigma, MaterialStateVariables const &material_state_variables) const =0
 
virtual ~MechanicsBase ()=default
 

Static Public Attributes

static int const KelvinVectorSize
 
static int const JacobianResidualSize
 

Private Member Functions

void calculateResidualBurgers (double const dt, const KelvinVector &strain_curr, const KelvinVector &strain_t, const KelvinVector &stress_curr, const KelvinVector &stress_t, const KelvinVector &strain_Kel_curr, const KelvinVector &strain_Kel_t, const KelvinVector &strain_Max_curr, const KelvinVector &strain_Max_t, ResidualVector &res, detail::LocalLubby2Properties< DisplacementDim > const &properties) const
 Calculates the 18x1 residual vector.
 
void calculateJacobianBurgers (double const t, ParameterLib::SpatialPosition const &x, double const dt, JacobianMatrix &Jac, double s_eff, const KelvinVector &sig_i, const KelvinVector &eps_K_i, detail::LocalLubby2Properties< DisplacementDim > const &properties) const
 Calculates the 18x18 Jacobian.
 

Private Attributes

NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters
 
Lubby2MaterialProperties _mp
 

Member Typedef Documentation

◆ JacobianMatrix

template<int DisplacementDim>
using MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::JacobianMatrix
Initial value:
Eigen::Matrix<double,
Eigen::RowMajor>
static int const JacobianResidualSize
Definition Lubby2.h:166

Definition at line 170 of file Lubby2.h.

◆ KelvinMatrix

template<int DisplacementDim>
using MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::KelvinMatrix
Initial value:
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType

Definition at line 163 of file Lubby2.h.

◆ KelvinVector

template<int DisplacementDim>
using MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::KelvinVector
Initial value:
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType

Definition at line 161 of file Lubby2.h.

◆ ResidualVector

template<int DisplacementDim>
using MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::ResidualVector = Eigen::Matrix<double, JacobianResidualSize, 1>

Definition at line 169 of file Lubby2.h.

Constructor & Destructor Documentation

◆ Lubby2()

template<int DisplacementDim>
MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::Lubby2 ( NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters,
Lubby2MaterialProperties & material_properties )
inlineexplicit

Definition at line 176 of file Lubby2.h.

179 : _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
180 _mp(material_properties)
181 {
182 }
NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters
Definition Lubby2.h:271
Lubby2MaterialProperties _mp
Definition Lubby2.h:272

Member Function Documentation

◆ calculateJacobianBurgers()

template<int DisplacementDim>
void MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::calculateJacobianBurgers ( double const t,
ParameterLib::SpatialPosition const & x,
double const dt,
JacobianMatrix & Jac,
double s_eff,
const KelvinVector & sig_i,
const KelvinVector & eps_K_i,
detail::LocalLubby2Properties< DisplacementDim > const & properties ) const
private

Calculates the 18x18 Jacobian.

Definition at line 246 of file Lubby2.cpp.

255{
256 Jac.setZero();
257
258 // build G_11
259 Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
260 .diagonal()
261 .setConstant(1.);
262
263 // build G_12
264 Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
265 .diagonal()
266 .setConstant(2.);
267
268 // build G_13
269 Jac.template block<KelvinVectorSize, KelvinVectorSize>(0,
271 .diagonal()
272 .setConstant(2.);
273
274 // build G_21
275 Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
276 .noalias() =
277 -0.5 * dt * properties.GM0 / properties.etaK * KelvinMatrix::Identity();
278 if (s_eff > 0.)
279 {
280 KelvinVector const eps_K_aid =
281 1. / (properties.etaK * properties.etaK) *
282 (properties.GM0 * sig_i - 2. * properties.GK * eps_K_i);
283
284 KelvinVector const dG_K = 1.5 * _mp.mK(t, x)[0] * properties.GK *
285 properties.GM0 / s_eff * sig_i;
286 KelvinVector const dmu_vK = 1.5 * _mp.mvK(t, x)[0] * properties.GM0 *
287 properties.etaK / s_eff * sig_i;
288 Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
289 0)
290 .noalias() += 0.5 * dt * eps_K_aid * dmu_vK.transpose() +
291 dt / properties.etaK * eps_K_i * dG_K.transpose();
292 }
293
294 // build G_22
295 Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
297 .diagonal()
298 .setConstant(1. + dt * properties.GK / properties.etaK);
299
300 // nothing to do for G_23
301
302 // build G_31
303 Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
304 0)
305 .noalias() =
306 -0.5 * dt * properties.GM0 / properties.etaM * KelvinMatrix::Identity();
307 if (s_eff > 0.)
308 {
309 KelvinVector const dmu_vM = 1.5 * _mp.mvM(t, x)[0] * properties.GM0 *
310 properties.etaM / s_eff * sig_i;
311 Jac.template block<KelvinVectorSize, KelvinVectorSize>(
312 2 * KelvinVectorSize, 0)
313 .noalias() += 0.5 * dt * properties.GM0 /
314 (properties.etaM * properties.etaM) * sig_i *
315 dmu_vM.transpose();
316 }
317
318 // nothing to do for G_32
319
320 // build G_33
321 Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
323 .diagonal()
324 .setConstant(1.);
325}
static int const KelvinVectorSize
Definition Lubby2.h:159
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
Definition Lubby2.h:161

References MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::etaK, MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::etaM, MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::GK, and MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::GM0.

◆ calculateResidualBurgers()

template<int DisplacementDim>
void MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::calculateResidualBurgers ( double const dt,
const KelvinVector & strain_curr,
const KelvinVector & strain_t,
const KelvinVector & stress_curr,
const KelvinVector & stress_t,
const KelvinVector & strain_Kel_curr,
const KelvinVector & strain_Kel_t,
const KelvinVector & strain_Max_curr,
const KelvinVector & strain_Max_t,
ResidualVector & res,
detail::LocalLubby2Properties< DisplacementDim > const & properties ) const
private

Calculates the 18x1 residual vector.

Definition at line 213 of file Lubby2.cpp.

225{
226 // calculate stress residual
227 res.template segment<KelvinVectorSize>(0).noalias() =
228 (stress_curr - stress_t) -
229 2. * ((strain_curr - strain_t) - (strain_Kel_curr - strain_Kel_t) -
230 (strain_Max_curr - strain_Max_t));
231
232 // calculate Kelvin strain residual
233 res.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
234 (strain_Kel_curr - strain_Kel_t) -
235 dt / (2. * properties.etaK) *
236 (properties.GM0 * stress_curr -
237 2. * properties.GK * strain_Kel_curr);
238
239 // calculate Maxwell strain residual
240 res.template segment<KelvinVectorSize>(2 * KelvinVectorSize).noalias() =
241 (strain_Max_curr - strain_Max_t) -
242 dt * 0.5 * properties.GM0 / properties.etaM * stress_curr;
243}

References MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::etaK, MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::etaM, MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::GK, and MaterialLib::Solids::Lubby2::detail::LocalLubby2Properties< DisplacementDim >::GM0.

◆ computeFreeEnergyDensity()

template<int DisplacementDim>
double MaterialLib::Solids::Lubby2::Lubby2< 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
inlineoverride

Definition at line 184 of file Lubby2.h.

192 {
193 assert(dynamic_cast<MaterialStateVariables const*>(
194 &material_state_variables) != nullptr);
195 MaterialStateVariables state(static_cast<MaterialStateVariables const&>(
196 material_state_variables));
197
198 auto const& eps_K = state.eps_K_j;
199 auto const& eps_K_prev = state.eps_K_t;
200
201 auto const& eps_M = state.eps_M_j;
202
203 auto local_lubby2_properties =
204 detail::LocalLubby2Properties<DisplacementDim>{t, x, _mp};
205
206 // calculation of deviatoric parts
208 auto const& P_dev = Invariants::deviatoric_projection;
209 KelvinVector const epsd_i = P_dev * eps;
210
211 // initial guess as elastic predictor
212 KelvinVector sigd_j = 2.0 * (epsd_i - state.eps_M_t - state.eps_K_t);
213
214 // Calculate effective stress and update material properties
215 double sig_eff = Invariants::equivalentStress(sigd_j);
216 local_lubby2_properties.update(sig_eff);
217
218 auto const& eta_K = local_lubby2_properties.etaK;
219
220 // This is specific to the backward Euler time scheme and needs to be
221 // updated if other time schemes are used.
222 return (eps - eps_K - eps_M).dot(sigma) / 2 +
223 eps_K.dot(sigma - eta_K * (eps_K - eps_K_prev) / dt) / 2;
224 }

References MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::_mp, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_K_j, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_K_t, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_M_j, and MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_M_t.

◆ createMaterialStateVariables()

template<int DisplacementDim>
std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables > MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::createMaterialStateVariables ( ) const
inlineoverridevirtual

Polymorphic creator for MaterialStateVariables objects specific for a material model.

Reimplemented from MaterialLib::Solids::MechanicsBase< DisplacementDim >.

Definition at line 151 of file Lubby2.h.

152 {
153 return std::unique_ptr<
154 typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
155 new MaterialStateVariables};
156 }

◆ getBulkModulus()

template<int DisplacementDim>
double MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::getBulkModulus ( double const t,
ParameterLib::SpatialPosition const & x,
KelvinMatrix const * const  ) const
inlineoverridevirtual

◆ integrateStress()

template<int DisplacementDim>
std::optional< std::tuple< typename Lubby2< DisplacementDim >::KelvinVector, std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables >, typename Lubby2< DisplacementDim >::KelvinMatrix > > MaterialLib::Solids::Lubby2::Lubby2< 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
override

Definition at line 75 of file Lubby2.cpp.

81{
82 auto const& eps_m = std::get<MPL::SymmetricTensor<DisplacementDim>>(
83 variable_array.mechanical_strain);
84 auto const& eps_m_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
85 variable_array_prev.mechanical_strain);
86 auto const& sigma_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
87 variable_array_prev.stress);
88
90
91 assert(dynamic_cast<MaterialStateVariables const*>(
92 &material_state_variables) != nullptr);
93 MaterialStateVariables state(
94 static_cast<MaterialStateVariables const&>(material_state_variables));
95 state.setInitialConditions();
96
97 auto local_lubby2_properties =
98 detail::LocalLubby2Properties<DisplacementDim>{t, x, _mp};
99
100 // calculation of deviatoric parts
101 auto const& P_dev = Invariants::deviatoric_projection;
102 KelvinVector const eps_m_d_i = P_dev * eps_m;
103 KelvinVector const eps_m_d_t = P_dev * eps_m_prev;
104
105 // initial guess as elastic predictor.
106 KelvinVector sigd_j = 2.0 * (eps_m_d_i - state.eps_M_t - state.eps_K_t);
107 // Note: sigd_t contains dimensionless stresses!
108 KelvinVector sigd_t = P_dev * sigma_prev / local_lubby2_properties.GM0;
109
110 // Calculate effective stress and update material properties
111 double sig_eff = Invariants::equivalentStress(sigd_j);
112 local_lubby2_properties.update(sig_eff);
113
114 using LocalJacobianMatrix =
115 Eigen::Matrix<double, KelvinVectorSize * 3, KelvinVectorSize * 3,
116 Eigen::RowMajor>;
117
118 // Linear solver for the newton loop is required after the loop with the
119 // same matrix. This saves one decomposition.
120 Eigen::FullPivLU<LocalJacobianMatrix> linear_solver;
121
122 // Different solvers are available for the solution of the local system.
123 // TODO Make the following choice of linear solvers available from the
124 // input file configuration:
125 // K_loc.partialPivLu().solve(-res_loc);
126 // K_loc.fullPivLu().solve(-res_loc);
127 // K_loc.householderQr().solve(-res_loc);
128 // K_loc.colPivHouseholderQr().solve(res_loc);
129 // K_loc.fullPivHouseholderQr().solve(-res_loc);
130 // K_loc.llt().solve(-res_loc);
131 // K_loc.ldlt().solve(-res_loc);
132
133 LocalJacobianMatrix K_loc;
134 { // Local Newton solver
135 using LocalResidualVector =
136 Eigen::Matrix<double, KelvinVectorSize * 3, 1>;
137
138 auto const update_residual = [&](LocalResidualVector& residual)
139 {
140 calculateResidualBurgers(dt, eps_m_d_i, eps_m_d_t, sigd_j, sigd_t,
141 state.eps_K_j, state.eps_K_t,
142 state.eps_M_j, state.eps_M_t, residual,
143 local_lubby2_properties);
144 };
145
146 auto const update_jacobian = [&](LocalJacobianMatrix& jacobian)
147 {
149 t, x, dt, jacobian, sig_eff, sigd_j, state.eps_K_j,
150 local_lubby2_properties); // for solution dependent Jacobians
151 };
152
153 auto const update_solution = [&](LocalResidualVector const& increment)
154 {
155 // increment solution vectors
156 sigd_j.noalias() += increment.template segment<KelvinVectorSize>(
157 KelvinVectorSize * 0);
158 state.eps_K_j.noalias() +=
159 increment.template segment<KelvinVectorSize>(KelvinVectorSize *
160 1);
161 state.eps_M_j.noalias() +=
162 increment.template segment<KelvinVectorSize>(KelvinVectorSize *
163 2);
164
165 // Calculate effective stress and update material properties
167 KelvinVectorSize>::equivalentStress(sigd_j);
168 local_lubby2_properties.update(sig_eff);
169 };
170
171 auto newton_solver = NumLib::NewtonRaphson(
172 linear_solver, update_jacobian, update_residual, update_solution,
174
175 auto const success_iterations = newton_solver.solve(K_loc);
176
177 if (!success_iterations)
178 {
179 return {};
180 }
181
182 // If the Newton loop didn't run, the linear solver will not be
183 // initialized.
184 // This happens usually for the first iteration of the first timestep.
185 if (*success_iterations == 0)
186 {
187 linear_solver.compute(K_loc);
188 }
189 }
190
192 tangentStiffnessA<DisplacementDim>(local_lubby2_properties.GM0,
193 local_lubby2_properties.KM0,
194 linear_solver);
195
196 // Hydrostatic part for the stress and the tangent.
197 double const delta_eps_m_trace = Invariants::trace(eps_m - eps_m_prev);
198 double const sigma_trace_prev = Invariants::trace(sigma_prev);
199 KelvinVector const sigma =
200 local_lubby2_properties.GM0 * sigd_j +
201 (local_lubby2_properties.KM0 * delta_eps_m_trace +
202 sigma_trace_prev / 3.) *
203 Invariants::identity2;
204 return {std::make_tuple(
205 sigma,
206 std::unique_ptr<
207 typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
208 new MaterialStateVariables{state}},
209 C)};
210}
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
Definition Lubby2.h:163
void calculateJacobianBurgers(double const t, ParameterLib::SpatialPosition const &x, double const dt, JacobianMatrix &Jac, double s_eff, const KelvinVector &sig_i, const KelvinVector &eps_K_i, detail::LocalLubby2Properties< DisplacementDim > const &properties) const
Calculates the 18x18 Jacobian.
Definition Lubby2.cpp:246
void calculateResidualBurgers(double const dt, const KelvinVector &strain_curr, const KelvinVector &strain_t, const KelvinVector &stress_curr, const KelvinVector &stress_t, const KelvinVector &strain_Kel_curr, const KelvinVector &strain_Kel_t, const KelvinVector &strain_Max_curr, const KelvinVector &strain_Max_t, ResidualVector &res, detail::LocalLubby2Properties< DisplacementDim > const &properties) const
Calculates the 18x1 residual vector.
Definition Lubby2.cpp:213
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > tangentStiffnessA(double const GM0, double const KM0, LinearSolver const &linear_solver)
Definition Lubby2.cpp:44

References MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_K_j, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_K_t, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_M_j, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::eps_M_t, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::setInitialConditions(), MaterialPropertyLib::VariableArray::stress, and MaterialLib::Solids::Lubby2::tangentStiffnessA().

Member Data Documentation

◆ _mp

◆ _nonlinear_solver_parameters

template<int DisplacementDim>
NumLib::NewtonRaphsonSolverParameters const MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::_nonlinear_solver_parameters
private

Definition at line 271 of file Lubby2.h.

◆ JacobianResidualSize

template<int DisplacementDim>
int const MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::JacobianResidualSize
static
Initial value:

Definition at line 166 of file Lubby2.h.

◆ KelvinVectorSize

template<int DisplacementDim>
int const MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::KelvinVectorSize
static
Initial value:
=
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.

Definition at line 159 of file Lubby2.h.

Referenced by MaterialLib::Solids::Lubby2::Lubby2< DisplacementDim >::MaterialStateVariables::MaterialStateVariables().


The documentation for this class was generated from the following files: