OGS
MaterialLib::Solids::Ehlers Namespace Reference

Classes

class  Damage
 
struct  DamageProperties
 
struct  DamagePropertiesParameters
 
struct  MaterialProperties
 
struct  MaterialPropertiesParameters
 
struct  OnePlusGamma_pTheta
 Holds powers of 1 + gamma_p*theta to base 0, m_p, and m_p-1. More...
 
struct  PhysicalStressWithInvariants
 
struct  PlasticStrain
 
class  SolidEhlers
 
struct  StateVariables
 

Enumerations

enum class  TangentType { Elastic , PlasticDamageSecant , Plastic }
 

Functions

std::unique_ptr< DamagePropertiesParameterscreateDamageProperties (std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, BaseLib::ConfigTree const &config)
 
template<int DisplacementDim>
std::unique_ptr< SolidEhlers< DisplacementDim > > createEhlers (std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, BaseLib::ConfigTree const &config)
 
template<int DisplacementDim>
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > sOdotS (MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &v)
 
template<int DisplacementDim>
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)
 
template<int DisplacementDim>
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)
 
template<int DisplacementDim>
double yieldFunction (MaterialProperties const &mp, PhysicalStressWithInvariants< DisplacementDim > const &s, double const k)
 
template<int DisplacementDim>
std::pair< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > thetaSigmaDerivatives (double theta, PhysicalStressWithInvariants< DisplacementDim > const &s)
 
template<int DisplacementDim>
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)
 
template<int DisplacementDim>
SolidEhlers< DisplacementDim >::JacobianMatrix calculatePlasticJacobian (double const dt, PhysicalStressWithInvariants< DisplacementDim > const &s, double const lambda, MaterialProperties const &mp)
 
template<int DisplacementDim>
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > calculateDResidualDEps (double const K, double const G)
 
double calculateIsotropicHardening (double const kappa, double const hardening_coefficient, double const eps_p_eff)
 
template<int DisplacementDim>
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)
 
template<typename ResidualVector , typename KelvinVector >
std::tuple< KelvinVector, PlasticStrain< KelvinVector >, double > splitSolutionVector (ResidualVector const &solution)
 
template<>
MathLib::KelvinVector::KelvinMatrixType< 3 > sOdotS< 3 > (MathLib::KelvinVector::KelvinVectorType< 3 > const &v)
 
template<>
MathLib::KelvinVector::KelvinMatrixType< 2 > sOdotS< 2 > (MathLib::KelvinVector::KelvinVectorType< 2 > const &v)
 
TangentType makeTangentType (std::string const &s)
 

Enumeration Type Documentation

◆ TangentType

Enumerator
Elastic 
PlasticDamageSecant 
Plastic 

Definition at line 37 of file Ehlers.h.

Function Documentation

◆ calculateDResidualDEps()

template<int DisplacementDim>
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > MaterialLib::Solids::Ehlers::calculateDResidualDEps ( double const K,
double const G )

Calculates the derivative of the residuals with respect to total strain. Implementation fully implicit only.

Definition at line 438 of file Ehlers.cpp.

440{
441 static int const KelvinVectorSize =
444
445 auto const& P_dev = Invariants::deviatoric_projection;
446 auto const& P_sph = Invariants::spherical_projection;
447 auto const& I =
449
450 return -2. * I * P_dev - 3. * K / G * I * P_sph;
451}
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

References MathLib::KelvinVector::kelvin_vector_dimensions().

◆ calculateIsotropicHardening()

double MaterialLib::Solids::Ehlers::calculateIsotropicHardening ( double const kappa,
double const hardening_coefficient,
double const eps_p_eff )
inline

Definition at line 453 of file Ehlers.cpp.

456{
457 return kappa * (1. + eps_p_eff * hardening_coefficient);
458}

Referenced by MaterialLib::Solids::Ehlers::SolidEhlers< DisplacementDim >::integrateStress().

◆ calculatePlasticJacobian()

