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 : #pragma once
11 :
12 : #include "StressUpdateBase.h"
13 :
14 : #include <array>
15 :
16 : /**
17 : * MultiParameterPlasticityStressUpdate performs the return-map
18 : * algorithm and associated stress updates for plastic
19 : * models where the yield function and flow directions
20 : * depend on multiple parameters (called "stress_params" in the
21 : * documentation and sp in the code) that
22 : * are themselves functions of stress.
23 : *
24 : * Let the stress_params be S = {S[0], S[1], ... S[N-1]}
25 : * and define _num_sp = N.
26 : *
27 : * For instance, CappedDruckerPrager plasticity has _num_sp = 2
28 : * and defines
29 : * S[0] = p = stress.trace()
30 : * S[1] = q = sqrt(stress.secondInvariant())
31 : *
32 : * The point of the stress_params is to reduce the number of
33 : * degrees of freedom involved in the return-map algorithm
34 : * (for computational efficiency as all the intensive computation
35 : * occurs in "stress_param space") and to allow users to
36 : * write their yield functions, etc, in forms that are
37 : * clear to them to avoid errors.
38 : *
39 : * Derived classes can describe plasticity via multiple yield
40 : * functions. The number of yield function is _num_yf.
41 : * The yield function(s) and flow potential(s)
42 : * must be functions of the stress_params.
43 : *
44 : * Derived classes can use any number of internal parameters.
45 : * The number of internal parameters is _num_intnl. The
46 : * derived classes must define the evolution of these internal
47 : * parameters, which are typically functions of certain
48 : * "distances" in the return-mapping procedure.
49 : *
50 : * For instance, CappedDruckerPrager plasticity has three
51 : * yield functions (a smoothed DruckerPrager cone, a tensile failure
52 : * envelope, and a compressive failure envelope) and two
53 : * internal parameters (plastic shear strain and plastic tensile strain).
54 : * The strength parameters (cohesion, friction angle,
55 : * dilation angle, tensile strength, compressive strength)
56 : * are functions of these internal parameters.
57 : *
58 : * A novel smoothing procedure is used to smooth the derived class's
59 : * yield functions into a single, smooth yield function. The
60 : * smoothing is also performed for the flow potential. All
61 : * return-mapping, etc, processes are performed on this single
62 : * yield function and flow potential.
63 : *
64 : * The return-map algorithm implements the updateState method of
65 : * StressUpdateBase. In particular, the following system of
66 : * equations is solved:
67 : * 0 = _rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j]
68 : * 0 = _rhs[1] = S[1] - S[1]^trial + ga * E[1, j] * dg/dS[j]
69 : * ...
70 : * 0 = _rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, j] * dg/dS[j]
71 : * 0 = _rhs[N] = f(S, intnl)
72 : * (here there is an implied sum over j)
73 : * for the _num_sp stress_params S, and the single scalar ga
74 : * (actually I solve for gaE = ga * _En, where _En is a normalising
75 : * factor so that gaE is of similar magnitude to the S variables).
76 : * In general E[a, b] = dS[a]/dstress(i, j) E(i, j, k, l) dS[b]/dstress(k, l),
77 : * and in the code below it is assumed that the E[a, b] are independent
78 : * of the stress_params, but adding a dependence would only require
79 : * some Jacobian terms to be modified.
80 : * There are N+1 rhs equations to solve. Here f is
81 : * the smoothed yield function, so the last equation is the admissibility
82 : * condition (the returned stress lies on the yield surface) and g
83 : * is the flow potential so the other conditions implement the normality
84 : * condition. It is up to the derived classes to defined E[i, j]
85 : * so that the rhs really represents the normality condition.
86 : *
87 : * For instance, CappedDruckerPrager has (with Eijkl being the elasticity
88 : * tensor):
89 : * E[0, 0] = _Epp = Eijkl.sum3x3()
90 : * E[0, 1] = 0
91 : * E[1, 0] = 0
92 : * E[1, 1] = Eijkl(0, 1, 0, 1)
93 : */
94 : class MultiParameterPlasticityStressUpdate : public StressUpdateBase
95 : {
96 : public:
97 : static InputParameters validParams();
98 :
99 : MultiParameterPlasticityStressUpdate(const InputParameters & parameters,
100 : unsigned num_sp,
101 : unsigned num_yf,
102 : unsigned num_intnl);
103 :
104 : protected:
105 : virtual void initQpStatefulProperties() override;
106 : using StressUpdateBase::updateState;
107 : virtual void updateState(RankTwoTensor & strain_increment,
108 : RankTwoTensor & inelastic_strain_increment,
109 : const RankTwoTensor & rotation_increment,
110 : RankTwoTensor & stress_new,
111 : const RankTwoTensor & stress_old,
112 : const RankFourTensor & elasticity_tensor,
113 : const RankTwoTensor & elastic_strain_old,
114 : bool compute_full_tangent_operator,
115 : RankFourTensor & tangent_operator) override;
116 :
117 : virtual void propagateQpStatefulProperties() override;
118 :
119 2568 : virtual TangentCalculationMethod getTangentCalculationMethod() override
120 : {
121 2568 : return TangentCalculationMethod::FULL;
122 : }
123 :
124 : /// Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics)
125 : constexpr static unsigned _tensor_dimensionality = 3;
126 :
127 : /// Number of stress parameters
128 : const unsigned _num_sp;
129 :
130 : /// An admissible value of stress_params at the initial time
131 : const std::vector<Real> _definitely_ok_sp;
132 :
133 : /// E[i, j] in the system of equations to be solved
134 : std::vector<std::vector<Real>> _Eij;
135 :
136 : /// normalising factor
137 : Real _En;
138 :
139 : /// _Cij[i, j] * _Eij[j, k] = 1 iff j == k
140 : std::vector<std::vector<Real>> _Cij;
141 :
142 : /// Number of yield functions
143 : const unsigned _num_yf;
144 :
145 : /// Number of internal parameters
146 : const unsigned _num_intnl;
147 :
148 : /// Maximum number of Newton-Raphson iterations allowed in the return-map process
149 : const unsigned _max_nr_its;
150 :
151 : /// Whether to perform finite-strain rotations
152 : const bool _perform_finite_strain_rotations;
153 :
154 : /// Smoothing tolerance: edges of the yield surface get smoothed by this amount
155 : const Real _smoothing_tol;
156 :
157 : /// Square of the smoothing tolerance
158 : const Real _smoothing_tol2;
159 :
160 : /// The yield-function tolerance
161 : const Real _f_tol;
162 :
163 : /// Square of the yield-function tolerance
164 : const Real _f_tol2;
165 :
166 : /**
167 : * In order to help the Newton-Raphson procedure, the applied
168 : * strain increment may be applied in sub-increments of size
169 : * greater than this value.
170 : */
171 : const Real _min_step_size;
172 :
173 : /// handles case of initial_stress that is inadmissible being supplied
174 : bool _step_one;
175 :
176 : /// Output a warning message if precision loss is encountered during the return-map process
177 : const bool _warn_about_precision_loss;
178 :
179 : /// plastic strain
180 : MaterialProperty<RankTwoTensor> & _plastic_strain;
181 :
182 : /// Old value of plastic strain
183 : const MaterialProperty<RankTwoTensor> & _plastic_strain_old;
184 :
185 : /// internal parameters
186 : MaterialProperty<std::vector<Real>> & _intnl;
187 :
188 : /// old values of internal parameters
189 : const MaterialProperty<std::vector<Real>> & _intnl_old;
190 :
191 : /// yield functions
192 : MaterialProperty<std::vector<Real>> & _yf;
193 :
194 : /// Number of Newton-Raphson iterations used in the return-map
195 : MaterialProperty<Real> & _iter;
196 :
197 : /// Maximum number of Newton-Raphson iterations used in the return-map during the course of the entire simulation
198 : MaterialProperty<Real> & _max_iter_used;
199 :
200 : /// Old value of maximum number of Newton-Raphson iterations used in the return-map during the course of the entire simulation
201 : const MaterialProperty<Real> & _max_iter_used_old;
202 :
203 : /// Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise)
204 : MaterialProperty<Real> & _linesearch_needed;
205 :
206 : /**
207 : * Struct designed to hold info about a single yield function
208 : * and its derivatives, as well as the flow directions
209 : */
210 : struct yieldAndFlow
211 : {
212 : Real f; // yield function value
213 : std::vector<Real> df; // df/d(stress_param[i])
214 : std::vector<Real> df_di; // df/d(intnl_variable[a])
215 : std::vector<Real> dg; // d(flow)/d(stress_param[i])
216 : std::vector<std::vector<Real>> d2g; // d^2(flow)/d(sp[i])/d(sp[j])
217 : std::vector<std::vector<Real>> d2g_di; // d^2(flow)/d(sp[i])/d(intnl[a])
218 :
219 625652 : yieldAndFlow() : yieldAndFlow(0, 0) {}
220 :
221 3135458 : yieldAndFlow(unsigned num_var, unsigned num_intnl)
222 3135458 : : f(0.0),
223 3135458 : df(num_var),
224 3135458 : df_di(num_intnl),
225 3135458 : dg(num_var),
226 3135458 : d2g(num_var, std::vector<Real>(num_var, 0.0)),
227 3135458 : d2g_di(num_var, std::vector<Real>(num_intnl, 0.0))
228 : {
229 3135458 : }
230 :
231 : // this may be involved in the smoothing of a group of yield functions
232 : bool operator<(const yieldAndFlow & fd) const { return f < fd.f; }
233 : };
234 :
235 : /**
236 : * Computes the smoothed yield function
237 : * @param stress_params The stress parameters (eg stress_params[0] = stress_zz and
238 : * stress_params[1] = sqrt(stress_zx^2 + stress_zy^2))
239 : * @param intnl The internal parameters (eg intnl[0] is shear, intnl[1] is tensile)
240 : * @return The smoothed yield function value
241 : */
242 : Real yieldF(const std::vector<Real> & stress_params, const std::vector<Real> & intnl) const;
243 :
244 : /**
245 : * Computes the smoothed yield function
246 : * @param yfs The values of the individual yield functions
247 : * @return The smoothed yield function value
248 : */
249 : Real yieldF(const std::vector<Real> & yfs) const;
250 :
251 : /**
252 : * Smooths yield functions. The returned value must be zero if abs(f_diff) >= _smoothing_tol
253 : * and otherwise must satisfy, over -_smoothing_tol <= f_diff <= _smoothing_tol:
254 : * (1) C2
255 : * (2) zero at f_diff = +/- _smoothing_tol
256 : * (3) derivative is +/-0.5 at f_diff = +/- _smoothing_tol
257 : * (4) derivative must be in [-0.5, 0.5]
258 : * (5) second derivative is zero at f_diff = +/- _smoothing_tol
259 : * (6) second derivative must be non-negative
260 : * in order to ensure C2 differentiability and convexity of the smoothed yield surface.
261 : */
262 : Real ismoother(Real f_diff) const;
263 :
264 : /**
265 : * Derivative of ismoother
266 : */
267 : Real smoother(Real f_diff) const;
268 :
269 : /**
270 : * Derivative of smoother
271 : */
272 : Real dsmoother(Real f_diff) const;
273 :
274 : /**
275 : * Calculates all yield functions and derivatives, and then performs the smoothing scheme
276 : * @param stress_params[in] The stress parameters (eg stress_params[0] = stress_zz and
277 : * stress_params[1] = sqrt(stress_zx^2 + stress_zy^2))
278 : * @param intnl[in] Internal parameters
279 : * @return The smoothed yield function and derivatives
280 : */
281 : yieldAndFlow smoothAllQuantities(const std::vector<Real> & stress_params,
282 : const std::vector<Real> & intnl) const;
283 :
284 : /**
285 : * Performs a line-search to find stress_params and gaE
286 : * Upon entry:
287 : * - rhs contains the *solution* to the Newton-Raphson (ie nrStep should have been called). If a
288 : * full
289 : * Newton step is used then stress_params[:] += rhs[0:_num_sp-1] and gaE += rhs[_num_sp]
290 : * - res2 contains the residual-squared before applying any of solution
291 : * - stress_params contains the stress_params before applying any of the solution
292 : * - gaE contains gaE before applying any of the solution (that is contained in rhs)
293 : * Upon exit:
294 : * - stress_params will be the stress_params after applying the solution
295 : * - gaE will be the stress_params after applying the solution
296 : * - rhs will contain the updated rhs values (after applying the solution) ready for the next
297 : * Newton-Raphson step,
298 : * - res2 will be the residual-squared after applying the solution
299 : * - intnl will contain the internal variables corresponding to the return from
300 : * trial_stress_params to stress_params
301 : * (and starting from intnl_ok)
302 : * - linesearch_needed will be 1.0 if a linesearch was needed
303 : * - smoothed_q will contain the value of the yield function and its derivatives, etc, at
304 : * (stress_params, intnl)
305 : * @param res2[in,out] the residual-squared, both as an input and output
306 : * @param stress_params[in,out] Upon input the value of the stress_params before the current
307 : * Newton-Raphson process was initiated.
308 : * Upon exit this will hold the values coming from the line search.
309 : * @param trial_stress_params[in] Trial value for the stress_params for this (sub)strain increment
310 : * @param gaE[in,out] Upon input the value of gaE before the current Newton-Raphson iteration was
311 : * initiated. Upon exit this will hold
312 : * the value coming from the line-search
313 : * @param smoothed_q[in,out] Upon input, the value of the smoothed yield function and derivatives
314 : * at the
315 : * prior-to-Newton configuration. Upon exit this is evaluated at the new (stress_params, intnl)
316 : * @param intnl_ok[in] The value of the internal parameters from the start of this (sub)strain
317 : * increment
318 : * @param intnl[in,out] The value of the internal parameters after the line-search has converged
319 : * @param rhs[in,out] Upon entry this contains the solution to the Newton-Raphson. Upon exit this
320 : * contains the updated rhs values
321 : * @return 0 if successful, 1 otherwise
322 : */
323 : int lineSearch(Real & res2,
324 : std::vector<Real> & stress_params,
325 : Real & gaE,
326 : const std::vector<Real> & trial_stress_params,
327 : yieldAndFlow & smoothed_q,
328 : const std::vector<Real> & intnl_ok,
329 : std::vector<Real> & intnl,
330 : std::vector<Real> & rhs,
331 : Real & linesearch_needed) const;
332 :
333 : /**
334 : * Performs a Newton-Raphson step to attempt to zero rhs
335 : * Upon return, rhs will contain the solution.
336 : * @param smoothed_q[in] The value of the smoothed yield function and derivatives prior to this
337 : * Newton-Raphson step
338 : * @param trial_stress_params[in] Trial value for the stress_params for this (sub)strain increment
339 : * @param stress_params[in] The current value of the stress_params
340 : * @param intnl[in] The current value of the internal parameters
341 : * @param gaE[in] The current value of gaE
342 : * @param rhs[in,out] Upon entry, the rhs to zero using Newton-Raphson. Upon exit, the solution
343 : * to the Newton-Raphson problem
344 : * @return 0 if successful, 1 otherwise
345 : */
346 : int nrStep(const yieldAndFlow & smoothed_q,
347 : const std::vector<Real> & trial_stress_params,
348 : const std::vector<Real> & stress_params,
349 : const std::vector<Real> & intnl,
350 : Real gaE,
351 : std::vector<Real> & rhs) const;
352 :
353 : /**
354 : * Calculates the residual-squared for the Newton-Raphson + line-search
355 : * @param rhs[in] The RHS vector
356 : * @return sum_i (rhs[i] * rhs[i])
357 : */
358 : Real calculateRes2(const std::vector<Real> & rhs) const;
359 :
360 : /**
361 : * Calculates the RHS in the following
362 : * 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j]
363 : * 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, j] * dg/dS[j]
364 : * ...
365 : * 0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, j] * dg/dS[j]
366 : * 0 = rhs[N] = f(S, intnl)
367 : * where N = _num_sp
368 : * @param trial_stress_params[in] The trial stress parameters for this (sub)strain increment,
369 : * S[:]^trial
370 : * @param stress_params[in] The current stress parameters, S[:]
371 : * @param gaE[in] ga*_En (the normalisation with _En is so that gaE is of similar magnitude to S)
372 : * @param smoothed_q[in] Holds the current value of yield function and derivatives evaluated at
373 : * the current stress parameters and the current internal parameters
374 : * @param rhs[out] The result
375 : */
376 : void calculateRHS(const std::vector<Real> & trial_stress_params,
377 : const std::vector<Real> & stress_params,
378 : Real gaE,
379 : const yieldAndFlow & smoothed_q,
380 : std::vector<Real> & rhs) const;
381 :
382 : /**
383 : * Derivative of -RHS with respect to the stress_params and gaE, placed
384 : * into an array ready for solving the linear system using
385 : * LAPACK gsev
386 : * @param smoothed_q[in] Holds the current value of yield function and derivatives evaluated at
387 : * the
388 : * current values of the stress_params and the internal parameters
389 : * @param dintnl[in] The derivatives of the internal parameters wrt the stress_params
390 : * @param stress_params[in] The current value of the stress_params during the Newton-Raphson
391 : * process
392 : * @param gaE[in] The current value of gaE
393 : * @param jac[out] The outputted derivatives
394 : */
395 : void dnRHSdVar(const yieldAndFlow & smoothed_q,
396 : const std::vector<std::vector<Real>> & dintnl,
397 : const std::vector<Real> & stress_params,
398 : Real gaE,
399 : std::vector<double> & jac) const;
400 :
401 : /**
402 : * Performs any necessary cleaning-up, then throw MooseException(message)
403 : * @param message The message to using in MooseException
404 : */
405 : virtual void errorHandler(const std::string & message) const;
406 :
407 : /**
408 : * Computes the values of the yield functions, given stress_params and intnl parameters.
409 : * Derived classes must override this, to provide the values of the yield functions
410 : * in yf.
411 : * @param stress_params[in] The stress parameters
412 : * @param intnl[in] The internal parameters
413 : * @param[out] yf The yield function values
414 : */
415 : virtual void yieldFunctionValuesV(const std::vector<Real> & stress_params,
416 : const std::vector<Real> & intnl,
417 : std::vector<Real> & yf) const = 0;
418 :
419 : /**
420 : * Completely fills all_q with correct values. These values are:
421 : * (1) the yield function values, yf[i]
422 : * (2) d(yf[i])/d(stress_params[j])
423 : * (3) d(yf[i])/d(intnl[j])
424 : * (4) d(flowPotential[i])/d(stress_params[j])
425 : * (5) d2(flowPotential[i])/d(stress_params[j])/d(stress_params[k])
426 : * (6) d2(flowPotential[i])/d(stress_params[j])/d(intnl[k])
427 : * @param stress_params[in] The stress parameters
428 : * @param intnl[in] The internal parameters
429 : * @param[out] all_q All the desired quantities
430 : */
431 : virtual void computeAllQV(const std::vector<Real> & stress_params,
432 : const std::vector<Real> & intnl,
433 : std::vector<yieldAndFlow> & all_q) const = 0;
434 :
435 : /**
436 : * Derived classes may employ this function to record stuff or do
437 : * other computations prior to the return-mapping algorithm. We
438 : * know that (trial_stress_params, intnl_old) is inadmissible when
439 : * this is called
440 : * @param trial_stress_params[in] The trial values of the stress parameters
441 : * @param stress_trial[in] Trial stress tensor
442 : * @param intnl_old[in] Old value of the internal parameters.
443 : * @param yf[in] The yield functions at (p_trial, q_trial, intnl_old)
444 : * @param Eijkl[in] The elasticity tensor
445 : */
446 : virtual void preReturnMapV(const std::vector<Real> & trial_stress_params,
447 : const RankTwoTensor & stress_trial,
448 : const std::vector<Real> & intnl_old,
449 : const std::vector<Real> & yf,
450 : const RankFourTensor & Eijkl);
451 :
452 : /**
453 : * Sets (stress_params, intnl) at "good guesses" of the solution to the Return-Map algorithm.
454 : * The purpose of these "good guesses" is to speed up the Newton-Raphson process by providing
455 : * it with a good initial guess.
456 : * Derived classes may choose to override this if their plastic models are easy
457 : * enough to solve approximately.
458 : * The default values, provided by this class, are simply gaE = 0, stress_params =
459 : * trial_stress_params,
460 : * that is, the "good guess" is just the trial point for this (sub)strain increment.
461 : * @param trial_stress_params[in] The stress_params at the trial point
462 : * @param intnl_old[in] The internal parameters before applying the (sub)strain increment
463 : * @param stress_params[out] The "good guess" value of the stress_params
464 : * @param gaE[out] The "good guess" value of gaE
465 : * @param intnl[out] The "good guess" value of the internal parameters
466 : */
467 : virtual void initializeVarsV(const std::vector<Real> & trial_stress_params,
468 : const std::vector<Real> & intnl_old,
469 : std::vector<Real> & stress_params,
470 : Real & gaE,
471 : std::vector<Real> & intnl) const;
472 :
473 : /**
474 : * Sets the internal parameters based on the trial values of
475 : * stress_params, their current values, and the old values of the
476 : * internal parameters.
477 : * Derived classes must override this.
478 : * @param trial_stress_params[in] The trial stress parameters (eg trial_p and trial_q)
479 : * @param current_stress_params[in] The current stress parameters (eg p and q)
480 : * @param intnl_old[out] Old value of internal parameters
481 : * @param intnl[out] The value of internal parameters to be set
482 : */
483 : virtual void setIntnlValuesV(const std::vector<Real> & trial_stress_params,
484 : const std::vector<Real> & current_stress_params,
485 : const std::vector<Real> & intnl_old,
486 : std::vector<Real> & intnl) const = 0;
487 :
488 : /**
489 : * Sets the derivatives of internal parameters, based on the trial values of
490 : * stress_params, their current values, and the current values of the
491 : * internal parameters.
492 : * Derived classes must override this.
493 : * @param trial_stress_params[in] The trial stress parameters
494 : * @param current_stress_params[in] The current stress parameters
495 : * @param intnl[in] The current value of the internal parameters
496 : * @param dintnl[out] The derivatives dintnl[i][j] = d(intnl[i])/d(stress_param j)
497 : */
498 : virtual void setIntnlDerivativesV(const std::vector<Real> & trial_stress_params,
499 : const std::vector<Real> & current_stress_params,
500 : const std::vector<Real> & intnl,
501 : std::vector<std::vector<Real>> & dintnl) const = 0;
502 :
503 : /**
504 : * Computes stress_params, given stress. Derived classes must
505 : * override this
506 : * @param stress[in] Stress tensor
507 : * @param stress_params[out] The compute stress_params
508 : */
509 : virtual void computeStressParams(const RankTwoTensor & stress,
510 : std::vector<Real> & stress_params) const = 0;
511 :
512 : /**
513 : * Derived classes may use this to perform calculations before
514 : * any return-map process is performed, for instance, to initialize
515 : * variables.
516 : * This is called at the very start of updateState, even before
517 : * any checking for admissible stresses, etc, is performed
518 : */
519 : virtual void initializeReturnProcess();
520 :
521 : /**
522 : * Derived classes may use this to perform calculations after the
523 : * return-map process has completed successfully in stress_param space
524 : * but before the returned stress tensor has been calculcated.
525 : * @param rotation_increment[in] The large-strain rotation increment
526 : */
527 : virtual void finalizeReturnProcess(const RankTwoTensor & rotation_increment);
528 :
529 : /**
530 : * Sets stress from the admissible parameters.
531 : * This is called after the return-map process has completed
532 : * successfully in stress_param space, just after finalizeReturnProcess
533 : * has been called.
534 : * Derived classes must override this function
535 : * @param stress_trial[in] The trial value of stress
536 : * @param stress_params[in] The value of the stress_params after the return-map process has
537 : * completed successfully
538 : * @param gaE[in] The value of gaE after the return-map process has completed successfully
539 : * @param intnl[in] The value of the internal parameters after the return-map process has
540 : * completed successfully
541 : * @param smoothed_q[in] Holds the current value of yield function and derivatives evaluated at
542 : * the returned state
543 : * @param Eijkl[in] The elasticity tensor
544 : * @param stress[out] The returned value of the stress tensor
545 : */
546 : virtual void setStressAfterReturnV(const RankTwoTensor & stress_trial,
547 : const std::vector<Real> & stress_params,
548 : Real gaE,
549 : const std::vector<Real> & intnl,
550 : const yieldAndFlow & smoothed_q,
551 : const RankFourTensor & Eijkl,
552 : RankTwoTensor & stress) const = 0;
553 :
554 : /**
555 : * Sets inelastic strain increment from the returned configuration
556 : * This is called after the return-map process has completed
557 : * successfully in stress_param space, just after finalizeReturnProcess
558 : * has been called.
559 : * Derived classes may override this function
560 : * @param stress_trial[in] The trial value of stress
561 : * @param gaE[in] The value of gaE after the return-map process has completed successfully
562 : * @param smoothed_q[in] Holds the current value of yield function and derivatives evaluated at
563 : * the returned configuration
564 : * @param elasticity_tensor[in] The elasticity tensor
565 : * @param returned_stress[in] The stress after the return-map process
566 : * @param inelastic_strain_increment[out] The inelastic strain increment resulting from this
567 : * return-map
568 : */
569 : virtual void
570 : setInelasticStrainIncrementAfterReturn(const RankTwoTensor & stress_trial,
571 : Real gaE,
572 : const yieldAndFlow & smoothed_q,
573 : const RankFourTensor & elasticity_tensor,
574 : const RankTwoTensor & returned_stress,
575 : RankTwoTensor & inelastic_strain_increment) const;
576 :
577 : /**
578 : * Calculates the consistent tangent operator.
579 : * Derived classes may choose to override this for computational
580 : * efficiency. The implementation in this class is quite expensive,
581 : * even though it looks compact and clean, because of all the
582 : * manipulations of RankFourTensors involved.
583 : * @param stress_trial[in] the trial value of the stress tensor for this strain increment
584 : * @param trial_stress_params[in] the trial values of the stress_params for this strain increment
585 : * @param stress[in] the returned value of the stress tensor for this strain increment
586 : * @param stress_params[in] the returned value of the stress_params for this strain increment
587 : * @param gaE[in] the total value of that came from this strain increment
588 : * @param smoothed_q[in] contains the yield function and derivatives evaluated at (p, q)
589 : * @param Eijkl[in] The elasticity tensor
590 : * @param compute_full_tangent_operator[in] true if the full consistent tangent operator is
591 : * needed,
592 : * otherwise false
593 : * @param dvar_dtrial[in] dvar_dtrial[i][j] = d({stress_param[i],gaE})/d(trial_stress_param[j])
594 : * for
595 : * this strain increment
596 : * @param[out] cto The consistent tangent operator
597 : */
598 : virtual void consistentTangentOperatorV(const RankTwoTensor & stress_trial,
599 : const std::vector<Real> & trial_stress_params,
600 : const RankTwoTensor & stress,
601 : const std::vector<Real> & stress_params,
602 : Real gaE,
603 : const yieldAndFlow & smoothed_q,
604 : const RankFourTensor & Eijkl,
605 : bool compute_full_tangent_operator,
606 : const std::vector<std::vector<Real>> & dvar_dtrial,
607 : RankFourTensor & cto);
608 :
609 : /**
610 : * d(stress_param[i])/d(stress) at given stress
611 : * @param stress stress tensor
612 : * @return d(stress_param[:])/d(stress)
613 : */
614 : virtual std::vector<RankTwoTensor> dstress_param_dstress(const RankTwoTensor & stress) const = 0;
615 :
616 : /**
617 : * d2(stress_param[i])/d(stress)/d(stress) at given stress
618 : * @param stress stress tensor
619 : * @return d2(stress_param[:])/d(stress)/d(stress)
620 : */
621 : virtual std::vector<RankFourTensor>
622 : d2stress_param_dstress(const RankTwoTensor & stress) const = 0;
623 :
624 : /// Sets _Eij and _En and _Cij
625 : virtual void setEffectiveElasticity(const RankFourTensor & Eijkl) = 0;
626 :
627 : /**
628 : * Calculates derivatives of the stress_params and gaE with repect to the
629 : * trial values of the stress_params for the (sub)strain increment.
630 : * After the strain increment has been fully applied, dvar_dtrial will
631 : * contain the result appropriate to the full strain increment. Before
632 : * that time (if applying in sub-strain increments) it will contain the
633 : * result appropriate to the amount of strain increment applied successfully.
634 : * @param elastic_only[in] whether this was an elastic step: if so then the updates to dvar_dtrial
635 : * are fairly trivial
636 : * @param trial_stress_params[in] Trial values of stress_params for this (sub)strain increment
637 : * @param stress_params[in] Returned values of stress_params for this (sub)strain increment
638 : * @param gaE[in] the value of gaE that came from this (sub)strain increment
639 : * @param intnl[in] the value of the internal parameters at the returned position
640 : * @param smoothed_q[in] contains the yield function and derivatives evaluated at (stress_params,
641 : * intnl)
642 : * @param step_size[in] size of this (sub)strain increment
643 : * @param compute_full_tangent_operator[in] true if the full consistent tangent operator is
644 : * needed,
645 : * otherwise false
646 : * @param dvar_dtrial[out] dvar_dtrial[i][j] = d({stress_param[i],gaE})/d(trial_stress_param[j])
647 : */
648 : void dVardTrial(bool elastic_only,
649 : const std::vector<Real> & trial_stress_params,
650 : const std::vector<Real> & stress_params,
651 : Real gaE,
652 : const std::vector<Real> & intnl,
653 : const yieldAndFlow & smoothed_q,
654 : Real step_size,
655 : bool compute_full_tangent_operator,
656 : std::vector<std::vector<Real>> & dvar_dtrial) const;
657 :
658 : /**
659 : * Check whether precision loss has occurred
660 : * @param[in] solution The solution to the Newton-Raphson system
661 : * @param[in] stress_params The currect values of the stress_params for this (sub)strain increment
662 : * @param[in] gaE The currenct value of gaE for this (sub)strain increment
663 : * @return true if precision loss has occurred
664 : */
665 : bool precisionLoss(const std::vector<Real> & solution,
666 : const std::vector<Real> & stress_params,
667 : Real gaE) const;
668 :
669 : private:
670 : /**
671 : * "Trial" value of stress_params that initializes the return-map process
672 : * This is derived from stress = stress_old + Eijkl * strain_increment.
673 : * However, since the return-map process can fail and be restarted by
674 : * applying strain_increment in multiple substeps, _trial_sp can vary
675 : * from substep to substep.
676 : */
677 : std::vector<Real> _trial_sp;
678 :
679 : /**
680 : * "Trial" value of stress that is set at the beginning of the
681 : * return-map process. It is fixed at stress_old + Eijkl * strain_increment
682 : * irrespective of any sub-stepping
683 : */
684 : RankTwoTensor _stress_trial;
685 :
686 : /**
687 : * 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, i] * dg/dS[i]
688 : * 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, i] * dg/dS[i]
689 : * ...
690 : * 0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, i] * dg/dS[i]
691 : * 0 = rhs[N] = f(S, intnl)
692 : * Here N = num_sp
693 : */
694 : std::vector<Real> _rhs;
695 :
696 : /**
697 : * d({stress_param[i], gaE})/d(trial_stress_param[j])
698 : */
699 : std::vector<std::vector<Real>> _dvar_dtrial;
700 :
701 : /**
702 : * The state (ok_sp, ok_intnl) is known to be admissible, so
703 : * ok_sp are stress_params that are "OK". If the strain_increment
704 : * is applied in substeps then ok_sp is updated after each
705 : * sub strain_increment is applied and the return-map is successful.
706 : * At the end of the entire return-map process _ok_sp will contain
707 : * the stress_params where (_ok_sp, _intnl) is admissible.
708 : */
709 : std::vector<Real> _ok_sp;
710 :
711 : /**
712 : * The state (ok_sp, ok_intnl) is known to be admissible
713 : */
714 : std::vector<Real> _ok_intnl;
715 :
716 : /**
717 : * _del_stress_params = trial_stress_params - ok_sp
718 : * This is fixed at the beginning of the return-map process,
719 : * irrespective of substepping. The return-map problem is:
720 : * apply del_stress_params to stress_prams, and then find
721 : * an admissible (returned) stress_params and gaE
722 : */
723 : std::vector<Real> _del_stress_params;
724 :
725 : /**
726 : * The current values of the stress params during the Newton-Raphson
727 : */
728 : std::vector<Real> _current_sp;
729 :
730 : /**
731 : * The current values of the internal params during the Newton-Raphson
732 : */
733 : std::vector<Real> _current_intnl;
734 :
735 : private:
736 : /**
737 : * The type of smoother function
738 : */
739 : enum class SmootherFunctionType
740 : {
741 : cos,
742 : poly1,
743 : poly2,
744 : poly3
745 : } _smoother_function_type;
746 : };
|
