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

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

#include <TensorMechanicsPlasticTensileMulti.h>

Inheritance diagram for TensorMechanicsPlasticTensileMulti:
[legend]

Public Member Functions

 TensorMechanicsPlasticTensileMulti (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 useCustomCTO () const override
 Returns false. You will want to override this in your derived class if you write a custom consistent tangent operator 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...
 
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 override
 Calculates a custom consistent tangent operator. 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...
 
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 tensile_strength (const Real internal_param) const
 tensile strength as a function of residual value, rate, and internal_param More...
 
virtual Real dtensile_strength (const Real internal_param) const
 d(tensile strength)/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 { tip = 0, edge = 1, plane = 2 }
 

Private Member Functions

Real dot (const std::vector< Real > &a, const std::vector< Real > &b) const
 dot product of two 3-dimensional vectors More...
 
Real triple (const std::vector< Real > &a, const std::vector< Real > &b, const std::vector< Real > &c) const
 triple product of three 3-dimensional vectors More...
 
bool returnTip (const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real initial_guess) const
 Tries to return-map to the Tensile tip. More...
 
bool returnEdge (const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real initial_guess) const
 Tries to return-map to the Tensile edge. More...
 
bool returnPlane (const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real initial_guess) const
 Tries to return-map to the Tensile plane The return value is true if the internal Newton-Raphson process has converged, otherwise it is false. More...
 
bool KuhnTuckerOK (const RankTwoTensor &returned_diagonal_stress, const std::vector< Real > &dpm, Real str, Real ep_plastic_tolerance) const
 Returns true if the Kuhn-Tucker conditions are satisfied. More...
 
virtual 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
 Just like returnMap, but a protected interface that definitely uses the algorithm, since returnMap itself does not use the algorithm if _use_returnMap=false. More...
 

Private Attributes

const TensorMechanicsHardeningModel_strength
 
const unsigned int _max_iters
 maximum iterations allowed in the custom return-map 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...
 
const bool _use_custom_cto
 Whether to use the custom consistent tangent operator calculation. More...
 

Detailed Description

FiniteStrainTensileMulti implements rate-independent associative tensile failure with hardening/softening in the finite-strain framework, using planar (non-smoothed) surfaces.

Definition at line 25 of file TensorMechanicsPlasticTensileMulti.h.

Member Enumeration Documentation

◆ return_type

Constructor & Destructor Documentation

◆ TensorMechanicsPlasticTensileMulti()

TensorMechanicsPlasticTensileMulti::TensorMechanicsPlasticTensileMulti ( const InputParameters &  parameters)

Definition at line 50 of file TensorMechanicsPlasticTensileMulti.C.

52  : TensorMechanicsPlasticModel(parameters),
53  _strength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
54  _max_iters(getParam<unsigned int>("max_iterations")),
55  _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
56  _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
57  _use_custom_cto(getParam<bool>("use_custom_cto"))
58 {
59  if (_shift < 0)
60  mooseError("Value of 'shift' in TensorMechanicsPlasticTensileMulti must not be negative\n");
61  if (_shift > _f_tol)
62  _console << "WARNING: value of 'shift' in TensorMechanicsPlasticTensileMulti is probably set "
63  "too high\n";
64  if (LIBMESH_DIM != 3)
65  mooseError("TensorMechanicsPlasticTensileMulti is only defined for LIBMESH_DIM=3");
66  MooseRandom::seed(0);
67 }
TensorMechanicsPlasticModel(const InputParameters &parameters)
const bool _use_custom_returnMap
Whether to use the custom return-map algorithm.
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
const unsigned int _max_iters
maximum iterations allowed in the custom return-map algorithm
const Real _f_tol
Tolerance on yield function.
const bool _use_custom_cto
Whether to use the custom consistent tangent operator calculation.
const TensorMechanicsHardeningModel & _strength

Member Function Documentation

◆ activeConstraints()

void TensorMechanicsPlasticTensileMulti::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 159 of file TensorMechanicsPlasticTensileMulti.C.

165 {
166  act.assign(3, false);
167 
168  if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol)
169  {
170  returned_stress = stress;
171  return;
172  }
173 
174  Real returned_intnl;
175  std::vector<Real> dpm(3);
176  RankTwoTensor delta_dp;
177  std::vector<Real> yf(3);
178  bool trial_stress_inadmissible;
179  doReturnMap(stress,
180  intnl,
181  Eijkl,
182  0.0,
183  returned_stress,
184  returned_intnl,
185  dpm,
186  delta_dp,
187  yf,
188  trial_stress_inadmissible);
189 
190  for (unsigned i = 0; i < 3; ++i)
191  act[i] = (dpm[i] > 0);
192 }
virtual 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
Just like returnMap, but a protected interface that definitely uses the algorithm, since returnMap itself does not use the algorithm if _use_returnMap=false.
const Real _f_tol
Tolerance on yield function.

◆ consistentTangentOperator()

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

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 from TensorMechanicsPlasticModel.

Definition at line 614 of file TensorMechanicsPlasticTensileMulti.C.

