www.mooseframework.org
Public Member Functions | Protected Member Functions | Protected Attributes | Static Protected Attributes | List of all members
CappedMohrCoulombCosseratStressUpdate Class Reference

CappedMohrCoulombCosseratStressUpdate implements rate-independent nonassociative Mohr-Coulomb plus tensile plus compressive plasticity with hardening/softening in the Cosserat setting. More...

#include <CappedMohrCoulombCosseratStressUpdate.h>

Inheritance diagram for CappedMohrCoulombCosseratStressUpdate:
[legend]

Public Member Functions

 CappedMohrCoulombCosseratStressUpdate (const InputParameters &parameters)
 
bool requiresIsotropicTensor () override
 The full elasticity tensor may be anisotropic, and usually is in the case of layered Cosserat. More...
 
void setQp (unsigned int qp)
 Sets the value of the global variable _qp for inheriting classes. More...
 
virtual Real computeTimeStepLimit ()
 
void resetQpProperties () final
 Retained as empty methods to avoid a warning from Material.C in framework. These methods are unused in all inheriting classes and should not be overwritten. More...
 
void resetProperties () final
 

Protected Member Functions

virtual void setStressAfterReturnV (const RankTwoTensor &stress_trial, const std::vector< Real > &stress_params, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, RankTwoTensor &stress) const override
 Sets stress from the admissible parameters. More...
 
virtual void preReturnMapV (const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress_trial, const std::vector< Real > &intnl_old, const std::vector< Real > &yf, const RankFourTensor &Eijkl) override
 Derived classes may employ this function to record stuff or do other computations prior to the return-mapping algorithm. More...
 
void setEffectiveElasticity (const RankFourTensor &Eijkl) override
 Sets _Eij and _En and _Cij. More...
 
virtual void consistentTangentOperatorV (const RankTwoTensor &stress_trial, const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, const std::vector< std::vector< Real >> &dvar_dtrial, RankFourTensor &cto) override
 Calculates the consistent tangent operator. More...
 
void computeStressParams (const RankTwoTensor &stress, std::vector< Real > &stress_params) const override
 Computes stress_params, given stress. More...
 
std::vector< RankTwoTensor > dstress_param_dstress (const RankTwoTensor &stress) const override
 d(stress_param[i])/d(stress) at given stress More...
 
std::vector< RankFourTensor > d2stress_param_dstress (const RankTwoTensor &stress) const override
 d2(stress_param[i])/d(stress)/d(stress) at given stress More...
 
void yieldFunctionValuesV (const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< Real > &yf) const override
 Computes the values of the yield functions, given stress_params and intnl parameters. More...
 
void computeAllQV (const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< yieldAndFlow > &all_q) const override
 Completely fills all_q with correct values. More...
 
void initializeVarsV (const std::vector< Real > &trial_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &stress_params, Real &gaE, std::vector< Real > &intnl) const override
 Sets (stress_params, intnl) at "good guesses" of the solution to the Return-Map algorithm. More...
 
void setIntnlValuesV (const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &intnl) const override
 Sets the internal parameters based on the trial values of stress_params, their current values, and the old values of the internal parameters. More...
 
void setIntnlDerivativesV (const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl, std::vector< std::vector< Real >> &dintnl) const override
 Sets the derivatives of internal parameters, based on the trial values of stress_params, their current values, and the current values of the internal parameters. More...
 
virtual void initQpStatefulProperties () override
 
virtual void updateState (RankTwoTensor &strain_increment, RankTwoTensor &inelastic_strain_increment, const RankTwoTensor &rotation_increment, RankTwoTensor &stress_new, const RankTwoTensor &stress_old, const RankFourTensor &elasticity_tensor, const RankTwoTensor &elastic_strain_old, bool compute_full_tangent_operator, RankFourTensor &tangent_operator) override
 Given a strain increment that results in a trial stress, perform some procedure (such as an iterative return-mapping process) to produce an admissible stress, an elastic strain increment and an inelastic strain increment. More...
 
virtual void propagateQpStatefulProperties () override
 If updateState is not called during a timestep, this will be. More...
 
virtual TangentCalculationMethod getTangentCalculationMethod () override
 
Real yieldF (const std::vector< Real > &stress_params, const std::vector< Real > &intnl) const
 Computes the smoothed yield function. More...
 
Real yieldF (const std::vector< Real > &yfs) const
 Computes the smoothed yield function. More...
 
Real ismoother (Real f_diff) const
 Smooths yield functions. More...
 
Real smoother (Real f_diff) const
 Derivative of ismoother. More...
 
Real dsmoother (Real f_diff) const
 Derivative of smoother. More...
 
yieldAndFlow smoothAllQuantities (const std::vector< Real > &stress_params, const std::vector< Real > &intnl) const
 Calculates all yield functions and derivatives, and then performs the smoothing scheme. More...
 
int lineSearch (Real &res2, std::vector< Real > &stress_params, Real &gaE, const std::vector< Real > &trial_stress_params, yieldAndFlow &smoothed_q, const std::vector< Real > &intnl_ok, std::vector< Real > &intnl, std::vector< Real > &rhs, Real &linesearch_needed) const
 Performs a line-search to find stress_params and gaE Upon entry: More...
 
int nrStep (const yieldAndFlow &smoothed_q, const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, const std::vector< Real > &intnl, Real gaE, std::vector< Real > &rhs) const
 Performs a Newton-Raphson step to attempt to zero rhs Upon return, rhs will contain the solution. More...
 
Real calculateRes2 (const std::vector< Real > &rhs) const
 Calculates the residual-squared for the Newton-Raphson + line-search. More...
 
void calculateRHS (const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, std::vector< Real > &rhs) const
 Calculates the RHS in the following 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j] 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, j] * dg/dS[j] ... More...
 
void dnRHSdVar (const yieldAndFlow &smoothed_q, const std::vector< std::vector< Real >> &dintnl, const std::vector< Real > &stress_params, Real gaE, std::vector< double > &jac) const
 Derivative of -RHS with respect to the stress_params and gaE, placed into an array ready for solving the linear system using LAPACK gsev. More...
 
virtual void errorHandler (const std::string &message) const
 Performs any necessary cleaning-up, then throw MooseException(message) More...
 
virtual void initializeReturnProcess ()
 Derived classes may use this to perform calculations before any return-map process is performed, for instance, to initialize variables. More...
 
virtual void finalizeReturnProcess (const RankTwoTensor &rotation_increment)
 Derived classes may use this to perform calculations after the return-map process has completed successfully in stress_param space but before the returned stress tensor has been calculcated. More...
 
virtual void setInelasticStrainIncrementAfterReturn (const RankTwoTensor &stress_trial, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &elasticity_tensor, const RankTwoTensor &returned_stress, RankTwoTensor &inelastic_strain_increment) const
 Sets inelastic strain increment from the returned configuration This is called after the return-map process has completed successfully in stress_param space, just after finalizeReturnProcess has been called. More...
 
void dVardTrial (bool elastic_only, const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, Real step_size, bool compute_full_tangent_operator, std::vector< std::vector< Real >> &dvar_dtrial) const
 Calculates derivatives of the stress_params and gaE with repect to the trial values of the stress_params for the (sub)strain increment. More...
 
bool precisionLoss (const std::vector< Real > &solution, const std::vector< Real > &stress_params, Real gaE) const
 Check whether precision loss has occurred. More...
 

Protected Attributes

const Real _host_young
 Young's modulus of the host material. More...
 
const Real _host_poisson
 Poisson's of the host material. More...
 
const Real _host_E0011
 E0011 = Lame lambda modulus of the host material. More...
 
const Real _host_E0000
 E0000 = Lame lambda + 2 * shear modulus of the host material. More...
 
const TensorMechanicsHardeningModel_tensile_strength
 Hardening model for tensile strength. More...
 
const TensorMechanicsHardeningModel_compressive_strength
 Hardening model for compressive strength. More...
 
const TensorMechanicsHardeningModel_cohesion
 Hardening model for cohesion. More...
 
const TensorMechanicsHardeningModel_phi
 Hardening model for friction angle. More...
 
const TensorMechanicsHardeningModel_psi
 Hardening model for dilation angle. More...
 
const bool _perfect_guess
 Whether to provide an estimate of the returned stress, based on perfect plasticity. More...
 
Real _poissons_ratio
 Poisson's ratio. More...
 
const Real _shifter
 When equal-eigenvalues are predicted from the stress initialization routine, shift them by this amount. More...
 
RankTwoTensor _eigvecs
 Eigenvectors of the trial stress as a RankTwoTensor, in order to rotate the returned stress back to stress space. More...
 
const unsigned _num_sp
 Number of stress parameters. More...
 
const std::vector< Real > _definitely_ok_sp
 An admissible value of stress_params at the initial time. More...
 
std::vector< std::vector< Real > > _Eij
 E[i, j] in the system of equations to be solved. More...
 
Real _En
 normalising factor More...
 
std::vector< std::vector< Real > > _Cij
 _Cij[i, j] * _Eij[j, k] = 1 iff j == k More...
 
const unsigned _num_yf
 Number of yield functions. More...
 
const unsigned _num_intnl
 Number of internal parameters. More...
 
const unsigned _max_nr_its
 Maximum number of Newton-Raphson iterations allowed in the return-map process. More...
 
const bool _perform_finite_strain_rotations
 Whether to perform finite-strain rotations. More...
 
const Real _smoothing_tol
 Smoothing tolerance: edges of the yield surface get smoothed by this amount. More...
 
const Real _f_tol
 The yield-function tolerance. More...
 
const Real _f_tol2
 Square of the yield-function tolerance. More...
 
const Real _min_step_size
 In order to help the Newton-Raphson procedure, the applied strain increment may be applied in sub-increments of size greater than this value. More...
 
bool _step_one
 handles case of initial_stress that is inadmissible being supplied More...
 
const bool _warn_about_precision_loss
 Output a warning message if precision loss is encountered during the return-map process. More...
 
MaterialProperty< RankTwoTensor > & _plastic_strain
 plastic strain More...
 
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
 Old value of plastic strain. More...
 
MaterialProperty< std::vector< Real > > & _intnl
 internal parameters More...
 
const MaterialProperty< std::vector< Real > > & _intnl_old
 old values of internal parameters More...
 
MaterialProperty< std::vector< Real > > & _yf
 yield functions More...
 
MaterialProperty< Real > & _iter
 Number of Newton-Raphson iterations used in the return-map. More...
 
MaterialProperty< Real > & _linesearch_needed
 Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) More...
 
const std::string _base_name
 Name used as a prefix for all material properties related to the stress update model. More...
 

Static Protected Attributes

static constexpr unsigned _tensor_dimensionality = 3
 Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics) More...
 

Detailed Description

CappedMohrCoulombCosseratStressUpdate implements rate-independent nonassociative Mohr-Coulomb plus tensile plus compressive plasticity with hardening/softening in the Cosserat setting.

The Mohr-Coulomb plasticity considers the symmetric part of the stress tensor only, and uses an isotropic elasticity tensor that is input by the user (the anti-symmetric parts of the stress tensor and the moment-stress tensor are not included in this plastic model, and any non-isometric parts of the elasticity tensor are ignored in the flow rule).

Definition at line 29 of file CappedMohrCoulombCosseratStressUpdate.h.

Constructor & Destructor Documentation

◆ CappedMohrCoulombCosseratStressUpdate()

CappedMohrCoulombCosseratStressUpdate::CappedMohrCoulombCosseratStressUpdate ( const InputParameters &  parameters)

Definition at line 34 of file CappedMohrCoulombCosseratStressUpdate.C.

36  : CappedMohrCoulombStressUpdate(parameters),
37  _host_young(getParam<Real>("host_youngs_modulus")),
38  _host_poisson(getParam<Real>("host_poissons_ratio")),
41 {
42 }
const Real _host_E0011
E0011 = Lame lambda modulus of the host material.
CappedMohrCoulombStressUpdate(const InputParameters &parameters)
const Real _host_poisson
Poisson&#39;s of the host material.
const Real _host_E0000
E0000 = Lame lambda + 2 * shear modulus of the host material.
const Real _host_young
Young&#39;s modulus of the host material.

Member Function Documentation

◆ calculateRes2()

Real MultiParameterPlasticityStressUpdate::calculateRes2 ( const std::vector< Real > &  rhs) const
protectedinherited

Calculates the residual-squared for the Newton-Raphson + line-search.

Parameters
rhs[in]The RHS vector
Returns
sum_i (rhs[i] * rhs[i])

Definition at line 730 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::lineSearch(), and MultiParameterPlasticityStressUpdate::updateState().

731 {
732  Real res2 = 0.0;
733  for (const auto & r : rhs)
734  res2 += r * r;
735  return res2;
736 }

◆ calculateRHS()

void MultiParameterPlasticityStressUpdate::calculateRHS ( const std::vector< Real > &  trial_stress_params,
const std::vector< Real > &  stress_params,
Real  gaE,
const yieldAndFlow smoothed_q,
std::vector< Real > &  rhs 
) const
protectedinherited

Calculates the RHS in the following 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j] 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, j] * dg/dS[j] ...

0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, j] * dg/dS[j] 0 = rhs[N] = f(S, intnl) where N = _num_sp

Parameters
trial_stress_params[in]The trial stress parameters for this (sub)strain increment, S[:]^trial
stress_params[in]The current stress parameters, S[:]
gaE[in]ga*_En (the normalisation with _En is so that gaE is of similar magnitude to S)
smoothed_q[in]Holds the current value of yield function and derivatives evaluated at the current stress parameters and the current internal parameters
rhs[out]The result

Definition at line 739 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::lineSearch(), and MultiParameterPlasticityStressUpdate::updateState().

744 {
745  const Real ga = gaE / _En;
746  for (unsigned i = 0; i < _num_sp; ++i)
747  {
748  rhs[i] = stress_params[i] - trial_stress_params[i];
749  for (unsigned j = 0; j < _num_sp; ++j)
750  rhs[i] += ga * _Eij[i][j] * smoothed_q.dg[j];
751  }
752  rhs[_num_sp] = smoothed_q.f;
753 }
std::vector< std::vector< Real > > _Eij
E[i, j] in the system of equations to be solved.
const unsigned _num_sp
Number of stress parameters.

◆ computeAllQV()

void CappedMohrCoulombStressUpdate::computeAllQV ( const std::vector< Real > &  stress_params,
const std::vector< Real > &  intnl,
std::vector< yieldAndFlow > &  all_q 
) const
overrideprotectedvirtualinherited

Completely fills all_q with correct values.

These values are: (1) the yield function values, yf[i] (2) d(yf[i])/d(stress_params[j]) (3) d(yf[i])/d(intnl[j]) (4) d(flowPotential[i])/d(stress_params[j]) (5) d2(flowPotential[i])/d(stress_params[j])/d(stress_params[k]) (6) d2(flowPotential[i])/d(stress_params[j])/d(intnl[k])

Parameters
stress_params[in]The stress parameters
intnl[in]The internal parameters
[out]all_qAll the desired quantities

Implements MultiParameterPlasticityStressUpdate.

Definition at line 149 of file CappedMohrCoulombStressUpdate.C.

