# General formulation for local integration in standard elastoplasticity with an arbitrary hardening model.

Anh Le Van, Géry de Saxcé, Philippe Le Grognec

► To cite this version:

Anh Le Van, Géry de Saxcé, Philippe Le Grognec. General formulation for local integration in standard elastoplasticity with an arbitrary hardening model. Computers & Structures, 2003, 81 (21), pp.2099-2109. 10.1016/S0045-7949(03)00256-6. hal-01006917

HAL Id: hal-01006917

https://hal.science/hal-01006917

Submitted on 6 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# General formulation for local integration in standard elastoplasticity with an arbitrary hardening model

Anh Le van<sup>a</sup>, Géry de Saxcé<sup>b</sup>, Philippe Le Grognec<sup>a</sup>

This paper describes a general method for deriving the plastic corrections and the consistent tangent modulus for a wide range of arbitrary non-linear hardening models within the framework of standard small strains elastoplasticity. The features of the proposed formulation are: (i) the local solution is obtained through an iterative procedure. The plastic corrections are given in closed forms exhibiting one scalar function denoted by $G_{\text{alg}}$ and three fourth-order tensors $\mathcal{D}_{\text{alg}}$, $\mathcal{G}_{\text{alg}}$, $\mathcal{L}_{\text{alg}}$, which are shown to be the algorithmic discrete counterparts of usual theoretical continuum quantities, (ii) the consistent tangent modulus has a symmetrical expression involving the same quantities. Finite element computations are performed using a particular non-linear kinematic hardening model and allow to exhibit the ratcheting phenomenon usually observed on mechanical components subjected to cyclic loadings.

**Keywords:** Elastoplasticity; Local integration; Plastic correction; Consistent tangent modulus; Ratcheting effect

## 1. Introduction

The numerical response of an elastoplastic structure is usually determined by an incremental method where the solution is computed at every time step by using the finite element method combined with a Newton-type iterative scheme. Within each iteration it is necessary to perform the two following important tasks characteristic of any elastoplastic algorithm: first, the local integration of the non-linear constitutive equations in order to compute the stress for a given strain increment and to build up the internal force vector, second, the computation of the consistent tangent modulus in order to form the structural tangent matrix. Since these tasks represent the main part of the elastoplastic calculation, a lot of work were devoted to their analysis.

In some cases an analytical local integration is possible. Thus, Krieg and Krieg [16] gave the exact local solution in perfect plasticity with von Mises yield function, by assuming a given constant and deviatoric increment of strain rate. Later, dealing with the plane stress situation, Simo and Govindjee [29] showed that the linear kinematic hardening leads to an algebraic equation of degree four which can be solved analytically. Recently, Alfano et al. [1] successfully generalized the treatment in [29] to include the significant case of isotropic hardening as well as elastoviscoplasticity of the Perzyna type, and also reduced the solution procedure to that of a quartic equation in the plastic multiplier.

In view of numerical purposes when analytical solutions are not available, the first local integration algorithm using Euler implicit scheme—referred to as radial return algorithm—was proposed by Wilkins [38] for perfect plasticity with von Mises yield criterion. Later, it was extended to linear isotropic and kinematic hardenings by Krieg and Key [15]. The further case of non-linear isotropic hardening combined with a special form

<footnote><sup>a</sup> Laboratoire de Génie Civil, Faculté des Sciences et des Techniques, 2 rue de la Houssinière, 44322 Nantes Cedex 3, France</footnote>

<footnote><sup>b</sup> Laboratoire de Mécanique de Lille, URA CNRS 1441, Université de Lille I, 59655 Villeneuve d'Ascq Cedex, France</footnote>
of non-linear kinematic hardening was investigated by Simo and Taylor [34]. In the slightly different frame of non-smooth mechanics, the catching up algorithm was proposed by Moreau [19,20] and used later on by Nguyen [22].

In order to handle other criteria than von Mises one, Ortiz et al. [24] proposed the closest point algorithm as a generalization of the radial return algorithm to arbitrary convex yield functions. The difficult case of plane stress was solved by Simo and Taylor [35]. Besides, other formulations using either the generalized trapezoidal or the generalized midpoint rule were presented by Ortiz and Popov [25], Ortiz and Martin [23], and Simo and Govindjee [30], with a view to encompassing the previous algorithms and combining numerical stability and accuracy. The cutting-plane algorithm successively proposed by Simo and Ortiz [33] in large deformations and Ortiz and Simo [26] in small strains is another stress update procedure lying on a different conception and may entail less computational efforts. After Hughes and Pister [12] introduced the notion of consistent linearization, the role of consistent tangent modulus was emphasized by Hughes and Taylor [13]. Later, a general local integration algorithm and the corresponding consistent tangent modulus for rate-independent plasticity were presented by Simo and Hughes [31,32], Ortiz and Martin [23], and Simo and Govindjee [30]. More recently, Hartmann and Haupt [11] considered a non-linear kinematic hardening model due to Armstrong and Frederick [3] and used the implicit Euler time integration procedure to solve the non-linear equation system for the stress computation. By the same backward scheme, Matzenmiller and Taylor [18] gave the stress update for general isotropic yield criteria containing all the invariants of the stress tensor. Completely explicit expressions for the plastic corrections and the consistent tangent modulus were given by Doghri [6–8] for $J_2$-plasticity possibly combined with damage. A first general approach to the derivation of the explicit expression for the consistent tangent modulus was provided by Alfano and Rosati [2] in rate-independent elastoplasticity with mixed hardening. More recently, Palazzo et al. [27] succeeded in obtaining explicit representations for the stress update and the consistent tangent operator for general isotropic elastoplastic models, by means of a novel linear combination of dyadic and square tensor products. Finally, it should be mentioned that one of the first solution for general inelastic analysis is the effective-stress-function algorithm proposed by Bathe and coworkers [4,14]. The complete stress state is obtained by solving a governing scalar equation called the effective stress function equation and the consistent tangent matrix is derived thereof. This method is actually applicable for more complex inelastic situations, like Drucker–Prager soil model, thermoelastoplasticity and creep solutions, as well as large strain elastoplasticity [9].

The above review makes no claim to completeness, the reader can find more details on the historical background to the local integration subject as well as more references in the quoted papers.

In this paper, we attempt to go further by proposing a more general, synthetic method to derive the plastic corrections and the consistent tangent modulus for a wide class of non-linear hardening models. The formulation is built for small strain elastoplasticity within the framework of the generalized standard materials theory described in [10,28]. The assumptions made are general enough for the proposed formulation to be valid with somewhat arbitrary yield functions and hardening energies. The local solution of the elastoplastic equations involves the one-step fully implicit integration scheme combined with a local Newton iterative procedure, as usually done in the literature. However, compared with Ref. [31] for instance, here the expressions for the plastic corrections and the consistent tangent modulus are obtained in closed (tensorial or matricial) forms ready for computer use, rather than in operator forms. The presented formulation highlights the role of one scalar function and three fourth-order tensors, which are shown to be the algorithmic discrete counterparts of usual theoretical continuum quantities. Also, consideration of general expressions for the hardening energies gives rise to some additional terms in the plastic corrections, which are not reported in the previous works.

The finite element computations carried out on the example of a particular non-linear kinematic hardening model allow to assess the validity of the proposed algorithm and to exhibit different types of ratcheting usually observed on mechanical components subjected to cyclic loadings.

## 2. Local integration

Let us consider a three-dimensional continuum body undergoing a quasi-static infinitesimal elastic–plastic transformation. The linearized strain tensor is decomposed into an elastic and plastic part, denoted by $\bar{\bar{\epsilon}}$ and $\bar{\bar{\epsilon}}^p$ respectively, according to the relationship

$$ \bar{\bar{\epsilon}} = \bar{\bar{\epsilon}}^c + \bar{\bar{\epsilon}}^p \quad (1) $$

The state and evolution equations for the elastoplastic problem are obtained according to the generalized standard materials theory which is based on thermodynamical foundations [10,28]. The following assumptions—not so restrictive—are made on the internal state of the body.

**Assumption 1.** The internal state at any point in the body and at any time is characterized by the plastic
strain $\bar{\mathbf{e}}^p$ and two other internal variables to account for
the deformation history, one second-order tensor harden-
ing variable denoted by $\bar{\boldsymbol{\alpha}}$ and one scalar hardening
variable denoted by $\alpha$.

**Assumption 2.** The hardening energy $w^{\text{hard}}$ splits up into two independent functions of the hardening variables

