M Brun

R Rezakhani
email: roozbeh.rezakhani@epfl.ch

J-F Molinari

Explicit dynamic approach for unbounded domains in frictional contact with Rate and State laws

Keywords: Finite Element Method, Hybrid Asynchronous Time Integrator, frictional contact, Rate and State friction law, unbounded domains, Perfectly Matched Layers

The paper explores the ability of an explicit time integration procedure to simulate the dynamics of shear rupture between unbounded elastic blocks on frictional interface, modeled with the finite element method. The behavior of the interface is governed by Rate and State (RS) friction laws, proposed to describe the rate dependent phenomena observed in experiments on rocks and many other materials in sliding contact. The method for the frictional contact between bodies is integrated in the framework of heterogeneous asynchronous time integrator (HATI). Dual formulation, requiring the introduction of Lagrange multipliers, is adopted by dealing in a weak way with the normal contact conditions as well as the tangential frictional conditions expressed in velocity. Taking advantage of the flexibility of the HATI framework, precise formulations of Perfectly Matched Layers (PML) are also incorporated, enabling us to treat semi-infinite media in frictional contact. Simulations concerning an unbounded elastic block on a rigid flat plane, initially compressed and sheared with remote loadings, are carried out. Advantages of the PML are underlined during the preliminary stages of loading as well as during the nucleation and propagation of the shear rupture along the interface. Numerical oscillations inherent from the high non-linearity of the RS friction law, characterized by significant variations of the friction coefficient with small changes in tangential velocity, are drastically reduced thanks to the HATI method.

Introduction

Frictional phenomena occur in various fields of mechanics at a broad range of scales, the engineering systems, such as gears and bearings, to the complex natural tectonic systems during earthquake processes. The mechanism of crustal earthquakes was extensively explored through experimental investigations for sliding rock blocks, concerning their stick-slip behavior and observed transient phenomena under abrupt changes in sliding velocity or normal stress [START_REF] Dieterich | Modeling of rock friction: 1. experimental results and constituve equations[END_REF][START_REF] Ruina | Slip instability and state variable friction law[END_REF][START_REF] Dieterich | Direct observation of frictional contacts: New insights for statedependent properties[END_REF]. The dependence of the observed friction coefficient on the sliding velocity was well characterized by slip rate and state (RS) constitutive formulations, such as the one proposed by Dieterich and Kilgore [START_REF] Dieterich | Direct observation of frictional contacts: New insights for statedependent properties[END_REF]: from an initial steady-state of sliding contact between two blocks, an abrupt change in sliding velocity results in an immediate jump in friction coefficient, followed by a displacement dependent decay and a stabilization at a new steady state of sliding contact. A state variable was also introduced in this formulation, representing a measure of the age of the continuously changing population of contacts between the sliding blocks, with a rate of change depending on the sliding velocity.

From a computational point of view, non-linear frictional constitutive laws, such as RS friction laws, were relatively recently employed in finite element simulations of shear rupture during earthquakes, but turn out to be very prone to numerical oscillations, in explicit integration dynamic framework as well as in implicit one. Needleman and co-authors [START_REF] Shi | Properties of dynamic rupture and energy partition in a solid with a frictional interface[END_REF][START_REF] Coker | Frictional sliding modes along an interface between identical elastic plates subject to shear impact loading[END_REF] developed an explicit Central Difference (CD) scheme to integrate the frictional law for two-dimensional sliding blocks under plane strain and small deformation assumption and noticed that the size of the time step should be kept very small: With respect to the classical CFL condition (Courant-Friedrichs-Lewy, [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]) for the stability of the explicit algorithm, small multiplicative coefficients were adopted, ranging from 10 -6 to 10 -2 . Recently, Tal and Hager [START_REF] Tal | Dynamic mortar finite element method for modeling of shear rupture on frictional rough surfaces[END_REF] employed an implicit time integration in a plane strain two-dimensional analysis in finite deformation for the integration of RS friction laws. They used a dual approach based on Lagrange multipliers for imposing the contact-normal conditions as well as the frictional-tangential conditions. Again, a high non-linearity was observed, leading to the increase of the number of iterations required for achieving the convergence at the end of the time step, as well as a strong decrease of the time step size during the propagation stage of the shear rupture. This finally results in a time step size close to the CFL condition, despite the adopted implicit time integration. In addition to the numerical oscillations due to the high non-linearity of the RS friction laws, wave reflections at the boundaries of the FE mesh can affect the behavior of the frictional interface. In order to ensure that no reflected elastic wave reaches the interface for a simulation time of interest, a simple solution is to consider a large dimension of the modeled solid in the direction normal to the interface [START_REF] Kammer | The existence of a critical length scale in regularised friction[END_REF], at the cost of increasing computation time. To avoid using an extended mesh, absorbing conditions have to be set up in order to damp out the main part of the incident waves, as done by DeDontney and co-authors [START_REF] Dedontney | Finite element modeling of branched ruptures including off-fault plasticity[END_REF], using the explicit dynamic FE package ABAQUS/Explicit for the investigation of branch activation in 2D plane strain model. Same strategy was employed by Ma and Beroza [START_REF] Ma | Rupture dynamics on a bimaterial interface for dipping faults[END_REF] for mimicking semi-infinite media so as to investigate dipping faults with bimaterial interface.

In this work, we propose to deal with unbounded domains with a RS frictional interface. For this purpose, the first main ingredient concerns the frictional contact, which is based on the explicit-implicit heterogeneous asynchronous time integrator (HATI) for nonsmooth transient dynamics with frictional contacts proposed by Fekak et al. [START_REF] Fekak | A new heterogeneous asynchronous explicitimplicit time integrator for nonsmooth dynamics[END_REF][START_REF] Di Stasio | Benchmark cases for robust explicit time integrators in non-smooth transient dynamics[END_REF][START_REF] Brun | Two FETI-based heterogeneous time step coupling methods for Newmark and α-schemes derived from the energy method[END_REF][START_REF] Gravouil | Heterogeneous asynchronous time integrators for computational structural dynamics[END_REF]. Thanks to the HATI formulation based on the introduction of Lagrange multipliers in the weak form of the equations of motion, explicit Central Difference (CD) scheme is employed for dealing with the region surrounding the sliding contact, whereas other parts of the problem can be handled with an implicit time integration scheme. In particular, in order to mimic unbounded domains, Hybrid Asynchronous Perfectly Matched Layers (HA-PML) at the boundaries of the truncated mesh were adopted and integrated using an implicit time integration scheme [START_REF] Brun | Hybrid asynchronous perfectly matched layer for seismic wave propagation in unbounded domains[END_REF], constituting the second main ingredient of our approach. Indeed, the PML is now considered as the most efficient way to absorb all types of incident waves, at all frequencies and incident angles [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF][START_REF] Basu | Erratum to: Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF][START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF][START_REF] Kucukcoban | Mixed Perfectly-Matched-Layers for direct transient analysis in 2D elastic heterogeneous media[END_REF][START_REF] Kucukcoban | A symmetric hybrid formulation for transient wave simulations in pml-truncated heterogeneous media[END_REF].

In order to assess the robustness of our approach, the test proposed by Kammer et al. [START_REF] Kammer | The existence of a critical length scale in regularised friction[END_REF], which concerns a 2D solid sliding on a rigid flat plane, is considered. By using the same nucleation procedure as in [START_REF] Kammer | The existence of a critical length scale in regularised friction[END_REF], that is by introducing an imposed fluid pressure varying progressively in time and space around a point belonging to the interface, it will be shown that a supershear single pulse in slip rate is created, followed with increasing numerical spurious vibrations, up to the creation of multiple peaks. Numerical oscillations were also observed in the literature, in particular in the case of 2D bilinear material simulation using Coulomb's law and slip weakening law, requiring the regularization of friction using, for example, the Prakash-Clifton friction law [START_REF] Prakash | Time resolved dynamic friction measurements in pressure-shear[END_REF] or artificial viscous damping such as Rayleigh viscous damping matrix [START_REF] Di Bartolomeo | Wave and rupture propagation at frictional bimaterial sliding interfaces: from local to global dynamics, from stick-slip to continuous sliding[END_REF][START_REF] Ma | Rupture dynamics on a bimaterial interface for dipping faults[END_REF]. Here, the HATI method is employed to mitigate these oscillations. In particular, the position of the HATI line, separating the contacting subdomain to the remaining of the medium, is located very close to the frictional interface to alleviate the numerical oscillations. It will be shown that the grid dependence of the results on the interface is significantly reduced, with a slip rate pulse converging with respect to the refinement in space and time.

