# A first order energy-dissipative momentum conserving scheme for elasto-plasticity using the variational updates formulation

Ludovic Noels, Laurent Stainier, Jean-Philippe Ponthot

► To cite this version:

Ludovic Noels, Laurent Stainier, Jean-Philippe Ponthot. A first order energy-dissipative momentum conserving scheme for elasto-plasticity using the variational updates formulation. Computer Methods in Applied Mechanics and Engineering, 2008, 197 (6-8), pp.706-726. 10.1016/j.cma.2007.09.008. hal-01004830

HAL Id: hal-01004830

https://hal.science/hal-01004830

Submitted on 19 Oct 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A first-order energy-dissipative momentum-conserving scheme for elasto-plasticity using the variational updates formulation

L. Noels¹, L. Stainier², J.-P. Ponthot

University of Liège, Continuum Mechanics & Thermomechanics, Chemin des Chevreuil 1, B-4000 Liège, Belgium

In a previous paper [L. Noels, L. Stainier, J.-P. Ponthot, An energy momentum conserving algorithm using the variational formulation of visco-plastic updates, Int. J. Numer. Methods Engrg. 65 (2006) 904–942] the authors demonstrated the efficiency of the variational formulation of elasto-plastic updates to develop energy–momentum conserving time integration algorithms. Indeed, within such a framework, the stress tensor always derives from an incremental potential, even when plastic behavior is considered. Therefore the verification of the conservation of energy in the non-linear range can easily be demonstrated: the sum of the reversible stored energy and irreversible dissipated energy exactly corresponds to the work of the external forces applied to the structure. Although this formulation was shown to be accurate and robust, the introduction of numerical dissipation for high-frequency numerical modes can be necessary to simulate complex phenomena. In this work, we propose a modification of the variational updates framework to introduce this numerical property, leading to a new energy-dissipative momentum-conserving time-integration algorithm for elasto-plasticity.

**Keywords:** Energy dissipation; Momentum conservation; Dynamics; Variational formulation; Elasto-plasticity; Finite-elements

## 1. Introduction

The time-integration of a finite-element space-discretization is commonly achieved by considering one-step finite difference schemes. These finite difference schemes can be solved explicitly or implicitly. Both families generally derive from the Newmark algorithm [1]. The stability of the Newmark algorithm can be demonstrated for linear models, while for non-linear models Belytschko and Schoeberle [2] proved that the discrete energy, computed from the work of internal forces and kinetic energy, is bounded

if it remains positive. Nevertheless, in the non-linear range, the Newmark algorithm leads to a work of internal forces that is different from the internal energy variation. In this context, Hughes et al. [3] have proved that the Newmark algorithm remains physically consistent only for small time step sizes.

In the past decade, interest in simulating non-linear dynamics has kept growing, which led to the development of new time integration algorithms that remain stable in the non-linear range. Since finite difference schemes are well suited to integrate finite-element space-discretizations – they can easily be implemented and are robust – and since the stability requirement can be achieved only with implicit algorithms – in the non-linear range even conditional stability can only be demonstrated for implicit schemes – this paper focuses on implicit finite difference time integration algorithms. For completeness, let us note that the properties of conservation can also be reached by using a Petrov–Galerkin time finite-element [4,5], by using an
approximation of the time Galerkin method [6,7], by using
a Runge–Kutta method (e.g. [8]) or even by considering an
approximation of the time discontinuous Galerkin method
[7,9]. Historically, the first finite difference algorithm verify-
ing the conserving properties in the non-linear range is the
“Energy Momentum Conserving Algorithm” or EMCA
proposed by Simo and Tarnow [10]. It consists of a mid-
point scheme with an adequate evaluation of internal
forces. This adequate evaluation was given for a Saint
Venant–Kirchhoff hyperelastic material. A general formu-
lation in term of the second Piola–Kirchhoff stress tensor
was proposed by Gonzalez [11], and is valid for general
hyperelastic materials. In this last formulation, the stress
tensor is decomposed into two parts: the first one is the
derivative of the free energy evaluated for a mid-configura-
tion strain tensor, and the second is a second-order correc-
tion that enforces the exact variation of energy. This
generalization was allowed by the form of the stress tensor
that derives from an internal energy. Actually, verification
of the conservation of the energy requires such a varia-
tional framework. For example, when considering
energy–momentum conserving algorithms applied to
contact interactions, the creation of a contact energy was
proposed by Armero and Petőcz [12,13] to ensure the
stability.

An issue when considering plasticity is the absence of an
energy potential from which to derive the stress tensor.
Energy-momentum conserving algorithms were then devel-
oped by Meng and Laursen [14,15] in such a way that the
work of internal forces is dissipative, consistently with plas-
tic behavior. They were also developed by Noels et al.
[16,17] and Armero [18] in such a way that for respectively
hypoelastic and hyperelastic based elasto-plastic materials,
the work of internal forces corresponds to the change of
reversible energy plus the plastically dissipated energy. In
general, the absence of an energy potential leads to a less
elegant mathematical formulation of the problem and can
lead to more complicated formulations of the internal
forces. This is particularly true when considering energy-
momentum conserving algorithms.

In a recent work, the authors [19] proposed to use the
variational updates framework to develop conserving algo-
rithms for elasto-plasticity. This variational constitutive
updates formulation, initially proposed by Radovitzky
and Ortiz [20] and developed by Ortiz and Stainier [21],
has the interesting property that the stress tensor always
derives from an incremental potential, even for elasto-plas-
tic models or visco-elastic models, see Fancello et al. [22].
This variational framework was further extended to the
general thermo-mechanical case by Yang et al. [23]. There-
fore, in such a framework we can use the formulation based
on the second Piola–Kirchhoff stress tensor as proposed by
Gonzalez [11] in a straightforward way. Moreover, the use
of the variational formulation does not lead to any *a priori*
restrictions on the material laws or parameters, even if in
this paper we focus on elasto-plasticity with isotropic
hardening.

Although recent developments in energy–momentum
conserving algorithms have demonstrated the accuracy of
the methods for non-linear structural dynamics, even for
elasto-plastic behaviors, the solution obtained can be pol-
luted by high-frequency numerical oscillations. These
non-physical modes are also present when time-integrating
with classical Newmark algorithms, and, in this particular
case, are commonly dissipated by introducing numerical
damping, leading to the generalized-α methods [24]. Never-
theless, the unconditional stability of these methods is
guaranteed only for linear systems or asymptotically for
high frequencies in the non-linear range [25]. It then
appears natural to introduce numerical dissipation in the
energy–momentum conserving algorithms, which were
therefore renamed Energy-Dissipative Momentum-Con-
serving (EDMC) algorithms by Armero and Romero
[26,27] who presented the first formulation of this kind
for hyperelastic materials. In the same way, the authors
[28] introduced dissipation for elasto-plastic hypoelastic-
based materials, which allowed the simulations of complex
phenomena such as blade loss in a turbo-engine. The nat-
ural path is then to introduce numerical dissipation into
the conserving formulation based on the variational consti-
tutive updates. Contrarily to previous work in the area, this
dissipation is introduced by adding a dissipative energy to
the potential resulting from the variational updates frame-
work. Therefore, the modification is not only introduced in
the second-order correction of Gonzalez’ stress formula-
tion, but also in the main part computed for the mid-con-
figuration strain tensor, ensuring that the correction of the
stress tensor remains second-order compared to the main
usual contribution, i.e. the stress tensor is equal to the der-
ivation of the energy (main part) corrected by a tensor
which is second order in the deformation increment.

This paper is organized as follow. In Section 2 the equa-
tions governing the continuum problem are described.
Attention is paid to the definition of an effective potential
when plasticity occurs, and how numerical dissipation
can be introduced in the resulting variational formulation.
These equations are then adapted to a finite-element dis-
cretization and to a finite time step increment in Section
3. It is then shown in Section 4 that the proposed varia-
tional formulation naturally leads to an Energy-Dissipative
Momentum-Conserving time integration, for elasto-plastic
behavior. This scheme is obtained by applying directly the
formulation of Gonzalez. The accuracy, robustness and
first-order property of the scheme are therefore demon-
strated in Section 5 by considering numerical examples.

## 2. Continuous dynamics

In this section we establish the variational formalism of
continuous dynamics, based on constitutive updates for
elasto-plastic behaviors. Moreover, attention is paid to
the introduction of numerical dissipation in this formula-
tion, conservation laws deduced from the resulting weak
formulation are also recalled.
## 2.1. Notations

Let $\mathcal{V}_0 \subset \mathbb{R}^3$ be the reference configuration (in the Euclidean space) of the body under study at time $t = t_0$. Given this convention, the material and spatial configurations of this representation are initially identical. From now on the subscript 0 will refer to the reference configuration. The deformation of this body in the time interval $\mathcal{T} = [t_0, t_f]$ is defined by the mapping $\vec{\phi} : \mathcal{V}_0 \times \mathcal{T} \to \mathcal{V}$, with

$$ \vec{x} \equiv \vec{\phi}(\vec{x}_0, t), \quad (1) $$

where $\vec{x}$ and $\vec{x}_0$ are respectively the current and reference positions. The two-point deformation gradient $\mathbf{F} : \mathcal{V}_0 \times \mathcal{T} \to \text{GL}_+(3, \mathbb{R})$, associated to the mapping, is restricted to $\text{GL}_+(3, \mathbb{R})$ the Lie group of invertible, orientation preserving linear transformations in $\mathbb{R}^3$ characterized by a positive determinant, and is defined by

$$ \mathbf{F} \equiv \frac{\partial \vec{x}}{\partial \vec{x}_0} = \vec{\nabla}_0 \vec{\phi}, \quad \text{with } J \equiv \det(\mathbf{F}) : \mathcal{V}_0 \times \mathcal{T} \to \mathbb{R}^+ \quad (2) $$

its Jacobian. The Jacobian allows to evaluate the time evolution of the density $\rho : \mathcal{V}_0 \times \mathcal{T} \to \mathbb{R}^+$ by the relation

$$ \rho J = \rho_0. \quad (3) $$

The reference boundary $\partial \mathcal{V}_0$ of the body $\mathcal{V}_0$ is oriented by the unit outward normal $\vec{n}_0$ and is decomposed into a Dirichlet part $\partial_D \mathcal{V}_0$ and a Neumann part $\partial_N \mathcal{V}_0$ where respectively the displacement and surface traction are constrained to $\vec{x}$ and $\vec{t}$, with

$$ \partial_N \mathcal{V}_0 \cap \partial_D \mathcal{V}_0 = \emptyset \quad \text{and} \quad \partial_N \mathcal{V}_0 \cup \partial_D \mathcal{V}_0 = \partial \mathcal{V}_0. \quad (4) $$

Given these definitions, the motion of the body $\vec{\phi}(t)$ is defined by the time evolution of the position $\vec{x} : \mathcal{V}_0 \times \mathcal{T} \to \mathbb{R}^3$ in the admissible manifold

$$ \mathcal{X} \equiv \{\vec{\phi}(t) : \mathcal{V}_0 \times \mathcal{T} \to \mathcal{V} \mid [J > 0, \vec{x}|_{\partial_D \mathcal{V}} = \vec{x}]\}. \quad (5) $$

## 2.2. Weak formulation

Let $\boldsymbol{\sigma}$ be the Cauchy stress tensor and $\vec{b}, \vec{t}$ respectively the volume forces and surface tractions applied to the body. Therefore, the continuum linear momentum equilibrium is stated (in the current configuration) by

$$ \vec{\nabla} \cdot \boldsymbol{\sigma} + \rho \vec{b} = \rho \ddot{\vec{\phi}} \quad \forall \vec{x}_0 \in \mathcal{V}_0, \quad (6) $$

$$ \boldsymbol{\sigma} \cdot \vec{n} = \vec{t} \quad \forall \vec{x}_0 \in \partial_N \mathcal{V}_0. \quad (7) $$

Let us consider the manifold of admissible virtual displacements