$$
w^{\text{hard}}(\bar{\boldsymbol{\alpha}}, \alpha) = w^{\boldsymbol{\beta}}(\bar{\boldsymbol{\alpha}}) + w^{\alpha}(\alpha) \quad (2)
$$

where $w^{\boldsymbol{\beta}}$ and $w^{\alpha}$ are strictly convex functions.

The state laws give us the stress $\bar{\sigma}$ and the thermo-
dynamic conjugates $\bar{a}$ and $a$, duals of the hardening
variables $\bar{\boldsymbol{\alpha}}$ and $\alpha$ respectively

$$
\begin{align}
\bar{\sigma} &= \mathcal{D} : \dot{\bar{\mathbf{e}}}^\mathbf{p} = \mathcal{D} : (\bar{\mathbf{e}} - \bar{\mathbf{e}}^\mathbf{p}) \\
\bar{a} &= \frac{\partial w^{\boldsymbol{\alpha}}}{\partial \bar{\boldsymbol{\alpha}}} (\bar{a}) \tag{3} \\
a &= \frac{\partial w^{\alpha}}{\partial \alpha}(\alpha)
\end{align}
$$

where $\mathcal{D}$ denotes the elasticity fourth-order tensor de-
fined as the second derivative of the elastic energy
$\frac{1}{2}\bar{\mathbf{e}}^\mathbf{p} : \mathcal{D} : \bar{\mathbf{e}}^\mathbf{p}$. The yield criterion expressed in terms of the
stress $\bar{\sigma}$ and the hardening variables $(\bar{a}, a)$ is satisfied at
any point in plastic loading

$$
f(\bar{\sigma}, \bar{a}, a) = 0 \quad (4)
$$

The evolution laws result from generalization of the
maximum plastic dissipation rule and provide the rates
of plastic strain and hardening variables by differenti-
ating the yield function $f$

$$
\begin{align}
\dot{\bar{\mathbf{e}}}^p &= \dot{\lambda} \frac{\partial f}{\partial \bar{\sigma}} \nonumber \\
- \dot{\bar{\boldsymbol{\alpha}}} &= \dot{\lambda} \frac{\partial f}{\partial \bar{a}} \tag{5} \\
- \dot{\alpha} &= \dot{\lambda} \frac{\partial f}{\partial a} \nonumber
\end{align}
$$

where $\dot{\lambda}$ is the plastic multiplier.

Given a strain increment $\Delta \bar{\mathbf{e}} = \bar{\mathbf{e}} - \bar{\mathbf{e}}_{n-1}$ at a point in
the body and at present time $t_n$, let us assume that there
is plastic loading and determine the elastoplastic solu-
tion in terms of $\Delta \bar{\mathbf{e}}$ and the other known quantities at
previous time $t_{n-1}$.

By differentiating the state laws (3) with respect to
time and eliminating the conjugate hardening variables
by means of the evolution laws (5), we get 14 equations
with 14 unknowns, viz. the stress $\bar{\sigma}$, the hardening
variables $(\bar{a}, a)$ and the plastic multiplier $\dot{\lambda}$

$$
\begin{align}
\dot{\bar{\boldsymbol{\alpha}}} &= \mathcal{D} : \dot{\bar{\mathbf{e}}} - \dot{\lambda} \mathcal{D} : \frac{\partial f}{\partial \bar{\boldsymbol{\sigma}}} (\bar{\boldsymbol{\sigma}}, \bar{\boldsymbol{\alpha}}, a) \notag \\
\dot{\bar{\boldsymbol{\alpha}}} &= \frac{\partial^2 w^{\boldsymbol{\alpha}}}{\partial \bar{\boldsymbol{\alpha}}^2} : \dot{\bar{\boldsymbol{\alpha}}} = - \dot{\lambda} \mathcal{G}(\bar{\boldsymbol{\alpha}}) : \frac{\partial f}{\partial \bar{\boldsymbol{\alpha}}} (\bar{\boldsymbol{\sigma}}, \bar{\boldsymbol{\alpha}}, a) \tag{6} \\
\dot{a} &= \frac{\partial^2 w^{\alpha}}{\partial \alpha^2} \dot{\alpha} = - \dot{\lambda} G(a) \frac{\partial f}{\partial a} (\bar{\boldsymbol{\sigma}}, \bar{\boldsymbol{\alpha}}, a) \notag \\
f(\bar{\boldsymbol{\sigma}}, \bar{\boldsymbol{\alpha}}, a) &= 0 \notag
\end{align}
$$

where $\mathcal{G}(\bar{\boldsymbol{\alpha}})$ is a symmetric tensor of order four and $G(a)$
a scalar function defined by

**Definition**

$$
\mathcal{G}(\bar{\boldsymbol{\alpha}}) \equiv \frac{\partial^2 w^{\boldsymbol{\alpha}}}{\partial \bar{\boldsymbol{\alpha}}^2}(\bar{\boldsymbol{\alpha}}(\bar{\boldsymbol{\alpha}))), \quad G(a) \equiv \frac{\partial^2 w^{\alpha}}{\partial \alpha^2}(\alpha(a)) \quad (7)
$$

In the sequel, $\mathcal{G}$ is referred to as the plastic (hard-
ening) tensor—as opposed to the elastic tensor $\mathcal{D}$—and
$G$ as the plastic (hardening) scalar function. In writing
(7), use has been made of the strict convexity of the
hardening energies $w^{\boldsymbol{\beta}}$ and $w^{\alpha}$ according to which the
state laws (3b) and (3c) can be inverted giving $(\bar{\boldsymbol{\alpha}}, \alpha)$ in
terms of $(\bar{\boldsymbol{\alpha}}, a)$.

Relations (6) are now discretized by means of the
implicit Euler scheme, giving rise to 14 scalar equations
with 14 unknowns $(\bar{\sigma}, \bar{a}, a, \Delta\lambda)$

$$
\begin{align}
\bar{\boldsymbol{\sigma}} - \bar{\boldsymbol{\sigma}}_{n-1} + \Delta\lambda \mathcal{D} : \frac{\partial f}{\partial \bar{\boldsymbol{\sigma}}} (\bar{\boldsymbol{\sigma}}, \bar{\boldsymbol{\alpha}}, a) - \mathcal{D} : \Delta\bar{\boldsymbol{\epsilon}} &= \bar{\boldsymbol{\emptyset}} \tag{8} \\
\bar{\boldsymbol{\alpha}} - \bar{\boldsymbol{\alpha}}_{n-1} + \Delta\lambda \mathcal{G}(\bar{\boldsymbol{\alpha}}) : \frac{\partial f}{\partial \bar{\boldsymbol{\alpha}}} (\bar{\boldsymbol{\sigma}}, \bar{\boldsymbol{\alpha}}, a) &= \bar{\boldsymbol{\emptyset}} \nonumber \\
a - a_{n-1} + \Delta\lambda G(a) \frac{\partial f}{\partial a} (\bar{\boldsymbol{\sigma}}, \bar{\boldsymbol{\alpha}}, a) &= 0 \nonumber \\
f(\bar{\boldsymbol{\sigma}}, \bar{\boldsymbol{\alpha}}, a) &= 0 \nonumber
\end{align}
$$

In the above equations, the subscripts for quantities at present time $t_n$ have been omitted for brevity, thus $\bar{\sigma}$ for instance designates $\bar{\sigma}_n$. The quantity $\Delta\lambda$ can be considered as an incremental plastic multiplier. The backward implicit Euler scheme ensures the numerical stability and the symmetry of the consistent tangent modulus, as has been shown in number of works. To solve the non-linear system (8) by means of Newton iterative method, we make the following hypothesis which is satisfied in practice.

**Assumption 3.** The yield function $f(\bar{\sigma}, \bar{a}, a)$ is decom-
posed as the sum of one term depending on $(\bar{\sigma}, \bar{a})$ and
another one depending on the scalar variable $a$

$$
f(\bar{\sigma}, \bar{a}, a) = f_1(\bar{\sigma}, \bar{a}) + f_2(a) \quad (9)
$$

The form (9) is not so restrictive in practice and is
general enough to enable one to deal with new plasticity
models. Note that variables $(\bar{\sigma}, \bar{a})$ do not necessarily
appear in the form of their difference $\bar{\sigma} - \bar{a}$, as we will see
below in the case of non-linear kinematic hardening.
Since the elastic and hardening energies are strictly convex, $\mathcal{D}$ and $\mathcal{G}(\bar{\boldsymbol{a}})$ are invertible. Taking into account (9) at each local iteration for solving (8), the plastic corrections $(\delta\bar{\boldsymbol{a}}, \delta\bar{\boldsymbol{a}}, \delta a, \delta\lambda)$ for $(\bar{\boldsymbol{a}}, \bar{\boldsymbol{a}}, a, \Delta\lambda)$ are given by