152 {
153  const Real ts = _tensile_strength.value(intnl[1]);
154  const Real dts = _tensile_strength.derivative(intnl[1]);
155  const Real cs = _compressive_strength.value(intnl[1]);
156  const Real dcs = _compressive_strength.derivative(intnl[1]);
157  const Real sinphi = std::sin(_phi.value(intnl[0]));
158  const Real dsinphi = std::cos(_phi.value(intnl[0])) * _phi.derivative(intnl[0]);
159  const Real sinpsi = std::sin(_psi.value(intnl[0]));
160  const Real dsinpsi = std::cos(_psi.value(intnl[0])) * _psi.derivative(intnl[0]);
161  const Real cohcos = _cohesion.value(intnl[0]) * std::cos(_phi.value(intnl[0]));
162  const Real dcohcos =
163  _cohesion.derivative(intnl[0]) * std::cos(_phi.value(intnl[0])) -
164  _cohesion.value(intnl[0]) * std::sin(_phi.value(intnl[0])) * _phi.derivative(intnl[0]);
165  const Real smax = stress_params[2]; // largest eigenvalue
166  const Real smid = stress_params[1];
167  const Real smin = stress_params[0]; // smallest eigenvalue
168 
169  // yield functions. See comment in yieldFunctionValuesV
170  all_q[0].f = smax - ts;
171  all_q[1].f = smid - ts;
172  all_q[2].f = smin - ts;
173  all_q[3].f = -smin - cs;
174  all_q[4].f = -smid - cs;
175  all_q[5].f = -smax - cs;
176  all_q[6].f = 0.5 * (smax - smin) + 0.5 * (smax + smin) * sinphi - cohcos;
177  all_q[7].f = 0.5 * (smid - smin) + 0.5 * (smid + smin) * sinphi - cohcos;
178  all_q[8].f = 0.5 * (smax - smid) + 0.5 * (smax + smid) * sinphi - cohcos;
179  all_q[9].f = 0.5 * (smid - smax) + 0.5 * (smid + smax) * sinphi - cohcos;
180  all_q[10].f = 0.5 * (smin - smid) + 0.5 * (smin + smid) * sinphi - cohcos;
181  all_q[11].f = 0.5 * (smin - smax) + 0.5 * (smin + smax) * sinphi - cohcos;
182 
183  // d(yield function)/d(stress_params)
184  for (unsigned yf = 0; yf < _num_yf; ++yf)
185  for (unsigned a = 0; a < _num_sp; ++a)
186  all_q[yf].df[a] = 0.0;
187  all_q[0].df[2] = 1.0;
188  all_q[1].df[1] = 1.0;
189  all_q[2].df[0] = 1.0;
190  all_q[3].df[0] = -1.0;
191  all_q[4].df[1] = -1.0;
192  all_q[5].df[2] = -1.0;
193  all_q[6].df[2] = 0.5 * (1.0 + sinphi);
194  all_q[6].df[0] = 0.5 * (-1.0 + sinphi);
195  all_q[7].df[1] = 0.5 * (1.0 + sinphi);
196  all_q[7].df[0] = 0.5 * (-1.0 + sinphi);
197  all_q[8].df[2] = 0.5 * (1.0 + sinphi);
198  all_q[8].df[1] = 0.5 * (-1.0 + sinphi);
199  all_q[9].df[1] = 0.5 * (1.0 + sinphi);
200  all_q[9].df[2] = 0.5 * (-1.0 + sinphi);
201  all_q[10].df[0] = 0.5 * (1.0 + sinphi);
202  all_q[10].df[1] = 0.5 * (-1.0 + sinphi);
203  all_q[11].df[0] = 0.5 * (1.0 + sinphi);
204  all_q[11].df[2] = 0.5 * (-1.0 + sinphi);
205 
206  // d(yield function)/d(intnl)
207  for (unsigned i = 0; i < 6; ++i)
208  all_q[i].df_di[0] = 0.0;
209  all_q[0].df_di[1] = all_q[1].df_di[1] = all_q[2].df_di[1] = -dts;
210  all_q[3].df_di[1] = all_q[4].df_di[1] = all_q[5].df_di[1] = -dcs;
211  for (unsigned i = 6; i < 12; ++i)
212  all_q[i].df_di[1] = 0.0;
213  all_q[6].df_di[0] = 0.5 * (smax + smin) * dsinphi - dcohcos;
214  all_q[7].df_di[0] = 0.5 * (smid + smin) * dsinphi - dcohcos;
215  all_q[8].df_di[0] = 0.5 * (smax + smid) * dsinphi - dcohcos;
216  all_q[9].df_di[0] = 0.5 * (smid + smax) * dsinphi - dcohcos;
217  all_q[10].df_di[0] = 0.5 * (smin + smid) * dsinphi - dcohcos;
218  all_q[11].df_di[0] = 0.5 * (smin + smax) * dsinphi - dcohcos;
219 
220  // the flow potential is just the yield function with phi->psi
221  // d(flow potential)/d(stress_params)
222  for (unsigned yf = 0; yf < 6; ++yf)
223  for (unsigned a = 0; a < _num_sp; ++a)
224  all_q[yf].dg[a] = all_q[yf].df[a];
225  all_q[6].dg[2] = all_q[7].dg[1] = all_q[8].dg[2] = all_q[9].dg[1] = all_q[10].dg[0] =
226  all_q[11].dg[0] = 0.5 * (1.0 + sinpsi);
227  all_q[6].dg[0] = all_q[7].dg[0] = all_q[8].dg[1] = all_q[9].dg[2] = all_q[10].dg[1] =
228  all_q[11].dg[2] = 0.5 * (-1.0 + sinpsi);
229 
230  // d(flow potential)/d(stress_params)/d(intnl)
231  for (unsigned yf = 0; yf < _num_yf; ++yf)
232  for (unsigned a = 0; a < _num_sp; ++a)
233  for (unsigned i = 0; i < _num_intnl; ++i)
234  all_q[yf].d2g_di[a][i] = 0.0;
235  all_q[6].d2g_di[2][0] = all_q[7].d2g_di[1][0] = all_q[8].d2g_di[2][0] = all_q[9].d2g_di[1][0] =
236  all_q[10].d2g_di[0][0] = all_q[11].d2g_di[0][0] = 0.5 * dsinpsi;
237  all_q[6].d2g_di[0][0] = all_q[7].d2g_di[0][0] = all_q[8].d2g_di[1][0] = all_q[9].d2g_di[2][0] =
238  all_q[10].d2g_di[1][0] = all_q[11].d2g_di[2][0] = 0.5 * dsinpsi;
239 
240  // d(flow potential)/d(stress_params)/d(stress_params)
241  for (unsigned yf = 0; yf < _num_yf; ++yf)
242  for (unsigned a = 0; a < _num_sp; ++a)
243  for (unsigned b = 0; b < _num_sp; ++b)
244  all_q[yf].d2g[a][b] = 0.0;
245 }
const TensorMechanicsHardeningModel & _psi
Hardening model for dilation angle.
const TensorMechanicsHardeningModel & _compressive_strength
Hardening model for compressive strength.
const unsigned _num_intnl
Number of internal parameters.
const TensorMechanicsHardeningModel & _tensile_strength
Hardening model for tensile strength.
const unsigned _num_yf
Number of yield functions.
virtual Real derivative(Real intnl) const
const TensorMechanicsHardeningModel & _phi
Hardening model for friction angle.
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
virtual Real value(Real intnl) const
const unsigned _num_sp
Number of stress parameters.

◆ computeStressParams()

void CappedMohrCoulombStressUpdate::computeStressParams ( const RankTwoTensor &  stress,
std::vector< Real > &  stress_params 
) const
overrideprotectedvirtualinherited

Computes stress_params, given stress.

Derived classes must override this

Parameters
stress[in]Stress tensor
stress_params[out]The compute stress_params

Implements MultiParameterPlasticityStressUpdate.

Definition at line 70 of file CappedMohrCoulombStressUpdate.C.

72 {
73  // stress_params[0] = smallest eigenvalue, stress_params[2] = largest eigenvalue
74  stress.symmetricEigenvalues(stress_params);
75 }

◆ computeTimeStepLimit()

Real StressUpdateBase::computeTimeStepLimit ( )
virtualinherited

Reimplemented in RadialReturnStressUpdate.

Definition at line 55 of file StressUpdateBase.C.

56 {
57  return std::numeric_limits<Real>::max();
58 }

◆ consistentTangentOperatorV()

void CappedMohrCoulombCosseratStressUpdate::consistentTangentOperatorV ( const RankTwoTensor &  stress_trial,
const std::vector< Real > &  trial_stress_params,
const RankTwoTensor &  stress,
const std::vector< Real > &  stress_params,
Real  gaE,
const yieldAndFlow smoothed_q,
const RankFourTensor &  Eijkl,
bool  compute_full_tangent_operator,
const std::vector< std::vector< Real >> &  dvar_dtrial,
RankFourTensor &  cto 
)
overrideprotectedvirtual

Calculates the consistent tangent operator.

Derived classes may choose to override this for computational efficiency. The implementation in this class is quite expensive, even though it looks compact and clean, because of all the manipulations of RankFourTensors involved.

Parameters
stress_trial[in]the trial value of the stress tensor for this strain increment
trial_stress_params[in]the trial values of the stress_params for this strain increment
stress[in]the returned value of the stress tensor for this strain increment
stress_params[in]the returned value of the stress_params for this strain increment
gaE[in]the total value of that came from this strain increment
smoothed_q[in]contains the yield function and derivatives evaluated at (p, q)
Eijkl[in]The elasticity tensor
compute_full_tangent_operator[in]true if the full consistent tangent operator is needed, otherwise false
dvar_dtrial[in]dvar_dtrial[i][j] = d({stress_param[i],gaE})/d(trial_stress_param[j]) for this strain increment
[out]ctoThe consistent tangent operator

Add the correction for the antisymmetric part of the elasticity tensor. CappedMohrCoulombStressUpdate computes cto(i, j, k, l) = dstress(i, j)/dstrain(k, l) and during the computations it explicitly performs certain contractions that result in symmetry between i and j, and k and l, viz, cto(i, j, k, l) = cto(j, i, k, l) = cto(i, j, l, k) That is correct because that plasticity model is only valid for symmetric stresses and strains. CappedMohrCoulombCosseratStressUpdate does not include contributions from the antisymmetric parts of stress (or strain), so the antisymmetric parts of cto are just the antisymmetric parts of the elasticity tensor, which must now get added to the cto computed by CappedMohrCoulombStressUpdate

Reimplemented from CappedMohrCoulombStressUpdate.

Definition at line 91 of file CappedMohrCoulombCosseratStressUpdate.C.

102 {
104  trial_stress_params,
105  stress,
106  stress_params,
107  gaE,
108  smoothed_q,
109  elasticity_tensor,
110  compute_full_tangent_operator,
111  dvar_dtrial,
112  cto);
113 
114  if (!compute_full_tangent_operator)
115  return;
116 
133  RankFourTensor anti;
134  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
135  for (unsigned j = 0; j < _tensor_dimensionality; ++j)
136  for (unsigned k = 0; k < _tensor_dimensionality; ++k)
137  for (unsigned l = 0; l < _tensor_dimensionality; ++l)
138  anti(i, j, k, l) = 0.5 * (elasticity_tensor(i, j, k, l) - elasticity_tensor(j, i, k, l));
139 
140  cto += anti;
141 }
virtual void consistentTangentOperatorV(const RankTwoTensor &stress_trial, const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, const std::vector< std::vector< Real >> &dvar_dtrial, RankFourTensor &cto) override
Calculates the consistent tangent operator.
static constexpr unsigned _tensor_dimensionality
Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics) ...

◆ d2stress_param_dstress()

std::vector< RankFourTensor > CappedMohrCoulombStressUpdate::d2stress_param_dstress ( const RankTwoTensor &  stress) const
overrideprotectedvirtualinherited

d2(stress_param[i])/d(stress)/d(stress) at given stress

Parameters
stressstress tensor
Returns
d2(stress_param[:])/d(stress)/d(stress)

Implements MultiParameterPlasticityStressUpdate.

Definition at line 87 of file CappedMohrCoulombStressUpdate.C.

88 {
89  std::vector<RankFourTensor> d2;
90  stress.d2symmetricEigenvalues(d2);
91  return d2;
92 }

◆ dnRHSdVar()

void MultiParameterPlasticityStressUpdate::dnRHSdVar ( const yieldAndFlow smoothed_q,
const std::vector< std::vector< Real >> &  dintnl,
const std::vector< Real > &  stress_params,
Real  gaE,
std::vector< double > &  jac 
) const
protectedinherited

Derivative of -RHS with respect to the stress_params and gaE, placed into an array ready for solving the linear system using LAPACK gsev.

Parameters
smoothed_q[in]Holds the current value of yield function and derivatives evaluated at the current values of the stress_params and the internal parameters
dintnl[in]The derivatives of the internal parameters wrt the stress_params
stress_params[in]The current value of the stress_params during the Newton-Raphson process
gaE[in]The current value of gaE
jac[out]The outputted derivatives

Definition at line 756 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::dVardTrial(), and MultiParameterPlasticityStressUpdate::nrStep().

761 {
762  for (auto & jac_entry : jac)
763  jac_entry = 0.0;
764 
765  const Real ga = gaE / _En;
766 
767  unsigned ind = 0;
768  for (unsigned var = 0; var < _num_sp; ++var)
769  {
770  for (unsigned rhs = 0; rhs < _num_sp; ++rhs)
771  {
772  if (var == rhs)
773  jac[ind] -= 1.0;
774  for (unsigned j = 0; j < _num_sp; ++j)
775  {
776  jac[ind] -= ga * _Eij[rhs][j] * smoothed_q.d2g[j][var];
777  for (unsigned k = 0; k < _num_intnl; ++k)
778  jac[ind] -= ga * _Eij[rhs][j] * smoothed_q.d2g_di[j][k] * dintnl[k][var];
779  }
780  ind++;
781  }
782  // now rhs = _num_sp (that is, the yield function)
783  jac[ind] -= smoothed_q.df[var];
784  for (unsigned k = 0; k < _num_intnl; ++k)
785  jac[ind] -= smoothed_q.df_di[k] * dintnl[k][var];
786  ind++;
787  }
788 
789  // now var = _num_sp (that is, gaE)
790  for (unsigned rhs = 0; rhs < _num_sp; ++rhs)
791  {
792  for (unsigned j = 0; j < _num_sp; ++j)
793  jac[ind] -= (1.0 / _En) * _Eij[rhs][j] * smoothed_q.dg[j];
794  ind++;
795  }
796  // now rhs = _num_sp (that is, the yield function)
797  jac[ind] = 0.0;
798 }
const unsigned _num_intnl
Number of internal parameters.
std::vector< std::vector< Real > > _Eij
E[i, j] in the system of equations to be solved.
const unsigned _num_sp
Number of stress parameters.

◆ dsmoother()

Real MultiParameterPlasticityStressUpdate::dsmoother ( Real  f_diff) const
protectedinherited

Derivative of smoother.

Definition at line 473 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::smoothAllQuantities().

474 {
475  if (std::abs(f_diff) >= _smoothing_tol)
476  return 0.0;
477  return 0.25 * M_PI / _smoothing_tol * std::cos(f_diff * M_PI * 0.5 / _smoothing_tol);
478 }
const Real _smoothing_tol
Smoothing tolerance: edges of the yield surface get smoothed by this amount.

◆ dstress_param_dstress()

std::vector< RankTwoTensor > CappedMohrCoulombStressUpdate::dstress_param_dstress ( const RankTwoTensor &  stress) const
overrideprotectedvirtualinherited

d(stress_param[i])/d(stress) at given stress