$$ \begin{aligned} \mathcal{D} \equiv \{ & \delta\vec{\phi}(t) : \mathcal{V}_0 \times \mathcal{T} \to \mathbb{R}^3 | [\delta\vec{\phi}(t)|_{\partial_D \mathcal{V}} = \vec{0} \quad \forall t \in \mathcal{T} \\ & \text{and } \delta\vec{\phi}(\vec{x}_0, t_0) = \delta\vec{\phi}(\vec{x}_0, t_f) = \vec{0} \quad \forall \vec{x}_0 \in \mathcal{V}_0 \} \end{aligned} \quad (8) $$

let us multiply Eq. (6) by such an admissible displacement and integrate on the body $\mathcal{V}$ and the time interval $\mathcal{T}$. When no confusion is possible, in order to simplify the notations, $\vec{\phi}(t)$ and $\delta\vec{\phi}(t)$ will be replaced in the equations

by, respectively, $\vec{\phi}$ and $\delta\vec{\phi}$. The weak form of the problem is then stated as finding $\vec{\phi}(t) \in \mathcal{X}$ such that

$$ \begin{aligned} & \int_{t_0}^{t_f} \left\{ \int_{\mathcal{V}_0} [\rho \ddot{\vec{\phi}} - \vec{\nabla} \cdot \boldsymbol{\sigma} - \rho \vec{b}] \cdot \delta\vec{\phi} \, d\mathcal{V} \right\} dt = 0 \quad \forall \delta\vec{\phi} \in \mathcal{D}, \\ & \forall t \in \mathcal{T}. \end{aligned} \quad (9) $$

Using Eqs. (2), (3), (5), (7) and (8), integration by parts of the stress divergence leads to a weak formulation of the problem stated in the reference configuration, which is finding $\vec{\phi}(t) \in \mathcal{X}$ such that

$$ \begin{aligned} & \int_{t_0}^{t_f} \left\{ \int_{\mathcal{V}_0} \rho_0 \ddot{\vec{\phi}} \cdot \delta\vec{\phi} \, d\mathcal{V}_0 + \int_{\mathcal{V}_0} [\boldsymbol{\sigma}\mathbf{F}^{-T}J] : \vec{\nabla}_0 \delta\vec{\phi} \, d\mathcal{V}_0 \right. \\ & \qquad \left. - \int_{\partial_N \mathcal{V}_0} \delta\vec{\phi} \cdot \vec{T} \, d\partial\mathcal{V}_0 - \int_{\mathcal{V}_0} \rho_0 \vec{b} \cdot \delta\vec{\phi} \, d\mathcal{V}_0 \right\} dt = 0 \\ & \forall \delta\vec{\phi} \in \mathcal{D}, \quad \forall t \in \mathcal{T}, \end{aligned} \quad (10) $$

where $\vec{T}$ is the surface traction $\vec{t}$ evaluated in the reference configuration.<fnref n="3" />

Since we are now referring to the reference configuration, it is natural to use the symmetric second Piola–Kirchhoff stress tensor

$$ \mathbf{S} \equiv J \mathbf{F}^{-1} \boldsymbol{\sigma} \mathbf{F}^{-T}. \quad (11) $$

Therefore, owing to the arbitrary nature of $\delta\vec{\phi}$, the weak form (10) reduces to finding $\vec{\phi}(t) \in \mathcal{X}$ such that

$$ \begin{aligned} & \int_{\mathcal{V}_0} \rho_0 \ddot{\vec{\phi}} \cdot \delta\vec{\phi} \, d\mathcal{V}_0 + \int_{\mathcal{V}_0} [\mathbf{F}\mathbf{S}] : \vec{\nabla}_0 \delta\vec{\phi} \, d\mathcal{V}_0 \\ & = \int_{\partial_N \mathcal{V}_0} \delta\vec{\phi} \cdot \vec{T} \, d\partial\mathcal{V}_0 + \int_{\mathcal{V}_0} \rho_0 \vec{b} \cdot \delta\vec{\phi} \, d\mathcal{V}_0 \quad \forall \delta\vec{\phi} \in \mathcal{D}, \\ & \forall t \in \mathcal{T}. \end{aligned} \quad (12) $$

This weak form will be useful to define the finite-elements formulation. But before proceeding the constitutive behavior linking the deformation gradient to the stress tensor needs to be defined in order to complete the formalism.

## 2.3. Variational formulation of the constitutive behavior

In order to relate the stress tensor to a functional, we use the variational framework of constitutive updates detailed in [20,21] particularized to elasto-plastic constitutive models. Toward this end, the strain tensor (2) is multiplicatively decomposed into a plastic part $\mathbf{F}^{\text{pl}} \in \text{SL}(3, \mathbb{R})$ and into an elastic part $\mathbf{F}^{\text{el}} \in \text{GL}_+(3, \mathbb{R})$ as

$$ \mathbf{F} = \mathbf{F}^{\text{el}} \mathbf{F}^{\text{pl}}, \quad (13) $$

with $\text{SL}(3, \mathbb{R}) \equiv \{\mathbf{F}^{\text{pl}} \in \text{GL}(3, \mathbb{R}) | [\det(\mathbf{F}^{\text{pl}}) = 1]\}$. The plastic deformation gradient can be coupled to internal variables by defining a flow rule. In the particular case of a von Mises flow rule, one has

$$ \dot{\mathbf{F}}^{\text{pl}} = \varepsilon^{\text{pl}} \mathbf{N} \mathbf{F}^{\text{pl}}, \quad (14) $$

<footnote n="3">³ One has $\vec{t} d\partial\mathcal{V} = \vec{T} d\partial\mathcal{V}_0$, while Eq. (7) becomes $[\boldsymbol{\sigma}\mathbf{F}^{-T}J] \cdot \vec{n}_0 = \vec{T}$, using Nanson's formula $\vec{n} d\partial\mathcal{V} = J\mathbf{F}^{-T}\vec{n}_0 d\partial\mathcal{V}_0$.</footnote>
where $\dot{\epsilon}^P \in \mathbb{R}_+$ is the equivalent plastic strain rate, and where $\mathbf{N} \in \mathcal{N}$ is the flow direction, with $\mathcal{N} = \{\mathbf{N} \in \text{Sym}(3, \mathbb{R}) | [\text{tr}(\mathbf{N}) = 0 \text{ and } \mathbf{N}: \mathbf{N} = 3/2]\}$.

It is also assumed that the Helmholtz free energy function $A: GL_+(3, \mathbb{R}) \times \mathbb{R}_+ \times \mathcal{N} \to \mathbb{R}_+$ can be additively decomposed into

$$A(\mathbf{F}, \epsilon^P, \mathbf{N}) \equiv \Phi^{\text{el}}(\mathbf{F}\mathbf{F}^{\text{pl}})^{-1}(\epsilon^P, \mathbf{N}) + \Phi^{\text{pl}}(\mathbf{F}^{\text{pl}}(\epsilon^P, \mathbf{N}), \epsilon^P), \quad (15)$$

where $\Phi^{\text{el}}: GL_+(3, \mathbb{R}) \to \mathbb{R}$ is the elastic potential and $\Phi^{\text{pl}}: SL(3, \mathbb{R}) \times \mathbb{R}_+ \to \mathbb{R}_+$ is the plastic potential. This additive decomposition translates the observation that, in most metals, the plastic work hardening does not influence the elastic (reversible) properties. In a general way, the stress tensor $\mathbf{S}: GL_+(3, \mathbb{R}) \times \mathbb{R}_+ \times \mathcal{N} \to \text{Sym}^2(3, \mathbb{R})$, the thermodynamic forces $\mathbf{T}: GL_+(3, \mathbb{R}) \times \mathbb{R}_+ \times \mathcal{N} \to GL(3, \mathbb{R})$ and $Y: GL_+(3, \mathbb{R}) \times \mathbb{R}_+ \times \mathcal{N} \to \mathbb{R}$ explicitly derive from this energy, with respectively

$$\mathbf{S} = 2 \frac{\partial A}{\partial \mathbf{C}}, \quad \mathbf{T} = - \frac{\partial A}{\partial \mathbf{F}^{\text{pl}}} \quad \text{and} \quad Y = - \frac{\partial A}{\partial \epsilon^P}, \qquad (16)$$

where the explicit dependency of $A$ on the right Cauchy strain tensor $\mathbf{C} = \mathbf{F}^T \mathbf{F}$ leads to the frame indifference. A dissipative pseudo-potential $\Psi: \mathbb{R}_+ \to \mathbb{R}_+$ can be associated to $\epsilon^P$ such that

$$\dot{\epsilon}^P = \frac{\partial \Psi(Y)}{\partial Y}, \qquad (17)$$

while its dual part $\Psi^*: \mathbb{R}_+ \to \mathbb{R}_+$ is obtained by the Legendre–Fenchel transform

$$\Psi^*(\dot{\epsilon}^P) = \sup_Y (Y \dot{\epsilon}^P - \Psi(Y)) \quad \text{with} \quad Y = \frac{\partial \Psi^*(\dot{\epsilon}^P)}{\partial \dot{\epsilon}^P}. \qquad (18)$$

Convexity of $\Psi$ (or $\Psi^*$) ensures the positiveness of the dissipation. The physical interpretation of $\mathbf{T}$ and $Y$ is given as follows. Combining Eqs. (13)–(16) leads to

$$\mathbf{S} = 2 \frac{\partial A}{\partial \mathbf{C}} = 2 \mathbf{F}^{\text{pl}} \frac{\partial \Phi^{\text{el}}(\mathbf{C}^{\text{el}})}{\partial \mathbf{C}^{\text{el}}} \mathbf{F}^{\text{pl}} \qquad (19)$$

$$\mathbf{T} = \mathbf{F}^{\text{el}^T} \mathbf{F} \mathbf{S} - \frac{\partial \Phi^{\text{pl}}(\mathbf{F}^{\text{pl}}, \epsilon^P)}{\partial \mathbf{F}^{\text{pl}}} = \mathbf{F}^{\text{el}^T} \mathbf{F} \mathbf{S} - \mathbf{T}_c \quad \text{and} \qquad (20)$$

$$Y = - \frac{\partial A}{\partial \mathbf{F}^{\text{pl}}} : \frac{\partial \mathbf{F}^{\text{pl}}}{\partial \epsilon^P} - \frac{\partial \Phi^{\text{pl}}(\epsilon^P)}{\partial \epsilon^P} = \mathbf{T} : [\mathbf{N} \mathbf{F}^{\text{pl}}] - \Phi^{\text{pl}}_{, \epsilon^P}, \qquad (21)$$

where $\mathbf{T}_c$ is a back-stress tensor and where $\Phi^{\text{pl}}_{, \epsilon^P}$ depends on the plastic potential.

From these general definitions, Ortiz and Stainier [21] proposed the new three-field functional $\dot{D}: GL_+(3, \mathbb{R}) \times \mathbb{R}_+ \times \mathcal{N} \to \mathbb{R}$ defined by

$$\dot{D}(\dot{\mathbf{F}}, \dot{\epsilon}^P, \mathbf{N}) \equiv \frac{\partial A}{\partial \dot{\mathbf{F}}} : \dot{\mathbf{F}} - Y \dot{\epsilon}^P + \Psi^*(\dot{\epsilon}^P) = \dot{A} + \Psi^*(\dot{\epsilon}^P). \qquad (22)$$

The minimization of the power $\dot{D}$ with respect to $\dot{\epsilon}^P$ gives back (18), while minimizing with respect to $\mathbf{N}$ corresponds to constrain $\mathbf{N}$ in the direction of the deviatoric stress $\mathbf{T}\mathbf{F}^{\text{pl}^T}$ (this corresponds to the principle of maximum plastic dissipation [29]). Moreover, the variation with respect to $\dot{\mathbf{F}}$ leads to

$$\frac{\partial \dot{D}(\dot{\mathbf{F}}, \dot{\epsilon}^P, \mathbf{N})}{\partial \dot{\mathbf{F}}} = \frac{\partial A}{\partial \dot{\mathbf{F}}} = \mathbf{F} \mathbf{S}. \qquad (23)$$

Since these variations of $\dot{D}$ do respect the prescribed conditions of the constitutive material, a one field functional can be defined by

$$\dot{D}^{\text{eff}}(\vec{x}_0, \dot{\mathbf{F}}) = \min_{\epsilon^P, \mathbf{N}} \dot{D}(\dot{\mathbf{F}}, \epsilon^P, \mathbf{N}), \quad \text{with} \quad \mathbf{S} = \mathbf{F}^{-1} \frac{\partial \dot{D}^{\text{eff}}(\dot{\mathbf{F}})}{\partial \dot{\mathbf{F}}}. \qquad (24)$$

The effective rate potential $\dot{D}^{\text{eff}}$ – it has the dimension of a power – corresponds to the stationary point of the functional $\dot{D}$ with respect to the internal variables.

The main idea of the present work is to add to this power a term that will introduce numerical dissipation during the time integration. In a general way, we define the numerical dissipative power $\Psi^d: GL_+(3, \mathbb{R}) \to \mathbb{R}_+$ as

$$\Psi^d \equiv \Psi^d(\dot{\mathbf{C}}^{\text{el}}), \qquad (25)$$

which depends on the rate of the elastic deformations. The plastic deformation used to compute this power is obtained by the minimization of $\dot{D}$. Indeed the minimization process yields the solution to the constitutive elasto-plastic problem, and the numerical dissipation should not modify this solution. This is the reason why the numerical dissipation is added to the functional after the minimization process instead of being added to the function to be minimized. Since the numerical dissipative power depends on the symmetric tensor $\mathbf{C}^{\text{el}}$, the second Piola–Kirchhoff stress tensor remains symmetric. So Eq. (24) is rewritten

$$\dot{D}^{\text{eff}}(\vec{x}_0, \dot{\mathbf{F}}) = \min_{\epsilon^P, \mathbf{N}} [\dot{D}(\dot{\mathbf{F}}, \epsilon^P, \mathbf{N})] + \Psi^d(\dot{\mathbf{C}}^{\text{el}}), \quad \text{with} \qquad (26)$$

$$\mathbf{S} = \mathbf{F}^{-1} \frac{\partial \dot{D}^{\text{eff}}(\dot{\mathbf{F}})}{\partial \dot{\mathbf{F}}}. \qquad (27)$$

These expressions define a stress tensor from a functional form. This resulting stress tensor can be seen as a combination of two terms: the first term is the derivation of the functional minimum and verifies the consistency with the original underlying elasto-dynamics problem (e.g. normality of the plastic flow, isotropy for isotropic laws, ...), while the second term is an artificial stress which purpose is to introduce numerical dissipation during the time increment scheme. The elastic tensor needed to compute this second term is obtained from the plastic flow determined during the minimization process. Consistency of the time stepping scheme will be ensured providing the form of this dissipation potential is chosen such that $\lim_{\Delta t \to 0^+} \frac{\partial \Psi^d}{\partial \dot{\mathbf{F}}} = 0$, i.e. the dissipation term should vanish for decreasing time steps. A specific expression of $\Psi^d$, verifying this condition, will be given below.

## 2.4. Laws of conservation

As it is shown in the forthcoming lines, assuming pure Neumann boundary conditions, the weak formulation
(12) verifies the conservation of the linear momentum defined by

$$ \vec{L} = \int_{\gamma_0} \rho_0 \dot{\vec{\phi}} d\gamma_0 : \mathcal{T} \rightarrow \mathbb{R}^3, \qquad (28) $$

of the angular momentum defined by

$$ \vec{J} = \int_{\gamma_0} \rho_0 \vec{\phi} \wedge \dot{\vec{\phi}} d\gamma_0 : \mathcal{T} \rightarrow \mathbb{R}^3, \qquad (29) $$

and of the total energy defined by

$$ E = K + U_{\text{el}} : \mathcal{T} \rightarrow \mathbb{R}_{+}, \qquad (30) $$

with respectively

$$ K = \int_{\gamma_0} \frac{1}{2} \rho_0 \dot{\vec{\phi}}^2 d\gamma_0, \quad \text{and} \qquad (31) $$

$$ U_{\text{el}} = \int_{\gamma_0} \Phi^{\text{el}} d\gamma_0, \qquad (32) $$

the kinetic and the reversible stored energy. Indeed, taking $\delta\vec{\phi} = \vec{c} \in \mathcal{D}$, where $\vec{c}$ is an arbitrary constant, in Eq. (12) leads to

$$ \dot{\vec{L}} = \int_{\gamma_0} \rho_0 \vec{b} d\gamma_0 + \int_{\partial_N \gamma_0} \vec{T} d\partial\gamma_0 \quad \forall t \in \mathcal{T}, \qquad (33) $$

which is the linear momentum conservation law. Similarly, the choice $\delta\vec{\phi} = \vec{c} \wedge \vec{\phi} \in \mathcal{D}$, leads to

$$ \dot{\vec{J}} = \int_{\gamma_0} \rho_0 \vec{\phi} \wedge \vec{b} d\gamma_0 + \int_{\partial_N \gamma_0} \vec{\phi} \wedge \vec{T} d\partial\gamma_0 \quad \forall t \in \mathcal{T}, \qquad (34) $$

since S is symmetric, which is the angular momentum conservation law. Eventually, the choice $\delta\vec{\phi} = \dot{\vec{\phi}} \in \mathcal{D}$ leads to

$$ \begin{aligned} \dot{K} + \int_{\gamma_0} D^{\text{eff}} d\gamma_0 &= \int_{\gamma_0} \rho_0 \dot{\vec{\phi}} \cdot \vec{b} d\gamma_0 + \int_{\partial_N \gamma_0} \dot{\vec{\phi}} \cdot \vec{T} d\partial\gamma_0 \\ &= \dot{W}_{\text{ext}} \quad \forall t \in \mathcal{T}, \end{aligned} \qquad (35) $$

where $\dot{W}_{\text{ext}}$ is the power of work the external forces applied to $\mathcal{V}$. Let us note that $\dot{D}^{\text{eff}}$ is equal to the power of internal forces (per unit undeformed volume) under the condition that $\Psi^*$ is homogeneous of order one in its arguments, which is only the case for rate-independent behavior.<fnref n="4" /> This power can be decomposed into

$$ \dot{W}_{\text{int}} = \int_{\gamma_0} \dot{D}^{\text{eff}} d\gamma_0 = \dot{U}_{\text{el}} + \dot{U}_i + \dot{W}_{\text{diss}} \quad \forall t \in \mathcal{T}, \qquad (36) $$

where

$$ \begin{aligned} \dot{W}_{\text{diss}} &= \int_{\gamma_0} \Psi^d d\gamma_0 \quad \text{and} \\ \dot{U}_i &= \int_{\gamma_0} \dot{D}^{\text{eff}} - \dot{\Phi}^{\text{el}} - \Psi^d d\gamma_0, \end{aligned} \qquad (37) $$

are respectively the numerically dissipated power and the irreversible power resulting from plastic behavior. Combining Eqs. (30), (35) and (37) leads to the conservation of energy

$$ \dot{K} + \dot{U}_{\text{el}} + \dot{U}_i + \dot{W}_{\text{diss}} = \dot{W}_{\text{ext}} \quad \forall t \in \mathcal{T}. \qquad (38) $$

### 3. Finite-elements discretization and time integration

In order to solve the weak formulation (12), a finite-element discretization is introduced. In order to perform the time integration, the constitutive updates formulation is adapted to the context of time increments, leading to an incremental functional. In addition we consider a formulation adapted to a quasi-incompressible constitutive behavior. It is then demonstrated that, in combination with the Energy-Dissipative Momentum-Conserving (EDMC) scheme, Gonzalez' expression [11] of the modified second Piola-Kirchhoff stress tensor $S^*$ can be used directly from this incremental functional in order to verify the conservation laws. The study of the spectral matrix associated to this time-integration method is also studied.

#### 3.1. Finite element discretization

The body $\mathcal{V}$ is now discretized into finite-elements $\mathcal{V}^e$, with $\mathcal{V} = \cup_e \mathcal{V}^e$, $\cap_e \mathcal{V}^e = \emptyset$. The approximation $\vec{\phi}_h(t)$ of the solution $\vec{\phi}(t)$ is a polynomial approximations $\mathbb{P}^k$ on the sub-domains $\mathcal{V}^e$. The solution is therefore restrained in the manifold

$$ \mathcal{X}_h = \{ \vec{\phi}_h(t) \in \mathcal{X} | [\vec{\phi}_h(t) \in C^0(\mathcal{V}) \text{ and } \vec{\phi}_h \in \mathbb{P}^k(\mathcal{V}^e)] \}. \qquad (39) $$

Commonly, nodes $\xi$ are associated to the finite elements, and shape functions $N^\xi: \gamma_0 \rightarrow \mathbb{R}$ are defined with the usual relations

$$ \begin{aligned} \vec{\phi}_h(\vec{x}_0, t) &= \sum_{\xi} N^{\xi}(\vec{x}_0) \vec{x}_0^{\xi}(t) \quad \text{and} \\ \delta\vec{\phi}_h(\vec{x}_0, t) &= \sum_{\xi} N^{\xi}(\vec{x}_0) \delta\vec{x}_0^{\xi}(t), \end{aligned} \qquad (40) $$

with $\vec{x}_0^\xi$ the positions of the approximate solution at node $\xi$. Introducing Eqs. (40) in Eq. (12), leads to the new finite elements weak form, which is finding $\vec{x}_0^\xi$ such that

$$ M^{\xi\mu} \ddot{\vec{x}}^{\mu} + \vec{f}_{\text{int}}^{\xi} = \vec{f}_{\text{ext}}^{\xi}, \quad \forall t \in \mathcal{T}, \quad \forall \xi, \qquad (41) $$

with

$$ M^{\xi\mu} = \int_{\gamma_0} \rho_0 N^{\xi} N^{\mu} d\gamma_0, \qquad (42) $$

$$ \vec{f}_{\text{int}}^{\xi} = \sum_{e} \int_{\gamma_0^e} \mathbf{F} \vec{D}^{\xi} d\gamma_0 \quad \text{and} \qquad (43) $$

$$ \vec{f}_{\text{ext}}^{\xi} = \int_{\partial_N \gamma_0} \vec{T} N^{\xi} d\partial\gamma_0 + \int_{\gamma_0} \rho_0 N^{\xi} \vec{b} d\gamma_0, \qquad (44) $$

respectively the mass matrix, the internal forces and the external forces, where $\vec{D}^\xi = \vec{\nabla}_0 N^\xi$ is the derivative of the shape functions.

While considering finite-elements discretization, locking in pressure is common, especially for linear elements. This issue can be solved by the method proposed by Simo and

<footnote n="4">⁴ Indeed, $\partial_{\vec{\phi}} \dot{\Psi}^* \dot{\epsilon}^p = \Psi^*$ only if $\Psi^*$ is linear with $\dot{\epsilon}^p$.</footnote>
Taylor [30] who considered a constant volume deformation on each element. So let us introduce the following definitions. Let $\theta^e \in \mathbb{R}_+$ be a constant value on $\mathcal{V}^e$, which allows to define the distortion gradient $\hat{\mathbf{F}} : \mathcal{V}_0 \times \mathcal{T} \to \text{SL}(3, \mathbb{R})$ and the modified deformation gradient $\bar{\mathbf{F}} : \mathcal{V}_0 \times \mathcal{T} \times \mathbb{R}_+ \to \text{SL}_+(3, \mathbb{R})$ respectively as

$$ \hat{\mathbf{F}} \equiv J^{-\frac{1}{3}} \mathbf{F} \quad \text{and} \quad \bar{\mathbf{F}} \equiv \theta^{\frac{1}{3}} \hat{\mathbf{F}} = \left[ \frac{\theta^e}{J} \right]^{\frac{1}{3}} \mathbf{F}, \quad \text{with} \quad (45) $$

$$ \hat{\mathbf{C}} \equiv \hat{\mathbf{F}}^T \hat{\mathbf{F}} = \left[ \frac{1}{J} \right]^{\frac{2}{3}} \mathbf{C} \quad \text{and} \quad \bar{\mathbf{C}} \equiv \bar{\mathbf{F}}^T \bar{\mathbf{F}} = \left[ \frac{\theta^e}{J} \right]^{\frac{2}{3}} \mathbf{C}. \quad (46) $$

## 3.2. Incremental updates of the variational elasto-plastic model

The time-integration of Eqs. (41) is accomplished via an incremental solution procedure in which the time interval of interest $\mathcal{T}$ is discretized into $n'$ time steps such that $\mathcal{T} = \bigcup_{n=0}^{n=n'+1} [t^n, t^{n+1}]$ and $\Delta t = t^{n+1} - t^n$ is the time step size. Therefore, the formulation described in Section 2.3 is modified in order to take into account the finite nature of the deformation gradients. The purpose is to adapt the variational power $D^{\text{eff}}$ defined in Eq. (24), leading to an incremental form $\Delta D^{\text{eff}}(\vec{x}_0, \bar{\mathbf{C}}(\vec{\phi}_h, \theta^e))$.

