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

HAL
open science

# On the finite element approximation and the solution by a penalty duality algorithm of an elasto-plastic problem

Bertrand Mercier

► To cite this version:

Bertrand Mercier. On the finite element approximation and the solution by a penalty duality algorithm of an elasto-plastic problem. [Research Report] IRIA-RR-149, IRIA. 1975, pp.35. hal-04716092

HAL Id: hal-04716092

https://inria.hal.science/hal-04716092v1

Submitted on 1 Oct 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
02167 RR149

RIA
laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
78150 - Le Chesnay
Tél. 954.90.20

laboratoire de recherche
en informatique
et automatique

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

On the finite element
approximation
and the solution
by a penalty-duality
algorithm
of an elasto-plastic
problem

B. Mercier

Rapport de Recherche n°149
Décembre 1975
ON THE FINITE ELEMENT APPROXIMATION AND THE SOLUTION
BY A PENALTY-DUALITY ALGORITHM OF AN ELASTO-PLASTIC PROBLEM

B. Mercier

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

RÉSUMÉ :

On étudie l'approximation par éléments finis d'un problème d'élastoplasticité qui provient d'un moindre de Henchy.

On donne pour un exemple la solution du problème approché par un nouvel algorithme de péna-
lisation-dualité, dont on démontre la convergence. La comparaison avec d'autres algorithmes montre
qu'il est bien adapté à ce type de problème.

ABSTRACT :

In this paper, we study a finite element approximation of the elasto-plastic problem which arises
from Hencky's law.

A numerical solution of the approximate problem is given for an example by a new algorithm of
a penalty-duality type, the convergence of which is proved. It is compared with other numerical
methods and seems to be better adapted to this type of problems.
<figure><img src="image_5.png" /></figure>

INTRODUCTION

A variational formulation of Hencky's law for elastic-perfectly-
plastic continua has been given by DUVAUT-LIONS [2].

This results in a mathematically well-posed stress problem.

By duality a displacement problem which we have used in [9] is
defined, in the case of the torsion of cylindrical beams, to obtain some
numerical results, using an approximation by piecewise linear finite
elements.

In this paper, it is shown that, even though the existence of a
solution for the dual problem is not proved, it is possible to prove the
convergence of the approximate stresses.

Then, we try to solve the approximate problem on a two-dimensional
.plane stress) example.

Numerically, two difficulties arise :

- the first one comes from the ill-conditionning of the stiffness
matrix, corresponding to the elasticity operator : if $\lambda_{\min}$ and
$\lambda_{max}$ are its smallest and greatest eigenvalue, $\frac{\lambda_{max}}{\lambda_{min}}$ may be as
large as $10^6$ depending on the shape of the domain.

- the second one comes from a strong non-linearity due to the
plasticity. This causes the functional to be minimized in the dual
problem to be only $C^1$.

On the other hand, we can consider two types of methods :

- algorithms minimizing the functional directly (S.O.R., conjugate gradient). They by-pass the second difficulty because each iteration is not costly.

- algorithms computing the solution of an elastic problem at each step.
  (auxiliary operator, penalty-duality). They circumvent the first
difficulty since the ill-conditionning effect occurs only in a linear
problem.
We introduce here an algorithm of a penalty-duality type which
surmounts the two difficulties, and we demonstrate its convergence.

The presentation is composed of the following four sections :

1. Description of the physical and mathematical problem.

2. Approximation and convergence.

3. Description and convergence of the penalty-duality algorithm.

4. Numerical results and comparison with other methods.

Conclusion.

--------- §---------
# 1. - DESCRIPTION OF THE CONTINUOUS PROBLEM.

## 1.1. Physical problem and Hencky's law.

We consider an elasto-plastic continuum with boundary $Γ$, the interior of which is $Ω ⊂ R^N$ ($N= 2$ or $3$).

The continuum is supposed to be fixed on one part $Γ_U$ of its boundary. It is subjected on the other part $Γ_F$ to a field of external pressures $F_i$, and inside $Ω$ to a field of body forces $f_i$ (Fig. 1).

The continuum is assumed to be isotropic and to satisfy Von Mises plasticity criterion.

We seek for :

- a displacement field $u$ statically admissible :

$$ u \in u_{ad} = \{v_i \in H^1(\Omega), i=1,...,N \mid v_i=0 \text{ on } \Gamma_U\} $$

- a stress tensor $\sigma$

- symmetrical :

$$ \sigma \in Y = \{ \phi_{ij} \in L^2(\Omega), i,j=1,...,N \mid \phi_{ij} = \phi_{ji} \} $$

- statically admissible :

$$ \sigma \in M = \{ \phi_{ij} \in Y \mid \phi_{ij} + f_i = 0 \quad (1)(2) \text{ in } \Omega \,;\, \phi_{ij} n_j = F_i \quad (3) \text{ on } \Gamma_F \} $$

- plastically admissible :

$$ \sigma \in K = \{ \phi_{ij} \in Y \mid \frac{1}{2} \phi_{ij}^D \phi_{ij}^D \le k^2, \text{ a.e. in } \Omega \} \quad (4) $$

- Finally, we suppose that the strain tensor $\epsilon(u)$ and $\sigma$ satisfy to Hencky's law : (see [2])

(1) Einstein convention.

(2) A comma stands for partial derivation.

(3) $n_j$: exterior, normal to $Γ$, outwards directed.

(4) $D$ stands for deviator: $\phi_{ij}^D = \frac{1}{N} \phi_{kk} \delta_{ij}$ where $\delta_{ij}$ is the Kronecker tensor.
$$
(1.1) \qquad (g^{-1}(\sigma), \varphi - \sigma) - (\varepsilon(u), \varphi - \sigma) \ge 0 \quad \forall \varphi \in K.
$$

We recall that $\varepsilon(v)$ is defined by $\varepsilon_{ij}(v) \equiv \frac{1}{2}(v_{i,j} + v_{j,i})$

$g : Y \to Y$ is defined from the compressibility and
shear modulus :

$$
(1.2) \quad (g(e))_{ij} \equiv K_0 e_{kk} \delta_{ij} + 2 \mu e_{ij}^D ,
$$

so that :

$$
(1.3) \quad (g^{-1}(\tau))_{ij} = \frac{1}{N^2 K_0} \tau_{kk} \delta_{ij} + \frac{1}{2\mu} e_{ij}^D .
$$

Furthermore

$$
\forall \tau, \varphi \in Y
$$

$$
(1.4) \qquad (\tau, \varphi) \equiv \int_{\Omega} \tau_{ij} \varphi_{ij} \, dX .
$$

We make the following hypothesis on the elasticity coefficients
(see [2]) :

There exists $\beta, \gamma > 0$ such that :

$$
(1.5) \qquad (g(e),e) \ge \beta|e|^2 \qquad \forall e \in Y
$$

$$
(1.6) \quad (g^{-1}(\tau), \tau) \ge \gamma |\tau|^2 \quad \forall \tau \in Y .
$$