Parameters
stressstress tensor
Returns
d(stress_param[:])/d(stress)

Implements MultiParameterPlasticityStressUpdate.

Definition at line 78 of file CappedMohrCoulombStressUpdate.C.

Referenced by CappedMohrCoulombStressUpdate::consistentTangentOperatorV().

79 {
80  std::vector<Real> sp;
81  std::vector<RankTwoTensor> dsp;
82  stress.dsymmetricEigenvalues(sp, dsp);
83  return dsp;
84 }

◆ dVardTrial()

void MultiParameterPlasticityStressUpdate::dVardTrial ( bool  elastic_only,
const std::vector< Real > &  trial_stress_params,
const std::vector< Real > &  stress_params,
Real  gaE,
const std::vector< Real > &  intnl,
const yieldAndFlow smoothed_q,
Real  step_size,
bool  compute_full_tangent_operator,
std::vector< std::vector< Real >> &  dvar_dtrial 
) const
protectedinherited

Calculates derivatives of the stress_params and gaE with repect to the trial values of the stress_params for the (sub)strain increment.

After the strain increment has been fully applied, dvar_dtrial will contain the result appropriate to the full strain increment. Before that time (if applying in sub-strain increments) it will contain the result appropriate to the amount of strain increment applied successfully.

Parameters
elastic_only[in]whether this was an elastic step: if so then the updates to dvar_dtrial are fairly trivial
trial_stress_params[in]Trial values of stress_params for this (sub)strain increment
stress_params[in]Returned values of stress_params for this (sub)strain increment
gaE[in]the value of gaE that came from this (sub)strain increment
intnl[in]the value of the internal parameters at the returned position
smoothed_q[in]contains the yield function and derivatives evaluated at (stress_params, intnl)
step_size[in]size of this (sub)strain increment
compute_full_tangent_operator[in]true if the full consistent tangent operator is needed, otherwise false
dvar_dtrial[out]dvar_dtrial[i][j] = d({stress_param[i],gaE})/d(trial_stress_param[j])

Definition at line 801 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

810 {
811  if (!_fe_problem.currentlyComputingJacobian())
812  return;
813 
814  if (!compute_full_tangent_operator)
815  return;
816 
817  if (elastic_only)
818  {
819  // no change to gaE, and all off-diag stuff remains unchanged from previous step
820  for (unsigned v = 0; v < _num_sp; ++v)
821  dvar_dtrial[v][v] += step_size;
822  return;
823  }
824 
825  const Real ga = gaE / _En;
826 
827  std::vector<std::vector<Real>> dintnl(_num_intnl, std::vector<Real>(_num_sp));
828  setIntnlDerivativesV(trial_stress_params, stress_params, intnl, dintnl);
829 
830  // rhs is described elsewhere, the following are changes in rhs wrt the trial_stress_param
831  // values
832  // In the following we use d(intnl)/d(trial variable) = - d(intnl)/d(variable)
833  std::vector<Real> rhs_cto((_num_sp + 1) * _num_sp);
834 
835  unsigned ind = 0;
836  for (unsigned a = 0; a < _num_sp; ++a)
837  {
838  // change in RHS[b] wrt changes in stress_param_trial[a]
839  for (unsigned b = 0; b < _num_sp; ++b)
840  {
841  if (a == b)
842  rhs_cto[ind] -= 1.0;
843  for (unsigned j = 0; j < _num_sp; ++j)
844  for (unsigned k = 0; k < _num_intnl; ++k)
845  rhs_cto[ind] -= ga * _Eij[b][j] * smoothed_q.d2g_di[j][k] * dintnl[k][a];
846  ind++;
847  }
848  // now b = _num_sp (that is, the yield function)
849  for (unsigned k = 0; k < _num_intnl; ++k)
850  rhs_cto[ind] -= smoothed_q.df_di[k] * dintnl[k][a];
851  ind++;
852  }
853 
854  // jac = d(-rhs)/d(var)
855  std::vector<double> jac((_num_sp + 1) * (_num_sp + 1));
856  dnRHSdVar(smoothed_q, dintnl, stress_params, gaE, jac);
857 
858  std::vector<int> ipiv(_num_sp + 1);
859  int info;
860  const int gesv_num_rhs = _num_sp + 1;
861  const int gesv_num_pq = _num_sp;
862  LAPACKgesv_(&gesv_num_rhs,
863  &gesv_num_pq,
864  &jac[0],
865  &gesv_num_rhs,
866  &ipiv[0],
867  &rhs_cto[0],
868  &gesv_num_rhs,
869  &info);
870  if (info != 0)
871  errorHandler("MultiParameterPlasticityStressUpdate: PETSC LAPACK gsev routine returned with "
872  "error code " +
873  Moose::stringify(info));
874 
875  ind = 0;
876  std::vector<std::vector<Real>> dvarn_dtrialn(_num_sp + 1, std::vector<Real>(_num_sp, 0.0));
877  for (unsigned spt = 0; spt < _num_sp; ++spt) // loop over trial stress-param variables
878  {
879  for (unsigned v = 0; v < _num_sp; ++v) // loop over variables in NR procedure
880  {
881  dvarn_dtrialn[v][spt] = rhs_cto[ind];
882  ind++;
883  }
884  // the final NR variable is gaE
885  dvarn_dtrialn[_num_sp][spt] = rhs_cto[ind];
886  ind++;
887  }
888 
889  const std::vector<std::vector<Real>> dvar_dtrial_old = dvar_dtrial;
890 
891  for (unsigned v = 0; v < _num_sp; ++v) // loop over variables in NR procedure
892  {
893  for (unsigned spt = 0; spt < _num_sp; ++spt) // loop over trial stress-param variables
894  {
895  dvar_dtrial[v][spt] = step_size * dvarn_dtrialn[v][spt];
896  for (unsigned a = 0; a < _num_sp; ++a)
897  dvar_dtrial[v][spt] += dvarn_dtrialn[v][a] * dvar_dtrial_old[a][spt];
898  }
899  }
900  // for gaE the formulae are a little different
901  const unsigned v = _num_sp;
902  for (unsigned spt = 0; spt < _num_sp; ++spt)
903  {
904  dvar_dtrial[v][spt] += step_size * dvarn_dtrialn[v][spt]; // note +=
905  for (unsigned a = 0; a < _num_sp; ++a)
906  dvar_dtrial[v][spt] += dvarn_dtrialn[v][a] * dvar_dtrial_old[a][spt];
907  }
908 }
void dnRHSdVar(const yieldAndFlow &smoothed_q, const std::vector< std::vector< Real >> &dintnl, const std::vector< Real > &stress_params, Real gaE, std::vector< double > &jac) const
Derivative of -RHS with respect to the stress_params and gaE, placed into an array ready for solving ...
virtual void errorHandler(const std::string &message) const
Performs any necessary cleaning-up, then throw MooseException(message)
const unsigned _num_intnl
Number of internal parameters.
std::vector< std::vector< Real > > _Eij
E[i, j] in the system of equations to be solved.
virtual void setIntnlDerivativesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl, std::vector< std::vector< Real >> &dintnl) const =0
Sets the derivatives of internal parameters, based on the trial values of stress_params, their current values, and the current values of the internal parameters.
const unsigned _num_sp
Number of stress parameters.

◆ errorHandler()

void MultiParameterPlasticityStressUpdate::errorHandler ( const std::string &  message) const
protectedvirtualinherited

Performs any necessary cleaning-up, then throw MooseException(message)

Parameters
messageThe message to using in MooseException

Definition at line 589 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::dVardTrial(), and MultiParameterPlasticityStressUpdate::updateState().

590 {
591  throw MooseException(message);
592 }

◆ finalizeReturnProcess()

void MultiParameterPlasticityStressUpdate::finalizeReturnProcess ( const RankTwoTensor &  rotation_increment)
protectedvirtualinherited

Derived classes may use this to perform calculations after the return-map process has completed successfully in stress_param space but before the returned stress tensor has been calculcated.

Parameters
rotation_increment[in]The large-strain rotation increment

Reimplemented in CappedDruckerPragerStressUpdate, CappedWeakPlaneStressUpdate, and CappedWeakInclinedPlaneStressUpdate.

Definition at line 600 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

602 {
603 }

◆ getTangentCalculationMethod()

virtual TangentCalculationMethod MultiParameterPlasticityStressUpdate::getTangentCalculationMethod ( )
inlineoverrideprotectedvirtualinherited

Reimplemented from StressUpdateBase.

Definition at line 122 of file MultiParameterPlasticityStressUpdate.h.

◆ initializeReturnProcess()

void MultiParameterPlasticityStressUpdate::initializeReturnProcess ( )
protectedvirtualinherited

Derived classes may use this to perform calculations before any return-map process is performed, for instance, to initialize variables.

This is called at the very start of updateState, even before any checking for admissible stresses, etc, is performed

Reimplemented in CappedDruckerPragerStressUpdate, CappedWeakPlaneStressUpdate, and CappedWeakInclinedPlaneStressUpdate.

Definition at line 595 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

596 {
597 }

◆ initializeVarsV()

void CappedMohrCoulombStressUpdate::initializeVarsV ( const std::vector< Real > &  trial_stress_params,
const std::vector< Real > &  intnl_old,
std::vector< Real > &  stress_params,
Real &  gaE,
std::vector< Real > &  intnl 
) const
overrideprotectedvirtualinherited

Sets (stress_params, intnl) at "good guesses" of the solution to the Return-Map algorithm.

The purpose of these "good guesses" is to speed up the Newton-Raphson process by providing it with a good initial guess. Derived classes may choose to override this if their plastic models are easy enough to solve approximately. The default values, provided by this class, are simply gaE = 0, stress_params = trial_stress_params, that is, the "good guess" is just the trial point for this (sub)strain increment.

Parameters
trial_stress_params[in]The stress_params at the trial point
intnl_old[in]The internal parameters before applying the (sub)strain increment
stress_params[out]The "good guess" value of the stress_params
gaE[out]The "good guess" value of gaE
intnl[out]The "good guess" value of the internal parameters

For CappedMohrCoulomb there are a number of possibilities for the returned stress configuration:

  • return to the Tensile yield surface to its plane
  • return to the Tensile yield surface to its max=mid edge
  • return to the Tensile yield surface to its tip
  • return to the Compressive yield surface to its plane
  • return to the Compressive yield surface to its max=mid edge
  • return to the Compressive yield surface to its tip
  • return to the Mohr-Coulomb yield surface to its plane
  • return to the Mohr-Coulomb yield surface to its max=mid edge
  • return to the Mohr-Coulomb yield surface to its mid=min edge
  • return to the Mohr-Coulomb yield surface to its tip
  • return to the edge between Tensile and Mohr-Coulomb
  • return to the edge between Tensile and Mohr-Coulomb, at max=mid point
  • return to the edge between Tensile and Mohr-Coulomb, at mid=min point
  • return to the edge between Compressive and Mohr-Coulomb
  • return to the edge between Compressive and Mohr-Coulomb, at max=mid point
  • return to the edge between Compressive and Mohr-Coulomb, at mid=min point Which of these is more prelevant depends on the parameters tensile strength, compressive strength, cohesion, etc. I simply check each possibility until i find one that works. _shifter is used to avoid equal eigenvalues

Reimplemented from MultiParameterPlasticityStressUpdate.

Definition at line 266 of file CappedMohrCoulombStressUpdate.C.

