# Classical and sequential limit analysis revisited

Jean-Baptiste Leblond, Djimedo Kondo, Léo Morin, Almahdi Remmal

► To cite this version:

Jean-Baptiste Leblond, Djimedo Kondo, Léo Morin, Almahdi Remmal. Classical and sequential limit analysis revisited. Comptes Rendus Mécanique, 2018, 346 (4), pp.336-349. 10.1016/j.crme.2017.12.015. hal-01768417

HAL Id: hal-01768417

https://hal.sorbonne-universite.fr/hal-01768417

Submitted on 17 Apr 2018

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
In order to answer questions (1) and (2), the principle of sequential limit analysis consists in considering a hardenable material in large strain as a discrete sequence of different, successive ideal plastic materials occupying different, successive geometrical configurations. In this approach, at a given discretized instant, the hardening and the geometry are considered as momentarily fixed. The hardenable material, considered within a fully general geometrical framework, then behaves like an ideal plastic material, considered within a linearized geometrical framework, the role of the previous hardening being merely to modify the instantaneous local yield criterion and the plastic flow rule, and that of the previous geometry changes to fix the present geometrical configuration. An "overall yield locus", supplemented by an "overall associated plastic flow rule", can then be determined using the standard theorems of classical limit analysis.

To answer question (3), the evolutions of the distribution of hardening parameters and the geometrical configuration are deduced a *posteriori*, by approximately updating these quantities using the trial velocity field used in the limit analysis, integrated within a time step.

The proposal has been received with some enthusiasm and sequential limit analysis has been used for various applications; see for instance the works [2–6], to quote just a few.

However, the mechanical foundations of sequential limit analysis have remained vague and to some large extent unclear, and this may prompt doubts about the overall validity of the method. How can one ignore, even momentarily, the evolutions of the hardening parameters and the geometry, whereas in reality they change constantly? Clearly, in order to settle such a question, it is necessary to scrutinize the connection between the evolution equations of the full elastic-plastic problem and the simpler equations of sequential limit analysis.

It will also be necessary to investigate the connection with the equations of classical limit analysis, because sequential limit analysis uses the results and methods of classical limit analysis, which is justified only provided the equations of both theories are analogous. Re-examining the equations of classical limit analysis will also be useful to re-assess, and remind the reader of the role played by elasticity in this theory. Indeed although Drucker et al. [7] have clearly shown that classical limit analysis perfectly tolerates the presence of elasticity, the memory of their work seems to have somewhat faded in time.

The paper is organized as follows.

* Section 2 is devoted to a short summary of the classical works the classical works of Drucker et al. [7] and Hill [8]. The discussion concentrates on those results of the theory actually used in its practical application. The discussion of the role of elasticity follows the work of Drucker et al. [7].

* Section 3 then discusses the foundations of sequential limit analysis, exploiting an analogy of its equations with those of classical limit analysis and using some, if not all, elements of Section 2. Special attention is again paid to the role of elasticity. The precise conditions of applicability of sequential limit analysis are made clear.

* Finally Section 4 considers, as a typical example, the case of porous ductile materials, which involves all possible complexities: elasticity, plasticity, strain hardening, and large strains. Using the results of Section 3, we discuss the applicability of sequential limit analysis to the derivation, through homogenization, of micromechanically-based models.

## 2. Classical limit analysis

### 2.1. *Hypotheses and notations – equations of the problem*

We consider a body $\Omega$ made of some *elastic–ideal plastic* material obeying the plastic flow rule associated with the yield criterion through Hill’s normality property. The elastic–plastic evolution problem is considered within a *linearized geometrical framework*, so that its equations are written using the initial position vector $\mathbf{X}$, the linearized strain tensor $\boldsymbol{\epsilon}$, and the Cauchy stress tensor $\boldsymbol{\sigma}$. With these hypotheses these equations read as follows (disregarding body forces for simplicity):

$$ 
\left\{
\begin{array}{ll}
\mathbf{div}_{\mathbf{X}}\boldsymbol{\sigma} = \mathbf{0} & \text{equilibrium} \\
\dot{\boldsymbol{\epsilon}} = \frac{1}{2}[\mathbf{grad}_{\mathbf{X}}\dot{\mathbf{u}} + (\mathbf{grad}_{\mathbf{X}}\dot{\mathbf{u}})^{\mathrm{T}}] & \text{definition of } \dot{\boldsymbol{\epsilon}} \\
\dot{\boldsymbol{\epsilon}} = \dot{\boldsymbol{\epsilon}}^{\mathrm{e}} + \dot{\boldsymbol{\epsilon}}^{\mathrm{P}} & \text{additive decomposition of } \dot{\boldsymbol{\epsilon}} \\
\dot{\boldsymbol{\epsilon}}^{\mathrm{e}} = \mathbf{S} : \dot{\boldsymbol{\sigma}} & \text{elasticity law} \\
\dot{\boldsymbol{\epsilon}}^{\mathrm{P}} = \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}}(\boldsymbol{\sigma}) & \text{plastic flow rule obeying the normality property} \\
f(\boldsymbol{\sigma}) \le 0, \quad \dot{\lambda} \ge 0, \quad f(\boldsymbol{\sigma})\dot{\lambda} = 0 & \text{Kuhn–Tucker’s complementarity conditions} \\
\text{B.C.} & \text{boundary conditions}
\end{array}
\right.
\quad (1)
$$

In these equations, $\dot{\mathbf{u}}$ denotes the velocity vector ($\mathbf{u}$ is the displacement vector), $\dot{\boldsymbol{\epsilon}}^{\mathrm{e}}$ the elastic strain rate tensor, $\dot{\boldsymbol{\epsilon}}^{\mathrm{P}}$ the plastic strain rate tensor, $\mathbf{S}$ the elastic compliance tensor, $f(\boldsymbol{\sigma})$ the yield function (depending only on $\boldsymbol{\sigma}$ in the absence of hardening), and $\dot{\lambda}$ the plastic multiplier. The boundary conditions, symbolized by the letters “B.C.”, need not be stated explicitly here; it will suffice to say that they depend linearly upon a finite number of load parameters $Q_1, Q_2, \dots, Q_N$, collectively defining a “load vector” $\mathbf{Q}$. The conjugate kinematic parameters $q_1, q_2, \dots, q_N$, collectively defining a “kinematic vector” $\mathbf{q}$, are defined in such a way that the virtual power of the external forces is $\sum_{i=1}^{N} Q_i \dot{q}_i = \mathbf{Q} \cdot \dot{\mathbf{q}}$.
Classical limit analysis studies loading paths ending at *limit loads* defined, following Drucker et al. [7], by the two conditions

$$ \dot{\mathbf{Q}}=\mathbf{0} ; \quad \dot{\mathbf{q}} \neq \mathbf{0} \qquad (2) $$

which mean that when the load is reached, it becomes stationary, while the structure keeps on deforming.<fnref n="1" />

## 2.2. A lemma of Drucker et al. [7]

Although Hill's [8] classical presentation of limit analysis was based on the assumption of the absence of elasticity, it has been clearly shown by Drucker et al. [7] that this hypothesis is by no means necessary. This stems from the following lemma.

**Lemma 2.1** (Drucker et al. [7]). When an arbitrary limit load **Q** is reached, the elastic strain rate $\dot{\epsilon}^e$ and the stress rate $\dot{\sigma}$ are zero at every point of the structure.

**Proof.** If a limit load is reached, combining condition (2)₁, the definition of conjugate parameters, and the principle of virtual work applied to the stress rate field {$\dot{\sigma}$} and the strain rate field {$\dot{\epsilon}$}, one gets

$$ 0 = \dot{\mathbf{Q}} \cdot \dot{\mathbf{q}} = \int_{\Omega} \dot{\sigma} : \dot{\epsilon} \, d\Omega $$

One gets from there, using the decomposition of the strain rate (1)₃, the elasticity law (1)₄ and the plastic flow rule (1)₅:

