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 1180 : CappedWeakPlaneStressUpdate::validParams()
18 : {
19 1180 : InputParameters params = TwoParameterPlasticityStressUpdate::validParams();
20 1180 : params.addClassDescription("Capped weak-plane plasticity stress calculator");
21 2360 : params.addRequiredParam<UserObjectName>(
22 : "cohesion",
23 : "A SolidMechanicsHardening UserObject that defines hardening of the cohesion. "
24 : "Physically the cohesion should not be negative.");
25 2360 : 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 2360 : 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 2360 : 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 2360 : 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 2360 : 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 2360 : params.addParam<bool>("perfect_guess",
49 2360 : 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 1180 : return params;
55 0 : }
56 :
57 888 : CappedWeakPlaneStressUpdate::CappedWeakPlaneStressUpdate(const InputParameters & parameters)
58 : : TwoParameterPlasticityStressUpdate(parameters, 3, 2),
59 888 : _cohesion(getUserObject<SolidMechanicsHardeningModel>("cohesion")),
60 888 : _tan_phi(getUserObject<SolidMechanicsHardeningModel>("tan_friction_angle")),
61 888 : _tan_psi(getUserObject<SolidMechanicsHardeningModel>("tan_dilation_angle")),
62 888 : _tstrength(getUserObject<SolidMechanicsHardeningModel>("tensile_strength")),
63 888 : _cstrength(getUserObject<SolidMechanicsHardeningModel>("compressive_strength")),
64 1776 : _small_smoother2(Utility::pow<2>(getParam<Real>("tip_smoother"))),
65 1776 : _perfect_guess(getParam<bool>("perfect_guess")),
66 888 : _stress_return_type(StressReturnType::nothing_special),
67 888 : _in_trial02(0.0),
68 888 : _in_trial12(0.0),
69 888 : _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 888 : 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 884 : 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 882 : if (_cohesion.value(0) < 0)
80 2 : mooseError("CappedWeakPlaneStressUpdate: Weak-plane cohesion must not be negative");
81 880 : 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 878 : }
85 :
86 : void
87 1225140 : CappedWeakPlaneStressUpdate::initializeReturnProcess()
88 : {
89 1225140 : _stress_return_type = StressReturnType::nothing_special;
90 1225140 : }
91 :
92 : void
93 435984 : CappedWeakPlaneStressUpdate::finalizeReturnProcess(const RankTwoTensor & /*rotation_increment*/)
94 : {
95 435984 : _stress_return_type = StressReturnType::nothing_special;
96 435984 : }
97 :
98 : void
99 435588 : 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 435588 : if (yf[1] >= 0)
108 236578 : _stress_return_type = StressReturnType::no_compression;
109 199010 : else if (yf[2] >= 0)
110 72288 : _stress_return_type = StressReturnType::no_tension;
111 :
112 : // The following are useful for the Cosserat case too
113 435588 : _in_trial02 = stress_trial(0, 2);
114 435588 : _in_trial12 = stress_trial(1, 2);
115 435588 : _in_q_trial = q_trial;
116 435588 : }
117 :
118 : void
119 1507012 : CappedWeakPlaneStressUpdate::computePQ(const RankTwoTensor & stress, Real & p, Real & q) const
120 : {
121 1507012 : p = stress(2, 2);
122 : // Because the following is not explicitly symmeterised, it is useful for the Cosserat case too
123 1507012 : q = std::sqrt(Utility::pow<2>(stress(0, 2)) + Utility::pow<2>(stress(1, 2)));
124 1507012 : }
125 :
126 : void
127 435588 : CappedWeakPlaneStressUpdate::setEppEqq(const RankFourTensor & Eijkl, Real & Epp, Real & Eqq) const
128 : {
129 435588 : Epp = Eijkl(2, 2, 2, 2);
130 435588 : Eqq = Eijkl(0, 2, 0, 2);
131 435588 : }
132 :
133 : void
134 429136 : 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 429136 : stress = stress_trial;
144 429136 : 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 429136 : stress(0, 0) -= Eijkl(2, 2, 0, 0) * gaE / _Epp * smoothed_q.dg[0];
148 429136 : stress(1, 1) -= Eijkl(2, 2, 1, 1) * gaE / _Epp * smoothed_q.dg[0];
149 429136 : if (_in_q_trial == 0.0)
150 10112 : stress(2, 0) = stress(2, 1) = stress(0, 2) = stress(1, 2) = 0.0;
151 : else
152 : {
153 419024 : stress(2, 0) = stress(0, 2) = _in_trial02 * q_ok / _in_q_trial;
154 419024 : stress(2, 1) = stress(1, 2) = _in_trial12 * q_ok / _in_q_trial;
155 : }
156 429136 : }
157 :
158 : void
159 1004502 : 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 1004502 : const Real tanpsi = _tan_psi.value(intnl[0]);
167 1004502 : dintnl[0][0] = 0.0;
168 1004502 : dintnl[0][1] = -1.0 / _Eqq;
169 1004502 : dintnl[1][0] = -1.0 / _Epp;
170 1004502 : dintnl[1][1] =
171 1004502 : tanpsi / _Eqq - (q_trial - q) * _tan_psi.derivative(intnl[0]) * dintnl[0][1] / _Eqq;
172 1004502 : }
173 :
174 : void
175 96096 : 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 96096 : cto = Eijkl;
188 96096 : if (!compute_full_tangent_operator)
189 : return;
190 :
191 96096 : const Real Ezzzz = Eijkl(2, 2, 2, 2);
192 96096 : const Real Ezxzx = Eijkl(2, 0, 2, 0);
193 96096 : const Real tanpsi = _tan_psi.value(_intnl[_qp][0]);
194 96096 : const Real dintnl0_dq = -1.0 / Ezxzx;
195 : const Real dintnl0_dqt = 1.0 / Ezxzx;
196 96096 : const Real dintnl1_dp = -1.0 / Ezzzz;
197 : const Real dintnl1_dpt = 1.0 / Ezzzz;
198 : const Real dintnl1_dq =
199 96096 : tanpsi / Ezxzx - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dq / Ezxzx;
200 : const Real dintnl1_dqt =
201 96096 : -tanpsi / Ezxzx - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dqt / Ezxzx;
202 :
203 384384 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
204 : {
205 288288 : const Real dpt_depii = Eijkl(2, 2, i, i);
206 288288 : cto(2, 2, i, i) = _dp_dpt * dpt_depii;
207 : const Real poisson_effect =
208 288288 : Eijkl(2, 2, 0, 0) / Ezzzz *
209 288288 : (_dgaE_dpt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dpt +
210 288288 : gaE * smoothed_q.d2g[0][1] * _dq_dpt +
211 288288 : gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dq * _dq_dpt) +
212 288288 : gaE * smoothed_q.d2g_di[0][1] *
213 288288 : (dintnl1_dpt + dintnl1_dp * _dp_dpt + dintnl1_dq * _dq_dpt)) *
214 288288 : dpt_depii;
215 288288 : cto(0, 0, i, i) -= poisson_effect;
216 288288 : cto(1, 1, i, i) -= poisson_effect;
217 288288 : if (q_trial > 0.0)
218 : {
219 281436 : cto(2, 0, i, i) = cto(0, 2, i, i) = _in_trial02 / q_trial * _dq_dpt * dpt_depii;
220 281436 : 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 96096 : -Eijkl(2, 2, 0, 0) / Ezzzz *
226 96096 : (_dgaE_dqt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dqt +
227 96096 : gaE * smoothed_q.d2g[0][1] * _dq_dqt +
228 96096 : gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dqt + dintnl0_dq * _dq_dqt) +
229 96096 : gaE * smoothed_q.d2g_di[0][1] * (dintnl1_dqt + dintnl1_dp * _dp_dqt + dintnl1_dq * _dq_dqt));
230 :
231 96096 : const Real dqt_dep20 = (q_trial == 0.0 ? 1.0 : _in_trial02 / q_trial) * Eijkl(2, 0, 2, 0);
232 96096 : cto(2, 2, 2, 0) = cto(2, 2, 0, 2) = _dp_dqt * dqt_dep20;
233 96096 : cto(0, 0, 2, 0) = cto(0, 0, 0, 2) = cto(1, 1, 2, 0) = cto(1, 1, 0, 2) =
234 96096 : poisson_effect * dqt_dep20;
235 96096 : if (q_trial > 0.0)
236 : {
237 : // for q_trial=0, Jacobian_mult is just given by the elastic case
238 93812 : cto(0, 2, 0, 2) = cto(2, 0, 0, 2) = cto(0, 2, 2, 0) = cto(2, 0, 2, 0) =
239 93812 : Eijkl(2, 0, 2, 0) * q / q_trial +
240 93812 : _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
241 93812 : cto(1, 2, 0, 2) = cto(2, 1, 0, 2) = cto(1, 2, 2, 0) = cto(2, 1, 2, 0) =
242 93812 : _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
243 : }
244 :
245 96096 : const Real dqt_dep21 = (q_trial == 0.0 ? 1.0 : _in_trial12 / q_trial) * Eijkl(2, 1, 2, 1);
246 96096 : cto(2, 2, 2, 1) = cto(2, 2, 1, 2) = _dp_dqt * dqt_dep21;
247 96096 : cto(0, 0, 2, 1) = cto(0, 0, 1, 2) = cto(1, 1, 2, 1) = cto(1, 1, 1, 2) =
248 96096 : poisson_effect * dqt_dep21;
249 96096 : if (q_trial > 0.0)
250 : {
251 : // for q_trial=0, Jacobian_mult is just given by the elastic case
252 93812 : cto(0, 2, 1, 2) = cto(2, 0, 1, 2) = cto(0, 2, 2, 1) = cto(2, 0, 2, 1) =
253 93812 : _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
254 93812 : cto(1, 2, 1, 2) = cto(2, 1, 1, 2) = cto(1, 2, 2, 1) = cto(2, 1, 2, 1) =
255 93812 : Eijkl(2, 1, 2, 1) * q / q_trial +
256 93812 : _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
257 : }
258 : }
259 :
260 : void
261 1661124 : CappedWeakPlaneStressUpdate::yieldFunctionValues(Real p,
262 : Real q,
263 : const std::vector<Real> & intnl,
264 : std::vector<Real> & yf) const
265 : {
266 3322248 : yf[0] = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * _tan_phi.value(intnl[0]) -
267 1661124 : _cohesion.value(intnl[0]);
268 :
269 1661124 : if (_stress_return_type == StressReturnType::no_tension)
270 0 : yf[1] = std::numeric_limits<Real>::lowest();
271 : else
272 1661124 : yf[1] = p - _tstrength.value(intnl[1]);
273 :
274 1661124 : if (_stress_return_type == StressReturnType::no_compression)
275 0 : yf[2] = std::numeric_limits<Real>::lowest();
276 : else
277 1661124 : yf[2] = -p - _cstrength.value(intnl[1]);
278 1661124 : }
279 :
280 : void
281 1365550 : 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 2731100 : all_q[0].f = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * _tan_phi.value(intnl[0]) -
288 1365550 : _cohesion.value(intnl[0]);
289 1365550 : if (_stress_return_type == StressReturnType::no_tension)
290 285312 : all_q[1].f = std::numeric_limits<Real>::lowest();
291 : else
292 1080238 : all_q[1].f = p - _tstrength.value(intnl[1]);
293 1365550 : if (_stress_return_type == StressReturnType::no_compression)
294 619460 : all_q[2].f = std::numeric_limits<Real>::lowest();
295 : else
296 746090 : all_q[2].f = -p - _cstrength.value(intnl[1]);
297 :
298 : // d(yield Function)/d(p, q)
299 : // derivatives wrt p
300 1365550 : all_q[0].df[0] = _tan_phi.value(intnl[0]);
301 1365550 : all_q[1].df[0] = 1.0;
302 1365550 : all_q[2].df[0] = -1.0;
303 :
304 : // derivatives wrt q
305 1365550 : if (_small_smoother2 == 0.0)
306 179482 : all_q[0].df[1] = 1.0;
307 : else
308 1186068 : all_q[0].df[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
309 1365550 : all_q[1].df[1] = 0.0;
310 1365550 : all_q[2].df[1] = 0.0;
311 :
312 : // d(yield Function)/d(intnl)
313 : // derivatives wrt intnl[0] (shear plastic strain)
314 1365550 : all_q[0].df_di[0] = p * _tan_phi.derivative(intnl[0]) - _cohesion.derivative(intnl[0]);
315 1365550 : all_q[1].df_di[0] = 0.0;
316 1365550 : all_q[2].df_di[0] = 0.0;
317 : // derivatives wrt intnl[q] (tensile plastic strain)
318 1365550 : all_q[0].df_di[1] = 0.0;
319 1365550 : all_q[1].df_di[1] = -_tstrength.derivative(intnl[1]);
320 1365550 : all_q[2].df_di[1] = -_cstrength.derivative(intnl[1]);
321 :
322 : // d(flowPotential)/d(p, q)
323 : // derivatives wrt p
324 1365550 : all_q[0].dg[0] = _tan_psi.value(intnl[0]);
325 1365550 : all_q[1].dg[0] = 1.0;
326 1365550 : all_q[2].dg[0] = -1.0;
327 : // derivatives wrt q
328 1365550 : if (_small_smoother2 == 0.0)
329 179482 : all_q[0].dg[1] = 1.0;
330 : else
331 1186068 : all_q[0].dg[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
332 1365550 : all_q[1].dg[1] = 0.0;
333 1365550 : all_q[2].dg[1] = 0.0;
334 :
335 : // d2(flowPotential)/d(p, q)/d(intnl)
336 : // d(dg/dp)/dintnl[0]
337 1365550 : all_q[0].d2g_di[0][0] = _tan_psi.derivative(intnl[0]);
338 1365550 : all_q[1].d2g_di[0][0] = 0.0;
339 1365550 : all_q[2].d2g_di[0][0] = 0.0;
340 : // d(dg/dp)/dintnl[1]
341 1365550 : all_q[0].d2g_di[0][1] = 0.0;
342 1365550 : all_q[1].d2g_di[0][1] = 0.0;
343 1365550 : all_q[2].d2g_di[0][1] = 0.0;
344 : // d(dg/dq)/dintnl[0]
345 1365550 : all_q[0].d2g_di[1][0] = 0.0;
346 1365550 : all_q[1].d2g_di[1][0] = 0.0;
347 1365550 : all_q[2].d2g_di[1][0] = 0.0;
348 : // d(dg/dq)/dintnl[1]
349 1365550 : all_q[0].d2g_di[1][1] = 0.0;
350 1365550 : all_q[1].d2g_di[1][1] = 0.0;
351 1365550 : 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 1365550 : all_q[0].d2g[0][0] = 0.0;
356 1365550 : all_q[1].d2g[0][0] = 0.0;
357 1365550 : all_q[2].d2g[0][0] = 0.0;
358 : // d(dg/dp)/dq
359 1365550 : all_q[0].d2g[0][1] = 0.0;
360 1365550 : all_q[1].d2g[0][1] = 0.0;
361 1365550 : all_q[2].d2g[0][1] = 0.0;
362 : // d(dg/dq)/dp
363 1365550 : all_q[0].d2g[1][0] = 0.0;
364 1365550 : all_q[1].d2g[1][0] = 0.0;
365 1365550 : all_q[2].d2g[1][0] = 0.0;
366 : // d(dg/dq)/dq
367 1365550 : if (_small_smoother2 == 0.0)
368 179482 : all_q[0].d2g[1][1] = 0.0;
369 : else
370 1186068 : all_q[0].d2g[1][1] = _small_smoother2 / std::pow(Utility::pow<2>(q) + _small_smoother2, 1.5);
371 1365550 : all_q[1].d2g[1][1] = 0.0;
372 1365550 : all_q[2].d2g[1][1] = 0.0;
373 1365550 : }
374 :
375 : void
376 436068 : 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 436068 : const Real tanpsi = _tan_psi.value(intnl_old[0]);
385 :
386 436068 : if (!_perfect_guess)
387 : {
388 33920 : p = p_trial;
389 33920 : q = q_trial;
390 33920 : gaE = 0.0;
391 : }
392 : else
393 : {
394 402148 : const Real coh = _cohesion.value(intnl_old[0]);
395 402148 : const Real tanphi = _tan_phi.value(intnl_old[0]);
396 402148 : const Real tens = _tstrength.value(intnl_old[1]);
397 402148 : const Real comp = _cstrength.value(intnl_old[1]);
398 402148 : const Real q_at_T = coh - tens * tanphi;
399 402148 : const Real q_at_C = coh + comp * tanphi;
400 :
401 402148 : if ((p_trial >= tens) && (q_trial <= q_at_T))
402 : {
403 : // pure tensile failure
404 156978 : p = tens;
405 156978 : q = q_trial;
406 156978 : gaE = p_trial - p;
407 : }
408 245170 : else if ((p_trial <= -comp) && (q_trial <= q_at_C))
409 : {
410 : // pure compressive failure
411 30672 : p = -comp;
412 30672 : q = q_trial;
413 30672 : gaE = p - p_trial;
414 : }
415 : else
416 : {
417 : // shear failure or a mixture
418 : // Calculate ga assuming a pure shear failure
419 214498 : const Real ga = (q_trial + p_trial * tanphi - coh) / (_Eqq + _Epp * tanphi * tanpsi);
420 214498 : 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 22032 : p = p_trial;
425 22032 : q = q_trial;
426 22032 : gaE = 0.0;
427 : }
428 : else
429 : {
430 192466 : q = q_trial - _Eqq * ga;
431 192466 : 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 5970 : q = 0.0;
436 5970 : p = coh / tanphi;
437 5970 : gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
438 : }
439 186496 : else if (q <= q_at_T)
440 : {
441 : // pure shear is incorrect: mixture of tension and shear is correct
442 51536 : q = q_at_T;
443 51536 : p = tens;
444 51536 : gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
445 : }
446 134960 : else if (q >= q_at_C)
447 : {
448 : // pure shear is incorrect: mixture of compression and shear is correct
449 26928 : q = q_at_C;
450 26928 : p = -comp;
451 26928 : if (p - p_trial < _Epp * tanpsi * (q_trial - q) / _Eqq)
452 : // trial point is sitting almost directly above corner
453 16128 : gaE = (q_trial - q) * _Epp / _Eqq;
454 : else
455 : // trial point is sitting to the left of the corner
456 10800 : gaE = (p - p_trial) / tanpsi;
457 : }
458 : else
459 : {
460 : // pure shear was correct
461 108032 : p = p_trial - _Epp * ga * tanpsi;
462 108032 : gaE = ga * _Epp;
463 : }
464 : }
465 : }
466 : }
467 436068 : setIntnlValues(p_trial, q_trial, p, q, intnl_old, intnl);
468 436068 : }
469 :
470 : void
471 1801534 : 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 1801534 : intnl[0] = intnl_old[0] + (q_trial - q) / _Eqq;
479 1801534 : const Real tanpsi = _tan_psi.value(intnl[0]);
480 1801534 : intnl[1] = intnl_old[1] + (p_trial - p) / _Epp - (q_trial - q) * tanpsi / _Eqq;
481 1801534 : }
482 :
483 : RankTwoTensor
484 435504 : CappedWeakPlaneStressUpdate::dpdstress(const RankTwoTensor & /*stress*/) const
485 : {
486 435504 : RankTwoTensor deriv = RankTwoTensor();
487 435504 : deriv(2, 2) = 1.0;
488 435504 : return deriv;
489 : }
490 :
491 : RankFourTensor
492 0 : CappedWeakPlaneStressUpdate::d2pdstress2(const RankTwoTensor & /*stress*/) const
493 : {
494 0 : return RankFourTensor();
495 : }
496 :
497 : RankTwoTensor
498 429616 : CappedWeakPlaneStressUpdate::dqdstress(const RankTwoTensor & stress) const
499 : {
500 429616 : RankTwoTensor deriv = RankTwoTensor();
501 429616 : const Real q = std::sqrt(Utility::pow<2>(stress(2, 0)) + Utility::pow<2>(stress(2, 1)));
502 429616 : if (q > 0.0)
503 : {
504 419504 : deriv(2, 0) = deriv(0, 2) = 0.5 * stress(2, 0) / q;
505 419504 : 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 10112 : deriv(2, 0) = deriv(0, 2) = 0.5;
511 10112 : deriv(2, 1) = deriv(1, 2) = 0.5;
512 : }
513 429616 : 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 : }
|