271 {
272  if (!_perfect_guess)
273  {
274  for (unsigned i = 0; i < _num_sp; ++i)
275  stress_params[i] = trial_stress_params[i];
276  gaE = 0.0;
277  }
278  else
279  {
280  const Real ts = _tensile_strength.value(intnl_old[1]);
281  const Real cs = _compressive_strength.value(intnl_old[1]);
282  const Real sinphi = std::sin(_phi.value(intnl_old[0]));
283  const Real cohcos = _cohesion.value(intnl_old[0]) * std::cos(_phi.value(intnl_old[0]));
284 
285  const Real trial_tensile_yf = trial_stress_params[2] - ts;
286  const Real trial_compressive_yf = -trial_stress_params[0] - cs;
287  const Real trial_mc_yf = 0.5 * (trial_stress_params[2] - trial_stress_params[0]) +
288  0.5 * (trial_stress_params[2] + trial_stress_params[0]) * sinphi -
289  cohcos;
290 
316  bool found_solution = false;
317 
318  if (trial_tensile_yf <= _f_tol && trial_compressive_yf <= _f_tol && trial_mc_yf <= _f_tol)
319  {
320  // this is needed because although we know the smoothed yield function is
321  // positive, the individual yield functions may not be
322  for (unsigned i = 0; i < _num_sp; ++i)
323  stress_params[i] = trial_stress_params[i];
324  gaE = 0.0;
325  found_solution = true;
326  }
327 
328  const bool tensile_possible = (ts < cohcos / sinphi); // tensile chops MC tip
329  const bool mc_tip_possible = (cohcos / sinphi < ts); // MC tip pokes through tensile
330  const bool mc_impossible = (0.5 * (ts + cs) + 0.5 * (ts - cs) * sinphi - cohcos <
331  _f_tol); // MC outside tensile and compressive
332 
333  const Real sinpsi = std::sin(_psi.value(intnl_old[0]));
334  const Real halfplus = 0.5 + 0.5 * sinpsi;
335  const Real neghalfplus = -0.5 + 0.5 * sinpsi;
336 
337  if (!found_solution && tensile_possible && trial_tensile_yf > _f_tol &&
338  (trial_compressive_yf <= _f_tol || (trial_compressive_yf > _f_tol && mc_impossible)))
339  {
340  // try pure tensile failure, return to the plane
341  // This involves solving yf[0] = 0 and the three flow-direction equations
342  // Don't try this if there is compressive failure, since returning to
343  // the tensile yield surface will only make compressive failure worse
344  const Real ga = (trial_stress_params[2] - ts) / _Eij[2][2];
345  stress_params[2] = ts; // largest eigenvalue
346  stress_params[1] = trial_stress_params[1] - ga * _Eij[1][2];
347  stress_params[0] = trial_stress_params[0] - ga * _Eij[0][2];
348 
349  // if we have to return to the edge, or tip, do that
350  Real dist_mod = 1.0;
351  const Real to_subtract1 = stress_params[1] - (ts - 0.5 * _shifter);
352  if (to_subtract1 > 0.0)
353  {
354  dist_mod += Utility::pow<2>(to_subtract1 / (trial_stress_params[2] - ts));
355  stress_params[1] -= to_subtract1;
356  }
357  const Real to_subtract0 = stress_params[0] - (ts - _shifter);
358  if (to_subtract0 > 0.0)
359  {
360  dist_mod += Utility::pow<2>(to_subtract0 / (trial_stress_params[2] - ts));
361  stress_params[0] -= to_subtract0;
362  }
363  if (mc_impossible) // might have to shift up to the compressive yield surface
364  {
365  const Real to_add0 = -stress_params[0] - cs;
366  if (to_add0 > 0.0)
367  {
368  dist_mod += Utility::pow<2>(to_add0 / (trial_stress_params[2] - ts));
369  stress_params[0] += to_add0;
370  }
371  const Real to_add1 = -cs + 0.5 * _shifter - stress_params[1];
372  if (to_add1 > 0.0)
373  {
374  dist_mod += Utility::pow<2>(to_add1 / (trial_stress_params[2] - ts));
375  stress_params[1] += to_add1;
376  }
377  }
378 
379  const Real new_compressive_yf = -stress_params[0] - cs;
380  const Real new_mc_yf = 0.5 * (stress_params[2] - stress_params[0]) +
381  0.5 * (stress_params[2] + stress_params[0]) * sinphi - cohcos;
382  if (new_mc_yf <= _f_tol && new_compressive_yf <= _f_tol)
383  {
384  gaE = std::sqrt(dist_mod) * (trial_stress_params[2] - stress_params[2]);
385  found_solution = true;
386  }
387  }
388  if (!found_solution && trial_compressive_yf > _f_tol &&
389  (trial_tensile_yf <= _f_tol || (trial_tensile_yf > _f_tol && mc_impossible)))
390  {
391  // try pure compressive failure
392  // Don't try this if there is tensile failure, since returning to
393  // the compressive yield surface will only make tensile failure worse
394  const Real ga = (trial_stress_params[0] + cs) / _Eij[0][0]; // this is negative
395  stress_params[0] = -cs;
396  stress_params[1] = trial_stress_params[1] - ga * _Eij[1][0];
397  stress_params[2] = trial_stress_params[2] - ga * _Eij[2][0];
398 
399  // if we have to return to the edge, or tip, do that
400  Real dist_mod = 1.0;
401  const Real to_add1 = -cs + 0.5 * _shifter - stress_params[1];
402  if (to_add1 > 0.0)
403  {
404  dist_mod += Utility::pow<2>(to_add1 / (trial_stress_params[0] + cs));
405  stress_params[1] += to_add1;
406  }
407  const Real to_add2 = -cs + _shifter - stress_params[2];
408  if (to_add2 > 0.0)
409  {
410  dist_mod += Utility::pow<2>(to_add2 / (trial_stress_params[0] + cs));
411  stress_params[2] += to_add2;
412  }
413  if (mc_impossible) // might have to shift down to the tensile yield surface
414  {
415  const Real to_subtract2 = stress_params[2] - ts;
416  if (to_subtract2 > 0.0)
417  {
418  dist_mod += Utility::pow<2>(to_subtract2 / (trial_stress_params[0] + cs));
419  stress_params[2] -= to_subtract2;
420  }
421  const Real to_subtract1 = stress_params[1] - (ts - 0.5 * _shifter);
422  if (to_subtract1 > 0.0)
423  {
424  dist_mod += Utility::pow<2>(to_subtract1 / (trial_stress_params[0] + cs));
425  stress_params[1] -= to_subtract1;
426  }
427  }
428 
429  const Real new_tensile_yf = stress_params[2] - ts;
430  const Real new_mc_yf = 0.5 * (stress_params[2] - stress_params[0]) +
431  0.5 * (stress_params[2] + stress_params[0]) * sinphi - cohcos;
432  if (new_mc_yf <= _f_tol && new_tensile_yf <= _f_tol)
433  {
434  gaE = std::sqrt(dist_mod) * (-trial_stress_params[0] + stress_params[0]);
435  found_solution = true;
436  }
437  }
438  if (!found_solution && !mc_impossible && trial_mc_yf > _f_tol)
439  {
440  // try pure shear failure, return to the plane
441  // This involves solving yf[6]=0 and the three flow-direction equations
442  const Real ga = ((trial_stress_params[2] - trial_stress_params[0]) +
443  (trial_stress_params[2] + trial_stress_params[0]) * sinphi - 2.0 * cohcos) /
444  (_Eij[2][2] - _Eij[0][2] + (_Eij[2][2] + _Eij[0][2]) * sinpsi * sinphi);
445  stress_params[2] =
446  trial_stress_params[2] - ga * (_Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus);
447  stress_params[1] = trial_stress_params[1] - ga * _Eij[1][0] * sinpsi;
448  stress_params[0] =
449  trial_stress_params[0] - ga * (_Eij[0][0] * neghalfplus + _Eij[0][2] * halfplus);
450  const Real f7 = 0.5 * (stress_params[1] - stress_params[0]) +
451  0.5 * (stress_params[1] + stress_params[0]) * sinphi - cohcos;
452  const Real f8 = 0.5 * (stress_params[2] - stress_params[1]) +
453  0.5 * (stress_params[2] + stress_params[1]) * sinphi - cohcos;
454  const Real new_tensile_yf = stress_params[2] - ts;
455  const Real new_compressive_yf = -stress_params[0] - cs;
456 
457  if (f7 <= _f_tol && f8 <= _f_tol && new_tensile_yf <= _f_tol && new_compressive_yf <= _f_tol)
458  {
459  gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
460  (stress_params[2] - stress_params[0])) /
461  (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2];
462  found_solution = true;
463  }
464  }
465  if (!found_solution && !mc_impossible && trial_mc_yf > _f_tol)
466  {
467  // Try return to the max=mid MC line.
468  // To return to the max=mid line, we need to solve f6 = 0 = f7 and
469  // the three flow-direction equations. In the flow-direction equations
470  // there are two plasticity multipliers, which i call ga6 and ga7,
471  // corresponding to the amounts of strain normal to the f6 and f7
472  // directions, respectively.
473  // So:
474  // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga7 Emax,a dg7/dSa
475  // = Smax^trial - ga6 cmax6 - ga7 cmax7 (with cmax6 and cmax7 evaluated below)
476  // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga7 Emid,a dg7/dSa
477  // = Smid^trial - ga6 cmid6 - ga7 cmid7
478  // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga7 Emin,a dg7/dSa
479  // = Smin^trial - ga6 cmin6 - ga7 cmin7
480  const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
481  const Real cmax7 = _Eij[2][1] * halfplus + _Eij[2][0] * neghalfplus;
482  // const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
483  // const Real cmid7 = _Eij[1][1] * halfplus + _Eij[1][0] * neghalfplus;
484  const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
485  const Real cmin7 = _Eij[0][1] * halfplus + _Eij[0][0] * neghalfplus;
486  // Substituting these into f6 = 0 yields
487  // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga7 (0.5(cmax6 -
488  // cmin6) + 0.5(cmax6 + cmin6)sinphi) = f6_trial - ga6 c6 - ga7 c7, where
489  const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
490  const Real c7 = 0.5 * (cmax7 - cmin7) + 0.5 * (cmax7 + cmin7) * sinphi;
491  // It isn't too hard to check that the other equation is
492  // 0 = f7_trial - ga6 c7 - ga7 c6
493  // These equations may be inverted to yield
494  if (c6 != c7)
495  {
496  const Real f6_trial = trial_mc_yf;
497  const Real f7_trial = 0.5 * (trial_stress_params[1] - trial_stress_params[0]) +
498  0.5 * (trial_stress_params[1] + trial_stress_params[0]) * sinphi -
499  cohcos;
500  const Real descr = Utility::pow<2>(c6) - Utility::pow<2>(c7);
501  Real ga6 = (c6 * f6_trial - c7 * f7_trial) / descr;
502  Real ga7 = (-c7 * f6_trial + c6 * f7_trial) / descr;
503  // and finally
504  stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga7 * cmax7;
505  stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga7 * cmin7;
506  stress_params[1] = stress_params[2] - 0.5 * _shifter;
507 
508  Real f8 = 0.5 * (stress_params[2] - stress_params[1]) +
509  0.5 * (stress_params[2] + stress_params[1]) * sinphi - cohcos;
510 
511  if (mc_tip_possible && f8 > _f_tol)
512  {
513  stress_params[2] = cohcos / sinphi;
514  stress_params[1] = stress_params[2] - 0.5 * _shifter;
515  stress_params[0] = stress_params[2] - _shifter;
516  f8 = 0.0;
517  ga6 = 1.0;
518  ga7 = 1.0;
519  }
520 
521  const Real new_tensile_yf = stress_params[2] - ts;
522  const Real new_compressive_yf = -stress_params[0] - cs;
523 
524  if (f8 <= _f_tol && new_tensile_yf <= _f_tol && new_compressive_yf <= _f_tol &&
525  ga6 >= 0.0 && ga7 >= 0.0)
526  {
527  gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
528  (stress_params[2] - stress_params[0])) /
529  (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2];
530  found_solution = true;
531  }
532  }
533  }
534  if (!found_solution && !mc_impossible && trial_mc_yf > _f_tol)
535  {
536  // Try return to the mid=min line.
537  // To return to the mid=min line, we need to solve f6 = 0 = f8 and
538  // the three flow-direction equations. In the flow-direction equations
539  // there are two plasticity multipliers, which i call ga6 and ga8,
540  // corresponding to the amounts of strain normal to the f6 and f8
541  // directions, respectively.
542  // So:
543  // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga8 Emax,a dg8/dSa
544  // = Smax^trial - ga6 cmax6 - ga8 cmax8 (with cmax6 and cmax8 evaluated below)
545  // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga8 Emid,a dg8/dSa
546  // = Smid^trial - ga6 cmid6 - ga8 cmid8
547  // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga8 Emin,a dg8/dSa
548  // = Smin^trial - ga6 cmin6 - ga8 cmin8
549  const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
550  const Real cmax8 = _Eij[2][2] * halfplus + _Eij[2][1] * neghalfplus;
551  // const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
552  // const Real cmid8 = _Eij[1][2] * halfplus + _Eij[1][1] * neghalfplus;
553  const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
554  const Real cmin8 = _Eij[0][2] * halfplus + _Eij[0][1] * neghalfplus;
555  // Substituting these into f6 = 0 yields
556  // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga8 (0.5(cmax6 -
557  // cmin6) + 0.5(cmax6 + cmin6)sinphi) = f6_trial - ga6 c6 - ga8 c8, where
558  const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
559  const Real c8 = 0.5 * (cmax8 - cmin8) + 0.5 * (cmax8 + cmin8) * sinphi;
560  // It isn't too hard to check that the other equation is
561  // 0 = f8_trial - ga6 c8 - ga8 c6
562  // These equations may be inverted to yield
563  if (c6 != c8)
564  {
565  const Real f6_trial = trial_mc_yf;
566  const Real f8_trial = 0.5 * (trial_stress_params[2] - trial_stress_params[1]) +
567  0.5 * (trial_stress_params[2] + trial_stress_params[1]) * sinphi -
568  cohcos;
569  const Real descr = Utility::pow<2>(c6) - Utility::pow<2>(c8);
570  Real ga6 = (c6 * f6_trial - c8 * f8_trial) / descr;
571  Real ga8 = (-c8 * f6_trial + c6 * f8_trial) / descr;
572  // and finally
573  stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga8 * cmax8;
574  stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga8 * cmin8;
575  stress_params[1] = stress_params[0] + 0.5 * _shifter;
576 
577  Real f7 = 0.5 * (stress_params[1] - stress_params[0]) +
578  0.5 * (stress_params[1] + stress_params[0]) * sinphi - cohcos;
579 
580  if (mc_tip_possible && f7 > _f_tol)
581  {
582  stress_params[2] = cohcos / sinphi;
583  stress_params[1] = stress_params[2] - 0.5 * _shifter;
584  stress_params[0] = stress_params[2] - _shifter;
585  f7 = 0.0;
586  ga6 = 1.0;
587  ga8 = 1.0;
588  }
589 
590  const Real new_tensile_yf = stress_params[2] - ts;
591  const Real new_compressive_yf = -stress_params[0] - cs;
592 
593  if (f7 <= _f_tol && new_tensile_yf <= _f_tol && new_compressive_yf <= _f_tol &&
594  ga6 >= 0.0 && ga8 >= 0.0)
595  {
596  gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
597  (stress_params[2] - stress_params[0])) /
598  (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2];
599  found_solution = true;
600  }
601  }
602  }
603  if (!found_solution && !mc_impossible && tensile_possible && trial_tensile_yf > _f_tol)
604  {
605  // Return to the line where yf[0] = 0 = yf[6].
606  // To return to this line, we need to solve f0 = 0 = f6 and
607  // the three flow-direction equations. In the flow-direction equations
608  // there are two plasticity multipliers, which i call ga0 and ga6
609  // corresponding to the amounts of strain normal to the f0 and f6
610  // directions, respectively.
611  // So:
612  // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga0 Emax,a dg0/dSa
613  // = Smax^trial - ga6 cmax6 - ga0 cmax0 (with cmax6 and cmax0 evaluated below)
614  // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga0 Emid,a dg0/dSa
615  // = Smid^trial - ga6 cmid6 - ga0 cmid0
616  // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga0 Emin,a dg0/dSa
617  // = Smin^trial - ga6 cmin6 - ga0 cmin0
618  const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
619  const Real cmax0 = _Eij[2][2];
620  const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
621  const Real cmid0 = _Eij[1][2];
622  const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
623  const Real cmin0 = _Eij[0][2];
624  // Substituting these into f6 = 0 yields
625  // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga0 (0.5(cmax0 -
626  // cmin0) + 0.5(cmax0 + cmin0)sinphi) = f6_trial - ga6 c6 - ga0 c0, where
627  const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
628  const Real c0 = 0.5 * (cmax0 - cmin0) + 0.5 * (cmax0 + cmin0) * sinphi;
629  // Substituting these into f0 = 0 yields
630  // 0 = f0_trial - ga6 cmax6 - ga0 cmax0
631  // These equations may be inverted to yield the following
632  const Real descr = c0 * cmax6 - c6 * cmax0;
633  if (descr != 0.0)
634  {
635  const Real ga0 = (-c6 * trial_tensile_yf + cmax6 * trial_mc_yf) / descr;
636  const Real ga6 = (c0 * trial_tensile_yf - cmax0 * trial_mc_yf) / descr;
637  stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga0 * cmax0;
638  stress_params[1] = trial_stress_params[1] - ga6 * cmid6 - ga0 * cmid0;
639  stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga0 * cmin0;
640 
641  // enforce physicality (go to corners if necessary)
642  stress_params[0] =
643  std::min(stress_params[0],
644  stress_params[2] - _shifter); // account for poor choice from user
645  stress_params[1] = std::max(std::min(stress_params[1], stress_params[2] - 0.5 * _shifter),
646  stress_params[0] + 0.5 * _shifter);
647 
648  const Real new_compressive_yf = -stress_params[0] - cs;
649  if (new_compressive_yf <= _f_tol && ga0 >= 0.0 && ga6 >= 0.0) // enforce ga>0!
650  {
651  gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
652  (stress_params[2] - stress_params[0])) /
653  (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2] +
654  (trial_stress_params[2] - stress_params[2]);
655  found_solution = true;
656  }
657  }
658  }
659  if (!found_solution && !mc_impossible)
660  {
661  // Return to the line where yf[3] = 0 = yf[6].
662  // To return to this line, we need to solve f3 = 0 = f6 and
663  // the three flow-direction equations. In the flow-direction equations
664  // there are two plasticity multipliers, which i call ga3 and ga6
665  // corresponding to the amounts of strain normal to the f3 and f6
666  // directions, respectively.
667  // So:
668  // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga3 Emax,a dg3/dSa
669  // = Smax^trial - ga6 cmax6 - ga3 cmax3 (with cmax6 and cmax3 evaluated below)
670  // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga3 Emid,a dg3/dSa
671  // = Smid^trial - ga6 cmid6 - ga3 cmid3
672  // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga3 Emin,a dg3/dSa
673  // = Smin^trial - ga6 cmin6 - ga3 cmin3
674  const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
675  const Real cmax3 = -_Eij[2][0];
676  const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
677  const Real cmid3 = -_Eij[1][0];
678  const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
679  const Real cmin3 = -_Eij[0][0];
680  // Substituting these into f6 = 0 yields
681  // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga3 (0.5(cmax3 -
682  // cmin3) + 0.5(cmax3 + cmin3)sinphi) = f6_trial - ga6 c6 - ga3 c3, where
683  const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
684  const Real c3 = 0.5 * (cmax3 - cmin3) + 0.5 * (cmax3 + cmin3) * sinphi;
685  // Substituting these into f3 = 0 yields
686  // 0 = - f3_trial - ga6 cmin6 - ga3 cmin3
687  // These equations may be inverted to yield the following
688  const Real descr = c3 * cmin6 - c6 * cmin3;
689  if (descr != 0.0)
690  {
691  const Real ga3 = (c6 * trial_compressive_yf + cmin6 * trial_mc_yf) / descr;
692  const Real ga6 = (-c3 * trial_compressive_yf - cmin3 * trial_mc_yf) / descr;
693  stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga3 * cmax3;
694  stress_params[1] = trial_stress_params[1] - ga6 * cmid6 - ga3 * cmid3;
695  stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga3 * cmin3;
696 
697  const Real new_tensile_yf = stress_params[2] - ts;
698  stress_params[2] =
699  std::max(stress_params[2],
700  stress_params[0] + _shifter); // account for poor choice from user
701  stress_params[1] = std::max(std::min(stress_params[1], stress_params[2] - 0.5 * _shifter),
702  stress_params[0] + 0.5 * _shifter);
703 
704  if (new_tensile_yf <= _f_tol && ga6 >= 0.0)
705  {
706  gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
707  (stress_params[2] - stress_params[0])) /
708  (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2] +
709  (-trial_stress_params[0] - stress_params[0]);
710 
711  found_solution = true;
712  }
713  }
714  }
715  if (!found_solution)
716  {
717  // Cannot find an acceptable initialisation
718  for (unsigned i = 0; i < _num_sp; ++i)
719  stress_params[i] = trial_stress_params[i];
720  gaE = 0.0;
721  mooseWarning("CappedMohrCoulombStressUpdate cannot initialize from max = ",
722  stress_params[2],
723  " mid = ",
724  stress_params[1],
725  " min = ",
726  stress_params[0]);
727  }
728  }
729  setIntnlValuesV(trial_stress_params, stress_params, intnl_old, intnl);
730 }
const TensorMechanicsHardeningModel & _psi
Hardening model for dilation angle.
const bool _perfect_guess
Whether to provide an estimate of the returned stress, based on perfect plasticity.
const TensorMechanicsHardeningModel & _compressive_strength
Hardening model for compressive strength.
const TensorMechanicsHardeningModel & _tensile_strength
Hardening model for tensile strength.
std::vector< std::vector< Real > > _Eij
E[i, j] in the system of equations to be solved.
const Real _f_tol
The yield-function tolerance.
void setIntnlValuesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &intnl) const override
Sets the internal parameters based on the trial values of stress_params, their current values...
const Real _shifter
When equal-eigenvalues are predicted from the stress initialization routine, shift them by this amoun...
const TensorMechanicsHardeningModel & _phi
Hardening model for friction angle.
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
virtual Real value(Real intnl) const
const unsigned _num_sp
Number of stress parameters.