The remainder of the paper is organized as follows: Section 2 introduces the problem with the strong form of the frictional contact between two bodies, the RS formulation and the PML in the time domain. Design of the PML is briefly reminded to achieve a target performance for its absorbing capabilities. The HATI framework is presented in section 3, by starting with the weak form of the equation of motion of unbounded domains: The velocity continuity between the PML region and the solid region as well as the contact and friction conditions are imposed in a weak way, by introducing Lagrange multipliers. Section 4 is devoted to the finite element discretization, in particular to the different terms in the weak form of the equation of motion related to the Lagrange multipliers, introduced at the PML interface and at the frictional contact interface. The normal contact and tangential frictional conditions are complemented with the RS friction law. 2D numerical applications are carried out in Section 5. First, the test proposed by Kammer et al. [START_REF] Kammer | The existence of a critical length scale in regularised friction[END_REF], concerning a 2D unbounded solid sliding on a rigid plane surface, with a constant friction coefficient is performed for different space and time discretisations, without regularisation. Then, the RS friction law is investigated and the numerical simulations show a supershear single pulse, propagating along the interface during a long time period of interest, whose intensity and shape are not altered by the space and time Right: Current configuration refinement, thanks to the HATI approach combined with the use of the PML.

Sliding contact between unbounded bodies

Strong form of the frictional contact between bodies

A two dimensional problem is considered under the small deformation assumption.

The two contacting bodies Ω 1 (slave body) and Ω 2 (master body) are illustrated in Figure 1. In Figure 1, Ω 1 and Ω 2 at time t are the current configurations of the two bodies with their interface Γ 1 and Γ 2 , assumed to be decomposed into three distinct parts, the Dirichlet (Γ u1 and Γ u2 ) and Neumann parts (Γ F1 and Γ F2 ) and the contact interface denoted by Γ C1 and Γ C2 .

In the normal direction, the contact condition has to be written for all the contact points at the current configuration:

         g N = [(X 1 + u 1 ) -(X 2 + u 2 )] • (-n 1 ) ≥ 0 τ N = σ i • n i • n i ≤ 0, i = 1, 2 g N τ N = 0 (1)
where X i , i = 1, 2, are the coordinates of nodes of the two bodies at the reference configuration and u i are the displacement fields at the current configuration; n 1 (resp. n 2 ) is the outer unit normal of the slave body (resp. master body) at contact points; σ i , i = 1, 2, are the stress tensors for the two bodies at contact points. The first inequality corresponds to the impenetrability condition with the definition of the gap function g N which can be either strictly positive when the bodies are separated and equal to zero when the bodies are in contact. The gap function is computed for all pairs of nodes potentially in contact, which are known a priori due to the small deformation assumption. The second inequality is called the intensility condition, expressing that the bodies press on each other with a compressive normal stress (τ N ≤ 0) or are inactive [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF][START_REF] Curnier | Unilateral contact[END_REF]. The third equation is called the complementarity equation, stating that the normal contact force does not work in the frictionless response. The three conditions correspond to Hertz-Signorini-Moreau conditions (HSM). The viability lemma proposed by Moreau [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF][START_REF] Moreau | Nonsmooth mechanics and applications, chapter Unilateral contact and dry friction in finite freedom dynamics[END_REF] allows to write the HSM conditions in terms of gap rate and contact impulse as given below:

   if g N > 0 then i N = 0 if g N = 0 then ġN = ( u1 -u2 ) • (-n 1 ) ≥ 0 , i N ≤ 0 and ġN i N = 0 (2)
where i N is the contact impulse and ui are the velocities. The heterogeneous asynchronous explicit-implicit time integrator proposed by [START_REF] Fekak | A new heterogeneous asynchronous explicitimplicit time integrator for nonsmooth dynamics[END_REF] is based on the same velocityimpulse formulation, inspired from Moreau's works.

In the tangential direction, in two dimension contact problems, the relative tangential velocity, also called slip rate, between the two nodes of a pair of nodes in contact is given by:

v T = ( u1 -u2 ) • t 1 (3)
where t 1 is the tangential vector at the interface relative to the slave side.

The frictional sliding is formulated according to Coulomb or Rate and State friction laws by introducing the slip function defined by ψ := |τ T | -µ |τ N |, in which τ T is the tangential stress and µ is the friction coefficient, constant in the case of Coulomb's law, and dependent on velocity and state variable in the case of the RS law. Then, the frictional sliding conditions are expressed as [START_REF] Popp | Dual mortar methods for computational contact mechanics -overview and recent developments[END_REF]:

ψ := |τ T | -µ |τ N | ≤ 0, v T + βτ T = 0, β ≥ 0, ψβ = 0 ( 4
)
where β is a scalar value. Finally, the behavior of the two contacting bodies is governed by the momentum equation in its strong form, complemented with the contact and stick-slip conditions. Let us note the whole domain Ω, including the two bodies:

Ω = Ω 1 ∪ Ω 2 .
It is assumed that both bodies Ω 1 and Ω 2 are linear elastic, characterized by the density ρ, Young's modulus E, and Poisson ratio ν. The time period of interest is denoted by

J = [0; T ].
The mechanical behaviour of the two sliding bodies is governed by the following system of equations:

                     ∇ • σ + b = ρü in Ω × J σ = C : ε in Ω × J ε = 1 2 ∇(u) + ∇(u) T in Ω × J Normal contact law (HSM) on Γ C × J Tangential friction law on Γ C × J (5)
where b is the body force per unit volume, σ and ε are the second-order stress and strain tensors, C is the fourth-order elastic constitutive tensor.

The Dirichlet and Neumann boundary conditions on Γ

D = Γ D 1 ∪Γ D 2 and Γ N = Γ N 1 ∪Γ N 2
are given by:

   u = u D on Γ D × J σ • n = t N on Γ N × J (6)
whereas u D is the Dirichlet prescribed displacement and t N the traction force.

Rate and State friction law

Experiments confirm that the friction coefficient depends on the slip rate and the state of interfacial contact areas. The dependence of friction on sliding history was observed for rocks as well as for many other materials. Dieterich [START_REF] Dieterich | Modeling of rock friction: 1. experimental results and constituve equations[END_REF][START_REF] Dieterich | Direct observation of frictional contacts: New insights for statedependent properties[END_REF] and Ruina [START_REF] Ruina | Slip instability and state variable friction law[END_REF] proposed an empirical rate and state-dependent constitutive formulation which has been extensively employed for analysis of fault slip in seismology. The coefficient of friction follows the empirical law:

µ = µ 0 + A ln v v + 1 + B ln θ θ + 1 ( 7
)
where v is the slip velocity, and θ is a state variable ; µ 0 , A and B are constant parameters which can be experimentally determined. The parameter v is the reference velocity and θ is a reference state value, verifying θ = D c /v , where D c is a characteristic sliding displacement required to stabilize friction coefficient after an abrupt change of sliding velocity. It can be interpreted as the required distance to renew surface contacts.

Indeed, it was established that this characteristic sliding distance is rather independent of slip velocity and normal stress, but would be well correlated with surface roughness.

The state variable describes the state of contact area between the contacting bodies, representing a measure of the age of continuously changing population of contacts. Ruina proposed the following aging law:

θ = 1 - θv D c (8)
The RS formula in Eq. ( 7) along with the aging law in Eq. ( 8), describe the experimental observations that, from a steady-state sliding characterized by an initial friction coefficient, a step increase in slip velocity at the interface results in an immediate jump in frictional force, followed by a friction decay depending on the displacement, and finally a stabilization at a new steady-state sliding friction at a new friction coefficient. There is velocity weakening, characterized by the new friction coefficient less than the initial one, if the parameters verify: A -B < 0.

Strong form of the PML

The sliding contact that has been investigated concerns the semi-infinite media. To appropriately mimic the sliding contact between semi-infinite media using a finite element procedure, the simplest way is to consider a large dimension in the normal direction to the interface. In such way, no reflected wave at the truncated limit of the mesh comes back towards the sliding interface within the duration of the simulation. Nonetheless, it requires very large mesh and leads to large computation time. Another way is to adopt appropriate absorbing conditions at the truncated boundaries of a significantly reduced mesh. Despite the computation time required by the introduced absorbing conditions, the computation time for the sliding contact simulation can be drastically reduced. Moreover, in terms of stability and regularization, the absorbing conditions can play a major role that we want to explore in the cases of Coulomb's friction and Rate and State friction laws.

Basu's formulation was adopted for the PML, expressed in the finite element framework [START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF][START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF]. Recently, this formulation was employed using Hybrid Asynchronous Time

Integration in order to propose appropriate PML with the capabilities to use implicit time integration and large time steps in comparison to the domain of interest in which the wave propagation is dealt using an explicit time integration scheme and a fine time step [START_REF] Brun | Hybrid asynchronous perfectly matched layer for seismic wave propagation in unbounded domains[END_REF]. First, the PML is formulated in a strong form in the frequency domain. It is defined by a complex-valued stretching functions λ i , replacing the real coordinates x i with the complex ones xi .

The complex-valued stretching function λ i is defined by:

∂ xi ∂x i = λ i (x i ) = 1 + f e i (x i ) -i f p i (x i ) bk s ( 9
)
where b denotes the characteristic length of the physical problem, k s = ω cs is the wavenumber, and c s is the S-wave velocity. The real-valued positive functions f e i and f p i vanish at the interface between the PML and the physical domain. The damping function f p i serves to attenuate the propagating waves in the x i direction, whereas the scaling function f e i attenuates the evanescent waves. In Eq. ( 9), the dependence of the complex term on the factor iω allows for an easy application of the inverse Fourier transform when expressing the PML in the time domain, resulting in a PML formulation independent on the frequency. In other words, all the frequencies are damped out in the same way.

The PML formulation is first obtained by modifying the governing equations defined in the frequency domain. For the sake of brevity, we consider in the following that the slave body Ω 1 is semi-infinite with an appropriate PML, wheareas the master body Ω 2 is finite. The formulation is easily extended to the case of two semi-infinite bodies involved in a sliding contact. The classical strong form of the equation of motion for an homogeneous isotropic medium under the plane strain assumption is written by substituting x i by xi as follows: [START_REF] Ma | Rupture dynamics on a bimaterial interface for dipping faults[END_REF] where C ijkl are the components of the fourth-order elastic constitutive tensor.

             j 1 λj (xj ) ∂σij ∂xj = -ω 2 ρu i σ ij = k,l C ijkl ε kl ε ij = 1 2 1 λj (xj ) ∂ui ∂xj + 1 λi(xi) ∂uj ∂xi
We introduce the following notations for the PML region: Ω PML is the region of the PML, bounded by the Γ PML = Γ D PML + Γ N PML , where Γ D PML ∩ Γ N PML = ∅, defining decomposition of the boundary between Dirichlet and Neumann conditions. In addition, g N denotes the prescribed tractions on Γ N PML . Thanks to the introduction of the stretching functions expressed in Eq. ( 9), the inverse Fourier transform can be easily applied to the previous frequency-domain equations, leading to the following time-domain equations [START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF][START_REF] Basu | Erratum to: Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF]:

               div(σ F e + Σ F p ) = ρf m ü + ρ cs b f c u + µ b 2 f k u in Ω PML × J σ = C : ε in Ω PML × J F eT ˙ F e + F pT ε F e + F eT ε F p + F pT E F p = ... 1 2 (∇ uT F e + F eT ∇ u) + 1 2 (∇u T F p + F pT ∇u) in Ω PML × J (11)
where the scalar-valued functions f m , f c and f k , and the second-order tensors F e , F p , F e and F p are all related to the PML stretching functions. Their expressions can be found in [START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF] and are summarized in this paper in Appendix A. The PML expressed for two-dimensional elastic waves in the time domain in Eq. ( 11) also involves time-integrals of the stress and strain tensors given by: Σ = t 0 σ dt and E = t 0 ε dt. For the boundary conditions, we assume that the displacements are equal to zero on the Dirichlet interface as well as zero traction force on the Neumann interface, given by:

   u = 0 on Γ D PML (σ F e + Σ F p ) • n = 0 on Γ N PML ( 12
)
The design of the PML is given by the thickness of the PML, denoted by L p , and the definition of the damping functions, f e i and f p i , which define the complex-valued stretching functions of the PML in Eq. ( 9). We consider polynomial functions in the thickness of the PML, with a degree equal to n, written as:

f α i = a α x i -x 0 L p n ; α = e, p; x 0 ≤ x i ≤ x 0 + L p ( 13
)
where a α is a positive integer, α being related to evanescent or propagating waves.

For a PML of thickness L p , Basu and Chopra [START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF] express the a α parameters as a function of the chosen target efficiency R related to the PML, by defining the efficiency as the ratio between the amplitude of the reflected waves and the amplitude of the incident waves:

a α = cp cs (n + 1) 2L p ln( 1 R ) ( 14
)
where c p is the P wave velocity. As a result, according to the above relationship, the efficiency of the PML is only controlled by the PML thickness L p and the degree n of the damping functions. Theoretically, in the continuous formulation, the reflected waves can be made arbitrary small, but it is not the case anymore after the space and time discretization built from the weak formulation presented in the next section.

3. Hybrid Asynchronous Time Integration of the semi-infinite bodies using PML

Weak form

The solution spaces and the test function spaces are defined for the two bodies: The solutions u 1 (t) and u 2 (t) in subdomains Ω 1 and Ω 2 are sought in the appropriate spaces

V 1 = u 1 ∈ [H 1 (Ω 1 )] d : u 1 = u D 1 on Γ D 1 and V 2 = u 2 ∈ [H 1 (Ω 2 )] d : u 2 = u D 2 on Γ D 2
, respectively, whereas the test space functions v 1 and v 2 belong to the spaces V 0 1 and V 0 2 , satisfying the zero value at the Dirichlet conditions. We consider that the slave body Ω 1 is semi-infinite, requiring the introduction of a PML defined in a third subdomain Ω 3 whose interface with the subdomain Ω 1 is noted as Γ p .

We introduce Lagrange multipliers, denoted by λ, in order to satisfy the normal contact and stick-slip conditions at the interface. The Lagrange multipliers belong to the appropriate spaces, corresponding to the adapted dual trace space denoted by M .

In addition, in order to connect the truncated semi-infinite slave body Ω 1 to its PML Ω 3 , the dual approach is again adopted, by introducing Lagrange multipliers λ p . Test functions related to the Lagrange multipliers are denoted by µ and µ p . Then, as proposed in [START_REF] Brun | Hybrid asynchronous perfectly matched layer for seismic wave propagation in unbounded domains[END_REF], it becomes possible to set up heterogeneous asynchronous strategy, enabling us to adopt a time integrator in the PML, independently from the choices made for the two contacting bodies. Indeed, we employ an explicit Central Difference scheme for the sliding contact region and an implicit constant acceleration average scheme for the PML. Now, the principle of virtual power for transient dynamics problems can be written for three subdomains: ∀t ∈ J, find the solution