template<int DisplacementDim>
SolidEhlers< DisplacementDim >::JacobianMatrix MaterialLib::Solids::Ehlers::calculatePlasticJacobian ( double const dt,
PhysicalStressWithInvariants< DisplacementDim > const & s,
double const lambda,
MaterialProperties const & mp )

Definition at line 257 of file Ehlers.cpp.

262{
263 static int const KelvinVectorSize =
266 using KelvinVector =
268 using KelvinMatrix =
270
271 auto const& P_dev = Invariants::deviatoric_projection;
272 auto const& identity2 = Invariants::identity2;
273
274 double const theta = s.J_3 / (s.J_2 * std::sqrt(s.J_2));
275 OnePlusGamma_pTheta const one_gt{mp.gamma_p, theta, mp.m_p};
276
277 auto const [dtheta_dsigma, d2theta_dsigma2] =
278 thetaSigmaDerivatives(theta, s);
279
280 // deviatoric flow
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));
284 KelvinVector const flow_D = plasticFlowDeviatoricPart(
285 s, one_gt, sqrtPhi, dtheta_dsigma, mp.gamma_p, mp.m_p);
286 KelvinVector const lambda_flow_D = lambda * flow_D;
287
288 typename SolidEhlers<DisplacementDim>::JacobianMatrix jacobian =
289 SolidEhlers<DisplacementDim>::JacobianMatrix::Zero();
290
291 // G_11
292 jacobian.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
293 .noalias() = KelvinMatrix::Identity();
294
295 // G_12
296 jacobian
297 .template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
298 .noalias() = 2 * KelvinMatrix::Identity();
299
300 // G_13
301 jacobian.template block<KelvinVectorSize, 1>(0, 2 * KelvinVectorSize)
302 .noalias() = mp.K / mp.G * identity2;
303
304 // G_14 and G_15 are zero
305
306 // G_21 -- derivative of deviatoric flow
307
308 double const gm_p = mp.gamma_p * mp.m_p;
309 // intermediate variable for derivative of deviatoric flow
310 KelvinVector const M0 = s.J_2 / one_gt.value * dtheta_dsigma;
311 // derivative of Phi w.r.t. sigma
312 KelvinVector const dPhi_dsigma =
313 one_gt.pow_m_p * (s.D + gm_p * M0) +
314 (mp.alpha_p * s.I_1 +
315 4 * boost::math::pow<2>(mp.delta_p) * boost::math::pow<3>(s.I_1)) *
316 identity2;
317
318 // intermediate variable for derivative of deviatoric flow
319 KelvinMatrix const M1 =
320 one_gt.pow_m_p *
321 (s.D * dPhi_dsigma.transpose() + gm_p * M0 * dPhi_dsigma.transpose());
322 // intermediate variable for derivative of deviatoric flow
323 KelvinMatrix const M2 =
324 one_gt.pow_m_p * (P_dev + s.D * gm_p * M0.transpose());
325
326 // intermediate variable for derivative of deviatoric flow
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);
331
332 // derivative of flow_D w.r.t. sigma
333 KelvinMatrix const dflow_D_dsigma =
334 (-M1 / (4 * boost::math::pow<3>(sqrtPhi)) + (M2 + M3) / (2 * sqrtPhi)) *
335 mp.G;
336 jacobian
337 .template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
338 .noalias() = -lambda * dflow_D_dsigma;
339
340 // G_22
341 jacobian
342 .template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
343 KelvinVectorSize)
344 .noalias() = KelvinMatrix::Identity() / dt;
345
346 // G_23 and G_24 are zero
347
348 // G_25
349 jacobian
350 .template block<KelvinVectorSize, 1>(KelvinVectorSize,
351 2 * KelvinVectorSize + 2)
352 .noalias() = -flow_D;
353
354 // G_31
355 {
356 // derivative of flow_V w.r.t. sigma
357 KelvinVector const dflow_V_dsigma =
358 3 * mp.G *
359 (-(mp.alpha_p * s.I_1 + 4 * boost::math::pow<2>(mp.delta_p) *
360 boost::math::pow<3>(s.I_1)) /
361 (4 * boost::math::pow<3>(sqrtPhi)) * dPhi_dsigma +
362 (mp.alpha_p * identity2 +
363 12 * boost::math::pow<2>(mp.delta_p * s.I_1) * identity2) /
364 (2 * sqrtPhi) +
365 2 * mp.epsilon_p * identity2);
366
367 jacobian.template block<1, KelvinVectorSize>(2 * KelvinVectorSize, 0)
368 .noalias() = -lambda * dflow_V_dsigma.transpose();
369 }
370
371 // G_32 is zero
372
373 // G_33
374 jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize) = 1. / dt;
375
376 // G_34 is zero
377
378 // G_35
379 {
380 double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
381 s, sqrtPhi, mp.alpha_p, mp.beta_p, mp.delta_p, mp.epsilon_p);
382 jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize + 2) = -flow_V;
383 }
384
385 // increment of effectiv plastic strain
386 double const eff_flow =
387 std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
388
389 if (eff_flow > 0)
390 {
391 // intermediate variable for derivative of plastic jacobian
392 KelvinVector const eff_flow23_lambda_flow_D =
393 -2 / 3. / eff_flow * lambda_flow_D;
394 // G_41
395 jacobian
396 .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 1, 0)
397 .noalias() = lambda * dflow_D_dsigma * eff_flow23_lambda_flow_D;
398 // G_45
399 jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 2) =
400 eff_flow23_lambda_flow_D.transpose() * flow_D;
401 }
402
403 // G_42 and G_43 are zero
404
405 // G_44
406 jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 1) = 1. / dt;
407
408 // G_51
409 {
410 double const one_gt_pow_m = std::pow(one_gt.value, mp.m);
411 double const gm = mp.gamma * mp.m;
412 // derivative of yield function w.r.t. sigma
413 KelvinVector const dF_dsigma =
414 mp.G *
415 (one_gt_pow_m * (s.D + gm * M0) +
416 (mp.alpha * s.I_1 + 4 * boost::math::pow<2>(mp.delta) *
417 boost::math::pow<3>(s.I_1)) *
418 identity2) /
419 (2. * sqrtPhi) +
420 mp.G * (mp.beta + 2 * mp.epsilon_p * s.I_1) * identity2;
421
422 jacobian
423 .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 2, 0)
424 .noalias() = dF_dsigma.transpose() / mp.G;
425 }
426
427 // G_54
428 jacobian(2 * KelvinVectorSize + 2, 2 * KelvinVectorSize + 1) =
429 -mp.kappa * mp.hardening_coefficient / mp.G;
430
431 // G_52, G_53, G_55 are zero
432 return jacobian;
433}
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
KV::KelvinVectorType< DisplacementDim > KelvinVector
Definition Base.h:27
Holds powers of 1 + gamma_p*theta to base 0, m_p, and m_p-1.
Definition Ehlers.cpp:96

