David Dureisseix

Paul Larousse

Anthony Gravouil

Jean Di Stasio

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

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

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 jump on the 2D interface surface with no thickness, and interfacial force distributions acting on the bodies apart the interface. The proposed test cases which can serve as benchmarks for simulation codes, 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.

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 [START_REF] Brogliato | Numerical simulation of finite dimensional multibody nonsmooth mechanical systems[END_REF][START_REF] Acary | On solving contact problems with Coulomb friction: Formulations and numerical comparisons[END_REF]. 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 [START_REF] Moreau | Indetermination due to dry friction in multibody dynamics[END_REF][START_REF] Moreau | Facing the plurality of solutions in nonsmooth mechanics[END_REF][START_REF] Alart | How to overcome indetermination and interpenetration in granular systems via nonsmooth contact dynamics. An exploratory investigation[END_REF]). 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 explicit scheme, namely the CD-Lagrange scheme [START_REF] Fekak | A new heterogeneous asynchronous explicit-implicit time integrator for nonsmooth dynamics[END_REF], embedding by construction the treatment of non-smooth evolutions, without the need for any penalization with artificial numerical parameters.

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 [START_REF] Jean | Non-smooth contact dynamics approach of cohesive materials[END_REF][START_REF] De Saxcé | A non incremental variational principle for brittle fracture[END_REF].

The CD-Lagrange scheme embedding the non-smooth Moreau framework [START_REF] Moreau | Une formulation du contact à frottement sec ; application au calcul numérique[END_REF][START_REF] Moreau | On Unilateral Constraints, Friction and Plasticity[END_REF][START_REF] Dubois | The Contact Dynamics method: A nonsmooth story[END_REF]) involves an interface behavior expressed not with a classical stress vs displacement jump function, but with an impulse-velocity jump relationship. It has been tested and used for different cases of non-regularized frictional contact [START_REF] Fekak | A new heterogeneous asynchronous explicit-implicit 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].

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 [START_REF] Meyer | Optimal control of perfect plasticity. Part I: stress tracking[END_REF]. Here, we focus on the plasticity developed in the interface model, for two phenomenological versions corresponding to two physical situations, and to derive benchmark 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 [START_REF] Francfort | Quasistatic brittle fracture seen as an energy minimizing movement[END_REF][START_REF] Kiener | 100 years after Griffith: From brittle bulk fracture to failure in 2D materials[END_REF], 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 [START_REF] Jean | Non-smooth contact dynamics approach of cohesive materials[END_REF]. The infinite stiffness in traction of the interface (before fracture) is a second one [START_REF] Collins-Craft | On the formulation and implementation of extrinsic cohesive zone models with contact[END_REF] 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, [START_REF] Raffa | A micromechanical model of a hard interface with micro-cracking damage[END_REF][START_REF] Lebon | On the emergence of adhesion in asymptotic analysis of piecewise linear anisotropic elastic bonded joints[END_REF]), 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 [START_REF] Wriggers | Computational Contact Mechanics[END_REF]). These cases will be described, implemented and tested in this article. The target cases considered herein deal with a constant time step Δ𝑡, 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 [START_REF] Ben Belgacem | The mortar finite element method for contact problems[END_REF][START_REF] Casadei | A mortar spectral/finite element method for complex 2D and 3D elastodynamic problems[END_REF][START_REF] Wohlmuth | A mortar finite element method using dual spaces for the Lagrange multiplier[END_REF][START_REF] Popp | Dual mortar methods for computational contact mechanicsoverview and recent developments[END_REF][START_REF] Pinto Carvalho | An efficient algorithm for rigid/deformable contact interaction based on the dual mortar method[END_REF]).

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 [START_REF] Park | Cohesive zone models: A critical review of traction-separation relationships across fracture surfaces[END_REF][START_REF] Zhang | Extrinsic cohesive modelling of dynamic fracture and microbranching instability in brittle materials[END_REF][START_REF] Collins-Craft | On the formulation and implementation of extrinsic cohesive zone models with contact[END_REF]; 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 symplectic CD-Lagrange scheme, with perfect frictionless unilateral impacts. Then we move to two new test cases, that could serve as benchmarks, involving perfect plasticity in the interface model: no elasticity is introduced to emphasize the non-smoothness feasibility of the approach.

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 𝑀 and 𝐾, submitted to an external force is driven by the second-order in time ordinary differential equation (ODE) of the form 𝑀 𝑈 + 𝐾𝑈 = 𝑓 ext , where 𝑈 (𝑡) is the displacement degreesof-freedom vector and 𝑓 ext (𝑡) 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 𝑛 and 𝑛 + 1/2 are related to the time steps 𝑡 𝑛 and 𝑡 𝑛+1/2 . Moreover, as a time-variational scheme, the order of the ODE is reduced by an order 1, and the time-discretized dynamic problem is

𝑀 (𝑉 𝑛+3/2 -𝑉 𝑛+1/2 ) + ℎ𝐾𝑈 𝑛+1 = ℎ𝑓 ext,𝑛+1 (1)
where ℎ is the time step (ℎ = Δ𝑡 assumed herein uniform), and 𝑉 is the velocity vector, such that 𝑈 𝑛+1 = 𝑈 𝑛 + ℎ𝑉 𝑛+1/2 is the explicit computation of the incremental displacement, provided that both 𝑈 𝑛 and 𝑉 𝑛+1/2 are known, and that 𝑈 𝑛+1 and 𝑉 𝑛+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 𝑀 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 𝐿 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 jump at contact interface is 𝑣 = 𝐿𝑉 , and the normal gap is 𝑔 = 𝑔 0 + 𝐿𝑈 where 𝑔 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 𝑟 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

𝑀 (𝑉 𝑛+3/2 -𝑉 𝑛+1/2 ) + ℎ𝐾𝑈 𝑛+1 = ℎ𝑓 ext,𝑛+1 + 𝐿 𝑇 𝑟 𝑛+3/2 (2)
or

𝑉 𝑛+3/2 = 𝑉 free + 𝑀 -1 𝐿 𝑇 𝑟 𝑛+3/2 (3)
where 𝑉 free = 𝑉 𝑛+1/2 + 𝑀 -1 ℎ(𝑓 ext,𝑛+1 -𝐾𝑈 𝑛+1 ) would be the velocity if no contact occurs, i.e. if 𝑟 𝑛+3/2 = 0. The reduced dynamics problem consists in expressing from (3) the relation between local unknown quantities at the contact interface 𝑣 𝑛+3/2 and 𝑟 𝑛+3/2 . This is done by pre-multiplying with 𝐿 to get

𝑣 𝑛+3/2 = 𝑣 free + 𝐻𝑟 𝑛+3/2 (4)
where 𝑣 free = 𝐿𝑉 free and 𝐻 = 𝐿𝑀 -1 𝐿 𝑇 is the so-called Delassus operator. When considering compatible finite element meshes between impacting bodies, or impact with a fixed rigid part, the operator 𝐿 is local to each impacting nodes. Due to the properties of the lumped mass matrix 𝑀; this leads to a Delassus operator which is diagonal, definite positive and spherical per node (spatial directions are equivalent, as in the mass matrix 𝑀 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 𝑔 and impulse 𝑟 , it leads to the KKT condition: 0 ≤ 𝑔 ⊥ 𝑟 ≥ 0. This fits within the thermodynamics framework, using a free energy potential Ψ(𝑔) = 𝐼 R + (𝑔) where 𝐼 R + is the indicator function of the convex positive cone R + (with a null value for a positive argument, and +∞ otherwise). Clearly, the non-smoothness arises from the non differential property of this potential, so that the state law involves a subdifferential operation as

-𝑟 ∈ 𝜕 𝑔 Ψ(𝑔) (5)
which is equivalent to the previous KKT condition. This condition can be replaced, thanks to Moreau viability lemma [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF][START_REF] Dubois | The Contact Dynamics method: A nonsmooth story[END_REF], with a relation between impulse and velocity jump, that can also be interpreted as a constitutive relation of the form: R (𝑣, 𝑟 ; 𝑔) = 0. This relation is also neither smooth nor univoque. It reads:

If 𝑔 > 0 then: 𝑟 = 0 (6) Else: 0 ≤ 𝑣 ⊥ 𝑟 ≥ 0 (7)
Finally the problem is (4) together with the behavior R (𝑣 𝑛+3/2 , 𝑟 𝑛+3/2 ; 𝑔 𝑛+1 ) = 0 and with diagonal Delassus operator, it provides an explicit solution as well. Once 𝑟 𝑛+3/2 is known, the dynamics (3) allows to get the full velocity vector 𝑉 𝑛+3/2 .

The classical Coulomb friction laws can also be taken easily into account [START_REF] Fekak | A new heterogeneous asynchronous explicit-implicit 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], with an underlying thermodynamics framework suited to this non-associated model thanks to the bi-potential method [START_REF] De Saxcé | The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms[END_REF][START_REF] Feng | The bi-potential method applied to the modeling of dynamic problems with friction[END_REF]. Indeed, there are other explicit formulations for contact problems, as the Nitsche approach in [START_REF] Chouly | Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems[END_REF]) but dedicated to deformable bodies.

