An acoustic Riemann solver for large strain computational contact dynamics

This article presents a vertex‐centered finite volume algorithm for the explicit dynamic analysis of large strain contact problems. The methodology exploits the use of a system of first order conservation equations written in terms of the linear momentum and a triplet of geometric deformation measures (comprising the deformation gradient tensor, its co‐factor, and its Jacobian) together with their associated jump conditions. The latter can be used to derive several dynamic contact models ensuring the preservation of hyperbolic characteristic structure across solution discontinuities at the contact interface, a clear advantage over the standard quasi‐static contact models where the influence of inertial effects at the contact interface is completely neglected. Taking advantage of the conservative nature of the formalism, both kinetic (traction) and kinematic (velocity) contact interface conditions are explicitly enforced at the fluxes through the use of appropriate jump conditions. Specifically, the kinetic condition is enforced in the usual linear momentum equation, whereas the kinematic condition can now be easily enforced in the geometric conservation equations without requiring a computationally demanding iterative algorithm. Additionally, a total variation diminishing shock capturing technique can be suitably incorporated in order to improve dramatically the performance of the algorithm at the vicinity of shocks. Moreover, and to guarantee stability from the spatial discretization standpoint, global entropy production is demonstrated through the satisfaction of semi‐discrete version of the classical Coleman–Noll procedure expressed in terms of the time rate of the so‐called Hamiltonian energy of the system. Finally, a series of numerical examples is examined in order to assess the performance and applicability of the algorithm suitably implemented in OpenFOAM. The knowledge of the potential contact loci between contact interfaces is assumed to be known a priori.


INTRODUCTION
In computational mechanics, the numerical modeling of contact and impact phenomena has been a major field of interest in industry, including numerous applications such as vehicle crash testing, 1 prototype testing, or manufacture processes. 2,3 These very complex problems are typically characterized by highly nonlinear deformation behavior with accompanying non-smooth response (or shocks) caused by transitions between various contact modes such as separation-to-contact, stick-to-slip, slip/stick-to-separation. Such problems must be solved by ensuring the satisfaction of linear momentum conservation equation (complemented by appropriate initial and boundary conditions) for each body individually, whilst at the same time enforcing the additional set of (kinematic and kinetic) contact interface conditions, that govern the interaction of these bodies. When considering a model with frictionless contact, these interface conditions act to prevent interpenetration of the bodies (kinematic condition) and to insure compressive interaction normal to the interface (kinetic condition). One challenging aspect, is that impenetrability cannot be expressed as an evolution (or algebraic) equation and so requires special numerical treatment. The most common techniques addressing this issue include penalty method, Lagrange multiplier method, or a combination of both. In the penalty method, 4 the impenetrability constraint is enforced as a penalty normal traction along the contact surface. The disadvantage of the penalty approach is that the enforcement of the impenetrability condition is only approximate and its effectiveness depends on the selection of the user-defined penalty parameters. If the value of the penalty parameter is too small, unpredictable amount of interpenetration would be observed. However, the penalty parameter cannot be arbitrarily large, as this can generate ill-conditioned systems that may require extremely small time steps for stability. 5,6 The correct choice for this parameter is key to success of the algorithm. For the Lagrange multiplier method, 7,8 the multipliers must be approximated and solved at the contact interface with the constraint such that the normal component of the traction must be compressive. The disadvantage of the method is that it requires the construction of a separate independent mesh and also requires the introduction of additional regularization techniques necessary in obtaining robust solutions. Such regularization procedures are usually ad-hoc and are not motivated by physical arguments. A popular example of regularization is the addition of von Neumann's artificial viscosity 9,10 to the Euler fluid equations to smear shocks over several computational cells. In absence of this artificial viscosity, central difference solutions to the Euler equation in the vicinity of shocks are oscillatory, eventually leading to the breakdown of a numerical scheme.
One of the earliest attempts at enforcing contact interface conditions via the physically-motivated jump conditions that derive from the linear momentum conservation equation and kinematic compatibility can be traced back to the work of Abedi and Haber. 11 In particular, a space-time discontinuous Galerkin finite element method for elasto-dynamic contact was presented. The presented examples were restricted to the case of small strain linear elasticity in two dimensions. Moreover, it is not yet clear if the overall finite element algorithm would satisfy the classical Coleman-Noll procedure in order to guarantee the production of non-negative entropy. On another front, some interesting works have also been explored using the computational fluid dynamics (CFD) platform "OpenFOAM" via the use of displacement-based finite volume discretization, 12,13 with special attention paid to the quasi-static simulation based on the use of penalty method 14,15 and lubricated contact models. 16 Aiming to resolve the shortcomings described above, the main goal of this article is to explore the solution of contact dynamics utilizing a set of first order conservation equations, [17][18][19][20] combined with the associated jump conditions across moving shocks * . [21][22][23][24][25][26][27][28][29][30][31][32][33] Building upon previous work developed by the authors, 34,35 a mixed methodology is presented in the form of a system of hyperbolic conservation laws, where the linear momentum and the minors of deformation (the deformation gradient tensor, its co-factor and its Jacobian) are regarded as the main conservation variables of this mixed approach. Taking advantage of this formalism, appropriate kinetic and kinematic contact interface conditions can be suitably enforced at the boundary fluxes of the underlying hyperbolic system by means of the Rankine-Hugoniot jump conditions. For instance, in the case of frictionless contact, the normal traction is enforced in the standard manner, that is, in the boundary fluxes of linear momentum conservation equation. One crucial advantage of solving the geometric conservation equations in this context is such that we now have the luxury of explicitly enforcing the normal component of the contact velocity in the associated boundary fluxes. On the other hand, upon the use of Rankine-Hugoniot jump conditions, 36 we can then naturally derive a series of dynamic contact models typically required in the simulation of contact problems. The objective of this article is to present a complete set of continuum Riemann-based type solutions for contact and separation (derived based on a system of first order conservation laws) assuming a priori knowledge of the potential contact loci. Physically, Riemann solutions describe correct fluxes in the form of traction and velocity at the contact interface. In linear elasticity, we can show that by enforcing appropriate boundary fluxes at the point of contact through the use of jump conditions would lead to exact energy transfer, provided the shock wave travels at the speed of sound.
From a spatial discretization point of view, a vertex based finite volume algorithm 34 in conjunction with (piecewise) linear reconstruction is employed. Additionally, a shock capturing technique 37 can also be incorporated in order to dramatically improve the resolution of the field variables at the vicinity of shocks. No ad-hoc algorithmic regularization procedures are needed. Insofar as contact-impact introduces discontinuities in the solution, the use of explicit time integrators is preferred (see Chapter 10 in Reference 5) as neither linearization nor a Newton's method is required. With this in mind, from a temporal discretization standpoint, we use the explicit type of two-stage Runge-Kutta time integrator. A crucial aspect that requires special attention is that of the stability of the overall algorithm. 38 This can be demonstrated by monitoring over time the Hamiltonian energy of the system, ensuring the production of entropy throughout the entire simulation. The overall methodology is shown to be capable of handling contact-impact problems without excessive spurious modes, even in the case of nearly incompressible elasticity and elasto-plasticity. Examples presented in the article are specifically chosen in order to illustrate the capability of the proposed framework addressing spurious oscillations in problems with shocks (or spatial jumps) without resorting to ad-hoc dissipative correction. Another contribution of the current work is to carry out its implementation into the OpenFOAM platform, widely accepted these days by industry.
The article is broken-down into the following sections. Section 2 starts by summarizing the total Lagrangian formulation of the conservation laws to be solved, comprising the linear momentum and the three geometric conservation laws. Section 3 provides the exact solution of the simple one dimensional two-bar impact derived from the associated jump conditions in linear elasticity. This leads to exact energy transfer from one bar to the other after contact without energy loss. Motivated by this, the section continues deriving a set of interface contact conditions (velocity and stress) applicable to multi-dimensional contact problems. Section 4 presents the second law of thermodynamics written in terms of the so-called Hamiltonian free energy. Section 5 describes the computational methodology of the vertex centered finite volume method. A proof of entropy production is included as a necessary condition for stability at the semi-discrete level. Section 6 includes the algorithmic flowchart of the resulting numerical scheme where special attention is paid to the procedure in addressing non-matching mesh interface. Section 7 presents a set of numerical examples to assess the accuracy and stability of the computational framework, with detailed comparison with other verified finite element code such as Abaqus. Section 8 presents some concluding remarks.