References MaterialLib::Solids::Ehlers::MaterialProperties::alpha, MaterialLib::Solids::Ehlers::MaterialProperties::alpha_p, MaterialLib::Solids::Ehlers::MaterialProperties::beta, MaterialLib::Solids::Ehlers::MaterialProperties::beta_p, MaterialLib::Solids::Ehlers::PhysicalStressWithInvariants< DisplacementDim >::D, MaterialLib::Solids::Ehlers::MaterialProperties::delta, MaterialLib::Solids::Ehlers::MaterialProperties::delta_p, MaterialLib::Solids::Ehlers::MaterialProperties::epsilon_p, MaterialLib::Solids::Ehlers::MaterialProperties::G, MaterialLib::Solids::Ehlers::MaterialProperties::gamma, MaterialLib::Solids::Ehlers::MaterialProperties::gamma_p, MaterialLib::Solids::Ehlers::MaterialProperties::hardening_coefficient, MaterialLib::Solids::Ehlers::PhysicalStressWithInvariants< DisplacementDim >::I_1, MaterialLib::Solids::Ehlers::PhysicalStressWithInvariants< DisplacementDim >::J_2, MaterialLib::Solids::Ehlers::PhysicalStressWithInvariants< DisplacementDim >::J_3, MaterialLib::Solids::Ehlers::MaterialProperties::K, MaterialLib::Solids::Ehlers::MaterialProperties::kappa, MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialLib::Solids::Ehlers::MaterialProperties::m, MaterialLib::Solids::Ehlers::MaterialProperties::m_p, plasticFlowDeviatoricPart(), and thetaSigmaDerivatives().

