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

Static Public Member Functions
static InputParameters	validParams ()

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

Member Enumeration Documentation

◆ return_type

enum TensorMechanicsPlasticTensileMulti::return_type

private

Enumerator
tip
edge
plane

Definition at line 219 of file TensorMechanicsPlasticTensileMulti.h.

   {
     tip = 0,
     edge = 1,
     plane = 2
   };

Constructor & Destructor Documentation

◆ TensorMechanicsPlasticTensileMulti()

TensorMechanicsPlasticTensileMulti::TensorMechanicsPlasticTensileMulti ( const InputParameters & parameters )

Definition at line 51 of file TensorMechanicsPlasticTensileMulti.C.

   : TensorMechanicsPlasticModel(parameters),
     _strength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
     _max_iters(getParam<unsigned int>("max_iterations")),
     _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
     _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
     _use_custom_cto(getParam<bool>("use_custom_cto"))
 {
   if (_shift < 0)
     mooseError("Value of 'shift' in TensorMechanicsPlasticTensileMulti must not be negative\n");
   if (_shift > _f_tol)
     _console << "WARNING: value of 'shift' in TensorMechanicsPlasticTensileMulti is probably set "
                 "too high\n";
   if (LIBMESH_DIM != 3)
     mooseError("TensorMechanicsPlasticTensileMulti is only defined for LIBMESH_DIM=3");
   MooseRandom::seed(0);
 }

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

	f	values of the yield functions
	stress	stress tensor
	intnl	internal parameter
	Eijkl	elasticity tensor (stress = Eijkl*strain)
[out]	act	act[i] = true if the i_th yield function is active
[out]	returned_stress	Approximate value of the returned stress

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 160 of file TensorMechanicsPlasticTensileMulti.C.

 {
   act.assign(3, false);
  
   if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol)
   {
     returned_stress = stress;
     return;
   }
  
   Real returned_intnl;
   std::vector<Real> dpm(3);
   RankTwoTensor delta_dp;
   std::vector<Real> yf(3);
   bool trial_stress_inadmissible;
   doReturnMap(stress,
               intnl,
               Eijkl,
               0.0,
               returned_stress,
               returned_intnl,
               dpm,
               delta_dp,
               yf,
               trial_stress_inadmissible);
  
   for (unsigned i = 0; i < 3; ++i)
     act[i] = (dpm[i] > 0);
 }

◆ 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_old	trial stress before returning
intnl_old	internal parameter before returning
stress	current returned stress state
intnl	internal parameter
E_ijkl	elasticity tensor
cumulative_pm	the cumulative plastic multipliers

