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 "CappedMohrCoulombStressUpdate.h"
11 : #include "libmesh/utility.h"
12 : #include "ElasticityTensorTools.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", CappedMohrCoulombStressUpdate);
15 :
16 : InputParameters
17 1696 : CappedMohrCoulombStressUpdate::validParams()
18 : {
19 1696 : InputParameters params = MultiParameterPlasticityStressUpdate::validParams();
20 3392 : params.addRequiredParam<UserObjectName>(
21 : "tensile_strength",
22 : "A SolidMechanicsHardening UserObject that defines hardening of the "
23 : "tensile strength. In physical situations this is positive (and always "
24 : "must be greater than negative compressive-strength.");
25 3392 : params.addRequiredParam<UserObjectName>(
26 : "compressive_strength",
27 : "A SolidMechanicsHardening UserObject that defines hardening of the "
28 : "compressive strength. In physical situations this is positive.");
29 3392 : params.addRequiredParam<UserObjectName>(
30 : "cohesion", "A SolidMechanicsHardening UserObject that defines hardening of the cohesion");
31 3392 : params.addRequiredParam<UserObjectName>("friction_angle",
32 : "A SolidMechanicsHardening UserObject "
33 : "that defines hardening of the "
34 : "friction angle (in radians)");
35 3392 : params.addRequiredParam<UserObjectName>(
36 : "dilation_angle",
37 : "A SolidMechanicsHardening UserObject that defines hardening of the "
38 : "dilation angle (in radians). Unless you are quite confident, this should "
39 : "be set positive and not greater than the friction angle.");
40 3392 : params.addParam<bool>("perfect_guess",
41 3392 : true,
42 : "Provide a guess to the Newton-Raphson procedure "
43 : "that is the result from perfect plasticity. With "
44 : "severe hardening/softening this may be "
45 : "suboptimal.");
46 1696 : params.addClassDescription("Nonassociative, smoothed, Mohr-Coulomb plasticity capped with "
47 : "tensile (Rankine) and compressive caps, with hardening/softening");
48 1696 : return params;
49 0 : }
50 :
51 1272 : CappedMohrCoulombStressUpdate::CappedMohrCoulombStressUpdate(const InputParameters & parameters)
52 : : MultiParameterPlasticityStressUpdate(parameters, 3, 12, 2),
53 1272 : _tensile_strength(getUserObject<SolidMechanicsHardeningModel>("tensile_strength")),
54 1272 : _compressive_strength(getUserObject<SolidMechanicsHardeningModel>("compressive_strength")),
55 1272 : _cohesion(getUserObject<SolidMechanicsHardeningModel>("cohesion")),
56 1272 : _phi(getUserObject<SolidMechanicsHardeningModel>("friction_angle")),
57 1272 : _psi(getUserObject<SolidMechanicsHardeningModel>("dilation_angle")),
58 2544 : _perfect_guess(getParam<bool>("perfect_guess")),
59 1272 : _poissons_ratio(0.0),
60 1272 : _shifter(_f_tol),
61 2544 : _eigvecs(RankTwoTensor())
62 : {
63 1272 : if (_psi.value(0.0) <= 0.0 || _psi.value(0.0) > _phi.value(0.0))
64 0 : mooseWarning("Usually the Mohr-Coulomb dilation angle is positive and not greater than the "
65 : "friction angle");
66 1272 : }
67 :
68 : void
69 188032 : CappedMohrCoulombStressUpdate::computeStressParams(const RankTwoTensor & stress,
70 : std::vector<Real> & stress_params) const
71 : {
72 : // stress_params[0] = smallest eigenvalue, stress_params[2] = largest eigenvalue
73 188032 : stress.symmetricEigenvalues(stress_params);
74 188032 : }
75 :
76 : std::vector<RankTwoTensor>
77 141464 : CappedMohrCoulombStressUpdate::dstress_param_dstress(const RankTwoTensor & stress) const
78 : {
79 : std::vector<Real> sp;
80 : std::vector<RankTwoTensor> dsp;
81 141464 : stress.dsymmetricEigenvalues(sp, dsp);
82 141464 : return dsp;
83 141464 : }
84 :
85 : std::vector<RankFourTensor>
86 0 : CappedMohrCoulombStressUpdate::d2stress_param_dstress(const RankTwoTensor & stress) const
87 : {
88 : std::vector<RankFourTensor> d2;
89 0 : stress.d2symmetricEigenvalues(d2);
90 0 : return d2;
91 0 : }
92 :
93 : void
94 128824 : CappedMohrCoulombStressUpdate::preReturnMapV(const std::vector<Real> & /*trial_stress_params*/,
95 : const RankTwoTensor & stress_trial,
96 : const std::vector<Real> & /*intnl_old*/,
97 : const std::vector<Real> & /*yf*/,
98 : const RankFourTensor & Eijkl)
99 : {
100 : std::vector<Real> eigvals;
101 128824 : stress_trial.symmetricEigenvaluesEigenvectors(eigvals, _eigvecs);
102 128824 : _poissons_ratio = ElasticityTensorTools::getIsotropicPoissonsRatio(Eijkl);
103 128824 : }
104 :
105 : void
106 128824 : CappedMohrCoulombStressUpdate::setStressAfterReturnV(const RankTwoTensor & /*stress_trial*/,
107 : const std::vector<Real> & stress_params,
108 : Real /*gaE*/,
109 : const std::vector<Real> & /*intnl*/,
110 : const yieldAndFlow & /*smoothed_q*/,
111 : const RankFourTensor & /*Eijkl*/,
112 : RankTwoTensor & stress) const
113 : {
114 : // form the diagonal stress
115 128824 : stress = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0);
116 : // rotate to the original frame
117 128824 : stress = _eigvecs * stress * (_eigvecs.transpose());
118 128824 : }
119 :
120 : void
121 299336 : CappedMohrCoulombStressUpdate::yieldFunctionValuesV(const std::vector<Real> & stress_params,
122 : const std::vector<Real> & intnl,
123 : std::vector<Real> & yf) const
124 : {
125 : // intnl[0] = shear, intnl[1] = tensile
126 299336 : const Real ts = _tensile_strength.value(intnl[1]);
127 299336 : const Real cs = _compressive_strength.value(intnl[1]);
128 299336 : const Real sinphi = std::sin(_phi.value(intnl[0]));
129 299336 : const Real cohcos = _cohesion.value(intnl[0]) * std::cos(_phi.value(intnl[0]));
130 299336 : const Real smax = stress_params[2]; // largest eigenvalue
131 299336 : const Real smid = stress_params[1];
132 299336 : const Real smin = stress_params[0]; // smallest eigenvalue
133 299336 : yf[0] = smax - ts;
134 299336 : yf[1] = smid - ts;
135 299336 : yf[2] = smin - ts;
136 299336 : yf[3] = -smin - cs;
137 299336 : yf[4] = -smid - cs;
138 299336 : yf[5] = -smax - cs;
139 299336 : yf[6] = 0.5 * (smax - smin) + 0.5 * (smax + smin) * sinphi - cohcos;
140 299336 : yf[7] = 0.5 * (smid - smin) + 0.5 * (smid + smin) * sinphi - cohcos;
141 299336 : yf[8] = 0.5 * (smax - smid) + 0.5 * (smax + smid) * sinphi - cohcos;
142 299336 : yf[9] = 0.5 * (smid - smax) + 0.5 * (smid + smax) * sinphi - cohcos;
143 299336 : yf[10] = 0.5 * (smin - smid) + 0.5 * (smin + smid) * sinphi - cohcos;
144 299336 : yf[11] = 0.5 * (smin - smax) + 0.5 * (smin + smax) * sinphi - cohcos;
145 299336 : }
146 :
147 : void
148 1130408 : CappedMohrCoulombStressUpdate::computeAllQV(const std::vector<Real> & stress_params,
149 : const std::vector<Real> & intnl,
150 : std::vector<yieldAndFlow> & all_q) const
151 : {
152 1130408 : const Real ts = _tensile_strength.value(intnl[1]);
153 1130408 : const Real dts = _tensile_strength.derivative(intnl[1]);
154 1130408 : const Real cs = _compressive_strength.value(intnl[1]);
155 1130408 : const Real dcs = _compressive_strength.derivative(intnl[1]);
156 1130408 : const Real sinphi = std::sin(_phi.value(intnl[0]));
157 1130408 : const Real dsinphi = std::cos(_phi.value(intnl[0])) * _phi.derivative(intnl[0]);
158 1130408 : const Real sinpsi = std::sin(_psi.value(intnl[0]));
159 1130408 : const Real dsinpsi = std::cos(_psi.value(intnl[0])) * _psi.derivative(intnl[0]);
160 1130408 : const Real cohcos = _cohesion.value(intnl[0]) * std::cos(_phi.value(intnl[0]));
161 : const Real dcohcos =
162 1130408 : _cohesion.derivative(intnl[0]) * std::cos(_phi.value(intnl[0])) -
163 1130408 : _cohesion.value(intnl[0]) * std::sin(_phi.value(intnl[0])) * _phi.derivative(intnl[0]);
164 1130408 : const Real smax = stress_params[2]; // largest eigenvalue
165 1130408 : const Real smid = stress_params[1];
166 1130408 : const Real smin = stress_params[0]; // smallest eigenvalue
167 :
168 : // yield functions. See comment in yieldFunctionValuesV
169 1130408 : all_q[0].f = smax - ts;
170 1130408 : all_q[1].f = smid - ts;
171 1130408 : all_q[2].f = smin - ts;
172 1130408 : all_q[3].f = -smin - cs;
173 1130408 : all_q[4].f = -smid - cs;
174 1130408 : all_q[5].f = -smax - cs;
175 1130408 : all_q[6].f = 0.5 * (smax - smin) + 0.5 * (smax + smin) * sinphi - cohcos;
176 1130408 : all_q[7].f = 0.5 * (smid - smin) + 0.5 * (smid + smin) * sinphi - cohcos;
177 1130408 : all_q[8].f = 0.5 * (smax - smid) + 0.5 * (smax + smid) * sinphi - cohcos;
178 1130408 : all_q[9].f = 0.5 * (smid - smax) + 0.5 * (smid + smax) * sinphi - cohcos;
179 1130408 : all_q[10].f = 0.5 * (smin - smid) + 0.5 * (smin + smid) * sinphi - cohcos;
180 1130408 : all_q[11].f = 0.5 * (smin - smax) + 0.5 * (smin + smax) * sinphi - cohcos;
181 :
182 : // d(yield function)/d(stress_params)
183 14695304 : for (unsigned yf = 0; yf < _num_yf; ++yf)
184 54259584 : for (unsigned a = 0; a < _num_sp; ++a)
185 40694688 : all_q[yf].df[a] = 0.0;
186 1130408 : all_q[0].df[2] = 1.0;
187 1130408 : all_q[1].df[1] = 1.0;
188 1130408 : all_q[2].df[0] = 1.0;
189 1130408 : all_q[3].df[0] = -1.0;
190 1130408 : all_q[4].df[1] = -1.0;
191 1130408 : all_q[5].df[2] = -1.0;
192 1130408 : all_q[6].df[2] = 0.5 * (1.0 + sinphi);
193 1130408 : all_q[6].df[0] = 0.5 * (-1.0 + sinphi);
194 1130408 : all_q[7].df[1] = 0.5 * (1.0 + sinphi);
195 1130408 : all_q[7].df[0] = 0.5 * (-1.0 + sinphi);
196 1130408 : all_q[8].df[2] = 0.5 * (1.0 + sinphi);
197 1130408 : all_q[8].df[1] = 0.5 * (-1.0 + sinphi);
198 1130408 : all_q[9].df[1] = 0.5 * (1.0 + sinphi);
199 1130408 : all_q[9].df[2] = 0.5 * (-1.0 + sinphi);
200 1130408 : all_q[10].df[0] = 0.5 * (1.0 + sinphi);
201 1130408 : all_q[10].df[1] = 0.5 * (-1.0 + sinphi);
202 1130408 : all_q[11].df[0] = 0.5 * (1.0 + sinphi);
203 1130408 : all_q[11].df[2] = 0.5 * (-1.0 + sinphi);
204 :
205 : // d(yield function)/d(intnl)
206 7912856 : for (unsigned i = 0; i < 6; ++i)
207 6782448 : all_q[i].df_di[0] = 0.0;
208 1130408 : all_q[0].df_di[1] = all_q[1].df_di[1] = all_q[2].df_di[1] = -dts;
209 1130408 : all_q[3].df_di[1] = all_q[4].df_di[1] = all_q[5].df_di[1] = -dcs;
210 7912856 : for (unsigned i = 6; i < 12; ++i)
211 6782448 : all_q[i].df_di[1] = 0.0;
212 1130408 : all_q[6].df_di[0] = 0.5 * (smax + smin) * dsinphi - dcohcos;
213 1130408 : all_q[7].df_di[0] = 0.5 * (smid + smin) * dsinphi - dcohcos;
214 1130408 : all_q[8].df_di[0] = 0.5 * (smax + smid) * dsinphi - dcohcos;
215 1130408 : all_q[9].df_di[0] = 0.5 * (smid + smax) * dsinphi - dcohcos;
216 1130408 : all_q[10].df_di[0] = 0.5 * (smin + smid) * dsinphi - dcohcos;
217 1130408 : all_q[11].df_di[0] = 0.5 * (smin + smax) * dsinphi - dcohcos;
218 :
219 : // the flow potential is just the yield function with phi->psi
220 : // d(flow potential)/d(stress_params)
221 7912856 : for (unsigned yf = 0; yf < 6; ++yf)
222 27129792 : for (unsigned a = 0; a < _num_sp; ++a)
223 20347344 : all_q[yf].dg[a] = all_q[yf].df[a];
224 1130408 : all_q[6].dg[2] = all_q[7].dg[1] = all_q[8].dg[2] = all_q[9].dg[1] = all_q[10].dg[0] =
225 1130408 : all_q[11].dg[0] = 0.5 * (1.0 + sinpsi);
226 1130408 : all_q[6].dg[0] = all_q[7].dg[0] = all_q[8].dg[1] = all_q[9].dg[2] = all_q[10].dg[1] =
227 1130408 : all_q[11].dg[2] = 0.5 * (-1.0 + sinpsi);
228 :
229 : // d(flow potential)/d(stress_params)/d(intnl)
230 14695304 : for (unsigned yf = 0; yf < _num_yf; ++yf)
231 54259584 : for (unsigned a = 0; a < _num_sp; ++a)
232 122084064 : for (unsigned i = 0; i < _num_intnl; ++i)
233 81389376 : all_q[yf].d2g_di[a][i] = 0.0;
234 1130408 : all_q[6].d2g_di[2][0] = all_q[7].d2g_di[1][0] = all_q[8].d2g_di[2][0] = all_q[9].d2g_di[1][0] =
235 1130408 : all_q[10].d2g_di[0][0] = all_q[11].d2g_di[0][0] = 0.5 * dsinpsi;
236 1130408 : all_q[6].d2g_di[0][0] = all_q[7].d2g_di[0][0] = all_q[8].d2g_di[1][0] = all_q[9].d2g_di[2][0] =
237 1130408 : all_q[10].d2g_di[1][0] = all_q[11].d2g_di[2][0] = 0.5 * dsinpsi;
238 :
239 : // d(flow potential)/d(stress_params)/d(stress_params)
240 14695304 : for (unsigned yf = 0; yf < _num_yf; ++yf)
241 54259584 : for (unsigned a = 0; a < _num_sp; ++a)
242 162778752 : for (unsigned b = 0; b < _num_sp; ++b)
243 122084064 : all_q[yf].d2g[a][b] = 0.0;
244 1130408 : }
245 :
246 : void
247 128824 : CappedMohrCoulombStressUpdate::setEffectiveElasticity(const RankFourTensor & Eijkl)
248 : {
249 : // Eijkl is required to be isotropic, so we can use the
250 : // frame where stress is diagonal
251 515296 : for (unsigned a = 0; a < _num_sp; ++a)
252 1545888 : for (unsigned b = 0; b < _num_sp; ++b)
253 1159416 : _Eij[a][b] = Eijkl(a, a, b, b);
254 128824 : _En = _Eij[2][2];
255 128824 : const Real denom = _Eij[0][0] * (_Eij[0][0] + _Eij[0][1]) - 2 * Utility::pow<2>(_Eij[0][1]);
256 515296 : for (unsigned a = 0; a < _num_sp; ++a)
257 : {
258 386472 : _Cij[a][a] = (_Eij[0][0] + _Eij[0][1]) / denom;
259 772944 : for (unsigned b = 0; b < a; ++b)
260 386472 : _Cij[a][b] = _Cij[b][a] = -_Eij[0][1] / denom;
261 : }
262 128824 : }
263 :
264 : void
265 140536 : CappedMohrCoulombStressUpdate::initializeVarsV(const std::vector<Real> & trial_stress_params,
266 : const std::vector<Real> & intnl_old,
267 : std::vector<Real> & stress_params,
268 : Real & gaE,
269 : std::vector<Real> & intnl) const
270 : {
271 140536 : if (!_perfect_guess)
272 : {
273 0 : for (unsigned i = 0; i < _num_sp; ++i)
274 0 : stress_params[i] = trial_stress_params[i];
275 0 : gaE = 0.0;
276 : }
277 : else
278 : {
279 140536 : const Real ts = _tensile_strength.value(intnl_old[1]);
280 140536 : const Real cs = _compressive_strength.value(intnl_old[1]);
281 140536 : const Real sinphi = std::sin(_phi.value(intnl_old[0]));
282 140536 : const Real cohcos = _cohesion.value(intnl_old[0]) * std::cos(_phi.value(intnl_old[0]));
283 :
284 140536 : const Real trial_tensile_yf = trial_stress_params[2] - ts;
285 140536 : const Real trial_compressive_yf = -trial_stress_params[0] - cs;
286 140536 : const Real trial_mc_yf = 0.5 * (trial_stress_params[2] - trial_stress_params[0]) +
287 140536 : 0.5 * (trial_stress_params[2] + trial_stress_params[0]) * sinphi -
288 140536 : cohcos;
289 :
290 : /**
291 : * For CappedMohrCoulomb there are a number of possibilities for the
292 : * returned stress configuration:
293 : * - return to the Tensile yield surface to its plane
294 : * - return to the Tensile yield surface to its max=mid edge
295 : * - return to the Tensile yield surface to its tip
296 : * - return to the Compressive yield surface to its plane
297 : * - return to the Compressive yield surface to its max=mid edge
298 : * - return to the Compressive yield surface to its tip
299 : * - return to the Mohr-Coulomb yield surface to its plane
300 : * - return to the Mohr-Coulomb yield surface to its max=mid edge
301 : * - return to the Mohr-Coulomb yield surface to its mid=min edge
302 : * - return to the Mohr-Coulomb yield surface to its tip
303 : * - return to the edge between Tensile and Mohr-Coulomb
304 : * - return to the edge between Tensile and Mohr-Coulomb, at max=mid point
305 : * - return to the edge between Tensile and Mohr-Coulomb, at mid=min point
306 : * - return to the edge between Compressive and Mohr-Coulomb
307 : * - return to the edge between Compressive and Mohr-Coulomb, at max=mid point
308 : * - return to the edge between Compressive and Mohr-Coulomb, at mid=min point
309 : * Which of these is more prelevant depends on the parameters
310 : * tensile strength, compressive strength, cohesion, etc.
311 : * I simply check each possibility until i find one that works.
312 : * _shifter is used to avoid equal eigenvalues
313 : */
314 :
315 : bool found_solution = false;
316 :
317 140536 : if (trial_tensile_yf <= _f_tol && trial_compressive_yf <= _f_tol && trial_mc_yf <= _f_tol)
318 : {
319 : // this is needed because although we know the smoothed yield function is
320 : // positive, the individual yield functions may not be
321 1152 : for (unsigned i = 0; i < _num_sp; ++i)
322 864 : stress_params[i] = trial_stress_params[i];
323 288 : gaE = 0.0;
324 : found_solution = true;
325 : }
326 :
327 140536 : const bool tensile_possible = (ts < cohcos / sinphi); // tensile chops MC tip
328 : const bool mc_tip_possible = (cohcos / sinphi < ts); // MC tip pokes through tensile
329 140536 : const bool mc_impossible = (0.5 * (ts + cs) + 0.5 * (ts - cs) * sinphi - cohcos <
330 140536 : _f_tol); // MC outside tensile and compressive
331 :
332 140536 : const Real sinpsi = std::sin(_psi.value(intnl_old[0]));
333 140536 : const Real halfplus = 0.5 + 0.5 * sinpsi;
334 140536 : const Real neghalfplus = -0.5 + 0.5 * sinpsi;
335 :
336 140536 : if (!found_solution && tensile_possible && trial_tensile_yf > _f_tol &&
337 11784 : (trial_compressive_yf <= _f_tol || (trial_compressive_yf > _f_tol && mc_impossible)))
338 : {
339 : // try pure tensile failure, return to the plane
340 : // This involves solving yf[0] = 0 and the three flow-direction equations
341 : // Don't try this if there is compressive failure, since returning to
342 : // the tensile yield surface will only make compressive failure worse
343 42552 : const Real ga = (trial_stress_params[2] - ts) / _Eij[2][2];
344 42552 : stress_params[2] = ts; // largest eigenvalue
345 42552 : stress_params[1] = trial_stress_params[1] - ga * _Eij[1][2];
346 42552 : stress_params[0] = trial_stress_params[0] - ga * _Eij[0][2];
347 :
348 : // if we have to return to the edge, or tip, do that
349 : Real dist_mod = 1.0;
350 42552 : const Real to_subtract1 = stress_params[1] - (ts - 0.5 * _shifter);
351 42552 : if (to_subtract1 > 0.0)
352 : {
353 12488 : dist_mod += Utility::pow<2>(to_subtract1 / (trial_stress_params[2] - ts));
354 12488 : stress_params[1] -= to_subtract1;
355 : }
356 42552 : const Real to_subtract0 = stress_params[0] - (ts - _shifter);
357 42552 : if (to_subtract0 > 0.0)
358 : {
359 6776 : dist_mod += Utility::pow<2>(to_subtract0 / (trial_stress_params[2] - ts));
360 6776 : stress_params[0] -= to_subtract0;
361 : }
362 42552 : if (mc_impossible) // might have to shift up to the compressive yield surface
363 : {
364 33920 : const Real to_add0 = -stress_params[0] - cs;
365 33920 : if (to_add0 > 0.0)
366 : {
367 0 : dist_mod += Utility::pow<2>(to_add0 / (trial_stress_params[2] - ts));
368 0 : stress_params[0] += to_add0;
369 : }
370 33920 : const Real to_add1 = -cs + 0.5 * _shifter - stress_params[1];
371 33920 : if (to_add1 > 0.0)
372 : {
373 0 : dist_mod += Utility::pow<2>(to_add1 / (trial_stress_params[2] - ts));
374 0 : stress_params[1] += to_add1;
375 : }
376 : }
377 :
378 42552 : const Real new_compressive_yf = -stress_params[0] - cs;
379 42552 : const Real new_mc_yf = 0.5 * (stress_params[2] - stress_params[0]) +
380 42552 : 0.5 * (stress_params[2] + stress_params[0]) * sinphi - cohcos;
381 42552 : if (new_mc_yf <= _f_tol && new_compressive_yf <= _f_tol)
382 : {
383 36552 : gaE = std::sqrt(dist_mod) * (trial_stress_params[2] - stress_params[2]);
384 : found_solution = true;
385 : }
386 : }
387 103984 : if (!found_solution && trial_compressive_yf > _f_tol &&
388 11784 : (trial_tensile_yf <= _f_tol || (trial_tensile_yf > _f_tol && mc_impossible)))
389 : {
390 : // try pure compressive failure
391 : // Don't try this if there is tensile failure, since returning to
392 : // the compressive yield surface will only make tensile failure worse
393 40264 : const Real ga = (trial_stress_params[0] + cs) / _Eij[0][0]; // this is negative
394 40264 : stress_params[0] = -cs;
395 40264 : stress_params[1] = trial_stress_params[1] - ga * _Eij[1][0];
396 40264 : stress_params[2] = trial_stress_params[2] - ga * _Eij[2][0];
397 :
398 : // if we have to return to the edge, or tip, do that
399 : Real dist_mod = 1.0;
400 40264 : const Real to_add1 = -cs + 0.5 * _shifter - stress_params[1];
401 40264 : if (to_add1 > 0.0)
402 : {
403 12192 : dist_mod += Utility::pow<2>(to_add1 / (trial_stress_params[0] + cs));
404 12192 : stress_params[1] += to_add1;
405 : }
406 40264 : const Real to_add2 = -cs + _shifter - stress_params[2];
407 40264 : if (to_add2 > 0.0)
408 : {
409 6928 : dist_mod += Utility::pow<2>(to_add2 / (trial_stress_params[0] + cs));
410 6928 : stress_params[2] += to_add2;
411 : }
412 40264 : if (mc_impossible) // might have to shift down to the tensile yield surface
413 : {
414 30600 : const Real to_subtract2 = stress_params[2] - ts;
415 30600 : if (to_subtract2 > 0.0)
416 : {
417 0 : dist_mod += Utility::pow<2>(to_subtract2 / (trial_stress_params[0] + cs));
418 0 : stress_params[2] -= to_subtract2;
419 : }
420 30600 : const Real to_subtract1 = stress_params[1] - (ts - 0.5 * _shifter);
421 30600 : if (to_subtract1 > 0.0)
422 : {
423 0 : dist_mod += Utility::pow<2>(to_subtract1 / (trial_stress_params[0] + cs));
424 0 : stress_params[1] -= to_subtract1;
425 : }
426 : }
427 :
428 40264 : const Real new_tensile_yf = stress_params[2] - ts;
429 40264 : const Real new_mc_yf = 0.5 * (stress_params[2] - stress_params[0]) +
430 40264 : 0.5 * (stress_params[2] + stress_params[0]) * sinphi - cohcos;
431 40264 : if (new_mc_yf <= _f_tol && new_tensile_yf <= _f_tol)
432 : {
433 35368 : gaE = std::sqrt(dist_mod) * (-trial_stress_params[0] + stress_params[0]);
434 : found_solution = true;
435 : }
436 : }
437 140536 : if (!found_solution && !mc_impossible && trial_mc_yf > _f_tol)
438 : {
439 : // try pure shear failure, return to the plane
440 : // This involves solving yf[6]=0 and the three flow-direction equations
441 68328 : const Real ga = ((trial_stress_params[2] - trial_stress_params[0]) +
442 68328 : (trial_stress_params[2] + trial_stress_params[0]) * sinphi - 2.0 * cohcos) /
443 68328 : (_Eij[2][2] - _Eij[0][2] + (_Eij[2][2] + _Eij[0][2]) * sinpsi * sinphi);
444 68328 : stress_params[2] =
445 68328 : trial_stress_params[2] - ga * (_Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus);
446 68328 : stress_params[1] = trial_stress_params[1] - ga * _Eij[1][0] * sinpsi;
447 68328 : stress_params[0] =
448 68328 : trial_stress_params[0] - ga * (_Eij[0][0] * neghalfplus + _Eij[0][2] * halfplus);
449 68328 : const Real f7 = 0.5 * (stress_params[1] - stress_params[0]) +
450 68328 : 0.5 * (stress_params[1] + stress_params[0]) * sinphi - cohcos;
451 68328 : const Real f8 = 0.5 * (stress_params[2] - stress_params[1]) +
452 68328 : 0.5 * (stress_params[2] + stress_params[1]) * sinphi - cohcos;
453 68328 : const Real new_tensile_yf = stress_params[2] - ts;
454 68328 : const Real new_compressive_yf = -stress_params[0] - cs;
455 :
456 68328 : if (f7 <= _f_tol && f8 <= _f_tol && new_tensile_yf <= _f_tol && new_compressive_yf <= _f_tol)
457 : {
458 21624 : gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
459 21624 : (stress_params[2] - stress_params[0])) /
460 21624 : (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2];
461 : found_solution = true;
462 : }
463 : }
464 140536 : if (!found_solution && !mc_impossible && trial_mc_yf > _f_tol)
465 : {
466 : // Try return to the max=mid MC line.
467 : // To return to the max=mid line, we need to solve f6 = 0 = f7 and
468 : // the three flow-direction equations. In the flow-direction equations
469 : // there are two plasticity multipliers, which i call ga6 and ga7,
470 : // corresponding to the amounts of strain normal to the f6 and f7
471 : // directions, respectively.
472 : // So:
473 : // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga7 Emax,a dg7/dSa
474 : // = Smax^trial - ga6 cmax6 - ga7 cmax7 (with cmax6 and cmax7 evaluated below)
475 : // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga7 Emid,a dg7/dSa
476 : // = Smid^trial - ga6 cmid6 - ga7 cmid7
477 : // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga7 Emin,a dg7/dSa
478 : // = Smin^trial - ga6 cmin6 - ga7 cmin7
479 46704 : const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
480 46704 : const Real cmax7 = _Eij[2][1] * halfplus + _Eij[2][0] * neghalfplus;
481 : // const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
482 : // const Real cmid7 = _Eij[1][1] * halfplus + _Eij[1][0] * neghalfplus;
483 46704 : const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
484 46704 : const Real cmin7 = _Eij[0][1] * halfplus + _Eij[0][0] * neghalfplus;
485 : // Substituting these into f6 = 0 yields
486 : // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga7 (0.5(cmax6 -
487 : // cmin6) + 0.5(cmax6 + cmin6)sinphi) = f6_trial - ga6 c6 - ga7 c7, where
488 46704 : const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
489 46704 : const Real c7 = 0.5 * (cmax7 - cmin7) + 0.5 * (cmax7 + cmin7) * sinphi;
490 : // It isn't too hard to check that the other equation is
491 : // 0 = f7_trial - ga6 c7 - ga7 c6
492 : // These equations may be inverted to yield
493 46704 : if (c6 != c7)
494 : {
495 : const Real f6_trial = trial_mc_yf;
496 46704 : const Real f7_trial = 0.5 * (trial_stress_params[1] - trial_stress_params[0]) +
497 46704 : 0.5 * (trial_stress_params[1] + trial_stress_params[0]) * sinphi -
498 46704 : cohcos;
499 46704 : const Real descr = Utility::pow<2>(c6) - Utility::pow<2>(c7);
500 46704 : Real ga6 = (c6 * f6_trial - c7 * f7_trial) / descr;
501 46704 : Real ga7 = (-c7 * f6_trial + c6 * f7_trial) / descr;
502 : // and finally
503 46704 : stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga7 * cmax7;
504 46704 : stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga7 * cmin7;
505 46704 : stress_params[1] = stress_params[2] - 0.5 * _shifter;
506 :
507 46704 : Real f8 = 0.5 * (stress_params[2] - stress_params[1]) +
508 46704 : 0.5 * (stress_params[2] + stress_params[1]) * sinphi - cohcos;
509 :
510 46704 : if (mc_tip_possible && f8 > _f_tol)
511 : {
512 7432 : stress_params[2] = cohcos / sinphi;
513 7432 : stress_params[1] = stress_params[2] - 0.5 * _shifter;
514 7432 : stress_params[0] = stress_params[2] - _shifter;
515 : f8 = 0.0;
516 : ga6 = 1.0;
517 : ga7 = 1.0;
518 : }
519 :
520 46704 : const Real new_tensile_yf = stress_params[2] - ts;
521 46704 : const Real new_compressive_yf = -stress_params[0] - cs;
522 :
523 46704 : if (f8 <= _f_tol && new_tensile_yf <= _f_tol && new_compressive_yf <= _f_tol &&
524 24216 : ga6 >= 0.0 && ga7 >= 0.0)
525 : {
526 12592 : gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
527 12592 : (stress_params[2] - stress_params[0])) /
528 12592 : (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2];
529 : found_solution = true;
530 : }
531 : }
532 : }
533 140536 : if (!found_solution && !mc_impossible && trial_mc_yf > _f_tol)
534 : {
535 : // Try return to the mid=min line.
536 : // To return to the mid=min line, we need to solve f6 = 0 = f8 and
537 : // the three flow-direction equations. In the flow-direction equations
538 : // there are two plasticity multipliers, which i call ga6 and ga8,
539 : // corresponding to the amounts of strain normal to the f6 and f8
540 : // directions, respectively.
541 : // So:
542 : // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga8 Emax,a dg8/dSa
543 : // = Smax^trial - ga6 cmax6 - ga8 cmax8 (with cmax6 and cmax8 evaluated below)
544 : // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga8 Emid,a dg8/dSa
545 : // = Smid^trial - ga6 cmid6 - ga8 cmid8
546 : // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga8 Emin,a dg8/dSa
547 : // = Smin^trial - ga6 cmin6 - ga8 cmin8
548 34112 : const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
549 34112 : const Real cmax8 = _Eij[2][2] * halfplus + _Eij[2][1] * neghalfplus;
550 : // const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
551 : // const Real cmid8 = _Eij[1][2] * halfplus + _Eij[1][1] * neghalfplus;
552 34112 : const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
553 34112 : const Real cmin8 = _Eij[0][2] * halfplus + _Eij[0][1] * neghalfplus;
554 : // Substituting these into f6 = 0 yields
555 : // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga8 (0.5(cmax6 -
556 : // cmin6) + 0.5(cmax6 + cmin6)sinphi) = f6_trial - ga6 c6 - ga8 c8, where
557 34112 : const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
558 34112 : const Real c8 = 0.5 * (cmax8 - cmin8) + 0.5 * (cmax8 + cmin8) * sinphi;
559 : // It isn't too hard to check that the other equation is
560 : // 0 = f8_trial - ga6 c8 - ga8 c6
561 : // These equations may be inverted to yield
562 34112 : if (c6 != c8)
563 : {
564 : const Real f6_trial = trial_mc_yf;
565 34112 : const Real f8_trial = 0.5 * (trial_stress_params[2] - trial_stress_params[1]) +
566 34112 : 0.5 * (trial_stress_params[2] + trial_stress_params[1]) * sinphi -
567 34112 : cohcos;
568 34112 : const Real descr = Utility::pow<2>(c6) - Utility::pow<2>(c8);
569 34112 : Real ga6 = (c6 * f6_trial - c8 * f8_trial) / descr;
570 34112 : Real ga8 = (-c8 * f6_trial + c6 * f8_trial) / descr;
571 : // and finally
572 34112 : stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga8 * cmax8;
573 34112 : stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga8 * cmin8;
574 34112 : stress_params[1] = stress_params[0] + 0.5 * _shifter;
575 :
576 34112 : Real f7 = 0.5 * (stress_params[1] - stress_params[0]) +
577 34112 : 0.5 * (stress_params[1] + stress_params[0]) * sinphi - cohcos;
578 :
579 34112 : if (mc_tip_possible && f7 > _f_tol)
580 : {
581 0 : stress_params[2] = cohcos / sinphi;
582 0 : stress_params[1] = stress_params[2] - 0.5 * _shifter;
583 0 : stress_params[0] = stress_params[2] - _shifter;
584 : f7 = 0.0;
585 : ga6 = 1.0;
586 : ga8 = 1.0;
587 : }
588 :
589 34112 : const Real new_tensile_yf = stress_params[2] - ts;
590 34112 : const Real new_compressive_yf = -stress_params[0] - cs;
591 :
592 34112 : if (f7 <= _f_tol && new_tensile_yf <= _f_tol && new_compressive_yf <= _f_tol &&
593 11936 : ga6 >= 0.0 && ga8 >= 0.0)
594 : {
595 11624 : gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
596 11624 : (stress_params[2] - stress_params[0])) /
597 11624 : (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2];
598 : found_solution = true;
599 : }
600 : }
601 : }
602 140536 : if (!found_solution && !mc_impossible && tensile_possible && trial_tensile_yf > _f_tol)
603 : {
604 : // Return to the line where yf[0] = 0 = yf[6].
605 : // To return to this line, we need to solve f0 = 0 = f6 and
606 : // the three flow-direction equations. In the flow-direction equations
607 : // there are two plasticity multipliers, which i call ga0 and ga6
608 : // corresponding to the amounts of strain normal to the f0 and f6
609 : // directions, respectively.
610 : // So:
611 : // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga0 Emax,a dg0/dSa
612 : // = Smax^trial - ga6 cmax6 - ga0 cmax0 (with cmax6 and cmax0 evaluated below)
613 : // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga0 Emid,a dg0/dSa
614 : // = Smid^trial - ga6 cmid6 - ga0 cmid0
615 : // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga0 Emin,a dg0/dSa
616 : // = Smin^trial - ga6 cmin6 - ga0 cmin0
617 17592 : const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
618 : const Real cmax0 = _Eij[2][2];
619 17592 : const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
620 : const Real cmid0 = _Eij[1][2];
621 17592 : const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
622 : const Real cmin0 = _Eij[0][2];
623 : // Substituting these into f6 = 0 yields
624 : // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga0 (0.5(cmax0 -
625 : // cmin0) + 0.5(cmax0 + cmin0)sinphi) = f6_trial - ga6 c6 - ga0 c0, where
626 17592 : const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
627 17592 : const Real c0 = 0.5 * (cmax0 - cmin0) + 0.5 * (cmax0 + cmin0) * sinphi;
628 : // Substituting these into f0 = 0 yields
629 : // 0 = f0_trial - ga6 cmax6 - ga0 cmax0
630 : // These equations may be inverted to yield the following
631 17592 : const Real descr = c0 * cmax6 - c6 * cmax0;
632 17592 : if (descr != 0.0)
633 : {
634 17592 : const Real ga0 = (-c6 * trial_tensile_yf + cmax6 * trial_mc_yf) / descr;
635 17592 : const Real ga6 = (c0 * trial_tensile_yf - cmax0 * trial_mc_yf) / descr;
636 17592 : stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga0 * cmax0;
637 17592 : stress_params[1] = trial_stress_params[1] - ga6 * cmid6 - ga0 * cmid0;
638 17592 : stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga0 * cmin0;
639 :
640 : // enforce physicality (go to corners if necessary)
641 17592 : stress_params[0] =
642 17592 : std::min(stress_params[0],
643 17592 : stress_params[2] - _shifter); // account for poor choice from user
644 : // if goto_corner then the max(min()) in the subsequent line will force the solution to lie
645 : // at the corner where max = mid = tensile. This means the signs of ga0 and ga6 become
646 : // irrelevant in the check below
647 17592 : const bool goto_corner = (stress_params[1] >= stress_params[2] - 0.5 * _shifter);
648 35184 : stress_params[1] = std::max(std::min(stress_params[1], stress_params[2] - 0.5 * _shifter),
649 17592 : stress_params[0] + 0.5 * _shifter);
650 :
651 17592 : const Real new_compressive_yf = -stress_params[0] - cs;
652 17592 : if (new_compressive_yf <= _f_tol &&
653 17592 : (goto_corner || (ga0 >= 0.0 && ga6 >= 0.0))) // enforce ga>=0 unless going to a corner
654 : {
655 6744 : gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
656 6744 : (stress_params[2] - stress_params[0])) /
657 6744 : (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2] +
658 6744 : (trial_stress_params[2] - stress_params[2]);
659 : found_solution = true;
660 : }
661 : }
662 : }
663 140536 : if (!found_solution && !mc_impossible)
664 : {
665 : // Return to the line where yf[3] = 0 = yf[6].
666 : // To return to this line, we need to solve f3 = 0 = f6 and
667 : // the three flow-direction equations. In the flow-direction equations
668 : // there are two plasticity multipliers, which i call ga3 and ga6
669 : // corresponding to the amounts of strain normal to the f3 and f6
670 : // directions, respectively.
671 : // So:
672 : // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga3 Emax,a dg3/dSa
673 : // = Smax^trial - ga6 cmax6 - ga3 cmax3 (with cmax6 and cmax3 evaluated below)
674 : // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga3 Emid,a dg3/dSa
675 : // = Smid^trial - ga6 cmid6 - ga3 cmid3
676 : // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga3 Emin,a dg3/dSa
677 : // = Smin^trial - ga6 cmin6 - ga3 cmin3
678 15744 : const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
679 15744 : const Real cmax3 = -_Eij[2][0];
680 15744 : const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
681 15744 : const Real cmid3 = -_Eij[1][0];
682 15744 : const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
683 15744 : const Real cmin3 = -_Eij[0][0];
684 : // Substituting these into f6 = 0 yields
685 : // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga3 (0.5(cmax3 -
686 : // cmin3) + 0.5(cmax3 + cmin3)sinphi) = f6_trial - ga6 c6 - ga3 c3, where
687 15744 : const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
688 15744 : const Real c3 = 0.5 * (cmax3 - cmin3) + 0.5 * (cmax3 + cmin3) * sinphi;
689 : // Substituting these into f3 = 0 yields
690 : // 0 = - f3_trial - ga6 cmin6 - ga3 cmin3
691 : // These equations may be inverted to yield the following
692 15744 : const Real descr = c3 * cmin6 - c6 * cmin3;
693 15744 : if (descr != 0.0)
694 : {
695 15744 : const Real ga3 = (c6 * trial_compressive_yf + cmin6 * trial_mc_yf) / descr;
696 15744 : const Real ga6 = (-c3 * trial_compressive_yf - cmin3 * trial_mc_yf) / descr;
697 15744 : stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga3 * cmax3;
698 15744 : stress_params[1] = trial_stress_params[1] - ga6 * cmid6 - ga3 * cmid3;
699 15744 : stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga3 * cmin3;
700 :
701 15744 : const Real new_tensile_yf = stress_params[2] - ts;
702 15744 : stress_params[2] =
703 15744 : std::max(stress_params[2],
704 15744 : stress_params[0] + _shifter); // account for poor choice from user
705 31488 : stress_params[1] = std::max(std::min(stress_params[1], stress_params[2] - 0.5 * _shifter),
706 15744 : stress_params[0] + 0.5 * _shifter);
707 :
708 15744 : if (new_tensile_yf <= _f_tol && ga6 >= 0.0)
709 : {
710 15744 : gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
711 15744 : (stress_params[2] - stress_params[0])) /
712 15744 : (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2] +
713 15744 : (-trial_stress_params[0] - stress_params[0]);
714 :
715 : found_solution = true;
716 : }
717 : }
718 : }
719 124792 : if (!found_solution)
720 : {
721 : // Cannot find an acceptable initialisation
722 0 : for (unsigned i = 0; i < _num_sp; ++i)
723 0 : stress_params[i] = trial_stress_params[i];
724 0 : gaE = 0.0;
725 0 : mooseWarning("CappedMohrCoulombStressUpdate cannot initialize from max = ",
726 : stress_params[2],
727 : " mid = ",
728 : stress_params[1],
729 : " min = ",
730 : stress_params[0]);
731 : }
732 : }
733 140536 : setIntnlValuesV(trial_stress_params, stress_params, intnl_old, intnl);
734 140536 : }
735 :
736 : void
737 1270944 : CappedMohrCoulombStressUpdate::setIntnlValuesV(const std::vector<Real> & trial_stress_params,
738 : const std::vector<Real> & current_stress_params,
739 : const std::vector<Real> & intnl_old,
740 : std::vector<Real> & intnl) const
741 : {
742 : // intnl[0] = shear, intnl[1] = tensile
743 1270944 : const Real smax = current_stress_params[2]; // largest eigenvalue
744 1270944 : const Real smin = current_stress_params[0]; // smallest eigenvalue
745 1270944 : const Real trial_smax = trial_stress_params[2]; // largest eigenvalue
746 1270944 : const Real trial_smin = trial_stress_params[0]; // smallest eigenvalue
747 1270944 : const Real ga_shear = ((trial_smax - trial_smin) - (smax - smin)) / (_Eij[2][2] - _Eij[0][2]);
748 1270944 : intnl[0] = intnl_old[0] + ga_shear;
749 1270944 : const Real sinpsi = std::sin(_psi.value(intnl[0]));
750 1270944 : const Real prefactor = (_Eij[2][2] + _Eij[0][2]) * sinpsi;
751 1270944 : const Real shear_correction = prefactor * ga_shear;
752 1270944 : const Real ga_tensile = (1.0 - _poissons_ratio) *
753 1270944 : ((trial_smax + trial_smin) - (smax + smin) - shear_correction) /
754 1270944 : _Eij[2][2];
755 1270944 : intnl[1] = intnl_old[1] + ga_tensile;
756 1270944 : }
757 :
758 : void
759 611884 : CappedMohrCoulombStressUpdate::setIntnlDerivativesV(const std::vector<Real> & trial_stress_params,
760 : const std::vector<Real> & current_stress_params,
761 : const std::vector<Real> & intnl,
762 : std::vector<std::vector<Real>> & dintnl) const
763 : {
764 : // intnl[0] = shear, intnl[1] = tensile
765 611884 : const Real smax = current_stress_params[2]; // largest eigenvalue
766 611884 : const Real smin = current_stress_params[0]; // smallest eigenvalue
767 611884 : const Real trial_smax = trial_stress_params[2]; // largest eigenvalue
768 611884 : const Real trial_smin = trial_stress_params[0]; // smallest eigenvalue
769 611884 : const Real ga_shear = ((trial_smax - trial_smin) - (smax - smin)) / (_Eij[2][2] - _Eij[0][2]);
770 : const std::vector<Real> dga_shear = {
771 611884 : 1.0 / (_Eij[2][2] - _Eij[0][2]), 0.0, -1.0 / (_Eij[2][2] - _Eij[0][2])};
772 : // intnl[0] = intnl_old[0] + ga_shear;
773 2447536 : for (std::size_t i = 0; i < _num_sp; ++i)
774 1835652 : dintnl[0][i] = dga_shear[i];
775 :
776 611884 : const Real sinpsi = std::sin(_psi.value(intnl[0]));
777 611884 : const Real dsinpsi_di0 = _psi.derivative(intnl[0]) * std::cos(_psi.value(intnl[0]));
778 :
779 611884 : const Real prefactor = (_Eij[2][2] + _Eij[0][2]) * sinpsi;
780 611884 : const Real dprefactor_di0 = (_Eij[2][2] + _Eij[0][2]) * dsinpsi_di0;
781 : // const Real shear_correction = prefactor * ga_shear;
782 611884 : std::vector<Real> dshear_correction(_num_sp);
783 2447536 : for (std::size_t i = 0; i < _num_sp; ++i)
784 1835652 : dshear_correction[i] = prefactor * dga_shear[i] + dprefactor_di0 * dintnl[0][i] * ga_shear;
785 : // const Real ga_tensile = (1 - _poissons_ratio) * ((trial_smax + trial_smin) - (smax + smin) -
786 : // shear_correction) /
787 : // _Eij[2][2];
788 : // intnl[1] = intnl_old[1] + ga_tensile;
789 2447536 : for (std::size_t i = 0; i < _num_sp; ++i)
790 1835652 : dintnl[1][i] = -(1.0 - _poissons_ratio) * dshear_correction[i] / _Eij[2][2];
791 611884 : dintnl[1][2] += -(1.0 - _poissons_ratio) / _Eij[2][2];
792 611884 : dintnl[1][0] += -(1.0 - _poissons_ratio) / _Eij[2][2];
793 611884 : }
794 :
795 : void
796 928 : CappedMohrCoulombStressUpdate::consistentTangentOperatorV(
797 : const RankTwoTensor & stress_trial,
798 : const std::vector<Real> & trial_stress_params,
799 : const RankTwoTensor & /*stress*/,
800 : const std::vector<Real> & stress_params,
801 : Real /*gaE*/,
802 : const yieldAndFlow & /*smoothed_q*/,
803 : const RankFourTensor & elasticity_tensor,
804 : bool compute_full_tangent_operator,
805 : const std::vector<std::vector<Real>> & dvar_dtrial,
806 : RankFourTensor & cto)
807 : {
808 928 : cto = elasticity_tensor;
809 928 : if (!compute_full_tangent_operator)
810 0 : return;
811 :
812 : // dvar_dtrial has been computed already, so
813 : // d(stress)/d(trial_stress) = d(eigvecs * stress_params * eigvecs.transpose())/d(trial_stress)
814 : // eigvecs is a rotation matrix, rot(i, j) = e_j(i) = i^th component of j^th eigenvector
815 : // d(rot_ij)/d(stress_kl) = d(e_j(i))/d(stress_kl)
816 : // = sum_a 0.5 * e_a(i) * (e_a(k)e_j(l) + e_a(l)e_j(k)) / (la_j - la_a)
817 : // = sum_a 0.5 * rot(i,a) * (rot(k,a)rot(l,j) + rot(l,a)*rot(k,j)) / (la_j - la_a)
818 928 : RankFourTensor drot_dstress;
819 3712 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
820 11136 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
821 33408 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
822 100224 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
823 300672 : for (unsigned a = 0; a < _num_sp; ++a)
824 : {
825 225504 : if (trial_stress_params[a] == trial_stress_params[j])
826 75168 : continue;
827 150336 : drot_dstress(i, j, k, l) +=
828 150336 : 0.5 * _eigvecs(i, a) *
829 150336 : (_eigvecs(k, a) * _eigvecs(l, j) + _eigvecs(l, a) * _eigvecs(k, j)) /
830 150336 : (trial_stress_params[j] - trial_stress_params[a]);
831 : }
832 :
833 928 : const RankTwoTensor eT = _eigvecs.transpose();
834 :
835 928 : RankFourTensor dstress_dtrial;
836 3712 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
837 11136 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
838 33408 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
839 100224 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
840 300672 : for (unsigned a = 0; a < _num_sp; ++a)
841 225504 : dstress_dtrial(i, j, k, l) +=
842 225504 : drot_dstress(i, a, k, l) * stress_params[a] * eT(a, j) +
843 225504 : _eigvecs(i, a) * stress_params[a] * drot_dstress(j, a, k, l);
844 :
845 928 : const std::vector<RankTwoTensor> dsp_trial = dstress_param_dstress(stress_trial);
846 3712 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
847 11136 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
848 33408 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
849 100224 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
850 300672 : for (unsigned a = 0; a < _num_sp; ++a)
851 902016 : for (unsigned b = 0; b < _num_sp; ++b)
852 676512 : dstress_dtrial(i, j, k, l) +=
853 676512 : _eigvecs(i, a) * dvar_dtrial[a][b] * dsp_trial[b](k, l) * eT(a, j);
854 :
855 928 : cto = dstress_dtrial * elasticity_tensor;
856 928 : }
|