◆ initQpStatefulProperties()

void MultiParameterPlasticityStressUpdate::initQpStatefulProperties ( )
overrideprotectedvirtualinherited

Reimplemented in CappedWeakInclinedPlaneStressUpdate.

Definition at line 126 of file MultiParameterPlasticityStressUpdate.C.

Referenced by CappedWeakInclinedPlaneStressUpdate::initQpStatefulProperties().

127 {
128  _plastic_strain[_qp].zero();
129  _intnl[_qp].assign(_num_intnl, 0);
130  _yf[_qp].assign(_num_yf, 0);
131  _iter[_qp] = 0.0;
132  _linesearch_needed[_qp] = 0.0;
133 }
MaterialProperty< std::vector< Real > > & _yf
yield functions
const unsigned _num_intnl
Number of internal parameters.
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
const unsigned _num_yf
Number of yield functions.
MaterialProperty< Real > & _iter
Number of Newton-Raphson iterations used in the return-map.
MaterialProperty< RankTwoTensor > & _plastic_strain
plastic strain
MaterialProperty< Real > & _linesearch_needed
Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) ...

◆ ismoother()

Real MultiParameterPlasticityStressUpdate::ismoother ( Real  f_diff) const
protectedinherited

Smooths yield functions.

The returned value must be zero if abs(f_diff) >= _smoothing_tol and otherwise must satisfy, over -_smoothing_tol <= f_diff <= _smoothing_tol: (1) C2 (2) zero at f_diff = +/- _smoothing_tol (3) derivative is +/-0.5 at f_diff = +/- _smoothing_tol (4) derivative must be in [-0.5, 0.5] (5) second derivative is zero at f_diff = +/- _smoothing_tol (6) second derivative must be non-negative in order to ensure C2 differentiability and convexity of the smoothed yield surface.

Definition at line 457 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::smoothAllQuantities(), and MultiParameterPlasticityStressUpdate::yieldF().

458 {
459  if (std::abs(f_diff) >= _smoothing_tol)
460  return 0.0;
461  return -_smoothing_tol / M_PI * std::cos(0.5 * M_PI * f_diff / _smoothing_tol);
462 }
const Real _smoothing_tol
Smoothing tolerance: edges of the yield surface get smoothed by this amount.

◆ lineSearch()

int MultiParameterPlasticityStressUpdate::lineSearch ( Real &  res2,
std::vector< Real > &  stress_params,
Real &  gaE,
const std::vector< Real > &  trial_stress_params,
yieldAndFlow smoothed_q,
const std::vector< Real > &  intnl_ok,
std::vector< Real > &  intnl,
std::vector< Real > &  rhs,
Real &  linesearch_needed 
) const
protectedinherited

Performs a line-search to find stress_params and gaE Upon entry:

  • rhs contains the solution to the Newton-Raphson (ie nrStep should have been called). If a full Newton step is used then stress_params[:] += rhs[0:_num_sp-1] and gaE += rhs[_num_sp]
  • res2 contains the residual-squared before applying any of solution
  • stress_params contains the stress_params before applying any of the solution
  • gaE contains gaE before applying any of the solution (that is contained in rhs) Upon exit:
  • stress_params will be the stress_params after applying the solution
  • gaE will be the stress_params after applying the solution
  • rhs will contain the updated rhs values (after applying the solution) ready for the next Newton-Raphson step,
  • res2 will be the residual-squared after applying the solution
  • intnl will contain the internal variables corresponding to the return from trial_stress_params to stress_params (and starting from intnl_ok)
  • linesearch_needed will be 1.0 if a linesearch was needed
  • smoothed_q will contain the value of the yield function and its derivatives, etc, at (stress_params, intnl)
    Parameters
    res2[in,out]the residual-squared, both as an input and output
    stress_params[in,out]Upon input the value of the stress_params before the current Newton-Raphson process was initiated. Upon exit this will hold the values coming from the line search.
    trial_stress_params[in]Trial value for the stress_params for this (sub)strain increment
    gaE[in,out]Upon input the value of gaE before the current Newton-Raphson iteration was initiated. Upon exit this will hold the value coming from the line-search
    smoothed_q[in,out]Upon input, the value of the smoothed yield function and derivatives at the prior-to-Newton configuration. Upon exit this is evaluated at the new (stress_params, intnl)
    intnl_ok[in]The value of the internal parameters from the start of this (sub)strain increment
    intnl[in,out]The value of the internal parameters after the line-search has converged
    rhs[in,out]Upon entry this contains the solution to the Newton-Raphson. Upon exit this contains the updated rhs values
    Returns
    0 if successful, 1 otherwise

Definition at line 481 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

490 {
491  const Real res2_old = res2;
492  const std::vector<Real> sp_params_old = stress_params;
493  const Real gaE_old = gaE;
494  const std::vector<Real> delta_nr_params = rhs;
495 
496  Real lam = 1.0; // line-search parameter
497  const Real lam_min = 1E-10; // minimum value of lam allowed
498  const Real slope = -2.0 * res2_old; // "Numerical Recipes" uses -b*A*x, in order to check for
499  // roundoff, but i hope the nrStep would warn if there were
500  // problems
501  Real tmp_lam; // cached value of lam used in quadratic & cubic line search
502  Real f2 = res2_old; // cached value of f = residual2 used in the cubic in the line search
503  Real lam2 = lam; // cached value of lam used in the cubic in the line search
504 
505  while (true)
506  {
507  // update variables using the current line-search parameter
508  for (unsigned i = 0; i < _num_sp; ++i)
509  stress_params[i] = sp_params_old[i] + lam * delta_nr_params[i];
510  gaE = gaE_old + lam * delta_nr_params[_num_sp];
511 
512  // and internal parameters
513  setIntnlValuesV(trial_stress_params, stress_params, intnl_ok, intnl);
514 
515  smoothed_q = smoothAllQuantities(stress_params, intnl);
516 
517  // update rhs for next-time through
518  calculateRHS(trial_stress_params, stress_params, gaE, smoothed_q, rhs);
519  res2 = calculateRes2(rhs);
520 
521  // do the line-search
522  if (res2 < res2_old + 1E-4 * lam * slope)
523  break;
524  else if (lam < lam_min)
525  return 1;
526  else if (lam == 1.0)
527  {
528  // model as a quadratic
529  tmp_lam = -0.5 * slope / (res2 - res2_old - slope);
530  }
531  else
532  {
533  // model as a cubic
534  const Real rhs1 = res2 - res2_old - lam * slope;
535  const Real rhs2 = f2 - res2_old - lam2 * slope;
536  const Real a = (rhs1 / Utility::pow<2>(lam) - rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
537  const Real b =
538  (-lam2 * rhs1 / Utility::pow<2>(lam) + lam * rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
539  if (a == 0.0)
540  tmp_lam = -slope / (2.0 * b);
541  else
542  {
543  const Real disc = Utility::pow<2>(b) - 3.0 * a * slope;
544  if (disc < 0)
545  tmp_lam = 0.5 * lam;
546  else if (b <= 0)
547  tmp_lam = (-b + std::sqrt(disc)) / (3.0 * a);
548  else
549  tmp_lam = -slope / (b + std::sqrt(disc));
550  }
551  if (tmp_lam > 0.5 * lam)
552  tmp_lam = 0.5 * lam;
553  }
554  lam2 = lam;
555  f2 = res2;
556  lam = std::max(tmp_lam, 0.1 * lam);
557  }
558 
559  if (lam < 1.0)
560  linesearch_needed = 1.0;
561  return 0;
562 }
virtual void setIntnlValuesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &intnl) const =0
Sets the internal parameters based on the trial values of stress_params, their current values...
void calculateRHS(const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, std::vector< Real > &rhs) const
Calculates the RHS in the following 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j] 0 = rhs[...
yieldAndFlow smoothAllQuantities(const std::vector< Real > &stress_params, const std::vector< Real > &intnl) const
Calculates all yield functions and derivatives, and then performs the smoothing scheme.
Real calculateRes2(const std::vector< Real > &rhs) const
Calculates the residual-squared for the Newton-Raphson + line-search.
const unsigned _num_sp
Number of stress parameters.

◆ nrStep()

int MultiParameterPlasticityStressUpdate::nrStep ( const yieldAndFlow smoothed_q,
const std::vector< Real > &  trial_stress_params,
const std::vector< Real > &  stress_params,
const std::vector< Real > &  intnl,
Real  gaE,
std::vector< Real > &  rhs 
) const
protectedinherited

Performs a Newton-Raphson step to attempt to zero rhs Upon return, rhs will contain the solution.

Parameters
smoothed_q[in]The value of the smoothed yield function and derivatives prior to this Newton-Raphson step
trial_stress_params[in]Trial value for the stress_params for this (sub)strain increment
stress_params[in]The current value of the stress_params
intnl[in]The current value of the internal parameters
gaE[in]The current value of gaE
rhs[in,out]Upon entry, the rhs to zero using Newton-Raphson. Upon exit, the solution to the Newton-Raphson problem
Returns
0 if successful, 1 otherwise

Definition at line 565 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

571 {
572  std::vector<std::vector<Real>> dintnl(_num_intnl, std::vector<Real>(_num_sp));
573  setIntnlDerivativesV(trial_stress_params, stress_params, intnl, dintnl);
574 
575  std::vector<double> jac((_num_sp + 1) * (_num_sp + 1));
576  dnRHSdVar(smoothed_q, dintnl, stress_params, gaE, jac);
577 
578  // use LAPACK to solve the linear system
579  const int nrhs = 1;
580  std::vector<int> ipiv(_num_sp + 1);
581  int info;
582  const int gesv_num_rhs = _num_sp + 1;
583  LAPACKgesv_(
584  &gesv_num_rhs, &nrhs, &jac[0], &gesv_num_rhs, &ipiv[0], &rhs[0], &gesv_num_rhs, &info);
585  return info;
586 }
void dnRHSdVar(const yieldAndFlow &smoothed_q, const std::vector< std::vector< Real >> &dintnl, const std::vector< Real > &stress_params, Real gaE, std::vector< double > &jac) const
Derivative of -RHS with respect to the stress_params and gaE, placed into an array ready for solving ...
const unsigned _num_intnl
Number of internal parameters.
virtual void setIntnlDerivativesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl, std::vector< std::vector< Real >> &dintnl) const =0
Sets the derivatives of internal parameters, based on the trial values of stress_params, their current values, and the current values of the internal parameters.
const unsigned _num_sp
Number of stress parameters.

◆ precisionLoss()

bool MultiParameterPlasticityStressUpdate::precisionLoss ( const std::vector< Real > &  solution,
const std::vector< Real > &  stress_params,
Real  gaE 
) const
protectedinherited

Check whether precision loss has occurred.

Parameters
[in]solutionThe solution to the Newton-Raphson system
[in]stress_paramsThe currect values of the stress_params for this (sub)strain increment
[in]gaEThe currenct value of gaE for this (sub)strain increment
Returns
true if precision loss has occurred

Definition at line 911 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

914 {
915  if (std::abs(solution[_num_sp]) > 1E-13 * std::abs(gaE))
916  return false;
917  for (unsigned i = 0; i < _num_sp; ++i)
918  if (std::abs(solution[i]) > 1E-13 * std::abs(stress_params[i]))
919  return false;
920  return true;
921 }
const unsigned _num_sp
Number of stress parameters.

◆ preReturnMapV()

void CappedMohrCoulombCosseratStressUpdate::preReturnMapV ( const std::vector< Real > &  trial_stress_params,
const RankTwoTensor &  stress_trial,
const std::vector< Real > &  intnl_old,
const std::vector< Real > &  yf,
const RankFourTensor &  Eijkl 
)
overrideprotectedvirtual

Derived classes may employ this function to record stuff or do other computations prior to the return-mapping algorithm.

We know that (trial_stress_params, intnl_old) is inadmissible when this is called

Parameters
trial_stress_params[in]The trial values of the stress parameters
stress_trial[in]Trial stress tensor
intnl_old[in]Old value of the internal parameters.
yf[in]The yield functions at (p_trial, q_trial, intnl_old)
Eijkl[in]The elasticity tensor

Reimplemented from CappedMohrCoulombStressUpdate.

Definition at line 45 of file CappedMohrCoulombCosseratStressUpdate.C.