$$
\begin{align}
\left(\mathcal{D}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}}^2}\right) : \delta\bar{\boldsymbol{a}} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}} \partial \bar{\boldsymbol{a}}} : \delta\bar{\boldsymbol{a}} &\nonumber \\
&= -\overline{\Psi}(\bar{\boldsymbol{a}}, \bar{\boldsymbol{a}}, \Delta\lambda) - \frac{\partial f}{\partial \bar{\boldsymbol{a}}} \delta\lambda \nonumber \\
\Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}} \partial \bar{\boldsymbol{a}}} : \delta\bar{\boldsymbol{a}} &+ \left[ \mathcal{G}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}}^2} \right. \nonumber \\
&\qquad \left. + \Delta\lambda \mathcal{G}^{-1} : \left( \frac{\partial f}{\partial \bar{\boldsymbol{a}}} : \frac{\partial^3 w^\lambda}{\partial \bar{\boldsymbol{a}}^3} \right) : \mathcal{G}^{-1} \right] : \delta\bar{\boldsymbol{a}} \tag{10} \\
&= -\overline{\Phi}(\bar{\boldsymbol{a}}, \bar{\boldsymbol{a}}, \Delta\lambda) - \frac{\partial f}{\partial \bar{\boldsymbol{a}}} \delta\lambda \nonumber \\
\left(\frac{1}{G(a)} + \Delta\lambda \frac{\partial^2 f}{\partial a^2} + \frac{\Delta\lambda}{G^2(a)} \frac{\partial^3 w^\lambda}{\partial \alpha^3} \frac{\partial f}{\partial a} \right) \delta a &\nonumber \\
&= -\Phi(a, \Delta\lambda) - \frac{\partial f}{\partial a} \delta\lambda \nonumber \\
\frac{\partial f}{\partial \bar{\boldsymbol{a}}} : \delta\bar{\boldsymbol{a}} + \frac{\partial f}{\partial \bar{\boldsymbol{a}}} : \delta\bar{\boldsymbol{a}} + \frac{\partial f}{\partial a} \delta a &= -f(\bar{\boldsymbol{a}}, \bar{\boldsymbol{a}}, a) \nonumber
\end{align}
$$

where $\overline{\Psi}$ and $\overline{\Phi}$ are second-order tensors and $\Phi$ a scalar function defined by

**Definition**

$$
\begin{align}
\overline{\Psi}(\bar{\boldsymbol{a}}, \bar{\boldsymbol{a}}, \Delta\lambda) &= \mathcal{D}^{-1} : (\bar{\boldsymbol{a}} - \bar{\boldsymbol{a}}_{n-1}) + \Delta\lambda \frac{\partial f}{\partial \bar{\boldsymbol{a}}} : (\bar{\boldsymbol{a}}, \bar{\boldsymbol{a}}) - \Delta\bar{\boldsymbol{e}} \nonumber \\
\overline{\Phi}(\bar{\boldsymbol{a}}, \bar{\boldsymbol{a}}, \Delta\lambda) &= \mathcal{G}^{-1}(\bar{\boldsymbol{a}}) : (\bar{\boldsymbol{a}} - \bar{\boldsymbol{a}}_{n-1}) + \Delta\lambda \frac{\partial f}{\partial \bar{\boldsymbol{a}}} : (\bar{\boldsymbol{a}}, \bar{\boldsymbol{a}}) \tag{11} \\
\Phi(a, \Delta\lambda) &= \frac{a - a_{n-1}}{G(a)} + \Delta\lambda \frac{\partial f}{\partial a}(a) \nonumber
\end{align}
$$

Note that the inner parentheses in the left-hand side of Eq. (10b) are essential since the colon products therein are not associative. Moreover, since the components of the fourth-order tensor $\frac{\partial f}{\partial \bar{\boldsymbol{a}}} : \frac{\partial^3 w^\lambda}{\partial \bar{\boldsymbol{a}}^3}$ in a fixed orthonormal basis read (with summation implied over repeated subscripts)

$$
\begin{align}
\left(\frac{\partial f}{\partial \bar{\boldsymbol{a}}} : \frac{\partial^3 w^\lambda}{\partial \bar{\boldsymbol{a}}^3}\right)_{ijkl} &= \left(\frac{\partial f}{\partial \bar{\boldsymbol{a}}}\right)_{mn} \left(\frac{\partial^3 w^\lambda}{\partial \bar{\boldsymbol{a}}^3}\right)_{nmijkl} \nonumber \\
&= \left(\frac{\partial f}{\partial \bar{\boldsymbol{a}}}\right)_{mn} \frac{\partial^3 w^\lambda}{\partial \bar{\boldsymbol{a}}_m \partial \bar{\boldsymbol{a}}_n \partial \bar{\boldsymbol{a}}_l \partial \bar{\boldsymbol{a}}_k} \nonumber \\
&= \left(\frac{\partial f}{\partial \bar{\boldsymbol{a}}}\right)_{klij} \tag{12}
\end{align}
$$

the fourth-order tensor $\frac{\partial f}{\partial \bar{\boldsymbol{a}}} : \frac{\partial^3 w^\lambda}{\partial \bar{\boldsymbol{a}}^3}$ is symmetric. Therefore, Eq. (10a) and (10b) form a *symmetric* system with unknowns $\delta\bar{\boldsymbol{a}}$ and $\delta\bar{\boldsymbol{a}}$. These equations, *uncoupled* from (10c) with unknown $\delta a$, take the following form

$$
\begin{align}
\tilde{\mathcal{A}} : \delta\bar{\boldsymbol{a}} + \tilde{\mathcal{B}} : \delta\bar{\boldsymbol{a}} &= \bar{\boldsymbol{p}} \nonumber \\
\tilde{\mathcal{B}}^T : \delta\bar{\boldsymbol{a}} + \tilde{\mathcal{C}} : \delta\bar{\boldsymbol{a}} &= \bar{\boldsymbol{q}} \tag{13}
\end{align}
$$

In practice, the incremental plastic multiplier $\Delta\lambda$ is small enough like the strain increment. Thus, the fourth-order tensors

$$
\tilde{\mathcal{C}} \equiv \mathcal{C}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}}^2} + \Delta\lambda \mathcal{G}^{-1} : \left( \frac{\partial f}{\partial \bar{\boldsymbol{a}}} : \frac{\partial^3 w^\lambda}{\partial \bar{\boldsymbol{a}}^3} \right) : \mathcal{G}^{-1}
$$

and

$$
\tilde{\mathcal{A}} \equiv \mathcal{D}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}}^2}
$$

are invertible like $\mathcal{G}^{-1}$ and $\mathcal{D}^{-1}$. In such conditions, the system (13) is solved as follows

$$
\begin{align}
\delta\bar{\boldsymbol{a}} &= \tilde{\mathcal{D}} : \bar{\boldsymbol{p}} + \tilde{\mathcal{L}} : \bar{\boldsymbol{q}} \nonumber \\
\delta\bar{\boldsymbol{a}} &= \tilde{\mathcal{L}}^T : \bar{\boldsymbol{p}} + \tilde{\mathcal{G}} : \bar{\boldsymbol{q}} \tag{14}
\end{align}
$$

$$
\begin{align*}
\text{with } \tilde{\mathcal{D}} &= (\tilde{\mathcal{A}} - \tilde{\mathcal{B}} : \tilde{\mathcal{C}}^{-1} : \tilde{\mathcal{B}}^T)^{-1}, \quad \tilde{\mathcal{L}} = -\tilde{\mathcal{A}}^{-1} : \tilde{\mathcal{B}} : \tilde{\mathcal{G}} \\
\text{and } \tilde{\mathcal{G}} &= (\tilde{\mathcal{C}} - \tilde{\mathcal{B}}^T : \tilde{\mathcal{A}}^{-1} : \tilde{\mathcal{B}})^{-1}.
\end{align*}
$$

From (14), the incremental plastic multiplier $\Delta\lambda$, the stress $\bar{\boldsymbol{a}}$ and the hardening variables $(\bar{\boldsymbol{a}}, a)$ are updated at each local iteration using the following plastic corrections [36,37]

