# An explicit dynamics framework suited to highly non-smooth interface behaviors

David Dureisseix, Paul Larousse, Anthony Gravouil, Jean Di Stasio

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

## To cite this version:

David Dureisseix, Paul Larousse, Anthony Gravouil, Jean Di Stasio. An explicit dynamics framework suited to highly non-smooth interface behaviors. 2023. hal-04146784v4

HAL Id: hal-04146784

https://hal.science/hal-04146784v4

Preprint submitted on 5 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# An explicit dynamics framework suited to highly non-smooth interface behaviors

David Dureisseix, Paul Larousse, Anthony Gravouil, Jean Di Stasio

► To cite this version:

David Dureisseix, Paul Larousse, Anthony Gravouil, Jean Di Stasio. An explicit dynamics framework suited to highly non-smooth interface behaviors. Journal of Theoretical, Computational and Applied Mechanics, 2024, pp.1-21. 10.46298/jtcam.11611. hal-04146784v5

HAL Id: hal-04146784

https://hal.science/hal-04146784v5

Submitted on 24 Sep 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# An explicit dynamics framework suited to highly non-smooth interface behaviors

David DUREISSEIX<fnref n="1" />, Paul LAROUSSE<fnref n="1" />,<fnref n="2" />, Anthony GRAVOUIL<fnref n="1" />, and Jean DI STASIO<fnref n="2" />

Dynamic systems, and in particular mechanical structures, may be subjected to non-smooth loadings such as impacts or shocks. Moreover, their behavior itself may exhibit more or less non-smooth evolutions, as when fracture occurs. Therefore, robust simulation models are of interest to capture such behaviors. A particular focus is made herein on time-stepping explicit dynamics schemes to allow efficient simulations, and non-smoothness is embedded within the discrete resolution model, so that robust simulations can be obtained, with a minimum number of numerical parameters. The original contributions of this article lie in the way the non-smooth behavior is formulated to be embedded in an explicit dynamics framework. This study focuses on the solver for dynamics with non-smooth interface behavior, rather than on the behavior models themselves. The applications concern non-smooth interface behaviors at macroscopic scale, between displacement *spatial* jump on the 2D interface surface with no thickness, and interfacial force distributions acting on the bodies apart the interface. The proposed test cases concern in a first step contact and perfectly plastic interface behavior (for illustrative purpose, on a 0D example). The last numerical test deals with contact, friction, fracture and adhesion for an extrinsic perfectly brittle interface behavior, to exemplify the feasibility on a full 3D finite element model.

Keywords matrix-free, perfect plasticity, brittle fracture, finite elements, time non-locality

## 1 Non-smoothness is difficult

Non-smooth behaviors in mechanical models often lead to difficulties in their numerical simulation, where the problem is often qualified as stiff. For instance, contact or impact problems, crack propagation, brittle behaviors... involve such difficulties (Brogliato et al. 2002; Acary et al. 2018). A quasi-static model in this case is usually difficult to tackle, while the use of dynamical behavior may induce some regularization, although rigid models persist in exhibiting multiple solutions (Moreau 2004; Moreau 2006; Alart 2014). In such cases, but for deformable models, explicit dynamics schemes are useful to reduce the computational cost. However, additional regularizations such as penalization are then not suited, since explicit schemes usually exhibit a critical time step which is driven to a possibly very small value with high penalization factors in the model.

In this study the focus is therefore on the use of a symplectic (when applied to Hamiltonian systems) and explicit scheme, namely the CD-Lagrange scheme (Fekak et al. 2017), embedding by construction the treatment of non-smooth evolutions, without the need for any penalization with artificial numerical parameters. Even for non-Hamiltonian cases (as those dealt with in this article), this family of integration schemes produced good energy preserving and balance equations properties as well.

Although non-smooth behaviors may also occur in the bulk of the structure, many different cases often arise also in the contact of an interface (assembling between mechanical parts, heterogeneous structure made by gluing of different parts, interaction between mechanical systems...) Therefore, this article focuses on these kinds of non-smooth behaviors. These stiff problems usually require particular formulations to overcome their resolution issue (Jean et al. 2001; de Saxcé 2022).

The CD-Lagrange scheme embedding the non-smooth Moreau framework (Moreau 1986; Moreau 2011; Dubois et al. 2018) involves an interface behavior expressed not with a classical stress vs displacement *spatial* jump function, but with an impulse-velocity *spatial* jump relationship. It

<footnote n="1">¹ INSA Lyon, CNRS, LaMCoS, UMR5259, F-69621 Villeurbanne, France</footnote>

<footnote n="2">² Centre de technologie de Ladoux, Manufacture française de pneumatiques Michelin, F-63000 Clermont-Ferrand, France</footnote>
**Identifiers**

DOI 10.46298/jtcam.11611

HAL hal-04146784v5

**History**

Received July 19, 2023

Accepted May 21, 2024

Published Sept 18, 2024

**Associate Editor**

Olivier THOMAS

**Reviewers**

Anonymous

Serge DUMONT

**Open Review**

HAL hal-04693281

**Supplementary Material**

See addendum

**Licence**

CC BY 4.0

©The Authors

# An explicit dynamics framework suited to highly non-smooth interface behaviors

David DUREISSEIX<fnref n="1" />, Paul LAROUSSE<fnref n="1" />,<fnref n="2" />, Anthony GRAVOUIL<fnref n="1" />, and Jean DI STASIO<fnref n="2" />

Dynamic systems, and in particular mechanical structures, may be subjected to non-smooth loadings such as impacts or shocks. Moreover, their behavior itself may exhibit more or less non-smooth evolutions, as when fracture occurs. Therefore, robust simulation models are of interest to capture such behaviors. A particular focus is made herein on time-stepping explicit dynamics schemes to allow efficient simulations, and non-smoothness is embedded within the discrete resolution model, so that robust simulations can be obtained, with a minimum number of numerical parameters. The original contributions of this article lie in the way the non-smooth behavior is formulated to be embedded in an explicit dynamics framework. This study focuses on the solver for dynamics with non-smooth interface behavior, rather than on the behavior models themselves. The applications concern non-smooth interface behaviors at macroscopic scale, between displacement spatial jump on the 2D interface surface with no thickness, and interfacial force distributions acting on the bodies apart the interface. The proposed test cases concern in a first step contact and perfectly plastic interface behavior (for illustrative purpose, on a 0D example). The last numerical test deals with contact, friction, fracture and adhesion for an extrinsic perfectly brittle interface behavior, to exemplify the feasibility on a full 3D finite element model.

**Keywords:** matrix-free, perfect plasticity, brittle fracture, finite elements, time non-locality

## 1 Non-smoothness is difficult

Non-smooth behaviors in mechanical models often lead to difficulties in their numerical simulation, where the problem is often qualified as stiff. For instance, contact or impact problems, crack propagation, brittle behaviors involve such difficulties (Brogliato et al. 2002; Acary et al. 2018). A quasi-static model in this case is usually difficult to tackle, while the use of dynamical behavior may induce some regularization, although rigid models persist in exhibiting multiple solutions (Moreau 2004; Moreau 2006; Alart 2014). In such cases, but for deformable models, explicit dynamics schemes are useful to reduce the computational cost. However, additional regularizations such as penalization are then not suited, since explicit schemes usually exhibit a critical time step which is driven to a possibly very small value with high penalization factors in the model.

In this study the focus is therefore on the use of a symplectic (when applied to Hamiltonian systems) and explicit scheme, namely the CD-Lagrange scheme (Fekak et al. 2017), embedding by construction the treatment of non-smooth evolutions, without the need for any penalization with artificial numerical parameters. Even for non-Hamiltonian cases (as those dealt with in this article), this family of integration schemes produced good energy preserving and balance equations properties as well.

Although non-smooth behaviors may also occur in the bulk of the structure, many different cases often arise also in the contact of an interface (assembling between mechanical parts, heterogeneous structure made by gluing of different parts, interaction between mechanical systems...) Therefore, this article focuses on these kinds of non-smooth behaviors. These stiff problems usually require particular formulations to overcome their resolution issue (Jean et al. 2001; De Saxcé 2022).

<footnote n="1">¹ INSA Lyon, CNRS, LaMCoS, UMR5259, Villeurbanne, France</footnote>

<footnote n="2">² Centre de technologie de Ladoux, Manufacture française de pneumatiques Michelin, Clermont-Ferrand, France</footnote>
has been tested and used for different cases of non-regularized frictional contact (Fekak et al. 2017; Di Stasio et al. 2019).

We restrict the study here to macroscopic phenomenological models of interfaces. There are many references dealing with lower scale interface phenomena (chemistry, bonds...) which are not under consideration here. Apart from the contact behaviors, there are still other non-smooth behaviors, for instance those arising from piecewise linear relationships, such as plasticity with a yield threshold and without viscosity. We then test the explicit scheme first on interfaces that exhibit such plastic behaviors. Indeed, the perfectly plastic behavior in the bulk of the material may be difficult to handle in standard finite-element codes. Indeed, the plastic threshold may lead algorithms to fail to converge, especially when non-admissible intermediate solutions are seek for during iterations (Meyer et al. 2022). Here, we focus on the plasticity developed in the interface model, for two phenomenological versions corresponding to two physical situations, and to derive challenging cases, we push the non-smoothness using perfect plasticity models. In a second step, a brittle interface with adhesion recovery is studied. The brittle behavior corresponds to a fracture of the interface which is a brutal decohesion. Although fracture models are usually derived from energy release when creating fracture surface (Francfort 2006; Kiener et al. 2022), the brittle fracture is usually associated with a low energy. The contact part of this interface behavior (in compression) is one source of non-smoothness (Jean et al. 2001). The infinite stiffness in traction of the interface (before fracture) is a second one (Collins-Craft et al. 2022) as it was also for the plastic interface. Finally the limit case of perfect brittle fracture, which is tested herein, exhibits a higher non-smooth response, corresponding to a zero energy dissipation. Although the physical interpretation of this limit case is not trivial (usually, looking at a smaller physical scale leads to some elastic behavior and some dissipation, (Raffa et al. 2022; Lebon et al. 2022)), this case is used herein to challenge the models, leading to difficult numerical non-smoothness issues.

It is noteworthy that these different cases do have some similarities with the non-smoothness involved in the frictional contact models: the rigid, perfect plastic behavior leads to a normal behavior looking similar to the tangential friction of Tresca model; then, the brittle fracture leads to a generalization of the perfect normal contact behavior, the so-called Hertz-Signorini-Moreau (HSM) or Karush-Kuhn-Tucker (KKT) conditions (Wriggers et al. 2007). These cases will be described, implemented and tested in this article. The target cases considered herein deal with a constant time step $\Delta t$, and the small displacement assumption. Moreover, compatible spatial discretizations on each side of an interface are assumed for sake of simplicity; this issue can be overcome for instance by considering the mortar approach (Ben Belgacem et al. 1998; Casadei et al. 2002; Wohlmuth 2000; Popp et al. 2014; Pinto Carvalho et al. 2022).

The main developments consist in keeping the main steps of the CD-Lagrange framework that is not a novelty for this paper (configuration update, free velocity, computation of the global velocity, and update of the displacement), but a central part is the innovative algorithmic issue of this article: the explicit local non-smoothness resolution (depending on the nature of the non-smoothness). This algorithmic part is therefore detailed, with the constraint to keep the resolution explicit and matrix-free.

Clearly, the different variants of the original Moreau schemes (Moreau 1988) (with the explicit gap prediction) are probably also suited for addressing the proposed non-smoothness issues. The classical version of this scheme has already been compared for benchmark cases in (Di Stasio et al. 2019) where symplecticity has been discussed. The comparison with the current proposed test cases are left as a direct perspective to this work.

## 2 A rigid, perfectly plastic interface behavior

We focus herein on perfectly plastic behavior of interfaces, as a macroscopic phenomenological model. To enforce non-smoothness, no elastic part is considered in the model, which is related to the so-called extrinsic case in fracture models (Park et al. 2011; Zhang et al. 2007; Collins-Craft et al. 2022); this is the limit case of an elastic interface when its stiffness is driven to infinity. To help resolution in standard finite element codes, these non-smoothness are often managed with regularization or penalization, i.e. using some numerical parameters similar to artificial stiffness.
The CD-Lagrange scheme embedding the non-smooth Moreau framework (Moreau 1986; Moreau 2011; Dubois et al. 2018) involves an interface behavior expressed not with a classical stress vs displacement spatial jump function, but with an impulse-velocity spatial jump relationship. It has been tested and used for different cases of non-regularized frictional contact (Fekak et al. 2017; Di Stasio et al. 2019).

We restrict the study here to macroscopic phenomenological models of interfaces. There are many references dealing with lower scale interface phenomena (chemistry, bonds...) which are not under consideration here. Apart from the contact behaviors, there are still other non-smooth behaviors, for instance those arising from piecewise linear relationships, such as plasticity with a yield threshold and without viscosity. We then test the explicit scheme first on interfaces that exhibit such plastic behaviors. Indeed, the perfectly plastic behavior in the bulk of the material may be difficult to handle in standard finite-element codes. Indeed, the plastic threshold may lead algorithms to fail to converge, especially when non-admissible intermediate solutions are seek for during iterations (Meyer et al. 2022). Here, we focus on the plasticity developed in the interface model, for two phenomenological versions corresponding to two physical situations, and to derive challenging cases, we push the non-smoothness using perfect plasticity models. In a second step, a brittle interface with adhesion recovery is studied. The brittle behavior corresponds to a fracture of the interface which is a brutal decohesion. Although fracture models are usually derived from energy release when creating fracture surface (Francfort 2006; Kiener et al. 2022), the brittle fracture is usually associated with a low energy. The contact part of this interface behavior (in compression) is one source of non-smoothness (Jean et al. 2001). The infinite stiffness in traction of the interface (before fracture) is a second one (Collins-Craft et al. 2022) as it was also for the plastic interface. Finally the limit case of perfect brittle fracture, which is tested herein, exhibits a higher non-smooth response, corresponding to a zero energy dissipation. Although the physical interpretation of this limit case is not trivial (usually, looking at a smaller physical scale leads to some elastic behavior and some dissipation (Raffa et al. 2022; Lebon et al. 2022)), this case is used herein to challenge the models, leading to difficult numerical non-smoothness issues.

It is noteworthy that these different cases do have some similarities with the non-smoothness involved in the frictional contact models: the rigid, perfect plastic behavior leads to a normal behavior looking similar to the tangential friction of Tresca model; then, the brittle fracture leads to a generalization of the perfect normal contact behavior, the so-called Hertz-Signorini-Moreau (HSM) or Karush-Kuhn-Tucker (KKT) conditions (Wriggers et al. 2007). These cases will be described, implemented and tested in this article. The target cases considered herein deal with a constant time step $\Delta t$, and the small displacement assumption. Moreover, compatible spatial discretizations on each side of an interface are assumed for sake of simplicity; this issue can be overcome for instance by considering the mortar approach (Belgacem et al. 1998; Casadei et al. 2002; Wohlmuth 2000; Popp et al. 2014; Carvalho et al. 2022).

The main developments consist in keeping the main steps of the CD-Lagrange framework that is not a novelty for this paper (configuration update, free velocity, computation of the global velocity, and update of the displacement), but a central part is the innovative algorithmic issue of this article: the explicit local non-smoothness resolution (depending on the nature of the non-smoothness). This algorithmic part is therefore detailed, with the constraint to keep the resolution explicit and matrix-free.

Clearly, the different variants of the original Moreau schemes (Moreau 1988) (with the explicit gap prediction) are probably also suited for addressing the proposed non-smoothness issues. The classical version of this scheme has already been compared for benchmark cases in (Di Stasio et al. 2019) where symplecticity has been discussed. The comparison with the current proposed test cases are left as a direct perspective to this work.

## 2 A rigid, perfectly plastic interface behavior

