Von Mises model

An implementation of classical J2 plasticity with isotropic hardening.

Index	Model Parameter	Description
0	\(E\)	Young’s modulus
1	\(\nu\)	Poisson’s ratio
2	\(f_\mathrm{y}^{0}\)	Initial yield stress
3	\(H_\mathrm{iso,lin}\)	First hardening parameter [1]
4	\(\Delta f_\mathrm{y}^{0,\infty}\)	Second hardening parameter [1]
5	\(\delta\)	Third hardening parameter [1]

State Variable	Name	Description
\(\kappa\)	kappa	Hardening variable [1]

[1] (1,2,3,4)

see Hardening Law.

Theory

Elastoplastic Constitutive Relations

The infinitesimal elastoplastic constitutive relations, which relate the stress and linearized strain rate tensors by means of the elastic continuum tangent operator, are given as

\[\dot{\sigma}_{ij} = C^\mathrm{e}_{ijkl}\dot{\varepsilon}^\mathrm{e}_{kl} = C^\mathrm{e}_{ijkl}\left(\dot{\varepsilon}_{kl}-\dot{\varepsilon}^\mathrm{p}_{kl}\right).\]

Therein, an additive split of the strain rate tensor in elastic and plastic parts is assumed.

The elastic continuum tangent operator is given in terms of Young’s modulus \(E\) and Poisson’s ratio \(\nu\) as

\[C^\mathrm{e}_{ijkl} = \frac{E\nu\,\delta_{ij}\delta_{kl}}{\left(1+\nu\right)\left(1-2\nu\right)} + \frac{E\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)}{2\left(1+\nu\right)}\]

Alternatively, the constitutive relations can be rewritten based on a volumetric-deviatoric split of the stress and strain rate tensors

\[\dot{\sigma}_{ij} = \dot{p}\delta_{ij}+\dot{s}_{ij} \qquad\text{with}\quad \dot{p}=\frac{\dot{\sigma}_{ii}}{3}\]

and

\[\dot{\varepsilon}_{ij} = \frac{\dot{\varepsilon}_\mathrm{vol}}{3}\delta_{ij}+\dot{e}_{ij} = \left( \dot{\varepsilon}^\mathrm{e}_\mathrm{vol}+\dot{\varepsilon}^\mathrm{p}_\mathrm{vol} \right)\frac{\delta_{ij}}{3} + \dot{e}^\mathrm{e}_{ij}+\dot{e}^\mathrm{p}_{ij} \qquad\text{with}\quad \dot{\varepsilon}_\mathrm{vol}=\dot{\varepsilon}_{ii}\]

as

\[\dot{\sigma}_{ij} = \dot{p}\delta_{ij}+\dot{s}_{ij} = K \dot{\varepsilon}^\mathrm{e}_\mathrm{vol} + 2G \dot{e}^\mathrm{e}_{ij}\]

Yield Function

\[s_{ij} = \sigma_{ij}-\frac{1}{3}\sigma_{kk},\]

The yield function, which is formulated in terms of the deviatoric part of the stress tensor is given as

\[f\left(\boldsymbol{\sigma},\kappa\right) = \sqrt{s_{ij}s_{ij}}-\sqrt{\frac{2}{3}}f_\mathrm{y}\left(\kappa\right),\]

where \(f_\mathrm{y}\left(\kappa\right)\) is the yield stress under uniaxial tension, which may depend on the hardening variable \(\kappa\).

Flow Rule

The model encompasses associated plastic flow, i.e., the plastic strain rate

\[\dot{\varepsilon}^\mathrm{p}_{ij} = \dot{\lambda}\frac{\partial{}f}{\partial{}\sigma_{ij}} = \dot{\lambda}\frac{s_{ij}}{\sqrt{s_{ij}s_{ij}}} = \dot{\lambda}\boldsymbol{n}\]

is proportional to the partial derivatives of the yield function with respect to the stress tensor. The plastic multiplier \(\dot{\lambda}\) is a proportionality constant.