621 {
622  if (!_use_custom_cto)
624  trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
625 
626  mooseAssert(cumulative_pm.size() == 3,
627  "TensorMechanicsPlasticTensileMulti size of cumulative_pm should be 3 but it is "
628  << cumulative_pm.size());
629 
630  if (cumulative_pm[2] <= 0) // All cumulative_pm are non-positive, so this is admissible
631  return E_ijkl;
632 
633  // Need the eigenvalues at the returned configuration
634  std::vector<Real> eigvals;
635  stress.symmetricEigenvalues(eigvals);
636 
637  // need to rotate to and from principal stress space
638  // using the eigenvectors of the trial configuration
639  // (not the returned configuration).
640  std::vector<Real> trial_eigvals;
641  RankTwoTensor trial_eigvecs;
642  trial_stress.symmetricEigenvaluesEigenvectors(trial_eigvals, trial_eigvecs);
643 
644  // The returnMap will have returned to the Tip, Edge or
645  // Plane. The consistentTangentOperator describes the
646  // change in stress for an arbitrary change in applied
647  // strain. I assume that the change in strain will not
648  // change the type of return (Tip remains Tip, Edge remains
649  // Edge, Plane remains Plane).
650  // I assume isotropic elasticity.
651  //
652  // The consistent tangent operator is a little different
653  // than cases where no rotation to principal stress space
654  // is made during the returnMap. Let S_ij be the stress
655  // in original coordinates, and s_ij be the stress in the
656  // principal stress coordinates, so that
657  // s_ij = diag(eigvals[0], eigvals[1], eigvals[2])
658  // We want dS_ij under an arbitrary change in strain (ep->ep+dep)
659  // dS = S(ep+dep) - S(ep)
660  // = R(ep+dep) s(ep+dep) R(ep+dep)^T - R(ep) s(ep) R(ep)^T
661  // Here R = the rotation to principal-stress space, ie
662  // R_ij = eigvecs[i][j] = i^th component of j^th eigenvector
663  // Expanding to first order in dep,
664  // dS = R(ep) (s(ep+dep) - s(ep)) R(ep)^T
665  // + dR/dep s(ep) R^T + R(ep) s(ep) dR^T/dep
666  // The first line is all that is usually calculated in the
667  // consistent tangent operator calculation, and is called
668  // cto below.
669  // The second line involves changes in the eigenvectors, and
670  // is called sec below.
671 
672  RankFourTensor cto;
673  const Real hard = dtensile_strength(intnl);
674  const Real la = E_ijkl(0, 0, 1, 1);
675  const Real mu = 0.5 * (E_ijkl(0, 0, 0, 0) - la);
676 
677  if (cumulative_pm[1] <= 0)
678  {
679  // only cumulative_pm[2] is positive, so this is return to the Plane
680  const Real denom = hard + la + 2 * mu;
681  const Real al = la * la / denom;
682  const Real be = la * (la + 2 * mu) / denom;
683  const Real ga = hard * (la + 2 * mu) / denom;
684  std::vector<Real> comps(9);
685  comps[0] = comps[4] = la + 2 * mu - al;
686  comps[1] = comps[3] = la - al;
687  comps[2] = comps[5] = comps[6] = comps[7] = la - be;
688  comps[8] = ga;
689  cto.fillFromInputVector(comps, RankFourTensor::principal);
690  }
691  else if (cumulative_pm[0] <= 0)
692  {
693  // both cumulative_pm[2] and cumulative_pm[1] are positive, so Edge
694  const Real denom = 2 * hard + 2 * la + 2 * mu;
695  const Real al = hard * 2 * la / denom;
696  const Real be = hard * (2 * la + 2 * mu) / denom;
697  std::vector<Real> comps(9);
698  comps[0] = la + 2 * mu - 2 * la * la / denom;
699  comps[1] = comps[2] = al;
700  comps[3] = comps[6] = al;
701  comps[4] = comps[5] = comps[7] = comps[8] = be;
702  cto.fillFromInputVector(comps, RankFourTensor::principal);
703  }
704  else
705  {
706  // all cumulative_pm are positive, so Tip
707  const Real denom = 3 * hard + 3 * la + 2 * mu;
708  std::vector<Real> comps(2);
709  comps[0] = hard * (3 * la + 2 * mu) / denom;
710  comps[1] = 0;
711  cto.fillFromInputVector(comps, RankFourTensor::symmetric_isotropic);
712  }
713 
714  cto.rotate(trial_eigvecs);
715 
716  // drdsig = change in eigenvectors under a small stress change
717  // drdsig(i,j,m,n) = dR(i,j)/dS_mn
718  // The formula below is fairly easily derived:
719  // S R = R s, so taking the variation
720  // dS R + S dR = dR s + R ds, and multiplying by R^T
721  // R^T dS R + R^T S dR = R^T dR s + ds .... (eqn 1)
722  // I demand that RR^T = 1 = R^T R, and also that
723  // (R+dR)(R+dR)^T = 1 = (R+dT)^T (R+dR), which means
724  // that dR = R*c, for some antisymmetric c, so Eqn1 reads
725  // R^T dS R + s c = c s + ds
726  // Grabbing the components of this gives ds/dS (already
727  // in RankTwoTensor), and c, which is:
728  // dR_ik/dS_mn = drdsig(i, k, m, n) = trial_eigvecs(m, b)*trial_eigvecs(n, k)*trial_eigvecs(i,
729  // b)/(trial_eigvals[k] - trial_eigvals[b]);
730  // (sum over b!=k).
731 
732  RankFourTensor drdsig;
733  for (unsigned k = 0; k < 3; ++k)
734  for (unsigned b = 0; b < 3; ++b)
735  {
736  if (b == k)
737  continue;
738  for (unsigned m = 0; m < 3; ++m)
739  for (unsigned n = 0; n < 3; ++n)
740  for (unsigned i = 0; i < 3; ++i)
741  drdsig(i, k, m, n) += trial_eigvecs(m, b) * trial_eigvecs(n, k) * trial_eigvecs(i, b) /
742  (trial_eigvals[k] - trial_eigvals[b]);
743  }
744 
745  // With diagla = diag(eigvals[0], eigvals[1], digvals[2])
746  // The following implements
747  // ans(i, j, a, b) += (drdsig(i, k, m, n)*trial_eigvecs(j, l)*diagla(k, l) + trial_eigvecs(i,
748  // k)*drdsig(j, l, m, n)*diagla(k, l))*E_ijkl(m, n, a, b);
749  // (sum over k, l, m and n)
750 
751  RankFourTensor ans;
752  for (unsigned i = 0; i < 3; ++i)
753  for (unsigned j = 0; j < 3; ++j)
754  for (unsigned a = 0; a < 3; ++a)
755  for (unsigned k = 0; k < 3; ++k)
756  for (unsigned m = 0; m < 3; ++m)
757  ans(i, j, a, a) += (drdsig(i, k, m, m) * trial_eigvecs(j, k) +
758  trial_eigvecs(i, k) * drdsig(j, k, m, m)) *
759  eigvals[k] * la; // E_ijkl(m, n, a, b) = la*(m==n)*(a==b);
760 
761  for (unsigned i = 0; i < 3; ++i)
762  for (unsigned j = 0; j < 3; ++j)
763  for (unsigned a = 0; a < 3; ++a)
764  for (unsigned b = 0; b < 3; ++b)
765  for (unsigned k = 0; k < 3; ++k)
766  {
767  ans(i, j, a, b) += (drdsig(i, k, a, b) * trial_eigvecs(j, k) +
768  trial_eigvecs(i, k) * drdsig(j, k, a, b)) *
769  eigvals[k] * mu; // E_ijkl(m, n, a, b) = mu*(m==a)*(n==b)
770  ans(i, j, a, b) += (drdsig(i, k, b, a) * trial_eigvecs(j, k) +
771  trial_eigvecs(i, k) * drdsig(j, k, b, a)) *
772  eigvals[k] * mu; // E_ijkl(m, n, a, b) = mu*(m==b)*(n==a)
773  }
774 
775  return cto + ans;
776 }
virtual Real dtensile_strength(const Real internal_param) const
d(tensile strength)/d(internal_param) as a function of residual value, rate, and internal_param ...
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.
const bool _use_custom_cto
Whether to use the custom consistent tangent operator calculation.