$$
\begin{align}
\delta\bar{\boldsymbol{a}} &= -\mathcal{D}_{\text{alg}} : \overline{\Psi} - \mathcal{L}_{\text{alg}} : \overline{\Phi} \nonumber \\
&\quad - \left( \mathcal{D}_{\text{alg}} : \frac{\partial f}{\partial \bar{\boldsymbol{a}}} + \mathcal{L}_{\text{alg}} : \frac{\partial f}{\partial \bar{\boldsymbol{a}}} \right) \delta\lambda \nonumber \\
\delta\bar{\boldsymbol{a}} &= -\mathcal{L}_{\text{alg}}^{\text{T}} : \overline{\Psi} - \mathcal{G}_{\text{alg}} : \overline{\Phi} \nonumber \\
&\quad - \left( \mathcal{L}_{\text{alg}}^{\text{T}} : \frac{\partial f}{\partial \bar{\boldsymbol{a}}} + \mathcal{G}_{\text{alg}} : \frac{\partial f}{\partial \bar{\boldsymbol{a}}} \right) \delta\lambda \tag{15} \\
\delta a &= -G_{\text{alg}} \left( \Phi + \frac{\partial f}{\partial a} \delta\lambda \right) \nonumber
\end{align}
$$

where the scalar function $G_{\text{alg}}$ and the fourth-order tensors $\mathcal{D}_{\text{alg}}$, $\mathcal{G}_{\text{alg}}$ and $\mathcal{L}_{\text{alg}}$ are defined by

**Definition**

$$
\begin{align}
\mathcal{D}_{\text{alg}} &= \left[ \mathcal{D}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}}^2} - (\Delta\lambda)^2 \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}} \partial \bar{\boldsymbol{a}}} : \left( \mathcal{G}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}}^2} \right. \right. \nonumber \\
&\qquad \left. \left. + \Delta\lambda \mathcal{G}^{-1} : \left( \frac{\partial f}{\partial \bar{\boldsymbol{a}}} : \frac{\partial^3 w^\lambda}{\partial \bar{\boldsymbol{a}}^3} \right) : \mathcal{G}^{-1} \right)^{-1} : \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}} \partial \bar{\boldsymbol{a}}} \right]^{-1} \nonumber \\
\mathcal{G}_{\text{alg}} &= \left[ \mathcal{G}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}}^2} + \Delta\lambda \mathcal{G}^{-1} : \left( \frac{\partial f}{\partial \bar{\boldsymbol{a}}} : \frac{\partial^3 w^\lambda}{\partial \bar{\boldsymbol{a}}^3} \right) : \mathcal{G}^{-1} \right. \nonumber \\
&\qquad \left. - (\Delta\lambda)^2 \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}} \partial \bar{\boldsymbol{a}}} : \left( \mathcal{D}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}}^2} \right)^{-1} : \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}} \partial \bar{\boldsymbol{a}}} \right]^{-1} \nonumber
\end{align}
$$

$$
\begin{align}
G_{\text{alg}} &= \frac{G(a)}{1 + \Delta\lambda G(a) \frac{\partial^2 f}{\partial a^2} + \frac{\Delta\lambda}{G(a)} \frac{\partial^3 w^\lambda}{\partial a^3} \frac{\partial f}{\partial a}} \nonumber \\
\mathcal{L}_{\text{alg}} &= - \left( \mathcal{D}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}}^2} \right)^{-1} : \Delta\lambda \frac{\partial^2 f}{\partial \bar{\boldsymbol{a}} \partial \bar{\boldsymbol{a}}} : \mathcal{G}_{\text{alg}} \tag{16}
\end{align}
$$

One can check that the transpose of $\mathcal{L}_{\text{alg}}$ can be put in the form
$$
\begin{align*}
\mathcal{L}_{\text{alg}}^{\text{T}} &= -\tilde{\mathcal{C}}^{-1} : \tilde{\mathcal{B}}^{\text{T}} : \tilde{\mathcal{D}} \\
&= -\left[ \mathcal{G}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\alpha}^2} + \Delta\lambda \mathcal{G}^{-1} : \left( \frac{\partial f}{\partial \bar{\alpha}} : \frac{\partial^3 w^{\bar{\alpha}}}{\partial \bar{\alpha}^3} \right) : \mathcal{G}^{-1} \right]^{-1} \\
&\quad : \Delta\lambda \frac{\partial^2 f}{\partial \bar{\alpha} \partial \bar{\alpha}} : \mathcal{G}_{\text{alg}}
\end{align*}
$$

Eventually, inserting relations (15) into (10d) provides
the correction for the plastic multiplier

$$
\delta\lambda = \left[ f - \left( \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{G}_{\text{alg}} + \frac{\partial f}{\partial \bar{\alpha}} : \mathcal{L}_{\text{alg}}^{\text{T}} \right) : \bar{\Psi} \right. \\
\left. - \left( \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{L}_{\text{alg}} + \frac{\partial f}{\partial \bar{\alpha}} : \mathcal{G}_{\text{alg}} \right) : \bar{\Phi} - G_{\text{alg}} \frac{\partial f}{\partial \alpha} \Phi \right] \\
\left/ \left[ \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{G}_{\text{alg}} : \frac{\partial f}{\partial \bar{\sigma}} + 2 \frac{\partial f}{\partial \bar{\alpha}} : \mathcal{L}_{\text{alg}} : \frac{\partial f}{\partial \bar{\alpha}} \right. \right. \\
\left. \left. + \frac{\partial f}{\partial \bar{\alpha}} : \mathcal{G}_{\text{alg}} : \frac{\partial f}{\partial \bar{\alpha}} + G_{\text{alg}} \left( \frac{\partial f}{\partial \alpha} \right)^2 \right] \right\} \quad (17)
$$

Recall that the right-hand sides of (15)–(17) are com-
puted with values obtained at the previous local itera-
tion. Since the fourth-order tensor $\frac{\partial f}{\partial \bar{\alpha}} : \frac{\partial^3 w^{\bar{\alpha}}}{\partial \bar{\alpha}^3}$ is symmetric,
$\mathcal{D}_{\text{alg}}$ and $\mathcal{G}_{\text{alg}}$ are symmetric too. To sum up, the local
solution of the elastoplastic problem is described by the
following algorithm.

Local integration loop

(a) Compute

$$
\left[ \mathcal{G}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\alpha}^2} + \Delta\lambda \mathcal{G}^{-1} : \left( \frac{\partial f}{\partial \bar{\alpha}} : \frac{\partial^3 w^{\bar{\alpha}}}{\partial \bar{\alpha}^3} \right) : \mathcal{G}^{-1} \right]^{-1}
$$

and derive $\mathcal{D}_{\text{alg}}$ from (16a).

(b) Compute

$$
\left( \mathcal{D}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\alpha}^2} \right)^{-1}
$$

and derive $\mathcal{G}_{\text{alg}}$ from (16b) and $\mathcal{L}_{\text{alg}}$ from (16d).

(c) Compute $G_{\text{alg}}$ from (16c).

(d) Compute the plastic corrections $\delta\lambda$ by (17), $(\delta\bar{\sigma}, \delta\bar{\alpha}, \delta\alpha)$ by (15), and update the solution $(\Delta\lambda, \bar{\sigma}, \bar{\alpha}, \alpha)$.

(e) Check the convergence, based on $\delta\lambda$ for instance, to exit the loop.

End of the loop

The proposed local integration procedure combines the
quadratic convergence rate of Newton scheme with the
unconditional stability of implicit Euler scheme. Expressions (15)–(17) for the plastic corrections show the major
role of tensors $\mathcal{D}_{\text{alg}}$, $\mathcal{G}_{\text{alg}}$, $\mathcal{L}_{\text{alg}}$ and scalar function $G_{\text{alg}}$ in
the inversion of the local equation set (8). The quantities $\mathcal{D}_{\text{alg}}$, $\mathcal{G}_{\text{alg}}$ and $G_{\text{alg}}$ represent the algorithmic or discrete
finite-step counterparts of the usual theoretical continuum
quantities: the elastic tensor $\mathcal{D}$, the plastic tensor $\mathcal{G}$ and the

plastic scalar function $G$, see definition (7), respectively. As
shown by definition (16), when the incremental plastic
multiplier $\Delta\lambda$ tends to zero, all the algorithmic quantities
tend to their theoretical counterparts.

