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

virtual void computeStress(double *stress, double *dStress_dStrain, const double *dStrain, const double *timeOld, const double dT, double &pNewDT) 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:

• stress – [inout] Cauchy stress

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

• dStrain – [in] linearized strain increment

• timeOld – [in] Old (pseudo-)time

• dt – [in] (Pseudo-)time increment from the old (pseudo-)time to the current (pseudo-)time

• pNewDT – [inout] Suggestion for a new time increment

virtual double getDensity() override

Get material density.

Throws:

std::runtime_error – if density is not defined.

Returns:

Density value.

inline virtual int getNumberOfRequiredStateVars() override

Returns:

Number of state variables required by the material.

virtual void assignStateVars(double *stateVars, int nStateVars) override

Assign state variable array to material.

Parameters:

• stateVars – [inout] Pointer to state variable array.

• nStateVars – [in] Number of state variables.

virtual StateView getStateView(const std::string &result) override

Access material state variables by name.

Parameters:

stateName – [in] Name of the requested state variable.

Returns:

A view into the state variable array.

Public Members

std::unique_ptr<VonMisesModelStateVarManager> managedStateVars

class VonMisesModelStateVarManager : public MarmotStateVarVectorManager

Inheritance diagram for Marmot::Materials::VonMisesModel::VonMisesModelStateVarManager:

Collaboration diagram for Marmot::Materials::VonMisesModel::VonMisesModelStateVarManager:

Public Functions

inline VonMisesModelStateVarManager(double *theStateVarVector)

Public Members

double &kappa

Hardening parameter.

Public Static Attributes

static const auto layout = makeLayout( {{ .name = "kappa", .length = 1 },} )