Sébastien Becker

Emilie Devijver

Rémi Molinier

Noël Jakse

Glass-forming ability of elemental zirconium

published or not. The documents may come

I. INTRODUCTION

The possibility to form a metallic glass (MG) from the melt was reported for the first time [1] in the 60's by splat quenching of a eutectic Au-Si alloy with cooling rates as high as 10 6 K/s. This has triggered fundamental questions regarding the possibility of forming a single-element metallic glass by cooling the melt at a sufficiently high cooling rate to avoid crystallization. Subsequent attempts for pure nickel by splat-quenching with cooling rate of 10 10 K/s showed an amorphous state with a crystallization temperature of 425 K higher than the room temperature [2], but it was suspected to be stabilized by impurities such as oxygen or carbon. Earlier studies using quench condensation on substrates [3] as well as other techniques [4] seemed to concur with the impossibility to form a metallic glass with a single element. Later, atomized iron droplets of about 30 nm in size were however produced [5,6] with mixed crystal and amorphous structure at large undercoolings with cooling rates of 10 7 K/s. It was further shown that a monoatomic glass is more easily formed from bodycentered cubic (bcc) than for face-centered cubic (fcc) or hexagonal closed-packed (hcp) metals [6]. This finding for bcc metals was further confirmed very recently on nanorods with cooling rates as large as 10 14 K/s [7]. This implies that homogeneous crystal nucleation for such large undercooling is impeded in the case of droplets while the growth might be suppressed in the case of the nanorods as they were attached to their crystalline ends [8].

As crystal nucleation is essentially initiated by heterogeneous nucleation from impurities, surfaces, or near grain boundaries [9], homogeneous crystal nucleation is driven by complex phenomena [10] difficult to capture experimentally [11] and thus remains a challenging issue.

Despite similarities of the liquid and the crystal in their density and coordination number, the seminal work by Turnbull [12] on undercooled metallic melts demonstrated that significant undercooling can be obtained. This implied that nucleation is an activated process in which the liquid transforms to the crystal by overcoming a free energy barrier. This is possible only if the liquid possesses an atomic structure substantially different from that of a simple crystal, and orders in the vicinity of a crystal surface [13]. It was then suggested by Frank [14] that the energy minimization in a monoatomic liquid favors the formation of local icosahedral structure. Such a five-fold symmetry, predominantly built from slightly distorted tetrahedral structures, is incompatible with simple closed-packed lattices such as the bcc, fcc and hcp crystals that satisfy a global energy minimum [START_REF] Nelson | Solid State Physics[END_REF].