We focus herein on perfectly plastic behavior of interfaces, as a macroscopic phenomenological model. To enforce non-smoothness, no elastic part is considered in the model, which is related to the so-called extrinsic case in fracture models (Park et al. 2011; Zhang et al. 2007; Collins-Craft
But the higher the artificial stiffness is, the stiffer the problem is. When one is concerned with explicit dynamics schemes, this cannot be a solution due to the classical conditional stability issue, reducing the critical time step, and rendering the problem ill-conditioned. A dedicated integration scheme to deal with non-smoothness may therefore be a suited solution.

In this section, we first recall the principles of the explicit time-stepping CD-Lagrange scheme, with perfect frictionless unilateral impacts. Then we move to two new test cases, involving perfect plasticity in the interface model: no elasticity is introduced to emphasize the non-smoothness feasibility of the approach.

## 2.1 Unilateral contact CD-Lagrange scheme

In a first step, the classical contact case is recalled to emphasize the different design stages of the algorithm.

For linear dynamics, a spatially discretized (by finite elements) system of mass and stiffness matrices $M$ and $K$, submitted to an external force is driven by the second-order in time ordinary differential equation (ODE) of the form $M\ddot{U} + KU = f_{\text{ext}}$, where $U(t)$ is the displacement degree-of-freedom vector and $f_{\text{ext}}(t)$ the external generalized nodal force vector. The CD-Lagrange scheme is based on the central difference time integration scheme, using a desynchronized time discretization where the indices $n$ and $n + 1/2$ are related to the time steps $t_n$ and $t_{n+1/2}$. This central difference variant (3-step version) is the one proposed in (Belytschko et al. 2014; Hairer et al. 2003), algorithmically equivalent to other versions, except when using non-smooth impact behaviors and to preserve symplecticity property for Hamiltonian systems of this kind. Moreover, as a time-variational scheme, the order of the ODE is reduced by an order 1, and the time-discretized dynamic problem is

$$M(V_{n+3/2} - V_{n+1/2}) + hKU_{n+1} = h f_{\text{ext},n+1} \quad (1)$$

where $h$ is the time step ($h = \Delta t$ assumed herein uniform), and $V$ is the velocity vector, such that $U_{n+1} = U_n + hV_{n+1/2}$ is the explicit computation of the incremental displacement, provided that both $U_n$ and $V_{n+1/2}$ are known, and that $U_{n+1}$ and $V_{n+3/2}$ are seek for. Contrary to event-driven schemes, the present time-stepping approach does not select the time steps to match with non-smooth event occurring times, allowing the use of pre-selected uniform time steps.

Moreover, to get an explicit and matrix-free resolution, a classical solution is to use the lumped mass matrix, so $M$ is diagonal, as considered herein. If part of the degrees of freedom (dof) are subjected to unilateral contact condition (without friction in a first step), one may use the localization matrix $L$ containing the dofs selection and projection onto the normal to the contact (which is constant for the small perturbation case studied herein). Therefore, the normal velocity spatial jump at contact interface is $v = LV$, and the normal gap is $g = g_0 + LU$ where $g_0$ is the initial gap. For impact problems, and moreover for spatially discretized ones, the velocity may exhibit jumps in time, and the acceleration, as well as the contact force are not always defined. The reduced-order system is therefore suited for such a situation, once one introduces the so-called impact impulse $r$ that replaces the contact force and is a convergent quantity when the time step is refined. The discretized non-smooth dynamics evolution equation is then written as

$$M(V_{n+3/2} - V_{n+1/2}) + hKU_{n+1} = h f_{\text{ext},n+1} + L^T r_{n+3/2} \quad (2)$$

or

$$V_{n+3/2} = V_{\text{free}} + M^{-1}L^T r_{n+3/2} \quad (3)$$

where $V_{\text{free}} = V_{n+1/2} + M^{-1}h(f_{\text{ext},n+1} - KU_{n+1})$ would be the velocity if no contact occurs, i.e. if $r_{n+3/2} = 0$.

The reduced dynamics problem consists in expressing from (3) the relation between local unknown quantities at the contact interface $v_{n+3/2}$ and $r_{n+3/2}$. This is done by pre-multiplying with $L$ to get

$$v_{n+3/2} = v_{\text{free}} + Hr_{n+3/2} \quad (4)$$
et al. 2022); this is the limit case of an elastic interface when its stiffness is driven to infinity. To help resolution in standard finite element codes, these non-smoothness are often managed with regularization or penalization, i.e. using some numerical parameters similar to artificial stiffness. But the higher the artificial stiffness is, the stiffer the problem is. When one is concerned with explicit dynamics schemes, this cannot be a solution due to the classical conditional stability issue, reducing the critical time step, and rendering the problem ill-conditioned. A dedicated integration scheme to deal with non-smoothness may therefore be a suited solution.

In this section, we first recall the principles of the explicit time-stepping CD-Lagrange scheme, with perfect frictionless unilateral impacts. Then we move to two new test cases, involving perfect plasticity in the interface model: no elasticity is introduced to emphasize the non-smoothness feasibility of the approach.

## 2.1 Unilateral contact CD-Lagrange scheme

In a first step, the classical contact case is recalled to emphasize the different design stages of the algorithm.

For linear dynamics, a spatially discretized (by finite elements) system of mass and stiffness matrices $M$ and $K$, submitted to an external force is driven by the second-order in time ordinary differential equation (ODE) of the form $M\ddot{U} + KU = f_{\text{ext}}$, where $U(t)$ is the displacement degrees-of-freedom vector and $f_{\text{ext}}(t)$ the external generalized nodal force vector. The CD-Lagrange scheme is based on the central difference time integration scheme, using a desynchronized time discretization where the indices $n$ and $n + 1/2$ are related to the time steps $t_n$ and $t_{n+1/2}$. This central difference variant (3-step version) is the one proposed in (Belytschko et al. 2014; Hairer et al. 2003), algorithmically equivalent to other versions, except when using non-smooth impact behaviors and to preserve symplecticity property for Hamiltonian systems of this kind. Moreover, as a time-variational scheme, the order of the ODE is reduced by an order 1, and the time-discretized dynamic problem is

$$M(V_{n+3/2} - V_{n+1/2}) + hKU_{n+1} = h f_{\text{ext},n+1} \quad (1)$$

where $h$ is the time step ($h = \Delta t$ assumed herein uniform), and $V$ is the velocity vector, such that $U_{n+1} = U_n + hV_{n+1/2}$ is the explicit computation of the incremental displacement, provided that both $U_n$ and $V_{n+1/2}$ are known, and that $U_{n+1}$ and $V_{n+3/2}$ are seek for. Contrary to event-driven schemes, the present time-stepping approach does not select the time steps to match with non-smooth event occurring times, allowing the use of pre-selected uniform time steps.

Moreover, to get an explicit and matrix-free resolution, a classical solution is to use the lumped mass matrix, so $M$ is diagonal, as considered herein. If part of the degrees of freedom (dof) are subjected to unilateral contact condition (without friction in a first step), one may use the localization matrix $L$ containing the dofs selection and projection onto the normal to the contact (which is constant for the small perturbation case studied herein). Therefore, the normal velocity spatial jump at contact interface is $v = LV$, and the normal gap is $g = g_0 + LU$ where $g_0$ is the initial gap. For impact problems, and moreover for spatially discretized ones, the velocity may exhibit jumps in time, and the acceleration, as well as the contact force are not always defined. The reduced-order system is therefore suited for such a situation, once one introduces the so-called impact impulse $r$ that replaces the contact force and is a convergent quantity when the time step is refined. The discretized non-smooth dynamics evolution equation is then written as

$$M(V_{n+3/2} - V_{n+1/2}) + hKU_{n+1} = h f_{\text{ext},n+1} + L^T r_{n+3/2} \quad (2)$$

or

$$V_{n+3/2} = V_{\text{free}} + M^{-1}L^T r_{n+3/2} \quad (3)$$

where $V_{\text{free}} = V_{n+1/2} + M^{-1}h(f_{\text{ext},n+1} - KU_{n+1})$ would be the velocity if no contact occurs, i.e. if $r_{n+3/2} = 0$.
where $v_{free} = LV_{free}$ and $H = LM^{-1}L^T$ is the so-called Delassus operator. When considering compatible finite element meshes between impacting bodies, or impact with a fixed rigid part, the operator $L$ is local to each impacting nodes. Due to the properties of the lumped mass matrix $M$; this leads to a Delassus operator which is diagonal, definite positive and spherical per node (spatial directions are equivalent, as in the mass matrix $M$ for volumic finite elements) as considered herein.

To close the problem, the local contact behavior should be added. When expressed with a relationship between gap $g$ and impulse $r$, it leads to the KKT condition: $0 \le g \perp r \ge 0$. This fits within the thermodynamics framework, using a free energy potential $\Psi(g) = I_{\mathbb{R}^+(g)}$ where $I_{\mathbb{R}^+}$ is the indicator function of the convex positive cone $\mathbb{R}^+$ (with a null value for a positive argument, and $+\infty$ otherwise). Clearly, the non-smoothness arises from the non differential property of this potential, so that the state law involves a subdifferential operation as

$$ -r \in \partial_g \Psi(g) \qquad (5) $$

which is equivalent to the previous KKT condition.

This condition can be replaced, thanks to Moreau viability lemma (Moreau 1999; Dubois et al. 2018), with a relation between impulse and velocity spatial jump, that can also be interpreted as a constitutive relation of the form: $\mathcal{R}(v, r; g) = 0$. This relation is also neither smooth nor univoque. It reads:

$$ \text{If } g > 0 \text{ then: } r = 0 \qquad (6) $$

$$ \text{Else: } 0 \le v \perp r \ge 0 \qquad (7) $$

Finally the problem is (4) together with the behavior $\mathcal{R}(v_{n+3/2}, r_{n+3/2}; g_{n+1}) = 0$ and with diagonal Delassus operator, it provides an explicit solution as well. Once $r_{n+3/2}$ is known, the dynamics (3) allows to get the full velocity vector $V_{n+3/2}$.

The classical Coulomb friction laws can also be taken easily into account (Fekak et al. 2017; Di Stasio et al. 2019), with an underlying thermodynamics framework suited to this non-associated model thanks to the bi-potential method (de Saxcé et al. 1998; Feng et al. 2005). Indeed, there are other explicit formulations for contact problems, as the Nitsche approach in (Chouly et al. 2018) but dedicated to deformable bodies.

## 2.2 A simple oD test case for normal rigid, perfectly plastic behavior

In this new test case, which is expected to be at least as non-smooth as the previous contact case, the only difference resides in the interface behavior, that should be expressed also as an impulse-velocity spatial jump relationship.

The classical phenomenological representation of the interface behavior, in the normal direction to the interface, is selected to be the parallel association of a perfect contact behavior (with index $c$) and a purely plastic behavior (with index $p$) with the stress threshold $\sigma_Y$, as depicted on Figure 1. This association leads to the following assembly relationships: $r = r_c + r_p$ for the impulses and $v_c = v_p = v$ for the velocity spatial jump (the gap is also the same for both).

<figure><img src="image_3.png" /><figcaption>Figure 1: Sketch of the phenomenological interface behavior</figcaption></figure>

Such a model may correspond to a physical situation where a thin layer of cohesive plastic material separates two volumic elastic bodies, see for instance some dedicated plasticity models (de los Ríos et al. 2009), elasto-plasticity (Spada et al. 2009; Brisotto et al. 2018), visco-plasticity (Monchiet et al. 2010), etc. for which the non-smoothness is nevertheless weaker than the one of perfectly plastic model.

Indeed, at a macroscopic scale (larger that the interphase thickness), the interphase is represented by a null-thickness interface with a dedicated plastic behavior. The contact part of the behavior is only a crude model of the fact that the possible compression of the interphase is limited.
The reduced dynamics problem consists in expressing from Equation (3) the relation between local unknown quantities at the contact interface $v_{n+3/2}$ and $r_{n+3/2}$. This is done by pre-multiplying with $L$ to get

$$v_{n+3/2} = v_{\text{free}} + Hr_{n+3/2} \quad (4)$$

where $v_{\text{free}} = LV_{\text{free}}$ and $H = LM^{-1}L^{\top}$ is the so-called Delassus operator. When considering compatible finite element meshes between impacting bodies, or impact with a fixed rigid part, the operator $L$ is local to each impacting nodes. Due to the properties of the lumped mass matrix $M$; this leads to a Delassus operator which is diagonal, definite positive and spherical per node (spatial directions are equivalent, as in the mass matrix $M$ for solid finite elements) as considered herein.

To close the problem, the local contact behavior should be added. When expressed with a relationship between gap $g$ and impulse $r$, it leads to the KKT condition $0 \le g \perp r \ge 0$. This fits within the thermodynamics framework, using a free energy potential $\Psi(g) = I_{\mathbb{R}^+(g)}$ where $I_{\mathbb{R}^+}$ is the indicator function of the convex positive cone $\mathbb{R}^+$ (with a null value for a positive argument, and $+\infty$ otherwise). Clearly, the non-smoothness arises from the non differential property of this potential, so that the state law involves a subdifferential operation as

$$-r \in \partial_g \Psi(g) \quad (5)$$

which is equivalent to the previous KKT condition.

This condition can be replaced, thanks to Moreau viability lemma (Moreau 1999; Dubois et al. 2018), with a relation between impulse and velocity spatial jump, that can also be interpreted as a constitutive relation of the form $\mathcal{R}(v, r; g) = 0$. This relation is also neither smooth nor single-valued. It reads

$$\text{If } g > 0 \text{ then } r = 0, \quad (6)$$

$$\text{Else } 0 \le v \perp r \ge 0. \quad (7)$$

Finally the problem is composed of Equation (4) together with the behavior $\mathcal{R}(v_{n+3/2}, r_{n+3/2}; g_{n+1}) = 0$ and with diagonal Delassus operator, it provides an explicit solution as well. Once $r_{n+3/2}$ is known, the dynamics in Equation (3) allows to get the full velocity vector $V_{n+3/2}$.

The classical Coulomb friction laws can also be taken easily into account (Fekak et al. 2017; Di Stasio et al. 2019), with an underlying thermodynamics framework suited to this non-associated model thanks to the bi-potential method (De Saxcé et al. 1998; Feng et al. 2005). Indeed, there are other explicit formulations for contact problems, as the Nitsche approach in (Chouly et al. 2018) but dedicated to deformable bodies.

## 2.2 A simple oD test case for normal rigid, perfectly plastic behavior

In this new test case, which is expected to be at least as non-smooth as the previous contact case, the only difference resides in the interface behavior, that should be expressed also as an impulse-velocity spatial jump relationship.

The classical phenomenological representation of the interface behavior, in the normal direction to the interface, is selected to be the parallel association of a perfect contact behavior (with index $c$) and a purely plastic behavior (with index $p$) with the stress threshold $\sigma_y$, as depicted in Figure 1. This association leads to the assembly relationships $r = r_c + r_p$, for the

<figure><img src="image_2.png" /><figcaption>**Figure 1** Sketch of the phenomenological interface behavior.</figcaption></figure>

impulses and $v_c = v_p = v$, for the velocity spatial jump (the gap is also the same for both).

Such a model may correspond to a physical situation where a thin layer of cohesive plastic material separates two solid elastic bodies, see for instance some dedicated plasticity models (Ríos
For the present model, the contact part of the behavior is similar to the one of the previous section. The plastic part of the behavior is usually modeled within the thermodynamic framework using a state potential for reversible part of the behavior, and a pseudo-potential of dissipation for the irreversible part (Lemaitre et al. 1994). Since the present model is perfectly plastic, no elastic part is involved and only a dissipation potential is used. The state variables are the normal stress $\sigma_p$ and the displacement spatial jump $u_p$ (the equivalent for an interface behavior to the strain for a volumic plasticity model). The non-smooth dissipation potential, for the oD case, is defined as

$$ \Phi^*(\sigma_p) = I_{\mathbb{R}^-}(f(\sigma_p)) \quad \text{with} \quad f(\sigma_p) = |\sigma_p| - \sigma_Y $$

where $f(\sigma_p)$ is the yield function, $\sigma_Y$ is the yield stress and $I_{\mathbb{R}^-}$ is the indicator function of the negative cone $\mathbb{R}^-$ (with a null value for a negative argument, and $+\infty$ otherwise). Clearly, the non-smoothness again arises from the non differential property of this potential (Moreau 2011), so that the evolution law involves also a subdifferential operation as

$$ \dot{u}_p \in \partial_{\sigma_p} \Phi^* $$

leading to

$$ \text{If } f(\sigma_p) < 0 \text{ then: } \dot{u}_p = 0 \qquad (8) $$

$$ \text{Else: } \dot{u}_p = \lambda \operatorname{sgn}(\sigma_p) \text{ with } \lambda \ge 0 \qquad (9) $$

Expressed with velocity spatial jump $v_p = \dot{u}_p$ and impulse $r_p = -h\sigma_p$ (the impulse corresponding to the normal plastic stress, this last one being positive in traction along the normal direction), this interface behavior is depicted on Figure 2, denoting $r_Y = h\sigma_Y$. Some corresponding stress vs velocity spatial jump answers are illustrated on Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 2: The considered constitutive behavior, as an impulse/velocity spatial jump relationship</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: The considered constitutive behavior, as a stress/displacement spatial jump relationship</figcaption></figure>

When the gap is strictly positive, the contact is disabled, so $r_c = 0$, $r = r_p$, and the reduced dynamics reads

$$ v = v_{\text{free}} + Hr_p \qquad (10) $$

With the diagonal and positive properties of $H$, the solution to the intersection of the linear reduced dynamics and non-smooth behavior relations is also explicit as illustrated on Figure 4 and detailed in Algorithm 1.

When the gap is not strictly positive, both the contact and the plastic behaviors are activated. This situation corresponds to the case of an hyperstatic association of two relationships, therefore leading not to a unique solution for each. Nevertheless, the solution is still unique for the sum of the impulses $r = r_c + r_p$ with the reduced dynamics

$$ v = v_{\text{free}} + Hr \qquad (11) $$

Indeed, consider the cases separation:
et al. 2009), elasto-plasticity (Spada et al. 2009; Brisotto et al. 2018), visco-plasticity (Monchiet et al. 2010), etc. for which the non-smoothness is nevertheless weaker than the one of perfectly plastic model.

Indeed, at a macroscopic scale (larger that the interphase thickness), the interphase is represented by a null-thickness interface with a dedicated plastic behavior. The contact part of the behavior is only a crude model of the fact that the possible compression of the interphase is limited.

For the present model, the contact part of the behavior is similar to the one of the previous section. The plastic part of the behavior is usually modeled within the thermodynamic framework using a state potential for reversible part of the behavior, and a pseudo-potential of dissipation for the irreversible part (Lemaitre et al. 1994). Since the present model is perfectly plastic, no elastic part is involved and only a dissipation potential is used. The state variables are the normal stress $\sigma_p$ and the displacement spatial jump $u_p$ (the equivalent for an interface behavior to the strain for a volumic plasticity model). The non-smooth dissipation potential, for the oD case, is defined as