The incremental form of Eq. (14) is

$$ \mathbf{F}^{\text{pl}n+1} = \exp([\varepsilon^{\text{pl}n+1} - \varepsilon^{\text{pl}}]\mathbf{N}) \mathbf{F}^{\text{pl}n} = \mathbf{A}(\Delta \varepsilon^{\text{pl}}) \mathbf{F}^{\text{pl}n}, \quad (47) $$

where tensor $\mathbf{A}(\Delta \varepsilon^{\text{pl}})$ verifies

$$ \det(\mathbf{A}(\Delta \varepsilon^{\text{pl}})) = \exp(\text{tr}(\Delta \varepsilon^{\text{pl}}\mathbf{N})) = 1 \quad \text{and} \\ \Delta \varepsilon^{\text{pl}} = \sqrt{\frac{2}{3}} \| \ln \mathbf{A}(\Delta \varepsilon^{\text{pl}}) \|, \quad (48) $$

with $\Delta \varepsilon^{\text{pl}} = \varepsilon^{\text{pl}n+1} - \varepsilon^{\text{pl}}$ the increment of equivalent plastic strain.

In order to use the quasi-incompressible technique defined in Section 3.3, we assume a split of the elastic potential $\Phi^{\text{el}}$ into a volume part $\Phi^{\text{vol}}_{\text{vol}} : \mathbb{R}_+ \to \mathbb{R}_+$ and into a deviatoric part $\Phi^{\text{dev}}_{\text{dev}} : \text{SL}(3, \mathbb{R}) \to \mathbb{R}_+$. Therefore the free Helmholtz energy (15) is rewritten

$$ A(\bar{\mathbf{F}}(\vec{\phi}_h, \theta^e), \varepsilon^{\text{pl}}, \mathbf{N}) = \Phi^{\text{vol}}_{\text{vol}}(\theta^e) + \Phi^{\text{vol}}_{\text{dev}}(\hat{\mathbf{C}}^{\text{vol}}) \\ + \Phi^{\text{pl}}(\mathbf{F}^{\text{pl}}(\varepsilon^{\text{pl}}, \mathbf{N}), \varepsilon^{\text{pl}}), \quad (49) $$

where it has been taken into account that the plastic deformation gradient $\mathbf{F}^{\text{pl}}$ is deviatoric, i.e. $\text{det}(\mathbf{F}^{\text{pl}}) = 1$. Therefore, $\Delta D : \text{GL}_+(3, \mathbb{R}) \times \mathbb{R}_+ \times \mathcal{N}$ the incremental form of Eq. (22) is stated as

$$ \Delta D(\bar{\mathbf{F}}(\vec{\phi}_h^{n+1}, \theta^{en+1}), \varepsilon^{\text{pl}n+1}, \mathbf{N}) \\ = \Phi^{\text{vol}}_{\text{vol}}(\theta^{en+1}) - \Phi^{\text{vol}}_{\text{vol}}(\theta^{en}) + \Delta D_{\text{dev}}(\hat{\mathbf{F}}, \varepsilon^{\text{pl}}, \mathbf{N}), \quad \text{with} \quad (50) $$

$$ \begin{aligned} \Delta D_{\text{dev}}(\hat{\mathbf{F}}, \varepsilon^{\text{pl}}, \mathbf{N}) ={}& \Phi^{\text{vol}}_{\text{dev}}(\hat{\mathbf{C}}^{\text{vol}n+1}(\hat{\mathbf{F}}^{n+1}, \varepsilon^{\text{vol}n+1}, \mathbf{N})) \\ & - \Phi^{\text{vol}}_{\text{dev}}(\hat{\mathbf{C}}^{\text{vol}n}) \\ & + \Phi^{\text{pl}}(\mathbf{F}^{\text{pl}n+1}(\varepsilon^{\text{vol}n+1}, \mathbf{N}), \varepsilon^{\text{vol}n+1}) \\ & - \Phi^{\text{pl}}(\mathbf{F}^{\text{pl}n}(\varepsilon^{\text{vol}n}), \varepsilon^{\text{vol}n}) + \Delta t \Psi^*(\frac{\Delta \varepsilon^{\text{pl}}}{\Delta t}). \end{aligned} \quad (51) $$

Similarly to what has been done in Section 2.3, since minimization of $\Delta D$ with respect to $\varepsilon^{\text{vol}n+1}$ leads to verify (18), and since minimization with respect to $\mathbf{N}$ leads to a radial return mapping scheme [19,21], the effective incremental potential $\Delta D^{\text{eff}} : \mathcal{V}_0 \times \text{GL}_+(3, \mathbb{R}) \to \mathbb{R}_+$ can be defined as the discrete form of Eq. (26), leading to

$$ \begin{aligned} \Delta D^{\text{eff}}(\vec{x}_0, \bar{\mathbf{F}}(\vec{\phi}_h^{n+1}, \theta^{en+1})) ={}& \min_{\varepsilon^{\text{vol}n+1}\mathbf{N}} [\Delta D] + \Delta t \Psi^{\text{vol}}_{\text{vol}}(\Delta \theta^e, \Delta t) \\ & + \Delta t \Psi^{\text{dev}}_{\text{dev}}(\Delta \hat{\mathbf{C}}^{\text{vol}}, \Delta t). \end{aligned} \quad (52) $$

Consistently with the continuous case (see Section 2.3), the numerical dissipation $\Psi^d$ is added after minimization of the energy functional increment $\Delta D$, in order to preserve the elasto-plastic solution of the constitutive problem. Once the plastic flow properties have been determined by the minimization process, the elastic tensor is known and the numerical dissipation can be evaluated directly. Nevertheless, when defining the explicit form of the dissipation, consistency of the method requires that $\lim_{\Delta t \to 0^+} \frac{\partial \Delta t}{\partial \hat{\mathbf{F}}} = 0$. For convenience, the numerical dissipation $\Psi^d$ has also been split into a volume part $\Psi^{\text{vol}}_{\text{vol}} : \mathbb{R}_+ \to \mathbb{R}_+$ and into a deviatoric part $\Psi^{\text{dev}}_{\text{dev}} : \text{SL}(3, \mathbb{R}) \to \mathbb{R}_+$. The dependence of the new terms being only to scalar and to symmetric tensors, the second Piola–Kirchhoff stress tensor remains symmetric. Since the minimization process is independent on $\theta^e$, this effective incremental potential can be split into the volume part $\Delta D^{\text{eff}}_{\text{vol}} : \mathbb{R}_+ \to \mathbb{R}_+$ and into the deviatoric part $\Delta D^{\text{eff}}_{\text{dev}} : \text{SL}(3, \mathbb{R}) \to \mathbb{R}_+$, such that

$$ \Delta D^{\text{eff}}(\vec{x}_0, \bar{\mathbf{F}}(\vec{\phi}_h^{n+1}, \theta^{en+1})) = \Delta D^{\text{eff}}_{\text{vol}}(\theta^{en+1}) + \Delta D^{\text{eff}}_{\text{dev}}(\hat{\mathbf{F}}^{n+1}), \quad (53) $$

with

$$ \begin{aligned} \Delta D^{\text{eff}}_{\text{vol}}(\theta^{en+1}) ={}& \Phi^{\text{vol}}_{\text{vol}}(\theta^{en+1}) - \Phi^{\text{vol}}_{\text{vol}}(\theta^{en}) \\ & + \Delta t \Psi^{\text{vol}}_{\text{vol}}(\Delta \theta^e, \Delta t), \quad \text{and} \end{aligned} \quad (54) $$

$$ \begin{aligned} \Delta D^{\text{eff}}_{\text{dev}}(\hat{\mathbf{F}}^{n+1}) ={}& \min_{\varepsilon^{\text{vol}n+1}\mathbf{N}} [\Delta D_{\text{dev}}(\hat{\mathbf{F}}, \varepsilon^{\text{pl}}, \mathbf{N})] \\ & + \Delta t \Psi^{\text{dev}}_{\text{dev}}(\Delta \hat{\mathbf{C}}^{\text{vol}}, \Delta t). \end{aligned} \quad (55) $$

## 3.3. Quasi-incompressible form of the weak formulation

Now that an effective internal energy increment has been defined, the volume deformation and the pressure can be computed as Simo and Taylor [30] did for an elastic material. A new three-field internal energy $W^e_{\text{int}} : \mathcal{X}_h \times \mathbb{R}_+ \times \mathbb{R}_+ \to \mathbb{R}$ is defined on the time step, such that
$$W_{\text{int}}^{e}(\vec{\phi}_{h}, \theta^{e}, p^{e}) \equiv \int_{\gamma^{e}} \{\Delta D^{\text{eff}}\vec{x}_{0}, \overline{\mathbf{C}}(\vec{\phi}_{h}, \theta^{e}) + p^{e}[J - \theta^{e}]\} d\mathcal{V}, \quad (56)$$

where $p^e$ is the constant thermodynamic force associated to $\theta^e$. It will be shown that this force is actually the constant pressure. Let us assume that $\delta\vec{\phi}_h(t^n) = \delta\vec{\phi}_h(t^{n+1})$. Since we are considering a Galerkin method, the same restrictions than for the displacements set (39) apply to the virtual displacement, with

$$\mathcal{D}_h \equiv \{\delta\vec{\phi}_h \in \mathcal{D} | [\delta\vec{\phi}_h \in C^0(\mathcal{V}), \delta\vec{\phi}_h \in \mathbb{P}^k(\mathcal{V}^e), \\ \text{and } \delta\vec{\phi}_h(t^n) = \delta\vec{\phi}_h(t^{n+1}) = 0]\}. \quad (57)$$

Therefore, the discretized form of the weak formulation (10) is the stationary point of the functional $I: \mathcal{X}_h \times \prod_e \mathbb{R}_+ \times \prod_e \mathbb{R}_+ \to \mathbb{R}$ such that

$$I(\vec{\phi}_h, (\theta^1, \dots), (p^1, \dots)) = \int_{t^n}^{t^{n+1}} \left\{ \sum_e \int_{\gamma_0^e} [\Delta D^{\text{eff}}(\vec{x}_0, \overline{\mathbf{C}}(\vec{\phi}_h, \theta^e)) + p^e[J - \theta^e]] d\mathcal{V}_0 - \int_{\gamma_0^e} \left[ \frac{1}{2} \rho_0 \dot{\vec{\phi}}_h^2 + \rho_0 \vec{b} \cdot \vec{\phi}_h \right] d\mathcal{V}_0 + \int_{\partial_N \mathcal{V}^e} \overline{\vec{T}} \cdot \vec{\phi}_h d\partial \mathcal{V}^e \right\} dt. \quad (58)$$

Since neither $K$, nor $W_{\text{ext}}$ depend on $p^e$, the variational principle leads to

$$\theta^e = \frac{1}{|\mathcal{V}_0^e|} \int_{\gamma_0^e} J d\mathcal{V}_0 \quad \text{and to} \quad (59)$$

$$p^e = \frac{1}{|\mathcal{V}_0^e|} \int_{\gamma_0^e} \frac{\partial D^{\text{eff}}(\vec{x}_0, \overline{\mathbf{F}}(\vec{\phi}, \theta^e))}{\partial \theta^e} d\mathcal{V}_0, \quad (60)$$

which gives an explicit dependency on $\vec{\phi}$ for $\theta^e$ and for $p^e$, where $|\mathcal{V}^e|$ is the volume measure of $\mathcal{V}^e$. The last variational principle leads to the new weak formulation by using the relations

$$\int_{t^n}^{t^{n+1}} \int_{\gamma_0^e} \rho_0 \ddot{\vec{\phi}}_h \cdot \delta \vec{\phi}_h \, d\mathcal{V}_0 \, dt = \underbrace{\left[ \int_{\gamma_0^e} \rho_0 \dot{\vec{\phi}}_h \cdot \delta \vec{\phi}_h \, d\mathcal{V}_0 \right]}_{=0}^{t^{n+1}} - \delta \int_{t^n}^{t^{n+1}} \int_{\gamma_0^e} \frac{1}{2} \rho_0 \dot{\vec{\phi}}_h^2 \, d\mathcal{V}_0 \, dt, \quad \text{and} \quad (61)$$

$$\int_{t^n}^{t^{n+1}} \int_{\gamma_0^e} [\mathbf{F}] : \vec{\nabla} \delta \vec{\phi}_h \, d\mathcal{V}_0 \, dt = \int_{t^n}^{t^{n+1}} \int_{\gamma_0^e} \frac{\partial \Delta D^{\text{eff}}}{\partial \mathbf{F}} : \delta \mathbf{F} \, d\mathcal{V}_0 \, dt = \delta \int_{t^n}^{t^{n+1}} \int_{\gamma_0^e} \Delta D^{\text{eff}} \, d\mathcal{V}_0 \, dt. \quad (62)$$

Indeed, since Gâteau derivatives lead to

$$\begin{align} \frac{\partial J}{\partial \vec{\phi}} \delta \vec{\phi} &= J \operatorname{tr}(\vec{\nabla} \delta \vec{\phi}) \quad \text{and to} \nonumber \\ \frac{\partial \overline{\mathbf{C}}}{\partial \vec{\phi}} \cdot \delta \vec{\phi} &= 2 \overline{\mathbf{F}}^T \vec{\nabla} \delta \vec{\phi}^T \overline{\mathbf{F}} - \frac{2}{3} \overline{\mathbf{C}} \operatorname{tr}(\vec{\nabla} \delta \vec{\phi}), \tag{63} \end{align}$$

the new weak formulation resulting from Eq. (58) now states as finding $\vec{\phi}_h \in \mathcal{X}_h$ such that

$$\begin{align} & \sum_e \int_{\gamma_0^e} \mathbf{F} \left[ 2J^{-\frac{3}{2}} \mathrm{DEV} \left( \frac{\partial \Delta D^{\mathrm{eff}}}{\partial \hat{\mathbf{C}}} \right) + p^e \mathbf{J} \mathbf{C}^{-1} \right] : \vec{\nabla}_0 \delta \vec{\phi}_h d\mathcal{V}_0 \nonumber \\ & \qquad + \int_{\gamma_0^e} \rho_0 \ddot{\vec{\phi}}_h \cdot \delta \vec{\phi}_h \, d\mathcal{V}_0 \nonumber \\ &= \int_{\partial_N \mathcal{V}^e} \delta \vec{\phi}_h \cdot \overline{\vec{T}} \, d\partial \mathcal{V}_0 + \int_{\gamma_0^e} \rho_0 \vec{b} \cdot \delta \vec{\phi}_h \, d\mathcal{V}_0 \nonumber \\ & \forall \delta \vec{\phi}_h \in \mathcal{D}_h, \quad \forall t \in [t^n, t^{n+1}], \tag{64} \end{align}$$

where $\mathrm{DEV}(\bullet) = \bullet - \frac{1}{3}\bullet : \mathbf{CC}^{-1}$ is the deviatoric part of a tensor and where $p^e$ appears clearly as the pressure field.

From these definitions, the stress tensor used to compute the internal forces (43) can be computed as

$$\mathbf{S}^{n+1} = \mathbf{S}_{\text{dev}}^{n+1} + p^{e n+1} J^{n+1} \mathbf{C}^{n+1-1} \quad \text{with} \quad (65)$$

$$p^e = \frac{\partial \Delta D^{\text{vol}}^{\text{eff}}}{\partial \theta}, \quad \text{and} \quad (66)$$

$$\mathbf{S}_{\text{dev}} = 2J^{-\frac{3}{2}}\text{DEV} \left( \mathbf{F}^{\text{pl}-1} \left[ \frac{\partial \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{\mathbf{C}}^{\text{el}}} + \Delta t \frac{\partial \Psi_{\text{dev}}^{\text{d}}}{\partial \hat{\mathbf{C}}^{\text{el}}} \right] \mathbf{F}^{\text{pl}-T} \right), \quad (67)$$

where it has been taken into account that, after the minimization process, the only explicit dependency on $\overline{\mathbf{F}}$ for $\Delta D_{\text{dev}}^{\text{eff}}$ is due to $\Phi_{\text{dev}}^{\text{el}}$ and to $\Psi_{\text{dev}}^{\text{d}}$.

### 3.4. Form of the potentials

Apart from convexity requirement, there is no a priori restriction on the potentials expressions. Nevertheless, in this work bi-logarithmic expressions are considered for the elastic energies, leading to