Hardening Law

The model includes non-linear isotropic hardening, which is controlled by the evolution of the yield stress as

\[f_\mathrm{y}\left(\kappa\right) = f_\mathrm{y}^0 + H_\mathrm{iso,lin}\kappa + \Delta f_\mathrm{y}^{0,\infty}\left(1-\exp\left(-\delta\kappa\right)\right).\]

Therein, \(f_\mathrm{y}^0\), \(H_\mathrm{iso,lin}\), \(\Delta{}f_\mathrm{y}^{0,\infty}\) and \(\delta\) are material parameters.

The evolution of the hardening variable is given as

\[\dot{\kappa} = \sqrt{\frac{2}{3}\dot{\varepsilon}^\mathrm{p}_{ij}\dot{\varepsilon}^\mathrm{p}_{ij}} = \dot{\lambda}\sqrt{\frac{2}{3}}.\]

Therein, the factor \(\sqrt{\frac{2}{3}}\) ensures that \(\dot{\kappa}\) is equal to \(\dot{\varepsilon}^\mathrm{p}_{11}\) for uniaxial loading in \(x_1\) direction.

Loading-Unloading Conditions, Consistency Condition

The plastic multiplier and the yield function obey the Karush–Kuhn–Tucker conditions

\[\dot{\lambda}f=0, \quad \dot{\lambda}\geq0, \quad f\leq0,\]

as well as the consistency condition

\[\dot{\lambda}\dot{f}=0.\]

Implementation

The model is implemented using an (implicit) backward Euler integration scheme. Given the stress tensor \(\boldsymbol{\sigma}^{n}\) and the value of the state variable \(\kappa^{n}\) at time \(t^{n}\) as well as an increment of the linearized strain tensor \(\Delta\boldsymbol{\varepsilon}=\boldsymbol{\varepsilon}^{n+1}-\boldsymbol{\varepsilon}^{n}\), an updated stress tensor \(\boldsymbol{\sigma}^{n+1}\) as well as an updated value of the state variable \(\kappa^{n+1}\) at time \(t^{n+1}\) are calculated by means of the stress update algorithm outlined below.

Stress Update Algorithm

Given the strain increment \(\Delta\boldsymbol{\varepsilon}\), an elastic trial stress is calculated as

\[\boldsymbol{\sigma}^\mathrm{trial} = \boldsymbol{\sigma}^n + \boldsymbol{C}^\mathrm{e}:\Delta\boldsymbol{\varepsilon}.\]

If the trial stress does not satisfy the yield condition, i.e., \(f\left(\boldsymbol{\boldsymbol{\sigma}^\mathrm{trial}},\kappa_n\right)<0\), the step is treated as elastic, which means that the updated stress is equal to the trial stress

\[\boldsymbol{\sigma}^{n+1} = \boldsymbol{\sigma}^\mathrm{trial}\]

and the hardening variable remains constant

\[\kappa_{n+1} = \kappa_n.\]

If the trial stress satisfies the yield condition, i.e., \(f\left(\boldsymbol{\boldsymbol{\sigma}^\mathrm{trial}},\kappa_n\right)\geq0\), the step is treated as elastoplastic, triggering the radial return mapping algorithm.

The updated deviatoric stress \(\boldsymbol{s}^{n+1}\) can be written as

\[\boldsymbol{s}^{n+1} = \boldsymbol{s}^{\mathrm{trial}} - 2G\Delta\boldsymbol{e}^\mathrm{p} = \boldsymbol{s}^{\mathrm{trial}} - 2G\Delta\lambda\boldsymbol{n}^{n+1}.\]

Taking the norm on both sides results in

\[\|\boldsymbol{s}^{n+1}\| = \|\boldsymbol{s}^{\mathrm{trial}}\| - 2G\Delta\lambda.\]

Inserting the resulting expression for \(\|\boldsymbol{s}^{n+1}\|\) into the yield function \(f^{n+1}\)

