www.mooseframework.org
Public Member Functions | Public Attributes | Protected Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
TensorMechanicsPlasticMohrCoulombMulti Class Reference

FiniteStrainMohrCoulombMulti implements rate-independent non-associative mohr-coulomb with hardening/softening in the finite-strain framework, using planar (non-smoothed) surfaces. More...

#include <TensorMechanicsPlasticMohrCoulombMulti.h>

Inheritance diagram for TensorMechanicsPlasticMohrCoulombMulti:
[legend]

Public Member Functions

 TensorMechanicsPlasticMohrCoulombMulti (const InputParameters &parameters)
 
virtual unsigned int numberSurfaces () const override
 The number of yield surfaces for this plasticity model. More...
 
virtual void yieldFunctionV (const RankTwoTensor &stress, Real intnl, std::vector< Real > &f) const override
 Calculates the yield functions. More...
 
virtual void dyieldFunction_dstressV (const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &df_dstress) const override
 The derivative of yield functions with respect to stress. More...
 
virtual void dyieldFunction_dintnlV (const RankTwoTensor &stress, Real intnl, std::vector< Real > &df_dintnl) const override
 The derivative of yield functions with respect to the internal parameter. More...
 
virtual void flowPotentialV (const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &r) const override
 The flow potentials. More...
 
virtual void dflowPotential_dstressV (const RankTwoTensor &stress, Real intnl, std::vector< RankFourTensor > &dr_dstress) const override
 The derivative of the flow potential with respect to stress. More...
 
virtual void dflowPotential_dintnlV (const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &dr_dintnl) const override
 The derivative of the flow potential with respect to the internal parameter. More...
 
virtual void activeConstraints (const std::vector< Real > &f, const RankTwoTensor &stress, Real intnl, const RankFourTensor &Eijkl, std::vector< bool > &act, RankTwoTensor &returned_stress) const override
 The active yield surfaces, given a vector of yield functions. More...
 
virtual std::string modelName () const override
 
virtual bool useCustomReturnMap () const override
 Returns false. You will want to override this in your derived class if you write a custom returnMap function. More...
 