$$\Phi_{\text{vol}}^{\text{el}}(\theta^e) = \frac{K_0}{2} \ln^2(\theta^e) \quad \text{and to} \quad \Phi_{\text{dev}}^{\text{el}}(\hat{\mathbf{C}}^e) = \frac{G_0}{4} \| \ln(\hat{\mathbf{C}}^e) \|^2. \quad (68)$$

In these expressions, $K_0$ and $G_0$ are respectively the initial bulk and shear modulus. The bi-logarithmic form of the potentials is particularly well suited to get a closed form of the minimization process [19,21]. This formalism is completed by defining the form of the power-law dissipation and isotropic plastic-hardening, which are respectively

$$\Psi^* = \begin{cases} \frac{m Y_0^* e^{\vec{p}_0}}{m+1} \left[ \frac{\dot{\vec{p}}}{e^{\vec{p}_0}} \right]^{\frac{m+1}{m}} & \text{if } \dot{\vec{p}} \ge 0 \\ \infty & \text{if } \dot{\vec{p}} < 0 \end{cases} \quad \text{and} \\ \Phi^{\text{pl}} = \sigma_0^{\text{vm}} e^{\vec{p}} + \frac{n Y_0 e_0^P}{n+1} \left[ \frac{e^{\vec{P}}}{e_0^P} \right]^{\frac{n+1}{n}}, \quad (69)$$

where $Y_0^*, Y_0, \dot{e}_0^P, e_0^P, m$ and $n$ are constants and where $\sigma_0^{\text{vm}}$ is the initial yield stress.
Finally the expressions of the numerical dissipation considered are

$$ \Psi_{\text{vol}}^{d}(\Delta\theta^{e}, \Delta t) = \frac{\chi}{2\Delta t} \frac{\partial^2 \Phi_{\text{vol}}^{\text{el}}(\theta^{\text{en}})}{\partial \theta^{e^2}} [\theta^{\text{en}+1} - \theta^{\text{en}}]^2, \quad \text{and} \quad (70) $$

$$ \begin{aligned} \Psi_{\text{dev}}^{d}(\Delta \hat{\mathbf{C}}^{\text{el}}, \Delta t) = & \frac{\chi}{2\Delta t} [\hat{\mathbf{C}}^{\text{el}n+1} - \hat{\mathbf{C}}^{\text{el}n}] : \frac{\partial^2 \Phi_{\text{dev}}^{\text{el}}(\hat{\mathbf{C}}^{\text{el}n})}{\partial \hat{\mathbf{C}}^{\text{el}} \partial \hat{\mathbf{C}}^{\text{el}}} \\ & : [\hat{\mathbf{C}}^{\text{el}n+1} - \hat{\mathbf{C}}^{\text{el}n}]. \end{aligned} \quad (71) $$

The proposed numerical scheme is thus similar in form to the addition of an artificial viscosity proportional to $\Delta t$ (or inversely proportional to the frequency $\omega$ in a continuous setting). These specific expressions will ensure verification of consistency conditions.

The choice of evaluating the second derivatives in the previous configuration simplifies the computation of the analytical stiffness matrix.

This particular choice will lead to a first-order time-integration scheme, as it will be shown. Parameter $\chi$ is a user parameter that controls the numerical dissipation. Similarly, a second-order time integration scheme can be obtained by evaluating the numerical dissipation at an intermediate configuration $\vec{\phi}^* = \vec{\phi}_h^n + \chi \Delta t (\dot{\vec{\phi}}^* - \dot{\vec{\phi}}_h^n)$. In a way similar to the one Armero and Romero [26,27] have proposed for elasto-dynamics, in this formalism, both $\vec{\phi}^*$ and $\vec{\phi}^*$ constitute new unknowns solved for each elements in order to enforce second-order accuracy. These developments are beyond the scope of this work and will be pursued in future work.

## 4. Time integration scheme

The finite-elements discretization combined with the variational updates formulation presented in Section 3 is now integrated in time with the EDMC algorithm. It is demonstrated that, in combination with the EDMC scheme, the Gonzalez' expression [11] of the modified second Piola-Kirchhoff stress tensor can be used directly from this incremental functional in order to verify the conservation laws. The study of the spectral matrix associated to this time-integration method is also considered.

### 4.1. The EDMC algorithm

When introducing numerical dissipation in the EMCA, Armero and Romero [26,27] have avoided bifurcation in the spectral analysis of the amplification matrix, by considering velocities dissipation $\vec{v}_{\text{diss}}$. Therefore, the relations between positions and velocities and accelerations are given by a mid-point approximation with

$$ \vec{x}^{\text{d}n+1} = \vec{x}^{\text{d}n} + \Delta t \left[ \frac{\dot{\vec{x}}^{\text{d}n+1} + \dot{\vec{x}}^{\text{d}n}}{2} + \vec{v}_{\text{diss}}^{\text{d}n+1} \right], \quad \text{and} \quad (72) $$

$$ \dot{\vec{x}}^{\text{d}n+1} = \dot{\vec{x}}^{\text{d}n} + \Delta t \frac{\ddot{\vec{x}}^{\text{d}n+1} + \ddot{\vec{x}}^{\text{d}n}}{2}. \quad (73) $$

This relation is a second-order approximation in $\Delta t$ if $\vec{v}_{\text{diss}} = \mathcal{O}(\Delta t^2)$ and is a first-order approximation if $\vec{v}_{\text{diss}} = \mathcal{O}(\Delta t)$. In this work, we use the first-order dissipation velocities used in [28], which are

$$ \vec{v}_{\text{diss}}^{\text{d}n+\frac{1}{2}} = \chi \frac{\|\dot{\vec{x}}_s^{\text{d}n+1}\| - \|\dot{\vec{x}}_s^{\text{d}n}\|}{\|\dot{\vec{x}}_s^{\text{d}n+1}\| + \|\dot{\vec{x}}_s^{\text{d}n}\|} \frac{\dot{\vec{x}}_s^{\text{d}n+1} + \dot{\vec{x}}_s^{\text{d}n}}{2}, \quad (74) $$

where $\chi$ is a user chosen control parameter. Its interval range is determined later owing a spectral analysis.

Balance Eq. (41) is discretized in time by

$$ M^{\vec{\epsilon}\mu} \frac{\ddot{\vec{x}}_s^{\mu n} + \dddot{\vec{x}}_s^{\mu n+1}}{2} = \vec{f}_{\text{ext}}^{\vec{\epsilon}n+\frac{1}{2}} - \vec{f}_{\text{int}}^{\vec{\epsilon}n+\frac{1}{2}}, \quad (75) $$

where $\vec{f}_{\text{ext}}^{\vec{\epsilon}n+\frac{1}{2}}$ and $\vec{f}_{\text{int}}^{\vec{\epsilon}n+\frac{1}{2}}$ are approximations of respectively the external and internal forces designed in a suitable way to respect the conservation laws. The set of Eqs. (72)-(75) is solved by a predictor-corrector algorithm, enhanced with a line-search method, see [28] for details.

The external forces formulation depends on the applied loads and is beyond the scope of this paper. In the context of Energy-Momentum Conserving Algorithm, Gonzalez [11] proposed a new expression of the internal forces. Given the incremental updates variational formulation developed in the previous section, Gonzalez expression can be directly applied, leading to

$$ \vec{f}_{\text{int}}^{\vec{\epsilon}n+\frac{1}{2}} = \int_{\gamma_0} \left\{ \frac{\mathbf{F}^{n+1} + \mathbf{F}^n}{2} [\mathbf{S}_{\text{dev}}^{*n+\frac{1}{2}} + 2p^{*n+\frac{1}{2}} \mathbf{dG}] \vec{D}^{\vec{\epsilon}} \right\} d\gamma_0, \quad (76) $$

where $\mathbf{S}_{\text{dev}}^*$ is the consistent deviatoric stress tensor, $p^*$ is the consistent pressure and $\mathbf{dG}$ is the modified differentiation of $J$ by $\mathbf{C}$. This last term is defined by

$$ \mathbf{dG} = \frac{\partial J^{n+\frac{1}{2}}}{\partial \mathbf{C}} + \left[ \frac{J^{n+1} - J^n - \frac{\partial J^{n+\frac{1}{2}}}{\partial \mathbf{C}} : \Delta \mathbf{C}}{\|\Delta \mathbf{C}\|^2} \right] \Delta \mathbf{C} \quad \text{with} \quad (77) $$

$$ \frac{\partial J^{n+\frac{1}{2}}}{\partial \mathbf{C}} = \frac{1}{2} \sqrt{\det \left( \frac{\mathbf{C}^{n+1} + \mathbf{C}^n}{2} \right) \left[ \frac{\mathbf{C}^{n+1} + \mathbf{C}^n}{2} \right]^{-1}}, \quad \text{and} \quad (78) $$

$$ \Delta \mathbf{C} = \mathbf{C}^{n+1} - \mathbf{C}^n. \quad (79) $$

Instead of using directly Eq. (66), the modified pressure is computed as

$$ \begin{aligned} p^{\text{d}n+\frac{1}{2}} = & \left[ \frac{\partial \Delta D_{\text{vol}}^{\text{eff}}}{\partial \theta^e} \right]^{n+\frac{1}{2}} \\ & + \frac{\Delta D_{\text{vol}}^{\text{eff}}(\theta^{e^{n+1}}, \theta^{e^n}) - \left[ \frac{\partial \Delta D_{\text{vol}}^{\text{eff}}}{\partial \theta^e} \right]^{n+\frac{1}{2}}}{\Delta \theta^{e^2}} \Delta \theta^e, \quad \text{with} \quad (80) \end{aligned} $$

$$ \left[ \frac{\partial \Delta D_{\text{vol}}^{\text{eff}}}{\partial \theta^e} \right]^{n+\frac{1}{2}} = \frac{\partial \Phi_{\text{vol}}^{\text{el}}}{\partial \theta^e} \left( \frac{\theta^{e^{n+1}} + \theta^{e^n}}{2} \right) + \Delta t \frac{\partial \Psi_{\text{vol}}^{\text{d}}}{\partial \Delta \theta^e} \left( \frac{\Delta \theta^e}{2}, \Delta t \right), $$

$$ \Delta \theta^e = \theta^{e^{n+1}} - \theta^{e^n}. \quad (82) $$

Similarly, Eq. (67) becomes
$$
S_{\text{dev}}^{n+\frac{1}{2}} = 2 \left[ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}}}{\partial \mathbf{C}} \right]^{n+\frac{1}{2}} \\
+ 2 \frac{\Delta D_{\text{dev}}^{\text{eff}} (\hat{\mathbf{C}}^{n+1}, \hat{\mathbf{C}}^n) - \left[ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}}}{\partial \mathbf{C}} \right]^{n+\frac{1}{2}}}{\| \Delta \mathbf{C} \|^2} : \Delta \mathbf{C}, \quad \text{with} \tag{83}
$$

$$
\left[ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}}}{\partial \mathbf{C}} \right]^{n+\frac{1}{2}} = (J^{n+\frac{1}{2}})^{-\frac{2}{3}} \text{DEV} \left( (\mathbf{F}^{\text{pl}n+\frac{1}{2}})^{-1} \left[ \frac{\partial \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{\mathbf{C}}^{\text{el}}} (\hat{\mathbf{C}}^{\text{el}n+\frac{1}{2}}) \right. \right. \\
\left. \left. + \Delta t \frac{\partial \mathbf{\Psi}_{\text{dev}}^{\text{d}}}{\partial \hat{\mathbf{C}}^{\text{el}}} (\hat{\mathbf{C}}^{\text{el}n+\frac{1}{2}} - \hat{\mathbf{C}}^{\text{el}n}, \Delta t) \right] (\mathbf{F}^{\text{pl}n+\frac{1}{2}})^{-T} \right), \quad (84)
$$

$$
\text{and } \Delta \mathbf{C} = \mathbf{C}^{n+1} - \mathbf{C}^n. \tag{85}
$$

In relation (84), $J^{n+\frac{1}{2}} = \det(\frac{\mathbf{C}^{n+1}+\mathbf{C}^n}{2})$, while terms $\mathbf{F}^{\text{pl}n+\frac{1}{2}}$, and $\hat{\mathbf{C}}^{\text{el}n+\frac{1}{2}}$ result from the minimization of $\Delta D$, Eq. (50), with the input tensor $\frac{\mathbf{C}^{n+1}+\mathbf{C}^n}{2}$.

Let us note that since the numerical dissipation is introduced into the effective potential, the approximations $p^* \simeq \left[ \frac{\partial \Delta D_{\text{vol}}^{\text{eff}}}{\partial \theta^e} \right]^{n+\frac{1}{2}}$ and $S_{\text{dev}}^* \simeq \left[ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}}}{\partial \bar{\mathbf{C}}^e} \right]^{n+\frac{1}{2}}$ are exact to the first-order as it was demonstrated by Gonzalez [11]. This would not be the case if the numerical dissipation was introduced only in the correcting terms. In such a case the correcting terms would be of the same order as the main terms.

Moreover, this formulation respects the conservation laws established in Section 2.4. Indeed, since $\sum_{\xi} \tilde{D}^{\xi} = 0$, applying a summation on the nodes in Eq. (75), and using Eq. (73) leads to

$$
\sum_{\xi} M^{\xi\mu} \dot{\mathbf{x}}^{\mu n+1} - \sum_{\xi} M^{\xi\mu} \dot{\mathbf{x}}^{\mu n} = \Delta t \sum_{\xi} \vec{f}_{\text{ext}}^{\xi^{n+1}}, \quad (86)
$$

which is the discrete expression of the law (33).

Using Eqs. (72)–(74), the vector product of Eq. (75) with $\vec{x}^{\xi^{n+1}/2} = \frac{\vec{x}^{\xi^{n+1}} + \vec{x}^n}{2}$ leads to

$$
\begin{align}
& M^{\xi\mu} \vec{x}^{\xi^{n+1}} \wedge \dot{\vec{x}}^{\mu n+1} - M^{\xi\mu} \vec{x}^{\xi^{n}} \wedge \dot{\vec{x}}^{\mu n} \nonumber \\
& \qquad = \Delta t \vec{x}^{\xi^{n+1}/2} \wedge [\vec{f}_{\text{ext}}^{\xi^{n+1}/2} - \vec{f}_{\text{int}}^{\xi^{n+1}/2}]. \tag{87}
\end{align}
$$

But if $\epsilon$ is the permutation tensor, and since $\mathbf{dG}$ and $\mathbf{S}_{\text{dev}}$ are symmetric, one has

$$
\vec{x}^{\xi^{n+1}/2} \wedge \vec{f}_{\text{int}}^{\xi^{n+1}/2} = \int_{\mathcal{Y}_0} \left\{ \epsilon : \left[ \frac{\mathbf{F}^{n+1} + \mathbf{F}^n}{2} [\mathbf{S}_{\text{dev}}^{n+1/2} + 2p^{n+1/2} \mathbf{dG}] \frac{\mathbf{F}^{n+1\mathrm{T}} + \mathbf{F}^{n\mathrm{T}}}{2} \right] \right\} d\mathcal{Y}_0 = \vec{0}, \quad (88)
$$

and Eq. (87) can be rewritten

$$
M^{\xi\mu} \vec{x}^{\xi^{n+1}} \wedge \dot{\vec{x}}^{\mu n+1} - M^{\xi\mu} \vec{x}^{\xi^{n}} \wedge \dot{\vec{x}}^{\mu n} = \Delta t \vec{x}^{\xi^{n+1}/2} \wedge \vec{f}_{\text{ext}}^{\xi^{n+1}/2}, \quad (89)
$$

which is the discrete expression of the law (34). Similarly, using Eqs. (72) and (73), the scalar product of Eq. (75) with $\vec{x}^{\xi^{n+1}} - \vec{x}^{\xi^n}$ leads to

$$
\begin{align}
& \frac{1}{2} M^{\xi\mu} \dot{\vec{x}}^{\xi^{n+1}} \cdot \dot{\vec{x}}^{\mu n+1} - \frac{1}{2} M^{\xi\mu} \dot{\vec{x}}^{\xi^n} \cdot \dot{\vec{x}}^{\mu n} + M^{\xi\mu} [\vec{x}^{\mu n+1} - \vec{x}^{\mu n}] \cdot \vec{v}_{\text{diss}}^{\xi^{n+1/2}} \nonumber \\
& = [\vec{f}_{\text{ext}}^{\xi^{n+1/2}} - \vec{f}_{\text{int}}^{\xi^{n+1/2}}] \cdot [\vec{x}^{\xi^{n+1}} - \vec{x}^{\xi^n}]. \tag{90}
\end{align}
$$

Regarding the dissipation velocities, Eq. (74) leads to

$$
\begin{align}
M^{\xi\mu}[\vec{x}^{\mu n+1} - \vec{x}^{\mu n}] \cdot \vec{v}_{\text{diss}}^{\xi^{n+1/2}} &\approx \frac{\chi}{2} m^{\xi} [\|\dot{\vec{x}}^{\xi^{n+1}}\| - \|\dot{\vec{x}}^{\xi^n}\|^2] \nonumber \\
&= \Delta W_{\text{diss}}^{K}, \tag{91}
\end{align}
$$

where we have assumed a lumped mass matrix $M^{\xi\mu} = m^\xi \delta_{\xi\mu}$ (no sum on $\xi$), and where $\Delta W_{\text{diss}}^K \ge 0$ is the part of the numerical dissipation $W_{\text{diss}}$ coming from the dissipation velocities.

On the other hand, combination of Eqs. (76)–(85) gives

