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 22 of file TensorMechanicsPlasticMohrCoulombMulti.h.

Member Enumeration Documentation

Constructor & Destructor Documentation

TensorMechanicsPlasticMohrCoulombMulti::TensorMechanicsPlasticMohrCoulombMulti ( const InputParameters &  parameters)

Definition at line 52 of file TensorMechanicsPlasticMohrCoulombMulti.C.

54  : TensorMechanicsPlasticModel(parameters),
55  _cohesion(getUserObject<TensorMechanicsHardeningModel>("cohesion")),
56  _phi(getUserObject<TensorMechanicsHardeningModel>("friction_angle")),
57  _psi(getUserObject<TensorMechanicsHardeningModel>("dilation_angle")),
58  _max_iters(getParam<unsigned int>("max_iterations")),
59  _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
60  _use_custom_returnMap(getParam<bool>("use_custom_returnMap"))
61 {
62  if (_shift < 0)
63  mooseError("Value of 'shift' in TensorMechanicsPlasticMohrCoulombMulti must not be negative\n");
64  if (_shift > _f_tol)
65  _console << "WARNING: value of 'shift' in TensorMechanicsPlasticMohrCoulombMulti is probably "
66  "set too high\n";
67  if (LIBMESH_DIM != 3)
68  mooseError("TensorMechanicsPlasticMohrCoulombMulti is only defined for LIBMESH_DIM=3");
69  MooseRandom::seed(0);
70 }
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

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 238 of file TensorMechanicsPlasticMohrCoulombMulti.C.

244 {
245  act.assign(6, false);
246 
247  if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol && f[3] <= _f_tol && f[4] <= _f_tol &&
248  f[5] <= _f_tol)
249  {
250  returned_stress = stress;
251  return;
252  }
253 
254  Real returned_intnl;
255  std::vector<Real> dpm(6);
256  RankTwoTensor delta_dp;
257  std::vector<Real> yf(6);
258  bool trial_stress_inadmissible;
259  doReturnMap(stress,
260  intnl,
261  Eijkl,
262  0.0,
263  returned_stress,
264  returned_intnl,
265  dpm,
266  delta_dp,
267  yf,
268  trial_stress_inadmissible);
269 
270  for (unsigned i = 0; i < 6; ++i)
271  act[i] = (dpm[i] > 0);
272 }
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.
Real TensorMechanicsPlasticMohrCoulombMulti::cohesion ( const Real  internal_param) const
protectedvirtual

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

Definition at line 275 of file TensorMechanicsPlasticMohrCoulombMulti.C.

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

276 {
277  return _cohesion.value(internal_param);
278 }
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
virtual Real value(Real intnl) const
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 249 of file TensorMechanicsPlasticModel.C.

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

256 {
257  return E_ijkl;
258 }
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 281 of file TensorMechanicsPlasticMohrCoulombMulti.C.

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

282 {
283  return _cohesion.derivative(internal_param);
284 }
virtual Real derivative(Real intnl) const
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
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 134 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by dyieldFunction_dstressV(), and flowPotentialV().

