The Prandtl-Reuss and Hencky model in perfect plasticity

The Prandtl-Reuss model in the theory of linear-elastic perfect plastic material in classical notation reads (cf. [14] and [7])

$$\sigma \epsilon \dot{u} \dot{\sigma} \lambda \Omega$$ (25)

$$\lambda \tau \sigma \leq \forall \tau \cal {F} \tau \leq \Omega \lambda \dot{\sigma} \Omega$$ (26)

$$\dot{u} \Gamma_{u}^{} \sigma \Gamma_{\sigma}^{}$$ (27)

This set of equations describes the deformation of an elasto-plastic body occupying a bounded domain $\Omega$ $\subset$ $\mathbb {R}$^d ( d = 2  or  3) which is fixed along a part  $\Gamma_{u}^{}$  of its boundary  $\partial$$\Omega$ , under the action of a body force with density f and a surface traction  g  along  $\Gamma_{\sigma}^{}$ = $\partial$$\Omega$ $\setminus$ $\Gamma_{u}^{}$ . Here, u denotes the displacement, $\epsilon$(u) = ${\frac{1}{2}}$($\nabla$u + $\nabla$u^T) the strain- and $\sigma$ the stress-tensor. A is the material tensor which is assumed to be symmetric and positive definite with inverse  C = A^-1. In addition, the plastic growth is denoted by $\lambda$, and $\cal {F}$ represents the continuous convex flow function defined by the von Mises yield function

$$\cal {F}$$($$\sigma$$) : = |$$\sigma^{D}_{}$$| - $$\sigma_{0}^{}$$ ,

for some  $\sigma_{0}^{}$ > 0 , depending on the deviatoric part of the stress tensor,  $\sigma^{D}_{}$ : = $\sigma$ - ${\tfrac{1}{3}}$tr($\sigma$)I. We consider the case of a quasi-static process, i.e., acceleration effects are neglected.

For applying a finite element method we have to reformulate problem (23) in a variational setting. Typically the behaviour of plastic material is a time dependent process on a time interval I : = [0, T] parametrizing the growing loading, i.e., the load functions are given in the form

f = f (t) = tf₀ ,        g = g(t) = tg₀ ,

with some prescribed load distributions  f₀(x)  and  g₀(x) . In order to give the variational formulation, we use the function spaces

$$\Omega \Omega \mathbb {R} \Omega \Omega \mathbb {R}$$

$$\tau \in \Omega \tau \in \Omega$$

and the sets of admissible stresses

$$\tau \in \tau \Omega \tau \Gamma_{\sigma}^{}$$ H^div_{f, g}

$$\Pi \Sigma \tau \in \Sigma \cal {F} \tau \leq \Omega$$

for any appropriate function space  $\Sigma$ . Furthermore, we introduce the velocities of the displacement, v = $\dot{u}$. We assume the process to start with a stress free state, i.e., $\sigma$(0) = 0. Using this notation the dual-mixed formulation seeks a pair {v,$\sigma$} : I$\to$L²($\Omega$)^d x $\Pi$H^div _{* , g}, such that  $\sigma$(0) = 0, and

$$\dot{\sigma} \tau \sigma \tau \sigma \ge \forall \tau \in \Pi$$ (28)

$$\sigma \varphi \varphi \forall \varphi \in \Omega$$ (29)

Integrating in (24) by parts, we obtain the primal-mixed variational formulation, where a pair {v,$\sigma$} : I$\to$V x $\Pi$L²($\Omega$)^{d x d}_sym is sought such that  $\sigma$(0) = 0, and

$$\dot{\sigma} \tau \sigma \epsilon \tau \sigma \ge \forall \tau \in \Pi \Omega$$ (30)

$$\sigma \epsilon \varphi \varphi \varphi \Gamma_{\sigma} \forall \varphi \in$$ (31)

Our computations for the elasto-plasticity problem are all based on this primal-mixed formulation.

If the rate dependence in the plastic law in (23),  $\epsilon$($\dot{u}$) = A$\dot{\sigma}$ + $\lambda$ , is neglected, we arrive at the Hencky model of (stationary) perfect plasticity (see, e.g., [7]). For this the primal-mixed variational formulation reads

$$\sigma \tau \sigma \epsilon \tau \sigma \ge \forall \tau \in \Pi \Omega$$ (32)

$$\sigma \epsilon \varphi \varphi \varphi \Gamma_{\sigma} \forall \varphi \in$$ (33)

The adventage of this formulation is that it can again be rewritten as a nonlinear variational equation analogously as for the Strang model problem,

$$\epsilon \epsilon \varphi \Pi \epsilon \epsilon \varphi \varphi \varphi \Gamma_{\sigma} \forall \varphi \in$$ (34)

