OGS
Lubby2.h
Go to the documentation of this file.
1
10#pragma once
11
15
16#include "MechanicsBase.h"
17
18namespace MaterialLib
19{
20namespace Solids
21{
22namespace Lubby2
23{
24//
25// Variables specific to the material model.
26//
28{
31 P const& GM0_,
32 P const& KM0_,
33 P const& etaK0_,
34 P const& etaM0_,
35 P const& mK_,
36 P const& mvK_,
37 P const& mvM_)
38 : GK0(GK0_),
39 GM0(GM0_),
40 KM0(KM0_),
41 etaK0(etaK0_),
42 etaM0(etaM0_),
43 mK(mK_),
44 mvK(mvK_),
45 mvM(mvM_)
46 {
47 }
48
49 // basic material parameters
50 P const& GK0;
51 P const& GM0;
52 P const& KM0;
53 P const& etaK0;
54 P const& etaM0;
55 P const& mK;
56 P const& mvK;
57 P const& mvM;
58};
59
60namespace detail
61{
62template <int DisplacementDim>
64{
65 LocalLubby2Properties(double const t,
68 : GM0(mp.GM0(t, x)[0]),
69 KM0(mp.KM0(t, x)[0]),
70 GK0(mp.GK0(t, x)[0]),
71 etaK0(mp.etaK0(t, x)[0]),
72 etaM0(mp.etaM0(t, x)[0]),
73 mK(mp.mK(t, x)[0]),
74 mvK(mp.mvK(t, x)[0]),
75 mvM(mp.mvM(t, x)[0])
76 {
77 }
78
79 void update(double const s_eff)
80 {
81 double const GM0_s_eff = GM0 * s_eff;
82 GK = GK0 * std::exp(mK * GM0_s_eff);
83 etaK = etaK0 * std::exp(mvK * GM0_s_eff);
84 etaM = etaM0 * std::exp(mvM * GM0_s_eff);
85 }
86
87 double const GM0;
88 double const KM0;
89 double const GK0;
90 double const etaK0;
91 double const etaM0;
92 double const mK;
93 double const mvK;
94 double const mvM;
95
96 // Solution dependent values.
97 double GK = std::numeric_limits<double>::quiet_NaN();
98 double etaK = std::numeric_limits<double>::quiet_NaN();
99 double etaM = std::numeric_limits<double>::quiet_NaN();
100};
101} // namespace detail
102
103template <int DisplacementDim>
104class Lubby2 final : public MechanicsBase<DisplacementDim>
105{
106public:
108 : public MechanicsBase<DisplacementDim>::MaterialStateVariables
109 {
111 {
112 // Previous time step values are not initialized and are set later.
115
116 // Initialize current time step values
117 eps_K_j.setZero(KelvinVectorSize);
118 eps_M_j.setZero(KelvinVectorSize);
119 }
120
122 {
125 }
126
127 void pushBackState() override
128 {
131 }
132
145
147 };
148
149 std::unique_ptr<
152 {
153 return std::unique_ptr<
156 }
157
158public:
159 static int const KelvinVectorSize =
165
166 static int const JacobianResidualSize =
167 3 * KelvinVectorSize; // Three is the number of components in the
168 // jacobian/residual, not the space dimension.
169 using ResidualVector = Eigen::Matrix<double, JacobianResidualSize, 1>;
170 using JacobianMatrix = Eigen::Matrix<double,
173 Eigen::RowMajor>;
174
175public:
176 explicit Lubby2(
177 NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters,
178 Lubby2MaterialProperties& material_properties)
179 : _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
180 _mp(material_properties)
181 {
182 }
183
185 double const t,
187 double const dt,
188 KelvinVector const& eps,
189 KelvinVector const& sigma,
191 material_state_variables) const override
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 =
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 }
225
226 double getBulkModulus(double const t,
228 KelvinMatrix const* const /*C*/) const override
229 {
230 return _mp.KM0(t, x)[0];
231 }
232
233 std::optional<std::tuple<KelvinVector,
234 std::unique_ptr<typename MechanicsBase<
235 DisplacementDim>::MaterialStateVariables>,
238 MaterialPropertyLib::VariableArray const& variable_array_prev,
239 MaterialPropertyLib::VariableArray const& variable_array,
240 double const t, ParameterLib::SpatialPosition const& x, double const dt,
242 material_state_variables) const override;
243
244private:
247 double const dt,
248 const KelvinVector& strain_curr,
249 const KelvinVector& strain_t,
250 const KelvinVector& stress_curr,
251 const KelvinVector& stress_t,
252 const KelvinVector& strain_Kel_curr,
253 const KelvinVector& strain_Kel_t,
254 const KelvinVector& strain_Max_curr,
255 const KelvinVector& strain_Max_t,
256 ResidualVector& res,
257 detail::LocalLubby2Properties<DisplacementDim> const& properties) const;
258
261 double const t,
263 double const dt,
264 JacobianMatrix& Jac,
265 double s_eff,
266 const KelvinVector& sig_i,
267 const KelvinVector& eps_K_i,
268 detail::LocalLubby2Properties<DisplacementDim> const& properties) const;
269
270private:
273};
274
275extern template class Lubby2<2>;
276extern template class Lubby2<3>;
277
278} // namespace Lubby2
279} // namespace Solids
280} // namespace MaterialLib
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
Definition Lubby2.cpp:75
NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters
Definition Lubby2.h:271
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
Definition Lubby2.h:184
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
Definition Lubby2.h:163
Lubby2(NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters, Lubby2MaterialProperties &material_properties)
Definition Lubby2.h:176
Eigen::Matrix< double, JacobianResidualSize, 1 > ResidualVector
Definition Lubby2.h:169
Lubby2MaterialProperties _mp
Definition Lubby2.h:272
Eigen::Matrix< double, JacobianResidualSize, JacobianResidualSize, Eigen::RowMajor > JacobianMatrix
Definition Lubby2.h:170
double getBulkModulus(double const t, ParameterLib::SpatialPosition const &x, KelvinMatrix const *const) const override
Definition Lubby2.h:226
static int const JacobianResidualSize
Definition Lubby2.h:166
static int const KelvinVectorSize
Definition Lubby2.h:159
std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables > createMaterialStateVariables() const override
Definition Lubby2.h:151
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
Definition Lubby2.h:161
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
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
Lubby2MaterialProperties(P const &GK0_, P const &GM0_, P const &KM0_, P const &etaK0_, P const &etaM0_, P const &mK_, P const &mvK_, P const &mvM_)
Definition Lubby2.h:30
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
Definition Lubby2.h:133
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
Definition Lubby2.h:135
LocalLubby2Properties(double const t, ParameterLib::SpatialPosition const &x, Lubby2MaterialProperties const &mp)
Definition Lubby2.h:65