$$ \Phi^*(\sigma_p) = I_{\mathbb{R}^-}(f(\sigma_p)) \quad \text{with} \quad f(\sigma_p) = |\sigma_p| - \sigma_Y \qquad (8) $$

where $f(\sigma_p)$ is the yield function, $\sigma_Y$ is the yield stress and $I_{\mathbb{R}^-}$ is the indicator function of the negative cone $\mathbb{R}^-$ (with a null value for a negative argument, and $+\infty$ otherwise). Clearly, the non-smoothness again arises from the non differential property of this potential (Moreau 2011), so that the evolution law involves also a subdifferential operation as $\dot{u}_p \in \partial_{\sigma_p} \Phi^*$ leading to

$$ \text{If } f(\sigma_p) < 0 \text{ then } \dot{u}_p = 0, \qquad (9) $$

$$ \text{Else } \dot{u}_p = \lambda \operatorname{sgn}(\sigma_p) \text{ with } \lambda \ge 0. \qquad (10) $$

Expressed with velocity spatial jump $v_p = \dot{u}_p$ and impulse $r_p = -h\sigma_p$ (the impulse corresponding to the normal plastic stress, this last one being positive in traction along the normal direction), this interface behavior is depicted in Figure 2, denoting $r_Y = h\sigma_Y$. Some corresponding stress vs velocity spatial jump answers are illustrated in Figure 3.

<figure><img src="image_3.png" /><figcaption>**Figure 3** Constitutive behavior as a stress/displacement spatial jump relationship.</figcaption></figure>

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

When the gap is strictly positive, the contact is disabled, so $r_c = 0$, $r = r_p$, and the reduced dynamics reads $v = v_{\text{free}} + Hr_p$. With the diagonal and positive properties of $H$, the solution to the intersection of the linear reduced dynamics and non-smooth behavior relations is also explicit as illustrated in Figure 4 and detailed in Algorithm 1. When the gap is not strictly positive, both the contact and the plastic behaviors are activated. This situation corresponds to the case of an hyperstatic association of two relationships, therefore leading not to a unique solution for each. Nevertheless, the solution is still unique for the sum of the impulses $r = r_c + r_p$ with the reduced dynamics

$$ v = v_{\text{free}} + Hr. \qquad (11) $$
<figure><img src="image_6.png" /><figcaption>**Figure 4:** The solution as the intersection of the behavior (plain line) and the reduced dynamics (dotted line)</figcaption></figure>

* if $v > 0$, then $r_c = 0$, so $r = r_p = -r_Y$ is computed as before. One therefore gets $v = v_{free} - Hr_Y := \tilde{v}$ that should be positive in this case;

* in the counter case, $v = 0$, and the reduced dynamics reads $0 = v_{free} + Hr = (v_{free} - Hr_Y) + H(r_c + r_p + r_Y)$, so $0 = \tilde{v} + H(r_c + r_Y - h\sigma_p)$ for which we should have $r_c \ge 0$ and $0 \le r_Y - h\sigma_p \le 2r_Y$. Now if $\tilde{v} > 0$ this is inconsistent and corresponds to the previous case. Therefore $\tilde{v} \le 0$ and $r = r_c - h\sigma_p = -H^{-1}v_{free} := \tilde{r}$

This cases separation can therefore be summarized as

* if $\tilde{v} > 0$ (which is equivalent to $\tilde{r} < -r_Y$), then $r = -r_Y$

* else $r = \tilde{r}$

This can be simplified also, as detailed in Algorithm 1.

**Input:** displacement $U_n$, velocity $V_{n+1/2}$

**Output:** displacement $U_{n+1}$, velocity $V_{n+3/2}$

1: Displacement $U_{n+1} \leftarrow U_n + hV_{n+1/2}$

► Explicit configuration

2: $g_{n+1} \leftarrow g_0 + LU_{n+1}$

3: $V_{free} \leftarrow V_{n+1/2} + M^{-1}h(f_{ext,n+1} - KU_{n+1})$

► Free velocity

4: $v_{free} \leftarrow LV_{free}$

5: $\tilde{r} \leftarrow -H^{-1}v_{free}$

6: $r_Y \leftarrow h\sigma_Y$

7: **if** $g_{n+1} > 0$ **then**

8:     **if** $\tilde{r} < -r_Y$ **then** $r_{n+3/2} \leftarrow -r_Y$

9:     **else if** $\tilde{r} > r_Y$ **then** $r_{n+3/2} \leftarrow r_Y$

10:    **else** $r_{n+3/2} \leftarrow \tilde{r}$

11: **end if**

12: **else** $r_{n+3/2} \leftarrow \max(\tilde{r}, -r_Y)$

13: **end if**

14: Velocity $V_{n+3/2} \leftarrow V_{free} + M^{-1}L^T r_{n+3/2}$

► Matrix-free dynamics

<caption>**Algorithm 1:** Time step increment for the perfectly plastic interface</caption>

To illustrate this resolution algorithm, a test case with a single dof problem is considered on Figure 5, with the non-dimensional parameters: $g_0 = 0$, $m = k = 2\sigma_Y = 1$ (the dimensional $m$ would be a mass per unit surface, the dimensional $k$ a stiffness homogeneous to a Young modulus). The characteristic time for the oscillations in the linear regime is $\tau = 2\pi\sqrt{k/m} = 2\pi$. Note that with a single degree-of-freedom problem, the distinction between implicit and explicit scheme is not significant, so there is not a real critical time step in this case and $h = \Delta t$ should be smaller than $\tau$ for accuracy criteria; here $h = 0.03$ is chosen.

The results are depicted on Figure 6. The loading consists in successive 3 forth and back increasing tractions followed by a highly dynamic compression whose maximal compression is maintained constant afterwards. The first traction is not high enough to produce plasticity, which can be seen on the displacement evolution: the interface remains perfectly glued. The second traction produces a plastic evolution: the interface plastically separates, and the interface impulse is constant during perfect plasticity. Finally the plasticity develops even more during the third traction, while the brutal compression allows to plastify in compression, up to the perfect contact activation with a velocity **spatial** jump and an impact impulse.

Since the model formulation acts on velocity and impulse, and since a time stepping approach is used, a residual displacement penetration may occur, related to a precision criteria (indeed, one cannot satisfy exactly simultaneously the displacement constraint, the velocity constraint, and
<figure><img src="image_5.png" /><figcaption>**Figure 4** Solution as the intersection of the behavior (solid line) and reduced dynamics (dotted line).</figcaption></figure>

Indeed, consider the cases separation:

* if $v > 0$, then $r_c = 0$, so $r = r_p = -r_Y$ is computed as before. One therefore gets $v = v_{free} - Hr_Y := \tilde{v}$ that should be positive in this case;

* otherwise, $v = 0$ and the reduced dynamics reads $0 = v_{free} + Hr = (v_{free} - Hr_Y) + H(r_c + r_p + r_Y)$, so $0 = \tilde{v} + H(r_c + r_Y - h\sigma_p)$ for which we should have $r_c \ge 0$ and $0 \le r_Y - h\sigma_p \le 2r_Y$. Now if $\tilde{v} > 0$ this is inconsistent and corresponds to the previous case. Therefore, $\tilde{v} \le 0$ and $r = r_c - h\sigma_p = -H^{-1}v_{free} := \tilde{r}$.

This separation can therefore be summarized as follows:

* if $\tilde{v} > 0$ (which is equivalent to $\tilde{r} < -r_Y$), then $r = -r_Y$,

* else $r = \tilde{r}$.

This can be simplified also, as detailed in Algorithm 1.

**Input** displacement $U_n$, velocity $V_{n+1/2}$

**Output** displacement $U_{n+1}$, velocity $V_{n+3/2}$

1. Displacement $U_{n+1} \leftarrow U_n + hV_{n+1/2}$

2. $g_{n+1} \leftarrow g_0 + LU_{n+1}$

3. $V_{free} \leftarrow V_{n+1/2} + M^{-1}h(f_{ext,n+1} - KU_{n+1})$

4. $v_{free} \leftarrow LV_{free}$

5. $\tilde{r} \leftarrow -H^{-1}v_{free}$

6. $r_Y \leftarrow h\sigma_Y$

7. **if** $g_{n+1} > 0$ **then**

8. **if** $\tilde{r} < -r_Y$ **then** $r_{n+3/2} \leftarrow -r_Y$

9. **else if** $\tilde{r} > r_Y$ **then** $r_{n+3/2} \leftarrow r_Y$

10. **else** $r_{n+3/2} \leftarrow \tilde{r}$

11. **end if**

12. **else** $r_{n+3/2} \leftarrow \max(\tilde{r}, -r_Y)$

13. **end if**

14. $V_{n+3/2} \leftarrow V_{free} + M^{-1}L^T r_{n+3/2}$

► Explicit configuration

► Free velocity

► Matrix-free dynamics

**Algorithm 1** Time step increment for the perfectly plastic interface.

To illustrate this resolution algorithm, a test case with a single-dof problem is considered in Figure 5, with the non-dimensional parameters: $g_0 = 0$, $m = k = 2\sigma_Y = 1$ (the dimensional $m$ would be a mass per unit surface, the dimensional $k$ a stiffness homogeneous to a Young modulus). The characteristic time for the oscillations in the linear regime is $\tau = 2\pi\sqrt{k/m} = 2\pi$.

<figure><img src="image_6.png" /><figcaption>**Figure 5** Sketch of the test case.</figcaption></figure>

Note that with a single degree-of-freedom problem, the distinction between implicit and explicit scheme is not significant, so there is not a real critical time step in this case and $h = \Delta t$ should be smaller than $\tau$ for accuracy criteria; here $h = 0.03$ is chosen.

The results are depicted in Figure 6. The loading consists in successive 3 forth and back increasing tractions followed by a highly dynamic compression whose maximal compression is
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Sketch of the test case</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6:** The numerical solution for the oD test case with the perfectly plastic interface</figcaption></figure>

the discretized relationship between them). Therefore, a relative error indicator can be easily computed with this residual as $\eta = - \min_t g(t) / \max_t |g(t)|$. On the present test case, it takes the value $\eta = 0.43\%$.

## 2.3 Another phenomenological association of contact and perfect plasticity

A serial association of the two previous elementary blocks for the oD test case of figure 7 would lead to a completely different macroscopic model. This last one now may correspond to a perfect contact interface between two bodies, taking into account some plasticity that may develop close to the interface but within the contacting bodies at a lower spatial scale (for instance roughness plasticity). In this case, the macroscopic model may embed the plasticity behavior not into the volumic bodies, but within the surfacic interface. Some examples of this association type may be found, for elasto-plasticity model (Liu et al. 2001; Ghaednia et al. 2017), for elasto-plastic contact at asperity scale (Kogut et al. 2003), or at a larger scale (Wriggers et al. 1993), etc.

This association leads to $r = r_c = r_p$, $v = v_c + v_p$ and $u = u_c + u_p$. The behavior of each block is the same as previously. The reduced dynamics still reads $v = v_{\text{free}} + Hr$.

<figure><img src="image_9.png" /><figcaption>**Figure 7:** Sketch of the test case for the plastic contact</figcaption></figure>

The same method as previously with case separation now leads to the following: if $u > u_p$
<figure><img src="image_7.png" /><figcaption>**Figure 6** Numerical solution for the oD test case with the perfectly plastic interface.</figcaption></figure>

maintained constant afterwards. The first traction is not high enough to produce plasticity, which can be seen on the displacement evolution: the interface remains perfectly glued. The second traction produces a plastic evolution: the interface plastically separates, and the interface impulse is constant during perfect plasticity. Finally the plasticity develops even more during the third traction, while the brutal compression allows to plastify in compression, up to the perfect contact activation with a velocity spatial jump and an impact impulse.

Since the model formulation acts on velocity and impulse, and since a time stepping approach is used, a residual displacement penetration may occur, related to a precision criteria (indeed, one cannot satisfy exactly simultaneously the displacement constraint, the velocity constraint, and the discretized relationship between them). Therefore, a relative error indicator can be easily computed with this residual as $\eta = - \min_t g(t)/\max_t |g(t)|$. On the present test case, it takes the value $\eta = 0.43 \%$.

## 2.3 Another phenomenological association of contact and perfect plasticity

A serial association of the two previous elementary blocks for the oD test case of Figure 7 would lead to a completely different macroscopic model. This last one now may corresponds to a perfect contact interface between two bodies, taking into account some plasticity that may develop close to the interface but within the contacting bodies at a lower spatial scale (for instance roughness plasticity). In this case, the macroscopic model may embed the plasticity behavior not into the solid bodies, but within the surfacic interface. Some examples of this association type may be found, for elasto-plasticity model (Liu et al. 2001; Ghaednia et al. 2017), for elasto-plastic contact at asperity scale (Kogut et al. 2003), or at a larger scale (Wriggers et al. 1993), etc.

This association leads to $r = r_c = r_p$, $v = v_c + v_p$ and $u = u_c + u_p$. The behavior of each block is the same as previously. The reduced dynamics still reads $v = v_{\text{free}} + Hr$. The same method as previously with case separation now leads to the following: if $u > u_p$ (i.e. $u_c > 0$, the contact is inactive) then $r_c = 0 = r_p = r$ (and $v_p = 0$, $v_c = v = v_{\text{free}}$); otherwise,

* if $v_{\text{free}} > 0$ (contact separation) then: $r_c = 0 = r_p = r$ (and $v_p = 0$, $v_c = v = v_{\text{free}}$)

* otherwise (active contact)

- if $v_{\text{free}} > -Hr_Y$ (no plasticity development) then: $r_c = r = r_p = -H^{-1}v_{\text{free}}$ (indeed less than $r_Y$) and $v_c = 0$, $v_p = 0$ so $v = 0$ (rigid interface)
(i.e. $u_c > 0$, the contact is inactive) then $r_c = 0 = r_p = r$ (and $v_p = 0$, $v_c = v = v_{\text{free}}$); otherwise,

* if $v_{\text{free}} > 0$ (contact separation) then: $r_c = 0 = r_p = r$ (and $v_p = 0$, $v_c = v = v_{\text{free}}$)

* otherwise (active contact)

- if $v_{\text{free}} > -Hr_Y$ (no plasticity development) then: $r_c = r = r_p = -H^{-1}v_{\text{free}}$ (indeed less than $r_Y$) and $v_c = 0$, $v_p = 0$ so $v = 0$ (rigid interface)

- else $v_{\text{free}} \le -Hr_Y$ (compressive plastification) so that: $r_p = r_Y$ (and $r_c = r = r_Y$ as well), $v = v_p = v_{\text{free}} + Hr_Y \le 0$ (and $v_c = 0$)

The numerical solution procedure is now described on algorithm 2 (one can note the need of the local plastic displacement $u_p$ as a state internal variable).

**Input:** displacement $U_n$ and $u_{p,n}$, velocity $V_{n+1/2}$ and $v_{p,n+1/2}$

**Output:** displacement $U_{n+1}$ and $u_{p,n+1}$, velocity $V_{n+3/2}$ and $v_{p,n+3/2}$

1: Displacement $U_{n+1} \leftarrow U_n + hV_{n+1/2}$ and $u_{p,n+1} \leftarrow u_{p,n} + hv_{p,n+1/2}$ ▶ Explicit configuration

2: $g_{n+1} \leftarrow g_0 + LU_{n+1}$

3: $V_{\text{free}} \leftarrow V_{n+1/2} + M^{-1}h(f_{\text{ext},n+1} - KU_{n+1})$ ▶ Free velocity

4: $v_{\text{free}} \leftarrow LV_{\text{free}}$

5: $r_Y \leftarrow h\sigma_Y$

6: **if** $g_{n+1} > u_{p,n+1}$ **then** ▶ Inactive contact

7: $r_{n+3/2} \leftarrow 0$

8: $v_{p,n+3/2} \leftarrow 0$

9: **else**

10: $r_{n+3/2} \leftarrow \max(0, \min(-H^{-1}v_{\text{free}}, r_Y))$

11: $v_{p,n+3/2} \leftarrow \min(0, v_{\text{free}} + Hr_Y)$

12: **end if**

13: Velocity $V_{n+3/2} \leftarrow V_{\text{free}} + M^{-1}L^T r_{n+3/2}$ ▶ Matrix-free dynamics

<caption>**Algorithm 2:** Time step increment for the perfectly plastic contact interface</caption>

The results of this test case are depicted in Figure 8. The loading is similar to the opposite of the previous one, but modulated to involve some plasticity in compression, and some contact release. When the traction load is positive, some contact release is obtained, while in compression for a sufficient value, the plastic flow is observed.

<figure><img src="image_10.png" /><figcaption>**Figure 8:** The numerical solution for the oD test case with the contact with perfect plasticity interface</figcaption></figure>

As previously, an error indicator can also be built on a possible residual contact displacement penetration. Now it reads $\eta = - \min_t u_c(t)/\max_t |u(t)|$ and takes the value $\eta = 1\%$.
<figure><img src="image_8.png" /><figcaption>**Figure 7** Sketch of the test case for the plastic contact.</figcaption></figure>

- else $v_{\text{free}} \le -Hr_Y$ (compressive plastification) so that: $r_p = r_Y$ (and $r_c = r = r_Y$ as well),
  $v = v_p = v_{\text{free}} + Hr_Y \le 0$ (and $v_c = 0$)

The numerical solution procedure is now described in Algorithm 2 (one can note the need of
the local plastic displacement $u_p$ as a state internal variable).

