Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
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 "CappedWeakPlaneStressUpdate.h"
11 :
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("TensorMechanicsApp", CappedWeakPlaneStressUpdate);
15 :
16 : InputParameters
17 590 : CappedWeakPlaneStressUpdate::validParams()
18 : {
19 590 : InputParameters params = TwoParameterPlasticityStressUpdate::validParams();
20 590 : params.addClassDescription("Capped weak-plane plasticity stress calculator");
21 1180 : params.addRequiredParam<UserObjectName>(
22 : "cohesion",
23 : "A TensorMechanicsHardening UserObject that defines hardening of the cohesion. "
24 : "Physically the cohesion should not be negative.");
25 1180 : params.addRequiredParam<UserObjectName>("tan_friction_angle",
26 : "A TensorMechanicsHardening UserObject that defines "
27 : "hardening of tan(friction angle). Physically the "
28 : "friction angle should be between 0 and 90deg.");
29 1180 : params.addRequiredParam<UserObjectName>(
30 : "tan_dilation_angle",
31 : "A TensorMechanicsHardening UserObject that defines hardening of the "
32 : "tan(dilation angle). Usually the dilation angle is not greater than "
33 : "the friction angle, and it is between 0 and 90deg.");
34 1180 : params.addRequiredParam<UserObjectName>(
35 : "tensile_strength",
36 : "A TensorMechanicsHardening UserObject that defines hardening of the "
37 : "weak-plane tensile strength. In physical situations this is positive "
38 : "(and always must be greater than negative compressive-strength.");
39 1180 : params.addRequiredParam<UserObjectName>("compressive_strength",
40 : "A TensorMechanicsHardening UserObject that defines "
41 : "hardening of the weak-plane compressive strength. In "
42 : "physical situations this is positive.");
43 1180 : params.addRequiredRangeCheckedParam<Real>(
44 : "tip_smoother",
45 : "tip_smoother>=0",
46 : "The cone vertex at shear-stress = 0 will be smoothed by "
47 : "the given amount. Typical value is 0.1*cohesion");
48 1180 : params.addParam<bool>("perfect_guess",
49 1180 : true,
50 : "Provide a guess to the Newton-Raphson procedure "
51 : "that is the result from perfect plasticity. With "
52 : "severe hardening/softening this may be "
53 : "suboptimal.");
54 590 : return params;
55 0 : }
56 :
57 444 : CappedWeakPlaneStressUpdate::CappedWeakPlaneStressUpdate(const InputParameters & parameters)
58 : : TwoParameterPlasticityStressUpdate(parameters, 3, 2),
59 444 : _cohesion(getUserObject<TensorMechanicsHardeningModel>("cohesion")),
60 444 : _tan_phi(getUserObject<TensorMechanicsHardeningModel>("tan_friction_angle")),
61 444 : _tan_psi(getUserObject<TensorMechanicsHardeningModel>("tan_dilation_angle")),
62 444 : _tstrength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
63 444 : _cstrength(getUserObject<TensorMechanicsHardeningModel>("compressive_strength")),
64 888 : _small_smoother2(Utility::pow<2>(getParam<Real>("tip_smoother"))),
65 888 : _perfect_guess(getParam<bool>("perfect_guess")),
66 444 : _stress_return_type(StressReturnType::nothing_special),
67 444 : _in_trial02(0.0),
68 444 : _in_trial12(0.0),
69 444 : _in_q_trial(0.0)
70 : {
71 : // With arbitary UserObjects, it is impossible to check everything,
72 : // but this will catch the common errors
73 444 : if (_tan_phi.value(0) < 0 || _tan_psi.value(0) < 0)
74 2 : mooseError("CappedWeakPlaneStressUpdate: Weak-plane friction and dilation angles must lie in "
75 : "[0, Pi/2]");
76 442 : if (_tan_phi.value(0) < _tan_psi.value(0))
77 1 : mooseError("CappedWeakPlaneStressUpdate: Weak-plane friction angle must not be less than "
78 : "dilation angle");
79 441 : if (_cohesion.value(0) < 0)
80 1 : mooseError("CappedWeakPlaneStressUpdate: Weak-plane cohesion must not be negative");
81 440 : if (_tstrength.value(0) + _cstrength.value(0) <= _smoothing_tol)
82 1 : mooseError("CappedWeakPlaneStressUpdate: Weak plane tensile strength plus compressive "
83 : "strength must be greater than smoothing_tol");
84 439 : }
85 :
86 : void
87 672874 : CappedWeakPlaneStressUpdate::initializeReturnProcess()
88 : {
89 672874 : _stress_return_type = StressReturnType::nothing_special;
90 672874 : }
91 :
92 : void
93 221056 : CappedWeakPlaneStressUpdate::finalizeReturnProcess(const RankTwoTensor & /*rotation_increment*/)
94 : {
95 221056 : _stress_return_type = StressReturnType::nothing_special;
96 221056 : }
97 :
98 : void
99 220858 : CappedWeakPlaneStressUpdate::preReturnMap(Real /*p_trial*/,
100 : Real q_trial,
101 : const RankTwoTensor & stress_trial,
102 : const std::vector<Real> & /*intnl_old*/,
103 : const std::vector<Real> & yf,
104 : const RankFourTensor & /*Eijkl*/)
105 : {
106 : // If it's obvious, then simplify the return-type
107 220858 : if (yf[1] >= 0)
108 120929 : _stress_return_type = StressReturnType::no_compression;
109 99929 : else if (yf[2] >= 0)
110 36208 : _stress_return_type = StressReturnType::no_tension;
111 :
112 : // The following are useful for the Cosserat case too
113 220858 : _in_trial02 = stress_trial(0, 2);
114 220858 : _in_trial12 = stress_trial(1, 2);
115 220858 : _in_q_trial = q_trial;
116 220858 : }
117 :
118 : void
119 813682 : CappedWeakPlaneStressUpdate::computePQ(const RankTwoTensor & stress, Real & p, Real & q) const
120 : {
121 813682 : p = stress(2, 2);
122 : // Because the following is not explicitly symmeterised, it is useful for the Cosserat case too
123 813682 : q = std::sqrt(Utility::pow<2>(stress(0, 2)) + Utility::pow<2>(stress(1, 2)));
124 813682 : }
125 :
126 : void
127 220858 : CappedWeakPlaneStressUpdate::setEppEqq(const RankFourTensor & Eijkl, Real & Epp, Real & Eqq) const
128 : {
129 220858 : Epp = Eijkl(2, 2, 2, 2);
130 220858 : Eqq = Eijkl(0, 2, 0, 2);
131 220858 : }
132 :
133 : void
134 217568 : CappedWeakPlaneStressUpdate::setStressAfterReturn(const RankTwoTensor & stress_trial,
135 : Real p_ok,
136 : Real q_ok,
137 : Real gaE,
138 : const std::vector<Real> & /*intnl*/,
139 : const yieldAndFlow & smoothed_q,
140 : const RankFourTensor & Eijkl,
141 : RankTwoTensor & stress) const
142 : {
143 217568 : stress = stress_trial;
144 217568 : stress(2, 2) = p_ok;
145 : // stress_xx and stress_yy are sitting at their trial-stress values
146 : // so need to bring them back via Poisson's ratio
147 217568 : stress(0, 0) -= Eijkl(2, 2, 0, 0) * gaE / _Epp * smoothed_q.dg[0];
148 217568 : stress(1, 1) -= Eijkl(2, 2, 1, 1) * gaE / _Epp * smoothed_q.dg[0];
149 217568 : if (_in_q_trial == 0.0)
150 7596 : stress(2, 0) = stress(2, 1) = stress(0, 2) = stress(1, 2) = 0.0;
151 : else
152 : {
153 209972 : stress(2, 0) = stress(0, 2) = _in_trial02 * q_ok / _in_q_trial;
154 209972 : stress(2, 1) = stress(1, 2) = _in_trial12 * q_ok / _in_q_trial;
155 : }
156 217568 : }
157 :
158 : void
159 507883 : CappedWeakPlaneStressUpdate::setIntnlDerivatives(Real /*p_trial*/,
160 : Real q_trial,
161 : Real /*p*/,
162 : Real q,
163 : const std::vector<Real> & intnl,
164 : std::vector<std::vector<Real>> & dintnl) const
165 : {
166 507883 : const Real tanpsi = _tan_psi.value(intnl[0]);
167 507883 : dintnl[0][0] = 0.0;
168 507883 : dintnl[0][1] = -1.0 / _Eqq;
169 507883 : dintnl[1][0] = -1.0 / _Epp;
170 507883 : dintnl[1][1] =
171 507883 : tanpsi / _Eqq - (q_trial - q) * _tan_psi.derivative(intnl[0]) * dintnl[0][1] / _Eqq;
172 507883 : }
173 :
174 : void
175 48048 : CappedWeakPlaneStressUpdate::consistentTangentOperator(const RankTwoTensor & /*stress_trial*/,
176 : Real /*p_trial*/,
177 : Real q_trial,
178 : const RankTwoTensor & /*stress*/,
179 : Real /*p*/,
180 : Real q,
181 : Real gaE,
182 : const yieldAndFlow & smoothed_q,
183 : const RankFourTensor & Eijkl,
184 : bool compute_full_tangent_operator,
185 : RankFourTensor & cto) const
186 : {
187 48048 : cto = Eijkl;
188 48048 : if (!compute_full_tangent_operator)
189 : return;
190 :
191 48048 : const Real Ezzzz = Eijkl(2, 2, 2, 2);
192 48048 : const Real Ezxzx = Eijkl(2, 0, 2, 0);
193 48048 : const Real tanpsi = _tan_psi.value(_intnl[_qp][0]);
194 48048 : const Real dintnl0_dq = -1.0 / Ezxzx;
195 : const Real dintnl0_dqt = 1.0 / Ezxzx;
196 48048 : const Real dintnl1_dp = -1.0 / Ezzzz;
197 : const Real dintnl1_dpt = 1.0 / Ezzzz;
198 : const Real dintnl1_dq =
199 48048 : tanpsi / Ezxzx - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dq / Ezxzx;
200 : const Real dintnl1_dqt =
201 48048 : -tanpsi / Ezxzx - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dqt / Ezxzx;
202 :
203 192192 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
204 : {
205 144144 : const Real dpt_depii = Eijkl(2, 2, i, i);
206 144144 : cto(2, 2, i, i) = _dp_dpt * dpt_depii;
207 : const Real poisson_effect =
208 144144 : Eijkl(2, 2, 0, 0) / Ezzzz *
209 144144 : (_dgaE_dpt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dpt +
210 144144 : gaE * smoothed_q.d2g[0][1] * _dq_dpt +
211 144144 : gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dq * _dq_dpt) +
212 144144 : gaE * smoothed_q.d2g_di[0][1] *
213 144144 : (dintnl1_dpt + dintnl1_dp * _dp_dpt + dintnl1_dq * _dq_dpt)) *
214 144144 : dpt_depii;
215 144144 : cto(0, 0, i, i) -= poisson_effect;
216 144144 : cto(1, 1, i, i) -= poisson_effect;
217 144144 : if (q_trial > 0.0)
218 : {
219 140700 : cto(2, 0, i, i) = cto(0, 2, i, i) = _in_trial02 / q_trial * _dq_dpt * dpt_depii;
220 140700 : cto(2, 1, i, i) = cto(1, 2, i, i) = _in_trial12 / q_trial * _dq_dpt * dpt_depii;
221 : }
222 : }
223 :
224 : const Real poisson_effect =
225 48048 : -Eijkl(2, 2, 0, 0) / Ezzzz *
226 48048 : (_dgaE_dqt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dqt +
227 48048 : gaE * smoothed_q.d2g[0][1] * _dq_dqt +
228 48048 : gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dqt + dintnl0_dq * _dq_dqt) +
229 48048 : gaE * smoothed_q.d2g_di[0][1] * (dintnl1_dqt + dintnl1_dp * _dp_dqt + dintnl1_dq * _dq_dqt));
230 :
231 48048 : const Real dqt_dep20 = (q_trial == 0.0 ? 1.0 : _in_trial02 / q_trial) * Eijkl(2, 0, 2, 0);
232 48048 : cto(2, 2, 2, 0) = cto(2, 2, 0, 2) = _dp_dqt * dqt_dep20;
233 48048 : cto(0, 0, 2, 0) = cto(0, 0, 0, 2) = cto(1, 1, 2, 0) = cto(1, 1, 0, 2) =
234 48048 : poisson_effect * dqt_dep20;
235 48048 : if (q_trial > 0.0)
236 : {
237 : // for q_trial=0, Jacobian_mult is just given by the elastic case
238 46900 : cto(0, 2, 0, 2) = cto(2, 0, 0, 2) = cto(0, 2, 2, 0) = cto(2, 0, 2, 0) =
239 46900 : Eijkl(2, 0, 2, 0) * q / q_trial +
240 46900 : _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
241 46900 : cto(1, 2, 0, 2) = cto(2, 1, 0, 2) = cto(1, 2, 2, 0) = cto(2, 1, 2, 0) =
242 46900 : _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
243 : }
244 :
245 48048 : const Real dqt_dep21 = (q_trial == 0.0 ? 1.0 : _in_trial12 / q_trial) * Eijkl(2, 1, 2, 1);
246 48048 : cto(2, 2, 2, 1) = cto(2, 2, 1, 2) = _dp_dqt * dqt_dep21;
247 48048 : cto(0, 0, 2, 1) = cto(0, 0, 1, 2) = cto(1, 1, 2, 1) = cto(1, 1, 1, 2) =
248 48048 : poisson_effect * dqt_dep21;
249 48048 : if (q_trial > 0.0)
250 : {
251 : // for q_trial=0, Jacobian_mult is just given by the elastic case
252 46900 : cto(0, 2, 1, 2) = cto(2, 0, 1, 2) = cto(0, 2, 2, 1) = cto(2, 0, 2, 1) =
253 46900 : _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
254 46900 : cto(1, 2, 1, 2) = cto(2, 1, 1, 2) = cto(1, 2, 2, 1) = cto(2, 1, 2, 1) =
255 46900 : Eijkl(2, 1, 2, 1) * q / q_trial +
256 46900 : _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
257 : }
258 : }
259 :
260 : void
261 893930 : CappedWeakPlaneStressUpdate::yieldFunctionValues(Real p,
262 : Real q,
263 : const std::vector<Real> & intnl,
264 : std::vector<Real> & yf) const
265 : {
266 1787860 : yf[0] = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * _tan_phi.value(intnl[0]) -
267 893930 : _cohesion.value(intnl[0]);
268 :
269 893930 : if (_stress_return_type == StressReturnType::no_tension)
270 0 : yf[1] = std::numeric_limits<Real>::lowest();
271 : else
272 893930 : yf[1] = p - _tstrength.value(intnl[1]);
273 :
274 893930 : if (_stress_return_type == StressReturnType::no_compression)
275 0 : yf[2] = std::numeric_limits<Real>::lowest();
276 : else
277 893930 : yf[2] = -p - _cstrength.value(intnl[1]);
278 893930 : }
279 :
280 : void
281 691567 : CappedWeakPlaneStressUpdate::computeAllQ(Real p,
282 : Real q,
283 : const std::vector<Real> & intnl,
284 : std::vector<yieldAndFlow> & all_q) const
285 : {
286 : // yield function values
287 1383134 : all_q[0].f = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * _tan_phi.value(intnl[0]) -
288 691567 : _cohesion.value(intnl[0]);
289 691567 : if (_stress_return_type == StressReturnType::no_tension)
290 143120 : all_q[1].f = std::numeric_limits<Real>::lowest();
291 : else
292 548447 : all_q[1].f = p - _tstrength.value(intnl[1]);
293 691567 : if (_stress_return_type == StressReturnType::no_compression)
294 317586 : all_q[2].f = std::numeric_limits<Real>::lowest();
295 : else
296 373981 : all_q[2].f = -p - _cstrength.value(intnl[1]);
297 :
298 : // d(yield Function)/d(p, q)
299 : // derivatives wrt p
300 691567 : all_q[0].df[0] = _tan_phi.value(intnl[0]);
301 691567 : all_q[1].df[0] = 1.0;
302 691567 : all_q[2].df[0] = -1.0;
303 :
304 : // derivatives wrt q
305 691567 : if (_small_smoother2 == 0.0)
306 90997 : all_q[0].df[1] = 1.0;
307 : else
308 600570 : all_q[0].df[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
309 691567 : all_q[1].df[1] = 0.0;
310 691567 : all_q[2].df[1] = 0.0;
311 :
312 : // d(yield Function)/d(intnl)
313 : // derivatives wrt intnl[0] (shear plastic strain)
314 691567 : all_q[0].df_di[0] = p * _tan_phi.derivative(intnl[0]) - _cohesion.derivative(intnl[0]);
315 691567 : all_q[1].df_di[0] = 0.0;
316 691567 : all_q[2].df_di[0] = 0.0;
317 : // derivatives wrt intnl[q] (tensile plastic strain)
318 691567 : all_q[0].df_di[1] = 0.0;
319 691567 : all_q[1].df_di[1] = -_tstrength.derivative(intnl[1]);
320 691567 : all_q[2].df_di[1] = -_cstrength.derivative(intnl[1]);
321 :
322 : // d(flowPotential)/d(p, q)
323 : // derivatives wrt p
324 691567 : all_q[0].dg[0] = _tan_psi.value(intnl[0]);
325 691567 : all_q[1].dg[0] = 1.0;
326 691567 : all_q[2].dg[0] = -1.0;
327 : // derivatives wrt q
328 691567 : if (_small_smoother2 == 0.0)
329 90997 : all_q[0].dg[1] = 1.0;
330 : else
331 600570 : all_q[0].dg[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
332 691567 : all_q[1].dg[1] = 0.0;
333 691567 : all_q[2].dg[1] = 0.0;
334 :
335 : // d2(flowPotential)/d(p, q)/d(intnl)
336 : // d(dg/dp)/dintnl[0]
337 691567 : all_q[0].d2g_di[0][0] = _tan_psi.derivative(intnl[0]);
338 691567 : all_q[1].d2g_di[0][0] = 0.0;
339 691567 : all_q[2].d2g_di[0][0] = 0.0;
340 : // d(dg/dp)/dintnl[1]
341 691567 : all_q[0].d2g_di[0][1] = 0.0;
342 691567 : all_q[1].d2g_di[0][1] = 0.0;
343 691567 : all_q[2].d2g_di[0][1] = 0.0;
344 : // d(dg/dq)/dintnl[0]
345 691567 : all_q[0].d2g_di[1][0] = 0.0;
346 691567 : all_q[1].d2g_di[1][0] = 0.0;
347 691567 : all_q[2].d2g_di[1][0] = 0.0;
348 : // d(dg/dq)/dintnl[1]
349 691567 : all_q[0].d2g_di[1][1] = 0.0;
350 691567 : all_q[1].d2g_di[1][1] = 0.0;
351 691567 : all_q[2].d2g_di[1][1] = 0.0;
352 :
353 : // d2(flowPotential)/d(p, q)/d(p, q)
354 : // d(dg/dp)/dp
355 691567 : all_q[0].d2g[0][0] = 0.0;
356 691567 : all_q[1].d2g[0][0] = 0.0;
357 691567 : all_q[2].d2g[0][0] = 0.0;
358 : // d(dg/dp)/dq
359 691567 : all_q[0].d2g[0][1] = 0.0;
360 691567 : all_q[1].d2g[0][1] = 0.0;
361 691567 : all_q[2].d2g[0][1] = 0.0;
362 : // d(dg/dq)/dp
363 691567 : all_q[0].d2g[1][0] = 0.0;
364 691567 : all_q[1].d2g[1][0] = 0.0;
365 691567 : all_q[2].d2g[1][0] = 0.0;
366 : // d(dg/dq)/dq
367 691567 : if (_small_smoother2 == 0.0)
368 90997 : all_q[0].d2g[1][1] = 0.0;
369 : else
370 600570 : all_q[0].d2g[1][1] = _small_smoother2 / std::pow(Utility::pow<2>(q) + _small_smoother2, 1.5);
371 691567 : all_q[1].d2g[1][1] = 0.0;
372 691567 : all_q[2].d2g[1][1] = 0.0;
373 691567 : }
374 :
375 : void
376 221098 : CappedWeakPlaneStressUpdate::initializeVars(Real p_trial,
377 : Real q_trial,
378 : const std::vector<Real> & intnl_old,
379 : Real & p,
380 : Real & q,
381 : Real & gaE,
382 : std::vector<Real> & intnl) const
383 : {
384 221098 : const Real tanpsi = _tan_psi.value(intnl_old[0]);
385 :
386 221098 : if (!_perfect_guess)
387 : {
388 17024 : p = p_trial;
389 17024 : q = q_trial;
390 17024 : gaE = 0.0;
391 : }
392 : else
393 : {
394 204074 : const Real coh = _cohesion.value(intnl_old[0]);
395 204074 : const Real tanphi = _tan_phi.value(intnl_old[0]);
396 204074 : const Real tens = _tstrength.value(intnl_old[1]);
397 204074 : const Real comp = _cstrength.value(intnl_old[1]);
398 204074 : const Real q_at_T = coh - tens * tanphi;
399 204074 : const Real q_at_C = coh + comp * tanphi;
400 :
401 204074 : if ((p_trial >= tens) && (q_trial <= q_at_T))
402 : {
403 : // pure tensile failure
404 81065 : p = tens;
405 81065 : q = q_trial;
406 81065 : gaE = p_trial - p;
407 : }
408 123009 : else if ((p_trial <= -comp) && (q_trial <= q_at_C))
409 : {
410 : // pure compressive failure
411 15336 : p = -comp;
412 15336 : q = q_trial;
413 15336 : gaE = p - p_trial;
414 : }
415 : else
416 : {
417 : // shear failure or a mixture
418 : // Calculate ga assuming a pure shear failure
419 107673 : const Real ga = (q_trial + p_trial * tanphi - coh) / (_Eqq + _Epp * tanphi * tanpsi);
420 107673 : if (ga <= 0 && p_trial <= tens && p_trial >= -comp)
421 : {
422 : // very close to one of the rounded corners: there is no advantage to guessing the
423 : // solution, so:
424 11016 : p = p_trial;
425 11016 : q = q_trial;
426 11016 : gaE = 0.0;
427 : }
428 : else
429 : {
430 96657 : q = q_trial - _Eqq * ga;
431 96657 : if (q <= 0.0 && q_at_T <= 0.0)
432 : {
433 : // user has set tensile strength so large that it is irrelevant: return to the tip of the
434 : // shear cone
435 2985 : q = 0.0;
436 2985 : p = coh / tanphi;
437 2985 : gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
438 : }
439 93672 : else if (q <= q_at_T)
440 : {
441 : // pure shear is incorrect: mixture of tension and shear is correct
442 25800 : q = q_at_T;
443 25800 : p = tens;
444 25800 : gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
445 : }
446 67872 : else if (q >= q_at_C)
447 : {
448 : // pure shear is incorrect: mixture of compression and shear is correct
449 13496 : q = q_at_C;
450 13496 : p = -comp;
451 13496 : if (p - p_trial < _Epp * tanpsi * (q_trial - q) / _Eqq)
452 : // trial point is sitting almost directly above corner
453 8080 : gaE = (q_trial - q) * _Epp / _Eqq;
454 : else
455 : // trial point is sitting to the left of the corner
456 5416 : gaE = (p - p_trial) / tanpsi;
457 : }
458 : else
459 : {
460 : // pure shear was correct
461 54376 : p = p_trial - _Epp * ga * tanpsi;
462 54376 : gaE = ga * _Epp;
463 : }
464 : }
465 : }
466 : }
467 221098 : setIntnlValues(p_trial, q_trial, p, q, intnl_old, intnl);
468 221098 : }
469 :
470 : void
471 912623 : CappedWeakPlaneStressUpdate::setIntnlValues(Real p_trial,
472 : Real q_trial,
473 : Real p,
474 : Real q,
475 : const std::vector<Real> & intnl_old,
476 : std::vector<Real> & intnl) const
477 : {
478 912623 : intnl[0] = intnl_old[0] + (q_trial - q) / _Eqq;
479 912623 : const Real tanpsi = _tan_psi.value(intnl[0]);
480 912623 : intnl[1] = intnl_old[1] + (p_trial - p) / _Epp - (q_trial - q) * tanpsi / _Eqq;
481 912623 : }
482 :
483 : RankTwoTensor
484 220816 : CappedWeakPlaneStressUpdate::dpdstress(const RankTwoTensor & /*stress*/) const
485 : {
486 220816 : RankTwoTensor deriv = RankTwoTensor();
487 220816 : deriv(2, 2) = 1.0;
488 220816 : return deriv;
489 : }
490 :
491 : RankFourTensor
492 0 : CappedWeakPlaneStressUpdate::d2pdstress2(const RankTwoTensor & /*stress*/) const
493 : {
494 0 : return RankFourTensor();
495 : }
496 :
497 : RankTwoTensor
498 217808 : CappedWeakPlaneStressUpdate::dqdstress(const RankTwoTensor & stress) const
499 : {
500 217808 : RankTwoTensor deriv = RankTwoTensor();
501 217808 : const Real q = std::sqrt(Utility::pow<2>(stress(2, 0)) + Utility::pow<2>(stress(2, 1)));
502 217808 : if (q > 0.0)
503 : {
504 210212 : deriv(2, 0) = deriv(0, 2) = 0.5 * stress(2, 0) / q;
505 210212 : deriv(2, 1) = deriv(1, 2) = 0.5 * stress(2, 1) / q;
506 : }
507 : else
508 : {
509 : // derivative is not defined here. For now i'll set:
510 7596 : deriv(2, 0) = deriv(0, 2) = 0.5;
511 7596 : deriv(2, 1) = deriv(1, 2) = 0.5;
512 : }
513 217808 : return deriv;
514 : }
515 :
516 : RankFourTensor
517 0 : CappedWeakPlaneStressUpdate::d2qdstress2(const RankTwoTensor & stress) const
518 : {
519 0 : RankFourTensor d2 = RankFourTensor();
520 :
521 0 : const Real q = std::sqrt(Utility::pow<2>(stress(2, 0)) + Utility::pow<2>(stress(2, 1)));
522 0 : if (q == 0.0)
523 : return d2;
524 :
525 0 : RankTwoTensor dq = RankTwoTensor();
526 0 : dq(2, 0) = dq(0, 2) = 0.25 * (stress(2, 0) + stress(0, 2)) / q;
527 0 : dq(2, 1) = dq(1, 2) = 0.25 * (stress(2, 1) + stress(1, 2)) / q;
528 :
529 0 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
530 0 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
531 0 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
532 0 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
533 0 : d2(i, j, k, l) = -dq(i, j) * dq(k, l) / q;
534 :
535 0 : d2(0, 2, 0, 2) += 0.25 / q;
536 0 : d2(0, 2, 2, 0) += 0.25 / q;
537 0 : d2(2, 0, 0, 2) += 0.25 / q;
538 0 : d2(2, 0, 2, 0) += 0.25 / q;
539 0 : d2(1, 2, 1, 2) += 0.25 / q;
540 0 : d2(1, 2, 2, 1) += 0.25 / q;
541 0 : d2(2, 1, 1, 2) += 0.25 / q;
542 0 : d2(2, 1, 2, 1) += 0.25 / q;
543 :
544 : return d2;
545 : }
|