u 1 (t) ∈ V 1 , u 2 (t) ∈ V 2 , u 3 (t) ∈ V 3 and λ(t) ∈ M , λ p (t) ∈ M p , for which the following weak form is satisfied ∀v 1 ∈ V 0 1 , ∀v 2 ∈ V 0 2 ,∀v 3 ∈ V 0 3 and ∀µ ∈ M , ∀µ p ∈ M p : Ω1 ρ 1 v 1 • ü1 dΩ + Ω1 ε(v 1 ) : σ 1 dΩ + Ω2 ρ 2 v 2 • ü2 dΩ + Ω2 ε(v 2 ) : σ 2 dΩ + ... Ω3 ρf m v 3 .ü 3 dΩ + Ω3 ρ cs b f c v 3 . u3 dΩ + Ω3 µ b 2 f k v 3 . u3 dΩ + ... Ω3 ε(v 3 ) e : σ 3 dΩ + Ω3 ε(v 3 ) p : Σ 3 dΩ = ... Ω1 v 1 • b 1 dΩ + Ω2 v 2 • b 2 dΩ + Γ N 1 v 1 • g N 1 dΓ + Γ N 2 v 2 • g N 2 dΓ + ... Γ C (v 1 -v 2 ) • λdΓ + Γ C µ • ( u1 -u2 )dΓ + Γp (v 1 -v 3 ) • λ p dΓ + Γp µ p • ( u1 -u3 )dΓ (15)
where the strain tensors ε(v 3 ) e and ε(v 3 ) p are derived from the strong form of the wave propagation in PML, in Eq. ( 11). They can be found in Appendix A.

The velocity continuity between the PML Ω 3 and the subdomain Ω 1 is enforced in a weak integral form:

Γp µ p .( u1 -u3 )dΓ = 0 (16)
In Eq. ( 15), the normal contact and tangential friction conditions are included in the In this paper, the weak form of the equilibrium equation in Eq. ( 15) is considered so as to build an explicit approach for the sliding contact with rate and state friction law for unbounded media separated by a frictional interface. As a consequence, no iterative method on the time step is required, but the time step size is limited to ensure the stability of the explicit time integration. Moreover, the hybrid approach using an explicit time integrator for the sliding contact and an implicit time integrator for the PML is adopted. Same time steps will be used for both time integrators in order to not alter the efficiency of the PML [START_REF] Brun | Hybrid asynchronous perfectly matched layer for seismic wave propagation in unbounded domains[END_REF].

In the following, the finite element discretization from the weak form in Eq. ( 15) will be summarized, by focusing on the constraint conditions concerning the PML and the sliding contacts as well as the interface power terms related to Lagrange multipliers λ p for the PML, λ N and λ T for the sliding contact.

The geometry of the three subodmains, their displacement and displacement time derivative, are discretized using the weak form of the equation of motion given in Eq. [START_REF] Brun | Hybrid asynchronous perfectly matched layer for seismic wave propagation in unbounded domains[END_REF]. We consider only matching meshes between the truncated slave body and its PML (interface Γ p ) as well as between the contacting bodies (interface Γ C ).

PML interface

The subdomains, including the elastic medium Ω 1 and the PML Ω 3 , are depicted in Figure 2 in the section devoted to the numerical application. As we consider matching meshes at the PML interface between the subdomains, the Lagrange multipliers can be assumed to be collocated at the common nodes between the subdomain Ω 1 and the PML.

As a result, contrarily to the mortar approach employed in [START_REF] Ehrl | A dual mortar approach for mesh tying within a variational multiscale method for incompressible flow[END_REF][START_REF] Klöppel | Fluid-structure interaction for nonconforming interfaces based on a dual mortar formulation[END_REF], the approximation of the Lagrange multiplier space does not involve shape functions defined at the interface, and the constraint matrices involved in the velocity continuity take the form of boolean matrices, with the only components equal to 1 for the degrees of freedom belonging to the interface. The weak constrain term Γp µ p • ( u1 -u3 )dΓ in Eq. ( 15) leads to the constraint relationship:

L (1) p U(1) + L (3) p U(3) = 0 on Γ p ( 17
)
where L (1) p and L (3) p are boolean matrices, operating on nodal vectors associated with the two subdomains Ω 1 and Ω 3 , and picks out the degrees of freedom lying on the interface Γ p in order to ensure the velocity continuity at the interface. The restricted velocity vectors on the interface are given by L (1) p U(1) and L (3) p U(3) , whereas the operators (L (1) p ) T and (L (3) p ) T represent the prolongation operators from the interface to the subdomains. The same relationships hold for the restricted velocity vectors on the interface related to the test functions v 1 and v 3 . According to the dual approach, Lagrange multipliers are introduced so as to ensure the velocity continuity at the PML interface. The interface power terms, related to the PML interface in Eq. ( 15), can now be expressed in a discrete format as:

Γp (v 1 -v 3 ) • λ p dΓ = (V (1) ) T (L (1) p ) T λ p + (V (3) ) T (L (3) p ) T λ p ( 18
)
where V (1) and V (3) denote the virtual nodal velocity vectors from the space discretization of the test functions v 1 and v 3 , and λ p is the Lagrange multiplier vector at the interface.

Contact frictional interface

The contact and frictional conditions are enforced in a weak form given by Eq. ( 15).

As the meshes are compatible at the contact interface Γ C , the discrete form of the normal conditions can be written for all the potential pairs of nodes in contact. After finite element discretisation of the two contacting bodies, we consider a set of pairs of nodes, with the index l = 1, • • • , N c , and N c being the total number of pairs of nodes potentially in contact. We have to introduce the constraint matrices L

C and L

C , selecting the nodes in contact such as we can write the gap vector:

g N = -L (1) C U (1) + L (2) C U (2) (19)
where the constant matrices L

C and L

C are boolean matrices. For a pair l of contacting nodes, the gap is scalar and is given by: g

l N = -X (1)
l + U

(1) l

• N

l + X

(2)

l + U (2) l • N (1)
l , where N

(1) l denotes the normal vector at the node l of the slave side, of Lagrangian coordinates X

(1) l and U (i)

l , i = 1, 2, are the displacement vectors of the two points in contact.

The contact conditions take the classical form of the Karush-Kuhn-Tucker conditions which can be directly related to the HSM condition [START_REF] Fekak | A new heterogeneous asynchronous explicitimplicit time integrator for nonsmooth dynamics[END_REF]. For all pairs of contacting nodes (∀l = 1, • • • , N c ), we have to check the HSM conditions:

               if g l N > 0 then λ l N = 0 if g l N ≤ 0 then          ġl N ≥ 0 λ l N ≥ 0 ġl N λ l N = 0 ( 20
)
with the discrete rate gap written for a given pair l of contacting nodes as: ġl

N = - U(1) l + U(2) l • N (1)
l . It remains to define the discrete tangential velocity for a pair l of contacting nodes:

v l T = U(1) l - U(2) l • T (1) l (21)
where T

(1) l

is the tangential vector related to the slave side.

For the sake of brevity, we gather the slave and master sides in order to express the gap vector, gape rate vector as well as the tangential velocity (slip rate) vector as:

         g N = L N U ġN = L N U v T = L T U ( 22
)
where the displacement vector U gathers all the displacements related to the two bodies in contact Ω 1 (slave) and Ω 2 (master), and the velocity vector U gathers all the velocities.

The constraint matrix L N associated with both subdomains is easily obtained from the constraint matrices L

C and L

C by using the previous expression of the gap in Eq. ( 19). In the same manner as for normal constraints, the constraint matrix L T associated with both subdomains is easily obtained from the pointwise expression of the tangential velocity in Eq. ( 21).

Time discretisation of the explicit-implicit HATI for frictional contact problems including PML