**Input** displacement $U_n$ and $u_{p,n}$, velocity $V_{n+1/2}$ and $v_{p,n+1/2}$

**Output** displacement $U_{n+1}$ and $u_{p,n+1}$, velocity $V_{n+3/2}$ and $v_{p,n+3/2}$

1 Displacement $U_{n+1} \leftarrow U_n + hV_{n+1/2}$ and $u_{p,n+1} \leftarrow u_{p,n} + hv_{p,n+1/2}$ ▶ Explicit configuration

2 $g_{n+1} \leftarrow g_0 + LU_{n+1}$
3 $V_{\text{free}} \leftarrow V_{n+1/2} + M^{-1}h(f_{\text{ext},n+1} - KU_{n+1})$ ▶ Free velocity

4 $v_{\text{free}} \leftarrow LV_{\text{free}}$

5 $r_Y \leftarrow h\sigma_Y$
6 **if** $g_{n+1} > u_{p,n+1}$ **then** ▶ Inactive contact

7 $r_{n+3/2} \leftarrow 0$
8 $v_{p,n+3/2} \leftarrow 0$

9 **else**

10 $r_{n+3/2} \leftarrow \max(0, \min(-H^{-1}v_{\text{free}}, r_Y))$

11 $v_{p,n+3/2} \leftarrow \min(0, v_{\text{free}} + Hr_Y)$

12 **end if**

13 $V_{n+3/2} \leftarrow V_{\text{free}} + M^{-1}L^T r_{n+3/2}$ ▶ Matrix-free dynamics

**Algorithm 2** Time step increment for the perfectly plastic contact interface.

The results of this test case are depicted in Figure 8. The loading is similar to the opposite of the previous one, but modulated to involve some plasticity in compression, and some contact release. When the traction load is positive, some contact release is obtained, while in compression for a sufficient value, the plastic flow is observed.

As previously, an error indicator can also be built on a possible residual contact displacement penetration. Now, it reads $\eta = - \min_t u_c(t) / \max_t |u(t)|$ and takes the value $\eta = 1\%$.

3 A rigid, brittle interface with adhesion recovery

An even more non-smooth case is designed herein for testing purpose as an interface behavior
that could be called “non-smooth post-it™”.

As for the previous example, a oD test (i.e. in normal direction of the interface only) is first
considered, and will after be generalized to a 2D interface.

The adhesion recovery is an other mechanism, allowing to recover some adhesion, when a
compression force has been applied. These kind of behavior has already been studied, but in the
case of progressive decohesion and damage evolution (Raous et al. 2006). A perfect brittle version
with non-smooth adhesion recovery is proposed in this section. The non-smoothness sources
reside in the perfect contact in compression, the infinite stiffness in traction before fracture, and
the null energy dissipation of the perfectly brittle fracture. Indeed this limit case is not well posed
for quasi-static problems, since it would lead to an instability leading to instantaneous debonding
of interfaces.
## 3 A rigid, brittle interface with adhesion recovery

An even more non-smooth case is designed herein for testing purpose as an interface behavior that could be called “non-smooth post-it™”.

As for the previous example, a 0D test (i.e. in normal direction of the interface only) is first considered, and will after be generalized to a 2D interface.

The adhesion recovery is an other mechanism, allowing to recover some adhesion, when a compression force has been applied. These kind of behavior has already been studied, but in the case of progressive decohesion and damage evolution (Raous et al. 2006). A perfect brittle version with non-smooth adhesion recovery is proposed in this section. The non-smoothness sources reside in the perfect contact in compression, the infinite stiffness in traction before fracture, and the null energy dissipation of the perfectly brittle fracture. Indeed this limit case is not well posed for quasi-static problems, since it would lead to an instability leading to instantaneous debonding of interfaces.

### 3.1 A normal behavior test with a 0D problem

The considered constitutive behavior is a compound of two different situations: a perfect contact or a gluing behavior with a yield traction load $f_Y$ (or impulse $r_Y = h f_Y$) above which a brittle rupture occurs. The sketch of this behavior is depicted in Figure 9. In compression, the unilateral contact prevents the interpenetration.

The reduced dynamics reads $v = v_{\text{free}} + Hr$. As before, with the diagonal and positivity properties of $H$, the solution is the intersection of the behavior and the reduced dynamics as illustrated on Figure 10, when the gap is null.

<figure><img src="image_11.png" /><figcaption>**Figure 9:** The considered brittle constitutive behavior, as an impulse/displacement spatial jump relationship (the dotted line denotes a brutal jump in the status)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10:** The brittle fracture solution as the intersection of the impulse/velocity behavior (plain line) and the reduced dynamics (dotted line)</figcaption></figure>

When the gap is strictly positive, the fracture has occurred, so $r = 0$.

When the gap is not strictly positive, three different cases may happen:

* either $\tilde{r} \le -r_Y$, then the solution is $r = 0$ (and $v \ge 0$, the fracture occurs);

* or $\tilde{r} > 0$, then the solution is $r = \tilde{r}$ (and $v = 0$, the contact is active);

* else, there are two possible solutions. In this case, the continuity of the previous status allows to decide which one is selected: either the previous status was inactive contact, so the status does not evolve and $r = 0$, or the previous status was active contact or glued, and the status is now glued and $r = \tilde{r}$.

This behavior does nevertheless not take into account that the yield stress $f_Y$ may evolve (and is therefore an internal variable). Indeed, the considered modified model is: when a rupture occurs, $f_Y$ is nullified; when the impulse is positive (compression), $f_Y$ may increase back (this is a recovery of adhesion). An induced simplification also resides in the previous status of the interface: it has been stored in the $f_Y$ internal variable, so that the $f_Y = 0$ case corresponds to the case where a fracture has occurred, and the constitutive behavior then collapses into the
<figure><img src="image_9.png" /><figcaption>**Figure 8** Numerical solution for the oD test case with the contact with perfect plasticity interface.</figcaption></figure>

### 3.1 A normal behavior test with a oD problem

The considered constitutive behavior is a compound of two different situations: a perfect contact or a gluing behavior with a yield traction load $f_Y$ (or impulse $r_Y = hf_Y$) above which a brittle rupture occurs. The sketch of this behavior is depicted in Figure 9. In compression, the unilateral contact prevents the interpenetration.

The reduced dynamics reads $v = v_{\text{free}} + Hr$. As before, with the diagonal and positivity properties of $H$, the solution is the intersection of the behavior and the reduced dynamics as illustrated in Figure 10, when the gap is null. When the gap is strictly positive, the fracture has

<figure><img src="image_10.png" /><figcaption>**Figure 10** Brittle fracture solution as the intersection of the impulse/velocity behavior (solid line) and the reduced dynamics (dotted line).</figcaption></figure>

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

occurred, so $r = 0$.

When the gap is not strictly positive, three different cases may happen:

* either $\tilde{r} \le -r_Y$, then the solution is $r = 0$ (and $v \ge 0$, the fracture occurs);
classical unilateral contact behavior; this allows to a more compact sub-cases selection (only two remaining): indeed, if rupture has occurred, one has $f_Y = r_Y = 0$, so that the first two previous cases cover all the possibilities.

An issue concerns the choice of the variable representing the threshold that should be a material intrinsic parameter. It was previously mentioned that either a yield stress $f_Y$ or an impulse distribution $r_Y$ are used, related with $r_Y = hf_Y$. Therefore, $h$ being a discretization parameter, $f_Y$ and $r_Y$ cannot both be intrinsic material parameters. An answer could a priori be obtained by looking at convergence of the solution when the time step $h$ goes to zero. If a oD problem is considered, loaded with an external force $f_{ext}(t)$, and a material parameter $r_Y$, the impulse solution looks like $r(t) = -hf_{ext}(t)$, therefore $r$ tends toward 0 (and $r/h$ has a finite limit), so whatever the load is, after a certain refinement in time discretization, one cannot have $r \le -r_Y$, so there is an inconsistency. Now consider a oD problem, loaded with an external impulse $r_{ext}(t)$ (so that the equivalent force is $f(t) = r_{ext}(t)/h$). If the material parameter is $f_Y$, since the solution looks like $r(t) = r_{ext}(t)$, after a certain refinement, and whatever the load is, one will always have $r > -r_Y = hf_Y$, so there is also an inconsistency.

Two physical fracture modes may be considered, that are physically different: either due to a smooth force, or due to a shock, so that the failure criteria is a compound: $r < -r_Y$ or $r < -hf_Y$. Herein, another physical model is considered: though the fracture happens brutally in time, the quantity that drives the fracture mechanism could be non-local in time. Indeed, at each interface point the quantity $R(t) = \int_{t-\tau_Y}^{t} f(t)dt$ may be defined, where $f$ is the interface normal stress, which is equivalent to the more general expression $R(t) = \sum_{i,t_i \in [t-\tau_Y, t]} r(t_i)$ homogeneous to an impulse. Its interpretation is that a smooth force may lead to a brittle fracture, provided that is has been applied during a sufficient time duration, and that shocks may also lead to a brittle fracture, provided that they have been repeated sufficiently rapidly or have been sufficiently intense. With this new quantity, the fracture criteria reads $R < -R_Y$, where $R_Y$ is a material parameter homogeneous to an impulse distribution, together with a second intrinsic parameter, namely a characteristic time $\tau_Y$. As mentioned previously, one may consider a variable-in-time $R_Y$. $R_Y$ is nullified when fracture occurs, and can be recovered by normal compression, i.e. when $r > 0$. A consistent criteria could then read: $R_Y(t) = \max(R_Y(t); R(t))$. And finally a maximum recovery could be added with a third material parameter, $R_{Ymax}$ leading to $R_Y(t) = \min(R_{Ymax}; \max(R_Y(t); R(t)))$. Such a modification with a non-locality in time is not related to a smoothing since no modification is made on the instantaneous fracture event, contrary to smoothing evolutions as in (Zhang et al. 2019).

A main implementation difference is due to the non-locality in time: the impulse field $r$ should be stored during a previous duration of $\tau_Y$, as a moving window in time. This indeed leads to a higher storage cost, but does not prevent the explicit matrix-free feature of the algorithm.

A bit more general model can be selected, such as one with an efficiency decrease of the adhesion recovery. Another internal variable is then used: $\alpha$, that measures the efficiency of the adhesion recovery: as an initial condition, $\alpha = 1$ (maximum efficiency). At each fracture event at time $t$, it is then decreased, for instance with the following evolution law: $\alpha(t^+) = \alpha(t^-)e^{-1/n_b}$, where $n_b$ is then an additional material parameter. This is a discrete evolution: the decrease of efficiency is not continuous in time, but occurs at fracture events. The Algorithm 3 implements the non-locality in time and the possible decrease in adhesion recovery efficiency.

It can also be noticed that when this non-smooth brittle fracture occurs, no energy is dissipated in the interface ($r = 0$, so $v^T r = 0$), this could therefore be called a perfectly brittle behavior.

The following test case is implemented: consider again the mechanical system of Figure 5, still with a dimensionless mass $m = 1$, but without stiffness $k = 0$ and with a brittle interface behavior. As a load, some shocks are produced, i.e. an external impulse $r_{ext}(t)$ is prescribed, Figure 11. Concerning the time discretization, 160 time steps of uniform duration $h$ are selected. The initial conditions (null displacement, negative initial velocity) are designed so that a periodic solution is obtained for the first following test.

Several tests are proposed:

* The test 1 consists of using as material parameters $R_{Ymax} = \infty$ (no adhesion limit), $n_b = \infty$
• or $\tilde{r} > 0$, then the solution is $r = \tilde{r}$ (and $v = 0$, the contact is active);

• else, there are two possible solutions. In this case, the continuity of the previous status allows to decide which one is selected: either the previous status was inactive contact, so the status does not evolve and $r = 0$, or the previous status was active contact or glued, and the status is now glued and $r = \tilde{r}$.

This behavior does nevertheless not take into account that the yield stress $f_Y$ may evolve (and is therefore an internal variable). Indeed, the considered modified model is: when a rupture occurs, $f_Y$ is nullified; when the impulse is positive (compression), $f_Y$ may increase back (this is a recovery of adhesion). An induced simplification also resides in the previous status of the interface: it has been stored in the $f_Y$ internal variable, so that the $f_Y = 0$ case corresponds to the case where a fracture has occurred, and the constitutive behavior then collapses into the classical unilateral contact behavior; this allows to a more compact sub-cases selection (only two remaining): indeed, if rupture has occurred, one has $f_Y = r_Y = 0$, so that the first two previous cases cover all the possibilities.

An issue concerns the choice of the variable representing the threshold that should be a material intrinsic parameter. It was previously mentioned that either a yield stress $f_Y$ or an impulse distribution $r_Y$ are used, related with $r_Y = hf_Y$. Therefore, $h$ being a discretization parameter, $f_Y$ and $r_Y$ cannot both be intrinsic material parameters. An answer could a priori be obtained by looking at convergence of the solution when the time step $h$ goes to zero. If a oD problem is considered, loaded with an external force $f_{ext}(t)$, and a material parameter $r_Y$, the impulse solution looks like $r(t) = -hf_{ext}(t)$, therefore $r$ tends toward 0 (and $r/h$ has a finite limit), so whatever the load is, after a certain refinement in time discretization, one cannot have $r \le -r_Y$, so there is an inconsistency. Now consider a oD problem, loaded with an external impulse $r_{ext}(t)$ (so that the equivalent force is $f(t) = r_{ext}(t)/h$). If the material parameter is $f_Y$, since the solution looks like $r(t) = r_{ext}(t)$, after a certain refinement, and whatever the load is, one will always have $r > -r_Y = hf_Y$, so there is also an inconsistency.

Two physical fracture modes may be considered, that are physically different: either due to a smooth force, or due to a shock, so that the failure criteria is a compound: $r < -r_Y$ or $r < -hf_Y$. Herein, another physical model is considered: though the fracture happens brutally in time, the quantity that drives the fracture mechanism could be non-local in time. Indeed, at each interface point the quantity $R(t) = \int_{t-\tau_Y}^{t} f(t)dt$ may be defined, where $f$ is the interface normal stress, which is equivalent to the more general expression $R(t) = \sum_{i,t_i \in [t-\tau_Y,t]} r(t_i)$ homogeneous to an impulse. Its interpretation is that a smooth force may lead to a brittle fracture, provided that is has been applied during a sufficient time duration, and that shocks may also lead to a brittle fracture, provided that they have been repeated sufficiently rapidly or have been sufficiently intense. With this new quantity, the fracture criteria reads $R < -R_Y$, where $R_Y$ is a material parameter homogeneous to an impulse distribution, together with a second intrinsic parameter, namely a characteristic time $\tau_Y$. As mentioned previously, one may consider a variable-in-time $R_Y$. $R_Y$ is nullified when fracture occurs, and can be recovered by normal compression, i.e. when $r > 0$. A consistent criteria could then read: $R_Y(t) = \max(R_Y(t); R(t))$. Finally, a maximum recovery could be added with the third material parameter $R_{Ymax}$ leading to $R_Y(t) = \min(R_{Ymax}; \max(R_Y(t); R(t)))$. Such a modification with a non-locality in time is not related to a smoothing since no modification is made on the instantaneous fracture event, contrary to smoothing evolutions as in (Zhang et al. 2019).

A main implementation difference is due to the non-locality in time: the impulse field $r$ should be stored during a previous duration of $\tau_Y$, as a moving window in time. This indeed leads to a higher storage cost, but does not prevent the explicit matrix-free feature of the algorithm.

A bit more general model can be selected, such as one with an efficiency decrease of the adhesion recovery. Another internal variable $\alpha$ is then used that measures the efficiency of the adhesion recovery: as an initial condition, $\alpha = 1$ (maximum efficiency). At each fracture event at time $t$, it is then decreased, for instance with the evolution law $\alpha(t^+) = \alpha(t^-)e^{-1/n_b}$ where $n_b$ is then an additional material parameter. This is a discrete evolution: the decrease of efficiency is not continuous in time, but occurs at fracture events. Algorithm 3 implements the non-locality in time and the possible decrease in adhesion recovery efficiency.

It can also be noticed that when this non-smooth brittle fracture occurs, no energy is
**Input:** $U_n$ and $V_{n+1/2}$ as well as interface status $R_{Y,n+1/2}$ and $\alpha_{n+1/2}$

**Output:** $U_{n+1}$ and $V_{n+3/2}$ as well as new interface status $R_{Y,n+3/2}$ and $\alpha_{n+3/2}$

1: Displacement $U_{n+1} \leftarrow U_n + hV_{n+1/2}$ ▶ Explicit configuration

2: $g_{n+1} \leftarrow g_0 + LU_{n+1}$

3: $V_{free} \leftarrow V_{n+1/2} + M^{-1}h(f_{ext,n+1} - KU_{n+1})$ ▶ Free velocity

4: $v_{free} \leftarrow LV_{free}$

5: **if** $g_{n+1} > 0$ **then** ▶ Inactive interaction

6: $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$

7: $R_{Y,n+3/2} \leftarrow 0$

8: $r_{n+3/2} \leftarrow 0$

9: **else**

10: $\tilde{r} \leftarrow -H^{-1}v_{free}$

11: $\tilde{R} \leftarrow \tilde{r} + \sum_{i,t_i \in [t_{n+3/2} - \tau_Y, t_{n+1/2}]} r_i$

12: $\tilde{R}_Y \leftarrow \min(R_{Ymax}; \max(R_{Y,n+1/2}; \alpha_{n+1/2}\tilde{R}))$