◆ calculatePlasticResidual()

template<int DisplacementDim>
SolidEhlers< DisplacementDim >::ResidualVectorType MaterialLib::Solids::Ehlers::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 )

Definition at line 196 of file Ehlers.cpp.

208{
209 static int const KelvinVectorSize =
212 using KelvinVector =
214
215 auto const& identity2 = Invariants::identity2;
216
217 double const theta = s.J_3 / (s.J_2 * std::sqrt(s.J_2));
218
220 // calculate stress residual
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;
224
225 auto const [dtheta_dsigma, d2theta_dsigma2] =
226 thetaSigmaDerivatives(theta, s);
227
228 OnePlusGamma_pTheta const one_gt{mp.gamma_p, theta, mp.m_p};
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));
232 KelvinVector const flow_D = plasticFlowDeviatoricPart(
233 s, one_gt, sqrtPhi, dtheta_dsigma, mp.gamma_p, mp.m_p);
234 KelvinVector const lambda_flow_D = lambda * flow_D;
235
236 residual.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
237 eps_p_D_dot - lambda_flow_D;
238
239 // plastic volume strain
240 {
241 double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
242 s, sqrtPhi, mp.alpha_p, mp.beta_p, mp.delta_p, mp.epsilon_p);
243 residual(2 * KelvinVectorSize, 0) = eps_p_V_dot - lambda * flow_V;
244 }
245
246 // evolution of plastic equivalent strain
247 residual(2 * KelvinVectorSize + 1) =
248 eps_p_eff_dot -
249 std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
250
251 // yield function (for plastic multiplier)
252 residual(2 * KelvinVectorSize + 2) = yieldFunction(mp, s, k) / mp.G;
253 return residual;
254}
Eigen::Matrix< double, JacobianResidualSize, 1 > ResidualVectorType
Definition Ehlers.h:282
std::pair< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > thetaSigmaDerivatives(double theta, PhysicalStressWithInvariants< DisplacementDim > const &s)
Definition Ehlers.cpp:155

References MaterialLib::Solids::Ehlers::MaterialProperties::alpha_p, MaterialLib::Solids::Ehlers::MaterialProperties::beta_p, MaterialLib::Solids::Ehlers::MaterialProperties::delta_p, MaterialLib::Solids::Ehlers::MaterialProperties::epsilon_p, MaterialLib::Solids::Ehlers::MaterialProperties::G, MaterialLib::Solids::Ehlers::MaterialProperties::gamma_p, MaterialLib::Solids::Ehlers::PhysicalStressWithInvariants< DisplacementDim >::I_1, MaterialLib::Solids::Ehlers::PhysicalStressWithInvariants< DisplacementDim >::J_2, MaterialLib::Solids::Ehlers::PhysicalStressWithInvariants< DisplacementDim >::J_3, MaterialLib::Solids::Ehlers::MaterialProperties::K, MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialLib::Solids::Ehlers::MaterialProperties::m_p, plasticFlowDeviatoricPart(), thetaSigmaDerivatives(), MaterialLib::Solids::Ehlers::PhysicalStressWithInvariants< DisplacementDim >::value, and yieldFunction().

◆ createDamageProperties()

std::unique_ptr< DamagePropertiesParameters > MaterialLib::Solids::Ehlers::createDamageProperties ( std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const & parameters,
BaseLib::ConfigTree const & config )
inline
Input File Parameter
material__solid__constitutive_relation__Ehlers__damage_properties__alpha_d
Input File Parameter
material__solid__constitutive_relation__Ehlers__damage_properties__beta_d
Input File Parameter
material__solid__constitutive_relation__Ehlers__damage_properties__h_d

Definition at line 23 of file CreateEhlers.h.

26{
28 auto& alpha_d =
29 ParameterLib::findParameter<double>(config, "alpha_d", parameters, 1);
30
31 DBUG("Use '{:s}' as alpha_d.", alpha_d.name);
32
34 auto& beta_d =
35 ParameterLib::findParameter<double>(config, "beta_d", parameters, 1);
36
37 DBUG("Use '{:s}' as beta_d.", beta_d.name);
38
40 auto& h_d =
41 ParameterLib::findParameter<double>(config, "h_d", parameters, 1);
42
43 DBUG("Use '{:s}' as h_d.", h_d.name);
44
45 return std::make_unique<DamagePropertiesParameters>(
46 DamagePropertiesParameters{alpha_d, beta_d, h_d});
47}
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30

References DBUG().

Referenced by createEhlers().

◆ createEhlers()

template<int DisplacementDim>
std::unique_ptr< SolidEhlers< DisplacementDim > > MaterialLib::Solids::Ehlers::createEhlers ( std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const & parameters,
BaseLib::ConfigTree const & config )
Input File Parameter
material__solid__constitutive_relation__type
Input File Parameter
material__solid__constitutive_relation__Ehlers__shear_modulus
Input File Parameter
material__solid__constitutive_relation__Ehlers__bulk_modulus
Input File Parameter
material__solid__constitutive_relation__Ehlers__kappa
Input File Parameter
material__solid__constitutive_relation__Ehlers__beta
Input File Parameter
material__solid__constitutive_relation__Ehlers__gamma
Input File Parameter
material__solid__constitutive_relation__Ehlers__hardening_modulus
Input File Parameter
material__solid__constitutive_relation__Ehlers__alpha
Input File Parameter
material__solid__constitutive_relation__Ehlers__delta
Input File Parameter
material__solid__constitutive_relation__Ehlers__eps
Input File Parameter
material__solid__constitutive_relation__Ehlers__m
Input File Parameter
material__solid__constitutive_relation__Ehlers__alphap
Input File Parameter
material__solid__constitutive_relation__Ehlers__deltap
Input File Parameter
material__solid__constitutive_relation__Ehlers__epsp
Input File Parameter
material__solid__constitutive_relation__Ehlers__mp
Input File Parameter
material__solid__constitutive_relation__Ehlers__betap
Input File Parameter
material__solid__constitutive_relation__Ehlers__gammap
Input File Parameter
material__solid__constitutive_relation__Ehlers__tangent_type
Input File Parameter
material__solid__constitutive_relation__Ehlers__damage_properties
Input File Parameter
material__solid__constitutive_relation__Ehlers__nonlinear_solver

Definition at line 50 of file CreateEhlers.h.

53{
55 config.checkConfigParameter("type", "Ehlers");
56 DBUG("Create Ehlers material");
57
59 auto& shear_modulus = ParameterLib::findParameter<double>(
60 config, "shear_modulus", parameters, 1);
61
62 DBUG("Use '{:s}' as shear modulus parameter.", shear_modulus.name);
63
65 auto& bulk_modulus = ParameterLib::findParameter<double>(
66 config, "bulk_modulus", parameters, 1);
67
68 DBUG("Use '{:s}' as bulk modulus parameter.", bulk_modulus.name);
69
71 auto& kappa =
72 ParameterLib::findParameter<double>(config, "kappa", parameters, 1);
73
74 DBUG("Use '{:s}' as kappa.", kappa.name);
75
77 auto& beta =
78 ParameterLib::findParameter<double>(config, "beta", parameters, 1);
79
80 DBUG("Use '{:s}' as beta.", beta.name);
81
83 auto& gamma =
84 ParameterLib::findParameter<double>(config, "gamma", parameters, 1);
85
86 DBUG("Use '{:s}' as gamma.", gamma.name);
87
89 auto& hardening_modulus = ParameterLib::findParameter<double>(
90 config, "hardening_modulus", parameters, 1);
91
92 DBUG("Use '{:s}' as hardening modulus parameter.", hardening_modulus.name);
93
95 auto& alpha =
96 ParameterLib::findParameter<double>(config, "alpha", parameters, 1);
97
98 DBUG("Use '{:s}' as alpha.", alpha.name);
99
101 auto& delta =
102 ParameterLib::findParameter<double>(config, "delta", parameters, 1);
103
104 DBUG("Use '{:s}' as delta.", delta.name);
105
107 auto& eps =
108 ParameterLib::findParameter<double>(config, "eps", parameters, 1);
109
110 DBUG("Use '{:s}' as eps.", eps.name);
111
113 auto& m = ParameterLib::findParameter<double>(config, "m", parameters, 1);
114
115 DBUG("Use '{:s}' as m.", m.name);
116
118 auto& alphap =
119 ParameterLib::findParameter<double>(config, "alphap", parameters, 1);
120
121 DBUG("Use '{:s}' as alphap.", alphap.name);
122
124 auto& deltap =
125 ParameterLib::findParameter<double>(config, "deltap", parameters, 1);
126
127 DBUG("Use '{:s}' as deltap.", deltap.name);
128
130 auto& epsp =
131 ParameterLib::findParameter<double>(config, "epsp", parameters, 1);
132
133 DBUG("Use '{:s}' as epsp.", epsp.name);
134
136 auto& paremeter_mp =
137 ParameterLib::findParameter<double>(config, "mp", parameters, 1);
138
139 DBUG("Use '{:s}' as mp.", paremeter_mp.name);
140
142 auto& betap =
143 ParameterLib::findParameter<double>(config, "betap", parameters, 1);
144
145 DBUG("Use '{:s}' as betap.", betap.name);
146
148 auto& gammap =
149 ParameterLib::findParameter<double>(config, "gammap", parameters, 1);
150
151 DBUG("Use '{:s}' as gammap.", gammap.name);
152
153 auto tangent_type =
155 makeTangentType(config.getConfigParameter<std::string>("tangent_type"));
156
158 shear_modulus, bulk_modulus, alpha, beta,
159 gamma, delta, eps, m,
160 alphap, betap, gammap, deltap,
161 epsp, paremeter_mp, kappa, hardening_modulus};
162
163 // Damage properties.
164 std::unique_ptr<DamagePropertiesParameters> ehlers_damage_properties;
165
166 auto const& ehlers_damage_config =
168 config.getConfigSubtreeOptional("damage_properties");
169 if (ehlers_damage_config)
170 {
171 ehlers_damage_properties =
172 createDamageProperties(parameters, *ehlers_damage_config);
173 }
174
175 auto const& nonlinear_solver_config =
177 config.getConfigSubtree("nonlinear_solver");
178 auto const nonlinear_solver_parameters =
179 NumLib::createNewtonRaphsonSolverParameters(nonlinear_solver_config);
180
181 return std::make_unique<SolidEhlers<DisplacementDim>>(
182 nonlinear_solver_parameters,
183 mp,
184 std::move(ehlers_damage_properties),
185 tangent_type);
186}
TangentType makeTangentType(std::string const &s)
Definition Ehlers.h:43
NewtonRaphsonSolverParameters createNewtonRaphsonSolverParameters(BaseLib::ConfigTree const &config)

References BaseLib::ConfigTree::checkConfigParameter(), createDamageProperties(), NumLib::createNewtonRaphsonSolverParameters(), DBUG(), BaseLib::ConfigTree::getConfigParameter(), BaseLib::ConfigTree::getConfigSubtree(), BaseLib::ConfigTree::getConfigSubtreeOptional(), and makeTangentType().

◆ makeTangentType()

TangentType MaterialLib::Solids::Ehlers::makeTangentType ( std::string const & s)
inline

Definition at line 43 of file Ehlers.h.

44{
45 if (s == "Elastic")
46 {
47 return TangentType::Elastic;
48 }
49 if (s == "PlasticDamageSecant")
50 {
51 return TangentType::PlasticDamageSecant;
52 }
53 if (s == "Plastic")
54 {
55 return TangentType::Plastic;
56 }
57 OGS_FATAL("Not valid string '{:s}' to create a tangent type from.", s);
58}
#define OGS_FATAL(...)
Definition Error.h:26

References Elastic, OGS_FATAL, Plastic, and PlasticDamageSecant.

Referenced by createEhlers().

◆ plasticFlowDeviatoricPart()

template<int DisplacementDim>
SolidEhlers< DisplacementDim >::KelvinVector MaterialLib::Solids::Ehlers::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 )

◆ plasticFlowVolumetricPart()

template<int DisplacementDim>
double MaterialLib::Solids::Ehlers::plasticFlowVolumetricPart ( PhysicalStressWithInvariants< DisplacementDim > const & s,
double const sqrtPhi,
double const alpha_p,
double const beta_p,
double const delta_p,
double const epsilon_p )

Definition at line 111 of file Ehlers.cpp.

115{
116 return 3 *
117 (alpha_p * s.I_1 +
118 4 * boost::math::pow<2>(delta_p) * boost::math::pow<3>(s.I_1)) /
119 (2 * sqrtPhi) +
120 3 * beta_p + 6 * epsilon_p * s.I_1;
121}

References MaterialLib::Solids::Ehlers::PhysicalStressWithInvariants< DisplacementDim >::I_1.

◆ predict_sigma()

template<int DisplacementDim>
SolidEhlers< DisplacementDim >::KelvinVector MaterialLib::Solids::Ehlers::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 )

Definition at line 461 of file Ehlers.cpp.

467{
468 static int const KelvinVectorSize =
471 auto const& P_dev = Invariants::deviatoric_projection;
472
473 // dimensionless initial hydrostatic pressure
474 double const pressure_prev = Invariants::trace(sigma_prev) / (-3. * G);
475 // initial strain invariant
476 double const e_prev = Invariants::trace(eps_prev);
477 // dimensioness hydrostatic stress increment
478 double const pressure = pressure_prev - K / G * (eps_V - e_prev);
479 // dimensionless deviatoric initial stress
480 typename SolidEhlers<DisplacementDim>::KelvinVector const sigma_D_prev =
481 P_dev * sigma_prev / G;
482 // dimensionless deviatoric stress
483 typename SolidEhlers<DisplacementDim>::KelvinVector const sigma_D =
484 sigma_D_prev + 2 * P_dev * (eps - eps_prev);
485 return sigma_D - pressure * Invariants::identity2;
486}
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
Definition Ehlers.h:295

References MathLib::KelvinVector::kelvin_vector_dimensions().

◆ sOdotS()

template<int DisplacementDim>
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > MaterialLib::Solids::Ehlers::sOdotS ( MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const & v)

Special product of v with itself: \(v \odot v\). The tensor v is given in Kelvin mapping.

Note
Implementation only for 2 and 3 dimensions.
Attention
Pay attention to the sign of the result, which normally would be negative, but the returned value is not negated. This has to do with \( d(A^{-1})/dA = -A^{-1} \odot A^{-1} \).

◆ sOdotS< 2 >()

Definition at line 869 of file Ehlers.cpp.

912{
914
915 result(0, 0) = v(0) * v(0);
916 result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
917 result(0, 2) = result(2, 0) = 0;
918 result(0, 3) = result(3, 0) = v(0) * v(3);
919
920 result(1, 1) = v(1) * v(1);
921 result(1, 2) = result(2, 1) = 0;
922 result(1, 3) = result(3, 1) = v(3) * v(1);
923
924 result(2, 2) = v(2) * v(2);
925 result(2, 3) = result(3, 2) = 0;
926
927 result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
928
929 return result;
930}

◆ sOdotS< 3 >()

Definition at line 869 of file Ehlers.cpp.

874{
876
877 result(0, 0) = v(0) * v(0);
878 result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
879 result(0, 2) = result(2, 0) = v(5) * v(5) / 2.;
880 result(0, 3) = result(3, 0) = v(0) * v(3);
881 result(0, 4) = result(4, 0) = v(3) * v(5) / std::sqrt(2.);
882 result(0, 5) = result(5, 0) = v(0) * v(5);
883
884 result(1, 1) = v(1) * v(1);
885 result(1, 2) = result(2, 1) = v(4) * v(4) / 2.;
886 result(1, 3) = result(3, 1) = v(3) * v(1);
887 result(1, 4) = result(4, 1) = v(1) * v(4);
888 result(1, 5) = result(5, 1) = v(3) * v(4) / std::sqrt(2.);
889
890 result(2, 2) = v(2) * v(2);
891 result(2, 3) = result(3, 2) = v(5) * v(4) / std::sqrt(2.);
892 result(2, 4) = result(4, 2) = v(4) * v(2);
893 result(2, 5) = result(5, 2) = v(5) * v(2);
894
895 result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
896 result(3, 4) = result(4, 3) =
897 v(3) * v(4) / 2. + v(5) * v(1) / std::sqrt(2.);
898 result(3, 5) = result(5, 3) =
899 v(0) * v(4) / std::sqrt(2.) + v(3) * v(5) / 2.;
900
901 result(4, 4) = v(1) * v(2) + v(4) * v(4) / 2.;
902 result(4, 5) = result(5, 4) =
903 v(3) * v(2) / std::sqrt(2.) + v(5) * v(4) / 2.;
904
905 result(5, 5) = v(0) * v(2) + v(5) * v(5) / 2.;
906 return result;
907}

◆ splitSolutionVector()

template<typename ResidualVector , typename KelvinVector >
std::tuple< KelvinVector, PlasticStrain< KelvinVector >, double > MaterialLib::Solids::Ehlers::splitSolutionVector ( ResidualVector const & solution)

Split the agglomerated solution vector in separate items. The arrangement must be the same as in the newton() function.

Definition at line 492 of file Ehlers.cpp.

493{
494 static auto const size = KelvinVector::SizeAtCompileTime;
495 return std::forward_as_tuple(
496 solution.template segment<size>(size * 0),
497 PlasticStrain<KelvinVector>{solution.template segment<size>(size * 1),
498 solution[size * 2], solution[size * 2 + 1]},
499 solution[size * 2 + 2]);
500}

◆ thetaSigmaDerivatives()

template<int DisplacementDim>
std::pair< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > MaterialLib::Solids::Ehlers::thetaSigmaDerivatives ( double theta,
PhysicalStressWithInvariants< DisplacementDim > const & s )

Definition at line 155 of file Ehlers.cpp.

157{
158 constexpr int KelvinVectorSize =
161 using KelvinVector =
163 using KelvinMatrix =
165
166 if (theta == 0)
167 {
168 return {KelvinVector::Zero(), KelvinMatrix::Zero()};
169 }
170
171 auto const& P_dev = Invariants::deviatoric_projection;
172
173 // inverse of deviatoric stress tensor
174 if (Invariants::determinant(s.D) == 0)
175 {
176 OGS_FATAL("Determinant is zero. Matrix is non-invertable.");
177 }
178 // inverse of sigma_D
179 KelvinVector const sigma_D_inverse = MathLib::KelvinVector::inverse(s.D);
180 KelvinVector const sigma_D_inverse_D = P_dev * sigma_D_inverse;
181
182 KelvinVector const dtheta_dsigma =
183 theta * sigma_D_inverse_D - 3. / 2. * theta / s.J_2 * s.D;
184 KelvinMatrix const d2theta_dsigma2 =
185 theta * P_dev * sOdotS<DisplacementDim>(sigma_D_inverse) * P_dev +
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();
190
191 return {dtheta_dsigma, d2theta_dsigma2};
192}
Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > inverse(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)

References MaterialLib::Solids::Ehlers::PhysicalStressWithInvariants< DisplacementDim >::D, MathLib::KelvinVector::inverse(), MaterialLib::Solids::Ehlers::PhysicalStressWithInvariants< DisplacementDim >::J_2, MathLib::KelvinVector::kelvin_vector_dimensions(), and OGS_FATAL.

Referenced by calculatePlasticJacobian(), and calculatePlasticResidual().

◆ yieldFunction()

template<int DisplacementDim>
double MaterialLib::Solids::Ehlers::yieldFunction ( MaterialProperties const & mp,
PhysicalStressWithInvariants< DisplacementDim > const & s,
double const k )