From the weak form in Eq. [START_REF] Brun | Hybrid asynchronous perfectly matched layer for seismic wave propagation in unbounded domains[END_REF], it is now possible to express the discrete equations of motion according to the Central Difference scheme over a time step ∆t = [t n ; t n+1 ]. The two bodies Ω 1 and Ω 2 are gathered in the same subdomain Ω to simplify the notations, whereas the PML subdomain is always denoted by Ω 3 .

                     M Ün+1 + C Un+1/2 + KU n+1 = f ext,n+1 + f cont,n+1 + (L p ) T λ p,n+1 M (3) Ü(3) n+1 + (C (3) + C (3)
) U(3) n+1 + (K (3) + K (3)
)U

(3)

n+1 + P (3) n = (L (3) p ) T λ p,n+1 L p Un+1 + L (3) p U(3)
n+1 = 0 Normal contact conditions HSM Tangential frictional conditions [START_REF] Di Bartolomeo | Wave and rupture propagation at frictional bimaterial sliding interfaces: from local to global dynamics, from stick-slip to continuous sliding[END_REF] In Eq. ( 23), the first equation represents the discrete equation of motion for the two contacting bodies, which will be integrated in time adopting the explicit Central Difference scheme. Following the classical explicit time integration procedure, the mass matrix related to the two bodies M is lumped and the viscous term C Un+1/2 is written at the mid-step [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]. At the right-hand side of the first equation, the contact force f cont,n+1 is deduced from the normal-contact and tangential-frictional conditions. The second equation is the discrete equation of motion associated with the PML whose matrices

M (3) , C (3) , C (3)
, K (3) and K [START_REF] Dieterich | Direct observation of frictional contacts: New insights for statedependent properties[END_REF] , as well as the vector P (3) n , depending only on strain and stress quantities known at the time t n , are defined in Appendix A. The PML region is integrated with an implicit time integrator (Constant Average Acceleration), thanks to the hybrid explicit/implicit approach [START_REF] Brun | Hybrid asynchronous perfectly matched layer for seismic wave propagation in unbounded domains[END_REF]. Finally, the third equation is the discrete velocity continuity at the PML interface, ensured by introducing the Lagrange multipliers λ p,n+1 on the right hand side of both equations of motion.

For the contacting bodies, the approach proposed by [START_REF] Fekak | A new heterogeneous asynchronous explicitimplicit time integrator for nonsmooth dynamics[END_REF], based on the velocityimpulse formulation for frictionless contacts, is extended to the case of frictional contacts. The approach is inspired from Moreau and Jean's works in nonsmooth dynamics [START_REF] Moreau | Nonsmooth mechanics and applications, chapter Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Moreau | Numerical aspects of the sweeping process[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF], setting up event-capturing methods which allow to take into account an accumulation of impacts in a given time step ∆t, contrarily to the family of event-tracking time integration methods which require the accurate detection of the impact times. For this purpose, the acceleration, which is not finite when the first impact occurs due to the velocity jump, is expressed in terms of mid-step velocity at time t n+1/2 , known on the current time step ∆t, and the mid-step velocity at time t n+3/2 related to the following time step: Ün+1 = Un+3/2 -Un+1/2 ∆t . Thus, after multiplying the first equation in Eq. [START_REF] Di Bartolomeo | Wave and rupture propagation at frictional bimaterial sliding interfaces: from local to global dynamics, from stick-slip to continuous sliding[END_REF] by the time step, we obtain : [START_REF] Curnier | Unilateral contact[END_REF] where I n+1 = ∆tf cont,n+1 is the impulse vector concerning all pairs of nodes of the interface potentially in contact. The above velocity-impulse formulation allows to deal with jump of velocities for the first impact between nodes as well as permanent contact between nodes during a prolonged contact stage. The impulse vector can be written in terms of transpose of previously defined constraint matrices and Lagrange multipliers introduced to ensure the normal contact conditions as well as the tangential frictional conditions:

M Un+3/2 = M Un+1/2 + ∆t f ext,n+1 -C Un+1/2 -K Un+1 + I n+1 + ∆t (L p ) T λ p,n+1
I n+1 = (L N,n+1 ) T λ N,n+3/2 + (L T,n+1 ) T λ T,n+3/2 (25)
First, the Lagrange multipliers λ N,n+3/2 related to the contact conditions is derived from the HSM conditions. The equation of motion is multiplied by L N,n+1 M -1 to obtain the reduced size interface equation involving only normal degrees of freedom at the interface:

H N λ N,n+3/2 = b N,n+3/2 (26)
where the operator H N is the Steklov-Poincaré (or Delassus) operator and the right-hand side vector are defined by:

   H N = L N,n+1 M -1 (L N,n+1 ) T b N,n+3/2 = -L N,n+1 Un+1/2 + ∆tM -1 f ext,n+1 -CU n+1/2 -K Un+1 (27)
After obtaining the Lagrange multipliers related to the normal conditions, the Lagrange multipliers related to the friction can be obtained as follows:

H T λ T,n+3/2 = b T,n+3/2 (28)
with

   H T = L T,n+1 M -1 (L T,n+1 ) T b T,n+3/2 = -L T,n+1 Un+1/2 + ∆tM -1 f ext,n+1 -CU n+1/2 -K Un+1 (29)
According to frictional conditions for Coulomb's law or RS law, the negative value of the slip function has to be checked. We give the expression of the slip function for a pair l of nodes, for which the normal and tangential Lagrange multipliers are scalar values:

We define

ψ l := λ l T,n+3/2 -µ(v l n+1/2 , θ l n+1 ) λ l N,n+3/2
, where the friction coefficient µ(v l n+1/2 , θ l n+1 ) depends on the local tangential velocity v l n+1/2 , known at the current time step, and the state variable θ l n+1 governed by the aging law. The aging law given in Eq. ( 8) is integrated implicitly in time according to a classical Euler forward scheme as follows:

θ l n+1 -θ l n ∆t = 1 - θ l n+1 v l n+1/2 Dc (30)
The local mid-step velocity v l n+1/2 at the current time step is employed to calculate the state variable θ l n+1 at the following time step as well as the estimated new friction coefficient µ(v l n+1/2 , θ l n+1 ). The initial value θ(t = 0) of the state variable is assumed to be equal to zero.

After calculating the new value of the friction coefficient, the frictional conditions µ(v l n+1/2 , θ l n+1 ) should be taken into account to determine whether the tangential state corresponds to a stick or a slip state for every pair of nodes l in contact:

     if ψ l ≤ 0 then v l T = 0 (stick state) if ψ l > 0 then v l T = 0 and λ T,n+3/2 = -b l T ,n+3/2 b l T ,n+3/2 µ(v l n+1/2 , θ l n+1 )λ N,n+3/2 (slip state) (31)
where the term

-b l T ,n+3/2 b l T ,n+3/2
gives the sign of the friction impulse.

Numerical examples

All the numerical applications presented in the following are developed in Matlab.

The time integration is explicit (Central Difference scheme) for the domain of interest, that is the elastic solid medium with the frictional interface, while the PML is integrated with an implicit time integration (Trapezoidal Rule), according to the HATI framework (Heterogeneous Asynchronous Time Integrator). Same time step size is used in both partitions.

Stick and slip test of a 2D unbounded solid sliding on a rigid flat surface

The case of a semi-infinite solid half-space sliding on a rigid flat surface is studied following a two dimensional analysis. The 2D solid is assumed to be linear elastic, under the plane strain assumption. This case was recently studied by Kammer et al.