A simple 0D 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 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 𝑐) and a purely plastic behavior (with index 𝑝) with the stress threshold 𝜎 𝑌 , as depicted on Figure 1. This association leads to the following assembly relationships: 𝑟 = 𝑟 𝑐 + 𝑟 𝑝 for the impulses and 𝑣 𝑐 = 𝑣 𝑝 = 𝑣 for the velocity 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 volumic elastic bodies, see for instance some dedicated plasticity models (de los [START_REF] De Los Ríos | Laminated materials with plastic interfaces: modeling and calculation[END_REF], elasto-plasticity [START_REF] Spada | Damage and plasticity at the interfaces in composite materials and structures[END_REF][START_REF] Brisotto | Simulation of interface behavior between FRP bars and concrete by an elastic-plastic theory via FEM[END_REF], visco-plasticity [START_REF] Monchiet | Interfacial models in viscoplastic composites materials[END_REF], 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 [START_REF] Lemaitre | Elasto-plastic Contact of Rough Surfaces[END_REF]. 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 𝜎 𝑝 and the displacement jump 𝑢 𝑝 (the equivalent for an interface behavior to the strain for a volumic plasticity model). The non-smooth dissipation potential, for the 0D case, is defined as

Φ ★ (𝜎 𝑝 ) = 𝐼 R -(𝑓 (𝜎 𝑝 )) with 𝑓 (𝜎 𝑝 ) = |𝜎 𝑝 | -𝜎 𝑌
where 𝑓 (𝜎 𝑝 ) is the yield function, 𝜎 𝑌 is the yield stress and 𝐼 R -is the indicator function of the negative cone R -(with a null value for a negative argument, and +∞ otherwise). Clearly, the non-smoothness again arises from the non differential property of this potential [START_REF] Moreau | On Unilateral Constraints, Friction and Plasticity[END_REF], so that the evolution law involves also a subdifferential operation as

𝑢 𝑝 ∈ 𝜕 𝜎 𝑝 Φ ★ leading to If 𝑓 (𝜎 𝑝 ) < 0 then: 𝑢 𝑝 = 0 (8) Else: 𝑢 𝑝 = 𝜆 sgn(𝜎 𝑝 ) with 𝜆 ≥ 0 (9)
Expressed with velocity jump 𝑣 𝑝 = 𝑢 𝑝 and impulse 𝑟 𝑝 = -ℎ𝜎 𝑝 (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 𝑟 𝑌 = ℎ𝜎 𝑌 . Some corresponding stress vs velocity jump answers are illustrated on Figure 3.

v p 0 hσ p r Y -r Y Figure 2:
The considered constitutive behavior, as an impulse/velocity jump relationship When the gap is strictly positive, the contact is disabled, so 𝑟 𝑐 = 0, 𝑟 = 𝑟 𝑝 , and the reduced dynamics reads

u p 0 σ p σ Y u p 0 σ p σ Y -σ Y
𝑣 = 𝑣 free + 𝐻𝑟 𝑝 (10)
With the diagonal and positive properties of 𝐻 , 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.

v p 0 r p r Y -r Y r Figure 4:
The solution as the intersection of the behavior (plain line) and the reduced dynamics (dotted line)

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 𝑟 = 𝑟 𝑐 + 𝑟 𝑝 with the reduced dynamics

𝑣 = 𝑣 free + 𝐻𝑟 (11)
Indeed, consider the cases separation:

• if 𝑣 > 0, then 𝑟 𝑐 = 0, so 𝑟 = 𝑟 𝑝 = -𝑟 𝑌 is computed as before. One therefore gets 𝑣 = 𝑣 free -𝐻𝑟 𝑌 := ṽ that should be positive in this case; • in the counter case, 𝑣 = 0, and the reduced dynamics reads 0 = 𝑣 free + 𝐻𝑟 = (𝑣 free -𝐻𝑟 𝑌 ) + 𝐻 (𝑟 𝑐 + 𝑟 𝑝 + 𝑟 𝑌 ), so 0 = ṽ + 𝐻 (𝑟 𝑐 + 𝑟 𝑌 -ℎ𝜎 𝑝 ) for which we should have 𝑟 𝑐 ≥ 0 and 0 ≤ 𝑟 𝑌 -ℎ𝜎 𝑝 ≤ 2𝑟 𝑌 . Now if ṽ > 0 this is inconsistent and corresponds to the previous case. Therefore ṽ ≤ 0 and

𝑟 = 𝑟 𝑐 -ℎ𝜎 𝑝 = -𝐻 -1 𝑣 free := r
This cases separation can therefore be summarized as

• if ṽ > 0 (which is equivalent to r < -𝑟 𝑌 ), then 𝑟 = -𝑟 𝑌 • else 𝑟 = r
This can be simplified also, as detailed in Algorithm 1.

Input: displacement 𝑈 𝑛 , velocity 𝑉 𝑛+1/2 Output: displacement 𝑈 𝑛+1 , velocity 𝑉 𝑛+3/2 1: Displacement 𝑈 𝑛+1 ← 𝑈 𝑛 + ℎ𝑉 𝑛+1/2 ⊲ Explicit configuration 2: 𝑔 𝑛+1 ← 𝑔 0 + 𝐿𝑈 𝑛+1 3: 𝑉 free ← 𝑉 𝑛+1/2 + 𝑀 -1 ℎ(𝑓 ext,𝑛+1 -𝐾𝑈 𝑛+1 ) ⊲ Free velocity 4: 𝑣 free ← 𝐿𝑉 free 5: r ← -𝐻 -1 𝑣 free 6: 𝑟 𝑌 ← ℎ𝜎 𝑌 7: if 𝑔 𝑛+1 > 0 then 8: if r < -𝑟 𝑌 then 𝑟 𝑛+3/2 ← -𝑟 𝑌 9: else if r > 𝑟 𝑌 then 𝑟 𝑛+3/2 ← 𝑟 𝑌 10: else 𝑟 𝑛+3/2 ← r 11: end if 12: else 𝑟 𝑛+3/2 ← max( r, -𝑟 𝑌 ) 13: end if 14: Velocity 𝑉 𝑛+3/2 ← 𝑉 free + 𝑀 -1 𝐿 𝑇 𝑟 𝑛+3/2
⊲ 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 on Figure 5, with the non-dimensional parameters: 𝑔 0 = 0, 𝑚 = 𝑘 = 2𝜎 𝑌 = 1 (the dimensional 𝑚 would be a mass per unit surface, the dimensional 𝑘 a stiffness homogeneous to a Young modulus). The characteristic time for the oscillations in the linear regime is 𝜏 = 2𝜋 √ 𝑘/𝑚 = 2𝜋. 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 ℎ = Δ𝑡 should be smaller than 𝜏 for accuracy criteria; here ℎ = 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 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 𝜂 =min 𝑡 𝑔(𝑡)/max 𝑡 |𝑔(𝑡)|. On the present test case, it takes the value 𝜂 = 0.43 %.

Another phenomenological association of contact and perfect plasticity

A serial association of the two previous elementary blocks for the 0D 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 volumic bodies, but within the surfacic interface. Some examples of this association type may be found, for elasto-plasticity model (Liu et al. 2001;[START_REF] Ghaednia | A review of elastic-plastic contact mechanics[END_REF], for elasto-plastic contact at asperity scale [START_REF] Kogut | A finite element based elastic-plastic model for the contact of rough surfaces[END_REF], or at a larger scale (Wriggers et al. 1993), etc. This association leads to 𝑟 = 𝑟 𝑐 = 𝑟 𝑝 , 𝑣 = 𝑣 𝑐 + 𝑣 𝑝 and 𝑢 = 𝑢 𝑐 + 𝑢 𝑝 . The behavior of each block is the same as previously. The reduced dynamics still reads 𝑣 = 𝑣 free + 𝐻𝑟 . The same method as previously with case separation now leads to the following: if 𝑢 > 𝑢 𝑝 (i.e. 𝑢 𝑐 > 0, the contact is inactive) then 𝑟 𝑐 = 0 = 𝑟 𝑝 = 𝑟 (and 𝑣 𝑝 = 0, 𝑣 𝑐 = 𝑣 = 𝑣 free ); otherwise,

• if 𝑣 free > 0 (contact separation) then: 𝑟 𝑐 = 0 = 𝑟 𝑝 = 𝑟 (and 𝑣 𝑝 = 0, 𝑣 𝑐 = 𝑣 = 𝑣 free )

• otherwise (active contact) if 𝑣 free > -𝐻𝑟 𝑌 (no plasticity development) then: 𝑟 𝑐 = 𝑟 = 𝑟 𝑝 = -𝐻 -1 𝑣 free (indeed less than 𝑟 𝑌 ) and 𝑣 𝑐 = 0, 𝑣 𝑝 = 0 so 𝑣 = 0 (rigid interface) else 𝑣 free ≤ -𝐻𝑟 𝑌 (compressive plastification) so that: 𝑟 𝑝 = 𝑟 𝑌 (and 𝑟 𝑐 = 𝑟 = 𝑟 𝑌 as well), 𝑣 = 𝑣 𝑝 = 𝑣 free + 𝐻𝑟 𝑌 ≤ 0 (and 𝑣 𝑐 = 0) The numerical solution procedure is now described on algorithm 2 (one can note the need of the local plastic displacement 𝑢 𝑝 as a state internal variable).

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

𝑟 𝑛+3/2 ← max 0 , min(-𝐻 -1 𝑣 free , 𝑟 𝑌 ) 11: 𝑣 𝑝,𝑛+3/2 ← min(0 , 𝑣 free + 𝐻𝑟 𝑌 ) 12: end if 13: Velocity 𝑉 𝑛+3/2 ← 𝑉 free + 𝑀 -1 𝐿 𝑇 𝑟 𝑛+3/2
⊲ Matrix-free dynamics

Algorithm 2: Time step increment for the perfectly plastic contact interface for a sufficient value, the plastic flow is observed.

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 TM ". 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 [START_REF] Raous | Restorable adhesion and friction[END_REF]. 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.

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 𝑓 𝑌 (or impulse 𝑟 𝑌 = ℎ𝑓 𝑌 ) 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 𝑣 = 𝑣 free + 𝐻𝑟 . As before, with the diagonal and positivity properties of 𝐻 , the solution is the intersection of the behavior and the reduced dynamics as illustrated on Figure 10, when the gap is null. When the gap is strictly positive, the fracture has occurred, so 𝑟 = 0. When the gap is not strictly positive, three different cases may happen: • either r ≤ -𝑟 𝑌 , then the solution is 𝑟 = 0 (and 𝑣 ≥ 0, the fracture occurs); • or r > 0, then the solution is 𝑟 = r (and 𝑣 = 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 𝑟 = 0, or the previous status was active contact or glued, and the status is now glued and 𝑟 = r . This behavior does nevertheless not take into account that the yield stress 𝑓 𝑌 may evolve (and is therefore an internal variable). Indeed, the considered modified model is: when a rupture occurs, 𝑓 𝑌 is nullified; when the impulse is positive (compression), 𝑓 𝑌 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 𝑓 𝑌 internal variable, so that the 𝑓 𝑌 = 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 𝑓 𝑌 = 𝑟 𝑌 = 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 𝑓 𝑌 or an impulse distribution 𝑟 𝑌 are used, related with 𝑟 𝑌 = ℎ𝑓 𝑌 . Therefore, ℎ being a discretization parameter, 𝑓 𝑌 and 𝑟 𝑌 cannot both be intrinsic material parameters. An answer could a priori be obtained by looking at convergence of the solution when the time step ℎ goes to zero. If a 0D problem is considered, loaded with an external force 𝑓 ext (𝑡), and a material parameter 𝑟 𝑌 , the impulse solution looks like 𝑟 (𝑡) = -ℎ𝑓 ext (𝑡), therefore 𝑟 tends toward 0 (and 𝑟 /ℎ has a finite limit), so whatever the load is, after a certain refinement in time discretization, one cannot have 𝑟 ≤ -𝑟 𝑌 , so there is an inconsistency. Now consider a 0D problem, loaded with an external impulse 𝑟 ext (𝑡) (so that the equivalent force is 𝑓 (𝑡) = 𝑟 ext (𝑡)/ℎ). If the material parameter is 𝑓 𝑌 , since the solution looks like 𝑟 (𝑡) = 𝑟 ext (𝑡), after a certain refinement, and whatever the load is, one will always have 𝑟 > -𝑟 𝑌 = ℎ𝑓 𝑌 , 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: 𝑟 < -𝑟 𝑌 or 𝑟 < -ℎ𝑓 𝑌 . 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 𝑅(𝑡) = ∫ 𝑡 𝑡 -𝜏 𝑌 𝑓 (𝑡)𝑑𝑡 may be defined, where 𝑓 is the interface normal stress, which is equivalent to the more general expression 𝑅(𝑡) = 𝑖,𝑡 𝑖 ∈ [𝑡 -𝜏 𝑌 ,𝑡 ] 𝑟 (𝑡 𝑖 ) 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 𝑅 < -𝑅 𝑌 , where 𝑅 𝑌 is a material parameter homogeneous to an impulse distribution, together with a second intrinsic parameter, namely a characteristic time 𝜏 𝑌 . As mentioned previously, one may consider a variable-in-time 𝑅 𝑌 . 𝑅 𝑌 is nullified when fracture occurs, and can be recovered by normal compression, i.e. when 𝑟 > 0. A consistent criteria could then read: 𝑅 𝑌 (𝑡) = max(𝑅 𝑌 (𝑡); 𝑅(𝑡)). And finally a maximum recovery could be added with a third material parameter, 𝑅 𝑌 max leading to 𝑅 𝑌 (𝑡) = min(𝑅 𝑌 max ; max(𝑅 𝑌 (𝑡); 𝑅(𝑡))). 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 [START_REF] Zhang | Model smoothing method of contact-impact dynamics in flexible multibody systems[END_REF].

A main implementation difference is due to the non-locality in time: the impulse field 𝑟 should be stored during a previous duration of 𝜏 𝑌 , 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: 𝛼, that measures the efficiency of the adhesion recovery: as an initial condition, 𝛼 = 1 (maximum efficiency). At each fracture event at time 𝑡, it is then decreased, for instance with the following evolution law: 𝛼 (𝑡 + ) = 𝛼 (𝑡 -)𝑒 -1/𝑛 𝑏 , where 𝑛 𝑏 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 (𝑟 = 0, so 𝑣 𝑇 𝑟 = 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 𝑚 = 1, but without stiffness 𝑘 = 0 and with a brittle interface behavior. As a load, some shocks are produced, i.e. an external impulse 𝑟 ext (𝑡) is prescribed, Figure 11. Concerning the time discretization, 160 time steps of uniform duration ℎ 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 𝑅 𝑌 max = ∞ (no adhesion limit), 𝑛 𝑏 = ∞ (no efficiency decrease), 𝜏 𝑌 = 0 (no delay effect). For these values, a cyclic response is produced, Figure 12. A first impulse 𝑟 (𝑡 = ℎ) > 0 is produced due to the initial condition that engenders a compressive shock and increases the yield impulse 𝑅 𝑌 ; then a negative impulse happens opposed to the external first positive shock 𝑟 ext which is not large enough to produce a fracture; the second external impulse is large enough for a fracture to happen (the impulse 𝑟 is therefore null, as the yield impulse 𝑅 𝑌 ), 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 𝑟 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 𝑛 𝑏 = 10. The solution is no more cyclic, Figure 13. Indeed the yield impulse reconstruction 𝑅 𝑌 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 𝜏 𝑌 = 0.08 (actually 64 time steps). Though the solution differs from test 1 for the yield 𝑅 𝑌 , Figure 14, due to the time delay that is a moving window in time, the associated material behavior is conservative, are allows to recover the same solution for the other quantities. To check time-convergence, an additional test is used, based on the previous test 3, but with an additional smooth loading 𝑓 ext (𝑡) = -𝐹 𝑑 sin( 45 𝜋𝑡/𝑇 𝑐 ) with a non-dimensional amplitude 𝐹 𝑑 = 2, 𝑇 𝑐 = 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 (ℎ ref ≈ 7.63 10 -8 ). Since for instance the velocity is discontinuous, a Hausdorff distance [START_REF] Moreau | Approximation en graphe d'une évolution discontinue[END_REF]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 𝑡/𝑇 , 𝑇 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.

3D structures and surfacic interfaces

To start with, the interface behavior is considered with a normal 𝑓 𝑁 and a tangential force distribution (or stress) 𝑓 𝑇 , together with a yield stress 𝜎 𝑌 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 equivalent stress is defined as

𝜎 (𝑓 𝑁 , 𝑓 𝑇 ) = √ < 𝑓 𝑁 > 2 + +∥𝑓 𝑇 ∥ 2 /𝜆 2 (12)
so only the tension normal part (hence the positive part < 𝑓 𝑁 > + is involved) acts on the fracture, and where 𝜆 is a weight for the influence of the tangential part, Figure 16. When 𝜎 (𝑓 𝑁 , 𝑓

𝑇

) < 𝜎 𝑌 no fracture occurs and the interface status is 'glued', i.e. the velocity jump is 𝑣 = 0; when 𝜎 (𝑓 𝑁 , 𝑓

𝑇

) > 𝜎 𝑌 a fracture has occurred and the interface status is 'contact', i.e. unilateral contact with Coulomb tangential friction, with a friction coefficient 𝜇, that may be either inactive (𝑣 𝑁 > 0, 𝑓 𝑁 = 0, 𝑓

𝑇 = 0) or sliding (𝑣 𝑇 ≠ 0, 𝑓 𝑁 > 0, ∥ 𝑓 𝑇 ∥ = 𝜇 𝑓 𝑁 ) or sticking (𝑣 𝑇 = 0, 𝑓 𝑁 > 0, ∥𝑓 𝑇 ∥ ≤ 𝜇 𝑓 𝑁 )
. 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 𝜎 𝑌 = 0, the threshold collapses to the positive semi-axis 𝑓 𝑁 < 0. When 𝜆 = ∞ the tangential direction has no influence on the fracture, the normal behavior is identical as in the 0D test case. When 𝑅 𝑌 max = 0, the behavior reduces to a classical contact with Coulomb friction. 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.

f N σ Y λσ Y σ = σ Y σ < σ Y σ > σ Y tension compression || f T ||
The resolution algorithm 4 is then built on the same framework as previously. In particular, the velocity-impulse formulation is used, with 𝑟 𝑌 = ℎ𝜎 𝑌 , 𝑟 𝑁 = -ℎ𝑓 𝑁 , 𝑟 𝑇 = -ℎ𝑓 𝑇 , 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 𝑘 is 𝑆 (𝑘 ) . To simplify the notations, their assembling on the full set of dofs will be denoted with the diagonal entry matrix 𝑆, 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

𝑉 𝑛+3/2 = 𝑉 free + 𝑀 -1 𝑆𝐿 𝑇 𝑟 𝑛+3/2 (13) = 𝑉 free + 𝑀 -1 𝑆𝐿 𝑇 𝑁 𝑟 𝑁 ,𝑛+3/2 + 𝑀 -1 𝑆𝐿 𝑇 𝑇 𝑟 𝑇 ,𝑛+3/2 (14)
where the mapping operators 𝐿 𝑇 𝑁 and 𝐿 𝑇 𝑇 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 𝑀, still valid in 3D, and similar properties for 𝑆, one has 𝐿 𝑁 𝑀 -1 𝑆𝐿 𝑇 𝑇 = 0 so that the reduced dynamics reads firstly: 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 [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF][START_REF] Tkachuk | Hybrid-mixed discretization of elastodynamic contact problems using consistent singular mass matrices[END_REF][START_REF] Dabaghi | A weighted finite element mass redistribution method for dynamic contact problems[END_REF][START_REF] Di Stasio | An explicit timeintegrator with singular mass for non-smooth dynamics[END_REF]. There is also a physical source of energy loss, due to the friction between fractured surfaces, driven by the friction coefficient 𝜇.

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 𝑅 1 = 1, the thickness 𝐿 1 = 0.5 and the length 𝐿 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 𝐸 = 1, a Poisson coefficient 𝜈 = 0.3 and a density of 𝜌 = 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 𝑅 𝑌 (𝑡 = 0) = 10 -4 . The studied time interval is 𝑇 = 16, with a time step ℎ = 𝑇 /8000 = 0.002, while the estimated critical time step is 𝑡 𝑐 ≈ 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 𝑔 0 = 0 and an initial velocity 𝑉 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. 𝑔 0 = -10 -10 , so that the initial state is active contact.

If 𝜆 = ∞ (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 𝜆 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 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 𝑐 (i.e. for 𝑐 3 times more finite elements) and a time step divided by 𝑐 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 𝑐 ∈ {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: where < • > -stands for the negative part. For the different refinement ratios 𝑐, 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 𝑐.

Conclusions and perspectives

The previously developed symplectic and 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 [START_REF] Collins-Craft | On the formulation and implementation of extrinsic cohesive zone models with contact[END_REF]. To avoid spatial localization a delay-effect model [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF][START_REF] Allix | Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading[END_REF][START_REF] Allix | A delay damage mesomodel of laminates under dynamic loading: basic aspects and identification issues[END_REF][START_REF] Desmorat | Delay-active damage versus nonlocal enhancement for anisotropic damage dynamics computations with alternated loading[END_REF]) will be more suited since non-locality in space [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Lasry | Localization limiters in transient problems[END_REF][START_REF] Kamensky | A Review of Nonlocality in Computational Contact Mechanics[END_REF] would prevent an explicit matrix-free approach. Generalization to asynchronous explicit schemes [START_REF] Gravouil | Heterogeneous asynchronous time integrators for computational structural dynamics[END_REF][START_REF] Niu | An asynchronous variational integrator for the phase field approach to dynamic fracture[END_REF]) 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.

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

Figure 1 :

 1 Figure 1: Sketch of the phenomenological interface behavior

Figure 3 :

 3 Figure 3: The considered constitutive behavior, as a stress/displacement jump relationship

Figure 5 :

 5 Figure 5: Sketch of the test case

Figure 6 :

 6 Figure 6: The numerical solution for the 0D test case with the perfectly plastic interface

Figure 7 :

 7 Figure 7: Sketch of the test case for the plastic contact

Figure 8 :

 8 Figure 8: The numerical solution for the 0D test case with the contact with perfect plasticity interface

Figure 9 :Figure 10 :

 910 Figure 9: The considered brittle constitutive behavior, as an impulse/displacement jump relationship (the dotted line denotes a brutal jump in the status)

Figure 11 :

 11 Figure 11: The non-smooth loading for the 0D test cases

Figure 12 :Figure 13 :Figure 14 :

 121314 Figure 12: The numerical solution for the first 0D brittle test case

Figure 15 :

 15 Figure15: The loading with forces and impulses, the numerical convergence of the solution with respect to time discretization and the reference solution for velocity

Figure 16 :

 16 Figure 16: Yield function for 2D fracture criteria

Figure 17 :

 17 Figure 17: Geometry of the 3D test case, with the vertical shear stress on the face end

Figure 18 :

 18 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 (𝑐 = 1); (e) test case 2: for space-time convergence (𝑐 = 3); (f) test case 2: for space-time convergence (𝑐 = 4)

  𝜂 =max 𝑡 (-mean 𝑥 𝑔(𝑥, 𝑡))max 𝑡 mean 𝑥 |𝑔(𝑥, 𝑡)| or 𝜂 ′ = mean 𝑥,𝑡 < 𝑔(𝑥, 𝑡) > - mean 𝑥,𝑡 |𝑔(𝑥, 𝑡)| (15)

  Input: displacement 𝑈 𝑛 and 𝑢 𝑝,𝑛 , velocity 𝑉 𝑛+1/2 and 𝑣 𝑝,𝑛+1/2 Output: displacement 𝑈 𝑛+1 and 𝑢 𝑝,𝑛+1 , velocity 𝑉 𝑛+3/2 and 𝑣 𝑝,𝑛+3/2 1: Displacement 𝑈 𝑛+1 ← 𝑈 𝑛 + ℎ𝑉 𝑛+1/2 and 𝑢 𝑝,𝑛+1 ← 𝑢 𝑝,𝑛 + ℎ𝑣 𝑝,𝑛+1/2 ⊲ Explicit configuration 2: 𝑔 𝑛+1 ← 𝑔 0 + 𝐿𝑈 𝑛+1 3: 𝑉 free ← 𝑉 𝑛+1/2 + 𝑀 -1 ℎ(𝑓 ext,𝑛+1 -𝐾𝑈 𝑛+1 ) ⊲ Free velocity 4: 𝑣 free ← 𝐿𝑉 free 5: 𝑟 𝑌 ← ℎ𝜎 𝑌 6: if 𝑔 𝑛+1 > 𝑢 𝑝,𝑛+1 then

		⊲ Inactive contact
	7:	𝑟 𝑛+3/2 ← 0
	8:	𝑣 𝑝,𝑛+3/2 ← 0
	9: else
	10:

Table 1 :

 1 𝑈 𝑛 and 𝑉 𝑛+1/2 as well as interface status 𝑅 𝑌 ,𝑛+1/2 and 𝛼 𝑛+1/2 Output: 𝑈 𝑛+1 and 𝑉 𝑛+3/2 as well as new interface status 𝑅 𝑌 ,𝑛+3/2 and 𝛼 𝑛+3/2 1: Displacement𝑈 𝑛+1 ← 𝑈 𝑛 + ℎ𝑉 𝑛+1/2 ⊲ Explicit configuration 2: 𝑔 𝑛+1 ← 𝑔 0 + 𝐿𝑈 𝑛+1 3: 𝑉 free ← 𝑉 𝑛+1/2 + 𝑀 -1 ℎ(𝑓 ext,𝑛+1 -𝐾𝑈 𝑛+1 ) 𝑖,𝑡 𝑖 ∈ [𝑡 𝑛+3/2 -𝜏 𝑌 ,𝑡 𝑛+1/2 ] 𝑟 𝑖 12:R𝑌 ← min(𝑅 𝑌 max ; max(𝑅 𝑌 ,𝑛+1/2 ; 𝛼 𝑛+1/2 R)) Interface non-dimensional parameters for the different test cases

	⊲ Free velocity

Input:

  𝑣 𝑁 = 𝑣 𝑁 free + 𝐻 𝑁 𝑟 𝑁 , where 𝐻 𝑁 = 𝐿 𝑁 𝑀 -1 𝑆𝐿 𝑇 𝑁 , 𝑣 𝑁 free = 𝐿 𝑁 𝑉 free , and secondly: 𝑣 𝑇 = 𝑣 𝑇 free + 𝐻 𝑇 𝑟 𝑇 , where 𝐻 𝑇 = 𝐿 𝑇 𝑀 -1 𝑆𝐿 𝑇 𝑇 , 𝑣 𝑇 free = 𝐿 𝑇 𝑉 free . The Delassus operators 𝐻 𝑁 and 𝐻 𝑇 are still definite positive, diagonal and spherical per node entry. 𝑈 𝑛 and 𝑉 𝑛+1/2 as well as interface status fields 𝑅 𝑌 ,𝑛+1/2 and 𝛼 𝑛+1/2 Output: 𝑈 𝑛+1 and 𝑉 𝑛+3/2 as well as new interface status fields 𝑅 𝑌 ,𝑛+3/2 and 𝛼 𝑛+3/2 1: Displacements 𝑈 𝑛+1 ← 𝑈 𝑛 + ℎ𝑉 𝑛+1/2 ⊲ Explicit configuration 2: 𝑔 𝑛+1 ← 𝑔 0 + 𝐿𝑈 𝑛+1 3: 𝑉 free ← 𝑉 𝑛+1/2 + 𝑀 -1 ℎ(𝑓 ext,𝑛+1 -𝐾𝑈 𝑛+1 ) ⊲ Free velocity 4: for each interface point 𝑘 do 𝑖,𝑡 𝑖 ∈ [𝑡 𝑛+3/2 -𝜏 𝑌 ,𝑡 𝑛+1/2 ] 𝑟 𝑁 ,𝑖 16: R𝑇 ← r𝑇 + 𝑖,𝑡 𝑖 ∈ [𝑡 𝑛+3/2 -𝜏 𝑌 ,𝑡 𝑛+1/2 ] 𝑟 𝑇 ,𝑖 17:R𝑌 ← min(𝑅 𝑌 max ; max(𝑅 𝑌 ,𝑛+1/2 ; 𝛼 𝑛+1/2 R𝑁 )) Velocity 𝑉 𝑛+3/2 ← 𝑉 free + 𝑀 -1 𝑆𝐿 𝑇 𝑟 𝑛+3/2

	5:	𝑣 𝑁 free ← 𝐿 (𝑘 ) 𝑁 𝑉 free	⊲ With the local to node 𝑘 operator
	6:	𝑣 𝑇 free ← 𝐿 (𝑘 ) 𝑇 𝑉 free
	7:	if 𝑔 𝑘 𝑛+1 > 0 then		⊲ Inactive interaction
	8:	𝛼 𝑛+3/2 ← 𝛼 𝑛+1/2
	9:	𝑅 𝑌 ,𝑛+3/2 ← 0
	10:	𝑟 𝑁 ,𝑛+3/2 ← 0
	11:	𝑟 𝑇 ,𝑛+3/2 ← 0
	12:	else
	13: 14:	r𝑁 ← -𝐻 -1 𝑁 𝑣 𝑁 free r𝑇 ← -𝐻 -1 𝑇 𝑣 𝑇 free
	15: R𝑁 ← r𝑁 + 18: if 𝜎 (-R𝑁 , -R𝑇 ) < R𝑌 then		⊲ Glued interaction
	19:	𝛼 𝑛+3/2 ← 𝛼 𝑛+1/2
	20:	𝑅 𝑌 ,𝑛+3/2 ← R𝑌
	21:	𝑟 𝑁 ,𝑛+3/2 ← r𝑁
	22:	𝑟 𝑇 ,𝑛+3/2 ← r𝑇
	23:	else
	24:	if 𝑅 𝑌 ,𝑛+1/2 > 0 then		⊲ Brittle fracture occurs
	25:	𝛼 𝑛+3/2 ← 𝛼 𝑛+1/2 𝑒 -1/𝑛 𝑏		⊲ Adhesion efficiency decrease
	26:	else
	27:	𝛼 𝑛+3/2 ← 𝛼 𝑛+1/2
	28:	end if
	29:	𝑅 𝑌 ,𝑛+3/2 ← 0
	30:	if r𝑁 > 0 then		⊲ Active contact
	31:	𝑟 𝑁 ,𝑛+3/2 ← r𝑁
	32:	if ∥ r𝑇 ∥ > 𝜇𝑟 𝑁 ,𝑛+3/2 then	⊲ Sliding
	33:	𝑟 𝑇 ,𝑛+3/2 ← 𝜇𝑟 𝑁 ,𝑛+3/2	r𝑇 ∥ r𝑇 ∥
	34:	else		⊲ Sticking
	35:	𝑟 𝑇 ,𝑛+3/2 ← r𝑇
	36:	end if
	37:	else		⊲ Released contact
	38:	𝑟 𝑁 ,𝑛+3/2 ← 0
	39:	𝑟 𝑇 ,𝑛+3/2 ← 0
	40:	end if
	41:	end if
	42:	end if
	43:	Assemble 𝑟 𝑁 ,𝑛+3/2 and 𝑟 𝑇 ,𝑛+3/2 in the right-hand-side 𝑟 𝑛+3/2
	44: end for
	45:

Input: ⊲ Matrix-free dynamics Algorithm 4: Time step increment for the 2D brittle fracture interface

Table 2 :

 2 Error indicators vs the refinement ratio

	𝑐	1	3	4
	𝜂 0.016 % 0.012 % 0.009 %
	𝜂 ′ 0.036 % 0.023 % 0.019 %

  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/


Antonio Castro
email: antonio.castro@eic.cefet-rj.br

Cefet / Rj

Cassio Souza
email: cassio.souza@eic.cefet-rj.br

Jorge Rodrigues
email: jorge.rodrigues@eic.cefet-rj.br

Fabio Porto

Esther Pacitti
email: esther.pacitti@lirmm.fr

Rafaelli Coutinho
email: rafaelli.coutinho@cefet-rj.br

Eduardo Ogasawara
email: eogasawara@ieee.org

GSTSM Package: Finding Frequent Sequences in Constrained Space and Time

Keywords: Data Mining, Spatial-Temporal, Time Series, Sequential Mining

Spatial time-stamped sequences have information about time and space where events occur. Mining such sequences can bring important insights. However, not all sequences are frequent over an entire dataset. Some are only common in subsets of time and space. This article explains the first tool for mining these sequences in constrained space and time: the GSTSM R package. It allows users to search for spatio-temporal patterns that are not frequent in the entire database, but are dense in restricted time-space intervals. Thus, making it possible to find non-trivial patterns that would not be found using common data mining tools.

Introduction

Data mining tools have been used to find interesting patterns in different areas of knowledge in various problems [START_REF] Alatrista-Salas | A knowledge discovery process for spatiotemporal data: Application to river water quality monitoring[END_REF]. The sequence mining knowledge area is a specialization of data mining, focused on finding sequences or series of events in datasets [START_REF] Parthasarathy | Incremental and interactive sequence mining[END_REF][START_REF] Zaki | Sequence mining in categorical domains: incorporating constraints[END_REF]. Such sequences may form patterns, a set of frequent attributes that appear persistently among the dataset. It means that its frequency exceeds a user-defined minimum threshold [START_REF] Aydin | Spatiotemporal event sequence mining from evolving regions[END_REF].

Several types of events involve both temporal and spatial data, such as financial to understand sales patterns over time and space [START_REF] Alatrista-Salas | A knowledge discovery process for spatiotemporal data: Application to river water quality monitoring[END_REF], and hydrological data for river water quality monitoring in different points over time [START_REF] Alatrista-Salas | A knowledge discovery process for spatiotemporal data: Application to river water quality monitoring[END_REF][START_REF] Alatrista-Salas | Spatio-sequential patterns mining: Beyond the boundaries[END_REF]. They correspond to Time-stamped Sequence (TS) events distributed in space [START_REF] Mooney | Sequential pattern mining -Approaches and algorithms[END_REF]. Mining sequences related to space and time enables to find knowledge related to phenomena that involve both spatial and temporal components, trying to find all sequences of significant, useful, interesting, and non-trivial events [START_REF] Aydin | Spatiotemporal event sequence mining from evolving regions[END_REF][START_REF] Sunitha | Mining frequent patterns from spatiotemporal data sets: A survey[END_REF][START_REF] Alatrista-Salas | A knowledge discovery process for spatiotemporal data: Application to river water quality monitoring[END_REF]. However, spatio-temporal sequential patterns may have low support if considered the entire dataset, but they can be frequent if considered only a period and region [START_REF] Huang | A framework for mining sequential patterns from spatio-temporal event data sets[END_REF]. The Generalized Spatial-Time Sequence Miner (GSTSM) package can find these patterns, being able to efficiently discover the region and period where they occur. This way, GSTSM would be the right tool to find time-localized patterns. This work describes the process, structure, and usage of the GSTSM package using a synthetic (but still complete) example. However, we also provide a glimpse of applicability in a real-world dataset. GSTSM was able to found sequential patterns in seismic data. They correspond to seismic horizons, which are important elements in the application domain.

preprint 2 Related Work

There are different methods for spatio-temporal data mining. Some use only data mining, searching for frequent patterns, considering only time [START_REF] Li | Prediction of human activity by discovering temporal sequence patterns[END_REF]. Others combine techniques by seeking in time and then grouping in space [START_REF] Flamand | Mining local climate data to assess spatiotemporal dengue fever epidemic patterns in French Guiana[END_REF]. Furthermore, there is a diversity in how constraints are handled. Some use global support, a value that is valid for the entire dataset [START_REF] Alatrista-Salas | Spatio-sequential patterns mining: Beyond the boundaries[END_REF]. Others consider local support, using predefined windows of time and space [START_REF] Koseoglu | ST Sequence Miner: visualization and mining of spatiotemporal event sequences[END_REF]. This work differs by seeking frequent sequences in time that occur in spatial groups. Instead of using predefined constraints for time and space, three density parameters are established: a minimum frequency to be achieved within the period, a maximum distance that a position can be from any other in the group, and a lower limit of distinct positions in the group. Thus, the formalization presented in this work can find different sizes of sequences, time intervals, and spatial regions where a sequence is frequent, based on the concepts of RG, KRG, and SRG introduced in [START_REF] Castro | Generalizac ¸ão de Minerac ¸ão de Sequências Restritas no Espac ¸o e no Tempo[END_REF].

As far as the conducted research has reached, the only work with a similar approach found in the literature is proposed by [START_REF] Campisano | Discovering tight space-time sequences[END_REF], which considers one-dimensional space. The present work is a generalization that presents a formalization considering space in its three-dimensional form.

Demonstration overview

GSTSM is a package which provides polymorphic functions that let the user extends its functionalities, as it is based on R [14] language S3 classes. The source code can be found at GitHub [START_REF] Castro | GSTSM Package[END_REF]. The package has a main class named GSTSM, that needs the parameters D, P, γ, β, and σ to instantiate an object, explained as follows:

• D and P represents the TS dataset with their respective positions. Each TS must be associated with one position. It means that the number of timestamped sequences (columns) in D must be equal to the number of positions (rows) in P. • for the user defined thresholds values in the range ]0, 1] for γ, values starting from 2 for β, and integer values starting from 1 for σ.

A GSTSM object is an instance of S3 Class built as a list with all this information. Furthermore, it generates an adjacency matrix that informs each position which other positions are at a maximum distance of σ.

The GSTSM package has the mine() method that implements the entire process of finding frequent sequences. It receives as input a GSTSM object and provides as output a list of the SRGs of all sizes found. The user does not need to call any other method to get the results. This method calls and passes all the necessary parameters to make the entire process transparent to the user.

The other methods used in each process step are polymorphic and can be extended by the user. It gives the user the ability to try its implementation. These are described as follows:

• find() has two input parameters: a GSTSM object and a set of candidate sequences of size k. It provides as output the KRGs for each candidate. • merge() has also two input parameters: a GSTSM object and a set of candidate sequences of size k containing information about the KRG of each one. The method returns the SRGs with the candidate sequences of size k. • generate candidates() has two input parameters: a GSTSM object and a set of SRGs of size k. There are no SRGs to pass to generate candidates of size one, a NULL value can be used. The method provides the candidate sequences of size k + 1.

An illustrative example shows the use of the GSTSM package functions. To start, the first action is installing and loading GSTSM package and then setting all the inputs for the package: D, P, γ, β, and σ. For D, we use a simple dataset. For P, positions in a row are used, with one unit distance. Each position is associated with a time series, such as p 1 to t 1 and p 2 to t 2 . The values for the user-defined thresholds are: γ = 0.8, β = 2, and σ = 1. After setting the input parameters, we can instantiate the GSTSM object and execute the mine() method. Listing 1 shows the code using the R command line.

Conclusion

GSTSM is the first tool for mining sequences in spatial time-stamped sequences datasets able to discover constrained patterns in time and space with all three dimensions. The package discovers patterns that may not be frequent over an entire dataset but are grouped in space and frequent in a time interval. It would not be easy to find these patterns without this tool. The results can differ from conventional data mining tools and give different about data behavior. The patterns are groups of positions and periods where the sequences are frequent according to the input parameters. The package is also extensible, enabling users to incorporate heuristics and optimizations to drive the discovery of patterns.

Listing 1 :

 1 R example of the use of GSTSM # l o a d i n g t h e GSTSM package > l i b r a r y ( " gstsm " )

Acknowledgements

The authors thank CNPq, CAPES, and FAPERJ for partially sponsoring this research.


Polarization Amplification in Dynamic Nuclear Polarization Magic-Angle Spinning Solid-State Nuclear Magnetic Resonance by Solubilizing Traditional Ionic Salts

Florian Ferrer, Marie Juramy, Ribal Jabbour, Samuel Cousin, Fabio Ziarelli, Giulia Mollica, Pierre Thureau,* and Stéphane Viel* ABSTRACT: Dynamic nuclear polarization can improve the sensitivity of magicangle spinning solid-state NMR experiments by one to two orders of magnitude. In aqueous media, experiments are usually performed using the so-called DNP juice, a glycerol-d8/D2O/H2O mixture (60/30/10, v/v/v) that can form a homogeneous glass at cryogenic temperatures. This acts as a cryoprotectant and prevents phase separation of the paramagnetic polarizing agents (PAs) that are added to the mixture to provide the source of electron spin polarization required for DNP. Here, we show that relatively high 1 H DNP enhancements (~60) can also be obtained in water without glycerol (or other glass forming agents) simply by dissolving high concentrations of electrolytes (such as NaCl or LiCl), which perturb the otherwise unavoidable ice crystallization observed upon cooling, hereby reducing PAs phase separation and restoring DNP efficiency.

Dynamic nuclear polarization has emerged as an efficient method for improving the sensitivity of magic-angle spinning (MAS) solid-state NMR (ssNMR) experiments. 1 Thanks to the pioneering works of R. G. Griffin and coworkers, 2-4 sensitivity enhancements over one to two orders of magnitude can now be readily achieved for a wide range of samples, [START_REF] Rossini | Dynamic nuclear polarization surface enhanced NMR spectroscopy[END_REF][START_REF] Lee | Is solid-state NMR enhanced by dynamic nuclear polarization? Solid State[END_REF][START_REF] Zhao | DNPenhanced solid-state NMR spectroscopy of active pharmaceutical ingredients[END_REF][START_REF] Mollica | Structural investigations of polymer materials by dynamic nuclear polarisation solid-state NMR[END_REF][START_REF] Rankin | Recent developments in MAS DNP-NMR of materials Solid State[END_REF][START_REF] Chow | Biomolecular and biological applications of solid-state NMR with dynamic nuclear polarization enhancement[END_REF][START_REF] Moroz | Dynamic nuclear polarization solid-state NMR spectroscopy for materials research Ann[END_REF][START_REF] Biedenbänder | Dynamic nuclear polarization for sensitivity enhancement in biomolecular solid-state NMR[END_REF] including at high magnetic fields (9.4 T and above). [START_REF] Corzilius | High-field dynamic nuclear polarization[END_REF][START_REF] Berruyer | Dynamic nuclear polarization enhancement of 200 at 21.15 T enabled by 65 kHz magic-angle spinning[END_REF] DNP is based on the transfer of electron spin polarization to nuclei promoted by microwave irradiation of the sample, [START_REF] Overhauser | Polarization of nuclei in metals[END_REF][START_REF] Carver | Polarization of nuclear spins in metals[END_REF] usually at cryogenic temperatures. For diamagnetic samples, the required source of electron spin polarization is provided by specific paramagnetic species (so-called polarizing agents or PAs) [START_REF] Menzildjian | Polarizing agents for efficient high field DNP solid-state NMR spectroscopy under magic-angle spinning: from design principles to formulation strategies[END_REF] that are added to the sample of interest. In aqueous media, as is typically the case for biological samples, one of the preferred strategies consists of mixing water with organic compounds such as glycerol, 3 DMSO, [START_REF] Hu | Dynamic nuclear polarization with biradicals[END_REF] or trehalose [START_REF] Kaushik | Trehalose matrices for high temperature dynamic nuclear polarization enhanced solid state NMR[END_REF] in predetermined relative volume fractions to form an amorphous glass at cryogenic temperatures (ca. 100 K). Indeed, water is known as a bad "glass former", 20 and only cooling rates >10 6 K.min -1 have been shown to promote the formation of glassy water at atmospheric pressure (i.e. hyperquenching). [START_REF] Bruggeller | Complete vitrification in pure liquid water and dilute aqueous-solutions[END_REF][START_REF] Mishima | Melting ice' I at 77 K and 10 kbar: a new method of making amorphous solids[END_REF] These cooling rates are far superior to those that can be achieved by inserting the sample into the pre-cooled low-temperature MAS ssNMR probehead used for DNP, which can be estimated in the 10 2 -10 3 K.min -1 range under typical implementation conditions. [START_REF] Ghosh | Cryogenic sample loading into a magic-angle spinning NMR spectrometer that preserves cellular viability[END_REF] Therefore, without the use of glassy agents, ice crystallization cannot be prevented. [START_REF] Bachler | Glass polymorphism and liquid-liquid phase transition in aqueous solutions: experiments and computer simulations[END_REF] This may lead to PA phase separation and self-aggregation (or crystallization), which in turn translate into a loss of DNP efficiency. As such, any DNP MAS ssNMR experiment performed in pure water typically ends up with hardly any sensitivity enhancement at all. In contrast, using glycerol or alternative "glass forming" compounds, good-quality glasses can be obtained in which PAs can be homogeneously dispersed, hereby allowing DNP efficiency to be restored. In addition, in some circumstances, these mixtures allow the analyte of interest (e.g. biomolecules) to be structurally preserved at low temperature (cryoprotection). One example of these binary, glass-forming mixtures involves glycerol and water. In particular, a 60/30/10 (v/v/v) glycerol-d8/D2O/H2O mixture has found widespread use in the related literature, where it has been referred to as 'DNP juice' for some time, [START_REF] Poschko | Spin noise detection of nuclear hyperpolarization at 1.2[END_REF] allowing large DNP sensitivity enhancements (>200 at 9.4 T) to be reached for a range of samples. Decreasing the amount of glycerol usually decreases the glass stability of the sample, [START_REF] Leavesley | Effect of water/glycerol polymorphism on dynamic nuclear polarization[END_REF] unless special rapid freezing techniques are used. [START_REF] Jeon | Application of millisecond time-resolved solid state NMR to the kinetics and mechanism of melittin self-assembly[END_REF] In the framework of our ongoing work related to the investigation of crystallization processes of glycine by DNP MAS ssNMR, [START_REF] Vioglio | Insights into the crystallization and structural evolution of glycine dihydrate by in situ solid-state NMR spectroscopy[END_REF][START_REF] Vioglio | A strategy for probing the evolution of crystallization processes by low-temperature solid-state NMR and dynamic nuclear polarization[END_REF][START_REF] Juramy | Monitoring crystallization processes in confined porous materials by dynamic nuclear polarization solid-state NMR[END_REF][START_REF] Cousin | Exploiting solidstate dynamic nuclear polarization NMR spectroscopy to establish the spatial distribution of polymorphic phases in a solid material[END_REF] we unexpectedly observed in one particular study a substantial 1 H DNP enhancement (x8) for

Glycine, H2O

with AMUPol

Glycine, H2O

with AMUPol

+ NaCl (4 M)

With a pinch of salt! glycine in water when NaCl was used for selecting a particular glycine polymorph (whereas no DNP enhancement could be observed in the absence of NaCl). [START_REF] Juramy | Monitoring the influence of additives on the crystallization processes of glycine with dynamic nuclear polarization solid-state NMR Solid State[END_REF] In the present work, we aim to comprehensively investigate this effect. As such, we first studied the impact of the NaCl concentration (from 0 M to 5 M) by recording 1 H® 13 C cross-polarization magic-angle spinning (CP MAS) spectra at 9.4 T and ca. 100 K (with and without microwave irradiation) for a series of aqueous solutions (100% H2O) of [1-13 C]glycine (0.5 M) at 3 different AMUPol concentrations (5, 10 and 20 mM). Further experimental details can be found in SI. We chose AMUPol because it can provide high 1 H DNP enhancements at 9.4 T (>200), [START_REF] Sauvée | Highly efficient, water-soluble polarizing agents for dynamic nuclear polarization at high frequency[END_REF] and its DNP properties (e.g. depolarization) have been thoroughly described in the literature. [START_REF] Mentink-Vigier | Nuclear depolarization and absolute sensitivity in magic-angle spinning cross effect dynamic nuclear polarization[END_REF] For all 0.5 M [1-13 C]glycine samples, we measured two quantities. First, we measured the 1 H DNP enhancement defined as the intensity ratio of the 1 H® [START_REF] Corzilius | High-field dynamic nuclear polarization[END_REF] C CP MAS resonances recorded with the microwave field on and off (Ion/Ioff). Second, because AMUPol is affected by depolarization, 33 a process that artificially improves the amplitude of the DNP enhancement by reducing the spectral intensities measured without microwave irradiation, [START_REF] Thurber | Perturbation of nuclear spin polarizations in solid state NMR of nitroxide-doped samples by magic-angle spinning without microwaves[END_REF] we also measured the relative sensitivity. The latter was herein expressed as the ratio of the integral of the [START_REF] Corzilius | High-field dynamic nuclear polarization[END_REF] CO resonance of a given glycine sample analyzed with microwaves to the integral of the same resonance for a reference glycine sample in pure water analyzed without microwaves (and without AMUPol). Moreover, to account for differences in recycle delays, this ratio of integrals was divided by the square root of the respective experimental time constants (see SI for details).

Both the 1 H DNP enhancement and the relative sensitivity are valuable in assessing the actual improvement in experimental sensitivity achieved through DNP. Results are reported in Figure 1.

As can be seen in Figure 1, substantial 1 H DNP enhancements were observed (Figure 1a), even for moderately high NaCl concentrations (1 M), for which a 1 H DNP enhancement of 20 could be measured. For comparison, as mentioned above, the 1 H DNP enhancement of glycine in water (without NaCl) using the same range of AMUPol concentrations is 1. This increase in 1 H DNP enhancement also translated into significant improvement in terms of relative sensitivity, as shown in Figure 1b. For instance, the relative sensitivity measured for an aqueous solution (100% H2O) with NaCl (4 M) and AMUPol (10 mM) is about 58% of that measured on the same glycine sample at the same concentration but with DNP juice (evidenced by the grey dashed line in Figure 1b). Interestingly, a 2 M NaCl solution in water (100% H2O) with 5 mM of AMUPol performs almost just as well. This increase in sensitivity allows 2D experiments to be readily recorded, as illustrated with a 2D 13 C-13 C INADEQUATE experiment [START_REF] Lesage | Determination of through-bond carbon-carbon connectivities in solid-state NMR using the INADEQUATE experiment[END_REF] performed in ca. 25 min on a 0.5 M U- [START_REF] Corzilius | High-field dynamic nuclear polarization[END_REF] C glycine sample containing LiCl and AMUPol at 4 M and 15 mM, respectively (Figure S1). Moreover, a similar gain in sensitivity (56% with respect to DNP juice) was obtained for another water-soluble compound, nicotinamide, which is a natural form of vitamin B3 that is widely used to form cocrystals with drugs to modify their properties [START_REF] Sanphui | Cocrystals of hydrochlorothiazide: solubility and diffusion/permeability enhancements through drug-coformer interactions[END_REF] (Figure S2).

Furthermore, note that the DNP juice only contains ca. 10% of protons, as it is composed of glycerol-d8/D2O/H2O (60/30/10, v/v/v). [START_REF] Biedenbänder | Dynamic nuclear polarization for sensitivity enhancement in biomolecular solid-state NMR[END_REF] Using the same relative concentration of protons, the relative sensitivity of a 4 M NaCl solution of glycine with 10 mM AMUPol in a D2O/H2O mixture (90/10, v/v) reached 75% of that obtained with the DNP juice (Figure 2). Interestingly, in Figure 2a, the [START_REF] Corzilius | High-field dynamic nuclear polarization[END_REF] C resonances due to the glycerol are clearly visible, which makes it difficult to observe other signals not only in the 50-80 ppm spectral region, but also in other regions due to the presence of multiple side bands of substantial intensity that may also cause spectral interference with signals of lower intensity. This is even more pronounced when considering proline as an example because the glycerol signals overlap almost entirely with the alpha-CH signal due to proline (Figure S3).

In addition, we investigated the effect of using different types of salts (Figure 3). In this case, the salt concentration was maintained constant (4 M) for comparison and we only used plain water (100% H2O) for all the electrolytic solutions. Figure 3 shows that substantial sensitivity Relative sensitivity enhancements could be obtained for electrolytic solutions containing Na + or Li + ions, with a relatively low impact of the respective counterion (Cl -, Br -, or I -), whereas solutions containing K + ions only yielded disappointing results. For comparison, recall that the 1 H DNP enhancement and the relative sensitivity of 0.5 M [1-13 C]glycine in pure water (100% H2O) without any added salt and without AMUPol are both equal to 1.

The fact that we do observe such sizeable DNP sensitivity enhancements suggests that, in these experimental conditions, ice crystallization may be hampered somehow (with a concomitant reduction in PA phase separation). In fact, we observed that 1 H spin-lattice relaxation times of aqueous solutions (100% H2O) of NaCl without AMUPol increased steadily with NaCl concentration, from T1[100 K] = 1 s without NaCl (i.e. a situation where water is likely to be present in the form of ice only) to T1[100 K] = 20 s at 6 M NaCl (near the solubility limit at 298 K). This indicates a possible change in the state of matter of water, although variations in the concentrations of dissolved paramagnetic O2 in these samples could also be contributing factors. [START_REF] Bok | The solubility of oxygen in water and saline solutions Front[END_REF] More precisely, data in Figure 3 showed that the best results were obtained for NaCl and LiCl, with 1 H DNP enhancements that reached 73 and 128, respectively. Still, the best relative sensitivity was observed for the NaCl solution due to a lower 1 H DNP build-up time constant (τ H,DNP 1 = 0.74 s and τ H,DNP 1 = 1.1 s for NaCl and LiCl, respectively), and because the integral of the glycine signal in the spectrum recorded on the LiCl solution without microwaves was 40% of that measured for the off spectrum of the NaCl solution. This would suggest that signal losses due to depolarization are more pronounced for LiCl solutions than for NaCl solutions.

Overall, the trends observed in Figure 3 are in line with the so-called "glass-forming" properties of these electrolytic solutions. It is indeed well established in the literature that concentrated aqueous solutions of electrolytes containing sodium, lithium, or potassium ions, behave differently upon cooling with respect to pure water, [START_REF] Angell | Liquid fragility and the glass transition in water and aqueous solutions[END_REF] allowing for the formation of glasses depending on experimental conditions such as pressure or cooling rates. For instance, for LiCl/H2O solutions, it has been shown that even low cooling rates such as 10 K.min -1 (amenable to our DNP setup) can lead to solution amorphization at ambient pressure -yielding a rather homogenous glass -as long as the LiCl concentration is close or above the eutectic (7.5 M). [START_REF] Bachler | Glass polymorphism and liquid-liquid phase transition in aqueous solutions: experiments and computer simulations[END_REF] In this respect, we also analyzed a highly concentrated LiCl/H2O solution (9 M), well above the eutectic point (Figure 3). A 1 H DNP enhancement of 166 was obtained in this case -the highest value of all -but the corresponding relative sensitivity was comparable to that obtained for LiCl at 4 M, mainly because = 1.1 s). In contrast, the cooling rates achievable in our DNP conditions should not be sufficiently fast to obtain glassy NaCl/water solutions, even at the eutectic concentration, because cooling rates larger than 10 5 K.min -1 are required at ambient pressure, [START_REF] Ludl | Quenching device for electrolytic aqueous solutions[END_REF] unless the aqueous NaCl solutions are confined in nanosized pores (which is not the case here). [START_REF] Zhao | Confinementinduced vitrification of aqueous sodium chloride solutions[END_REF] Moreover, in that latter study, [START_REF] Zhao | Confinementinduced vitrification of aqueous sodium chloride solutions[END_REF] it was observed that nano-confinement could vitrify NaCl/H2O solutions but not KCl/H2O solutions, indicating that glasses from KCl/H2O solutions are even more difficult to obtain than those from NaCl/H2O solutions. [START_REF] Kanno | Homogeneous nucleation and glass formation in aqueous alkali-halide solutions at high-pressures[END_REF] In other words, while it is reasonable to assume that our experimental DNP setup could lead to the formation of glasses for the Li + solutions, this should not be the case for the Na + solutions (and even more so for the K + solutions). Clearly, one could argue that the pressure inside the spinning rotor is not exactly equal to the ambient pressure because of the fast-spinning speed of the rotor that exerts an additional centrifugal pressure on its inner wall (which can reach up to ca. 2 MPa when spinning at room temperature a 3.2 mm MAS rotor at 8 kHz). [START_REF] Bartalucci | Disentangling the effect of pressure and mixing on a mechanochemical bromination reaction by solid-state NMR spectroscopy[END_REF] Still, the pressures necessary to suppress crystallization for concentrated NaCl solutions (4.5 M) were shown to be superior to 100 MPa. [START_REF] Kanno | Homogeneous nucleation and glass formation in aqueous alkali-halide solutions at high-pressures[END_REF] Remarkably, however, the [START_REF] Ghosh | Cryogenic sample loading into a magic-angle spinning NMR spectrometer that preserves cellular viability[END_REF] Na MAS ssNMR spectra recorded at ca. 100 K on all the aqueous solutions of NaCl only showed a single and broad resonance at ca. -6 ppm (with a linewidth at half height of ca. 1400 Hz/13 ppm), which indicated that the sodium ions were not in a well-defined crystalline state but experienced a variation of local chemical environments that is symptomatic of a microscopically disordered state. This suggests that the presence of NaCl, even at low concentrations, affects ice crystallization in our experimental conditions, leading to a diminished or absent PA phase separation typically observed in water. Consequently, the restoration of DNP efficiency becomes feasible and reliable through an approach referred to as PASTIS (Polarization Amplification by Solubilizing Traditional Ionic Salts). We emphasize here that this approach is different from the work by Kubicki et al. where they added solid (insoluble) dielectric particles to the solution, resulting in an amplification of the microwave field in the sample and a concurrent increase in the DNP enhancement. [START_REF] Kubicki | Amplifying dynamic nuclear polarization of frozen solutions by incorporating dielectric particles[END_REF] PASTIS represents a possible alternative to the use of the DNP juice whenever glycerol should be avoided, [START_REF] Biedenbänder | Dynamic nuclear polarization for sensitivity enhancement in biomolecular solid-state NMR[END_REF] but it is likely to be more suitable for analyzing water-soluble molecules of small to moderate molecular weights rather than proteins or other large biological assemblies, since the structure of the latter can be altered in the presence of large amounts of electrolytes. Still, using these experimental conditions, it could be of interest to study other water-soluble PAs, 17 including at higher magnetic fields and faster MAS rates. Applications in dissolution DNP could also be explored. [START_REF] Jannin | Application and methodology of dissolution dynamic nuclear polarization in physical, chemical and biological contexts[END_REF] Investigations along these lines will be addressed in our lab.

n ASSOCIATED CONTENT

Experimental details regarding sample preparation, NMR measurements and data analysis, are available (PDF). §n AUTHOR INFORMATION

Corresponding Authors

Pierre Thureau -Aix-Marseille Univ, CNRS, ICR, 13013, Marseille, France; orcid.org/0000-0002-9157-256X Email: pierre.thureau@univ-amu.fr Stéphane Viel -Aix-Marseille Univ, CNRS, ICR, 13013, Marseille, France; orcid.org/0000-0003-3338-0559 Email: s.viel@univ-amu.fr

Figure 1 .

 1 Figure 1. Evolution of the (a) 1 H DNP enhancement and (b) relative sensitivity measured at 9.4 T and ca. 100 K on aqueous solutions (100% H2O) of 0.5 M [1-13 C]glycine with increasing concentration of NaCl for 3 distinct concentrations of AMUPol (5, 10, and 20 mM). In (b), the grey dashed line indicates the relative sensitivity obtained on a comparable glycine solution but in DNP juice (glycerol-d8/D2O/H2O, 60/30/10, v/v/v) with 10 mM of AMUPol (and without NaCl).

Figure 2 .

 2 Figure 2. DNP-enhanced 1 H® 13 C CP MAS spectra of [1-13 C]glycine (0.5 M) with AMUPol (10 mM) recorded at 9.4 T and ca. 100 K in (a) DNP juice and (b) a 4 M NaCl solution in D2O/H2O (90/10, v/v). The noise levels of both spectra are illustrated in the corresponding insets. These spectra have been normalized to the square root of their 1 H DNP build-up time constants (3.4 s and 1.6 s, respectively, see SI for details). The * and x symbols indicate spinning side bands due to glycine and glycerol, respectively. The 13 CO and CH2 resonances due to [1-13 C]glycine and the CH and CH2 resonances due to glycerol are shown with green and orange labels, respectively.

Figure 3 .

 3 Figure 3. Evolution of (a) the 1 H DNP enhancement and (b) the relative sensitivity measured in aqueous solutions (100% H2O) of 0.5 M [1-13 C]glycine in different salts (at 9.4 T and ca 100 K). In each case, the AMUPol concentration is 15 mM and the salt concentration is 4 M, except for the 3 rd case from the left (LiCl(*)) for which the LiCl concentration is 9.5 M.

= 4 . 4 s

 44 higher 1 H DNP build-up time constant (τ H,DNP 1 instead of τ H,DNP 1