13: **if** $\tilde{R} > -\tilde{R}_Y$ **then** ▶ Glued interaction

14: $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$

15: $R_{Y,n+3/2} \leftarrow \tilde{R}_Y$

16: $r_{n+3/2} \leftarrow \tilde{r}$

17: **else**

18: **if** $\tilde{R}_Y > 0$ **then** ▶ Brittle fracture occurs

19: $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}e^{-1/n_b}$ ▶ Adhesion efficiency decrease

20: **else**

21: $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$

22: **end if**

23: $R_{Y,n+3/2} \leftarrow 0$

24: $r_{n+3/2} \leftarrow <\tilde{r}>_+$

25: **end if**

26: **end if**

27: Velocity $V_{n+3/2} \leftarrow V_{free} + M^{-1}L^T r_{n+3/2}$ ▶ Matrix-free dynamics

<figure><img src="image_13.png" /><figcaption>**Figure 11:** The non-smooth loading for the oD test cases</figcaption></figure>

<table><caption>**Table 1:** Interface non-dimensional parameters for the different test cases</caption><thead><tr><th>test</th><th>1</th><th>2</th><th>3</th><th>3D</th></tr></thead><tbody><tr><td>efficiency decrease n<sub>b</sub></td><td>&infin;</td><td>10</td><td>&infin;</td><td>&infin;</td></tr><tr><td>adhesion limit R<sub>Ymax</sub></td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>10<sup>-4</sup></td></tr><tr><td>delay effect &tau;<sub>Y</sub></td><td>0</td><td>0</td><td>0.08</td><td>0.2</td></tr><tr><td>normal to tangent influence &lambda;</td><td>-</td><td>-</td><td>-</td><td>0.8</td></tr><tr><td>friction coefficient &mu;</td><td>-</td><td>-</td><td>-</td><td>0.1</td></tr></tbody></table>
**Input** $U_n$ and $V_{n+1/2}$ as well as interface status $R_{Y,n+1/2}$ and $\alpha_{n+1/2}$

**Output** $U_{n+1}$ and $V_{n+3/2}$ as well as new interface status $R_{Y,n+3/2}$ and $\alpha_{n+3/2}$

1. Displacement $U_{n+1} \leftarrow U_n + hV_{n+1/2}$

2. $g_{n+1} \leftarrow g_0 + LU_{n+1}$

3. $V_{\text{free}} \leftarrow V_{n+1/2} + M^{-1}h(f_{\text{ext},n+1} - KU_{n+1})$

4. $v_{\text{free}} \leftarrow LV_{\text{free}}$

5. **if** $g_{n+1} > 0$ **then**

6. $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$

7. $R_{Y,n+3/2} \leftarrow 0$

8. $r_{n+3/2} \leftarrow 0$

9. **else**

10. $\tilde{r} \leftarrow -H^{-1}v_{\text{free}}$

11. $\tilde{R} \leftarrow \tilde{r} + \sum_{i,t_i \in [t_{n+3/2} - \tau_Y, t_{n+1/2}]} r_i$

12. $\tilde{R}_Y \leftarrow \min(R_{Y,\text{max}}; \max(R_{Y,n+1/2}; \alpha_{n+1/2}\tilde{R}))$

13. **if** $\tilde{R} > -\tilde{R}_Y$ **then**

14. $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$

15. $R_{Y,n+3/2} \leftarrow \tilde{R}_Y$

16. $r_{n+3/2} \leftarrow \tilde{r}$

17. **else**

18. **if** $\tilde{R}_Y > 0$ **then**

19. $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}e^{-1/n_b}$

20. **else**

21. $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$

22. **end if**

23. $R_{Y,n+3/2} \leftarrow 0$

24. $r_{n+3/2} \leftarrow \langle \tilde{r} \rangle_+$

25. **end if**

26. **end if**

27. Velocity $V_{n+3/2} \leftarrow V_{\text{free}} + M^{-1}L^T r_{n+3/2}$

► Explicit configuration

► Free velocity

► Inactive interaction

► Glued interaction

► Brittle fracture occurs

► Adhesion efficiency decrease

► Matrix-free dynamics

**Algorithm 3** Time step increment for the brittle fracture interface (normal behavior only).

dissipated in the interface ($r = 0$, so $v^T r = 0$), this could therefore be called a perfectly brittle behavior.

The following test case is implemented: consider again the mechanical system of Figure 5, still with a dimensionless mass $m = 1$, but without stiffness $k = 0$ and with a brittle interface behavior. As a load, some shocks are produced, i.e. an external impulse $r_{\text{ext}}(t)$ is prescribed, Figure 11. Concerning the time discretization, 160 time steps of uniform duration $h$ are selected.

<figure><img src="image_12.png" /><figcaption>**Figure 11** Non-smooth loading for the oD test cases.</figcaption></figure>

As detailed below, various tests are proposed and the initial conditions (vanishing displacement, negative initial velocity) are designed so that a periodic solution is obtained for Test 1.

* Test 1 consists of using as material parameters $R_{Y,\text{max}} = \infty$ (no adhesion limit), $n_b = \infty$ (no efficiency decrease), $\tau_Y = 0$ (no delay effect). For these values, a cyclic response is produced, see Figure 12. A first impulse $r(t = h) > 0$ is produced due to the initial condition that engenders a compressive shock and increases the yield impulse $R_Y$; then a negative impulse happens opposed to the external first positive shock $r_{\text{ext}}$ which is not large enough to produce a fracture; the second external impulse is large enough for a fracture to happen (the impulse $r$ is therefore null, as the yield impulse $R_Y$), and to induce kinetic energy of detachment. The third external impulse
(no efficiency decrease), $\tau_Y = 0$ (no delay effect). For these values, a cyclic response is produced, Figure 12. A first impulse $r(t = h) > 0$ is produced due to the initial condition that engenders a compressive shock and increases the yield impulse $R_Y$; then a negative impulse happens opposed to the external first positive shock $r_{ext}$ which is not large enough to produce a fracture; the second external impulse is large enough for a fracture to happen (the impulse $r$ is therefore null, as the yield impulse $R_Y$), and to induce kinetic energy of detachment. The third external impulse in negative and has been adjusted to get back the mass into contact with a positive impulse $r$ identical to the initial one, so cycling is obtained afterwards. Thanks to good discrete energy conservation of the integration scheme, cycles are preserved at long time duration.

* For the test 2, only the efficiency decrease is changed to $n_b = 10$. The solution is no more cyclic, Figure 13. Indeed the yield impulse reconstruction $R_Y$ decreases up to the point when the first type of external impulse can produce the fracture, so the impulse and velocity solutions are no more cyclic.

* The test 3 uses the same parameters as for test 1, except that a delay effect is used with $\tau_Y = 0.08$ (actually 64 time steps). Though the solution differs from test 1 for the yield $R_Y$, Figure 14, due to the time delay that is a moving window in time, the associated material behavior is conservative, allows to recover the same solution for the other quantities.

<figure><img src="image_14.png" /><figcaption>Figure 12: The numerical solution for the first oD brittle test case</figcaption></figure>

To check time-convergence, an additional test is used, based on the previous test 3, but with an additional smooth loading $f_{ext}(t) = -F_d \sin(\frac{4}{5}\pi t/T_c)$ with a non-dimensional amplitude $F_d = 2$, $T_c = 0.04$ being the non-dimensional duration of one cycle. Different time steps are used, and an over-killing one is used to obtain a reference solution ($h_{ref} \approx 7.63 \ 10^{-8}$). Since for instance the velocity is discontinuous, a Hausdorff distance (Moreau 1978; Acary 2012) is used to get an error with respect to this reference, both for the displacement and the velocity solutions. Note that normalizations of evolution graphs are needed for the metric of such a norm: the time scale is rescaled as $t/T$, $T$ being the total studied time duration and both displacement and velocity are rescaled with respect to their maximum absolute values of the reference solutions. In such a way the Hausdorff distance is an absolute error measure and is reported on Figure 15.

Apart from checking the convergence, and as expected for this central-difference-like scheme for non-smooth problems, a convergence rate of $O(h)$ is obtained, see e.g. (Chen et al. 2013).

## 3.2 3D structures and surfacic interfaces

To start with, the interface behavior is considered with a normal $f_N$ and a tangential force distribution (or stress) $f_T$, together with a yield stress $\sigma_Y$ at an interface point. No non-locality in
<table><caption>**Table 1** Interface non-dimensional parameters for the different test cases.</caption><thead><tr><th>test</th><th>1</th><th>2</th><th>3</th><th>3D</th></tr></thead><tbody><tr><td>efficiency decrease n<sub>b</sub></td><td>∞</td><td>10</td><td>∞</td><td>∞</td></tr><tr><td>adhesion limit R<sub>Y max</sub></td><td>∞</td><td>∞</td><td>∞</td><td>10<sup>-4</sup></td></tr><tr><td>delay effect τ<sub>Y</sub></td><td>0</td><td>0</td><td>0.08</td><td>0.2</td></tr><tr><td>normal to tangent influence λ</td><td>—</td><td>—</td><td>—</td><td>0.8</td></tr><tr><td>friction coefficient μ</td><td>—</td><td>—</td><td>—</td><td>0.1</td></tr></tbody></table>

in negative and has been adjusted to get back the mass into contact with a positive impulse *r* identical to the initial one, so cycling is obtained afterwards. Thanks to good discrete energy conservation of the integration scheme, cycles are preserved at long time duration.

* For Test 2, only the efficiency decrease is changed to $n_b = 10$. The solution is no more cyclic, see Figure 13. Indeed the yield impulse reconstruction $R_Y$ decreases up to the point when the first type of external impulse can produce the fracture, so the impulse and velocity solutions are no more cyclic.

* Test 3 uses the same parameters as for test 1, except that a delay effect is used with $τ_Y = 0.08$ (actually 64 time steps). Though the solution differs from test 1 for the yield $R_Y$, see Figure 14, due to the time delay that is a moving window in time, the associated material behavior is conservative, allows to recover the same solution for the other quantities.

<figure><img src="image_13.png" /><figcaption>**Figure 12** Numerical solution for the first 0D brittle test case.</figcaption></figure>

To check time-convergence, an additional test is used, based on the previous test 3, but with an additional smooth loading $f_{ext}(t) = -F_d \sin(\frac{4}{5}\pi t/T_c)$ with a non-dimensional amplitude $F_d = 2$, $T_c = 0.04$ being the non-dimensional duration of one cycle. Different time steps are used, and an over-killing one is used to obtain a reference solution ($h_{ref} \approx 7.63 \ 10^{-8}$). Since for instance the velocity is discontinuous, a Hausdorff distance (Moreau 1978; Acary 2012) is used to get an error with respect to this reference, both for the displacement and the velocity solutions. Note that normalizations of evolution graphs are needed for the metric of such a norm: the time scale is rescaled as $t/T$, $T$ being the total studied time duration and both displacement and velocity are rescaled with respect to their maximum absolute values of the reference solutions. In such a way the Hausdorff distance is an absolute error measure and is reported in Figure 15.

Apart from checking the convergence, and as expected for this central-difference-like scheme
<figure><img src="image_15.png" /><figcaption>Figure 13: The numerical solution for the second oD brittle test case</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: The numerical solution for the third oD brittle test case</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 13** Numerical solution for the second oD brittle test case.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14** Numerical solution for the third oD brittle test case.</figcaption></figure>

for non-smooth problems, a convergence rate of $O(h)$ is obtained, see e.g. (Chen et al. 2013).

## 3.2 3D structures and surfacic interfaces

To start with, the interface behavior is considered with a normal $f_N$ and a tangential force distribution (or stress) $f_T$, together with a yield stress $\sigma_Y$ at an interface point. No non-locality in time are considered in a first step.

Since fracture may happen with a shear stress, a threshold criteria should be used. A proposed
<figure><img src="image_17.png" /><figcaption>Figure 15: The loading with forces and impulses, the numerical convergence of the solution with respect to time discretization and the reference solution for velocity</figcaption></figure>

time are considered in a first step.

Since fracture may happen with a shear stress, a threshold criteria should be used. A proposed equivalent stress is defined as

$$ \sigma(f_N, \underline{f}_T) = \sqrt{\langle f_N \rangle_+^2 + \|\underline{f}_T\|^2 / \lambda^2} \quad (12) $$

so only the tension normal part (hence the positive part $\langle f_N \rangle_+$ is involved) acts on the fracture, and where $\lambda$ is a weight for the influence of the tangential part, Figure 16. When $\sigma(f_N, \underline{f}_T) < \sigma_Y$ no fracture occurs and the interface status is ‘glued’, i.e. the velocity spatial jump is $v = 0$; when $\sigma(f_N, \underline{f}_T) > \sigma_Y$ a fracture has occurred and the interface status is ‘contact’, i.e. unilateral contact with Coulomb tangential friction, with a friction coefficient $\mu$, that may be either inactive ($v_N > 0, f_N = 0, \underline{f}_T = 0$) or sliding ($v_T \neq 0, f_N > 0, \|\underline{f}_T\| = \mu f_N$) or sticking ($v_T = 0, f_N > 0, \|\underline{f}_T\| \le \mu f_N$). In this considered model, no smooth transition occurs between the glued and contact status, to keep the behavior as non smooth as possible for testing the proposed resolution algorithms.

It can be noticed that when $\sigma_Y = 0$, the threshold collapses to the positive semi-axis $f_N < 0$. When $\lambda = \infty$ the tangential direction has no influence on the fracture, the normal behavior is identical as in the oD test case. When $R_{Ymax} = 0$, the behavior reduces to a classical contact with Coulomb friction.

<figure><img src="image_18.png" /><figcaption>Figure 16: Yield function for 2D fracture criteria</figcaption></figure>

Since with a 2D interface, tangential components are involved as vector quantities, underlying has been used to differentiate them from scalar quantities. Note however that in the following, when dealing with discretized quantities as column vectors, no underlying will be used anymore.

The resolution algorithm 4 is then built on the same framework as previously. In particular, the velocity-impulse formulation is used, with $r_Y = h\sigma_Y, r_N = -hf_N, r_T = -hf_T$, but as 3D cases discretized by finite elements are concerned, one has to add a loop on the nodes of the interface.
<figure><img src="image_16.png" /><figcaption>**Figure 15** Loading with forces and impulses: numerical convergence of the solution with respect to time discretization and the reference solution for velocity.</figcaption></figure>

equivalent stress is defined as

$$ \sigma(f_N, f_T) = \sqrt{\langle f_N \rangle_+^2 + \|f_T\|^2 / \lambda^2} \qquad (12) $$

so only the tension normal part (hence the positive part $\langle f_N \rangle_+$ is involved) acts on the fracture, and where $\lambda$ is a weight for the influence of the tangential part, Figure 16. When $\sigma(f_N, f_T) < \sigma_Y$ no fracture occurs and the interface status is ‘glued’, i.e. the velocity spatial jump is $v = 0$; when $\sigma(f_N, f_T) > \sigma_Y$ a fracture has occurred and the interface status is ‘contact’, i.e. unilateral contact with Coulomb tangential friction, with a friction coefficient $\mu$, that may be either inactive ($v_N > 0$, $f_N = 0$, $f_T = 0$) or sliding ($v_T \neq 0$, $f_N > 0$, $\|f_T\| = \mu f_N$) or sticking ($v_T = 0$, $f_N > 0$, $\|f_T\| \le \mu f_N$). In this considered model, no smooth transition occurs between the glued and contact status, to keep the behavior as non smooth as possible for testing the proposed resolution algorithms.

It can be noticed that when $\sigma_Y = 0$, the threshold collapses to the positive semi-axis $f_N < 0$. When $\lambda = \infty$ the tangential direction has no influence on the fracture, the normal behavior is identical as in the oD test case. When $R_{Ymax} = 0$, the behavior reduces to a classical contact with Coulomb friction. Since with a 2D interface, tangential components are involved as vector

<figure><img src="image_17.png" /><figcaption>**Figure 16** Yield function for 2D fracture criteria.</figcaption></figure>

quantities, underlying has been used to differentiate them from scalar quantities. Note however that in the following, when dealing with discretized quantities as column vectors, no underlying will be used anymore.

Algorithm 4 is then built on the same framework as previously. In particular, the velocity-impulse formulation is used, with $r_Y = h\sigma_Y$, $r_N = -hf_N$, and $r_T = -hf_T$, but as 3D cases discretized by finite elements are concerned, one has to add a loop on the nodes of the interface.
The non-locality in time and the adhesion recovery decrease are also trivially extended to these cases.

Moreover, since the local quantities on the interface are impulse distributions, an additional integration on the interface has to be performed. To keep these quantities at nodal position on the interface, a nodal integration is therefore selected, the integration weight at node k is $S^{(k)}$. To simplify the notations, their assembling on the full set of dofs will be denoted with the diagonal entry matrix $S$, with duplicated entries on the 3 spatial directions and 3 zeros entries on the nodes that are not located on the interface. For 3D cases, the dynamics therefore reads

$$V_{n+3/2} = V_{\text{free}} + M^{-1}SL^T r_{n+3/2} \quad (13)$$

$$= V_{\text{free}} + M^{-1}SL_N^T r_{N,n+3/2} + M^{-1}SL_T^T r_{T,n+3/2} \quad (14)$$