$$ 0 = \int_{\Omega} \dot{\sigma} : \dot{\epsilon}^e \, d\Omega + \int_{\Omega} \dot{\sigma} : \dot{\epsilon}^P \, d\Omega = \int_{\Omega} \dot{\sigma} : S : \dot{\sigma} \, d\Omega + \int_{\Omega} \dot{\sigma} : \dot{\lambda} \frac{\partial f}{\partial \dot{\sigma}}(\sigma) \, d\Omega $$

But, in the absence of any hardening parameter, $f$ depends only on $\sigma$; thus $\dot{\sigma} : \frac{\partial f}{\partial \sigma}(\sigma) = \frac{d}{dt}[f(\sigma)] = \dot{f}$, so that the last integral may be rewritten as $\int_{\Omega} \dot{\lambda} \dot{f} \, d\Omega$. Now at every point there are two possibilities: either $\dot{\lambda} = 0$ (elastic behavior),<fnref n="2" /> or $\dot{\lambda} > 0$ (plastic behavior) but then $\dot{f} = 0$ (the yield condition is satisfied and remains so). In both cases, $\dot{\lambda} \dot{f} = 0$. It follows that the integral $\int_{\Omega} \dot{\sigma} : \dot{\lambda} \frac{\partial f}{\partial \dot{\sigma}}(\sigma) \, d\Omega = \int_{\Omega} \dot{\lambda} \dot{f} \, d\Omega$ is zero, so that the preceding equation yields

$$ \int_{\Omega} \dot{\sigma} : S : \dot{\sigma} \, d\Omega = 0 $$

But the compliance tensor **S** is positive-definite, so that the quantity $\dot{\sigma} : S : \dot{\sigma}$ is non-negative at every point. Therefore, its integral over $\Omega$ can only be zero if it is uniformly zero, which implies that $\dot{\sigma} = \mathbf{0}$ and therefore $\dot{\epsilon}^e = \mathbf{0}$ at every point. □

With Lemma 2.1, at a limit load, equations (1) of the problem may be re-written in the following form:

$$ \left\{ \begin{array}{l} \dot{\mathbf{q}} \neq \mathbf{0} \\ \mathrm{div}_{\chi} \sigma = \mathbf{0} \\ \dot{\epsilon} \equiv \frac{1}{2} [\mathrm{grad}_{\chi} \dot{\mathbf{u}} + (\mathrm{grad}_{\chi} \dot{\mathbf{u}})^{\mathrm{T}}] = \dot{\epsilon}^{\mathrm{P}} \\ \dot{\epsilon}^{\mathrm{P}} = \dot{\lambda} \frac{\partial f}{\partial \dot{\sigma}}(\sigma) \\ f(\sigma) \le 0, \quad \dot{\lambda} \ge 0, \quad f(\sigma) \dot{\lambda} = 0 \\ \text{B.C.} \end{array} \right. \qquad (3) $$

where the condition $\dot{\mathbf{Q}}=\mathbf{0}$ has been left out because it does not play any role beyond the justification of the elimination of $\dot{\epsilon}^e$ from the equations.

<footnote n="1">¹ Condition (2)₁ applies to most loading paths but rules out special ones. For instance, for some paths **Q** moves along the set of limit loads; then $\dot{\mathbf{Q}}$ never vanishes. A limit load may also be reached only asymptotically, so that $\dot{\mathbf{Q}}$ never truly vanishes, although it continuously decreases in magnitude.</footnote>

<footnote n="2">² This occurrence is perfectly possible at some points; a structure does not necessarily have to be entirely plastic for a limit load to be reached. Simple systems of bars provide elementary examples. Another, more elaborate one is supplied by the coalescence of voids in porous ductile solids, wherein the plastic strain rate suddenly concentrates within thin layers, the regions between these layers undergoing elastic unloading.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of the proof of consequence 1 of Hill's [8] lemma.</figcaption></figure>

## 2.3. A lemma of Hill [8]

We are now in a position to establish a classical lemma of Hill [8] without resorting to his unnecessary hypothesis of absence of elasticity.

The set K of potentially sustainable loads is defined as consisting of those loads $Q^*$ for which there exists a stress field $\{\sigma^*\}$ which is both statically admissible with $Q^*$ (it obeys the equilibrium equations and the boundary conditions with the given value of $Q^*$) and plastically admissible (it satisfies the condition $f(\sigma^*) \le 0$ at every point). Then we obtain the following lemma.

**Lemma 2.2** (Hill [8]). Let $Q$ denote a limit load, $\dot{q}$ the rate of the kinematic vector it generates along the loading path considered, and $Q^* \in K$ an arbitrary potentially sustainable load. Then

$$ (Q - Q^*) \cdot \dot{q} \ge 0 \qquad (4) $$

**Proof.** Let $\{\sigma\}$ and $\{\dot{\epsilon}\}$ denote the stress and strain rate fields generated by the load $Q$, and $\{\sigma^*\}$ the stress field associated with the load $Q^*$ in its definition as a potentially sustainable load. Then, by the principle of virtual work applied to the stress field $\{\sigma\} - \{\sigma^*\}$ and the strain rate field $\{\dot{\epsilon}\}$, plus Lemma 2.1,

$$ (Q - Q^*) \cdot \dot{q} = \int_{\Omega} (\sigma - \sigma^*) : \dot{\epsilon} \, d\Omega = \int_{\Omega} (\sigma - \sigma^*) : \dot{\epsilon}^P \, d\Omega $$

But $(\sigma - \sigma^*) : \dot{\epsilon}^P \ge 0$ at every point by Hill's principle (equivalent to the normality law). The result follows. □

**Remarks.** 2.1. Inequality (4) stands as a kind of “overall equivalent” of Hill's local principle. There is however a significant difference arising from the presence of elasticity. Hill's principle applies to any plastically admissible stress tensor $\sigma$, even if it gives rise to a purely elastic behavior; indeed if $f(\sigma) < 0$, then $\dot{\epsilon}^P = 0$, so that $(\sigma - \sigma^*) : \dot{\epsilon}^P = 0$. In contrast inequality (4) stands for loads $Q$ that cannot be exceeded, but says nothing about loads that can; indeed if $Q$ can be exceeded, $\dot{\epsilon}^e$ has no reason to be zero at every point, so that the quantity $(Q - Q^*) \cdot \dot{q}$ differs from $\int_{\Omega} (\sigma - \sigma^*) : \dot{\epsilon}^P \, d\Omega$, and nothing can be said about its sign. This is an illustration of the subtleties brought by elasticity in limit analysis.

2.2. The only difference between the proof of inequality (4) given here and that provided by Hill [8] lies in the justification of the equality $\dot{\epsilon} = \dot{\epsilon}^P$, based on Drucker's Lemma 2.1 here, and on the mere neglect of elasticity in [8].

## 2.4. Consequences of Hill's [8] lemma

**Consequence 2.1.** All limit loads necessarily lie on the boundary of the set K of potentially sustainable loads.

**Proof.** Fig. 1 illustrates the elementary geometric reasoning leading to this conclusion, in the simple case where the vector $Q$ of load parameters consists of two components only. Obviously, all limit loads lie in $K$. Assume that the above statement is wrong and that one of them, $Q$, does not lie on the boundary of this set; it then lies in its interior. Therefore, there exists an open ball $B(Q, \eta)$ of center $Q$ and radius $\eta > 0$ entirely contained inside $K$. By Lemma 2.2, for every $Q^* \in B(Q, \eta)$, $(Q - Q^*) \cdot \dot{q} \ge 0$, where $\dot{q}$ denotes the rate of the kinematic vector generated by $Q$. But when $Q^*$ spans $B(Q, \eta)$, the vector $Q - Q^*$ spans all the directions of the $N$-dimensional space; therefore, the inequality $X \cdot \dot{q} \ge 0$ holds for all vectors $X$ in this space. Applying it to the pair $(X, -X)$, one sees that in fact $X \cdot \dot{q} = 0$ for every $X$, which implies that necessarily $\dot{q} = 0$. But this contradicts condition (2)$_2$. □
This result, of fundamental importance in limit analysis, ascribes to the set K of potentially sustainable loads, defined a priori as a purely mathematical object, the essential physical role of being the set of actually sustainable loads: all load paths governed by the equations of the elastic–plastic problem must end on its boundary, not at some interior point.

The clearest statement of the second consequence is probably that provided by Mandel [9]:

**Consequence 2.2.** (i) The set K is convex; (ii) if **Q** is a limit load (necessarily lying on the boundary of K), the corresponding rate **q̇** of the kinematic vector is directed along some local exterior normal to K.

The geometrical proof is very classical (see, e.g., Salençon's [10] textbook on plasticity and limit analysis) and need not be recalled. Note that property (i) does not in fact bring anything new, since the convexity of K was obvious from the start from its definition, the linearity of the equilibrium equations and the boundary conditions, and the convexity of the local reversibility domain {$\sigma^*|f(\sigma^*) \le 0$}.

Assume that an explicit analytic expression of the set K has been found in the form of an inequality of the form $F(\mathbf{Q}) \le 0$ for some “overall yield function” F. If this function is regular at the point $\mathbf{Q}$ (the boundary of K admits a local tangent hyperplane and therefore a single local normal), property (ii) above means that there exists a scalar $\dot{\Lambda}$ – the “overall plastic multiplier” – such that

$$ \dot{q} = \dot{\Lambda} \frac{\partial F}{\partial \mathbf{Q}}(\mathbf{Q}), \quad \dot{\Lambda} \ge 0 \qquad (5) $$

If F is not regular at the point $\mathbf{Q}$ (the boundary of K admits a local cone of normals), this equation must be replaced by the more general one

$$ \dot{q} \in \partial F(\mathbf{Q}) \qquad (6) $$

where $\partial F(\mathbf{Q})$ is the sub-differential of F at the point $\mathbf{Q}$ (Moreau [11]).

## 2.5. Practical procedures of classical limit analysis

The results of Subsection 2.4 show that, in order to study limit loads and the corresponding overall directions of the flow, it suffices to determine the set K of potentially sustainable loads, which depends upon the constitutive law of the material only through its yield criterion.

The static or interior approach (Hill [8]) provides the most intuitive, if not the most practically efficient, way of doing so. It simply consists in determining some elements of K by exhibiting stress fields statically admissible with them and plastically admissible. If $\mathbf{Q}$ denotes such an element, the convexity of K implies that the whole segment $[0, \mathbf{Q}]$ is contained inside this set, so that the limit load lying on the straight half-line originating from the origin and passing through $\mathbf{Q}$ must lie beyond this point (whence the terminology “interior approach”).

The less intuitive but more efficient kinematic or exterior approach relies on another result of Hill [8]. This result makes use of the function $\pi(\dot{\epsilon})$ defined for any $\dot{\epsilon}$ by the formula

$$ \pi(\dot{\epsilon}) = \sup_{\sigma^*|f(\sigma^*) \le 0} \sigma^* : \dot{\epsilon} \qquad (7) $$

This function is the support function (the Legendre–Fenchel transform of the indicator function) of the convex set {$\sigma^*|f(\sigma^*) \le 0$}. (It is often called the plastic dissipation. This terminology is slightly inadequate in that $\dot{\epsilon}$ in equation (7) is a priori just an arbitrary symmetric second-rank tensor, not the plastic part of an actual strain rate tensor. The inadequacy will be seen to become more pronounced in sequential limit analysis, see Section 3 below.)

**The fundamental inequality of limit analysis** (Hill [8]). If $\mathbf{Q}^*$ denotes an arbitrary element of K, $\dot{\mathbf{q}}$ an arbitrary rate of kinematic parameters (unconnected to $\mathbf{Q}^*$), and $[\mathbf{v}]$ an arbitrary velocity field kinematically admissible with $\dot{\mathbf{q}}$ (compatible with the given value of this rate), generating a strain rate field {$\dot{\epsilon}$}, then

$$ \mathbf{Q}^* \cdot \dot{\mathbf{q}} \le \int_{\Omega} \pi(\dot{\epsilon}) \, d\Omega \qquad (8) $$

**Proof.** The result is a simple consequence of the definition of K, the principle of virtual work and the fact that $\sigma^* : \dot{\epsilon} \le \pi(\dot{\epsilon})$ for every $\sigma^*$ respecting the inequality $f(\sigma^*) \le 0$. $\square$

Inequality (8) provides “upper estimates”, in a vectorial sense, of the limit loads (whence the terminology “exterior approach”), provided that one is able to exhibit suitable velocity fields. (This generally requires more than simply satisfying the requirement of kinematic admissibility, because the function $\pi(\dot{\epsilon})$ is often infinite – making inequality (8) void of information – except for some specific strain rates; for instance, for yield criteria independent of the mean stress, $\dot{\epsilon}$ must be traceless, that is the velocity field must be incompressible, in order for $\pi(\dot{\epsilon})$ to be finite.)
### 3. Sequential limit analysis

#### 3.1. New hypotheses, notations, and evolution equations

We now come to the more general problem of a body $\Omega$ made of some elastic-plastic material subject to strain hardening, and envisaged within a general geometrical framework involving large displacements and strains. The collection of local (scalar, vectorial or tensorial) hardening parameters is collectively denoted by $\chi$, and the present position vector by $\mathbf{x}$. The equations of the evolution problem are written using a Eulerian setting, on the current configuration, and involve the Eulerian strain rate $\mathbf{d}$ and the Cauchy stress tensor $\boldsymbol{\sigma}$:

$$ 
\left\{
\begin{array}{ll}
\mathbf{div}_x \boldsymbol{\sigma} = \mathbf{0} & \text{equilibrium} \\
\mathbf{d} \equiv \frac{1}{2} [\mathbf{grad}_x \mathbf{v} + (\mathbf{grad}_x \mathbf{v})^\mathrm{T}] & \text{definition of } \mathbf{d} \\
\mathbf{d} = \mathbf{d}^\mathrm{e} + \mathbf{d}^\mathrm{P} & \text{additive decomposition of } \mathbf{d} \\
\mathbf{d}^\mathrm{e} = \mathbf{S} : \frac{\mathrm{D}\boldsymbol{\sigma}}{\mathrm{D}t} & \text{hypoelasticity law} \\
\mathbf{d}^\mathrm{P} = \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}}(\boldsymbol{\sigma}, \chi) & \text{plastic flow rule obeying normality} \\
f(\boldsymbol{\sigma}, \chi) \le 0, \quad \dot{\lambda} \ge 0, \quad f(\boldsymbol{\sigma}, \chi) \dot{\lambda} = 0 & \text{Kuhn-Tucker's complementarity conditions} \\
\frac{\mathrm{D}\chi}{\mathrm{D}t} = \dots & \text{evolution of hardening parameters} \\
\dot{\mathbf{x}} = \mathbf{v} & \text{evolution of geometry} \\
\text{B.C.} & \text{boundary conditions}
\end{array}
\right.
\qquad (9) $$

In these equations, $\mathbf{v}$ denotes the velocity, $\mathbf{d}^\mathrm{e}$ the elastic strain rate tensor, $\mathbf{d}^\mathrm{P}$ the plastic strain rate tensor, and the symbol $\frac{\mathrm{D}}{\mathrm{D}t}$ represents some objective time-derivative (for instance that of Jaumann). In comparison to system (1), system (9) encompasses two new equations, (9)$_7$ and (9)$_8$, providing the evolutions of the hardening parameters and the geometry. The precise expression of $\frac{\mathrm{D}}{\mathrm{D}t}$ in equation (9)$_7$ is not specified because it is immaterial here, but in practice the value of this derivative is determined by those of the plastic strain rate $\mathbf{d}^\mathrm{P}$ and the vector $\chi$ itself.

#### 3.2. Conditions for reduction to the equations of limit analysis

In general, there is no way to reduce the general evolution equations (9) to the simpler system of equations (3) of limit analysis, for at least two reasons:

* in the presence of strain hardening and/or geometry changes, the load vector $\mathbf{Q}$ has no reason to ever reach a stationary value, so that the condition $\dot{\mathbf{Q}} = \mathbf{0}$ forming the basis of Drucker et al.'s [7] proof of Lemma 2.1 of Subsection 2.2, is never satisfied;

* even if this condition were satisfied, Drucker et al.'s [7] proof would break down, because in the presence of hardening $\dot{\boldsymbol{\sigma}} : \frac{\partial f}{\partial \boldsymbol{\sigma}}(\boldsymbol{\sigma}, \chi)$ differs from $\dot{f}$ by the term $\dot{\boldsymbol{\chi}} : \frac{\partial f}{\partial \boldsymbol{\chi}}(\boldsymbol{\sigma}, \chi)$; hence the integral $\int_{\Omega} \dot{\boldsymbol{\sigma}} : \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}}(\boldsymbol{\sigma}) \mathrm{d}\Omega$ cannot be expressed as $\int_{\Omega} \dot{\lambda} \dot{f} \mathrm{d}\Omega$.

However, in some cases it is physically reasonable to completely disregard elasticity. (It is difficult to state general conditions ensuring the possibility of such a simplification, since its soundness simultaneously depends on the material, the geometry, and the loading; each individual case warrants a specific discussion. An example will be considered in detail in Section 4.) This means considering the limit $\mathbf{S} \to \mathbf{0}$ in equations (9); the system then reduces to

$$ 
\left\{
\begin{array}{ll}
\dot{\mathbf{q}} \neq \mathbf{0} & \\
\mathbf{div}_x \boldsymbol{\sigma} = \mathbf{0} & \\
\mathbf{d} \equiv \frac{1}{2} [\mathbf{grad}_x \mathbf{v} + (\mathbf{grad}_x \mathbf{v})^\mathrm{T}] = \mathbf{d}^\mathrm{P} & \\
\mathbf{d}^\mathrm{P} = \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}}(\boldsymbol{\sigma}, \chi) & \\
f(\boldsymbol{\sigma}, \chi) \le 0, \quad \dot{\lambda} \ge 0, \quad f(\boldsymbol{\sigma}, \chi) \dot{\lambda} = 0 & \\
\frac{\mathrm{D}\chi}{\mathrm{D}t} = \dots & \\
\dot{\mathbf{x}} = \mathbf{v} & \\
\text{B.C.} &
\end{array}
\right.
\qquad (10) $$

where the additional condition $\dot{\mathbf{q}} \neq \mathbf{0}$ simply states that we are interested only in loads $\mathbf{Q}$ generating an actual deformation of the structure. Now consider, within system (10), the subsystem consisting only of the first five equations plus the last one:
$$
\left\{
\begin{array}{l}
\dot{\mathbf{q}} \neq \mathbf{0} \\
\mathbf{div}_{\mathbf{x}} \sigma = \mathbf{0} \\
\mathbf{d} = \frac{1}{2} [\mathbf{grad}_{\mathbf{x}} \mathbf{v} + (\mathbf{grad}_{\mathbf{x}} \mathbf{v})^{\mathrm{T}}] = \mathbf{d}^{\mathrm{P}} \\
\mathbf{d}^{\mathrm{P}} = \dot{\lambda} \frac{\partial f}{\partial \sigma}(\sigma, \chi) \\
f(\sigma, \chi) \le 0, \quad \dot{\lambda} \ge 0, \quad f(\sigma, \chi) \dot{\lambda} = 0 \\
\text{B.C.}
\end{array}
\right.
\qquad (11)
$$

If considered over a finite time interval, this subsystem is incomplete, because it does not specify how the hardening parameters and the geometry evolve in time, whereas such features have an explicit impact (for the hardening parameters), or an implicit one (for the geometry), upon the equations. However, *assume that the full distribution of hardening parameters and the full present configuration are known at some given instant, and write system (11) at this sole instant*. This system is then complete in itself and does not require the expressions of $\frac{\partial \chi}{\partial \mathbf{D}}$ and $\dot{\mathbf{x}}$, since these derivatives do not enter the equations.

Moreover system (11), viewed in this light, is completely analogous to the system (3) of equations of classical limit analysis, with the correspondences $\mathbf{x} \leftrightarrow \mathbf{X}$, $\mathbf{v} \leftrightarrow \dot{\mathbf{u}}$, $\mathbf{d} \leftrightarrow \dot{\boldsymbol{\epsilon}}$, $\mathbf{d}^{\mathrm{P}} \leftrightarrow \dot{\boldsymbol{\epsilon}}^{\mathrm{P}}$; in the former system the hardening parameters and the geometry play the role of mere parameters. *This justifies – if the restrictive hypothesis of absence of elasticity is acceptable – the theory of sequential limit analysis:* that is the use of the results and methods of classical limit analysis, the distribution of hardening parameters and the geometry being considered as given and fixed, their role being merely to modify the yield criterion and prescribe the present configuration.

It is worth stressing, however, that the wording "sequential limit analysis", although employed in all papers using this theory – see the selection cited in the Introduction –, is particularly unhappy in that it gives the erroneous impression that it focuses on limit loads. Of course, in reality such loads do not exist in the presence of hardening and/or geometric changes. *The loads* $\mathbf{Q}$ *considered and studied by sequential limit analysis are simply those which – in the absence of elasticity – generate an actual deformation of the structure, $\dot{\mathbf{q}} \neq \mathbf{0}$, and have nothing to do with non-existent limit loads.*

### 3.3. Procedures of sequential limit analysis

The analogy between the equations of sequential and classical limit analysis just mentioned justifies – in the absence of elasticity – the following definitions and results:

* one defines a set $K_{G,\{\chi\}}$ consisting of those loads $\mathbf{Q}^*$ for which there exists a stress field $\{\boldsymbol{\sigma}^*\}$ statically admissible with $\mathbf{Q}^*$ and satisfying the condition $f(\boldsymbol{\sigma}, \chi) \le 0$ at every point. The subscripts $G$ and $\{\chi\}$ here underline the dependence of this set upon the full geometry symbolized by the letter $G$, and the full field $\{\chi\}$ of internal parameters;

* the loads $\mathbf{Q}$ generating an actual deformation of the structure at the instant considered, $\dot{\mathbf{q}} \neq \mathbf{0}$, necessarily lie on the boundary of the set $K_{G,\{\chi\}}$;

* assume that an explicit inequality $F_{G,\{\chi\}}(\mathbf{Q}) \le 0$ defining the set $K_{G,\{\chi\}}$ has been found. The rate $\dot{\mathbf{q}}$ of the kinematic vector generated by $\mathbf{Q}$ is then of the form

$$ \dot{\mathbf{q}} = \dot{\lambda} \frac{\partial F_{G,\{\chi\}}}{\partial \mathbf{Q}}(\mathbf{Q}) , \quad \dot{\lambda} \ge 0 \qquad (12) $$

or

$$ \dot{\mathbf{q}} \in \partial F_{G,\{\chi\}}(\mathbf{Q}) \qquad (13) $$

depending on whether the function $F_{G,\{\chi\}}$ is regular or not at the point $\mathbf{Q}$.

The set $K_{G,\{\chi\}}$ may be determined using both static and kinematic approaches. In the static approach, one merely determines elements of this set by exhibiting stress fields $\{\boldsymbol{\sigma}^*\}$ satisfying the conditions requested by its definition.

In the kinematic approach, one defines a function $\pi_{\chi}$ by the formula

$$ \pi_{\chi}(\mathbf{d}) = \sup_{\boldsymbol{\sigma}^* | f(\boldsymbol{\sigma}^*, \chi) \le 0} \boldsymbol{\sigma}^* : \mathbf{d} \qquad (14) $$

where the subscript $\chi$ underlines the dependence upon the vector of hardening parameters. This function is the support function (the Legendre–Fenchel transform of the indicator function) of the convex set $\{\boldsymbol{\sigma}^* | f(\boldsymbol{\sigma}^*, \chi) \le 0\}$, where $\chi$ acts as a parameter. (The wording *plastic dissipation* should be avoided for this function even more than for the function $\pi$ of classical limit analysis, since the true plastic dissipation should include some contribution from the evolution of the hardening variable $\chi$, absent in the right-hand side of the definition (14).) Then, the inequality

$$ \mathbf{Q}^* \cdot \dot{\mathbf{q}} \le \int_{\Omega} \pi_{\chi}(\mathbf{d}) \, d\Omega \qquad (15) $$
holds for all elements $Q^*$ of $K_{G,\chi}$, all rates $\dot{q}$ of the kinematic vector (unconnected to $Q^*$), and all velocity fields $v$ kinematically admissible with $\dot{q}$, generating a strain rate field $d$.

### 3.4. Update of hardening parameters and geometry

Even when the absence of elasticity warrants applicability of sequential limit analysis, this theory, at the present stage of its presentation, remains incomplete in that it cannot, by definition, say anything about the evolution of the hardening parameters and the geometry, since the expressions $(10)_6$ and $(10)_7$ of $\frac{DX}{Dt}$ and $\dot{x}$ are left out in the set of equations (11) defining it. (Note that it says nothing either about the possible reduction of the full field $\{\chi\}$ of hardening parameters, defined everywhere in $\Omega$, to a finite number of quantities approximately “summarizing” this field; this very important point is the object of active researches, see, e.g., Michel and Suquet [12].)

Sequential limit analysis must therefore be completed with some *ad hoc* procedure aimed at determining such evolutions. Within a kinematic approach, a natural way of doing so is to assume that the trial velocity fields used define an “acceptable approximation” of the true one. These fields can then be used in equation (10)$_7$ to define the evolution of the geometry, in a direct manner, and in equation (10)$_6$ to define the evolution of the hardening parameters, upon calculation of the corresponding strain rate field.

It must be stressed that although this procedure is quite natural, it does not rely on any variational setting, so that the judgment on the quality of the trial velocity fields is an essentially subjective matter and there is no objective control on the quality of the approximate evolutions of the hardening parameters and the geometry obtained. This marks an important difference with the results of limit analysis, the quality of which is guaranteed, or at least controlled, by virtue of the variational setting.

### 3.5. Assessment of sequential limit analysis

The advantages of sequential limit analysis are two-fold:

* it extends the results and methods of classical limit analysis so as to incorporate the effects of strain hardening and/or large displacements and strains. This is important because these two effects have a systematic (for the former) or frequent (for the latter) impact upon plasticity problems of practical significance;

* its use, with a mere transposition, of the well-established procedures of classical limit analysis, makes it easy to use and powerful. The kinematic approach, in particular, is versatile and efficient, thanks to the minimum restrictions put on the possible trial velocity fields.

But the theory also has shortcomings:

* in contrast to classical limit analysis, it is, strictly speaking, applicable only in the complete absence of elasticity. Such a hypothesis is of course never rigorously satisfied. In the absence of general conditions ensuring its soundness, the approximate validity, or complete invalidity of the theory must thus be judged in each specific case, and depends on whether the elastic-plastic coupling plays a minor, or major role in the situation considered;

* its variational setting is incomplete and its procedures and results must be supplemented with *ad hoc* estimates of the evolutions of the hardening parameters and the geometry. The quality of these estimates cannot be guaranteed by variational arguments.

With these advantages and drawbacks, sequential limit analysis provides complete answers to all three questions raised in the Introduction, thus permitting to predict the full evolution in time of *rigid* plastic structures:

* the answer to question (1), pertaining to the mechanical conditions promoting overall plastic flow, lies in the overall yield criterion, depending on the instantaneous distribution of hardening parameters and instantaneous geometrical configuration;

* the answer to question (2), pertaining to the deformation mode corresponding to some load obeying this criterion, lies in the overall flow rule obeying the normality property;

* the answer to question (3), pertaining to the evolution in time of the distribution of hardening parameters and the geometrical configuration, lies in the complementary procedure based on the trial velocity fields of the limit analysis.

In particular, when combined with homogenization of some “representative volume element”, the theory is able to provide a full set of constitutive equations for heterogeneous, *rigid* plastic materials; the specific example of porous plastic solids will now be considered.

## 4. Example: the case of porous ductile materials

In problems of large strain plasticity, elastic strains are often much smaller than plastic ones, so that it is reasonable to disregard elasticity and apply sequential limit analysis. But there are exceptions; probably the most notable of these is
the problem of cavitation in an elastic-plastic medium subjected to some highly hydrostatic loading. In such a situation, a spherical cavity may spontaneously appear and subsequently grow through the plastic flow of the surrounding material, as a result of the relaxation of the elastic energy stored there; see, e.g., Tvergaard et al. [13]. In this case the elastic-plastic coupling is of paramount importance, and sequential limit analysis is clearly inapplicable.

In the sequel we shall consider the case of the quasi-static ductile rupture of metals, which is especially interesting in that (i) it involves all possible complexities: elasticity, plasticity, strain hardening and large displacements and strains; but (ii) the importance of the elastic-plastic coupling is eminently variable and warrants a detailed discussion.

## 4.1. Gurson's [14] model

The most classical constitutive model for materials subjected to ductile rupture is that of Gurson [14]. When discussing it, one must carefully distinguish between the derivation proposed by Gurson himself and alternative, possibly improved ones.

* Gurson's [14] original derivation consisted of two steps.

In a first step, he performed an approximate (classical) limit analysis of a typical representative volume element in a porous plastic medium: namely a hollow sphere, made of an ideal plastic material obeying the von Mises yield criterion and the Prandtl-Reuss associated plastic flow rule, and subjected to some arbitrary loading through conditions of homogeneous boundary strain rate (Hill [15], Mandel [16]). The approximate overall yield criterion he obtained read

$$F(\Sigma, f) = \frac{\Sigma_{eq}^2}{\bar{\sigma}^2} + 2f \cosh\left(\frac{3}{2} \frac{\Sigma_m}{\bar{\sigma}}\right) - 1 - f^2 \le 0 \quad (16)$$

where $\Sigma$ denotes the overall stress tensor, $f$ the porosity (void volume fraction), $\Sigma_{eq}$ the overall von Mises equivalent stress (von Mises norm of the deviator of $\Sigma$), $\Sigma_m$ the overall mean stress (1/3 of the trace of $\Sigma$), and $\bar{\sigma}$ the yield stress of the sound (unvoided) material in uniaxial tension. Also, the overall flow rule followed the normality rule, in agreement with the general property recalled in Subsection 2.4:

$$\mathbf{D}^p = \dot{\Lambda} \frac{\partial F}{\partial \Sigma}(\Sigma, f) , \quad \dot{\Lambda} \ge 0 \quad (17)$$

where $\mathbf{D}^p$ denotes the overall plastic strain rate and $\dot{\Lambda}$ the plastic multiplier.

In a second step, Gurson proposed a complete constitutive model for porous plastic materials by combining several elements. Two of these were the overall yield criterion (16) and the plastic flow rule (17) derived from (classical) limit analysis. Other elements included unconnected, heuristic modelings of the effects of elasticity, isotropic hardening and geometry changes.

* Elasticity was incorporated by supplementing the overall plastic strain rate with some overall elastic strain rate given by some hypoelasticity law (disregarding the influence of porosity).

* Isotropic hardening was introduced by assuming that its effect could be summarized through a single internal variable. Gurson thus heuristically replaced, in his criterion (16) for ideal plastic materials, the constant yield stress $\bar{\sigma}$ by some “average yield stress”, still denoted $\bar{\sigma}$; he assumed this average yield stress to be connected to some “average equivalent strain” $\bar{\epsilon}$ through the relation $\bar{\sigma} = \sigma(\bar{\epsilon})$, where $\sigma(\epsilon)$ denotes the function providing the strain-dependent yield stress of the sound material; and he proposed the following evolution equation for $\bar{\epsilon}$:

$$(1-f)\bar{\sigma} \frac{d\bar{\epsilon}}{dt} = \Sigma : \mathbf{D}^p \quad (18)$$

The basis of this equation was a heuristic identification of the plastic dissipations in the real heterogeneous material, $\Sigma : \mathbf{D}^p$, and in a fictitious “equivalent”, homogeneous, porous material with equivalent strain $\bar{\epsilon}$ and yield stress $\bar{\sigma}$, $(1-f)\bar{\sigma} \frac{d\bar{\epsilon}}{dt}$.

* Finally geometric changes were introduced essentially through the following approximate evolution equation for the porosity:

$$\dot{f} = (1 - f) \operatorname{tr} \mathbf{D}^p \quad (19)$$

resulting from the neglect of elasticity plus the plastic incompressibility of the sound matrix.

* Whether this derivation can be improved, by using the elements expounded in Sections 2 and 3, will now be debated. The first remark is that classical limit analysis can be of no use here. Indeed, its use of a linearized geometrical framework rules out geometry updates and therefore any evolution of the porosity – since this evolution is related to changes in the volume of the voids. Thus it could only justify a “small-displacement, small-strain” version of Gurson's [14] model of no practical interest, since the whole purpose of the model is to depict the gradual degradation and loss of stress-carrying capacity of the material arising from the progressive increase of the porosity.
However, in most problems of ductile rupture under monotone loadings, elasticity plays a minor role, so that sequential limit analysis may be invoked to justify some aspects of the model considered in its rigid-plastic limit.<fnref n="3" /> Such aspects include the overall yield criterion (16), the plastic flow rule (17) obeying the normality property, and the evolution equation (19) of the porosity. (Note, however, that in the absence of elasticity, this evolution equation becomes a rigorous consequence of plastic incompressibility of the matrix, so that the help of sequential limit analysis is not really needed here.)

On the other hand, sequential limit analysis cannot be used to justify Gurson's [14] approach of strain hardening effects in porous ductile materials, and in particular his evolution equation (18) for the internal variable $\bar{\epsilon}$. Indeed, the identification of plastic dissipations it rested upon had nothing to do with, and in fact was incompatible with sequential limit analysis. Gurson never tried, in the spirit of sequential limit analysis, to account for the effect of the detailed heterogeneous distribution of hardening upon the overall criterion and plastic flow rule; nor did he use the trial velocity fields of his limit analysis to specify the evolution in time of this distribution.

## 4.2. Later models

An improved variant of Gurson's [14] modeling of strain hardening effects in porous plastic materials was proposed by Leblond et al. [20]. These authors first evidenced various shortcomings of Gurson's model in the hardenable case – notably its incompatibility with the exact solution to the problem of a hollow sphere made of rigid hardenable material and subjected to hydrostatic tension. They showed that the main drawback in this model was the use of a single internal variable $\bar{\sigma}$ to describe (isotropic) hardening effects. They then proposed the following slight modification of the criterion:

$$F(\Sigma, f) = \frac{\Sigma_{\text{eq}}^2}{\Sigma_1^2} + 2f \cosh\left(\frac{3}{2} \frac{\Sigma_m}{\Sigma_2}\right) - 1 - f^2 \le 0 \quad (20)$$

where $\Sigma_1$ and $\Sigma_2$ denote distinct internal variables connected to the overall yield stresses under purely deviatoric and purely hydrostatic loadings, respectively. Using again the solution to the hollow sphere problem as a basis, they proposed the following expressions for these quantities:

$$\left\{
\begin{aligned}
\Sigma_1 &= \frac{1}{b^3 - a^3} \int_{a^3}^{b^3} \sigma (\langle \epsilon \rangle_{S(r)}) d(r^3) \\
\Sigma_2 &= \frac{1}{\ln(b^3/a^3)} \int_{a^3}^{b^3} \sigma (\langle \epsilon \rangle_{S(r)}) \frac{d(r^3)}{r^3}
\end{aligned}
\right.
\quad (21)$$

where a and b are the inner and outer radii of the sphere (connected through the relation $f = a^3/b^3$), $\sigma(\epsilon)$ denotes the function providing the strain-dependent yield stress of the sound material like above, and $\langle \epsilon \rangle_{S(r)}$ is the average value of the equivalent cumulated strain $\epsilon$ over the spherical surface $S(r)$ of radius $r$. A complex expression of this average value completing the model, which needs not be given here, was also provided, assuming proportional straining at the global scale.

The connections of this work with sequential limit analysis were as follows:

* Leblond et al.'s [20] expressions of the limit loads under purely deviatoric and purely hydrostatic loadings, implicitly contained in equations (20) and (21), were based on homogenization through consideration of the problem of the hollow sphere, considered as a typical representative volume element in the porous material. The derivation of these equations did not, however, resort to sequential limit analysis;

* on the other hand, Leblond et al.'s [20] approximate expression of the average value $\langle \epsilon \rangle_{S(r)}$ was derived, in the absence of an exact solution to the hollow sphere problem for a general, non-hydrostatic loading, by using the trial velocity fields of Gurson's [14] limit analysis. This approach was quite in line with the update of internal parameters completing the procedures of sequential limit analysis (see Subsection 3.4), although this was not mentioned by the authors.

Leblond et al.'s [20] model was recently refined by Lacroix et al. [21] through introduction and radial discretization, at each "macroscopic material point", of an underlying hollow microsphere, and numerical calculation of the integrals in equation (21) based on this discretization. The interval of integration $[a^3, b^3]$ in these integrals was thus divided into $N$

<footnote n="3">³ This limit has been explicitly considered in some works devoted to analytically-based applications of the model: see, e.g., [17–19].</footnote>
sub-intervals, $[a^3 = r_0^3, r_1^3]$, $[r_1^3, r_2^3]$, ..., $[r_{N-1}^3, r_N^3 = b^3]$, and the average values $\langle \epsilon \rangle_i$ of $\epsilon$ in the various intervals $[r_{i-1}^3, r_i^3]$ were calculated (and stored as internal variables) using the evolution equations

$$ \frac{d}{dt} \langle \epsilon \rangle_i = \sqrt{4 \frac{b^6}{\bar{r}_i^6} (D_m^p)^2 + (D_{eq}^p)^2} \quad , \quad \bar{r}_i \equiv \frac{1}{2} (r_{i-1} + r_i) \qquad (22) $$

where $D_m^p$ denotes the overall mean plastic strain rate (1/3 of the trace of $D^p$) and $D_{eq}^p$ the overall equivalent plastic strain rate (von Mises norm of the deviator of $D^p$).<fnref n="4" /> By permitting to relax Leblond et al.'s [20] hypothesis of proportional global straining, this approach led to a much improved description of the distribution of hardening around voids, and thus to more accurate evaluations of the parameters $\Sigma_1$ and $\Sigma_2$.

The relations between Lacroix et al.'s [21] model and sequential limit analysis were the same as for Leblond et al.'s [20] model: the justifications of the criterion (20) and the expressions (21) of the parameters $\Sigma_1$ and $\Sigma_2$ were no different from those in the earlier work and thus made no reference to sequential limit analysis; but the evolution equations (22) of the average values $\langle \epsilon \rangle_i$ were again obtained from Gurson's [14] trial velocity fields, that is using, even if implicitly, the procedure of update of internal parameters completing sequential limit analysis.

Even more recently, Morin et al. [22] used Lacroix et al.'s [21] idea of radial discretization of an underlying hollow microsphere in explicit conjunction – for the first time – with sequential limit analysis, to define extensions of Lacroix et al.'s [21] model to various types of hardening: linear kinematic, nonlinear kinematic (Armstrong and Frederick's [23] model), mixed isotropic/kinematic – a complex model of cyclic plasticity of Chaboche [24] including a memory of the maximum amplitude of previous cycles is even currently considered, with applications to ductile rupture under cyclic loadings in view.<fnref n="5" />

Consider, for instance, the case of mixed isotropic/kinematic hardening, with a local backstress (center of the local reversibility domain) denoted $\alpha$. Morin et al. [22] assumed the distribution of this backstress in each spherical layer $[r_{i-1}^3, r_i^3]$ to be approximately of the form (with classical notations for spherical coordinates)

$$ \alpha^i(r, \theta, \varphi) = A_1^i + A_2^i (-2\mathbf{e}_r \otimes \mathbf{e}_r + \mathbf{e}_\theta \otimes \mathbf{e}_\theta + \mathbf{e}_\varphi \otimes \mathbf{e}_\varphi) \qquad (23) $$

where $A_1^i$ and $A_2^i$ denote a second-order traceless tensor and a scalar, respectively, both uniform within the layer $[r_{i-1}^3, r_i^3]$. With this hypothesis, and using the kinematic approach of sequential limit analysis with Gurson's [14] trial velocity fields, they obtained the following approximate overall yield criterion:

$$ F(\Sigma, f) \equiv \frac{(\Sigma - \mathbf{A}_1)_{\mathrm{eq}}^2}{\Sigma_1^2} + 2f \cosh \left( \frac{3}{2} \frac{\Sigma_m - A_2}{\Sigma_2} \right) - 1 - f^2 \le 0 \qquad (24) $$

where $\mathbf{A}_1$ and $\mathbf{A}_2$ are a second-order traceless tensor and a scalar given by

$$ \left\{ \begin{array}{l} \mathbf{A}_1 = \sum_{i=1}^{N} \mathbf{A}_1^i \frac{r_i^3 - r_{i-1}^3}{b^3} \\ \\ \mathbf{A}_2 = \sum_{i=1}^{N} 2\mathbf{A}_2^i \ln \left( \frac{r_i^3}{r_{i-1}^3} \right) \end{array} \right. \qquad (25) $$

Morin et al. [22] also obtained the evolution equations of the quantities $A_1^i$ and $A_2^i$ characterizing kinematic hardening in the layers, by using the procedure of update of internal parameters completing sequential limit analysis, again with Gurson's [14] trial velocity fields. In the simplest case of linear kinematic hardening, governed at the local scale by the evolution equation

$$ \frac{D\alpha}{Dt} = h\mathbf{d}^p \qquad (26) $$

<footnote n="4">⁴ Note that the implementation of the model did not introduce unreasonable demands with respect to CPU time and storage capacity, since the integrals to be calculated numerically were purely radial, no numerical integration over the spherical variables $\theta$ and $\phi$ being requested.</footnote>

<footnote n="5">⁵ It is also worth mentioning here the very recent paper of Paux et al. [25], in which the authors also used sequential-limit analysis as a help to define a constitutive model, in the slightly different context of plasticity of porous single crystals. These authors' approach rested on a combination of ideas borrowed from the work of Leblond et al. [20], on the one hand, and sequential limit analysis of a hollow sphere, on the other hand, without resorting to Lacroix et al.'s [21] radial discretization of this sphere followed by numerical integration.</footnote>
where $\frac{D}{Dt}$ denotes some objective time-derivative and $h$ a material parameter, the evolution equations of the $A_1^i$ and $A_2^i$ read

$$ \left\{ \begin{aligned} \frac{DA_1^i}{Dt} &= h(D^P)' \\ \dot{A}_2^i &= h \frac{b^3}{\bar{r}_i^3} D_m^P \end{aligned} \right. \qquad (27) $$

where $(D^P)'$ denotes the deviator of $D^P$. These equations complete the model.

Globally, the combination of the concept of the underlying microsphere and the powerful methods of sequential limit analysis revealed very efficient to derive models for porous plastic materials having internal variables obeying various evolution laws. Finite element micromechanical simulations of the behavior of a hollow sphere made of rigid-hardenable materials obeying various hardening laws, and subjected to various cyclic loadings, are currently under progress in order to assess the quality of the predictions of Morin et al.'s [22] model in the absence of elasticity.

However, this model, being explicitly derived from sequential limit analysis, is subjected to the basic restriction of this theory to materials devoid of elasticity. Unfortunately, it so happens that for problems of ductile rupture under cyclic loadings, for which the model was specifically designed, elasticity effects are not negligible. This was noted for the first time by Devaux et al. [26], who identified elasticity and strain hardening as the two factors responsible for the effect of ratcheting of the porosity under cyclic loadings – that is the gradual increase with the number of cycles, under conditions of constant absolute value of the triaxiality, of the porosity averaged over one cycle. The impact of elasticity upon the ratcheting of the porosity was fully confirmed later, notably in the works of Mbiakop et al. [27] and Cheng et al. [28], elasticity being estimated to be responsible for roughly 50% of the effect, the other 50% being attributable to strain hardening.

Thus, in spite of its qualities, Morin et al.'s [22] model is still incomplete in that, by its very construction, it cannot claim to describe the effect of elasticity upon the ratcheting of the porosity under cyclic loadings. Refined homogenization theories under development for nonlinear heterogeneous solids (Lahellec and Suquet [29]) may, in the future, lead to rigorous descriptions of the elastic-plastic coupling in porous plastic solids, and in particular of this effect. In the meantime, however, there is no other solution than to introduce this coupling in a heuristic way in the model considered – be it that of Morin et al. [22] derived from sequential limit analysis, or some other one.

It is worth noting that an interesting first step in this direction has been very recently made by Cheng et al. [28] (independently of Morin et al.'s [22] work and model). Their treatment was based on a hypothesis of “separation of phases”, purely elastic on the one hand, purely plastic (devoid of elasticity) on the other hand. The rates of the porosity in these two types of phases, $\dot{f}^e$ in the elastic ones, $\dot{f}^P$ in the plastic ones, were evaluated separately. This is approximate because in reality presence of plasticity does not mean absence of elasticity; thus $\dot{f}^e$ has no reason to vanish during the plastic phases. However, in spite of the errors introduced by this hypothesis, the model seems promising in that it seems to at least qualitatively capture the mechanism responsible for the influence of elasticity upon the ratcheting of the porosity under cyclic loadings. This mechanism consists of a slight asymmetry introduced by elasticity between the tensile and compressive parts of each cycle: during the tensile part, the volume of the voids first increases a bit reversibly during the initial elastic phase, which facilitates their subsequent irreversible growth through the plastic flow of the surrounding matrix, whereas during the compressive part, this volume first decreases a bit elastically, making the subsequent plastic shrinkage of the voids more difficult.

## 5. Concluding summary

The main purpose of this paper was to examine and precisely state the conditions of applicability of Yang's [1] sequential limit analysis, which stands as a heuristic extension of the classical theory of limit analysis developed by Drucker et al. [7] and Hill [8], incorporating the effects of strain hardening and geometry changes. The treatment was based on a study of the connections between the equations of the full elastic-plastic evolution problem and the simpler ones of both classical and sequential limit analysis. Particular emphasis was put on the role of elasticity.

In a first step, we briefly re-examined the derivation of the results and methods of classical limit analysis, starting from the equations of the elastic-plastic evolution problem and the definition of limit loads. Although not original, this short re-exposition was deemed necessary in view of the use made by sequential limit analysis, our main object of study, of the results and methods of the more classical theory. It was also an occasion to re-state that classical limit analysis, although intrinsically limited to ideal plastic materials considered within a linearized geometrical framework, perfectly tolerates the presence of elasticity. This was deemed useful in view of the fact that Hill's [8] treatment of the subject, based on neglect of elasticity from the start, seems to have in time overshadowed, to some extent, that of Drucker et al. [7], incorporating elasticity.

We then examined if the equations of sequential limit analysis could be deduced from those of the elastic-plastic evolution problem, now considered in the most general context of a hardenable material and a large displacements/large strains framework. Such a deduction proved impossible without any additional hypotheses, showing that the methods of sequential limit analysis are not justified in general. However, in the absence of elasticity, those equations of the evolution problem that do not pertain to the evolution of the hardening parameters and the geometry precisely reduce to those of sequential limit
analysis, which are themselves equivalent to those of classical limit analysis considering the full distribution of hardening parameters and the geometrical configuration as known and fixed. This fully justifies the use of the methods of sequential limit analysis, derived from those of classical limit analysis, for rigid plastic materials, even if subject to strain hardening and geometrical changes. The results derived from these methods must, however, be completed with evolution equations for the hardening parameters and the geometry; these evolution equations can approximately be deduced from the trial velocity fields used in the limit analysis, but without the safeguard of a variational setting guaranteeing the quality of the results. With this additional procedure, the theory is complete and fully answers the three questions raised in the Introduction, for rigid plastic materials.

Of course, in reality, materials always involve some elasticity, but in large-strain plasticity problems, plastic strains are often much larger than elastic ones; following the preceding discussion, for such problems, the use of sequential limit analysis is either totally unjustified or approximately justified, depending on whether the elastic-plastic coupling is important or minor. A discussion of this question was finally presented, in the context of the homogenization-based derivation of models for porous ductile materials. For monotone loadings, elastic-plastic coupling is generally minor, so that elasticity may reasonably be neglected. With this hypothesis, sequential limit analysis is applicable, and justifies Gurson's [14] classical model – except for his independent, purely heuristic modeling of strain hardening effects. For cyclic loadings, the elastic-plastic coupling is important, being responsible for roughly half of the effect of ratcheting of the porosity under conditions of constant absolute value of the triaxiality [26-28]. A very recent model of Morin et al. [22], explicitly based on the use of sequential limit analysis, permits to describe the impact of strain hardening upon this phenomenon accurately, for various isotropic, kinematic or mixed hardening laws of the matrix. This model, however, suffers from the basic restriction of sequential limit analysis to rigid plastic materials, and therefore needs completion in the future by some more or less heuristic modeling of the impact of elasticity upon the ratcheting of the porosity.

## Compliance with ethical standards

The authors declare they have no conflict of interest.

## Acknowledgements

Jean-Baptiste Leblond acknowledges financial support from the "Institut universitaire de France" under grant éOTP L10K018.

## References

[1] W.H. Yang, Large deformation of structures by sequential limit analysis, Int. J. Solids Struct. 30 (1993) 1001-1013.

[2] L. Corradi, N. Panzeri, A triangular finite element for sequential limit analysis of shells, Adv. Eng. Softw. 35 (2004) 633-643.

[3] S.Y. Leu, Analytical and numerical investigation of strain-hardening viscoplastic thick walled-cylinders under internal pressure by using sequential limit analysis, Comput. Methods Appl. Mech. Eng. 196 (2007) 2713-2722.

[4] S.Y. Leu, R.S. Li, Exact solutions of sequential limit analysis of pressurized cylinders with combined hardening based on a generalized Holder inequality: formulation and validation, Int. J. Mech. Sci. 64 (2012) 47-53.

[5] D. Kong, C.M. Martin, B.W. Byrne, Modelling large plastic deformations of cohesive soils using sequential limit analysis, Int. J. Numer. Anal. Methods Geomech. 41 (2017) 1781-1806, https://doi.org/10.1002/nag.2700.

[6] X.P. Yuan, B. Maillot, Y.M. Leroy, Deformation pattern during normal faulting: a sequential limit analysis, J. Geophys. Res., Solid Earth 122 (2017) 1496-1516, https://doi.org/10.1002/2016JB013430.

[7] D.C. Drucker, W. Prager, M.J. Greenberg, Extended limit analysis theorems for continuous media, Q. Appl. Math. 9 (1952) 381-389.

[8] R. Hill, On the state of stress in a plastic-rigid body at the yield point, Philos. Mag. 42 (1951) 868-875.

[9] J. Mandel, Cours de mécanique des milieux continus, Gauthier-Villars, Paris, 1966.

[10] J. Salençon, Calcul à la rupture et analyse limite, Presses de l'École nationale des Ponts et Chaussées, Paris, 1983.

[11] J.J. Moreau, Fonctionnelles Convexes, publication of the Consiglio Nazionale delle Richerche, Roma and the Facolta di Ingegneria di Roma "Torr Vergata", 2003.

[12] J.-C. Michel, M. Suquet, A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations, J. Mech. Phys. Solids 90 (2016) 254-285.

[13] V. Tvergaard, Y. Huang, J.W. Hutchinson, Cavitation instabilities in a power hardening elastic-plastic solid, Eur. J. Mech. A, Solids 11 (1992) 215-231.

[14] A.L. Gurson, Continuum theory of ductile rupture by void nucleation and growth: Part I - yield criteria and flow rules for porous ductile media, ASME J. Eng. Mater. Technol. 99 (1977) 2-15.

[15] R. Hill, The essential structure of constitutive laws for metal composites and polycrystals, J. Mech. Phys. Solids 15 (1967) 79-95.

[16] J. Mandel, Contribution théorique à l'étude de l'écrouissage et des lois d'écoulement plastique, in: Proceedings of the 11th International Congress on Applied Mechanics, Springer, Munich, FRG, 1964, pp. 502-509.

[17] G. Perrin, J.-B. Leblond, Analytical study of a hollow sphere made of porous plastic material and subjected to hydrostatic tension - application to some problems in ductile fracture of metals, Int. J. Plast. 6 (1990) 677-699.

[18] J.-B. Leblond, G. Perrin, A self-consistent approach to coalescence of cavities in inhomogeneously voided ductile solids, J. Mech. Phys. Solids 47 (1999) 1823-1841.

[19] G. Perrin, J.-B. Leblond, Accelerated void growth in porous ductile solids containing two populations of cavities, Int. J. Plast. 16 (2000) 91-120.

[20] J.-B. Leblond, G. Perrin, J. Devaux, An improved Gurson-type model for hardenable ductile metals, Eur. J. Mech. A, Solids 14 (1995) 499-527.

[21] R. Lacroix, J.-B. Leblond, G. Perrin, Numerical study and theoretical modelling of void growth in porous ductile materials subjected to cyclic loadings, Eur. J. Mech. A, Solids 55 (2016) 100-109.

[22] L. Morin, J.-C. Michel, J.-B. Leblond, A Gurson-type layer model for ductile porous solids with isotropic and kinematic hardening, Int. J. Solids Struct. 118-119 (2017) 167-178.
[23] P. Armstrong, C. Frederick, A mathematical representation of the multiaxial Bauschinger effect, Mater. High Temp. 24 (1966) 11–26.

[24] J.-L. Chaboche, Constitutive equations for cyclic plasticity and cyclic viscoplasticity, Int. J. Plast. 5 (1989) 247–302.

[25] J. Paux, R. Brenner, D. Kondo, Plastic yield criterion and hardening of porous single crystals, Int. J. Solids Struct. (2017), https://doi.org/10.1016/
j.ijsolstr.2017.08.030.

[26] J. Devaux, M. Gologanu, J.-B. Leblond, G. Perrin, On continued void growth in ductile metals subjected to cyclic loadings, in: J. Willis (Ed.), Proceedings of the IUTAM Symposium on Nonlinear Analysis of Fracture, Kluwer, Cambridge, GB, 1997, pp. 299–310.

[27] A. Mbiakop, A. Constantinescu, K. Danas, On void shape effects of periodic elasto-plastic materials subjected to cyclic loading, Eur. J. Mech. A, Solids 49 (2014) 481–499.

[28] L. Cheng, K. Danas, A. Constantinescu, D. Kondo, A homogenization model for porous ductile solids under cyclic loads comprising a matrix with isotropic and linear kinematic hardening, Int. J. Solids Struct. 121 (2017) 174–190.

[29] N. Lahellec, P. Suquet, Effective response and field statistics in elasto-plastic and elasto-viscoplastic composites under radial and non-radial loadings, Int. J. Plast. 42 (2013) 1–30.