We can see that Hencky's law (1.1) is equivalent to saying that $\sigma$ is
the projection on $K$ of $g(\varepsilon(v))$. We write :

$$
(1.7) \qquad \sigma = f'(\varepsilon(v)) \qquad \text{a.e. in } \Omega
$$

where we set $f': Y \to K$.

$$
(1.8) \quad f'(e) \equiv K_0 e_{kk} \delta_{ij} + \min(2\mu, \frac{k\sqrt{2}}{\sqrt{e_{ij}^D e_{ij}^D}}) e_{ij}^D \quad \text{a.e. in } \Omega.
$$
## 1.2. Primal problem.

Apply Green's formula to (1.1):

$$ (1.9) \qquad (\sigma, \epsilon(v)) = \int_{\Omega} \sigma_{ij} n_j v_i d\Gamma - \int_{\Omega} \sigma_{ij,j} v_i dX $$

we obtain, as $\sigma \in K \cap M$ and $u \in u_{ad}$:

$$ (1.10) \qquad (g^{-1}(\sigma), \varphi - \sigma) - \int_{\Gamma_F} (F_i - \varphi_{ij} n_j) u_i d\sigma + (f_i + \varphi_{ij,j}) u_i dX \ge 0 \\ \forall \varphi \in K \cap M, $$

then choosing $\varphi \in M$ results in the elimination of $u$:

$$ (1.11) \qquad (g^{-1}(\sigma), \varphi - \sigma) \ge 0 \qquad \forall \varphi \in K \cap M $$

thus $\sigma$ is the solution of the **primal problem**:

$$ (P) \qquad \begin{array}{l} \text{Find } \sigma \text{ minimizing} \\ \frac{1}{2}(g^{-1}(\varphi), \varphi) \text{ on } K \cap M \end{array} $$

$K$ and $M$ are closed and convex. The hypothesis (1.6) shows that the bilinear form $(g^{-1}(\tau), \varphi)$ is coercive, which implies:

### THEOREM 1.1

If $K \cap M \neq \emptyset$, then there exists a unique solution $\sigma$ of $(P)$. ■
### 1.3. Dual problem.

As in [2][9] we may introduce a dual problem, the solution of which is the displacement field $u$, if it does exist.

We notice first that $\varphi \in M$ is equivalent to

$$ (1.12) \quad \int_{\Omega} (\varphi_{ij,j} + f_i) v_i dX + \int_{\Gamma_F} (F_i - \varphi_{ij,n_j}) v_i d\Gamma = 0 \quad \forall v \in u_{ad} $$

Let us set $L: V \to R$

$$ (1.13) \quad L(v) = \int_{\Omega} f_i v_i dX + \int_{\Gamma_F} F_i v_i d\Gamma $$

A transformation of (1.12) by Green's formula (1.9) gives :

$$ (1.14) \quad (\varphi, \varepsilon(v)) = L(v) \quad \forall v \in u_{ad} $$

Thus we may write

$$ (1.15) \quad M = \{ \varphi \in Y \mid (\varphi, \varepsilon(v)) = L(v) \quad \forall v \in u_{ad} \} $$

Noticing now that

$$ (1.16) \quad \sup_{v \in u_{ad}} \left[ L(v) - (\varepsilon(v), \varphi) \right] = \begin{cases} 0 & \text{if } \varphi \in M \\ +\infty & \text{otherwise} \end{cases} $$

the primal problem (P) is equivalent to

$$ (1.17) \quad \inf_{\varphi \in K} \sup_{v \in u_{ad}} \left[ \frac{1}{2}(g(\varphi), \varphi) - (\varepsilon(v), \varphi) + L(v) \right] $$

Thus we can define the dual problem by permutation of $\inf$ and $\sup$ over $\varphi$

$$ (1.18) \quad \sup_{v \in u_{ad}} \inf_{\varphi \in K} \left\{ \frac{1}{2} (g^{-1}(\varphi), \varphi) - (\varepsilon(v), \varphi) + L(v) \right\} $$

taking into account the definition (1.8) of $f'$, we see that for each $v$ the optimal $\varphi$ minimizing $\frac{1}{2}(g^{-1}(\varphi), \varphi) - (\varepsilon(v), \varphi)$ on $K$ is :

$$ (1.19) \quad \varphi = f'(\varepsilon(v)) \quad \text{a.e. in } \Omega. $$
the dual problem is then :

( ρ* )