[8] in order to highlight the existence of a critical length when Prakash-Clifton law [START_REF] Prakash | Time resolved dynamic friction measurements in pressure-shear[END_REF] is employed to regularize Coulomb's friction law. This test was also investigated by Radiguet et al. [START_REF] Radiguet | Survival of heterogeneous stress distributions created by precursory slip at frictional interfaces[END_REF] using a slip weakening friction law on the frictional interface to study the dynamics of successive slip events. Prakash-Clifton regularization was widely employed to alleviate the grid size dependence observed when dissimilar materials were in contact [START_REF] Cochard | Fault rupture between dissimilar materials: Ill-posedness, regularization, and slip-pulse response[END_REF][START_REF] Ranjith | Slip dynamics at an interface between dissimilar materials[END_REF][START_REF] Di Bartolomeo | Wave and rupture propagation at frictional bimaterial sliding interfaces: from local to global dynamics, from stick-slip to continuous sliding[END_REF][START_REF] Coker | Frictional sliding modes along an interface between identical elastic plates subject to shear impact loading[END_REF][START_REF] Dedontney | Finite element modeling of branched ruptures including off-fault plasticity[END_REF]. Here, as the effect of PML on the grid size dependence is investigated, the Prakash-Clifton regularization is not incorporated into the following simulations. The original problem proposed by Kammer et al. The initial remote loads are applied on the top of the solid so as to generate an initial uniform stress state at the interface with a normal pressure p equal to 0.0150 Pa and a shear load t equal to 0.0105 Pa, corresponding to shear-pressure ratio equal to 0.7. The solid is in equilibrium at the initial stage, with a sticking state along all the interface points, by considering Coulomb's law with a constant friction coefficient µ equal to 0.75.

As a result, the initial shear state at the interface corresponds to a ratio t/p between the frictional force and the normal force equal to 0.7, that is slightly less than the Coulomb's frictional strength. It enables us to propagate an interface shear rupture, nucleated from a fixed area around the point located at X = 10 m. The nucleation procedure is the same as the one employed by different authors [START_REF] Kammer | The existence of a critical length scale in regularised friction[END_REF][START_REF] Cochard | Fault rupture between dissimilar materials: Ill-posedness, regularization, and slip-pulse response[END_REF][START_REF] Andrews | Wrinkle-like slip pulse on a fault between different materials[END_REF]: It is introduced during the dynamics simulation by decreasing locally the normal contact pressure on an area around the nucleation point, whose size and intensity are progressively increased in time. It aims at reproducing the smooth increase in intensity and in space of a fluid pressure source, acting at the interface. The space and time increase of this fluid pressure is described by an elliptical region in the x -t plane. It has been determined that a maximum decrease of 50% with respect to the initial pressure value p is sufficient to induce a slip pulse from the nucleation area in the case of Coulomb's law.

The solid is discretized using regular square finite elements, with linear shape function and four Gauss points. The coarsest mesh is characterized by a uniform finite element size L FE equal to 0.5 m, with a stable time step at about 50% of the CFL critical time step (∆t = 0.10 s). The PML has a thickness L p , corresponding to 10 finite elements in the thickness. We study the effect of the refinement of the mesh, characterized by the parameter nb, varying from 1, for the coarsest mesh, to 10, for the finest mesh, the associated time step being divided by the same parameter in order to ensure the stability of the CD explicit time integration.

The set of equations in Eq. ( 23) is integrated in time using the same time step for the body in contact with the rigid flat plane and its PML. Using the previously presented HATI approach, the PML is integrated using an implicit time integration scheme (Constant Average Acceleration) whereas the contacting truncated mesh is integrated explicitly using the Central Difference scheme.

Preliminary loading stages using PML

Two preliminary stages are considered before starting the nucleation procedure: First we apply dynamically the normal pressure, increasing progressively over the duration T p = 20s, and then, we let the remote pressure equal to the target constant level p = 0.0150P a over the same duration T p = 20s; in the second stage, we apply the same procedure with the remote shear loading t over the duration Tt = 20s and then, we impose a constant target shear load ( t = 0.0105P a) over the same duration Tt = 20s.