Investigations of the interplay between homogeneous nucleation and local polymorphism [START_REF] Ten Wolde | [END_REF] remain experimentally challenging to observe [11,[17][18][19][20][21][22], especially at large undercoolings ∆T = T -T M , T being the temperature of the studied liquid and T M its melting or liquidus temperature. They were often carried out using direct particle based simulation methods like molecular dynamics (MD) or Monte-Carlo (MC) [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Smit | Understanding molecular simulations[END_REF]. For Al, a prototype fcc metal, nucleation might occurs primarily by formation of fcc nuclei along with hcp defects [START_REF] Mahata | [END_REF]26]. The situation is less clear for Cu since former MD simulations [27] showed an initial stage of nucleation with a fcc local ordering while very recent ones [28] involve metastable bcc polymorph pre-critical nuclei. For Fe and Mg, early stages of nucleation involve formation of respectively bcc and hcp phase corresponding to their stable phase.

In Na [29], another bcc metal, it was shown that the supercooled liquid, with significant icosahedral and defective icosahedral local order, transforms during the nucleation process into the stable bcc structure via defective bcc local structure. Such a competition between the bcc and icosahedral local orderings in the undercooled liquid was observed in Zr [30,31] and Ta [32] by means of ab initio molecular dynamics (AIMD) simulations. However, for these two metals, it was shown that the structural ordering can be interpreted on the basis of a competition between a polytetrahedral local order, more complex that the simple icosahedral one, and a bcc-type ordering. For Zr, the bcc-type ordering increases with the degree of undercooling together with more numerous Frank-Kasper (FK) polyhedra such as Z14, Z15, and Z16 [33]. A similar situation occurs for Ta but with a local order close to that of the A15 phase.

Complex polytetrahedral orderings found in some bcc elemental undercooled liquid metals [30][31][32] might therefore play a significant role in the frustration against crystallization scenario suggested by Tanaka [34,35] when approaching the glass transition temperature T G . Interestingly, for pure Ta, formation of a stable glass was observed [7] that could be a result of this concept of frustration. However, the situation is less clear for Zr since attempts to form glassy spheres using atomized droplets [5,6] failed and formation of a metallic glass at high pressure [36] was subject to misinterpretation of the experiments [36,37]. Therefore, the question of the existence of a stable amorphous state for Zr remains open.

The aim of the present work is therefore to study the glass formation as well as the onset of crystallization features of pure Zr as a function of the degree of undercooling ∆T in the range between T M and T G , by means of large-scale MD simulations [38] up to one million atoms. For this purpose we have designed a semi-empirical potential within the framework of the Modified Embedded Atom Model (MEAM) [39] based on a bcc reference structure, which reproduce the structural and thermodynamic features of the liquid state and solid β-phase in the vicinity of the melting point. We determined the glass transition temperature and timetemperature transformation (TTT) diagram at large undercoolings. Our findings indicate that elemental Zr might form a metallic glass with cooling rates higher than 10 12 K/s, now being achievable experimentally [7], and results from a strong competition between bcc and fivefold FK polytetrahedral local structural orderings, impeding homogeneous nucleation.

The layout of the paper is the following. Section II is devoted to the technical aspects of the simulations and the design of the potential. Section III presents and discusses the optimization and tests of the potential as well as the results of the simulations. Finally in Section IV we draw our conclusions.

II. SIMULATION BACKGROUND

A. Modified Embedded Atom Model

Several semi-empirical interatomic potentials have been developed for atomistic simulations for pure Zr including second-moment tight-binding method [40], Embedded Atom Model (EAM) [41], MEAM [39,42,43], Reax Force Field [44] as well as Machine-Learning technique within the Kernel-Ridge Regression [45]. However, only some of then were designed to take into account the martensitic hcp-bcc transformation at ambient pressure [40,41,45]. Moreover, none of them were designed to take into account the liquid state properties explicitly and specifically local structural properties, self-diffusion and the melting temperature. In the context of the present work, i.e. the study of the liquid and undercooled properties of Zr, our strategy is then to develop a MEAM potential using a bcc reference structure which corresponds to the underlying high temperature crystalline structure (β-phase). In the original MEAM potential, many-body and directional bonding aspects are taken into account in the first neighborhood, but here we use its second nearest-neighbor (2NN) extension [46] in order to take larger varying atomic environment into account. Considering a system of N atoms, the potential energy functional in the MEAM formalism can be expressed in a general form as

E = N i=1   F i (ρ i ) + 1 2 N j =i=1 Φ(r ij )   , (1)
where F i (ρ i ) is the embedding energy function depending on the background electronic density ρ i , Φ(r ij ) is the pair potential interaction as a function of interatomic distance r ij between atoms i and j. The MEAM formalism is well documented in the literature, therefore we refer the reader to Refs. [46] and [43] for a detailed description.

The electron density is calculated first for each atomic site from an analytic expression taking into account the length and directionality in bonding of the neighbors [46], and the background electron density is taken in the form 2ρ i /(1 + exp(Γ)), where Γ is electronic density expansion over the t (k) (k = 0, . . . , 3) electronic structure parameters. Then from a specific form of the embedding function and the total energy, the pair potential Φ(r) in Eq.

(1) is estimated using the equation of state of Rose [48]

E Rose = -E coh 1 + a * + da * 3 e -a * , (2)
where a * = α(r/r e -1) where r e is the nearest-neighbor distance in the bcc reference structure and α = (9BΩ/E coh ) 1/2 , B, being the bulk modulus, E coh the cohesive energy, and Ω the atomic volume. The 2NN part is controlled through the screening parameters C min and C max , with C min taking values below 2, and we follow strictly the original paper by Lee et al.

[46] which is implemented in LAMMPS [38]. Accordingly, the cutoff radius of the potential takes a value of 7.60 Å.

From a previous attempt [47], we kept the electronic structure parameters β (k) and t (k) (k = 0, . . . , 3) and optimized the parameters r e , d, C min , and C max only. The value of r e was adjusted to reproduce at best the density of the liquid and bcc crystalline phase, and d was determined to get at best the experimental melting temperature of Zr. It is worth mentioning that d is usually fitted on the pressure variation of the bulk modulus in the solid phase (∂B/∂P ). However, the amplitude of variation of d considered here should not alter much this quantity [43] while having a significant impact on melting, as it affects essentially on the well depth of the pair potential φ(r ij ) (see Eq. 2). The screening parameters C min , C max were fitted to get the better representation of at least the two first peaks of the pair correlation function obtained from the AIMD simulation (see Sec. II B). et al. [43] C min , and C max correspond to the best compromise. It should be also noted that the structural energy difference ∆E hcp→bcc = 0.017 eV is consistent with the value obtained by them [43] and similar values of the mechanical properties than those reported previously [47] are found.

Interestingly enough, the values of the hcp version of the MEAM potential obtained by Kim

E coh (eV) r e ( Å) α A β (0) β (1) β (2) β (3) t (1) t (2) t (

B. Ab initio molecular dynamics

In order to provide reference values for the pair-correlation function and self-diffusion coefficient for the MEAM potential fitting, AIMD simulations have been performed using the VASP code [49,50]. The projected augmented-wave (PAW) method was used to describe the electron-ion interaction [51,52] and the exchange-correlation energy was taken from the generalized gradient approximation in the PBE form [53]. In the present PAW potential, the 5s, 4d orbitals as well as the semi-core 4s and 4p orbitals are treated as valence orbitals with a plane-wave cutoff of 300 eV. Only the Γ-point has been considered to sample the supercell Brillouin zone.

The simulations were carried out in the NVT ensemble (constant volume V , temperature T , and number of particles N ), with a Nosé thermostat [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Smit | Understanding molecular simulations[END_REF], similarly than in our previous work [30] but using N = 256 atoms in order to have a more refined description of the structure and dynamics in the liquid state. Equations of motions were solved using Verlet's algorithm in the velocity form with a time step of 1 fs. The volume V of the cell was fixed to reproduce the experimental densities [61] and was subject to standard periodic boundary conditions. We have considered T = 2500 K, well above the experimental melting temperature T M,exp = 2128 K. The simulation was started from an initial configuration with random positions and equilibrated at T = 2500 K for 10 ps. The run was continued for 30 ps to calculate the pair-correlation function, and the velocity auto-correlation function.

C. Classical molecular dynamics

All MD simulations were performed with the LAMMPS code [38] in the isobaricisothermal ensemble (constant temperature, pressure and number of particles) with ambient pressure, using the Nose-Hoover thermostat and barostat [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Smit | Understanding molecular simulations[END_REF] with periodic boundary conditions (PBC) in the three directions of space. As for AIMD, Verlet's algorithm in the velocity form was used for the numerical integration of the equations of movements, with a time step of 2 fs. We conducted two types of simulations whose characteristics are described in Table II: (i) MD with N 128000 atoms of liquid-solid interfaces (LSI) for the purpose of determining the melting temperature by the two-phase coexistence technique with different crystalline orientations; (ii) melt-quenching simulations of bulk systems with N = 5488, 16000 or 1024000 atoms down to the glass transition temperature as well as simulation of isotherms in the undercooled region to determine the time-temperature-transformation curves from homogeneous nucleation events.

The liquid-solid interface simulations are firstly carried out for the determination of the melting temperature T M of the MEAM potential. The procedure is similar to the two-phase coexistence approach formerly introduced by Morris et al. [54,55] and slightly improved subsequently by Sun et al. [56]. A simulation cell containing around N = 128000 atoms was set up with an initial crystalline configuration with PBC applied to the three directions of space. The crystal is heated and equilibrated at constant pressure to a temperature 100 K below a guess of T M . We observed that the system remained fully crystalline. Half of the simulation cell is subsequently heated and maintained at a temperature approximately 1500 K above T M until a complete melting was observed. The liquid part was then cooled down and equilibrated at a temperature 100 K above the estimated T M . This procedure creates a coexisting solid-liquid simulation cell containing two crystal-melt interfaces due to the PBC.

The simulation of the entire system is pursued in the isobaric-isoenthalpic ensemble in which the simulation box was allowed to change its length only along the x direction in order to keep a zero pressure. As a result, the temperature of the coexisting liquid-solid interface evolved toward a steady state corresponding to the thermodynamic melting temperature.

The simulation is continued for 1 ns, and the average melting temperature is determined on the last 100 ps if a steady position of the two interfaces is observed. If a complete melting or solidification occurs, the procedure is started over again with a refined guess of the melting temperature. This procedure is applied for three crystalline orientations, i.e.

(100), (111), and (110) whose characteristics are described in Table II. Two box shapes are also considered for the (100) interface for which a snapshot obtained with the OVITO software [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO -the Open Visualization Tool[END_REF] is shown in Figure 1.

For the melt quenching of the bulk systems, simulations were started at a temperature T = 2500 K from an initial configuration with random atomic positions. After equilibration at this temperature, the system was cooled down to 300 K with a cooling rate of 10 12 K/s. During this process, thermodynamic properties are calculated and configurations are recorded each 1000 time steps, for the purpose of starting subsequent MD simulations along chosen isotherms in the undercooled and glassy states. It is worth mentioning that for Bulk 2 system, various cooling rates were explored, namely 10 13 , 10 12 , and 10 11 K/s. While for the two faster cooling rates, homogeneous nucleation can be avoided, for 10 11 K/s crystallisation during cooling was observed, guiding our choice for the cooling rate with Bulk 3 system to study the homogeneous nucleation process.

The average position location of each solid-liquid interface and the local structural anal- ysis were determined using, as a descriptor, the common-neighbor analysis (CNA) [58] with the indexing of Faken and Jonsson [59] and a bond-based algorithm as implemented in the OVITO software [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO -the Open Visualization Tool[END_REF] where a uniform cut-off radius corresponding to the first minimum of should be mentioned that a better adjustment could have been obtained, but d has also a correlated effect together with r e on the density and a reasonable compromise between the two physical quantities had to be made.

N Shape T M (K) ρ L M (K) ( Å) -3 ρ S M (K) ( Å) -3
The liquid and solid densities at melting, extracted respectively from the liquid and solid parts of each LSI simulations, are given in Table II. Again, essentially no shape and orientation effect is seen. The various measurements of the liquid densities and their temperature evolution span over a period of 7 decades and were reviewed recently [60,61].

They display a significant departure of 13% between 0.0363 Å-3 and 0.0412 Å-3 . Our gives the probability of finding a particle j at distances r relative to a particle i located at the origin, and reads: n(r) represents the mean number of particles j in a spherical shell of radius r and thickness ∆r centered on particle i. Integrating g(r) up to their first minimum gives access to the partial coordination numbers. The comparison with AIMD simulations and existing experimental [64] data are shown in Fig. 2(a). For T = 2200 K and T = 1850 K the AIMD curves are taken from our previous work [30]. A good agreement is found with AIMD simulations indicating that the functional form given by Eq. ( 1) in the framework of the 2NN-MEAM [46] is able to reproduce microscopic structural aspects of Zr in the liquid state and their temperature evolution in the undercooled region. Both classical and ab initio MD simulations show a reasonable agreement with the experimental data [64], especially the position of the first peak and subsequent oscillations indicating that the bond length is well reproduced.

g(r) = N V n(r) 4πr 2 ∆r . (3)
Taking a cutoff radius of 4.31 Å. the resulting coordination number is 12.4 at T = 2500 K which is a bit higher than experimental value of 11.9 at 2290 K but still slightly lower than the ab initio value of 13. The departure being less than one atom, the coordination number obtained for the MEAM potential is deemed reasonable.

B. Tests of the MEAM potential

The reliability of the potential can be tested on the dynamics as well as some thermodynamic quantities. As for the short-time dynamics, we have considered the velocity auto-correlation function (VACF) :

Ψ(t) = 1 N N i=1 [v i (t + t 0 ) • v i (t 0 )] 2 t 0 , (4)
where v i (t) denotes the velocity of atom i at time t. The angular brackets correspond to an averaging over time origins t 0 . Fig. 2(b) displays the velocity auto-correlation function normalized to their value at time t = 0, Ψ(t = 0), determined from the MEAM potential and compared to the AIMD results at T = 2500 K. Both curves show a negative well at short time, characteristic of a back-scattering effect due to the first neighbor atoms shell, so-called cage effect, followed by a rapid damping towards zero, attesting the relaxation of the system on the sub-picosecond time scale. The cage effect is more pronounced for the MEAM, leading to a lower diffusion, however its characteristic time is similar to the AIMD.

The self-diffusion coefficient D is determined from the linear behavior at long times of the mean-square displacement written as

R 2 (t) = 1 N N l=1 [r l (t + t 0 ) -r l (t 0 )] 2 t 0 , (5)
where r l (t) denotes the position of atom l at time t. The long time behavior of R 2 (t) is linear, which is a characteristic feature of the liquid state, as shown in the Inset of Fig. 2(b).

For the MEAM potential D = 0.34 ± 0.07 Å2 /ps, which underestimates the AIMD value, i.e., D = 0.53 ± 0.1 Å2 /ps but still remains in reasonable agreement.

Turning now to the thermodynamics, we consider the enthalpy which is directly calculated from positions and velocities during the MD simulation by

H(P, T ) = 1 2 N i=1 mv 2 i + 1 2 N i =j=1 u(r ij ) + P ext V, (6)
where P ext is the pressure imposed to the simulation box, and the volume V is a dynamical variable, m being the mass of the atoms. The heat capacity at constant pressure C P is then determined from a numerical derivative of the enthalpy using its standard definition:

C P (T, P ) = ∂H(T, P ) ∂T . ( 7
)
The latent heat of fusion, also called the enthalpy of melting is an important quantity for solidification phenomena. The latter can be determined from the enthalpy difference between the liquid and solid branches at the melting temperature T M . The enthalpy curves of Zr were determined at ambient pressure for the solid and liquid branches of the MEAM potential using MD simulation with Bulk 1 system (see Table II), and are drawn in Fig. 3.

The simulation is started at T = 300 K with a perfect bcc crystal. A simulation is performed during 100 ps (50 ps equilibration and 50 ps production) and followed by a step-wise heating with a temperature step of 100 K which amounts to an average heating rate of 10 12 K/s.

The average value of the enthalpy is calculated from the production stage and error bars are evaluated from the standard deviation of the Gaussian distribution. The production time along the isotherms is not long enough to observe a transition bcc-hcp transition below the experimental martensitic transition at 1139 K, which occurs on the 10 second time scale as shown in recent electrostatic levitation measurements [65]. The procedure is repeated until a dynamic melting is observed in the simulation box, which occurred at T = 2400 K, significantly higher than the thermodynamic melting temperature T M determined from LSI simulations due to overheating effects. For the liquid branch, the simulations are started at T = 3400 K with a disordered initial configuration.

The same procedure as for the solid branch is followed but with a step-wise cooling down to 300 K. Above T = 2400 K, the difference in the enthalpy from the heating and cooling processes is negligible, indicating that there is no reminiscence of the crystalline state. The cooling procedure for the liquid branch is repeated with Bulk 3 system and a cooling rate of 10 12 K/s, and as can be seen on Fig. 3, no size effect can be detected. From the liquid and solid branches we find an enthalpy of melting of 18.2 kJ/mol which compares quite well with the assessed experimental value [66] of 15.6 kJ/mol. Taking the numerical derivative of solid branch yields a value of C P at the melting point of 34.4 J/mol/K which is in very good agreement with the experimental value of 35.7 J/mol/K. For the liquid a value of 34.5 J/mol/K is obtained which understimate the experimental data of 39.7 J/mol/K [66] and of 39.9 J/mol/K [63].

All these results bring us to the conclusion that the present optimized MEAM potential is reliable to treat the structural, dynamic and thermodynamic properties of Zr at ambient pressure.

FIG. 3: Enthalpy as a function of temperature for the solid and liquid branches (see text for details). The red line marks the slope of the glassy states as a guide for the eyes for the crossover between the liquid and glassy regimes.

C. Glass transition temperature

The melt-quenching simulations described in Section III B and shown in Fig. 3 reveal that with a cooling rate of 10 12 K/s crystallization can be avoided. The enthalpy H(T ) exhibit a change of slope around 1000 K that signals that the glass transition sets in. It is worth mentioning that similarly Mo et al. [67], using an EAM potential [68] from a highly optimized ab initio procedure [69], could form a glass at a cooling-rate of 5 × 10 12 with a glass transition in the same range of temperatures as in the present work. Due to the fact that the glass transition depends on the cooling rate [70,71], the temperature at which it is observed, so-called fictive temperature, is generally higher that the genuine value of the glass transition temperature, T G . A more refined value of T G can be obtained based on the concept of Potential Energy Landscape (PEL) [72] from the evolution of the inherent structure energy (ISE) [73] as a function of temperature, which is displayed in Fig. 4. For this purpose, atomic configurations were recorded at the interval of 100 K during the quench run of Bulk 3 system at 10 12 K/s (from 4000 K to 300 K). At each temperature, separate production runs were performed and five independent atomic configurations were chosen in order to extract their inherent structure (IS). The IS are obtained by minimizing their energy by means of a conjugate gradient algorithm to bring the system to its local minimum in the PES, and an average value of the ISE is taken over the 5 configurations for each temperature.

In the high-temperature range corresponding to the free diffusion regime, an almost linear behavior is found. At about T D = 2500 K, the ISE starts to decrease more rapidly in the landscape influence regime. Note that the crossover T D occurs above the melting point T M , which seems to be a universal feature for monoatomic liquids [77]. Finally, it is almost constant as the system, having found a deep minimum, can be considered as undergone the glass transition. From the crossover between the landscape influenced and glassy regimes [76],

we estimate the glass transition temperature to be T G = 890 K which is significantly lower than that obtained from the H(T ) curve shown in Fig. 3. The ratio between the melting temperature and the glass transition temperature T rg = T G /T M , representing also one of the simplest glass forming ability (GFA) criterion [85], takes a value T rg = 0.42 from our simulations.

The microscopic structural features of this glass forming ability is depicted in Fig. 5 using the common-neighbor analysis. The CNA bonds were determined from the same inherent structures from the ISE curve shown Fig. 4 using a cutoff being the first minimum of the corresponding pair-correlation function. At T = 2500 K the proportion of the pairs is consistent with the AIMD results [30]. The main feature emerging from Fig. 5 is an increase of the 555, 444 and 666 pairs and a decrease of the 544 and 433, representing defective 555 pairs, with decreasing temperature up to the glass transition, while the 421 and 422 pairs, correspond to the fcc and hcp structures, remain always small. The simultaneous increase of the 555 together with the decrease of the 544 and 433 indicates that the five-fold symmetry becomes better defined. Given the evolution of the amount of 433 pairs, it cannot be excluded that their decrease participates also to the increase of the 444 ones. The significant amount of 444 and 666 pairs signals the presence of bcc-like ordering. In the perfect bcc crystalline structure the ratio of the proportions of 444 to the 666 pairs is 0.75, while the ratio observed in the simulations done here is around 0.60. This implies that a non negligible amount of 666 pairs can be associated to the 555 pairs to form Frank-Kasper (FK) polyhedra as described in details in Ref. [START_REF] Nelson | Solid State Physics[END_REF]. More precisely, FK polyhedra Z14, Z15 and Z16, with 14, 15 and 16 neighbors have a specific signature in terms of the CNA: 12 pairs 555 with 2, 3 and 4 pairs 666, respectively, which can be easily detected in a simulation. A similar evolution of the basic pairs have been obtained recently by Su et al. [74] using an EAM potential extracted from the binary Cu-Zr system [75]. This is also the case for the results of Mo et al. [67] mentioned above.

The fact that the coordination number of the IS evolves from 13.24 to 13.64 (from 2500 K to 500 K) is consistent with this view. Moreover, a detail inspection of full polyhedra in a typical glass configuration (T = 300 K) indicates the presence of 20186 Z14, 11373 Z15, and 1793 Z16, as well as 59894 bcc polyhedra. It is worth mentioning that these polyhedra are separated and homogeneously distributed in the configuration. The simultaneous increase of FK and bcc polyhedra during cooling, and especially Z14 with the same coordination number as the bcc highlights the competition between bcc and Frank-Kasper polytetrahedral fivefold symmetry orderings and concur with the concept of geometrical frustration propounded by Tanaka [34,35].

D. Time-temperature-transformation curve

To evaluate more precisely the conditions for which Zr can form a glass or crystallizes, the crystal nucleation time as a function of degree of undercooling ∆T is quantified. The resulting time-temperature-transformation curve is shown in Fig. 6. The melt quenching simulation Bulk 3 system at a cooling rate of 10 12 K/s was considered. For temperatures at intervals of 50 K in the range from 1400 K to 1000 K (as well as additional temperatures i.e. 1125, 1175 et 1275 K), a configuration is taken during the quench and considered as an initial configuration for an isothermal simulation at constant pressure. This range of temperatures is well below the melting temperature and just above T G , and corresponds to deep undercooling conditions with fast homogeneous nucleation rates. At lower ∆T , the nucleation time becomes large and beyond the reach for standard molecular dynamics used here. It would require advanced techniques like umbrella sampling [78], metadynamics [START_REF] Laio | Proc. Natl. Acad. Sci[END_REF], various forward-flux sampling methods [START_REF] Allen | [END_REF][81][82], or seeding approaches [83] that will be considered elsewhere.

For each temperature, the crystal nucleation is detected by a significant drop in the potential energy with time and a threshold of 40% of the atoms having a crystalline environment using the CNA. The mean nucleation time for each ∆T is evaluated as the average over five independent observations, each of them being carried out with the same initial configuration but having a different random velocity distribution. The resulting TTT curve exhibits a typical nose shape with a minimum nucleation time of 0.39 ns around T N = 1250 K. As an important feature, the million-atom simulation cells are large enough to observe numerous nucleation events, as shown in Fig. 7. This is not always the case for Bulk 2 system with N = 16000 atoms resulting in significant size effects, as was also pointed out very recently [26] for the study of homogeneous nucleation phenomena. As a matter of fact, repeating the simulation procedure for N = 16000 (Bulk 2 system) leads to a TTT curve with a nose occurring at a significantly lower ∆T , as can be seen in Fig. 6, and the simulation cell contains hardly more than one critical nucleus at the onset of nucleation (see Fig. 7). It is worth mentioning that, in principle, the fraction of crystalline structure follows Avrami's equation involving multiple nucleation events, which should in turn be fitted to infer the nucleation time. However, aiming to compare with Bulk 2 system, the chosen methods remains a good approximation for single nucleation events, given the fact that for such monoatomic systems nucleation is extremely rapid as can be seen in the inset of Fig. 6. This is an illustration that nucleation phenomena can hardly be treated completely and accurately with simulation box size of the order of 10 4 atoms only.

E. Discussion

Renewed interests in glass formation in pure metals was recently triggered by the possibility to achieve experimentally cooling rates comparable to what is currently carried out by MD simulations [7]. The ability to form a glass is favored in systems where the melting temperature T M is closer to the glass transition one, T G , as for the case of alloys with eutectic composition [1] for which it is the closest. Therefore, high values of the reduced glass T rg = T G /T M might favor glass formation. From Turnbull's observations [85], the nucleation rate maximum occurring between T M and T G is strongly correlated with T rg . Thus it might be easier to form a glass by suppressing homogeneous nucleation for metals having high values of T rg . For the case of Zr we are concerned about here, MD simulations yield a value T rg = 0.42 significantly higher that the standard rule for pure metals T G ≈ 0.3T M [11] but obviously still lower than T rg = 2/3 representing a good GFA that is nearly achieved for For all temperatures 40% corresponds roughly to the half of the potential energy drop during crystallization and is taken as the measure of the nucleation time.

instance for Cu-Zr binary alloys [86].

In recent experiments, Zhong et al. [7] showed the possibility to form a glass for pure bcc metals (Mo, V, W, and Ta) by nanosecond electrical pulses on nano-bridges achieving cooling rates as high as 10 14 K/s. The fact that both ends of the nanobridges found to be glassy were crystalline suggests primarily that the crystal growth was suppressed. However, due to the high cooling rate and high purity of the samples, the possibility of suppressing homogeneous nucleation on the short time scale cannot either be excluded in an early stage.

Focussing on pure Ta, the comparative analysis of Orava and Greer [84] reports a relatively experimentally.

IV. CONCLUSIONS

In summary, a semi-empirical potential for elemental Zr was designed in the framework of modified embedded atom model within the second nearest neighbor formulation [46], for the purpose of performing large-scale molecular dynamics simulation of stable liquid, metastable undercooled and amorphous phases as well as homogeneous nucleation phenomena of the high temperature bcc crystalline phase. Using liquid-solid interface MD simulations, it is

shown that the interaction model reproduces accurately experimental densities of the liquid and solid states in the vicinity of the melting point as well as a thermodynamic melting temperature T M = 2119 K close to the experimental value.

The MEAM potential is able to predict dynamics properties such as the self-diffusion coefficient and the velocity auto-correlation functions as compared with ab initio molecular dynamics simulations. Furthermore, it gives a good representation of the thermodynamics properties near melting such as the enthalpy of melting of 18.2 kJ/mol and the constant pressure specific heat of 34.4 and 34.5 J/mol/K, respectively for the solid and the liquid, in good agreement with assessed experimental data.

From a melt-quenching simulation with one million atoms at a quenching rate of 10 12 K/s, a glass transition temperature T G = 890 K is observed from the evolution of the inherent structure energy, namely the crossover between the landscape-influenced and the glassy regimes of the potential-energy lanscape. It yields a glass forming ability in the sense of Turnbull T rg = T G /T M = 0.42 which is higher than the common rule T rg 0.3 and similar to other bcc elemental metals [84]. The time-temperature transformation curve, determined on isotherms from initial configuration taken during the quench, displays a nose in the deep undercooling regime around T N = 1250 K with a minimal average nucleation time of 0.39 ns.

This finding implies the possibility for elemental Zr to form a metallic glass with cooling rates higher than 10 12 K/s, achievable experimentally [7], and indicates that the formation of the glass might come from the avoidance of homogeneous nucleation by geometrical frustration [34,35] resulting from a strong competition between bcc and fivefold FK polytetrahedral structural orderings. How this competition impacts the crystallisation kinetics

still remains an open question. The present work bring the relevant computational context

FIG. 1 :

 1 FIG. 1: Upper panel: snapshot of the liquid-solid (100) interface with 640 × 10 × 10 shape; Lower panel: snapshot of the liquid-solid (100) interface with 80 × 40 × 20 shape.

  the pair-correlation function of the liquid is applied to create bonds between pairs of particles. The CNA classifies pairs around each atom by sets of three indices: the first index represents the number of near neighbors common to this pair, the second index corresponds to the number of nearest-neighbor bonds among the shared neighbors, and the third index indicates the longest chain of bonded atoms among them. For instance, 421 and 422 bonded pairs are characteristic of close packed structures fcc and hcp, respectively. The occurrence of 444 and 666 pairs, with specific proportions, signals the presence of bcc ordering. The degree of five-fold symmetry is obtained from the proportion of 555, 554 and 433 pairs, which represent perfect (555) and distorted FFS based motifs. The system at the interface is considered as solid at a given distance if at least 50% of the atoms have a local environment characterized as bcc according to the CNA.III. RESULTS AND DISCUSSIONA. Optimization of the MEAM potentialIn a first stage, the optimization of the parameters d, r e , C min , and C max was carried out to reproduce simultaneously at best the experimental melting temperature, the density of the liquid and bcc solid phase, and the two first peaks of the pair-correlation function, respectively. TableIIgathers the values of melting temperature determined from the liquidsolid interface for the three crystalline orientations and the two box shapes (see Fig.1). All the values are mutually compatible indicating that shape and crystalline orientation effect are essentially negligible. Finally, T M = 2119 K was accepted as the average of the 4 values determined, which is in very good agreement with the experimental one, namely 2128 K. It

  simulated liquid densities at melting are close to the most recent experimental value of 0.0396 Å-3 measured by means of electrostatic levitation technique by Wang et al. [61], with a departure of 4.5%. The density of the bcc solid at melting compare well with the experimental data of Lu et al. [62] with underestimation of 4% only. Interestingly, solidliquid density difference in the present simulations is ∆ρ S-L = 0.007 Å-3 , close to the only existing consistent measurements of Paradis and Rhim [63] giving 0.006 Å-3 . The values of the screening parameters C min and C max are optimized at the microscopic level on the pair-correlation function, g(r), in the liquid state. Taking the AIMD curve obtained in the present work at 2500 K as a reference, they have been obtained by minimizing the least-square deviation. The function g(r)

FIG. 2 :

 2 FIG. 2: (a) Pair-correlation functions at T = 2500 K, T = 2200 K, and T = 1850 K, from the MEAM potential, AIMD of the present work as well as from the previous one [30], and compared to the experiment [64]; (b) Velocity-autocorrelation at function at T = 2500 K from the MEAM potential, AIMD of the present work.

FIG. 4 :

 4 FIG. 4: Inherent-structure energy as a function of temperature for Bulk 3 system with N = 1024000 (see text). Open symbols correspond to the ISE of the perfect crystalline phases.

FIG. 5 :

 5 FIG. 5: Temperature evolution of the abundance of CNA bonds during the melt quenching. The vertical yellow bar indicates the glass transition region.

FIG. 6 :

 6 FIG. 6: Time-Temperature-transformation curve for Bulk 3 system with N = 1024000 and Bulk 2 system with N = 16000 in the temperature range between T = 1400 K and T = 1000 K (see text). Inset: Potential energy as a function of time for all the temperatures in the same temperature range. For each temperature, the chosen curve corresponds to the simulation in which the crystallization event occurs at a time closest to the average nucleation time. Arrows indicate the location on the potential energy curve where 15% and 40% of the atoms have a full bcc structural environment.

TABLE I :

 I 3) C min C max d 6.29 3.10 4.27 0.72 4.10 1.00 1.00 1.00 5.80 -0.350 -1.30 1.00 1.44 -0.069 Parameters of MEAM potentials (see text for their definition).

TABLE II :

 II Parameters of classical molecular dynamics simulations for bulk and solid -liquid interface. N stands for the number of atoms and the shape of the simulation box is given in units of primitive bcc cells in the x, y, and z direction. T M , ρ L M and ρ S M are respectively the calculated melting temperature, liquid and solid densities for the different interfaces. The uncertainty of the melting temperature is taken as the standard deviation with a confidence interval of 95% over the four estimations of the melting temperature.

	Bulk 1	5488	14 × 14 × 14
	Bulk 2	16000	20 × 20 × 20
	Bulk 3	1024000	80 × 80 × 80
	Interface 1 (100) 128000	640 × 10 × 10 2119 ± 15	0.0371 ± 0.03	0.0386 ± 0.01
	Interface 2 (100) 128000	80 × 40 × 20 2111 ± 15	0.0377 ± 0.01	0.0388 ± 0.01
	Interface 1 (110) 127960	640 × 10 × 10 2132 ± 15	0.0373 ± 0.01	0.0386 ± 0.04
	Interface 1 (111) 128016	640 × 10 × 10 2114 ± 15	0.0377 ± 0.03	0.0381 ± 0.01

Acknowledgement

We acknowledge the CINES and IDRIS under Project N INP2227/72914, as well as CI-MENT/GRICAD for computational resources. This work was performed within the framework of the Centre of Excellence of Multifunctional Architectured Materials "CEMAM" ANR-10-LABX-44-01 funded by the "Investments for the Future" Program. This work has been partially supported by MIAI @ Grenoble Alpes, (ANR-19-P3IA-0003). Fruitful discussions within the French collaborative network in high temperature thermodynamics GDR CNRS 3584 (TherMatHT) are also acknowledged.

low value of T rg = 0.34. Nevertheless, pure Ta still shows the ability to form a glass [7] with a cooling rate of the order of a 10 14 K/s. Interestingly, with a value of T rg = 0.42, higher than Ta, it is shown here that Zr can form a glass with an even lower cooling rate of 10 12 K/s.

Given the fact that the nose of the TTT curve, corresponding to the lowest nucleation time takes the value of 0.39 ns, this cooling rate allows clearly to avoid homogeneous nucleation.

Interestingly, from the nose position of the TTT curve of Bulk 3 system we obtain a critical cooling rate of 2.2 × 10 12 K/s, consistent with the fact that Zr crystallize during cooling with a cooling rate of 10 11 K/s. The first experimental attempt to form a metallic glass for Zr was done by Zhang and Zhao [36], but they retracted their paper later since what they interpreted as amorphous was in fact rapid crystal growth, and was later confirmed [37]. In their experiments, the cooling rate surely lower than the one used here might lead inevitably to crystal nucleation even if they were conducted under high pressure. Our results might also explain why Kim et al. [5,6] were not able to form an amorphous with Zr, as their cooling rate were limited to 10 7 K/s. Nevertheless, with the possibility to reach experimentally very high cooling rates [7], the results presented here indicated that obtaining Zr metallic glass may be reachable


Michel Raous
email: raous@lma.cnrs-mrs.fr

FRICTION AND ADHESION

Keywords: Unilateral contact, friction, adhesion, non smooth dynamics

The studies carried out on adhesion by the group "Modeling in Contact Mechanics" at the LMA are reviewed in this paper and recent applications are presented. Based on the introduction of the adhesion intensity variable developed by M. Frémond, different forms of a model coupling adhesion to unilateral contact and friction have been developed. The formulations are given either under the form of implicit variational inequalities or the one of complementarity problems. Both quasi-static and dynamic formulations are considered. The model is non smooth because we do not use any regularization for the unilateral conditions and for the friction, i.e. Signorini conditions and strict Coulomb law are written. In the thermodynamics analysis, the state and the complementarity laws are then written using subdifferentials and differential inclusions because of the non convexity and non differentiability of the potentials. For the dynamics, the formulation is given in term of differential measures in order to deal with the non continuity of the velocities that may occur in the solutions. This work therefore owes much to the theories and the numerical scheme developed by J.-J. Moreau and M.

Introduction

In order to describe the smooth transition from a completely adhesive contact to a usual unilateral contact (Signorini conditions) with Coulomb friction, a model based on interface damage has been first developed for quasi-static problems in [START_REF] Raous | Un modèle couplant adhérence et frottement pour le contact unilatéral entre deux solides déformables[END_REF] [START_REF] Cangémi | Frottement et adhérence : modèle, traitement numérique et application à l'interface fibre/matrice[END_REF] 1 [START_REF] Raous | Consistent model coupling adhesion, friction and unilateral contact[END_REF] [START_REF] Raous | Quasistatic Signorini problem with Coulomb friction and coupling to adhesion[END_REF]. Using a dynamic formulation, the model was then extended to account for the brittle behaviour occurring when a crack interacts with fiber-matrix interfaces in composite materials in (Raous-Monerie, 2002) [START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF]. More recently, the model has been used to study metal/concrete interfaces in reinforced concrete in civil engineering (Karray et al, submitted), delamination of coated bodies (Raous et al, 2002), delamination of glued assembling in civil engineering [START_REF] Raous | Modélisation de l'adhésion par collage[END_REF], cohesive masonry [START_REF] Jean | Non-smooth contact dynamics approach of cohesive materials[END_REF]) [START_REF] Acary | Contribution la modélisation mécanique et numérique des édifices maçonnés[END_REF] and production of wear particles in bio-engineering [START_REF] Baudriller | Approche micro-mécanique et modélisation de l'usure : application aux biomatéraiux des prothéses totales de hanche[END_REF]. The quasi-static formulation was extended to deal with hyperelasticity in [START_REF] Bretelle | Unilateral contact with adhesion and friction between two hyperelastic bodies[END_REF]. Mathematical results about the existence of the solutions were given in [START_REF] Cocou | Existence results for unilateral quasistatic contact problems with friction and adhesion[END_REF] without using any regularization on the contact conditions.

THE MODEL

The RCC model (Raous-Cangémi-Cocou) has been first given in [START_REF] Raous | Un modèle couplant adhérence et frottement pour le contact unilatéral entre deux solides déformables[END_REF] [START_REF] Cangémi | Frottement et adhérence : modèle, traitement numérique et application à l'interface fibre/matrice[END_REF] and then extensively presented in [START_REF] Raous | Quasistatic Signorini problem with Coulomb friction and coupling to adhesion[END_REF]. It has been extended to the present form including progressive friction with the term (1 -β) in [START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF]) (Raous-Monerie, 2002).

Adhesion is characterized in this model by the internal variable β, introduced by Frémond [START_REF] Frémond | Adhérence des solides[END_REF] [START_REF] Frémond | Contact with adhesion[END_REF], which denotes the intensity of adhesion. The introduction of a damageable stiffness of the interface ensures a good continuity between the two models during the competition between friction and adhesion. The behaviour of the interface is described by the following relations, where (1) gives the unilateral contact with adhesion, (2) gives the Coulomb friction with adhesion and

(3) gives the evolution of the adhesion intensity β. Initially, when the adhesion is complete, the interface is elastic as long as the energy threshold w is not reached. After that, damage of the interface occurs and consequently, on the one hand, the adhesion intensity β and the apparent stiffness β 2 C N and β 2 C T decrease, and on the other hand, friction begins to operate. When the adhesion is completely broken (β = 0), we get the classical Signorini problem with Coulomb friction.

-R r N + β 2 C N u N ≥ 0 , u N ≥ 0 , -R r N + β 2 C N u N . u N = 0, (1) R r T = β 2 C T u T , R r N = R N , R T -R r T ≤ µ(1 -β) R N -β 2 C N u N , (2)
with if R T -R r T < µ(1 -β) R N -β 2 C N u N ⇒ uT = 0, if R T -R r T = µ(1 -β) R N -β 2 C N u N ⇒ ∃λ ≥ 0, uT = λ(R T -R r T ), β = -w -β (C N u 2 N + C T u T 2 ) -k∆β - b 1/p if β ∈ [0, 1[, β ≤ -w -β (C N u 2 N + C T u T 2 ) -k∆β - b 1/p if β = 1 .
(3) R N and u N are the algebraic values of the normal components of the contact force and those of the relative displacement between the two bodies (occupying domains Ω 1 and Ω 2 ) defined on the contact boundary Γ c , and R T and u T are the tangential components of this contact force and those of this relative displacement. Subscript r denotes the reversible parts. The constitutive parameters of this interface law are as follows :

C N and C T are the initial stiffnesses of the interface (full adhesion), µ is the friction coefficient, w is the decohesion energy, p is a power coefficient (p = 1 in what follows) and k = 0 in what follows.

= C T = C, u 0 = ω/C and R 0 = √ ωC) .
It should be noted that the Signorini conditions are strictly imposed when compression occurs. References on other models can be found in [START_REF] Raous | Quasistatic Signorini problem with Coulomb friction and coupling to adhesion[END_REF] [START_REF] Raous | Consistent model coupling adhesion, friction and unilateral contact[END_REF]) and a comparison between some of them is made in [START_REF] Monerie | Comparaison de lois d'interface fibre/matrice sur la base d'un modèle uni-axial d'expérience de micro-indentation[END_REF] (models developed by Tvergaard-Needleman, Girard-Feyel-Chaboche, Michel-Suquet, Allix-Ladevèze). Using penalization and augmented Lagrangian on a similar model to the RCC one, Talon-Curnier have solved the quasistatic problem using generalized Newton method [START_REF] Talon | A model of adhesion coupled to contact and friction[END_REF].

THE THERMODYNAMICS

In the framework of continuum thermodynamics, the contact zone is assumed to be a material surface and the local constitutive laws are obtained by choosing two specific forms of the free energy and the dissipation potential associated to the surface. The following thermodynamic variables are introduced : the relative displacements (u N n,u T ) and the adhesion intensity β are chosen as the state variables, and the contact force R and a decohesion force G, as the associated thermodynamic forces. The thermodynamic analysis given for the RCC model in [START_REF] Raous | Consistent model coupling adhesion, friction and unilateral contact[END_REF] has been extended to the present model in [START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF]) (Raous-Monerie, 2002) in order to obtain relations (1) to (3), where friction is progressively introduced in the form of the term (1 -β)µ into (2). Expressions (4) and ( 5) are adopted for the free energy Ψ(u N , u T , β) and the potential of dissipation Φ uT , β . In (4), the indicator function I K (where K = {v / v ≥ 0} ) imposes the unilateral condition u N ≥ 0 and the indicator function I P (where P = {γ / 0 ≤ γ ≤ 1}) imposes the condition β ∈ [0, 1]. In (5), the indicator function I C -( β) (where C -= {γ / γ ≤ 0}) imposes that β ≤ 0 : the adhesion can only decrease and cannot be regenerated (it is irreversible) in the present model.

Ψ(u N , u T , β) = 1/2 β 2 C N u 2 N +1/2 β 2 C T u T 2 -w β+I K (u N )+I P (β) (4) Φ uT , β = µ(1-β) R N -β 2 C N u N uT + b/(p+1) β p+1 + I C -( β) (5)
Ψ has a part which is convex but not differentiable and another part which is differentiable but not convex with respect to the pair (u, β). Φ is convex but has a part which is not differentiable. The state laws and the complementarity laws are then written as follows in order to obtain the contact behaviour laws given in section 1 [START_REF] Raous | Consistent model coupling adhesion, friction and unilateral contact[END_REF].

           R r T = ∂Ψ d ∂[u T ] R r N - ∂Ψ d ∂[u N ] ∈ ∂I K ([u N ]) -G β - ∂Ψ d ∂β ∈ ∂I [0,1] (β) (6) R N = R r N (R ir T , G β ) ∈ ∂Φ( [u T ], β) (7)

THE QUASI-STATIC FORMULATION

The formulation and the approximation of quasi-static frictional problems given in [START_REF] Cocou | Formulation and approximation of quasistatic frictional contact[END_REF] have been extended to adhesion problems in [START_REF] Raous | Un modèle couplant adhérence et frottement pour le contact unilatéral entre deux solides déformables[END_REF] [START_REF] Raous | Consistent model coupling adhesion, friction and unilateral contact[END_REF] [START_REF] Raous | Quasistatic Signorini problem with Coulomb friction and coupling to adhesion[END_REF]. The problem can be here set as the coupling between two variational inequalities (one of which is implicit) and a differential equation.

Problem (P 1 ) : Find ( u, β) ∈ W 1,2 (0, T ; V ) × W 1,2 (0, T ; H) such that u(0) = u 0 ∈ K, β(0) = β 0 ∈ H [0, 1[ and for ∀t ∈ [0, T ], u(t) ∈ K, and ∀v ∈ V a( u, v -˙ u) + j(β, u N , v T ) -j(β, u N , uT ) + Γ C β 2 C T u T . (v T -uT )ds ≥ ( F , v -˙ u) -R N , v N -uN , ( 8
) -R N , z -u N + Γ C β 2 C N u N .(z -u N )ds ≥ 0 ∀z ∈ K, (9)
β = -1/b w -(C N u 2 N + C T u T 2 )β - a.e. on Γ C , (10)
where :

-u = (u 1 , u 2 ) where u 1 and u 2 define the displacements in Ω 1 and Ω 2 , -

V = (V 1 , V 2 ), V α = v α ∈ H 1 (Ω α ) 3 ; v α = 0 a.e. on Γ α U , α = 1, 2, -H = L ∞ (Γ c ), -K = v = (v 1 , v 2 ) ∈ V 1 × V 2 ; v N ≥ 0 a.e. on Γ c
, where Γ c is the contact boundary between the two solids Ω 1 and Ω 2 , -a(. , .) is the bilinear form associated to the elasticity mapping,

-j(β, u N , v T ) = Γ C µ(1 -β) R N -β 2 C N u N v T ds , -F = (F 1 , F 2 )
are the given force densities applied to solid 1 and to solid 2 respectively. By using an incremental approximation, it has been established in [START_REF] Raous | Un modèle couplant adhérence et frottement pour le contact unilatéral entre deux solides déformables[END_REF] [START_REF] Raous | Consistent model coupling adhesion, friction and unilateral contact[END_REF]) that the numerical solutions can be obtained by adapting the methods that we have developed for dealing with classical frictional unilateral contacts [START_REF] Raous | Quasistatic Signorini problem with Coulomb friction and coupling to adhesion[END_REF]. The solutions have mainly been obtained as follows :

-either by taking a fixed point on the sliding threshold and solving a sequence of minimization problems with the choice between a projected Gauss-Seidel method (accelerated by relaxation or Aitken processes) and a projected conjugate gradient method [START_REF] Raous | Preconditioned conjugate gradient method for a unilateral problem with friction[END_REF], -or by taking a complementarity formulation, using a mathematical programming method (Lemke's method). With the adhesion model, these solvers are coupled with the numerical integration of the differential equation on β by using θ-methods (θ = 1 is often chosen). Implementation of the algorithms has been conducted in our finite element code GYPTIS90 (Latil-Raous, 1991).

THE DYNAMICS

Depending on the characteristics of the interface, especially when b tends towards zero (no viscosity for the evolution of the intensity of adhesion), brittle behaviour can be obtained and the inertia effects have to be taken into account. A 3D dynamic formulation has been developed (Raous-Monerie, 2002) (Monerie-Raous, 1999) (Monerie-Acary, 2001). Because of the non smooth character of these interface laws, the dynamics is written in terms of differential measures as follows (where q denotes the displacement and r the contact force) :

M (q, t)d q = F (q, q, t)dt -rdν , ( 11
)
where d q is a differential measure associated with q(t) :

t 2 t 1 d q = q(t + 2 ) -q(t - 1 ) ∀t 2 > t 1 , (12)
and hence :

t 2 t 1 M (q, t)d q = t 2 t 1 F (q, q, t)dt - ]t 1 , t 2 ] rdν, (13)
q(t 2 ) = q(t 1 ) + t 2 t 1 q(τ )dτ . ( 14
)
The Non Smooth Contact Dynamics method developed by (Moreau, 1998) [START_REF] Moreau | Some numerical methods in multi-body dynamics: application to granular media[END_REF] [START_REF] Moreau | Some basics of unilateral dynamics[END_REF] [START_REF] Jean | The Non Smooth Contact Dynamics method[END_REF] has been extended to the treatment of the RCC model. A solver for dealing with adhesion has been implemented in the finite element code LMGC [START_REF] Jean | The Non Smooth Contact Dynamics method[END_REF]) (Monerie-Acary, 2001). Another solver based on complementarity formulation (Lemke's method) dealing with non smooth dynamics problems has been implemented in the finite element code SIMEM3 at the LMA [START_REF] Vola | Consistent time discretization for dynamical frictional contact problems and complementarity techniques[END_REF].

APPLICATIONS Delamination benchmarks

Various benchmarks for simulating delamination have been developed in the framework of a joint project with the Laboratoire Central des Ponts et Chaussées [START_REF] Raous | Modélisation de l'adhésion par collage[END_REF] focusing on adhesion and gluing in civil engineering. Micro-indentation of a single fiber in a composite material This model was first developed in order to describe the behaviour of the fiber-matrix interface of composite materials [START_REF] Cangémi | Frottement et adhérence : modèle, traitement numérique et application à l'interface fibre/matrice[END_REF]) [START_REF] Raous | Consistent model coupling adhesion, friction and unilateral contact[END_REF]. The parameters of the model were identified on micro-indentation experiments performed at the ONERA on single fibers. Physical and mechanical considerations are taken into account in the identification procedure : for example, the values of the initial interface stiffness C N , C T are taken to be in the range corresponding to the elastic properties of the oxides located at the interfaces. The viscosity parameter b is particularly difficult to identify : experiments with different loading velocities are required for this purpose. The validity of the model was then confirmed by taking various kinds of loadings (especially cyclic loadings) and comparing experimental results with those obtained in the numerical simulation whith the same interface parameters. Details of these studies can be found in [START_REF] Cangémi | Frottement et adhérence : modèle, traitement numérique et application à l'interface fibre/matrice[END_REF] [START_REF] Raous | Consistent model coupling adhesion, friction and unilateral contact[END_REF] and later in [START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF].

Interaction between cracks and fiber/matrix interfaces in composite materials

Again in collaboration with ONERA, during the thesis [START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF], the RCC model has been used to investigate the different ways in which cracks propagate through a composite material and how they depend on the interface properties. The adhesion model has been used to account both the crack propagation (decohesive crack with no viscosity, i.e. b = 0) and the interface behaviour. Crack bridging, crack trapping and fiber breaking can be observed depending on the interface characteristics (Raous-Monerie, 2002) [START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF]. On Fig. 1.5 adhesion is broken in the black zones.

Metal/concrete interfaces in reinforced concrete

In civil engineering, we are now testing the RCC model for the adhesive contact between steel and concrete in reinforced concrete. Pull-out tests of a steel shaft are being simulated in the framework of a joint project between LMA and ENIT in Tunisia (Karray et al, submitted). In that case, a variable friction coefficient (depending on the sliding displacement) has been used in order to take into account the wear of the surface which occurs during the sliding which seems to be quite significant with concrete. Another version of the way used to introduce friction has been used : (1 -β 2 ) instead of (1 -β). Results with C T = C N =16N/mm 3 , µ varying from 0.45 to 0.3 and f(β) = (1 -β 2 ) are given in Fig. 1.6.

Delamination of a coated body

A simplified approach to the delamination of a coated body has been conducted to simulate the indentation of a thin Chrome layer adhering to a metal body (Raous et al, 2002).

Cohesive masonry

The RCC model has been used in [START_REF] Jean | Non-smooth contact dynamics approach of cohesive materials[END_REF]) [START_REF] Acary | Contribution la modélisation mécanique et numérique des édifices maçonnés[END_REF] to simulate the behaviour of a cohesive dome. In Fig. 1.8, the deformation of a dome let on pillars and submitted to gravity is shown.

Wear in biomechanics

In order to simulate the wear occurring in biomaterials and especially the production of wear particles, H. Baudriller [START_REF] Baudriller | Approche micro-mécanique et modélisation de l'usure : application aux biomatéraiux des prothéses totales de hanche[END_REF] has used the RCC model in a volumic sense. The adhesive model is used on all the interfaces between the finite elements.

Figure 1 . 2 .

 12 Figure 1.1. Normal behaviour of the interface

  Figure 1.3. Delamination of two layers submitted to vertical traction

Figure 1 . 5 .

 15 Figure 1.5. Interaction of a crack and a fiber/matrix interface

Figure 1 . 6 .

 16 Figure 1.6. Simulation of a pull-out test on reinforced concrete

Figure 1 . 7 .

 17 Figure 1.7. Indentation of a coated surface

Figure 1 . 8 .

 18 Figure 1.8. Behaviour of a cohesive dome

  Fig 1.9 gives an example of production of fragments.

Figure 1 . 9 .

 19 Figure 1.9. Fragmentation of a block

Acknowledgments

M. Cocou, J.-C. Latil, I. Rosu, in the group MMC at the LMA, L. Cangémi, Y. Monerie, during the preparation of their thesis, and M.A. Karray (ENIT) have contributed to the works presented here. Thanks to V. Acary and H. Baudrier for providing figures 1.8 and 1.9.