virtual bool returnMap (const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const override
 Performs a custom return-map. More...
 
void initialize ()
 
void execute ()
 
void finalize ()
 
virtual void hardPotentialV (const RankTwoTensor &stress, Real intnl, std::vector< Real > &h) const
 The hardening potential. More...
 
virtual void dhardPotential_dstressV (const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &dh_dstress) const
 The derivative of the hardening potential with respect to stress. More...
 
virtual void dhardPotential_dintnlV (const RankTwoTensor &stress, Real intnl, std::vector< Real > &dh_dintnl) const
 The derivative of the hardening potential with respect to the internal parameter. More...
 
virtual bool useCustomCTO () const
 Returns false. You will want to override this in your derived class if you write a custom consistent tangent operator function. More...
 
virtual RankFourTensor consistentTangentOperator (const RankTwoTensor &trial_stress, Real intnl_old, const RankTwoTensor &stress, Real intnl, const RankFourTensor &E_ijkl, const std::vector< Real > &cumulative_pm) const
 Calculates a custom consistent tangent operator. More...
 
bool KuhnTuckerSingleSurface (Real yf, Real dpm, Real dpm_tol) const
 Returns true if the Kuhn-Tucker conditions for the single surface are satisfied. More...
 

Public Attributes

const Real _f_tol
 Tolerance on yield function. More...
 
const Real _ic_tol
 Tolerance on internal constraint. More...
 

Protected Member Functions

virtual Real cohesion (const Real internal_param) const
 cohesion as a function of residual value, rate, and internal_param More...
 
virtual Real dcohesion (const Real internal_param) const
 d(cohesion)/d(internal_param) as a function of residual value, rate, and internal_param More...
 
virtual Real phi (const Real internal_param) const
 phi as a function of residual value, rate, and internal_param More...
 
virtual Real dphi (const Real internal_param) const
 d(phi)/d(internal_param) as a function of residual value, rate, and internal_param More...
 
virtual Real psi (const Real internal_param) const
 psi as a function of residual value, rate, and internal_param More...
 
virtual Real dpsi (const Real internal_param) const
 d(psi)/d(internal_param) as a function of residual value, rate, and internal_param More...
 
virtual Real yieldFunction (const RankTwoTensor &stress, Real intnl) const
 The following functions are what you should override when building single-plasticity models. More...
 
virtual RankTwoTensor dyieldFunction_dstress (const RankTwoTensor &stress, Real intnl) const
 The derivative of yield function with respect to stress. More...
 
virtual Real dyieldFunction_dintnl (const RankTwoTensor &stress, Real intnl) const
 The derivative of yield function with respect to the internal parameter. More...
 
virtual RankTwoTensor flowPotential (const RankTwoTensor &stress, Real intnl) const
 The flow potential. More...
 
virtual RankFourTensor dflowPotential_dstress (const RankTwoTensor &stress, Real intnl) const
 The derivative of the flow potential with respect to stress. More...
 
virtual RankTwoTensor dflowPotential_dintnl (const RankTwoTensor &stress, Real intnl) const
 The derivative of the flow potential with respect to the internal parameter. More...
 
virtual Real hardPotential (const RankTwoTensor &stress, Real intnl) const
 The hardening potential. More...
 
virtual RankTwoTensor dhardPotential_dstress (const RankTwoTensor &stress, Real intnl) const
 The derivative of the hardening potential with respect to stress. More...
 
virtual Real dhardPotential_dintnl (const RankTwoTensor &stress, Real intnl) const
 The derivative of the hardening potential with respect to the internal parameter. More...
 

Private Types

enum  return_type {
  tip110100 = 0, tip010101 = 1, edge010100 = 2, edge000101 = 3,
  plane000100 = 4
}
 

Private Member Functions

void yieldFunctionEigvals (Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector< Real > &f) const
 Calculates the yield functions given the eigenvalues of stress. More...
 
void df_dsig (const RankTwoTensor &stress, Real sin_angle, std::vector< RankTwoTensor > &df) const
 this is exactly dyieldFunction_dstress, or flowPotential, depending on whether sin_angle = sin(phi), or sin_angle = sin(psi), respectively More...
 
void perturbStress (const RankTwoTensor &stress, std::vector< Real > &eigvals, std::vector< RankTwoTensor > &deigvals) const
 perturbs the stress tensor in the case of almost-equal eigenvalues. More...
 
bool KuhnTuckerOK (const std::vector< Real > &yf, const std::vector< Real > &dpm, Real ep_plastic_tolerance) const
 Returns true if the Kuhn-Tucker conditions are satisfied. More...
 
bool doReturnMap (const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const
 See doco for returnMap function. More...
 
bool returnTip (const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
 Tries to return-map to the MC tip using the THREE directions given in n, and THREE dpm values are returned. More...
 
bool returnPlane (const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
 Tries to return-map to the MC plane using the n[3] direction The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk. More...
 
bool returnEdge000101 (const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, Real mag_E, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
 Tries to return-map to the MC edge using the n[4] and n[6] directions The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk. More...
 
bool returnEdge010100 (const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, Real mag_E, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
 Tries to return-map to the MC edge using the n[1] and n[3] directions The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk. More...
 

Private Attributes

const TensorMechanicsHardeningModel_cohesion
 Hardening model for cohesion. More...
 
const TensorMechanicsHardeningModel_phi
 Hardening model for phi. More...
 
const TensorMechanicsHardeningModel_psi
 Hardening model for psi. More...
 
const unsigned int _max_iters
 Maximum Newton-Raphison iterations in the custom returnMap algorithm. More...
 
const Real _shift
 yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalues More...
 
const bool _use_custom_returnMap
 Whether to use the custom return-map algorithm. More...
 

Detailed Description

FiniteStrainMohrCoulombMulti implements rate-independent non-associative mohr-coulomb with hardening/softening in the finite-strain framework, using planar (non-smoothed) surfaces.

Definition at line 25 of file TensorMechanicsPlasticMohrCoulombMulti.h.

Member Enumeration Documentation

◆ return_type

Constructor & Destructor Documentation

◆ TensorMechanicsPlasticMohrCoulombMulti()

TensorMechanicsPlasticMohrCoulombMulti::TensorMechanicsPlasticMohrCoulombMulti ( const InputParameters &  parameters)

Definition at line 58 of file TensorMechanicsPlasticMohrCoulombMulti.C.

60  : TensorMechanicsPlasticModel(parameters),
61  _cohesion(getUserObject<TensorMechanicsHardeningModel>("cohesion")),
62  _phi(getUserObject<TensorMechanicsHardeningModel>("friction_angle")),
63  _psi(getUserObject<TensorMechanicsHardeningModel>("dilation_angle")),
64  _max_iters(getParam<unsigned int>("max_iterations")),
65  _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
66  _use_custom_returnMap(getParam<bool>("use_custom_returnMap"))
67 {
68  if (_shift < 0)
69  mooseError("Value of 'shift' in TensorMechanicsPlasticMohrCoulombMulti must not be negative\n");
70  if (_shift > _f_tol)
71  _console << "WARNING: value of 'shift' in TensorMechanicsPlasticMohrCoulombMulti is probably "
72  "set too high\n";
73  if (LIBMESH_DIM != 3)
74  mooseError("TensorMechanicsPlasticMohrCoulombMulti is only defined for LIBMESH_DIM=3");
75  MooseRandom::seed(0);
76 }
const unsigned int _max_iters
Maximum Newton-Raphison iterations in the custom returnMap algorithm.
TensorMechanicsPlasticModel(const InputParameters &parameters)
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
const TensorMechanicsHardeningModel & _phi
Hardening model for phi.
const TensorMechanicsHardeningModel & _psi
Hardening model for psi.
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
const Real _f_tol
Tolerance on yield function.
const bool _use_custom_returnMap
Whether to use the custom return-map algorithm.

Member Function Documentation

◆ activeConstraints()

void TensorMechanicsPlasticMohrCoulombMulti::activeConstraints ( const std::vector< Real > &  f,
const RankTwoTensor &  stress,
Real  intnl,
const RankFourTensor &  Eijkl,
std::vector< bool > &  act,
RankTwoTensor &  returned_stress 
) const
overridevirtual

The active yield surfaces, given a vector of yield functions.

This is used by FiniteStrainMultiPlasticity to determine the initial set of active constraints at the trial (stress, intnl) configuration. It is up to you (the coder) to determine how accurate you want the returned_stress to be. Currently it is only used by FiniteStrainMultiPlasticity to estimate a good starting value for the Newton-Rahson procedure, so currently it may not need to be super perfect.

Parameters
fvalues of the yield functions
stressstress tensor
intnlinternal parameter
Eijklelasticity tensor (stress = Eijkl*strain)
[out]actact[i] = true if the i_th yield function is active
[out]returned_stressApproximate value of the returned stress

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 244 of file TensorMechanicsPlasticMohrCoulombMulti.C.

250 {
251  act.assign(6, false);
252 
253  if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol && f[3] <= _f_tol && f[4] <= _f_tol &&
254  f[5] <= _f_tol)
255  {
256  returned_stress = stress;
257  return;
258  }
259 
260  Real returned_intnl;
261  std::vector<Real> dpm(6);
262  RankTwoTensor delta_dp;
263  std::vector<Real> yf(6);
264  bool trial_stress_inadmissible;
265  doReturnMap(stress,
266  intnl,
267  Eijkl,
268  0.0,
269  returned_stress,
270  returned_intnl,
271  dpm,
272  delta_dp,
273  yf,
274  trial_stress_inadmissible);
275 
276  for (unsigned i = 0; i < 6; ++i)
277  act[i] = (dpm[i] > 0);
278 }
bool doReturnMap(const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const
See doco for returnMap function.
const Real _f_tol
Tolerance on yield function.

◆ cohesion()

Real TensorMechanicsPlasticMohrCoulombMulti::cohesion ( const Real  internal_param) const
protectedvirtual

cohesion as a function of residual value, rate, and internal_param

Definition at line 281 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap(), dyieldFunction_dintnlV(), returnEdge000101(), returnEdge010100(), returnPlane(), returnTip(), and yieldFunctionV().

282 {
283  return _cohesion.value(internal_param);
284 }
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
virtual Real value(Real intnl) const

◆ consistentTangentOperator()

RankFourTensor TensorMechanicsPlasticModel::consistentTangentOperator ( const RankTwoTensor &  trial_stress,
Real  intnl_old,
const RankTwoTensor &  stress,
Real  intnl,
const RankFourTensor &  E_ijkl,
const std::vector< Real > &  cumulative_pm 
) const
virtualinherited

Calculates a custom consistent tangent operator.

You may choose to over-ride this in your derived TensorMechanicsPlasticXXXX class.

(Note, if you over-ride returnMap, you will probably want to override consistentTangentOpertor too, otherwise it will default to E_ijkl.)

Parameters
stress_oldtrial stress before returning
intnl_oldinternal parameter before returning
stresscurrent returned stress state
intnlinternal parameter
E_ijklelasticity tensor
cumulative_pmthe cumulative plastic multipliers
Returns
the consistent tangent operator: E_ijkl if not over-ridden

Reimplemented in TensorMechanicsPlasticTensileMulti, TensorMechanicsPlasticDruckerPragerHyperbolic, TensorMechanicsPlasticMeanCapTC, and TensorMechanicsPlasticJ2.

Definition at line 253 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticJ2::consistentTangentOperator(), TensorMechanicsPlasticDruckerPragerHyperbolic::consistentTangentOperator(), TensorMechanicsPlasticMeanCapTC::consistentTangentOperator(), and TensorMechanicsPlasticTensileMulti::consistentTangentOperator().

260 {
261  return E_ijkl;
262 }

◆ dcohesion()

Real TensorMechanicsPlasticMohrCoulombMulti::dcohesion ( const Real  internal_param) const
protectedvirtual

d(cohesion)/d(internal_param) as a function of residual value, rate, and internal_param

Definition at line 287 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by dyieldFunction_dintnlV(), returnEdge000101(), returnEdge010100(), returnPlane(), and returnTip().

288 {
289  return _cohesion.derivative(internal_param);
290 }
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
virtual Real derivative(Real intnl) const

◆ df_dsig()

void TensorMechanicsPlasticMohrCoulombMulti::df_dsig ( const RankTwoTensor &  stress,
Real  sin_angle,
std::vector< RankTwoTensor > &  df 
) const
private

this is exactly dyieldFunction_dstress, or flowPotential, depending on whether sin_angle = sin(phi), or sin_angle = sin(psi), respectively

Definition at line 140 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by dyieldFunction_dstressV(), and flowPotentialV().

143 {
144  std::vector<Real> eigvals;
145  std::vector<RankTwoTensor> deigvals;
146  stress.dsymmetricEigenvalues(eigvals, deigvals);
147 
148  if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
149  perturbStress(stress, eigvals, deigvals);
150 
151  df.resize(6);
152  df[0] = 0.5 * (deigvals[0] - deigvals[1]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
153  df[1] = 0.5 * (deigvals[1] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
154  df[2] = 0.5 * (deigvals[0] - deigvals[2]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
155  df[3] = 0.5 * (deigvals[2] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
156  df[4] = 0.5 * (deigvals[1] - deigvals[2]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
157  df[5] = 0.5 * (deigvals[2] - deigvals[1]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
158 }
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
void perturbStress(const RankTwoTensor &stress, std::vector< Real > &eigvals, std::vector< RankTwoTensor > &deigvals) const
perturbs the stress tensor in the case of almost-equal eigenvalues.

◆ dflowPotential_dintnl()

RankTwoTensor TensorMechanicsPlasticModel::dflowPotential_dintnl ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited

The derivative of the flow potential with respect to the internal parameter.

Parameters
stressthe stress at which to calculate the flow potential
intnlinternal parameter
Returns
dr_dintnl(i, j) = dr(i, j)/dintnl

Reimplemented in TensorMechanicsPlasticMeanCapTC, TensorMechanicsPlasticDruckerPrager, TensorMechanicsPlasticJ2, TensorMechanicsPlasticWeakPlaneShear, TensorMechanicsPlasticWeakPlaneTensile, TensorMechanicsPlasticMohrCoulomb, TensorMechanicsPlasticTensile, TensorMechanicsPlasticMeanCap, TensorMechanicsPlasticSimpleTester, and TensorMechanicsPlasticWeakPlaneTensileN.

Definition at line 132 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::dflowPotential_dintnlV().

134 {
135  return RankTwoTensor();
136 }

◆ dflowPotential_dintnlV()

void TensorMechanicsPlasticMohrCoulombMulti::dflowPotential_dintnlV ( const RankTwoTensor &  stress,
Real  intnl,
std::vector< RankTwoTensor > &  dr_dintnl 
) const
overridevirtual

The derivative of the flow potential with respect to the internal parameter.

Parameters
stressthe stress at which to calculate the flow potential
intnlinternal parameter
[out]dr_dintnldr_dintnl[alpha](i, j) = dr[alpha](i, j)/dintnl

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 224 of file TensorMechanicsPlasticMohrCoulombMulti.C.

226 {
227  const Real cos_angle = std::cos(psi(intnl));
228  const Real dsin_angle = cos_angle * dpsi(intnl);
229 
230  std::vector<Real> eigvals;
231  std::vector<RankTwoTensor> deigvals;
232  stress.dsymmetricEigenvalues(eigvals, deigvals);
233 
234  if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
235  perturbStress(stress, eigvals, deigvals);
236 
237  dr_dintnl.resize(6);
238  dr_dintnl[0] = dr_dintnl[1] = 0.5 * (deigvals[0] + deigvals[1]) * dsin_angle;
239  dr_dintnl[2] = dr_dintnl[3] = 0.5 * (deigvals[0] + deigvals[2]) * dsin_angle;
240  dr_dintnl[4] = dr_dintnl[5] = 0.5 * (deigvals[1] + deigvals[2]) * dsin_angle;
241 }
virtual Real dpsi(const Real internal_param) const
d(psi)/d(internal_param) as a function of residual value, rate, and internal_param ...
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
void perturbStress(const RankTwoTensor &stress, std::vector< Real > &eigvals, std::vector< RankTwoTensor > &deigvals) const
perturbs the stress tensor in the case of almost-equal eigenvalues.
virtual Real psi(const Real internal_param) const
psi as a function of residual value, rate, and internal_param

◆ dflowPotential_dstress()

RankFourTensor TensorMechanicsPlasticModel::dflowPotential_dstress ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited

◆ dflowPotential_dstressV()

void TensorMechanicsPlasticMohrCoulombMulti::dflowPotential_dstressV ( const RankTwoTensor &  stress,
Real  intnl,
std::vector< RankFourTensor > &  dr_dstress 
) const
overridevirtual

The derivative of the flow potential with respect to stress.

Parameters
stressthe stress at which to calculate the flow potential
intnlinternal parameter
[out]dr_dstressdr_dstress[alpha](i, j, k, l) = dr[alpha](i, j)/dstress(k, l)

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 200 of file TensorMechanicsPlasticMohrCoulombMulti.C.

202 {
203  std::vector<RankFourTensor> d2eigvals;
204  stress.d2symmetricEigenvalues(d2eigvals);
205 
206  const Real sinpsi = std::sin(psi(intnl));
207 
208  dr_dstress.resize(6);
209  dr_dstress[0] =
210  0.5 * (d2eigvals[0] - d2eigvals[1]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
211  dr_dstress[1] =
212  0.5 * (d2eigvals[1] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
213  dr_dstress[2] =
214  0.5 * (d2eigvals[0] - d2eigvals[2]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
215  dr_dstress[3] =
216  0.5 * (d2eigvals[2] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
217  dr_dstress[4] =
218  0.5 * (d2eigvals[1] - d2eigvals[2]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
219  dr_dstress[5] =
220  0.5 * (d2eigvals[2] - d2eigvals[1]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
221 }
virtual Real psi(const Real internal_param) const
psi as a function of residual value, rate, and internal_param

◆ dhardPotential_dintnl()

Real TensorMechanicsPlasticModel::dhardPotential_dintnl ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited

The derivative of the hardening potential with respect to the internal parameter.

Parameters
stressthe stress at which to calculate the hardening potentials
intnlinternal parameter
Returns
the derivative

Reimplemented in TensorMechanicsPlasticMeanCapTC.

Definition at line 173 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::dhardPotential_dintnlV().

175 {
176  return 0.0;
177 }

◆ dhardPotential_dintnlV()

void TensorMechanicsPlasticModel::dhardPotential_dintnlV ( const RankTwoTensor &  stress,
Real  intnl,
std::vector< Real > &  dh_dintnl 
) const
virtualinherited

The derivative of the hardening potential with respect to the internal parameter.

Parameters
stressthe stress at which to calculate the hardening potentials
intnlinternal parameter
[out]dh_dintnldh_dintnl[alpha] = dh[alpha]/dintnl

Definition at line 179 of file TensorMechanicsPlasticModel.C.

182 {
183  dh_dintnl.resize(numberSurfaces(), dhardPotential_dintnl(stress, intnl));
184 }
virtual unsigned int numberSurfaces() const
The number of yield surfaces for this plasticity model.
virtual Real dhardPotential_dintnl(const RankTwoTensor &stress, Real intnl) const
The derivative of the hardening potential with respect to the internal parameter. ...

◆ dhardPotential_dstress()

RankTwoTensor TensorMechanicsPlasticModel::dhardPotential_dstress ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited

The derivative of the hardening potential with respect to stress.

Parameters
stressthe stress at which to calculate the hardening potentials
intnlinternal parameter
Returns
dh_dstress(i, j) = dh/dstress(i, j)

Reimplemented in TensorMechanicsPlasticMeanCapTC.

Definition at line 159 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::dhardPotential_dstressV().

161 {
162  return RankTwoTensor();
163 }

◆ dhardPotential_dstressV()

void TensorMechanicsPlasticModel::dhardPotential_dstressV ( const RankTwoTensor &  stress,
Real  intnl,
std::vector< RankTwoTensor > &  dh_dstress 
) const
virtualinherited

The derivative of the hardening potential with respect to stress.

Parameters
stressthe stress at which to calculate the hardening potentials
intnlinternal parameter
[out]dh_dstressdh_dstress[alpha](i, j) = dh[alpha]/dstress(i, j)

Definition at line 165 of file TensorMechanicsPlasticModel.C.

168 {
169  dh_dstress.assign(numberSurfaces(), dhardPotential_dstress(stress, intnl));
170 }
virtual unsigned int numberSurfaces() const
The number of yield surfaces for this plasticity model.
virtual RankTwoTensor dhardPotential_dstress(const RankTwoTensor &stress, Real intnl) const
The derivative of the hardening potential with respect to stress.

◆ doReturnMap()

bool TensorMechanicsPlasticMohrCoulombMulti::doReturnMap ( const RankTwoTensor &  trial_stress,
Real  intnl_old,
const RankFourTensor &  E_ijkl,
Real  ep_plastic_tolerance,
RankTwoTensor &  returned_stress,
Real &  returned_intnl,
std::vector< Real > &  dpm,
RankTwoTensor &  delta_dp,
std::vector< Real > &  yf,
bool &  trial_stress_inadmissible 
) const
private

See doco for returnMap function.

The interface is identical to this one. This one can be called internally regardless of the value of _use_custom_returnMap

Definition at line 359 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by activeConstraints(), and returnMap().

369 {
370  mooseAssert(dpm.size() == 6,
371  "TensorMechanicsPlasticMohrCoulombMulti size of dpm should be 6 but it is "
372  << dpm.size());
373 
374  std::vector<Real> eigvals;
375  RankTwoTensor eigvecs;
376  trial_stress.symmetricEigenvaluesEigenvectors(eigvals, eigvecs);
377  eigvals[0] += _shift;
378  eigvals[2] -= _shift;
379 
380  Real sinphi = std::sin(phi(intnl_old));
381  Real cosphi = std::cos(phi(intnl_old));
382  Real coh = cohesion(intnl_old);
383  Real cohcos = coh * cosphi;
384 
385  yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, yf);
386 
387  if (yf[0] <= _f_tol && yf[1] <= _f_tol && yf[2] <= _f_tol && yf[3] <= _f_tol && yf[4] <= _f_tol &&
388  yf[5] <= _f_tol)
389  {
390  // purely elastic (trial_stress, intnl_old)
391  trial_stress_inadmissible = false;
392  return true;
393  }
394 
395  trial_stress_inadmissible = true;
396  delta_dp.zero();
397  returned_stress = RankTwoTensor();
398 
399  // these are the normals to the 6 yield surfaces, which are const because of the assumption of no
400  // psi hardening
401  std::vector<RealVectorValue> norm(6);
402  const Real sinpsi = std::sin(psi(intnl_old));
403  const Real oneminus = 0.5 * (1 - sinpsi);
404  const Real oneplus = 0.5 * (1 + sinpsi);
405  norm[0](0) = oneplus;
406  norm[0](1) = -oneminus;
407  norm[0](2) = 0;
408  norm[1](0) = -oneminus;
409  norm[1](1) = oneplus;
410  norm[1](2) = 0;
411  norm[2](0) = oneplus;
412  norm[2](1) = 0;
413  norm[2](2) = -oneminus;
414  norm[3](0) = -oneminus;
415  norm[3](1) = 0;
416  norm[3](2) = oneplus;
417  norm[4](0) = 0;
418  norm[4](1) = oneplus;
419  norm[4](2) = -oneminus;
420  norm[5](0) = 0;
421  norm[5](1) = -oneminus;
422  norm[5](2) = oneplus;
423 
424  // the flow directions are these norm multiplied by Eijkl.
425  // I call the flow directions "n".
426  // In the following I assume that the Eijkl is
427  // for an isotropic situation. Then I don't have to
428  // rotate to the principal-stress frame, and i don't
429  // have to worry about strange off-diagonal things
430  std::vector<RealVectorValue> n(6);
431  for (unsigned ys = 0; ys < 6; ++ys)
432  for (unsigned i = 0; i < 3; ++i)
433  for (unsigned j = 0; j < 3; ++j)
434  n[ys](i) += E_ijkl(i, i, j, j) * norm[ys](j);
435  const Real mag_E = E_ijkl(0, 0, 0, 0);
436 
437  // With non-zero Poisson's ratio and hardening
438  // it is not computationally cheap to know whether
439  // the trial stress will return to the tip, edge,
440  // or plane. The following at least
441  // gives a not-completely-stupid guess
442  // trial_order[0] = type of return to try first
443  // trial_order[1] = type of return to try second
444  // trial_order[2] = type of return to try third
445  // trial_order[3] = type of return to try fourth
446  // trial_order[4] = type of return to try fifth
447  // In the following the "binary" stuff indicates the
448  // deactive (0) and active (1) surfaces, eg
449  // 110100 means that surfaces 0, 1 and 3 are active
450  // and 2, 4 and 5 are deactive
451  const unsigned int number_of_return_paths = 5;
452  std::vector<int> trial_order(number_of_return_paths);
453  if (yf[1] > _f_tol && yf[3] > _f_tol && yf[5] > _f_tol)
454  {
455  trial_order[0] = tip110100;
456  trial_order[1] = edge010100;
457  trial_order[2] = plane000100;
458  trial_order[3] = edge000101;
459  trial_order[4] = tip010101;
460  }
461  else if (yf[1] <= _f_tol && yf[3] > _f_tol && yf[5] > _f_tol)
462  {
463  trial_order[0] = edge000101;
464  trial_order[1] = plane000100;
465  trial_order[2] = tip110100;
466  trial_order[3] = tip010101;
467  trial_order[4] = edge010100;
468  }
469  else if (yf[1] <= _f_tol && yf[3] > _f_tol && yf[5] <= _f_tol)
470  {
471  trial_order[0] = plane000100;
472  trial_order[1] = edge000101;
473  trial_order[2] = edge010100;
474  trial_order[3] = tip110100;
475  trial_order[4] = tip010101;
476  }
477  else
478  {
479  trial_order[0] = edge010100;
480  trial_order[1] = plane000100;
481  trial_order[2] = edge000101;
482  trial_order[3] = tip110100;
483  trial_order[4] = tip010101;
484  }
485 
486  unsigned trial;
487  bool nr_converged = false;
488  bool kt_success = false;
489  std::vector<RealVectorValue> ntip(3);
490  std::vector<Real> dpmtip(3);
491 
492  for (trial = 0; trial < number_of_return_paths; ++trial)
493  {
494  switch (trial_order[trial])
495  {
496  case tip110100:
497  for (unsigned int i = 0; i < 3; ++i)
498  {
499  ntip[0](i) = n[0](i);
500  ntip[1](i) = n[1](i);
501  ntip[2](i) = n[3](i);
502  }
503  kt_success = returnTip(eigvals,
504  ntip,
505  dpmtip,
506  returned_stress,
507  intnl_old,
508  sinphi,
509  cohcos,
510  0,
511  nr_converged,
512  ep_plastic_tolerance,
513  yf);
514  if (nr_converged && kt_success)
515  {
516  dpm[0] = dpmtip[0];
517  dpm[1] = dpmtip[1];
518  dpm[3] = dpmtip[2];
519  dpm[2] = dpm[4] = dpm[5] = 0;
520  }
521  break;
522 
523  case tip010101:
524  for (unsigned int i = 0; i < 3; ++i)
525  {
526  ntip[0](i) = n[1](i);
527  ntip[1](i) = n[3](i);
528  ntip[2](i) = n[5](i);
529  }
530  kt_success = returnTip(eigvals,
531  ntip,
532  dpmtip,
533  returned_stress,
534  intnl_old,
535  sinphi,
536  cohcos,
537  0,
538  nr_converged,
539  ep_plastic_tolerance,
540  yf);
541  if (nr_converged && kt_success)
542  {
543  dpm[1] = dpmtip[0];
544  dpm[3] = dpmtip[1];
545  dpm[5] = dpmtip[2];
546  dpm[0] = dpm[2] = dpm[4] = 0;
547  }
548  break;
549 
550  case edge000101:
551  kt_success = returnEdge000101(eigvals,
552  n,
553  dpm,
554  returned_stress,
555  intnl_old,
556  sinphi,
557  cohcos,
558  0,
559  mag_E,
560  nr_converged,
561  ep_plastic_tolerance,
562  yf);
563  break;
564 
565  case edge010100:
566  kt_success = returnEdge010100(eigvals,
567  n,
568  dpm,
569  returned_stress,
570  intnl_old,
571  sinphi,
572  cohcos,
573  0,
574  mag_E,
575  nr_converged,
576  ep_plastic_tolerance,
577  yf);
578  break;
579 
580  case plane000100:
581  kt_success = returnPlane(eigvals,
582  n,
583  dpm,
584  returned_stress,
585  intnl_old,
586  sinphi,
587  cohcos,
588  0,
589  nr_converged,
590  ep_plastic_tolerance,
591  yf);
592  break;
593  }
594 
595  if (nr_converged && kt_success)
596  break;
597  }
598 
599  if (trial == number_of_return_paths)
600  {
601  sinphi = std::sin(phi(intnl_old));
602  cosphi = std::cos(phi(intnl_old));
603  coh = cohesion(intnl_old);
604  cohcos = coh * cosphi;
605  yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, yf);
606  Moose::err << "Trial stress = \n";
607  trial_stress.print(Moose::err);
608  Moose::err << "which has eigenvalues = " << eigvals[0] << " " << eigvals[1] << " " << eigvals[2]
609  << "\n";
610  Moose::err << "and yield functions = " << yf[0] << " " << yf[1] << " " << yf[2] << " " << yf[3]
611  << " " << yf[4] << " " << yf[5] << "\n";
612  Moose::err << "Internal parameter = " << intnl_old << "\n";
613  mooseError("TensorMechanicsPlasticMohrCoulombMulti: FAILURE! You probably need to implement a "
614  "line search if your hardening is too severe, or you need to tune your tolerances "
615  "(eg, yield_function_tolerance should be a little smaller than (young "
616  "modulus)*ep_plastic_tolerance).\n");
617  return false;
618  }
619 
620  // success
621 
622  returned_intnl = intnl_old;
623  for (unsigned i = 0; i < 6; ++i)
624  returned_intnl += dpm[i];
625  for (unsigned i = 0; i < 6; ++i)
626  for (unsigned j = 0; j < 3; ++j)
627  delta_dp(j, j) += dpm[i] * norm[i](j);
628  returned_stress = eigvecs * returned_stress * (eigvecs.transpose());
629  delta_dp = eigvecs * delta_dp * (eigvecs.transpose());
630  return true;
631 }
void yieldFunctionEigvals(Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector< Real > &f) const
Calculates the yield functions given the eigenvalues of stress.
virtual Real phi(const Real internal_param) const
phi as a function of residual value, rate, and internal_param
bool returnEdge000101(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, Real mag_E, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC edge using the n[4] and n[6] directions The return value is true if the...
bool returnPlane(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC plane using the n[3] direction The return value is true if the internal...
bool returnEdge010100(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, Real mag_E, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC edge using the n[1] and n[3] directions The return value is true if the...
virtual Real cohesion(const Real internal_param) const
cohesion as a function of residual value, rate, and internal_param
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
bool returnTip(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC tip using the THREE directions given in n, and THREE dpm values are ret...
const Real _f_tol
Tolerance on yield function.
virtual Real psi(const Real internal_param) const
psi as a function of residual value, rate, and internal_param

◆ dphi()

Real TensorMechanicsPlasticMohrCoulombMulti::dphi ( const Real  internal_param) const
protectedvirtual

d(phi)/d(internal_param) as a function of residual value, rate, and internal_param

Definition at line 299 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by dyieldFunction_dintnlV(), returnEdge000101(), returnEdge010100(), returnPlane(), and returnTip().

300 {
301  return _phi.derivative(internal_param);
302 }
const TensorMechanicsHardeningModel & _phi
Hardening model for phi.
virtual Real derivative(Real intnl) const

◆ dpsi()

Real TensorMechanicsPlasticMohrCoulombMulti::dpsi ( const Real  internal_param) const
protectedvirtual

d(psi)/d(internal_param) as a function of residual value, rate, and internal_param

Definition at line 311 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by dflowPotential_dintnlV().

312 {
313  return _psi.derivative(internal_param);
314 }
const TensorMechanicsHardeningModel & _psi
Hardening model for psi.
virtual Real derivative(Real intnl) const

◆ dyieldFunction_dintnl()

Real TensorMechanicsPlasticModel::dyieldFunction_dintnl ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited

The derivative of yield function with respect to the internal parameter.

Parameters
stressthe stress at which to calculate the yield function
intnlinternal parameter
Returns
the derivative

Reimplemented in TensorMechanicsPlasticMeanCapTC, TensorMechanicsPlasticDruckerPrager, TensorMechanicsPlasticJ2, TensorMechanicsPlasticWeakPlaneShear, TensorMechanicsPlasticWeakPlaneTensile, TensorMechanicsPlasticMohrCoulomb, TensorMechanicsPlasticTensile, TensorMechanicsPlasticMeanCap, TensorMechanicsPlasticSimpleTester, and TensorMechanicsPlasticWeakPlaneTensileN.

Definition at line 91 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::dyieldFunction_dintnlV().

93 {
94  return 0.0;
95 }

◆ dyieldFunction_dintnlV()

void TensorMechanicsPlasticMohrCoulombMulti::dyieldFunction_dintnlV ( const RankTwoTensor &  stress,
Real  intnl,
std::vector< Real > &  df_dintnl 
) const
overridevirtual

The derivative of yield functions with respect to the internal parameter.

Parameters
stressthe stress at which to calculate the yield function
intnlinternal parameter
[out]df_dintnldf_dintnl[alpha] = df[alpha]/dintnl

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 169 of file TensorMechanicsPlasticMohrCoulombMulti.C.

172 {
173  std::vector<Real> eigvals;
174  stress.symmetricEigenvalues(eigvals);
175  eigvals[0] += _shift;
176  eigvals[2] -= _shift;
177 
178  const Real sin_angle = std::sin(phi(intnl));
179  const Real cos_angle = std::cos(phi(intnl));
180  const Real dsin_angle = cos_angle * dphi(intnl);
181  const Real dcos_angle = -sin_angle * dphi(intnl);
182  const Real dcohcos = dcohesion(intnl) * cos_angle + cohesion(intnl) * dcos_angle;
183 
184  df_dintnl.resize(6);
185  df_dintnl[0] = df_dintnl[1] = 0.5 * (eigvals[0] + eigvals[1]) * dsin_angle - dcohcos;
186  df_dintnl[2] = df_dintnl[3] = 0.5 * (eigvals[0] + eigvals[2]) * dsin_angle - dcohcos;
187  df_dintnl[4] = df_dintnl[5] = 0.5 * (eigvals[1] + eigvals[2]) * dsin_angle - dcohcos;
188 }
virtual Real phi(const Real internal_param) const
phi as a function of residual value, rate, and internal_param
virtual Real dphi(const Real internal_param) const
d(phi)/d(internal_param) as a function of residual value, rate, and internal_param ...
virtual Real cohesion(const Real internal_param) const
cohesion as a function of residual value, rate, and internal_param
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
virtual Real dcohesion(const Real internal_param) const
d(cohesion)/d(internal_param) as a function of residual value, rate, and internal_param ...

◆ dyieldFunction_dstress()

RankTwoTensor TensorMechanicsPlasticModel::dyieldFunction_dstress ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited

The derivative of yield function with respect to stress.

Parameters
stressthe stress at which to calculate the yield function
intnlinternal parameter
Returns
df_dstress(i, j) = dyieldFunction/dstress(i, j)

Reimplemented in TensorMechanicsPlasticMeanCapTC, TensorMechanicsPlasticIsotropicSD, TensorMechanicsPlasticDruckerPrager, TensorMechanicsPlasticJ2, TensorMechanicsPlasticOrthotropic, TensorMechanicsPlasticWeakPlaneShear, TensorMechanicsPlasticWeakPlaneTensile, TensorMechanicsPlasticMohrCoulomb, TensorMechanicsPlasticTensile, TensorMechanicsPlasticMeanCap, TensorMechanicsPlasticSimpleTester, and TensorMechanicsPlasticWeakPlaneTensileN.

Definition at line 76 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::dyieldFunction_dstressV().

78 {
79  return RankTwoTensor();
80 }

◆ dyieldFunction_dstressV()

void TensorMechanicsPlasticMohrCoulombMulti::dyieldFunction_dstressV ( const RankTwoTensor &  stress,
Real  intnl,
std::vector< RankTwoTensor > &  df_dstress 
) const
overridevirtual

The derivative of yield functions with respect to stress.

Parameters
stressthe stress at which to calculate the yield function
intnlinternal parameter
[out]df_dstressdf_dstress[alpha](i, j) = dyieldFunction[alpha]/dstress(i, j)

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 161 of file TensorMechanicsPlasticMohrCoulombMulti.C.

163 {
164  const Real sinphi = std::sin(phi(intnl));
165  df_dsig(stress, sinphi, df_dstress);
166 }
void df_dsig(const RankTwoTensor &stress, Real sin_angle, std::vector< RankTwoTensor > &df) const
this is exactly dyieldFunction_dstress, or flowPotential, depending on whether sin_angle = sin(phi)...
virtual Real phi(const Real internal_param) const
phi as a function of residual value, rate, and internal_param

◆ execute()

void TensorMechanicsPlasticModel::execute ( )
inherited

Definition at line 46 of file TensorMechanicsPlasticModel.C.

47 {
48 }

◆ finalize()

void TensorMechanicsPlasticModel::finalize ( )
inherited

Definition at line 51 of file TensorMechanicsPlasticModel.C.

52 {
53 }

◆ flowPotential()

RankTwoTensor TensorMechanicsPlasticModel::flowPotential ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited

◆ flowPotentialV()

void TensorMechanicsPlasticMohrCoulombMulti::flowPotentialV ( const RankTwoTensor &  stress,
Real  intnl,
std::vector< RankTwoTensor > &  r 
) const
overridevirtual

The flow potentials.

Parameters
stressthe stress at which to calculate the flow potential
intnlinternal parameter
[out]rr[alpha] is the flow potential for the "alpha" yield function

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 191 of file TensorMechanicsPlasticMohrCoulombMulti.C.

194 {
195  const Real sinpsi = std::sin(psi(intnl));
196  df_dsig(stress, sinpsi, r);
197 }
void df_dsig(const RankTwoTensor &stress, Real sin_angle, std::vector< RankTwoTensor > &df) const
this is exactly dyieldFunction_dstress, or flowPotential, depending on whether sin_angle = sin(phi)...
virtual Real psi(const Real internal_param) const
psi as a function of residual value, rate, and internal_param

◆ hardPotential()

Real TensorMechanicsPlasticModel::hardPotential ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited

The hardening potential.

Parameters
stressthe stress at which to calculate the hardening potential
intnlinternal parameter
Returns
the hardening potential

Reimplemented in TensorMechanicsPlasticMeanCapTC.

Definition at line 146 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::hardPotentialV().

147 {
148  return -1.0;
149 }

◆ hardPotentialV()

void TensorMechanicsPlasticModel::hardPotentialV ( const RankTwoTensor &  stress,
Real  intnl,
std::vector< Real > &  h 
) const
virtualinherited

The hardening potential.

Parameters
stressthe stress at which to calculate the hardening potential
intnlinternal parameter
[out]hh[alpha] is the hardening potential for the "alpha" yield function

Definition at line 151 of file TensorMechanicsPlasticModel.C.

154 {
155  h.assign(numberSurfaces(), hardPotential(stress, intnl));
156 }
virtual Real hardPotential(const RankTwoTensor &stress, Real intnl) const
The hardening potential.
virtual unsigned int numberSurfaces() const
The number of yield surfaces for this plasticity model.

◆ initialize()

void TensorMechanicsPlasticModel::initialize ( )
inherited

Definition at line 41 of file TensorMechanicsPlasticModel.C.

42 {
43 }

◆ KuhnTuckerOK()

bool TensorMechanicsPlasticMohrCoulombMulti::KuhnTuckerOK ( const std::vector< Real > &  yf,
const std::vector< Real > &  dpm,
Real  ep_plastic_tolerance 
) const
private

Returns true if the Kuhn-Tucker conditions are satisfied.

Parameters
yfThe six yield function values
dpmThe six plastic multipliers
ep_plastic_toleranceThe tolerance on the plastic strain (if dpm>-ep_plastic_tolerance then it is grouped as "non-negative" in the Kuhn-Tucker conditions).

Definition at line 1141 of file TensorMechanicsPlasticMohrCoulombMulti.C.

1144 {
1145  mooseAssert(
1146  yf.size() == 6,
1147  "TensorMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires yf to be size 6, but your is "
1148  << yf.size());
1149  mooseAssert(
1150  dpm.size() == 6,
1151  "TensorMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires dpm to be size 6, but your is "
1152  << dpm.size());
1153  for (unsigned i = 0; i < 6; ++i)
1154  if (!TensorMechanicsPlasticModel::KuhnTuckerSingleSurface(yf[i], dpm[i], ep_plastic_tolerance))
1155  return false;
1156  return true;
1157 }
bool KuhnTuckerSingleSurface(Real yf, Real dpm, Real dpm_tol) const
Returns true if the Kuhn-Tucker conditions for the single surface are satisfied.

◆ KuhnTuckerSingleSurface()

bool TensorMechanicsPlasticModel::KuhnTuckerSingleSurface ( Real  yf,
Real  dpm,
Real  dpm_tol 
) const
inherited

Returns true if the Kuhn-Tucker conditions for the single surface are satisfied.

Parameters
yfYield function value
dpmplastic multiplier
dpm_toltolerance on plastic multiplier: viz dpm>-dpm_tol means "dpm is non-negative"

Definition at line 247 of file TensorMechanicsPlasticModel.C.

Referenced by KuhnTuckerOK(), TensorMechanicsPlasticTensileMulti::KuhnTuckerOK(), and TensorMechanicsPlasticModel::returnMap().

248 {
249  return (dpm == 0 && yf <= _f_tol) || (dpm > -dpm_tol && yf <= _f_tol && yf >= -_f_tol);
250 }
const Real _f_tol
Tolerance on yield function.

◆ modelName()

std::string TensorMechanicsPlasticMohrCoulombMulti::modelName ( ) const
overridevirtual

Implements TensorMechanicsPlasticModel.

Definition at line 317 of file TensorMechanicsPlasticMohrCoulombMulti.C.

318 {
319  return "MohrCoulombMulti";
320 }

◆ numberSurfaces()

unsigned int TensorMechanicsPlasticMohrCoulombMulti::numberSurfaces ( ) const
overridevirtual

The number of yield surfaces for this plasticity model.

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 79 of file TensorMechanicsPlasticMohrCoulombMulti.C.

80 {
81  return 6;
82 }

◆ perturbStress()

void TensorMechanicsPlasticMohrCoulombMulti::perturbStress ( const RankTwoTensor &  stress,
std::vector< Real > &  eigvals,
std::vector< RankTwoTensor > &  deigvals 
) const
private

perturbs the stress tensor in the case of almost-equal eigenvalues.

Note that, upon entry, this algorithm assumes that eigvals are the eigvenvalues of stress

Parameters
stressinput stress
[in,out]eigvalseigenvalues after perturbing.
[in,out]deigvalsd(eigenvalues)/dstress_ij after perturbing.

Definition at line 120 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by df_dsig(), and dflowPotential_dintnlV().

123 {
124  Real small_perturbation;
125  RankTwoTensor shifted_stress = stress;
126  while (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
127  {
128  for (unsigned i = 0; i < 3; ++i)
129  for (unsigned j = 0; j <= i; ++j)
130  {
131  small_perturbation = 0.1 * _shift * 2 * (MooseRandom::rand() - 0.5);
132  shifted_stress(i, j) += small_perturbation;
133  shifted_stress(j, i) += small_perturbation;
134  }
135  shifted_stress.dsymmetricEigenvalues(eigvals, deigvals);
136  }
137 }
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...

◆ phi()

Real TensorMechanicsPlasticMohrCoulombMulti::phi ( const Real  internal_param) const
protectedvirtual

phi as a function of residual value, rate, and internal_param

Definition at line 293 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap(), dyieldFunction_dintnlV(), dyieldFunction_dstressV(), returnEdge000101(), returnEdge010100(), returnPlane(), returnTip(), and yieldFunctionV().

294 {
295  return _phi.value(internal_param);
296 }
const TensorMechanicsHardeningModel & _phi
Hardening model for phi.
virtual Real value(Real intnl) const

◆ psi()

Real TensorMechanicsPlasticMohrCoulombMulti::psi ( const Real  internal_param) const
protectedvirtual

psi as a function of residual value, rate, and internal_param

Definition at line 305 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by dflowPotential_dintnlV(), dflowPotential_dstressV(), doReturnMap(), and flowPotentialV().

306 {
307  return _psi.value(internal_param);
308 }
const TensorMechanicsHardeningModel & _psi
Hardening model for psi.
virtual Real value(Real intnl) const

◆ returnEdge000101()

bool TensorMechanicsPlasticMohrCoulombMulti::returnEdge000101 ( const std::vector< Real > &  eigvals,
const std::vector< RealVectorValue > &  n,
std::vector< Real > &  dpm,
RankTwoTensor &  returned_stress,
Real  intnl_old,
Real &  sinphi,
Real &  cohcos,
Real  initial_guess,
Real  mag_E,
bool &  nr_converged,
Real  ep_plastic_tolerance,
std::vector< Real > &  yf 
) const
private

Tries to return-map to the MC edge using the n[4] and n[6] directions The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk.

Parameters
eigvalsThe three stress eigenvalues, sorted in ascending order
nThe six return directions, n=E_ijkl*r. Note this algorithm assumes isotropic elasticity, so these are vectors in principal stress space
dpm[out]The six plastic multipliers resulting from the return-map to the edge
returned_stress[out]The returned stress. This will be diagonal, with the return-mapped eigenvalues in the diagonal positions, sorted in ascending order
intnl_oldThe internal parameter at stress=eigvals. This algorithm doesn't form the plastic strain, so you will have to use intnl=intnl_old+sum(dpm) if you need the new internal-parameter value at the returned point.
sinphi[out]The new value of sin(friction angle) after returning.
cohcos[out]The new value of cohesion*cos(friction angle) after returning.
mag_EAn approximate value for the magnitude of the Young's modulus. This is used to set appropriate tolerances in the Newton-Raphson procedure
initial_guessA guess for sum(dpm)
nr_converged[out]Whether the internal Newton-Raphson process converged
ep_plastic_toleranceThe user-set tolerance on the plastic strain.
yf[out]The yield functions after return

Definition at line 907 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap().

919 {
920  // This returns to the Mohr-Coulomb edge
921  // with 000101 being active. This means that
922  // f3=f5=0. Denoting the returned stress by "a"
923  // (in principal stress space), this means that
924  // a0=a1 = (2Ccosphi + a2(1+sinphi))/(sinphi-1)
925  //
926  // Also, a = eigvals - dpm3*n3 - dpm5*n5,
927  // where the dpm are the plastic multipliers
928  // and the n are the flow directions.
929  //
930  // Hence we have 5 equations and 5 unknowns (a,
931  // dpm3 and dpm5). Eliminating all unknowns
932  // except for x = dpm3+dpm5 gives the residual below
933  // (I used mathematica)
934 
935  Real x = initial_guess;
936  sinphi = std::sin(phi(intnl_old + x));
937  Real cosphi = std::cos(phi(intnl_old + x));
938  Real coh = cohesion(intnl_old + x);
939  cohcos = coh * cosphi;
940 
941  const Real numer_const =
942  -n[3](2) * eigvals[0] - n[5](1) * eigvals[0] + n[5](2) * eigvals[0] + n[3](2) * eigvals[1] +
943  n[5](0) * eigvals[1] - n[5](2) * eigvals[1] - n[5](0) * eigvals[2] + n[5](1) * eigvals[2] +
944  n[3](0) * (-eigvals[1] + eigvals[2]) - n[3](1) * (-eigvals[0] + eigvals[2]);
945  const Real numer_coeff1 = 2 * (-n[3](0) + n[3](1) + n[5](0) - n[5](1));
946  const Real numer_coeff2 =
947  n[5](2) * (eigvals[0] - eigvals[1]) + n[3](2) * (-eigvals[0] + eigvals[1]) +
948  n[5](1) * (eigvals[0] + eigvals[2]) + (n[3](0) - n[5](0)) * (eigvals[1] + eigvals[2]) -
949  n[3](1) * (eigvals[0] + eigvals[2]);
950  Real numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
951  const Real denom_const = -n[3](2) * (n[5](0) - n[5](1)) - n[3](1) * (-n[5](0) + n[5](2)) +
952  n[3](0) * (-n[5](1) + n[5](2));
953  const Real denom_coeff = -n[3](2) * (n[5](0) - n[5](1)) - n[3](1) * (n[5](0) + n[5](2)) +
954  n[3](0) * (n[5](1) + n[5](2));
955  Real denom = denom_const + denom_coeff * sinphi;
956  Real residual = -x + numer / denom;
957 
958  Real deriv_phi;
959  Real deriv_coh;
960  Real jacobian;
961  const Real tol = Utility::pow<2>(_f_tol / (mag_E * 10.0));
962  unsigned int iter = 0;
963  do
964  {
965  do
966  {
967  deriv_phi = dphi(intnl_old + dpm[3]);
968  deriv_coh = dcohesion(intnl_old + dpm[3]);
969  jacobian = -1 +
970  (numer_coeff1 * deriv_coh * cosphi - numer_coeff1 * coh * sinphi * deriv_phi +
971  numer_coeff2 * cosphi * deriv_phi) /
972  denom -
973  numer * denom_coeff * cosphi * deriv_phi / denom / denom;
974  x -= residual / jacobian;
975  if (iter > _max_iters) // not converging
976  {
977  nr_converged = false;
978  return false;
979  }
980 
981  sinphi = std::sin(phi(intnl_old + x));
982  cosphi = std::cos(phi(intnl_old + x));
983  coh = cohesion(intnl_old + x);
984  cohcos = coh * cosphi;
985  numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
986  denom = denom_const + denom_coeff * sinphi;
987  residual = -x + numer / denom;
988  iter++;
989  } while (residual * residual > tol);
990 
991  // now must ensure that yf[3] and yf[5] are both "zero"
992  const Real dpm3minusdpm5 =
993  (2.0 * (eigvals[0] - eigvals[1]) + x * (n[3](1) - n[3](0) + n[5](1) - n[5](0))) /
994  (n[3](0) - n[3](1) + n[5](1) - n[5](0));
995  dpm[3] = (x + dpm3minusdpm5) / 2.0;
996  dpm[5] = (x - dpm3minusdpm5) / 2.0;
997 
998  for (unsigned i = 0; i < 3; ++i)
999  returned_stress(i, i) = eigvals[i] - dpm[3] * n[3](i) - dpm[5] * n[5](i);
1001  returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
1002  } while (yf[3] * yf[3] > _f_tol * _f_tol && yf[5] * yf[5] > _f_tol * _f_tol);
1003 
1004  // so the NR process converged, but we must
1005  // check Kuhn-Tucker
1006  nr_converged = true;
1007 
1008  if (dpm[3] < -ep_plastic_tolerance || dpm[5] < -ep_plastic_tolerance)
1009  // Kuhn-Tucker failure. This is a potential weak-point: if the user has set _f_tol quite
1010  // large, and ep_plastic_tolerance quite small, the above NR process will quickly converge, but
1011  // the solution may be wrong and violate Kuhn-Tucker.
1012  return false;
1013 
1014  if (yf[0] > _f_tol || yf[1] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol)
1015  // Kuhn-Tucker failure
1016  return false;
1017 
1018  // success
1019  dpm[0] = dpm[1] = dpm[2] = dpm[4] = 0;
1020  return true;
1021 }
const unsigned int _max_iters
Maximum Newton-Raphison iterations in the custom returnMap algorithm.
void yieldFunctionEigvals(Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector< Real > &f) const
Calculates the yield functions given the eigenvalues of stress.
const double tol
Definition: Setup.h:19
virtual Real phi(const Real internal_param) const
phi as a function of residual value, rate, and internal_param
virtual Real dphi(const Real internal_param) const
d(phi)/d(internal_param) as a function of residual value, rate, and internal_param ...
virtual Real cohesion(const Real internal_param) const
cohesion as a function of residual value, rate, and internal_param
virtual Real dcohesion(const Real internal_param) const
d(cohesion)/d(internal_param) as a function of residual value, rate, and internal_param ...
const Real _f_tol
Tolerance on yield function.

◆ returnEdge010100()

bool TensorMechanicsPlasticMohrCoulombMulti::returnEdge010100 ( const std::vector< Real > &  eigvals,
const std::vector< RealVectorValue > &  n,
std::vector< Real > &  dpm,
RankTwoTensor &  returned_stress,
Real  intnl_old,
Real &  sinphi,
Real &  cohcos,
Real  initial_guess,
Real  mag_E,
bool &  nr_converged,
Real  ep_plastic_tolerance,
std::vector< Real > &  yf 
) const
private

Tries to return-map to the MC edge using the n[1] and n[3] directions The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk.

Parameters
eigvalsThe three stress eigenvalues, sorted in ascending order
nThe six return directions, n=E_ijkl*r. Note this algorithm assumes isotropic elasticity, so these are vectors in principal stress space
dpm[out]The six plastic multipliers resulting from the return-map to the edge
returned_stress[out]The returned stress. This will be diagonal, with the return-mapped eigenvalues in the diagonal positions, sorted in ascending order
intnl_oldThe internal parameter at stress=eigvals. This algorithm doesn't form the plastic strain, so you will have to use intnl=intnl_old+sum(dpm) if you need the new internal-parameter value at the returned point.
sinphi[out]The new value of sin(friction angle) after returning.
cohcos[out]The new value of cohesion*cos(friction angle) after returning.
mag_EAn approximate value for the magnitude of the Young's modulus. This is used to set appropriate tolerances in the Newton-Raphson procedure
initial_guessA guess for sum(dpm)
nr_converged[out]Whether the internal Newton-Raphson process converged
ep_plastic_toleranceThe user-set tolerance on the plastic strain.
yf[out]The yield functions after return

Definition at line 1024 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap().

1036 {
1037  // This returns to the Mohr-Coulomb edge
1038  // with 010100 being active. This means that
1039  // f1=f3=0. Denoting the returned stress by "a"
1040  // (in principal stress space), this means that
1041  // a1=a2 = (2Ccosphi + a0(1-sinphi))/(sinphi+1)
1042  //
1043  // Also, a = eigvals - dpm1*n1 - dpm3*n3,
1044  // where the dpm are the plastic multipliers
1045  // and the n are the flow directions.
1046  //
1047  // Hence we have 5 equations and 5 unknowns (a,
1048  // dpm1 and dpm3). Eliminating all unknowns
1049  // except for x = dpm1+dpm3 gives the residual below
1050  // (I used mathematica)
1051 
1052  Real x = initial_guess;
1053  sinphi = std::sin(phi(intnl_old + x));
1054  Real cosphi = std::cos(phi(intnl_old + x));
1055  Real coh = cohesion(intnl_old + x);
1056  cohcos = coh * cosphi;
1057 
1058  const Real numer_const = -n[1](2) * eigvals[0] - n[3](1) * eigvals[0] + n[3](2) * eigvals[0] -
1059  n[1](0) * eigvals[1] + n[1](2) * eigvals[1] + n[3](0) * eigvals[1] -
1060  n[3](2) * eigvals[1] + n[1](0) * eigvals[2] - n[3](0) * eigvals[2] +
1061  n[3](1) * eigvals[2] - n[1](1) * (-eigvals[0] + eigvals[2]);
1062  const Real numer_coeff1 = 2 * (n[1](1) - n[1](2) - n[3](1) + n[3](2));
1063  const Real numer_coeff2 =
1064  n[3](2) * (-eigvals[0] - eigvals[1]) + n[1](2) * (eigvals[0] + eigvals[1]) +
1065  n[3](1) * (eigvals[0] + eigvals[2]) + (n[1](0) - n[3](0)) * (eigvals[1] - eigvals[2]) -
1066  n[1](1) * (eigvals[0] + eigvals[2]);
1067  Real numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
1068  const Real denom_const = -n[1](0) * (n[3](1) - n[3](2)) + n[1](2) * (-n[3](0) + n[3](1)) +
1069  n[1](1) * (-n[3](2) + n[3](0));
1070  const Real denom_coeff =
1071  n[1](0) * (n[3](1) - n[3](2)) + n[1](2) * (n[3](0) + n[3](1)) - n[1](1) * (n[3](0) + n[3](2));
1072  Real denom = denom_const + denom_coeff * sinphi;
1073  Real residual = -x + numer / denom;
1074 
1075  Real deriv_phi;
1076  Real deriv_coh;
1077  Real jacobian;
1078  const Real tol = Utility::pow<2>(_f_tol / (mag_E * 10.0));
1079  unsigned int iter = 0;
1080  do
1081  {
1082  do
1083  {
1084  deriv_phi = dphi(intnl_old + dpm[3]);
1085  deriv_coh = dcohesion(intnl_old + dpm[3]);
1086  jacobian = -1 +
1087  (numer_coeff1 * deriv_coh * cosphi - numer_coeff1 * coh * sinphi * deriv_phi +
1088  numer_coeff2 * cosphi * deriv_phi) /
1089  denom -
1090  numer * denom_coeff * cosphi * deriv_phi / denom / denom;
1091  x -= residual / jacobian;
1092  if (iter > _max_iters) // not converging
1093  {
1094  nr_converged = false;
1095  return false;
1096  }
1097 
1098  sinphi = std::sin(phi(intnl_old + x));
1099  cosphi = std::cos(phi(intnl_old + x));
1100  coh = cohesion(intnl_old + x);
1101  cohcos = coh * cosphi;
1102  numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
1103  denom = denom_const + denom_coeff * sinphi;
1104  residual = -x + numer / denom;
1105  iter++;
1106  } while (residual * residual > tol);
1107 
1108  // now must ensure that yf[1] and yf[3] are both "zero"
1109  Real dpm1minusdpm3 =
1110  (2 * (eigvals[1] - eigvals[2]) + x * (n[1](2) - n[1](1) + n[3](2) - n[3](1))) /
1111  (n[1](1) - n[1](2) + n[3](2) - n[3](1));
1112  dpm[1] = (x + dpm1minusdpm3) / 2.0;
1113  dpm[3] = (x - dpm1minusdpm3) / 2.0;
1114 
1115  for (unsigned i = 0; i < 3; ++i)
1116  returned_stress(i, i) = eigvals[i] - dpm[1] * n[1](i) - dpm[3] * n[3](i);
1118  returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
1119  } while (yf[1] * yf[1] > _f_tol * _f_tol && yf[3] * yf[3] > _f_tol * _f_tol);
1120 
1121  // so the NR process converged, but we must
1122  // check Kuhn-Tucker
1123  nr_converged = true;
1124 
1125  if (dpm[1] < -ep_plastic_tolerance || dpm[3] < -ep_plastic_tolerance)
1126  // Kuhn-Tucker failure. This is a potential weak-point: if the user has set _f_tol quite
1127  // large, and ep_plastic_tolerance quite small, the above NR process will quickly converge, but
1128  // the solution may be wrong and violate Kuhn-Tucker
1129  return false;
1130 
1131  if (yf[0] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol || yf[5] > _f_tol)
1132  // Kuhn-Tucker failure
1133  return false;
1134 
1135  // success
1136  dpm[0] = dpm[2] = dpm[4] = dpm[5] = 0;
1137  return true;
1138 }
const unsigned int _max_iters
Maximum Newton-Raphison iterations in the custom returnMap algorithm.
void yieldFunctionEigvals(Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector< Real > &f) const
Calculates the yield functions given the eigenvalues of stress.
const double tol
Definition: Setup.h:19
virtual Real phi(const Real internal_param) const
phi as a function of residual value, rate, and internal_param
virtual Real dphi(const Real internal_param) const
d(phi)/d(internal_param) as a function of residual value, rate, and internal_param ...
virtual Real cohesion(const Real internal_param) const
cohesion as a function of residual value, rate, and internal_param
virtual Real dcohesion(const Real internal_param) const
d(cohesion)/d(internal_param) as a function of residual value, rate, and internal_param ...
const Real _f_tol
Tolerance on yield function.

◆ returnMap()

bool TensorMechanicsPlasticMohrCoulombMulti::returnMap ( const RankTwoTensor &  trial_stress,
Real  intnl_old,
const RankFourTensor &  E_ijkl,
Real  ep_plastic_tolerance,
RankTwoTensor &  returned_stress,
Real &  returned_intnl,
std::vector< Real > &  dpm,
RankTwoTensor &  delta_dp,
std::vector< Real > &  yf,
bool &  trial_stress_inadmissible 
) const
overridevirtual

Performs a custom return-map.

You may choose to over-ride this in your derived TensorMechanicsPlasticXXXX class, and you may implement the return-map algorithm in any way that suits you. Eg, using a Newton-Raphson approach, or a radial-return, etc. This may also be used as a quick way of ascertaining whether (trial_stress, intnl_old) is in fact admissible.

For over-riding this function, please note the following.

(1) Denoting the return value of the function by "successful_return", the only possible output values should be: (A) trial_stress_inadmissible=false, successful_return=true. That is, (trial_stress, intnl_old) is in fact admissible (in the elastic domain). (B) trial_stress_inadmissible=true, successful_return=false. That is (trial_stress, intnl_old) is inadmissible (outside the yield surface), and you didn't return to the yield surface. (C) trial_stress_inadmissible=true, successful_return=true. That is (trial_stress, intnl_old) is inadmissible (outside the yield surface), but you did return to the yield surface. The default implementation only handles case (A) and (B): it does not attempt to do a return-map algorithm.

(2) you must correctly signal "successful_return" using the return value of this function. Don't assume the calling function will do Kuhn-Tucker checking and so forth!

(3) In cases (A) and (B) you needn't set returned_stress, returned_intnl, delta_dp, or dpm. This is for computational efficiency.

(4) In cases (A) and (B), you MUST place the yield function values at (trial_stress, intnl_old) into yf so the calling function can use this information optimally. You will have already calculated these yield function values, which can be quite expensive, and it's not very optimal for the calling function to have to re-calculate them.

(5) In case (C), you need to set: returned_stress (the returned value of stress) returned_intnl (the returned value of the internal variable) delta_dp (the change in plastic strain) dpm (the plastic multipliers needed to bring about the return) yf (yield function values at the returned configuration)

(Note, if you over-ride returnMap, you will probably want to override consistentTangentOpertor too, otherwise it will default to E_ijkl.)

Parameters
trial_stressThe trial stress
intnl_oldValue of the internal parameter
E_ijklElasticity tensor
ep_plastic_toleranceTolerance defined by the user for the plastic strain
[out]returned_stressIn case (C): lies on the yield surface after returning and produces the correct plastic strain (normality condition). Otherwise: not defined
[out]returned_intnlIn case (C): the value of the internal parameter after returning. Otherwise: not defined
[out]dpmIn case (C): the plastic multipliers needed to bring about the return. Otherwise: not defined
[out]delta_dpIn case (C): The change in plastic strain induced by the return process. Otherwise: not defined
[out]yfIn case (C): the yield function at (returned_stress, returned_intnl). Otherwise: the yield function at (trial_stress, intnl_old)
[out]trial_stress_inadmissibleShould be set to false if the trial_stress is admissible, and true if the trial_stress is inadmissible. This can be used by the calling prorgram
Returns
true if a successful return (or a return-map not needed), false if the trial_stress is inadmissible but the return process failed

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 323 of file TensorMechanicsPlasticMohrCoulombMulti.C.

333 {
335  return TensorMechanicsPlasticModel::returnMap(trial_stress,
336  intnl_old,
337  E_ijkl,
338  ep_plastic_tolerance,
339  returned_stress,
340  returned_intnl,
341  dpm,
342  delta_dp,
343  yf,
344  trial_stress_inadmissible);
345 
346  return doReturnMap(trial_stress,
347  intnl_old,
348  E_ijkl,
349  ep_plastic_tolerance,
350  returned_stress,
351  returned_intnl,
352  dpm,
353  delta_dp,
354  yf,
355  trial_stress_inadmissible);
356 }
bool doReturnMap(const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const
See doco for returnMap function.
virtual bool returnMap(const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const
Performs a custom return-map.
const bool _use_custom_returnMap
Whether to use the custom return-map algorithm.

◆ returnPlane()

bool TensorMechanicsPlasticMohrCoulombMulti::returnPlane ( const std::vector< Real > &  eigvals,
const std::vector< RealVectorValue > &  n,
std::vector< Real > &  dpm,
RankTwoTensor &  returned_stress,
Real  intnl_old,
Real &  sinphi,
Real &  cohcos,
Real  initial_guess,
bool &  nr_converged,
Real  ep_plastic_tolerance,
std::vector< Real > &  yf 
) const
private

Tries to return-map to the MC plane using the n[3] direction The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk.

Parameters
eigvalsThe three stress eigenvalues, sorted in ascending order
nThe six return directions, n=E_ijkl*r. Note this algorithm assumes isotropic elasticity, so these are vectors in principal stress space
dpm[out]The six plastic multipliers resulting from the return-map to the plane
returned_stress[out]The returned stress. This will be diagonal, with the return-mapped eigenvalues in the diagonal positions, sorted in ascending order
intnl_oldThe internal parameter at stress=eigvals. This algorithm doesn't form the plastic strain, so you will have to use intnl=intnl_old+sum(dpm) if you need the new internal-parameter value at the returned point.
sinphi[out]The new value of sin(friction angle) after returning.
cohcos[out]The new value of cohesion*cos(friction angle) after returning.
initial_guessA guess for sum(dpm)
nr_converged[out]Whether the internal Newton-Raphson process converged
ep_plastic_toleranceThe user-set tolerance on the plastic strain.
yf[out]The yield functions after return

Definition at line 802 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap().

813 {
814  // This returns to the Mohr-Coulomb plane using n[3] (ie 000100)
815  //
816  // The returned point is defined by the f[3]=0 and
817  // a = eigvals - dpm[3]*n[3]
818  // where "a" is the returned point and dpm[3] is the plastic multiplier.
819  // This equation is a vector equation in principal stress space.
820  // (Kuhn-Tucker also demands that dpm[3]>=0, but we leave checking
821  // that condition for the end.)
822  // Since f[3]=0, we must have
823  // a[2]*(1+sinphi) + a[0]*(-1+sinphi) - 2*coh*cosphi = 0
824  // which gives dpm[3] as the solution of
825  // alpha*dpm[3] + eigvals[2] - eigvals[0] + beta*sinphi - 2*coh*cosphi = 0
826  // with alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0))*sinphi
827  // beta = eigvals[2] + eigvals[0]
828 
829  mooseAssert(n.size() == 6,
830  "TensorMechanicsPlasticMohrCoulombMulti: Custom plane-return algorithm must be "
831  "supplied with n of size 6, whereas yours is "
832  << n.size());
833  mooseAssert(dpm.size() == 6,
834  "TensorMechanicsPlasticMohrCoulombMulti: Custom plane-return algorithm must be "
835  "supplied with dpm of size 6, whereas yours is "
836  << dpm.size());
837  mooseAssert(yf.size() == 6,
838  "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
839  "supplied with yf of size 6, whereas yours is "
840  << yf.size());
841 
842  dpm[3] = initial_guess;
843  sinphi = std::sin(phi(intnl_old + dpm[3]));
844  Real cosphi = std::cos(phi(intnl_old + dpm[3]));
845  Real coh = cohesion(intnl_old + dpm[3]);
846  cohcos = coh * cosphi;
847 
848  Real alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0)) * sinphi;
849  Real deriv_phi;
850  Real dalpha;
851  const Real beta = eigvals[2] + eigvals[0];
852  Real deriv_coh;
853 
854  Real residual =
855  alpha * dpm[3] + eigvals[2] - eigvals[0] + beta * sinphi - 2.0 * cohcos; // this is 2*yf[3]
856  Real jacobian;
857 
858  const Real f_tol2 = Utility::pow<2>(_f_tol);
859  unsigned int iter = 0;
860  do
861  {
862  deriv_phi = dphi(intnl_old + dpm[3]);
863  dalpha = -(n[3](2) + n[3](0)) * cosphi * deriv_phi;
864  deriv_coh = dcohesion(intnl_old + dpm[3]);
865  jacobian = alpha + dalpha * dpm[3] + beta * cosphi * deriv_phi - 2.0 * deriv_coh * cosphi +
866  2.0 * coh * sinphi * deriv_phi;
867 
868  dpm[3] -= residual / jacobian;
869  if (iter > _max_iters) // not converging
870  {
871  nr_converged = false;
872  return false;
873  }
874 
875  sinphi = std::sin(phi(intnl_old + dpm[3]));
876  cosphi = std::cos(phi(intnl_old + dpm[3]));
877  coh = cohesion(intnl_old + dpm[3]);
878  cohcos = coh * cosphi;
879  alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0)) * sinphi;
880  residual = alpha * dpm[3] + eigvals[2] - eigvals[0] + beta * sinphi - 2.0 * cohcos;
881  iter++;
882  } while (residual * residual > f_tol2);
883 
884  // so the NR process converged, but we must
885  // check Kuhn-Tucker
886  nr_converged = true;
887  if (dpm[3] < -ep_plastic_tolerance)
888  // Kuhn-Tucker failure
889  return false;
890 
891  for (unsigned i = 0; i < 3; ++i)
892  returned_stress(i, i) = eigvals[i] - dpm[3] * n[3](i);
894  returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
895 
896  // by construction abs(yf[3]) = abs(residual/2) < _f_tol/2
897  if (yf[0] > _f_tol || yf[1] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol || yf[5] > _f_tol)
898  // Kuhn-Tucker failure
899  return false;
900 
901  // success!
902  dpm[0] = dpm[1] = dpm[2] = dpm[4] = dpm[5] = 0;
903  return true;
904 }
const unsigned int _max_iters
Maximum Newton-Raphison iterations in the custom returnMap algorithm.
void yieldFunctionEigvals(Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector< Real > &f) const
Calculates the yield functions given the eigenvalues of stress.
virtual Real phi(const Real internal_param) const
phi as a function of residual value, rate, and internal_param
virtual Real dphi(const Real internal_param) const
d(phi)/d(internal_param) as a function of residual value, rate, and internal_param ...
virtual Real cohesion(const Real internal_param) const
cohesion as a function of residual value, rate, and internal_param
virtual Real dcohesion(const Real internal_param) const
d(cohesion)/d(internal_param) as a function of residual value, rate, and internal_param ...
const Real _f_tol
Tolerance on yield function.

◆ returnTip()

bool TensorMechanicsPlasticMohrCoulombMulti::returnTip ( const std::vector< Real > &  eigvals,
const std::vector< RealVectorValue > &  n,
std::vector< Real > &  dpm,
RankTwoTensor &  returned_stress,
Real  intnl_old,
Real &  sinphi,
Real &  cohcos,
Real  initial_guess,
bool &  nr_converged,
Real  ep_plastic_tolerance,
std::vector< Real > &  yf 
) const
private

Tries to return-map to the MC tip using the THREE directions given in n, and THREE dpm values are returned.

Note that you must supply THREE suitale n vectors out of the total of SIX flow directions, and then interpret the THREE dpm values appropriately. The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress, dpm) will be junk.

Parameters
eigvalsThe three stress eigenvalues, sorted in ascending order
nThe three return directions, n=E_ijkl*r. Note this algorithm assumes isotropic elasticity, so these are 3 vectors in principal stress space
dpm[out]The three plastic multipliers resulting from the return-map to the tip.
returned_stress[out]The returned stress. This will be diagonal, with the return-mapped eigenvalues in the diagonal positions, sorted in ascending order
intnl_oldThe internal parameter at stress=eigvals. This algorithm doesn't form the plastic strain, so you will have to use intnl=intnl_old+sum(dpm) if you need the new internal-parameter value at the returned point.
sinphi[out]The new value of sin(friction angle) after returning.
cohcos[out]The new value of cohesion*cos(friction angle) after returning.
initial_guessA guess for sum(dpm)
nr_converged[out]Whether the internal Newton-Raphson process converged
ep_plastic_toleranceThe user-set tolerance on the plastic strain.
yf[out]The yield functions after return

Definition at line 634 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap().

645 {
646  // This returns to the Mohr-Coulomb tip using the THREE directions
647  // given in n, and yields the THREE dpm values. Note that you
648  // must supply THREE suitable n vectors out of the total of SIX
649  // flow directions, and then interpret the THREE dpm values appropriately.
650  //
651  // Eg1. You supply the flow directions n[0], n[1] and n[3] as
652  // the "n" vectors. This is return-to-the-tip via 110100.
653  // Then the three returned dpm values will be dpm[0], dpm[1] and dpm[3].
654 
655  // Eg2. You supply the flow directions n[1], n[3] and n[5] as
656  // the "n" vectors. This is return-to-the-tip via 010101.
657  // Then the three returned dpm values will be dpm[1], dpm[3] and dpm[5].
658 
659  // The returned point is defined by the three yield functions (corresonding
660  // to the three supplied flow directions) all being zero.
661  // that is, returned_stress = diag(cohcot, cohcot, cohcot), where
662  // cohcot = cohesion*cosphi/sinphi
663  // where intnl = intnl_old + dpm[0] + dpm[1] + dpm[2]
664  // The 3 plastic multipliers, dpm, are defiend by the normality condition
665  // eigvals - cohcot = dpm[0]*n[0] + dpm[1]*n[1] + dpm[2]*n[2]
666  // (Kuhn-Tucker demands that all dpm are non-negative, but we leave
667  // that checking for the end.)
668  // (Remember that these "n" vectors and "dpm" values must be interpreted
669  // differently depending on what you pass into this function.)
670  // This is a vector equation with solution (A):
671  // dpm[0] = triple(eigvals - cohcot, n[1], n[2])/trip;
672  // dpm[1] = triple(eigvals - cohcot, n[2], n[0])/trip;
673  // dpm[2] = triple(eigvals - cohcot, n[0], n[1])/trip;
674  // where trip = triple(n[0], n[1], n[2]).
675  // By adding the three components of that solution together
676  // we can get an equation for x = dpm[0] + dpm[1] + dpm[2],
677  // and then our Newton-Raphson only involves one variable (x).
678  // In the following, i specialise to the isotropic situation.
679 
680  mooseAssert(n.size() == 3,
681  "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
682  "supplied with n of size 3, whereas yours is "
683  << n.size());
684  mooseAssert(dpm.size() == 3,
685  "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
686  "supplied with dpm of size 3, whereas yours is "
687  << dpm.size());
688  mooseAssert(yf.size() == 6,
689  "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
690  "supplied with yf of size 6, whereas yours is "
691  << yf.size());
692 
693  Real x = initial_guess;
694  const Real trip = triple_product(n[0], n[1], n[2]);
695  sinphi = std::sin(phi(intnl_old + x));
696  Real cosphi = std::cos(phi(intnl_old + x));
697  Real coh = cohesion(intnl_old + x);
698  cohcos = coh * cosphi;
699  Real cohcot = cohcos / sinphi;
700 
701  if (_cohesion.modelName().compare("Constant") != 0 || _phi.modelName().compare("Constant") != 0)
702  {
703  // Finding x is expensive. Therefore
704  // although x!=0 for Constant Hardening, solution (A)
705  // demonstrates that we don't
706  // actually need to know x to find the dpm for
707  // Constant Hardening.
708  //
709  // However, for nontrivial Hardening, the following
710  // is necessary
711  // cohcot_coeff = [1,1,1].(Cross[n[1], n[2]] + Cross[n[2], n[0]] + Cross[n[0], n[1]])/trip
712  Real cohcot_coeff =
713  (n[0](0) * (n[1](1) - n[1](2) - n[2](1)) + (n[1](2) - n[1](1)) * n[2](0) +
714  (n[1](0) - n[1](2)) * n[2](1) + n[0](2) * (n[1](0) - n[1](1) - n[2](0) + n[2](1)) +
715  n[0](1) * (n[1](2) - n[1](0) + n[2](0) - n[2](2)) +
716  (n[0](0) - n[1](0) + n[1](1)) * n[2](2)) /
717  trip;
718  // eig_term = eigvals.(Cross[n[1], n[2]] + Cross[n[2], n[0]] + Cross[n[0], n[1]])/trip
719  Real eig_term = eigvals[0] *
720  (-n[0](2) * n[1](1) + n[0](1) * n[1](2) + n[0](2) * n[2](1) -
721  n[1](2) * n[2](1) - n[0](1) * n[2](2) + n[1](1) * n[2](2)) /
722  trip;
723  eig_term += eigvals[1] *
724  (n[0](2) * n[1](0) - n[0](0) * n[1](2) - n[0](2) * n[2](0) + n[1](2) * n[2](0) +
725  n[0](0) * n[2](2) - n[1](0) * n[2](2)) /
726  trip;
727  eig_term += eigvals[2] *
728  (n[0](0) * n[1](1) - n[1](1) * n[2](0) + n[0](1) * n[2](0) - n[0](1) * n[1](0) -
729  n[0](0) * n[2](1) + n[1](0) * n[2](1)) /
730  trip;
731  // and finally, the equation we want to solve is:
732  // x - eig_term + cohcot*cohcot_coeff = 0
733  // but i divide by cohcot_coeff so the result has the units of
734  // stress, so using _f_tol as a convergence check is reasonable
735  eig_term /= cohcot_coeff;
736  Real residual = x / cohcot_coeff - eig_term + cohcot;
737  Real jacobian;
738  Real deriv_phi;
739  Real deriv_coh;
740  unsigned int iter = 0;
741  do
742  {
743  deriv_phi = dphi(intnl_old + x);
744  deriv_coh = dcohesion(intnl_old + x);
745  jacobian = 1.0 / cohcot_coeff + deriv_coh * cosphi / sinphi -
746  coh * deriv_phi / Utility::pow<2>(sinphi);
747  x += -residual / jacobian;
748 
749  if (iter > _max_iters) // not converging
750  {
751  nr_converged = false;
752  return false;
753  }
754 
755  sinphi = std::sin(phi(intnl_old + x));
756  cosphi = std::cos(phi(intnl_old + x));
757  coh = cohesion(intnl_old + x);
758  cohcos = coh * cosphi;
759  cohcot = cohcos / sinphi;
760  residual = x / cohcot_coeff - eig_term + cohcot;
761  iter++;
762  } while (residual * residual > _f_tol * _f_tol / 100);
763  }
764 
765  // so the NR process converged, but we must
766  // calculate the individual dpm values and
767  // check Kuhn-Tucker
768  nr_converged = true;
769  if (x < -3 * ep_plastic_tolerance)
770  // obviously at least one of the dpm are < -ep_plastic_tolerance. No point in proceeding. This
771  // is a potential weak-point: if the user has set _f_tol quite large, and ep_plastic_tolerance
772  // quite small, the above NR process will quickly converge, but the solution may be wrong and
773  // violate Kuhn-Tucker.
774  return false;
775 
776  // The following is the solution (A) written above
777  // (dpm[0] = triple(eigvals - cohcot, n[1], n[2])/trip, etc)
778  // in the isotropic situation
779  RealVectorValue v;
780  v(0) = eigvals[0] - cohcot;
781  v(1) = eigvals[1] - cohcot;
782  v(2) = eigvals[2] - cohcot;
783  dpm[0] = triple_product(v, n[1], n[2]) / trip;
784  dpm[1] = triple_product(v, n[2], n[0]) / trip;
785  dpm[2] = triple_product(v, n[0], n[1]) / trip;
786 
787  if (dpm[0] < -ep_plastic_tolerance || dpm[1] < -ep_plastic_tolerance ||
788  dpm[2] < -ep_plastic_tolerance)
789  // Kuhn-Tucker failure. No point in proceeding
790  return false;
791 
792  // Kuhn-Tucker has succeeded: just need returned_stress and yf values
793  // I do not use the dpm to calculate returned_stress, because that
794  // might add error (and computational effort), simply:
795  returned_stress(0, 0) = returned_stress(1, 1) = returned_stress(2, 2) = cohcot;
796  // So by construction the yield functions are all zero
797  yf[0] = yf[1] = yf[2] = yf[3] = yf[4] = yf[5] = 0;
798  return true;
799 }
const unsigned int _max_iters
Maximum Newton-Raphison iterations in the custom returnMap algorithm.
virtual Real phi(const Real internal_param) const
phi as a function of residual value, rate, and internal_param
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
const TensorMechanicsHardeningModel & _phi
Hardening model for phi.
virtual Real dphi(const Real internal_param) const
d(phi)/d(internal_param) as a function of residual value, rate, and internal_param ...
virtual std::string modelName() const =0
virtual Real cohesion(const Real internal_param) const
cohesion as a function of residual value, rate, and internal_param
virtual Real dcohesion(const Real internal_param) const
d(cohesion)/d(internal_param) as a function of residual value, rate, and internal_param ...
const Real _f_tol
Tolerance on yield function.

◆ useCustomCTO()

bool TensorMechanicsPlasticModel::useCustomCTO ( ) const
virtualinherited

Returns false. You will want to override this in your derived class if you write a custom consistent tangent operator function.

Reimplemented in TensorMechanicsPlasticTensileMulti, TensorMechanicsPlasticMeanCapTC, TensorMechanicsPlasticDruckerPragerHyperbolic, and TensorMechanicsPlasticJ2.

Definition at line 214 of file TensorMechanicsPlasticModel.C.

215 {
216  return false;
217 }

◆ useCustomReturnMap()

bool TensorMechanicsPlasticMohrCoulombMulti::useCustomReturnMap ( ) const
overridevirtual

Returns false. You will want to override this in your derived class if you write a custom returnMap function.

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 1160 of file TensorMechanicsPlasticMohrCoulombMulti.C.

1161 {
1162  return _use_custom_returnMap;
1163 }
const bool _use_custom_returnMap
Whether to use the custom return-map algorithm.

◆ yieldFunction()

Real TensorMechanicsPlasticModel::yieldFunction ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited

The following functions are what you should override when building single-plasticity models.

The yield function

Parameters
stressthe stress at which to calculate the yield function
intnlinternal parameter
Returns
the yield function

Reimplemented in TensorMechanicsPlasticMeanCapTC, TensorMechanicsPlasticIsotropicSD, TensorMechanicsPlasticDruckerPrager, TensorMechanicsPlasticJ2, TensorMechanicsPlasticOrthotropic, TensorMechanicsPlasticWeakPlaneShear, TensorMechanicsPlasticWeakPlaneTensile, TensorMechanicsPlasticMohrCoulomb, TensorMechanicsPlasticTensile, TensorMechanicsPlasticDruckerPragerHyperbolic, TensorMechanicsPlasticMeanCap, TensorMechanicsPlasticSimpleTester, and TensorMechanicsPlasticWeakPlaneTensileN.

Definition at line 62 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::yieldFunctionV().

63 {
64  return 0.0;
65 }

◆ yieldFunctionEigvals()

void TensorMechanicsPlasticMohrCoulombMulti::yieldFunctionEigvals ( Real  e0,
Real  e1,
Real  e2,
Real  sinphi,
Real  cohcos,
std::vector< Real > &  f 
) const
private

Calculates the yield functions given the eigenvalues of stress.

Parameters
e0Smallest eigenvalue
e1Middle eigenvalue
e2Largest eigenvalue
sinphisin(friction angle)
cohcoscohesion*cos(friction angle)
[out]fthe yield functions

Definition at line 102 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap(), returnEdge000101(), returnEdge010100(), returnPlane(), and yieldFunctionV().

104 {
105  // Naively it seems a shame to have 6 yield functions active instead of just
106  // 3. But 3 won't do. Eg, think of a loading with eigvals[0]=eigvals[1]=eigvals[2]
107  // Then to return to the yield surface would require 2 positive plastic multipliers
108  // and one negative one. Boo hoo.
109 
110  f.resize(6);
111  f[0] = 0.5 * (e0 - e1) + 0.5 * (e0 + e1) * sinphi - cohcos;
112  f[1] = 0.5 * (e1 - e0) + 0.5 * (e0 + e1) * sinphi - cohcos;
113  f[2] = 0.5 * (e0 - e2) + 0.5 * (e0 + e2) * sinphi - cohcos;
114  f[3] = 0.5 * (e2 - e0) + 0.5 * (e0 + e2) * sinphi - cohcos;
115  f[4] = 0.5 * (e1 - e2) + 0.5 * (e1 + e2) * sinphi - cohcos;
116  f[5] = 0.5 * (e2 - e1) + 0.5 * (e1 + e2) * sinphi - cohcos;
117 }

◆ yieldFunctionV()

void TensorMechanicsPlasticMohrCoulombMulti::yieldFunctionV ( const RankTwoTensor &  stress,
Real  intnl,
std::vector< Real > &  f 
) const
overridevirtual

Calculates the yield functions.

Note that for single-surface plasticity you don't want to override this - override the private yieldFunction below

Parameters
stressthe stress at which to calculate the yield function
intnlinternal parameter
[out]fthe yield functions

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 85 of file TensorMechanicsPlasticMohrCoulombMulti.C.

88 {
89  std::vector<Real> eigvals;
90  stress.symmetricEigenvalues(eigvals);
91  eigvals[0] += _shift;
92  eigvals[2] -= _shift;
93 
94  const Real sinphi = std::sin(phi(intnl));
95  const Real cosphi = std::cos(phi(intnl));
96  const Real cohcos = cohesion(intnl) * cosphi;
97 
98  yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, f);
99 }
void yieldFunctionEigvals(Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector< Real > &f) const
Calculates the yield functions given the eigenvalues of stress.
virtual Real phi(const Real internal_param) const
phi as a function of residual value, rate, and internal_param
virtual Real cohesion(const Real internal_param) const
cohesion as a function of residual value, rate, and internal_param
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...

Member Data Documentation

◆ _cohesion

const TensorMechanicsHardeningModel& TensorMechanicsPlasticMohrCoulombMulti::_cohesion
private

Hardening model for cohesion.

Definition at line 99 of file TensorMechanicsPlasticMohrCoulombMulti.h.

Referenced by cohesion(), dcohesion(), and returnTip().

◆ _f_tol

const Real TensorMechanicsPlasticModel::_f_tol
inherited

◆ _ic_tol

const Real TensorMechanicsPlasticModel::_ic_tol
inherited

Tolerance on internal constraint.

Definition at line 177 of file TensorMechanicsPlasticModel.h.

◆ _max_iters

const unsigned int TensorMechanicsPlasticMohrCoulombMulti::_max_iters
private

Maximum Newton-Raphison iterations in the custom returnMap algorithm.

Definition at line 108 of file TensorMechanicsPlasticMohrCoulombMulti.h.

Referenced by returnEdge000101(), returnEdge010100(), returnPlane(), and returnTip().

◆ _phi

const TensorMechanicsHardeningModel& TensorMechanicsPlasticMohrCoulombMulti::_phi
private

Hardening model for phi.

Definition at line 102 of file TensorMechanicsPlasticMohrCoulombMulti.h.

Referenced by dphi(), phi(), and returnTip().

◆ _psi

const TensorMechanicsHardeningModel& TensorMechanicsPlasticMohrCoulombMulti::_psi
private

Hardening model for psi.

Definition at line 105 of file TensorMechanicsPlasticMohrCoulombMulti.h.

Referenced by dpsi(), and psi().

◆ _shift

const Real TensorMechanicsPlasticMohrCoulombMulti::_shift
private

yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalues

Definition at line 111 of file TensorMechanicsPlasticMohrCoulombMulti.h.

Referenced by df_dsig(), dflowPotential_dintnlV(), doReturnMap(), dyieldFunction_dintnlV(), perturbStress(), TensorMechanicsPlasticMohrCoulombMulti(), and yieldFunctionV().

◆ _use_custom_returnMap

const bool TensorMechanicsPlasticMohrCoulombMulti::_use_custom_returnMap
private

Whether to use the custom return-map algorithm.

Definition at line 114 of file TensorMechanicsPlasticMohrCoulombMulti.h.

Referenced by returnMap(), and useCustomReturnMap().


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