In this study, we assume that the hardening energies $w^\alpha(\bar{\alpha})$ and $w^\alpha(\alpha)$ in Eq. (2) can be given general expressions. Thus, the plastic (hardening) tensor $\mathcal{G}$ and the plastic (hardening) scalar function $G$, defined by relation (7), are not constant a priori. In particular, this leads to the term $\frac{\partial f}{\partial \bar{\alpha}} : \frac{\partial^3 w^{\bar{\alpha}}}{\partial \bar{\alpha}^3}$ in Eq. (16a) and (16b), which is not reported in the literature, see e.g. Box 3.5 in [31]. Although to date usual models involve quadratic expressions for $w^\alpha$ so that $\frac{\partial f}{\partial \bar{\alpha}} : \frac{\partial^3 w^{\bar{\alpha}}}{\partial \bar{\alpha}^3}$ is zero, theoretical considerations give rise to the latter term which must be taken into account in general.

In the case of arbitrary yield functions and hardening
energies, fourth-order tensors $\mathcal{D}_{\text{alg}}$, $\mathcal{G}_{\text{alg}}$ and $\mathcal{L}_{\text{alg}}$ must
be computed numerically from definition (16), and this
requires inversions of 6×6 symmetric matrices.

**Remark.** In some cases, these numerical computations
can be avoided as recently shown by Palazzo et al. [27].
These authors, by using an original linear combination
of dyadic and square tensor products, have successfully
derived explicit representations for the stress update and
the consistent tangent operator, for general isotropic
yield functions depending upon all the three invariants
of the relative stress $\bar{\sigma} - \bar{\alpha}$.

Furthermore, using the symmetry of $\mathcal{D}_{\text{alg}}$ and $\mathcal{G}_{\text{alg}}$,
relations (15a), (15b) and (17) can be recast as follows

$$
\delta\lambda = \left[ f - \left( \frac{\partial f}{\partial \bar{\sigma}} \right)_{\text{alg}} : \mathcal{D}_{\text{alg}} : \bar{\Psi} - \left( \frac{\partial f}{\partial \bar{\alpha}} \right)_{\text{alg}} : \mathcal{G}_{\text{alg}} : \bar{\Phi} \right. \\
\left. - G_{\text{alg}} \frac{\partial f}{\partial \alpha} \Phi \right] \bigg/ \left[ \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{D}_{\text{alg}} : \left( \frac{\partial f}{\partial \bar{\sigma}} \right)_{\text{alg}} \right. \\
\left. + \frac{\partial f}{\partial \bar{\alpha}} : \mathcal{G}_{\text{alg}} : \left( \frac{\partial f}{\partial \bar{\alpha}} \right)_{\text{alg}} + G_{\text{alg}} \left( \frac{\partial f}{\partial \alpha} \right)^2 \right] \quad (18)
$$

and

$$
\begin{equation}
\begin{split}
\delta \bar{\sigma} = \mathcal{D}_{\text{alg}} : \left[ \right. & - \bar{\Psi} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\alpha} \partial \bar{\sigma}} : (\mathcal{G}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\alpha}^2} \\
& \left. + \Delta\lambda \mathcal{G}^{-1} : \left( \frac{\partial f}{\partial \bar{\alpha}} : \frac{\partial^3 w^{\bar{\alpha}}}{\partial \bar{\alpha}^3} \right) : \mathcal{G}^{-1} \right]^{-1} : \bar{\Phi} - \left( \frac{\partial f}{\partial \bar{\sigma}} \right)_{\text{alg}} \delta\lambda \\
& \left. \left. \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{G}_{\text{alg}} : \left[ \right. \right. \right. \\
& \left. \left. \left. - \bar{\Phi} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\alpha} \partial \bar{\sigma}} : (\mathcal{D}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\alpha}^2}) \right. \right. \\
& \left. \left. \left. : \bar{\Psi} - \left( \frac{\partial f}{\partial \bar{\alpha}} \right)_{\text{alg}} \delta\lambda \right] \right] \right\}
\end{split}
\tag{19}
\end{equation}
$$

where $(\frac{\partial f}{\partial \bar{\sigma}})_{\text{alg}}$ and $(\frac{\partial f}{\partial \bar{\alpha}})_{\text{alg}}$ are second-order tensors defined by
**Definition**

$$
\begin{align}
\left(\frac{\partial f}{\partial \bar{\sigma}}\right)_{\text{alg}} &= \frac{\partial f}{\partial \bar{\sigma}} - \Delta\lambda \frac{\partial f}{\partial \bar{a}} : \left[ \mathcal{G}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{a}^2} \right. \nonumber \\
&\qquad \left. + \Delta\lambda \mathcal{G}^{-1} : \left( \frac{\partial f}{\partial \bar{a}} : \frac{\partial^3 w^{\bar{a}}}{\partial \bar{a}^3} \right) : \mathcal{G}^{-1} \right]^{-1} : \frac{\partial^2 f}{\partial \bar{\sigma} \partial \bar{a}} \nonumber \\
\left(\frac{\partial f}{\partial \bar{a}}\right)_{\text{alg}} &= \frac{\partial f}{\partial \bar{a}} - \Delta\lambda \frac{\partial f}{\partial \bar{\sigma}} : \left( \mathcal{D}^{-1} + \Delta\lambda \frac{\partial^2 f}{\partial \bar{\sigma}^2} \right)^{-1} : \frac{\partial^2 f}{\partial \bar{a} \partial \bar{\sigma}} \tag{20}
\end{align}
$$

Expression (18) for $\delta\lambda$, strictly equivalent to (17), is quite
analogous to the theoretical expression for $\hat{\lambda}$ derived
from the continuous problem:

$$
\dot{\lambda} = \frac{\frac{\partial f}{\partial \bar{\sigma}} : \mathcal{D} : \dot{\bar{\varepsilon}}}{\frac{\partial f}{\partial \bar{\sigma}} : \mathcal{D} : \frac{\partial f}{\partial \bar{a}} + \frac{\partial f}{\partial \bar{a}} : \mathcal{G}(\bar{a}) : \dot{\bar{\varepsilon}} + G(a) \left(\frac{\partial f}{\partial a}\right)^2} \quad (21)
$$

where the sum appearing in the denominator

$$
\bar{h} \equiv \frac{\partial f}{\partial \bar{a}} : \mathcal{G}(\bar{a}) : \frac{\partial f}{\partial \bar{a}} + G(a) \left(\frac{\partial f}{\partial a}\right)^2
$$

is usually referred to as the plastic modulus. When the
incremental plastic multiplier $\Delta\lambda$ tends to zero, tensors
$(\frac{\partial f}{\partial \bar{\sigma}})_{\text{alg}}$ and $(\frac{\partial f}{\partial \bar{a}})_{\text{alg}}$ tend to their theoretical continuum
counterparts $\frac{\partial f}{\partial \bar{\sigma}}$ and $\frac{\partial f}{\partial \bar{a}}$, respectively.

3. Consistent tangent modulus

The local integration is now assumed to be achieved
so that the local solution $(\bar{\sigma}, \bar{a}, a, \Delta\lambda)$ at present time $t_n$ of
the system (3)–(5) is available and a formal implicit re-
lation has been established between stress and strain,
$\bar{\sigma} = \dot{\bar{\sigma}}(\bar{e}, \bar{a}, \text{state at } t_{n-1})$, where the expression “state at $t_{n-1}$”
designates the known values of all the variables at pre-
vious time $t_{n-1}$. The concepts developed in the previous
section will allow us to readily derive the consistent
tangent modulus $(\frac{\partial \bar{\sigma}}{\partial \bar{e}})_{\text{alg}} = \frac{\partial \bar{\sigma}}{\partial \bar{e}}$ as the algorithmic coun-
terpart of the theoretical tangent modulus [13,21,34].

Differentiating relations (8) with respect to time yields

$$
\begin{align}
\dot{\bar{\sigma}} &= \mathcal{D} : \dot{\bar{\varepsilon}} - \dot{\Delta}\lambda \mathcal{D} : \frac{\partial f}{\partial \bar{\sigma}} - \Delta\lambda \mathcal{D} : \left( \frac{\partial^2 f}{\partial \bar{\sigma}^2} : \dot{\bar{\varepsilon}} + \frac{\partial^2 f}{\partial \bar{a}\bar{a}} : \dot{\bar{\varepsilon}} \right) \nonumber \\
\dot{\bar{a}} &= -\dot{\Delta}\lambda \mathcal{G} : \frac{\partial f}{\partial \bar{a}} - \Delta\lambda \left( \mathcal{G} : \frac{\partial^2 f}{\partial \bar{\sigma}\bar{a}} : \dot{\bar{\varepsilon}} \right. \nonumber \\
&\qquad \left. + \mathcal{G} : \frac{\partial^2 f}{\partial \bar{a}^2} : \dot{\bar{a}} + \frac{\partial f}{\partial \bar{a}} : \frac{\partial \mathcal{G}}{\partial \bar{a}} : \dot{\bar{\varepsilon}} \right) \nonumber \\
\dot{a} &= -\dot{\Delta}\lambda G \frac{\partial f}{\partial a} - \Delta\lambda \left( G \frac{\partial^2 f}{\partial a^2} + \frac{\partial f}{\partial a} \frac{\partial G}{\partial a} \right) \dot{a} \nonumber \\
\dot{f} &= \frac{\partial f}{\partial \bar{\sigma}} : \dot{\bar{\varepsilon}} + \frac{\partial f}{\partial \bar{a}} : \dot{\bar{a}} + \frac{\partial f}{\partial a} \dot{a} = 0 \tag{22}
\end{align}
$$

