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