\[f^{n+1} = \|\boldsymbol{s}^{n+1}\| -\sqrt{\frac{2}{3}}f_\mathrm{y}\left(\kappa^{n+1}\right) = 0\]

and substituting \(\Delta\lambda\) with \(\sqrt{\frac{3}{2}}\Delta\kappa\) results in the scalar expression

\[g\left(\Delta\kappa\right) = \|\boldsymbol{s}^\mathrm{trial}\| -\sqrt{6}G\Delta\kappa -\sqrt{\frac{2}{3}}f_\mathrm{y}\left(\Delta\kappa\right) = 0,\]

which can be solved numerically to obtain \(\Delta\kappa\).

If \(\Delta\kappa\) is known, the state variable is updated as

\[\kappa^{n+1} = \kappa^{n} + \Delta\kappa\]

and the stress update is performed by exploiting \(\Delta\lambda=\sqrt{\frac{3}{2}}\Delta\kappa\)

\[\boldsymbol{\sigma}^{n+1} = \boldsymbol{\sigma}^\mathrm{trial} - 2G\Delta\lambda\frac{\boldsymbol{s}^{\mathrm{trial}}}{\|\boldsymbol{s}^{\mathrm{trial}}\|}\]

Consistent Algorithmic Tangent Operator

The consistent tangent operator is given as

\[\frac{\partial{\boldsymbol{s}^{n+1}}}{\partial\boldsymbol{\varepsilon}^{n+1}} = \mathbf{C}^\mathrm{e} - 2G\left( \left( 1+\frac{\partial f(\kappa^{n+1})}{\partial\kappa} \right)^{-1} -\frac{2G\Delta\lambda^{n+1}}{\|\boldsymbol{s}^\mathrm{trial}\|} \right) \boldsymbol{n}^{\mathrm{trial}}\otimes\boldsymbol{n}^{\mathrm{trial}} - \frac{4G^2\Delta\lambda}{\|\boldsymbol{s}^\mathrm{trial}\|}\boldsymbol{I}^{\mathrm{dev}}.\]

class VonMisesModel : public MarmotMaterialHypoElastic

Inheritance diagram for Marmot::Materials::VonMisesModel:

Collaboration diagram for Marmot::Materials::VonMisesModel:

An implementation of classical J2 plasticity with isotropic hardening.

Public Functions

VonMisesModel(const double *materialProperties, const int nMaterialProperties, const int materialLabel)

virtual void computeStress(state3D &state, Marmot::Matrix6d &dStressDDStrain, const Marmot::Vector6d &dStrain, const timeInfo &timeInfo) const override

For a given linearized strain increment \(\Delta\boldsymbol{\varepsilon}\) at the old and the current time, compute the Cauchy stress and the algorithmic tangent \(\frac{\partial\boldsymbol{\sigma}^{(n+1)}}{\partial\boldsymbol{\varepsilon}^{(n+1)}}\).

Parameters:

• state – [inout] A state3D instance carrying stress, strain energy, and state variables

• dStress_dStrain – [inout] Algorithmic tangent representing the derivative of the Cauchy stress tensor with respect to the linearized strain

• dStrain – [in] linearized strain increment

• timeInfo – Structure carrying the current (pseudo-)time and the (pseudo-)time increment

virtual void computeStressExplicit(state3D &state, const Marmot::Vector6d &dStrain, const timeInfo &timeInfo) const override

Explicit version of computeStress for use in explicit time integration schemes. The algorithmic tangent is not needed in explicit schemes and will therefore not be computed.

Note

The default implementation calls computeStress and ignores the algorithmic tangent.

Note

Derived classes may override this method for efficiency reasons.

Parameters:

• state – [inout] A state3D instance carrying stress, strain energy, and state variables

• dStrain – [in] linearized strain increment

• timeInfo – Structure carrying time information

virtual double getDensity(const double *stateVars) const override

Get the mass density of the material.

Parameters:

stateVars – Pointer to the state variable array

Returns:

Mass density of the material