In the above relations, computations are made with the
values (Δλ, σ, a, a) deduced from the local integration.

By comparison with the theoretical tangent tensor $\frac{\partial \sigma}{\partial \bar{\varepsilon}}$
derived from system (6) for the continuous problem, the
algorithmic or consistent tangent tensor $(\frac{\partial \sigma}{\partial \bar{\varepsilon}})_{\text{alg}}$ must be
different because of additional terms in (22).

By noting that system (22a) and (22b) is analogous to
system (10a) and (10b) (replace tensor $\bar{\Psi}$ in (10a) by $-\dot{\bar{\varepsilon}}$
and make $\bar{\Phi} = \bar{\mathbf{0}}$ in (10b)), it can be inverted using (14)
as follows

$$
\begin{align}
\dot{\bar{\sigma}} &= \mathcal{D}_{\text{alg}} : \dot{\bar{\varepsilon}} - \left( \mathcal{D}_{\text{alg}} : \frac{\partial f}{\partial \bar{\sigma}} + \mathcal{L}_{\text{alg}} : \frac{\partial f}{\partial \bar{a}} \right) \dot{\Delta}\lambda \nonumber \\
\dot{\bar{a}} &= \mathcal{L}_{\text{alg}}^{\text{T}} : \dot{\bar{\varepsilon}} - \left( \mathcal{L}_{\text{alg}}^{\text{T}} : \frac{\partial f}{\partial \bar{\sigma}} + \mathcal{G}_{\text{alg}} : \frac{\partial f}{\partial \bar{a}} \right) \dot{\Delta}\lambda \tag{23}
\end{align}
$$

where tensors $\mathcal{D}_{\text{alg}}$, $\mathcal{G}_{\text{alg}}$ and $\mathcal{L}_{\text{alg}}$ are defined in (16).
Inserting relations (22c) and (23) into (22d) leads to

$$
\begin{equation}
\begin{split}
\dot{\Delta}\lambda = \left[ \left( \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{D}_{\text{alg}} + \frac{\partial f}{\partial \bar{a}} : \mathcal{L}_{\text{alg}}^{\text{T}} \right) / \left( \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{D}_{\text{alg}} : \frac{\partial f}{\partial \bar{\sigma}} \right) \right. \\
\left. \qquad + 2 \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{L}_{\text{alg}} : \frac{\partial f}{\partial \bar{a}} + \frac{\partial f}{\partial \bar{a}} : \mathcal{G}_{\text{alg}} : \frac{\partial f}{\partial \bar{a}} \right. \\
\left. \qquad \left. + G_{\text{alg}} \left( \frac{\partial f}{\partial a} \right)^2 \right) \right] : \dot{\bar{\varepsilon}}
\end{split}
\tag{24}
\end{equation}
$$

Eventually, by putting (24) in (23a) we obtain the con-
sistent tangent modulus

$$
\begin{equation}
\begin{split}
\left(\frac{\partial \bar{\sigma}}{\partial \bar{\varepsilon}}\right)_{\text{alg}} = \mathcal{D}_{\text{alg}} - & \left[ \left( \mathcal{D}_{\text{alg}} : \frac{\partial f}{\partial \bar{\sigma}} + \mathcal{L}_{\text{alg}} : \frac{\partial f}{\partial \bar{a}} \right) \otimes \left( \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{D}_{\text{alg}} + \frac{\partial f}{\partial \bar{a}} : \mathcal{L}_{\text{alg}}^{\text{T}} \right) \right] \\
& / \left[ \left( \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{D}_{\text{alg}} : \frac{\partial f}{\partial \bar{a}} + 2 \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{L}_{\text{alg}} : \frac{\partial f}{\partial \bar{a}} \right. \right. \\
& \left. \left. + \frac{\partial f}{\partial \bar{a}} : \mathcal{G}_{\text{alg}} : \frac{\partial f}{\partial \bar{a}} + G_{\text{alg}} \left( \frac{\partial f}{\partial a} \right)^2 \right) \right]
\end{split}
\tag{25}
\end{equation}
$$

where symbol $\otimes$ denotes dyadic tensor product defined
by $(\bar{e} \otimes \bar{a})_{ijkl} = c_{ij} d_{kl}$, $\forall$ second-order tensors $\bar{e}$ and $\bar{a}$.
Relation (25) can be recast into a simpler expression
using algorithmic tensors $(\frac{\partial f}{\partial \bar{\sigma}})_{\text{alg}}$ and $(\frac{\partial f}{\partial \bar{a}})_{\text{alg}}$ introduced
in (20)

$$
\left(\frac{\partial \bar{\sigma}}{\partial \bar{\varepsilon}}\right)_{\text{alg}} = \mathcal{D}_{\text{alg}} - \left[ \mathcal{D}_{\text{alg}} : \left(\frac{\partial f}{\partial \bar{\sigma}}\right)_{\text{alg}} \otimes \left(\frac{\partial f}{\partial \bar{\sigma}}\right)_{\text{alg}} : \mathcal{D}_{\text{alg}} \right] \\
\left/ \left[ \left[ \frac{\partial f}{\partial \bar{\sigma}} : \mathcal{D}_{\text{alg}} : \left(\frac{\partial f}{\partial \bar{\sigma}}\right)_{\text{alg}} + \frac{\partial f}{\partial \bar{a}} : \mathcal{G}_{\text{alg}} : \left(\frac{\partial f}{\partial \bar{a}}\right)_{\text{alg}} \right. \right. \right. \\
\left. \left. + G_{\text{alg}} \left(\frac{\partial f}{\partial a}\right)^2 \right] \right\} \tag{26}
$$

The consistent tangent modulus is symmetric since $\mathcal{D}_{\text{alg}}$
is so. Expression (26) shows strong analogy with the
continuum theoretical tangent tensor
$$ \frac{\partial \bar{\sigma}}{\partial \bar{\varepsilon}} = \mathcal{D} - \frac{\mathcal{D} : \frac{\partial f}{\partial \bar{\sigma}} \otimes \frac{\partial f}{\partial \bar{\varepsilon}} : \mathcal{D}}{\frac{\partial f}{\partial \bar{\varepsilon}} : \mathcal{D} : \frac{\partial f}{\partial \bar{\sigma}} + \frac{\partial f}{\partial \bar{a}} : \mathcal{G}(\bar{a}) : \frac{\partial f}{\partial \bar{a}} + G(a) \left(\frac{\partial f}{\partial a}\right)^2} \quad (27) $$

Also, the consistency of tensor $(\partial\bar{\sigma}/\partial\bar{\varepsilon})_{\text{alg}}$ is directly confirmed according to the limit properties discussed in the previous section when $\Delta\lambda$ tends to zero.

## Remarks

(i) The proposed formulation can be extended to non-associated plasticity by using a potential $F$ different from the yield function $f$ in the evolution laws (5). Then one observes that the new tensors $\mathcal{G}_{\text{alg}}$ and $\mathcal{G}_{\text{alg}}$—equivalent to (16a) and (16b)—remain symmetric but the new tangent modulus—equivalent to (25) or (26)—does not, as expected.

(ii) Only two hardening variables have been considered as usually done in plasticity, one second-order tensor and one scalar function. In more elaborate plasticity models or in damage mechanics where a larger number of internal variables is required, it is more difficult to derive closed-form expressions like (18), (19) and (26) for the local iteration corrections and the consistent tangent modulus.

## 4. Numerical results

The above local integration procedure is now incorporated in a three-dimensional finite element program. The description of a given hardening model must include the first and second derivatives of the yield function $f$, Eq. (9), as well as the derivatives of the hardening energy $w^{\text{hard}}$, Eq. (2), in order to define the plastic quantities (7) and the algorithmic tensors (16) and (20).

Elasticity is assumed to be isotropic, characterized by the Young modulus $E$ and the Poisson ratio $\nu$, and the hardening energy and the yield function are given by [3,5]

