---
abstract: 'We devise and evaluate numerically a Hybrid High-Order (HHO) method for finite plasticity within a logarithmic strain framework. The HHO method uses as discrete unknowns piecewise polynomials of order $k\ge1$ on the mesh skeleton, together with cell-based polynomials that can be eliminated locally by static condensation. The HHO method leads to a primal formulation, supports polyhedral meshes with non-matching interfaces, is free of volumetric locking, the integration of the behavior law is performed only at cell-based quadrature nodes, and the tangent matrix in Newton’s method is symmetric. Moreover, the principle of virtual work is satisfied locally with equilibrated tractions. Various two- and three-dimensional benchmarks are presented, as well as comparisons against known solutions obtained with an industrial software using conforming and mixed finite elements.'
address:
- ', '
- ' , '
- ' , '
- ' , '
author:
- Mickaël Abbas
- Alexandre Ern
- 'Nicolas Pignet\*'
- 'Mickaël Abbas$^1,2$, Alexandre Ern$^{3,4}$ and Nicolas Pignet$^{1,2,3,4}$'
bibliography:
- 'Bibliographie.bib'
title:
- 'A Hybrid High-Order method for finite elastoplastic deformations within a logarithmic strain framework'
- 'A Hybrid High-Order method for finite elastoplastic deformations within a logarithmic strain framework'
---

[**Keywords:**]{} Finite strain plasticity – Hybrid High-Order methods – Polyhedral meshes – Locking-free\

Introduction
============

Modelling plasticity, particularly in the regime of finite deformations, is of a major importance in industrial applications since this is one of the main nonlinearites that can be encountered in nonlinear solid mechanics. Moreover, finite elastoplastic deformations have a major influence on the life time of a mechanical structure. The present contribution is an extension to the finite strain regime of the Hybrid High-Order (HHO) method for incremental associative plasticity with small deformations [@AbErPi2018a]. This extension hinges on a logarithmic strain framework [@Miehe2002] for anisotropic finite elastoplasticity. This framework provides a natural extension of small elastoplastic deformations to finite elastoplastic deformations by means of purely geometric transformations. Indeed, the weak form of the plasticity problem is derived from the minimization of an energy functional based on an incremental pseudo-energy density.

The present work aims at addressing the following important issues. Firstly, the incompressibility of the plastic deformations generally leads to volumetric locking when employing a continuous Galerkin (cG) approximation based on low-order $H^1$-conforming finite elements. In these methods, only the displacement field is approximated globally, whereas the variables associated with the plastic behavior are computed locally in each mesh cell (typically at quadrature nodes). A way to circumvent the volumetric locking is to use high-order $H^1$-conforming finite elements or NURBS for small [@Elguedj2014] and finite [@Elguedj2008] elastoplastic deformations. Therein, the displacement is still the only field which is approximated globally. However, the resulting discrete problem is more costly to solve because of the larger support of the basis functions. Another possible way to prevent volumetric locking is to introduce additional global unknowns as in the Enhanced Assumed Strain (EAS) methods [@Simo1992a] and in mixed methods [@AlAkhrass2014; @Saracibar2006; @SouzaNeto2005; @Simo1985] on simplicial or hexahedral meshes (the variables associated with the plastic behavior are still computed locally). However, the introduction of additional globally coupled unknowns generally increases the cost of building and solving the discrete problem. Moreover, devising mixed methods on polyhedral meshes with non-matching interfaces is a delicate question. On the positive side, cG methods as well as EAS and mixed methods require to perform the integration of the behavior law only at the quadrature nodes in the mesh cells. Another class of methods free of volumetric locking are discontinuous Galerkin (dG) methods. We mention in particular [@Eyck2006; @Eyck2008; @Eyck2008a; @Noels2006] for hyperelasticity. Interior penalty dG methods have been developed for classical plasticity with small [@Hansbo2010; @Liu2010] and finite [@Liu2013] deformations, and for gradient plasticity with small [@Djoko2007; @Djoko2007a] and finite [@McBride2009] deformations. However, dG methods from the literature generally require to perform the integration of the behavior law also at additional quadrature nodes located at the mesh faces. Moreover, if the plasticity problem is solved using a Newton’s method, which is often the case, the tangent matrix from the dG formulation is generally non-symmetric owing to the nonlinear nature of the consistency term. Thus, the solving cost can increase significantly, particularly with iterative solvers (the memory requirements can become important for direct solvers). We also mention the lowest-order Virtual Element Method (VEM) for inelastic problems with small deformations [@BeiraodaVeiga2015] (and its two-dimensional higher-order extension [@Artioli2017]), whereas the case of finite deformations is treated in [@Wriggers2017a; @Hudobivnik2018], still in the lowest-order case. We also mention the recent study of low-order hybrid dG methods with conforming traces [@Wulfinghoff2017] and the hybridizable weakly conforming Galerkin method with nonconforming traces in the context of linear [@Kramer2016] and nonlinear [@Bayat2018] solid mechanics. Moreover, finite volume methods for plasticity problems have been developped for small deformations [@Taylor2003] and for large deformations [@Barton2010; @Cardiff2014].

In the present work, we devise and evaluate numerically a HHO method for finite plasticity within a logarithmic strain framework. HHO methods have been introduced a few years ago for diffusion problems [@DiPEL:2014] and for linear elasticity problems [@DiPEr:2015]. Since then, the development of HHO methods has received a vigorous interest. Examples include in solids mechanics Biot’s problem [@BoBoD:2016], nonlinear elasticity [@BoDPS:2017] and associative plasticity [@AbErPi2018a] with small deformations, and hyperelasticity with finite deformations [@AbErPi2018], and in fluid mechanics, the incompressible Stokes equations [@DPELS:16], the steady incompressible Navier–Stokes equations [@DiPKr:2017], and viscoplatic flows with yield stress [@Cascavita2018]. The discrete unknowns in HHO methods in computational mechanics are face-based vector-valued polynomials of arbitrary order $k \geq 1$ on the mesh skeleton. Cell-based vector-valued polynomials are also introduced for the stability and approximation properties of the method. These cell-based vector-valued polynomials are eliminated locally by using the well-known static condensation technique (based on a local Schur complement).

The devising of HHO methods hinges on two key ideas: a local reconstruction operator acting on the face and cell unknowns that builds a tensor-valued polynomial representing the displacement gradient in the polynomial space $\Pkd(T; \Rdd)$, where $T$ is a generic mesh cell and $d$ is the space dimension [@BoDPS:2017; @AbErPi2018]; a local stabilization operator that weakly enforces on each mesh face the consistency between the local face unknowns and the trace of the cell unknowns [@DiPEL:2014; @DiPEr:2015]. HHO methods offer several advantages: general meshes (including fairly general polyhedral mesh cells and non-matching interfaces) are supported; a local formulation using equilibrated fluxes is available; computational benefits owing to the static condensation of the cell unknowns and the higher-order convergence rates, and the construction is dimension-independent. Moreover, an open-source implementation of HHO methods, the `DiSk++` library, is available using generic programming tools [@CicDE:2018]. In computational mechanics, other salient features of HHO methods are: a displacement-based formulation avoiding the need to introduce additional globally coupled unknowns; absence of volumetric locking; the integration of the behavior law only at the cell quadrature nodes; and the tangent matrix arising in the Newton’s method is symmetric. Furthermore, HHO methods have been bridged [@CoDPE:2016] to Hybridizable Discontinuous Galerkin (HDG) methods [@CoGoL:09] and to nonconforming Virtual Element Methods (ncVEM) [@Ayuso2016]. The essential difference with HDG methods is that the HHO stabilization is different so as to deliver $O(h^{k+1})$ energy-error estimates for linear model problems with smooth solutions on general meshes, where $h$ is the mesh-size. Concerning ncVEM, the devising viewpoint is different (ncVEM considers the computable projection of virtual functions instead of a reconstruction operator), whereas the stabilization achieves similar convergence rates as HHO but is written differently. We also notice that, to our knowledge, HDG methods have not yet been devised for finite elastoplasticity problems (in contrast to hyperelasticity problems [@KaLeC:2015; @Nguyen2012a]). Owing to the close links between HHO and HDG methods, this work can thus be seen as the first HDG-like method for plasticity problems in finite deformations.

This paper is organized as follows: in Section 2, we present the plasticity model within a logarithmic strain framework and the weak formulation of the governing equations. In Section 3, we devise the HHO method and highlight some of its theoretical aspects. In Section 4, we investigate numerically the HHO method on two- and three-dimensional benchmarks from the literature, and we compare our results to analytical solutions whenever available and to numerical results obtained using established cG and mixed methods implemented in the open-source industrial software [@CodeAster].