◆ 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 TensorMechanicsPlasticTensileMulti::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 140 of file TensorMechanicsPlasticTensileMulti.C.

142 {
143  dr_dintnl.assign(3, RankTwoTensor());
144 }

◆ dflowPotential_dstress()

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

◆ dflowPotential_dstressV()

void TensorMechanicsPlasticTensileMulti::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 133 of file TensorMechanicsPlasticTensileMulti.C.

135 {
136  stress.d2symmetricEigenvalues(dr_dstress);
137 }

◆ 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 TensorMechanicsPlasticTensileMulti::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
privatevirtual

Just like returnMap, but a protected interface that definitely uses the algorithm, since returnMap itself does not use the algorithm if _use_returnMap=false.

Definition at line 253 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by activeConstraints(), and returnMap().

263 {
264  mooseAssert(dpm.size() == 3,
265  "TensorMechanicsPlasticTensileMulti size of dpm should be 3 but it is "
266  << dpm.size());
267 
268  std::vector<Real> eigvals;
269  RankTwoTensor eigvecs;
270  trial_stress.symmetricEigenvaluesEigenvectors(eigvals, eigvecs);
271  eigvals[0] += _shift;
272  eigvals[2] -= _shift;
273 
274  Real str = tensile_strength(intnl_old);
275 
276  yf.resize(3);
277  yf[0] = eigvals[0] - str;
278  yf[1] = eigvals[1] - str;
279  yf[2] = eigvals[2] - str;
280 
281  if (yf[0] <= _f_tol && yf[1] <= _f_tol && yf[2] <= _f_tol)
282  {
283  // purely elastic (trial_stress, intnl_old)
284  trial_stress_inadmissible = false;
285  return true;
286  }
287 
288  trial_stress_inadmissible = true;
289  delta_dp.zero();
290  returned_stress.zero();
291 
292  // In the following i often assume that E_ijkl is
293  // for an isotropic situation. This reduces FLOPS
294  // substantially which is important since the returnMap
295  // is potentially the most compute-intensive function
296  // of a simulation.
297  // In many comments i write the general expression, and
298  // i hope that might guide future coders if they are
299  // generalising to a non-istropic E_ijkl
300 
301  // n[alpha] = E_ijkl*r[alpha]_kl expressed in principal stress space
302  // (alpha = 0, 1, 2, corresponding to the three surfaces)
303  // Note that in principal stress space, the flow
304  // directions are, expressed in 'vector' form,
305  // r[0] = (1,0,0), r[1] = (0,1,0), r[2] = (0,0,1).
306  // Similar for _n:
307  // so _n[0] = E_ij00*r[0], _n[1] = E_ij11*r[1], _n[2] = E_ij22*r[2]
308  // In the following I assume that the E_ijkl is
309  // for an isotropic situation.
310  // In the anisotropic situation, we couldn't express
311  // the flow directions as vectors in the same principal
312  // stress space as the stress: they'd be full rank-2 tensors
313  std::vector<RealVectorValue> n(3);
314  n[0](0) = E_ijkl(0, 0, 0, 0);
315  n[0](1) = E_ijkl(1, 1, 0, 0);
316  n[0](2) = E_ijkl(2, 2, 0, 0);
317  n[1](0) = E_ijkl(0, 0, 1, 1);
318  n[1](1) = E_ijkl(1, 1, 1, 1);
319  n[1](2) = E_ijkl(2, 2, 1, 1);
320  n[2](0) = E_ijkl(0, 0, 2, 2);
321  n[2](1) = E_ijkl(1, 1, 2, 2);
322  n[2](2) = E_ijkl(2, 2, 2, 2);
323 
324  // With non-zero Poisson's ratio and hardening
325  // it is not computationally cheap to know whether
326  // the trial stress will return to the tip, edge,
327  // or plane. The following is correct for zero
328  // Poisson's ratio and no hardening, and at least
329  // gives a not-completely-stupid guess in the
330  // more general case.
331  // trial_order[0] = type of return to try first
332  // trial_order[1] = type of return to try second
333  // trial_order[2] = type of return to try third
334  const unsigned int number_of_return_paths = 3;
335  std::vector<int> trial_order(number_of_return_paths);
336  if (yf[0] > _f_tol) // all the yield functions are positive, since eigvals are ordered eigvals[0]
337  // <= eigvals[1] <= eigvals[2]
338  {
339  trial_order[0] = tip;
340  trial_order[1] = edge;
341  trial_order[2] = plane;
342  }
343  else if (yf[1] > _f_tol) // two yield functions are positive
344  {
345  trial_order[0] = edge;
346  trial_order[1] = tip;
347  trial_order[2] = plane;
348  }
349  else
350  {
351  trial_order[0] = plane;
352  trial_order[1] = edge;
353  trial_order[2] = tip;
354  }
355 
356  unsigned trial;
357  bool nr_converged = false;
358  for (trial = 0; trial < number_of_return_paths; ++trial)
359  {
360  switch (trial_order[trial])
361  {
362  case tip:
363  nr_converged = returnTip(eigvals, n, dpm, returned_stress, intnl_old, 0);
364  break;
365  case edge:
366  nr_converged = returnEdge(eigvals, n, dpm, returned_stress, intnl_old, 0);
367  break;
368  case plane:
369  nr_converged = returnPlane(eigvals, n, dpm, returned_stress, intnl_old, 0);
370  break;
371  }
372 
373  str = tensile_strength(intnl_old + dpm[0] + dpm[1] + dpm[2]);
374 
375  if (nr_converged && KuhnTuckerOK(returned_stress, dpm, str, ep_plastic_tolerance))
376  break;
377  }
378 
379  if (trial == number_of_return_paths)
380  {
381  Moose::err << "Trial stress = \n";
382  trial_stress.print(Moose::err);
383  Moose::err << "Internal parameter = " << intnl_old << "\n";
384  mooseError("TensorMechanicsPlasticTensileMulti: FAILURE! You probably need to implement a "
385  "line search\n");
386  // failure - must place yield function values at trial stress into yf
387  str = tensile_strength(intnl_old);
388  yf[0] = eigvals[0] - str;
389  yf[1] = eigvals[1] - str;
390  yf[2] = eigvals[2] - str;
391  return false;
392  }
393 
394  // success
395 
396  returned_intnl = intnl_old;
397  for (unsigned i = 0; i < 3; ++i)
398  {
399  yf[i] = returned_stress(i, i) - str;
400  delta_dp(i, i) = dpm[i];
401  returned_intnl += dpm[i];
402  }
403  returned_stress = eigvecs * returned_stress * (eigvecs.transpose());
404  delta_dp = eigvecs * delta_dp * (eigvecs.transpose());
405  return true;
406 }
bool returnTip(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real initial_guess) const
Tries to return-map to the Tensile tip.
bool KuhnTuckerOK(const RankTwoTensor &returned_diagonal_stress, const std::vector< Real > &dpm, Real str, Real ep_plastic_tolerance) const
Returns true if the Kuhn-Tucker conditions are satisfied.
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
bool returnEdge(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real initial_guess) const
Tries to return-map to the Tensile edge.
virtual Real tensile_strength(const Real internal_param) const
tensile strength as a function of residual value, rate, and internal_param
const Real _f_tol
Tolerance on yield function.
bool returnPlane(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real initial_guess) const
Tries to return-map to the Tensile plane The return value is true if the internal Newton-Raphson proc...

