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 "TensorMechanicsPlasticDruckerPragerHyperbolic.h"
11 : #include "RankFourTensor.h"
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("TensorMechanicsApp", TensorMechanicsPlasticDruckerPragerHyperbolic);
15 :
16 : InputParameters
17 138 : TensorMechanicsPlasticDruckerPragerHyperbolic::validParams()
18 : {
19 138 : InputParameters params = TensorMechanicsPlasticDruckerPrager::validParams();
20 276 : params.addParam<bool>("use_custom_returnMap",
21 276 : true,
22 : "Whether to use the custom returnMap "
23 : "algorithm. Set to true if you are using "
24 : "isotropic elasticity.");
25 276 : params.addParam<bool>("use_custom_cto",
26 276 : true,
27 : "Whether to use the custom consistent tangent "
28 : "operator computations. Set to true if you are "
29 : "using isotropic elasticity.");
30 138 : params.addClassDescription("J2 plasticity, associative, with hardening");
31 414 : params.addRangeCheckedParam<Real>("smoother",
32 276 : 0.0,
33 : "smoother>=0",
34 : "The cone vertex at J2=0 is smoothed. The maximum mean "
35 : "stress possible, which is Cohesion*Cot(friction_angle) for "
36 : "smoother=0, becomes (Cohesion - "
37 : "smoother/3)*Cot(friction_angle). This is a non-hardening "
38 : "parameter");
39 414 : params.addRangeCheckedParam<unsigned>(
40 : "max_iterations",
41 276 : 40,
42 : "max_iterations>0",
43 : "Maximum iterations to use in the custom return map function");
44 138 : params.addClassDescription(
45 : "Non-associative Drucker Prager plasticity with hyperbolic smoothing of the cone tip.");
46 138 : return params;
47 0 : }
48 :
49 69 : TensorMechanicsPlasticDruckerPragerHyperbolic::TensorMechanicsPlasticDruckerPragerHyperbolic(
50 69 : const InputParameters & parameters)
51 : : TensorMechanicsPlasticDruckerPrager(parameters),
52 69 : _smoother2(Utility::pow<2>(getParam<Real>("smoother"))),
53 138 : _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
54 138 : _use_custom_cto(getParam<bool>("use_custom_cto")),
55 207 : _max_iters(getParam<unsigned>("max_iterations"))
56 : {
57 69 : }
58 :
59 : Real
60 49472 : TensorMechanicsPlasticDruckerPragerHyperbolic::yieldFunction(const RankTwoTensor & stress,
61 : Real intnl) const
62 : {
63 : Real aaa;
64 : Real bbb;
65 49472 : bothAB(intnl, aaa, bbb);
66 49472 : return std::sqrt(stress.secondInvariant() + _smoother2) + stress.trace() * bbb - aaa;
67 : }
68 :
69 : RankTwoTensor
70 0 : TensorMechanicsPlasticDruckerPragerHyperbolic::df_dsig(const RankTwoTensor & stress, Real bbb) const
71 : {
72 0 : return 0.5 * stress.dsecondInvariant() / std::sqrt(stress.secondInvariant() + _smoother2) +
73 0 : stress.dtrace() * bbb;
74 : }
75 :
76 : RankFourTensor
77 0 : TensorMechanicsPlasticDruckerPragerHyperbolic::dflowPotential_dstress(const RankTwoTensor & stress,
78 : Real /*intnl*/) const
79 : {
80 0 : RankFourTensor dr_dstress;
81 0 : dr_dstress = 0.5 * stress.d2secondInvariant() / std::sqrt(stress.secondInvariant() + _smoother2);
82 0 : dr_dstress += -0.5 * 0.5 * stress.dsecondInvariant().outerProduct(stress.dsecondInvariant()) /
83 0 : std::pow(stress.secondInvariant() + _smoother2, 1.5);
84 0 : return dr_dstress;
85 : }
86 :
87 : std::string
88 0 : TensorMechanicsPlasticDruckerPragerHyperbolic::modelName() const
89 : {
90 0 : return "DruckerPragerHyperbolic";
91 : }
92 :
93 : bool
94 49472 : TensorMechanicsPlasticDruckerPragerHyperbolic::returnMap(const RankTwoTensor & trial_stress,
95 : Real intnl_old,
96 : const RankFourTensor & E_ijkl,
97 : Real ep_plastic_tolerance,
98 : RankTwoTensor & returned_stress,
99 : Real & returned_intnl,
100 : std::vector<Real> & dpm,
101 : RankTwoTensor & delta_dp,
102 : std::vector<Real> & yf,
103 : bool & trial_stress_inadmissible) const
104 : {
105 49472 : if (!(_use_custom_returnMap))
106 0 : return TensorMechanicsPlasticModel::returnMap(trial_stress,
107 : intnl_old,
108 : E_ijkl,
109 : ep_plastic_tolerance,
110 : returned_stress,
111 : returned_intnl,
112 : dpm,
113 : delta_dp,
114 : yf,
115 : trial_stress_inadmissible);
116 :
117 49472 : yf.resize(1);
118 :
119 49472 : yf[0] = yieldFunction(trial_stress, intnl_old);
120 :
121 49472 : if (yf[0] < _f_tol)
122 : {
123 : // the trial_stress is admissible
124 19072 : trial_stress_inadmissible = false;
125 19072 : return true;
126 : }
127 :
128 30400 : trial_stress_inadmissible = true;
129 30400 : const Real mu = E_ijkl(0, 1, 0, 1);
130 30400 : const Real lambda = E_ijkl(0, 0, 0, 0) - 2.0 * mu;
131 30400 : const Real bulky = 3.0 * lambda + 2.0 * mu;
132 30400 : const Real Tr_trial = trial_stress.trace();
133 30400 : const Real J2trial = trial_stress.secondInvariant();
134 :
135 : // Perform a Newton-Raphson to find dpm when
136 : // residual = (1 + dpm*mu/ll)sqrt(ll^2 - s^2) - sqrt(J2trial) = 0, with s=smoother
137 : // with ll = sqrt(J2 + s^2) = aaa - bbb*Tr(stress) = aaa - bbb(Tr_trial - p*3*bulky*bbb_flow)
138 : Real aaa;
139 : Real daaa;
140 : Real bbb;
141 : Real dbbb;
142 : Real bbb_flow;
143 : Real dbbb_flow;
144 : Real ll;
145 : Real dll;
146 : Real residual;
147 : Real jac;
148 30400 : dpm[0] = 0;
149 : unsigned int iter = 0;
150 : do
151 : {
152 208664 : bothAB(intnl_old + dpm[0], aaa, bbb);
153 208664 : dbothAB(intnl_old + dpm[0], daaa, dbbb);
154 208664 : onlyB(intnl_old + dpm[0], dilation, bbb_flow);
155 208664 : donlyB(intnl_old + dpm[0], dilation, dbbb_flow);
156 208664 : ll = aaa - bbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow);
157 208664 : dll = daaa - dbbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow) +
158 208664 : bbb * bulky * 3 * (bbb_flow + dpm[0] * dbbb_flow);
159 208664 : residual = bbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow) - aaa +
160 208664 : std::sqrt(J2trial / Utility::pow<2>(1 + dpm[0] * mu / ll) + _smoother2);
161 208664 : jac = dbbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow) -
162 208664 : bbb * bulky * 3 * (bbb_flow + dpm[0] * dbbb_flow) - daaa +
163 208664 : 0.5 * J2trial * (-2.0) * (mu / ll - dpm[0] * mu * dll / ll / ll) /
164 208664 : Utility::pow<3>(1 + dpm[0] * mu / ll) /
165 208664 : std::sqrt(J2trial / Utility::pow<2>(1.0 + dpm[0] * mu / ll) + _smoother2);
166 208664 : dpm[0] += -residual / jac;
167 208664 : if (iter > _max_iters) // not converging
168 : return false;
169 208664 : iter++;
170 208664 : } while (residual * residual > _f_tol * _f_tol);
171 :
172 : // set the returned values
173 30400 : yf[0] = 0;
174 30400 : returned_intnl = intnl_old + dpm[0];
175 :
176 30400 : bothAB(returned_intnl, aaa, bbb);
177 30400 : onlyB(returned_intnl, dilation, bbb_flow);
178 30400 : ll = aaa - bbb * (Tr_trial - dpm[0] * bulky * 3.0 * bbb_flow);
179 : returned_stress =
180 30400 : trial_stress.deviatoric() / (1.0 + dpm[0] * mu / ll); // this is the deviatoric part only
181 :
182 30400 : RankTwoTensor rij = 0.5 * returned_stress.deviatoric() /
183 : ll; // this is the derivatoric part the flow potential only
184 :
185 : // form the returned stress and the full flow potential
186 30400 : const Real returned_trace_over_3 = (aaa - ll) / bbb / 3.0;
187 121600 : for (unsigned i = 0; i < 3; ++i)
188 : {
189 91200 : returned_stress(i, i) += returned_trace_over_3;
190 91200 : rij(i, i) += bbb_flow;
191 : }
192 :
193 30400 : delta_dp = rij * dpm[0];
194 :
195 : return true;
196 : }
197 :
198 : RankFourTensor
199 30400 : TensorMechanicsPlasticDruckerPragerHyperbolic::consistentTangentOperator(
200 : const RankTwoTensor & trial_stress,
201 : Real intnl_old,
202 : const RankTwoTensor & stress,
203 : Real intnl,
204 : const RankFourTensor & E_ijkl,
205 : const std::vector<Real> & cumulative_pm) const
206 : {
207 30400 : if (!_use_custom_cto)
208 0 : return TensorMechanicsPlasticModel::consistentTangentOperator(
209 : trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
210 :
211 30400 : const Real mu = E_ijkl(0, 1, 0, 1);
212 30400 : const Real la = E_ijkl(0, 0, 0, 0) - 2.0 * mu;
213 30400 : const Real bulky = 3.0 * la + 2.0 * mu;
214 : Real bbb;
215 30400 : onlyB(intnl, friction, bbb);
216 : Real bbb_flow;
217 30400 : onlyB(intnl, dilation, bbb_flow);
218 : Real dbbb_flow;
219 30400 : donlyB(intnl, dilation, dbbb_flow);
220 30400 : const Real bbb_flow_mod = bbb_flow + cumulative_pm[0] * dbbb_flow;
221 30400 : const Real J2 = stress.secondInvariant();
222 30400 : const RankTwoTensor sij = stress.deviatoric();
223 30400 : const Real sq = std::sqrt(J2 + _smoother2);
224 :
225 : const Real one_over_h =
226 30400 : 1.0 / (-dyieldFunction_dintnl(stress, intnl) + mu * J2 / Utility::pow<2>(sq) +
227 30400 : 3.0 * bbb * bbb_flow_mod * bulky); // correct up to hard
228 :
229 : const RankTwoTensor df_dsig_timesE =
230 30400 : mu * sij / sq + bulky * bbb * RankTwoTensor(RankTwoTensor::initIdentity); // correct
231 : const RankTwoTensor rmod_timesE =
232 30400 : mu * sij / sq +
233 30400 : bulky * bbb_flow_mod * RankTwoTensor(RankTwoTensor::initIdentity); // correct up to hard
234 :
235 : const RankFourTensor coeff_ep =
236 30400 : E_ijkl - one_over_h * rmod_timesE.outerProduct(df_dsig_timesE); // correct
237 :
238 30400 : const RankFourTensor dr_dsig_timesE = -0.5 * mu * sij.outerProduct(sij) / Utility::pow<3>(sq) +
239 60800 : mu * stress.d2secondInvariant() / sq; // correct
240 : const RankTwoTensor df_dsig_E_dr_dsig =
241 30400 : 0.5 * mu * _smoother2 * sij / Utility::pow<4>(sq); // correct
242 :
243 : const RankFourTensor coeff_si =
244 30400 : RankFourTensor(RankFourTensor::initIdentitySymmetricFour) +
245 30400 : cumulative_pm[0] *
246 60800 : (dr_dsig_timesE - one_over_h * rmod_timesE.outerProduct(df_dsig_E_dr_dsig));
247 :
248 30400 : RankFourTensor s_inv;
249 : try
250 : {
251 30400 : s_inv = coeff_si.invSymm();
252 : }
253 0 : catch (MooseException & e)
254 : {
255 0 : return coeff_ep; // when coeff_si is singular return the "linear" tangent operator
256 0 : }
257 :
258 30400 : return s_inv * coeff_ep;
259 : }
260 :
261 : bool
262 0 : TensorMechanicsPlasticDruckerPragerHyperbolic::useCustomReturnMap() const
263 : {
264 0 : return _use_custom_returnMap;
265 : }
266 :
267 : bool
268 30400 : TensorMechanicsPlasticDruckerPragerHyperbolic::useCustomCTO() const
269 : {
270 30400 : return _use_custom_cto;
271 : }
|