Returns: the consistent tangent operator: E_ijkl if not over-ridden

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 615 of file TensorMechanicsPlasticTensileMulti.C.

 {
   if (!_use_custom_cto)
     return TensorMechanicsPlasticModel::consistentTangentOperator(
         trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
  
   mooseAssert(cumulative_pm.size() == 3,
               "TensorMechanicsPlasticTensileMulti size of cumulative_pm should be 3 but it is "
                   << cumulative_pm.size());
  
   if (cumulative_pm[2] <= 0) // All cumulative_pm are non-positive, so this is admissible
     return E_ijkl;
  
   // Need the eigenvalues at the returned configuration
   std::vector<Real> eigvals;
   stress.symmetricEigenvalues(eigvals);
  
   // need to rotate to and from principal stress space
   // using the eigenvectors of the trial configuration
   // (not the returned configuration).
   std::vector<Real> trial_eigvals;
   RankTwoTensor trial_eigvecs;
   trial_stress.symmetricEigenvaluesEigenvectors(trial_eigvals, trial_eigvecs);
  
   // The returnMap will have returned to the Tip, Edge or
   // Plane.  The consistentTangentOperator describes the
   // change in stress for an arbitrary change in applied
   // strain.  I assume that the change in strain will not
   // change the type of return (Tip remains Tip, Edge remains
   // Edge, Plane remains Plane).
   // I assume isotropic elasticity.
   //
   // The consistent tangent operator is a little different
   // than cases where no rotation to principal stress space
   // is made during the returnMap.  Let S_ij be the stress
   // in original coordinates, and s_ij be the stress in the
   // principal stress coordinates, so that
   // s_ij = diag(eigvals[0], eigvals[1], eigvals[2])
   // We want dS_ij under an arbitrary change in strain (ep->ep+dep)
   // dS = S(ep+dep) - S(ep)
   //    = R(ep+dep) s(ep+dep) R(ep+dep)^T - R(ep) s(ep) R(ep)^T
   // Here R = the rotation to principal-stress space, ie
   // R_ij = eigvecs[i][j] = i^th component of j^th eigenvector
   // Expanding to first order in dep,
   // dS = R(ep) (s(ep+dep) - s(ep)) R(ep)^T
   //      + dR/dep s(ep) R^T + R(ep) s(ep) dR^T/dep
   // The first line is all that is usually calculated in the
   // consistent tangent operator calculation, and is called
   // cto below.
   // The second line involves changes in the eigenvectors, and
   // is called sec below.
  
   RankFourTensor cto;
   const Real hard = dtensile_strength(intnl);
   const Real la = E_ijkl(0, 0, 1, 1);
   const Real mu = 0.5 * (E_ijkl(0, 0, 0, 0) - la);
  
   if (cumulative_pm[1] <= 0)
   {
     // only cumulative_pm[2] is positive, so this is return to the Plane
     const Real denom = hard + la + 2 * mu;
     const Real al = la * la / denom;
     const Real be = la * (la + 2 * mu) / denom;
     const Real ga = hard * (la + 2 * mu) / denom;
     std::vector<Real> comps(9);
     comps[0] = comps[4] = la + 2 * mu - al;
     comps[1] = comps[3] = la - al;
     comps[2] = comps[5] = comps[6] = comps[7] = la - be;
     comps[8] = ga;
     cto.fillFromInputVector(comps, RankFourTensor::principal);
   }
   else if (cumulative_pm[0] <= 0)
   {
     // both cumulative_pm[2] and cumulative_pm[1] are positive, so Edge
     const Real denom = 2 * hard + 2 * la + 2 * mu;
     const Real al = hard * 2 * la / denom;
     const Real be = hard * (2 * la + 2 * mu) / denom;
     std::vector<Real> comps(9);
     comps[0] = la + 2 * mu - 2 * la * la / denom;
     comps[1] = comps[2] = al;
     comps[3] = comps[6] = al;
     comps[4] = comps[5] = comps[7] = comps[8] = be;
     cto.fillFromInputVector(comps, RankFourTensor::principal);
   }
   else
   {
     // all cumulative_pm are positive, so Tip
     const Real denom = 3 * hard + 3 * la + 2 * mu;
     std::vector<Real> comps(2);
     comps[0] = hard * (3 * la + 2 * mu) / denom;
     comps[1] = 0;
     cto.fillFromInputVector(comps, RankFourTensor::symmetric_isotropic);
   }
  
   cto.rotate(trial_eigvecs);
  
   // drdsig = change in eigenvectors under a small stress change
   // drdsig(i,j,m,n) = dR(i,j)/dS_mn
   // The formula below is fairly easily derived:
   // S R = R s, so taking the variation
   // dS R + S dR = dR s + R ds, and multiplying by R^T
   // R^T dS R + R^T S dR = R^T dR s + ds .... (eqn 1)
   // I demand that RR^T = 1 = R^T R, and also that
   // (R+dR)(R+dR)^T = 1 = (R+dT)^T (R+dR), which means
   // that dR = R*c, for some antisymmetric c, so Eqn1 reads
   // R^T dS R + s c = c s + ds
   // Grabbing the components of this gives ds/dS (already
   // in RankTwoTensor), and c, which is:
   //   dR_ik/dS_mn = drdsig(i, k, m, n) = trial_eigvecs(m, b)*trial_eigvecs(n, k)*trial_eigvecs(i,
   //   b)/(trial_eigvals[k] - trial_eigvals[b]);
   // (sum over b!=k).
  
   RankFourTensor drdsig;
   for (unsigned k = 0; k < 3; ++k)
     for (unsigned b = 0; b < 3; ++b)
     {
       if (b == k)
         continue;
       for (unsigned m = 0; m < 3; ++m)
         for (unsigned n = 0; n < 3; ++n)
           for (unsigned i = 0; i < 3; ++i)
             drdsig(i, k, m, n) += trial_eigvecs(m, b) * trial_eigvecs(n, k) * trial_eigvecs(i, b) /
                                   (trial_eigvals[k] - trial_eigvals[b]);
     }
  
   // With diagla = diag(eigvals[0], eigvals[1], digvals[2])
   // The following implements
   // ans(i, j, a, b) += (drdsig(i, k, m, n)*trial_eigvecs(j, l)*diagla(k, l) + trial_eigvecs(i,
   // k)*drdsig(j, l, m, n)*diagla(k, l))*E_ijkl(m, n, a, b);
   // (sum over k, l, m and n)
  
   RankFourTensor ans;
   for (unsigned i = 0; i < 3; ++i)
     for (unsigned j = 0; j < 3; ++j)
       for (unsigned a = 0; a < 3; ++a)
         for (unsigned k = 0; k < 3; ++k)
           for (unsigned m = 0; m < 3; ++m)
             ans(i, j, a, a) += (drdsig(i, k, m, m) * trial_eigvecs(j, k) +
                                 trial_eigvecs(i, k) * drdsig(j, k, m, m)) *
                                eigvals[k] * la; // E_ijkl(m, n, a, b) = la*(m==n)*(a==b);
  
   for (unsigned i = 0; i < 3; ++i)
     for (unsigned j = 0; j < 3; ++j)
       for (unsigned a = 0; a < 3; ++a)
         for (unsigned b = 0; b < 3; ++b)
           for (unsigned k = 0; k < 3; ++k)
           {
             ans(i, j, a, b) += (drdsig(i, k, a, b) * trial_eigvecs(j, k) +
                                 trial_eigvecs(i, k) * drdsig(j, k, a, b)) *
                                eigvals[k] * mu; // E_ijkl(m, n, a, b) = mu*(m==a)*(n==b)
             ans(i, j, a, b) += (drdsig(i, k, b, a) * trial_eigvecs(j, k) +
                                 trial_eigvecs(i, k) * drdsig(j, k, b, a)) *
                                eigvals[k] * mu; // E_ijkl(m, n, a, b) = mu*(m==b)*(n==a)
           }
  
   return cto + ans;
 }

◆ 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

stress	the stress at which to calculate the flow potential
intnl	internal 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 133 of file TensorMechanicsPlasticModel.C.

 {
   return RankTwoTensor();
 }

Referenced by TensorMechanicsPlasticModel::dflowPotential_dintnlV().

◆ 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

	stress	the stress at which to calculate the flow potential
	intnl	internal parameter
[out]	dr_dintnl	dr_dintnl[alpha](i, j) = dr[alpha](i, j)/dintnl

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 141 of file TensorMechanicsPlasticTensileMulti.C.

 {
   dr_dintnl.assign(3, RankTwoTensor());
 }

◆ dflowPotential_dstress()

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

protectedvirtualinherited

The derivative of the flow potential with respect to stress.

Parameters

stress	the stress at which to calculate the flow potential
intnl	internal parameter

Returns: dr_dstress(i, j, k, l) = dr(i, j)/dstress(k, l)

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

Definition at line 119 of file TensorMechanicsPlasticModel.C.

 {
   return RankFourTensor();
 }

Referenced by TensorMechanicsPlasticModel::dflowPotential_dstressV().

◆ 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

	stress	the stress at which to calculate the flow potential
	intnl	internal parameter
[out]	dr_dstress	dr_dstress[alpha](i, j, k, l) = dr[alpha](i, j)/dstress(k, l)

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 134 of file TensorMechanicsPlasticTensileMulti.C.

 {
   stress.d2symmetricEigenvalues(dr_dstress);
 }

◆ 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

stress	the stress at which to calculate the hardening potentials
intnl	internal parameter

Returns: the derivative

Reimplemented in TensorMechanicsPlasticMeanCapTC.

Definition at line 174 of file TensorMechanicsPlasticModel.C.

 {
   return 0.0;
 }

Referenced by TensorMechanicsPlasticModel::dhardPotential_dintnlV().

◆ 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

	stress	the stress at which to calculate the hardening potentials
	intnl	internal parameter
[out]	dh_dintnl	dh_dintnl[alpha] = dh[alpha]/dintnl

Definition at line 180 of file TensorMechanicsPlasticModel.C.

 {
   dh_dintnl.resize(numberSurfaces(), dhardPotential_dintnl(stress, intnl));
 }

◆ dhardPotential_dstress()

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

protectedvirtualinherited

The derivative of the hardening potential with respect to stress.

Parameters

stress	the stress at which to calculate the hardening potentials
intnl	internal parameter

Returns: dh_dstress(i, j) = dh/dstress(i, j)

Reimplemented in TensorMechanicsPlasticMeanCapTC.

Definition at line 160 of file TensorMechanicsPlasticModel.C.

 {
   return RankTwoTensor();
 }

Referenced by TensorMechanicsPlasticModel::dhardPotential_dstressV().

◆ 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

	stress	the stress at which to calculate the hardening potentials
	intnl	internal parameter
[out]	dh_dstress	dh_dstress[alpha](i, j) = dh[alpha]/dstress(i, j)

Definition at line 166 of file TensorMechanicsPlasticModel.C.

 {
   dh_dstress.assign(numberSurfaces(), dhardPotential_dstress(stress, intnl));
 }

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

 {
   mooseAssert(dpm.size() == 3,
               "TensorMechanicsPlasticTensileMulti size of dpm should be 3 but it is "
                   << dpm.size());
  
   std::vector<Real> eigvals;
   RankTwoTensor eigvecs;
   trial_stress.symmetricEigenvaluesEigenvectors(eigvals, eigvecs);
   eigvals[0] += _shift;
   eigvals[2] -= _shift;
  
   Real str = tensile_strength(intnl_old);
  
   yf.resize(3);
   yf[0] = eigvals[0] - str;
   yf[1] = eigvals[1] - str;
   yf[2] = eigvals[2] - str;
  
   if (yf[0] <= _f_tol && yf[1] <= _f_tol && yf[2] <= _f_tol)
   {
     // purely elastic (trial_stress, intnl_old)
     trial_stress_inadmissible = false;
     return true;
   }
  
   trial_stress_inadmissible = true;
   delta_dp.zero();
   returned_stress.zero();
  
   // In the following i often assume that E_ijkl is
   // for an isotropic situation.  This reduces FLOPS
   // substantially which is important since the returnMap
   // is potentially the most compute-intensive function
   // of a simulation.
   // In many comments i write the general expression, and
   // i hope that might guide future coders if they are
   // generalising to a non-istropic E_ijkl
  
   // n[alpha] = E_ijkl*r[alpha]_kl expressed in principal stress space
   // (alpha = 0, 1, 2, corresponding to the three surfaces)
   // Note that in principal stress space, the flow
   // directions are, expressed in 'vector' form,
   // r[0] = (1,0,0), r[1] = (0,1,0), r[2] = (0,0,1).
   // Similar for _n:
   // so _n[0] = E_ij00*r[0], _n[1] = E_ij11*r[1], _n[2] = E_ij22*r[2]
   // In the following I assume that the E_ijkl is
   // for an isotropic situation.
   // In the anisotropic situation, we couldn't express
   // the flow directions as vectors in the same principal
   // stress space as the stress: they'd be full rank-2 tensors
   std::vector<RealVectorValue> n(3);
   n[0](0) = E_ijkl(0, 0, 0, 0);
   n[0](1) = E_ijkl(1, 1, 0, 0);
   n[0](2) = E_ijkl(2, 2, 0, 0);
   n[1](0) = E_ijkl(0, 0, 1, 1);
   n[1](1) = E_ijkl(1, 1, 1, 1);
   n[1](2) = E_ijkl(2, 2, 1, 1);
   n[2](0) = E_ijkl(0, 0, 2, 2);
   n[2](1) = E_ijkl(1, 1, 2, 2);
   n[2](2) = E_ijkl(2, 2, 2, 2);
  
   // With non-zero Poisson's ratio and hardening
   // it is not computationally cheap to know whether
   // the trial stress will return to the tip, edge,
   // or plane.  The following is correct for zero
   // Poisson's ratio and no hardening, and at least
   // gives a not-completely-stupid guess in the
   // more general case.
   // trial_order[0] = type of return to try first
   // trial_order[1] = type of return to try second
   // trial_order[2] = type of return to try third
   const unsigned int number_of_return_paths = 3;
   std::vector<int> trial_order(number_of_return_paths);
   if (yf[0] > _f_tol) // all the yield functions are positive, since eigvals are ordered eigvals[0]
                       // <= eigvals[1] <= eigvals[2]
   {
     trial_order[0] = tip;
     trial_order[1] = edge;
     trial_order[2] = plane;
   }
   else if (yf[1] > _f_tol) // two yield functions are positive
   {
     trial_order[0] = edge;
     trial_order[1] = tip;
     trial_order[2] = plane;
   }
   else
   {
     trial_order[0] = plane;
     trial_order[1] = edge;
     trial_order[2] = tip;
   }
  
   unsigned trial;
   bool nr_converged = false;
   for (trial = 0; trial < number_of_return_paths; ++trial)
   {
     switch (trial_order[trial])
     {
       case tip:
         nr_converged = returnTip(eigvals, n, dpm, returned_stress, intnl_old, 0);
         break;
       case edge:
         nr_converged = returnEdge(eigvals, n, dpm, returned_stress, intnl_old, 0);
         break;
       case plane:
         nr_converged = returnPlane(eigvals, n, dpm, returned_stress, intnl_old, 0);
         break;
     }
  
     str = tensile_strength(intnl_old + dpm[0] + dpm[1] + dpm[2]);
  
     if (nr_converged && KuhnTuckerOK(returned_stress, dpm, str, ep_plastic_tolerance))
       break;
   }
  
   if (trial == number_of_return_paths)
   {
     Moose::err << "Trial stress = \n";
     trial_stress.print(Moose::err);
     Moose::err << "Internal parameter = " << intnl_old << "\n";
     mooseError("TensorMechanicsPlasticTensileMulti: FAILURE!  You probably need to implement a "
                "line search\n");
     // failure - must place yield function values at trial stress into yf
     str = tensile_strength(intnl_old);
     yf[0] = eigvals[0] - str;
     yf[1] = eigvals[1] - str;
     yf[2] = eigvals[2] - str;
     return false;
   }
  
   // success
  
   returned_intnl = intnl_old;
   for (unsigned i = 0; i < 3; ++i)
   {
     yf[i] = returned_stress(i, i) - str;
     delta_dp(i, i) = dpm[i];
     returned_intnl += dpm[i];
   }
   returned_stress = eigvecs * returned_stress * (eigvecs.transpose());
   delta_dp = eigvecs * delta_dp * (eigvecs.transpose());
   return true;
 }

Referenced by activeConstraints(), and returnMap().

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

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

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

 {
   return _strength.derivative(internal_param);
 }

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

◆ 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

stress	the stress at which to calculate the yield function
intnl	internal parameter

Returns: the derivative

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

Definition at line 92 of file TensorMechanicsPlasticModel.C.

 {
   return 0.0;
 }

Referenced by TensorMechanicsPlasticModel::dyieldFunction_dintnlV().

◆ 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

	stress	the stress at which to calculate the yield function
	intnl	internal parameter
[out]	df_dintnl	df_dintnl[alpha] = df[alpha]/dintnl

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 117 of file TensorMechanicsPlasticTensileMulti.C.

 {
   df_dintnl.assign(3, -dtensile_strength(intnl));
 }

◆ dyieldFunction_dstress()

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

protectedvirtualinherited

The derivative of yield function with respect to stress.

Parameters

stress	the stress at which to calculate the yield function
intnl	internal 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 77 of file TensorMechanicsPlasticModel.C.

 {
   return RankTwoTensor();
 }

Referenced by TensorMechanicsPlasticModel::dyieldFunction_dstressV().

◆ 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

	stress	the stress at which to calculate the yield function
	intnl	internal parameter
[out]	df_dstress	df_dstress[alpha](i, j) = dyieldFunction[alpha]/dstress(i, j)

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 92 of file TensorMechanicsPlasticTensileMulti.C.

 {
   std::vector<Real> eigvals;
   stress.dsymmetricEigenvalues(eigvals, df_dstress);
  
   if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
   {
     Real small_perturbation;
     RankTwoTensor shifted_stress = stress;
     while (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
     {
       for (unsigned i = 0; i < 3; ++i)
         for (unsigned j = 0; j <= i; ++j)
         {
           small_perturbation = 0.1 * _shift * 2 * (MooseRandom::rand() - 0.5);
           shifted_stress(i, j) += small_perturbation;
           shifted_stress(j, i) += small_perturbation;
         }
       shifted_stress.dsymmetricEigenvalues(eigvals, df_dstress);
     }
   }
 }

Referenced by flowPotentialV().

◆ execute()

void TensorMechanicsPlasticModel::execute ( )

inherited

Definition at line 47 of file TensorMechanicsPlasticModel.C.

48 {

49 }

◆ finalize()

void TensorMechanicsPlasticModel::finalize ( )

inherited

Definition at line 52 of file TensorMechanicsPlasticModel.C.

53 {

54 }

◆ flowPotential()

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

protectedvirtualinherited

The flow potential.

Parameters

stress	the stress at which to calculate the flow potential
intnl	internal parameter

Returns: the flow potential

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

Definition at line 106 of file TensorMechanicsPlasticModel.C.

 {
   return RankTwoTensor();
 }

Referenced by TensorMechanicsPlasticModel::flowPotentialV().

◆ flowPotentialV()

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

overridevirtual

The flow potentials.

Parameters

	stress	the stress at which to calculate the flow potential
	intnl	internal parameter
[out]	r	r[alpha] is the flow potential for the "alpha" yield function

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 125 of file TensorMechanicsPlasticTensileMulti.C.

 {
   // This plasticity is associative so
   dyieldFunction_dstressV(stress, intnl, r);
 }

◆ hardPotential()

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

protectedvirtualinherited

The hardening potential.

Parameters

stress	the stress at which to calculate the hardening potential
intnl	internal parameter

Returns: the hardening potential

Reimplemented in TensorMechanicsPlasticMeanCapTC.

Definition at line 147 of file TensorMechanicsPlasticModel.C.

 {
   return -1.0;
 }

Referenced by TensorMechanicsPlasticModel::hardPotentialV().

◆ hardPotentialV()

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

virtualinherited

The hardening potential.

Parameters

	stress	the stress at which to calculate the hardening potential
	intnl	internal parameter
[out]	h	h[alpha] is the hardening potential for the "alpha" yield function

Definition at line 152 of file TensorMechanicsPlasticModel.C.

 {
   h.assign(numberSurfaces(), hardPotential(stress, intnl));
 }

◆ initialize()

void TensorMechanicsPlasticModel::initialize ( )

inherited

Definition at line 42 of file TensorMechanicsPlasticModel.C.

43 {

44 }

◆ 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_stress	The eigenvalues (sorted in ascending order as is standard in this Class) are stored in the diagonal components
dpm	The three plastic multipliers
str	The yield strength
ep_plastic_tolerance	The 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 602 of file TensorMechanicsPlasticTensileMulti.C.

 {
   for (unsigned i = 0; i < 3; ++i)
     if (!TensorMechanicsPlasticModel::KuhnTuckerSingleSurface(
             returned_diagonal_stress(i, i) - str, dpm[i], ep_plastic_tolerance))
       return false;
   return true;
 }

Referenced by doReturnMap().

◆ 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

yf	Yield function value
dpm	plastic multiplier
dpm_tol	tolerance on plastic multiplier: viz dpm>-dpm_tol means "dpm is non-negative"

Definition at line 248 of file TensorMechanicsPlasticModel.C.

 {
   return (dpm == 0 && yf <= _f_tol) || (dpm > -dpm_tol && yf <= _f_tol && yf >= -_f_tol);
 }

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

◆ modelName()

std::string TensorMechanicsPlasticTensileMulti::modelName ( ) const

overridevirtual

Implements TensorMechanicsPlasticModel.

Definition at line 212 of file TensorMechanicsPlasticTensileMulti.C.

 {
   return "TensileMulti";
 }

◆ numberSurfaces()

unsigned int TensorMechanicsPlasticTensileMulti::numberSurfaces ( ) const

overridevirtual

The number of yield surfaces for this plasticity model.

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 71 of file TensorMechanicsPlasticTensileMulti.C.

 {
   return 3;
 }

◆ 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

eigvals	The three stress eigenvalues, sorted in ascending order
n	The three return directions, n=E_ijkl*r. Note this algorithm assumes isotropic elasticity, so these are 3 vectors in principal stress space
dpm	The 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_stress	The returned stress. This will be diagonal, with the return-mapped eigenvalues in the diagonal positions, sorted in ascending order
intnl_old	The 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_guess	A guess of dpm[1]+dpm[2]

Definition at line 487 of file TensorMechanicsPlasticTensileMulti.C.

 {
   // work out the point to which we would return, "a".  It is defined by
   // f1 = 0 = f2, and the normality condition:
   //   (eigvals - a).(n1 x n2) = 0,
   // where eigvals is the starting position
   // (it is a vector in principal stress space).
   // To get f1=0=f2, we need a = (a0, str, str), and a0 is found
   // by expanding the normality condition to yield:
   //   a0 = (-str*n1xn2[1] - str*n1xn2[2] + edotn1xn2)/n1xn2[0];
   // where edotn1xn2 = eigvals.(n1 x n2)
   //
   // We need to find the plastic multipliers, dpm, defined by
   //   eigvals - a = dpm[1]*n1 + dpm[2]*n2
   // For the isotropic case, and defining eminusa = eigvals - a,
   // the solution is easy:
   //   dpm[0] = 0;
   //   dpm[1] = (eminusa[1] - eminusa[0])/(n[1][1] - n[1][0]);
   //   dpm[2] = (eminusa[2] - eminusa[0])/(n[2][2] - n[2][0]);
   //
   // Now specialise to the isotropic case.  Define
   //   x = dpm[1] + dpm[2] = (eigvals[1] + eigvals[2] - 2*str)/(n[0][0] + n[0][1])
   // Notice that the RHS is a function of x, so we solve using
   // Newton-Raphson starting with x=initial_guess
   Real x = initial_guess;
   const Real denom = n[0](0) + n[0](1);
   Real str = tensile_strength(intnl_old + x);
  
   if (_strength.modelName().compare("Constant") != 0)
   {
     // Finding x is expensive.  Therefore
     // although x!=0 for Constant Hardening, solution
     // for dpm above demonstrates that we don't
     // actually need to know x to find the dpm for
     // Constant Hardening.
     //
     // However, for nontrivial Hardening, the following
     // is necessary
     const Real eig = eigvals[1] + eigvals[2];
     Real residual = denom * x - eig + 2 * str;
     Real jacobian;
     unsigned int iter = 0;
     do
     {
       jacobian = denom + 2 * dtensile_strength(intnl_old + x);
       x += -residual / jacobian;
       if (iter > _max_iters) // not converging
         return false;
       str = tensile_strength(intnl_old + x);
       residual = denom * x - eig + 2 * str;
       iter++;
     } while (residual * residual > _f_tol * _f_tol);
   }
  
   dpm[0] = 0;
   dpm[1] = ((eigvals[1] * n[0](0) - eigvals[2] * n[0](1)) / (n[0](0) - n[0](1)) - str) / denom;
   dpm[2] = ((eigvals[2] * n[0](0) - eigvals[1] * n[0](1)) / (n[0](0) - n[0](1)) - str) / denom;
  
   returned_stress(0, 0) = eigvals[0] - n[0](1) * (dpm[1] + dpm[2]);
   returned_stress(1, 1) = returned_stress(2, 2) = str;
   return true;
 }

Referenced by doReturnMap().

◆ 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_stress	The trial stress
	intnl_old	Value of the internal parameter
	E_ijkl	Elasticity tensor
	ep_plastic_tolerance	Tolerance defined by the user for the plastic strain
[out]	returned_stress	In case (C): lies on the yield surface after returning and produces the correct plastic strain (normality condition). Otherwise: not defined
[out]	returned_intnl	In case (C): the value of the internal parameter after returning. Otherwise: not defined
[out]	dpm	In case (C): the plastic multipliers needed to bring about the return. Otherwise: not defined
[out]	delta_dp	In case (C): The change in plastic strain induced by the return process. Otherwise: not defined
[out]	yf	In case (C): the yield function at (returned_stress, returned_intnl). Otherwise: the yield function at (trial_stress, intnl_old)
[out]	trial_stress_inadmissible	Should 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 218 of file TensorMechanicsPlasticTensileMulti.C.

 {
   if (!_use_custom_returnMap)
     return TensorMechanicsPlasticModel::returnMap(trial_stress,
                                                   intnl_old,
                                                   E_ijkl,
                                                   ep_plastic_tolerance,
                                                   returned_stress,
                                                   returned_intnl,
                                                   dpm,
                                                   delta_dp,
                                                   yf,
                                                   trial_stress_inadmissible);
  
   return doReturnMap(trial_stress,
                      intnl_old,
                      E_ijkl,
                      ep_plastic_tolerance,
                      returned_stress,
                      returned_intnl,
                      dpm,
                      delta_dp,
                      yf,
                      trial_stress_inadmissible);
 }

◆ 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

eigvals	The three stress eigenvalues, sorted in ascending order
n	The three return directions, n=E_ijkl*r. Note this algorithm assumes isotropic elasticity, so these are 3 vectors in principal stress space
dpm	The 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_stress	The returned stress. This will be diagonal, with the return-mapped eigenvalues in the diagonal positions, sorted in ascending order
intnl_old	The 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_guess	A guess of dpm[2]

Definition at line 556 of file TensorMechanicsPlasticTensileMulti.C.

 {
   // the returned point, "a", is defined by f2=0 and
   // a = p - dpm[2]*n2.
   // This is a vector equation in
   // principal stress space, and dpm[2] is the third
   // plasticity multiplier (dpm[0]=0=dpm[1] for return
   // to the plane) and "p" is the starting
   // position (p=eigvals).
   // (Kuhn-Tucker demands that dpm[2]>=0, but we leave checking
   // that condition for later.)
   // Since f2=0, we must have a[2]=tensile_strength,
   // so we can just look at the [2] component of the
   // equation, which yields
   // n[2][2]*dpm[2] - eigvals[2] + str = 0
   // For hardening, str=tensile_strength(intnl_old+dpm[2]),
   // and we want to solve for dpm[2].
   // Use Newton-Raphson with initial guess dpm[2] = initial_guess
   dpm[2] = initial_guess;
   Real residual = n[2](2) * dpm[2] - eigvals[2] + tensile_strength(intnl_old + dpm[2]);
   Real jacobian;
   unsigned int iter = 0;
   do
   {
     jacobian = n[2](2) + dtensile_strength(intnl_old + dpm[2]);
     dpm[2] += -residual / jacobian;
     if (iter > _max_iters) // not converging
       return false;
     residual = n[2](2) * dpm[2] - eigvals[2] + tensile_strength(intnl_old + dpm[2]);
     iter++;
   } while (residual * residual > _f_tol * _f_tol);
  
   dpm[0] = 0;
   dpm[1] = 0;
   returned_stress(0, 0) = eigvals[0] - dpm[2] * n[2](0);
   returned_stress(1, 1) = eigvals[1] - dpm[2] * n[2](1);
   returned_stress(2, 2) = eigvals[2] - dpm[2] * n[2](2);
   return true;
 }

Referenced by doReturnMap().

◆ 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

eigvals	The three stress eigenvalues, sorted in ascending order
n	The three return directions, n=E_ijkl*r. Note this algorithm assumes isotropic elasticity, so these are 3 vectors in principal stress space
dpm	The 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_stress	The returned stress. This will be diagonal, with the return-mapped eigenvalues in the diagonal positions, sorted in ascending order
intnl_old	The 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_guess	A guess of dpm[0]+dpm[1]+dpm[2]

Definition at line 410 of file TensorMechanicsPlasticTensileMulti.C.

 {
   // The returned point is defined by f0=f1=f2=0.
   // that is, returned_stress = diag(str, str, str), where
   // str = tensile_strength(intnl),
   // where intnl = intnl_old + dpm[0] + dpm[1] + dpm[2]
   // The 3 plastic multipliers, dpm, are defiend by the normality condition
   //   eigvals - str = dpm[0]*n[0] + dpm[1]*n[1] + dpm[2]*n[2]
   // (Kuhn-Tucker demands that all dpm are non-negative, but we leave
   // that checking for later.)
   // This is a vector equation with solution (A):
   //   dpm[0] = triple(eigvals - str, n[1], n[2])/trip;
   //   dpm[1] = triple(eigvals - str, n[2], n[0])/trip;
   //   dpm[2] = triple(eigvals - str, n[0], n[1])/trip;
   // where trip = triple(n[0], n[1], n[2]).
   // By adding the three components of that solution together
   // we can get an equation for x = dpm[0] + dpm[1] + dpm[2],
   // and then our Newton-Raphson only involves one variable (x).
   // In the following, i specialise to the isotropic situation.
  
   Real x = initial_guess;
   const Real denom = (n[0](0) - n[0](1)) * (n[0](0) + 2 * n[0](1));
   Real str = tensile_strength(intnl_old + x);
  
   if (_strength.modelName().compare("Constant") != 0)
   {
     // Finding x is expensive.  Therefore
     // although x!=0 for Constant Hardening, solution (A)
     // demonstrates that we don't
     // actually need to know x to find the dpm for
     // Constant Hardening.
     //
     // However, for nontrivial Hardening, the following
     // is necessary
     const Real eig = eigvals[0] + eigvals[1] + eigvals[2];
     const Real bul = (n[0](0) + 2 * n[0](1));
  
     // and finally, the equation we want to solve is:
     // bul*x - eig + 3*str = 0
     // where str=tensile_strength(intnl_old + x)
     // and x = dpm[0] + dpm[1] + dpm[2]
     // (Note this has units of stress, so using _f_tol as a convergence check is reasonable.)
     // Use Netwon-Raphson with initial guess x = 0
     Real residual = bul * x - eig + 3 * str;
     Real jacobian;
     unsigned int iter = 0;
     do
     {
       jacobian = bul + 3 * dtensile_strength(intnl_old + x);
       x += -residual / jacobian;
       if (iter > _max_iters) // not converging
         return false;
       str = tensile_strength(intnl_old + x);
       residual = bul * x - eig + 3 * str;
       iter++;
     } while (residual * residual > _f_tol * _f_tol);
   }
  
   // The following is the solution (A) written above
   // (dpm[0] = triple(eigvals - str, n[1], n[2])/trip, etc)
   // in the isotropic situation
   dpm[0] = (n[0](0) * (eigvals[0] - str) + n[0](1) * (eigvals[0] - eigvals[1] - eigvals[2] + str)) /
            denom;
   dpm[1] = (n[0](0) * (eigvals[1] - str) + n[0](1) * (eigvals[1] - eigvals[2] - eigvals[0] + str)) /
            denom;
   dpm[2] = (n[0](0) * (eigvals[2] - str) + n[0](1) * (eigvals[2] - eigvals[0] - eigvals[1] + str)) /
            denom;
   returned_stress(0, 0) = returned_stress(1, 1) = returned_stress(2, 2) = str;
   return true;
 }

Referenced by doReturnMap().

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

 {
   return _strength.value(internal_param);
 }

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

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

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

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

 {
   return _use_custom_cto;
 }

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

 {
   return _use_custom_returnMap;
 }

◆ validParams()

InputParameters TensorMechanicsPlasticTensileMulti::validParams ( )

static

Definition at line 21 of file TensorMechanicsPlasticTensileMulti.C.

 {
   InputParameters params = TensorMechanicsPlasticModel::validParams();
   params.addClassDescription("Associative tensile plasticity with hardening/softening");
   params.addRequiredParam<UserObjectName>(
       "tensile_strength",
       "A TensorMechanicsHardening UserObject that defines hardening of the tensile strength");
   params.addParam<Real>("shift",
                         "Yield surface is shifted by this amount to avoid problems with "
                         "defining derivatives when eigenvalues are equal.  If this is "
                         "larger than f_tol, a warning will be issued.  Default = f_tol.");
   params.addParam<unsigned int>("max_iterations",
                                 10,
                                 "Maximum number of Newton-Raphson iterations "
                                 "allowed in the custom return-map algorithm. "
                                 " For highly nonlinear hardening this may "
                                 "need to be higher than 10.");
   params.addParam<bool>("use_custom_returnMap",
                         true,
                         "Whether to use the custom returnMap "
                         "algorithm.  Set to true if you are using "
                         "isotropic elasticity.");
   params.addParam<bool>("use_custom_cto",
                         true,
                         "Whether to use the custom consistent tangent "
                         "operator computations.  Set to true if you are "
                         "using isotropic elasticity.");
   return params;
 }

◆ 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

stress	the stress at which to calculate the yield function
intnl	internal parameter

Returns: the yield function

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

Definition at line 63 of file TensorMechanicsPlasticModel.C.

 {
   return 0.0;
 }

Referenced by TensorMechanicsPlasticModel::yieldFunctionV().

◆ 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

	stress	the stress at which to calculate the yield function
	intnl	internal parameter
[out]	f	the yield functions

Reimplemented from TensorMechanicsPlasticModel.

Definition at line 77 of file TensorMechanicsPlasticTensileMulti.C.

 {
   std::vector<Real> eigvals;
   stress.symmetricEigenvalues(eigvals);
   const Real str = tensile_strength(intnl);
  
   f.resize(3);
   f[0] = eigvals[0] + _shift - str;
   f[1] = eigvals[1] - str;
   f[2] = eigvals[2] - _shift - str;
 }

Member Data Documentation

◆ _f_tol

const Real TensorMechanicsPlasticModel::_f_tol

inherited

Tolerance on yield function.

Definition at line 175 of file TensorMechanicsPlasticModel.h.

◆ _ic_tol

const Real TensorMechanicsPlasticModel::_ic_tol

inherited

Tolerance on internal constraint.

Definition at line 178 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 99 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 102 of file TensorMechanicsPlasticTensileMulti.h.

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

◆ _strength

const TensorMechanicsHardeningModel& TensorMechanicsPlasticTensileMulti::_strength

private

Definition at line 96 of file TensorMechanicsPlasticTensileMulti.h.

Referenced by dtensile_strength(), returnEdge(), returnTip(), and tensile_strength().

◆ _use_custom_cto

const bool TensorMechanicsPlasticTensileMulti::_use_custom_cto

private

Whether to use the custom consistent tangent operator calculation.

Definition at line 108 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 105 of file TensorMechanicsPlasticTensileMulti.h.

Referenced by returnMap(), and useCustomReturnMap().

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

tensor_mechanics/include/userobjects/TensorMechanicsPlasticTensileMulti.h
tensor_mechanics/src/userobjects/TensorMechanicsPlasticTensileMulti.C

Public Member Functions

Static Public Member Functions

Public Attributes

Protected Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Enumeration Documentation

◆ return_type

Constructor & Destructor Documentation

◆ TensorMechanicsPlasticTensileMulti()

Member Function Documentation

◆ activeConstraints()

◆ consistentTangentOperator()

◆ dflowPotential_dintnl()

◆ dflowPotential_dintnlV()

◆ dflowPotential_dstress()

◆ dflowPotential_dstressV()

◆ dhardPotential_dintnl()

◆ dhardPotential_dintnlV()

◆ dhardPotential_dstress()

◆ dhardPotential_dstressV()

◆ doReturnMap()

◆ dot()

◆ dtensile_strength()

◆ dyieldFunction_dintnl()

◆ dyieldFunction_dintnlV()

◆ dyieldFunction_dstress()

◆ dyieldFunction_dstressV()

◆ execute()

◆ finalize()

◆ flowPotential()

◆ flowPotentialV()

◆ hardPotential()

◆ hardPotentialV()

◆ initialize()

◆ KuhnTuckerOK()

◆ KuhnTuckerSingleSurface()

◆ modelName()

◆ numberSurfaces()

◆ returnEdge()

◆ returnMap()

◆ returnPlane()

◆ returnTip()

◆ tensile_strength()

◆ triple()

◆ useCustomCTO()

◆ useCustomReturnMap()

◆ validParams()

◆ yieldFunction()

◆ yieldFunctionV()

Member Data Documentation

◆ _f_tol

◆ _ic_tol

◆ _max_iters

◆ _shift

◆ _strength

◆ _use_custom_cto

◆ _use_custom_returnMap