FIRST ORDER HYPERBOLIC SYSTEM FOR SOLID DYNAMICS
Consider the three dimensional deformation of an isothermal body of material density R moving from its initial undeformed configuration Ω V , with boundary Ω V defined by an outward unit normal N, to a current deformed configuration Ω v (t) at time t, with boundary Ω v (t) defined by an outward unit normal n. The time dependent motion (X, t) of the body can be described by the following system of total Lagrangian global conservation laws 35,39-52 with the surface flux vector being defined as  N = ∑ 3 I=1  I N I . Here,  is the vector of conservation variables,  I is the flux vector at Ith material direction and  is the source term. Their components are with the Cartesian material coordinate basis being defined In terms of notations, p = R v is the linear momentum per unit material volume, v represents the velocity field, f R is the body force per unit material volume, {F, H, J} are the triplet deformation measures representing deformation gradient tensor, its co-factor and its Jacobian, P represents the first Piola-Kirchhoff stress tensor. The symbol × represents the tensor cross product between vectors and/or second order tensors, see References 41,53,and 54. Given the fact that the above system (1) has more equations than needed, suitable compatibility conditions (also known as involutions 45,46,55 ) are necessary, namely CURL and DIV represent the material curl and divergence operators carried out with respect to the material configuration. For smooth functions, expression (1) is equivalent to the set of first order local differential equations described as The above local form implies that the variables describing the state of the solid in motion (such as, velocity v and stresses P) are continuous functions throughout the solid. In other words, it is always possible to find their spatial derivatives as required by the divergence operators that appear in Equation (5). This is indeed usually the case but situations may arise when these variables experience sudden jumps in value, that is, they become discontinuous across surfaces which move across the body. These jumps are known as shocks and are the result of sudden physical phenomena such as contact-impact problems.
To account for shock phenomena, the integral Equation (1) also leads to the following jump conditions across a discontinuity surface with normal N propagating with speed U, that is These jump conditions are sometimes referred to as Rankine-Hugoniot equations 56,57 describing the behavior of a material across a shock. These conditions can then be particularized for the set of conservation variables considered in this article, namely the linear momentum and the triplet of deformation measures Here, ⟦•⟧ = [•] + − [•] − denotes the jump operator between the right and left states of a discontinuous surface. For the particular case of a reversible process, the closure of system (1) requires the introduction of a suitable constitutive law relating the stress tensor P with the triplet of geometric strain measures {F, H, J}, obeying the principle of objectivity 5,58,59 and thermodynamic consistency (via the Coleman-Noll procedure 60 ). In this work, a Mooney-Rivlin model is employed and is summarized in Remark 2 for completeness. Finally, for a complete definition of the initial boundary value problem, initial and boundary (essential and natural) conditions must be specified as appropriate.
Remark 1. For the conservation of mass, the material density at the reference configuration cannot be a function of time, that is R t = 0. 61 This implies that R is given by the initial conditions of the solids and it remains constant throughout the motion, and does not need to be considered as part of the unknowns in system (1) to be solved in time. Since there is no possible flow of mass across the physical interface (meaning, the associated normal flux vector vanishes), the jump condition associated with the mass conservation becomes U⟦ R ⟧ = 0. This intrinsically implies that the material density at the reference configuration must be continuous across a shock. Equation (7a) thus reduces to where , , and (bulk modulus) are positive material parameters. Appropriate values for the material parameters and can be defined in terms of the shear modulus , that is, 2 + 3 It is now possible to express the first Piola Kirchhoff stress tensor as 41 where the conjugate stresses { F , H , Σ J } are defined by and For = 0, the Mooney-Rivlin model described above (9) degenerates into the so-called nearly incompressible neo-Hookean model. In order to model irrecoverable plastic behavior, the standard rate-independent von-Mises plasticity model 58 with isotropic hardening is used and the basic structure was already summarized in Algorithm 1 in Reference 45.
Remark 3. It is often necessary to obtain expressions for the symmetric Kirchhoff (or Cauchy) stress tensor since it is needed either to express plasticity models or to display the solution results. Such expressions can be easily obtained from the following standard relationship between these tensors 58 To achieve this, substitution of (10) into (13) for P, and following the procedure described in Reference 41, gives the resulting expression for the Kirchhoff (or Cauchy) stress where

