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