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

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.

Notation The notation used are quite classical and are not defined in an exhaustive way. We follow the notation introduced in the standard textbooks. For vector and tensors, we choose the following notation: The set T d s is the set of symmetric tensors of order s in dimension d. For a convex set C ⊂ IR n , the normal cone to C at x is defined by N C (x) = {s ∈ IR n | s ⊤ (y − x) ⩽ 0, for all y ∈ C}. (2) For a convex function φ, its subdifferential is denoted by ∂φ. The indicator function of C is denoted i C and we have N C = ∂i C . For a function F : IR n → IR n , the inclusion defines a (finite dimensional) variational inequality.

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;Majzoobi 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 onesided 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 (1948Hill ( , 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 (1970Moreau ( , 1971Moreau ( , 1974Moreau ( , 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 (1958Ziegler ( , 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 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, 2003Mielke, , 2005Mielke 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).
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 detailed in the Simo and Hughes (1998) uncontested book. 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 andMoreau (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 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(Brüls et al., , 2018Capobianco 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 finitedimensional 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.

Elasto-plasticity dynamics with unilateral contacts
The formulation of a model for the elasto-plasticity dynamics is based on the standard textbooks (Maugin, 1992;Nguyen, 2000). 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 (1970Moreau ( , 1974Moreau ( , 1986 in the context of plasticity, friction and contact, introducing super-potentials and pseudopotentials of dissipation. In this work, the notion of Generalized Standard Materials (GSM) pioneered by Halphen and Nguyen (1975) is used for the constitutive models with the Moreau-Ziegler assumption of normal dissipation. The presentation of the elasto-plastic model with hardening is made in the GSM framework in order to satisfy the principles of thermo-mechanics.

Principle of virtual powers
The material body is assumed to be an open set Ω ⊂ IR d , d ∈ 1, 3 for t ∈ [t 0 , T ] with a smooth boundary ∂Ω. 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 ρ(x), the density of the material, that is dm(x) = ρ(x)dv(x). The displacement is denoted by u(x, t) : Ω × [0, T ] → IR d and the velocity is denoted by v(x, t) : Ω × [0, T ] → IR d . We assume that we are in the framework of the smallperturbation hypothesis (SPH) and the strain rate is given byε(x, t), where the second order symmetric strain tensor ε is given by (∇ denotes the gradient) Principle of virtual power The principle of virtual power in the setting of classical mechanics postulates 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 : where 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 where σ ∈ T d 2 is the Cauchy stress tensor, the virtual power of volume forces f and the power of applied contact forces t(x, t) on ∂Ω F ⊂ ∂Ω as and the unknown contact forces r(x, t) on ∂Ω C ⊂ ∂Ω as Energy balance The kinetic energy is defined as and its time derivative, when it is defined, gives the power of inertial forces as d dt Choosing the virtual field v ⋆ equal to the actual velocity field v, we obtain the energy balance from the principle of virtual power: d dt 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 IR 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.

Elasto-plasticity in the Generalized Standard Material (GSM) setting
As usual under the SPH for the elasto-plastic solids, the strain tensor is decomposed in an additive way into an elastic part ε e , i.e. the reversible part of the strain, and a plastic part ε p as follows Free energy and state laws We start by the definition of the free energy per unit of volume ψ(ε e , α) in the isothermal setting (the free energy does not depend on the temperature T ) assuming a separation in ε e and α: where α is a tensor of hardening parameters, ψ e is the elastic potential energy and ψ p is the plastic potential energy due to the hardening process. The state laws govern the generalized (driving) forces σ(x, t) ∈ T d and a(x, t) as the derivative of the free energy: Following the simplest and classical case presented in Maugin (1992), we choose quadratic potentials, ψ e (ε e ) = 1 2 ε e : E : ε e and ψ p (α) = 1 2 α · D · α with E a symmetric 4-th order elastic tensor and D a symmetric second order tensor, the state laws are then linear : 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 must be positive. Developingψ(ε e , α), we get With the state laws defined above, this can be simplified in The total intrinsic dissipation (or dissipated power) D(t) and the total free energy Ψ are defined as Then, the power of internal forces is expressed as that is 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 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 ε p and α 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 σ and forces a, and the complementary laws are defined by the following differential inclusion The complementary laws (24) imply the Clausius-Duhem inequality (23) if the origin belongs to C. Using the definition of the normal cone (2), the complementary laws are equivalent to the following variational inequalitẏ also know as the Hill's principle of maximum dissipatioṅ 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 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.

Unilateral contact
Ω Ω n(x) Figure 1: A contact pair and a local contact frame.
In this work, the unilateral contact is considered as perfect unilateral constraint. For a contact point C A of position x ∈ ∂Ω C candidate to contact, we assume that we are able to associate a unique contact point C B of positionx = Γ(x) ∈ ∂Ω C such that the normal gap function is defined as where n(x) is the outward normal vector at x (see Figure 1). When the contact occurs with a fixed obstacle, the pointx is possibly defined on the boundary of this obstacle. For a detailed description of the practical computation of Γ, 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 ψ s of the system is defined with the indicator function of IR + , The Signorini condition is then obtained by sub-differentiating 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: It is also possible to write a stronger condition at the velocity level that implies the standard Signorini condition as where T IR+ (g N ) is the tangent cone of IR + (see Moreau (1988b)) and the normal relative velocity v N is defined as the time derivative of the gap function Equivalently, the complementarity condition at the velocity level is The contact force r(x, t) are related to the normal contact reaction by duality using the virtual power of contact forces as Power of contact forces. The power done by the contact forces vanishes almost everywhere due to (34), and does not dissipate or release energy: 3 Space and Time discretizations

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 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 σ e . A discretization of the stresses σ 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 Hence, the space-discretized equation of motion (37) reads as 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: 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 α ∈ 1 . . . m , the gap function is written as a function of the nodal displacement u as g α N (u). The normal relative velocity is then given by The vector v N is defined by collecting all the relative velocities at the contacts points, v N := col(v α N , α ∈ 1, m ), and the matrix H is defined as By duality, the contact forces are expressed as The Signorini condition at the velocity level is then written as −r N ∈ N T IR m 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 di 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 is extended as measure equation as If we define r N as the density of di N with respect to the Lebesgue measure, the equation (45) appears as the equation of motion, valid almost everywhere. At any time point t i , the measure equation amounts to solving the impact equation where the percussion p i,N is the normal percussion, i.e. the density of di 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 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 (44) are formulated as a measure inclusion as 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 ∈ [t 0 , T ] read as with the initial conditions v(t 0 ) = v 0 , u(t 0 ) = u 0 and ε p (t 0 ) = 0 for some t 0 .
In order to prepare the time discretization of (50), the differential equation are separated from the algebraic equations and the normal cone inclusion introducing two slack variables y = −α and z = −ε p . Assuming that the elastic law is valid at time t 0 , that is, σ(t 0 ) = Eε(t 0 ), the elasticity can be expressed asσ or equivalently, with the compliance matrix S = E −1 as With the slack variable y and assuming that a(t 0 ) = −Dα(t 0 ), the linear hardening law reads Finally, with the slack variables y and z, the normal cone inclusion for the plasticity flow rule is given by The complete system of equations for the elasto-plastic evolution with contact reduces to with the initial conditions v(t 0 ) = v 0 , u(t 0 ) = u 0 and ε p (t 0 ) = 0. Let us recall that we assume that a(t 0 ) = −Dα(t 0 ) and σ(t 0 ) = Eε(t 0 ).

Time discretization with an extended Moreau-Jean time-stepping scheme
The time-discretization of the system (55) is given by an extension of the Moreau-Jean scheme for nonsmooth dynamics (Jean and Moreau, 1987) and by the catching-up algorithm for plasticity (Moreau, 1976). The dynamics equation in the first line of (55) is discretized using the definition of the differential measure over a time interval (t k , t k+1 ] as 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: 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 (139) are approximated with a θ−method, for θ > 0 as where the following notation is used x k+θ = θx k+1 + (1 − θ)x k for any variable x. The non impulsive dynamical equations in (55) are also discretized using a θ−method The normal cone inclusion governing the plastic flow rule is discretized as follows 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 Hence, the impact law is written as where v α N,k+1 = H α,⊤ v k+1 . By collecting only the contact variables for the index set I k as v N := col(v α N , α ∈ I k ) and p N := col(p α N , α ∈ I k ), the contact impact law is given by where m is the cardinal of I k . Using the reduced matrix Summary of the space and time discretized equations Altogether, the space and time discretized equations are given by where u k+1 = u k + hv k+θ can be solved afterwards, and Using the relation x k+θ − x k = θ(x k+1 − x k ), the discretized equations can be formulated in terms of variables approximated at t k+θ

Well-posedness
To study the well posedness to the problem (66), 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 Proposition 1. Using the first order optimality conditions of (67), the solutions (v,ε, σ, a, v N ) of (67) are solutions of (66) In order to derive the optimality conditions for (67), we consider the Lagrangian function The optimality conditions of the problem (67) is given by 0 ∈ ∂L(v,ε, σ, a, v N , λ, µ). Computing the subgradients, we obtain for the optimality conditions Note that hθσ appears as the Lagrange multiplier that enforces the condition Bv =ε. Introducing the variables z, y, p N such that and simplifying the equations, we obtain This completes the proof.
Let us introduce the following assumptions Assumption 1. The matrices M , S and D are symmetric definite positive matrices.
To prove that the min-max problem has solutions, i.e. saddle points exist, we consider a reduced version of the problem (67) where the equality constraints Bv =ε and θv N = H ⊤ v − (1 − θ)v N,k are substituted. This yields a reduced saddle point problem with its associated Lagrangian function of the form Assuming that M , S and D are symmetric definite positive matrices, the Lagrangian functionL is a convex function in v and a concave function of (σ, a). With the assumption (72) , the following coercivity properties hold: From Theorem 4.3.1 in (Hiriart-Urruty and Lemaréchal, 1993), the set of solutions of the saddle point problem (73) is a nonempty compact and convex set. The strict convexity ofL with respect to v for all (σ, a) ∈ C and the strict concavity with respect to σ, a for all v such that H ⊤ v 0 + (θ(1 + e) − 1)v N,k ⩾ 0 imply the uniqueness of the solution of (73). To obtain a complete solution of (67) it suffices to computė ε = Bv and θv N = Proposition 3. Under Assumptions 1 and 2, the problem (66) has a solution for (v k+θ , σ k+θ , a k+θ , v N,k+1 ) and (z k+θ , y k+θ , p N,k+1 ). Furthermore, the solution is unique for (v k+θ , σ k+θ , a k+θ , v N,k+1 , z k+θ , y k+θ ).
As it is noted in (Hiriart-Urruty and Lemaréchal, 1993, Remark 4.3.2), the coercivity properties (75) imply that the following necessary optimality conditions has a solution With our setting we get Let us introduce θv N such that and z, y, p Then, the optimality conditions can be written as A solution of (80) is a solution of the problem (66). The uniqueness of (v k+θ , σ k+θ , a k+θ , v N,k+1 , z k+θ , y k+θ ) 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 (72) 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 (65) or the min-max problem (67) 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 (65) 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 (66), we obtain that can be simplified in In order to get a symmetric problem, we multiply the previous equation by θ 2 with r = θ 2 ev N,k + H v k + θM −1 hf ext,k+θ . The last line of the inclusion is finally rewritten as Since the inclusion involves a cone, the problem (65) amounts to solving the following affine variational inequality Equivalent convex quadratic optimization problem The goal of the section is to show that the affine variational inequality (86) is equivalent to a convex quadratic problem if Assumptions 1 hold. The matrix Q can be written as 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 (86) is equivalent to solve the following convex quadratic optimization problem The well-posedness of the problem (89) 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 (89) has a unique solution (σ, a, p N ) if the set C is a non empty convex set.
Under Assumptions 1 and 3, the matrix W = θ 2 H ⊤ M −1 H is a symmetric positive definite matrix. The matrix Q can be written as For a sufficiently small time step h, the matrix Q is then a positive definite matrix. Since the optimization problem (89) is strictly convex, a unique solution exists if the set C is a non-empty convex set.
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.

Numerical methods
Since the problem (89) 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: where f is a smooth vector-valued function with non-vanishing gradients at the solution. The optimization problem (89) takes the form min σ,a,pN 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;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 σ.
where s is the deviatoric part of the tensor σ, and the hydrostatic stress σ h defined by The Huber-Von Mises criterion is With an alternative expression of the J 2 invariant, another expression of the Huber-Von Mises criterion is in the original stress tensor: Using the Voigt notation for the stress σ = σ 11 σ 22 σ 33 σ 23 σ 13 σ 12 ⊤ , the J 2 invariant can be expressed as a quadratic function 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 min σ,a,pN 3σ ⊤ g Jσ g ⩾ κ 2 for g ∈ 1, n e,g p N ⩾ 0.

(99)
The case of Drucker-Prager plasticity. The Drucker-Prager plasticity involves the following definition of the convex set C as where ∥s∥ = J 2 , k d is a constant that depends on the dimension d, c is the cohesion and φ the friction angle. It is well know that the set K σ is a translated second order cone that cannot be straightforwardly written with a smooth function f as in (91). 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 min σ,a,pN 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 (σ h , s) g ∈ K σ defines the following flow rule where e = ε − 1 3 ε h I and ε h = tr ε. 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 where θ is the dilatancy angle ranging from 0 to φ. 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 (89) or (92) amounts to finding a solution of the optimality conditions given by the variational inequality (86). 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): In the case of (86), we obtain the following equation for ρ > 0 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 Saxcé, 1994). The method is robust by rather slow, and a second order method is preferred. To this aim, the equation R(x) = 0 in (104) is solved with a semismooth Newton method as 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.

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: where E(t) is total mechanical energy E(t) = T (t) + Ψ(t). After the space discretization by a finite element method, the energy balance takes the form where the space-discretized energies are given by with the convention that the kinetic energy is also a right continuous function of bounded variations, i.e, T(t) = T + (t). The dissipation and the power of external forces are Since in a discrete system, impacts may occur, the power of the reaction impulse is given by Let us now compute the variation of the total mechanical energy over a time-step, starting by the kinetic energy Using the relation ) and the first equation in (65), we obtain For the potential energy, we obtain Using the approximation of works of external forces and the dissipated work by the θ-method as and W k+1 p k and finally, an approximation of the work dissipated by the percussion (see Acary (2016)) the increment of total energy is then given by The following proposition summarizes the properties of the discrete energy balance (118).
Proposition 5. The main conclusions that can be drawn for the energy dissipation of the scheme are 1. For θ = 1 2 , the time-stepping scheme satisfies the approximation of the discrete energy balance without introducing artificial dissipation: (119) 2. The dissipated work due to plasticity is positive and for θ ⩽ 1 1 + e , the dissipated work due to impact is also positive.

For
A result in (Acary, 2016) shows that Furthermore, the inclusion implies that W p ⩾ 0. The inequality (120) ensures the practical unconditional stability of the scheme.
Remark 7. Naturally, the θ-method has some shortcomings such as excessive numerical dissipation in dynamics for θ > 1/2 which makes the scheme more suitable for quasi-static evolutions. For θ = 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-α scheme (Brüls et al., 2014(Brüls et al., , 2018Chen, Acary, et al., 2013) is much better suited. In the article, we have decided to keep θ-method for simplicity, but nothing prevents the use of more advanced schemes.

Numerical applications
In this section, a serie 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). 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 : 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.

Piecewise linear yield surface
The case of a feasible set C defined by linear relations When the set C is finite represented as in (91) and assuming that f is sufficiently smooth with nonvanishing gradients, the normal cone to C can expressed as Then the plastic flow rule is expressed as the classical flow rule for plasticity with the complementarity condition where λ is the vector of plastic multipliers. If the yield function f is linear, the gradients are given by constant matrices denoted by and we have With linear yield functions, the variational inequality (86) can be reduced to a Linear Complementarity Problem (LCP). To this aim, we express the stress tensor from (82) as Substituting (129) in (84) and (128) yields and The problem in terms of λ and p N,k+1 can be defined in the form of the follwing LCP 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). 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 ⊤ σ , J ⊤ a , κ Y are defined in the present case by 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: Both kinematic and isotropic hardenings are considered leading to the definition of the force vector a g , expressed at each Gauss point as 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 ⩾ 0 is associated with isotropic hardening. The vector of plastic multipliers λ g at Gauss point g is denoted by 4.2 Example 1 : Cantilever beam subjected to quasi-static displacements 1 2 3 4 5 6 7 8ū x u y L elements end points elements mid points Gauss points y Figure 4: Example 1: Cantilever beam discretized in 4 elements subjected to quasi-static displacements u x andū y at its free end.
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ū x andū y at its free end. The properties of the structure are summarized in Table 1.
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 θ = 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 the beam free end. The axial and transversal components of the velocity are 1×10 −4 m/s and 4×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 × 10 −3 m/s until t = 2.5 s. Then, the transversal velocity is cancelled and a constant axial velocity, equal to 2 × 10 −4 m/s, is applied (see Figure 5). 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.84 s. Once the yield criterion is reached, the plastic multiplier λ 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.
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 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 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.
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.4 s (Figure 14a), the plastic multipliers λ 1,3 and λ 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.5 s - Figure 12), the loading path is modified. Thus, only one inequality constraint is active and the multiplier λ 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). : 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.   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.  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. λ 4,3 λ 3,3 λ 2,3 λ 1,3 Figure 11: Example 1, case 2: time evolutions of the plastic multipliers λ 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. 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.  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.

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 θ parameter).
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, θ was set at 1/2 to prevent artificial energy dissipation and the time step was set at h = 10 −4 s.

Cross section area
A (m 2 ) 3.4 × 10 −3 Mass of the sphere m sphere (kg) 250 Initial height of the sphere q y sphere (t = 0) (m) 1

Moment of inertia
Restitution coefficient e 0 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 k+1 ext 0 , i.e. the sum of the kinetic energy T k+1 and the free energy Ψ k+1 , is always equal to the sum of elastic mechanical energy E k+1 , plastic dissipation W k+1 p 0 and contact dissipation W k+1 c 0 . Indeed, setting θ = 1/2 prevents artificial energy dissipation. Just before impact, the work of external forces W k+1 ext 0 is equal to the kinetic energy of the sphere. During the simulation, W k+1 ext 0 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 k+1 p 0 ) and to contact (W k+1 c 0 ), 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 k+1 c 0 remains constant while W k+1 p 0 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 θ = 1/2. 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. n el = 10 n el = 50 n el = 100 n el = 500 n el = 1000 n el = 5000 n el = 10000 n el = 50000 n el = 100000

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 weigth of 804kg). The structure is clamped at its two extremities. The projectile (mass : m p = 1000kg) impacts the structure at 0.67m from the ground with a horizontal velocity v i = 6 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 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 θ = 1/2 to prevent artificial energy dissipation and e = 0 to maximize energy dissipation at contact.   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 −3 s (Figure 24a), the projectile bounces on the structure with a velocity after rebound equal to −4.24 m s −1 (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 −1 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 −1 , 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 −1 , 6 m s −1 and 8 m s −1 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    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. 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.

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's 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 θ parameter. For θ = 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 (67) to the work of Mielke (2003Mielke ( , 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 nonlinear 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.
for the consistent mass matrix, and for the nodal internal, external and contact forces, we obtain the contribution to equation of motion of the element e as M ev (t) + f int,e (t) = f ext,e (t) + r e (t).
The discrete equation are obtained by summing the contributions of the elements yielding √ ω e,g B ⊤ e (x g ) √ ω e,g σ e (x g ) where x e,g ∈ IR d , g ∈ 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 ε e,g = B e,g u g , with ε e,g = ε e (x g ) and B e,g = B e (x g ) The strains at Gauss points in an element and in the structure are collected scaled by their weights as follows ε e = col( √ w e,g ε e,g , g ∈ 1, n e,g ) and ε = col(ε e , e ∈ 1, n el ) and, doing in the same way for the matrix B e,g , we get ε e = B e u and ε = Bu, with B e = col( √ w e,g B e,g , g ∈ 1, n e,g ) and B = col(B e , e ∈ 1, n el ).
In the same way, the stresses, evaluated at Gauss points of an element are gathered as follows σ e,g = σ e (x g ), σ e = col( √ w e,g σ e,g , g ∈ 1, n e,g ) and σ = col(σ e , e ∈ 1, n el ) The matrix B appears as the discrete divergence operator, also called the equilibrium matrix in the classical force method.

C Beam element
The elementary form function matrix is The derivative matrix of N is then: The elementary consistent mass matrix associated with the element is : where L e is the element length, A is the cross section area and ρ is the material density. The elementary elasticity matrix C e , which relates stresses σ e and elastic strains ε e e , is : The elementary hardening matrix D e , which relates internal forces a e and hardening parameters α e , is: where H ki and H is are the kinematic and isotropic hardening moduli, respectively.