Motivation: The local one-dimensional contact solution
In order to motivate the more complex contact-impact solutions (e.g., stick-slip-separation transition) developed in this section, consider first a simple one-dimensional case comprising two bars, as illustrated in Figure 1A, where the left bar is travelling with a given velocity v 0 impacts the right bar which is at rest. When the contact between two bars takes place, the resulting contact-impact motion is governed by a reduced set of (one-dimensional) jump conditions described in system (7) as is the gap separation between two bar at t 0 . Left column represents the velocity profile v x and right column represents the stress profile P xX (not the traction).
For ease of understanding, both bars are assumed to be made of the same linear elastic material defined as P xX = ( + 2 ) (F xX − 1). 58 With this linear model, and noting that ⟦P xX ⟧ = ( + 2 ) ⟦F xX ⟧, we can then obtain the shock wave speed c p by substituting expression (16b) into (16a). Expression (16a) after some algebra becomes and which after rearranging gives † This corresponds to the speed of the sound wave in the bar which can be obtained considering the classical wave propagation theory. 61 Once the shock wave speed is determined, attention is now focused on the evaluation of the velocity (kinematic) and traction (kinetic) at the contact-impact scenario. The same evaluation procedure would also be repeated when considering the separation process. When contact is made between two points of the bar, shock waves are generated and travel in opposite directions along each bar (from the contact point to the free end of the bar) as shown in Figure 1A-C. When such a compressive stress wave reaches the end of a bar, wave reflection occurs (see Figure 1D). The reflective wave varies depending on the actual physical boundary of the problem under consideration. In the case of a free end (i.e., traction is zero and particle velocity is doubled), the reflected wave becomes a tensile stress wave which is an inverted image of the compressive stress wave. The frequency and amplitude of the velocity wave however in this case remains unchanged in reflection. Finally, as soon as the tensile stress wave arrives at the contact point, both bars would undergo the separation process. The above procedure describing the wave evolution for a two-bar impact in one dimension is graphically represented in Figure 1 for clarity.
Let us now focus on the mathematical solutions of contact-impact state between two bars. Upon contact, both instantaneous velocity v C x and traction t C x are obtained by applying the jump condition (16a) between the values of variables before and after the impact through appropriate initial conditions, to give The first equation (19a) corresponds to the jump relation between the left bar (travelling with a given speed v − x = v 0 ) and the contact point, whereas the second equation (19b) refers to the jump between the right bar (at rest v + x = 0) and the contact point. Additionally, both ends of the bars are traction-free right before contact, which implies that t − x = t + x = 0. Solving the above system (19) analytically gives the common (or continuous) velocity and traction at the contact point for both bars This is generally known as contact-stick mode in one dimension. Following the same procedure described above, it is also possible to determine the release velocity for each of the bars right after separation. In this separation mode, the release traction must be zero ensuring the traction-free compatibility condition. Focusing first on the left bar, the release velocity can be achieved by introducing the jump condition (16a) between the values of variables before and after separation to give Suitable conditions for velocity and traction right before separation must be enforced. In this specific demonstration, and referring to Figure 1E, we use the values of v − x = t − x = 0 at the contact point prior to separation. With these at hand, the associated release velocity in expression (21) becomes null, shown as below Analogously, we can now repeat the same demonstration for the right bar. The corresponding jump condition in relation to the right bar follows Using the exact same value of the contact traction as for the left bar (i.e., t + x = 0) and the velocity to be v + x = v 0 , the above expression yields v C,+ Observing the fact that the release velocity for the left bar is null (22) and for the right bar is v 0 (24), this represents a complete transfer of kinetic energy (and also total energy) from the left bar to the right bar at the point of contact without energy loss. This is illustrated in Figure 1A,F. From the point of view of hyperbolic differential equations, the sudden impact between two bars results in a Riemann problem in each bar with a simple analytical solution. The boundary fluxes at the point of impact, namely the velocity and traction, must be compatible with the jump conditions. In linear elasticity, this compatibility leads to exact energy balance ‡ since the shocks wave travels at the speed of sound. From a mathematical standpoint, the elastodynamic contact problem is well posed without the need to introduce artificial ad-hoc dissipative effects, provided that the first order conservation equations (1) together with the associated jump conditions (7) are used. The sections that follow will conceptually extend the above simple local contact solution to three dimensional local solution with a possibility to consider bi-material between contact.

Extension to general contact procedure
It is worthwhile recalling the general solution process for the contact algorithm (for instance, stick-slip-separation transition) that would ensure the satisfaction of Karush-Kuhn-Tucker condition. 5,11 To begin with, it is instructive to first determine the trial solution assuming contact-stick mode. Such trial solution is used as a criteria in the subsequent development, to check if two bodies are in contact or separate. The contact between two bodies will only take place if the following two conditions hold The first (kinetic) condition ensures that the normal component of the trial contact traction t C,trial n must be in compression, whereas the second (kinematic) condition ensures two bodies are in contact, that is the normal separation between points of contact is zero. We next need to examine the nature of the contact motion, whether they are in stick mode or slip mode, depending on the tangential friction introduced in the model. To achieve this, it is possible to introduce a slip criterion Φ accounting for the difference between the value of a trial tangential traction vector t C,trial t and a tangential friction arising from isotropic Coulomb friction, 5 that is with the magnitude (or norm) of a vector being defined as || in the above equation represents the positive part of the scalar value. The value of the slip criterion determines the transition between contact-stick and contact-slip modes. If Φ ≤ 0, contact-stick conditions hold given the fact that the computed tangential force does not exceed the Coulomb limit. Otherwise, when the value of Φ > 0, we then accept the solution to be in contact-slip mode.
Finally, the transition from either contact-slip or contact-stick to separation would happen when t C,trial n ≥ 0 (that is, the violation of the kinetic condition), even with the value of = 0. The overall procedure described above is summarized in Algorithm 1. This however would require the evaluation of contact traction and velocity associated with various dynamic contact models involved, and which will be presented in the following sections.

