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

HAL
open science

# Variational approach for nonsmooth elasto-plastic dynamics with contact and impacts

Vincent Acary, Franck Bourrier, Benoit Viano

► To cite this version:

Vincent Acary, Franck Bourrier, Benoit Viano. Variational approach for nonsmooth elasto-plastic dynamics with contact and impacts. Computer Methods in Applied Mechanics and Engineering, 2023, 414, pp.116156. 10.1016/j.cma.2023.116156. hal-03978387v4

HAL Id: hal-03978387

https://inria.hal.science/hal-03978387v4

Submitted on 18 Sep 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Variational approach for nonsmooth elasto-plastic dynamics
with contact and impacts

Vincent Acary*, Franck Bourrier*,† Benoit Viano*

May 25, 2023

**Abstract** The objective of this article is the modelling and the numerical simulation of the response of elastoplastic structures to impacts. To this end, a numerical method is proposed that takes into account one-sided contact (Signorini condition) and impact phenomena together with plasticity in a monolithic solver, while accounting for the non-smooth character of the dynamics. The formulation of the plasticity and the contact laws are based on inclusions into normal cones of convex sets, or equivalently, variational inequalities following the pioneering work of Moreau (1974) and Halphen and Nguyen (1975), who introduced the assumptions of normal dissipation and of generalised standard materials (GSM) in the framework of associated plasticity with strain hardening. The proposed time-stepping method is an extension of the Jean and Moreau (1987) scheme for nonsmooth dynamics. The discrete energy balance shows that spurious numerical damping can be suppressed and that the scheme is in practice unconditionally stable. Furthermore, the finite-dimensional variational inequality at each time-step is well-posed and can be solved by optimization methods for convex quadratic programs, providing an interesting alternative to the return mapping algorithm coupled with a dedicated frictional contact method. The paper is completed by illustrative numerical examples of impacts on steel structures consisting of beams.

**Keywords:** contact, impact, plasticity, nonsmooth dynamics, variational inequality, normal cone, generalised standard materials, saddle point, time-stepping methods, numerical optimisation.

**Highlights.**

1. A formulation of all the dynamical equations of the elastoplastic law and contact law in a single variational inequality, respecting the principles of the thermodynamics

2. A time-stepping method, on the basis of the Moreau-Jean time stepping scheme is developed enabling the integration of the nonsmooth dynamics consistently with the impulsive motion.

3. Formulation of a convex quadratic program with well-posedness results (existence and uniqueness) that can be solved efficiently with large time/load increments. This ensures that the scheme is well-posed

4. A discrete energy balance is given that ensures that the scheme is practically stable.

5. Simulation of multi-criteria elasto-plastic flow rules with contact constraints. Application to impact on structures made of beams.

<footnote>*Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, Institute of Engineering, LJK, 38000, Grenoble</footnote>

<footnote>†Univ. Grenoble Alpes, INRAE, ETNA, 38000, Grenoble, France</footnote>
**Notation** For vector and tensors, we choose the following notation:

$$ \|\mathbf{x}\|^2 = \|\mathbf{x}\|^2 = \underbrace{x_i x^i}_{\text{indicial notation}} = \underbrace{\mathbf{x} \cdot \mathbf{x}}_{\text{tensor notation}} = \underbrace{x^\top x}_{\text{vector notation}} \quad (1) $$

The set $T_s^d$ is the set of symmetric tensors of order $s$ in dimension $d$. For a convex set $C \subset \mathbb{R}^n$, the normal cone to $C$ at $x$ is defined by

$$ N_C(x) = \{s \in \mathbb{R}^n \mid s^\top(y-x) \le 0, \text{ for all } y \in C\}. \quad (2) $$

For a convex function $\varphi$, its subdifferential is denoted by $\partial\varphi$. The indicator function of $C$ is denoted $i_C$ and we have $N_C = \partial i_C$. For a function $F: \mathbb{R}^n \to \mathbb{R}^n$, the inclusion

$$ -F(x) \in N_C(x) \text{ or equivalently, } F^\top(x)(x-y) \ge 0, \text{ for all } y \in C \quad (3) $$

defines a (finite dimensional) variational inequality. When the set $C$ is a cone, the variational inequality (3) is a complementarity problem denoted as

$$ C^* \ni F(x) \perp x \in C, \quad (4) $$

where $C^*$ is the dual cone of $C$. The complementarity symbol is defined by $x \perp y \iff x^\top y = 0$. When $C$ is equal to the positive orthant, i.e. $C = \mathbb{R}_+^n$, we write

$$ 0 \le F(x) \perp x \ge 0, \quad (5) $$

where the inequalities have to be understood component-wise for vectors.

For the thermo-mechanics of plasticity, we follow the notation introduced in Maugin (1992). This notation used are quite classical and are not defined in an exhaustive way, referring to the textbook when needed.

# 1 Introduction

The impact on elastoplastic solids and structures is an important research topic due to its numerous applications. In solid mechanics, the modelling of impact and plasticity phenomena is required for the reliability of metallic mechanical parts to impacts or in shot peening processes (Levers and Prior, 1998; Majzooobi et al., 2005; Nouguier-Lehon et al., 2013), for example. Generally, impacts, even at low speed, generate plastic deformations, possibly localised in the materials (Johnson, 1985), which can have strong consequences on their dimensioning. Examples include civil engineering applications such as vehicle impacts on buildings and structures (Chen, Wu, et al., 2021; Heng, Jia, et al., 2022; Heng, Hjiaj, et al., 2016, 2017), or ship impacts on harbour infrastructures (Guo et al., 2020; Sha et al., 2021). In the field of natural hazard mitigation, many applications require the joint modeling and simulation of plasticity and impact. Among them, we can cite building pounding (Langlade et al., 2021), the impact with soil foundation during earthquakes, the impact of rock blocks in the mountains against protective structures (walls, nets, or trees)(Bertrand et al., 2012; Di Giacinto et al., 2020; Dupire et al., 2016a,b), and finally impacts on rail, road or electrical transport infrastructures (Kaewunruen et al., 2018; Zeng et al., 2018).

The objective of this article is the modelling and the numerical simulation of the response of elastoplastic structures to impacts. To this end, a numerical method is proposed that takes into account one-sided contact (Signorini condition) and impact phenomena together with plasticity in a monolithic solver, while accounting for the non-smooth character of the dynamics. Furthermore, the goal is to ensure that the method is stable by formulating a discrete energy balance and yields an optimisation problem that is well-posed at each time step.

**Variational approach to elastoplastic systems.** One of the ideas of the proposed contribution is to return to the basis of the plasticity formulation, as variational principles, and to incorporate contact and impact conditions. In doing so, we will also be able to take advantage of a large set of numerical methods, coming from optimisation and mathematical programming and to be able to prove well-posedness results.
The formulation, through variational inequality of the plastic flow in solids is not new, and well known. Hill (1948, 1950) proposed a variational formulation for the plasticity problem in the case of perfect plasticity from the principle of maximum plastic work. Maier (1968a) formulated the plasticity constitutive laws as a mathematical problem which is also a variational inequality. The formulation as a mathematical programming problem, more precisely of complementarity type, and the maximum dissipation principle problem are the two sides of the same coin. Moreau (1970, 1971, 1974, 1976), one the father of convex analysis, recognised this structure of the problem as an inclusion in a normal cone to a convex set. Extending the works of Ziegler (1958, 1962) on normality for nonsmooth potentials, he postulated the existence of a pseudo-potential of dissipation given by a lower semi-continuous convex function. With the normality rule and the convexity of the yield criterion, the constitutive equations of plasticity are formulated as a variational inequality, or equivalently, a normal cone inclusion into a convex set. With the notion of subdifferential, the principle can easily be extended to non-smooth yield surfaces described by multiple yield criteria. Following this pioneering research, Halphen and Nguyen (1975) introduced the notion of generalised standard materials (GSM). The plasticity law with internal variables follows the normality rule between the plastic strain rate and the yield surface, by assuming pseudo-potential of dissipation. A recent account on this subject can be found in Houlsby (2019). Variational formulations enabled to prove the first mathematical results in plasticity (Duvaut and Lions, 1976; Moreau, 1976; Suquet, 1981). De Saxcé (1992) introduced the notions of bi-potential function and implicit standard materials (ISM), that extend the framework to non-associated plasticity (see Cheng et al. (2015) for a recent application of the method).

To complete the picture of variational approaches to plasticity, one can consider work that has sought to establish full variational principles. i.e., not only restricted to the principle of maximum dissipation. These variational principles, which seek to express the plasticity problem as an extremum (minimum, maximum or saddle point) were first written on the incremental problem, i.e. discretised in time (Comi and Perego, 1995; Maier et al., 1991; Ortiz and Stainier, 1999; Romano et al., 1993; Simo and Honein, 1990). The principle of maximum dissipation provides the basic building block for writing this type of problem. In this article, these principles will be recovered for the problem adding contact, identifying them from the discrete variational inequalities in the case where the assumptions of normal dissipativity, positive hardening and convexity of the elastic set are respected. Based on this work, variational principles have been developed in continuous time to provide a complete variational approach to the problem (Carstensen et al., 2002; Miehe, Apel, et al., 2002; Miehe, Schotte, et al., 2002; Miehe, 2002; Mielke, 2003, 2005; Mielke et al., 2002). The main advantage of this approach is that it makes it possible to use a wide range of numerical methods developed in the optimisation and mathematical programming community, which are robust, experienced and with global convergence (Boyd and Vandenberghe, 2004).

**Mathematical programming and optimization techniques.** Once the continuous problem is formulated in the GSM or ISM framework, the discretization in space, and possibly in time, leads to problems of finite-dimensional variational inequalities and complementarity problems (Facchinei and Pang, 2003). With some further assumptions such as the normality rule and the convexity of admissible sets, these problems can be recast into optimisation problems. Historically, Maier (1968a), with its coworkers (Capurso and Maier, 1970; Maier, 1968b, 1969) first considered mathematical programming techniques to formulate the plastic constitutive laws. Many subsequent research results have been done using these techniques (to cite a few, see Corradi, 1990; Donato and Maier, 1972; Feijoo and Zouain, 1988; Maier et al., 1991; Martin et al., 1987; Pereira et al., 1988; Romano et al., 1993; Smith, 1990 and references therein). More recently, several methods of modern numerical optimization have been successfully used to solve elasto-plastic problems: sequential quadratic programming method (Wieners, 2007), complementarity problem solvers (Tangaramvong et al., 2012; Zheng et al., 2020), accelerated gradient algorithms (Kanno, 2016; Shimizu and Kanno, 2020). Limit and shakedown analyses are problems for which an optimization approach of plasticity is natural, using second order cone programming technique and quadratic optimisation (Bisbos et al., 2005; Delbecq et al., 1977; Makrodimopoulos and Martin, 2005a,b; Mercier, 1976; Pastor, Thoré, et al., 2008; Pastor, Pastor, et al., 2015). In (Krabbenhøft, Lyamin, Sloan, and Wriggers, 2007; Krabbenhøft, Lyamin, Hjiaj, et al., 2005), developments are carried out for Mohr-Coulomb yield criterion in the context of cone variational inequality, extending the pioneering works of Berga and De Saxcé (1994).
and Hjiaj et al. (2003). In Zhang (2014) and subsequent work pieces (Meng et al., 2020; Zhang, Krabbenhøft, Pedroso, et al., 2013; Zhang, Sheng, et al., 2017; Zhang, Krabbenhøft, Sheng, et al., 2015; Zhao et al., 2022; Zhou et al., 2023), the mathematical programming approach of Krabbenhøft, Lyamin, Sloan, and Wriggers (2007) based on a saddle-point problem is further extended to new space-discretization techniques such as PFEM and non-associative flow rule for applications in geotechnics.

**Comparison with standard computational approach to elastoplastic systems: the return mapping method.** The numerical solution of elastoplasticity problems is generally based on the return mapping method, attributed to Wilkins (1963) and exposed in standard textbooks such as Simo and Hughes (1998) or Souza Neto et al. (2011). This method is very efficient and robust in many cases provided that the initial iterate is sufficiently close to the solution (see Scherzinger (2017) for recent advances). The local quadratic convergence is a major asset of the method, but it requires a rigorous work of construction and development of consistent tangent operators, which contain conditional statements that are difficult to deal with. Furthermore, as it is noted in Zheng et al. (2020), convergence of the return mapping algorithm is only local, meaning that we need to start from a good initial point preventing large load increments. Some globalisation methods exist, but they generally required a potential to minimize.

We want to stress that return mapping algorithm is not completely distinct from mathematical programming or optimisation methods. As explained by Simo and Honein (1990), the discrete flow and the hardening rules are the mathematical expression of the so-called closest-point-projection algorithm (Nguyen, 1977; Wilkins, 1963) and they recall that “from a mathematical standpoint, these algorithms appear to have been first studied in Moreau (1977), who coined the expression catching-up algorithms”. More generally, Christensen (2002b) shown that this method is a special case, for particular parameter choices, of a semismooth Newton method applied to a projection-based reformulation of the variational inequality that describes plasticity. A similar approach is followed in Bruno et al. (2020) using conic optimization. In conclusion, mathematical programming methods and the return mapping algorithm are not opposed. It seems clear today that the return mapping algorithm is a special case of mathematical programming methods, the semi-smooth Newton methods. Recognizing this point, it is possible to take advantage of the efficiency of Newton’s methods, to globalize them, to add other constraints such as the contact and to simplify the calculation of linearized operators and the numerical implementation.

**Nonsmooth dynamics with unilateral contact and impacts.** Another fundamental contribution of the nonsmooth mechanics community under the seminal impulse of Moreau is the dynamics of mechanical systems with contact, Coulomb friction and impact. Again, these constitutive laws can be written as variational inequalities and complementarity problems. Unilateral contact in mechanical systems implies jumps in the velocities. For discrete systems, measure differential equations are necessary to formulate rigorously the problem in order to design robust and efficient time-stepping schemes. Without entering into details, our work will be based on the nonsmooth contact dynamics method developed by Jean and Moreau (1987, 1992) (see also Jean (1999) and Moreau (1988b) for further developments and Acary and Brogliato (2008) and Dubois et al. (2018) for reviews of the method). In this time-stepping scheme, the velocity and the impulses are the main discrete unknown variables, yielding to a robust and consistent integration scheme. When only unilateral contact and impact are involved, the discrete problem can also be written as an optimisation problem. It seems therefore natural to propose an unified formulation for the problem of elasto-plastic structure subjected to contact conditions and impact, that enables the use of the powerful algorithms developed in the numerical optimisation community.

Of course, the NSCD method is not the only approach to model and simulate dynamics with contact, friction and impacts. First of all, we mention that one way of approaching non-smooth dynamics is to smooth the problem by penalisation (Belytschko and Neal, 1991; Fuente and Felippa, 1991; Meier et al., 2016; Yang, 2006), by singular mass method (Dabaghi et al., 2016; Di Stasio et al., 2021; Khenous et al., 2008; Renard, 2010) and by Nitsche method (Chouly, Hild, et al., 2015; Chouly and Renard, 2018). These approaches are effective and simple to implement. We have decided to stay in the framework of non-smooth dynamics for three main reasons: a) the consistency of the method in the case where the dynamics presents percussions (Dirac atoms in the forces, see Remark 1), b) the simplicity of writing the discrete variational inequality with plasticity and c) the possibility to extend the scheme to plastic
rigid models (Delbecq et al., 1977; Frémond et al., 1975). In the framework of non-smooth dynamics, the Moreau-Jean scheme is not the only scheme for the approximation of dynamic systems with impact. We can in particular quote the works of Paoli and Schatzman (2002a,b, 2007). Finally, we mention recent works that improve the performances of these historical schemes in nonsmooth dynamics (Time Discontinuous Galerkin, nonsmooth Newmark and nonsmooth generalized-α, time finite element) but which have not been used here for simplicity (Brüls et al., 2014, 2018; Capobianco and Eugster, 2018; Chen, Acary, et al., 2013; Dumont and Paoli, 2006; Schindler and Acary, 2013; Schindler, Rezaei, et al., 2015).

**Dynamics of elastoplastic structures with contact** There are finally very few works that consider the dynamics of elastoplastic problems with contact by formulating the problem as a variational inequality and using optimization algorithms to solve the problem. In Heng, Hjiaj, et al. (2016, 2017), Khan, Ahmad, et al. (2021), and Khan, Smith, et al. (2013), the effect of impact in elasto-plastic structures with the help of mathematical programming techniques is studied, but there is no monolithic algorithm to solve the contact and impact multipliers with the plasticity flow rule. In Christensen (2002b), a monolithic semi-smooth algorithm is developed but in the quasi-static case. Finally, the most advanced work in this direction is (Meng et al., 2020) but a) the non-smoothness of the dynamics is not taken into account and b) there are no results that guarantee that the system is well-posed.

**Contribution and outline of the article** The outline of the article is as follows. In Section 2, the equations of the dynamics of an elastoplastic system with unilateral contact are recalled. The goal of this section is, first, to formulate all the equations in a variational framework, respecting the principles of the thermodynamics and second, to write the elastoplastic law and contact law as a variational inequality. These two features allow one to state some energy principles. The formulation starts from the well-known framework of generalised standard materials and is extended to the case of unilateral contact in dynamics. Starting from the continuous time and space formulation, the model is discretized in Section 3. The space discretization is based on a standard iso-parametric finite element application to keep the presentation as simple as possible. We end up with a finite-dimensional differential variational inequality that takes into account the plasticity and the contact constraints at the end of Section 3.1. Since we then deal with finite-dimensional systems with finite masses, the discontinuities in the velocity imply impulsive forces. This is the reason why the system is recast with the help of differential measures and an impact law is introduced in the formulation to close the system equations. In Section 3.2, a time-stepping method is developed on the basis of the Moreau-Jean time stepping scheme. This scheme integrates the nonsmooth dynamics dealing consistently with the impulsive motion. Furthermore, it allows one to give some well-posedness results and a discrete energy balance that ensures that the scheme is well-posed and stable. Some of elements of the literature are relatively closed from our approach, but we differ mainly in the facts that a) we consider explicitly dynamics with impacts in a nonsmooth setting, b) we prove the well-posedness of the one-step problem and c) we give a discrete energy balance that ensures the stability of the scheme. Finally, in Section 4, some numerical illustrations are developed to show the interest of the proposed approach. The method is able to simulate multi-criteria elasto-plastic flow rules with contact constraints. The discrete energy balance highlights which physical processes dissipate the energy. Especially, the question of the coefficient of restitution and the energy dissipated by the impact is discussed with respect to the mesh size and the time-step.

## 2 Elasto-plasticity dynamics with unilateral contacts

The constitutive equations for such elasto-plastic systems involves constraints on the state variables and their derivatives. Typically, the stresses and the forces associated to the hardening parameters are constrained to be in an admissible set defined by the yield criteria. The most convenient framework for dealing with constraints is Convex Analysis as it has been pioneered by Moreau (1970, 1974, 1986) in the context of plasticity, friction and contact, introducing super-potentials and pseudo-potentials of dissipation. With Moreau's developments as a solid basis, Halphen and Nguyen (1975) invented the notion of Generalized Standard Materials (GSM), which makes it possible to write constitutive laws with internal variables on a rigorous thermodynamical basis. In this article, the notion of GSM is chosen for the constitutive models of
elasto-plasticity with hardening, with the Moreau-Ziegler assumption of normal dissipation in order to satisfy the principles of thermo-mechanics.

## 2.1 Principle of virtual powers

The material body is assumed to be an open set $\Omega \subset \mathbb{R}^d$, $d \in [1, 3]$ for $t \in [t_0, T]$ with a smooth boundary $\partial\Omega$. The volume Lebesgue measure is denoted by $dv(x)$ and the mass measure is denoted by $dm(x)$. We assume that mass measure has only a density with respect to $dv(x)$, denoted by $\rho(x)$, the density of the material, that is $dm(x) = \rho(x)dv(x)$. The displacement is denoted by $\mathbf{u}(x, t) : \Omega \times [0, T] \to \mathbb{R}^d$ and the velocity is denoted by $\mathbf{v}(x, t) : \Omega \times [0, T] \to \mathbb{R}^d$. We assume that we are in the framework of the small-perturbation hypothesis and the strain rate is given by $\dot{\varepsilon}(x, t)$, where the second order symmetric strain tensor $\varepsilon$ is given by $(\nabla$ denotes the gradient)

$$ \varepsilon(\mathbf{u}) = \nabla_s(\mathbf{u}) = \frac{1}{2}(\nabla^\top(\mathbf{u}) + \nabla(\mathbf{u})). \qquad (6) $$