51 {
52  std::vector<Real> eigvals;
53  stress_trial.symmetricEigenvaluesEigenvectors(eigvals, _eigvecs);
55 }
const Real _host_poisson
Poisson&#39;s of the host material.
RankTwoTensor _eigvecs
Eigenvectors of the trial stress as a RankTwoTensor, in order to rotate the returned stress back to s...

◆ propagateQpStatefulProperties()

void MultiParameterPlasticityStressUpdate::propagateQpStatefulProperties ( )
overrideprotectedvirtualinherited

If updateState is not called during a timestep, this will be.

This method allows derived classes to set internal parameters from their Old values, for instance

Reimplemented from StressUpdateBase.

Definition at line 136 of file MultiParameterPlasticityStressUpdate.C.

137 {
139  for (unsigned i = 0; i < _num_intnl; ++i)
140  _intnl[_qp][i] = _intnl_old[_qp][i];
141 }
const unsigned _num_intnl
Number of internal parameters.
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
const MaterialProperty< std::vector< Real > > & _intnl_old
old values of internal parameters
MaterialProperty< RankTwoTensor > & _plastic_strain
plastic strain
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
Old value of plastic strain.

◆ requiresIsotropicTensor()

bool CappedMohrCoulombCosseratStressUpdate::requiresIsotropicTensor ( )
inlineoverridevirtual

The full elasticity tensor may be anisotropic, and usually is in the case of layered Cosserat.

However, this class only uses the isotropic parts of it (corresponding to the "host" material) that are encoded in _host_young and _host_poisson

Implements StressUpdateBase.

Definition at line 40 of file CappedMohrCoulombCosseratStressUpdate.h.

40 { return false; }

◆ resetProperties()

void StressUpdateBase::resetProperties ( )
inlinefinalinherited

Definition at line 117 of file StressUpdateBase.h.

117 {}

◆ resetQpProperties()

void StressUpdateBase::resetQpProperties ( )
inlinefinalinherited

Retained as empty methods to avoid a warning from Material.C in framework. These methods are unused in all inheriting classes and should not be overwritten.

Definition at line 116 of file StressUpdateBase.h.

116 {}

◆ setEffectiveElasticity()

void CappedMohrCoulombCosseratStressUpdate::setEffectiveElasticity ( const RankFourTensor &  Eijkl)
overrideprotectedvirtual

Sets _Eij and _En and _Cij.

Implements MultiParameterPlasticityStressUpdate.

Definition at line 58 of file CappedMohrCoulombCosseratStressUpdate.C.

59 {
60  _Eij[0][0] = _Eij[1][1] = _Eij[2][2] = _host_E0000;
61  _Eij[0][1] = _Eij[1][0] = _Eij[0][2] = _Eij[2][0] = _Eij[1][2] = _Eij[2][1] = _host_E0011;
62  _En = _Eij[2][2];
63  const Real denom = _Eij[0][0] * (_Eij[0][0] + _Eij[0][1]) - 2 * Utility::pow<2>(_Eij[0][1]);
64  for (unsigned a = 0; a < _num_sp; ++a)
65  {
66  _Cij[a][a] = (_Eij[0][0] + _Eij[0][1]) / denom;
67  for (unsigned b = 0; b < a; ++b)
68  _Cij[a][b] = _Cij[b][a] = -_Eij[0][1] / denom;
69  }
70 }
const Real _host_E0011
E0011 = Lame lambda modulus of the host material.
std::vector< std::vector< Real > > _Cij
_Cij[i, j] * _Eij[j, k] = 1 iff j == k
std::vector< std::vector< Real > > _Eij
E[i, j] in the system of equations to be solved.
const Real _host_E0000
E0000 = Lame lambda + 2 * shear modulus of the host material.
const unsigned _num_sp
Number of stress parameters.

◆ setInelasticStrainIncrementAfterReturn()

void MultiParameterPlasticityStressUpdate::setInelasticStrainIncrementAfterReturn ( const RankTwoTensor &  stress_trial,
Real  gaE,
const yieldAndFlow smoothed_q,
const RankFourTensor &  elasticity_tensor,
const RankTwoTensor &  returned_stress,
RankTwoTensor &  inelastic_strain_increment 
) const
protectedvirtualinherited

Sets inelastic strain increment from the returned configuration This is called after the return-map process has completed successfully in stress_param space, just after finalizeReturnProcess has been called.

Derived classes may override this function

Parameters
stress_trial[in]The trial value of stress
gaE[in]The value of gaE after the return-map process has completed successfully
smoothed_q[in]Holds the current value of yield function and derivatives evaluated at the returned configuration
elasticity_tensor[in]The elasticity tensor
returned_stress[in]The stress after the return-map process
inelastic_strain_increment[out]The inelastic strain increment resulting from this return-map

Reimplemented in TwoParameterPlasticityStressUpdate.

Definition at line 714 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

721 {
722  const std::vector<RankTwoTensor> dsp_dstress = dstress_param_dstress(returned_stress);
723  inelastic_strain_increment = RankTwoTensor();
724  for (unsigned i = 0; i < _num_sp; ++i)
725  inelastic_strain_increment += smoothed_q.dg[i] * dsp_dstress[i];
726  inelastic_strain_increment *= gaE / _En;
727 }
virtual std::vector< RankTwoTensor > dstress_param_dstress(const RankTwoTensor &stress) const =0
d(stress_param[i])/d(stress) at given stress
const unsigned _num_sp
Number of stress parameters.

◆ setIntnlDerivativesV()

void CappedMohrCoulombStressUpdate::setIntnlDerivativesV ( const std::vector< Real > &  trial_stress_params,
const std::vector< Real > &  current_stress_params,
const std::vector< Real > &  intnl,
std::vector< std::vector< Real >> &  dintnl 
) const
overrideprotectedvirtualinherited

Sets the derivatives of internal parameters, based on the trial values of stress_params, their current values, and the current values of the internal parameters.

Derived classes must override this.

Parameters
trial_stress_params[in]The trial stress parameters
current_stress_params[in]The current stress parameters
intnl[in]The current value of the internal parameters
dintnl[out]The derivatives dintnl[i][j] = d(intnl[i])/d(stress_param j)

Implements MultiParameterPlasticityStressUpdate.

Definition at line 755 of file CappedMohrCoulombStressUpdate.C.

759 {
760  // intnl[0] = shear, intnl[1] = tensile
761  const Real smax = current_stress_params[2]; // largest eigenvalue
762  const Real smin = current_stress_params[0]; // smallest eigenvalue
763  const Real trial_smax = trial_stress_params[2]; // largest eigenvalue
764  const Real trial_smin = trial_stress_params[0]; // smallest eigenvalue
765  const Real ga_shear = ((trial_smax - trial_smin) - (smax - smin)) / (_Eij[2][2] - _Eij[0][2]);
766  const std::vector<Real> dga_shear = {
767  1.0 / (_Eij[2][2] - _Eij[0][2]), 0.0, -1.0 / (_Eij[2][2] - _Eij[0][2])};
768  // intnl[0] = intnl_old[0] + ga_shear;
769  for (std::size_t i = 0; i < _num_sp; ++i)
770  dintnl[0][i] = dga_shear[i];
771 
772  const Real sinpsi = std::sin(_psi.value(intnl[0]));
773  const Real dsinpsi_di0 = _psi.derivative(intnl[0]) * std::cos(_psi.value(intnl[0]));
774 
775  const Real prefactor = (_Eij[2][2] + _Eij[0][2]) * sinpsi;
776  const Real dprefactor_di0 = (_Eij[2][2] + _Eij[0][2]) * dsinpsi_di0;
777  // const Real shear_correction = prefactor * ga_shear;
778  std::vector<Real> dshear_correction(_num_sp);
779  for (std::size_t i = 0; i < _num_sp; ++i)
780  dshear_correction[i] = prefactor * dga_shear[i] + dprefactor_di0 * dintnl[0][i] * ga_shear;
781  // const Real ga_tensile = (1 - _poissons_ratio) * ((trial_smax + trial_smin) - (smax + smin) -
782  // shear_correction) /
783  // _Eij[2][2];
784  // intnl[1] = intnl_old[1] + ga_tensile;
785  for (std::size_t i = 0; i < _num_sp; ++i)
786  dintnl[1][i] = -(1.0 - _poissons_ratio) * dshear_correction[i] / _Eij[2][2];
787  dintnl[1][2] += -(1.0 - _poissons_ratio) / _Eij[2][2];
788  dintnl[1][0] += -(1.0 - _poissons_ratio) / _Eij[2][2];
789 }
const TensorMechanicsHardeningModel & _psi
Hardening model for dilation angle.
std::vector< std::vector< Real > > _Eij
E[i, j] in the system of equations to be solved.
virtual Real derivative(Real intnl) const
virtual Real value(Real intnl) const
const unsigned _num_sp
Number of stress parameters.

◆ setIntnlValuesV()

void CappedMohrCoulombStressUpdate::setIntnlValuesV ( const std::vector< Real > &  trial_stress_params,
const std::vector< Real > &  current_stress_params,
const std::vector< Real > &  intnl_old,
std::vector< Real > &  intnl 
) const
overrideprotectedvirtualinherited

Sets the internal parameters based on the trial values of stress_params, their current values, and the old values of the internal parameters.

Derived classes must override this.

Parameters
trial_stress_params[in]The trial stress parameters (eg trial_p and trial_q)
current_stress_params[in]The current stress parameters (eg p and q)
intnl_old[out]Old value of internal parameters
intnl[out]The value of internal parameters to be set

Implements MultiParameterPlasticityStressUpdate.

Definition at line 733 of file CappedMohrCoulombStressUpdate.C.

Referenced by CappedMohrCoulombStressUpdate::initializeVarsV().

737 {
738  // intnl[0] = shear, intnl[1] = tensile
739  const Real smax = current_stress_params[2]; // largest eigenvalue
740  const Real smin = current_stress_params[0]; // smallest eigenvalue
741  const Real trial_smax = trial_stress_params[2]; // largest eigenvalue
742  const Real trial_smin = trial_stress_params[0]; // smallest eigenvalue
743  const Real ga_shear = ((trial_smax - trial_smin) - (smax - smin)) / (_Eij[2][2] - _Eij[0][2]);
744  intnl[0] = intnl_old[0] + ga_shear;
745  const Real sinpsi = std::sin(_psi.value(intnl[0]));
746  const Real prefactor = (_Eij[2][2] + _Eij[0][2]) * sinpsi;
747  const Real shear_correction = prefactor * ga_shear;
748  const Real ga_tensile = (1.0 - _poissons_ratio) *
749  ((trial_smax + trial_smin) - (smax + smin) - shear_correction) /
750  _Eij[2][2];
751  intnl[1] = intnl_old[1] + ga_tensile;
752 }
const TensorMechanicsHardeningModel & _psi
Hardening model for dilation angle.
std::vector< std::vector< Real > > _Eij
E[i, j] in the system of equations to be solved.
virtual Real value(Real intnl) const

◆ setQp()

void StressUpdateBase::setQp ( unsigned int  qp)
inherited

Sets the value of the global variable _qp for inheriting classes.

Definition at line 42 of file StressUpdateBase.C.

43 {
44  _qp = qp;
45 }

◆ setStressAfterReturnV()

void CappedMohrCoulombCosseratStressUpdate::setStressAfterReturnV ( const RankTwoTensor &  stress_trial,
const std::vector< Real > &  stress_params,
Real  gaE,
const std::vector< Real > &  intnl,
const yieldAndFlow smoothed_q,
const RankFourTensor &  Eijkl,
RankTwoTensor &  stress 
) const
overrideprotectedvirtual

Sets stress from the admissible parameters.

This is called after the return-map process has completed successfully in stress_param space, just after finalizeReturnProcess has been called. Derived classes must override this function

Parameters
stress_trial[in]The trial value of stress
stress_params[in]The value of the stress_params after the return-map process has completed successfully
gaE[in]The value of gaE after the return-map process has completed successfully
intnl[in]The value of the internal parameters after the return-map process has completed successfully
smoothed_q[in]Holds the current value of yield function and derivatives evaluated at the returned state
Eijkl[in]The elasticity tensor
stress[out]The returned value of the stress tensor

Reimplemented from CappedMohrCoulombStressUpdate.

Definition at line 73 of file CappedMohrCoulombCosseratStressUpdate.C.

81 {
82  // form the diagonal stress
83  stress = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0);
84  // rotate to the original frame, to give the symmetric part of the stress
85  stress = _eigvecs * stress * (_eigvecs.transpose());
86  // add the non-symmetric parts
87  stress += 0.5 * (stress_trial - stress_trial.transpose());
88 }
RankTwoTensor _eigvecs
Eigenvectors of the trial stress as a RankTwoTensor, in order to rotate the returned stress back to s...

◆ smoothAllQuantities()

MultiParameterPlasticityStressUpdate::yieldAndFlow MultiParameterPlasticityStressUpdate::smoothAllQuantities ( const std::vector< Real > &  stress_params,
const std::vector< Real > &  intnl 
) const
protectedinherited

Calculates all yield functions and derivatives, and then performs the smoothing scheme.

Parameters
stress_params[in]The stress parameters (eg stress_params[0] = stress_zz and stress_params[1] = sqrt(stress_zx^2 + stress_zy^2))
intnl[in]Internal parameters
Returns
The smoothed yield function and derivatives

Definition at line 373 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::lineSearch(), and MultiParameterPlasticityStressUpdate::updateState().