Plasticity model {#sec::model}
================

In what follows, we write $v$ or $V$ for scalar-valued fields, $\vv$ or $\vecteur{V}$ for vector-valued fields, $\matrice{v}$ or $\matrice{V}$ for second-order tensor-valued fields, and $\tenseur{4}{V}$ for fourth-order tensor-valued fields. Contrary to the hyperelastic model, the elastoplastic model is based on the assumption that the deformations are no longer reversible.

Kinematics
----------

Let $\mathcal{B}_{0}$ be an elastoplastic material body that occupies the domain $\Omega_0$ in the reference configuration. Here, $\Omega_0 \subset \Rd$, $d \in \{2,3\}$, is a bounded connected polyhedral domain with Lipschitz boundary $\Gamma := \partial \Omega_0$ decomposed in the two relatively open subsets $\Bn$ and $\Bd$, where a Neumann and a Dirichlet condition is enforced respectively, and such that $\overline{\Bn} \cup \overline{\Bd} = \Gamma$, $\Bn \cap \Bd = \emptyset $, and $\Bd$ has positive Hausdorff-measure (so as to prevent rigid-body motions). Due to the deformation, a point $\Xo \in \Omega_0$ is mapped to a point $\xc(t) = \Xo + \vu(\Xo,t)$ in the equilibrium configuration, where $\vu : \Omega_0 \times I \rightarrow \Rd$ is the displacement mapping. The deformation gradient $\Fdef(\vu) = \matrice{I} +\gradX \, \vu$ takes values in $\Rdd_{+}$ which is the set of $\Rdd$-matrices with positive determinant. In what follows, the gradient and divergence operators are taken with respect to the coordinate $\Xo$ of the reference configuration (we use the subscript $X$ to indicate it).

We use the logarithmic strain framework [@Miehe2002] developed for anisotropic finite elastoplasticity. Hence, it allows us to define the logarithmic strain tensor $\Elog \in \Msym$ as $$\Elog := \frac{1}{2} \ln(\Fdef^T \Fdef).$$ This measure of the deformations $\Elog$ is objective. Moreover, if the eigenvectors of $\Elog$ do not change with time (the eigenvalues may change in time), then $\dot{\Elog} = \dot{\matrice{U}}\matrice{U}^{-1}$, where $\matrice{U} \in \Msym $ is the right stretch tensor from the polar decomposition $\Fdef = \matrice{R} \matrice{U}$. The plastic deformations are measured by means of the plastic logarithmic strain tensor $\pElog \in \Msym$. We assume the following additive decomposition of the logarithmic strain tensor $\Elog$: $$\eElog := \Elog - \pElog,$$ where $\eElog \in \Msym$ is the elastic logarithmic strain tensor. Finally, the plastic strains are assumed to be incompressible, i.e. $$\trace \pElog = 0.$$

Constitutive logarithmic strain model {#ss:helmhotz}
-------------------------------------

In what follows, we place ourselves within the framework of generalized standard materials initially introduced in [@Halphen1975] and further developed in [@Lemaitre1994]. Moreover, the plasticity model is assumed to be strain-hardening (or perfect) and rate-independent, i.e., the time and the speed of the deformations have no influence on the strains. For this reason, only the incremental plasticity problem is considered. The local material state is described by the logarithmic strain tensor $\Elog \in \Msym$, the plastic logarithmic strain tensor $\pElog \in \Msym$, and a finite collection of internal variables $\IV:= ( \alpha_1, \cdots, \alpha_m) \in \Reel^m$, which contain typically at least the equivalent plastic strain $p$, see Sect. \[ss:nl\_hardening\] for a simple example or [@Lemaitre1994] for more detailed examples. For simplicity, we denote $\IVG:=( \pElog, \IV ) \in \IVGS$ the generalized internal variables, where the space of the generalized internal variables is $$\label{eq:def_space_X}
\IVGS := \left\lbrace \IVG=(\pElog, \IV) \in  \Msym \times \Reel^{m}  \: | \: \trace(\pElog)= 0 \right\rbrace.$$ Moreover, we assume that there exists a Helmholtz free energy $\psi: \Msym \times \Reel^m \rightarrow \Reel$ acting on the pair $(\eElog, \IV)$ and satisfying the following hypothesis.

\[hyp\_free\_ener\] $\psi$ can be decomposed additively into an elastic and a plastic part as follows: $$\psi(\eElog, \IV) = \psi^e( \eElog) + \psi^p( \IV).$$ where the function $\psi^p$ is strictly convex and the function $\psi^e$ is polyconvex.

Following the second principle of thermodynamics, the logarithmic stress tensor $\stressT \in \Msym$ and the thermodynamic forces $\IF$ are derived from $\psi$ as follows: $$\label{eq_IF}
\stressT = \partial_{ \eElog} \psi^e(\eElog)  \quad \textrm{and} \quad
\IF =   \partial_{ \IV} \psi^p(\IV).$$ The criterion to determine whether the deformations become plastic hinges on the scalar-valued yield function $\Phi: \Msym \times \Reel^m \rightarrow \Reel$, which is a continuous and convex function of the logarithmic stress tensor $\stressT$ and the thermodynamic forces $\IF$. Letting $\mathcal{A} := \left\lbrace (\stressT, \IF) \in \Msym \times \Reel^m \, \vert \, \Phi(\stressT, \IF) \leq 0 \right\rbrace$ be the convex set of admissible states, the elastic domain $\mathcal{A}^e$ is composed of all the pairs $ (\stressT, \IF)$ such that $\Phi(\stressT, \IF) < 0 $, and the yield surface $\partial \mathcal{A}$ of all the pairs $ (\stressT, \IF)$ such that $\Phi(\stressT, \IF) = 0 $.

\[hyp\_yield\] The yield function $\Phi: \Msym \times \Reel^m \rightarrow \Reel$ satisfies the following properties: $\Phi$ is a piecewise analytical function; the point $(\matrice{0}, \vecteur{0})$ lies in the elastic domain, i.e., $\Phi(\matrice{0}, \vecteur{0}) < 0$; and $\Phi$ is differentiable at all points on the yield surface $\partial\mathcal{A}$.

Finally, the incremental dissipation function $D :\IVGS \rightarrow \Reel$ is a convex function which is positively homogeneous of degree one and is defined as follows: $$\label{eq::dissipation}
D(\textrm{d} \IVG) = \sup_{(\stressT, \IF) \in \mathcal{A}} \left( \stressT : \textrm{d} \pElog - \IF \SCAL \textrm{d} \IV  \right),$$ where $\textrm{d} \IVG, \textrm{d} \pElog$, and $\textrm{d} \IV $ are the finite increments of $ \IVG, \pElog$, and $ \IV $, respectively.

Finite elastoplasticity model problem in incremental form {#ss:mech_model}
---------------------------------------------------------

We are interested in finding the quasi-static evolution in the pseudo-time interval $I=[0,t_F]$, $t_F > 0$, of the elastoplastic material body $\mathcal{B}_0$. The pseudo-time interval $I$ is discretized into $N$ subintervals such that $t^0=0 < t^1 < \cdots < t^N = t_F$. The evolution occurs under the action of a body force $\loadext:\Omega_0 \times I \rightarrow \Rd $, a traction force $\Tn: \Bn \times I \rightarrow \Rd$ exerted on the Neumann boundary $\Bn$, and a displacement $\vu{}_{\textrm{D}}: \Bd \times I \rightarrow \Rd$ imposed on the Dirichlet boundary $\Bd$. We denote by $V_{\textrm{D}}^n$, resp. $V_0$, the set of all kinematically admissible displacements which satisfy the Dirichlet conditions, resp. homogeneous Dirichlet conditions on $\Bd$: $$V_{\textrm{D}}^n = \left\lbrace \vv \in H^1(\Omega_0; \Rd) \: | \: \vv = \vu{}_{\textrm{D}}(t^n) \: \mbox{ on } \Bd \right\rbrace, \; V_0 = \left\lbrace \vv \in H^1(\Omega_0; \Rd) \: | \: \vv = \vecteur{0} \: \mbox{ on } \Bd \right\rbrace.$$ Following [@Miehe2002], we define for any pseudo-time step $1 \leq n \leq N$, the incremental pseudo-energy density $\Psi: \Rdd_{+} \times \IVGS \rightarrow \Reel$ acting on the pair $ (\Fdef, \IVG)$ such that $$\label{eq::pseudo-energy}
\Psi(\Fdef, \IVG) = \left\lbrace \big( \psi^e(\eElog) + \psi^p(\IV) \big) - \big( \psi^e(\eElog(\vu^{n-1}))  + \psi^p(\IV^{n-1}) \big) \right\rbrace + D(\IVG -\IVG^{n-1}),$$ where $\vu^{n-1} \in V_{\textrm{D}}^{n-1}$ and $\IVG^{n-1} \in  L^2(\Omega_0; \IVGS)$ are given from the previous pseudo-time step or the initial condition. Note that the second term in , which is evaluated at $t^{n-1}$, is irrelevant for minimization purposes. It is added so that the pseudo-energy is in a time-incremental form. This allows us to define the energy functional $\mathcal{E}^n: V_{\textrm{D}}^n \times L^2(\Omega_0; \IVGS)  \rightarrow \Reel$ such that, for all kinematically admissible displacement fields $\vv \in V_{\textrm{D}}^n$ and all generalized internal variables $\IVG \in L^2(\Omega_0; \IVGS)$. $$\label{eq:def_calE}
\mathcal{E}^n (\vv, \IVG) = \int_{\Omega_0} \Psi(\Fdef(\vv), \IVG) \dV - \int_{\Omega_0} \loadext^n\SCAL\vv \dV - \int_{\Bn} \Tn^n \SCAL\vv \dA.$$ The quasi-static equilibrium of the elastoplastic body $\mathcal{B}_0$ is determined at each pseudo-time step $1 \leq n \leq N$ by finding a displacement field $\vu^n \in V_{\textrm{D}}^n$ and generalized internal variables $\IVG^n \in L^2(\Omega_0; \IVGS)$ which minimize the energy functional $\mathcal{E}^n$ in , i.e., $$\label{argmin}
(\vu^n, \IVG^n) \in \argmin_{\vv \in V_{\textrm{D}}^n, \, \IVG \in L^2(\Omega_0; \IVGS) }\mathcal{E}^n (\vv, \IVG).$$ On the one hand, the first variation of $\mathcal{E}^n$ with respect to the displacement field leads to: $$\label{weak_form}
0 = D \mathcal{E}^n(\vu^n, \IVG^n)[\vv] =  \int_{\Omega_0} \PK^n : \gradX\vv \,d\Omega_0 - \int_{\Omega_0} \loadext^n \SCAL \vv \,d\Omega_0 - \int_{\Bn} \Tn^n \SCAL \vv \,d\Gamma, \quad  \textrm{for all } \vv\in V_0,$$ where $\PK^n := \partial_{\Fdef} \Psi(\Fdef(\vu^n), \IVG^n)$ is the first Piola–Kirchhoff stress tensor. On the other hand, the first variation of $\mathcal{E}^n$ with respect to the generalized internal variables leads to the following incremental nonlinear equations (see [@Miehe2002; @Simo1992b] for example): $$\label{normal_rules}
\pElogn - \pElogpn = \Lambda(\stressT^n, \IF^n) \, \partial_{\stressT} \Phi(\stressT^n, \IF^n), \quad \textrm{and} \quad  \IV^n - \IV^{n-1} = -  \Lambda(\stressT^n, \IF^n) \, \partial_{\IF} \Phi(\stressT^n, \IF^n),$$ where the plastic multiplier $\Lambda(\stressT^n, \IF^n)$ verifies the Kuhn–Tucker conditions $$\label{consistency_cond}
\Lambda(\stressT^n, \IF^n) \geq 0, \quad \Phi(\stressT^n, \IF^n) \leq 0, \quad \textrm{and} \quad \Lambda(\stressT^n, \IF^n) \,\Phi(\stressT^n, \IF^n) = 0.$$ We assume additionally that the incremental pseudo-energy density $\Psi$ is polyconvex so that local minimizers of the energy functional exist (cf. e.g [@Miehe2002; @Ball1976]). Thus, the minimization problem can be reformulated, in a more classical way as follows: For all $1 \leq n \leq N$, given $\vu^{n-1} \in V_{\textrm{D}}^{n-1}$ and $\IVG^{n-1} \in  L^2(\Omega_0; \IVGS)$ from the previous pseudo-time step or the initial condition, find $\vu^n \in V_{\textrm{D}}^{n}$ and $\IVG^n \in  L^2(\Omega_0; \IVGS)$ such that

\[weak\_form\_cont\] $$\label{weak_equil}
  \int_{\Omega_0} \PK^n : \gradX\vv \,d\Omega_0 = \int_{\Omega_0} \loadext^n \SCAL \vv \,d\Omega_0 + \int_{\Bn} \Tn^n \SCAL \vv \,d\Gamma, \quad  \textrm{for all } \vv\in V_0,$$ $$\label{weak_plas}
( \IVG^n, \PK^n, \epmodulePKn) =  \textrm{FINITE\_PLASTICITY}(\IVG^{n-1}, \Fdef(\vu^{n-1}), \Fdef(\vu^n) ).$$

The procedure $\textrm{FINITE\_PLASTICITY}$ allows one to compute the new values of the generalized internal variables $\IVG$, the first Piola–Kirchhoff stress tensor $\PK$ and the consistent nominal elastoplastic tangent modulus $\epmodulePK$ at each pseudo-time step. This procedure is detailed in Section \[ss:flow\_rules\].

Algorithmic aspects {#ss:flow_rules}
-------------------

The incremental elastoplasticity problem that has to be solved is to find the new value, after incrementation, of the generalized internal variables $\IVG^{\textrm{new}} = (\Elog^{p,\textrm{new}}, \IV^{\textrm{new}}) \in \IVGS$, the first Piola–Kirchhoff stress tensor $\PK^{\textrm{new}} \in \Msym$, and the consistent nominal elastoplastic tangent modulus $\epmodulePKnew$, given the generalized internal variables $\IVG \in \IVGS$, the deformation gradient $\Fdef \in \Rdd_{+}$, and the new value of the deformation gradient $\Fdef^{\textrm{new}} \in \Rdd_{+}$. Solving this problem is denoted as previously $$(\IVG^{\textrm{new}}, \PK^{\textrm{new}}, \epmodulePKnew)= \textrm{FINITE\_PLASTICITY}(\IVG, \Fdef, \Fdef^{\textrm{new}}).$$ The procedure to compute $(\IVG^{\textrm{new}}, \PK^{\textrm{new}}, \epmodulePKnew)$ is described in Algorithm \[algo::plasticity\] and is composed of three different steps. Firstly, a geometric pre-processing is applied in order to compute the logarithmic strain tensors $\Elog$ and $\Elog^{\textrm{new}}$. Secondly, the procedure $\textrm{SMALL\_PLASTICITY}$ is used to solve the nonlinear incremental problem - so as to compute $(\IVG^{\textrm{new}}, \stressT^{\textrm{new}}, \epmoduleTnew)$. The resolution of - inside the procedure $\textrm{SMALL\_PLASTICITY}$ requires to solve a constrained nonlinear problem which is the same as in the case of plasticity with small deformations and thus makes it possible to extend the procedures already developed for small deformations to finite deformations without modifications (further details about the procedure $\textrm{SMALL\_PLASTICITY}$ can be found in [@AbErPi2018a Sect. 2.3]). One significant example of such a procedure is the standard radial return mapping [@Simo1992b; @Simo1998]. Finally, a geometric post-processing step is applied to compute the new values of the first Piola–Kirchhoff stress tensor $\PK^{\textrm{new}}$ and the consistent nominal elastoplastic tangent modulus $\epmodulePKnew$ from the logarithmic stress tensor $\stressT^{\textrm{new}}$ and the consistent logarithmic elastoplastic tangent modulus $\epmoduleTnew$. Detailed explanations to compute the pre- and post-processing steps are given in [@Miehe2002 Box. 4].

Set $\Elog = \frac{1}{2}\ln(\Fdef^T \Fdef)$, $\Elog^{\textrm{new}} = \frac{1}{2}\ln(\Fdef^{\textrm{new},T} \Fdef^{\textrm{new}})$ and $\mathrm{d}\Elog := \Elog^{\textrm{new}}- \Elog$ Compute $(\IVG^{\textrm{new}}, \stressT^{\textrm{new}}, \epmoduleTnew)= \textrm{SMALL\_PLASTICITY}(\IVG, \Elog, \mathrm{d}{\Elog})$. Compute $\PK^{\textrm{new}} =  \stressT^{\textrm{new}} : (\partial_{\Fdef} \Elog)^{\textrm{new}}$ and $ \epmodulePKnew = (\partial_{\Fdef} \Elog)^{\textrm{new},T} : \epmoduleTnew : (\partial_{\Fdef} \Elog)^{\textrm{new}} + \stressT^{\textrm{new}} : (\partial_{\Fdef \Fdef} \Elog)^{\textrm{new}}$ $(\IVG^{\textrm{new}} , \PK^{\textrm{new}}, \epmodulePKnew)$

Note that $\epmodulePKnew := \partial_{\Fdef \Fdef} \Psi(\Fdef^{\textrm{new}}, \IVG^{\textrm{new}})$ is the consistent elastoplastic tangent modulus and is a fourth-order tensor having only the major symmetries contrary to $\epmoduleTnew$ which has the major and minor symmetries. For a finite incremental strain, the consistent elastoplastic tangent modulus generally differs from the so-called continuous elastoplastic tangent modulus which is obtained by letting the incremental strain tend to zero [@Simo1985a].

Example: nonlinear isotropic hardening with a von Mises yield criterion {#ss:nl_hardening}
-----------------------------------------------------------------------

An illustration of the plasticity model defined above is the nonlinear isotropic hardening model with a von Mises criterion. The elastic part of the free energy is such that $$\psi^e(\eElog) = \frac12 \eElog : \elasticmodule:\eElog,$$ where the elastic modulus is $\elasticmodule = 2\mu \tenseur{4}{I}^{s} + \lambda \matrice{I} \otimes \matrice{I}$, with $\mu>0$, $3\lambda + 2 \mu > 0$, $(\tenseur{4}{I}^{s})_{ij,kl}= \frac{1}{2}(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})$, and $(\matrice{I} \otimes \matrice{I})_{ij,kl}=\delta_{ij}\delta_{kl}$. The internal variable is $\IV := p$, where $p \geq 0$ is the equivalent plastic strain. The plastic part of the free energy is such that $$\psi^p(\IV) = \sigma_{y,0}p +   \frac{H}{2} p^2 + (\sigma_{y,\infty} -\sigma_{y,0})(p - \frac{1 - e^{-\delta p}}{\delta}),$$ where $H \geq 0$ is the isotropic hardening modulus, $ \sigma_{y,0} >0$, resp. $ \sigma_{y,\infty} \geq 0$, is the initial, resp. infinite, yield stress and $\delta \geq 0$ is the saturation parameter. The associated thermodynamic force $\IF = \sigma_{y,0} + Hp + (\sigma_{y,\infty} -\sigma_{y,0})(1 - e^{-\delta p})$ is called the internal stress. Concerning the yield function, we consider a $J_2$-plasticity model with a von Mises criterion: $$\Phi(\stressT, \IF)  = \sqrt{\frac{3}{2}} \Vert \dev( \stressT) \Vert_{\matrice{\ell}^2} - \IF,$$ where $\dev(\matrice{\tau}) := \matrice{\tau} - \frac{1}{d} \trace(\matrice{\tau})\matrice{I}$ is the deviatoric operator, and the Frobenius norm is defined as $\Vert\matrice{\tau}\Vert_{\matrice{\ell}^2} = \sqrt{ \matrice{\tau}  : \matrice{\tau} }$, for all $\matrice{\tau} \in \Reel^{d \times d}$. Moreover, the above model describes with a reasonable accuracy the behaviour of metals [@Lemaitre1994]. This model is used for the numerical examples presented in Section \[sec::sec\_numexp\].

The Hybrid High-Order method {#sec:HHO}
============================

Discrete setting {#ss::discrete_setting}
----------------

We consider a mesh sequence $(\Th)_{h>0}$, where for each $h>0$, the mesh $\Th$ is composed of nonempty disjoint open polyhedra with planar faces such that $\overline{\Omega}_0 = \bigcup_{T\in\Th} \overline{T}$. The mesh-size is $h = \max_{T\in\Th} h_T$, where $h_T$ stands for the diameter of the cell $T$. A closed subset $F$ of $\overline{\Omega}_0$ is called a mesh face if it is a subset with nonempty relative interior of some affine hyperplane $H_F$ and if either there exist two distinct mesh cells $T_-, \: T_+ \in \Th$ such that $F = \partial T_- \cap \partial T_+ \cap H_F$ (and $F$ is called an interface) or there exists one mesh cell $T\in\Th$ such that $F = \partial T \cap \Gamma \cap H_F$ (and $F$ is called a boundary face). The mesh faces are collected in the set $\Fh$ which is further partitioned into the subset $\Fhi$ which is the collection of the interfaces and the subset $\Fhb$ which is the collection of the boundary faces. We assume that the mesh is compatible with the partition of the boundary $\Gamma$ into $\Bd$ and $\Bn$, so that we can further split the set $\Fhb$ into the disjoint subsets $\Fhbd$ and $\Fhbn$ with obvious notation. For all $T \in \Th$, $\FT$ is the collection of the mesh faces that are subsets of $\partial T$ and $\nT$ is the unit outward normal to $T$. We assume that the mesh sequence $(\Th)_{h>0}$ is shape-regular in the sense specified in [@DiPEr:2015], i.e., there is a matching simplicial submesh of $\Th$ that belongs to a shape-regular family of simplicial meshes in the usual sense of Ciarlet [@Ciarlet1978] and such that each mesh cell $T \in \Th$ (resp., mesh face $F \in \Fh$) can be decomposed in a finite number of sub-cells (resp., sub-faces) which belong to only one mesh cell (resp., to only one mesh face or to the interior of a mesh cell) with uniformly comparable diameter.

Let $k \geq 1$ be a fixed polynomial degree and $l \in \lbrace k, k+1 \rbrace$. In each mesh cell $T\in\Th$, the local HHO unknowns consist of a pair $(\vT,\vdT)$, where the cell unknown $\vT\in \Pld(T; \Rd)$ is a vector-valued $d$-variate polynomial of degree at most $l$ in the mesh cell $T$, and $\vdT \in \PkF(\FT; \Rd) = \bigtimes_{F \in\FT} \PkF(F; \Rd)$ is a piecewise, vector-valued $(d-1)$-variate polynomial of degree at most $k$ on each face $F\in\FT$. We write more concisely that $$\label{local_dof}
(\vT,\vdT) \in \UklT := \Pld(T;\Rd) \times \PkF(\FT;\Rd).$$ We write the superscript $k$ first since $k$ is the value that determines the convergence rates of the approximation. The degrees of freedom are illustrated in Fig. \[fig\_HHO\_dofs\], where a dot indicates one degree of freedom (which is not necessarily computed as a point evaluation) and the geometric shape of the cell is only illustrative. We equip the space $\UklT$ with the following local discrete strain semi-norm: $$\label{eq:snorme}
\snorme[1,T]{( \vT, \vdT)}^2 := \normem[T]{\gradX \vT}^2 + \normev[\dT]{\gamma_{\dT}^{\frac12}(\vT-\vdT)}^2,$$ with the piecewise constant function $\gamma_{\dT}$ such that $\gamma_{\dT|F}=h_F^{-1}$ for all $F\in\FT$, where $h_F$ is the diameter of $F$. We notice that $\snorme[1,T]{( \vT, \vdT)}=0$ implies that $\vT$ is a constant and that $\vdT$ is the trace of $\vT$ on $\partial T$.

     

Local gradient reconstruction and stabilization {#sec:grad_rec}
-----------------------------------------------

The first key ingredient in the devising of the HHO method is a local gradient reconstruction in each mesh cell $T \in \Th$. This reconstruction is materialized by an operator $\GkT : \UklT \rightarrow \Pkd(T;\Rdd)$ mapping onto the space composed of $\Rdd$-valued polynomials in $T$. The main reason for reconstructing the gradient in a larger space than the space $\gradX \Pkpd(T; \Rd)$ (as for the linear elasticity problem [@DiPEr:2015] ) is that the reconstructed gradient of a test function acts against a discrete stress tensor which is not in gradient form, see [@DiPietro2018 Section 4] for further insight. For all $(\vT,\vdT) \in \UklT$, the reconstructed gradient $\GkT(\vT,\vdT) \in \Pkd(T;\Rdd)$ is obtained by solving the following local problem: For all $\matrice{\tau} \in \Pkd(T;\Rdd)$, $$\label{eq_reconstruction_grad}
\psm[T]{\GkT( \vT,\vdT)}{\matrice{\tau}} = \psm[T]{\gradX{\vT}}{\matrice{\tau}} + \psv[\dT]{\vdT - \vT{}_{|\dT}}{\matrice{\tau} \SCAL \nT}.$$ Solving this problem entails choosing a basis of the polynomial space $\Pkd(T;\Reel)$ and inverting the associated mass matrix for each component of the tensor $\GkT(\vT,\vdT)$. The second key ingredient in the HHO method is a local stabilization operator that enforces weakly the matching between the faces unknowns and the trace of the cell unknowns. Following [@DiPEL:2014; @DiPEr:2015], the stabilization operator $\SdTk : \PlF(\FT;\Rd) \rightarrow \PkF(\FT;\Rd)$ acts on the difference $\vecteur{\theta} = \vdT - \vT{}_{|\dT} \in \PlF(\FT;\Rd)$, and in the mixed-order case $l=k+1$ is such that, for all $\vecteur{\theta}\in \PkpF(\FT;\Rd)$, $$\label{eq:stabHDG}
\SdTk(\vecteur{\theta}) = \PikdTv\big(\vecteur{\theta}\big),$$ where $\PikdTv$ denotes the $L^2$-orthogonal projectors onto $\PkF(\FT;\Rd)$, and in the equal-order case $l=k$ is such that, for all $\vecteur{\theta}\in \PkF(\FT;\Rd)$, $$\label{eq:stabHHO}
\SdTk(\vecteur{\theta}) = \PikdTv\big(\vecteur{\theta}-(\matrice{I}-\PikTv)\DTkp(\vecteur{0},\vecteur{\theta})_{|\dT}\big),$$ where $\PikTv$ denotes the $L^2$-orthogonal projectors onto $\Pkd(T;\Rd)$. The local displacement reconstruction operator $\DTkp : \UklT \rightarrow \Poly_d^{k+1} (T; \Rd)$ is such that, for all $(\vT,\vdT) \in \UklT$, $\DTkp(\vT,\vdT) \in \Poly_d^{k+1}(T;\Rd)$ is obtained by solving the following local Neumann problem: For all $\vecteur{w} \in \Poly_d^{k+1} (T;\Rd)$, $$\label{eq_reconstruction_depl}
\psm[T]{\gradX{\DTkp(\vT,\vdT)}}{\gradX{\vecteur{w}}} = \psm[T]{\gradX{\vT}}{\gradX{\vecteur{w}}} + \psv[\dT]{\vdT - \vT{}_{|\dT}}{\gradX{\vecteur{w}} \SCAL \nT}.$$ together with the mean-value conditions $\int_T \DTkp(\vT,\vdT) \:dT= \int_T \vT \:dT$. Comparing with , one readily sees that $\gradX{\DTkp(\vT,\vdT)}$ is the $L^2$-orthogonal projection of $\GkT(\vT,\vdT)$ onto the subspace $\gradX \Poly^{k+1}_d(T;\Rd)$. Note also that the right-hand side of  can be rewritten as $\PikdTv(\vdT - \vT{}_{|\dT}-(\matrice{I}-\PikTv)\DTkp(\vT,\vdT)_{|\dT})$. Adapting [@DiPEr:2015 Lemma 4], it is straightforward to establish the following stability and boundedness properties (the proof is omitted for brevity).

\[lemma\_stability\_stab\] Let the gradient reconstruction operator be defined by  and the stabilization operator be defined by  or . Let $\gamma_{\dT}$ be defined below . Then, we have the following properties: *(i)* Boundedness: there exists $\alpha_\sharp< +\infty$, uniform w.r.t. $h$, such that, for all $T \in \Th$ and for $(\vT,\vdT) \in \UklT$, $$\bigg(\normem[T]{\GkT (\vT, \vdT)}^2 + \normev[\dT]{\gamma_{\dT}^{\frac12}\SdTk(\vdT - \vT{}_{|\dT})}^2\bigg)^{\frac12} \leq \alpha_\sharp \snorme[{1,T}]{(\vT, \vdT)}.$$ *(ii)* Stability: there exists $\alpha_\flat > 0$, uniform w.r.t. $h$, such that, for all $T \in \Th$ and all $(\vT,\vdT) \in \UklT$, $$\alpha_\flat \snorme[{1,T}]{(\vT, \vdT)} \leq \bigg(\normem[T]{\GkT (\vT, \vdT)}^2 + \normev[\dT]{\gamma_{\dT}^{\frac12}\SdTk(\vdT - \vT{}_{|\dT})}^2\bigg)^{\frac12}.$$

As shown in [@DiPEr:2015], the following important commuting property holds true: $$\label{eq:commuting}
\GkT(\vecteur{I}_{T,\dT}(\vecteur{v})) = \matrice{\Pi}^k_T(\gradX\vecteur{v}),\qquad \forall \vecteur{v}\in H^1(T;\Rd),$$ with the reduction operator $\vecteur{I}_{T,\dT} : H^1(T;\Rd)\rightarrow \UklT$ defined as $\vecteur{I}_{T,\dT}(\vecteur{v})=(\PilTv(\vecteur{v}),\PikdTv(\vecteur{v}{}_{|\dT}))$. Taking the trace in , we infer that $$\label{eq:trace_stab}
\trace\big(\GkT(\vecteur{I}_{T,\dT}(\vecteur{v}))\big) = \Pi_T^k(\divergenceX{\vecteur{v}}), \qquad \forall \vecteur{v}\in H^1(T;\Rd),$$ which is the key commuting property to prove robustness for quasi-incompressible linear elasticity, see [@DiPEr:2015]. This absence of volumetric locking is confirmed in the numerical experiments performed in Section \[sec::sec\_numexp\] in the nonlinear setting of finite elastoplasticity. Finally, proceeding as in [@DiPEr:2015 Thm. 8], one can show that for the linear elasticity problem and smooth solutions, the energy error converges as $h^{k+1}|\vecteur{u}|_{\vecteur{H}^{k+2}(\Omega_0)}$.

The stabilization operator is essential to prove the above mentioned convergence rates in the equal-order case for linear problems and smooth solutions on general meshes. In general, HDG methods use the stabilization operator $\SdTk(\vecteur{\theta}) = \vecteur{\theta}$ in the equal-order case which differs from the stabilization operator and allows one to show only that the energy error converges as $h^{k}|\vecteur{u}|_{\vecteur{H}^{k+1}(\Omega_0)}$ for linear problems and smooth solutions on general meshes. In the mixed-order case, the stabilization operator has been initially introduced in [@Lehrenfeld:10] and the same convergence rates as for the HHO method are obtained.

Global discrete problem
-----------------------

Let us now devise the global discrete problem. We set $\Pld(\Th;\Rd) :=\bigtimes_{T \in\Th} \Pld(T; \Rd)$ and $\PkF(\Fh;\Rd) := \bigtimes_{F \in\Fh} \PkF(F; \Rd)$. The global space of discrete HHO unknowns is defined as $$\label{eq:def_Ukh}
\Uklh := \Pld(\Th;\Rd) \times \PkF(\Fh;\Rd).$$ For an element $\vh \in \Uklh$, we use the generic notation $\vh = (\vTh,\vFh)$. For any mesh cell $T \in \Th$, we denote by $(\vT,\vdT)\in \UklT$ the local components of $\vh$ attached to the mesh cell $T$ and to the faces composing its boundary $\dT$, and for any mesh face $F\in\Fh$, we denote by $\vF$ the component of $\vh$ attached to the face $F$. The Dirichlet boundary condition on the displacement field can be enforced explicitly on the discrete unknowns attached to the boundary faces in $\Fhbd$. Letting $\PikFv$ denote the $L^2$-orthogonal projector onto $\PkF(F;\Rd)$, we set

$$\begin{aligned}
\Uklnhd &:= \left\lbrace (\vTh, \vFh ) \in \Uklh \: \vert \: \vF = \PikFv(\vu_{\textrm{D}}(t^n)), \; \forall F \in \Fhbd   \right \rbrace, \\
\Uklhz &:= \left\lbrace (\vTh, \vFh ) \in \Uklh \: \vert \: \vF = \vecteur{0}, \; \forall F \in \Fhbd   \right \rbrace.\end{aligned}$$

Note that the map $\vh\mapsto (\sum_{T\in\Th} \snorme[1,T]{( \vT, \vdT)}^2)^{\frac12}$ defines a norm on $\Uklhz$.

A key feature of the present HHO method is that the discrete generalized internal variables are computed only at the quadrature points in each mesh cell. We introduce for all $T \in \Th$, the quadrature points $\Qp_T= (\Qp_{T,j})_{1 \leq j \leq m_Q}$, with $\Qp_{T,j} \in T$ for all $1 \leq j \leq m_Q$, and the quadrature weights $\Wp_T= (\Wp_{T,j})_{1 \leq j \leq m_Q}$, with $\Wp_{T,j} \in \Reel$ for all $1 \leq j \leq m_Q$. We denote by $k_Q$ the order of the quadrature. Then, the discrete internal variables are sought in the space $$\IVGTh := \bigtimes_{T \in\Th}  \IVGS^{m_Q},$$ that is, for all $T\in\Th$, the internal variables attached to $T$ form a vector $\IVGT=(\IVGT(\Qp_{T,j}))_{1\le j\le m_Q}$ with $\IVGT(\Qp_{T,j})\in \IVGS$ for all $1\le j\le m_Q$.

We can now formulate the global discrete problem. We use the following notation for two tensor-valued functions $(\matrice{p}, \matrice{q})$ defined on $T$: $$\label{eq_prods}
\psmd[T]{\matrice{p}}{\matrice{q}} := \sum_{j=1}^{m_Q} \Wp_{T,j} \, \matrice{p}(\Qp_{T,j}) : \matrice{q}(\Qp_{T,j}).$$ We also need to consider the case where we know the tensor $\matrice{\tilde p}$ only at the quadrature nodes (we use a tilde to indicate this situation), i.e., we have $\matrice{\tilde p} = (\matrice{\tilde p}(\Qp_{T,j}))_{1\le j\le m_Q} \in (\Rdd)^{m_Q}$. In this case, we slightly abuse the notation by denoting again by $\psmd[T]{\matrice{\tilde p}}{\matrice{q}}$ the quantity equal to the right-hand side of . The discrete energy functional $\mathcal{E}_h^n: \Uklnhd \times \IVGTh \rightarrow \Reel$ is defined for any pseudo-time step $ 1 \leq n \leq N$ by $$\begin{aligned}
\label{eq::Eh}
\mathcal{E}_h^n\big((\vTh, \vFh), \IVGh\big) =& \sum_{T\in\Th}  \left\lbrace \psd[T]{\dQP{ \stressPotW}(\FkT(\vT, \vdT), \IVGT)}{1} - \psv[T]{\loadext^n}{\vT} \right\rbrace - \sum_{F\in\Fhbn} \psv[F]{ \Tn^n}{\vF} \nonumber \\
&+  \sum_{T\in\Th} \frac{\beta}{2} \normev[\dT]{\gamma_{\dT}^{\frac12}\SdTk(\vdT- \vT{}_{|\dT})}^2,\end{aligned}$$ for all $(\vTh, \vFh) \in \Uklnhd$ and $\IVGh \in \IVGTh$, with the local deformation gradient operator $\FkT:\UklT \rightarrow \Pkd(T;\Rdd)$ such that $\FkT(\vT, \vdT) := \matrice{I} + \GkT(\vT, \vdT)$. Moreover, in the second line of , the stabilization employs a weight of the form $\beta=2\mu \beta_0$ with $\beta_0>0$. In the original HHO method for linear elasticity [@DiPEr:2015], the choice $\beta_0=1$ is considered. In the present setting, the choice for $\beta_0$ is further discussed in Section \[sec:newton\] and in Section \[ss:stab\]. The global discrete problem consists in seeking for any pseudo-time step $ 1 \leq n \leq N$, a stationary point of the discrete energy functional: Find the pair of discrete displacements $(\uThn, \uFhn) \in \Uklnhd$ and the discrete internal variables $\IVGh^n \in \IVGTh$ such that, for all $(\dvTh, \dvFh) \in \Uklhz$, $$\begin{aligned}
\label{discrete_problem_ptv}
&\sum_{T\in\Th} \psmd[T]{\tilde{\PK}^{n}}{\GkT(\dvT,\dvdT)}
+  \sum_{T\in\Th} \beta \psv[\dT]{\gamma_{\dT}\SdTk(\udTn- \uTn{}_{|\dT})}{\SdTk(\dvdT - \dvT{}_{|\dT})} \nonumber \\
&= \sum_{T\in\Th}  \psv[T]{ \loadext^n}{\dvT} + \sum_{F\in\Fhbn} \psv[F]{\Tn^n}{\dvF},\end{aligned}$$ where for all $T \in \Th$ and all $1 \leq j \leq m_Q$, $$\label{discrete_problem_plast}
( \IVGT^n(\Qp_{T,j}), \tilde{\PK}^n(\Qp_{T,j}), \depmodulePKn(\Qp_{T,j}))  = 
\textrm{FINITE\_PLASTICITY}( \IVGT^{n-1}(\Qp_{T,j}), \FkT(\uT^{n-1},\udT^{n-1})(\Qp_{T,j}),\FkT(\uT^{n},\udT^{n})(\Qp_{T,j})),$$$$\begin{gathered}
\label{discrete_problem_plast}
( \IVGT^n(\Qp_{T,j}), \tilde{\PK}^n(\Qp_{T,j}), \depmodulePKn(\Qp_{T,j}))  = \\
\textrm{FINITE\_PLASTICITY}( \IVGT^{n-1}(\Qp_{T,j}), \FkT(\uT^{n-1},\udT^{n-1})(\Qp_{T,j}),\FkT(\uT^{n},\udT^{n})(\Qp_{T,j})),\end{gathered}$$with $(\uTh^{n-1}, \uFh^{n-1}) \in \Uklpnhd$ and $\IVGh^{n-1} \in \IVGTh$ given either from the previous pseudo-time step or the initial condition.

Discrete principle of virtual work
----------------------------------

The discrete problem  expresses the principle of virtual work at the global level, and adapting the ideas introduced in [@CoDPE:2016] (see also [@AbErPi2018a; @BoDPS:2017]), it is possible to infer a local principle of virtual work in terms of face-based discrete tractions that comply with the law of action and reaction.

Let $\SdTks : \PkF(\FT;\Rd) \rightarrow \PkF(\FT;\Rd)$ be the adjoint operator of $\SdTk$ with respect to the $L^2(\dT;\Rd)$-inner product so that we have $\psv[\dT]{\gamma_{\dT}\SdTk(\vecteur{\theta})}{\SdTk(\vecteur{\zeta})} = \psv[\dT]{\SdTks(\gamma_{\dT}\SdTk(\vecteur{\theta}))}{\vecteur{\zeta}}$ (recall that the weight $\gamma_{\dT}$ is piecewise constant on $\dT$). Let $\PikTmd: (\Rdd)^{m_Q} \rightarrow \Pkd(T;\Rdd )$ denote the $L^2_Q$-orthogonal projector such that for all $\tilde{\matrice{p}} \in (\Rdd)^{m_Q}$, $\psm[T]{\PikTmd(\tilde{\matrice{p}})}{\matrice{q}} = \psmd[T]{\tilde{\matrice{p}}}{\matrice{q}}$ for all $\matrice{q} \in \Pkd(T; \Rdd)$. Finally, for any pseudo-time step $1\le n\le N$ and all $T\in\Th$, let us define the discrete traction: $$\label{eq:traction_stab}
\vecteur{T}_{T}^{n} := \PikTmd(\dPK_T^{n})\SCAL\nT
+ \beta \SdTks(\gamma_{\dT}\SdTk(\udTn- \uTn{}_{|\dT}))\in \PkF(\FT; \Rd),$$ where $\dPK_T^n =  (\dPK_{T}^n(\Qp_{T,j}))_{1 \leq j \leq m_Q} \in( \Rdd)^{m_Q}$ and $(\uTn, \udTn) \in \UklT$.

\[lem:equil\_sHHO\] Assume that $k_Q\ge 2k$. Then, for any pseudo-time step $1 \leq n \leq N$, the following local principle of virtual work holds true for all $T\in\Th$: $$\label{eq:pwk}
\psmd[T]{\dPK^n_T}{\gradX{\dvT}}-
\psv[\dT]{\vecteur{T}_{T}^n}{\dvT{}_{|\dT}} = \psv[T]{ \loadext^n}{\dvT},
\qquad \forall \dvT \in \Pld(T;\Rd),$$ where the discrete tractions $\vecteur{T}_{T}^n$ defined by  satisfy the following law of action and reaction for all $F\in\Fhi\cup\Fhbn$:

\[eq:balance\] $$\begin{aligned}
{2}
&\vecteur{T}_{T_-|F}^n + \vecteur{T}_{T_+|F}^n = \vecteur{0},
&\quad&\text{if $F\in\Fhi$ with $F=\partial T_- \cap \partial T_+ \cap H_F$},\\
&\vecteur{T}_{T|F}^n  = \PikFv(\Tn^n),&\quad&\text{if $F\in\Fhbn$ with $F=\partial T\cap \Bn\cap H_F$}.\end{aligned}$$

Nonlinear solver and implementation {#sec:newton}
-----------------------------------

The nonlinear problem - arising at any pseudo-time step $1 \leq n \leq N$ is solved using a Newton’s method. Given $(\uTh^{n-1}, \uFh^{n-1}) \in \Uklpnhd$ and $\IVGh^{n-1} \in \IVGTh$ from the previous pseudo-time step or the initial condition, the Newton’s method is initialized by setting $(\uTh^{n,0}, \uFh^{n,0})=(\uTh^{n-1}, \uFh^{n-1})$, up to the update of the Dirichlet condition, and $\IVGh^{n,0}= \IVGh^{n-1}$. Then, at each Newton’s step $i\ge0$, one computes the incremental displacement $(\duTh^{n,i}, \duFh^{n,i}) \in \Uklhz$ and updates the discrete displacement as $(\uTh^{n,i+1} \uFh^{n,i+1}) = (\uTh^{n,i}, \uFh^{n,i})+(\duTh^{n,i}, \duFh^{n,i})$. The linear system of equations to be solved is $$\begin{aligned}
\label{eq_stiffness_matrix}
& \hphantom{+} \sum_{T\in\Th} \psmd[T]{\depmodulePKni: \GkT(\duT^{n,i},\dudT^{n,i})}{\GkT(\dvT,\dvdT)} \nonumber \\ 
&+  \sum_{T\in\Th} \beta \psv[\dT]{\gamma_{\dT}\SdTk(\dudT^{n,i} - 
\duT^{n,i}{}_{|\dT})}{\SdTk(\dvdT - \dvT{}_{|\dT})} \nonumber \\ 
&= -R_h^{n,i}(\dvTh,\dvFh),\end{aligned}$$ for all $(\dvT,\dvdT)\in\Uklhz$, where for all $T \in \Th$ and all $1 \leq j \leq m_Q$, $$(\IVGT^{n,i}(\Qp_{T,j}), \dPK^{n,i}(\Qp_{T,j}), \depmodulePKni(\Qp_{T,j})) = \textrm{FINITE\_PLASTICITY}(\IVG_{T,j}^{n-1},\Fdef_{T,j}^{n-1},\Fdef_{T,j}^{n,i}),$$ with $\IVG_{T,j}^{n-1} := \IVGT^{n-1}(\Qp_{T,j})$, $\Fdef_{T,j}^{n,i} := \FkT(\uT^{n,i},\udT^{n,i})(\Qp_{T,j})$, $\Fdef_{T,j}^{n-1} := \FkT(\uT^{n-1},\udT^{n-1})(\Qp_{T,j})$, and the residual term $$\begin{aligned}
R_h^{n,i}(\dvTh, \dvFh) 
={}& \sum_{T\in\Th} \psmd[T]{\dPK^{n,i}}{\GkT(\dvT,\dvdT)} - \sum_{T\in\Th}  \psv[T]{ \loadext^n}{\dvT} - \sum_{F\in\Fhbn} \psv[F]{\Tn^n}{\dvF}  \nonumber \\
&+  \sum_{T\in\Th} \beta \psv[\dT]{\gamma_{\dT}\SdTk(\udT^{n,i} - \uT^{n,i}{}_{|\dT})}{\SdTk(\dvdT - \dvT{}_{|\dT})}.\end{aligned}$$ The assembling of the stiffness matrix resulting from the left-hand side of  is local (and thus fully parallelizable). The discrete internal variables $\IVGh^{n} \in \IVGTh$ are updated at the end of each pseudo-time step. Moreover, since the consistent elastoplastic tangent modulus $\epmodulePK$ has major symmetries, its eigenvalues are real. Let us set $\theta_{\Th, Q} := \min_{(T,j) \in \Th\times \{1,\ldots,m_Q\}} \theta^{\min}(\depmodulePK(\Qp_{T,j}))$, where $\theta^{\min}(\tenseur{4}{{\mathbb{M}}})$ denotes the smallest eigenvalue of a symmetric fourth-order tensor $\tenseur{4}{{\mathbb{M}}}$. The following result shows that the linear system  arising at each Newton’s step is coercive under the simple choice $\beta_0>0$ on the stabilization parameter if $\depmodulePK$ is positive-definite. Notice that strain-hardening plasticity is not a sufficient condition for positive-definiteness of $\depmodulePK$ (only for $\depmoduleT$) since in finite elastoplasticity, geometrical nonlinearities also exist.

\[th::coer\_newton\] Assume that $k_Q\ge 2k$ and that all the quadrature weights are positive. Moreover, assume that the consistent elastoplastic tangent modulus $\epmodulePK$ is positive-definite, i.e., $\theta_{\Th, Q} >0$. Then, the linear system  in each Newton’s step is coercive for all $\beta_0 > 0$, i.e., there exists $C_{\rm ell}>0$, independent of $h$, such that for all $(\vTh, \vFh) \in \Uklhz$, $$\begin{gathered}
\label{eq_stable}
\sum_{T\in\Th} \psmd[T]{\depmodulePK: \GkT(\vT,\vdT)}{\GkT(\vT,\vdT)}+  \sum_{T\in\Th} \beta \normev[\dT]{\gamma_{\dT}^{\frac12}\SdTk(\vdT- \vT{}_{|\dT})}^2 \\
\geq C_{\rm ell} \min\left(\beta_0, \frac{\theta_{\Th, Q}}{2\mu} \right) 2\mu \sum_{T\in\Th} \snorme[{1,T}]{(\vT, \vdT)}^2.\end{gathered}$$

The proof follows directly from [@AbErPi2018a Theorem 6]. Note that Theorem \[th::coer\_newton\] states that one iteration of the Newton’s process is stable under a positive-definiteness assumption but does not state that a solution to the whole Newton’s process exists. The existence of such a solution has been showed in the context of nonlinear elliptic equations for dG methods [@Yadav2013]. Moreover, if the consistent elastoplastic tangent modulus $\epmodulePK$ is no longer positive-definite for at least one quadrature point which is a possibility in finite elastoplasticity since geometrical nonlinearities also exist, then Theorem \[th::coer\_newton\] is no longer valid (even if $\depmoduleT$ remains positive-definite). Moreover, a reasonable choice of the stabilization parameter appears to be $\beta_0 \geq \max(1,\frac{\theta_{\Th, Q}}{2\mu})$ because $\beta_0 = 1$ is a natural choice for linear elasticity [@DiPEr:2015] and the choice $\beta_0 \geq  \frac{\theta_{\Th, Q}}{2\mu}$ allows one to adjust the stabilization parameter if the evolution is plastic. We investigate numerically the choice of $\beta_0$ in Section \[ss:stab\].

From a numerical point of view, as is classical with HHO methods [@DiPEr:2015; @DiPEL:2014], and more generally with hybrid approximation methods, the cell unknowns $\duT^{n,i}$ in  can be eliminated locally by using a static condensation (or Schur complement) technique. This allows one to reduce  to a linear system in terms of the face unknowns only. The reduced system is of size $N_{\Fh} \times  d{k+d-1\choose d-1}$, where $N_{\Fh}$ denotes the number of mesh faces (unknowns attached to Dirichlet boundary faces can be eliminated by enforcing the boundary condition explicitly). The implementation of HHO methods is realized using the open-source library `DiSk++` which provides generic programming tools for the implementation of HHO methods and is available online[^1]. We refer the reader to [@CicDE:2018] and [@AbErPi2018a Section 3.6] for further aspects about the implementation.

Numerical examples {#sec::sec_numexp}
==================

The goal of this section is to evaluate the proposed HHO method on two- and three-dimensional benchmarks from the literature: (*i*) a necking of a 2D rectangular bar subjected to uniaxial extension, (*ii*) a Cook’s membrane subjected to bending, (*iii*) a torsion of a square-section bar, and (*iv*) a quasi-incompressible sphere under internal pressure. We compare our results to the analytical solution whenever available or to numerical results obtained using the industrial open-source FEM software [@CodeAster]. In this case, we consider a linear, resp. quadratic, cG formulation, referred to as Q1, resp. T2 or Q2 when full integration is used, or, Q2\_RI when reduced integration is used, depending on the mesh, and a three-field mixed formulation in which the unknowns are the displacement, the pressure and the volumetric strain fields referred to as UPG [@AlAkhrass2014]; in the UPG method, the displacement field is quadratic, whereas both the pressure and the volumetric strain fields are linear. The conforming Q1, T2 and Q2 methods with full integration, contrary to Q2\_RI method with reduced integration in most of the situations, are known to present volumetric locking due to plastic incompressibility, whereas the UPG method is known to be robust but costly. Numerical results obtained using the UPG method are used as a reference solution whenever an analytical solution is not available.

The nonlinear isotropic plasticity model with a von Mises yield criterion described in Section \[ss:nl\_hardening\] is used for the test cases. For the first three test cases, strain-hardening plasticity is considered with the following material parameters: Young modulus $E =206.9~\GPa$, Poisson ratio $\nu=0.29$, hardening parameter $H = 129.2~\MPa$, initial yield stress $\sigma_{y,0} = 450~\MPa$, infinite yield stress $\sigma_{y,\infty} = 715~\MPa$, and saturation parameter $\delta = 16.93$. For the fourth case, perfect plasticity is considered with the following material parameters: Young modulus $E =28.85~\MPa$, Poisson ratio $\nu=0.499$, hardening parameter $H =0~\MPa$, initial and infinite yield stresses $\sigma_{y,0} = \sigma_{y,\infty} = 6~\MPa$, and saturation parameter $\delta = 0$. Moreover, for the two-dimensional test cases (*i*) and (*ii*), we assume additionally a plane strain condition. In the numerical experiments reported in this section, the stabilization parameter is taken to be $\beta = 2 \mu$ ($\beta_0=1$), and all the quadratures use positive weights. In particular, for the HHO method, we employ a quadrature of order $k_Q=2k$ for the behavior cell integration. We employ the notation HHO($k$, $l$) when using face polynomials of order $k$ and cell polynomials of order $l$.

In Section \[ss::num\_inv\], we perform further numerical investigations to test other aspects of HHO methods such as the support of general meshes with possibly non-conforming interfaces, the possibility of considering the lowest-order case $k=0$, and the dependence on the stabilization parameter $\beta$.

Necking of a 2D rectangular bar {#ss::bar}
-------------------------------

In this first benchmark, we consider a 2D rectangular bar with an initial imperfection. The bar is subjected to uniaxial extension. This example has been studied previously by many authors as a necking problem [@Elguedj2014; @Simo1992a; @SouzaNeto2005; @Saracibar2006; @Wriggers2017a] and can be used to test the robustness of the different methods. The bar has a length of $53.334~\mm$ and a variable width from an initial width value of $12.826~\mm$ at the top to a width of $12.595~\mm$ at the center of the bar to create a geometric imperfection. A vertical displacement $u_y = 5~\mm$ is imposed at both ends, as shown in Fig. \[fig::necking\_geom\]. For symmetry reasons, only one-quarter of the bar is discretized, and the mesh is composed of 400 quadrangles, see Fig. \[fig::necking\_mesh\]. The load-displacement curve is plotted in Fig. \[fig::necking\_load\]. We observe that except for Q1, all the other methods give very similar results. Moreover, the equivalent plastic strain $p$, respectively the trace of the Cauchy stress tensor $\stress$, are shown in Fig. \[fig::necking\_p\], resp. in Fig. \[fig::necking\_trace\], at the quadrature points on the final configuration. A sign of locking is the presence of strong oscillations in the trace of the Cauchy stress tensor $\stress$. We notice that the cG formulations Q1 and Q2 lock, contrary to the HHO, Q2\_RI, and UPG methods which deliver similar results. We remark however that the results for HHO(1;1), HHO(1;2), and Q2\_RI are slightly less smooth than for HHO(2;2), HHO(2;3), and UPG. The reason is that on a fixed mesh, the three former methods have less quadrature points than the three latter ones, see Table \[tab::necking\] (HHO(2;2), HHO(2;3), and UPG have the same number of quadrature points). Therefore, the stress is evaluated using less points in HHO(1;1), HHO(1;2), and Q2\_RI. It is sufficient to refine the mesh or to increase the order of the quadrature by two in HHO(1;1) and HHO(1;2) to retrieve similar results to those for the three other methods (not shown for brevity).

       Method        Q1     Q2    [[Q2\_RI]{}]{}   UPG    HHO(1;1)   HHO(1;2)   HHO(2;2)   HHO(2;3)
  ---------------- ------ ------ ---------------- ------ ---------- ---------- ---------- ----------
   Number of dofs   884    2566    [[2566]{}]{}    3450     3364       3364       5046       5046
   Number of QPs    1600   3600    [[1600]{}]{}    3600     1600       1600       3600       3600

  : Necking of a 2D rectangular bar: number of globally coupled degrees of freedom (dofs) and quadrature points (QPs) for the different methods.[]{data-label="tab::necking"}

   

     

     

     

     

Cook’s membrane {#ss::cook}
---------------

We consider the Cook’s membrane problem which is a well known bending-dominated test case [@Simo1992a; @AlAkhrass2014; @Elguedj2008]. It consists of a tapered panel, clamped on one side, and subjected to a total vertical load $F_y = 5~\kN$ applied uniformly along all the opposite side, as shown in Fig. \[fig::cook\_geom\]. The simulation is performed on a sequence of refined quadrangular meshes such that each side contains $2^N$ edges with $0\leq N \leq 6$. The vertical displacement of the point $A$ versus the number of degrees of freedom is plotted in Fig. \[fig::cook\_depl\] for the different methods. As expected when comparing the number of degrees of freedom, the linear cG formulation Q1 has the slower convergence, HHO(1;2), Q2\_RI, and UPG converge slightly faster than HHO(1;1), Q2, whereas HHO(2;2) and HHO(2;3) outperform all the other methods and give almost the same results. Moreover, we show in Fig. \[fig::cook\_trace\] the trace of the Cauchy stress tensor $\stress$ at the quadrature points on the final configuration. The cG formulations Q1 and Q2 present oscillations that confirm the presence of volumetric locking, contrary to the HHO, Q2\_RI, and UPG methods which deliver similar and smooth results (even if the cG formulations Q1 and Q2 present volumetric locking in terms of stress, they converge in terms of displacement with mesh refinement). However, if we compare the trace of the Cauchy stress tensor $\stress$ for HHO(1;1) and HHO(1;2), we remark that the trace is slightly smoother near the upper-right corner for HHO(1;2) than for HHO(1;1). This can be explained by the presence of non-physical vertical localization bands of plastic deformations for HHO(1;1) and not for HHO(1;2). Localization bands constitute a well-known problem when the plasticity model is local. Computational practice with cG approximations indicates that increasing the order of the finite elements mitigates this issue. The same effect is observed here by increasing the degree of the cell unknowns (further numerical investigations are performed in Sect. \[ss::poly\_mesh\] for HHO methods). An alternative is to use a non-local plasticity model [@McBride2009].

 

   

   

 

Torsion of a square-section bar {#ss:torsion}
-------------------------------

This third benchmark [@Hudobivnik2018] allows one to test the robustness of HHO methods under large torsion. The bar has a square-section of length $L=1~\mm$ and a height of $H=5~\mm$ along the $z$-direction. The bottom end is clamped and the top end is subjected to a planar rotation of an angle $\Theta $ around its center along the $z$-direction and remains plane (the displacement of the top end along the $z$-direction is blocked), see Fig. \[fig::bar\_geom\]. The mesh is composed of 1900 hexahedra, see Fig. \[fig::bar\_mesh\]. The equivalent plastic strain $p$ is plotted in Fig. \[fig::bar\_p\] for HHO(1;1) and for different rotation angles $\Theta$. There is no sign of localization of the plastic deformations even for large rotations whatever the HHO variant is used. Moreover, the trace of the Cauchy stress tensor $\stress$ is plotted on the final configuration for $\Theta=360^{\circ}$ and for the Q2, UPG, and HHO methods in Fig. \[fig::bar\_trace\]. As expected, there is no sign of volumetric locking for the HHO and UPG methods which give similar results contrary to Q2. The small oscillations observed at both ends are due to the imposed conditions on the displacement.

 

         

         

Quasi-incompressible sphere under internal pressure {#ss::sphere}
---------------------------------------------------

This last benchmark [@AlAkhrass2014] consists of a quasi-incompressible sphere under internal pressure for which an analytical solution is known when the entire sphere has reached a plastic state. This benchmark is particularly challenging compared to the previous ones since we consider here perfect plasticity. The sphere has an inner radius $R_{in}=0.8~\mm$ and an outer radius $R_{out}=1~\mm$. An internal radial pressure $P$ is imposed. For symmetry reasons, only one-eighth of the sphere is discretized, and the mesh is composed of 1580 tetrahedra, see Fig. \[fig::sphere\_mesh\]. The simulation is performed until the limit load corresponding to an internal pressure $P_{lim} \simeq 2.54~\MPa$ is reached.

 

The equivalent plastic strain $p$ is plotted for HHO(1;2) in Fig. \[fig::sphere\_p\], and the trace of the Cauchy stress tensor $\stress$ is compared for HHO, UPG and T2 methods in Fig. \[fig::sphere\_trace\] at all the quadrature points on the final configuration for the limit load. We notice that the quadratic element T2 locks, whereas HHO and UPG do not present any sign of locking and produce results that are very close to the analytical solution. However, the trace of the Cauchy stress tensor $\stress$ is slightly more dispersed around the analytical solution for HHO(2;2) and HHO(2;3) than for HHO(1;1) and HHO(1;2) near the outer boundary. For this test case, we do not expect that HHO(2;2) and HHO(2;3) will deliver more accurate solutions than HHO(1;1) and HHO(1;2) since the geometry is discretized using tetrahedra with planar faces.

We next investigate the influence of the quadrature order $k_Q$ on the accuracy of the solution. The trace of the Cauchy stress tensor $\stress$ is compared for HHO(1;1), HHO(2;2), and UPG methods in Fig. \[fig::sphere\_trace\_kq\] at all the quadrature points on the final configuration for the limit load, and for a quadrature order $k_Q$ higher than the one employed in Fig. \[fig::sphere\_trace\] (HHO(1;2) and HHO(2;3) give similar results and are not shown for brevity). We remark that when we increase the quadrature order, UPG locks for quasi-incompressible finite deformations, whereas HHO does not lock, and the results are (only) a bit more dispersed around the analytical solution. Moreover, HHO(2;2) is less sensitive than HHO(1;1) to the choice of the quadrature order $k_Q$. Note that this problem is not present for HHO methods with small deformations [@AbErPi2018a]. Furthermore, this sensitivity to the quadrature order seems to be absent for finite deformations when the elastic deformations are compressible (the plastic deformations are still incompressible). To illustrate this claim, we perform the same simulations as before but for a compressible material. The Poisson ratio is taken now as $\nu=0.3$ (recall that we used $\nu=0.499$ in the quasi-incompressible case) whereas the other material parameters are unchanged. Unfortunately, an analytical solution is no longer available in the compressible case. We compare again the trace of the Cauchy stress tensor $\stress$ for HHO(1;1), HHO(2;2), and UPG methods in Fig. \[fig::sphere\_trace\_comp\] at all the quadrature points on the final configuration and for different quadrature orders $k_Q$. We observe a quite marginal dependence on the quadrature order for HHO methods (as in the quasi-incompressible case); whereas the UPG method still locks if the order of the quadrature is increased. Moreover, in the compressible case, HHO(2;2) gives a more accurate solution than HHO(1;1).\

 

 

 

 

 

Summary of the above results
----------------------------

The proposed HHO method has been tested successfully on four benchmarks in two and three dimensions. A first conclusion is that the proposed HHO method is robust for large elastoplastic deformations and is locking-free as mixed methods but without the need to introduce additional globally coupled unknowns. HHO(2;2) and HHO(2;3) give generally more accurate results both for the displacement and the Cauchy stress tensor than HHO(1;1), HHO(1;2), and UPG on a fixed mesh (cG methods lock). Moreover, contrary to the UPG method, HHO methods are not very sensitive to the choice of the quadrature order (particularly for $k=2$). Finally, HHO(1;1) appears to be more prone to the localization of the plastic deformations, contrary to the other variants HHO(1;2), HHO(2;2), and HHO(2;3).

Further numerical investigations {#ss::num_inv}
================================

In this section, we perform further numerical investigations to test other capacities of HHO methods such as the support of general meshes with possibly non-conforming interfaces, the possibility of considering the lowest-order case $k=0$, and the dependence on the stabilization parameter $\beta$.

Polygonal meshes {#ss::poly_mesh}
----------------

In the previous sections, the proposed HHO method has been tested on simplicial and hexahedral meshes so as to be able to compare it to the UPG method which only supports this type of meshes. Our goal is now to illustrate that the HHO method supports general meshes with possibly non-matching interfaces. For our test cases, the polygonal meshes are generated from quadrangular meshes by removing the common face for some pairs of neighbouring cells and then merging the two cells in question (about 30% of the cells are merged) thereby producing non-matching interfaces materialized by hanging nodes for a significant portion of the mesh cells. We consider the Cook’s membrane problem from Section \[ss::cook\] and we use a mesh composed of 719 polygonal cells including quads, pentagons (quads with one hanging node) and hexagons (quads with two hanging nodes). The trace of the Cauchy stress tensor $\stress$ is shown in Fig. \[fig::cook\_poly\] at the quadrature points on the reference configuration; and we compare the results with a reference solution computed with HHO(2;3) on a $32\times 32$ quadrangular mesh. The results agree very well except for HHO(1;1) where the trace is not smooth due to the localization of the plastic deformations (as for the quadrangular mesh, see Section \[ss::cook\]). A reason for the localization of the plastic deformations is the loss of coercivity of the consistent elastoplastic tangent modulus $\depmodulePK$. The evolution of the magnitude of the smallest eigenvalue $\theta_{\Th, Q}$ of $\depmodulePK$ during the loading is plotted in Fig. \[fig::cook\_evol\_vp\]. We remark that the magnitude of $\theta_{\Th, Q}$ decreases quickly when the plastic evolution begins (around $F_y = 1~\kN$), then $\theta_{\Th, Q}$ continues to decrease more slowly and remains positive for all HHO methods except for HHO(1;1) where it decreases more quickly and becomes negative (so that Theorem \[th::coer\_newton\] is no longer valid). As mentioned above, this loss of positive-definiteness of $\depmodulePK$ for HHO(1;1) can explain the presence of nonphysical plastic localization. Moreover, the distribution of the smallest eigenvalue $\theta^{\min}(\depmodulePK)$ of $\depmodulePK$ at the quadrature points at the end of the loading $F_y = 5~\kN $ for the different HHO methods is summarized in Table \[tab::cook\_vp\]. Only HHO(1;1) has negative eigenvalues. We remark that the number of quadrature points where $\depmodulePK$ has negative eigenvalues for HHO(1;1) or close to 0 (&lt;0.5) for the others HHO methods is small compared to the total number of quadrature points. This confirms that the loss of positive-definiteness and the presence of plastic localization are local and confined to few quadrature points. Note that for a total vertical load $F_y > 6.1~\kN $, $\theta_{\Th, Q}$ is negative for all HHO methods.

   

 

 

Lowest-order variant
--------------------

The main reason to take $k \geq 1$ in the HHO method applied to the linear elasticity problem is that the rigid-body motions $\vecteur{RM}(T)$ are then a subset of $\vecteur{U}_{T}^{1,1}$. The lowest-order case $k=0$ and $l=1$ is interesting since there are only $d$ unknowns per face, i.e, two in 2D and three in 3D; and we could expect that the energy error, resp. the $L^2$-error, converges as $h|\vecteur{u}|_{\vecteur{H}^{2}(\Omega_0)}$, resp. as $h^2|\vecteur{u}|_{\vecteur{H}^{2}(\Omega_0)}$, for the linear elasticity problem. The difficulty with this lowest-order case is to deal with the rigid-body motions on the faces since unfortunately $\vecteur{RM}(T){}_{|\dT} \nsubseteq \Poly^{0}_{d-1}(\FT; \Rd)$. Therefore, at the theoretical level, it is not clear that Lemma \[lemma\_stability\_stab\] still holds true. Nevertheless, we observed numerically that for the linear elasticity problem, the energy error, resp. the $L^2$-error, converges as $h|\vecteur{u}|_{\vecteur{H}^{2}(\Omega_0)}$, resp. as $h^2|\vecteur{u}|_{\vecteur{H}^{2}(\Omega_0)}$, (the expected optimal rates) if all the cells have at least $2d$ faces, i.e, four faces in 2D and six faces in 3D. This observation seems to be confirmed for small elastoplasticity. However, the conclusions are less clear for finite elastoplasticity. The equivalent plastic strain $p$ and the trace of the Cauchy stress tensor $\stress$ are plotted in Fig. \[fig::low-order\] at the quadrature points on the final configuration for the necking of a rectangular bar (see Section \[ss::bar\]) approximated using the HHO(0;1) variant. We observe the absence of volumetric locking and that the results are close to those obtained for $k\geq 1$ (see Fig. \[fig::necking\_p\] and Fig. \[fig::necking\_trace\]). However, for the Cook’s membrane problem (see Section \[ss::cook\]), the displacement is not correct (not shown for brevity).

In[@Fu2015], a theoretical study of an HDG method for the linear elasticity problem in the equal-order case $k=0$ is performed where the main difference with the present HHO method is the stabilization weight which is of the form $O(1)$ and no longer $O(h^{-1})$ as here. In this case, for the linear elasticity problem, the the energy error, the $L^2$-error, and the stress error converge as $h|\vecteur{u}|_{\vecteur{H}^{2}(\Omega_0)}$, $h|\vecteur{u}|_{\vecteur{H}^{2}(\Omega_0)}$, as $h^{\frac12}|\vecteur{u}|_{\vecteur{H}^{2}(\Omega_0)}$, respectively, on general meshes. Moreover, recent numerical results still for the linear elasticity problem [@Sevilla2019] indicate that the stress error can converge as $h|\vecteur{u}|_{\vecteur{H}^{2}(\Omega_0)}$. Nevertheless, the $L^2$-error converges slower for this HDG variant than what we could expect for the lowest-order HHO method.

 

Influence of the stabilization parameter {#ss:stab}
----------------------------------------

To evaluate the influence of the stabilization parameter $\beta_0$, we compare the total number of Newton’s iterations needed to solve the nonlinear problem versus the magnitude of the stabilization parameter $\beta_0$. The Newton’s iterations are stopped under the relative criterion $\norme[\Th]{R_h(\dvTh,\dvFh)} \leq 10^{-6} \norme[\Th]{F_{\textrm{int}}(\dvTh,\dvFh)}$ where $F_{\textrm{int}} $ are the internal forces. We perform this comparison on two of the previous benchmarks, the Cook’s membrane problem (see Section \[ss::cook\]) and the quasi-incompressible sphere under internal pressure (see Section \[ss::sphere\]). In Fig. \[fig:cook\_beta\], we report the total number of Newton’s iterations for the Cook’s membrane problem with strain-hardening plasticity. We use a $32 \times 32$ quadrangular mesh, and 15 load increments of equal size are considered. On the one hand, we remark that the different HHO variants need almost the same total number of Newton’s iterations (around 78 compared to 75 for UPG) if $\beta_0 \geq 0.1$; On the other hand, if $\beta_0  < 0.01$, the Newton’s method stops converging whatever the HHO variant and the number of load increments. For the quasi-incompressible sphere under internal pressure, the pressure is applied in 15 increments of equal size. Recall that this experiment is particularly challenging since we are considering here perfect plasticity for which the stability result from Theorem \[th::coer\_newton\] is not applicable. In Fig. \[fig:sphere\_beta\], we plot the total number of Newton’s iterations to perform the simulation. On the one hand, if $\beta_0 \geq 10$, all the HHO variants need almost to the same total number of Newton’s iterations (around 57 compared to 55 for UPG). On the other hand, if $\beta_0 \leq 1$, the HHO variants with $k =2$ need more Newton’s iterations than the HHO variants with $k=1$. As previously, if $\beta_0 < 0.1$, the Newton’s method stops converging.

A first conclusion is therefore that the proposed HHO methods are stable for a large range of values of the stabilization parameter $\beta_0$. A second conclusion is that it seems reasonable to take $\beta_0 \in [1, 100]$ since the number of Newton’s iterations is lower and close to the value for UPG, and the condition number does not increase too much. Note that for extremely large values of $\beta_0$, HHO methods reduce to $H^1$-conforming methods due to the matching of the face unknowns with the trace of the cell unknowns, and volumetric-locking can appear (not shown here for brevity).

 

Conclusions
===========

We have devised and evaluated numerically a Hybrid High-Order (HHO) method to approximate finite elastoplastic deformations within a logarithmic strain framework. This framework allows one to re-use behavior laws developed originally for small deformations in the context of finite deformations.The HHO method exhibits a robust behavior for strain-hardening plasticity as well as for perfect plasticity, and produces accurate solutions with a moderate number of degrees of freedom for various benchmarks from the literature. In particular, as mixed methods, the HHO method avoids volumetric locking due to plastic incompressiblity, but with less unknowns than mixed methods for the same accuracy. Moreover, the HHO method supports general meshes with non-matching interfaces.

This work can be pursued in several directions. One could use a non-local plasticity model, as for example a strain-gradient plasticity model, to take into account scale-dependent effects [@McBride2009] and possibly prevent unphysical localization of the plastic deformations. Furthermore, the extension of the present HHO method to contact and friction problems is the subject of ongoing work.

[^1]: `https://github.com/wareHHOuse/diskpp`