$$
\begin{align}
\vec{f}_{\text{int}}^{\xi^{n+1/2}} \cdot [\vec{x}^{\xi^{n+1}} - \vec{x}^{\xi^n}] &= \int_{\mathcal{Y}_0} \left\{ [\mathbf{S}_{\text{dev}}^{*n+1/2} + 2p^{*n+1/2} \mathbf{dG}] : \frac{\mathbf{\Delta C}}{2} \right\} d\mathcal{Y}_0 \nonumber \\
&= \int_{\mathcal{Y}_0} \left\{ \Delta D_{\text{dev}}^{\text{eff}} (\hat{\mathbf{C}}^{n+1}, \hat{\mathbf{C}}^n) \right. \nonumber \\
&\qquad \left. + p^{*n+1/2} [J^{n+1} - J^n] \right\} d\mathcal{Y}_0. \tag{92}
\end{align}
$$

Introducing the definition (59) in these latest results and remembering that $p^*$ is constant on $\mathcal{Y}^e$ allows to rewrite this last equation as

$$
\begin{align}
W_{\text{int}}^{n+1} - W_{\text{int}}^{n} &= \int_{\mathcal{Y}_0} \Delta D_{\text{dev}}^{\text{eff}} (\hat{\mathbf{C}}^{n+1}, \hat{\mathbf{C}}^n) d\mathcal{Y}_0 \nonumber \\
&\quad + \sum_e p^{*n+1/2} \int_{\mathcal{Y}_0} [J^{n+1} - J^n] d\mathcal{Y}_0 \nonumber \\
&= \int_{\mathcal{Y}_0} \left\{ \Delta D_{\text{dev}}^{\text{eff}} (\hat{\mathbf{C}}^{n+1}, \hat{\mathbf{C}}^n) \right. \nonumber \\
&\qquad \left. + \Delta D_{\text{vol}}^{\text{eff}} (\theta^{e n+1}, \theta^{e n}) \right\} d\mathcal{Y}_0. \tag{93}
\end{align}
$$

Time discretization of Eqs. (32), (36) and (37) reads

$$
\begin{equation}
\begin{split}
U_{\text{el}}^{n+1} - U_{\text{el}}^{n} = \int_{\mathcal{Y}_0} & \{ \boldsymbol{\Phi}_{\text{vol}}^{\text{el}}(\theta^{e n+1}) + \boldsymbol{\Phi}_{\text{dev}}^{\text{el}}(\hat{\boldsymbol{C}}^{\text{el} n+1}) \\
& - \boldsymbol{\Phi}_{\text{vol}}^{\text{el}}(\theta^{e n}) - \boldsymbol{\Phi}_{\text{dev}}^{\text{el}}(\hat{\boldsymbol{C}}^{\text{el} n}) \} d\mathcal{Y}_0,
\end{split}
\tag{94}
\end{equation}
$$

$$
W_{\text{diss}}^{W^{n+1}} - W_{\text{diss}}^{W^n} = \Delta t \int_{\mathcal{Y}_0} \left\{ \Psi_{\text{vol}}^{\text{d}}(\theta^{e n+1} - \theta^{e n}, \Delta t) \right. \\
\left. \qquad + \Psi_{\text{dev}}^{\text{d}}(\hat{\mathbf{C}}^{\text{el} n+1} - \hat{\mathbf{C}}^{\text{el} n}, \Delta t) \right\} d\mathcal{Y}_0, \quad \text{and} \tag{95}
$$

$$
U_i^{n+1} - U_i^n = \int_{\mathcal{Y}_0} \left\{ \begin{aligned} & \boldsymbol{\Phi}^{\text{pl}}(\mathbf{F}^{\text{pl}n+1}, \varepsilon^{\text{p}n+1}) - \boldsymbol{\Phi}^{\text{pl}}(\mathbf{F}^{\text{pl}n}, \varepsilon^{\text{p}n}) \\ & + \Delta t \boldsymbol{\Psi}^* \left( \frac{\Delta \varepsilon^{\text{p}}}{\Delta t} \right) \end{aligned} \right\} d\mathcal{Y}_0. \quad (96)
$$

Using these last relations and Eqs. (53)–(55), the work increment of the internal forces (93) can be rewritten

$$
W_{\text{int}}^{n+1} - W_{\text{int}}^{n} = U_{\text{el}}^{n+1} - U_{\text{el}}^{n} + W_{\text{diss}}^{W^{n+1}} - W_{\text{diss}}^{W^n} + U_{i}^{n+1} - U_{i}^{n}, \quad (97)
$$

and, therefore, Eq. (90) becomes

$$
\Delta K + \Delta U_{\text{el}} = \Delta W_{\text{ext}} - \Delta U_{i} - \Delta W_{\text{diss}}, \quad (98)
$$
where $W_{\text{diss}} = W_{\text{diss}}^K + W_{\text{diss}}^W$ is the numerical dissipation and where $K = \frac{1}{2} M \xi^\mu \dot{x}^\mu \dot{x}^\mu$ is the discretized kinetic energy. Relation (98) is the discrete form of the energy conservation (38), which demonstrates the consistency of the proposed formulation.

Let us notice that the irreversible energy (96) dissipated during the time step does not exactly correspond to the work accomplished by the irreversible part of the stress tensor. Indeed, the Perzyna dissipation $\Psi^*$ (69) leads to

$$ \int_{t^n}^{t^{n+1}} \frac{\partial \Psi^* (\varepsilon^{\text{P}})}{\partial \varepsilon^{\text{P}}} dt \neq \Delta t \Psi^* \left( \frac{\Delta \varepsilon^{\text{P}}}{\Delta t} \right), \quad (99) $$

unless $m \to \infty$, see Eq. (96). Moreover, the correcting term of $S_{\text{dev}}^*$ (83) is no longer second order compared to the main derivation $[\frac{\partial \Delta D}{\partial C}]^{n+\frac{1}{2}}$. Nevertheless, since this dissipation of energy is positive, the G-stability of the method is still ensured (the energy of the system $K + U_{\text{el}}$ is either constant or decreasing). In the present work we will assume that there is no viscous dissipation, and since the numerical dissipation $\Psi^d$ does not suffer from this behavior, the formulation remains consistent and the correcting term is still second order compared to the main derivation.

To solve efficiently a time increment with the predictor-corrector scheme developed for EDMC algorithms [28], the consistent stiffness matrix associated to the internal forces formulation (76) has to be defined in a closed form. This is done in Appendix A.

## 4.2. Numerical properties

