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 "FiniteStrainPlasticMaterial.h"
11 : #include "libmesh/utility.h"
12 :
13 : registerMooseObject("SolidMechanicsApp", FiniteStrainPlasticMaterial);
14 :
15 : InputParameters
16 72 : FiniteStrainPlasticMaterial::validParams()
17 : {
18 72 : InputParameters params = ComputeStressBase::validParams();
19 :
20 144 : params.addRequiredParam<std::vector<Real>>(
21 : "yield_stress",
22 : "Input data as pairs of equivalent plastic strain and yield stress: Should "
23 : "start with equivalent plastic strain 0");
24 144 : params.addParam<Real>("rtol", 1e-8, "Plastic strain NR tolerance");
25 144 : params.addParam<Real>("ftol", 1e-4, "Consistency condition NR tolerance");
26 144 : params.addParam<Real>("eptol", 1e-7, "Equivalent plastic strain NR tolerance");
27 72 : params.addClassDescription("Associative J2 plasticity with isotropic hardening.");
28 :
29 72 : return params;
30 0 : }
31 :
32 54 : FiniteStrainPlasticMaterial::FiniteStrainPlasticMaterial(const InputParameters & parameters)
33 : : ComputeStressBase(parameters),
34 54 : _yield_stress_vector(getParam<std::vector<Real>>("yield_stress")), // Read from input file
35 54 : _plastic_strain(declareProperty<RankTwoTensor>(_base_name + "plastic_strain")),
36 108 : _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "plastic_strain")),
37 54 : _eqv_plastic_strain(declareProperty<Real>(_base_name + "eqv_plastic_strain")),
38 108 : _eqv_plastic_strain_old(getMaterialPropertyOld<Real>(_base_name + "eqv_plastic_strain")),
39 108 : _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
40 108 : _strain_increment(getMaterialProperty<RankTwoTensor>(_base_name + "strain_increment")),
41 108 : _rotation_increment(getMaterialProperty<RankTwoTensor>(_base_name + "rotation_increment")),
42 54 : _elasticity_tensor_name(_base_name + "elasticity_tensor"),
43 54 : _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
44 108 : _rtol(getParam<Real>("rtol")),
45 108 : _ftol(getParam<Real>("ftol")),
46 108 : _eptol(getParam<Real>("eptol")),
47 54 : _deltaOuter(RankTwoTensor::Identity().times<i_, j_, k_, l_>(RankTwoTensor::Identity())),
48 108 : _deltaMixed(RankTwoTensor::Identity().times<i_, k_, j_, l_>(RankTwoTensor::Identity()))
49 : {
50 54 : }
51 :
52 : void
53 192 : FiniteStrainPlasticMaterial::initQpStatefulProperties()
54 : {
55 192 : ComputeStressBase::initQpStatefulProperties();
56 192 : _plastic_strain[_qp].zero();
57 192 : _eqv_plastic_strain[_qp] = 0.0;
58 192 : }
59 :
60 : void
61 33856 : FiniteStrainPlasticMaterial::computeQpStress()
62 : {
63 :
64 : // perform the return-mapping algorithm
65 33856 : returnMap(_stress_old[_qp],
66 33856 : _eqv_plastic_strain_old[_qp],
67 33856 : _plastic_strain_old[_qp],
68 33856 : _strain_increment[_qp],
69 33856 : _elasticity_tensor[_qp],
70 33856 : _stress[_qp],
71 33856 : _eqv_plastic_strain[_qp],
72 33856 : _plastic_strain[_qp]);
73 :
74 : // Rotate the stress tensor to the current configuration
75 33856 : _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
76 :
77 : // Rotate plastic strain tensor to the current configuration
78 33856 : _plastic_strain[_qp] =
79 33856 : _rotation_increment[_qp] * _plastic_strain[_qp] * _rotation_increment[_qp].transpose();
80 :
81 : // Calculate the elastic strain_increment
82 33856 : _elastic_strain[_qp] = _mechanical_strain[_qp] - _plastic_strain[_qp];
83 :
84 33856 : _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
85 33856 : }
86 :
87 : /**
88 : * Implements the return-map algorithm via a Newton-Raphson process.
89 : * This idea is fully explained in Simo and Hughes "Computational
90 : * Inelasticity" Springer 1997, for instance, as well as many other
91 : * books on plasticity.
92 : * The basic idea is as follows.
93 : * Given: sig_old - the stress at the start of the "time step"
94 : * plastic_strain_old - the plastic strain at the start of the "time step"
95 : * eqvpstrain_old - equivalent plastic strain at the start of the "time step"
96 : * (In general we would be given some number of internal
97 : * parameters that the yield function depends upon.)
98 : * delta_d - the prescribed strain increment for this "time step"
99 : * we want to determine the following parameters at the end of this "time step":
100 : * sig - the stress
101 : * plastic_strain - the plastic strain
102 : * eqvpstrain - the equivalent plastic strain (again, in general, we would
103 : * have an arbitrary number of internal parameters).
104 : *
105 : * To determine these parameters, introduce
106 : * the "yield function", f
107 : * the "consistency parameter", flow_incr
108 : * the "flow potential", flow_dirn_ij
109 : * the "internal potential", internalPotential (in general there are as many internalPotential
110 : * functions as there are internal parameters).
111 : * All three of f, flow_dirn_ij, and internalPotential, are functions of
112 : * sig and eqvpstrain.
113 : * To find sig, plastic_strain and eqvpstrain, we need to solve the following
114 : * resid_ij = 0
115 : * f = 0
116 : * rep = 0
117 : * This is done by using Newton-Raphson.
118 : * There are 8 equations here: six from resid_ij=0 (more generally there are nine
119 : * but in this case resid is symmetric); one from f=0; one from rep=0 (more generally, for N
120 : * internal parameters there are N of these equations).
121 : *
122 : * resid_ij = flow_incr*flow_dirn_ij - (plastic_strain - plastic_strain_old)_ij
123 : * = flow_incr*flow_dirn_ij - (E^{-1}(trial_stress - sig))_ij
124 : * Here trial_stress = E*(strain - plastic_strain_old)
125 : * sig = E*(strain - plastic_strain)
126 : * Note: flow_dirn_ij is evaluated at sig and eqvpstrain (not the old values).
127 : *
128 : * f is the yield function, evaluated at sig and eqvpstrain
129 : *
130 : * rep = -flow_incr*internalPotential - (eqvpstrain - eqvpstrain_old)
131 : * Here internalPotential are evaluated at sig and eqvpstrain (not the old values).
132 : *
133 : * The Newton-Raphson procedure has sig, flow_incr, and eqvpstrain as its
134 : * variables. Therefore we need the derivatives of resid_ij, f, and rep
135 : * with respect to these parameters
136 : *
137 : * In this associative J2 with isotropic hardening, things are a little more specialised.
138 : * (1) f = sqrt(3*s_ij*s_ij/2) - K(eqvpstrain) (this is called "isotropic hardening")
139 : * (2) associativity means that flow_dirn_ij = df/d(sig_ij) = s_ij*sqrt(3/2/(s_ij*s_ij)), and
140 : * this means that flow_dirn_ij*flow_dirn_ij = 3/2, so when resid_ij=0, we get
141 : * (plastic_strain_dot)_ij*(plastic_strain_dot)_ij = (3/2)*flow_incr^2, where
142 : * plastic_strain_dot = plastic_strain - plastic_strain_old
143 : * (3) The definition of equivalent plastic strain is through
144 : * eqvpstrain_dot = sqrt(2*plastic_strain_dot_ij*plastic_strain_dot_ij/3), so
145 : * using (2), we obtain eqvpstrain_dot = flow_incr, and this yields
146 : * internalPotential = -1 in the "rep" equation.
147 : */
148 : void
149 33856 : FiniteStrainPlasticMaterial::returnMap(const RankTwoTensor & sig_old,
150 : const Real eqvpstrain_old,
151 : const RankTwoTensor & plastic_strain_old,
152 : const RankTwoTensor & delta_d,
153 : const RankFourTensor & E_ijkl,
154 : RankTwoTensor & sig,
155 : Real & eqvpstrain,
156 : RankTwoTensor & plastic_strain)
157 : {
158 : // the yield function, must be non-positive
159 : // Newton-Raphson sets this to zero if trial stress enters inadmissible region
160 : Real f;
161 :
162 : // the consistency parameter, must be non-negative
163 : // change in plastic strain in this timestep = flow_incr*flow_potential
164 33856 : Real flow_incr = 0.0;
165 :
166 : // direction of flow defined by the potential
167 33856 : RankTwoTensor flow_dirn;
168 :
169 : // Newton-Raphson sets this zero
170 : // resid_ij = flow_incr*flow_dirn_ij - (plastic_strain - plastic_strain_old)
171 33856 : RankTwoTensor resid;
172 :
173 : // Newton-Raphson sets this zero
174 : // rep = -flow_incr*internalPotential - (eqvpstrain - eqvpstrain_old)
175 : Real rep;
176 :
177 : // change in the stress (sig) in a Newton-Raphson iteration
178 33856 : RankTwoTensor ddsig;
179 :
180 : // change in the consistency parameter in a Newton-Raphson iteration
181 33856 : Real dflow_incr = 0.0;
182 :
183 : // change in equivalent plastic strain in one Newton-Raphson iteration
184 : Real deqvpstrain = 0.0;
185 :
186 : // convenience variable that holds the change in plastic strain incurred during the return
187 : // delta_dp = plastic_strain - plastic_strain_old
188 : // delta_dp = E^{-1}*(trial_stress - sig), where trial_stress = E*(strain - plastic_strain_old)
189 33856 : RankTwoTensor delta_dp;
190 :
191 : // d(yieldFunction)/d(stress)
192 33856 : RankTwoTensor df_dsig;
193 :
194 : // d(resid_ij)/d(sigma_kl)
195 33856 : RankFourTensor dr_dsig;
196 :
197 : // dr_dsig_inv_ijkl*dr_dsig_klmn = 0.5*(de_ij de_jn + de_ij + de_jm), where de_ij = 1 if i=j, but
198 : // zero otherwise
199 33856 : RankFourTensor dr_dsig_inv;
200 :
201 : // d(yieldFunction)/d(eqvpstrain)
202 : Real fq;
203 :
204 : // yield stress at the start of this "time step" (ie, evaluated with
205 : // eqvpstrain_old). It is held fixed during the Newton-Raphson return,
206 : // even if eqvpstrain != eqvpstrain_old.
207 : Real yield_stress;
208 :
209 : // measures of whether the Newton-Raphson process has converged
210 : Real err1, err2, err3;
211 :
212 : // number of Newton-Raphson iterations performed
213 : unsigned int iter = 0;
214 :
215 : // maximum number of Newton-Raphson iterations allowed
216 : unsigned int maxiter = 100;
217 :
218 : // plastic loading occurs if yieldFunction > toly
219 : Real toly = 1.0e-8;
220 :
221 : // Assume this strain increment does not induce any plasticity
222 : // This is the elastic-predictor
223 33856 : sig = sig_old + E_ijkl * delta_d; // the trial stress
224 33856 : eqvpstrain = eqvpstrain_old;
225 33856 : plastic_strain = plastic_strain_old;
226 :
227 33856 : yield_stress = getYieldStress(eqvpstrain); // yield stress at this equivalent plastic strain
228 33856 : if (yieldFunction(sig, yield_stress) > toly)
229 : {
230 : // the sig just calculated is inadmissable. We must return to the yield surface.
231 : // This is done iteratively, using a Newton-Raphson process.
232 :
233 : delta_dp.zero();
234 :
235 33856 : sig = sig_old + E_ijkl * delta_d; // this is the elastic predictor
236 :
237 33856 : flow_dirn = flowPotential(sig);
238 :
239 33856 : resid = flow_dirn * flow_incr - delta_dp; // Residual 1 - refer Hughes Simo
240 33856 : f = yieldFunction(sig, yield_stress);
241 33856 : rep = -eqvpstrain + eqvpstrain_old - flow_incr * internalPotential(); // Residual 3 rep=0
242 :
243 33856 : err1 = resid.L2norm();
244 : err2 = std::abs(f);
245 : err3 = std::abs(rep);
246 :
247 160832 : while ((err1 > _rtol || err2 > _ftol || err3 > _eptol) &&
248 : iter < maxiter) // Stress update iteration (hardness fixed)
249 : {
250 126976 : iter++;
251 :
252 126976 : df_dsig = dyieldFunction_dstress(sig);
253 126976 : getJac(sig, E_ijkl, flow_incr, dr_dsig); // gets dr_dsig = d(resid_ij)/d(sig_kl)
254 126976 : fq = dyieldFunction_dinternal(eqvpstrain); // d(f)/d(eqvpstrain)
255 :
256 : /**
257 : * The linear system is
258 : * ( dr_dsig flow_dirn 0 )( ddsig ) ( - resid )
259 : * ( df_dsig 0 fq )( dflow_incr ) = ( - f )
260 : * ( 0 1 -1 )( deqvpstrain ) ( - rep )
261 : * The zeroes are: d(resid_ij)/d(eqvpstrain) = flow_dirn*d(df/d(sig_ij))/d(eqvpstrain) = 0
262 : * and df/d(flow_dirn) = 0 (this is always true, even for general hardening and
263 : * non-associative)
264 : * and d(rep)/d(sig_ij) = -flow_incr*d(internalPotential)/d(sig_ij) = 0
265 : */
266 :
267 126976 : dr_dsig_inv = dr_dsig.invSymm();
268 :
269 : /**
270 : * Because of the zeroes and ones, the linear system is not impossible to
271 : * solve by hand.
272 : * NOTE: andy believes there was originally a sign-error in the next line. The
273 : * next line is unchanged, however andy's definition of fq is negative of
274 : * the original definition of fq. andy can't see any difference in any tests!
275 : */
276 126976 : dflow_incr = (f - df_dsig.doubleContraction(dr_dsig_inv * resid) + fq * rep) /
277 126976 : (df_dsig.doubleContraction(dr_dsig_inv * flow_dirn) - fq);
278 126976 : ddsig =
279 253952 : dr_dsig_inv *
280 126976 : (-resid -
281 126976 : flow_dirn * dflow_incr); // from solving the top row of linear system, given dflow_incr
282 126976 : deqvpstrain =
283 : rep + dflow_incr; // from solving the bottom row of linear system, given dflow_incr
284 :
285 : // update the variables
286 126976 : flow_incr += dflow_incr;
287 126976 : delta_dp -= E_ijkl.invSymm() * ddsig;
288 126976 : sig += ddsig;
289 126976 : eqvpstrain += deqvpstrain;
290 :
291 : // evaluate the RHS equations ready for next Newton-Raphson iteration
292 126976 : flow_dirn = flowPotential(sig);
293 126976 : resid = flow_dirn * flow_incr - delta_dp;
294 126976 : f = yieldFunction(sig, yield_stress);
295 126976 : rep = -eqvpstrain + eqvpstrain_old - flow_incr * internalPotential();
296 :
297 126976 : err1 = resid.L2norm();
298 : err2 = std::abs(f);
299 : err3 = std::abs(rep);
300 : }
301 :
302 33856 : if (iter >= maxiter)
303 0 : mooseError("Constitutive failure");
304 :
305 33856 : plastic_strain += delta_dp;
306 : }
307 33856 : }
308 :
309 : Real
310 194688 : FiniteStrainPlasticMaterial::yieldFunction(const RankTwoTensor & stress, const Real yield_stress)
311 : {
312 194688 : return getSigEqv(stress) - yield_stress;
313 : }
314 :
315 : RankTwoTensor
316 541760 : FiniteStrainPlasticMaterial::dyieldFunction_dstress(const RankTwoTensor & sig)
317 : {
318 541760 : RankTwoTensor deriv = sig.dsecondInvariant();
319 541760 : deriv *= std::sqrt(3.0 / sig.secondInvariant()) / 2.0;
320 541760 : return deriv;
321 : }
322 :
323 : Real
324 126976 : FiniteStrainPlasticMaterial::dyieldFunction_dinternal(const Real equivalent_plastic_strain)
325 : {
326 126976 : return -getdYieldStressdPlasticStrain(equivalent_plastic_strain);
327 : }
328 :
329 : RankTwoTensor
330 287808 : FiniteStrainPlasticMaterial::flowPotential(const RankTwoTensor & sig)
331 : {
332 287808 : return dyieldFunction_dstress(sig); // this plasticity model assumes associative flow
333 : }
334 :
335 : Real
336 160832 : FiniteStrainPlasticMaterial::internalPotential()
337 : {
338 160832 : return -1;
339 : }
340 :
341 : Real
342 321664 : FiniteStrainPlasticMaterial::getSigEqv(const RankTwoTensor & stress)
343 : {
344 321664 : return std::sqrt(3 * stress.secondInvariant());
345 : }
346 :
347 : // Jacobian for stress update algorithm
348 : void
349 126976 : FiniteStrainPlasticMaterial::getJac(const RankTwoTensor & sig,
350 : const RankFourTensor & E_ijkl,
351 : Real flow_incr,
352 : RankFourTensor & dresid_dsig)
353 : {
354 126976 : RankTwoTensor sig_dev, df_dsig, flow_dirn;
355 126976 : RankTwoTensor dfi_dft, dfi_dsig;
356 126976 : RankFourTensor dft_dsig, dfd_dft, dfd_dsig;
357 : Real sig_eqv;
358 : Real f1, f2, f3;
359 126976 : RankFourTensor temp;
360 :
361 126976 : sig_dev = sig.deviatoric();
362 126976 : sig_eqv = getSigEqv(sig);
363 126976 : df_dsig = dyieldFunction_dstress(sig);
364 126976 : flow_dirn = flowPotential(sig);
365 :
366 126976 : f1 = 3.0 / (2.0 * sig_eqv);
367 126976 : f2 = f1 / 3.0;
368 126976 : f3 = 9.0 / (4.0 * Utility::pow<3>(sig_eqv));
369 :
370 126976 : dft_dsig = f1 * _deltaMixed - f2 * _deltaOuter - f3 * sig_dev.outerProduct(sig_dev);
371 :
372 126976 : dfd_dsig = dft_dsig;
373 126976 : dresid_dsig = E_ijkl.invSymm() + dfd_dsig * flow_incr;
374 126976 : }
375 :
376 : // Obtain yield stress for a given equivalent plastic strain (input)
377 : Real
378 33856 : FiniteStrainPlasticMaterial::getYieldStress(const Real eqpe)
379 : {
380 : unsigned nsize;
381 :
382 33856 : nsize = _yield_stress_vector.size();
383 :
384 33856 : if (_yield_stress_vector[0] > 0.0 || nsize % 2 > 0) // Error check for input inconsitency
385 0 : mooseError("Error in yield stress input: Should be a vector with eqv plastic strain and yield "
386 : "stress pair values.\n");
387 :
388 : unsigned int ind = 0;
389 : Real tol = 1e-8;
390 :
391 111872 : while (ind < nsize)
392 : {
393 111872 : if (std::abs(eqpe - _yield_stress_vector[ind]) < tol)
394 3424 : return _yield_stress_vector[ind + 1];
395 :
396 108448 : if (ind + 2 < nsize)
397 : {
398 108448 : if (eqpe > _yield_stress_vector[ind] && eqpe < _yield_stress_vector[ind + 2])
399 30432 : return _yield_stress_vector[ind + 1] +
400 30432 : (eqpe - _yield_stress_vector[ind]) /
401 30432 : (_yield_stress_vector[ind + 2] - _yield_stress_vector[ind]) *
402 30432 : (_yield_stress_vector[ind + 3] - _yield_stress_vector[ind + 1]);
403 : }
404 : else
405 0 : return _yield_stress_vector[nsize - 1];
406 :
407 : ind += 2;
408 : }
409 :
410 : return 0.0;
411 : }
412 :
413 : Real
414 126976 : FiniteStrainPlasticMaterial::getdYieldStressdPlasticStrain(const Real eqpe)
415 : {
416 : unsigned nsize;
417 :
418 126976 : nsize = _yield_stress_vector.size();
419 :
420 126976 : if (_yield_stress_vector[0] > 0.0 || nsize % 2 > 0) // Error check for input inconsitency
421 0 : mooseError("Error in yield stress input: Should be a vector with eqv plastic strain and yield "
422 : "stress pair values.\n");
423 :
424 : unsigned int ind = 0;
425 :
426 440000 : while (ind < nsize)
427 : {
428 440000 : if (ind + 2 < nsize)
429 : {
430 440000 : if (eqpe >= _yield_stress_vector[ind] && eqpe < _yield_stress_vector[ind + 2])
431 126976 : return (_yield_stress_vector[ind + 3] - _yield_stress_vector[ind + 1]) /
432 126976 : (_yield_stress_vector[ind + 2] - _yield_stress_vector[ind]);
433 : }
434 : else
435 : return 0.0;
436 :
437 : ind += 2;
438 : }
439 :
440 : return 0.0;
441 : }
|