137 {
138  std::vector<Real> eigvals;
139  std::vector<RankTwoTensor> deigvals;
140  stress.dsymmetricEigenvalues(eigvals, deigvals);
141 
142  if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
143  perturbStress(stress, eigvals, deigvals);
144 
145  df.resize(6);
146  df[0] = 0.5 * (deigvals[0] - deigvals[1]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
147  df[1] = 0.5 * (deigvals[1] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
148  df[2] = 0.5 * (deigvals[0] - deigvals[2]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
149  df[3] = 0.5 * (deigvals[2] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
150  df[4] = 0.5 * (deigvals[1] - deigvals[2]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
151  df[5] = 0.5 * (deigvals[2] - deigvals[1]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
152 }
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.
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 128 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::dflowPotential_dintnlV().

130 {
131  return RankTwoTensor();
132 }
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 218 of file TensorMechanicsPlasticMohrCoulombMulti.C.

220 {
221  const Real cos_angle = std::cos(psi(intnl));
222  const Real dsin_angle = cos_angle * dpsi(intnl);
223 
224  std::vector<Real> eigvals;
225  std::vector<RankTwoTensor> deigvals;
226  stress.dsymmetricEigenvalues(eigvals, deigvals);
227 
228  if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
229  perturbStress(stress, eigvals, deigvals);
230 
231  dr_dintnl.resize(6);
232  dr_dintnl[0] = dr_dintnl[1] = 0.5 * (deigvals[0] + deigvals[1]) * dsin_angle;
233  dr_dintnl[2] = dr_dintnl[3] = 0.5 * (deigvals[0] + deigvals[2]) * dsin_angle;
234  dr_dintnl[4] = dr_dintnl[5] = 0.5 * (deigvals[1] + deigvals[2]) * dsin_angle;
235 }
virtual Real psi(const Real internal_param) const
psi 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 dpsi(const Real internal_param) const
d(psi)/d(internal_param) as a function of residual value, rate, and internal_param ...
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.
RankFourTensor TensorMechanicsPlasticModel::dflowPotential_dstress ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited
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 194 of file TensorMechanicsPlasticMohrCoulombMulti.C.

196 {
197  std::vector<RankFourTensor> d2eigvals;
198  stress.d2symmetricEigenvalues(d2eigvals);
199 
200  const Real sinpsi = std::sin(psi(intnl));
201 
202  dr_dstress.resize(6);
203  dr_dstress[0] =
204  0.5 * (d2eigvals[0] - d2eigvals[1]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
205  dr_dstress[1] =
206  0.5 * (d2eigvals[1] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
207  dr_dstress[2] =
208  0.5 * (d2eigvals[0] - d2eigvals[2]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
209  dr_dstress[3] =
210  0.5 * (d2eigvals[2] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
211  dr_dstress[4] =
212  0.5 * (d2eigvals[1] - d2eigvals[2]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
213  dr_dstress[5] =
214  0.5 * (d2eigvals[2] - d2eigvals[1]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
215 }
virtual Real psi(const Real internal_param) const
psi as a function of residual value, rate, and internal_param
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 169 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::dhardPotential_dintnlV().

171 {
172  return 0.0;
173 }
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 175 of file TensorMechanicsPlasticModel.C.

178 {
179  dh_dintnl.resize(numberSurfaces(), dhardPotential_dintnl(stress, intnl));
180 }
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. ...
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 155 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::dhardPotential_dstressV().

157 {
158  return RankTwoTensor();
159 }
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 161 of file TensorMechanicsPlasticModel.C.

164 {
165  dh_dstress.assign(numberSurfaces(), dhardPotential_dstress(stress, intnl));
166 }
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.
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 353 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by activeConstraints(), and returnMap().

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

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

294 {
295  return _phi.derivative(internal_param);
296 }
virtual Real derivative(Real intnl) const
const TensorMechanicsHardeningModel & _phi
Hardening model for phi.
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 305 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by dflowPotential_dintnlV().

306 {
307  return _psi.derivative(internal_param);
308 }
virtual Real derivative(Real intnl) const
const TensorMechanicsHardeningModel & _psi
Hardening model for psi.
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 87 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::dyieldFunction_dintnlV().

89 {
90  return 0.0;
91 }
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 163 of file TensorMechanicsPlasticMohrCoulombMulti.C.

166 {
167  std::vector<Real> eigvals;
168  stress.symmetricEigenvalues(eigvals);
169  eigvals[0] += _shift;
170  eigvals[2] -= _shift;
171 
172  const Real sin_angle = std::sin(phi(intnl));
173  const Real cos_angle = std::cos(phi(intnl));
174  const Real dsin_angle = cos_angle * dphi(intnl);
175  const Real dcos_angle = -sin_angle * dphi(intnl);
176  const Real dcohcos = dcohesion(intnl) * cos_angle + cohesion(intnl) * dcos_angle;
177 
178  df_dintnl.resize(6);
179  df_dintnl[0] = df_dintnl[1] = 0.5 * (eigvals[0] + eigvals[1]) * dsin_angle - dcohcos;
180  df_dintnl[2] = df_dintnl[3] = 0.5 * (eigvals[0] + eigvals[2]) * dsin_angle - dcohcos;
181  df_dintnl[4] = df_dintnl[5] = 0.5 * (eigvals[1] + eigvals[2]) * dsin_angle - dcohcos;
182 }
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 ...
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 ...
virtual Real cohesion(const Real internal_param) const
cohesion as a function of residual value, rate, and internal_param
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 72 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::dyieldFunction_dstressV().

74 {
75  return RankTwoTensor();
76 }
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 155 of file TensorMechanicsPlasticMohrCoulombMulti.C.

157 {
158  const Real sinphi = std::sin(phi(intnl));
159  df_dsig(stress, sinphi, df_dstress);
160 }
virtual Real phi(const Real internal_param) const
phi as a function of residual value, rate, and internal_param
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)...
void TensorMechanicsPlasticModel::execute ( )
inherited

Definition at line 42 of file TensorMechanicsPlasticModel.C.

43 {
44 }
void TensorMechanicsPlasticModel::finalize ( )
inherited

Definition at line 47 of file TensorMechanicsPlasticModel.C.

48 {
49 }
RankTwoTensor TensorMechanicsPlasticModel::flowPotential ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited
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 185 of file TensorMechanicsPlasticMohrCoulombMulti.C.

188 {
189  const Real sinpsi = std::sin(psi(intnl));
190  df_dsig(stress, sinpsi, r);
191 }
virtual Real psi(const Real internal_param) const
psi as a function of residual value, rate, and internal_param
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)...
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 142 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::hardPotentialV().

143 {
144  return -1.0;
145 }
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 147 of file TensorMechanicsPlasticModel.C.

150 {
151  h.assign(numberSurfaces(), hardPotential(stress, intnl));
152 }
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.
void TensorMechanicsPlasticModel::initialize ( )
inherited

Definition at line 37 of file TensorMechanicsPlasticModel.C.

38 {
39 }
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 1132 of file TensorMechanicsPlasticMohrCoulombMulti.C.

1135 {
1136  mooseAssert(
1137  yf.size() == 6,
1138  "TensorMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires yf to be size 6, but your is "
1139  << yf.size());
1140  mooseAssert(
1141  dpm.size() == 6,
1142  "TensorMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires dpm to be size 6, but your is "
1143  << dpm.size());
1144  for (unsigned i = 0; i < 6; ++i)
1145  if (!TensorMechanicsPlasticModel::KuhnTuckerSingleSurface(yf[i], dpm[i], ep_plastic_tolerance))
1146  return false;
1147  return true;
1148 }
bool KuhnTuckerSingleSurface(Real yf, Real dpm, Real dpm_tol) const
Returns true if the Kuhn-Tucker conditions for the single surface are satisfied.
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 243 of file TensorMechanicsPlasticModel.C.

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

244 {
245  return (dpm == 0 && yf <= _f_tol) || (dpm > -dpm_tol && yf <= _f_tol && yf >= -_f_tol);
246 }
const Real _f_tol
Tolerance on yield function.
std::string TensorMechanicsPlasticMohrCoulombMulti::modelName ( ) const
overridevirtual

Implements TensorMechanicsPlasticModel.

Definition at line 311 of file TensorMechanicsPlasticMohrCoulombMulti.C.

312 {
313  return "MohrCoulombMulti";
314 }
unsigned int TensorMechanicsPlasticMohrCoulombMulti::numberSurfaces ( ) const
overridevirtual

The number of yield surfaces for this plasticity model.

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 73 of file TensorMechanicsPlasticMohrCoulombMulti.C.

74 {
75  return 6;
76 }
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 114 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by df_dsig(), and dflowPotential_dintnlV().

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

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

Definition at line 287 of file TensorMechanicsPlasticMohrCoulombMulti.C.

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

288 {
289  return _phi.value(internal_param);
290 }
const TensorMechanicsHardeningModel & _phi
Hardening model for phi.
virtual Real value(Real intnl) const
Real TensorMechanicsPlasticMohrCoulombMulti::psi ( const Real  internal_param) const
protectedvirtual

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

Definition at line 299 of file TensorMechanicsPlasticMohrCoulombMulti.C.

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

300 {
301  return _psi.value(internal_param);
302 }
const TensorMechanicsHardeningModel & _psi
Hardening model for psi.
virtual Real value(Real intnl) const
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 898 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap().

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

Referenced by doReturnMap().

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

327 {
329  return TensorMechanicsPlasticModel::returnMap(trial_stress,
330  intnl_old,
331  E_ijkl,
332  ep_plastic_tolerance,
333  returned_stress,
334  returned_intnl,
335  dpm,
336  delta_dp,
337  yf,
338  trial_stress_inadmissible);
339 
340  return doReturnMap(trial_stress,
341  intnl_old,
342  E_ijkl,
343  ep_plastic_tolerance,
344  returned_stress,
345  returned_intnl,
346  dpm,
347  delta_dp,
348  yf,
349  trial_stress_inadmissible);
350 }
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.
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 bool _use_custom_returnMap
Whether to use the custom return-map algorithm.
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 793 of file TensorMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap().

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

Referenced by doReturnMap().

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

211 {
212  return false;
213 }
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 1151 of file TensorMechanicsPlasticMohrCoulombMulti.C.

1152 {
1153  return _use_custom_returnMap;
1154 }
const bool _use_custom_returnMap
Whether to use the custom return-map algorithm.
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, TensorMechanicsPlasticDruckerPrager, TensorMechanicsPlasticIsotropicSD, TensorMechanicsPlasticJ2, TensorMechanicsPlasticOrthotropic, TensorMechanicsPlasticWeakPlaneShear, TensorMechanicsPlasticWeakPlaneTensile, TensorMechanicsPlasticMohrCoulomb, TensorMechanicsPlasticTensile, TensorMechanicsPlasticDruckerPragerHyperbolic, TensorMechanicsPlasticMeanCap, TensorMechanicsPlasticSimpleTester, and TensorMechanicsPlasticWeakPlaneTensileN.

Definition at line 58 of file TensorMechanicsPlasticModel.C.

Referenced by TensorMechanicsPlasticModel::yieldFunctionV().

59 {
60  return 0.0;
61 }
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 96 of file TensorMechanicsPlasticMohrCoulombMulti.C.

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

98 {
99  // Naively it seems a shame to have 6 yield functions active instead of just
100  // 3. But 3 won't do. Eg, think of a loading with eigvals[0]=eigvals[1]=eigvals[2]
101  // Then to return to the yield surface would require 2 positive plastic multipliers
102  // and one negative one. Boo hoo.
103 
104  f.resize(6);
105  f[0] = 0.5 * (e0 - e1) + 0.5 * (e0 + e1) * sinphi - cohcos;
106  f[1] = 0.5 * (e1 - e0) + 0.5 * (e0 + e1) * sinphi - cohcos;
107  f[2] = 0.5 * (e0 - e2) + 0.5 * (e0 + e2) * sinphi - cohcos;
108  f[3] = 0.5 * (e2 - e0) + 0.5 * (e0 + e2) * sinphi - cohcos;
109  f[4] = 0.5 * (e1 - e2) + 0.5 * (e1 + e2) * sinphi - cohcos;
110  f[5] = 0.5 * (e2 - e1) + 0.5 * (e1 + e2) * sinphi - cohcos;
111 }
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 79 of file TensorMechanicsPlasticMohrCoulombMulti.C.

82 {
83  std::vector<Real> eigvals;
84  stress.symmetricEigenvalues(eigvals);
85  eigvals[0] += _shift;
86  eigvals[2] -= _shift;
87 
88  const Real sinphi = std::sin(phi(intnl));
89  const Real cosphi = std::cos(phi(intnl));
90  const Real cohcos = cohesion(intnl) * cosphi;
91 
92  yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, f);
93 }
virtual Real phi(const Real internal_param) const
phi 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 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 cohesion(const Real internal_param) const
cohesion as a function of residual value, rate, and internal_param

Member Data Documentation

const TensorMechanicsHardeningModel& TensorMechanicsPlasticMohrCoulombMulti::_cohesion
private

Hardening model for cohesion.

Definition at line 96 of file TensorMechanicsPlasticMohrCoulombMulti.h.

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

const Real TensorMechanicsPlasticModel::_f_tol
inherited
const Real TensorMechanicsPlasticModel::_ic_tol
inherited

Tolerance on internal constraint.

Definition at line 174 of file TensorMechanicsPlasticModel.h.

const unsigned int TensorMechanicsPlasticMohrCoulombMulti::_max_iters
private

Maximum Newton-Raphison iterations in the custom returnMap algorithm.

Definition at line 105 of file TensorMechanicsPlasticMohrCoulombMulti.h.

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

const TensorMechanicsHardeningModel& TensorMechanicsPlasticMohrCoulombMulti::_phi
private

Hardening model for phi.

Definition at line 99 of file TensorMechanicsPlasticMohrCoulombMulti.h.

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

const TensorMechanicsHardeningModel& TensorMechanicsPlasticMohrCoulombMulti::_psi
private

Hardening model for psi.

Definition at line 102 of file TensorMechanicsPlasticMohrCoulombMulti.h.

Referenced by dpsi(), and psi().

const Real TensorMechanicsPlasticMohrCoulombMulti::_shift
private

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

Definition at line 108 of file TensorMechanicsPlasticMohrCoulombMulti.h.

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

const bool TensorMechanicsPlasticMohrCoulombMulti::_use_custom_returnMap
private

Whether to use the custom return-map algorithm.

Definition at line 111 of file TensorMechanicsPlasticMohrCoulombMulti.h.

Referenced by returnMap(), and useCustomReturnMap().


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