375 {
376  std::vector<yieldAndFlow> all_q(_num_yf, yieldAndFlow(_num_sp, _num_intnl));
377  computeAllQV(stress_params, intnl, all_q);
378 
379  /* This routine holds the key to my smoothing strategy. It
380  * may be proved that this smoothing strategy produces a
381  * yield surface that is both C2 differentiable and convex,
382  * assuming the individual yield functions are C2 and
383  * convex too.
384  * Of course all the derivatives must also be smoothed.
385  * Also, I assume that d(flow potential)/dstress gets smoothed
386  * by the Yield Function (which produces a C2 flow potential).
387  * See the line identified in the loop below.
388  * Only time will tell whether this is a good strategy, but it
389  * works well in all tests so far. Convexity is irrelevant
390  * for the non-associated case, but at least the return-map
391  * problem should always have a unique solution.
392  * For two yield functions+flows, labelled 1 and 2, we
393  * should have
394  * d(g1 - g2) . d(f1 - f2) >= 0
395  * If not then the return-map problem for even the
396  * multi-surface plasticity with no smoothing won't have a
397  * unique solution. If the multi-surface plasticity has
398  * a unique solution then the smoothed version defined
399  * below will too.
400  */
401 
402  // res_f is the index that contains the smoothed yieldAndFlow
403  std::size_t res_f = 0;
404 
405  for (std::size_t a = 1; a < all_q.size(); ++a)
406  {
407  if (all_q[res_f].f >= all_q[a].f + _smoothing_tol)
408  // no smoothing is needed: res_f is already indexes the largest yield function
409  continue;
410  else if (all_q[a].f >= all_q[res_f].f + _smoothing_tol)
411  {
412  // no smoothing is needed, and res_f needs to index to all_q[a]
413  res_f = a;
414  continue;
415  }
416  else
417  {
418  // smoothing is required
419  const Real f_diff = all_q[res_f].f - all_q[a].f;
420  const Real ism = ismoother(f_diff);
421  const Real sm = smoother(f_diff);
422  const Real dsm = dsmoother(f_diff);
423  // we want: all_q[res_f].f = 0.5 * all_q[res_f].f + all_q[a].f + _smoothing_tol) + ism,
424  // but we have to do the derivatives first
425  for (unsigned i = 0; i < _num_sp; ++i)
426  {
427  for (unsigned j = 0; j < _num_sp; ++j)
428  all_q[res_f].d2g[i][j] =
429  0.5 * (all_q[res_f].d2g[i][j] + all_q[a].d2g[i][j]) +
430  dsm * (all_q[res_f].df[j] - all_q[a].df[j]) * (all_q[res_f].dg[i] - all_q[a].dg[i]) +
431  sm * (all_q[res_f].d2g[i][j] - all_q[a].d2g[i][j]);
432  for (unsigned j = 0; j < _num_intnl; ++j)
433  all_q[res_f].d2g_di[i][j] = 0.5 * (all_q[res_f].d2g_di[i][j] + all_q[a].d2g_di[i][j]) +
434  dsm * (all_q[res_f].df_di[j] - all_q[a].df_di[j]) *
435  (all_q[res_f].dg[i] - all_q[a].dg[i]) +
436  sm * (all_q[res_f].d2g_di[i][j] - all_q[a].d2g_di[i][j]);
437  }
438  for (unsigned i = 0; i < _num_sp; ++i)
439  {
440  all_q[res_f].df[i] = 0.5 * (all_q[res_f].df[i] + all_q[a].df[i]) +
441  sm * (all_q[res_f].df[i] - all_q[a].df[i]);
442  // whether the following (smoothing g with f's smoother) is a good strategy remains to be
443  // seen...
444  all_q[res_f].dg[i] = 0.5 * (all_q[res_f].dg[i] + all_q[a].dg[i]) +
445  sm * (all_q[res_f].dg[i] - all_q[a].dg[i]);
446  }
447  for (unsigned i = 0; i < _num_intnl; ++i)
448  all_q[res_f].df_di[i] = 0.5 * (all_q[res_f].df_di[i] + all_q[a].df_di[i]) +
449  sm * (all_q[res_f].df_di[i] - all_q[a].df_di[i]);
450  all_q[res_f].f = 0.5 * (all_q[res_f].f + all_q[a].f + _smoothing_tol) + ism;
451  }
452  }
453  return all_q[res_f];
454 }
const unsigned _num_intnl
Number of internal parameters.
Real smoother(Real f_diff) const
Derivative of ismoother.
const unsigned _num_yf
Number of yield functions.
virtual void computeAllQV(const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< yieldAndFlow > &all_q) const =0
Completely fills all_q with correct values.
Real dsmoother(Real f_diff) const
Derivative of smoother.
const Real _smoothing_tol
Smoothing tolerance: edges of the yield surface get smoothed by this amount.
const unsigned _num_sp
Number of stress parameters.
Real ismoother(Real f_diff) const
Smooths yield functions.

◆ smoother()

Real MultiParameterPlasticityStressUpdate::smoother ( Real  f_diff) const
protectedinherited

Derivative of ismoother.

Definition at line 465 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::smoothAllQuantities().

466 {
467  if (std::abs(f_diff) >= _smoothing_tol)
468  return 0.0;
469  return 0.5 * std::sin(f_diff * M_PI * 0.5 / _smoothing_tol);
470 }
const Real _smoothing_tol
Smoothing tolerance: edges of the yield surface get smoothed by this amount.

◆ updateState()

void MultiParameterPlasticityStressUpdate::updateState ( RankTwoTensor &  strain_increment,
RankTwoTensor &  inelastic_strain_increment,
const RankTwoTensor &  rotation_increment,
RankTwoTensor &  stress_new,
const RankTwoTensor &  stress_old,
const RankFourTensor &  elasticity_tensor,
const RankTwoTensor &  elastic_strain_old,
bool  compute_full_tangent_operator,
RankFourTensor &  tangent_operator 
)
overrideprotectedvirtualinherited

Given a strain increment that results in a trial stress, perform some procedure (such as an iterative return-mapping process) to produce an admissible stress, an elastic strain increment and an inelastic strain increment.

If _fe_problem.currentlyComputingJacobian() = true, then updateState also computes d(stress)/d(strain) (or some approximation to it).

This method is called by ComputeMultipleInelasticStress. This method is pure virutal: all inheriting classes must overwrite this method.

Parameters
strain_incrementUpon input: the strain increment. Upon output: the elastic strain increment
inelastic_strain_incrementThe inelastic_strain resulting from the interative procedure
rotation_incrementThe finite-strain rotation increment
stress_newUpon input: the trial stress that results from applying strain_increment as an elastic strain. Upon output: the admissible stress
stress_oldThe old value of stress
elasticity_tensorThe elasticity tensor
compute_full_tangent_operatorThe calling routine would like the full consistent tangent operator to be placed in tangent_operator, if possible. This is irrelevant if _fe_problem.currentlyComputingJacobian() = false
tangent_operatord(stress)/d(strain), or some approximation to it If compute_full_tangent_operator=false, then tangent_operator=elasticity_tensor is an appropriate choice. tangent_operator is only computed if _fe_problem.currentlyComputingJacobian() = true

Implements StressUpdateBase.

Definition at line 144 of file MultiParameterPlasticityStressUpdate.C.

153 {
155 
156  if (_t_step >= 2)
157  _step_one = false;
158 
159  // initially assume an elastic deformation
160  std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _intnl[_qp].begin());
161  _iter[_qp] = 0.0;
162  _linesearch_needed[_qp] = 0.0;
163 
164  computeStressParams(stress_new, _trial_sp);
166 
167  if (yieldF(_yf[_qp]) <= _f_tol)
168  {
170  inelastic_strain_increment.zero();
171  if (_fe_problem.currentlyComputingJacobian())
172  tangent_operator = elasticity_tensor;
173  return;
174  }
175 
176  _stress_trial = stress_new;
177  /* The trial stress must be inadmissible
178  * so we need to return to the yield surface. The following
179  * equations must be satisfied.
180  *
181  * 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, i] * dg/dS[i]
182  * 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, i] * dg/dS[i]
183  * ...
184  * 0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, i] * dg/dS[i]
185  * 0 = rhs[N] = f(S, intnl)
186  *
187  * as well as equations defining intnl parameters as functions of
188  * stress_params, trial_stress_params and intnl_old
189  *
190  * The unknowns are S[0], ..., S[N-1], gaE, and the intnl parameters.
191  * Here gaE = ga * _En (the _En serves to make gaE similar magnitude to S)
192  * I find it convenient to solve the first N+1 equations for p, q and gaE,
193  * while substituting the "intnl parameters" equations into these during the solve process
194  */
195 
196  for (auto & deriv : _dvar_dtrial)
197  deriv.assign(_num_sp, 0.0);
198 
199  preReturnMapV(_trial_sp, stress_new, _intnl_old[_qp], _yf[_qp], elasticity_tensor);
200 
201  setEffectiveElasticity(elasticity_tensor);
202 
203  if (_step_one)
204  std::copy(_definitely_ok_sp.begin(), _definitely_ok_sp.end(), _ok_sp.begin());
205  else
206  computeStressParams(stress_old, _ok_sp);
207  std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _ok_intnl.begin());
208 
209  // Return-map problem: must apply the following changes in stress_params,
210  // and find the returned stress_params and gaE
211  for (unsigned i = 0; i < _num_sp; ++i)
212  _del_stress_params[i] = _trial_sp[i] - _ok_sp[i];
213 
214  Real step_taken = 0.0; // amount of del_stress_params that we've applied and the return-map
215  // problem has succeeded
216  Real step_size = 1.0; // potentially can apply del_stress_params in substeps
217  Real gaE_total = 0.0;
218 
219  // current values of the yield function, derivatives, etc
220  yieldAndFlow smoothed_q;
221 
222  // In the following sub-stepping procedure it is possible that
223  // the last step is an elastic step, and therefore smoothed_q won't
224  // be computed on the last step, so we have to compute it.
225  bool smoothed_q_calculated = false;
226 
227  while (step_taken < 1.0 && step_size >= _min_step_size)
228  {
229  if (1.0 - step_taken < step_size)
230  // prevent over-shoots of substepping
231  step_size = 1.0 - step_taken;
232 
233  // trial variables in terms of admissible variables
234  for (unsigned i = 0; i < _num_sp; ++i)
235  _trial_sp[i] = _ok_sp[i] + step_size * _del_stress_params[i];
236 
237  // initialize variables that are to be found via Newton-Raphson
239  Real gaE = 0.0;
240 
241  // flags indicating failure of Newton-Raphson and line-search
242  int nr_failure = 0;
243  int ls_failure = 0;
244 
245  // NR iterations taken in this substep
246  unsigned step_iter = 0;
247 
248  // The residual-squared for the line-search
249  Real res2 = 0.0;
250 
251  if (step_size < 1.0 && yieldF(_trial_sp, _ok_intnl) <= _f_tol)
252  // This is an elastic step
253  // The "step_size < 1.0" in above condition is for efficiency: we definitely
254  // know that this is a plastic step if step_size = 1.0
255  smoothed_q_calculated = false;
256  else
257  {
258  // this is a plastic step
259 
260  // initialize current_sp, gaE and current_intnl based on the non-smoothed situation
262  // and find the smoothed yield function, flow potential and derivatives
264  smoothed_q_calculated = true;
265  calculateRHS(_trial_sp, _current_sp, gaE, smoothed_q, _rhs);
266  res2 = calculateRes2(_rhs);
267 
268  // Perform a Newton-Raphson with linesearch to get current_sp, gaE, and also smoothed_q
269  while (res2 > _f_tol2 && step_iter < _max_nr_its && nr_failure == 0 && ls_failure == 0)
270  {
271  // solve the linear system and store the answer (the "updates") in rhs
272  nr_failure = nrStep(smoothed_q, _trial_sp, _current_sp, _current_intnl, gaE, _rhs);
273  if (nr_failure != 0)
274  break;
275 
276  // handle precision loss
277  if (precisionLoss(_rhs, _current_sp, gaE))
278  {
280  {
281  Moose::err << "MultiParameterPlasticityStressUpdate: precision-loss in element "
282  << _current_elem->id() << " quadpoint=" << _qp << ". Stress_params =";
283  for (unsigned i = 0; i < _num_sp; ++i)
284  Moose::err << " " << _current_sp[i];
285  Moose::err << " gaE = " << gaE << "\n";
286  }
287  res2 = 0.0;
288  break;
289  }
290 
291  // apply (parts of) the updates, re-calculate smoothed_q, and res2
292  ls_failure = lineSearch(res2,
293  _current_sp,
294  gaE,
295  _trial_sp,
296  smoothed_q,
297  _ok_intnl,
299  _rhs,
300  _linesearch_needed[_qp]);
301  step_iter++;
302  }
303  }
304 
305  if (res2 <= _f_tol2 && step_iter < _max_nr_its && nr_failure == 0 && ls_failure == 0 &&
306  gaE >= 0.0)
307  {
308  // this Newton-Raphson worked fine, or this was an elastic step
309  std::copy(_current_sp.begin(), _current_sp.end(), _ok_sp.begin());
310  gaE_total += gaE;
311  step_taken += step_size;
313  std::copy(_intnl[_qp].begin(), _intnl[_qp].end(), _ok_intnl.begin());
314  // calculate dp/dp_trial, dp/dq_trial, etc, for Jacobian
315  dVardTrial(!smoothed_q_calculated,
316  _trial_sp,
317  _ok_sp,
318  gaE,
319  _ok_intnl,
320  smoothed_q,
321  step_size,
322  compute_full_tangent_operator,
323  _dvar_dtrial);
324  if (static_cast<Real>(step_iter) > _iter[_qp])
325  _iter[_qp] = static_cast<Real>(step_iter);
326  step_size *= 1.1;
327  }
328  else
329  {
330  // Newton-Raphson + line-search process failed
331  step_size *= 0.5;
332  }
333  }
334 
335  if (step_size < _min_step_size)
336  errorHandler("MultiParameterPlasticityStressUpdate: Minimum step-size violated");
337 
338  // success!
339  finalizeReturnProcess(rotation_increment);
340  yieldFunctionValuesV(_ok_sp, _intnl[_qp], _yf[_qp]);
341 
342  if (!smoothed_q_calculated)
343  smoothed_q = smoothAllQuantities(_ok_sp, _intnl[_qp]);
344 
346  _stress_trial, _ok_sp, gaE_total, _intnl[_qp], smoothed_q, elasticity_tensor, stress_new);
347 
349  gaE_total,
350  smoothed_q,
351  elasticity_tensor,
352  stress_new,
353  inelastic_strain_increment);
354 
355  strain_increment = strain_increment - inelastic_strain_increment;
356  _plastic_strain[_qp] = _plastic_strain_old[_qp] + inelastic_strain_increment;
357 
358  if (_fe_problem.currentlyComputingJacobian())
359  // for efficiency, do not compute the tangent operator if not currently computing Jacobian
361  _trial_sp,
362  stress_new,
363  _ok_sp,
364  gaE_total,
365  smoothed_q,
366  elasticity_tensor,
367  compute_full_tangent_operator,
368  _dvar_dtrial,
369  tangent_operator);
370 }
MaterialProperty< std::vector< Real > > & _yf
yield functions
std::vector< Real > _trial_sp
"Trial" value of stress_params that initializes the return-map process This is derived from stress = ...
virtual void setEffectiveElasticity(const RankFourTensor &Eijkl)=0
Sets _Eij and _En and _Cij.
const bool _warn_about_precision_loss
Output a warning message if precision loss is encountered during the return-map process.
virtual void errorHandler(const std::string &message) const
Performs any necessary cleaning-up, then throw MooseException(message)
virtual void initializeVarsV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &stress_params, Real &gaE, std::vector< Real > &intnl) const
Sets (stress_params, intnl) at "good guesses" of the solution to the Return-Map algorithm.
int lineSearch(Real &res2, std::vector< Real > &stress_params, Real &gaE, const std::vector< Real > &trial_stress_params, yieldAndFlow &smoothed_q, const std::vector< Real > &intnl_ok, std::vector< Real > &intnl, std::vector< Real > &rhs, Real &linesearch_needed) const
Performs a line-search to find stress_params and gaE Upon entry:
virtual void consistentTangentOperatorV(const RankTwoTensor &stress_trial, const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, const std::vector< std::vector< Real >> &dvar_dtrial, RankFourTensor &cto)
Calculates the consistent tangent operator.
virtual void setIntnlValuesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &intnl) const =0
Sets the internal parameters based on the trial values of stress_params, their current values...
std::vector< Real > _ok_intnl
The state (ok_sp, ok_intnl) is known to be admissible.
const Real _f_tol2
Square of the yield-function tolerance.
const Real _min_step_size
In order to help the Newton-Raphson procedure, the applied strain increment may be applied in sub-inc...
int nrStep(const yieldAndFlow &smoothed_q, const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, const std::vector< Real > &intnl, Real gaE, std::vector< Real > &rhs) const
Performs a Newton-Raphson step to attempt to zero rhs Upon return, rhs will contain the solution...
void calculateRHS(const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, std::vector< Real > &rhs) const
Calculates the RHS in the following 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j] 0 = rhs[...
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
yieldAndFlow smoothAllQuantities(const std::vector< Real > &stress_params, const std::vector< Real > &intnl) const
Calculates all yield functions and derivatives, and then performs the smoothing scheme.
const Real _f_tol
The yield-function tolerance.
virtual void setInelasticStrainIncrementAfterReturn(const RankTwoTensor &stress_trial, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &elasticity_tensor, const RankTwoTensor &returned_stress, RankTwoTensor &inelastic_strain_increment) const
Sets inelastic strain increment from the returned configuration This is called after the return-map p...
const std::vector< Real > _definitely_ok_sp
An admissible value of stress_params at the initial time.
const MaterialProperty< std::vector< Real > > & _intnl_old
old values of internal parameters
RankTwoTensor _stress_trial
"Trial" value of stress that is set at the beginning of the return-map process.
virtual void computeStressParams(const RankTwoTensor &stress, std::vector< Real > &stress_params) const =0
Computes stress_params, given stress.
Real calculateRes2(const std::vector< Real > &rhs) const
Calculates the residual-squared for the Newton-Raphson + line-search.
virtual void yieldFunctionValuesV(const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< Real > &yf) const =0
Computes the values of the yield functions, given stress_params and intnl parameters.
MaterialProperty< Real > & _iter
Number of Newton-Raphson iterations used in the return-map.
MaterialProperty< RankTwoTensor > & _plastic_strain
plastic strain
std::vector< Real > _current_intnl
The current values of the internal params during the Newton-Raphson.
virtual void finalizeReturnProcess(const RankTwoTensor &rotation_increment)
Derived classes may use this to perform calculations after the return-map process has completed succe...
MaterialProperty< Real > & _linesearch_needed
Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) ...
std::vector< std::vector< Real > > _dvar_dtrial
d({stress_param[i], gaE})/d(trial_stress_param[j])
std::vector< Real > _current_sp
The current values of the stress params during the Newton-Raphson.
bool _step_one
handles case of initial_stress that is inadmissible being supplied
void dVardTrial(bool elastic_only, const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, Real step_size, bool compute_full_tangent_operator, std::vector< std::vector< Real >> &dvar_dtrial) const
Calculates derivatives of the stress_params and gaE with repect to the trial values of the stress_par...
virtual void setStressAfterReturnV(const RankTwoTensor &stress_trial, const std::vector< Real > &stress_params, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, RankTwoTensor &stress) const =0
Sets stress from the admissible parameters.
virtual void preReturnMapV(const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress_trial, const std::vector< Real > &intnl_old, const std::vector< Real > &yf, const RankFourTensor &Eijkl)
Derived classes may employ this function to record stuff or do other computations prior to the return...
bool precisionLoss(const std::vector< Real > &solution, const std::vector< Real > &stress_params, Real gaE) const
Check whether precision loss has occurred.
std::vector< Real > _del_stress_params
_del_stress_params = trial_stress_params - ok_sp This is fixed at the beginning of the return-map pro...
const unsigned _max_nr_its
Maximum number of Newton-Raphson iterations allowed in the return-map process.
virtual void initializeReturnProcess()
Derived classes may use this to perform calculations before any return-map process is performed...
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
Old value of plastic strain.
Real yieldF(const std::vector< Real > &stress_params, const std::vector< Real > &intnl) const
Computes the smoothed yield function.
std::vector< Real > _rhs
0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, i] * dg/dS[i] 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1...
std::vector< Real > _ok_sp
The state (ok_sp, ok_intnl) is known to be admissible, so ok_sp are stress_params that are "OK"...
const unsigned _num_sp
Number of stress parameters.