In this section we demonstrate the first-order convergence of the time integration. Toward this end, the tensile test of elongation $x$ of an elastic beam (length $l$, area $A$, Young's modulus $E$ and Poisson ratio $\nu$) is considered. The deformation gradient associated to this test is

$$ \mathbf{F} = \begin{pmatrix} 1 + \frac{x}{l} & 0 & 0 \\ 0 & 1 - \nu \frac{x}{l} & 0 \\ 0 & 0 & 1 - \nu \frac{x}{l} \end{pmatrix}. \quad (100) $$

As it is shown in Appendix B, the internal forces (76) corresponding to a deformation from $x^n$ to $x^{n+1}$ reads

$$ \vec{f}_{\text{int}}^{n+\frac{1}{2}} = \frac{EA}{l} \begin{pmatrix} \left[ 1 + \chi \frac{\chi^{n+1} - \chi^n}{\chi^{n+1} + \chi^n} \right] \frac{\chi^n + \chi^{n+1}}{2} \\ 0 \\ 0 \end{pmatrix}. \quad (101) $$

Assuming that the velocity keeps the same direction during a time step, the dissipation velocity (74) is rewritten

$$ v_{\text{diss}}^{n+\frac{1}{2}} = \chi \frac{\dot{x}^{n+1} - \dot{x}^n}{\dot{x}^{n+1} + \dot{x}^n} \frac{\dot{x}^{n+1} + \dot{x}^n}{2}. \quad (102) $$

If the mass associated to the degree of freedom is $m = \frac{\rho A l}{2}$ the frequency of the system is $\omega^2 = \frac{2E}{\rho A l^2}$ and the set of Eqs. (72)–(75) are rewritten

$$ x^{n+1} = x^n + \frac{\Delta t}{2} \left[ 1 + \chi \frac{\dot{x}^{n+1} - \dot{x}^n}{\dot{x}^{n+1} + \dot{x}^n} \right] [\dot{x}^{n+1} + \dot{x}^n], \quad (103) $$

$$ \dot{x}^{n+1} = \dot{x}^n + \frac{\Delta t}{2} [\dot{x}^{n+1} + \dot{x}^n], \quad (104) $$

$$ \ddot{x}^{n+1} = -\ddot{x}^n - \omega^2 \left[ 1 + \chi \frac{\dot{x}^{n+1} - \dot{x}^n}{\dot{x}^{n+1} + \dot{x}^n} \right] [\dot{x}^{n+1} + \dot{x}^n]. \quad (105) $$

This system of equations is equivalent to the one obtained by Armero and Romero [26] for a spring model and leads to an error $\Omega_d$ on the adimensional frequency $\Omega = \omega \Delta t$ and to an error $\xi_d$ on the damping ratio, that are respectively

$$ \Omega_d = \Omega - \frac{1}{4} \left[ \chi^2 + \frac{1}{3} \right] \Omega^3 + \mathcal{O}(\Omega^4), \quad \text{and} \quad (106) $$

$$ \xi_d = \frac{\chi}{2} \Omega + \mathcal{O}(\Omega^2). \quad (107) $$

Due to the presence of the dissipation parameter $\chi$, the scheme is first-order accurate. Moreover the spectral radius $\rho_d$ of the amplification matrix associated to this scheme is

$$ \rho_d(\Omega) = \frac{1}{1 + \frac{\Omega^2}{4} [1 + \chi]^2} \sqrt{\left[1 - \frac{\Omega^2}{4} [1 - \chi^2]\right]^2 + \Omega^2}. \quad (108) $$

If $0 \le \chi \le 1$ the spectral radius is always lower than unity and is decreasing with $\Omega$, demonstrating the absence of bifurcation. It is close to unity for the low frequencies and tends to its minimal value $\rho_\infty = \frac{1-\chi}{1+\chi}$ for the high frequencies.

# 5. Numerical examples

The accuracy of the developed scheme is demonstrated by considering several examples involving highly non-linear behavior. These non-linearities are either geometrical, or material due to the large deformation or the elasto-plastic behavior.

## 5.1. Tumbling L-shaped block

This example is particularly well suited to verify the conservation of the angular momentum and energy. It consists into an elasto-plastic L-shaped block subjected to an initial loading. Its geometry is described in Fig. 1 and the material properties are reported in Table 1. On face ‘a’, the forces, depending on time $t$, applied at every nodes, are

$$ \vec{f}_{\text{ext}}^{a}(t) = \begin{pmatrix} 4 \\ 8 \\ 12 \end{pmatrix} \text{ N s}^{-1} \times \begin{cases} t & 0 \le t \le 2.5 \text{ s} \\ (5-t) & 2.5 < t \le 5 \text{ s} \end{cases}, \quad (109) $$

while the opposite forces are applied on face ‘b’ ($\vec{f}_{\text{ext}}^{b}(t) = -\vec{f}_{\text{ext}}^{a}(t)$). After 5 s, the forces are relaxed and the block starts tumbling.

Let us solve this problem with the developed EDMC algorithm. A spectral radius $\rho_\infty = 0.7$ (leading to $\chi \approx 0.17$) is considered, while a first simulation with a time step $\Delta t = 0.25 \text{ s}$ is performed. Fig. 2 shows some snapshots
of the simulation. It can be seen that the block is plastically and elastically deformed under the initial loads, and that it enters into rotation. Particularly, the plastic deformations

<figure><img src="image_2.png" /><figcaption>Fig. 1. Geometry of the L-shaped block (m).</figcaption></figure>

do not increase indefinitely with time, as it would be the case if instabilities occured. Time evolution of the angular momenta is illustrated in Fig. 3a and it appears that they are preserved after the initial loads are released, in accordance with the theory. In the same way, the energy conservation is illustrated in Fig. 3b. It can be seen that the total energy of the system ($K + U_{el} + U_i$) decreases from the initial energy introduced in the system by the external forces, to a constant corresponding to the sum of the kinetic

<table><caption>Table 1
Material properties for the tumbling of a L-shaped block</caption><thead><tr><th>Property</th><th>Value</th></tr></thead><tbody><tr><td>Density</td><td>ρ = 100 kg m<sup>-3</sup></td></tr><tr><td>Bulk modulus</td><td>K<sub>0</sub> = 5000 N m<sup>-2</sup></td></tr><tr><td>Shear modulus</td><td>G<sub>0</sub> = 1000 N m<sup>-2</sup></td></tr><tr><td>Yield stress</td><td>σ<sub>0</sub><sup>vm</sup> = 300 N m<sup>-2</sup></td></tr><tr><td>Hardening exponent</td><td>n = 1</td></tr><tr><td>Linear hardening</td><td>h = Y<sub>0</sub>/ε<sub>0</sub><sup>p</sup> = 400 N m<sup>-2</sup></td></tr><tr><td>Perzyna exponent</td><td>m = ∞</td></tr><tr><td>Perzyna hardening</td><td>Y<sub>0</sub><sup>*</sup> = 0</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Snapshots of deformed geometry and equivalent plastic strain for the tumbling L-shaped block.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Verification of the conservation laws for the tumbling of a L-shaped block.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Influence of the time step size on the energy of the system for the tumbling of a L-shaped block.</figcaption></figure>

energy of the rigid rotation with the energy plastically dissipated. Similarly, time evolution of the kinetic energy $K$ and of the elastic energy $U_{el}$ exhibit a decrease towards the rigid motion energy and towards zero respectively. The annihilation of the oscillations with time is clearly illustrated in this picture. Influence of the time step size is shown in Fig. 4. The time step considered increased from 0.25 s to 1 s. For this last case, the Newmark time-integration scheme leading to the same spectral radius $\rho_\infty$ is also considered. Regarding the EDMC algorithm, if the time step size increases, the elastic energy of the system is dissipated faster (Fig. 4a), and the total energy of the system tends toward the same rigid body motion energy. This behavior is in accordance with the theory, since a larger time step size correspond to a larger adimensional frequency $\Omega$ and therefore to a larger damping ratio (107). The asymptotic behavior of the Newmark simulation is not the rigid rotation of the system but corresponds to a total dissipation of the energy. Moreover, this dissipation is much more important than the dissipation of the EDMC algorithm. The resulting error in the plastically dissipated energy (Fig. 4b) is therefore more important for the Newmark algorithm (>400%) than for the EDMC scheme (25%), when a time step size of 1 s instead of 0.25 s is considered.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Geometry of the bar for the Taylor's impact test (mm).</figcaption></figure>
<table><caption>Table 2
Material properties for the Taylor's impact test</caption><thead><tr><th>Property</th><th>Value</th></tr></thead><tbody><tr><td>Density</td><td>ρ = 8930 kg m<sup>-3</sup></td></tr><tr><td>Bulk modulus</td><td>K<sub>0</sub> = 130,000 N mm<sup>-2</sup></td></tr><tr><td>Shear modulus</td><td>G<sub>0</sub> = 43,333 N mm<sup>-2</sup></td></tr><tr><td>Yield stress</td><td>σ<sub>0</sub><sup>ym</sup> = 400 N mm<sup>-2</sup></td></tr><tr><td>Hardening exponent</td><td>n = 1</td></tr><tr><td>Hardening</td><td>h = Y<sub>0</sub>/c<sub>0</sub><sup>P</sup> = 100 N mm<sup>-2</sup></td></tr><tr><td>Perzyna exponent</td><td>m = ∞</td></tr><tr><td>Perzyna hardening</td><td>Y<sub>0</sub><sup>*</sup> = 0</td></tr></tbody></table>

## 5.2. Taylor's bar impact

This test has for purpose to demonstrate the accuracy of the scheme when large plastic deformations are involved.

The bar described in Fig. 5 impacts a rigid wall with an initial velocity of 227 m s⁻¹. Material properties are reported in Table 2.

First, this simulation is solved with a time step $Δt = 0.4$ µs, and solutions obtained with the EDMC algorithm for different values of the spectral radius are compared to the solutions obtained with the Newmark [1] and the Chung-Hulbert [24] time integration algorithms. The importance of considering the elastic deformation when computing the dissipation potential $Ψ_{dev}^d$ (71) is also illustrated by solving this test while considering $Ψ_{dev}^d = \frac{\gamma}{2Δt}ΔC : \frac{\partial^2 φ_{dev}^e}{\partial C^e \partial C^e} : ΔC$. This last case is referenced to as modified EDMC. Final deformations are illustrated in Fig. 6, and final geometries are reported in Table 3. All the simulations give the same solution within 5% error, except the

<figure><img src="image_7.png" /><figcaption>Fig. 6. Comparison of the final equivalent plastic strain ($t = 80$ µs, $Δt = 0.4$ µs) for Taylor's bar impact using: (a)–(c) EDMC scheme with a spectral radius evolving from $ρ_∞ = 0.9$ to $ρ_∞ = 0.65$; (d) EDMC scheme but with the numerical dissipation $Ψ_{dev}^d$ considering the total deformations instead of the elastic ones as in Eq. (71); (e) second-order Chung-Hulbert time integration scheme; (f) first-order Newmark time-integration scheme.</figcaption></figure>
<table><caption>Table 3
Final configuration ($t = 80$ µs, $\Delta t = 0.4$ µs) for the Taylor's impact test</caption><thead><tr><th>Scheme</th><th>Max<br>e<sup>P</sup></th><th>Radius<br>(mm)</th><th>Length<br>(mm)</th></tr></thead><tbody><tr><td>EDMC; Δt = 0.4; ρ<sub>∞</sub> = 0.9</td><td>2.69</td><td>6.87</td><td>21.46</td></tr><tr><td>EDMC; Δt = 0.4; ρ<sub>∞</sub> = 0.7</td><td>2.76</td><td>6.91</td><td>21.50</td></tr><tr><td>EDMC; Δt = 0.4; ρ<sub>∞</sub> = 0.65</td><td>2.77</td><td>6.91</td><td>21.51</td></tr><tr><td>Modified EDMC; Δt = 0.4;<br>ρ<sub>∞</sub> = 0.7</td><td>1.63</td><td>5.78</td><td>31.45</td></tr><tr><td>CH; Δt = 0.4; ρ<sub>∞</sub> = 0.7</td><td>2.79</td><td>6.91</td><td>21.51</td></tr><tr><td>Newmark; Δt = 0.4; ρ<sub>∞</sub> = 0.7</td><td>2.79</td><td>6.85</td><td>21.56</td></tr></tbody></table>

modified EDMC one, which leads to a different solution. This demonstrates the accuracy of the EDMC algorithms when the dissipation potential (71) is considered.

The order of convergence is studied by considering the plastically dissipated energy $U_i$ at the end of simulation, for different time step sizes. These energies are reported in Fig. 7a, showing that almost the whole initial kinetic energy has been dissipated by plasticity. All the schemes lead to a similar solution, whatever the time step size. The corresponding errors are reported in Fig. 7b. All the

schemes are first-order accurate. This is in accordance with the theory for the EDMC and for the Newmark algorithms, but in contradiction to the theory for the Chung-Hulbert scheme. This contradiction is justified by the fact that the theoretical analysis was assuming a linear behavior, which is not the case in the present simulation.

### 5.3. Lateral impact of a square tube

The purpose of this example is to emphasize the robustness of the scheme by considering a more elaborated benchmark. It consists of the side impact of the square tube illustrated in Fig. 8, on a rigid cylindrical punch. Taking advantage of the symmetry of the problem, only one fourth of the structure is simulated, with appropriate boundary conditions on the symmetry sides. The tube is discretized with 30 elements on its half length, with a non-uniform distribution (elements close to the punch are three times smaller than elements on the other extremity), with 15 elements on its height (out of the curved corners) and with 5 elements on each corner curves. The tube is made of Aluminum (properties reported in Table 4), modeled by a

<figure><img src="image_8.png" /><figcaption>Fig. 7. Plastically dissipated energy for the Taylor's impact test.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Geometry of the square tube lateral impact test (mm).</figcaption></figure>
<table><caption>Table 4
Material properties for the square-tube lateral impact test</caption><thead><tr><th>Property</th><th>Value</th></tr></thead><tbody><tr><td>Density</td><td>ρ = 2700 kg m<sup>-3</sup></td></tr><tr><td>Bulk modulus</td><td>K<sub>0</sub> = 59167 N mm<sup>-2</sup></td></tr><tr><td>Shear modulus</td><td>G<sub>0</sub> = 27308 N mm<sup>-2</sup></td></tr><tr><td>Yield stress</td><td>σ<sub>0</sub><sup>vm</sup> = 200 N mm<sup>-2</sup></td></tr><tr><td>Saturated yield stress</td><td>σ<sub>∞</sub> = 211.64 N mm<sup>-2</sup></td></tr><tr><td>Exponential hardening parameter</td><td>h<sub>e</sub> = 100</td></tr><tr><td>Linear hardening parameter</td><td>h = 209</td></tr><tr><td>Perzyna exponent</td><td>m = ∞</td></tr><tr><td>Perzyna hardening</td><td>Y<sub>0</sub><sup>*</sup> = 0</td></tr></tbody></table>

saturated isotropic hardening law, $\sigma^{vm} = \sigma_0^{vm} + [\sigma_\infty - \sigma_0^{vm}] [1 - \exp^{-h_e \epsilon^P}] + h_e^P$, which corresponds to the plastic dissipation potential

$$ \Phi^{pl} = \sigma_0^{vm} \epsilon^P + [\sigma_\infty - \sigma_0^{vm}] \left[ \epsilon^P + \frac{\exp^{-h_e \epsilon^P}}{h_e} \right] + \frac{h}{2} \epsilon^{P2}. \quad (110) $$

The mass of the tube can be increased by adding point masses at each extremity in order to represent the inertia

of the whole structure (vehicle, e.g.). Two cases are considered, without this addition of mass and with an addition of 0.02 kg at each extremity of the complete tube. The tube and the masses have an initial velocity of 50 m s⁻¹. The EDMC algorithm developed is used with the parameter $\rho_\infty = 0.8$ while contact interactions are treated using the method proposed by Armero and Petőcz [12], see also [31] for numerical implementation. Time step size is automatically computed using the method proposed in [32,33].

Snapshots of these simulations at time $t = 1$ ms can be found in Fig. 9. When additional mass is considered, the initial kinetic energy is multiplied by two and the resulting deformations are more important. In particular the two sides of the tube collapse and enter into contact together, contrarily to the case without additional mass.

Let us now detail the simulation without additional mass. At time $t = 0.5$ ms the tube is completely crushed against the punch, and the remaining kinetic energy $K$ is close to zero, see Fig. 10a. Between time $t = 0.5$ ms and time $t = 1$ ms the tube rebounds. This rebound occurs since

<figure><img src="image_10.png" /><figcaption>Fig. 9. Snapshots of deformed geometry and equivalent plastic strain for the side impact of square tube simulation at time $t = 1$ ms: (a) when no additional masses are considered and (b) with an additional mass of 0.02 kg added on each extremity (of the full square tube).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Time evolution of the energies in play when no additional masses are considered.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 11. Time evolution of the energies in play when additional masses are considered.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Snapshots and equivalent plastic strain for the side impact of square tube simulation at time $t = 2$ ms with an additional mass of 0.02 kg added on each extremity (of the full square tube): (a) 1/4 of the tube is represented and (b) full view of the tube is represented.</figcaption></figure>

the reversible stored energy $U_{el}$ can be transformed into kinetic energy $K$. This transfer can be observed in Fig. 10a. It is also obvious that the numerical oscillations are annihilated by the numerical dissipation, which allows the simulation to simulate the rebound. The contact energy, which is defined by the work of the contact forces, remains negative during the simulation (Fig. 10b). This ensures the stability of the time integration since no additional energy is introduced in the system. Since this contact energy remains negligible compared to the initial kinetic energy (less than 1%), accuracy is guaranteed.

Regarding the case where additional masses are considered, the time evolutions of these energies are represented in Fig. 11. It can be seen that after 1 ms the remaining kinetic energy $K$ of the system is equal to about 15% of the initial kinetic energy. The tube is wrapped around the punch in such a way that the central part has already rebound while the extremity is still moving forward. Contrarily to the previous case, the tube continues deforming

after 1 ms. Snapshots of this simulation at time $t = 2$ ms, when the tube has rebounded are reported in Fig. 12.

## 6. Conclusions

In this paper a new Energy-Dissipative Momentum-Conserving time integration algorithm is proposed for non-linear elasto-plastic dynamics.

This algorithm is based on an incremental variational updates formulation of the elasto-plastic behavior. Owing to this formulation, an effective stress potential can always be defined, even when irreversible deformations occur. The main originality of this paper is to add a dissipative function to this incremental potential, which leads to a new definition of a potential for the system.

From this potential, the internal forces expression developed by Gonzalez [11] to integrate in time the linear momentum equations, in such a way that the energy of the system is preserved, can directly be applied. The
introduction of the dissipative function in this potential does not only lead to a conserving algorithm for elasto-plastic behavior, but also introduce numerical dissipation, leading to an Energy-Dissipative Momentum-Conserving algorithm. As proposed by Armero and Romero [26,27], a dissipation velocity is also introduced to avoid bifurcation in the spectral properties. The resulting algorithm is first-order accurate, verifies the conservation of linear and angular momentum and ensures that the numerical dissipation is always positive.

All these properties are verified on numerical examples. The dynamics of a tumbling block is studied, and it is shown that the angular momenta are conserved, contrarily to Newmark-based time integration algorithms, leading to a better accuracy on the evaluation of the plastic deformation. An impact test also demonstrates the accuracy and robustness of the algorithm when large plastic deformations occur.

## Appendix A. Stiffness matrix

Let the consistent tangent stiffness matrix K associated to the internal forces (76) be decomposed in the following way

$$ \mathbf{K}^{\tilde{\xi}\mu} = \frac{\partial \vec{J}^{\tilde{\xi}}_{\text{int}}}{\partial \vec{x}^{\mu}} = \mathbf{K}_{\text{geo}}^{\tilde{\xi}\mu} + \mathbf{K}_{\text{vol}}^{\tilde{\xi}\mu} + \mathbf{K}_{\text{dev}}^{\tilde{\xi}\mu}, \quad (A.1) $$

with respectively

$$ \mathbf{K}_{\text{geo}}^{\tilde{\xi}\mu} = \int_{\gamma_0} \left\{ \frac{1}{2} \frac{\partial \mathbf{F}^{n+1}}{\partial \vec{x}^{\mu}} [\mathbf{S}_{\text{dev}}^* + 2p^* \mathbf{dG}] \vec{D}^{\tilde{\xi}} \right\} d\gamma_0, \quad (A.2) $$

$$ \mathbf{K}_{\text{dev}}^{\tilde{\xi}\mu} = \int_{\gamma_0} \left\{ \frac{\mathbf{F}^{n+1} + \mathbf{F}^n}{2} \frac{\partial \mathbf{S}_{\text{dev}}^*}{\partial \vec{x}^{\mu}} \vec{D}^{\tilde{\xi}} \right\} d\gamma_0, \quad (A.3) $$

$$ \mathbf{K}_{\text{vol}}^{\tilde{\xi}\mu} = \int_{\gamma_0} \left\{ [\mathbf{F}^{n+1} + \mathbf{F}^n] \frac{\partial p^*}{\partial \vec{x}^{\mu}} \mathbf{dG} \vec{D}^{\tilde{\xi}} \right\} d\gamma_0, \quad (A.4) $$

the geometrical part, the deviatoric part and the volume part. In these expressions the derivation is with respect to $\vec{x}^{n+1}$ but the superscript $n+1$ has been omitted for clarity purpose.

The following straightforward results will be used to establish the stiffness component:

$$ \frac{\partial \mathbf{F}_{ij}}{\partial \vec{x}_k^{\mu}} = \vec{D}_j^{\mu} \delta_{ik}, \quad \frac{\partial \mathbf{C}_{ij}}{\partial \vec{x}_k^{\mu}} = [\delta_{li} \mathbf{F}_{kj} + \mathbf{F}_{ki} \delta_{lj}] \vec{D}_l^{\mu} \frac{\partial \|\mathbf{\Delta C}\|^2}{\partial \mathbf{C}} = 2\mathbf{\Delta C} \\ \text{and} \quad \frac{\partial \sqrt{\det(\mathbf{C})}}{\partial \mathbf{C}} = \frac{1}{2} \sqrt{\det(\mathbf{C})} \mathbf{C}^{-T}. \quad (A.5) $$

### A.1. Geometrical part

Using relations (A.5), the geometrical part (A.2) can be computed as

$$ \mathbf{K}_{\text{geo}}^{\tilde{\xi}\mu} = \int_{\gamma_0} \left\{ \frac{\mathbf{I}}{2} \vec{D}^{\mu} \cdot [\mathbf{S}_{\text{dev}}^* + 2p^* \mathbf{dG}] \cdot \vec{D}^{\tilde{\xi}} \right\} d\mathcal{V}_0, \quad (A.6) $$

with **I** the unity tensor.

### A.2. Volume part

The volume part is decomposed into two terms: the first one $[\mathbf{K}_{\text{vol}}]^{1}$ results from the differentiation of the constant pressure (over the element), and the second one $[\mathbf{K}_{\text{vol}}]^{2}$ results from the differentiation of $\mathbf{dG}$.

Using (59), the first contribution states

$$ \begin{aligned} [\mathbf{K}_{\text{vol}}^{\tilde{\xi}\mu}]_{ik}^{1} &= \int_{\gamma_0} \left\{ \left[ \frac{\mathbf{F}^{n+1} + \mathbf{F}^n}{2} \mathbf{dG} \right]_{ij} \vec{D}_j^{\tilde{\xi}} \frac{\partial 2p^*}{\partial \theta^e} \frac{\partial \theta^e}{\partial \vec{x}_k^{\mu}} \right\} d\gamma_0 \\ &= \int_{\gamma_0} \left\{ \vec{D}_j^{\tilde{\xi}} \left[ \frac{\mathbf{F}^{n+1} + \mathbf{F}^n}{2} \mathbf{dG} \right]_{ij} \right. \\ & \qquad \left. \times \delta_{kl} \frac{\partial 2p^*}{\partial \theta^e} \frac{1}{\gamma_0^{\tilde{\xi}e}} \int_{\gamma_0} \{ \mathbf{J} \mathbf{F}^{n+1-T} \vec{D}^{\mu} \} \mathbf{l} d\gamma_0 \right\} d\gamma_0, \end{aligned} \quad (A.7) $$

with, using Eq. (80),

$$ \frac{\partial 2p^*}{\partial \theta^e} = \begin{cases} \frac{\partial^2 [\phi_{\text{vol}}^{\text{el}} + \Delta t \psi^{\text{d}}_{\text{vol}}]}{\partial \theta^e \partial \theta^e} \left( \frac{\theta^{e n+1} + \theta^{e n}}{2} \right) & \text{if } \Delta \theta^e \to 0 \\ \frac{2}{\Delta \theta^e} \frac{\partial [\phi_{\text{vol}}^{\text{el}} + \Delta t \psi^{\text{d}}_{\text{vol}}]}{\partial \theta^e} (\theta^{e n+1}) - \frac{2}{\Delta \theta^e} \frac{\phi_{\text{vol}}^{\text{el}}(\theta^{e n+1}) - \phi_{\text{vol}}^{\text{el}}(\theta^{e n}) - \Delta t \psi^{\text{d}}_{\text{vol}}(\Delta \theta^e)}{\Delta \theta^e} & \text{if } \Delta \theta^e \neq 0. \end{cases} \quad (A.8) $$

Using Eq. (A.5), the second contribution is

$$ \begin{aligned} [\mathbf{K}_{\text{vol}}^{\tilde{\xi}\mu}]_{ik}^{2} &= \int_{\gamma_0} \left\{ \frac{\mathbf{F}_{ir}^{n+1} + \mathbf{F}_{ir}^n}{2} 2p^* \frac{\partial \mathbf{dG}_{rf}}{\partial \vec{x}_k^{\mu}} \vec{D}_j^{\tilde{\xi}} \right\} d\gamma_0 \\ &= \int_{\gamma_0} \left\{ \vec{D}_j^{\tilde{\xi}} [\mathcal{H}_{\text{vol}}]_{ijkl} \vec{D}_l^{\mu} \right\} d\gamma_0, \end{aligned} \quad (A.9) $$

with the fourth order tensor term

$$ [\mathcal{H}_{\text{vol}}]_{ijkl} = 4p^* \frac{\mathbf{F}_{ir}^{n+1} + \mathbf{F}_{ir}^n}{2} \frac{\partial \mathbf{dG}_{rf}}{\partial \mathbf{C}_{ml}} \mathbf{F}_{km}^{n+1}. \quad (A.10) $$

The missing term can be computed by using Eqs. (A.5) and (77). Indeed, one has

$$ \begin{aligned} \frac{\partial \mathbf{dG}}{\partial \mathbf{C}} = \frac{1}{2} \left[ \mathcal{I} - \frac{\Delta \mathbf{C} \otimes \Delta \mathbf{C}}{\|\Delta \mathbf{C}\|^2} \right] : \frac{\partial^2 J^{n+\frac{1}{2}}}{\partial \mathbf{C}^2} \\ + \frac{J^{n+1}}{2} \frac{\Delta \mathbf{C} \otimes \mathbf{C}^{n+1-1}}{\|\Delta \mathbf{C}\|^2} - \frac{\Delta \mathbf{C} \otimes \frac{\partial \vec{J}^{n+\frac{1}{2}}}{\partial \mathbf{C}}}{\|\Delta \mathbf{C}\|^2} \\ + \left[ \frac{\Delta J - \frac{\partial J^{n+\frac{1}{2}}}{\partial \mathbf{C}} : \Delta \mathbf{C}}{\|\Delta \mathbf{C}\|^2} \right] \left[ \mathcal{I} - 2 \frac{\Delta \mathbf{C} \otimes \Delta \mathbf{C}}{\|\Delta \mathbf{C}\|^2} \right], \end{aligned} \quad (A.11) $$

with

$$ \frac{\partial^2 J^{n+\frac{1}{2}}}{\partial \mathbf{C}^2} = \frac{1}{2} \sqrt{\det(\mathbf{C}^{n+\frac{1}{2}})} \left[ \frac{1}{2} \mathbf{C}^{n+\frac{1}{2}-1} \otimes \mathbf{C}^{n+\frac{1}{2}-1} - \mathcal{I} \mathbf{C}^{n+\frac{1}{2}-1} \right]. \quad (A.12) $$

In this last expression $\mathbf{C}^{n+\frac{1}{2}} = \frac{\mathbf{C}^{n+1} + \mathbf{C}^n}{2}$ and $[\mathcal{I}\mathbf{A}]_{ijkl} = \frac{1}{2}\mathbf{A}_{ik}\mathbf{A}_{jl} + \frac{1}{2}\mathbf{A}_{il}\mathbf{A}_{jk}$.
### A.3. Deviatoric part

Owing to Eq. (A.5), the deviatoric term becomes

$$ egin{aligned} [\mathbf{K}_{\text{dev}}^{\xi\mu}]_{ik}^{1} &= \int_{\gamma_0} \left\{ \frac{\mathbf{F}_{ir}^{n+1} + \mathbf{F}_{ir}^{n}}{2} \frac{\partial \mathbf{S}_{\text{dev}rj}^{*}}{\partial \vec{x}_{k}^{\mu}} \vec{D}_{j}^{\xi} \right\} d\gamma_0 \\ &= \int_{\gamma_0} \left\{ \vec{D}_{j}^{\xi} [\mathcal{H}_{\text{dev}}]_{ijkl} \vec{D}_{l}^{\mu} \right\} d\gamma_0, \end{aligned} \quad (A.13) $$

with

$$ [\mathcal{H}_{\text{dev}}]_{ijkl} = \frac{\mathbf{F}_{ir}^{n+1} + \mathbf{F}_{ir}^{n}}{2} 2 \frac{\partial \mathbf{S}_{\text{dev}rj}^{*}}{\mathbf{C}_{ml}} \mathbf{F}_{km}^{n+1}. \quad (A.14) $$

Using relations (83), the derivation of the consistent deviatoric stress can be evaluated as

$$ egin{aligned} 2 \frac{\partial \mathbf{S}_{\text{dev}}^{*}}{\partial \mathbf{C}} ={}& \frac{1}{2} \left[ \mathcal{I} - \frac{\Delta \mathbf{C} \otimes \Delta \mathbf{C}}{\|\Delta \mathbf{C}\|^2} \right] : 4 \frac{d}{d\mathbf{C}^{n+\frac{1}{2}}} \left[ \frac{\partial \Delta \mathbf{D}^{\text{eff}}_{\text{dev}}}{\partial \mathbf{C}} \right]^{n+\frac{1}{2}} \\ & + \frac{\Delta \mathbf{C}}{\|\Delta \mathbf{C}\|^2} \otimes 4 \frac{\partial \Delta \mathbf{D}^{\text{eff}}_{\text{dev}}(\hat{\mathbf{C}}^{n+1}, \hat{\mathbf{C}}^n)}{\mathbf{d}\mathbf{C}^{n+1}} \\ & - \frac{\Delta \mathbf{C}}{\|\Delta \mathbf{C}\|^2} \otimes 4 \left[ \frac{\partial \Delta \mathbf{D}^{\text{eff}}_{\text{dev}}}{\partial \mathbf{C}} \right]^{n+\frac{1}{2}} \\ & + 4 \left[ \frac{\Delta \mathbf{D}^{\text{eff}}_{\text{dev}}(\hat{\mathbf{C}}^{n+1}, \hat{\mathbf{C}}^n) - \left[ \frac{\partial \Delta \mathbf{D}^{\text{eff}}_{\text{dev}}}{\partial \mathbf{C}} \right]^{n+\frac{1}{2}}}{\|\Delta \mathbf{C}\|^2} : \Delta \mathbf{C} \right] \\ & \times \left[ \mathcal{I} - 2 \frac{\Delta \mathbf{C} \otimes \Delta \mathbf{C}}{\|\Delta \mathbf{C}\|^2} \right]. \end{aligned} \quad (A.15) $$

In this expression, differentiation with symbol *d* is used instead of $\partial$ because the minimum value of $\Delta D$ depends on C. Moreover, exponent $n + \frac{1}{2}$ refers to values computed for $\frac{C^{n+1}+C^n}{2}$. Material tensor $M = 4 \frac{d}{dC^{n+\frac{1}{2}}} \left[ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}}}{\partial C} \right]^{n+\frac{1}{2}}$ can be decomposed into