Thus, the complete duration of the preliminary stages is equal to T 0 = 80s, assuming during this time period a constant friction coefficient (Coulomb's law). Both pressure and shear loadings are applied at a line of the mesh just below the PML interface (see Figure 2).

Global normal and tangential forces are plotted versus time in Figure 3. In Figure 3(a), the time history of the global forces is obtained with the PML. It can be remarked that the equilibrium can be reached very quickly thanks to the dissipative effect provided by the PML.

For the purpose of comparison, the stabilization in normal and tangential forces is also investigated by adopting a strong viscous damping in the bulk, without PML. Thus, for this second strategy, a Rayleigh viscous matrix C is introduced into the equation of motion of the solid without PML, by constructing it from the mass matrix (C = αM). Indeed, it is known that this choice does not decrease the stable time step of the Central Difference scheme, when the viscous term is taken into account at mid-step (at

t n+1/2
) into the equilibrium equation [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF], as done in this study and detailed in previous sections. The parameter α is computed from a representative duration corresponding to the duration T P -wave needed by the P-waves to propagate back and forth into the solid (from the line, at which forces are imposed, to the interface and back to the force line). So, we take α = ξω 0 , where ω 0 is the pulsation associated with the representative duration:

ω 0 = 2π
T P -wave . In Figure 3(b), the stabilization in time of the global forces during the two preliminary stages is plotted, including a Rayleigh viscous matrix associated with a strong damping coefficient (ξ = 0.3). It can be observed that the PML efficiently damps out the pressure and shear waves, reaching more quickly the target constant values than in the case of the Rayleigh viscous damping matrix, despite the fact that no damping is taken into account in the solid bulk when PML is considered. propagates along the interface in the positive direction of X, at a super-shear velocity, until it stops naturally. Similar results were obtained in [START_REF] Kammer | The existence of a critical length scale in regularised friction[END_REF]. Finally, the slip rate at a point of the interface located at 6m from the nucleation point is plotted for the different meshes in Figure 6. It can be seen that, the finer is the mesh, more spikes appear with increasing peak values, indicating that the sliding problem needs to be regularized [START_REF] Kammer | The existence of a critical length scale in regularised friction[END_REF]. Indeed, it is well known that Coulomb's friction is ill-posed, generating the grid-dependence, illustrated by Cochard and Rice in [START_REF] Cochard | Fault rupture between dissimilar materials: Ill-posedness, regularization, and slip-pulse response[END_REF]. The RS friction law is now considered with the same meshes and preliminary loading stages as in previous section. The RS friction law is given in Eq. ( 7), with the parameters:

µ 0 = 0.5, v = 0.1, D c = 10 -6 , A = 0.005, B = 0.075µ 0 and θ = D c /v . Moreover,
to prevent singularities at low velocities, the value of the friction coefficient is taken as equal to 0.75 for velocities less than 10 -4 m/s. The steady state curve related to the RS friction law is shown in Figure 7. The slip velocity after nucleation is between 10 -2 m/s and 10 -3 m/s, indicating that the rate and state regime is a velocity-weakening regime. on different degrees of refinement (mesh parameter nb equal to 4 and 7) as well as the time history of the slip rate pulse at three different locations, X = 20m, X = 25m and X = 30m. It can be seen that a slip rate pulse is generated from the nucleation area in the region around the point X = 10m, and propagates with a rupture velocity, denoted by V r , assessed at about 1.4m/s, which is comprised between the S wave velocity (1m/s) and the P wave velocity (1.73m/s). Thus, as in Coulomb's case, the rupture propagates with a super-shear velocity, but, this time, it does not decelerate: The slip rate pulse propagates in a self-sustaining manner on all the frictional interface, up to the right periodic condition, at about the time 102s. The self-sustained dynamic shear is due to the destabilizing effect of the RS friction law in a velocity-weakening regime. In addition, it can be observed that the intensity of the pulse slightly increases as the pulse propagates. Finally, numerical disturbances after the passage of the main wave are clearly highlighted: After this first peak, additional peaks are created, which propagate with the same velocity along the interface. Similar observations were made in the literature, noting the high sensitivity of the RS friction law. It is particularly the case in explicit simulations requiring very small time steps compared to the CFL condition [START_REF] Shi | Properties of dynamic rupture and energy partition in a solid with a frictional interface[END_REF][START_REF] Ben-Zion | Dynamic simulations of slip on a smooth fault in an elastic solid[END_REF][START_REF] Lapusta | Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate-and state-dependent friction[END_REF],

but also in implicit simulations [START_REF] Tal | Dynamic mortar finite element method for modeling of shear rupture on frictional rough surfaces[END_REF] with a surprisingly short time step size, meeting to be of the order of the CFL condition. In Figure 10, snapshots of the slip rate in the bulk are displayed for three different times. It can be observed that the curved shape in the bulk, displayed by velocity isovalues, is triggered by the nucleation procedure at the time 85s. The shape of the velocity pulse in the bulk tends to be more straight at the two following times 90s and 94s. It is of particular interest that there is no visible spurious reflections at the PML interface, enabling us to model the shear rupture propagation with a relatively small dimension of the solid in the normal direction to the frictional interface. Then, the energies during the simulation are investigated in order to assess the energy of the shear rupture at the interface compared with the elastic energy in the solid, as done by Shi and co-authors [START_REF] Shi | Properties of dynamic rupture and energy partition in a solid with a frictional interface[END_REF]. The ratio between the kinetic energy in the bulk with the elastic energy is plotted in time in Figure 11 for the two preliminary stages, followed by the shear rupture simulation with the RS friction law. It is referred in the literature as the seismic efficiency. During the two preliminary stages, it can be confirmed that the remote pressure and shear loadings have been applied dynamically in a very progressive way, with the kinetic energy very low (less than 1% of the elastic energy) during the loading stages (0s ≤ t ≤ 20s and 40s ≤ t ≤ 60s), and becoming close to zero during the plateau stages for the stabilization (20s ≤ t ≤ 40s and 60s ≤ t ≤ 80s). In the shear rupture stage, triggered by the nucleation procedure, it is underlined that the pulse-like rupture is associated with very low level of kinetic energy, assessed at 0.25% at the time t = 102s when the pulse reaches the right periodic boundary condition. Such as very low seismic efficiency (less than 1%) has been observed in the literature for pulse-like rupture at the interface with a single pulse [START_REF] Shi | Properties of dynamic rupture and energy partition in a solid with a frictional interface[END_REF]. The ratio between the interface energy with the elastic energy is also plotted in Figure 11: its value, computed from the Lagrange multipliers, is equal to zero during the pressure and shear preliminary stages and increases during the rupture stage, reaching at about 4% at the time t = 102s. In Figure 14, it can be seen that the slip rate pulse propagates along the interface without spurious oscillations up to reach the right periodic boundary condition at time t = 102s. The time history of the slip rate pulse exhibits an unchanged pulse, propagating at a constant rupture velocity at about 1.4m/s, with the same intensity and without secondary pulses. The rupture wave is displayed in terms of nodal normal and tangential contact forces in Figure 15 for three different times. The shape of the wave is not altered anymore by the numerical noise coming from the RS friction law. Finally, the time history of the tangential shear force is given in Figure 16 at two different points (X = 20m and X = 30m) belonging to the interface. In comparison to the previous test where the PML is located at 4m from the frictional interface, this test shows that the numerical disturbances inherent to the friction law are filtered out.

PML

We consider the same test but, this time, without the PML which is replaced by an elastic medium. The size of the elastic medium is sufficient to ensure that there are no waves reflected back from the top boundary of the elastic subdomain during a period of interest. After this period, reflected waves trigger secondary slips and pulses at the interface. The fine mesh is adopted with the mesh parameter nb = 7. We start with the HATI line located to eight Finite Elements from the frictional interface. In other words, the contacting subdomain is composed of eight rows of Finite Elements, coupled with a sufficiently large elastic subdomain on the top. Without PML, the Rayleigh damping, proportional to the mass matrix, is employed in the preliminary stages to achieve a stabilized equilibrium at the interface as explained in section 4.2, after which, the Rayleigh damping is deactivated and the nucleation procedure is introduced. The time history of the slip rate at the interface is displayed in Figure 18. Again we observe that the numerical disturbances appear after the passage of the main pulse in Figure18(a)

and (b). In addition, in comparison to the previous case with PML, we can see in Figure The shape of the wave is not altered anymore by the numerical noise coming from the RS friction law.

As a consequence, it can be concluded that the HATI coupling, located here very close to the interface, provides a very desirable filtering effect for the RS friction law in explicit dynamics. It constitutes a useful method to filter this type of friction law, at the expense of decomposing the problem at hand in two subdomains, one subdomain composed of one row of elements along the frictional interface and the other subdomain for the remaining of the medium. Obviously, for unbounded domains, the PML is needed, otherwise reflected waves from the boundaries disturb the interface behaviour. In the case where we only focus on the interface behaviour, the PML can be located at one element from the frictional interface as presented previously.

Conclusion

We extend the new heterogeneous asynchronous time integrator for nonsmooth dynamics in 2D small deformation to the case of frictional interfaces governed by Rate and State friction laws. The time integration of the frictional contact adopts a dual formulation, using Lagrange multipliers derived from the normal and tangential conditions expressed in terms of velocity. The area concerning the sliding contact is integrated explicitly using the explicit Central Difference (CD) scheme, whereas unbounded domains are taken into account thanks to the introduction of Hybrid Asynchronous Perfectly Matched Layers (HA-PML) at the boundaries of the truncated mesh, integrated using the implicit Constant Average Acceleration (CAA) scheme.

We apply our method to the case of a 2D sliding block on a flat rigid plane, initially compressed and sheared by remote loadings, and subsequently subjected to a nucleation procedure at the interface, triggering a shear rupture along the interface. First, during the preliminary pressure and shear loadings to reach the initial state before the nucleation, the advantages of the PML, with respect to classical Rayleigh viscous damping introduced into the bulk, are highlighted. Then, it has been been shown that the nucleation procedure, in the case of RS friction law, generates a self-sustaining single pulse, propagating along the interface with a supershear constant rupture velocity. As already observed in the literature in simulations of shear rupture with RS friction laws, numerical oscillations following the first main pulse have been highlighted, becoming secondary pulses with increasing intensity as the rupture wave propagates. We show that these numerical oscillations are filtered out by a PML close to the RS interface with the HATI coupling. The rupture wave in terms of slip rate, normal and shear stresses, propagates along the interface without changing its shape and the obtained quantities related to the interface, describing the propagation of the supershear wave rupture, converge with respect to the space and time refinement. In addition to this optimal test when we only focus on the interface behaviour, the filtering effect of the HATI coupling has been highlighted. As a result, the HATI approach is particularly useful when adopting nonlinear RS friction formulation, which is highly dependent on interfacial slip velocity and very prone to the generation of numerical disturbances. The perspective of this work is to apply this method for unbounded domains in contact with different material parameters.

where N I is the shape function related to the node I and I d is the identity matrix.

Taking into account Eq. ( 34), the internal force term p e n+1 , derived from the energy terms Ω εe : σ dΩ + Ω εp : Σ dΩ in Eq. ( 34), can be written as:

For the time stepping, Basu and Chopra assume the relationship between stress and time-integral stress vectors:

Σn+1 = Σn + dt σn+1 (40)
It follows:

p e n+1 = Ωe BT σn+1 + Ωe BpT Σn (41)
where:

BT = BeT + dt BpT (42)
Additional assumptions are required for the time stepping procedure:

   εn+1 = ε n+1 -εn dt E n+1 = E n + dt ε n+1 (43)
Using the assumptions given in Eq. ( 43), the third equation of the system in Eq. [START_REF] Fekak | A new heterogeneous asynchronous explicitimplicit time integrator for nonsmooth dynamics[END_REF] leads to the expression of the strain εn+1 at the end of the time step [START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF]:

εn+1 = 1 dt B ε Un+1 + B Q U n+1 + 1 dt Fε εn -FQ Ên ( 44
)
where B ε , B Q , Fε and FQ are matrices given in the following. The matrix B is given in terms of the submatrix for the node I by:

B I =      F 11 N l I1 F 21 N l I1 F 12 N l I2 F 22 N l I2 F 11 N l I2 + F 12 N l I1 F 21 N l I2 + F 22 N l I1      (45)
with:

F l = F p + F e dt -1 F = F e F l F Q = F p F l (46)
and

N l Ii = F l ij N I,j (47)
The matrix B Q is defined, similarly, by replacing F with F Q . Finally, the matrix F is given by: Likewise, the matrix FQ is obtained by replacing F with F Q .

F =      (F
The stress σn+1 at the end of the time step is computed from the strain εn+1 thanks to the elastic constitutive matrix. The element-wise internal force p e n+1 can now be written in terms of the element vectors of velocity and displacement ( Ue n+1 and U e n+1 ) as well as a term, denoted P(ε e n , E e n , Σ e n ), depending only on strains and stresses known at the beginning of the time step t n = ndt. The element-wise internal force p e n+1 is written as: where the element matrices c e and k e are defined as:

Figure 1 :

 1 Figure 1: Configurations of the slave body Ω 1 against the master body Ω 2 . Left: Reference configuration;

  term Γ C µ • ( u1 -u2 )dΓ. It can be decomposed into two constraint equations related to the contact and friction conditions, involving the test function µ N related to the normal Lagrange multipliers λ N and the test function µ T related to the tangential Lagrange multipliers λ T . Both Lagrange multipliers are scalar-valued functions defined on the interface Γ C .

[ 8 ] is displayed in Figure 2 :

 82 The longitudinal dimension L is equal to 40 m and the truncated height H equal to 4 m, where the PML will be incorporated to mimic the semi-infinite medium. It has to be noted that Kammer et al.[START_REF] Kammer | The existence of a critical length scale in regularised friction[END_REF] employed an extended mesh in the vertical direction with a large height H equal to 20 m to prevent reflected waves from influencing the interface behaviour. The PML thickness, denoted by L p , has to be chosen depending on the target performance R of the PML, as given in the design relationship of the PML in Eq. (14). The elastic parameters are: Young's modulus E = 2.5 Pa, Poisson's ratio ν = 0.25 and density ρ = 0.25 kg/m 3 , leading to the P wave velocity c P = 1.73 m/s and the S wave velocity c S = 1 m/s. On the left and right sides of the solid, we assume periodic boundary conditions.

Figure 2 :

 2 Figure 2: Two-dimensional test of a semi-infinite linear elastic solid of length L = 40 m, initially subjected to remote pressure and shear loadings, and subsequently nucleated artificially at X = 10 m (at the star), generating a propagating slip along the interface. The PML is introduced at an height H and Periodic Boundary Conditions are considered on the left and right sides.

Figure 3 : 4 . 3 .

 343 Figure 3: Time-history of the global normal contact force and tangential frictional force at the interface during the first two preliminary loading stages with the remote pressure p and shear t; (a) With PML (without Rayleigh viscous damping matrix); (b) Without PML, but with Rayleigh viscous damping matrix

Figure 4 :

 4 Figure 4: Time history of the slip rate along the interface (0m ≤ X ≤ 40m)

Figure 5 :

 5 Figure 5: Isovalues of the velocities into the solid and PML at three different time instants.

Figure 6 : 4 . 4 .

 644 Figure 6: Time history of the slip rate at a point belonging to the interface governed by Coulomb friction law (at X = 16m, that is 6m, on the right from the nucleation point) as function of the refinement in space (nb equal to 1, 4, 7, 10 corresponding to finite element size equal to 0.5m, 0.125m, 0.071m, 0.05m, respectively) and time (CFL condition for the time step size)

  law RS law with threshold at low velocities

Figure 7 :

 7 Figure 7: Rate and state law and threshold at the low velocities

Figure 8 :Figure 9 :

 89 Figure 8: Time history of the slip rate along the interface (0m ≤ X ≤ 40m) and time history slip rate at three different locations along the interface (X = 20m, X = 25m, X = 30m) for the mesh parameter nb = 4

Figure 10 :

 10 Figure 10: Isovalues of the slip rate into the solid and PML (edges in black) at three different times

Figure 11 :

 11 Figure 11: Percentage of kinetic energy and sliding contact energy with respect to the maximum of the internal energy

Figure 12 :

 12 Figure 12: Contact (in red) and tangential (in blue) forces along the interface at three different times

Figure 13 :

 13 Figure 13: Time history of the tangential force at two points of the interface (X = 20m and X = 30m)

Figure 14 :

 14 Figure 14: Time history of the slip rate along the interface (0m ≤ X ≤ 40m) and time history of the slip rate at three different locations along the interface (X = 20m, X = 25m, X = 30m) for the mesh parameter nb = 7 with a PML close to the interface

Figure 15 :

 15 Figure 15: Contact (in red) and tangential (in blue) forces along the interface at three different times

Figure 16 :Figure 17 : 4 . 6 .

 161746 Figure 16: Time history of the tangential force at two points of the interface (X = 20m and X = 30m)

  18(c) that reflected waves from the top of the elastic medium trigger new slips and pulses at the interface at around the time t = T 0 + 16s. Then these new pulses propagate along the interface as shown in Figure18(a) and (b). The development of the secondary peaks, characteristic of the RS friction law, is clearly visible on the shear rupture wave, plotted at time t = T 0 + 15s, shown in 18(d).

Figure 18 :

 18 Figure 18: HATI line located to eight Finite Elements from the frictional interface (mesh parameter nb = 7): (a) Time history of the slip rate along the interface (0m ≤ X ≤ 40m) with red arrows indicating the secondary pulses ; (b) Time history of the slip rate at three different locations along the interface (X = 20m, X = 25m, X = 30m) with red arrows indicating the secondary pulses ; (c) Slip rate isovalues at time T 0 + 16s with red circle indicating the propagation in the solid of the secondary pulses ; (d) Contact (in red) and tangential (in blue) forces along the interface at time t = 95s

Figure 19 :

 19 Figure 19: HATI line located to one Finite Element from the frictional interface (mesh parameter nb = 7): (a) Time history of the slip rate along the interface (0m ≤ X ≤ 40m) ; (b) Time history of the slip rate at three different locations along the interface (X = 20m, X = 25m, X = 30m) ; (c) Slip rate isovalues at time T 0 + 16s ; (d) Contact (in red) and tangential (in blue) forces along the interface at time t = 95s

Both

  Voigt notation is adopted for the stress and time-integral stress tensors, giving the following vectors: σ = Be and Bp matrices depend on the derivatives of the shape functions and damping functions expressed in terms of their nodal submatrices as: e ji N I,j et Ñ p Ii = F p ji N I,j

  p e n+1 = c e U e n+1 + k e U e n+1 + P(ε e n , E e n ,

  22 2F 11 F 12 2F 21 F 22 F 11 F 22 + F 12 F 21

	
	   	(48)

11 ) 2 (F 21 ) 2 F 11 F 21 (F 12 ) 2 (F 22 ) 2 F 12 F

Appendix A

In the strong form of the wave equation within a PML, the following tensors, depending on the scalar-valued scaling f e and damping f p functions, are involved:

with:

Let v be the test function belonging to an appropriate space, the weak formulation is obtained by integrating over the computational domain Ω:

Γ = ∂Ω is the boundary of Ω, and n the outer unit normal vector related to Γ. The strain tensors εe and εp depend on the scaling and damping functions of the PML as follows [START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF]:

After the finite element disceretization, we define the mass, the damping and the stiffness matrices, given by blocks of size 2 × 2:

)

D is the material constitutive matrix under the plane strain assumption expressed as:

where µ L is the shear modulus (or Lamé's second coefficient) and κ is the bulk modulus.

Finally, on the whole PML, the discrete form of the equation of motion at the end time t n+1 of the time step h = [t n t n+1 ] can be expressed as:

in which the global matrices C and K, related to velocities and displacements, respectively, are derived from the classical assembly procedure of the previous elementary matrices c e and k e ; P(ε e n , E e n , Σ e n ) is the known part of the internal force at the beginning of the time step; M, C and K are the global matrices resulting from the assembly of the element matrices given in Eq. [START_REF] Lapusta | Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate-and state-dependent friction[END_REF].