where

$$\Pi \tau \begin{cases} \quad \tau \,, &\text{if} \;\; \vert\tau^D\ver... ... \,, &\text{if} \;\; \vert\tau^D\vert>\sigma _0 \,. \end{cases}$$ (35)

Splitting the continuous loading process  f (t) , g(t) in the Prandtl-Reuss model into a sequence of incremental load steps,

 fⁿ = f^{n - 1} + k_nf₀ ,        gⁿ = g^{n - 1} + k_ng₀ ,

with a pseudo-time step length  k_n : = t_n - t_{n - 1}, we obtain a sequence of Hencky-type problems,

$$\sigma^{n}_{} \tau \sigma^{n}_{} \sigma^{n-1}_{} \epsilon \tau \sigma^{n}_{} \ge \forall \tau \in \Pi \Omega$$ (36)

$$\sigma^{n}_{} \epsilon \varphi \varphi \varphi \Gamma_{\sigma} \forall \varphi \in$$ (37)

with the initial value  $\sigma^{0}_{}$ = 0. According to the above discussion, each such load step is equivalent to a nonlinear problem of the form

$$\epsilon \epsilon \varphi \Pi \sigma^{n-1}_{} \epsilon \epsilon \varphi \varphi \varphi \Gamma_{\sigma} \forall \varphi \in$$ (38)

Solving this equation gives us an update for the deformation velocity  vⁿ $\in$ V and then the new stress and deformation as

$$\sigma^{n}_{}$$ : = $$\Pi$$($$\sigma^{n-1}_{}$$ + k_nC$$\epsilon$$(vⁿ))   $$\in$$  $$\Pi$$L²($$\Omega$$)^{d x d}_sym ,    uⁿ : = u^{n - 1} + k_nvⁿ   $$\in$$  V .

This process is equivalent to applying the backward Euler time-stepping scheme to the nonstationary problem (25). We will assume that the incremental step  k_nf₀ , k_ng₀  is chosen small enough that the resulting time-discretization error can be neglected compared to the error resulting from the finite element approximation. Our numerical tests show that this is realistic, at least for the type of problems considered. The stationary Hencky model (26) may be viewed as the approximation of the Prandtl-Reuss model by one time step of length  k = 1  starting form the initial state  $\sigma^{0}_{}$ $\equiv$ 0 , u⁰ $\equiv$ 0.

Finally, the nonlinear problems (30) are approximated by a damped Newton iteration. Starting from the result at the preceding load level,  v^{n, 0} : = v^{n - 1}, the iteration step  v^{n, i - 1} $\rightarrow$ v^{n, i}  reads

$$\epsilon \epsilon \delta \epsilon \varphi \varphi \varphi \Gamma_{\sigma} \epsilon \epsilon \varphi \forall \varphi \in$$ (39)

with the Jacobian  C'( ^. )  of the function  C( ^. ) , followed by the update

v^{n, i} = v^{n, i - 1} + $$\lambda_{i}^{}$$$$\delta$$v^{n, i} .

The damping parameter  0 < $\lambda_{i}^{}$$\le$1  is determined in the form  $\lambda_{i}^{}$ = 2^-r  such that the residual norm is decreased. Hence, we have reduced the nonlinear plasticity problem (23) on the continuous level to a sequence of linear problems of the form

$$\epsilon \epsilon \epsilon \varphi \varphi \forall \varphi \in$$ (40)

with certain functionals  l ( ^. )  defined on  V. These may be solved by the CR-method with multigrid acceleration. The use of multigrid is rather natural as in the course of the mesh refinement process a sequence of nested meshes is automatically generated. Since the emphasis of this paper is on the aspect of a posteriori error estimation and mesh design, we do not go here into the details of the algebraic solution techniques and instead refer to [33] and [34]. Again the error estimation uses the general approach described above applied to the nonlinear problem (30), i.e., the meshes are kept fixed during the Newton iteration, but are dynamically adapted in the coarse of the time stepping. The corresponding linearized dual problems

$$\epsilon \epsilon \varphi \epsilon \varphi \forall \varphi \in$$ (41)

are solved by the same technique as used for solving the Newton steps (31).

We remark that, in contrast to really time dependent problems involving acceleration terms, the control of the time-stepping error in the quasi-stationary Prandtl-Reuss model is less critical. Since in each incremental load step the deformation velocity  vⁿ = $\dot{u}^{n}_{}$  is updated, we do not expect much accumulation of these local errors over the loading path. Hence it seems justified to estimate the global spatial error in this process by simply treating it as a sequence of stationary problems. This approach, of course, would be disastrous in the case of a parabolic problem like, e.g., the heat equation.