$$ M = \det(\mathbf{C}^{n+\frac{1}{2}})^{-\frac{2}{3}} \left[ \begin{aligned} & \mathbf{C}^{n+\frac{1}{2}} + \frac{4}{3} \left[ \mathcal{I} \hat{\mathbf{C}}^{-1} - \frac{1}{3} \hat{\mathbf{C}}^{-1} \otimes \hat{\mathbf{C}}^{-1} \right]^{n+\frac{1}{2}} \\ & \times \left[ \mathbf{F}^{\text{pl}-1} \frac{\partial[\mathbf{\Phi}_{\text{dev}}^{\text{el}} + \Delta t \mathbf{\Psi}_{\text{dev}}^{\text{d}}]}{\partial \hat{\mathbf{C}}^{\text{el}}} \mathbf{F}^{\text{pl}-T} : \hat{\mathbf{C}} \right]^{n+\frac{1}{2}} \\ & - \frac{4}{3} \left[ \hat{\mathbf{C}}^{-1} \otimes \left[ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}}}{\partial \mathbf{C}} \right]^{n+\frac{1}{2}} + \left[ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}}}{\partial \mathbf{C}} \right]^{n+\frac{1}{2}} \otimes \hat{\mathbf{C}}^{-1} \right]^{n+\frac{1}{2}} \end{aligned} \right], \quad (A.16) $$

with

$$ \begin{aligned} \mathcal{C} ={}& \frac{d\partial 4\Delta D_{\text{dev}}^{\text{eff}}}{d\hat{\mathbf{C}}\partial\hat{\mathbf{C}}} - \frac{1}{3}\hat{\mathbf{C}}^{-1} \otimes \left[ \hat{\mathbf{C}} : \frac{d\partial 4\Delta D_{\text{dev}}^{\text{eff}}}{d\hat{\mathbf{C}}\partial\hat{\mathbf{C}}} \right] \\ & - \frac{1}{3} \left[ \frac{d\partial 4\Delta D_{\text{dev}}^{\text{eff}}}{d\hat{\mathbf{C}}\partial\hat{\mathbf{C}}} : \hat{\mathbf{C}} \right] \otimes \hat{\mathbf{C}}^{-1} \\ & + \frac{1}{9} \left[ \hat{\mathbf{C}} : \frac{d\partial 4\Delta D_{\text{dev}}^{\text{eff}}}{d\hat{\mathbf{C}}\partial\hat{\mathbf{C}}} : \hat{\mathbf{C}} \right] \hat{\mathbf{C}}^{-1} \otimes \hat{\mathbf{C}}^{-1}. \end{aligned} \quad (A.17) $$

$$ \text{The first total derivative } 4 \frac{d\Delta D_{\text{dev}}^{\text{eff}}(\hat{\mathbf{C}}^{n+1}, \hat{\mathbf{C}}^n)}{d\mathbf{C}^{n+1}} \text{ is given by } \frac{d\Delta D_{\text{dev}}^{\text{eff}}(\hat{\mathbf{C}}^{n+1}, \hat{\mathbf{C}}^n)}{d\mathbf{C}^{n+1}} = J^{n+1-\frac{2}{3}} \text{DEV} \frac{d\Delta D_{\text{dev}}^{\text{eff}}(\hat{\mathbf{C}}^{n+1}, \hat{\mathbf{C}}^n)}{d\hat{\mathbf{C}}}, \quad (A.18) $$

with

$$ \begin{aligned} \frac{d\Delta D_{\text{dev}}^{\text{eff}}}{d\hat{\mathbf{C}}} ={}& \frac{\partial \Delta D_{\text{dev}}}{\partial \hat{\mathbf{C}}} + \frac{\partial \Delta D_{\text{dev}}}{\partial \varepsilon^P} \frac{\partial \varepsilon^P}{\partial \hat{\mathbf{C}}} \\ & + \frac{\partial \Delta D_{\text{dev}}}{\partial \mathbf{N}} : \frac{\partial \mathbf{N}}{\partial \hat{\mathbf{C}}} + \frac{\partial \Delta t \mathbf{\Psi}_{\text{dev}}^{\text{d}}}{\partial \hat{\mathbf{C}}^{\text{el}}} : \frac{d\hat{\mathbf{C}}^{\text{el}}}{d\hat{\mathbf{C}}}. \end{aligned} \quad (A.19) $$

Some terms of this relation depend on the choice of the potentials and will be given in the logarithm case. Nevertheless, in the general case, one has

$$ \frac{d\hat{\mathbf{C}}^{\text{el}}}{d\hat{\mathbf{C}}} = \frac{d\mathbf{F}^{\text{pl}-T}}{d\hat{\mathbf{C}}} \hat{\mathbf{C}} \mathbf{F}^{\text{pl}-1} + \mathcal{I}_{\mathbf{F}^{\text{pl}-T}} + \mathbf{F}^{\text{pl}-T} \frac{d\mathbf{F}^{\text{pl}-1}}{d\hat{\mathbf{C}}}, \quad \text{with} \quad (A.20) $$

$$ \begin{aligned} \frac{d\mathbf{F}_{ij}^{\text{pl}-1}}{d\hat{\mathbf{C}}_{kl}} ={}& \frac{\partial \mathbf{F}_{ij}^{\text{pl}-1}}{\partial \varepsilon^P} \frac{\partial \varepsilon^P}{\partial \hat{\mathbf{C}}_{kl}} + \frac{\partial \mathbf{F}_{ij}^{\text{pl}-1}}{\partial \mathbf{N}_{mn}} \frac{\partial \mathbf{N}_{mn}}{\partial \hat{\mathbf{C}}_{kl}} \\ ={}& - \mathbf{F}_{im}^{\text{pl}-1} \mathbf{N}_{mj} \frac{\partial \varepsilon^P}{\partial \hat{\mathbf{C}}_{kl}} - \mathbf{F}_{ip}^{\text{pl}-1} \mathcal{G}_{pqmn}^{\text{exp}} [\mathbf{F}_{ip}^{\text{pl}n} \mathbf{F}_{ip}^{\text{pl}-1}]_{qj} \frac{\partial \mathbf{N}_{mn}}{\partial \hat{\mathbf{C}}_{kl}}, \end{aligned} \quad (A.21) $$

where Eq. (47) has been used, and where $\mathcal{G}^{\text{exp}} = \partial \exp(\Delta \varepsilon^P \mathbf{N}) / \partial \mathbf{N}$ is computed using the spectral decomposition. Therefore Eq. (A.19) is rewritten

$$ \frac{d\Delta D_{\text{dev}}^{\text{eff}}}{d\hat{\mathbf{C}}} = \mathbf{F}^{\text{pl}-1} \frac{\partial \mathbf{\Phi}_{\text{dev}}^{\text{el}}}{\partial \hat{\mathbf{C}}^{\text{el}}} \mathbf{F}^{\text{pl}-T} \\ + \frac{\partial \Delta D_{\text{dev}}}{\partial \mathbf{N}} : \frac{\partial \mathbf{N}}{\partial \hat{\mathbf{C}}} + \frac{\partial \Delta t \mathbf{\Psi}_{\text{dev}}^{\text{d}}}{\partial \hat{\mathbf{C}}^{\text{el}}} : \frac{d\hat{\mathbf{C}}^{\text{el}}}{d\hat{\mathbf{C}}}, \quad (A.22) $$

$$ \frac{\partial \Delta D_{\text{dev}}}{\partial \mathbf{N}} = - \left[ \mathbf{F}^{\text{pl}-1} \frac{\partial \mathbf{\Phi}_{\text{dev}}^{\text{el}}}{\partial \hat{\mathbf{C}}^{\text{el}}} \mathbf{F}^{\text{pl}-T} \right] : \left[ \mathbf{F}^{\text{pl}T} \hat{\mathbf{C}}^{\text{el}} \mathcal{G}^{\text{exp}} \mathbf{F}^{\text{pl}nT} \right] \\ - \left[ \mathbf{F}^{\text{pl}-1} \frac{\partial \mathbf{\Phi}_{\text{dev}}^{\text{el}}}{\partial \hat{\mathbf{C}}^{\text{el}}} \mathbf{F}^{\text{pl}-T} \right] : \left[ \mathbf{F}^{\text{pl}T} \hat{\mathbf{C}}^{\text{el}} \mathcal{G}^{\text{exp}} \mathbf{F}^{\text{pl}nT} \right]^{\text{T}}, \quad (A.23) $$

since $\frac{\partial \Delta D_{\text{dev}}}{\partial \varepsilon^P} = 0$ (due to the minimization process, but $\frac{\partial \Delta D_{\text{dev}}}{\partial \mathbf{N}} \neq 0$, since there are constraints on N).

The second total derivative $\frac{d\Delta 4\Delta D_{\text{dev}}^{\text{eff}}}{d\hat{\mathbf{C}}\partial\hat{\mathbf{C}}}$ is obtained by using the same method, leading to

$$ \begin{aligned} \frac{d\partial \Delta D_{\text{dev}}^{\text{eff}}}{d\hat{\mathbf{C}}\partial\hat{\mathbf{C}}} ={}& \frac{\partial^2 \mathbf{\Phi}_{\text{dev}}^{\text{el}}}{\partial\hat{\mathbf{C}}\partial\hat{\mathbf{C}}} + \frac{1}{2} \frac{\partial^2 \mathbf{\Phi}_{\text{dev}}^{\text{el}}}{\partial\varepsilon^P\partial\hat{\mathbf{C}}} \otimes \frac{\partial\varepsilon^P}{\partial\hat{\mathbf{C}}} + \frac{1}{2} \frac{\partial\varepsilon^P}{\partial\hat{\mathbf{C}}} \otimes \frac{\partial^2 \mathbf{\Phi}_{\text{dev}}^{\text{el}}}{\partial\varepsilon^P\partial\hat{\mathbf{C}}} \\ & + \frac{1}{2} \frac{\partial^2 \mathbf{\Phi}_{\text{dev}}^{\text{el}}}{\partial\mathbf{N}\partial\hat{\mathbf{C}}} : \frac{\partial\mathbf{N}}{\partial\hat{\mathbf{C}}} + \frac{1}{2} \left[ \frac{\partial\mathbf{N}}{\partial\hat{\mathbf{C}}} \right]^{\text{T}} : \left[ \frac{\partial^2 \mathbf{\Phi}_{\text{dev}}^{\text{el}}}{\partial\mathbf{N}\partial\hat{\mathbf{C}}} \right]^{\text{T}} \\ & + \mathbf{F}^{\text{pl}-1} \frac{\partial^2 \mathbf{\Psi}_{\text{dev}}^{\text{d}}}{\partial\hat{\mathbf{C}}^{\text{el}}\partial\hat{\mathbf{C}}^{\text{el}}} \mathbf{F}^{\text{pl}-T} : \frac{\partial\hat{\mathbf{C}}^{\text{el}}}{\partial\hat{\mathbf{C}}} \\ & + \Delta t \frac{\partial\mathbf{F}^{\text{pl}-1}}{\partial\hat{\mathbf{C}}} \frac{\partial\mathbf{\Psi}_{\text{dev}}^{\text{d}}}{\partial\hat{\mathbf{C}}^{\text{el}}} \mathbf{F}^{\text{pl}-T} + \Delta t \mathbf{F}^{\text{pl}-1} \frac{\partial\mathbf{\Psi}_{\text{dev}}^{\text{d}}}{\partial\hat{\mathbf{C}}^{\text{el}}} \frac{\partial\mathbf{F}^{\text{pl}-T}}{\partial\hat{\mathbf{C}}}, \end{aligned} \quad (A.24) $$
where $\mathcal{H}_{ijkl}^{TT} = \mathcal{H}_{klij}$, and where Eqs. (A.20) and (A.21) define the derivation of elastic and plastic tensors. The explicit derivations of the elastic potential can be evaluated using (47), as being

$$ \frac{\partial^2 \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C} \partial \hat{C}} = \mathbf{F}^{\text{pl}-1} \mathbf{F}^{\text{pl}-1} \frac{\partial^2 \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C}^{\text{el}} \partial \hat{C}^{\text{el}}} \mathbf{F}^{\text{pl}-1} \mathbf{F}^{\text{pl}-1}, \quad (\text{A.25}) $$

$$ \frac{\partial^2 \Phi_{\text{dev}}^{\text{el}}}{\partial \varepsilon^p \partial \hat{C}} = -\mathbf{F}^{\text{pl}-1} \mathbf{N} \frac{\partial \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C}^{\text{el}}} \mathbf{F}^{\text{pl}-T} - \mathbf{F}^{\text{pl}-1} \frac{\partial \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C}^{\text{el}}} \mathbf{N} \mathbf{F}^{\text{pl}-T} \\ - \frac{\partial^2 \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C} \partial \hat{C}} : [\mathbf{F}^{\text{pl}T} \hat{C}^{\text{el}} \mathbf{N} \mathbf{F}^{\text{pl}} + \mathbf{F}^{\text{pl}T} \mathbf{N} \hat{C}^{\text{el}} \mathbf{F}^{\text{pl}}], \quad (\text{A.26}) $$

$$ \frac{\partial^2 \Phi_{\text{dev}}^{\text{el}}}{\partial \mathbf{N} \partial \hat{C}} = -\mathbf{F}^{\text{pl}-1} \mathcal{D}^{\text{exp}} \mathbf{F}^{\text{pl}T} \mathbf{F}^{\text{pl}-1} \frac{\partial \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C}^{\text{el}}} \mathbf{F}^{\text{pl}-T} \\ - \mathbf{F}^{\text{pl}-1} \frac{\partial \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C}^{\text{el}}} \mathbf{F}^{\text{pl}-T} [\mathbf{F}^{\text{pl}-1} \mathcal{D}^{\text{exp}} \mathbf{F}^{\text{pl}T}]^T \\ - \frac{\partial^2 \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C} \partial \hat{C}} : [\mathbf{F}^{\text{pl}T} \hat{C}^{\text{el}} \mathcal{D}^{\text{exp}} \mathbf{F}^{\text{pl}T} + [\mathbf{F}^{\text{pl}T} \hat{C}^{\text{el}} \mathcal{D}^{\text{exp}} \mathbf{F}^{\text{pl}T}]^T]. \quad (\text{A.27}) $$