$$
\left\{
\begin{array}{l}
\text{Find } u \in \mathcal{U}_{\text{ad}} \text{ minimizing} \\
\text{J}(v) = f(\varepsilon(v)) - L(v) \\
\text{on } \mathcal{U}_{\text{ad}}.
\end{array}
\right.
$$

where $f:Y \rightarrow R$ is defined as follows :

$$
(1.20) \quad f(e) \equiv K_o (e_{kk})^2 + \begin{cases}
\mu e_{ij}^D e_{ij}^D & \text{if } e_{ij}^D e_{ij}^D \le \frac{k^2}{2\mu^2} \\
k \left( \left(2e_{ij}^D e_{ij}^D\right)^{\frac{1}{2}} - \frac{k}{2\mu} \right) & \text{otherwise}
\end{cases}
$$

In fact, $f$ is $C^1$-Gâteaux-differentiable and its gradient is $f'$
defined in (1.8). Thus, $(\rho^*)$ is equivalent to the variational equation :

$$
(1.21) \quad (f'(\epsilon(v)), \epsilon(w)) - L(w) = 0 \quad \forall w \in \mathcal{U}_{\text{ad}}.
$$

1.4. Existence.

. The functional J to be minimized is infinite when v is infinite
only in (W^1,1(Ω))^N when the following hypothesis is satisfied : (see [2])

$$
(1.22) \quad \exists \sigma_1 K \mathbb{P} M \text{ such that } \frac{1}{2} \sigma_{ij}^{1D} \sigma_{ij}^{1D} \le k_1^2 < k^2 \text{ a.e.}
$$

Since $W^{1,1}(\Omega)$ is not reflexive, the existence of a solution for $(\rho^*)$
may be expected in its bidual $(W^{1,1}(\Omega))$".

Note that the hypothesis (1.22) implies in particular $K \cap M \neq \emptyset$,
and it is satisfied only if the loads $f_i$ and $F_i$ are sufficiently
small.
## 1.5. Extremality solution.

If the problem $(\rho^*)$ has a solution $u$, it is related to the
solution $\sigma$ of the primal problem $(\mathcal{P})$ by the extremality relation (1.19)
(see [3]):

$$
(1.23) \qquad \sigma = f'(\epsilon(u)) \quad \text{a.e. in } \Omega
$$

This is precisely Hencky's law, relating displacements and stresses.

In the continuum there will be elastic and plastic zones, separated
by a free boundary.

The elastic zones will be characterized by :

$$
\frac{1}{2} \sigma_{ij}^{D} \sigma_{ij}^{D} < k^2 \Leftrightarrow \epsilon_{ij}^{D}(v) \epsilon_{ij}^{D}(v) < \frac{k^2}{2\mu^2} .
$$

The plastic zones will be characterized by :

$$
\frac{1}{2} \sigma_{ij}^{D} \sigma_{ij}^{D} = k^2 \Leftrightarrow \epsilon_{ij}^{D}(v) \epsilon_{ij}^{D}(v) \ge \frac{k^2}{2\mu^2} .
$$

Let $\Omega_e$ and $\Omega_p$ be these zones, they are unique because $\sigma$ is unique.

1.6. Uniqueness.

If there exist two solutions $u$ and $\bar{u}$, they must both satisfy
(1.23). This shows that $\varepsilon(u) = \varepsilon(\bar{u})$ on $\Omega_e$.

All we can say is that if $\Omega_e$ is connected and $\partial\Omega_e \cap \Gamma_U$ has a non zero measure, then $u = \bar{u}$ on $\Omega_e$, if it is smooth.

On $\Omega_p$ we know only that $\epsilon_{ij}(\bar{u}) = \lambda \epsilon_{ij}^D(u) + \frac{1}{N} \epsilon_{kk}(u) \delta_{ij}$
where $\lambda \in L^\infty(\Omega_p)$.
## 2. APPROXIMATION.

### 2.1. Approximating spaces.

We suppose, for the sake of simplicity, that $N=2$ and that $\Omega$ is polyhedric.

We have seen that only the primal problem is mathematically well-posed. But on the other hand, the dual problem seems easier to solve numerically because there are no constraints.

In order to approximate both the primal and the dual problems, we approximate $Y$ by $Y_h \subset Y$ and $V$ by $V_h \subset U_{ad}$.

We suppose that:

$$ (2.1) \qquad \forall v \in V_h \quad \text{then} \quad \varepsilon(v) \in Y_h. $$

In the following, we shall choose a "regular" family of triangulations of $\Omega$ of $\mathcal{L}_h$ (see [1]), such that the size of the triangles in $\mathcal{L}_h$ be smaller than $h$.

We take:

$$ (2.2) \qquad V_h = \{ v_i \in C^0(\Omega), i=1,2 \mid v_i \mid_{T_i} \in P_i^{(1)} \cup \mathcal{L}_h \mid v_i = 0 \text{ on } \Gamma_i \} $$

$$ (2.3) \qquad Y_h = \{ \tau \in Y \mid \tau_{ij} \in P_o \text{ on each } T_i \in \mathcal{L}_h \} $$

We check immediately that $Y_h$ and $V_h$ satisfy (2.1).

We know that $V_h$ is an internal approximation of $U_{ad}$ that is (see [1]):

$$ (2.4) \qquad \forall v \in U_{ad} \quad \text{we may find} \quad v_h \in V_h \quad \text{such that} \quad v_h \to v \quad \text{strongly in } U_{ad}. $$

On the other hand, let us define $r_h$ as the projector $Y \to Y_h$. That is, we have, for any $\varphi \in Y$, $r_h \varphi \in Y_h$ and

$$ (2.5) \qquad (\varphi - r_h \varphi, z_h)_Y = 0 \quad \forall z_h \in Y_h. $$

<footnote>(1) $P_k$ is the set of polynomials of degree $\le k$.</footnote>
Obviously

$$
(2.6) \qquad r_h \varphi \to \varphi \text{ in } Y \text{ strongly when } h \to 0.
$$

2.2. <ins>Approximation of the primal problem.</ins>

We approximate M by

$$
(2.7) \quad M_h = \{ \varphi \in Y_h \mid (\varphi, \varepsilon(v)) = L(v) \quad \forall v \in V_h \}
$$

and
$$
K \cap M \quad \text{by} \quad K \cap M_h.
$$

We shall prove the following properties for K ∩ M_h :

LEMMA 2.1. Let $Y_h$ and $V_h$ be defined by (2.2) and (2.3) and $M_h$ by (2.7)
Then :

$$
(2.8) \quad K \cap M_h \text{ is a closed convex subset of } Y_h
$$

(i) $\forall \varphi \in K \cap M$ we can find $\varphi_h \in K \cap M_h$ with $\varphi_h \to \varphi$ in $Y$ **strongly**.

(ii) If $\tau_h \in K \cap M_h$, $\tau_h \to \tau$ in $Y$ **weakly**, then $\tau \in K \cap M$.

Proof.

Since $M_h$ is an affine linear variety in a finite dimensional space,
then obviously (2.8) holds. Let us prove (2.9) now.

(i) Let $\varphi \in K \cap M$, we choose $\varphi_h = r_h \varphi$.

Apply (2.5) with $z_h = \varepsilon(v_h)$ for any $v_h \in V_h$.

Then :

$$
(2.10) \qquad (r_h \varphi, \epsilon(v_h)) = (\varphi, \epsilon(v_h)) .
$$

And $\varphi \in M$ implies now :

$$
(2.11) \qquad (r_h \varphi, \varepsilon(v_h)) = L(v_h)
$$

which proves that $r_h \varphi \in M_h$.
On the other hand, T being any triangle of $\ell_h$, let us apply (2.5)
with

$$
z_h = \begin{cases} r_h \varphi^{\mathrm{D}} & \text{on } T \\ 0 & \text{otherwise.} \end{cases}
$$

we deduce that

$$
\int_T \varphi_{ij}^D \cdot r_h \varphi_{ij}^D \, dX = \int_T r_h \varphi_{ij}^D \cdot \iota_h \varphi_{ij}^D \, dX \quad \forall T \in \mathcal{L}_h
$$

$r_h \varphi$ being constant on each triangle, we may write :

$$
(2.12) \quad |T| |r_h \varphi|^2 = r_h \varphi_{ij}^D \cdot \int_T \varphi_{ij}^D \le |r_h \varphi|^2 \cdot \int_T |\varphi|^2 .
$$

applying Schwartz' inequality.

Noticing now that $\varphi \in K \Rightarrow |\varphi|^D = \sqrt{\varphi_{ij}^D \varphi_{ij}^D} \le k\sqrt{2}$ (2.12) gives

$$
(2.13) \qquad |r_h \varphi|^D \le k \sqrt{2}
$$

this shows that $r_h \varphi \in K$

■ (i)

(ii) If $\tau_h \to \tau$ in $Y$ weakly, then for any $z \in Y$ with $z_h \to z$ in $Y$ strongly :

$$
(2.14) \quad (z_h, \tau_h) \to (z, \tau) \quad \text{when} \quad h \to 0.
$$

For any $v \in u_{ad}$ it is possible to find $v_h \to v$ in $u_{ad}$ strongly.

Thus $\varepsilon(v_h) \to \varepsilon(v)$ in $Y$ strongly.

Then, applying (2.14) with $\tau_h = \varepsilon(v_h)$ and $\tau = \varepsilon(v)$, we get

$$
(2.15) \quad (\tau_h, \varepsilon(v_h)) \to (\tau, \varepsilon(v)).
$$
Moreover, since $\tau_h \in M_h$, $(\tau_h, \varepsilon(v_h)) = L(v_h) \rightarrow L(v)$ when $h \rightarrow 0$

Then, in the limit :

$$ (\tau, \varepsilon(v)) = L(v) \quad \forall v \in u_{ad} $$

which shows that $\tau \in M$.

On the other hand $K$ being closed and convex, it is weakly closed, and $\tau \in K$.

Then, we approximate the primal problem (P) by

$$ (\mathcal{P}_h) \qquad \begin{array}{l} \text{Find } \sigma_h \text{ minimizing} \\ \frac{1}{2} (g^{-1}(\varphi_h), \varphi_h) \text{ on } K \cap M_h \end{array} $$

(2.9)(i) implying that $K \cap M_h \neq \emptyset$, there exists a unique solution $\sigma_h$ to $(\mathcal{P}_h)$.

In fact, $(\mathcal{P}_h)$ is equivalent to the variational inequality :

$$ (2.16) \qquad (g^{-1}(\sigma_h), \varphi - \sigma_h) \ge 0 \quad \forall \varphi \in K \cap M_h. $$

With the results of Lemma 2.1., we satisfy the hypothesis of a theorem of GLOWINSKI-LIONS-TREMOLIERES [6] for the convergence of an internal approximation. Thus, we can state the following theorem :

THEOREM 2.1.

If $Y_h$ and $V_h$ are chosen as in (2.2)(2.3) and $M_h$ as in (2.7) then the solution $\sigma_h$ of the approximate problem $(\mathcal{P}_h)$ converges strongly in $Y$ to $\sigma$ solution of (P).
### 2.3. Approximation of the dual problem.

We notice first that for $\varphi \in Y_h$

$$ (2.17) \qquad \sup_{v \in V_h} \left[ L(v) - (\varphi, \varepsilon(v)) \right] = \begin{cases} 0 & \text{if } \varphi \in M_h \\ +\infty & \text{otherwise} \end{cases} $$

Thus $(\rho_h^*)$ may be written :

$$ (2.18) \qquad \inf_{\varphi \in Y_h^*} \sup_{v \in V_h} \left[ \frac{1}{2}(g^{-1}(\varphi, \varphi) - (\varphi, \varepsilon(v)) + L(v)) \right] $$

Once again, we can define a dual problem by permutation of

Inf and Sup :

$$ (2.19) \qquad \sup_{v \in V_h} \left\{ \inf_{\varphi \in Y_h^*} \sup_{v \in V_h} \left[ \frac{1}{2}(g^{-1}(\varphi, \varphi) - (\varphi, \varepsilon(v))) + L(v) \right] \right\}. $$

Since (2.1) is satisfied, the evaluation of $\inf_{\varphi \in Y_h^*}$ gives the same result as in (1.19). Thus the dual problem of $(\rho_h^*)$ may be written:

$$ \begin{array}{ll} \text{(}\rho_h^*\text{)} & \text{Find } u_h \in V_h \text{ minimizing} \\ & J(v) = f(\varepsilon(v)) - L(v) \text{ on } V_h \end{array} $$

This in turn can be interpreted as the approximate problem of $(\rho_h^*)$ corresponding to the internal approximation of $V$ by $V_h$.

Similarly, the extremality relation (1.23) will still hold for $\sigma_h$ and $u_h$.

$V_h$ being finite dimensional, we can see that $J(v_h) \to +\infty$ when $\|v_h\| \to \infty$. Then, there exists at least a solution $u_h$ to $(\rho_h^*)$.
## 2.4. Uniqueness in the approximate problem.

We have a uniqueness result analogous to the one of § 1.6.
That is, the approximate elastic and plastic zones $\Omega_e^h$ and $\Omega_p^h$ are unique.

If we suppose in addition (as in §1.6), that

$$ (2.20) \qquad \left\{ \begin{array}{l} \Omega_e^h \text{ is connected} \\ \operatorname{mes}(\partial \Omega_e^h \cap \Gamma_U) \neq 0 \end{array} \right. $$

(that is there is at least a triangle of $\Omega_e^h$ having one side on $\Gamma_U$), then $u_h$ is unique on $\Omega_e^h$.

Note that as in § 1.6 (2.20) can occur only if the external loads $F_i$ and $f_i$ are sufficiently small.

# 3: DESCRIPTION AND CONVERGENCE OF A PENALTY-DUALITY ALGORITHM TO SOLVE THE APPROXIMATE PROBLEM.

## 3.1. Idea.

The dual problem ($\Gamma_h^p$) is equivalent to the unconstrained minimization of a function of v, when only v and $\epsilon(v)$ appear.

As in [5], we introduce an artificial constraint $e = \epsilon(v)$, and we minimize the new functional in e and v subject to this constraint.

Then, as in [5], to solve the new problem, we penalize and dualize simultaneously the artificial constraint $e = \epsilon(v)$ following an idea of HESTENES [7].
### 3.2. An equivalent problem.

Following the idea of §3.1., we notice that the minimization problem $(\mathcal{P}_h^*)$ is equivalent to :

$$ (3.1) \qquad \begin{array}{ll} \underset{e \in Y_h, v \in V_h}{\operatorname{Inf}} & \left[f(e) - L(v)\right] \\ & e = \sigma(v) \end{array} $$

Let $\mathcal{L}_\alpha$ be the lagrangian defined on $(V_h \times_{\mathcal{Y}_h} \mathcal{Y}_h) \times_{\mathcal{Y}_h}$ for $\alpha > 0$:

$$ (3.2) \qquad \mathcal{L}_{\alpha}(v, e; \sigma) \equiv f(e) - L(v) + (\sigma, \varepsilon(v) - e) + \frac{1}{2\alpha}(g(e - \varepsilon(v)), e - \varepsilon(v)). $$

We introduce the problem :

$(\mathcal{P}_\alpha)$

Find $(u, e; \sigma)$ saddle-point of $\mathcal{L}_\alpha$ on $(V_h \times_{\mathcal{Y}_h} \mathcal{Y}_h) \times_{\mathcal{Y}_h}$.

We can prove the following result :

#### PROPOSITION 3.1.

If $u_h$ and $\sigma_h$ are solutions of $(\mathcal{P}_h^*)$ and $(\mathcal{P}_h)$, then $(u_h, \varepsilon(u_h); \sigma_h)$ is a solution of $(\mathcal{P}_\alpha)$.

Reciprocally if $(v, e; \sigma)$ is solution of $(\mathcal{P}_\alpha)$, then $v$ is a solution of $(\mathcal{P}_h^*)$, while $e = \varepsilon(v)$ and $\sigma = \sigma_h$ solution of $(\mathcal{P}_h)$.

#### Proof.

As $\mathcal{L}_\alpha$ is $C^1$-Gateaux-differentiable, $(\mathcal{P}_\alpha)$ is equivalent to the system of Kuhn and Tucker conditions :
$$
(3.3) \qquad \left( \frac{\partial \ell_{\alpha}}{\partial v} , w \right) = 0 \qquad \forall w \in V_h
$$

$$
(3.4) \qquad \left( \frac{\partial \ell_{\alpha}}{\partial \epsilon} , \varphi \right) = 0 \qquad \forall \varphi \in Y_h
$$

$$
(3.5) \qquad \left( \frac{\partial \ell_{\alpha}}{\partial \sigma}, \tau \right) = 0 \qquad \forall \tau \in Y_h,
$$

where

$$
(3.6) \quad \left( \frac{\partial \ell_{\alpha}}{\partial v}, w \right) = \left( \sigma, \varepsilon(w) \right) - L(w) + \frac{1}{\alpha} \left( g(\varepsilon(v)-e), \varepsilon(w) \right)
$$

$$
(3.7) \qquad \frac{\partial \ell_{\alpha}}{\partial e} = f'(e) - \sigma + \frac{1}{\alpha} g(e - \varepsilon(v)) \quad \text{on each } T \in \mathcal{T}_{h}
$$

$$
(3.8) \qquad \frac{\partial \ell_{\alpha}}{\partial \sigma} = \varepsilon(v) - e \qquad \text{on each } T \in \mathcal{T}_{h}
$$

We can see that if $u_h$ and $\sigma_h$ are solutions of $(\rho_h^*)$ and $(\rho_h)$,
the extremality relation (1.23) gives us :

$$
\sigma_h = f'(\varepsilon(u_h)).
$$

Then obviously the conditions (3.3)(3.4)(3.5) hold.

Reciprocally if (u,e; σ) is solution of (ρ) then (3.5) implies :

$$
(3.9) \qquad e = \varepsilon(v).
$$

Inserting in (3.4), we get

$$
(3.10) \qquad f'(e) = \sigma
$$

so that (3.3) gives (f'(ε(v)), ε(w)) = L(w) ∀ w ∈ V_h
thus v is solution of (P_h^*).

3.3. Algorithm.

In order to solve (P_α) we use the following algorithm which is a
variation of Uzawa's algorithm.
$e^o$ and $\sigma^o$ are given,

by induction, $e^n$ and $\sigma^n$ being given, compute

• $v^{n+1}$ solution of the variational equation

$$ (3.11) \qquad (g(\varepsilon(v^{n+1})), \varepsilon(w)) = (g(e^n) - \alpha\sigma^n, \varepsilon(w)) + \alpha L(w) \quad \forall w \in V_h $$

• $e^{n+1}$ solution of the non-linear equation :

$$ (3.12) \qquad f'(e^{n+1}) + \frac{1}{\alpha} g(e^{n+1}) = \frac{1}{\alpha} g(\varepsilon(v^{n+1})) + \sigma^n $$

• and $\sigma^{n+1}$ :

$$ (3.13) \qquad \sigma^{n+1} = \sigma^n + \rho g(e^{n+1} - \varepsilon(v^{n+1})) \\ \qquad \qquad \qquad \qquad \qquad \text{(penalty-duality algorithm)} $$

we can show the following theorem :

**THEOREM 3.1.**

Let $Y_h$ and $V_h$ be chosen as in (2.2) (2.3).

If $\rho = \frac{1}{\alpha}$ then the penalty-duality algorithm (3.17) (3.18) (3.19) converges in the sense that $\sigma^{n+1} \to \sigma_h$, a solution of $\rho_h$.

If hypothesis (2.20) holds, then $v^{n+1} \to u_h$ on $\Omega_e^h$ and $u_h$ is a solution of $(\rho_h^*)$.

**Proof** : a) **preliminaries** : we need the notion of "g-projection".

Let us consider first the following two subspaces of $Y_h$ :

$$ M_o = \{ \varphi \in Y_h \mid (\varphi, \varepsilon(v)) = 0 \quad \forall v \in V_h \} $$

$$ M_l = \{ \varphi \in Y_h \mid \exists v \in V_h \text{ such that } \varphi = g(\varepsilon(v)) \} $$
$M_o$ and $M_l$ have the following properties :

1°) $M_o$ and $M_l$ are "g-orthogonal", that is

$$ (3.14) \qquad (\varphi_o, \varphi_l)_g \equiv (\varphi_o, g^{-1}(\varphi_l)) = 0 \quad \forall \varphi_o \in M_o \text{ and } \varphi_l \in M_l. $$
(This follows from the definition of $M_o$ and $M_l$).

$$ 2°) \quad M_o \oplus M_l = Y_h. $$

In fact : $\forall \varphi \in Y_h$ we have $\varphi = \varphi_o + \varphi_l$ where

$$ (3.15) \qquad \left\{ \begin{array}{ll} \varphi_o = \varphi - g(\varepsilon(v)) & \in M_o \\ \varphi_l = g(\varepsilon(v)) & \in M_l \end{array} \right. $$

if $v \in V_h$ is the unique solution of the variational problem

$$ (3.16) \qquad (g(\varepsilon(v)), \varepsilon(w)) = (\varphi_g \llcorner (w)) \quad \forall w \in V_h $$

and (3.14) implies :

$$ (3.17) \qquad |\varphi|^2_g = |\varphi_o|^2_g + |\varphi_l|^2_g $$

Then $\varphi_o$ is the "g-projection" of $\varphi$ on $M_o$ and $\varphi_l$ on $M_l$

$$ (3.18) \qquad \left\{ \begin{array}{l} \varphi_o = P_{M_o} (\varphi) \\ \varphi_l = P_{M_l}^o (\varphi) \end{array} \right. $$

where $P_{M_o}$ and $P_{M_l}$ are "g-projectors", that is satisfy all the properties of a projector but for the scalar product $(, )_g$ defined in (3.14).
b) estimates : subtracting (3.13) from (3.12) we deduce that if $\rho = \frac{1}{\alpha}$

$$
(3.19) \qquad f'(e^{n+1}) = \sigma^{n+1} \qquad (\forall n \ge 0)
$$

We use this fact subsequently :

Let $(v,e; \sigma_h)$ be solution of $(\rho_d)$ :

Subtracting (3.4) from (3.12), and (3.5) from (3.11), we get

$$
(3.20) \qquad (\sigma^{n+1} - \sigma_h) + \frac{1}{\alpha} g(e^{n+1}-e) = \frac{1}{\alpha} g(\varepsilon(v^{n+1}-v)) + (\sigma^n - \sigma_h)
$$

and

$$
(3.21) \qquad \left( g(\varepsilon(v^{n+1}-v)) - g(e^n-e) + \alpha(\sigma^n - \sigma_h) \right), \varepsilon(w) = 0 \quad \forall w \in V_h
$$

(3.21) along with (3.15)(3.16) shows that :

$$
(3.22) \qquad g(\varepsilon(v^{n+1})) = P_{M_1} (g(e^n-e)) - \alpha P_{M_1} (\sigma^n - \sigma_h) .
$$

Thus (3.20) can be written ;

$$
(3.23) \qquad (\sigma^{n+1} - \sigma_h) + \frac{1}{\alpha} g(e^{n+1}-e) = \frac{1}{\alpha} P_{M_1}(g(e^n-e)) - P_{M_0}(\sigma^n - \sigma_h)
$$

Taking the *g*-scalar product of (3.24) with itself gives :

$$
(3.24) \qquad \begin{aligned}[t]
|\sigma^{n+1} - \sigma_h|^2_g + \frac{2}{\alpha}(\sigma^{n+1} - \sigma_h, e^{n+1}-e) + \frac{1}{\alpha^2}|e^{n+1}-e|^2_g &\le \\
&\le |\sigma^n - \sigma_h|^2_g + \frac{1}{\alpha^2}|e^n-e|^2_g
\end{aligned}
.
$$

(3.19) and the convexity of $f$ implies that $(\sigma^{n+1}-\sigma_h, e^{n+1}-e) \ge 0$
thus the sequence $\{\tau^n\}$ where

$$
(3.25) \qquad \tau^n = |\sigma^n - \sigma_h|^2_g + \frac{1}{\alpha^2} |e^n - e|^2_g
$$

is a decreasing sequence. Since $\tau^n \ge 0$, this sequence has a limit
as $n \to \infty$ and $(\tau^{n+1} - \tau^n) \to 0$.
Returning now to (3.24), we see that :

$$
(3.16) \qquad ( \sigma^{n+1} - \sigma_h, e^{n+1} - e ) \to 0 .
$$

Notice that (3.10)(3.19) is Hencky's law (1.7) for ($\sigma^{n+1}$, $e^{n+1}$) on the one hand, and for ($\sigma_h$, $e$) on the other hand. Then (1.1) implies :

$$
(3.27) \qquad (g^{-1}(\sigma^{n+1}) - e^{n+1}, \tau - \sigma^{n+1}) \ge 0 \qquad \forall \tau \in K
$$

and

$$
(3.28) \qquad (g^{-1}(\sigma_h) - e, \tau - \sigma_h) \ge 0 \qquad \forall \tau \in K
$$

$$
\text{Since } e^{n+1}-e = e^{n+1}-g^{-1}(\sigma^{n+1})+g^{-1}(\sigma^{n+1})-g^{-1}(\sigma_h)+g^{-1}(\sigma_h)-e
$$

we deduce from (3.26)(3.27)(3.28) that

$$
|\sigma^{n+1} - \sigma_h|^2_g = (\sigma^{n+1} - \sigma_h, g^{-1}(\sigma^{n+1} - \sigma_h)) \to 0.
$$

Thus $\sigma^{n+1} \to \sigma_h$ the solution of ($P_h^*$).

With an argument similar to the one used in §1.6, we see that $|e^{n+1}-e|\to 0$
at least on each triangle $T \in \Omega_e^h$ (approximate elastic zone), and
on the other hand (3.20) now implies that $(e^{n+1}-\epsilon(v^{n+1})) \to 0$ on $\Omega_e^h$.

The conclusion is that if (2.20) is satisfied, then
$v^{n+1} \to v = u_h$ a solution of ($P_h^*$) which is unique on $\Omega_e^h$.

3.4. Notes and comments.

In fact, it is proved in an analogous case that the algorithm
converges for $0 < \rho < \frac{2}{\alpha}$. (see [4]).

For the programming considerations, the penalty-duality algorithm
consists of solving at each step :

* an elasticity problem, with a fixed matrix, and a variable second member,

* a non-linear problem which can be solved explicitly triangle by triangle by the following formula :
$$ e_{ij}^{n+1} = \frac{\alpha}{2\mu} \max \left( \frac{1}{1+\alpha'}, \frac{1-\sqrt{\frac{k\sqrt{2}}{\psi_{ij}^D \psi_{ij}^D}}}{\psi_{ij}^D} \right) \psi_{ij}^D \frac{1}{N^{2K_o}} \frac{\alpha}{1+\alpha} \psi_{kk} \delta_{ij} $$

where $\psi = \frac{1}{\alpha} g(\varepsilon(v^{n+1})) + \sigma^n$

* an iteration on $\sigma^{n+1}$ (cf. (3.19)) which has the same property.

We can see that at each step the real work will be to solve a linear system always with the same matrix. Then we can compute the Cholesky factorization of the matrix once and for all.
## 4. - NUMERICAL RESULTS.

### 4.1. Physical problem solved.

All the numerical experiments have been performed on the following two-dimensional problems.

Let us consider a beam of length $l$, and thickness $2a$, which is cantilevered at one end. At the other end, the beam is subjected to a load $F$ directed downwards (Fig. 2). We consider this as a plane stress problem.

In the following we call $U$ the vertical displacement at the free end ($x_2 = l$).

<figure><img src="image_6.png" /><figcaption>Fig. 2 Cantilevered beam considered</figcaption></figure>

This is a particular case of the problem considered in 1.1.
with

$$ \Omega = ]-a,a[ \times ]0,l[ \subset \mathbb{R}^2 $$

$$ \Gamma_U = \Gamma_1 $$

$$ f_i = 0 \quad \text{on } \Omega $$

$$ \Gamma_F = \Gamma_2 \cup \Gamma_3 \cup \Gamma_4 $$

$$ F_i = 0 \quad \text{on } \Gamma_2 \text{ and } \Gamma_3 $$

$$ F_i = \begin{cases} F & \text{on } \Gamma_4 \\ 0 & \end{cases} $$
## 4.2. Numerical data.

We take

$$a = .2 \qquad l = 20$$

$$\{ E = 15 \text{ (elasticity modulus)} $$

$$ \gamma = 0.33 \text{ (Poisson's coefficient)} $$

so that $K_o = E/2(1-\gamma)$ and $\mu = E/2(1+\gamma)$

$$ \frac{k^2}{2\mu} = 6.10^{-6} \text{ arbitrarily} $$

Finally we take $|F_o| = 0.01$

and $F = s \cdot F_o$ where s is a scale factor.

## 4.3. Discretization.

Since $\Omega$ is rectangular we choose a regular mesh of finite elements (see Fig. 4) with two discretization steps :

$h_1 = \frac{4}{n_1}$ in the $x_1$ direction

and $h_2 = \frac{20}{n_2}$ in the $x_2$ direction.

Thus there are $2 \cdot n_1 \cdot n_2$ finite elements and $2(n_1+1)n_2$ independent variables in $V_h$ ($6 n_1 n_2$ in $Y_h$).

## 4.4. Solution of the elastic problem.

Since at each step of the penalty-duality algorithm we have to solve an elasticity problem, we tried first to solve it alone.

The elasticity problem is the particular case of our problem where $k = +\infty$, so that $(P_h^*)$ is the minimization of a quadratic functional.

We compare several methods, that we present in increasing order of difficulty, from the point of view of programming.
- successive over-relaxation (S.O.R.)

with $h_1 = h_2 = \frac{1}{3}$ (24 independent variables), the algorithm failed to converge in 400 iterations. This seems to be due to a strong ill-conditioning in the "stiffness" matrix :

if $\lambda_{\min}$ and $\lambda_{\max}$ are its smallest and greatest eigenvalue, $\lambda_{\max}^*$ may be of order $10^6$, depending on the shape of $\Omega$ !

- steepest descent.

under the same conditions the pure gradient method, with
optimal stepsize gave only slight improvements.

- conjugate gradient.

on the other hand the method of conjugate gradient worked very
well, and in each case we obtained convergence in a number of
iterations less than the number of variables. However double
precision was required in the computations (see Fig. 3 and 5).

- direct method.

since the functional to be minimized in $(\rho_h^*)$ is quadratic in
the elastic case, the numerical problem is equivalent to solving
a linear system with a stiffness matrix, which is symmetric posi-
tive definite, and moreover sparse. Its band width is $2(1+\min(n_p n_2))$
when the unknowns are judiciously ordered.

We employed a Cholesky band factorization with good results.

4.5. Solution of the elasto-plastic problem by non-linear conjugate gradient

In § 4.4. we saw that among the iterative methods with which we
tried to solve the elastic problem, only the conjugate gradient gave satis-
faction.

This method in turn, can be adapted easily to the elasto-plastic case to solve the dual problem $(\rho_h^*)$ directly. However a few problems arise :

- frequency of reinitialization : (see [6])

that is : if $p$ is the number of variables, there is no justification to look for more than $p$ conjugate directions. One should
restart in the gradient direction. However practically it
is better to reinitialize every 2p or 3p iteration.

- minimization in one variable.

J is only $C^1$ Gateaux-differentiable, but we can compute
locally its Hessian J''. So, if we start from $v^0$ and try
to minimize J in the direction $w$, we can proceed by
induction as follows : (Newton algorithm)

$$
v^{m+1} = v^m + \rho_m w
$$

where

$$
\rho_m = \frac{(J'(v^m), w)}{(J''(v^m), w, w)}.
$$

In fact, with the choice made for $v_h$, $J''(v)$ is constant
on each triangle and has two determinations only where
$\epsilon_{ij}^D(v) \epsilon_{ij}^D(v) = \frac{k^2}{2\mu^2}$ : the "elastic" one and the "plastic" one.

We take by convention the elastic one, and try the algorithm.

If it begins to diverge, we correct it by dichotomy.

Practically we took a stopping test for this algorithm :

$$
\rho_m \le \epsilon = 10^{-10} .
$$

It converges in less than 4 iterations, and most of the time,
the second iteration checked that the first iteration was
enough.

Finally (see results on Fig. 3 and 5) the number of iterations
necessary to the non linear conjugate gradient to converge to the
solution of $(\rho_h^*)$ with a stopping test : $|J'(v^n)|^2 \le \epsilon$ is smoothly
increasing with the area of the plastic zones, as the load F increases.

For a typical problem it can require 4 or 5 times more itera-
tions than in the elastic case.
### 4.6. Solution of the elasto-plastic problem by the gradient method with auxiliary operator.

This is the algorithm (see [6]):

$v^0$ being given

by induction $v^n$ being given

determine $v^{n+1}$ solution of the variational equation

$$ (4.1) \qquad \left( g(\varepsilon(v^{n+1})), \varepsilon(w) \right) = \left( g(\varepsilon(v^n)), \varepsilon(w) \right) - \rho \left[ f'(\varepsilon(v^{n+1})), \varepsilon(w) \right] - L(w) $$

Note that the convergence of this method is proved (see [6]) for $\rho$ sufficiently small when the non linear operator $\varepsilon^* \cdot f \cdot \varepsilon$ is coercive, which is not the case here.

Practically the method converges correctly when plasticity is not important, but the number of iterations increases much faster than in penalty-duality algorithm (see Fig. 5).

The program is very similar to the one of the penalty-duality algorithm since both need mainly to solve an elastic problem at each step, with always the same matrix.

The best choice for $\rho$ seems $\rho = 1$.

### 4.7. Solution of the elasto-plastic problem by the penalty-duality algorithm.

The algorithm used is the one presented in §3.4.

The stopping test was:

$$ \sum_T |e^{n+1} - \varepsilon(v^{n+1})| \le \varepsilon_1 $$

with $\varepsilon_1 = 10^{-8}$.

#### Choice of $\alpha$:

When $F$ is small enough, the problem is fully elastic, and it is possible to show that for $\alpha = 1$, the sequence $v^n$ converges in one
iteration, even though $e^n$ and $\sigma^n$ converge only linearly.

When the problem is elasto-plastic, $\alpha$ must certainly be chosen greater than 1. A typical value may be $\alpha = 2$ or 3 (see Fig. 11). **But this choice is not crucial.**

<table><caption>Fig. 3: Results for $n_1 = 6$ and $n_2 = 10$. (140 independent variables).</caption><thead><tr><th>Nb.</th><th>F/F<sub>o</sub></th><th>U</th><th>Nb. of iterations<br>for conjugate gradient</th><th>Nb. of iterations<br>for penalty-duality</th></tr></thead><tbody><tr><td>0</td><td>0.2</td><td>0.101</td><td>133 (elastic)</td><td>--</td></tr><tr><td>1</td><td>0.4</td><td>0.203</td><td>245</td><td>--</td></tr><tr><td>2</td><td>0.5</td><td>0.270</td><td>411</td><td>--</td></tr><tr><td>3</td><td>0.55</td><td>0.322</td><td>547 (CPU=42")<sup>(1)</sup></td><td>35 (CPU= 3")</td></tr><tr><td>4</td><td>0.6</td><td>0.404</td><td>674</td><td>--</td></tr><tr><td>5</td><td>0.65</td><td>0.542</td><td>645</td><td>--</td></tr><tr><td>6</td><td>0.76</td><td>0.92</td><td>800</td><td>--</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 4: Order of appearance of the plastic triangles.</figcaption></figure>

The numbers are referring to the loads of Fig. 3 : e.g. for $F = 0.55 F_o$ the plastic zones observed were all the triangles noted 1, 2, and 3.

<footnote>(1) CPU is the execution time of the program on IBM 370/168.</footnote>
<table><caption>Fig.5 - Results for $n_1$ = 8 and $n_2$ = 14 (252 independent variables)</caption><thead><tr><th colspan="4">Conjugate gradient.Auxiliary operator</th><th colspan="4">Penalty-duality</th></tr><tr><th>Nb</th><th>F/F<sub>o</sub></th><th>U</th><th>Number of iterations</th><th>CPU</th><th>Number of iterations</th><th>CPU</th><th>Number of iterations</th><th>α</th><th>CPU</th></tr></thead><tbody><tr><td>0</td><td>0.2</td><td>0.116</td><td>181</td><td>30"</td><td>1 (elastic)</td><td>--</td><td>1 (elastic)</td><td>--</td><td>--</td></tr><tr><td>1</td><td>0.4</td><td>0.236</td><td>377</td><td>58"</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td></tr><tr><td>2</td><td>0.5</td><td>0.346</td><td>645</td><td>80"</td><td>82</td><td>10"</td><td>51</td><td>2</td><td>6"</td></tr><tr><td>3</td><td>0.55</td><td>0.460</td><td>&gt; 800<sup>(+)</sup></td><td>107"</td><td>--</td><td>--</td><td>57</td><td>2</td><td>7"</td></tr><tr><td>4</td><td>0.6</td><td>0.680</td><td>&gt;1200<sup>(+)</sup></td><td>2'30</td><td>215</td><td>24"</td><td>79</td><td>2.5</td><td>9"</td></tr><tr><td>5</td><td>0.65</td><td>1.00</td><td>--</td><td>--</td><td>&gt; 250<sup>(+)</sup></td><td>28"</td><td>88</td><td>3</td><td>10"</td></tr></tbody></table>

(the cost of one iteration is about the same time for each method).

<figure><img src="image_8.png" /><figcaption>Fig.6 - Order of appearance of the plastic triangles ($n_1$ = 8, $n_2$ = 14).</figcaption></figure>

(as in Fig. 4 : the numbers are referring to the load of Fig. 5).
<figure><img src="image_9.png" /><figcaption>Fig. 7 - Variation of the terminal displacement U in function of the load F.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8 Appearance of the plastic zones (after smoothing).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9 Displacement of the beam (magnified by 2.85) for F = 0.6 × F₀, with (n₁ = 8, n₂ = 14).</figcaption></figure>

Each element corresponds to 4 triangles; the pecked lines show the initial position.
<figure><img src="image_12.png" /><figcaption>Fig. 10 - Evolution of the stresses in one of the corner element (plastic zone) ($n_3 = 6, n_2 = 10$)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11 - Choice of $\alpha$ in the penalty-duality algorithm.</figcaption></figure>
### 4.8. Comments.

Figures 2 and 4 show the performances of the three methods. As one iteration costs about the same in each method, we can see that penalty-duality is about 2 or 3 times faster than auxiliary operator, and 10 or 20 times faster than conjugate gradient in these examples.

Figures 3 and 5 show the order of appearance of the plastic triangles.

Figure 6 shows the variation of the displacement U in function of the load F.

Figure 7 is a smoothing of Figures 3 and 5 to get an idea of how the plastic zones spread for the continuous problem.

Figure 8 shows the straining of the whole beam for $F = 0.6 F_o$.

Figure 9 shows the evolution of $\sigma$ in a plastic triangle on the convex K, when the load F is increasing. The graphic is made in the deviator space, that is taking $\sigma_{12}^D$ as first variable and $\sigma_{11}^D = - \sigma_{12}^D$ as second variable. The plot of the convex $K(\sigma_{ij}^D, \sigma_{ij}^D < 1)$ is then a disk centered at the origin.

We can check on Fig. 9 that Hencky's law is satisfied, that is, for each load F, $\sigma_{ij}^D$ is well the projection on K of $2\mu\varepsilon_{ij}^D(u)$.

As there is no unloading we can consider that Hencky's law is a good approximation of Prandtl and Reuss Law which says that for each increasing of load $\Delta F$, $\Delta\sigma_{ij}^D$ is the projection of K of $2\mu\varepsilon_{ij}^D(\Delta u)$.

At last, Figure 10 shows how to choose $\alpha$ for $F = F_o$, and ($n_1 = n_2 = 3$), (4 elements on 18 were plastic).

**Note 4.1.** : A variation of the penalty-duality algorithm.

We have tried a variation of the penalty-duality algorithm which consists in penalizing and dualizing only $e_{ij}^D = \varepsilon_{ij}^D$.

It is possible because $f'(e)$ (see (1.8)) is non linear only in $e^D$.

It gives approximatively the same numerical results, but the choice of $\alpha$ may be a little bit more difficult. ■
CONCLUSION.

For the elasto-plastic problem arising from Hencky's law,
we are looking for a stress field $\sigma$, and a displacement field $u$.

We saw that $\sigma$ and $u$ are the solutions of two problems in duality : ($\rho$) and ($\rho^*$), ($\rho$) being coercive. If the applied loads are sufficiently small, we have the existence and uniqueness of $\sigma$. If the elastic zones $\Omega_e$ satisfy some hypothesis, then we may have the uniqueness of $u$ on $\Omega_e$.

We chose approximation of $u$ and $\sigma$ by finite elements such that their approximations $\sigma_h$ and $u_h$ be solutions of problems ($\rho_h$) and ($\rho_h^*$) themselves in duality. Using then the coercivity of ($\rho$), we showed that $\sigma_h \to \sigma$ if $h \to 0$.

To solve the approximate problem, we used a penalty-duality algorithm which iterates simultaneously on $\sigma^n$ and $u^n$.

Thus, it can be shown that at least $\sigma^n \to \sigma^h$ when $n \to \infty$ and $u^n \to u^h$ on $\Omega_e$, when $n \to \infty$.

Finally, the numerical results show that this algorithm gives very good results when compared to the non-linear conjugate gradient method, or to the gradient method with auxiliary operator.

-------
BIBLIOGRAPHY.

[1] P.G. CIARLET, P.A. RAVIART: "General Lagrange and Hermite interpolation in $\mathbb{R}^n$ with application to finite element methods". Arch. Rat. Mech. Anal. 46, pp. 177-189 (1972).

[2] G. DUVAUT, J.L. LIONS "Les Inéquations en Mécanique et en Physique" Dunod (1972).

[3] I. EKELAND, R. TEMAM : "Analyse convexe et. Problèmes variationnels" Dunod (1973).

[4] D. GABAY, B. MERCIER :( To appear).

[5] R. GLOWINSKI, A. MARROCCO : "Sur l'approximation, par éléments finis d'ordre 1, et la résolution par pénalisation-dualité, d'une classe de problèmes de Dirichlet non linéaires". C.R.A.S. : t. 278 (24/6/1974).

[6] R. GLOWINSKI, J.L. LIONS, R. TREMOLIERES :" Méthodes numériques de résolution des inéquations en mécanique et en physique" Dunod (1975).

[7] M. HESTENES : "Multiplier and gradient methods". J.O.T.A. 4, N° 5 (1969).

[8] H. LANCHON : Thèse d'Etat, Paris (1972).

[9] B. MERCIER : "Recherche numérique d'un champ de déplacements en élasto-plasticité". Rapport LABORIA N° 71 (1974).

[10] B. MERCIER: "Approximation par éléments finis et résolution par un algorithme de pénalisation-dualité d'un problème d'élasto-plasticité". C.R.A.S(to appear in 1975).

[11] O.C. ZIENKIEWICZ : "The finite element method in engineering science" McGrawhill (1970).