$$ \begin{aligned} w^{\text{hard}} &= w^{\bar{\sigma}}(\bar{\alpha}) = \frac{1}{2} \frac{2H}{3} \bar{\alpha} : \bar{\alpha} \\ f(\bar{\alpha}, \bar{a}) &= \sqrt{\frac{3}{2} \|\bar{s} - \bar{a}\|} + \frac{1}{2} \frac{3K}{2H} \bar{a} : \bar{a} - \sigma_0 \end{aligned} \quad (28) $$

where $\bar{s}$ denotes the deviator of stress $\bar{\sigma}$, $\sigma_0$ the initial yield stress, $H$ the kinematic hardening modulus and $K$ the saturation velocity. To investigate cyclically loaded structures, the model (28) is commonly used in the frame of non-associative plasticity, that is, the evolution laws

(5) are derived from potential $f$ given by (28), whereas the yield function is given by (28) with $K=0$. Here, the model (28) is used in the frame of standard associative plasticity which means that the potential coincide with the yield function given by (28). Although any combination of non-linear isotropic and non-linear kinematic hardenings is possible, the numerical results presented in this paper only involve a pure kinematic non-linear hardening, as it is enough to describe the ratcheting phenomenon.

The material properties have the same values for all the numerical examples, see Table 1.

Throughout, the finite element computations are carried out using three-dimensional 20-node hexahedral elements with $3 \times 3 \times 3 = 27$ Gauss integration points.

### 4.1. Uniaxial ratcheting effect

The first numerical computation is performed on a bar of constant section under uniaxial loading. The uniaxial stress is applied cyclically from $\sigma_m = -2 \times 10^8$ Pa to $\sigma_M = 4.5 \times 10^8$ Pa corresponding to an asymmetric loading process. The stress–strain curve resulting from the tension-compression cyclic loading is shown in Fig. 1.

The non-zero mean load has led to accumulation of the plastic strain which is known as the uniaxial ratcheting. The ratcheting step, corresponding to the accumulation of plastic strain during one cycle, only depends on the material properties and loading process, and is given by

$$ \Delta\varepsilon^p = \frac{1}{K} \ln \left[ \frac{\left(\frac{H}{K}\right)^2 - a_m^2}{\left(\frac{H}{K}\right)^2 - a_M^2} \right] \quad (29) $$

where

$$ a_m = \frac{-1 + \sqrt{1 + 2(\sigma_m + \sigma_0)\frac{K}{H}}}{\frac{K}{H}} \quad (30) $$

$$ a_M = \frac{1 - \sqrt{1 - 2(\sigma_M - \sigma_0)\frac{K}{H}}}{\frac{K}{H}} $$

The computed value $\Delta\varepsilon^p = 1.88 \times 10^{-3}$ is identical with the analytical solution within 1%.

### 4.2. Multiaxial ratcheting effect

Now consider a bar of rectangular section subjected to a constant tension $\sigma_{11} = \sigma$ and symmetrically alternated

<table><caption>Table 1
Material properties</caption><thead><tr><th>Young modulus <i>E</i></th><th>Poisson ratio <i>ν</i></th><th>Initial yield stress <i>σ</i><sub>0</sub></th><th>Kinematic hardening modulus <i>H</i></th><th>Saturation velocity <i>K</i></th></tr></thead><tbody><tr><td>2×10<sup>11</sup> Pa</td><td>0.3</td><td>3×10<sup>8</sup> Pa</td><td>3×10<sup>10</sup> Pa</td><td>60</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Uniaxial ratcheting, axial stress σ versus axial strain ε.</figcaption></figure>

shear stress σ₁₂, going from −τ to τ, which represents the coupling between traction and cyclic twisting. Here the stresses σ and τ take the same value 2×10⁸ Pa.

The stress–strain curve, Fig. 2, resulting from the cyclic loading shows that combining a constant tensile stress with a symmetric alternated shear yields accumulation of the axial plastic strain. The multiaxial ratcheting is essentially due to the non-zero value of the axial stress.

In this case, a constant ratcheting step is reached after a few cycles corresponding to a transitional stage, subsequent to the progressive accommodation encountered in the case of linear kinematic hardening.

Since it is difficult to obtain an analytical relation like Eq. (29), we shall confine ourselves to an approximate expression for the limit ratcheting step. This can be done

by calculating the evolution laws (5) with potential f given by (28), whereas the yield function is given by (28) with K = 0 (this is what is done in non-associated plasticity). One obtains the following relation between the variations of the two plastic strains ε̇₁₁ and ε̇₁₂

$$ \Delta\epsilon_{11}^{p} = \frac{4}{\sqrt{3}} \frac{\sigma}{\sqrt{(\frac{H}{K} + \sigma_{0})^{2} - \sigma^{2}}} \Delta\epsilon_{12}^{p} = 0.596\Delta\epsilon_{12}^{p} \quad (31) $$

The numerical factor 0.612 obtained is 3% higher than the value given by (31).

### 4.3. Cylinder under internal pressure

Let us consider a thick-walled cylinder of internal and external radii Rᵢ = R = 1 m and Rₑ = 2R = 2 m. Boundary conditions are applied so as to realize the plane strain conditions. The mesh of one quarter of the cylinder contains 503 nodes and 60 elements (10 elements along the radial direction by 6 elements along the circumferential direction). The cylinder is subjected to a cyclic internal pressure, going from p = 0 to a given maximum pressure p = pmax = 430 MPa. Fig. 3 depicts the variation of pressure p and loop stress σ₀₀ versus the evolution of inner radial displacement.

Since pmax is sufficiently high to induce plasticity during the unloading, a ratcheting phenomenon occurs with a nearly constant ratcheting step, similar to Fig. 1. However, the situation here is different from the previous ones, as the stress state is not homogeneous in the cylinder and the plastic zone spreads out progressively in the radial direction. The peak in σ₀₀ in Fig. 3 is not proper to non-linear kinematic hardening, it also occurs for the same cylinder in perfect plasticity.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Multiaxial ratcheting, shear σ₁₂ versus axial strain ε₁₁.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Thick-walled cylinder, internal pressure p and loop stress σ₀₀ versus inner radial displacement u(R).</figcaption></figure>
### 4.4. Beam under bending moment

The last numerical example is related to a beam of square section (the width and the height are equal, $b = h = 2$ cm, and the length is $l = 10$ cm), subjected to a bending moment $M$ applied at each end so that the beam undergoes a pure bending. The beam is discretized as a three-dimensional body with 20-node hexahedral elements. The mesh is made of eight elements in length and six elements in height to take account of the non-homogeneous stress state in the section.

Loading and unloading cycles are computed between $M_m = -600$ Nm and $M_M = 800$ Nm. The load-displacement curve, Fig. 4, where $M_0$ is the moment corresponding to the onset of plastification, clearly shows that ratcheting occurs with a nearly constant ratcheting step. The result is similar to that obtained by Mahbadi and Eslami [17] using a one-dimensional beam model. The maximum deflection at the center of the beam increases regularly and always takes positive values (except for the first cycle), see Table 2. Thus, the beam remains in the upper half-plane during the bending process.

All the numerical results obtained above have been compared to analytical or numerical models only, when

<figure><img src="image_5.png" /><figcaption>Fig. 4. Beam under pure bending, load-displacement curve: moment $M$ versus maximum deflection $v$.</figcaption></figure>

<table><caption>Table 2
Beam under pure bending, extremum deflections at extremum moments</caption><tr><td>M/M<sub>0</sub></td><td>2</td><td>-1.5</td><td>2</td></tr><tr><td>v/l</td><td>1.24 × 10<sup>-2</sup></td><td>-1.53 × 10<sup>-3</sup></td><td>1.43 × 10<sup>-2</sup></td></tr><tr><td>M/M<sub>0</sub></td><td>-1.5</td><td>2</td><td>-1.5</td></tr><tr><td>v/l</td><td>3.65 × 10<sup>-4</sup></td><td>1.61 × 10<sup>-2</sup></td><td>2.16 × 10<sup>-3</sup></td></tr></table>

these are available. It is also necessary to experimentally identify the mechanical properties of the theoretical model in order to assess its validity. This important work is not considered in this paper.

## 5. Conclusions

In this paper, a general local integration procedure in small strain elastoplasticity has been analyzed together with the corresponding consistent tangent modulus. The proposed formulation—established within the framework of the generalized standard materials theory—utilizes a minimal number of assumptions, and is general enough to be applied to arbitrary hardening models. The salient features of the present work are summarized below.