At this point, the general form of the stiffness matrix is defined. The missing terms ($\frac{\partial \Phi_{\text{vol}}^{\text{el}}}{\partial \theta^e}$, $\frac{\partial^2 \Phi_{\text{vol}}^{\text{el}}}{\partial \theta^e \partial \theta^e}$, $\Delta t \frac{\partial \Psi^{\text{d}}}{\partial \theta^e \partial \theta^e}$, $\Delta t \frac{\partial^2 \Psi^{\text{d}}}{\partial \theta^e \partial \theta^e}$ for the volume part and $\frac{\partial \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C}^{\text{el}}}$, $\frac{\partial^2 \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C}^{\text{el}} \partial \hat{C}^{\text{el}}}$, $\Delta t \frac{\partial \Psi^{\text{d}}}{\partial \hat{C}^{\text{el}}}$, $\Delta t \frac{\partial^2 \Psi^{\text{d}}}{\partial \hat{C}^{\text{el}} \partial \hat{C}^{\text{el}}}$, $\frac{\partial \varepsilon^p}{\partial \hat{C}}$ and $\frac{\partial \mathbf{N}}{\partial \hat{C}}$ for the deviatoric part) depend on the choice of the potential.

### A.4. Bi-logarithmic potential

In this particular case, the derivation of the elastic potential are directly deduced from Eq. (68), with

$$ \frac{\partial \Phi_{\text{vol}}^{\text{el}}}{\partial \theta^e} = K_0 \frac{\ln(\theta^e)}{\theta^e}, \quad (\text{A.28}) $$

$$ \frac{\partial^2 \Phi_{\text{vol}}^{\text{el}}}{\partial \theta^{e2}} = K_0 \frac{1 - \ln(\theta^e)}{\theta^{e2}} \quad (\text{A.29}) $$

and $\frac{\partial \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C}^{\text{el}}}$, $\frac{\partial^2 \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C}^{\text{el}} \partial \hat{C}^{\text{el}}}$ are computed using the spectral decomposition.

Terms resulting from the numerical dissipation are computed from (70), leading to

$$ \Delta t \frac{\partial \Psi_{\text{vol}}^{\text{d}}}{\partial \theta^e} = \chi \frac{\partial^2 \Phi_{\text{vol}}^{\text{el}}(\theta^{e'n})}{\partial \theta^{e2}} \Delta \theta^e, \quad (\text{A.30}) $$

$$ \Delta t \frac{\partial^2 \Psi_{\text{vol}}^{\text{d}}}{\partial \theta^{e2}} = \chi \frac{\partial^2 \Phi_{\text{vol}}^{\text{el}}(\theta^{e'n})}{\partial \theta^{e2}}, \quad (\text{A.31}) $$

$$ \Delta t \frac{\partial \Psi_{\text{dev}}^{\text{d}}}{\partial \hat{C}^{\text{el}}} = \chi \frac{\partial^2 \Phi_{\text{dev}}^{\text{el}}(\hat{C}^{\text{el}''})}{\partial \hat{C}^{\text{el}} \partial \hat{C}^{\text{el}}} : \Delta \hat{C}^{\text{el}}, \quad \text{and to} \quad (\text{A.32}) $$

$$ \Delta t \frac{\partial^2 \Psi_{\text{dev}}^{\text{d}}}{\partial \hat{C}^{\text{el}} \partial \hat{C}^{\text{el}}} = \chi \frac{\partial^2 \Phi_{\text{dev}}^{\text{el}}(\hat{C}^{\text{el}''})}{\partial \hat{C}^{\text{el}} \partial \hat{C}^{\text{el}}}. \quad (\text{A.33}) $$

Using developments done in [19], the last remaining terms can be evaluated in the closed-form

$$ \frac{\partial \varepsilon^p}{\partial \hat{C}} = \frac{3}{2} \frac{G_0}{(3G_0 + h) \sqrt{\ln \hat{C}^{\text{pr}} : \mathcal{D}^{\ln, \text{pr}} : \mathcal{J}_{\mathbf{F}^{\text{pl}T} \partial \hat{C}^{\text{el}}}}} \ln \hat{C}^{\text{pr}} : \mathcal{D}^{\ln, \text{pr}} : \mathcal{J}_{\mathbf{F}^{\text{pl}T} \partial \hat{C}^{\text{el}}}, \quad (\text{A.34}) $$

$$ \frac{\partial \mathbf{N}}{\partial \hat{C}} = \sqrt{\frac{3}{2 \ln \hat{C}^{\text{pr}} : \ln \hat{C}^{\text{pr}}}} \mathcal{D}^{\ln, \text{pr}} : \mathcal{J}_{\mathbf{F}^{\text{pl}T} \partial \hat{C}^{\text{el}}} - \sqrt{\frac{3}{2 \ln \hat{C}^{\text{pr}} : \ln \hat{C}^{\text{pr}}}} \ln \hat{C}^{\text{pr}} \otimes \ln \hat{C}^{\text{pr}} : \mathcal{D}^{\ln, \text{pr}} : \mathcal{J}_{\mathbf{F}^{\text{pl}T} \partial \hat{C}^{\text{el}}}, \quad (\text{A.35}) $$

where $h$ is the hardening value, $\hat{C}^{\text{pr}}$ is the elastic predictor [19], and where $\mathcal{D}^{\ln, \text{pr}} = \frac{\partial \ln \hat{C}^{\text{pr}}}{\partial \hat{C}^{\text{pr}}}$ is computed using the spectral decomposition.

## Appendix B. Linearization

Considering the deformation gradient (100), a linearization leads to the following deformations,

$$ \begin{aligned} \mathbf{C} &= 2\mathbf{F} - \mathbf{I} + \mathcal{O}\left(\frac{x^2}{l^2}\right) \\ &= \begin{pmatrix} 1 + 2\frac{x}{l} & 0 & 0 \\ 0 & 1 - 2v\frac{x}{l} & 0 \\ 0 & 0 & 1 - 2v\frac{x}{l} \end{pmatrix} + \mathcal{O}\left(\frac{x^2}{l^2}\right), \end{aligned} \quad (\text{B.1}) $$

$$ J = 1 + \frac{x}{l}(1 - 2v) + \mathcal{O}\left(\frac{x^2}{l^2}\right), \quad \text{and} \quad (\text{B.2}) $$

$$ \hat{\mathbf{C}} = \begin{pmatrix} 1 + \frac{4x}{3l}(1 + v) & 0 & 0 \\ 0 & 1 - \frac{2x}{3l}(1 + v) & 0 \\ 0 & 0 & 1 - \frac{2x}{3l}(1 + v) \end{pmatrix} + \mathcal{O}\left(\frac{x^2}{l^2}\right). \quad (\text{B.3}) $$

The modified pressure (80) of the internal forces is therefore rewritten

$$ \Delta \theta^e = (1 - 2v) \frac{x^{n+1} - x^n}{l} + \mathcal{O}\left(\frac{x^2}{l^2}\right), \quad (\text{B.4}) $$

$$ \Delta D_{\text{vol}}^{\text{eff}} = \frac{K_0 (1 - 2v)^2 \Delta x}{2l} \left( \frac{x^{n+1} + x^n}{l} + \chi \frac{\Delta x}{l} \right) + \mathcal{O}\left(\frac{x^3}{l^3}\right), \quad (\text{B.5}) $$

$$ p^* = \frac{K_0 (1 - 2v)}{2} \left( \frac{x^{n+1} + x^n}{l} + \chi \frac{\Delta x}{l} \right) + \mathcal{O}\left(\frac{x^2}{l^2}\right), \quad (\text{B.6}) $$

while the pressure base (77) is rewritten

$$ \Delta \mathbf{C} = \frac{2\Delta x}{l} \begin{pmatrix} 1 & 0 & 0 \\ 0 & -v & 0 \\ 0 & 0 & -v \end{pmatrix} + \mathcal{O}\left(\frac{x^2}{l^2}\right), \quad (\text{B.7}) $$

$$ \frac{\partial J^{n+\frac{1}{2}}}{\partial \mathbf{C}} = \frac{1}{2} \begin{pmatrix} 1 - \frac{x^{n+1}+x^n}{l} \left(\frac{1}{2} + v\right) & 0 & 0 \\ 0 & 1 + \frac{x^{n+1}+x^n}{2l} & 0 \\ 0 & 0 & 1 + \frac{x^{n+1}+x^n}{2l} \end{pmatrix} + \mathcal{O}\left(\frac{x^2}{l^2}\right), \quad (\text{B.8}) $$
$$ \mathbf{dG} = \frac{1}{2}\mathbf{I} + \mathcal{O}\left(\frac{\mathbf{x}}{l}\right), \qquad (B.9) $$

$$ \text{since } \frac{\Delta J - \frac{\partial J^{n+\frac{1}{2}}}{\partial C} \cdot \Delta C}{\Delta C \cdot \Delta C} \Delta \mathbf{C} = \mathcal{O}\left(\frac{\mathbf{x}}{l}\right). $$

Concerning the deviatoric part (83), it comes

$$ \Delta D_{\text{dev}}^{\text{eff}} = \frac{2G_0(1+v)^2}{3} \frac{x^{n+12} - x^n^2}{l^2} + \chi \frac{(x^{n+1} - x^n)^2}{l^2} + \mathcal{O}\left(\frac{x^3}{l^3}\right), \qquad (B.10) $$

$$ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}^{n+\frac{1}{2}}}}{\partial \hat{\mathbf{C}}} = \frac{G_0(1+v)}{3} \left( \frac{x^{n+1} + x^n}{l} + \chi \frac{x^{n+1} - x^n}{l} \right) \times \begin{pmatrix} 1 & 0 & 0 \\ 0 & -\frac{1}{2} & 0 \\ 0 & 0 & -\frac{1}{2} \end{pmatrix} + \mathcal{O}\left(\frac{x^2}{l^2}\right), \qquad (B.11) $$

$$ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}^{n+\frac{1}{2}}}}{\partial \mathbf{C}} = \frac{G_0(1+v)}{3} \left( \frac{x^{n+1} + x^n}{l} + \chi \frac{x^{n+1} - x^n}{l} \right) \times \begin{pmatrix} 1 & 0 & 0 \\ 0 & -\frac{1}{2} & 0 \\ 0 & 0 & -\frac{1}{2} \end{pmatrix} + \mathcal{O}\left(\frac{x^2}{l^2}\right), \qquad (B.12) $$

$$ \mathbf{S}_{\text{dev}}^{*} = \frac{G_0(1+v)}{3} \left( \frac{x^{n+1} + x^n}{l} + \chi \frac{x^{n+1} - x^n}{l} \right) \begin{pmatrix} 1 & 0 & 0 \\ 0 & -\frac{1}{2} & 0 \\ 0 & 0 & -\frac{1}{2} \end{pmatrix} + \mathcal{O}\left(\frac{x^2}{l^2}\right), \qquad (B.13) $$

$$ \text{since } \frac{\Delta D_{\text{dev}}^{\text{eff}} - \frac{\partial \Delta D_{\text{dev}}^{\text{eff}}}{\partial \mathbf{C}} \cdot \Delta \mathbf{C}}{\Delta \mathbf{C} \cdot \Delta \mathbf{C}} \Delta \mathbf{C} = \mathcal{O}\left(\frac{x^2}{l^2}\right). $$

Using the relations $G_0 = \frac{E}{2+2\nu}$ and $K_0 = \frac{E}{3-6\nu}$, and combining Eqs. (B.6), (B.9) and (B.13) leads to

$$ \frac{\mathbf{F}^{n+1} + \mathbf{F}^n}{2} (\mathbf{S}_{\text{dev}}^* + 2p^* \mathbf{dG}) = \left( \frac{x^{n+1} + x^n}{l} + \chi \frac{x^{n+1} - x^n}{l} \right) \begin{pmatrix} \frac{E}{2} \\ 0 \\ 0 \end{pmatrix} \qquad (B.14) $$

which allows to evaluate the internal forces (101).

## References

[1] N. Newmark, A method of computation for structural dynamics, J. Engrg. Mech. Div. ASCE 85 (EM3) (1959) 67–94.

[2] T. Belytschko, D. Schoeberle, On the unconditional stability of an implicit algorithm for non-linear structural dynamics, J. Appl. Mech. 42 (1975) 865–869.

[3] T. Hughes, A note on the stability of Newmark’s algorithm in nonlinear structural dynamics, Int. J. Numer. Methods Engrg. 11 (2) (1977) 383–386.

[4] P. Betsch, P. Steinmann, Conservation properties of a time fe method. part I: Time-stepping schemes for n-body problems, Int. J. Numer. Methods Engrg. 49 (2000) 599–638.

[5] P. Betsch, P. Steinmann, Conservation properties of a time FE method. part II: Time-stepping schemes for non-linear elastodynamics, Int. J. Numer. Methods Engrg. 50 (2001) 1931–1955.

[6] O. Bauchau, N. Theron, Energy decaying scheme for non-linear beam models, Comput. Methods Appl. Mech. Engrg. 134 (1996) 37–56.

[7] O. Bauchau, T. Joo, Computational schemes for non-linear elastodynamics, Int. J. Numer. Methods Engrg. 45 (1999) 693–719.

[8] C. Bottasso, M. Borri, Integrating finite rotation, Comput. Methods Appl. Mech. Engrg. 164 (1998) 307–331.

[9] C. Bottasso, M. Borri, Energy preserving/decaying schemes for non-linear beam dynamics using the helicoidal approximation, Comput. Methods Appl. Mech. Engrg. 143 (1997) 393–415.

[10] J. Simo, N. Tarnow, The discrete energy–momentum method. Conserving algorithms for nonlinear elastodynamics, ZAMP 43 (1992) 757–792.

[11] O. Gonzalez, Exact energy and momentum conserving algorithms for general models in nonlinear elasticity, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1763–1783.

[12] F. Armero, E. Petőcz, Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems, Comput. Methods Appl. Mech. Engrg. 158 (1998) 269–300.

[13] F. Armero, E. Petőcz, A new dissipative time-stepping algorithm for frictional contact problems: formulation and analysis, Comput. Methods Appl. Mech. Engrg. 179 (1999) 151–178.

[14] X. Meng, T. Laursen, Energy consistent algorithms for dynamic finite deformation plasticity, Comput. Methods Appl. Mech. Engrg. 191 (2001) 1639–1675.

[15] X. Meng, T. Laursen, On energy consistency of large deformation plasticity models, with application to the design of unconditionally stable time integrators, Finite Elem. Anal. Des. 38 (2002) 949–963.

[16] L. Noels, L. Stainier, J.-P. Ponthot, Energy–momentum conserving algorithm for non-linear hypoelastic constitutive models, Int. J. Numer. Methods Engrg. 59 (2004) 83–114.

[17] L. Noels, L. Stainier, J.-P. Ponthot, On the use of large time steps with an energy–momentum conserving algorithm for non-linear hypoelastic constitutive models, Int. J. Solids Struct. 41 (2004) 663–693.

[18] F. Armero, Energy-dissipative momentum-conserving time-stepping algorithms for finite strain multiplicative plasticity, Comput. Methods Appl. Mech. Engrg. 195 (2006) 4862–4889.

[19] L. Noels, L. Stainier, J.-P. Ponthot, An energy momentum conserving algorithm using the variational formulation of visco-plastic updates, Int. J. Numer. Methods Engrg. 65 (2006) 904–942.

[20] R. Radovitzky, M. Ortiz, Error estimation and adaptative meshing in strongly nonlinear dynamics problems, Comput. Methods Appl. Mech. Engrg. 172 (1999) 203–240.

[21] M. Ortiz, L. Stainier, The variational formulation of viscoplastic updates, Comput. Methods Appl. Mech. Engrg. 171 (1999) 419–444.

[22] E. Fancello, J.-P. Ponthot, L. Stainier, A variational formulation of constitutive models and updates in non-linear finite viscoelasticity, Int. J. Numer. Methods Engrg. 65 (2006) 1831–1864.

[23] Q. Yang, L. Stainier, M. Ortiz, A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids, J. Mech. Phys. Solids 54 (2006) 401–424.

[24] J. Chung, G. Hulbert, A time integration algorithms for structural dynamics with improved numerical dissipations: the generalized-α method, J. Appl. Mech. 60 (1993) 371–375.

[25] S. Erlicher, L. Bonaventura, O. Bursi, The analysis of the α-generalized method for non-linear dynamic problems, Comput. Mech. 28 (2002) 83–104.

[26] F. Armero, I. Romero, On the formulation of high-frequency dissipative time-stepping algorithms for non-linear dynamics. Part I: low-order methods for two model problems and nonlinear elastodynamics, Comput. Methods Appl. Mech. Engrg. 190 (2001) 2603–2649.

[27] F. Armero, I. Romero, On the formulation of high-frequency dissipative time-stepping algorithms for non-linear dynamics. Part II: second-order methods, Comput. Methods Appl. Mech. Engrg. 190 (2001) 6783–6824.

[28] L. Noels, L. Stainier, J.-P. Ponthot, Simulation of complex impact problems with implicit time algorithm, Application to a blade-loss problem, Int. J. Impact Engrg. 32 (2005) 358–386.
[29] J. Lubliner, Plasticity Theory, Macmillan, 1990.

[30] J. Simo, R. Taylor, Quasi-incompressible finite elasticity in principal stretches continuum basis and numerical algorithms, Comput. Methods Appl. Mech. Engrg. 85 (1991) 273–310.

[31] L. Noels, L. Stainier, J.-P. Ponthot, Simulation of crashworthiness problems with improved contact algorithms for implicit time integration, Int. J. Impact Engrg. 32 (2006) 799–825.

[32] L. Noels, L. Stainier, J.-P. Ponthot, J. Bonini, Automatic time stepping algorithms for implicit numerical simulations of non-linear dynamics, Adv. Engrg. Soft. 33 (10) (2002) 581–595.

[33] L. Noels, L. Stainier, J.-P. Ponthot, Self-adapting time integration management in crash-worthiness and sheet metal forming computations, Int. J. Vehicle Des. 30 (2) (2002) 67–114.