**Principle of virtual power** The principle of virtual power in the setting of classical mechanics pos- tulates that the virtual power of inertial forces balances the virtual power of all other forces, volume, internal, external and contact applied to the system for all virtual velocity fields :

$$ \mathcal{P}_{\text{inertia}}^{*}(t) = \mathcal{P}_{\text{v}}^{*}(t) + \mathcal{P}_{\text{int}}^{*}(t) + \mathcal{P}_{\text{ext}}^{*}(t) + \mathcal{P}_{\text{contact}}^{*}(t), \quad \forall \mathbf{v}^{*} \qquad (7) $$

where $\mathbf{v}^*$ is the virtual velocity and the virtual powers are defined as follows: the virtual power of inertial forces and the virtual power of internal forces as

$$ \mathcal{P}_{\text{inertia}}^{*}(t) = \int_{\Omega} \mathbf{v}^{*}(x, t) \cdot \dot{\mathbf{v}}(x, t) \, dm(x), \quad \mathcal{P}_{\text{int}}^{*}(t) = - \int_{\Omega} \dot{\mathbf{\varepsilon}}^{*}(x, t) : \boldsymbol{\sigma}(x, t) \, dv(x), \qquad (8) $$

where $\boldsymbol{\sigma} \in T_2^d$ is the Cauchy stress tensor, the virtual power of volume forces $\mathbf{f}$ and the power of applied contact forces $\mathbf{t}(x, t)$ on $\partial\Omega_F \subset \partial\Omega$ as

$$ \mathcal{P}_{\text{v}}^{*}(t) = \int_{\Omega} \mathbf{v}^{*}(x, t) \cdot \mathbf{f}(x, t) \, dv(x), \quad \mathcal{P}_{\text{ext}}^{*}(t) = - \int_{\partial\Omega_F} \mathbf{v}^{*}(x, t) \cdot \mathbf{t}(x, t) \, ds(x), \qquad (9) $$

and the unknown contact forces $\mathbf{r}(x, t)$ on $\partial\Omega_C \subset \partial\Omega$ as

$$ \mathcal{P}_{\text{contact}}^{*}(t) = \int_{\partial\Omega_C} \mathbf{v}^{*}(x, t) \cdot \mathbf{r}(x, t) \, ds(x). \qquad (10) $$

**Energy balance** The kinetic energy is defined as

$$ \mathcal{T}(t) = \frac{1}{2} \int_{\Omega} \mathbf{v}(x, t) \cdot \mathbf{v}(x, t) \, dm(x), \qquad (11) $$

and its time derivative, when it is defined, gives the power of inertial forces as

$$ \frac{d}{dt} \mathcal{T}(t) = \mathcal{P}_{\text{inertia}}(t) = \int_{\Omega} \mathbf{v}(x, t) \cdot \dot{\mathbf{v}}(x, t) \, dm(x). \qquad (12) $$

Choosing the virtual field $\mathbf{v}^*$ equal to the actual velocity field $\mathbf{v}$, we obtain the energy balance from the principle of virtual power:

$$ \frac{d}{dt} \mathcal{T}(t) = \mathcal{P}_{\text{v}}(t) + \mathcal{P}_{\text{int}}(t) + \mathcal{P}_{\text{ext}}(t) + \mathcal{P}_{\text{contact}}(t). \qquad (13) $$

**Remark 1.** We have implicitly made the assumption in this section that the continuous dynamics in space and time has no impulsive forces (Dirac atoms). This assumption is usually found in the mechanical literature on continuum mechanics with discontinuities, especially, with shocks waves in elasto-plastic solids, for
instance in (Davison, 2008; Germain and Lee, 1973; Lee and Liu, 1964; Mandel, 1964). In a seminal paper, Lebeau and Schatzman (1984) show the existence and uniqueness of solutions and prove the conservation of energy for an elastic half-space with Signorini condition, where there is no impulsive forces, or impulsive stresses. In dimension one ($d = 1$) it corresponds to the equation of an elastic bar and in dimension two ($d = 2$) to that of a membrane. For a more general object, in particular a manifold of codimension 1 in $\mathbb{R}^d$ (beams, plates, shells), the question is open and simple examples show that impulsive forces are possible.

The assumption that the set of discontinuities is of zero measure with respect to $dm$ is also questionable, even if the mass has only a density with respect to the volume measure (or surface measure). When one-sided constraints are imposed on elastic bodies in dynamics, the question of the existence and regularity of solutions remains open.

## 2.2 Elasto-plasticity in the Generalized Standard Material (GSM) setting

As usual under the small perturbation hypothesis for the elasto-plastic solids, the strain tensor is decomposed in an additive way into an elastic part $\varepsilon^e$, i.e. the reversible part of the strain, and a plastic part $\varepsilon^p$ as follows

$$ \varepsilon = \varepsilon^e + \varepsilon^p. \qquad (14) $$

**Free energy and state laws** We start by the definition of the free energy per unit of volume $\psi(\varepsilon^e, \alpha)$ in the isothermal setting (the free energy does not depend on the temperature $T$) assuming a separation in $\varepsilon^e$ and $\alpha$:

$$ \psi(\varepsilon^e, \alpha) = \psi^e(\varepsilon^e) + \psi^p(\alpha), \qquad (15) $$

where $\alpha$ is a tensor of hardening parameters, $\psi^e$ is the elastic potential energy and $\psi^p$ is the plastic potential energy due to the hardening process. The state laws govern the generalized (driving) forces $\sigma(x, t) \in T^d$ and $a(x, t)$ as the derivative of the free energy:

$$ \sigma = \frac{\partial \psi(\varepsilon^e, \alpha)}{\partial \varepsilon^e}, \quad \text{and} \quad a = -\frac{\partial \psi(\varepsilon^e, \alpha)}{\partial \alpha}. \qquad (16) $$

Following the simplest case, we choose quadratic potentials,

$$ \psi^e(\varepsilon^e) = \frac{1}{2}\varepsilon^e : \mathbf{E} : \varepsilon^e \quad \text{and} \quad \psi^p(\alpha) = \frac{1}{2}\alpha \cdot \mathbf{D} \cdot \alpha \qquad (17) $$

with $\mathbf{E}$ a symmetric 4-th order elastic tensor and $\mathbf{D}$ a symmetric second order tensor, the state laws are then linear:

$$ \sigma = \mathbf{E} : \varepsilon^e \quad \text{and} \quad a = -\mathbf{D} \cdot \alpha. \qquad (18) $$

It is worth noting that, although the choice of quadratic potentials simplifies the presentation, this choice is not compulsory and taking into account more general non-linear potentials does not pose specific problems.

**Intrinsic dissipation** From the first and second principles of thermodynamics, the intrinsic dissipation per unit of volume $d$ in a isothermal setting, defined by

$$ d = \sigma : \dot{\varepsilon} - \dot{\psi}(\varepsilon^e, \alpha), \qquad (19) $$

must be positive. Developing $\dot{\psi}(\varepsilon^e, \alpha)$, we get

$$ d = \sigma : \dot{\varepsilon} - \left[ \frac{\partial \psi(\varepsilon^e, \alpha)}{\partial \varepsilon^e} : \dot{\varepsilon}^e + \frac{\partial \psi(\varepsilon^e, \alpha)}{\partial \alpha} \cdot \dot{\alpha} \right]. \qquad (20) $$

With the state laws defined above, this can be simplified in

$$ d = \sigma : \dot{\varepsilon} - [\sigma : \dot{\varepsilon}^e - a \cdot \dot{\alpha}] = \sigma : [\dot{\varepsilon} - \dot{\varepsilon}^e] + a \cdot \dot{\alpha} = \sigma : \dot{\varepsilon}^p + a \cdot \dot{\alpha}. \qquad (21) $$
The total intrinsic dissipation (or dissipated power) $D(t)$ and the total free energy $\Psi$ are defined as

$$D(t) = \int_{\Omega} d(x, t) \, dx, \quad \text{and} \quad \Psi(t) = \int_{\Omega} \psi(\epsilon^e(x, t), \alpha(x, t)) \, dx. \qquad (22)$$

Then, the power of internal forces is expressed as

$$P_{\text{int}}(t) = - \int_{\Omega} d(x, t) \, dx - \int_{\Omega} \dot{\psi}(\epsilon^e(x, t), \alpha(x, t)) \, dx, \qquad (23)$$

that is

$$P_{\text{int}}(t) = -D(t) - \dot{\Psi}(t). \qquad (24)$$

For a detailed energy balance in the context of elastoplasticity with hardening, we refer to Chrysochoos et al. (1989).

**Normality rule and complementary laws** From the second principle of the thermodynamics, the intrinsic dissipation must satisfy the Clausius-Duhem inequality, that is, in an isothermal setting

$$d = \sigma : \dot{\epsilon}^p + a \cdot \dot{\alpha} \ge 0. \qquad (25)$$

A convenient way to ensure that complementary laws will imply this inequality is to postulate the existence of a pseudo potential of dissipation from which the rate of changes of the $\epsilon^p$ and $\alpha$ derive (Moreau, 1974). In the case of elasto-plastic models that obey to the normality rule (associated plasticity), we introduce a convex set $C$ as the set of admissible stresses $\sigma$ and forces $a$, and the complementary laws are defined by the following differential inclusion

$$\left(\dot{\epsilon}^p \quad \dot{\alpha}\right) \in N_C\left(\frac{\sigma}{a}\right). \qquad (26)$$

The complementary laws (26) imply the Clausius-Duhem inequality (25) if the origin belongs to $C$. Using the definition of the normal cone (2), the complementary laws are equivalent to the following variational inequality

$$\dot{\epsilon}^p : (\sigma - \hat{\sigma}) + \dot{\alpha} \cdot (a - \hat{a}) \ge 0, \quad \forall (\hat{\sigma}, \hat{a}) \in C, \qquad (27)$$

also known as the Hill principle of maximum dissipation

$$\dot{\epsilon}^p : \sigma + \dot{\alpha} \cdot a \ge \dot{\epsilon}^p : \hat{\sigma} + \dot{\alpha} \cdot \hat{a}, \quad \forall (\hat{\sigma}, \hat{a}) \in C. \qquad (28)$$

**Remark 2.** We would like to stress that $C$ is a general convex set. Naturally, in the numerical practice, we need a handable description in order to evaluate the model. A usual description is obtained when the set $C$ is finite representable with a set of smooth functions $f$ as

$$C = \{(\sigma, a) | f(\sigma, a) \ge 0\}. \qquad (29)$$

Some qualification conditions are generally added to ensure that the normal cone can be coherently computed (see for instance Brogliato et al., 2006 in the context of dynamical systems). Note that the multi-criteria case does not pose any particular problems. The set $C$ is not a priori a polyhedral set and the yield functions $f$ can be non-linear.

## 2.3 Unilateral contact

In this work, the unilateral contact is considered as perfect unilateral constraint. For a contact point $C_A$ of position $x \in \partial\Omega_C$ candidate to contact, we assume that we are able to associate a unique contact point $C_B$ of position $\tilde{x} = \Gamma(x) \in \partial\Omega_C$ such that the normal gap function is defined as

$$g_n(x) = (\tilde{x} - x) \cdot n(x) = (\Gamma(x) - x) \cdot n(x), \qquad (30)$$
<figure><img src="image_3.png" /><figcaption>Figure 1: A contact pair and a local contact frame.</figcaption></figure>

where $n(x)$ is the outward normal vector at $x$ (see Figure 1). When the contact occurs with a fixed obstacle, the point $\tilde{x}$ is possibly defined on the boundary of this obstacle. For a detailed description of the practical computation of $\Gamma$, we refer to Wriggers (2006, Section 4.1) and Konyukhov and Schweizerhof (2012). The Signorini condition can be viewed as an additional state law associated with the additional state variable $g_N$. The surface free energy $\psi^s$ of the system is defined with the indicator function of $\mathbb{R}_+$, $i_{\mathbb{R}_+}$, taken at $g_N$

$$ \psi^s(g_N) = i_{\mathbb{R}_+}(g_N). \qquad (31) $$

The Signorini condition is then obtained by sub-differentiating

$$ -r_N \in \partial_{g_N} \psi^s(g_N) = \partial i_{\mathbb{R}_+}(g_N) \iff -r_N \in N_{\mathbb{R}_+}(g_N). \qquad (32) $$

where $r_N$ the normal reaction at contact. The definition of the subdifferential of the indicator function leads to the Signorini condition in terms of complementarity:

$$ 0 \le g_N \perp r_N \ge 0. \qquad (33) $$

It is also possible to write a stronger condition at the velocity level that implies the standard Signorini condition as

$$ -r_N \in N_{T_{\mathbb{R}_+}(g_N)}(v_N), \qquad (34) $$

where $T_{\mathbb{R}_+}(g_N)$ is the tangent cone of $\mathbb{R}_+$ (see Moreau (1988b)) and the normal relative velocity $v_N$ is defined as the time derivative of the gap function

$$ v_N(x(t)) = \frac{d}{dt} g_N(x(t)). \qquad (35) $$

Equivalently, the complementarity condition at the velocity level is

$$ 0 \le r_N \perp v_N \ge 0 \text{ if } g_N = 0, \text{ else } r_N = 0. \qquad (36) $$

The contact force $r(x, t)$ are related to the normal contact reaction by duality using the virtual power of contact forces as

$$ \mathcal{P}_{\text{contact}}^{*}(t) = \int_{\partial \Omega_{C}} v^{*}(x, t) \cdot r(x, t) \, ds(x) = \int_{\partial \Omega_{C}} v_{N}^{*}(x, t) r_{N}(x, t) \, ds(x). \qquad (37) $$

**Power of contact forces.** The power done by the contact forces vanishes almost everywhere due to (36), and does not dissipate or release energy:

$$ \mathcal{P}_{\text{contact}}(t) = \int_{\partial \Omega} r(x, t) \cdot v(x, t) \, ds(x) = \int_{\partial \Omega} r_N(x, t) v_N(x, t) \, ds(x) = 0. \qquad (38) $$
# 3 Space and Time discretizations

## 3.1 Standard finite element discretization

In this section, we use a standard finite element discretization of the elasto-plastic bodies. We have chosen to leave it as simple as possible since this is not the main objective of the paper. For a comprehensive review of the specificity of finite element application to plasticity, we refer to Nodargi (2019) where enhanced mixed FE discretization can be found. The discrete equation are obtained by summing the contributions of the elements yielding

$$M\dot{v}(t) + f_{\text{int}}(t) = f_{\text{ext}}(t) + r(t). \qquad (39)$$

Details of the notation are given in Appendix A. The Dirichlet boundary condition are taken into account by modifying the mass matrix and the vector of forces accordingly, following the method described in Zienkiewicz et al. (2005, Chapter 2, Section 2.5).

**Equilibrium matrix** In order to obtain a convenient matrix notation for the internal forces, we note that the internal forces are linear in $\sigma_e$. A discretization of the stresses $\sigma_e(x, t)$ in the element should result in a linear formula expressed as a matrix product with the introduction of an equilibrium matrix. In the sequel, we use the discretization directly based in the Gauss quadrature rule for evaluating the integrals. More enhanced approaches have been discussed in (Corradi, 1985). Using the discrete equilibrium matrix $B$ defined in Appendix B, we get the matrix of the internal forces as

$$f_{\text{int},e} = B_e^T \sigma_e, \text{ and } f_{\text{int}} = \sum_{e=1}^{n_{\text{el}}} f_{\text{int},e} = B^T \sigma. \qquad (40)$$

Hence, the space-discretized equation of motion (39) reads as

$$M\dot{v}(t) + B^T \sigma(t) = f_{\text{ext}}(t) + r(t). \qquad (41)$$

**Remark 3.** In Appendix B, the strain and stress discrete values are gathered into column vector using a scaling by the square root of the Gauss weights. This allows a convenient matrix notation, where the structure equilibrium matrix is directly the adjoint operator of the strain-displacement matrix (see for instance (Hughes, 1980; Pellegrino, 1993; Reddy and Martin, 1991)). It would also have been possible to introduce another diagonal matrix containing the Gauss weight leading to a direct gathering without scaling of the strain and stress values at Gauss points.

**Discretization of the elasto-plastic equations** Since the stress and strain vectors are evaluated at Gauss points, we choose to apply directly the elastoplastic model at the Gauss points. It results in the following system of equation in matrix vector notation:

$$\left\{ \begin{array}{l} \sigma = E\varepsilon^e = E(\varepsilon - \varepsilon^p) \\ a = -D\alpha \\ \left(\dot{\varepsilon}^p\right) \in N_C\left(\frac{\sigma}{a}\right), \end{array} \right. \qquad (42)$$

where $E$ is the elasticity matrix obtained by reformulating the elasticity tensor $E$ in matrix with Voigt notation for instance. The same applies for $D$. The only point of attention is the definition of the convex set $C$ that is a Cartesian product of the set $C$ for each Gauss points and taking into account the scaling due to the Gauss weights introduced in appendix B.

**Discretization of the contact conditions** The discretization of the contact conditions is a delicate subject. An accurate discretization of the contact traction $r(x, t)$ call for the use of mortar finite element techniques (see for instance (Belgacem et al., 1998; Popp et al., 2009)). Since it is beyond the scope of the article and to keep the presentation as simple as possible, the discretization is made by choosing to apply
the Signorini condition to a set of discrete points on the contact boundary associated with point load contact forces. These discrete points can be simply a node of the mesh, or a point for which the position is expressed in terms of the node displacements in a iso-parametric setting, with the shape function matrix. For a contact point labeled by $\alpha \in [1... m]$, the gap function is written as a function of the nodal displacement $u$ as $g_N^\alpha(u)$. The normal relative velocity is then given by

$$v_N^\alpha = \frac{d}{dt} g_N(u) = \nabla_u^\top (g_N^\alpha(u)) v := H^{\alpha, \top}(u)v. \quad (43)$$

The vector $v_N$ is defined by collecting all the relative velocities at the contacts points, $v_N := \text{col}(v_N^\alpha, \alpha \in [1, m])$, and the matrix $H$ is defined as

$$v_N = H^\top(u)v. \quad (44)$$

By duality, the contact forces are expressed as

$$r^\alpha = H^\alpha(u)r_N^\alpha. \quad (45)$$

The Signorini condition at the velocity level is then written as

$$-r_N \in N_{T_{\mathbb{R}^m_+}(g_N)}(v_N). \quad (46)$$

**Nonsmooth dynamics for finite-freedom mechanical systems**

Once the space-discretization is performed, we end up with a finite freedom mechanical system that is subjected to unilateral constraints. In standard finite element discretization, all degrees of freedom possess a finite mass. Since the system may encounter some jumps in the velocities, the fact that we have finite masses leads to the possibility to have percussion in the dynamics. In other words, the assumption that we made in Section 2.1 is no longer valid and we need to take into account the velocity jump into the equation of motion. To this aim, we consider that the velocity is a right continuous function of bounded variations with respect to time and the acceleration is therefore considered as a differential measure $dv$ associated with $v$ (Moreau, 1988a,b). We will further assume that the reaction due to the unilateral constraint is also a measure denoted by $d\dot{i}_N$ that has the same support as $dv$. Assuming that the internal forces have only a density with the Lebesgue measure on the real line, the equation of motion

$$M\dot{v}(t) + B^T \sigma(t) = f_{\text{ext}}(t) + H(u(t))r_N(t), \quad (47)$$

is extended as measure equation as

$$M dv + B^\top \sigma(t) dt = f_{\text{ext}}(t) dt + H(u(t)) d\dot{i}_N. \quad (48)$$

If we define $r_N$ as the density of $d\dot{i}_N$ with respect to the Lebesgue measure, the equation (47) appears as the equation of motion, valid almost everywhere. At any time point $t_i$, the measure equation amounts to solving the impact equation

$$M(v(t_i) - v^-(t_i)) = H(u(t_i))p_{i,N}, \quad (49)$$

where the percussion $p_{i,N}$ is the normal percussion, i.e. the density of $d\dot{i}_N$ with respect to the Dirac atom at $t_i$. In order to solve the system with a new unknown $p_{i,N}$, we need to introduce an impact law. The simplest choice is the Newton impact law

$$v_N^\alpha(t_i) = -ev_N^{\alpha, -}(t_i), \quad \text{if } g_N^\alpha(u(t_i)) = 0 \text{ and } v_N^{\alpha, -}(t_i) \le 0, \quad (50)$$

where $e$ is the Newton coefficient of restitution. Following the seminal work of Moreau (1988b), the Newton impact law and the contact law expressed in (46) are formulated as a measure inclusion as

$$-d\dot{i}_N \in N_{T_{\mathbb{R}^m_+}(g_N(t))}(v_N(t) + ev_N^-(t)). \quad (51)$$
**Summary of the space-discretized equations** To sum-up, after the space discretization, the equations of the elasto-plastic system with unilateral contact for $t \in [t_0, T]$ read as

$$ egin{cases} 
Mdv + B^\top \sigma(t)dt = f_{\text{ext}}(t)dt + H(u(t))di_N \\ 
\dot{u}(t) = v(t) \\ 
\sigma(t) = E(\varepsilon(t) - \varepsilon^p(t)) \\ 
a(t) = -D\alpha(t) \\ 
\begin{pmatrix} \dot{\varepsilon}^p(t) \\ \dot{\alpha}(t) \end{pmatrix} \in N_C \begin{pmatrix} \sigma(t) \\ a(t) \end{pmatrix} \\ 
v_N(t) = H^\top(u(t))v(t) \\ 
-di_N \in N_{T_{\mathbb{R}^m_+}(g_N(t))} (v_N(t) + ev_N^-(t)), 
\end{cases} 
\quad (52) $$

with the initial conditions $v(t_0) = v_0$, $u(t_0) = u_0$ and $\varepsilon^p(t_0) = 0$ for some $t_0$.

In order to prepare the time discretization of (52), the differential equation are separated from the algebraic equations and the normal cone inclusion introducing two slack variables $y = -\dot{\alpha}$ and $z = -\dot{\varepsilon}^p$. Assuming that the elastic law is valid at time $t_0$, that is, $\sigma(t_0) = E\varepsilon(t_0)$, the elasticity can be expressed as

$$ \dot{\sigma}(t) = E(\dot{\varepsilon}(t) - \dot{\varepsilon}^p(t)) = E(Bv(t) + z(t)), \quad (53) $$

or equivalently, with the compliance matrix $S = E^{-1}$ as

$$ S\dot{\sigma}(t) = Bv(t) + z(t). \quad (54) $$

With the slack variable $y$ and assuming that $a(t_0) = -D\alpha(t_0)$, the linear hardening law reads

$$ D^{-1}\dot{a}(t) = -\dot{\alpha}(t) = y(t). \quad (55) $$

Finally, with the slack variables $y$ and $z$, the normal cone inclusion for the plasticity flow rule is given by

$$ - \begin{pmatrix} z \\ y \end{pmatrix} \in N_C \begin{pmatrix} \sigma \\ a \end{pmatrix}. \quad (56) $$

The complete system of equations for the elasto-plastic evolution with contact reduces to

$$ egin{cases} 
Mdv + B^\top \sigma(t)dt = f_{\text{ext}}(t)dt + H(u(t))di_N \\ 
\dot{u}(t) = v(t) \\ 
S\dot{\sigma}(t) = Bv(t) + z(t) \\ 
D^{-1}\dot{a}(t) = y(t) \\ 
- \begin{pmatrix} z(t) \\ y(t) \end{pmatrix} \in N_C \begin{pmatrix} \sigma(t) \\ a(t) \end{pmatrix} \\ 
v_N(t) = H^\top(u(t))v(t) \\ 
-di_N \in N_{T_{\mathbb{R}^m_+}(g_N(t))} (v_N(t) + ev_N^-(t)), 
\end{cases} 
\quad (57) $$

with the initial conditions $v(t_0) = v_0$, $u(t_0) = u_0$ and $\varepsilon^p(t_0) = 0$. Let us recall that we assume that $a(t_0) = -D\alpha(t_0)$ and $\sigma(t_0) = E\varepsilon(t_0)$.

## 3.2 Time discretization with an extended Moreau-Jean time-stepping scheme

The time-discretization of the system (57) is given by an extension of the Moreau-Jean scheme for non-smooth dynamics (Jean and Moreau, 1987) and by the catching-up algorithm for plasticity (Moreau, 1976). The dynamics equation in the first line of (57) is discretized using the definition of the differential measure over a time interval $(t_k, t_{k+1}]$ as

$$ \int_{(t_k, t_{k+1}]} Mdv = M(v^+(t_{k+1}) - v^+(t_k)) = - \int_{t_k}^{t_{k+1}} B^\top \sigma(t)dt + \int_{t_k}^{t_{k+1}} f_{\text{ext}}(t)dt + \int_{(t_k, t_{k+1}]} H(u(t))di_N. $$
(58)

As in the original Moreau–Jean scheme, the impulse of the reaction measure is kept as primary variable to ensure consistency of the scheme when an impact occurs:

$$ \int_{(t_k, t_{k+1}]} H(u(t))di_N \approx H(u_k)p_{N,k+1}. \quad (59) $$

The explicit approximation of $H(u(t))$ is justified by the fact that we assume to be under the small perturbation hypothesis. Using the following standard notation $x_k$ and $x_{k+1}$ for the approximations of the function $x(t)$ at time $t_k$ and $t_{k+1}$, the remaining nonimpulsive terms of (141) are approximated with a $\theta$-method, for $\theta > 0$ as

$$ M(v_{k+1} - v_k) + hB^\top \sigma_{k+\theta} = h f_{\text{ext},k+\theta} + H(u_k) p_{N,k+1}, \quad (60) $$

where the following notation is used $x_{k+\theta} = \theta x_{k+1} + (1-\theta)x_k$ for any variable $x$. The non impulsive dynamical equations in (57) are also discretized using a $\theta$-method

$$ \begin{aligned} u_{k+1} &= u_k + h v_{k+\theta}, \\ S(\sigma_{k+1} - \sigma_k) - h B v_{k+\theta} &= h z_{k+\theta}, \\ D^{-1}(a_{k+1} - a_k) &= h y_{k+\theta}. \end{aligned} \quad (61) $$

The normal cone inclusion governing the plastic flow rule is discretized as follows

$$ - \begin{pmatrix} z_{k+\theta} \\ y_{k+\theta} \end{pmatrix} \in N_C \begin{pmatrix} \sigma_{k+\theta} \\ a_{k+\theta} \end{pmatrix}. \quad (62) $$

For the contact impact conditions, the active constraints at the velocity level must be selected. To this aim, the following index set of contact is introduced

$$ \mathcal{I}_k = \{ \alpha \mid g_{N,k}^{\alpha} \le 0, v_{N,k}^{\alpha} \le 0 \} \quad (63) $$

with $v_{N,k}^{\alpha} = H^{\alpha,\top} v_k$. Hence, the impact law is written as

$$ \left\{ \begin{array}{ll} -p_{N,k+1}^{\alpha} \in N_{\mathbb{R}^+}(v_{N,k+1}^{\alpha} + e v_{N,k}^{\alpha}), & \text{for } \alpha \in \mathcal{I}_k \\ p_{N,k+1}^{\alpha} = 0, & \text{otherwise.} \end{array} \right. \quad (64) $$

where $v_{N,k+1}^{\alpha} = H^{\alpha,\top} v_{k+1}$. By collecting only the contact variables for the index set $\mathcal{I}_k$ as $v_N := \text{col}(v_N^{\alpha}, \alpha \in \mathcal{I}_k)$ and $p_N := \text{col}(p_N^{\alpha}, \alpha \in \mathcal{I}_k)$, the contact impact law is given by

$$ -p_{N,k+1} \in N_{\mathbb{R}^+}(v_{N,k+1} + e v_{N,k}), \quad (65) $$

where $m$ is the cardinal of $\mathcal{I}_k$. Using the reduced matrix $H = \text{col}(H^{\alpha}, \alpha \in \mathcal{I}_k)$, we also have $v_{N,k} = H^{\top} v_k$ and

$$ M(v_{k+1} - v_k) + h B^\top \sigma_{k+\theta} = h f_{\text{ext},k+\theta} + H p_{N,k+1}. \quad (66) $$

**Summary of the space and time discretized equations** Altogether, the space and time discretized equations are given by

$$ \left\{ \begin{array}{l} M(v_{k+1} - v_k) + h B^\top \sigma_{k+\theta} = h f_{\text{ext},k+\theta} + H p_{N,k+1} \\ S(\sigma_{k+1} - \sigma_k) - h B v_{k+\theta} = h z_{k+\theta}, \\ D^{-1}(a_{k+1} - a_k) = h y_{k+\theta} \\ v_{N,k+1} = H^\top v_{k+1} \\ - \begin{pmatrix} z_{k+\theta} \\ y_{k+\theta} \\ p_{N,k+1} \end{pmatrix} \in N_{C \times \mathbb{R}_+^m} \begin{pmatrix} \sigma_{k+\theta} \\ a_{k+\theta} \\ v_{N,k+1} + e v_{N,k} \end{pmatrix} \end{array} \right. \quad (67) $$
where $u_{k+1} = u_k + h v_{k+\theta}$ can be solved afterwards, and $u_0 = u(t_0)$, $v_0 = v(t_0)$, $z_0 = 0$.

Using the relation $x_{k+\theta} - x_k = \theta(x_{k+1} - x_k)$, the discretized equations can be formulated in terms of variables approximated at $t_{k+\theta}$.

$$ egin{cases} 
M(v_{k+\theta} - v_k) + h\theta B^\top \sigma_{k+\theta} = h\theta f_{\text{ext},k+\theta} + \theta H p_{N,k+1} \\ 
S(\sigma_{k+\theta} - \sigma_k) - h\theta B v_{k+\theta} = h\theta z_{k+\theta}, \\ 
D^{-1}(a_{k+\theta} - a_k) = h\theta y_{k+\theta} \\ 
\theta v_{N,k+1} = H^\top v_{k+\theta} - (1-\theta)v_{N,k} \\ 
- \begin{pmatrix} z_{k+\theta} \\ y_{k+\theta} \\ v_{N,k+1} + e v_{N,k} \end{pmatrix} \in N_{C \times \mathbb{R}_+^m} \begin{pmatrix} \sigma_{k+\theta} \\ a_{k+\theta} \\ p_{N,k+1} \end{pmatrix} 
\end{cases} 
\quad (68) $$

### 3.3 Well-posedness

To study the well posedness to the problem (68), we propose to show that it corresponds to the optimality conditions of a saddle point problem, and then to show that this saddle point problem admits a unique solution.

**A saddle point problem.** Let us consider the following saddle point problem

$$ \min_{v, \dot{\epsilon}} \max_{\sigma, a} \frac{1}{2} (v - v_k)^\top M (v - v_k) - \frac{1}{2} (\sigma - \sigma_k)^\top S (\sigma - \sigma_k) - \frac{1}{2} (a - a_k)^\top D^{-1} (a - a_k) + h \theta \sigma^\top \dot{\epsilon} - h \theta f_{\text{ext},k+1}^\top v $$

$$ \begin{array}{ll} \text{s.t.} & Bv = \dot{\epsilon} \\ & \theta v_N = H^\top v - (1-\theta)v_{N,k} \\ & \begin{pmatrix} \sigma \\ a \\ v_N + e v_{N,k} \end{pmatrix} \in C \times \mathbb{R}_+^m. \end{array} \quad (69) $$

**Proposition 1.** Using the first order optimality conditions of (69), the solutions $(v, \dot{\epsilon}, \sigma, a, v_N)$ of (69) are solutions of (68) with $(v_{k+\theta}, z_{k+\theta}, \sigma_{k+\theta}, a_{k+\theta}, v_{N,k+1}) = (v, -\dot{\epsilon}, \sigma, a, v_N)$.

In order to derive the optimality conditions for (69), we consider the Lagrangian function

$$ \begin{aligned} \mathcal{L}(v, \dot{\epsilon}, v_N, \sigma, a, \lambda, \mu) = & \frac{1}{2}(v - v_k)^\top M(v - v_k) - \frac{1}{2}(\sigma - \sigma_k)^\top S(\sigma - \sigma_k) - \frac{1}{2}(a - a_k)^\top D^{-1}(a - a_k) \\ & + h\theta\sigma^\top\dot{\epsilon} - h\theta f_{\text{ext},k+1}^\top v + \lambda^\top(Bv - \dot{\epsilon}) + \mu^\top(\theta v_N - H^\top v + (1-\theta)v_{N,k}) \\ & - h\theta\Psi_C\begin{pmatrix} \sigma \\ a \end{pmatrix} + \theta^2 \Psi_{\mathbb{R}_+^m}(v_N + e v_{N,k}). \end{aligned} \quad (70) $$

The optimality conditions of the problem (69) is given by $0 \in \partial \mathcal{L}(v, \dot{\epsilon}, \sigma, a, v_N, \lambda, \mu)$. Computing the
subgradients, we obtain for the optimality conditions

$$ (\nabla_v \mathcal{L}): \quad 0 = M(v - v_k) - h\theta f_{\text{ext},k+1}^{\top} + B^{\top}\lambda - H\mu $$

$$ (\nabla_{\dot{\epsilon}} \mathcal{L}): \quad 0 = h\theta\sigma - \lambda $$

$$ (\partial_{\sigma} \mathcal{L}): \quad 0 \in -S(\sigma - \sigma_k) + h\theta\dot{\epsilon} - h\theta\partial_{\sigma} \Psi_C \begin{pmatrix} \sigma \\ a \end{pmatrix} $$

$$ (\partial_a \mathcal{L}): \quad 0 \in -D^{-1}(a - a_k) - h\theta\partial_a \Psi_C \begin{pmatrix} \sigma \\ a \end{pmatrix} \qquad (71) $$

$$ (\partial_{v_N} \mathcal{L}): \quad 0 \in \theta\mu + \theta^2\partial\Psi(v_N + ev_{N,k}) $$

$$ (\nabla_{\lambda} \mathcal{L}): \quad 0 = Bv - \dot{\epsilon} $$

$$ (\nabla_{\mu} \mathcal{L}): \quad 0 = \theta v_N - H^{\top}v + (1 - \theta)v_{N,k} $$

Note that $h\theta\sigma$ appears as the Lagrange multiplier that enforces the condition $Bv = \dot{\epsilon}$. Introducing the variables $z, y, p_N$ such that

$$ - \begin{pmatrix} z \\ y \\ p_N \end{pmatrix} \in \partial \Psi_{C \times \mathbb{R}_+^m} \begin{pmatrix} \sigma \\ a \\ v_N + ev_{N,k} \end{pmatrix}, \qquad (72) $$

and simplifying the equations, we obtain

$$ \left\{ \begin{array}{l} M(v - v_k) + h\theta B^{\top} \sigma = h\theta f_{\text{ext},k+1}^{\top} + \theta H p_N \\ S(\sigma - \sigma_k) - h\theta Bv = h\theta z \\ D^{-1}(a - a_k) = h\theta y \\ \theta v_N = H^{\top}v - (1 - \theta)v_{N,k} \\ - \begin{pmatrix} z \\ y \\ p_N \end{pmatrix} \in \partial \Psi_{C \times \mathbb{R}_+^m} \begin{pmatrix} \sigma \\ a \\ v_N + ev_{N,k} \end{pmatrix} \end{array} \right. \qquad (73) $$

This completes the proof.

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

Let us introduce the following assumptions

**Assumption 1.** The matrices $M, S$ and $D$ are symmetric definite positive matrices.

**Assumption 2.** It exists $v^0, \sigma^0, a^0$ such that

$$ \left\{ \begin{array}{l} \begin{pmatrix} \sigma^0 \\ a^0 \end{pmatrix} \in C \\ H^{\top}v^0 + (\theta(1 + e) - 1)v_{N,k} \ge 0 \end{array} \right. \qquad (74) $$

**Proposition 2.** Under Assumptions 1 and 2, the problem (69) has a unique solution $(v, \dot{\epsilon}, \sigma, a, v_N)$.

To prove that the min-max problem has solutions, i.e. saddle points exist, we consider a reduced version of the problem (69) where the equality constraints $Bv = \dot{\epsilon}$ and $\theta v_N = H^{\top}v - (1 - \theta)v_{N,k}$ are
substituted. This yields a reduced saddle point problem

$$
\begin{equation}
\begin{aligned}
\min_{v} \max_{\sigma, a} & \quad \frac{1}{2}(v - v_k)^{\top} M(v - v_k) - \frac{1}{2}(\sigma - \sigma_k)^{\top} S(\sigma - \sigma_k) - \frac{1}{2}(a - a_k)^{\top} D^{-1}(a - a_k) \\
& \quad + h \theta \sigma^{\top} B v - h \theta f_{\text{ext}, k+1}^{\top} v \\
\text{s.t.} & \quad \left( \begin{array}{c} \sigma \\ a \\ H^{\top} v + (\theta(1+e) - 1) v_{N,k} \end{array} \right) \in C \times \mathbb{R}_{+}^{m}.
\end{aligned}
\tag{75}
\end{equation}
$$

$$
\text{s.t.} \quad \begin{pmatrix} \sigma \\ a \\ H^\top v + (\theta(1+e) - 1)v_{N,k} \end{pmatrix} \in C \times \mathbb{R}_+^m.
$$

with its associated Lagrangian function of the form

$$
\begin{equation}
\begin{aligned}
\bar{\mathcal{L}}(v, \sigma, a) = & \frac{1}{2}(v - v_k)^{\top} M(v - v_k) - \frac{1}{2}(\sigma - \sigma_k)^{\top} S(\sigma - \sigma_k) - \frac{1}{2}(a - a_k)^{\top} D^{-1}(a - a_k) \\
& + h\theta\sigma^{\top}Bv - h\theta f_{\text{ext},k+1}^{\top}v - h\theta\Psi_C \begin{pmatrix} \sigma \\ a \end{pmatrix} + \theta^2 \Psi_{\mathbb{R}_+^m} (H^{\top}v + (\theta(1+e)-1)v_{N,k}).
\end{aligned}
\tag{76}
\end{equation}
$$

Assuming that *M*, *S* and *D* are symmetric definite positive matrices, the Lagrangian function $\bar{\mathcal{L}}$ is a convex function in *v* and a concave function of $(\sigma, a)$. With the assumption (74), the following coercivity properties hold:

$$
\lim_{\|v\| \to +\infty} \bar{\mathcal{L}}(v, \sigma^0, a^0) = +\infty, \quad \lim_{\|(\sigma, a)^\top\| \to +\infty} \bar{\mathcal{L}}(v^0, \sigma, a) = -\infty \qquad (77)
$$

From Theorem 4.3.1 in (Hiriart-Urruty and Lemaréchal, 1993), the set of solutions of the saddle point problem (75) is a nonempty compact and convex set. The strict convexity of $\bar{\mathcal{L}}$ with respect to $v$ for all $(\sigma, a) \in C$ and the strict concavity with respect to $\sigma, a$ for all $v$ such that $H^\top v^0 + (\theta(1+e)-1)v_{N,k} \ge 0$ imply the uniqueness of the solution of (75). To obtain a complete solution of (69) it suffices to compute $\dot{\epsilon} = Bv$ and $\theta v_N = H^\top v - (1-\theta)v_{N,k}$. $\square$

**Proposition 3.** *Under Assumptions 1 and 2, the problem (68) has a solution for $(v_{k+\theta}, \sigma_{k+\theta}, a_{k+\theta}, v_{N,k+1})$ and $(z_{k+\theta}, y_{k+\theta}, p_{N,k+1})$. Furthermore, the solution is unique for $(v_{k+\theta}, \sigma_{k+\theta}, a_{k+\theta}, v_{N,k+1}, z_{k+\theta}, y_{k+\theta})$.*

As it is noted in (Hiriart-Urruty and Lemaréchal, 1993, Remark 4.3.2), the coercivity properties (77)
imply that the following necessary optimality conditions has a solution

$$
\left\{
\begin{array}{l}
0 \in \partial_v \bar{\mathcal{L}}(v, \sigma, a) \\
0 \in \partial_{\sigma,a} \bar{\mathcal{L}}(v, \sigma, a)
\end{array}
\right.
\qquad (78)
$$

With our setting we get

$$
(\partial_v \bar{\mathcal{L}}:) \quad 0 \in M(v - v_k) - h\theta f_{\text{ext},k+1}^{\top} + h\theta B^{\top} \sigma + \theta^2 H \partial \Psi_{\mathbb{R}_+^m} (H^{\top} v + (\theta(1+e)-1)v_{N,k})
$$

$$
(\partial_{\sigma} \bar{\mathcal{L}}:) \quad 0 \in -S(\sigma - \sigma_k) + h\theta \dot{\epsilon} - h\theta \partial_{\sigma} \Psi_C \begin{pmatrix} \sigma \\ a \end{pmatrix}
$$

$$
(\partial_a \bar{\mathcal{L}}:) \quad 0 \in -D^{-1}(a - a_k) - h\theta \partial_a \Psi_C \begin{pmatrix} \sigma \\ a \end{pmatrix} \tag{79}
$$

Let us introduce $\theta v_N$ such that

$$
\theta v_N = H^\top v + (\theta(1+e) - 1)v_{N,k} \quad (80)
$$
and $z, y, p_N$

$$ - \begin{pmatrix} z \\ y \\ p_N \end{pmatrix} \in \partial \Psi_{C \times \mathbb{R}_+^m} \begin{pmatrix} \sigma \\ a \\ v_N \end{pmatrix} \qquad (81) $$

Then, the optimality conditions can be written as

$$ \left\{ \begin{array}{l} M(v - v_k) + h\theta B^\top \sigma = h\theta f_{\text{ext},k+1}^\top + \theta H p_N \\ S(\sigma - \sigma_k) - h\theta Bv = h\theta z \\ D^{-1}(a - a_k) = h\theta y \\ \theta v_N = H^\top v - (1 - \theta)v_{N,k} \\ - \begin{pmatrix} z \\ y \\ p_N \end{pmatrix} \in \partial \Psi_{C \times \mathbb{R}_+^m} \begin{pmatrix} \sigma \\ a \\ v_N + ev_{N,k} \end{pmatrix} \end{array} \right. \qquad (82) $$

A solution of (82) is a solution of the problem (68). The uniqueness of $(v_{k+\theta}, \sigma_{k+\theta}, a_{k+\theta}, v_{N,k+1}, z_{k+\theta}, y_{k+\theta})$ is directly obtained from the Proposition 2. □

**Remark 4.** The assumption on the positivity of the matrix $D$ implies a positive hardening that prevents from softening behavior. In the case of perfect plasticity with $D = 0$, the variable $a$ and the associated terms can be removed from the problem and we get the same results for existence and uniqueness of solutions of the discrete problem. The assumption on the velocity $v^0$ in (74) is a standard assumption for the feasibility of the problem that can be found in a more general context of contact and friction in (Acary, Cadoux, et al., 2011). Note that if $H$ has full rank that condition is satisfied.

**Reduction of the linear equations** Solving directly the problem (67) or the min-max problem (69) is possible with numerical methods of optimization, but, in the sequel, we prefer reduce the linear equations to reduce the number of unknown variables. The first three lines of (67) are linear equations and the discrete velocity $v_{k+1}$ can be substituted to obtain a reduced variational inequality. Using the two first relations in (68), we obtain

$$ h\theta z_{k+\theta} = S(\sigma_{k+\theta} - \sigma_k) - h\theta B \left( v_k + \theta M^{-1} \left( h f_{\text{ext},k+\theta} + H p_{N,k+1} - h B^\top \sigma_{k+\theta} \right) \right), \qquad (83) $$

that can be simplified in

$$ h\theta z_{k+\theta} = \underbrace{[S + h^2 \theta^2 B M^{-1} B^\top]}_{U} \sigma_{k+\theta} - \underbrace{[h\theta^2 B M^{-1} H]}_{V} p_{N,k+1} + s, \qquad (84) $$

with $s = -S\sigma_k - h\theta B(v_k + \theta h M^{-1} f_{\text{ext},k+\theta})$. A similar expression for $v_{N,k+1} + ev_{N,k}$ is obtained by multiplying by $H^\top$:

$$ v_{N,k+1} + ev_{N,k} = ev_{N,k} + H^\top (v_k + h\theta M^{-1} f_{\text{ext},k+\theta}) + H^\top M^{-1} H p_{N,k+1} - hH^\top M^{-1} B^\top \sigma_{k+\theta}. \qquad (85) $$

In order to get a symmetric problem, we multiply the previous equation by $\theta^2$

$$ \theta^2 (v_{N,k+1} + ev_{N,k}) = \underbrace{\theta^2 H^\top M^{-1} H}_{W} p_{N,k+1} - \underbrace{h\theta^2 H^\top M^{-1} B^\top}_{V^\top} \sigma_{k+\theta} + r \qquad (86) $$
with $r = \theta^2 (ev_{N,k} + H(v_k + \theta M^{-1} (hf_{\text{ext},k+\theta})))$. The last line of the inclusion is finally rewritten as

$$h\theta y_{k+\theta} = D^{-1}(a_{k+\theta} - a_k). \qquad (87)$$

Since the inclusion involves a cone, the problem (67) amounts to solving the following affine variational inequality

$$- \left( Q \begin{pmatrix} \sigma_{k+\theta} \\ a_{k+\theta} \\ p_{N,k+1} \end{pmatrix} + p \right) \in N_{C \times \mathbb{R}_+^m} \begin{pmatrix} \sigma_{k+\theta} \\ a_{k+\theta} \\ p_{N,k+1} \end{pmatrix}, \qquad (88)$$

with

$$Q = \begin{pmatrix} U & 0 & -V \\ 0 & D^{-1} & 0 \\ -V^{\top} & 0 & W \end{pmatrix} \quad \text{and} \quad p = \begin{pmatrix} s \\ D^{-1}a_k \\ r \end{pmatrix}. \qquad (89)$$

**Equivalent convex quadratic optimization problem** The goal of the section is to show that the affine variational inequality (88) is equivalent to a convex quadratic problem if Assumptions 1 hold. The matrix $Q$ can be written as

$$\begin{align}
Q &= \begin{pmatrix} S & 0 & 0 \\ 0 & D^{-1} & 0 \\ 0 & 0 & 0 \end{pmatrix} + \begin{pmatrix} h^2 \theta^2 B M^{-1} B^{\top} & 0 & -h \theta^2 B M^{-1} H \\ 0 & 0 & 0 \\ -h \theta^2 H^{\top} M^{-1} B^{\top} & 0 & \theta^2 H^{\top} M^{-1} H \end{pmatrix} \tag{90} \\
&= \begin{pmatrix} S & 0 & 0 \\ 0 & D^{-1} & 0 \\ 0 & 0 & 0 \end{pmatrix} + \begin{pmatrix} -h \theta B^{\top} & 0 & \theta H \end{pmatrix}^{\top} M^{-1} \begin{pmatrix} -h \theta B^{\top} & 0 & \theta H \end{pmatrix}. \nonumber
\end{align}$$

In this form, the matrix $Q$ appears to be a symmetric semi-definite matrix as a sum of two symmetric semi-definite matrices. With this property, solving the affine variational inequality (88) is equivalent to solve the following convex quadratic optimization problem

$$\begin{align}
\min_{\sigma, a, p_N} & \quad \frac{1}{2} \begin{pmatrix} \sigma \\ a \\ p_N \end{pmatrix}^{\top} Q \begin{pmatrix} \sigma \\ a \\ p_N \end{pmatrix} + p^{\top} \begin{pmatrix} \sigma \\ a \\ p_N \end{pmatrix} \tag{91} \\
\text{s.t.} & \quad \begin{pmatrix} \sigma \\ a \\ p_N \end{pmatrix} \in C \times \mathbb{R}_+^m. \nonumber
\end{align}$$

The well-posedness of the problem (91) is then equivalent of the well-posedness of the optimization problem. We already ensure that the solution exists and is unique except for $p_N$. In the convex case without strict convexity, results can be quite technical and needs assumption on the matrix $H$ such as Slater condition. In the following paragraph, we prove it in the simpler case of strict convexity of the cost function adding a rank condition on $H$.

**Existence and uniqueness in the strictly convex case.** To prove the existence and uniqueness of solutions, let us define the following assumption.

**Assumption 3.** The matrix $H$ has full rank.

**Proposition 4.** Under Assumptions 1 and 3 and for a sufficiently small time step, the problem (91) has a unique solution $(\sigma, a, p_N)$ if the set $C$ is a non empty convex set.
Under Assumptions 1 and 3, the matrix $W = \theta^2 H^\top M^{-1} H$ is a symmetric positive definite matrix. The matrix $Q$ can be written as

$$Q = \begin{pmatrix} S & 0 & 0 \\ 0 & D^{-1} & 0 \\ 0 & 0 & W \end{pmatrix} + h \begin{pmatrix} h\theta^2 B M^{-1} B^\top & 0 & -\theta^2 B M^{-1} H \\ 0 & 0 & 0 \\ -\theta^2 H^\top M^{-1} B^\top & 0 & 0 \end{pmatrix} \quad (92)$$

For a sufficiently small time step $h$, the matrix $Q$ is then a positive definite matrix. Since the optimization problem (91) is strictly convex, a unique solution exists if the set $C$ is a non-empty convex set. $\square$

**Remark 5.** This assumption 3 deserves some comments. The full rank property of the matrix $H$ is standard in finite element application where the constraints are chosen to be linearly independent.

## 3.4 Numerical methods

Since the problem (91) appears to be a convex quadratic optimization problem, numerous numerical methods can be applied to solve it. In practice, the set $C$ is generally finitely represented and given by a set of inequalities that describes the yield function:

$$C = \{(\sigma, a) \mid f(\sigma, a) \ge 0\}, \quad (93)$$

where $f$ is a smooth vector-valued function with non-vanishing gradients at the solution. The optimization problem (91) takes the form

$$\min_{\sigma, a, p_N} \quad \frac{1}{2} \begin{pmatrix} \sigma \\ p_N \\ a \end{pmatrix}^\top Q \begin{pmatrix} \sigma \\ p_N \\ a \end{pmatrix} + p^\top \begin{pmatrix} \sigma \\ p_N \\ a \end{pmatrix} \quad (94)$$

$$\text{s.t.} \quad \begin{array}{l} f(\sigma, a) \ge 0 \\ p_N \ge 0 \end{array}$$

that can be solved by any kind of numerical optimization techniques for solving quadratic optimization (Bonnans et al., 2006), interior point methods (Wright, 1996), semi-smooth Newton methods, proximal algorithms (Parikh, 2014) such as Alternating Direction Methods of Multipliers (ADMM). In the context of quasistatic plasticity, some of these methods have already been validated, see for instance (Kanno, 2020; Krabbenhøft, Lyamin, Sloan, and Wriggers, 2007; Krabbenhøft, Lyamin, and Sloan, 2007; Shimizu and Kanno, 2018) and also in the coupled problem of contact and plasticity (Christensen, 2002a).

**The case of $J_2$ perfect plasticity.** The $J_2$ plasticity is based on the Huber-Von Mises criterion amounts to defining the yield function $f$ with the $J_2$ invariant of the stress $\sigma$.

$$J_2 = \frac{1}{2} s : s = \frac{1}{2} s_{ij} s_{ji} = \frac{1}{2} \operatorname{tr}(s \cdot s) \quad (95)$$

where $s$ is the deviatoric part of the tensor $\sigma$, and the hydrostatic stress $\sigma^h$ defined by

$$s = \sigma - \sigma^h \mathbf{I} = \sigma - \frac{1}{3} \operatorname{tr}(\sigma) \mathbf{I} = \sigma - \frac{1}{3} I_1 \mathbf{I} \quad \text{with } \sigma^h = \frac{1}{3} I_1, I_1 = \sigma : \mathbf{I} = \sigma_{ii} = \operatorname{tr}(\sigma). \quad (96)$$

The Huber-Von Mises criterion is

$$\sqrt{3J_2} \le \kappa, \text{ or } 3J_2 \le \kappa^2. \quad (97)$$

With an alternative expression of the $J_2$ invariant, another expression of the Huber-Von Mises criterion is in the original stress tensor:

$$\frac{1}{2} \left( (\sigma_{11} - \sigma_{22})^2 + (\sigma_{22} - \sigma_{33})^2 + (\sigma_{33} - \sigma_{11})^2 + 6(\sigma_{12}^2 + \sigma_{13}^2 + \sigma_{23}^2) \right) \le \kappa^2. \quad (98)$$
Using the Voigt notation for the stress $\sigma = (\sigma_{11} \ \sigma_{22} \ \sigma_{33} \ \sigma_{23} \ \sigma_{13} \ \sigma_{12})^\top$, the $J_2$ invariant can be expressed as a quadratic function

$$s^\top s = \sigma^\top J \sigma, \qquad (99)$$

with

$$J = \frac{1}{3} \begin{pmatrix} 2 & -1 & -1 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 & 0 \\ -1 & -1 & 2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 6 & 0 & 0 \\ 0 & 0 & 0 & 0 & 6 & 0 \\ 0 & 0 & 0 & 0 & 0 & 6 \end{pmatrix}. \qquad (100)$$

The matrix $J$ is semi-definite positive which shows that $C$ is convex. In that case, we have to solve a quadratic program with quadratic constraints that can be solved efficiently with SQP solvers

$$\begin{array}{ll} \min_{\sigma,a,p_N} & \frac{1}{2} \begin{pmatrix} \sigma \\ p_N \\ a \end{pmatrix}^\top Q \begin{pmatrix} \sigma \\ p_N \\ a \end{pmatrix} + p^\top \begin{pmatrix} \sigma \\ p_N \\ a \end{pmatrix} \\ \text{s.t.} & 3\sigma_g^\top J \sigma_g \ge \kappa^2 \text{ for } g \in [1, n_{e,g}] \\ & p_N \ge 0. \end{array} \qquad (101)$$

**The case of Drucker-Prager plasticity.** The Drucker-Prager plasticity involves the following definition of the convex set $C$ as

$$C := K_{\sigma} = \left\{ (\sigma^h, s) \mid \frac{1}{k_d} \|s\| + \sigma^h \tan(\varphi) \le c \right\} \qquad (102)$$

where $\|s\| = \sqrt{J_2}$, $k_d$ is a constant that depends on the dimension $d$, $c$ is the cohesion and $\varphi$ the friction angle. It is well known that the set $K_\sigma$ is a translated second order cone that cannot be straightforwardly written with a smooth function $f$ as in (93). The fact that $C$ has the structure of a pointed convex cone renders difficult to express the set $C$ with a smooth vector-valued function, but it can be formulated on symmetric cone of semi-definite type or second order types (Berga and De Saxcé, 1994; Bisbos et al., 2005; Hjiaj et al., 2003). Then the tools from second-order cone or semi-definite programming can be used to solve

$$\begin{array}{ll} \min_{\sigma,a,p_N} & \frac{1}{2} \begin{pmatrix} \sigma \\ p_N \\ a \end{pmatrix}^\top Q \begin{pmatrix} \sigma \\ p_N \\ a \end{pmatrix} + p^\top \begin{pmatrix} \sigma \\ p_N \\ a \end{pmatrix} \\ \text{s.t.} & \sigma^h = \frac{1}{3} \operatorname{tr}(\sigma), s = \sigma - \sigma^h I \\ & (\sigma^h, s)_g \in K_\sigma \text{ for } g \in [1, n_{e,g}] \\ & p_N \ge 0. \end{array} \qquad (103)$$

**Remark 6.** It is important to note that the solution of the convex optimisation problem leads to a model of associated plasticity, which is unusual for geomaterial plasticity related to the Drucker-Prager criteria. Indeed, the optimality conditions related to the constraints $(\sigma^h, s)_g \in K_\sigma$ defines the following flow rule

$$\left( \begin{array}{c} \dot{\varepsilon}^{p,h} \\ \dot{e}^p \end{array} \right) \in N_{K_\sigma} \left( \begin{array}{c} \sigma^h \\ s \end{array} \right), \qquad (104)$$

where $e = \varepsilon - \frac{1}{3}\varepsilon^h I$ and $\varepsilon^h = \operatorname{tr}\varepsilon$. As it is introduced and explained in Berga and De Saxcé, 1994, the nonassociated plasticity needs the introduction of Implicit Standard Materials (ISM) yielding the flow rule

$$\left( \begin{array}{c} \dot{\varepsilon}^{p,h} + k_d (\tan(\varphi) - \tan(\theta)) \|\dot{e}^p\| \\ \dot{e}^p \end{array} \right) \in N_{K_\sigma} \left( \begin{array}{c} \sigma^h \\ s \end{array} \right), \qquad (105)$$
where $\theta$ is the dilatancy angle ranging from 0 to $\phi$. In the latter case, there is no direct associated convex optimization problem. Hence, the complementarity problem over cones must be directly solved.

**Nonsmooth Newton algorithm, projection formulation and return mapping.** Solving the optimization problems (91) or (94) amounts to finding a solution of the optimality conditions given by the variational inequality (88). For a standard variational inequality on a convex set $X$, it is well known that it is possible to reformulate it with the natural map operator (see Facchinei and Pang (2003) for details):

$$ -F(x) \in N_X(x) \iff R(x) := x - \operatorname{proj}_X(x - \rho F(x)) = 0, \rho > 0. \quad (106) $$

In the case of (88), we obtain the following equation for $\rho > 0$

$$ R \left( \begin{pmatrix} \sigma \\ a \\ p_N \end{pmatrix} \right) = \begin{pmatrix} \sigma \\ a \\ p_N \end{pmatrix} - \operatorname{proj}_{C \times \mathbb{R}_+^m} \left( \begin{pmatrix} \sigma \\ a \\ p_N \end{pmatrix} - \rho \left( Q \begin{pmatrix} \sigma \\ a \\ p_N \end{pmatrix} + p \right) \right) = 0. \quad (107) $$

In practical plasticity problems, the projection on $C$ is known. Then, this equation can be solved by fixed point techniques with projection, also called Uzawa method in computational mechanics (Berga and De Saçcé, 1994). The method is robust but rather slow, and a second order method is preferred. To this aim, the equation $R(x) = 0$ in (106) is solved with a semismooth Newton method as

$$ x_{k+1} = x_k + H_k^{-1}(-R(x_k)), \quad \text{with } H_k \in \partial R(x_k). \quad (108) $$

When the formulation is based on the projection operator, it is rather easy to compute $H_k$, an element of the generalized Jacobian of $R$. The semi-smooth Newton method contains the return-mapping algorithm as a special case as it has been shown in Christensen (2002a). If the problem is equivalent to an optimization problem, then the method can also be globalized.

## 3.5 Structure preserving properties

In this section, we are interested in the structure preserving properties of the proposed time-stepping scheme, especially, the discrete energy balance of the scheme when plastic flow and unilateral contact occur. In the continuous time and space setting the energy balance reads as:

$$ \frac{d}{dt} \mathcal{E}(t) = \frac{d}{dt} (\mathcal{T}(t) + \mathcal{P}(t)) = -\mathcal{D}(t) + \mathcal{P}_v(t) + \mathcal{P}_{\text{ext}}(t), \quad (109) $$

where $\mathcal{E}(t)$ is total mechanical energy $\mathcal{E}(t) = \mathcal{T}(t) + \mathcal{P}(t)$. After the space discretization by a finite element method, the energy balance takes the form

$$ d(\mathcal{T}(t) + \mathcal{P}(t)) = -\mathcal{D}(t)dt + \mathcal{P}_{\text{ext}}(t)dt + d\mathcal{P}_{\text{impact}}, \quad (110) $$

where the space-discretized energies are given by

$$ \mathcal{E}(t) = \mathcal{T}(t) + \mathcal{P}(t) \quad \text{with} \quad \mathcal{T}(t) = \frac{1}{2} v^{\top}(t) M v(t), \quad \text{and} \quad \mathcal{P}(t) = \frac{1}{2} \varepsilon^{\top}(t) E \varepsilon^{\top}(t) + \frac{1}{2} \alpha^{\top}(t) D \alpha(t) \quad (111) $$

with the convention that the kinetic energy is also a right continuous function of bounded variations, i.e, $\mathcal{T}(t) = \mathcal{T}^{+}(t)$. The dissipation and the power of external forces are

$$ \mathcal{D}(t) = \sigma^{\top}(t) \dot{\varepsilon}^{p}(t) + a^{\top}(t) \dot{\alpha}(t) \quad \text{and} \quad \mathcal{P}_{\text{ext}}(t) = f_{\text{ext}}^{\top}(t) v(t). \quad (112) $$

Since in a discrete system, impacts may occur, the power of the reaction impulse is given by

$$ d\mathcal{P}_{\text{impact}} = \frac{1}{2}(v_N^+ + v_N^-)d\dot{i}_N. \quad (113) $$
Let us now compute the variation of the total mechanical energy over a time-step, starting by the kinetic energy

$$ \Delta T_k^{k+1} = T(t_{k+1}) - T(t_k) = \frac{1}{2}(v_{k+1} + v_k)^\top M(v_{k+1} - v_k). \quad (114) $$

Using the relation $\frac{1}{2}(v_{k+1} + v_k) = v_{k+\theta} + (\frac{1}{2} - \theta)(v_{k+1} - v_k)$ and the first equation in (67), we obtain

$$ \begin{aligned} \Delta T_k^{k+1} &= (\frac{1}{2} - \theta)\|v_{k+1} - v_k\|_M^2 + v_{k+\theta}^\top (h f_{\text{ext},k+\theta} - h B^\top \sigma_{k+\theta} + H p_{N,k+1}) \\ &= (\frac{1}{2} - \theta)\|v_{k+1} - v_k\|_M^2 + h v_{k+\theta}^\top f_{\text{ext},k+\theta} - h \sigma_{k+\theta}^\top \dot{\varepsilon}_{k+\theta} + v_{N,k+\theta}^\top p_{N,k+1} \end{aligned} \quad (115) $$

with the notation $\|x\|_M^2 = x^\top M x$ for a symmetric definite positive matrix $M$. For the potential energy, we obtain

$$ \begin{aligned} \Delta \Psi_k^{k+1} &= \Psi(t_{k+1}) - \Psi(t_k) = \frac{1}{2}(\varepsilon_{k+1}^e + \varepsilon_k^e)^\top E (\varepsilon_{k+1}^e - \varepsilon_k^e) + \frac{1}{2}(\alpha_{k+1} + \alpha_k)^\top D (\alpha_{k+1} - \alpha_k) \\ &= (\frac{1}{2} - \theta)\|\varepsilon_{k+1}^e - \varepsilon_k^e\|_E^2 + h \sigma_{k+\theta}^\top \dot{\varepsilon}_{k+\theta} + (\frac{1}{2} - \theta)\|\alpha_{k+1} - \alpha_k\|_D^2 + h y_{k+\theta}^\top \alpha_{k+\theta}, \end{aligned} \quad (116) $$

with the notation $\|x\|_E^2 = x^\top E x$ and $\|x\|_D^2 = x^\top D x$ for symmetric definite positive matrices $E$ and $D$. Using the approximation of works of external forces and the dissipated work by the $\theta$-method as

$$ W_{\text{ext } k}^{k+1} := h v_{k+\theta}^{\top} f_{\text{ext},k+\theta} \approx \int_{t_k}^{t_{k+1}} P_{\text{ext}}(t) dt, \quad (117) $$

and

$$ W_{p \ k}^{k+1} := h \sigma_{k+\theta}^{\top} \dot{\varepsilon}_{k+\theta}^{p} - h a_{k+\theta}^{\top} y_{k+\theta} \approx \int_{t_k}^{t_{k+1}} D(t) dt, \quad (118) $$

and finally, an approximation of the work dissipated by the percussion (see Acary (2016))

$$ W_{c \ k}^{k+1} := v_{N,k+\theta}^{\top} p_{N,k+1} = (1 - \theta(1 + e)) v_{N,k}^{\top} p_{N,k+1}, \quad (119) $$

the increment of total energy is then given by

$$ \Delta E_k^{k+1} - W_{\text{ext } k}^{k+1} + W_{p \ k}^{k+1} - W_{c \ k}^{k+1} = (\frac{1}{2} - \theta) (\|v_{k+1} - v_k\|_M^2 + \|\varepsilon_{k+1}^e - \varepsilon_k^e\|_E^2 + \|\alpha_{k+1} - \alpha_k\|_D^2). \quad (120) $$

The following proposition summarizes the properties of the discrete energy balance (120).

**Proposition 5.** Under Assumption 1, the main conclusions that can be drawn for the energy dissipation of the scheme are

1. When $\theta = \frac{1}{2}$, the time-stepping scheme satisfies the approximation of the discrete energy balance without introducing artificial dissipation:

$$ \Delta E_k^{k+1} - W_{\text{ext } k}^{k+1} = -W_{p \ k}^{k+1} + W_{c \ k}^{k+1}. \quad (121) $$

2. The dissipated work due to plasticity is positive and for $\theta \le \frac{1}{1+e}$, the dissipated work due to impact is also positive.
3. When $\frac{1}{2} \le \theta \le \frac{1}{1+e} \le 1$, we have the following dissipation inequality

$$ \Delta E_k^{k+1} - W_{\text{ext}} k^{k+1} \le 0. \qquad (122) $$

A result in (Acary, 2016) shows that

$$ W_{c k}^{k+1} \le 0 \text{ for } \theta \le \frac{1}{1+e}. \qquad (123) $$

Furthermore, the inclusion

$$ - \begin{pmatrix} z_{k+\theta} \\ y_{k+\theta} \end{pmatrix} \in N_{C \times \mathbb{R}_+^m} \begin{pmatrix} \sigma_{k+\theta} \\ a_{k+\theta} \end{pmatrix}, \qquad (124) $$

implies that $W_p \ge 0$.

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

The inequality (122) ensures the practical unconditional stability of the scheme.

**Remark 7.** Naturally, the $\theta$-method has some shortcomings such as excessive numerical dissipation in dynamics for $\theta > 1/2$ which makes the scheme more suitable for quasi-static evolutions. For $\theta = 1/2$, the scheme has very good energy conservation properties that do not filter out spurious oscillations due to the discretization in space. In this case, it is clear that a non-smooth Newmark or non-smooth generalized-$\alpha$ scheme (Brüls et al., 2014, 2018; Chen, Acary, et al., 2013) is much better suited. In the article, we have decided to keep $\theta$-method for simplicity, but nothing prevents the use of more advanced schemes.

# 4 Numerical applications

In this section, a series of examples demonstrates the capabilities of the proposed formulation. Even if the approach can be applied to complex structures and nonlinear yield surfaces, the choice was made to propose simple examples for illustration purposes. We focused on beam structures able to sustain both tensile and bending forces. Assuming that plastic strains occur from the combined effects of axial and bending forces, we considered a piecewise linear yield surface expressed in terms of tensile and bending stresses. Associated flow rule as well as isotropic and kinematic hardening were considered. Both statics and dynamics responses of the structures were analyzed. Dedicated Python scripts were developed using the LCPs solvers implemented in Siconos code (Acary, Bremond, et al., 2023).

## 4.1 Finite element model of elastoplastic beam structures

### 4.1.1 Finite element method for Elasto-plastic Euler-Bernoulli beams with piecewise linear yield functions - a simple illustrative case

<figure><img src="image_6.png" /><figcaption>Figure 2: Description of the beam element used</figcaption></figure>
In the examples, specific beam elements able to sustain bending and tensile loadings are used. Beam elements based on Euler-Bernoulli beam theory are classically formulated using third order shape functions and numerical integration techniques based on Gaussian quadrature that involve two Gauss points per element. In order to consider two Gauss points for the axial displacements, we introduced elements with three nodes and quadratic shape functions in the axial direction. The nodal displacements and Gauss point locations of the element used are presented in Figure 2. For each element, the generalized stresses are expressed at Gauss points as :

$$ \sigma_e = [n_1 \ m_1 \ n_2 \ m_2]^\top, \qquad (125) $$

where $n_g$ is the normal force at Gauss point g and $m_g$ is the bending moment at Gauss point g. Standard elementary matrices are recalled in Appendix C.

### 4.1.2 Piecewise linear yield surface

**The case of a feasible set C defined by linear relations** When the set C is finite represented as in (93) and assuming that f is sufficiently smooth with nonvanishing gradients, the normal cone to C can expressed as

$$ N_C \begin{pmatrix} \sigma \\ a \end{pmatrix} = \left\{ (z, y) \mid \begin{pmatrix} z \\ y \end{pmatrix} = \nabla_{\sigma,a}^{\top} f(\sigma, a)\lambda, \lambda \ge 0, \lambda^{\top} f(\sigma, a) = 0 \right\}. \qquad (126) $$

Then the plastic flow rule

$$ \left( \begin{array}{c} \dot{\varepsilon}^p \\ \dot{\alpha} \end{array} \right) \in N_C \left( \begin{array}{c} \sigma \\ a \end{array} \right) \qquad (127) $$

is expressed as the classical flow rule for plasticity with the complementarity condition

$$ \left\{ \begin{array}{l} \dot{\varepsilon}^p = -\nabla_{\sigma}^{\top} f(\sigma, a)\lambda \\ \dot{\alpha} = -\nabla_{a}^{\top} f(\sigma, a)\lambda \\ 0 \le \lambda \perp f(\sigma, a) \ge 0, \end{array} \right. \qquad (128) $$

where $\lambda$ is the vector of plastic multipliers. If the yield function f is linear, the gradients are given by constant matrices denoted by

$$ J_{\sigma} = \nabla_{\sigma}^{\top} f(\sigma, a) \quad \text{and} \quad J_{a} = \nabla_{a}^{\top} f(\sigma, a), \qquad (129) $$

and we have

$$ f(\sigma, a) = J_{\sigma}^{\top} \sigma + J_{a}^{\top} a + \kappa_{Y}. \qquad (130) $$

With linear yield functions, the variational inequality (88) can be reduced to a Linear Complementarity Problem (LCP). To this aim, we express the stress tensor from (84) as

$$ \sigma_{k+\theta} = h\theta U^{-1} J_{\sigma} \lambda - U^{-1}s + U^{-1}V p_{N,k+1}. \qquad (131) $$

Substituting (131) in (86) and (130) yields

$$ \theta^2(v_{N,k+1} + ev_{N,k}) = (W - V^\top U^{-1}V)p_{N,k+1} - h\theta V^\top U^{-1}J_\sigma\lambda + V^\top U^{-1}s + r, \qquad (132) $$

and

$$ \begin{aligned} f_{k+\theta} &= J_{\sigma}^{\top} \sigma_{k+\theta} + J_{a}^{\top} a_{k+\theta} + \kappa_{Y} \\ &= h\theta (J_{\sigma}^{\top} U^{-1} J_{\sigma} + J_{a}^{\top} D J_{a}) \lambda + J_{\sigma}^{\top} U^{-1} V p_{N,k+1} - J_{\sigma}^{\top} U^{-1} s + J_{a}^{\top} a_{k} + \kappa_{Y}. \end{aligned} \qquad (133) $$
The problem in terms of $\lambda$ and $p_{N,k+1}$ can be defined in the form of the following LCP

$$ \left\{ \begin{array}{l} \left( \begin{array}{c} f_{k+\theta} \\ v_{N,k+1} + e v_{N,k} \end{array} \right) = L \left( \begin{array}{c} h\theta\lambda \\ p_{N,k+1} \end{array} \right) + q \\ 0 \le h\theta\lambda \perp f_{k+\theta} \ge 0 \\ 0 \le v_{N,k+1} + e v_{N,k} \perp p_{N,k+1}, \end{array} \right. \qquad (134) $$

with

$$ L = \begin{pmatrix} J_{\sigma}^{\top} U^{-1} J_{\sigma} + J_{a}^{\top} D J_{a} & J_{\sigma}^{\top} U^{-1} V \\ -V^{\top} U^{-1} J_{\sigma} & W - V^{\top} U^{-1} V \end{pmatrix} \quad \text{and} \quad q = \begin{pmatrix} -J_{\sigma}^{\top} U^{-1} s + J_{a}^{\top} a_{k} + \kappa_{Y} \\ V^{\top} U^{-1} s + r \end{pmatrix}. \quad (135) $$

In the case of quasistatic plasticity without contact conditions, a very similar formulation of the problem as a LCP can be found in (Maier, 1984).

<figure><img src="image_7.png" /><figcaption>Figure 3: Yield surface in the case of perfect plasticity (a), kinematic hardening (b), and isotropic hardening (c) in the space of generalized stresses $(n, m)$.</figcaption></figure>

The chosen plastic yield surface takes into account the combined effects of axial and bending forces when a plastic flow occurs. A yield criterion based on the approach proposed in Heng, Hjiaj, et al. (2016) and Duan and Chen (1990) is chosen, in the specific case of linear yield functions (see Figure 3). Considering a set of linear yield criteria at a given Gauss point $g$, the matrices $J_{\sigma}^{\top}$, $J_{a}^{\top}$, $\kappa_{Y}$ are defined in the present case by

$$ J_{\sigma}^{\top} = \begin{bmatrix} \frac{1}{n_p} & \frac{1}{m_p} \\ \frac{1}{n_p} & \frac{1}{m_p} \\ \frac{1}{n_p} & \frac{1}{m_p} \\ \frac{1}{n_p} & \frac{1}{m_p} \end{bmatrix}, \quad J_{a}^{\top} = \begin{bmatrix} \frac{1}{m_p} & \frac{1}{m_p} & \frac{1}{m_p} \\ \frac{1}{m_p} & \frac{1}{m_p} & \frac{1}{m_p} \\ \frac{1}{m_p} & \frac{1}{m_p} & \frac{1}{m_p} \\ \frac{1}{m_p} & \frac{1}{m_p} & \frac{1}{m_p} \end{bmatrix}, \quad \kappa_{Y} = \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix}, \qquad (136) $$

where $n_p$ is the limit admissible tensile force and $m_p$ is the limit admissible bending moment. At Gauss point $g$, the generalized stresses are noted by:

$$ \sigma_{e,g} = [n_g \quad m_g]^{\top}, \qquad (137) $$

Both kinematic and isotropic hardenings are considered leading to the definition of the force vector $a_g$, expressed at each Gauss point as

$$ a_g = [a_{ki,n,g} \quad a_{ki,m,g} \quad a_{is,g}]^{\top}, \qquad (138) $$
where $a_{ki,n,g}$ and $a_{ki,m,g}$ are associated with kinematic hardening in terms of normal force $n_g$ and bending moment $m_g$, respectively and $a_{is,g} \ge 0$ is associated with isotropic hardening. The vector of plastic multipliers $\lambda_g$ at Gauss point $g$ is denoted by

$$ \lambda_g = [\lambda_{g,1} \ \lambda_{g,2} \ \lambda_{g,3} \ \lambda_{g,4}]^\top. \qquad (139) $$

## 4.2 Example 1: Cantilever beam subjected to quasi-static displacements

<figure><img src="image_8.png" /><figcaption>Figure 4: Example 1: Cantilever beam discretized in 4 elements subjected to quasi-static displacements $\bar{u}_x$ and $\bar{u}_y$ at its free end.</figcaption></figure>

The first example is dedicated to evaluate the capability of the approach to simulate the response of elastoplastic structures subjected to various quasi-static loading scenarios. This example also illustrates the evolution of the different variables, in particular of the plastic multipliers, during the onset of plastic strains. The structure is an elastoplastic cantilever beam of length $L$, discretized in 4 equal length beam elements to keep it as simple as possible, and subjected to quasi-static displacements $\bar{u}_x$ and $\bar{u}_y$ at its free end. The properties of the structure are summarized in Table 1.

<table><caption>Table 1: Properties of the elastoplastic cantilever beam used in example 1.</caption><tr><td>Length</td><td>L (m)</td><td>0.8</td></tr><tr><td>Moment of inertia</td><td>I (m<sup>4</sup>)</td><td>0.0283 &times; 10<sup>-6</sup></td></tr><tr><td>Cross section area</td><td>A (m<sup>2</sup>)</td><td>0.22 &times; 10<sup>-3</sup></td></tr><tr><td>Young modulus</td><td>E (MPa)</td><td>210 000</td></tr><tr><td>Kinematic hardening modulus</td><td>H<sub>ki</sub> (MPa)</td><td>E/10 = 21 000 &times; 10<sup>6</sup></td></tr><tr><td>Isotropic hardening modulus</td><td>H<sub>is</sub> (MPa)</td><td>E/100 = 2100 &times; 10<sup>6</sup></td></tr><tr><td>Limit admissible tensile force</td><td>n<sub>P</sub> (N)</td><td>33 000</td></tr><tr><td>Limit admissible bending moment</td><td>m<sub>P</sub> (Nm)</td><td>283.5</td></tr><tr><td>Material density</td><td>&rho; (kg/m<sup>3</sup>)</td><td>7701</td></tr></table>

As the formulation explicitly integrates dynamical effects, we can directly model the quasi-static loadings as low velocity dynamical loadings. In these simulations, the time step was set at $h = 10^{-2}$s. In addition, to ensure a maximal dissipation of the kinetic energy, we used $\theta = 1$. Three cases are considered with different combinations of applied velocities and plastic behavior in order to reach different regimes for the plastic flows (see Figure 5 for a definition of the applied velocities).

* **Case 1** We consider a plastic law without hardening and a constant velocity is applied at the beam free end. The axial and transversal components of the velocity are $1 \times 10^{-4}$m/s and $4 \times 10^{-3}$m/s, respectively (see Figure 5).

* **Case 2** We consider a plastic law with isotropic and kinematic hardening with the same prescribed velocities.
* **Case 3** In this case, we consider a plastic law with isotropic and kinematic hardening and we apply first a constant transversal velocity, set at $8 \times 10^{-3}$ m/s until $t = 2.5$ s. Then, the transversal velocity is cancelled and a constant axial velocity, equal to $2 \times 10^{-4}$ m/s, is applied (see Figure 5).

<figure><img src="image_9.png" /><figcaption>Figure 5: Example 1: Time evolution of the applied velocities at the free extremity for the different loading cases.</figcaption></figure>

**Case 1: Transversal and axial displacement - perfect plasticity** As shown in Figure 6, the axial force $n_g$ and the bending moment $m_g$ first linearly increase until the yield criterion is reached for Gauss point 1 at time 1.84s. Once the yield criterion is reached, the plastic multiplier $\lambda_{1,3}$ becomes positive (Figure 8). Hence, the generalized stresses at Gauss point 1 evolve along the yield surface for increasing displacements (Figure 7). The onset of plastic strains at Gauss point 1 entails changes in the evolution of the stresses at all Gauss points. In particular, the bending moments decrease due to the formation of a plastic hinge at Gauss point 1.

<figure><img src="image_10.png" /><figcaption>Figure 6: Example 1, case1: time evolutions of the axial forces $n_g$ (a) and bending moments $m_g$ (b) at the Gauss points $g$ of a cantilever beam subjected to axial and transversal displacements simultaneously in the case of perfect plasticity.</figcaption></figure>

**Case 2: Transversal and axial displacement - isotropic and kinematic hardening** The introduction of hardening entails both translation of the yield surface and increases in this surface once the yield
<figure><img src="image_11.png" /><figcaption>Figure 7: Example 1, case 1: stress path at the Gauss points of the cantilever beam subjected to axial and transversal displacements simultaneously in the case of perfect plasticity. The bending moments $m_g$ at the Gauss points $g$ are plotted against the axial force $n_g$. The diamond shape is the yield surface $C$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 8: Example 1, case 1: time evolution of the plastic multipliers $\lambda_{1,3}$ at Gauss point 1 of the cantilever beam subjected to axial and transversal displacements simultaneously in the case of perfect plasticity.</figcaption></figure>
criterion is reached (Figure 10). Contrary to the case of perfect plasticity, the bending moments do not decrease when the yield criterion is reached (Figure 9). In this case, we do not observe a plastic hinge located only at Gauss point 1. The yield criterion is successively reached at Gauss point 2, and then at Gauss point 3, for increasing displacements. The plastic multipliers successively become positive (Figure 11) when the yield criterion is reached at a given Gauss point.

<figure><img src="image_13.png" /><figcaption>Figure 9: Example 1, case 2: Time evolution of the axial forces $n_g$ (a) and bending moments $m_g$ (b) at the Gauss points $g$ of the cantilever beam subjected to axial and transversal displacements simultaneously in the case of isotropic and kinematic hardenings.</figcaption></figure>

**Case 3: Transversal followed by axial displacement - perfect plasticity** The third case illustrates the specific situation when the stress at a Gauss point reaches a corner of the admissible set, i.e. an intersection of two linear limits of the yield surface. The node displacements are depicted in Figure 12. When the yield criterion is reached at Gauss point 1 for $t = 1.4s$ (Figure 14a), the plastic multipliers $\lambda_{1,3}$ and $\lambda_{1,4}$ both become positive (Figure 15). Indeed, at the corner, the yield criterion fullfillment requires the activation of two conditions. The onset of plastic strains results in the formation of a plastic hinge at Gauss point 1. As classically for perfectly plastic cantilever beams subjected to bending loading only, the bending moments remain constant when the yield criterion is reached (Figure 13). When the transversal velocity is cancelled and the constant axial velocity is applied ($t = 2.5s$ - Figure 12), the loading path is modified. Thus, only one inequality constraint is active and the multiplier $\lambda_{1,4}$ vanishes (Figure 15). To fullfill the yield criterion under this loading path, the bending moments decrease while the axial forces increase (Figure 13).
<figure><img src="image_14.png" /><figcaption>Figure 10: Example 1, case 2: stress path at the Gauss points of the cantilever beam subjected to axial and transversal displacements simultaneously in the case of isotropic and kinematic hardenings. The bending moments $m_g$ at the Gauss points $g$ are plotted against the axial force $n_g$. The grey diamond shape is the initial yield surface and the black diamond shape is the final one.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 11: Example 1, case 2: time evolutions of the plastic multipliers $\lambda_{g,3}$ at the first 4 Gauss points $g$ of the cantilever beam subjected to axial and transversal displacements simultaneously in the case of isotropic and kinematic hardening.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 12: Example 1, case 3: time evolution of the axial displacements $u_x$ and transversal displacements $u_y$ at the 4 free nodes of the cantilever beam successively subjected to transversal and axial displacements in the case of perfect plasticity.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 13: Example 1, case 3: Time evolution of the axial forces $n_g$ (a) and bending moments $m_g$ (b) at the Gauss points of the cantilever beam successively subjected to transversal and axial displacements in the case of perfect plasticity.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 14: Example 1, case 3: stress path at the Gauss points of the cantilever beam successively subjected to transversal and then axial displacements. The bending moments $m_g$ at the Gauss points $g$ are plotted against the axial force $n_g$. The diamond shape is the yield surface.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 15: Example 1, case 3: time evolution of the plastic multipliers $\lambda_{1,3}$ (a) and $\lambda_{1,4}$ (b) at Gauss point 1 of the cantilever beam successively subjected to transversal and then axial displacements in the case of perfect plasticity.</figcaption></figure>
## 4.3 Example 2: Cantilever beam subjected to impact

The second example illustrates the relevance of the formulation for the joint assessment of plasticity and impact in a monolithic algorithm. The analysis focuses on the energy balance of the system, with specific interest on the evolution of energy dissipation depending on the parameters of the numerical scheme (time discretization, space discretization, and $\theta$ parameter).

<figure><img src="image_20.png" /><figcaption>Figure 16: Example 2: impacted cantilever beam.</figcaption></figure>

In this example, an elasto-plastic cantilever beam of length $L = 3m$, with isotropic and kinematic hardening, discretized with $n_{el} = 100$ equal length elements, is impacted at its free end by a spherical projectile of mass $m_{sphere}$, initially located at a distance $q_{y sphere} = g_{N}(t = 0) = 1m$ (Figure 16). The properties of the beam and of the projectile are summarized in Table 2. The restitution coefficient was set at $e = 0$ to emphasize energy dissipation due to impact, $\theta$ was set at $1/2$ to prevent artificial energy dissipation and the time step was set at $h = 10^{-4}s$.

<table><caption>Table 2: Example 2: properties of the cantilever beam and of the projectile.</caption><tr><td>Cross section area</td><td>A (m<sup>2</sup>)</td><td>3.4 &times; 10<sup>-3</sup></td></tr><tr><td>Moment of inertia</td><td>I (m<sup>4</sup>)</td><td>8.643 &times; 10<sup>-6</sup></td></tr><tr><td>Young modulus</td><td>E (MPa)</td><td>210 000</td></tr><tr><td>Material density</td><td>&rho; (kg/m<sup>3</sup>)</td><td>7701</td></tr><tr><td>Kinematic hardening modulus</td><td>H<sub>ki</sub> (MPa)</td><td>E/10 = 21 000</td></tr><tr><td>Isotropic hardening modulus</td><td>H<sub>is</sub> (MPa)</td><td>E/100 = 2100</td></tr><tr><td>Limit admissible tensile force</td><td>n<sub>P</sub> (N)</td><td>799 000</td></tr><tr><td>Limit admissible bending moment</td><td>m<sub>P</sub> (Nm)</td><td>33 840</td></tr><tr><td>Mass of the sphere</td><td>m<sub>sphere</sub> (kg)</td><td>250</td></tr><tr><td>Initial height of the sphere</td><td>q<sub>y sphere</sub>(t = 0) (m)</td><td>1</td></tr><tr><td>Restitution coefficient</td><td>e</td><td>0</td></tr></table>

The impact is first characterized by a succession of short interactions between the projectile and the beam (Figure 20 - light blue zone). As the restitution coefficient $e$ is equal to zero, each interaction lasts few timesteps. The existence and the duration of this transient phase is also related with the stiffness of the beam, and the high frequency dynamics generated by the FEM discretization. This phase is followed by a longer phase of permanent contact between the beam and the sphere (Figure 20 - blue zone). During this phase, the relative velocity vanishes while the impulse is always positive (Figure 20), as specified in the complementarity condition.

The impulses applied to the beam and to the projectile during the interaction induce the bouncing of the sphere and the bending of the beam (Figure 17). After the bouncing of the sphere, the undamped vibration of the beam is observed until a second impact. The beam oscillates around a negative equilibrium position since plastic strains developed during the first interaction phase.

**Energy transfer** In Figure 18, the discrete work of external forces $W_{ext 0}^{k+1}$, i.e. the sum of the kinetic energy $T_{k+1}$ and the free energy $\Psi_{k+1}$, is always equal to the sum of elastic mechanical energy $E_{k+1}$,
plastic dissipation $W_p^k$ and contact dissipation $W_c^k$. Indeed, setting $\theta = 1/2$ prevents artificial energy dissipation. Just before impact, the work of external forces $W_{ext0}^{k+1}$ is equal to the kinetic energy of the sphere. During the simulation, $W_{ext0}^{k+1}$ varies depending on the vertical position of the sphere only (for the sake of clarity, gravity was not applied to the beam in this example). The first interaction between the sphere and the beam results in substantial energy dissipation (Figure 19), due to plastic flow ($W_p^{k+1}$) and to contact ($W_c^{k+1}$), to a lesser extent, until a permanent contact establishes between the sphere and the beam (Figure 17 - blue colored area). Then, as the relative velocities of the beam and the sphere at the contact point are equal, $W_c^{k+1}$ remains constant while $W_p^{k+1}$ still increases until the sphere reaches its lower position (Figure 19). Just before the end of the first interaction, a second transition period occurs, with activation and desactivation of the contact, associated to very small energy dissipation at contact as the relative velocities are very small. The same energy dissipation process occurs for the second interaction with substantially smaller energy dissipation. Between the two interaction phases, no energy dissipation occurs, as $\theta = 1/2$.

<figure><img src="image_21.png" /><figcaption>Figure 17: Example 2: vertical positions of the beam free end and sphere along time.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 18: Time evolutions of the different energies and works involved in the impact of a cantilever beam by a sphere. The work of external forces $W_{ext0}^{k+1}$ is equal to the sum of mechanical energy $E_{k+1}$, plastic dissipation $W_p^{k+1}$ and contact dissipation $W_c^{k+1}$. And the elastic energy $E_{k+1}$ is the sum of the kinetic energy $T_{k+1}$ and the free energy $\Psi_{k+1}$.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 19: Time evolutions of the energies dissipated by plastic flow ($W_p^{k+1}$) and contact ($W_c^{k+1}$) during the impact of a sphere on a cantilever beam.</figcaption></figure>

**Time and space discretization** To illustrate the influence of the time and space discretization on the discrete energy balance, we consider the first period from $t = 0$ s to $t = 0.457$ s. As illustrated for the dissipation (Figure 21), energy dissipations due to plastic flow and contact tend to converge to a given value for refined time and space discretizations, i.e. increasing numbers of elements and decreasing time steps. A focus on the first percussion between the sphere and the beam (Figure 22) shows that convergence towards a threshold value of $p_N$ is reached for any spatial resolution. In addition, smaller time steps are necessary to reach convergence for finer spatial discretizations. The threshold value of $p_N$ depends on the spatial discretization: it decreases for increasing numbers of elements. This observation can be related to the open scientific question of the existence of percussion in this configuration for continuous beam (see Chatterjee (2004)). In this particular case, it seems that the percussion $p_n$ vanishes with the time-step and the mesh size. It also questions the relevance of the energy dissipation at contact in practice on such an application.
<figure><img src="image_24.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(c)</figcaption></figure>

Figure 20: Example 2: Time evolutions of the reaction impulse ($p_N$), gap function ($g_N$) and relative normal velocity ($v_N$) during the impact of a sphere on a cantilever beam. We can observe a contact activation period where several contact activations are observed as in the first colored area (light blue), before a permanent contact occurs (second colored area - blue).
<figure><img src="image_27.png" /><figcaption>Figure 21: Example 2: Convergence of the energy dissipated by plastic flow (a) and contact (b) after the end of the first phase of intermittent interactions between the sphere and the cantilever beam (end of the light blue area on Figure 20). The dissipated energies are plotted as a function of the time step $h$ for different numbers of elements $n_{el}$.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Figure 22: Example 2: First percussion as a function of the time step $h$ when the projectile impacts the free extremity of the cantilever beam.</figcaption></figure>

## 4.4 Example 3 : Frame made of beams

The last example illustrates the interest of the approach for the analysis of the response of more complex structures subjected to impact, in terms of deformation and energy dissipation.

The chosen structure is a rectangular frame impacted by a projectile (see Figure 23). The frame is 2 meters high and 3 meters large, made of HEB300 in S355 steel (total weigh of 804kg). The structure is clamped at its two extremities. The projectile (mass : $m_p = 1000\text{kg}$) impacts the structure at $0.67\text{m}$ from the ground with a horizontal velocity $v_i = 6\text{ m/s}$. The characteristics of the structure are summarized in Table 3.

The vertical beams are discretized with 30 elements and the horizontal one with 60 elements. The time step was set at $h = 10^{-5}\text{s}$, which guarantees the accurate approximation of the impact phenomenon and of the structure vibration, as the first natural frequency of the structure is 71Hz. We also set $\theta = 1/2$ to prevent artificial energy dissipation and $e = 0$ to maximize energy dissipation at contact.

<figure><img src="image_29.png" /><figcaption>Figure 23: Example 3: Schematic representation of the impacted frame.</figcaption></figure>
<table><caption>Table 3: Example 3: structure and simulation characteristics.</caption><tr><td>Cross section area</td><td>A (m²)</td><td>14.9 × 10⁻³</td></tr><tr><td>Moment of inertia</td><td>I (m⁴)</td><td>256 × 10⁻⁶</td></tr><tr><td>Young modulus</td><td>E (MPa)</td><td>210 000</td></tr><tr><td>Material density</td><td>ρ (kg/m³)</td><td>7701</td></tr><tr><td>Limit admissible tensile force</td><td>n<sub>P</sub> (MN)</td><td>5.3</td></tr><tr><td>Limit admissible bending moment</td><td>m<sub>P</sub> (MNm)</td><td>0.60</td></tr><tr><td>Kinematic hardening modulus</td><td>H<sub>kin</sub> (MPa)</td><td>E/10 = 21 000</td></tr><tr><td>Isotropic hardening modulus</td><td>H<sub>is</sub> (MPa)</td><td>E/100 = 2100</td></tr><tr><td>Number of elements per vertical beam</td><td>n<sub>el,v</sub></td><td>30</td></tr><tr><td>Number of elements per horizontal beam</td><td>n<sub>el,h</sub></td><td>60</td></tr><tr><td>Projectile mass</td><td>m<sub>p</sub> (kg)</td><td>1000</td></tr><tr><td>Projectile impact velocity</td><td>v<sub>i</sub> (m/s)</td><td>6</td></tr></table>

**Influence of elastoplasticity modelling** The analysis focused on the influence of the modelling of the material behavior on the energy dissipation and on the plastic strains location, based on comparisons between elasticity, perfect plasticity and plasticity with isotropic and kinematic hardening.

If an elastic material is considered, the impact results in a substantial kinetic energy transfer from the projectile to the structure. After a contact phase lasting around 0.25 × 10⁻³ s (Figure 24a), the projectile bounces on the structure with a velocity after rebound equal to -4.24 m s⁻¹ (Figure 24b), which corresponds to 71% of the impact velocity. Only a small part of the energy is dissipated at the contact point (1.85% of the impact energy - Figure 25). Consequently, the energy transferred to the structure (48.2% of the impact energy) mainly results in undamped oscillations of the frame.

In the case of perfect plasticity, most of the projectile energy is transferred to the structure and the velocity of the projectile is almost nil after impact (-0.12 m s⁻¹ i.e. 2% of the impact velocity). The impact energy is dissipated by plastic flow in the structure (19.3% of the impact energy) while a large amount results in structure oscillation (79.3% of the impact energy). The energy dissipated at the contact is similar to the energy dissipated for an elastic structure (1.4% of the impact energy).

The integration of hardening only slightly modifies the energy transfers and dissipation compared to perfect plasticity. A slightly smaller amount of energy is dissipated by plastic flow (12.1% of the impact energy) which results in a larger velocity of the projectile after impact -2.28 m s⁻¹, i.e. 38% of the impact velocity). The mechanical energy of the structure after impact and the energy dissipated at the contact (71.9% and 1.5% of the impact energy, respectively) are in the same range as for perfect plasticity.

On the contrary, the locations of the plastic strains in the structure are also significantly different (Figure 26). In the cases of perfect plasticity, local plastic hinges develop near the contact point, the extremities and the connections between beams while, for hardening plasticity, the plastic strains are more broadly distributed around these specific points of the structure. This result seems to be coherent with the results of (Khan, Ahmad, et al., 2021) even if the system is not exactly the same. The creation of plastic hinges for perfect plasticity drastically reduces structure resistance. This may explain the significantly larger duration of the interaction between the structure and the projectile for perfect plasticity (Figure 24b).

**Influence of the impact energy** Several simulations for increasing values of the impact energy have been done considering either three fixed projectile mass 250 kg, 600 kg and 1000 kg and increasing velocities or a fixed impact velocity 4 m s⁻¹, 6 m s⁻¹ and 8 m s⁻¹ and increasing masses. Table 4 shows the values of velocities and masses associated to each kinetic energy level for both cases.

Minor differences in terms of energy dissipation due to plastic flow are observed both for fixed masses and velocities (Figure 27a). More significant differences are observed for the energy dissipated at the contact point (Figure 27b). For increasing velocities and fixed masses, the energy dissipated linearly increases as a function of impact energy with larger increase for larger velocities. On the contrary, for increasing masses, i.e. fixed velocities, the dissipated energy first increases for increasing impact energy
<figure><img src="image_30.png" /><figcaption>Figure 24: Example 3: Projectile and structure impacted point horizontal positions (a) and projectile velocity (b) as a function of time for different material constitutive and evolution laws (elasticity, perfect plasticity and hardening plasticity).</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 25: Example 3: Distribution of the initial kinetic energy of the projectile at the end of the simulation. The bar length is the initial projectile energy.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 26: Example 3: Deformed structure 0.0265 seconds after the beginning of the contact (displacements multiplied 100 times) with maximum plastic curvature at each Gauss point.</figcaption></figure>
<table><caption>Table 4: Example 3: Analysis of the influence of the impact energy on the energy dissipation by plasticity and contact. The impact energy was varied using either fixed masses or fixed velocities. The values of the associated velocities and masses, respectively, are also given.</caption><thead><tr><th>Impact energy (kJ)</th><th>0.5</th><th>2</th><th>4.5</th><th>8</th><th>12.5</th><th>18</th><th>24.5</th><th>32</th></tr></thead><tbody><tr><td>v<sub>i</sub> (m/s) (m<sub>p</sub> = 1000kg)</td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td><td>8</td></tr><tr><td>v<sub>i</sub> (m/s) (m<sub>p</sub> = 600kg)</td><td>1.29</td><td>2.58</td><td>3.87</td><td>5.16</td><td>6.45</td><td>7.75</td><td>9.04</td><td>10.33</td></tr><tr><td>v<sub>i</sub> (m/s) (m<sub>p</sub> = 250kg)</td><td>2</td><td>4</td><td>6</td><td>8</td><td>10</td><td>12</td><td>14</td><td>16</td></tr><tr><td>m<sub>p</sub> (kg) (v<sub>i</sub> = 8m/s)</td><td>15.63</td><td>62.5</td><td>140.63</td><td>250</td><td>390.63</td><td>562.5</td><td>765.63</td><td>1000</td></tr><tr><td>m<sub>p</sub> (kg) (v<sub>i</sub> = 6m/s)</td><td>27.78</td><td>111.11</td><td>250</td><td>444.44</td><td>694.44</td><td>1000</td><td>1361.11</td><td>1777.78</td></tr><tr><td>m<sub>p</sub> (kg) (v<sub>i</sub> = 4 m/s)</td><td>62.5</td><td>250</td><td>562.5</td><td>1000</td><td>1562.5</td><td>2250</td><td>3062.5</td><td>4000</td></tr></tbody></table>

and, then, reaches threshold values that are larger for larger impact velocities. This difference is due to the use of a kinematic impact law that depends only on the relative velocities at the impact point and not on the masses of the interacting bodies. Different evolutions of the energies dissipated at the contact point could be observed for different impact laws, based on energy restitution coefficient for example.

<figure><img src="image_33.png" /><figcaption>Figure 27: Energy dissipated by plasticity (a) and by contact (b) in function of the impact energy. The impact energy was varied by using either fixed masses or fixed velocities.</figcaption></figure>

# 5 Conclusion

This work presents a monolithic algorithm for solving elastoplasticity and contact. The constitutive and evolutionary laws of elastoplasticity are based on the assumption of normal dissipativity. This can be expressed as a variational inequality, or equivalently, as a normal cone inclusion. This assumption is also equivalent to the Hill principle of maximum dissipation, and corresponds to the well-known framework of generalized standard materials. As usual in the nonsmooth contact dynamics method, the Signorini contact is also expressed as a variational inequality. In total, the dynamics of an elastoplastic system is written as a single variational inequality.

After spatial discretization using a standard iso-parametric finite element method, we end up with a finite-dimensional differential variational inequality that takes into account both the plasticity and the contact constraints. A time-stepping method is developed on the basis of the Moreau-Jean time stepping scheme in order to get a consistent scheme when impact occurs. Especially, some well-posedness results and a discrete energy balance are given which ensure that the scheme is well-posed and stable. The scheme benefits from excellent energy properties under conditions on the $\theta$ parameter. For $\theta = 1/2$, the
discrete energy balance shows that the numerical dissipation is removed.

The discrete one-step system can be recast into a constrained quadratic optimisation problem which has a unique solution under conditions of positive strain hardening, including perfect plasticity, and a sufficiently small time step. It is then possible to use numerous solving methods for mathematical programming to determine the solution.

Finally, numerical examples illustrate the possibilities of the method with the use of a LCP solver to solve elastoplastic dynamical systems with contact without resorting to the return mapping algorithm. The one-step system is solved a the machine accuracy with a single call of a LCP pivoting method. The singular points of the yield function are not a problem, and the dynamics under impact are validated by the numerical preservation of the energy balance of the system for different configurations.

The framework presented in this paper is general enough to be used in the context of perfect plasticity, or with positive strain hardening in 2D and 3D. The simulation of nonlinear plasticity yield criteria combining several mechanisms should not pose particular problems. The passage to finite strains is one of the objectives of future works using the work of Lee, 1968; Miehe, Apel, et al., 2002; Simo, Taylor, et al., 1985, as well as the use of different spatial discretization and formulation (Arbitrary Lagrangian-Eulerian (ALE) method, Particle Finite Element Method (PFEM), or Material Point Method (MPM)).

An immediate perspective is to complete the variational approach by linking the saddle point problem (69) to the work of Mielke (2003, 2005) and Mielke et al. (2002) integrating the inertia terms. The objective is to obtain a complete variational principle including dynamical effects and contact.

Another perspective is the extension of the proposed approach to non-associated plasticity and Coulomb friction using the De Saxcé (1992) bi-potential approach. This should lead to non-convex optimization problems for which it will be necessary to guarantee the existence of solutions and to find robust non-linear programming algorithms. Finally, a coupling with cohesive zone models for fracture seems quite simple, based on the recent work of Collins-Craft et al. (2022). The question of softening behaviors is more difficult because it would require using second gradient models or Cosserat media, but this crucial point could deserve a dedicated development.

## References

Acary, V., F. Cadoux, et al. (2011). "A formulation of the linear discrete Coulomb friction problem via convex optimization". In: *ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik* 91.2, pp. 155-175. ISSN: 1521-4001. [DOI], [HAL].

Acary, V. (2016). "Energy conservation and dissipation properties of time-integration methods for non-smooth elastodynamics with contact". In: *ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik* 96.5, pp. 585-603. [DOI], [ARXIV].

Acary, V., M. Bremond, et al. (Apr. 2023). *Siconos*. Version 4.4.0. URL: https://hal.science/hal-04056972.

Acary, V. and B. Brogliato (Jan. 30, 2008). *Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics*. Springer. 529 pp. ISBN: 978-3-540-75392-6.

Bathe, K. (1996). *Finite Element Procedures*. Second. Englewood Cliffs, New Jersey: Prentice-Hall.

Belgacem, F., P. Hild, and P. Laborde (1998). "The mortar finite element method for contact problems". In: *Mathematical and Computer Modelling* 28.4-8. Recent Advances in Contact Mechanics, pp. 263-271. ISSN: 0895-7177. [DOI], [OA].

Belytschko, T. and M. O. Neal (1991). "Contact-impact by the pinball algorithm with penalty and Lagrangian methods". In: *International Journal for Numerical Methods in Engineering* 31.3, pp. 547-572. [DOI].

Berga, A. and G. De Saxcé (1994). "Elastoplastic finite element analysis of soil problems with implicit standard material constitutive laws". In: *Revue Européenne des Éléments Finis* 3.3, pp. 411-456. [DOI].

Bertrand, D. et al. (2012). "Full-Scale Dynamic Analysis of an Innovative Rockfall Fence Under Impact Using the Discrete Element Method: From the Local Scale to the Structure Scale". In: *Rock Mechanics and Rock Engineering* 45.5, pp. 885-900. ISSN: 1434-453X. [DOI].
Bisbos, C. D., A. Makrodimopoulos, and P. M. Pardalos (Feb. 1, 2005). "Second-Order Cone Programming Approaches to Static Shakedown Analysis in Steel Plasticity". In: *Optimization Methods and Software* 20.1, pp. 25–52. ISSN: 1055-6788. (Visited on 09/27/2021). [DOI].

Bonnans, J.-F. et al. (2006). *Numerical optimization: theoretical and practical aspects*. Springer Science & Business Media.

Boyd, S. and L. Vandenberghe (2004). *Convex optimization*. Cambridge University Press. [DOI].

Brogliato, B. et al. (2006). "On the equivalence between complementarity systems, projected systems and differential inclusions". In: *Systems & Control Letters* 55.1, pp. 45–51. [DOI], [HAL].

Brüls, O., V. Acary, and A. Cardona (2014). "Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-α scheme". In: *Computer Methods in Applied Mechanics and Engineering* 281, pp. 131–161. [DOI], [HAL].

— (2018). "On the Constraints Formulation in the Nonsmooth Generalized-α Method". In: *Advanced Topics in Nonsmooth Dynamics*. Ed. by S. I. Publishing. Springer International Publishing, pp. 335–374.

Bruno, H. et al. (2020). "Return-mapping algorithms for associative isotropic hardening plasticity using conic optimization". In: *Applied Mathematical Modelling* 78, pp. 724–748. ISSN: 0307-904X. [DOI].

Capobianco, G. and S. R. Eugster (2018). "Time finite element based Moreau-type integrators". In: *International Journal for Numerical Methods in Engineering* 114.3, pp. 215–231. [DOI].

Capurso, M. and G. Maier (June 1, 1970). "Incremental Elastoplastic Analysis and Quadratic Optimization". In: *Meccanica* 5.2, pp. 107–116. ISSN: 1572-9648. (Visited on 09/24/2021). [DOI].

Carstensen, C., K. Hackl, and A. Mielke (2002). "Non-convex potentials and microstructures in finite-strain plasticity". In: *Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences* 458.2018, pp. 299–317. [DOI].

Chatterjee, A. (May 5, 2004). "The Short-Time Impulse Response of Euler-Bernoulli Beams". In: *Journal of Applied Mechanics* 71.2, pp. 208–218. ISSN: 0021-8936. (Visited on 08/28/2022). [DOI].

Chen, L., H. Wu, and T. Liu (2021). "Vehicle Collision with Bridge Piers: A State-of-the-Art Review". In: *Advances in Structural Engineering* 24.2, pp. 385–400. ISSN: 1369-4332. [DOI].

Chen, Q.-Z., V. Acary, et al. (2013). "A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints". In: *International Journal for Numerical Methods in Engineering* 96.8, pp. 487–511. [DOI], [HAL].

Cheng, L. et al. (2015). "A bipotential-based limit analysis and homogenization of ductile porous materials with non-associated Drucker-Prager matrix". In: *Journal of the Mechanics and Physics of Solids* 77, pp. 1–26. ISSN: 0022-5096. [DOI], [HAL].

Chouly, F., P. Hild, and Y. Renard (2015). "A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 49.2, pp. 481–502. [DOI], [OA].

Chouly, F. and Y. Renard (2018). "Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems". In: *Advanced Modeling and Simulation in Engineering Sciences* 5.1, pp. 1–38. [DOI], [OA].

Christensen, P. W. (2002a). "A nonsmooth Newton method for elastoplastic problems". In: *Computer Methods in Applied Mechanics and Engineering* 191.11-12, pp. 1189–1219. [DOI].

— (2002b). "A semi-smooth Newton method for elasto-plastic contact problems". In: *International Journal of Solids and Structures* 39.8, pp. 2323–2341. [DOI].

Chrysochoos, A. et al. (1989). "Plastic and dissipated work and stored energy". In: *Nuclear Engineering and Design* 114.3, pp. 323–333. [DOI], [HAL].

Collins-Craft, N., F. Bourrier, and V. Acary (2022). "On the formulation and implementation of extrinsic cohesive zone models with contact". In: *Computer Methods in Applied Mechanics and Engineering* 400, p. 115545. [DOI], [HAL].

Comi, C. and U. Perego (1995). "A unified approach for variationally consistent finite elements in elasto-plasticity". In: *Computer Methods in Applied Mechanics and Engineering* 121.1-4, pp. 323–344. ISSN: 0045-7825. [DOI].

Corradi, L. (1990). "Variational Statements and Mathematical Programming Formulations in Elastic-Plastic Analysis". In: *Mathematical Programming Methods in Structural Plasticity*. CISM International Centre for Mechanical Sciences. 978-3-211-82191-6: Springer Vienna, pp. 231–253. ISBN: 978-3-211-82191-6. [DOI].
Corradi, L. (1985). "Finite-element formulation of some extremum theorems of incremental plasticity". In: *Engineering Fracture Mechanics* 21.4, pp. 807-816. [DOI].

Dabaghi, F. et al. (2016). "A robust finite element redistribution approach for elastodynamic contact problems". In: *Applied Numerical Mathematics* 103, pp. 48-71. ISSN: 0168-9274. [DOI], [HAL].

Davison, L. (2008). *Fundamentals of shock wave propagation in solids*. Springer Berlin Heidelberg. [DOI], [OA].

De Saxcé, G. (Jan. 1, 1992). "Une Generalisation de l'inegalite de Fenchel et Ses Applications Aux Lois Constitutives". In: *C.R. Acad. Sci. Paris* 314.

Delbecq, J. et al. (1977). "Eléments finis en plasticité et visco-plasticité". In: *J. Mécanique Appliquée* 1.3, pp. 267-304.

Di Giacinto, D. et al. (2020). "A Novel Steel Damping System for Rockfall Protection Galleries". In: *Journal of Constructional Steel Research* 175, p. 106360. ISSN: 0143-974X. [DOI].

Di Stasio, J. et al. (2021). "An explicit time-integrator with singular mass for non-smooth dynamics". In: *Computational Mechanics* 68.1, pp. 97-112. [DOI], [HAL].

Donato, O. de and G. Maier (1972). "Mathematical Programming Methods for the Inelastic Analysis of Reinforced Concrete Frames Allowing for Limited Rotation Capacity". In: *International Journal for Numerical Methods in Engineering* 7.1, pp. 42-43. ISSN: 1097-0207. (Visited on 09/30/2021). [DOI].

Duan, L. and W.-F. Chen (Apr. 1, 1990). "A Yield Surface Equation for Doubly Symmetrical Sections". In: *Engineering Structures* 12.2, pp. 114-119. ISSN: 0141-0296. (Visited on 02/07/2022). [DOI].

Dubois, F., V. Acary, and M. Jean (2018). "The Contact Dynamics method: A nonsmooth story". In: *Comptes Rendus Mécanique* 346.3, pp. 247-262. [DOI], [OA].

Dumont, Y. and L. Paoli (2006). "Vibrations of a beam between obstacles. Convergence of a fully discretized approximation." In: *ESAIM: Mathematical Modelling and Numerical Analysis* 40.4, pp. 705-734. [DOI], [OA].

Dupire, S. et al. (2016a). "Novel Quantitative Indicators to Characterize the Protective Effect of Mountain Forests against Rockfall". In: *Ecological Indicators* 67, pp. 98-107. ISSN: 1470-160X. [DOI], [HAL].

— (2016b). "The protective effect of forests against rockfalls across the French Alps: Influence of forest diversity". In: *Forest Ecology and Management* 382, pp. 269-279. [DOI], [HAL].

Duvaut, G. and J. L. Lions (1976). *Inequalities in mechanics and physics*. Vol. 219. Springer Berlin Heidelberg. [DOI], [OA].

Facchinei, F. and J. S. Pang (2003). *Finite-dimensional Variational Inequalities and Complementarity Problems*. Vol. I & II. Springer Series in Operations Research. Springer New York. [DOI], [OA].

Feijoo, R. A. and N. Zouain (1988). "Formulations in rates and increments for elastic-plastic analysis". In: *International Journal for Numerical Methods in Engineering* 26.9, pp. 2031-2048. [DOI].

Felippa, C. A. (2004). "Introduction to finite element methods". In: *University of Colorado* 885.

Frémond, M., A. Pecker, and J. Salençon (1975). *Méthode variationnelle pour le matériau rigide-plastique*. Tech. rep. Rapport interne. Ecole Polytechnique, Laboratoire de Mécanique des Solides.

Fuente, H. M. de la and C. A. Felippa (1991). "Ephemeral penalty functions for contact-impact dynamics". In: *Finite Elements in Analysis and Design* 9.3, pp. 177-191. ISSN: 0168-874X. [DOI].

Germain, P. and E. Lee (1973). "On shock waves in elastic-plastic solids". In: *Journal of the Mechanics and Physics of Solids* 21.6, pp. 359-382. ISSN: 0022-5096. [DOI].

Guo, X., C. Zhang, and Z. Chen (2020). "Dynamic Performance and Damage Evaluation of a Scoured Double-Pylon Cable-Stayed Bridge under Ship Impact". In: *Engineering Structures* 216, p. 110772. ISSN: 0141-0296. [DOI].

Halphen, B. and Q. S. Nguyen (1975). "Sur les matériaux standard généralisés". In: *Journal de Mécanique* 14, pp. 39-63.

Heng, K., P. Jia, et al. (2022). "Vehicular impact resistance of highway bridge with seismically-designed UHPC pier". In: *Engineering Structures* 252, p. 113635. [DOI].

Heng, P., M. Hjiaj, et al. (2016). "A Simplified Model for Nonlinear Dynamic Analysis of Steel Column Subjected to Impact". In: *International Journal of Non-Linear Mechanics* 86, pp. 37-54. ISSN: 0020-7462. [DOI].

— (2017). "An Enhanced SDOF Model to Predict the Behaviour of a Steel Column Impacted by a Rigid Body". In: *Engineering Structures* 152, pp. 771-789. ISSN: 0141-0296. [DOI].
Hill, R. (Jan. 1, 1948). "A Variational Principle of Maximum Plastic Work in Classical Plasticity". In: *The Quarterly Journal of Mechanics and Applied Mathematics* 1.1, pp. 18-28. ISSN: 0033-5614. (Visited on 09/24/2021). [DOI].

- (1950). *The Mathematical Theory of Plasticity*. Oxford University Press. Oxford Classic Texts in the Physical Sciences. Oxford. ISBN: 0-19-850367-9.

Hiriart-Urruty, J.-B. and C. Lemaréchal (1993). *Convex Analysis and Minimization Algorithms*. Vol. I and II. Springer Berlin Heidelberg. [DOI], [OA].

Hjiaj, M., J. Fortin, and G. de Saxcé (2003). "A complete stress update algorithm for the non-associated Drucker-Prager model including treatment of the apex". In: *International Journal of Engineering Science* 41.10, pp. 1109-1143. [DOI].

Houlsby, G. T. (2019). "Frictional Plasticity in a Convex Analytical Setting". en. In: *Open Geomechanics* 1, 3, pp. 1-10. [DOI], [OA].

Hughes, T. J. R. (1980). "Generalization of selective integration procedures to anisotropic and nonlinear media". In: *International Journal for Numerical Methods in Engineering* 15.9, pp. 1413-1418. [DOI].

Jean, M. (July 20, 1999). "The Non-Smooth Contact Dynamics Method". In: *Computer Methods in Applied Mechanics and Engineering* 177.3-4, pp. 235-257. ISSN: 0045-7825. (Visited on 10/28/2021). [DOI], [HAL].

Jean, M. and J. J. Moreau (1987). "Dynamics in the presence of unilateral contacts and dry friction: a numerical approach". In: *Unilateral Problems in Structural Analysis — 2*. Ed. by G. Del Pietro and F. Maceri. Springer Vienna, pp. 151-196. [DOI], [HAL].

- (1992). "Unilaterality and dry friction in the dynamics of rigid body collections". In: *1st Contact Mechanics International Symposium*, pp. 31-48.

Johnson, K. L. (1985). *Contact mechanics*. Cambridge University Press. [DOI].

Kaewunruen, S., C. Ngamkhanong, and C. H. Lim (2018). "Damage and Failure Modes of Railway Pre-stressed Concrete Sleepers with Holes/Web Openings Subject to Impact Loading Conditions". In: *Engineering Structures* 176, pp. 840-848. ISSN: 0141-0296. [DOI], [OA].

Kanno, Y. (Dec. 1, 2016). "A Fast First-Order Optimization Approach to Elastoplastic Analysis of Skeletal Structures". In: *Optimization and Engineering* 17.4, pp. 861-896. ISSN: 1573-2924. (Visited on 10/27/2021). [DOI], [ARXIV].

- (2020). "A note on a family of proximal gradient methods for quasi-static incremental problems in elastoplastic analysis". In: *Theoretical and Applied Mechanics Letters* 10.5, pp. 315-320. [DOI], [OA].

Khan, A., I. Ahmad, et al. (Nov. 1, 2021). "A Modified Lemke Algorithm for Dynamic Rigid Plastic Response of Skeletal Structures". In: *Computers & Structures* 256, p. 106638. ISSN: 0045-7949. (Visited on 10/12/2021). [DOI].

Khan, A., D. L. Smith, and B. A. Izzuddin (May 1, 2013). "Investigation of Rigid-Plastic Beams Subjected to Impact Using Linear Complementarity". In: *Engineering Structures*. Engineering Structures: Modelling and Computations (Special Issue IASS-IACM 2012) 50, pp. 137-148. ISSN: 0141-0296. (Visited on 10/12/2021). [DOI].

Khenous, H. B., P. Laborde, and Y. Renard (2008). "Mass redistribution method for finite element contact problems in elastodynamics". In: *European Journal of Mechanics - A/Solids* 27.5, pp. 918-932. ISSN: 0997-7538. [DOI], [HAL].

Konyukhov, A. and K. Schweizerhof (2012). *Computational contact mechanics: geometrically exact theory for arbitrary shaped bodies*. Vol. 67. Springer Berlin Heidelberg. [DOI].

Krabbenhøft, K., A. V. Lyamin, S. W. Sloan, and P. Wriggers (2007). "An interior-point algorithm for elastoplasticity". In: *International Journal for Numerical Methods in Engineering* 69.3, pp. 592-626. [DOI].

Krabbenhøft, K., A. Lyamin, and S. Sloan (Mar. 1, 2007). "Formulation and Solution of Some Plasticity Problems as Conic Programs". In: *International Journal of Solids and Structures* 44.5, pp. 1533-1549. ISSN: 0020-7683. (Visited on 09/22/2021). [DOI], [OA].

Krabbenhøft, K., A. V. Lyamin, M. Hjiaj, et al. (2005). "A new discontinuous upper bound limit analysis formulation". In: *International Journal for Numerical Methods in Engineering* 63.7, pp. 1069-1088. [DOI].

Langlade, T. et al. (2021). "Modelling of Earthquake-Induced Pounding between Adjacent Structures with a Non-Smooth Contact Dynamics Method". In: *Engineering Structures* 241, p. 112426. ISSN: 0141-0296. [DOI].

Lebeau, G. and M. Schatzman (1984). "A wave problem in a half-space with a unilateral constraint at the boundary". In: *Journal of Differential Equations* 53.3, pp. 309-361. [DOI], [OA].
Lee, E. H. (1968). *Elastic-plastic deformation at finite strains*. Tech. rep. 183. Division of Engineering Mechanics. University of Stanford. [DOI].

Lee, E. H. and D. T. Liu (1964). "An Example of the Influence of Yield on High Pressure Wave Propagation". In: *Stress Waves in Anelastic Solids*. Springer Berlin Heidelberg, pp. 239–254. [DOI].

Levers, A. and A. Prior (1998). "Finite element analysis of shot peening". In: *Journal of Materials Processing Technology* 80-81, pp. 304–308. [DOI].

Maier, G. (Jan. 1, 1984). "Mathematical Programming Applications to Structural Mechanics: Some Introductory Thoughts". In: *Engineering Structures* 6.1, pp. 2–6. ISSN: 0141-0296. (Visited on 09/23/2021). [DOI].

Maier, G., C. Comi, and A. Corigliano (1991). "Extremum Properties of Finite-Step Solutions in Elastoplasticity with Nonlinear Mixed Hardening". In: pp. 99–102. [DOI].

Maier, G. (June 1, 1968a). "A Quadratic Programming Approach for Certain Classes of Non Linear Structural Problems". In: *Meccanica* 3.2, pp. 121–130. ISSN: 1572-9648. (Visited on 09/29/2021). [DOI].

— (Dec. 1, 1968b). "Quadratic Programming and Theory of Elastic-Perfectly Plastic Structures". In: *Meccanica* 3.4, pp. 265–273. ISSN: 1572-9648. (Visited on 09/29/2021). [DOI].

— (Sept. 1, 1969). "Shakedown Theory in Perfect Elastoplasticity with Associated and Nonassociated Flow-Laws: A Finite Element, Linear Programming Approach". In: *Meccanica* 4.3, pp. 250–260. ISSN: 1572-9648. (Visited on 09/29/2021). [DOI].

Majzoobi, G., R. Azizi, and A. Alavi Nia (2005). "A three-dimensional simulation of shot peening process using multiple shot impacts". In: *Journal of Materials Processing Technology* 164-165, pp. 1226–1234. [DOI].

Makrodimopoulos, A. and C. Martin (Jan. 1, 2005a). "A Novel Formulation of Upper Bound Limit Analysis as a Second-Order Cone Programming Problem". In.

— (Jan. 1, 2005b). "Limit Analysis Using Large-Scale SOCP Optimization". In.

Mandel, J. (1964). "Propagation des surfaces de discontinuité dans un milieu élastoplastique". In: *Stress Waves in Anelastic Solids*. Springer Berlin Heidelberg, pp. 331–340. [DOI].

Martin, J. et al. (1987). "Applications of mathematical programming concepts to incremental elastic-plastic analysis". In: *Engineering Structures* 9.3, pp. 171–176. ISSN: 0141-0296. [DOI].

Maugin, G. A. (1992). *The thermomechanics of plasticity and fracture*. Vol. 7. Cambridge University Press.

Meier, C., A. Popp, and W. A. Wall (2016). "A finite element approach for the line-to-line contact interaction of thin beams with arbitrary orientation". In: *Computer Methods in Applied Mechanics and Engineering* 308, pp. 377–413. [DOI], [ARXIV].

Meng, J. et al. (2020). "A smoothed finite element method using second-order cone programming". In: *Computers and Geotechnics* 123, p. 103547. ISSN: 0266-352X. [DOI].

Mercier, B. (1976). "Régularisation, approximation et résolution du problème des charges limites". fre. In: *Publications mathématiques et informatique de Rennes* S5, pp. 1–30.

Miehe, C., N. Apel, and M. Lambrecht (2002). "Anisotropic additive plasticity in the logarithmic strain space: modular kinematic formulation and implementation based on incremental minimization principles for standard materials". In: *Computer Methods in Applied Mechanics and Engineering* 191.47-48, pp. 5383–5425. [DOI].

Miehe, C., J. Schotte, and M. Lambrecht (2002). "Homogenization of inelastic solid materials at finite strains based on incremental minimization principles. Application to the texture analysis of polycrystals". In: *Journal of the Mechanics and Physics of Solids* 50.10, pp. 2123–2167. [DOI].

Miehe, C. (2002). "Strain-driven homogenization of inelastic microstructures and composites based on an incremental variational formulation". In: *International Journal for Numerical Methods in Engineering* 55.11, pp. 1285–1322. [DOI].

Mielke, A. (2003). "Energetic formulation of multiplicative elasto-plasticity using dissipation distances". In: *Continuum Mechanics and Thermodynamics* 15.4, pp. 351–382. [DOI].

— (2005). "Handbook of differential equations. Evolutionary equations". In: vol. 2. Elsevier. Chap. 6. Evolution of rate-independent systems.

Mielke, A., F. Theil, and V. I. Levitas (2002). "A Variational Formulation of Rate-Independent Phase Transformations Using an Extremum Principle". In: *Archive for Rational Mechanics and Analysis* 162.2, pp. 137–177. [DOI].
Moreau, J. J. (1970). "Sur les lois de frottement, de plasticité et de viscosité". In: *Comptes rendus de l'Académie des sciences. Série A - Sciences mathématiques* 271, pp. 608-611.

— (1971). *Fonctions de Résistance et Fonctions de Dissipation*. (Visited on 09/24/2021).

— (1974). "On unilateral constraints, friction and plasticity". In: *New Variational Techniques in Mathematical Physics*. Ed. by G. Capriz and G. Stampacchia. Springer Berlin Heidelberg, pp. 171-322. [DOI], [HAL].

— (1976). "Applications of convex analysis to the treatment of elasto-plastic systems". In: *Applications of Methods of Functional Analysis to Problems in Mechanics*. Ed. by P. Germain and B. Nayroles. Vol. 503. Lecture Notes in Mathematics. Springer Berlin Heidelberg, pp. 56-89. [DOI].

— (1977). "Evolution Problem Associated with a Moving Convex Set in a Hilbert Space". In: *Journal of Differential Equations* 26.3, pp. 347-374. [DOI], [HAL].

— (1986). "Une formulation du contact à frottement sec; application au calcul numérique". In: *Comptes rendus de l'Académie des sciences. Série 2, Mécanique, Physique, Chimie, Sciences de l'univers, Sciences de la Terre* 302.13, pp. 799-801. ISSN: 0764-4450.

— (1988a). "Bounded Variation in Time". In: *Topics in Nonsmooth Mechanics*. Ed. by J. J. Moreau, P. Panagiotopoulos, and G. Strang. Basel: Birkhäuser, pp. 1-74.

— (1988b). "Unilateral contact and dry friction in finite freedom dynamics". In: *Nonsmooth Mechanics and Applications*. Ed. by J. Moreau and P. P.D. CISM, Courses and lectures 302. Formulation mathematiques tire du livre Contacts mechanics. Wien- New York: Springer Vienna, pp. 1-82. [DOI], [HAL].

Nguyen, Q. S. (1977). "On the elastic-plastic initial-boundary value problem and its numerical integration". In: *International Journal for Numerical Methods in Engineering* 11.5, pp. 817-832. [DOI].

Nodargi, N. A. (2019). "An overview of mixed finite elements for the analysis of inelastic bidimensional structures". In: *Archives of Computational Methods in Engineering* 26.4, pp. 1117-1151. [DOI], [ARXIV].

Nouguier-Lehon, C. et al. (2013). "Surface impact analysis in shot peening process". In: *Wear* 302.1-2, pp. 1058-1063. [DOI].

Ortiz, M. and L. Stainier (1999). "The variational formulation of viscoplastic constitutive updates". In: *Computer Methods in Applied Mechanics and Engineering* 171.3-4, pp. 419-444. ISSN: 0045-7825. [DOI], [HAL].

Paoli, L. and M. Schatzman (2002a). "A Numerical Scheme for Impact Problems I: The One-Dimensional Case". In: *SIAM Journal on Numerical Analysis* 40.2, pp. 702-733. [DOI], [HAL].

— (2002b). "A Numerical Scheme for Impact Problems II: The Multi-Dimensional Case". In: *SIAM Journal on Numerical Analysis* 40.2, pp. 734-768. [DOI].

— (2007). "Numerical simulation of the dynamics of an impacting bar." English. In: *Computer Methods in Applied Mechanics and Engineering* 196.29-30, pp. 2839-2851. [DOI], [HAL].

Parikh, N. (2014). "Proximal algorithms". In: *Foundations and Trends® in Optimization* 1.3, pp. 127-239. [DOI].

Pastor, F., P. Thoré, et al. (2008). "Convex optimization and limit analysis: Application to Gurson and porous Drucker-Prager materials". In: *Engineering Fracture Mechanics* 75.6, pp. 1367-1383. ISSN: 0013-7944. [DOI].

Pastor, F., J. Pastor, and D. Kondo (2015). "Numerical limit analysis and plasticity criterion of a porous Coulomb material with elliptic cylindrical voids". In: *Comptes Rendus Mécanique* 343.3, pp. 199-209. ISSN: 1631-0721. [DOI].

Pellegrino, S. (1993). "Structural computations with the singular value decomposition of the equilibrium matrix". In: *International Journal of Solids and Structures* 30.21, pp. 3025-3035. ISSN: 0020-7683. [DOI].

Pereira, N. Z., L. A. Borges, and M. B. Hecke (1988). "A force method for elastic-plastic analysis of frames by quadratic optimization". In: *International Journal of Solids and Structures* 24.2, pp. 211-230. ISSN: 0020-7683. [DOI].

Popp, A., M. W. Gee, and W. A. Wall (2009). "A finite deformation mortar contact formulation using a primal-dual active set strategy". In: *International Journal for Numerical Methods in Engineering* 79.11, pp. 1354-1391. [DOI].

Reddy, B. and J. Martin (1991). "Algorithms for the solution of internal variable problems in plasticity". In: *Computer Methods in Applied Mechanics and Engineering* 93.2, pp. 253-273. ISSN: 0045-7825. [DOI].
Renard, Y. (2010). "The singular dynamic method for constrained second order hyperbolic equations: Application to dynamic contact problems". In: *Journal of Computational and Applied Mathematics* 234.3, pp. 906-923. ISSN: 0377-0427. [DOI], [OA].

Romano, G., L. Rosati, and F. de Sciarra (1993). "Variational principles for a class of finite step elastoplastic problems with non-linear mixed hardening". In: *Computer Methods in Applied Mechanics and Engineering* 109.3-4, pp. 293-314. [DOI].

Scherzinger, W. (2017). "A return mapping algorithm for isotropic and anisotropic plasticity models using a line search method". In: *Computer Methods in Applied Mechanics and Engineering* 317, pp. 526-553. [DOI], [OA].

Schindler, T. and V. Acary (2013). "Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: Definition and outlook". In: *Mathematics and Computers in Simulation* 95, pp. 180-199. ISSN: 0378-4754. [DOI], [HAL].

Schindler, T., S. Rezaei, et al. (2015). "Half-explicit timestepping schemes on velocity level based on time-discontinuous Galerkin methods". In: *Computer Methods in Applied Mechanics and Engineering* 290.15, pp. 250-276. [DOI], [HAL].

Sha, Y., J. Amdahl, and C. Dørum (2021). "Numerical and Analytical Studies of Ship Deckhouse Impact with Steel and RC Bridge Girders". In: *Engineering Structures* 234, p. 111868. ISSN: 0141-0296. [DOI], [OA].

Shimizu, W. and Y. Kanno (2018). "Accelerated proximal gradient method for elastoplastic analysis with von Mises yield criterion". In: *Japan Journal of Industrial and Applied Mathematics* 35.1, pp. 1-32. [DOI].

- (2020). "A note on accelerated proximal gradient method for elastoplastic analysis with Tresca yield criterion". In: *Journal of the Operations Research Society of Japan* 63.3, pp. 78-92. [DOI], [OA].

Simo, J. C. and T. Honein (1990). "Variational Formulation, Discrete Conservation Laws, and Path-Domain Independent Integrals for Elasto-Viscoplasticity". In: *Journal of Applied Mechanics* 57.3, pp. 488-497. ISSN: 0021-8936. [DOI].

Simo, J. C. and T. J. R. Hughes (1998). *Computational Inelasticity*. Springer-Verlag. 405 pp. ISBN: 978-0-387-22763-4. [DOI], [OA].

Simo, J., R. Taylor, and K. Pister (1985). "Variational and projection methods for the volume constraint in finite deformation elasto-plasticity". In: *Computer Methods in Applied Mechanics and Engineering* 51.1-3, pp. 177-208. ISSN: 0045-7825. [DOI].

Smith, D. L. (1990). *Mathematical programming methods in structural plasticity*. Springer.

Souza Neto, E. A. de, D. Peric, and D. R. Owen (2011). *Computational methods for plasticity: theory and applications*. John Wiley & Sons.

Suquet, P.-M. (1981). "Sur les équations de la plasticité: existence et régularité des solutions". In: *J. Mécanique* 20.1, pp. 3-39.

Tangaramvong, S., F. Tin-Loi, and C. Song (2012). "A Direct Complementarity Approach for the Elastoplastic Analysis of Plane Stress and Plane Strain Structures". In: *International Journal for Numerical Methods in Engineering* 90.7, pp. 838-866. ISSN: 1097-0207. (Visited on 10/27/2021). [DOI].

Wieners, C. (2007). "Nonlinear solution methods for infinitesimal perfect plasticity". In: *ZAMM* 87.8-9, pp. 643-660. [DOI], [OA].

Wilkins, M. L. (1963). *Calculation of elastic-plastic flow*. Tech. rep. California Univ Livermore Radiation Lab.

Wriggers, P. (2006). *Computational contact mechanics*. 2nd. Springer Berlin Heidelberg. [DOI], [OA].

Wright, S. J. (1996). *Primal-Dual Interior-Point Methods*. Philadelphia: Society for Industrial and Applied Mathematics. [DOI].

Yang, B. (2006). *Mortar finite element methods for large deformation contact mechanics*. Vol. 68. 02.

Zeng, Q., C. D. Stoura, and E. G. Dimitrakopoulos (2018). "A Localized Lagrange Multipliers Approach for the Problem of Vehicle-Bridge-Interaction". In: *Engineering Structures* 168, pp. 82-92. ISSN: 0141-0296. [DOI].

Zhang, X., K. Krabbenhøft, D. Pedroso, et al. (2013). "Particle finite element analysis of large deformation and granular flow problems". In: *Computers and Geotechnics* 54, pp. 133-142. [DOI].
Zhang, X., D. Sheng, et al. (2017). “Lagrangian modelling of large deformation induced by progressive failure of sensitive clays with elastoviscoplasticity”. In: *International Journal for Numerical Methods in Engineering* 112.8, pp. 963–989. [DOI], [HAL].

Zhang, X. (2014). “Particle finite element method in geomechanics”. PhD thesis. University of Newcastle Australia.

Zhang, X., K. Krabbenhøft, D. Sheng, et al. (2015). “Numerical simulation of a flow-like landslide using the particle finite element method”. In: *Computational Mechanics* 55.1, pp. 167–177. [DOI].

Zhao, R. et al. (2022). “A sequential linear complementarity problem for multisurface plasticity”. In: *Applied Mathematical Modelling* 103, pp. 557–579. ISSN: 0307-904X. [DOI].

Zheng, H., T. Zhang, and Q. Wang (2020). “The mixed complementarity problem arising from non-associative plasticity with non-smooth yield surfaces”. In: *Computer Methods in Applied Mechanics and Engineering* 361, p. 112756. ISSN: 0045-7825. [DOI].

Zhou, X.-W. et al. (2023). “A mixed selective edge-based smoothed PFEM with second-order cone programming for geotechnical large deformation analysis”. In: *Computers and Geotechnics* 153, p. 105047. ISSN: 0266-352X. [DOI], [OA].

Ziegler, H. (Mar. 1, 1958). “An Attempt to Generalize Onsager’s Principle, and Its Significance for Rheological Problems”. In: *Zeitschrift für angewandte Mathematik und Physik* ZAMP 9.5-6, pp. 748–763. ISSN: 1420-9039. (Visited on 09/28/2021). [DOI].

— (1962). *Some Extremum Principles in Irreversible Thermodynamics, with Application to Continuum Mechanics*. Swiss Federal Institute of Technology. 242 pp.

Zienkiewicz, O. C., R. L. Taylor, and J. Z. Zhu (2005). *The finite element method: its basis and fundamentals*. Elsevier.

## A Finite element notation

The FE discretization of the displacement in an element labeled by $e \in [1, n_{el}]$ is given by

$$u_e(x) = \bar{N}_e(x)u_e = N_e(x)u, \qquad (140)$$

where $u_e \in \mathbb{R}^{d_u}$ is the vector composed of the nodal displacements of the element and $u \in \mathbb{R}^n$ the nodal displacement vector of the whole structure. The matrix $N_e$ denotes the shape function matrix of the element $e$.

Note that $u$ is not the simple concatenation of $u_e$, in other words $u \neq \text{col}(u_e, e \in [1, n_{el}])$ but is related by a local-to-global mapping $A_e$ such that $u_e = A_e u$. This mapping is also used for the assembly of the structural matrices. To avoid too complex notation related to the assembly process, we prefer to formulate the problem directly with $N_e(x)$ (see Bathe, 1996).

The strain, in vector notation, $\varepsilon_e \in \mathbb{R}^{d_e}$, using Voigt notation in an element is given by

$$\varepsilon_e(x) = B_e(x)u, \qquad (141)$$

where $B_e(x)$ is the strain-displacement matrix for the element $e$ obtained by applying (140) in the definition of the strain tensor from the displacement. The nodal velocity vector $v_e$ is assumed to be given by the time derivative of the nodal displacement vector leading to

$$v_e(x) = N_e(x)v. \qquad (142)$$

Using the same Galerkin approximation for virtual velocities, the principle of virtual power yields the following space discretized equation for the element $e$

$$\int_{\Omega_e} N_e^{\top}(x) N_e(x) dm(x, t) \dot{v}(t) + \int_{\Omega_e} B_e^{\top}(x) \sigma(x, t) dv(x) = \int_{\Omega_e} N_e^{\top}(x) f(x, t) dv(x) + \int_{\partial \Omega_e} N_e^{\top}(x) \tau(x, t) ds(x) + \int_{\partial \Omega_e} N_e^{\top}(x) r(x, t) ds(x). \qquad (143)$$

Using the following standard notation,

$$M_e = \int_{\Omega_e} N_e^{\top}(x) N_e(x) dm(x), \qquad (144)$$
for the consistent mass matrix, and

$$
\begin{align}
f_{\text{int},e}(t) &= \int_{\Omega_e} B_e^\top(x) \sigma(x,t) dv(x) \nonumber \\
f_{\text{ext},e}(t) &= \int_{\Omega_e} N_e^\top(x) f(x,t) dv(x) + \int_{\partial\Omega_e} N_e^\top(x) \tau(x,t) ds(x) \tag{145} \\
r_e(t) &= \int_{\partial\Omega_e} N_e^\top(x) r(x,t) ds(x), \nonumber
\end{align}
$$

for the nodal internal, external and contact forces, we obtain the contribution to equation of motion of
the element e as

$$
M_e \dot{v}(t) + f_{\text{int},e}(t) = f_{\text{ext},e}(t) + r_e(t). \quad (146)
$$

The discrete equation are obtained by summing the contributions of the elements yielding

$$
M \dot{v}(t) + f_{\text{int}}(t) = f_{\text{ext}}(t) + r(t). \tag{147}
$$

where

$$
M = \sum_{e=1}^{n_{\text{el}}} M_e(t), \quad f_{\text{int}}(t) = \sum_{e=1}^{n_{\text{el}}} f_{\text{int},e}, \quad f_{\text{ext}}(t) = \sum_{e=1}^{n_{\text{el}}} f_{\text{ext},e}, \quad \text{and} \quad r(t) = \sum_{e=1}^{n_{\text{el}}} r_e. \tag{148}
$$

**B Discrete equilibrium matrix**

The evaluation of the internal forces using the Gauss quadrature rule is

$$
f_{\text{int},e} = \sum_{g=1}^{n_{e,g}} \omega_g B_e^\top(x_g) \sigma_e(x_g) = \sum_{g=1}^{n_{e,g}} \sqrt{\omega_{e,g}} B_e^\top(x_g) \sqrt{\omega_{e,g}} \sigma_e(x_g) \quad (149)
$$

where $x_{e,g} \in \mathbb{R}^d$, $g \in [1, n_{e,g}]$ denotes the Gauss points associated with their weights $w_{e,g}$ that are positive. Using the standard Gauss quadrature rule, the strain is also evaluated at the Gauss points. We define the notation

$$
\epsilon_{e,g} = B_{e,g} u_g, \quad \text{with} \quad \epsilon_{e,g} = \epsilon_e(x_g) \text{ and } B_{e,g} = B_e(x_g) \tag{150}
$$

The strains at Gauss points in an element and in the structure are collected scaled by their weights as
follows

$$
\epsilon_e = \operatorname{col}(\sqrt{w_{e,g}} \epsilon_{e,g}, g \in [1, n_{e,g}]) \quad \text{and} \quad \epsilon = \operatorname{col}(\epsilon_e, e \in [1, n_{el}]) \tag{151}
$$

and, doing in the same way for the matrix $B_{e,g}$, we get

$$
\epsilon_e = B_e u \quad \text{and} \quad \epsilon = B u, \tag{152}
$$

with

$$
B_e = \operatorname{col}(\sqrt{w_{e,g}} B_{e,g}, g \in [1, n_{e,g}]) \quad \text{and} \quad B = \operatorname{col}(B_e, e \in [1, n_{el}]). \tag{153}
$$

In the same way, the stresses, evaluated at Gauss points of an element are gathered as follows

$$
\sigma_{e,g} = \sigma_e(x_g), \quad \sigma_e = \operatorname{col}(\sqrt{w_{e,g}} \sigma_{e,g}, g \in [1, n_{e,g}]) \quad \text{and} \quad \sigma = \operatorname{col}(\sigma_e, e \in [1, n_{el}]) \quad (154)
$$

The matrix *B* appears as the discrete divergence operator, also called the equilibrium matrix in the
classical force method.
C Beam element

For the sake of completeness, we recall in this appendix the isoparametric beam element used in the numerical applications. Third order shape functions are used for the bending and second order one for the tension is order to match the location of the Gauss points. Apart from that, the development follows the classical approach that can be found in (Felippa, 2004). The elementary form function matrix is

$$
\bar{N}_e = \begin{bmatrix} \frac{\xi}{2}(\xi - 1) & 0 & 0 & 1 - \xi^2 & \frac{\xi}{2}(\xi + 1) & 0 & 0 \\ 0 & \frac{1}{4}(1 - \xi)^2(2 + \xi) & \frac{L_e}{8}(1 - \xi)^2(1 + \xi) & 0 & 0 & \frac{1}{4}(1 + \xi)^2(2 - \xi) & -\frac{L_e}{8}(1 + \xi)^2(1 - \xi) \end{bmatrix} \tag{155}
$$

The derivative matrix of N is then:

$$
B_e = \begin{bmatrix} \frac{2\xi - 1}{L_e} & 0 & 0 & \frac{-4\xi}{L_e} & \frac{2\xi + 1}{L_e} & 0 & 0 \\ 0 & \frac{6\xi}{L_e^2} & \frac{3\xi - 1}{L_e} & 0 & 0 & \frac{-6\xi}{L_e^2} & \frac{3\xi + 1}{L_e} \end{bmatrix} \quad (156)
$$

The elementary consistent mass matrix associated with the element is :

$$
M_e = \frac{\rho A L_e}{420} \begin{bmatrix} 56 & 0 & 0 & 28 & -14 & 0 & 0 \\ 0 & 156 & 22L_e & 0 & 0 & 54 & -13L_e \\ 0 & 22L_e & 4L_e^2 & 0 & 0 & 13L_e & -3L_e^2 \\ 28 & 0 & 0 & 224 & 28 & 0 & 0 \\ -14 & 0 & 0 & 28 & 56 & 0 & 0 \\ 0 & 54L_e & 13L_e & 0 & 0 & 156 & -22L_e \\ 0 & -13L_e & -3L_e^2 & 0 & 0 & -22L_e & 4L_e^2 \end{bmatrix} \quad (157)
$$

where $L_e$ is the element length, $A$ is the cross section area and $\rho$ is the material density.
The elementary elasticity matrix $C_e$, which relates stresses $\sigma_e$ and elastic strains $\epsilon_e$, is :

$$
C_e = S_e^{-1} = E \begin{bmatrix} A & 0 & 0 & 0 \\ 0 & I & 0 & 0 \\ 0 & 0 & A & 0 \\ 0 & 0 & 0 & I \end{bmatrix} \quad (158)
$$

The elementary hardening matrix $D_e$, which relates internal forces $a_e$ and hardening parameters $\alpha_e$, is:

$$
D_e = \begin{bmatrix} H_{ki}A & 0 & 0 & 0 & 0 & 0 \\ 0 & H_{ki}I & 0 & 0 & 0 & 0 \\ 0 & 0 & H_{is}I & 0 & 0 & 0 \\ 0 & 0 & 0 & H_{ki}A & 0 & 0 \\ 0 & 0 & 0 & 0 & H_{ki}I & 0 \\ 0 & 0 & 0 & 0 & 0 & H_{is}I \end{bmatrix} \quad (159)
$$

where $H_{ki}$ and $H_{is}$ are the kinematic and isotropic hardening moduli, respectively.