The yield function has been decomposed as the sum of one term depending on the stress and the back-stress tensor, and another one depending on a scalar. Here the stress and the back stress do not necessarily appear in the form of their difference, as seen in the case of non-linear kinematic hardening.

When the hardening energy $w^\alpha(\bar{\alpha})$ in relation (2) take general expressions other than a quadratic one—similarly to $w^\alpha(\alpha)$ in non-linear isotropic hardening—it has been shown that there appears the additional term $\frac{\partial f}{\partial \bar{\alpha}} : \frac{\partial^3 w^\alpha}{\partial \bar{\alpha}^3}$ in Eqs. (16a) and (16b), which does not exist in the literature. This term is a *symmetric* fourth-order tensor which should be taken into account in general theoretical considerations.

The plastic corrections have been given in closed forms, Eqs. (15)–(19). The consistent tangent modulus has a symmetrical expression, Eqs. (25) or (26). All these expressions have been given in tensorial forms which are ready to be converted to matrical expressions for computational use.

The expressions for the plastic corrections and the consistent tangent modulus involve the scalar function $G_{alg}$ and the fourth-order tensors $D_{alg}$, $G_{alg}$ and $L_{alg}$, the first two being *symmetric*. These algorithmic quantities represent the discrete finite-step counterparts of the usual theoretical continuum quantities: the plastic scalar function $G$, the elastic tensor $D$ and the plastic tensor $G$, respectively. When the time step tends to zero, all the algorithmic quantities tend to their theoretical counterparts.

The numerical computations using a particular non-linear kinematic hardening model have provided results which are in good agreement with analytical or other numerical solutions available in the literature. They also show that non-linear kinematic hardening models are effective in reproducing the important ratcheting phenomenon, responsible for one of the structural failure cases.
## References

[1] Alfano G, De Angelis F, Rosati L. Closed-form solution of the finite-step constitutive problem in plane stress $j_2$-elastoviscoplasticity. ECCM'99, European Conference on Computational Mechanics, München, Germany, 31 August-3 September, 1999.

[2] Alfano G, Rosati L. A general approach to the evaluation of consistent tangent operators for rate-independent elastoplasticity. Comput Meth Appl Mech Eng 1998;167:75-89.

[3] Armstrong PJ, Frederick CO. A mathematical representation of the multiaxial Bauschinger effect. Technical Report RD/B/N731, General Electricity Generating Board, Berkeley Nuclear Laboratories, 1966.

[4] Bathe KJ. Finite element procedures. Englewood Cliffs, NJ: Prentice Hall; 1996.

[5] Chaboche JL. Unified cyclic viscoplastic constitutive equations: development, capabilities, and thermodynamic framework. In: Krausz AS, Krausz K, editors. Unified constitutive laws of plastic deformation. Academic Press; 1996. p. 1-68.

[6] Doghri I. Fully implicit integration and consistent tangent modulus in elastoplasticity. Int J Numer Meth Eng 1993;36:3915-32.

[7] Doghri I. Numerical implementation and analysis of a class of metal plasticity models coupled with ductile damage. Int J Numer Meth Eng 1995;38:3403-31.

[8] Doghri I. Mechanics of deformable solids: linear and nonlinear, analytical and computational aspects. Springer; 2000.

[9] Eterovic AL, Bathe KJ. A hyperelastic-based large strain elasto-plastic constitutive formulation with combined isotropic-kinematic hardening using the logarithmic stress and strain measures. Int J Numer Meth Eng 1990;30:1099-114.

[10] Halphen B, Nguyen QS. Sur les matériaux standard généralisés. J Mécanique 1975;14(1):39-63.

[11] Hartmann S, Haupt P. Stress computation and consistent tangent operator using non-linear kinematic hardening models. Int J Numer Meth Eng 1993;36:3801-14.

[12] Hughes TJR, Pister KS. Consistent linearization in mechanics of solids and structures. Comput Struct 1978; 9:391-7.

[13] Hughes TJR, Taylor RL. Unconditionally stable algorithms for quasi-static elasto-viscoplastic finite element analysis. Comput Struct 1978;8:169-73.

[14] Kojic M, Bathe KJ. The “effective-stress-function” algorithm for thermo-elasto-plastic and creep. Int J Numer Meth Eng 1987;24:1509-32.

[15] Krieg RD, Key SW. Implementation of a time dependent plasticity theory into structural computer programs. In: Stricklin JA, Saczalski KJ, editors. Constitutive equations in viscoplasticity: computational and engineering aspects, AMD-20. New York: ASME; 1976.

[16] Krieg RD, Krieg DB. Accuracies of numerical solution methods for the elastic-perfectly plastic model. J Press Vess Technol, ASME 1977;99:510-5.

[17] Mahbadi H, Eslami MR. Cyclic loading of beams based on the Prager and Frederick-Armstrong kinematic hardening models. Int J Mech Sci 2002;44:859-79.

[18] Matzenmiller A, Taylor RL. A return mapping algorithm for isotropic elastoplasticity. Int J Numer Meth Eng 1994; 37:813-26.

[19] Moreau JJ. Proximité et dualité dans un espace hilbertien. Bull Soc Math France 1965;93:273-99.

[20] Moreau JJ. On unilateral constraints, friction and plasticity. Lecture notes, CIME, Bressanone, Cremonese, Roma, 1973.

[21] Nagtegaal JC. On the implementation of inelastic constitutive equations with special reference to large deformation problems. Comput Meth Appl Mech Eng 1982;33:469-84.

[22] Nguyen QS. On the elastic-plastic initial-boundary value problem and its numerical integration. Int J Numer Meth Eng 1977;11:817-32.

[23] Ortiz M, Martin JB. Symmetry-preserving return mapping algorithms and incrementally extremal paths: a unification of concepts. Int J Numer Meth Eng 1989;28:1839-53.

[24] Ortiz M, Pinski PM, Taylor RL. Operator split methods for the numerical solution of the elastoplastic dynamic problem. Comput Meth Appl Mech Eng 1983;39:137-57.

[25] Ortiz M, Popov EP. Accuracy and stability of integration algorithms for elastoplastic constitutive relations. Int J Numer Meth Eng 1985;21:1561-76.

[26] Ortiz M, Simo JC. An analysis of a new class of integration algorithms for elastoplastic constitutive equations. Int J Numer Meth Eng 1986;23:353-66.

[27] Palazzo V, Rosati L, Valoroso N. Computational issues of general isotropic elastoplastic models. In: VI International Conference on Computational Plasticity (COMPLAS VI), Barcelona, Spain, 11-14 September, 2000.

[28] Reddy BD, Martin JB. Internal variable formulations of problems in elastoplasticity: constitutive and algorithmic aspect. Appl Mech Rev 1994;47(9):429-56.

[29] Simo JC, Govindjee S. Exact closed-formed solution of the return mapping algorithm for plane stress elasto-viscoplasticity. Eng Comput 1988;3:254-8.

[30] Simo JC, Govindjee S. Non-linear b-stability and symmetry preserving return mapping algorithms for plasticity and viscoplasticity. Int J Numer Meth Eng 1991;31:151-76.

[31] Simo JC, Hughes TJR. Computational inelasticity. Springer; 1998.

[32] Simo JC, Hughes TJR. General return mapping algorithms for rate-independent plasticity. In: Desai CS et al., editors. 2nd Int Conf on Constitutive Laws for Eng Materials: Theory and Applications, vol. I. Tucson, Arizona, USA, 5-8 January, 1987. p. 221-31.

[33] Simo JC, Ortiz M. A unified approach to finite deformation elastoplasticity based on the use of hyperelastic constitutive equations. Comput Meth Appl Mech Eng 1985;49:221-45.

[34] Simo JC, Taylor RL. Consistent tangent operators for rate-independent elastoplasticity. Comput Meth Appl Mech Eng 1985;48:101-18.

[35] Simo JC, Taylor RL. A return mapping algorithm for plane stress elastoplasticity. Int J Numer Meth Eng 1986; 22:649-70.

[36] Le van A, de Saxcé G. Exact expressions for the plastic corrections and the consistent tangent modulus in rate-independent elastoplasticity. In: VI International Conference on Computational Plasticity (COMPLAS VI), Barcelona, Spain, 11-14 September, 2000.
[37] Le van A. Une formulation lagrangienne totale en élastoplasticité finie. Habilitation Thesis. Ecole Centrale de Nantes, Université de Nantes, France, 1999.

[38] Wilkins ML. Calculation of elastic-plastic flow. In: Adler B et al., editors. Methods of computational physics 3. New York: Academic Press; 1964.