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 "CappedWeakPlaneStressUpdate.h"
11 :
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", CappedWeakPlaneStressUpdate);
15 :
16 : InputParameters
17 1280 : CappedWeakPlaneStressUpdate::validParams()
18 : {
19 1280 : InputParameters params = TwoParameterPlasticityStressUpdate::validParams();
20 1280 : params.addClassDescription("Capped weak-plane plasticity stress calculator");
21 2560 : params.addRequiredParam<UserObjectName>(
22 : "cohesion",
23 : "A SolidMechanicsHardening UserObject that defines hardening of the cohesion. "
24 : "Physically the cohesion should not be negative.");
25 2560 : params.addRequiredParam<UserObjectName>("tan_friction_angle",
26 : "A SolidMechanicsHardening UserObject that defines "
27 : "hardening of tan(friction angle). Physically the "
28 : "friction angle should be between 0 and 90deg.");
29 2560 : params.addRequiredParam<UserObjectName>(
30 : "tan_dilation_angle",
31 : "A SolidMechanicsHardening 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 2560 : params.addRequiredParam<UserObjectName>(
35 : "tensile_strength",
36 : "A SolidMechanicsHardening 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 2560 : params.addRequiredParam<UserObjectName>("compressive_strength",
40 : "A SolidMechanicsHardening UserObject that defines "
41 : "hardening of the weak-plane compressive strength. In "
42 : "physical situations this is positive.");
43 2560 : 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 2560 : params.addParam<bool>("perfect_guess",
49 2560 : 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 1280 : return params;
55 0 : }
56 :
57 963 : CappedWeakPlaneStressUpdate::CappedWeakPlaneStressUpdate(const InputParameters & parameters)
58 : : TwoParameterPlasticityStressUpdate(parameters, 3, 2),
59 963 : _cohesion(getUserObject<SolidMechanicsHardeningModel>("cohesion")),
60 963 : _tan_phi(getUserObject<SolidMechanicsHardeningModel>("tan_friction_angle")),
61 963 : _tan_psi(getUserObject<SolidMechanicsHardeningModel>("tan_dilation_angle")),
62 963 : _tstrength(getUserObject<SolidMechanicsHardeningModel>("tensile_strength")),
63 963 : _cstrength(getUserObject<SolidMechanicsHardeningModel>("compressive_strength")),
64 1926 : _small_smoother2(Utility::pow<2>(getParam<Real>("tip_smoother"))),
65 1926 : _perfect_guess(getParam<bool>("perfect_guess")),
66 963 : _stress_return_type(StressReturnType::nothing_special),
67 963 : _in_trial02(0.0),
68 963 : _in_trial12(0.0),
69 963 : _in_q_trial(0.0)
70 : {
71 : // With arbitrary UserObjects, it is impossible to check everything,
72 : // but this will catch the common errors
73 963 : if (_tan_phi.value(0) < 0 || _tan_psi.value(0) < 0)
74 4 : mooseError("CappedWeakPlaneStressUpdate: Weak-plane friction and dilation angles must lie in "
75 : "[0, Pi/2]");
76 959 : if (_tan_phi.value(0) < _tan_psi.value(0))
77 2 : mooseError("CappedWeakPlaneStressUpdate: Weak-plane friction angle must not be less than "
78 : "dilation angle");
79 957 : if (_cohesion.value(0) < 0)
80 2 : mooseError("CappedWeakPlaneStressUpdate: Weak-plane cohesion must not be negative");
81 955 : if (_tstrength.value(0) + _cstrength.value(0) <= _smoothing_tol)
82 2 : mooseError("CappedWeakPlaneStressUpdate: Weak plane tensile strength plus compressive "
83 : "strength must be greater than smoothing_tol");
84 953 : }
85 :
86 : void
87 1550876 : CappedWeakPlaneStressUpdate::initializeReturnProcess()
88 : {
89 1550876 : _stress_return_type = StressReturnType::nothing_special;
90 1550876 : }
91 :
92 : void
93 544084 : CappedWeakPlaneStressUpdate::finalizeReturnProcess(const RankTwoTensor & /*rotation_increment*/)
94 : {
95 544084 : _stress_return_type = StressReturnType::nothing_special;
96 544084 : }
97 :
98 : void
99 543496 : 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 543496 : if (yf[1] >= 0)
108 297646 : _stress_return_type = StressReturnType::no_compression;
109 245850 : else if (yf[2] >= 0)
110 90232 : _stress_return_type = StressReturnType::no_tension;
111 :
112 : // The following are useful for the Cosserat case too
113 543496 : _in_trial02 = stress_trial(0, 2);
114 543496 : _in_trial12 = stress_trial(1, 2);
115 543496 : _in_q_trial = q_trial;
116 543496 : }
117 :
118 : void
119 1909696 : CappedWeakPlaneStressUpdate::computePQ(const RankTwoTensor & stress, Real & p, Real & q) const
120 : {
121 1909696 : p = stress(2, 2);
122 : // Because the following is not explicitly symmeterised, it is useful for the Cosserat case too
123 1909696 : q = std::sqrt(Utility::pow<2>(stress(0, 2)) + Utility::pow<2>(stress(1, 2)));
124 1909696 : }
125 :
126 : void
127 543496 : CappedWeakPlaneStressUpdate::setEppEqq(const RankFourTensor & Eijkl, Real & Epp, Real & Eqq) const
128 : {
129 543496 : Epp = Eijkl(2, 2, 2, 2);
130 543496 : Eqq = Eijkl(0, 2, 0, 2);
131 543496 : }
132 :
133 : void
134 536756 : 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 536756 : stress = stress_trial;
144 536756 : 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 536756 : stress(0, 0) -= Eijkl(2, 2, 0, 0) * gaE / _Epp * smoothed_q.dg[0];
148 536756 : stress(1, 1) -= Eijkl(2, 2, 1, 1) * gaE / _Epp * smoothed_q.dg[0];
149 536756 : if (_in_q_trial == 0.0)
150 13265 : stress(2, 0) = stress(2, 1) = stress(0, 2) = stress(1, 2) = 0.0;
151 : else
152 : {
153 523491 : stress(2, 0) = stress(0, 2) = _in_trial02 * q_ok / _in_q_trial;
154 523491 : stress(2, 1) = stress(1, 2) = _in_trial12 * q_ok / _in_q_trial;
155 : }
156 536756 : }
157 :
158 : void
159 1260606 : 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 1260606 : const Real tanpsi = _tan_psi.value(intnl[0]);
167 1260606 : dintnl[0][0] = 0.0;
168 1260606 : dintnl[0][1] = -1.0 / _Eqq;
169 1260606 : dintnl[1][0] = -1.0 / _Epp;
170 1260606 : dintnl[1][1] =
171 1260606 : tanpsi / _Eqq - (q_trial - q) * _tan_psi.derivative(intnl[0]) * dintnl[0][1] / _Eqq;
172 1260606 : }
173 :
174 : void
175 121448 : 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 121448 : cto = Eijkl;
188 121448 : if (!compute_full_tangent_operator)
189 : return;
190 :
191 121448 : const Real Ezzzz = Eijkl(2, 2, 2, 2);
192 121448 : const Real Ezxzx = Eijkl(2, 0, 2, 0);
193 121448 : const Real tanpsi = _tan_psi.value(_intnl[_qp][0]);
194 121448 : const Real dintnl0_dq = -1.0 / Ezxzx;
195 : const Real dintnl0_dqt = 1.0 / Ezxzx;
196 121448 : const Real dintnl1_dp = -1.0 / Ezzzz;
197 : const Real dintnl1_dpt = 1.0 / Ezzzz;
198 : const Real dintnl1_dq =
199 121448 : tanpsi / Ezxzx - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dq / Ezxzx;
200 : const Real dintnl1_dqt =
201 121448 : -tanpsi / Ezxzx - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dqt / Ezxzx;
202 :
203 485792 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
204 : {
205 364344 : const Real dpt_depii = Eijkl(2, 2, i, i);
206 364344 : cto(2, 2, i, i) = _dp_dpt * dpt_depii;
207 : const Real poisson_effect =
208 364344 : Eijkl(2, 2, 0, 0) / Ezzzz *
209 364344 : (_dgaE_dpt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dpt +
210 364344 : gaE * smoothed_q.d2g[0][1] * _dq_dpt +
211 364344 : gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dq * _dq_dpt) +
212 364344 : gaE * smoothed_q.d2g_di[0][1] *
213 364344 : (dintnl1_dpt + dintnl1_dp * _dp_dpt + dintnl1_dq * _dq_dpt)) *
214 364344 : dpt_depii;
215 364344 : cto(0, 0, i, i) -= poisson_effect;
216 364344 : cto(1, 1, i, i) -= poisson_effect;
217 364344 : if (q_trial > 0.0)
218 : {
219 355785 : cto(2, 0, i, i) = cto(0, 2, i, i) = _in_trial02 / q_trial * _dq_dpt * dpt_depii;
220 355785 : 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 121448 : -Eijkl(2, 2, 0, 0) / Ezzzz *
226 121448 : (_dgaE_dqt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dqt +
227 121448 : gaE * smoothed_q.d2g[0][1] * _dq_dqt +
228 121448 : gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dqt + dintnl0_dq * _dq_dqt) +
229 121448 : gaE * smoothed_q.d2g_di[0][1] * (dintnl1_dqt + dintnl1_dp * _dp_dqt + dintnl1_dq * _dq_dqt));
230 :
231 121448 : const Real dqt_dep20 = (q_trial == 0.0 ? 1.0 : _in_trial02 / q_trial) * Eijkl(2, 0, 2, 0);
232 121448 : cto(2, 2, 2, 0) = cto(2, 2, 0, 2) = _dp_dqt * dqt_dep20;
233 121448 : cto(0, 0, 2, 0) = cto(0, 0, 0, 2) = cto(1, 1, 2, 0) = cto(1, 1, 0, 2) =
234 121448 : poisson_effect * dqt_dep20;
235 121448 : if (q_trial > 0.0)
236 : {
237 : // for q_trial=0, Jacobian_mult is just given by the elastic case
238 118595 : cto(0, 2, 0, 2) = cto(2, 0, 0, 2) = cto(0, 2, 2, 0) = cto(2, 0, 2, 0) =
239 118595 : Eijkl(2, 0, 2, 0) * q / q_trial +
240 118595 : _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
241 118595 : cto(1, 2, 0, 2) = cto(2, 1, 0, 2) = cto(1, 2, 2, 0) = cto(2, 1, 2, 0) =
242 118595 : _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
243 : }
244 :
245 121448 : const Real dqt_dep21 = (q_trial == 0.0 ? 1.0 : _in_trial12 / q_trial) * Eijkl(2, 1, 2, 1);
246 121448 : cto(2, 2, 2, 1) = cto(2, 2, 1, 2) = _dp_dqt * dqt_dep21;
247 121448 : cto(0, 0, 2, 1) = cto(0, 0, 1, 2) = cto(1, 1, 2, 1) = cto(1, 1, 1, 2) =
248 121448 : poisson_effect * dqt_dep21;
249 121448 : if (q_trial > 0.0)
250 : {
251 : // for q_trial=0, Jacobian_mult is just given by the elastic case
252 118595 : cto(0, 2, 1, 2) = cto(2, 0, 1, 2) = cto(0, 2, 2, 1) = cto(2, 0, 2, 1) =
253 118595 : _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
254 118595 : cto(1, 2, 1, 2) = cto(2, 1, 1, 2) = cto(1, 2, 2, 1) = cto(2, 1, 2, 1) =
255 118595 : Eijkl(2, 1, 2, 1) * q / q_trial +
256 118595 : _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
257 : }
258 : }
259 :
260 : void
261 2094960 : CappedWeakPlaneStressUpdate::yieldFunctionValues(Real p,
262 : Real q,
263 : const std::vector<Real> & intnl,
264 : std::vector<Real> & yf) const
265 : {
266 4189920 : yf[0] = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * _tan_phi.value(intnl[0]) -
267 2094960 : _cohesion.value(intnl[0]);
268 :
269 2094960 : if (_stress_return_type == StressReturnType::no_tension)
270 0 : yf[1] = std::numeric_limits<Real>::lowest();
271 : else
272 2094960 : yf[1] = p - _tstrength.value(intnl[1]);
273 :
274 2094960 : if (_stress_return_type == StressReturnType::no_compression)
275 0 : yf[2] = std::numeric_limits<Real>::lowest();
276 : else
277 2094960 : yf[2] = -p - _cstrength.value(intnl[1]);
278 2094960 : }
279 :
280 : void
281 1709170 : 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 3418340 : all_q[0].f = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * _tan_phi.value(intnl[0]) -
288 1709170 : _cohesion.value(intnl[0]);
289 1709170 : if (_stress_return_type == StressReturnType::no_tension)
290 352792 : all_q[1].f = std::numeric_limits<Real>::lowest();
291 : else
292 1356378 : all_q[1].f = p - _tstrength.value(intnl[1]);
293 1709170 : if (_stress_return_type == StressReturnType::no_compression)
294 779616 : all_q[2].f = std::numeric_limits<Real>::lowest();
295 : else
296 929554 : all_q[2].f = -p - _cstrength.value(intnl[1]);
297 :
298 : // d(yield Function)/d(p, q)
299 : // derivatives wrt p
300 1709170 : all_q[0].df[0] = _tan_phi.value(intnl[0]);
301 1709170 : all_q[1].df[0] = 1.0;
302 1709170 : all_q[2].df[0] = -1.0;
303 :
304 : // derivatives wrt q
305 1709170 : if (_small_smoother2 == 0.0)
306 211594 : all_q[0].df[1] = 1.0;
307 : else
308 1497576 : all_q[0].df[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
309 1709170 : all_q[1].df[1] = 0.0;
310 1709170 : all_q[2].df[1] = 0.0;
311 :
312 : // d(yield Function)/d(intnl)
313 : // derivatives wrt intnl[0] (shear plastic strain)
314 1709170 : all_q[0].df_di[0] = p * _tan_phi.derivative(intnl[0]) - _cohesion.derivative(intnl[0]);
315 1709170 : all_q[1].df_di[0] = 0.0;
316 1709170 : all_q[2].df_di[0] = 0.0;
317 : // derivatives wrt intnl[q] (tensile plastic strain)
318 1709170 : all_q[0].df_di[1] = 0.0;
319 1709170 : all_q[1].df_di[1] = -_tstrength.derivative(intnl[1]);
320 1709170 : all_q[2].df_di[1] = -_cstrength.derivative(intnl[1]);
321 :
322 : // d(flowPotential)/d(p, q)
323 : // derivatives wrt p
324 1709170 : all_q[0].dg[0] = _tan_psi.value(intnl[0]);
325 1709170 : all_q[1].dg[0] = 1.0;
326 1709170 : all_q[2].dg[0] = -1.0;
327 : // derivatives wrt q
328 1709170 : if (_small_smoother2 == 0.0)
329 211594 : all_q[0].dg[1] = 1.0;
330 : else
331 1497576 : all_q[0].dg[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
332 1709170 : all_q[1].dg[1] = 0.0;
333 1709170 : all_q[2].dg[1] = 0.0;
334 :
335 : // d2(flowPotential)/d(p, q)/d(intnl)
336 : // d(dg/dp)/dintnl[0]
337 1709170 : all_q[0].d2g_di[0][0] = _tan_psi.derivative(intnl[0]);
338 1709170 : all_q[1].d2g_di[0][0] = 0.0;
339 1709170 : all_q[2].d2g_di[0][0] = 0.0;
340 : // d(dg/dp)/dintnl[1]
341 1709170 : all_q[0].d2g_di[0][1] = 0.0;
342 1709170 : all_q[1].d2g_di[0][1] = 0.0;
343 1709170 : all_q[2].d2g_di[0][1] = 0.0;
344 : // d(dg/dq)/dintnl[0]
345 1709170 : all_q[0].d2g_di[1][0] = 0.0;
346 1709170 : all_q[1].d2g_di[1][0] = 0.0;
347 1709170 : all_q[2].d2g_di[1][0] = 0.0;
348 : // d(dg/dq)/dintnl[1]
349 1709170 : all_q[0].d2g_di[1][1] = 0.0;
350 1709170 : all_q[1].d2g_di[1][1] = 0.0;
351 1709170 : 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 1709170 : all_q[0].d2g[0][0] = 0.0;
356 1709170 : all_q[1].d2g[0][0] = 0.0;
357 1709170 : all_q[2].d2g[0][0] = 0.0;
358 : // d(dg/dp)/dq
359 1709170 : all_q[0].d2g[0][1] = 0.0;
360 1709170 : all_q[1].d2g[0][1] = 0.0;
361 1709170 : all_q[2].d2g[0][1] = 0.0;
362 : // d(dg/dq)/dp
363 1709170 : all_q[0].d2g[1][0] = 0.0;
364 1709170 : all_q[1].d2g[1][0] = 0.0;
365 1709170 : all_q[2].d2g[1][0] = 0.0;
366 : // d(dg/dq)/dq
367 1709170 : if (_small_smoother2 == 0.0)
368 211594 : all_q[0].d2g[1][1] = 0.0;
369 : else
370 1497576 : all_q[0].d2g[1][1] = _small_smoother2 / std::pow(Utility::pow<2>(q) + _small_smoother2, 1.5);
371 1709170 : all_q[1].d2g[1][1] = 0.0;
372 1709170 : all_q[2].d2g[1][1] = 0.0;
373 1709170 : }
374 :
375 : void
376 544168 : 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 544168 : const Real tanpsi = _tan_psi.value(intnl_old[0]);
385 :
386 544168 : if (!_perfect_guess)
387 : {
388 42200 : p = p_trial;
389 42200 : q = q_trial;
390 42200 : gaE = 0.0;
391 : }
392 : else
393 : {
394 501968 : const Real coh = _cohesion.value(intnl_old[0]);
395 501968 : const Real tanphi = _tan_phi.value(intnl_old[0]);
396 501968 : const Real tens = _tstrength.value(intnl_old[1]);
397 501968 : const Real comp = _cstrength.value(intnl_old[1]);
398 501968 : const Real q_at_T = coh - tens * tanphi;
399 501968 : const Real q_at_C = coh + comp * tanphi;
400 :
401 501968 : if ((p_trial >= tens) && (q_trial <= q_at_T))
402 : {
403 : // pure tensile failure
404 197362 : p = tens;
405 197362 : q = q_trial;
406 197362 : gaE = p_trial - p;
407 : }
408 304606 : else if ((p_trial <= -comp) && (q_trial <= q_at_C))
409 : {
410 : // pure compressive failure
411 38820 : p = -comp;
412 38820 : q = q_trial;
413 38820 : gaE = p - p_trial;
414 : }
415 : else
416 : {
417 : // shear failure or a mixture
418 : // Calculate ga assuming a pure shear failure
419 265786 : const Real ga = (q_trial + p_trial * tanphi - coh) / (_Eqq + _Epp * tanphi * tanpsi);
420 265786 : 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 27844 : p = p_trial;
425 27844 : q = q_trial;
426 27844 : gaE = 0.0;
427 : }
428 : else
429 : {
430 237942 : q = q_trial - _Eqq * ga;
431 237942 : 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 7442 : q = 0.0;
436 7442 : p = coh / tanphi;
437 7442 : gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
438 : }
439 230500 : else if (q <= q_at_T)
440 : {
441 : // pure shear is incorrect: mixture of tension and shear is correct
442 64952 : q = q_at_T;
443 64952 : p = tens;
444 64952 : gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
445 : }
446 165548 : else if (q >= q_at_C)
447 : {
448 : // pure shear is incorrect: mixture of compression and shear is correct
449 33316 : q = q_at_C;
450 33316 : p = -comp;
451 33316 : if (p - p_trial < _Epp * tanpsi * (q_trial - q) / _Eqq)
452 : // trial point is sitting almost directly above corner
453 19960 : gaE = (q_trial - q) * _Epp / _Eqq;
454 : else
455 : // trial point is sitting to the left of the corner
456 13356 : gaE = (p - p_trial) / tanpsi;
457 : }
458 : else
459 : {
460 : // pure shear was correct
461 132232 : p = p_trial - _Epp * ga * tanpsi;
462 132232 : gaE = ga * _Epp;
463 : }
464 : }
465 : }
466 : }
467 544168 : setIntnlValues(p_trial, q_trial, p, q, intnl_old, intnl);
468 544168 : }
469 :
470 : void
471 2253254 : 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 2253254 : intnl[0] = intnl_old[0] + (q_trial - q) / _Eqq;
479 2253254 : const Real tanpsi = _tan_psi.value(intnl[0]);
480 2253254 : intnl[1] = intnl_old[1] + (p_trial - p) / _Epp - (q_trial - q) * tanpsi / _Eqq;
481 2253254 : }
482 :
483 : RankTwoTensor
484 543412 : CappedWeakPlaneStressUpdate::dpdstress(const RankTwoTensor & /*stress*/) const
485 : {
486 543412 : RankTwoTensor deriv = RankTwoTensor();
487 543412 : deriv(2, 2) = 1.0;
488 543412 : return deriv;
489 : }
490 :
491 : RankFourTensor
492 0 : CappedWeakPlaneStressUpdate::d2pdstress2(const RankTwoTensor & /*stress*/) const
493 : {
494 0 : return RankFourTensor();
495 : }
496 :
497 : RankTwoTensor
498 537428 : CappedWeakPlaneStressUpdate::dqdstress(const RankTwoTensor & stress) const
499 : {
500 537428 : RankTwoTensor deriv = RankTwoTensor();
501 537428 : const Real q = std::sqrt(Utility::pow<2>(stress(2, 0)) + Utility::pow<2>(stress(2, 1)));
502 537428 : if (q > 0.0)
503 : {
504 524163 : deriv(2, 0) = deriv(0, 2) = 0.5 * stress(2, 0) / q;
505 524163 : 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 13265 : deriv(2, 0) = deriv(0, 2) = 0.5;
511 13265 : deriv(2, 1) = deriv(1, 2) = 0.5;
512 : }
513 537428 : 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 : }
|