where the mapping operators $L_N^T$ and $L_T^T$ contain respectively the normal to the interface, two tangent directions at the interface, and the mapping to the global dofs. With the previously mentioned properties of matrix $M$, still valid in 3D, and similar properties for $S$, one has $L_N M^{-1} S L_T^T = 0$ so that the reduced dynamics reads firstly: $v_N = v_{Nfree} + H_N r_N$, where $H_N = L_N M^{-1} S L_N^T$, $v_{Nfree} = L_N V_{\text{free}}$, and secondly: $v_T = v_{Tfree} + H_T r_T$, where $H_T = L_T M^{-1} S L_T^T$, $v_{Tfree} = L_T V_{\text{free}}$. The Delassus operators $H_N$ and $H_T$ are still definite positive, diagonal and spherical per node entry.

Concerning the energy conservation, no energy is lost during fracture, but kinetic energy of back impacting node is lost at impact (no Newton restitution coefficient is used). This is a drawback of the discretized problem (the energy loss is driven to zero when mesh size is refined), but can be dealt with using an advanced modeling of mass matrix: the so-called redistributed or singular mass technique (Khenous et al. 2008; Tkachuk et al. 2013; Dabaghi et al. 2019; Di Stasio et al. 2021). There is also a physical source of energy loss, due to the friction between fractured surfaces, driven by the friction coefficient $\mu$.

### 3.3 3D test case

The 3D test case is modeled with classical finite elements, the geometry is given on Figure 17. It consists on a one-fourth of a tube slice (the internal radius is $R_1 = 1$, the thickness $L_1 = 0.5$ and the length $L_2 = 0.5$) meshed with 5712 cubic elements with 4 nodes each, leading to 20475 dof. The non-dimensional volumic parameters are: a Young modulus $E = 1$, a Poisson coefficient $\nu = 0.3$ and a density of $\rho = 1$. Two interfaces with the rigid foundation are considered, on the two rectangular ends of the tube quarter. Both have the same behavior, whose parameters are given in Table 1. With the previous spatial discretization, there are 390 interface nodes.

The boundary conditions are a symmetry plane, and a shear stress on the other end face (a quadratic radial distribution, whose resultant value is depicted on Figures 18 (a) and (c)) engendering a non-uniform solicitation both with normal and tangential components on the interfaces. Two time evolutions of the loading are considered, Figure 18 (a) and (c), the first one being less smooth than the second and involving more dynamical effects.

The initial conditions are a null displacement and velocity, and an initial adhesion threshold $R_Y(t = 0) = 10^{-4}$. The studied time interval is $T = 16$, with a time step $h = T/8000 = 0.002$, while the estimated critical time step is $t_c \approx 0.047$.

One can note that there could be some unstable initial determination of the evolution of the state variables, subjected to numerical precision issues when one uses at initial condition both an initial gap $g_0 = 0$ and an initial velocity $V_0 = 0$. Indeed, this corresponds to a grazing contact for which the distinction between active and inactive contact is undetermined. To fix this issue, a simple solution is to initialize the initial gap with a 'small' value (larger than the numerical precision), e.g. $g_0 = -10^{-10}$, so that the initial state is active contact.

If $\lambda = \infty$ (no influence of tangential component), no failure occurs for this test case, since the interfaces are subjected to a normal compression. On the other hand, with the selected value of $\lambda$ of Table 1, shearing is taken into account and the evolution of the fraction of cracked surface of the interfaces are depicted on Figure 18 (b) and (d), for the two loading time evolutions. Indeed, the smooth loading leads to a delayed fracture, since the loading increases slower than for the first case. Moreover, the fracture curves are smoother as well, since there are less forward and
The non-locality in time and the adhesion recovery decrease are also trivially extended to these cases.

Moreover, since the local quantities on the interface are impulse distributions, an additional integration on the interface has to be performed. To keep these quantities at nodal position on the interface, a nodal integration is therefore selected, the integration weight at node k is $S^{(k)}$. To simplify the notations, their assembling on the full set of dofs will be denoted with the diagonal entry matrix S, with duplicated entries on the 3 spatial directions and 3 zeros entries on the nodes that are not located on the interface. For 3D cases, the dynamics therefore reads

$$V_{n+3/2} = V_{\text{free}} + M^{-1}SL^{\top}r_{n+3/2} \quad (13)$$

$$= V_{\text{free}} + M^{-1}SL_{N}^{\top}r_{N,n+3/2} + M^{-1}SL_{T}^{\top}r_{T,n+3/2} \quad (14)$$

where the mapping operators $L_N^\top$ and $L_T^\top$ contain respectively the normal to the interface, two tangent directions at the interface, and the mapping to the global dofs. With the previously mentioned properties of matrix M, still valid in 3D, and similar properties for S, one has $L_N M^{-1} S L_T^\top = 0$ so that the reduced dynamics reads firstly: $v_N = v_{Nfree} + H_N r_N$, where $H_N = L_N M^{-1} S L_N^\top$, $v_{Nfree} = L_N V_{\text{free}}$, and secondly: $v_T = v_{Tfree} + H_T r_T$, where $H_T = L_T M^{-1} S L_T^\top$, $v_{Tfree} = L_T V_{\text{free}}$. The Delassus operators $H_N$ and $H_T$ are still definite positive, diagonal and spherical per node entry.

Concerning the energy conservation, no energy is lost during fracture, but kinetic energy of back impacting node is lost at impact (no Newton restitution coefficient is used). This is a drawback of the discretized problem (the energy loss is driven to zero when mesh size is refined), but can be dealt with using an advanced modeling of mass matrix: the so-called redistributed or singular mass technique (Khenous et al. 2008; Tkachuk et al. 2013; Dabaghi et al. 2019; Di Stasio et al. 2021). There is also a physical source of energy loss, due to the friction between fractured surfaces, driven by the friction coefficient $\mu$.

### 3.3 3D test case

The 3D test case is modeled with classical finite elements, the geometry is given in Figure 17. It consists on a one-fourth of a tube slice (the internal radius is $R_1 = 1$, the thickness $L_1 = 0.5$ and the length $L_2 = 0.5$) meshed with 5712 cubic elements with 4 nodes each, leading to 20475 dof. The non-dimensional parameters are: a Young modulus $E = 1$, a Poisson coefficient $\nu = 0.3$ and a density of $\rho = 1$. Two interfaces with the rigid foundation are considered, on the two rectangular ends of the tube quarter. Both have the same behavior, whose parameters are given in Table 1. With the previous spatial discretization, there are 390 interface nodes. The boundary conditions

<figure><img src="image_18.png" /><figcaption>**Figure 17** Geometry of the 3D test case with the vertical shear stress on the face end.</figcaption></figure>

are a symmetry plane, and a shear stress on the other end face (a quadratic radial distribution, whose resultant value is depicted in Figures 18(a,c)) engendering a non-uniform solicitation both with normal and tangential components on the interfaces. Two time histories of the loading are
**Input:** $U_n$ and $V_{n+1/2}$ as well as interface status fields $R_{Y,n+1/2}$ and $\alpha_{n+1/2}$

**Output:** $U_{n+1}$ and $V_{n+3/2}$ as well as new interface status fields $R_{Y,n+3/2}$ and $\alpha_{n+3/2}$

1: Displacements $U_{n+1} \leftarrow U_n + hV_{n+1/2}$ ▶ Explicit configuration

2: $g_{n+1} \leftarrow g_0 + LU_{n+1}$

3: $V_{free} \leftarrow V_{n+1/2} + M^{-1}h(f_{ext,n+1} - KU_{n+1})$ ▶ Free velocity

4: **for each interface point k do**
5: &nbsp;&nbsp;&nbsp;&nbsp;$v_{Nfree} \leftarrow L_N^{(k)} V_{free}$ ▶ With the local to node k operator

6: &nbsp;&nbsp;&nbsp;&nbsp;$v_{Tfree} \leftarrow L_T^{(k)} V_{free}$
7: &nbsp;&nbsp;&nbsp;&nbsp;**if** $g_{n+1}^k > 0$ **then** ▶ Inactive interaction

8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$

9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$R_{Y,n+3/2} \leftarrow 0$

10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$r_{N,n+3/2} \leftarrow 0$

11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$r_{T,n+3/2} \leftarrow 0$

12: **else**
13: &nbsp;&nbsp;&nbsp;&nbsp;$\tilde{r}_N \leftarrow -H_N^{-1}v_{Nfree}$
14: &nbsp;&nbsp;&nbsp;&nbsp;$\tilde{r}_T \leftarrow -H_T^{-1}v_{Tfree}$
15: &nbsp;&nbsp;&nbsp;&nbsp;$\tilde{R}_N \leftarrow \tilde{r}_N + \sum_{i,t_i \in [t_{n+3/2}-\tau_Y, t_{n+1/2}]} r_{N,i}$
16: &nbsp;&nbsp;&nbsp;&nbsp;$\tilde{R}_T \leftarrow \tilde{r}_T + \sum_{i,t_i \in [t_{n+3/2}-\tau_Y, t_{n+1/2}]} r_{T,i}$
17: &nbsp;&nbsp;&nbsp;&nbsp;$\tilde{R}_Y \leftarrow \min(R_{Y,max}; \max(R_{Y,n+1/2}; \alpha_{n+1/2}\tilde{R}_N))$

18: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\sigma(-\tilde{R}_N, -\tilde{R}_T) < \tilde{R}_Y$ **then** ▶ Glued interaction
19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$R_{Y,n+3/2} \leftarrow \tilde{R}_Y$
21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$r_{N,n+3/2} \leftarrow \tilde{r}_N$
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$r_{T,n+3/2} \leftarrow \tilde{r}_T$

23: **else**
24: &nbsp;&nbsp;&nbsp;&nbsp;**if** $R_{Y,n+1/2} > 0$ **then** ▶ Brittle fracture occurs
25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}e^{-1/n_b}$ ▶ Adhesion efficiency decrease

26: **else**
27: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$

28: **end if**
29: &nbsp;&nbsp;&nbsp;&nbsp;$R_{Y,n+3/2} \leftarrow 0$
30: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\tilde{r}_N > 0$ **then** ▶ Active contact

31: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$r_{N,n+3/2} \leftarrow \tilde{r}_N$
32: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\|\tilde{r}_T\| > \mu r_{N,n+3/2}$ **then** ▶ Sliding

33: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$r_{T,n+3/2} \leftarrow \mu r_{N,n+3/2} \frac{\tilde{r}_T}{\|\tilde{r}_T\|}$
34: **else** ▶ Sticking

35: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$r_{T,n+3/2} \leftarrow \tilde{r}_T$

36: **end if**
37: **else** ▶ Released contact
38: &nbsp;&nbsp;&nbsp;&nbsp;$r_{N,n+3/2} \leftarrow 0$

39: &nbsp;&nbsp;&nbsp;&nbsp;$r_{T,n+3/2} \leftarrow 0$
40: **end if**
41: **end if**
42: **end if**
43: Assemble $r_{N,n+3/2}$ and $r_{T,n+3/2}$ in the right-hand-side $r_{n+3/2}$

44: **end for**
45: Velocity $V_{n+3/2} \leftarrow V_{free} + M^{-1}SL^T r_{n+3/2}$ ▶ Matrix-free dynamics

<caption>**Algorithm 4:** Time step increment for the 2D brittle fracture interface</caption>
Input $U_n$ and $V_{n+1/2}$ as well as interface status fields $R_{Y,n+1/2}$ and $\alpha_{n+1/2}$

Output $U_{n+1}$ and $V_{n+3/2}$ as well as new interface status fields $R_{Y,n+3/2}$ and $\alpha_{n+3/2}$

1 Displacements $U_{n+1} \leftarrow U_n + hV_{n+1/2}$ ▶ Explicit configuration

2 $g_{n+1} \leftarrow g_0 + LU_{n+1}$
3 $V_{\text{free}} \leftarrow V_{n+1/2} + M^{-1}h(f_{\text{ext},n+1} - KU_{n+1})$ ▶ Free velocity

4 for each interface point $k$ do
5 $v_{N\text{free}} \leftarrow L_N^{(k)} V_{\text{free}}$ ▶ With the local to node $k$ operator

6 $v_{T\text{free}} \leftarrow L_T^{(k)} V_{\text{free}}$
7 if $g_{n+1}^k > 0$ then ▶ Inactive interaction

8 $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$
9 $R_{Y,n+3/2} \leftarrow 0$
10 $r_{N,n+3/2} \leftarrow 0$
11 $r_{T,n+3/2} \leftarrow 0$
12 else

13 $\tilde{r}_N \leftarrow -H_N^{-1}v_{N\text{free}}$
14 $\tilde{r}_T \leftarrow -H_T^{-1}v_{T\text{free}}$
15 $\tilde{R}_N \leftarrow \tilde{r}_N + \sum_{i,t_i \in [t_{n+3/2}-\tau_Y, t_{n+1/2}]} r_{N,i}$
16 $\tilde{R}_T \leftarrow \tilde{r}_T + \sum_{i,t_i \in [t_{n+3/2}-\tau_Y, t_{n+1/2}]} r_{T,i}$
17 $\tilde{R}_Y \leftarrow \min(R_{Y,\text{max}}; \max(R_{Y,n+1/2}; \alpha_{n+1/2}\tilde{R}_N))$

18 if $\sigma(-\tilde{R}_N, -\tilde{R}_T) < \tilde{R}_Y$ then ▶ Glued interaction
19 $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$
20 $R_{Y,n+3/2} \leftarrow \tilde{R}_Y$
21 $r_{N,n+3/2} \leftarrow \tilde{r}_N$
22 $r_{T,n+3/2} \leftarrow \tilde{r}_T$

23 else
24 if $R_{Y,n+1/2} > 0$ then ▶ Brittle fracture occurs
25 $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}e^{-1/n_b}$ ▶ Adhesion efficiency decrease

26 else
27 $\alpha_{n+3/2} \leftarrow \alpha_{n+1/2}$

28 end if
29 $R_{Y,n+3/2} \leftarrow 0$
30 if $\tilde{r}_N > 0$ then ▶ Active contact

31 $\tilde{r}_{N,n+3/2} \leftarrow \tilde{r}_N$
32 if $\|\tilde{r}_T\| > \mu r_{N,n+3/2}$ then ▶ Sliding

33 $r_{T,n+3/2} \leftarrow \mu r_{N,n+3/2} \frac{\tilde{r}_T}{\|\tilde{r}_T\|}$
34 else ▶ Sticking

35 $r_{T,n+3/2} \leftarrow \tilde{r}_T$
36 end if

37 else
38 $r_{N,n+3/2} \leftarrow 0$
39 $r_{T,n+3/2} \leftarrow 0$

40 end if
41 end if
42 end if

43 Assemble $r_{N,n+3/2}$ and $r_{T,n+3/2}$ in the right-hand-side $r_{n+3/2}$

44 end for
45 Velocity $V_{n+3/2} \leftarrow V_{\text{free}} + M^{-1}SL^T r_{n+3/2}$ ▶ Matrix-free dynamics

**Algorithm 4** Time step increment for the 2D brittle fracture interface.

considered, Figures 18(a,c), the first one being less smooth than the second and involving more dynamical effects.

The initial conditions are a null displacement and velocity, and an initial adhesion threshold $R_Y(t=0) = 10^{-4}$. The studied time interval is $T=16$, with a time step $h = T/8000 = 0.002$, while the estimated critical time step is $t_c \approx 0.047$.

One can note that there could be some unstable initial determination of the evolution of the state variables, subjected to numerical precision issues when one uses at initial condition both an initial gap $g_0 = 0$ and an initial velocity $V_0 = 0$. Indeed, this corresponds to a grazing contact for which the distinction between active and inactive contact is undetermined. To fix this issue, a simple solution is to initialize the initial gap with a 'small' value (larger than the numerical precision), e.g. $g_0 = -10^{-10}$, so that the initial state is active contact.
<figure><img src="image_19.png" /><figcaption>**Figure 17:** Geometry of the 3D test case, with the vertical shear stress on the face end</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18:** (a) test case 1, and (c) test case 2: time loading evolutions; (b) test case 1, and (d) test case 2: evolutions of the fraction of cracked surfaces ($c = 1$); (e) test case 2: for space-time convergence ($c = 3$); (f) test case 2: for space-time convergence ($c = 4$)</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 18** Solutions for the 3D problem. (a) test case 1, and (c) test case 2: time loading; (b) test case 1, and (d) test case 2: evolution of the fraction of cracked surfaces (c = 1); (e) test case 2: for space-time convergence (c = 3); (f) test case 2: for space-time convergence (c = 4).</figcaption></figure>

If $\lambda = \infty$ (no influence of tangential component), no failure occurs for this test case, since the interfaces are subjected to a normal compression. On the other hand, with the selected value of $\lambda$ of Table 1, shearing is taken into account and the evolution of the fraction of cracked surface of the interfaces are depicted in Figures 18(b,d), for the two loads. Indeed, the smooth loading leads to a delayed fracture, since the loading increases slower than for the first case. Moreover, the fracture curves are smoother as well, since there are less forward and backward elastic waves traveling in the structure. In each case, we can observe that the fracture, induced mainly herein by shearing of the interface, is first developing in the interface number 2, Figure 17, more subjected to shearing.

To question the space-time convergence, the same test is performed again, with a mesh size divided by a factor c (i.e. for $c^3$ times more finite elements) and a time step divided by c also (the ratio of the time step over the critical ones is therefore constant). The corresponding evolutions of the fraction of cracked surfaces are depicted in Figures 18(d-f) for $c \in \{1, 3, 4\}$. The qualitative comparison is consistent; one can notice a kind of smoothing of the evolution curves for these macroscopic quantities, that may be due to the higher resolution.