Contact-stick condition
Motivated by the one-dimensional problem illustrated in Section 3.1, we now extend the concept to multi-dimensions by postulating that impact between two bodies travelling at different speeds leads to a common velocity and traction at the point of contact as shown in Figure 2. The normal components of the velocity and traction at the point of contact are defined as Both of these values are likely to be different for the left and right bodies right before contact and are therefore denoted as v − n and t − n for the left body and v + n and t + n for the right body. The common values after contact are denoted by v C n and t C n . Note first that the impact will generate two types of shock waves travelling from the contact point into each of the two bodies.
In the case of frictionless contact, the generated shock waves will travel with volumetric speed U p . The evaluation of the common contact velocity and traction vectors is governed by the jump conditions across the two shocks, obtained by applying Equation (8) on each body as follows Algorithm 1. General procedure for stick-to-slip-to-separation transition if = 0 then Obtain trial contact-stick traction: t C,trial = t C (34a) (see Section 3.2.1) Determine the normal contact traction: t C,trial Note that the negative sign in front of N − and N + are necessary as the shocks propagate into the body in directions opposite to N − and N + . Multiplying the above expressions by a unique normal vector gives The difference in sign between expressions (29a) and (29b) is because n is normal to the surface of the left body and hence t + n = −n ⋅ (P + N + ) and t C n = −n ⋅ P C N + , whereas t − n = n ⋅ (P − N − ) and t C n = n ⋅ (P C N − ). Expressions (29) represent a system of two equations for four unknowns, namely t C n and v C n (expressed in terms of the left and right normal tractions and velocity before the impact) and also {U − p , U + p } the speeds of the shocks after impact. Unfortunately, the shock speeds in Equation (29) are in general also a function of the unknowns t C n and v C n rendering the system of equations highly nonlinear. A much simpler case emerges when the speed of the shock after contact is assumed to be equal to the speed of sound c p (18) (derived on the basis of the linear elastic model), which only depends on the material properties under consideration. This is usually referred to as an acoustic Riemann solver 56,57 widely known in the field of CFD. Doing this will give closed-form expressions for the normal components of velocity v C n and traction t C n as In the situation where friction is present to a sufficient degree to prevent relative sliding, similar common tangential components of the velocity and traction can be derived. Consequently shear shocks are also generated and are, again, assumed to be identical to the simple shear wave speed obtained via linear elasticity 45 If v t and t t are defined as a similar derivation for the common tangential traction and velocity vectors gives where c − s and c + s are the left and right body shear shock speeds. Finally, the complete common velocity and traction vectors at the contact point can be combined This is typically known as contact-stick mode. Numerically, expressions (34) can be viewed as the summation of the average states (unstable) and the associated upwinding stabilization terms depending on the jumps. This has been extensively exploited by the authors in developing stabilized methods with the objective to improve the numerical solutions by alleviating unwanted spurious hour-glassing and pressure instabilities. 34,43,45,[47][48][49] Provided the interface conditions (25) and also the slip function Φ ≤ 0 (26) hold, we accept the contact-stick solution (34) as the actual local solution for contact. Otherwise, we must investigate other possible contact models such as contact-slip or separation. This will be presented in the next section.
Remark 4. It is also useful to consider the case where the jump in traction (or the first Piola Kirchhoff stress P (10)) is dominated by the jump in the pressure component of the stress (which in this case is related to Σ J ), whilst the rest of the components of the stress { F , H } can be neglected. This is indeed the case when attempting to model problems with predominant nearly incompressible behavior. Use of (10) in conjunction with the Nanson's rule H Ave N = Λ H n (where Λ H is the ratio between the current area and the undeformed area), enables the jump in the traction vector to be Neglecting the stress components F and H implies that ⟦ F + H × F⟧N = 0. The jump in the traction vector in the normal direction can now be derived by multiplying (35c) with a normal vector n to yield It is also interesting to notice that the jump in the tangential component of the traction vanishes. This is easily shown as below by making use of expressions (32b), (36), and (35c).
Remark 5. When considering the exact same material properties on the left and right sides at a point of contact, the density and the shock wave speeds are identical and constant for both sides, namely − Enforcing these conditions in (30) and (33) especially in the case of nearly incompressible materials (see . (38b) Remark 6. Consider a contact system where the material on the right side is much stiffer than the material on the left side. Under this circumstance, both pressure and shear shock wave speeds on the stiffer material are approximated to be c + p ≈ ∞ and c + s ≈ ∞, which, upon substitution into Equations (30) and (33), gives Observe the fact that only the velocity of the stiffer side, that is v + , enters the solutions. These solutions indeed coincide with the Dirichlet boundary conditions already discussed in Reference 45, where the velocity v + is prescribed on the boundary of the domain. For instance, when considering no-slip wall boundary condition, the values of v + t and v + n are set to zero.

Contact-slip conditions
When the magnitude of the tangential traction described in (33) exceeds the Coulomb friction limit, that is Φ > 0 (26), relative sliding between two surfaces is then allowed. This phenomenon is known as contact-slip mode. In this mode, only normal component of the velocity is continuous across the contact surface between the left and right bodies. The tangential components of the velocity however may suffer jumps. Mathematically, the complete velocity field associated with slip mode for both the left and right surfaces are postulated as with v C n being defined in (30a). The remaining variables to be determined are the respective tangential velocity components {v C,− t , v C,+ t }. In order to achieve this, it is instructive to consider the slip traction vector to be With regard to the first term on the right hand side of the above equation, the normal component of the contact traction t C n must be in compression (i.e., its value must be strictly negative) and its expression remains exactly the same as the one described in contact-stick mode (30b). The second term t B t represents the tangential frictional traction arising from the Coulomb model of friction and is in the direction n ⊥ opposing the motion predicted using the tangential traction (33) in stick mode. Use of (41) in conjunction with (28), enables the tangential components of slip velocity vectors to become Re-arranging the above expressions render It is useful to notice that frictionless contact condition (also known as symmetric condition) can be easily recovered by enforcing the value of frictional coefficient k = 0 in (41), which in turn implies that t B t = 0.

Separation conditions
In the current work, we assume homogeneous prescribed traction (that is, traction-free conditions) in separation mode, but non-vanishing prescribed tractions due to viscous fluid loading are also possible. This reveals the fact that the traction after separation must ensure traction-free compatibility conditions, namely With this, we are now in a position to proceed with the evaluation of release (or post-separation) velocity on respective surfaces by re-applying Equation (8) between the values of variables before and after separation, repeated below again for convenience The normal components of the release velocities are obtained by multiplying expressions above with the normal vector n and which, after rearranging, becomes In line with linear elasticity theory, above expressions coincide with the expressions shown in (22) and (24). A similar derivation for the tangential velocity vectors can now follow Combining (47) and (48) enables the release velocities to be expressed as It is worth noticing that all the velocity components described in (49) are generally distinct and independent on opposing sides of the contact surface.

SECOND LAW OF THERMODYNAMICS
In order to pave the way for the proof of entropy production presented in a subsequent section, it is useful to introduce the notion of Hamiltonian (X, t). 61,63 For the isothermal case, this indeed can be understood as a generalized convex entropy function of the system of conservation equations (1), coinciding with the definition of total energy per unit of undeformed volume. Specifically, the Hamiltonian  is defined by which represents the summation of the kinetic energy per unit of undeformed volume (i.e., the first term on the right hand side of (50)) and the internal energy  expressed in terms of the triplet deformation measures {F, H, J} and a set of state variables [64][65][66] (such as plastic deformation or similar) collected in the form of a tensor . Note here that (X, t) and (p, F, H, J, ) represent alternative functional representations of the same quantity. It is instructive to revisit the second law of thermodynamics when written in terms of the Hamiltonian. Taking the derivatives of (50) with respect to its arguments, the time rate of the Hamiltonian for one of the bodies involving contact is obtained via the chain rule as follows where, Equations (5) and (10) have been substituted in the third and fifth lines of (51), respectively. Subsequently, we can substitute the linear momentum conservation Equation (1) into (51) to give Recalling By performing integration by parts of the DIV term in expression (53), and after some re-arrangement, it renders whereΠ ext denotes the power introduced by external forces, defined aṡ Here, Γ represents the boundary faces on contact region and Ω V ⧵ Γ represents the remaining boundary faces that are not in contact. In the above expression, the first term on the right hand side represents external force acting on a body, the second term represents the non-contact boundary forces obtained via the enforcement of standard Neumann or Dirichlet boundary conditions, and the third term represents the contact boundary forces describing appropriately the contact-impact phenomenon. Such contact boundary contributions are suitably enforced by solving a Riemann-like problem as already presented in Section 3.
Consider the case of elasto-plasticity 47,52,58 where the elastic energy is expressed in terms of elastic left Cauchy-Green tensor b e = FC −1 p F T , thus in this case the internal state variable is indeed the inverse of the plastic right Cauchy Green tensor, that is = C −1 p . With this, the rate of plastic dissipationḊ described in (54) becomeṡ In view of the fact that the rate of plastic straiṅp has been defined as the work conjugate to the von Mises equivalent stress , 58 equation above can be alternatively expressed as 58 where ′ represents the deviatoric component of the Kirchhoff stress. Noticing that in the above expression the rate of dissipation is always non-negative, that isḊ ≥ 0, Equation (54) can be transformed into the following inequality which represents a valid expression for the second law of thermodynamics. 59 Satisfaction of inequality (58) is a necessary ab initio condition to ensure stability, otherwise referred to as the Coleman-Noll procedure. 34 This key concept will be further exploited in Section 5.2 at a semi-discrete level.

Semi-discrete formulation for dynamic contact
The vertex centered finite volume spatial discretization presented in this work requires the generation of a median dual mesh 34,35,50 for the definition of control volumes (see Figure 3). With this in mind, expression (1) can now be spatially discretized over an undeformed control volume Ω a V , to give Here,  a and  a are the average values of both the conservation variables and source term vector within the control volume, respectively. Moreover, the surface integral of (59) is approximated by means of appropriate numerical fluxes, resulting in § where b ∈ Λ a represents the set of neighboring control volumes b associated with the control volume a and C , = A , 3 N , represents the (tributary) boundary area vector. For a given edge connecting nodes a and b, the mean undeformed area F I G U R E 3 Dual mesh of (A) an interior node and (B) a boundary node using the medial dual approach in two dimensional triangular mesh.
vector C ab satisfies the reciprocal relation, that is C ab = −C ba . The terms within the parenthesis in (61) correspond to the evaluation of the control volume internal interface flux  I N ab , non-contact boundary fluxes  B a and contact boundary fluxes  C a . This evaluation is comprised of a summation over edges (first term in the parenthesis), a summation over non-contact boundary faces (second term in the parenthesis) and a summation over contact boundary faces (third term in the parenthesis). The internal interface flux  I N ab =  I N ab ( − ab ,  + ab , N ab ) must be evaluated on the basis of the contact-stick condition (see Section 3.2.1), which depends on the reconstructed states at both sides of the mid-edge of ab, namely  − ab and  + ab . The non-contact boundary flux  B a is enforced through either Neumann or Dirichlet boundaries and the contact flux  C a is determined following strictly the contact procedure presented in Algorithm 1 (obeying appropriate contact-impact physics). Notice that in (61), C = 0 when the boundary face is not in contact.
It is worth noticing that Equation (61) would only lead to a first order solution in space 45 provided that  − ab and  + ab are modeled following a piecewise constant representation. For instance,  − ab =  a and  + ab =  b , thus leading to excessive numerical dissipation in the solution. The physics of the problem can no longer be captured accurately unless excessively fine meshes are used, which is clearly undesirable especially for large scale problems in practice. To overcome this drawback, and to guarantee second order accuracy in space, a suitable linear reconstruction procedure is used. A detailed discussion of this reconstruction procedure has already been presented in References 45,48,and 47. Expression (61) is now particularized for each individual component of  , yielding Here In ensuring discrete satisfaction of the involutions (4), and following the work of Reference 34, one viable option is to approximate the updates of F (62b) and H (62c) using central difference approximations by neglecting the jump in traction in (30) and (33) (or jump in Σ J (38a)) for v I in (62b) and (62c). Additionally, in order to ensure the triplet deformation measures to be solved in a consistent manner, the average strain variables F Ave and H Ave appearing in expressions (62c) and (62d) will be replaced by F a and H a . With this, and assuming the jump in traction is dominated by jump in pressure (see Remark 4), the geometric conservation Equations (62b)-(62d) reduce to Here, the strictly positive parameter S Σ J ab and the weighted average velocity field are defined as with their components being described as It is important to emphasize that strong compatibility between the different kinematic fields {F, H, J} is lost at the semi-discrete level. However, weak compatibility is maintained due to the coupled nature of the semi-discrete system of conservation equations. For visualization purposes, the current deformed geometry is recovered by integrating in time the discrete nodal velocity field obtained using (62a) With respect to the time integration of the above system (62a), (63a)-(63c) along with the geometry x (66), and keeping in mind a fast and efficient algorithm, we advocate for an explicit one-step two-stage total variation diminishing Runge-Kutta (TVD-RK) method, thoroughly reported by the authors in Reference 48 and references therein.

Entropy production
In this section, inequality (58) is assessed for the above set of semi-discrete equations (62a), (63a), (63b), (63c). For illustrative purposes, the body under consideration is said to be homogeneous. Additionally, and in line with the Godunov's theorem, 56,57 we assume piecewise constant approximation (first order in space) for variables across each control volume. Making use of expression (57), the semi-discrete counterpart of (51) is where, Equations (63a)-(63c) and (10) have been substituted in the second and third lines of (67), respectively. Subsequently, we can substitute the linear momentum conservation Equation (62a), the deformation gradient conservation Equation (63a) and, after some algebra, gives Here,Π ext denotes the semi-discrete power contribution, expressed aṡ and the positive definite matrices are Noticing that the nested summation is carried out over control volumes in (68) and the anti-symmetric nature of the first line of the right hand side, we can conclude that these terms cancel and thus (68) reduces to It is worth pointing out that the first two terms on the right hand side can be equivalently written by swapping indices a and b. Simple averaging the first line and the second line of the equation above, and noticing that S v ab = S v ba and S Indeed, the first two terms in the square bracket of (72) are always non-negative. Moreover, in the case of elasto-plasticity, the third term representing the rate of plastic dissipation is also non-negative.

ALGORITHMIC DESCRIPTION
For ease of understanding, Algorithm 2 summarizes the complete algorithmic description of the proposed finite volume methodology for large strain contact dynamics. This algorithm ¶ is implemented in modern CFD code "OpenFOAM," with an eye on large scale contact simulation in future works. Remark 7. Numerical simulations of contact problems often involve modeling the interaction of multiple bodies across a non-conforming (or non-matching) interface mesh. In order to address these scenarios, a pre-existing OpenFOAM library was employed by the proposed method. This OpenFOAM library, known as Arbitrary Mesh Interface (AMI), is based on the conservative local Galerkin projection procedure presented by Farrell and Maddison in Reference 68. By harnessing this library a projection weighting is calculated based on the overlapping face area of the two contact surfaces in the reference configuration. Since this AMI procedure is based on a surface to surface projection, additional piecewise surface-to-vertex reconstruction algorithm is then required. For illustrative purposes, Algorithm 3 summarizes the non-conforming mapping procedure in two dimensions, projecting variables from the "−" contact surface to the "+" contact surface. Its graphical representation is also depicted in Figure 4.

NUMERICAL EXAMPLES
In this section, a wide variety of numerical examples are presented in order to assess the robustness, applicability, and performance of the proposed formulation presented above. In the following sections, it is important to demonstrate the overall algorithm • ensures consistency and accuracy of the field variables at the contact interface for both conforming and non-conforming interface meshes, • guarantees long term stability by satisfying the discrete version of the second law of thermodynamics (72), and • circumvents hour-glassing and pressure instabilities even in the case of nearly incompressible material and elasto-plasticity.
In the following numerical computations, we consider only the frictionless contact where the value of friction coefficient k in (41) is set to zero. We also assume that for simplicity the contact points between potential contacting interfaces are known a priori. In terms of the temporal stability of the algorithm, the Courant-Friedrichs-Lewy number of CFL = 0.3 has been chosen. In addition, comparisons are also carried out against the results simulated using the commercial software package Abaqus. 69 From the spatial discretization standpoint, the standard linear finite element method (triangular element in two dimension and tetrahedral element in three dimension) and mean dilatation approach (quadrilateral element in two dimension and hexahedral element in three dimension) are employed, in conjunction with a set of built-in artificial viscosity parameters in order to dissipate high frequency oscillations. One-dimensional two-bar impact for similar bars The first example corresponds to the impact of two bars having equal length with an initial gap of = 0.01 m, as shown in Figure 5. Bar 1 (on the left), is travelling with a given velocity v 1 0 = 0.1 m∕s towards bar 2 (on the right) which is at rest. Material properties for both bars are exactly the same and are summarized in Table 1. The main objective of this classical benchmarked problem is to examine the robustness and reliability of the proposed algorithm in capturing contact mode transition. As reported in literature, 6,70-75 most of the methods still exhibit severe non-physical oscillations in the velocity resolution throughout the duration of contact and also post separation. Specific ad-hoc regularization procedure is generally required to limit these numerical artefacts.
In this example, a linear elastic model is considered and four different levels of mesh refinements are used. First, we demonstrate the proposed algorithm is capable of satisfying the second law of thermodynamics, and hence ensuring long-term stability. Figure 6A shows the time history of global total energy of the two bars. Its resolution after the contact at time t = 0.1 s is better represented by refining the mesh. Another interesting variable of interest is the accumulated numerical entropy (dissipation) present in the algorithm. This is achieved by integrating the Hamiltonian energy of the system described in (72) over time # , which decreases over time for the entire simulation. This is seen in Figure 6B. The total numerical dissipation is reduced when successively increasing the mesh density. In addition, Figure 6C,D illustrate the time histories of different forms of energy for bar 1 and bar 2, respectively. These include kinetic energy  dΩ V , and the total energy being defined as the summation of kinetic energy and elastic strain energy. At time t = 0, bar 1 (travelling at a given velocity v 1 0 ) has the total energy fully dominated by kinetic energy, whereas bar 2 has zero total energy (no movement and stress-free) since it is at rest. When impact takes place at time t = 0.1 s, some of the total energy of bar 1 is transferred into kinetic energy of bar 2 (both bars travel together whilst in contact) and some into elastic strain energy of bar 2 (as both bar elastically deform during the period of contact). The elastic strain energy of the bars peak at time t = 0.2 s due to the fact that the compressive stress wave arrives at the free end of the bars. The stress wave is then reflected back from the free end (now the stress wave becomes tensile stress wave) into the contact point which result in separation at time t = 0.3 s. After separation, as expected, nearly all total energy is transferred from bar 1 to bar 2, with only approximately 0.4 % of the total energy (via Mesh IV) dissipated numerically due to the use of a Riemann-based algorithm (refer to the first two terms on the right hand side of (72)).
Second, we highlight the importance of incorporating a slope limiter into the proposed algorithm in order to improve the resolution of field variables in the vicinity of shocks. Figure 7 shows the time histories of velocity and stress for bar 1.  (9) is used. Their corresponding material parameters are summarized in Table 1.
The exact (analytical) solution is also provided for verification purposes. As it can be observed, using piecewise constant representation (first order accurate in space) for field interpolation, the solutions are fairly dissipative over time for both velocity v x and axial stress xx unless excessively fine mesh is used. To enhance the accuracy, we introduce a piecewise linear reconstruction. Such enhancement, as seen in Figure 7B, gives much better resolution in stress but fails prior to separation, where non-physical oscillations are generated. In order to control these spurious modes, the classical Barth and Jespersen limiter 57 is implemented. A great improvement is observed in Figure 7B. As compared to the classical finite element method, no post-separation velocity oscillations are observed in the proposed algorithm. Finally, for qualitative comparison purposes, we monitor the velocity v x , displacement u x and stress xx evolutions at two locations, namely point A in bar 1 and point B in bar 2. Our results are in very good agreement with the given exact solutions (see Figure 8), without showing undershoots/overshoots near a discontinuity.

7.1.2
One-dimensional two-bar impact for dissimilar bars We next analyze the same impact problem but now considering a softer material (bar 1) impacts against a stiffer material (bar 2). This is achieved by reducing the value of Young's modulus E of bar 1 from 100 to 49 N∕m 2 . The remaining to position X = 10 m. A neo-Hookean constitutive model as described in (9) is used. Their corresponding material parameters are summarized in Table 1. material properties for bars 1 and 2 are exactly the same as those reported in Section 7.1.1, and are summarized again in Table 2 for completeness. The purpose of this test case is to examine the applicability of the algorithm in addressing impact between two different materials. Figure 9 illustrates the time evolutions of velocity v x , displacement u x and stress xx monitored at points A and B. The proposed algorithm can accurately capture the solutions during impact process and release process, displaying extremely good agreement with the closed-form solutions. No specific ad-hoc regularization procedure is required.

7.2
Objective 2: Spurious mechanism 7.2.1 Two dimensional compressible ring collision As previously explored in References 76 and 77, the main aim of this classical benchmark problem is to examine the capability of the proposed algorithm in alleviating unwanted spurious modes that may potentially arise in the contact (or shock) interface. The problem consists of simulating the collision of two rubber rings, with an initial gap of 8 mm, coming together at a relative speed of 1.18 m∕s. The geometry of the problem is displayed in Figure 10. In this example, a neo-Hookean model presented in (9) is considered. The values of all the relevant material parameters used can be found in Table 3. Aiming to show mesh convergence || , four successively refined meshes (see Figure 11) are used. These include (Mesh I) 2480, (Mesh II) 10,080, (Mesh III) 81,280, and (Mesh IV) 163,200 number of linear triangular elements for each ring. In order to ensure the algorithm correctly reproduces the second law of thermodynamics, the global entropy and total energy are monitored (see Figure 12A,B). Indeed, for all four meshes, both the global entropy and total energy of the system decrease over time, whereby the irreversibility is caused by numerical stabilizations introduced into the algorithm (which is precisely the square bracket term of (72)). Before contact takes place, the total energy of the system is completely dominated by kinetic energy. When time t > 6.8 ms, that is after collision takes place, the kinetic energy of the system decreases and transforms into elastic strain energy. Additionally, a very small amount of the kinetic energy also converts to monotonic decreasing numerical dissipation during a deformation process. This is seen in Figure 12C,D.
For comparison purposes, we also simulate the same problem discretized using the standard linear finite element method (using 163,840 number of linear triangular meshes with 82,944 nodes) and the mean dilatation approach (using  (9) is used. Their corresponding material parameters are summarized in Table 1.
TA B L E 2 Two-bar impact: Material parameters used in the simulation.

Bar 1 Bar 2
Young's modulus  Figure 13 illustrates the deformation process of the two rings at time t = {10, 20, 30, 40} ms, displaying how the two rings collide, bounce off and then oscillate. No spurious modes are observed. In comparison to the mean dilatation technique, very similar results in terms of deformed shape and pressure field are observed. For completeness, the time evolutions of the components of velocity v x , displacement u x and stress xx are also monitored. As shown in Figure 14, the obtained solutions converge to the results of mean dilatation technique by refining the mesh. In comparison to the standard linear finite element method, the proposed algorithm clearly outperforms it by accurately capturing stress discontinuities (in this case, xx ) without any spurious oscillations. This is seen in Figure 15C,D.

Two dimensional nearly incompressible bar impact
Similar to the objectives described in Section 7.2.1, another standard benchmark problem previously adopted in Reference 47 is considered. As shown in Figure 16, the example presents the impact of two nearly incompressible rectangular bars travelling at each other with a relative velocity of v 0 = [100, 0] T m∕s. For each bar, its width is of w = 6.4 mm and its length is of L = 32.4 mm. The normal separation between two bars is 8 mm. A neo-Hookean model is used. The values of all the simulation parameters are summarized in Table 4. For completeness, we discretize the problem using four different levels of mesh refinement, including (Mesh I) 640, (Mesh II) 2560, (Mesh III) 10,240, and (Mesh IV) 40,960 number of linear triangular elements per bar. First, a mesh refinement study for the proposed algorithm is carried out. In Figure 17, the deformation pattern of the structure predicted using a small number of elements (Mesh I) agrees well with the results obtained using finer discretizations (Mesh II-Mesh IV). As for the latter, improved pressure resolution is observed. For qualitative comparison purposes, time evolution of velocity v x and displacement u x are monitored and compared in Figure 18. Interestingly, double contact occurs between 80 and 100 μs. It is well-known that pressure checker-boarding is commonly encountered in standard linear finite elements when attempting to model materials with predominant nearly incompressible behavior. This numerical artefact can be completely resolved by the algorithm proposed in the current paper, without resorting to any ad-hoc regularization procedure. Comparing with mean dilatation technique (see Figure 19), smoother version in pressure profile is observed. Figure 20 shows the time evolution of the deformation behavior of two bars come into contact. Again, very smooth pressure profile is seen throughout the entire contact-impact process. Neither hour-glassing nor pressure checker-boarding are observed.

Objective 3: Non-matching contact interface
We now extend the above two dimensional bar impact problem to three dimension as displayed in Figure 21. This example serves the purpose to examine the accuracy and reliability of the proposed algorithm when considering non-matching meshes at the contact interface. A neo-Hookean model is chosen and the corresponding material properties remain exactly the same as the one listed in Table 4. Aiming to show mesh independent convergence for this problem, we begin this example by performing a series of non-conforming mesh refinement analysis. The (bulk) mesh information for the two bodies are presented in Table 5 and their respective non-matching interface meshes are depicted in Figure 22. As shown in Figures 23 and  24, it is remarkable that the deformation pattern together with pressure profile converge even with the use of a A neo-Hookean model as described in (9) is used. Their corresponding material parameters are summarized in Table 1.  Figure 26. The solutions indeed converge with a progressive level of mesh refinement. In order for the overall algorithm to ensure long term stability, the global entropy is monitored (see Figure 27). As expected, the global total entropy of the system decreases over time throughout the entire simulation duration.
Observe that our solution is slightly more dissipative than that of the mean dilatation technique via tri-linear hexahedral elements (see Figure 28). No spurious modes are seen comparing with the linear tetrahedral finite element method.  (9) is used. Their corresponding material parameters are summarized in Table 3.

Objective 4: Highly nonlinear impact
In the last example, we consider the rebound of a torus of outer radius R o = 40 mm, inner radius r i = 30 mm and diameter d 0 = 1 mm. The torus impacts against a rigid frictionless wall with an initial velocity of v 0 = [1.18, 0, 0] T m∕s where the separation distance between the torus and the wall is of = 4 mm. This is illustrated in Figure 29. A neo-Hookean model is first chosen with the material properties summarized in the third column of Table 6. Aiming to demonstrate the consistency of the algorithm, the domain is discretized using four different levels of refinement, namely (Mesh I) 4643, (Mesh II) 12,439, (Mesh III) 29,748, and (Mesh IV) 56,955 number of unstructured linear tetrahedral meshes. Figure 30B shows the reduction of total numerical dissipation when successively increasing the mesh density. More crucially, the global (numerical) entropy is non-positive and reduces over time for the entire simulation. In Figure 30D, the evolution in time of the kinetic energy, elastic strain energy, and total  (9) is used. Their corresponding material parameters are summarized in Table 3. A neo-Hookean model as described in (9) is used. Their corresponding material parameters are summarized in Table 3.

F I G U R E 17
Nearly incompressible bar impact: Comparison of bar impact at time t = 90 μs using various mesh refinements. In each subfigure, the first row depicts the current deformed state discretized with linear triangular mesh and the second row illustrates pressure contour. A neo-Hookean model described in (9) is used. Their corresponding material parameters are summarized in Table 4 Table 4. A neo-Hookean model described in (9) is used. Their corresponding material parameters are summarized in Table 4. energy (being the summation of both kinetic and elastic strain energy) are monitored. At the start of the simulation, the total energy of the system is completely dominated by kinetic energy as the torus is moving with the initial velocity until approximately 3.4 ms where impact occurs. In the elastic case, the kinetic energy is mostly transferred into elastic strain energy. When time is approximately 30 ms, that is when the separation begins to occur, the elastic strain energy is then converted back to kinetic energy as the torus bounces off the rigid surface. Comparing with the standard linear finite element method (101,045 number of linear tetrahedra and 20,427 number of nodes), the proposed algorithm can be used without experiencing any spurious mechanism (see Figure 31). By making use of the proposed method, a series of deformed states are shown in Figure 32, where the color contour plot indicates the pressure distribution. Moreover, the same problem is further analyzed by employing a Hencky-based von Mises plasticity (with isotropic hardening) model. Its associated materials properties are summarized in the fourth column of Table 6. When the torus impacts against a frictionless wall, the total kinetic energy is partially converted into elastic strain energy whilst most of the kinetic energy in this case is transferred into irrecoverable plastic energy dissipation. Indeed, the amount of physical plastic dissipation introduced into this model can then be monitored by integrating in time the term related to the internal dissipationḊ appearing in the Hamiltonian energy of the system described in (72). Again, as seen in Figures 33 and  34, the proposed algorithm can effectively alleviate non-physical pressure instabilities. Our solutions match well with the results obtained via mean dilatation technique (which is discretized with 29,441 number of tri-linear hexahedra and 33,793 number of nodes). For visualization purposes, Figure 35 shows the time evolution of the plastic deformation of a torus with very smooth pressure field.

F I G U R E 20
(A) (B) F I G U R E 22 Impact with non-matching interface mesh: For ease of explication, we choose to display a series of non-conforming mesh refinements on X-Y plane view. In (A), the first four rows represent the non-conforming mesh discretizations for both bars (Mesh I to Mesh IV). Their respective close-up view of the interface between two bodies can be seen in the first four columns of (B). For completeness, we also discretize the problem using conforming mesh. This is illustrated in the last row of (A) and its associated close-up view is displayed in the last column of (B).

F I G U R E 23
Impact with non-matching interface mesh: Comparison of three-dimensional bar impact at time t = 120 μs using (A) non-conforming mesh discretizations (Mesh I-Mesh IV, from first row to fourth row) and (B) conforming mesh discretization (fifth row). Color contour plot indicates the pressure profile. A neo-Hookean model (9) is used. Their corresponding material parameters are summarized in Table 4.

F I G U R E 24
Impact with non-matching interface mesh: Comparison of three-dimensional bar impact at time t = 260 μs using (A) non-conforming mesh discretizations (Mesh I-Mesh IV, from first row to fourth row) and (B) conforming mesh discretization (fifth row). Color contour plot indicates the pressure profile. A neo-Hookean model (9) is used. Their corresponding material parameters are summarized in Table 4.  Table 4.  Table 4.  (9) is used. Their corresponding material parameters are summarized in Table 4.

F I G U R E 28
Impact with non-matching interface mesh: Time evolution of global total energy using the proposed method (Mesh IV), mean dilatation and linear finite element method. Comparison of deformed shapes at time t = 205 μs is shown, where the color plot indicates pressure distribution. A neo-Hookean constitutive model as described in (9) is used. Their corresponding material parameters are summarized in Table 4. The results of mean dilatation is obtained using 5000 number of trilinear hexahedra with 6171 number of nodes for bar 1 and 8640 number of trilinear hexahedra with 10,309 number of nodes for bar 2. The results of linear finite element method is obtained using 30,778 number of linear tetrahedra with 6267 number of nodes for bar 1 and 50,261 number of linear tetrahedra with 9976 number of nodes for bar 2.

F I G U R E 29
Torus impact: Geometry and problem setup.

TA B L E 6
Torus impact: Material parameters used in the simulation.

neo-Hookean von-Mises plasticity
Young's modulus    (9) is used and the material parameters are summarized in third column of Table 6.

F I G U R E 32
Elastic torus impact: A sequence of deformed structures with pressure resolution at times t = {2.5, 5, 7.5, … , 50} ms (from left to right and top to bottom). Results obtained using the proposed algorithm discretized with linear tetrahedra (Mesh IV). A neo-Hookean model is used and the corresponding material parameters are summarized in third column of Table 6.  Table 6.

F I G U R E 35
Elasto-plastic torus impact: A sequence of deformed structures experiencing plasticity with pressure resolution at times t = {2.5, 5, 7.5, … , 50} ms (from left to right and top to bottom). Results obtained using the proposed algorithm discretized with linear tetrahedra (Mesh IV). A Hencky-based von Mises plasticity model (with isotropic hardening) is used and the corresponding material parameters are summarized in fourth column of Table 6.

CONCLUSIONS
The article presents an explicit vertex centered finite volume method for the dynamic solution of non-smooth contact problems, where a mixed system of first order conservation equations together with the associated jump conditions is used. Using the specific jump equation for the conservation of linear momentum, several dynamic contact models are derived ensuring the preservation of hyperbolic characteristic structure across contact interface. The formulation has been implemented within the modern CFD code "OpenFOAM," aiming to bridge the gap between CFD and Computational Solid Dynamics. Through the examples presented in this article, the proposed algorithm proves to perform extremely well in dynamic contact-impact problems without resorting to ad-hoc algorithmic regularization correction. Specifically, the proposed algorithm by construction overcomes a number of persistent numerical drawbacks commonly found in the literature. No spurious hour-glassing is observed and correct (smooth) pressure pattern are obtained in contrast to alternative finite element approaches, such as the well known linear tetrahedral element technology. Crucially, the overall algorithm ensures long-term stability by monitoring the global entropy production via the Hamiltonian energy of the system. The consideration of nonlinear shock wave speeds in the contact-impact conditions within the current computational framework is the next step of our work.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. . This can be achieved by simply setting the value of Poisson's ratio equal to zero. ‡ Exact energy transfer at the contact interface holds when an exact Riemann solver is used. § In the second and third terms in the parenthesis of (61), the weighted average stencil proposed by Löhner et al. 67 is used by computing the boundary flux over a boundary face (and ) in three dimensions as

ENDNOTES
where b, c are the other two nodes that together with node a define boundary face (and ). ¶ Moreover, we have also implemented the algorithm using an in-house software for the proof-of-concept one-dimensional and two-dimensional examples presented in the article. # Insofar as the linear elastic model is used in this case, we thus neglect the physical dissipation introduced by the model, that is the rate of plastic dissipation. || It is important to emphasize that the simulation of rubber ring collision could have been performed with only one rubber ring via appropriate symmetry condition. However, in our case, we decided to consider the complete two-ring impact in order to check the resolution at shock interface between two bodies.