◆ dot()

Real TensorMechanicsPlasticTensileMulti::dot ( const std::vector< Real > &  a,
const std::vector< Real > &  b 
) const
private

dot product of two 3-dimensional vectors

Definition at line 195 of file TensorMechanicsPlasticTensileMulti.C.

197 {
198  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
199 }

◆ dtensile_strength()

Real TensorMechanicsPlasticTensileMulti::dtensile_strength ( const Real  internal_param) const
protectedvirtual

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

Definition at line 153 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by consistentTangentOperator(), dyieldFunction_dintnlV(), returnEdge(), returnPlane(), and returnTip().

154 {
155  return _strength.derivative(internal_param);
156 }
virtual Real derivative(Real intnl) const
const TensorMechanicsHardeningModel & _strength

◆ 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 TensorMechanicsPlasticTensileMulti::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 116 of file TensorMechanicsPlasticTensileMulti.C.

119 {
120  df_dintnl.assign(3, -dtensile_strength(intnl));
121 }
virtual Real dtensile_strength(const Real internal_param) const
d(tensile strength)/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 TensorMechanicsPlasticTensileMulti::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 91 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by flowPotentialV().

93 {
94  std::vector<Real> eigvals;
95  stress.dsymmetricEigenvalues(eigvals, df_dstress);
96 
97  if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
98  {
99  Real small_perturbation;
100  RankTwoTensor shifted_stress = stress;
101  while (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
102  {
103  for (unsigned i = 0; i < 3; ++i)
104  for (unsigned j = 0; j <= i; ++j)
105  {
106  small_perturbation = 0.1 * _shift * 2 * (MooseRandom::rand() - 0.5);
107  shifted_stress(i, j) += small_perturbation;
108  shifted_stress(j, i) += small_perturbation;
109  }
110  shifted_stress.dsymmetricEigenvalues(eigvals, df_dstress);
111  }
112  }
113 }
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...

◆ 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 TensorMechanicsPlasticTensileMulti::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 124 of file TensorMechanicsPlasticTensileMulti.C.

127 {
128  // This plasticity is associative so
129  dyieldFunction_dstressV(stress, intnl, r);
130 }
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.

◆ 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 TensorMechanicsPlasticTensileMulti::KuhnTuckerOK ( const RankTwoTensor &  returned_diagonal_stress,
const std::vector< Real > &  dpm,
Real  str,
Real  ep_plastic_tolerance 
) const
private

Returns true if the Kuhn-Tucker conditions are satisfied.

Parameters
returned_diagonal_stressThe eigenvalues (sorted in ascending order as is standard in this Class) are stored in the diagonal components
dpmThe three plastic multipliers
strThe yield strength
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 601 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by doReturnMap().

605 {
606  for (unsigned i = 0; i < 3; ++i)
608  returned_diagonal_stress(i, i) - str, dpm[i], ep_plastic_tolerance))
609  return false;
610  return true;
611 }
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 TensorMechanicsPlasticMohrCoulombMulti::KuhnTuckerOK(), 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 TensorMechanicsPlasticTensileMulti::modelName ( ) const
overridevirtual

Implements TensorMechanicsPlasticModel.

Definition at line 211 of file TensorMechanicsPlasticTensileMulti.C.

212 {
213  return "TensileMulti";
214 }

◆ numberSurfaces()

unsigned int TensorMechanicsPlasticTensileMulti::numberSurfaces ( ) const
overridevirtual

The number of yield surfaces for this plasticity model.

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 70 of file TensorMechanicsPlasticTensileMulti.C.

71 {
72  return 3;
73 }

◆ returnEdge()

bool TensorMechanicsPlasticTensileMulti::returnEdge ( const std::vector< Real > &  eigvals,
const std::vector< RealVectorValue > &  n,
std::vector< Real > &  dpm,
RankTwoTensor &  returned_stress,
Real  intnl_old,
Real  initial_guess 
) const
private

Tries to return-map to the Tensile edge.

The return value is true if the internal Newton-Raphson process has converged, otherwise it is false

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
dpmThe three plastic multipliers resulting from the return-map to the edge. This algorithm doesn't do Kuhn-Tucker checking, so these could be positive or negative or zero (dpm[0]=0 always for Edge return).
returned_stressThe 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.
initial_guessA guess of dpm[1]+dpm[2]

Definition at line 486 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by doReturnMap().

492 {
493  // work out the point to which we would return, "a". It is defined by
494  // f1 = 0 = f2, and the normality condition:
495  // (eigvals - a).(n1 x n2) = 0,
496  // where eigvals is the starting position
497  // (it is a vector in principal stress space).
498  // To get f1=0=f2, we need a = (a0, str, str), and a0 is found
499  // by expanding the normality condition to yield:
500  // a0 = (-str*n1xn2[1] - str*n1xn2[2] + edotn1xn2)/n1xn2[0];
501  // where edotn1xn2 = eigvals.(n1 x n2)
502  //
503  // We need to find the plastic multipliers, dpm, defined by
504  // eigvals - a = dpm[1]*n1 + dpm[2]*n2
505  // For the isotropic case, and defining eminusa = eigvals - a,
506  // the solution is easy:
507  // dpm[0] = 0;
508  // dpm[1] = (eminusa[1] - eminusa[0])/(n[1][1] - n[1][0]);
509  // dpm[2] = (eminusa[2] - eminusa[0])/(n[2][2] - n[2][0]);
510  //
511  // Now specialise to the isotropic case. Define
512  // x = dpm[1] + dpm[2] = (eigvals[1] + eigvals[2] - 2*str)/(n[0][0] + n[0][1])
513  // Notice that the RHS is a function of x, so we solve using
514  // Newton-Raphson starting with x=initial_guess
515  Real x = initial_guess;
516  const Real denom = n[0](0) + n[0](1);
517  Real str = tensile_strength(intnl_old + x);
518 
519  if (_strength.modelName().compare("Constant") != 0)
520  {
521  // Finding x is expensive. Therefore
522  // although x!=0 for Constant Hardening, solution
523  // for dpm above demonstrates that we don't
524  // actually need to know x to find the dpm for
525  // Constant Hardening.
526  //
527  // However, for nontrivial Hardening, the following
528  // is necessary
529  const Real eig = eigvals[1] + eigvals[2];
530  Real residual = denom * x - eig + 2 * str;
531  Real jacobian;
532  unsigned int iter = 0;
533  do
534  {
535  jacobian = denom + 2 * dtensile_strength(intnl_old + x);
536  x += -residual / jacobian;
537  if (iter > _max_iters) // not converging
538  return false;
539  str = tensile_strength(intnl_old + x);
540  residual = denom * x - eig + 2 * str;
541  iter++;
542  } while (residual * residual > _f_tol * _f_tol);
543  }
544 
545  dpm[0] = 0;
546  dpm[1] = ((eigvals[1] * n[0](0) - eigvals[2] * n[0](1)) / (n[0](0) - n[0](1)) - str) / denom;
547  dpm[2] = ((eigvals[2] * n[0](0) - eigvals[1] * n[0](1)) / (n[0](0) - n[0](1)) - str) / denom;
548 
549  returned_stress(0, 0) = eigvals[0] - n[0](1) * (dpm[1] + dpm[2]);
550  returned_stress(1, 1) = returned_stress(2, 2) = str;
551  return true;
552 }
virtual Real dtensile_strength(const Real internal_param) const
d(tensile strength)/d(internal_param) as a function of residual value, rate, and internal_param ...
virtual std::string modelName() const =0
const unsigned int _max_iters
maximum iterations allowed in the custom return-map algorithm
virtual Real tensile_strength(const Real internal_param) const
tensile strength as a function of residual value, rate, and internal_param
const Real _f_tol
Tolerance on yield function.
const TensorMechanicsHardeningModel & _strength

◆ returnMap()

bool TensorMechanicsPlasticTensileMulti::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 217 of file TensorMechanicsPlasticTensileMulti.C.

227 {
229  return TensorMechanicsPlasticModel::returnMap(trial_stress,
230  intnl_old,
231  E_ijkl,
232  ep_plastic_tolerance,
233  returned_stress,
234  returned_intnl,
235  dpm,
236  delta_dp,
237  yf,
238  trial_stress_inadmissible);
239 
240  return doReturnMap(trial_stress,
241  intnl_old,
242  E_ijkl,
243  ep_plastic_tolerance,
244  returned_stress,
245  returned_intnl,
246  dpm,
247  delta_dp,
248  yf,
249  trial_stress_inadmissible);
250 }
const bool _use_custom_returnMap
Whether to use the custom return-map algorithm.
virtual 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
Just like returnMap, but a protected interface that definitely uses the algorithm, since returnMap itself does not use the algorithm if _use_returnMap=false.
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.

◆ returnPlane()

bool TensorMechanicsPlasticTensileMulti::returnPlane ( const std::vector< Real > &  eigvals,
const std::vector< RealVectorValue > &  n,
std::vector< Real > &  dpm,
RankTwoTensor &  returned_stress,
Real  intnl_old,
Real  initial_guess 
) const
private

Tries to return-map to the Tensile plane The return value is true if the internal Newton-Raphson process has converged, otherwise it is false.

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
dpmThe three plastic multipliers resulting from the return-map to the plane. This algorithm doesn't do Kuhn-Tucker checking, so dpm[2] could be positive or negative or zero (dpm[0]=dpm[1]=0 always for Plane return).
returned_stressThe 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.
initial_guessA guess of dpm[2]

Definition at line 555 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by doReturnMap().

561 {
562  // the returned point, "a", is defined by f2=0 and
563  // a = p - dpm[2]*n2.
564  // This is a vector equation in
565  // principal stress space, and dpm[2] is the third
566  // plasticity multiplier (dpm[0]=0=dpm[1] for return
567  // to the plane) and "p" is the starting
568  // position (p=eigvals).
569  // (Kuhn-Tucker demands that dpm[2]>=0, but we leave checking
570  // that condition for later.)
571  // Since f2=0, we must have a[2]=tensile_strength,
572  // so we can just look at the [2] component of the
573  // equation, which yields
574  // n[2][2]*dpm[2] - eigvals[2] + str = 0
575  // For hardening, str=tensile_strength(intnl_old+dpm[2]),
576  // and we want to solve for dpm[2].
577  // Use Newton-Raphson with initial guess dpm[2] = initial_guess
578  dpm[2] = initial_guess;
579  Real residual = n[2](2) * dpm[2] - eigvals[2] + tensile_strength(intnl_old + dpm[2]);
580  Real jacobian;
581  unsigned int iter = 0;
582  do
583  {
584  jacobian = n[2](2) + dtensile_strength(intnl_old + dpm[2]);
585  dpm[2] += -residual / jacobian;
586  if (iter > _max_iters) // not converging
587  return false;
588  residual = n[2](2) * dpm[2] - eigvals[2] + tensile_strength(intnl_old + dpm[2]);
589  iter++;
590  } while (residual * residual > _f_tol * _f_tol);
591 
592  dpm[0] = 0;
593  dpm[1] = 0;
594  returned_stress(0, 0) = eigvals[0] - dpm[2] * n[2](0);
595  returned_stress(1, 1) = eigvals[1] - dpm[2] * n[2](1);
596  returned_stress(2, 2) = eigvals[2] - dpm[2] * n[2](2);
597  return true;
598 }
virtual Real dtensile_strength(const Real internal_param) const
d(tensile strength)/d(internal_param) as a function of residual value, rate, and internal_param ...
const unsigned int _max_iters
maximum iterations allowed in the custom return-map algorithm
virtual Real tensile_strength(const Real internal_param) const
tensile strength as a function of residual value, rate, and internal_param
const Real _f_tol
Tolerance on yield function.

◆ returnTip()

bool TensorMechanicsPlasticTensileMulti::returnTip ( const std::vector< Real > &  eigvals,
const std::vector< RealVectorValue > &  n,
std::vector< Real > &  dpm,
RankTwoTensor &  returned_stress,
Real  intnl_old,
Real  initial_guess 
) const
private

Tries to return-map to the Tensile tip.

The return value is true if the internal Newton-Raphson process has converged, otherwise it is false

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
dpmThe three plastic multipliers resulting from the return-map to the tip. This algorithm doesn't do Kuhn-Tucker checking, so these could be positive or negative or zero
returned_stressThe 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.
initial_guessA guess of dpm[0]+dpm[1]+dpm[2]

Definition at line 409 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by doReturnMap().

415 {
416  // The returned point is defined by f0=f1=f2=0.
417  // that is, returned_stress = diag(str, str, str), where
418  // str = tensile_strength(intnl),
419  // where intnl = intnl_old + dpm[0] + dpm[1] + dpm[2]
420  // The 3 plastic multipliers, dpm, are defiend by the normality condition
421  // eigvals - str = dpm[0]*n[0] + dpm[1]*n[1] + dpm[2]*n[2]
422  // (Kuhn-Tucker demands that all dpm are non-negative, but we leave
423  // that checking for later.)
424  // This is a vector equation with solution (A):
425  // dpm[0] = triple(eigvals - str, n[1], n[2])/trip;
426  // dpm[1] = triple(eigvals - str, n[2], n[0])/trip;
427  // dpm[2] = triple(eigvals - str, n[0], n[1])/trip;
428  // where trip = triple(n[0], n[1], n[2]).
429  // By adding the three components of that solution together
430  // we can get an equation for x = dpm[0] + dpm[1] + dpm[2],
431  // and then our Newton-Raphson only involves one variable (x).
432  // In the following, i specialise to the isotropic situation.
433 
434  Real x = initial_guess;
435  const Real denom = (n[0](0) - n[0](1)) * (n[0](0) + 2 * n[0](1));
436  Real str = tensile_strength(intnl_old + x);
437 
438  if (_strength.modelName().compare("Constant") != 0)
439  {
440  // Finding x is expensive. Therefore
441  // although x!=0 for Constant Hardening, solution (A)
442  // demonstrates that we don't
443  // actually need to know x to find the dpm for
444  // Constant Hardening.
445  //
446  // However, for nontrivial Hardening, the following
447  // is necessary
448  const Real eig = eigvals[0] + eigvals[1] + eigvals[2];
449  const Real bul = (n[0](0) + 2 * n[0](1));
450 
451  // and finally, the equation we want to solve is:
452  // bul*x - eig + 3*str = 0
453  // where str=tensile_strength(intnl_old + x)
454  // and x = dpm[0] + dpm[1] + dpm[2]
455  // (Note this has units of stress, so using _f_tol as a convergence check is reasonable.)
456  // Use Netwon-Raphson with initial guess x = 0
457  Real residual = bul * x - eig + 3 * str;
458  Real jacobian;
459  unsigned int iter = 0;
460  do
461  {
462  jacobian = bul + 3 * dtensile_strength(intnl_old + x);
463  x += -residual / jacobian;
464  if (iter > _max_iters) // not converging
465  return false;
466  str = tensile_strength(intnl_old + x);
467  residual = bul * x - eig + 3 * str;
468  iter++;
469  } while (residual * residual > _f_tol * _f_tol);
470  }
471 
472  // The following is the solution (A) written above
473  // (dpm[0] = triple(eigvals - str, n[1], n[2])/trip, etc)
474  // in the isotropic situation
475  dpm[0] = (n[0](0) * (eigvals[0] - str) + n[0](1) * (eigvals[0] - eigvals[1] - eigvals[2] + str)) /
476  denom;
477  dpm[1] = (n[0](0) * (eigvals[1] - str) + n[0](1) * (eigvals[1] - eigvals[2] - eigvals[0] + str)) /
478  denom;
479  dpm[2] = (n[0](0) * (eigvals[2] - str) + n[0](1) * (eigvals[2] - eigvals[0] - eigvals[1] + str)) /
480  denom;
481  returned_stress(0, 0) = returned_stress(1, 1) = returned_stress(2, 2) = str;
482  return true;
483 }
virtual Real dtensile_strength(const Real internal_param) const
d(tensile strength)/d(internal_param) as a function of residual value, rate, and internal_param ...
virtual std::string modelName() const =0
const unsigned int _max_iters
maximum iterations allowed in the custom return-map algorithm
virtual Real tensile_strength(const Real internal_param) const
tensile strength as a function of residual value, rate, and internal_param
const Real _f_tol
Tolerance on yield function.
const TensorMechanicsHardeningModel & _strength