Since we got a solution which is a space-time field, the norm on the residual penetration to get an error indicator is questionable. Indeed, we can use an spatial average and a maximum
backward elastic waves traveling in the structure. In each case, we can observe that the fracture, induced mainly herein by shearing of the interface, is first developing in the interface number 2, Figure 17, more subjected to shearing.

To question the space-time convergence, the same test is performed again, with a mesh size divided by a factor c (i.e. for c³ times more finite elements) and a time step divided by c also (the ratio of the time step over the critical ones is therefore constant). The corresponding evolutions of the fraction of cracked surfaces are depicted on Figures 18(d-f) for c ∈ {1, 3, 4}. The qualitative comparison is consistent; one can notice a kind of smoothing of the evolution curves for these macroscopic quantities, that may be due to the higher resolution.

Since we got a solution which is a space-time field, the norm on the residual penetration to get an error indicator is questionable. Indeed, we can use an spatial average and a maximum norm in time, or a full averaged one (here using the mean value), for instance:

$$ \eta = \frac{\max_t(-\text{mean}_x g(x, t))}{\max_t \text{mean}_x |g(x, t)|} \quad \text{or} \quad \eta' = \frac{-\text{mean}_{x,t} < g(x, t) >_-}{\text{mean}_{x,t} |g(x, t)|} \quad (15) $$

where <•> stands for the negative part. For the different refinement ratios c, simultaneously for the mesh size and the time steps, the obtained values are reported in Table 2. These errors are quite small, and the tendency is a decrease with the refinement factor c.

<table><caption>**Table 2:** Error indicators vs the refinement ratio</caption><thead><tr><th>c</th><th>1</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>η</td><td>0.016 %</td><td>0.012 %</td><td>0.009 %</td></tr><tr><td>η'</td><td>0.036 %</td><td>0.023 %</td><td>0.019 %</td></tr></tbody></table>

## 4 Conclusions and perspectives

The previously developed explicit dynamic scheme CD-Lagrange has been tested on different highly non-smooth behaviors localized in mechanical interfaces. It provides a framework for the implementation of non-smooth (and not regularized) relationships, provided that the constitutive behaviors are expressed between impulse and velocity. Perfect plasticity with rigid interface and perfectly brittle fracture have been tested for a single degree-of-freedom problem, and validated on a full 3D finite element peeling-like problem.

Concerning perspectives, testing if the non-smooth (and non-univoque) constitutive relations may lead to the so-called deterministic chaos, i.e. an a priori deterministic problem, but whose solution may strongly depend on initial conditions, is under concern. Some less non-smooth cases can also be studied, as ductile damage (cohesive zone models), for which an extrinsic case is feasible with no initial compliance (Collins-Craft et al. 2022). To avoid spatial localization a delay-effect model (Needleman 1988; Allix et al. 1997; Allix et al. 2003; Desmorat et al. 2010) will be more suited since non-locality in space (Pijaudier-Cabot et al. 1987; Lasry et al. 1988; Kamensky et al. 2022) would prevent an explicit matrix-free approach. Generalization to asynchronous explicit schemes (Gravouil et al. 2014; Niu et al. 2022) is also under concern. Comparing the present test cases, in particular the 3D one, with other time integration schemes is also a direct perspective to this work. Indeed, the proposed test cases could be used as benchmark cases, as additional ones to those already proposed in (Di Stasio et al. 2019), to test the implementation on other codes or with other time integration schemes, as an other direct perspective.

**Funding** We gratefully acknowledge the French National Association for Research and Technology (ANRT, CIFRE grant number 2021/0957). This work was supported by the "Manufacture Française de Pneumatiques Michelin".

## Bibliographic references

Acary, V. (2012). "Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts". *Applied Numerical Mathematics* 62.10, pp. 1259–1275. DOI: 10.1016/j.apnum.2012.06.026. URL: https://hal.science/inria-00476398/
norm in time, or a full averaged one (here using the mean value), for instance

$$ \eta = \frac{\max_t(-\text{mean}_x g(x, t))}{\max_t \text{mean}_x |g(x, t)|} \quad \text{or} \quad \eta' = \frac{-\text{mean}_{x,t} \langle g(x, t) \rangle_-}{\text{mean}_{x,t} |g(x, t)|} \qquad (15) $$

where $\langle \cdot \rangle_-$ stands for the negative part. For the different refinement ratios c, simultaneously for the mesh size and the time steps, the obtained values are reported in Table 2. These errors are quite small, and the tendency is a decrease with the refinement factor c.

<table><caption>**Table 2** Error indicators versus refinement ratio.</caption><thead><tr><th>c</th><th>1</th><th>3</th><th>4</th></tr></thead><tbody><tr><th>η</th><td>0.016 %</td><td>0.012 %</td><td>0.009 %</td></tr><tr><th>η'</th><td>0.036 %</td><td>0.023 %</td><td>0.019 %</td></tr></tbody></table>

# 4 Conclusions and perspectives

The previously developed explicit dynamic scheme CD-Lagrange has been tested on different highly non-smooth behaviors localized in mechanical interfaces. It provides a framework for the implementation of non-smooth (and not regularized) relationships, provided that the constitutive behaviors are expressed between impulse and velocity. Perfect plasticity with rigid interface and perfectly brittle fracture have been tested for a single degree-of-freedom problem, and validated on a full 3D finite element peeling-like problem.

Concerning perspectives, testing if the non-smooth (and non single-valued) constitutive relations may lead to the so-called deterministic chaos, i.e. an a priori deterministic problem, but whose solution may strongly depend on initial conditions, is under concern. Some less non-smooth cases can also be studied, as ductile damage (cohesive zone models), for which an extrinsic case is feasible with no initial compliance (Collins-Craft et al. 2022). To avoid spatial localization a delay-effect model (Needleman 1988; Allix et al. 1997; Allix et al. 2003; Desmorat et al. 2010) will be more suited since non-locality in space (Pijaudier-Cabot et al. 1987; Lasry et al. 1988; Kamensky et al. 2022) would prevent an explicit matrix-free approach. Generalization to asynchronous explicit schemes (Gravouil et al. 2014; Niu et al. 2022) is also under concern. Comparing the present test cases, in particular the 3D one, with other time integration schemes is also a direct perspective to this work. Indeed, the proposed test cases could be used as benchmark cases, as additional ones to those already proposed in (Di Stasio et al. 2019), to test the implementation on other codes or with other time integration schemes, as an other direct perspective.

# References

* Acary, V. (2012). Higher order event capturing time-stepping schemes for nonsmooth multi-body systems with unilateral constraints and impacts. *Applied Numerical Mathematics* 62(10):1259-1275. [DOI], [HAL].

* Acary, V., M. Brémond, and O. Huber (2018). On solving contact problems with Coulomb friction: Formulations and numerical comparisons. *Advanced Topics in Nonsmooth Dynamics*. Springer, 375-457. [DOI], [HAL].

* Alart, P. (2014). How to overcome indetermination and interpenetration in granular systems via nonsmooth contact dynamics. An exploratory investigation. *Computer Methods in Applied Mechanics and Engineering* 270:37-56. [DOI], [HAL].

* Allix, O., P. Feissel, and P. Thévenet (2003). A delay damage mesomodel of laminates under dynamic loading: basic aspects and identification issues. *Computers & Structures* 81(12):1177-1191. [DOI], [HAL].

* Allix, O. and J.-F. Deü (1997). Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading. *Engineering Transactions* 45(1):29-46. [DOI], [OA].

* Belgacem, F., P. Hild, and P. Laborde (1998). The mortar finite element method for contact problems. *Mathematical and Computer Modelling* 28(4-8):263-271. [DOI], [OA].
Acary, V., M. Brémond, and O. Huber (2018). "On solving contact problems with Coulomb friction: Formulations and numerical comparisons". *Advanced Topics in Nonsmooth Dynamics*. Springer International Publishing, pp. 375-457. DOI: 10.1007/978-3-319-75972-2_10. URL: https://hal.science/hal-01878539

Alart, P. (2014). "How to overcome indetermination and interpenetration in granular systems via nonsmooth contact dynamics. An exploratory investigation". *Computer Methods in Applied Mechanics and Engineering* 270, pp. 37-56. DOI: 10.1016/j.cma.2013.11.020. URL: https://hal.science/hal-02075883

Allix, O. and J.-F. Deü (1997). "Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading". *Engineering Transactions* 45.1, pp. 29-46. DOI: 10.24423/engtrans.680.1997. URL: https://et.ippt.gov.pl/index.php/et/article/view/680/501

Allix, O., P. Feissel, and P. Thévenet (2003). "A delay damage mesomodel of laminates under dynamic loading: basic aspects and identification issues". *Computers & Structures* 81.12, pp. 1177-1191. DOI: 10.1016/S0045-7949(03)00035-X. URL: https://hal.science/hal-01745211

Belytschko, T., W. K. Liu, B. Moran, and K. Elkhodary (2014). *Nonlinear finite elements for continua and structures*. Wiley

Ben Belgacem, F., P. Hild, and P. Laborde (1998). "The mortar finite element method for contact problems". *Mathematical and Computer Modelling* 28.4-8, pp. 263-271. DOI: 10.1016/S0895-7177(98)00121-6

Brisotto, D. d. S., E. Bittencourt, and V. M. R. D. Bessa (2018). "Simulation of interface behavior between FRP bars and concrete by an elastic-plastic theory via FEM". *Revista IBRACON de Estruturas e Materiais* 11.6, pp. 1381-1390. DOI: 10.1590/S1983-41952018000600011

Brogliato, B., T. ten Dam, L. Paoli, F. Génot, and M. Abadie (2002). "Numerical simulation of finite dimensional multibody nonsmooth mechanical systems". *Applied Mechanics Reviews* 55.2, pp. 107-150. DOI: 10.1115/1.1454112. URL: https://hal.science/hal-01349852/

Casadei, F., E. Gabellini, G. Fotia, F. Maggio, and A. Quarteroni (2002). "A mortar spectral/finite element method for complex 2D and 3D elastodynamic problems". *Computer Methods in Applied Mechanics and Engineering* 191.45, pp. 5119-5148. DOI: 10.1016/S0045-7825(02)00294-3

Chen, Q. Z., V. Acary, G. Virlez, and O. Brüls (2013). "A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints". *International Journal for Numerical Methods in Engineering* 96.8, pp. 487-511. DOI: 10.1002/nme.4563. URL: https://hal.science/hal-00858066

Chouly, F. and Y. Renard (2018). "Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems". *Advanced Modeling and Simulation in Engineering Science* 5.1, pp. 31-68. DOI: 10.1186/S40323-018-0124-5. URL: https://hal.science/hal-01814774

Collins-Craft, N. A., F. Bourrier, and V. Acary (2022). "On the formulation and implementation of extrinsic cohesive zone models with contact". *Computer Methods in Applied Mechanics and Engineering* 400, p. 115545. DOI: 10.1016/j.cma.2022.115545. URL: https://hal.science/hal-03371667

Dabaghi, F., P. Krejčí, A. Petrov, J. Pousin, and Y. Renard (2019). "A weighted finite element mass redistribution method for dynamic contact problems". *Journal of Computational and Applied Mathematics* 345, pp. 338-356. DOI: 10.1016/j.cam.2018.06.030. URL: https://hal.science/hal-01237748

de los Ríos, G. S. A., R. C. Balderas, V. A. Duong, S. Chataigner, J.-F. Caron, A. Ehrlacher, G. Foret, and A. D. Diaz (2009). "Laminated materials with plastic interfaces: modeling and calculation". *Modelling and Simulation in Materials Science and Engineering* 17.2, p. 025008. DOI: 10.1088/0965-0393/17/2/025008

de Saxcé, G. (2022). "A non incremental variational principle for brittle fracture". *International Journal of Solids and Structures* 252, p. 111761. DOI: 10.1016/j.ijsolstr.2022.111761. URL: https://arxiv.org/abs/2112.03065

de Saxcé, G. and Z.-Q. Feng (1998). "The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms". *Mathematical and Computer Modelling* 28.4-8, pp. 225-245. DOI: 10.1016/S0895-7177(98)00119-8

Desmorat, R., M. Chambart, F. Gatuingt, and D. Guilbaud (2010). "Delay-active damage versus non-local enhancement for anisotropic damage dynamics computations with alternated loading".
Belytschko, T., W. K. Liu, B. Moran, and K. Elkhodary (2014). *Nonlinear Finite Elements for Continua and Structures*. 2nd ed. Wiley. ISBN: 9781118632703.

Brisotto, D. S., E. Bittencourt, and V. M. R. D. Bessa (2018). Simulation of interface behavior between FRP bars and concrete by an elastic-plastic theory via FEM. *Revista IBRACON de Estruturas e Materiais* 11(6):1381-1390. [DOI], [OA].

Brogliato, B., A. ten Dam, L. Paoli, F. Génot, and M. Abadie (2002). Numerical simulation of finite dimensional multibody nonsmooth mechanical systems. *Applied Mechanics Reviews* 55(2):107-150. [DOI], [HAL].

Carvalho, R. P., A. M. C. Carneiro, F. M. A. Pires, and A. Popp (2022). An efficient algorithm for rigid/deformable contact interaction based on the dual mortar method. *Computational Mechanics* 71(1):143-167. [DOI], [OA].

Casadei, F., E. Gabellini, G. Fotia, F. Maggio, and A. Quarteroni (2002). A mortar spectral/finite element method for complex 2D and 3D elastodynamic problems. *Computer Methods in Applied Mechanics and Engineering* 191(45):5119-5148. [DOI].

Chen, Q., V. Acary, G. Virlez, and O. Brüls (2013). A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints. *International Journal for Numerical Methods in Engineering* 96(8):487-511. [DOI], [HAL].

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

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

Dabaghi, F., P. Krejčí, A. Petrov, J. Pousin, and Y. Renard (2019). A weighted finite element mass redistribution method for dynamic contact problems. *Journal of Computational and Applied Mathematics* 345:338-356. [DOI], [OA].

De Saxcé, G. and Z.-Q. Feng (1998). The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms. *Mathematical and Computer Modelling* 28(4-8):225-245. [DOI], [OA].

De Saxcé, G. (2022). A non incremental variational principle for brittle fracture. *International Journal of Solids and Structures* 252:111761. [DOI], [OA].

Desmorat, R., M. Chambart, F. Gatuingt, and D. Guilbaud (2010). Delay-active damage versus non-local enhancement for anisotropic damage dynamics computations with alternated loading. *Engineering Fracture Mechanics* 77(12):2294-2315. [DOI], [HAL].

Di Stasio, J., D. Dureisseix, G. Georges, A. Gravouil, and T. Homolle (2021). An explicit time-integrator with singular mass for non-smooth dynamics. *Computational Mechanics* 68(1):97-112. [DOI], [HAL].

Di Stasio, J., D. Dureisseix, A. Gravouil, G. Georges, and T. Homolle (2019). Benchmark cases for robust explicit time integrators in non-smooth transient dynamics. *Advanced Modeling and Simulation in Engineering Sciences* 6(1). [DOI], [OA].

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

Fekak, F.-E., M. Brun, A. Gravouil, and B. Depale (2017). A new heterogeneous asynchronous explicit-implicit time integrator for nonsmooth dynamics. *Computational Mechanics* 60(1):1-21. [DOI], [HAL].

Feng, Z.-Q., P. Joli, J.-M. Cros, and B. Magnain (2005). The bi-potential method applied to the modeling of dynamic problems with friction. *Computational Mechanics* 36(5):375-383. [DOI], [HAL].

Francfort, G. (2006). Quasistatic brittle fracture seen as an energy minimizing movement. *GAMM-Mitteilungen* 29(2):172-191. [DOI].

Ghaednia, H., X. Wang, S. Saha, Y. Xu, A. Sharma, and R. L. Jackson (2017). A review of elastic-plastic contact mechanics. *Applied Mechanics Reviews* 69(6). [DOI], [OA].

Gravouil, A., A. Combescure, and M. Brun (2014). Heterogeneous asynchronous time integrators for computational structural dynamics. *International Journal for Numerical Methods in Engineering* 102(3-4):202-232. [DOI], [HAL].
*Engineering Fracture Mechanics* 77.12, pp. 2294-2315. DOI: 10.1016/j.engfracmech.2010.04.006.
URL: https://hal.science/hal-00994292

Di Stasio, J., D. Dureisseix, G. Georges, A. Gravouil, and T. Homolle (2021). "An explicit time-integrator with singular mass for non-smooth dynamics". *Computational Mechanics* 68.1. Ed. by S. Verlag, pp. 97-112. DOI: 10.1007/s00466-021-02021-5. URL: https://hal.science/hal-03212569/

Di Stasio, J., D. Dureisseix, A. Gravouil, G. Georges, and T. Homolle (2019). "Benchmark cases for robust explicit time integrators in non-smooth transient dynamics". *Advanced Modeling and Simulation in Engineering Sciences* 6.1. DOI: 10.1186/s40323-019-0126-y. URL: https://amses-journal.springeropen.com/counter/pdf/10.1186/s40323-019-0126-y

Dubois, F., V. Acary, and M. Jean (2018). "The Contact Dynamics method: A nonsmooth story". *Comptes Rendus Mécanique* 346.3, pp. 247-262. DOI: 10.1016/j.crme.2017.12.009. URL: https://hal.science/hal-01676287

Fekak, F.-E., M. Brun, A. Gravouil, and B. Depale (2017). "A new heterogeneous asynchronous explicit-implicit time integrator for nonsmooth dynamics". *Computational Mechanics* 60.1, pp. 1-21. DOI: 10.1007/s00466-017-1397-0. URL: https://hal.science/hal-01852099

