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

Member Enumeration Documentation

Constructor & Destructor Documentation

TensorMechanicsPlasticTensileMulti::TensorMechanicsPlasticTensileMulti ( const InputParameters &  parameters)

Definition at line 44 of file TensorMechanicsPlasticTensileMulti.C.

46  : TensorMechanicsPlasticModel(parameters),
47  _strength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
48  _max_iters(getParam<unsigned int>("max_iterations")),
49  _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
50  _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
51  _use_custom_cto(getParam<bool>("use_custom_cto"))
52 {
53  if (_shift < 0)
54  mooseError("Value of 'shift' in TensorMechanicsPlasticTensileMulti must not be negative\n");
55  if (_shift > _f_tol)
56  _console << "WARNING: value of 'shift' in TensorMechanicsPlasticTensileMulti is probably set "
57  "too high\n";
58  if (LIBMESH_DIM != 3)
59  mooseError("TensorMechanicsPlasticTensileMulti is only defined for LIBMESH_DIM=3");
60  MooseRandom::seed(0);
61 }
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

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 153 of file TensorMechanicsPlasticTensileMulti.C.

159 {
160  act.assign(3, false);
161 
162  if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol)
163  {
164  returned_stress = stress;
165  return;
166  }
167 
168  Real returned_intnl;
169  std::vector<Real> dpm(3);
170  RankTwoTensor delta_dp;
171  std::vector<Real> yf(3);
172  bool trial_stress_inadmissible;
173  doReturnMap(stress,
174  intnl,
175  Eijkl,
176  0.0,
177  returned_stress,
178  returned_intnl,
179  dpm,
180  delta_dp,
181  yf,
182  trial_stress_inadmissible);
183 
184  for (unsigned i = 0; i < 3; ++i)
185  act[i] = (dpm[i] > 0);
186 }
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.
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 608 of file TensorMechanicsPlasticTensileMulti.C.

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

136 {
137  dr_dintnl.assign(3, RankTwoTensor());
138 }
RankFourTensor TensorMechanicsPlasticModel::dflowPotential_dstress ( const RankTwoTensor &  stress,
Real  intnl 
) const
protectedvirtualinherited
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 127 of file TensorMechanicsPlasticTensileMulti.C.

129 {
130  stress.d2symmetricEigenvalues(dr_dstress);
131 }
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 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 247 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by activeConstraints(), and returnMap().

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

191 {
192  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
193 }
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 147 of file TensorMechanicsPlasticTensileMulti.C.

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

148 {
149  return _strength.derivative(internal_param);
150 }
virtual Real derivative(Real intnl) const
const TensorMechanicsHardeningModel & _strength
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 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 110 of file TensorMechanicsPlasticTensileMulti.C.

113 {
114  df_dintnl.assign(3, -dtensile_strength(intnl));
115 }
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 ...
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 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 85 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by flowPotentialV().

87 {
88  std::vector<Real> eigvals;
89  stress.dsymmetricEigenvalues(eigvals, df_dstress);
90 
91  if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
92  {
93  Real small_perturbation;
94  RankTwoTensor shifted_stress = stress;
95  while (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
96  {
97  for (unsigned i = 0; i < 3; ++i)
98  for (unsigned j = 0; j <= i; ++j)
99  {
100  small_perturbation = 0.1 * _shift * 2 * (MooseRandom::rand() - 0.5);
101  shifted_stress(i, j) += small_perturbation;
102  shifted_stress(j, i) += small_perturbation;
103  }
104  shifted_stress.dsymmetricEigenvalues(eigvals, df_dstress);
105  }
106  }
107 }
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
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 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 118 of file TensorMechanicsPlasticTensileMulti.C.

121 {
122  // This plasticity is associative so
123  dyieldFunction_dstressV(stress, intnl, r);
124 }
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.
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 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 595 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by doReturnMap().

599 {
600  for (unsigned i = 0; i < 3; ++i)
602  returned_diagonal_stress(i, i) - str, dpm[i], ep_plastic_tolerance))
603  return false;
604  return true;
605 }
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 TensorMechanicsPlasticMohrCoulombMulti::KuhnTuckerOK(), 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 TensorMechanicsPlasticTensileMulti::modelName ( ) const
overridevirtual

Implements TensorMechanicsPlasticModel.

Definition at line 205 of file TensorMechanicsPlasticTensileMulti.C.

206 {
207  return "TensileMulti";
208 }
unsigned int TensorMechanicsPlasticTensileMulti::numberSurfaces ( ) const
overridevirtual

The number of yield surfaces for this plasticity model.

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 64 of file TensorMechanicsPlasticTensileMulti.C.

65 {
66  return 3;
67 }
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 480 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by doReturnMap().

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

221 {
223  return TensorMechanicsPlasticModel::returnMap(trial_stress,
224  intnl_old,
225  E_ijkl,
226  ep_plastic_tolerance,
227  returned_stress,
228  returned_intnl,
229  dpm,
230  delta_dp,
231  yf,
232  trial_stress_inadmissible);
233 
234  return doReturnMap(trial_stress,
235  intnl_old,
236  E_ijkl,
237  ep_plastic_tolerance,
238  returned_stress,
239  returned_intnl,
240  dpm,
241  delta_dp,
242  yf,
243  trial_stress_inadmissible);
244 }
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.
const bool _use_custom_returnMap
Whether to use the custom return-map algorithm.
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 549 of file TensorMechanicsPlasticTensileMulti.C.

Referenced by doReturnMap().

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

Referenced by doReturnMap().

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

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

142 {
143  return _strength.value(internal_param);
144 }
virtual Real value(Real intnl) const
const TensorMechanicsHardeningModel & _strength
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 196 of file TensorMechanicsPlasticTensileMulti.C.

199 {
200  return a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) +
201  a[2] * (b[0] * c[1] - b[1] * c[0]);
202 }
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 779 of file TensorMechanicsPlasticTensileMulti.C.

780 {
781  return _use_custom_cto;
782 }
const bool _use_custom_cto
Whether to use the custom consistent tangent operator calculation.
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 773 of file TensorMechanicsPlasticTensileMulti.C.

774 {
775  return _use_custom_returnMap;
776 }
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 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 70 of file TensorMechanicsPlasticTensileMulti.C.

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

Member Data Documentation

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

maximum iterations allowed in the custom return-map algorithm

Definition at line 95 of file TensorMechanicsPlasticTensileMulti.h.

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

const Real TensorMechanicsPlasticTensileMulti::_shift
private

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

Definition at line 98 of file TensorMechanicsPlasticTensileMulti.h.

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

const TensorMechanicsHardeningModel& TensorMechanicsPlasticTensileMulti::_strength
private
const bool TensorMechanicsPlasticTensileMulti::_use_custom_cto
private

Whether to use the custom consistent tangent operator calculation.

Definition at line 104 of file TensorMechanicsPlasticTensileMulti.h.

Referenced by consistentTangentOperator(), and useCustomCTO().

const bool TensorMechanicsPlasticTensileMulti::_use_custom_returnMap
private

Whether to use the custom return-map algorithm.

Definition at line 101 of file TensorMechanicsPlasticTensileMulti.h.

Referenced by returnMap(), and useCustomReturnMap().


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