◆ yieldF() [1/2]

Real MultiParameterPlasticityStressUpdate::yieldF ( const std::vector< Real > &  stress_params,
const std::vector< Real > &  intnl 
) const
protectedinherited

Computes the smoothed yield function.

Parameters
stress_paramsThe stress parameters (eg stress_params[0] = stress_zz and stress_params[1] = sqrt(stress_zx^2 + stress_zy^2))
intnlThe internal parameters (eg intnl[0] is shear, intnl[1] is tensile)
Returns
The smoothed yield function value

Definition at line 616 of file MultiParameterPlasticityStressUpdate.C.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

618 {
619  std::vector<Real> yfs(_num_yf);
620  yieldFunctionValuesV(stress_params, intnl, yfs);
621  return yieldF(yfs);
622 }
const unsigned _num_yf
Number of yield functions.
virtual void yieldFunctionValuesV(const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< Real > &yf) const =0
Computes the values of the yield functions, given stress_params and intnl parameters.
Real yieldF(const std::vector< Real > &stress_params, const std::vector< Real > &intnl) const
Computes the smoothed yield function.

◆ yieldF() [2/2]

Real MultiParameterPlasticityStressUpdate::yieldF ( const std::vector< Real > &  yfs) const
protectedinherited

Computes the smoothed yield function.

Parameters
yfsThe values of the individual yield functions
Returns
The smoothed yield function value

Definition at line 625 of file MultiParameterPlasticityStressUpdate.C.

626 {
627  Real yf = yfs[0];
628  for (std::size_t i = 1; i < yfs.size(); ++i)
629  if (yf >= yfs[i] + _smoothing_tol)
630  // no smoothing is needed, and yf is the biggest yield function
631  continue;
632  else if (yfs[i] >= yf + _smoothing_tol)
633  // no smoothing is needed, and yfs[i] is the biggest yield function
634  yf = yfs[i];
635  else
636  yf = 0.5 * (yf + yfs[i] + _smoothing_tol) + ismoother(yf - yfs[i]);
637  return yf;
638 }
const Real _smoothing_tol
Smoothing tolerance: edges of the yield surface get smoothed by this amount.
Real ismoother(Real f_diff) const
Smooths yield functions.

◆ yieldFunctionValuesV()

void CappedMohrCoulombStressUpdate::yieldFunctionValuesV ( const std::vector< Real > &  stress_params,
const std::vector< Real > &  intnl,
std::vector< Real > &  yf 
) const
overrideprotectedvirtualinherited

Computes the values of the yield functions, given stress_params and intnl parameters.

Derived classes must override this, to provide the values of the yield functions in yf.

Parameters
stress_params[in]The stress parameters
intnl[in]The internal parameters
[out]yfThe yield function values

Implements MultiParameterPlasticityStressUpdate.

Definition at line 122 of file CappedMohrCoulombStressUpdate.C.

125 {
126  // intnl[0] = shear, intnl[1] = tensile
127  const Real ts = _tensile_strength.value(intnl[1]);
128  const Real cs = _compressive_strength.value(intnl[1]);
129  const Real sinphi = std::sin(_phi.value(intnl[0]));
130  const Real cohcos = _cohesion.value(intnl[0]) * std::cos(_phi.value(intnl[0]));
131  const Real smax = stress_params[2]; // largest eigenvalue
132  const Real smid = stress_params[1];
133  const Real smin = stress_params[0]; // smallest eigenvalue
134  yf[0] = smax - ts;
135  yf[1] = smid - ts;
136  yf[2] = smin - ts;
137  yf[3] = -smin - cs;
138  yf[4] = -smid - cs;
139  yf[5] = -smax - cs;
140  yf[6] = 0.5 * (smax - smin) + 0.5 * (smax + smin) * sinphi - cohcos;
141  yf[7] = 0.5 * (smid - smin) + 0.5 * (smid + smin) * sinphi - cohcos;
142  yf[8] = 0.5 * (smax - smid) + 0.5 * (smax + smid) * sinphi - cohcos;
143  yf[9] = 0.5 * (smid - smax) + 0.5 * (smid + smax) * sinphi - cohcos;
144  yf[10] = 0.5 * (smin - smid) + 0.5 * (smin + smid) * sinphi - cohcos;
145  yf[11] = 0.5 * (smin - smax) + 0.5 * (smin + smax) * sinphi - cohcos;
146 }
const TensorMechanicsHardeningModel & _compressive_strength
Hardening model for compressive strength.
const TensorMechanicsHardeningModel & _tensile_strength
Hardening model for tensile strength.
const TensorMechanicsHardeningModel & _phi
Hardening model for friction angle.
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
virtual Real value(Real intnl) const

Member Data Documentation

◆ _base_name

const std::string StressUpdateBase::_base_name
protectedinherited

Name used as a prefix for all material properties related to the stress update model.

Definition at line 122 of file StressUpdateBase.h.

◆ _Cij

std::vector<std::vector<Real> > MultiParameterPlasticityStressUpdate::_Cij
protectedinherited

◆ _cohesion

const TensorMechanicsHardeningModel& CappedMohrCoulombStressUpdate::_cohesion
protectedinherited

◆ _compressive_strength

const TensorMechanicsHardeningModel& CappedMohrCoulombStressUpdate::_compressive_strength
protectedinherited

◆ _definitely_ok_sp

const std::vector<Real> MultiParameterPlasticityStressUpdate::_definitely_ok_sp
protectedinherited

◆ _eigvecs

RankTwoTensor CappedMohrCoulombStressUpdate::_eigvecs
protectedinherited

Eigenvectors of the trial stress as a RankTwoTensor, in order to rotate the returned stress back to stress space.

Definition at line 65 of file CappedMohrCoulombStressUpdate.h.

Referenced by CappedMohrCoulombStressUpdate::consistentTangentOperatorV(), preReturnMapV(), CappedMohrCoulombStressUpdate::preReturnMapV(), setStressAfterReturnV(), and CappedMohrCoulombStressUpdate::setStressAfterReturnV().

◆ _Eij

std::vector<std::vector<Real> > MultiParameterPlasticityStressUpdate::_Eij
protectedinherited

◆ _En

Real MultiParameterPlasticityStressUpdate::_En
protectedinherited

◆ _f_tol

const Real MultiParameterPlasticityStressUpdate::_f_tol
protectedinherited

◆ _f_tol2

const Real MultiParameterPlasticityStressUpdate::_f_tol2
protectedinherited

Square of the yield-function tolerance.

Definition at line 164 of file MultiParameterPlasticityStressUpdate.h.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

◆ _host_E0000

const Real CappedMohrCoulombCosseratStressUpdate::_host_E0000
protected

E0000 = Lame lambda + 2 * shear modulus of the host material.

Definition at line 53 of file CappedMohrCoulombCosseratStressUpdate.h.

Referenced by setEffectiveElasticity().

◆ _host_E0011

const Real CappedMohrCoulombCosseratStressUpdate::_host_E0011
protected

E0011 = Lame lambda modulus of the host material.

Definition at line 50 of file CappedMohrCoulombCosseratStressUpdate.h.

Referenced by setEffectiveElasticity().

◆ _host_poisson

const Real CappedMohrCoulombCosseratStressUpdate::_host_poisson
protected

Poisson's of the host material.

Definition at line 47 of file CappedMohrCoulombCosseratStressUpdate.h.

Referenced by preReturnMapV().

◆ _host_young

const Real CappedMohrCoulombCosseratStressUpdate::_host_young
protected

Young's modulus of the host material.

Definition at line 44 of file CappedMohrCoulombCosseratStressUpdate.h.

◆ _intnl

MaterialProperty<std::vector<Real> >& MultiParameterPlasticityStressUpdate::_intnl
protectedinherited

◆ _intnl_old

const MaterialProperty<std::vector<Real> >& MultiParameterPlasticityStressUpdate::_intnl_old
protectedinherited

◆ _iter

MaterialProperty<Real>& MultiParameterPlasticityStressUpdate::_iter
protectedinherited

Number of Newton-Raphson iterations used in the return-map.

Definition at line 195 of file MultiParameterPlasticityStressUpdate.h.

Referenced by MultiParameterPlasticityStressUpdate::initQpStatefulProperties(), and MultiParameterPlasticityStressUpdate::updateState().

◆ _linesearch_needed

MaterialProperty<Real>& MultiParameterPlasticityStressUpdate::_linesearch_needed
protectedinherited

Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise)

Definition at line 198 of file MultiParameterPlasticityStressUpdate.h.

Referenced by MultiParameterPlasticityStressUpdate::initQpStatefulProperties(), and MultiParameterPlasticityStressUpdate::updateState().

◆ _max_nr_its

const unsigned MultiParameterPlasticityStressUpdate::_max_nr_its
protectedinherited

Maximum number of Newton-Raphson iterations allowed in the return-map process.

Definition at line 152 of file MultiParameterPlasticityStressUpdate.h.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

◆ _min_step_size

const Real MultiParameterPlasticityStressUpdate::_min_step_size
protectedinherited

In order to help the Newton-Raphson procedure, the applied strain increment may be applied in sub-increments of size greater than this value.

Definition at line 171 of file MultiParameterPlasticityStressUpdate.h.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

◆ _num_intnl

const unsigned MultiParameterPlasticityStressUpdate::_num_intnl
protectedinherited

◆ _num_sp

const unsigned MultiParameterPlasticityStressUpdate::_num_sp
protectedinherited

◆ _num_yf

const unsigned MultiParameterPlasticityStressUpdate::_num_yf
protectedinherited

◆ _perfect_guess

const bool CappedMohrCoulombStressUpdate::_perfect_guess
protectedinherited

Whether to provide an estimate of the returned stress, based on perfect plasticity.

Definition at line 52 of file CappedMohrCoulombStressUpdate.h.

Referenced by CappedMohrCoulombStressUpdate::initializeVarsV().

◆ _perform_finite_strain_rotations

const bool MultiParameterPlasticityStressUpdate::_perform_finite_strain_rotations
protectedinherited

◆ _phi

const TensorMechanicsHardeningModel& CappedMohrCoulombStressUpdate::_phi
protectedinherited

◆ _plastic_strain

MaterialProperty<RankTwoTensor>& MultiParameterPlasticityStressUpdate::_plastic_strain
protectedinherited

◆ _plastic_strain_old

const MaterialProperty<RankTwoTensor>& MultiParameterPlasticityStressUpdate::_plastic_strain_old
protectedinherited

◆ _poissons_ratio

Real CappedMohrCoulombStressUpdate::_poissons_ratio
protectedinherited

◆ _psi

const TensorMechanicsHardeningModel& CappedMohrCoulombStressUpdate::_psi
protectedinherited

◆ _shifter

const Real CappedMohrCoulombStressUpdate::_shifter
protectedinherited

When equal-eigenvalues are predicted from the stress initialization routine, shift them by this amount.

This avoids equal-eigenvalue problems, but also accounts for the smoothing of the yield surface

Definition at line 62 of file CappedMohrCoulombStressUpdate.h.

Referenced by CappedMohrCoulombStressUpdate::initializeVarsV().

◆ _smoothing_tol

const Real MultiParameterPlasticityStressUpdate::_smoothing_tol
protectedinherited

◆ _step_one

bool MultiParameterPlasticityStressUpdate::_step_one
protectedinherited

handles case of initial_stress that is inadmissible being supplied

Definition at line 174 of file MultiParameterPlasticityStressUpdate.h.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

◆ _tensile_strength

const TensorMechanicsHardeningModel& CappedMohrCoulombStressUpdate::_tensile_strength
protectedinherited

◆ _tensor_dimensionality

constexpr unsigned MultiParameterPlasticityStressUpdate::_tensor_dimensionality = 3
staticprotectedinherited

◆ _warn_about_precision_loss

const bool MultiParameterPlasticityStressUpdate::_warn_about_precision_loss
protectedinherited

Output a warning message if precision loss is encountered during the return-map process.

Definition at line 177 of file MultiParameterPlasticityStressUpdate.h.

Referenced by MultiParameterPlasticityStressUpdate::updateState().

◆ _yf

MaterialProperty<std::vector<Real> >& MultiParameterPlasticityStressUpdate::_yf
protectedinherited

The documentation for this class was generated from the following files: