Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "CappedMohrCoulombCosseratStressUpdate.h"
11 : #include "libmesh/utility.h"
12 :
13 : registerMooseObject("SolidMechanicsApp", CappedMohrCoulombCosseratStressUpdate);
14 :
15 : InputParameters
16 360 : CappedMohrCoulombCosseratStressUpdate::validParams()
17 : {
18 360 : InputParameters params = CappedMohrCoulombStressUpdate::validParams();
19 360 : params.addClassDescription("Capped Mohr-Coulomb plasticity stress calculator for the Cosserat "
20 : "situation where the host medium (ie, the limit where all Cosserat "
21 : "effects are zero) is isotropic. Note that the return-map flow rule "
22 : "uses an isotropic elasticity tensor built with the 'host' properties "
23 : "defined by the user.");
24 720 : params.addRequiredRangeCheckedParam<Real>("host_youngs_modulus",
25 : "host_youngs_modulus>0",
26 : "Young's modulus for the isotropic host medium");
27 720 : params.addRequiredRangeCheckedParam<Real>("host_poissons_ratio",
28 : "host_poissons_ratio>=0 & host_poissons_ratio<0.5",
29 : "Poisson's ratio for the isotropic host medium");
30 360 : return params;
31 0 : }
32 :
33 270 : CappedMohrCoulombCosseratStressUpdate::CappedMohrCoulombCosseratStressUpdate(
34 270 : const InputParameters & parameters)
35 : : CappedMohrCoulombStressUpdate(parameters),
36 270 : _host_young(getParam<Real>("host_youngs_modulus")),
37 540 : _host_poisson(getParam<Real>("host_poissons_ratio")),
38 270 : _host_E0011(_host_young * _host_poisson / (1.0 + _host_poisson) / (1.0 - 2.0 * _host_poisson)),
39 270 : _host_E0000(_host_E0011 + _host_young / (1.0 + _host_poisson))
40 : {
41 270 : }
42 :
43 : void
44 11616 : CappedMohrCoulombCosseratStressUpdate::preReturnMapV(
45 : const std::vector<Real> & /*trial_stress_params*/,
46 : const RankTwoTensor & stress_trial,
47 : const std::vector<Real> & /*intnl_old*/,
48 : const std::vector<Real> & /*yf*/,
49 : const RankFourTensor & /*Eijkl*/)
50 : {
51 : std::vector<Real> eigvals;
52 11616 : stress_trial.symmetricEigenvaluesEigenvectors(eigvals, _eigvecs);
53 11616 : _poissons_ratio = _host_poisson;
54 11616 : }
55 :
56 : void
57 11616 : CappedMohrCoulombCosseratStressUpdate::setEffectiveElasticity(const RankFourTensor & /*Eijkl*/)
58 : {
59 11616 : _Eij[0][0] = _Eij[1][1] = _Eij[2][2] = _host_E0000;
60 11616 : _Eij[0][1] = _Eij[1][0] = _Eij[0][2] = _Eij[2][0] = _Eij[1][2] = _Eij[2][1] = _host_E0011;
61 11616 : _En = _Eij[2][2];
62 11616 : const Real denom = _Eij[0][0] * (_Eij[0][0] + _Eij[0][1]) - 2 * Utility::pow<2>(_Eij[0][1]);
63 46464 : for (unsigned a = 0; a < _num_sp; ++a)
64 : {
65 34848 : _Cij[a][a] = (_Eij[0][0] + _Eij[0][1]) / denom;
66 69696 : for (unsigned b = 0; b < a; ++b)
67 34848 : _Cij[a][b] = _Cij[b][a] = -_Eij[0][1] / denom;
68 : }
69 11616 : }
70 :
71 : void
72 11616 : CappedMohrCoulombCosseratStressUpdate::setStressAfterReturnV(
73 : const RankTwoTensor & stress_trial,
74 : const std::vector<Real> & stress_params,
75 : Real /*gaE*/,
76 : const std::vector<Real> & /*intnl*/,
77 : const yieldAndFlow & /*smoothed_q*/,
78 : const RankFourTensor & /*Eijkl*/,
79 : RankTwoTensor & stress) const
80 : {
81 : // form the diagonal stress
82 11616 : stress = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0);
83 : // rotate to the original frame, to give the symmetric part of the stress
84 11616 : stress = _eigvecs * stress * (_eigvecs.transpose());
85 : // add the non-symmetric parts
86 11616 : stress += 0.5 * (stress_trial - stress_trial.transpose());
87 11616 : }
88 :
89 : void
90 256 : CappedMohrCoulombCosseratStressUpdate::consistentTangentOperatorV(
91 : const RankTwoTensor & stress_trial,
92 : const std::vector<Real> & trial_stress_params,
93 : const RankTwoTensor & stress,
94 : const std::vector<Real> & stress_params,
95 : Real gaE,
96 : const yieldAndFlow & smoothed_q,
97 : const RankFourTensor & elasticity_tensor,
98 : bool compute_full_tangent_operator,
99 : const std::vector<std::vector<Real>> & dvar_dtrial,
100 : RankFourTensor & cto)
101 : {
102 256 : CappedMohrCoulombStressUpdate::consistentTangentOperatorV(stress_trial,
103 : trial_stress_params,
104 : stress,
105 : stress_params,
106 : gaE,
107 : smoothed_q,
108 : elasticity_tensor,
109 : compute_full_tangent_operator,
110 : dvar_dtrial,
111 : cto);
112 :
113 256 : if (!compute_full_tangent_operator)
114 0 : return;
115 :
116 : /**
117 : * Add the correction for the antisymmetric part of the elasticity
118 : * tensor.
119 : * CappedMohrCoulombStressUpdate computes
120 : * cto(i, j, k, l) = dstress(i, j)/dstrain(k, l)
121 : * and during the computations it explicitly performs certain
122 : * contractions that result in symmetry between i and j, and k and l,
123 : * viz, cto(i, j, k, l) = cto(j, i, k, l) = cto(i, j, l, k)
124 : * That is correct because that plasticity model is only valid for
125 : * symmetric stresses and strains.
126 : * CappedMohrCoulombCosseratStressUpdate does not include contributions from the
127 : * antisymmetric parts of stress (or strain), so the antisymmetric
128 : * parts of cto are just the antisymmetric parts of the elasticity
129 : * tensor, which must now get added to the cto computed by
130 : * CappedMohrCoulombStressUpdate
131 : */
132 256 : RankFourTensor anti;
133 1024 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
134 3072 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
135 9216 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
136 27648 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
137 20736 : anti(i, j, k, l) = 0.5 * (elasticity_tensor(i, j, k, l) - elasticity_tensor(j, i, k, l));
138 :
139 256 : cto += anti;
140 : }
|