◆ tensile_strength()

Real TensorMechanicsPlasticTensileMulti::tensile_strength ( const Real  internal_param) const
protectedvirtual

tensile strength as a function of residual value, rate, and internal_param

Definition at line 147 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by doReturnMap(), returnEdge(), returnPlane(), returnTip(), and yieldFunctionV().

148 {
149  return _strength.value(internal_param);
150 }
virtual Real value(Real intnl) const
const TensorMechanicsHardeningModel & _strength

◆ triple()

Real TensorMechanicsPlasticTensileMulti::triple ( const std::vector< Real > &  a,
const std::vector< Real > &  b,
const std::vector< Real > &  c 
) const
private

triple product of three 3-dimensional vectors

Definition at line 202 of file TensorMechanicsPlasticTensileMulti.C.

205 {
206  return a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) +
207  a[2] * (b[0] * c[1] - b[1] * c[0]);
208 }

◆ useCustomCTO()

bool TensorMechanicsPlasticTensileMulti::useCustomCTO ( ) const
overridevirtual

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

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 785 of file TensorMechanicsPlasticTensileMulti.C.

786 {
787  return _use_custom_cto;
788 }
const bool _use_custom_cto
Whether to use the custom consistent tangent operator calculation.

◆ useCustomReturnMap()

bool TensorMechanicsPlasticTensileMulti::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 779 of file TensorMechanicsPlasticTensileMulti.C.

780 {
781  return _use_custom_returnMap;
782 }
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 }

◆ yieldFunctionV()

void TensorMechanicsPlasticTensileMulti::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 76 of file TensorMechanicsPlasticTensileMulti.C.

79 {
80  std::vector<Real> eigvals;
81  stress.symmetricEigenvalues(eigvals);
82  const Real str = tensile_strength(intnl);
83 
84  f.resize(3);
85  f[0] = eigvals[0] + _shift - str;
86  f[1] = eigvals[1] - str;
87  f[2] = eigvals[2] - _shift - str;
88 }
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
virtual Real tensile_strength(const Real internal_param) const
tensile strength as a function of residual value, rate, and internal_param

Member Data Documentation

◆ _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 TensorMechanicsPlasticTensileMulti::_max_iters
private

maximum iterations allowed in the custom return-map algorithm

Definition at line 98 of file TensorMechanicsPlasticTensileMulti.h.

Referenced by returnEdge(), returnPlane(), and returnTip().

◆ _shift

const Real TensorMechanicsPlasticTensileMulti::_shift
private

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

Definition at line 101 of file TensorMechanicsPlasticTensileMulti.h.

Referenced by doReturnMap(), dyieldFunction_dstressV(), TensorMechanicsPlasticTensileMulti(), and yieldFunctionV().

◆ _strength

const TensorMechanicsHardeningModel& TensorMechanicsPlasticTensileMulti::_strength
private

◆ _use_custom_cto

const bool TensorMechanicsPlasticTensileMulti::_use_custom_cto
private

Whether to use the custom consistent tangent operator calculation.

Definition at line 107 of file TensorMechanicsPlasticTensileMulti.h.

Referenced by consistentTangentOperator(), and useCustomCTO().

◆ _use_custom_returnMap

const bool TensorMechanicsPlasticTensileMulti::_use_custom_returnMap
private

Whether to use the custom return-map algorithm.

Definition at line 104 of file TensorMechanicsPlasticTensileMulti.h.

Referenced by returnMap(), and useCustomReturnMap().


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