Feng, Z.-Q., P. Joli, J.-M. Cros, and B. Magnain (2005). "The bi-potential method applied to the modeling of dynamic problems with friction". *Computational Mechanics* 36.5, pp. 375-383. DOI: 10.1007/s00466-005-0663-8. URL: https://hal.science/hal-00342921

Francfort, G. A. (2006). "Quasistatic brittle fracture seen as an energy minimizing movement". *GAMM-Mitteilungen* 29.2, pp. 172-191. DOI: 10.1002/gamm.201490029

Ghaednia, H., X. Wang, S. Saha, Y. Xu, A. Sharma, and R. L. Jackson (2017). "A review of elastic-plastic contact mechanics". *Applied Mechanics Reviews* 69.6. DOI: 10.1115/1.4038187

Gravouil, A., A. Combescure, and M. Brun (2014). "Heterogeneous asynchronous time integrators for computational structural dynamics". *International Journal for Numerical Methods in Engineering* 102.3-4, pp. 202-232. DOI: 10.1002/nme.4818

Hairer, E., C. Lubich, and G. Wanner (2003). "Geometric numerical integration illustrated by the Störmer-Verlet method". *Acta Numerica* 12, pp. 399-450. DOI: 10.1017/s0962492902000144

Jean, M., V. Acary, and Y. Monerie (2001). "Non-smooth contact dynamics approach of cohesive materials". *Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences* 359.1789. Ed. by F. G. Pfeiffer, pp. 2497-2518. DOI: 10.1098/rsta.2001.0906. URL: https://hal.science/inria-00423556

Kamensky, D., M. D. Alaydin, and Y. Bazilevs (2022). "A Review of Nonlocality in Computational Contact Mechanics". *Current Trends and Open Problems in Computational Mechanics*. Springer International Publishing, pp. 239-246. DOI: 10.1007/978-3-030-87312-7_23

Khenous, H. B., P. Laborde, and Y. Renard (2008). "Mass redistribution method for finite element contact problems in elastodynamics". *European Journal of Mechanics, A/Solids* 27.5, pp. 918-932. DOI: 10.1016/j.euromechsol.2008.01.001. URL: https://hal.science/hal-00582045

Kiener, D. and S. M. Han (2022). "100 years after Griffith: From brittle bulk fracture to failure in 2D materials". *MRS Bulletin* 47.8, pp. 792-799. DOI: 10.1557/s43577-022-00379-2

Kogut, L. and I. Etsion (2003). "A finite element based elastic-plastic model for the contact of rough surfaces". *Tribology Transactions* 46.3, pp. 383-390. DOI: 10.1080/10402000308982641

Lasry, D. and T. Belytschko (1988). "Localization limiters in transient problems". *International Journal of Solids and Structures* 24.6, pp. 581-597. DOI: 10.1016/0020-7683(88)90059-5

Lebon, F. and R. Rizzoni (2022). "On the emergence of adhesion in asymptotic analysis of piecewise linear anisotropic elastic bonded joints". *European Journal of Mechanics - A/Solids* 93, p. 104512. DOI: 10.1016/j.euromechsol.2022.104512. URL: https://hal.science/hal-03546068

Lemaitre, J. and J.-L. Chaboche (1994). *Mechanics of Solid Materials*. Cambridge University Pres

Liu, G., J. Zhu, L. Yu, and Q. J. Wang (2001). "Elasto-plastic Contact of Rough Surfaces". *Tribology Transactions* 44.3, pp. 437-443. DOI: 10.1080/10402000108982478

Meyer, C. and S. Walther (2022). "Optimal control of perfect plasticity. Part I: stress tracking". *Mathematical control and related fields* 12.2, pp. 275-301. DOI: 10.3934/mcrf.2021022. URL: https://arxiv.org/abs/2001.02969

Monchiet, V. and G. Bonnet (2010). "Interfacial models in viscoplastic composites materials". *International Journal of Engineering Science* 48.12, pp. 1762-1768. DOI: 10.1016/j.ijengsci.2010.09.024. URL: https://hal.science/hal-00687820
Hairer, E., C. Lubich, and G. Wanner (2003). Geometric numerical integration illustrated by the Störmer–Verlet method. *Acta Numerica* 12:399–450. [DOI], [OA].

Jean, M., V. Acary, and Y. Monerie (2001). Non-smooth contact dynamics approach of cohesive materials. *Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences* 359(1789):2497–2518. [DOI], [HAL].

Kamensky, D., M. D. Alaydin, and Y. Bazilevs (2022). A Review of Nonlocality in Computational Contact Mechanics. *Current Trends and Open Problems in Computational Mechanics*. Springer, 239–246. [DOI].

Khenous, H. B., P. Laborde, and Y. Renard (2008). Mass redistribution method for finite element contact problems in elastodynamics. *European Journal of Mechanics - A/Solids* 27(5):918–932. [DOI], [HAL].

Kiener, D. and S. M. Han (2022). 100 years after Griffith: From brittle bulk fracture to failure in 2D materials. *MRS Bulletin* 47(8):792–799. [DOI], [OA].

Kogut, L. and I. Etsion (2003). A finite element based elastic-plastic model for the contact of rough surfaces. *Tribology Transactions* 46(3):383–390. [DOI].

Lasry, D. and T. Belytschko (1988). Localization limiters in transient problems. *International Journal of Solids and Structures* 24(6):581–597. [DOI], [HAL].

Lebon, F. and R. Rizzoni (2022). On the emergence of adhesion in asymptotic analysis of piecewise linear anisotropic elastic bonded joints. *European Journal of Mechanics - A/Solids* 93:104512. [DOI], [HAL].

Lemaitre, J. and J.-L. Chaboche (1994). *Mechanics of Solid Materials*. Cambridge University Press. [DOI].

Liu, G., J. Zhu, L. Yu, and Q. J. Wang (2001). Elasto-plastic Contact of Rough Surfaces. *Tribology Transactions* 44(3):437–443. [DOI].

Meyer, C. and S. Walther (2022). Optimal control of perfect plasticity. Part I: stress tracking. *Mathematical Control & Related Fields* 12(2):275. [DOI], [OA].

Monchiet, V. and G. Bonnet (2010). Interfacial models in viscoplastic composites materials. *International Journal of Engineering Science* 48(12):1762–1768. [DOI], [HAL].

Moreau, J. J. (1988). Unilateral Contact and Dry Friction in Finite Freedom Dynamics. *Nonsmooth Mechanics and Applications*. Vol. 302. CISM International Centre for Mechanical Sciences. Springer. Chap. 1, 1–82. [DOI], [HAL].

Moreau, J. J. (1999). Numerical aspects of the sweeping process. *Computer Methods in Applied Mechanics and Engineering* 177(3–4):329–349. [DOI], [HAL].

Moreau, J. J. (1978). Approximation en graphe d'une évolution discontinue. *RAIRO. Analyse numérique* 12(1):75–84. [DOI], [OA].

Moreau, J. J. (1986). Une formulation du contact à frottement sec ; application au calcul numérique. *Comptes Rendus de l'Académie des Sciences Série II* 302(13):799–801. [HAL].

Moreau, J. J. (2004). Indetermination due to dry friction in multibody dynamics. *European Congress on Computational Methods in Applied Sciences and Engineering* (Jyväskylä, Finland, July 24–24, 2004). [HAL].

Moreau, J. J. (2006). Facing the plurality of solutions in nonsmooth mechanics. *Second International Conference on Nonsmooth/Nonconvex Mechanics with Applications in Engineering* (Thessaloniki, Greece, July 7–8, 2006), pp 3–12. [HAL].

Moreau, J. J. (2011). On Unilateral Constraints, Friction and Plasticity. *New Variational Techniques in Mathematical Physics*. Springer, 171–322. [DOI], [HAL].

Needleman, A. (1988). Material rate dependence and mesh sensitivity in localization problems. *Computer Methods in Applied Mechanics and Engineering* 67(1):69–85. [DOI].

Niu, Z., V. Ziaei-Rad, Z. Wu, and Y. Shen (2022). An asynchronous variational integrator for the phase field approach to dynamic fracture. *International Journal for Numerical Methods in Engineering* 124(2):434–457. [DOI], [OA].

Park, K. and G. H. Paulino (2011). Cohesive zone models: A critical review of traction-separation relationships across fracture surfaces. *Applied Mechanics Reviews* 64(6). [DOI], [HAL].

Pijaudier-Cabot, G. and Z. P. Bažant (1987). Nonlocal damage theory. *Journal of Engineering Mechanics* 113(10):1512–1533. [DOI], [HAL].

Popp, A. and W. A. Wall (2014). Dual mortar methods for computational contact mechanics –
Moreau, J. J. (1978). "Approximation en graphe d'une évolution discontinue". *RAIRO. Analyse numérique* 12.1, pp. 75-84. URL: https://hal.science/hal-01788586

Moreau, J. J. (1986). "Une formulation du contact à frottement sec ; application au calcul numérique". *Comptes Rendus de l'Académie des Sciences Serie II* 302.13, pp. 799-801. URL: https://hal.science/hal-01349850

Moreau, J. J. (1988). "Unilateral Contact and Dry Friction in Finite Freedom Dynamics". *Nonsmooth Mechanics and Applications*. Vol. 302. CISM International Centre for Mechanical Sciences. Springer Vienna. Chap. 1, pp. 1-82. DOI: 10.1007/978-3-7091-2624-0_1. URL: https://hal.science/hal-01713847

Moreau, J. J. (1999). "Numerical aspects of the sweeping process". *Computer Methods in Applied Mechanics and Engineering* 177.3-4, pp. 329-349. DOI: 10.1016/S0045-7825(98)00387-9. URL: https://hal.science/hal-01349847

Moreau, J. J. (2004). "Indetermination due to dry friction in multibody dynamics". *European Congress on Computational Methods in Applied Sciences and Engineering*. ECCOMAS 2004 proceedings. Jyväskylä, Finland. URL: https://hal.science/hal-01825251

Moreau, J. J. (2006). "Facing the plurality of solutions in nonsmooth mechanics". *Nonsmooth/Nonconvex Mechanics with Applications in Engineering*. Ed. by C. C. Baniotopoulos. Thessalonique, Greece, pp. 3-12. URL: https://hal.science/hal-01824588

Moreau, J. J. (2011). "On Unilateral Constraints, Friction and Plasticity". *New Variational Techniques in Mathematical Physics*. Springer Berlin Heidelberg, pp. 171-322. DOI: 10.1007/978-3-642-10960-7_7

Needleman, A. (1988). "Material rate dependence and mesh sensitivity in localization problems". *Computer Methods in Applied Mechanics and Engineering* 67.1, pp. 69-85. DOI: 10.1016/0045-7825(88)90069-2

Niu, Z., V. Ziaei-Rad, Z. Wu, and Y. Shen (2022). "An asynchronous variational integrator for the phase field approach to dynamic fracture". *International Journal for Numerical Methods in Engineering* 124.2, pp. 434-457. DOI: 10.1002/nme.7127. URL: https://arxiv.org/abs/2203.08797

Park, K. and G. H. Paulino (2011). "Cohesive zone models: A critical review of traction-separation relationships across fracture surfaces". *Applied Mechanics Reviews* 64.6. DOI: 10.1115/1.4023110

Pijaudier-Cabot, G. and Z. P. Bažant (1987). "Nonlocal damage theory". *Journal of Engineering Mechanics* 113.10, pp. 1512-1533. DOI: 10.1061/(ASCE)0733-9399(1987)113:10(1512)

Pinto Carvalho, R., A. M. Couto Carneiro, F. M. Andrade Pires, and A. Popp (2022). "An efficient algorithm for rigid/deformable contact interaction based on the dual mortar method". *Computational Mechanics* 71.1, pp. 143-167. DOI: 10.1007/s00466-022-02226-2. URL: https://arxiv.org/abs/2201.01165

Popp, A. and W. A. Wall (2014). "Dual mortar methods for computational contact mechanics – overview and recent developments". *GAMM-Mitteilungen* 37.1, pp. 66-84. DOI: 10.1002/gamm.201410004

Raffa, M. L., F. Lebon, and R. Rizzoni (2022). "A micromechanical model of a hard interface with micro-cracking damage". *International Journal of Mechanical Sciences* 216, p. 106974. DOI: 10.1016/j.ijmecsci.2021.106974. URL: https://hal.science/hal-03404654

Raous, M., M. Schryve, and M. Cocou (2006). "Restorable adhesion and friction". *Nonsmooth/Nonconvex Mechanics with applications in Engineering*. ZITI Publisher, pp. 165-172. URL: https://hal.science/hal-00476784

Spada, A., G. Giambanco, and P. Rizzo (2009). "Damage and plasticity at the interfaces in composite materials and structures". *Computer Methods in Applied Mechanics and Engineering* 198.49-52, pp. 3884-3901. DOI: 10.1016/j.cma.2009.08.024

Tkachuk, A., B. I. Wohlmuth, and M. Bischoff (2013). "Hybrid-mixed discretization of elasto-dynamic contact problems using consistent singular mass matrices". *International Journal for Numerical Methods in Engineering* 94.5, pp. 473-493. DOI: 10.1002/nme.4457

Wohlmuth, B. I. (2000). "A mortar finite element method using dual spaces for the Lagrange multiplier". *SIAM Journal on Numerical Analysis* 38.3, pp. 989-1012. DOI: 10.1137/S0036142999350929. URL: https://hal.science/hal-03747534
overview and recent developments. *GAMM-Mitteilungen* 37(1):66–84. [DOI].

Raffa, M. L., F. Lebon, and R. Rizzoni (2022). A micromechanical model of a hard interface with micro-cracking damage. *International Journal of Mechanical Sciences* 216:106974. [DOI], [HAL].

Raous, M., M. Schryve, and M. Cocou (2006). Recoverable adhesion and friction. *Second International Conference on Nonsmooth/Nonconvex Mechanics with Applications in Engineering* (Thessaloniki, Greece, July 7–8, 2006), pp 165–172. [HAL].

Ríos, G. S. Aquino de los, R. C. Balderas, V. A. Duong, S. Chataigner, J.-F. Caron, A. Ehrlacher, G. Foret, and A. Diaz Diaz (2009). Laminated materials with plastic interfaces: modeling and calculation. *Modelling and Simulation in Materials Science and Engineering* 17(2):025008. [DOI], [HAL].

Spada, A., G. Giambanco, and P. Rizzo (2009). Damage and plasticity at the interfaces in composite materials and structures. *Computer Methods in Applied Mechanics and Engineering* 198(49–52):3884–3901. [DOI].

Tkachuk, A., B. I. Wohlmuth, and M. Bischoff (2013). Hybrid-mixed discretization of elasto-dynamic contact problems using consistent singular mass matrices. *International Journal for Numerical Methods in Engineering* 94(5):473–493. [DOI], [HAL].

Wohlmuth, B. I. (2000). A mortar finite element method using dual spaces for the Lagrange multiplier. *SIAM Journal on Numerical Analysis* 38(3):989–1012. [DOI], [HAL].

Wriggers, P. and G. Zavarise (1993). Application of augmented Lagrangian techniques for non-linear constitutive laws in contact interfaces. *Communications in Numerical Methods in Engineering* 9(10):815–824. [DOI], [HAL].

Wriggers, P. and T. A. Laursen (2007). *Computational Contact Mechanics*. Springer. [DOI].

Zhang, X., Z. Qi, G. Wang, and S. Guo (2019). Model smoothing method of contact-impact dynamics in flexible multibody systems. *Mechanism and Machine Theory* 138:124–148. [DOI].

Zhang, Z. J., G. H. Paulino, and W. Celes (2007). Extrinsic cohesive modelling of dynamic fracture and microbranching instability in brittle materials. *International Journal for Numerical Methods in Engineering* 72(8):893–923. [DOI].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** DD carried out most of the study, performed numerical simulations, and drafted the manuscript. PL helped with implementation and numerical issues. All authors developed the methodology, conceived of the study, and participated in its design, coordination, and critical review of the manuscript. All authors read and approved the final manuscript.

**Supplementary Material** The 3D results exposed in the present paper can be reproduced using the datasets and software solutions archived at the [permalink](https://permalink.10.5281/zenodo.13778789).

**Funding** We gratefully acknowledge the French National Association for Research and Technology (ANRT, CIFRE grant number 2021/0957). This work was supported by the "Manufacture Française de Pneumatiques Michelin".

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.
Wriggers, P. and G. Zavarise (1993). "Application of augmented Lagrangian techniques for non-linear constitutive laws in contact interfaces". *Communications in Numerical Methods in Engineering* 9.10, pp. 815-824. DOI: 10.1002/cnm.1640091005

Wriggers, P. and T. A. Laursen (2007). *Computational Contact Mechanics*. Springer Vienna. DOI: 10.1007/978-3-211-77298-0

Zhang, X., Z. Qi, G. Wang, and S. Guo (2019). "Model smoothing method of contact-impact dynamics in flexible multibody systems". *Mechanism and Machine Theory* 138, pp. 124-148. DOI: 10.1016/j.mechmachtheory.2019.03.039

Zhang, Z. J., G. H. Paulino, and W. Celes (2007). "Extrinsic cohesive modelling of dynamic fracture and microbranching instability in brittle materials". *International Journal for Numerical Methods in Engineering* 72.8, pp. 893-923. DOI: 10.1002/nme.2030