Katharina Immel

Thang X Duong

Vu-Hieu Nguyen

Guillaume Haïat
email: guillaume.haiat@univ-paris-est.fr

Roger A Sauer

A modified Coulomb's law for the tangential debonding of osseointegrated implants

Keywords: sticking and sliding friction, adhesion, computational contact mechanics, boneimplant interface, state variable friction, finite element method

. Overall, the results show close agreement with the experimental data, especially the peak and the softening part of the torque curve with a relative error of less than 2.25 %. In addition, better estimates of the bone's shear modulus and the adhesion energy are obtained. The proposed model is particularly suitable to account for partial osseointegration, as is also shown.

Introduction

Endosseous implants are widely used in orthopedic and dental surgery, for instance in total hip and knee arthoplasty or tooth replacement. Implant stability is a strong determinant of implant success, and it is established through bone growing around and into the rough surface of the implant, a process called osseointegration [START_REF] Albrektsson | Osseointegrated titanium implants: requirements for ensuring a long-lasting, direct bone-to-implant anchorage in man[END_REF]. However, debonding-related implant failures still occur and are difficult to anticipate, as bone remodeling and osseointegration phenomena are complex and remain poorly understood. Although these phenomena include biological, biochemical, and biomechanical factors, the biomechanical properties of the bone-implant interface (BII) are the governing factors of implant stability [START_REF] Haïat | Effects of biomechanical properties of the boneimplant interface on dental implant stability: from in silico approaches to the patient's mouth[END_REF][START_REF] Mathieu | Ultrasonic evaluation of dental implant biomechanical stability: an in vitro study[END_REF][START_REF] Mathieu | Numerical simulation of ultrasonic wave propagation for the evaluation of dental implant biomechanical stability[END_REF][START_REF] Mathieu | Biomechanical determinants of the stability of dental implants: Influence of the boneimplant interface properties[END_REF].

Two kinds of implant stability can be distinguished: primary (or initial) stability and secondary (or long term) stability (see for example Refs. [START_REF] Kim | Removal torque values of surface-treated mini-implants after loading[END_REF][START_REF] Sul | Oxidized implants and their influence on the bone response[END_REF]). Primary stability is achieved during surgery and is mainly governed by interlocking phenomena and bone quality, whereas secondary stability is achieved several weeks or months after surgery, through the formation and maturation of newly formed bone tissue at the BII. While the evolution of secondary implant stability is governed by complex biochemical processes, the mechanical behavior of the BII remains crucial for the surgical outcome [START_REF] Gao | Biomechanical behaviours of the boneimplant interface: a review[END_REF]. The difficulty of predicting stability arises from the complex nature of the BII, related to i) the implant surface roughness, ii) the non-homogeneous contact between bone and implant, iii) adhesion and friction phenomena, iv) the time-dependent change of periprosthetic bone properties [START_REF] Mathieu | Biomechanical determinants of the stability of dental implants: Influence of the boneimplant interface properties[END_REF] and v) the widely varying load cases during the implant life cycle. While the implant material, roughness and surface coating are important factors determining implant stability, friction phenomena also play a major role, and in turn depend on the surface roughness and bone quality [START_REF] Shirazi-Adl | Experimental determination of friction characteristics at the trabecular bone and porous-coated metal interface in cementless implants[END_REF].

Most studies on bone attachment to implants have used push-in or pull-out in vitro tests (see for example Refs. [START_REF] Berahmani | An experimental study to investigate biomechanical aspects of the initial stability of press-fit implants[END_REF][START_REF] Bishop | The influence of bone damage on press-fit mechanics[END_REF][START_REF] Damm | Friction coefficient and effective interference at the implant-bone interface[END_REF][START_REF] Wennerberg | Nanostructures and hydrophilicity influence osseointegration: a biomechanical study in the rabbit tibia[END_REF]). As the implant geometry influences the test results [START_REF] Brånemark | Biomechanical characterization of osseointegration: an experimental in vivo investigation in the beagle dog[END_REF] and leads to spatially complex, non-uniform, multiaxial stress fields [START_REF] Shirazi-Adl | Finite element stress analysis of a push-out test part 1: fixed interface using stress compatible elements[END_REF] and instantaneous crack propagation, using realistic implant geometries makes it difficult to estimate a physically meaningful value for the interfacial mechanical strength. Therefore, models with a planar BII were designed to minimize the effects of friction and mechanical forces introduced by the geometry [START_REF] Skripitz | Attachment of PMMA cement to bone: force measurements in rats[END_REF], see e.g. the experiments on coin-shaped implants by Rønold et al. [START_REF] Rønold | The use of a coin shaped implant for direct in situ measurement of attachment strength for osseointegrating biomaterial surfaces[END_REF][START_REF] Rønold | Analysing the optimal value for titanium implant roughness in bone attachment using a tensile test[END_REF]. However, the measured pull-out force in these experiments cannot be used to retrieve information about the strength of the interface. More recently, mode III cleavage experiments applied to coin-shaped implants have been proposed by Mathieu et al. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF]. A rotation of the implant with respect to the bone was imposed and the resulting moment was recorded. This resulted in stable crack propagation, which allows to assess the adhesion energy. However, the agreement between the analytical and experimental results was only moderate and significant discrepancies were obtained.

The tangential load-displacement behavior at the BII was found to be nonlinear [START_REF] Dammak | Analysis of cementless implants using interface nonlinear friction-experimental and finite element studies[END_REF][START_REF] Rancourt | Friction properties of the interface between porous-surfaced metals and tibial cancellous bone[END_REF][START_REF] Shirazi-Adl | Experimental determination of friction characteristics at the trabecular bone and porous-coated metal interface in cementless implants[END_REF]. However, numerical analyses of implant stability still only consider either fully bonded, frictionless contact or pure Coulomb's friction (see e.g. Refs. [START_REF] Chong | Analysis of bone-prosthesis interface micromotion for cementless tibial prosthesis fixation and the influence of loading conditions[END_REF][START_REF] Ghosh | Bone remodelling around cementless composite acetabular components: The effects of implant geometry and implant-bone interfacial conditions[END_REF][START_REF] Pettersen | Subject specific finite element analysis of implant stability for a cementless femoral stem[END_REF]). While such assumptions may be valid for the assessment of initial stability, prediction of failure of osseointegrated implants requires a more accurate description of the contact behavior and the inclusion of frictional and adhesive effects [START_REF] Castellani | Bone-implant interface strength and osseointegration: Biodegradable magnesium alloy versus standard titanium control[END_REF][START_REF] Tschegg | Characterization methods of bone-implantinterfaces of bioresorbable and titanium implants by fracture mechanical means[END_REF].

A simple way to model macroscopic friction phenomena that are related to different states of the BII are state variable friction laws, introduced by Rice and Ruina [START_REF] Rice | Stability of steady frictional slipping[END_REF][START_REF] Ruina | Slip instability and state variable friction laws[END_REF], which were motivated by the experimental findings of Dieterich [START_REF] Dieterich | Time-dependent friction and the mechanics of stick-slip[END_REF][START_REF] Dieterich | Modeling of rock friction: 1. experimental results and constitutive equations[END_REF]. These laws focus on the observed phenomena of i) fading memory and steady state, ii) positive dependence on the instantaneous slip rate, and iii) negative dependence on past slip rates. In general, it is assumed that at any given time, the contact surface has a certain state and the frictional stress only depends on that state, the slip rate and the contact pressure. Similarly, the rate of change for the state only depends on the current state, the slip rate, and the pressure at the analyzed point. Although these laws have so far been mainly applied in geology and geophysics, one can also interpret the state variable as the osseointegration degree of the BII.

The aim of this work is to propose a phenomenological model for the frictional contact behavior of debonding osseointegrated implants. The classical Coulomb's law is extended from a constant to a varying friction coefficient, that models the transition from an unbroken to a broken state, based on a state variable depending on the total sliding distance of the implant. While the unbroken state denotes osseointegration and thus the presence of adhesive bonds and a higher friction coefficient, the broken state denotes pure frictional contact behavior of the interface with a lower friction coefficient. Thus, this model can account for the higher tangential forces observed in osseointegrated implants compared to unbonded implants.

The remainder of this work is structured as follows: Sec. 2 describes the governing equations and the contact formulation, including the new friction model. In Sec. 3, the setup of the mode III cleavage experiment of osseointegrated implant is presented. Additionally, the corresponding analytical model developed by Mathieu et al. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] and the analytical and numerical model based on the proposed friction model of Sec. 2.2.3 are introduced. The new analytical and numerical models are calibrated with experimental data and compared to the analytical model from Ref. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] in Sec. 4. Finally, Sec. 5 gives a conclusion and an outlook to possible extensions.

Models and Methods

This section discusses the governing equations describing the contact behavior, as well as the new friction model. The resulting equations are discretized within a finite element framework to obtain a numerical solution, which is briefly summarized in Appendix A. The readers are referred to Duong and Sauer [START_REF] Duong | A concise frictional contact formulation based on surface potentials and isogeometric discretization[END_REF] and Sauer and De Lorenzis [START_REF] Sauer | A computational contact formulation based on surface potentials[END_REF][START_REF] Sauer | An unbiased computational contact formulation for 3d friction[END_REF] for a more detailed derivation of the considered contact formulation.

Material laws

Throughout this work, a hyper-elastic Neo-Hookean material model for both, bone and implant, is used. According to Zienkiewicz and Taylor [START_REF] Zienkiewicz | The finite element method for solid and structural mechanics[END_REF] the stress-strain relation for the Cauchy stress σ of this model is given by

σ = Λ J (lnJ) I + G J (b -I) , (1)
with the volume change J, the identity tensor I, and the left Cauchy-Green tensor b. The Lam parameters G (shear modulus) and Λ can be expressed in terms of Young's modulus E and Poisson ration ν, by

G = E 2(1 + ν) and Λ = 2Gν 1 -2ν . (2)

Friction laws

Given two 3D bodies B 1 and B 2 and their contact surface ∂ c B, the contact traction t c can be decomposed into a normal and tangential component, i.e.

t c = t n -t t . (3)
The magnitude of the normal traction is equal to the normal pressure, i.e.

t n = p. (4)
For frictional contact, the tangential traction is determined by the behavior during sticking and sliding, and the distinction between these two states is based on a slip criterion of the form f s < 0, for sticking, = 0, for sliding.

(

It can be formulated as

f s = t t -t max t , t t ≤ t max t , (6)
where t max t > 0 is the maximal tangential contact traction. For sticking, the tangential traction is defined by the constraint that no relative tangential motion occurs, while for sliding, the tangential traction is defined by a sliding law.

Coulomb's law

For classical Coulomb's friction, the slip function f s can be defined as

f s = t t -µp, (7)
where µ denotes the (constant) friction coefficient. Then, the tangential sliding traction is given by

t t = -µp ġs ġs , (8)
with ġs denoting the relative sliding velocity between the two bodies in contact, see Appendix A.

State variable friction laws

In so called rate-state or state variable friction laws (see, e.g. Rice and Ruina [START_REF] Rice | Stability of steady frictional slipping[END_REF][START_REF] Ruina | Slip instability and state variable friction laws[END_REF]) it is assumed, that at any given time, any point x on the contact surface has a state φ = φ(x, t).

The tangential contact traction then depends in general on the contact pressure p, the sliding rate ġs = ġs , and the state variable φ, i.e.

t t = t t (p, ġs , φ). (9)
For a given x, the state φ is assumed to be a continuous function of ġs (t). Moreover, the rate of change of φ at x generally only depends on the pressure, the sliding rate and the instantaneous state of this point, i.e. φ = F(p, ġs , φ).

Since φ does not depend on the state at other points, Eq. ( 10) is a local evolution law. State variable friction laws are able to model a change in the frictional contact traction due to past states (referred to as memory [START_REF] Rice | Stability of steady frictional slipping[END_REF][START_REF] Ruina | Slip instability and state variable friction laws[END_REF]). They are also able to model an asymptotic approach to steady state sliding when ġs becomes constant [START_REF] Rice | Stability of steady frictional slipping[END_REF].

Modified Coulomb's friction law

In contrast to a constant friction coefficient as is used in the classical Coulomb's law [START_REF] Brånemark | Biomechanical characterization of osseointegration: an experimental in vivo investigation in the beagle dog[END_REF], we propose to model µ as a function of the scalar state variable φ, as

µ := µ(φ) = φ µ ub + (1 -φ) µ b , (11)
where µ ub and µ b are the friction coefficient for the unbroken (initial) and broken state, respectively, that are weighted according to a state variable. As this is a local friction model, where φ and thus µ can change pointwise, it allows for the description of locally varying bonding states, such as occur in crack propagation and partial osseointegration.

According to Eq. ( 11), the state variable φ determines whether a point is in an unbroken, partially broken or fully broken state. Here, φ is considered to depend on the accumulated sliding distance

g s = t 0 ġs dt ( 12
)
at a certain point x, according to the smooth function

1 φ(g s ) = φ 0 •        1 if gs as < 1, 1 2 -1 2 sin π 2bs gs as -b s -1 if 1 ≤ gs as ≤ 1 + 2b s , 0 if gs as > 1 + 2b s . (13)
This function was designed such that it captures the experimental behavior shown in Sec. 4.2 (see Fig. 4). Equation ( 13) depends on the initial bonding state φ 0 = [0, 1], where φ 0 (x) = 0 denotes no initial bonding and φ 0 (x) = 1 represents full initial bonding. 2 The bonding state variable indicates that for φ(x) > 0 tangential contact at x is governed by the proposed friction law [START_REF] Chong | Analysis of bone-prosthesis interface micromotion for cementless tibial prosthesis fixation and the influence of loading conditions[END_REF], while for φ(x) = 0 it is governed by classical Coulomb's law [START_REF] Brånemark | Biomechanical characterization of osseointegration: an experimental in vivo investigation in the beagle dog[END_REF] with µ = µ b . This definition results in three possible states for every point on the contact surface: fully bonded (φ = 1), debonding (0 < φ < 1) and fully debonded/sliding (φ = 0), which is illustrated in Figures 1a and1b. The parameter a s represents the sliding threshold up to which tangential adhesion takes effect (µ = µ ub for g s ≤ a s ), while b s defines the size of the transition zone between the two friction coefficients. This implies, that up to a sliding length of g s = a s we have a higher resistance to tangential displacement, similar to the effect of adhesion. After the sliding distance a s is reached, the friction coefficient starts to decrease to µ = µ b , corresponding to the sliding of a fully debonded body. The sliding distance needed for a certain point on the contact surface to fully debond is then controlled by the parameter b s . 0 0.1 0. For simple cases, this model can be solved analytically, as presented in Sec. 3.3. However, for the complex geometries of endosseous implants, the nonlinear and anisotropic behavior of bone tissue and loading conditions inside the human body may require a numerical solution. This becomes particularly important when considering inhomogeneous initial bonding (where φ 0 (x) = const ∀x), due to imperfect or partial osseointegration, as presented in Sec. 4.4.

Mode III debonding of osseointegrated coin-shaped implants

To calibrate the friction model described above, it is applied to the mode III debonding of an osseointegrated implant. The test case is originally presented by Mathieu et al. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] (in the following also referred to as the reference study). First, a short summary of the experimental and analytical results from the reference study are given. Second, a new analytical model based on the proposed friction model of Sec. 2.2.3 is introduced. Third, the numerical setup and parameters are presented.

Experimental setup

In Mathieu et al. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF], two coin-shaped implants made of titanium alloy (Ti-6Al-4V), with a radius R = 2.5 mm and a height H = 3 mm, were implanted into the tibiae of a rabbit and left in vivo during seven weeks. Polytetrafluoroethylene (PTFE) caps were placed around the implants, to ensure that bone in-growth only occurred at the bottom of the cylindrical implants. After seven weeks, the rabbit was sacrificed and the bone samples including the implants were extracted and conserved. Then, mode III cleavage experiments were carried out. The bone specimen was rigidly fixed to minimize the remaining normal force. The implant was fixed by a chuck screwed to a torque meter. Then, a 10 • rotation with a constant rotation speed of 0.01 • s -1 was imposed. Finally, the torque and rotation angle as a function of time were extracted via post processing.

Analytical model of Mathieu et al. [32]

To explain the experimental results, Mathieu et al. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] developed an analytical model that couples crack propagation with friction, based on Chateauminois et al. [START_REF] Chateauminois | Friction and shear fracture of an adhesive contact under torsion[END_REF]. The model assumes that the crack starts from the outside and propagates purely in a circular manner and in mode III to the center of the contact surface, using a Dugdale crack model, without any normal force. In this description, the interfacial forces are supposed to be constant up to a given separation distance between the surfaces. Let R be the radius of the implant and the initial contact area, then c defines the radius of the current crack (see Fig. 2a), corresponding to the twisting angle

θ = π wadh 4cG + τ 0 2G cosh -1 R c , (14)
where wadh denotes the adhesion energy. c separates the contact area into an adhering/sticking region for r < c and a debonded/sliding region for c < r < R. In the sliding region, the orthoradial shear stress σ θz is assumed to be constant, according to

σ θz (r) = τ 0 for c < r < R, (15)
with constant τ 0 = 3T ∞ /2πR 3 , where T ∞ corresponds to the torque for a rotation angle equal to infinity (i.e., where the surface is fully debonded), and z denotes the axial direction. In the sticking region, the orthoradial stress becomes

σ θz (r) = 2 π   πG wadh c r √ c 2 -r 2 + τ 0 sin -1   R c 2 -1 R r 2 -1     for r < c. (16)

New analytical model with modified Coulomb's friction

Given the modified Coulomb's law in Sec. 2.2.3, a new analytical model can be derived. Here, only the key aspects of the new analytical model are presented. The complete derivation can be found in Appendix B.

Due to symmetry, the tangential traction component σ θz is distributed radially symmetric along the BII, while the radial traction component σ zr is zero (as in the model presented in Sec. 3.2). We assume both bodies to be linear elastic and the normal contact pressure p to be distributed homogeneously along the contact surface. We define the hyperbola

c(θ) := θ lin θ R (17)
to be the function of the critical radius for the stick/slip transition, such that r < c denotes the sticking region and c ≥ r ≥ R denotes the sliding region. Furthermore, we define θ lin as the limit for which the tangential traction will be a linear function of the implant radius and θ max is the rotation angle for which the whole contact surface starts sliding (see Fig. 12). It also marks the location of the maximal torque T max in the torque-angle-curve (see e.g. Fig. 4a).

The tangential contact traction as a function of the rotation angle can then be expressed as

σ θz (θ, r) = t max t r c(θ) , for r < c (sticking), µp, for r ≥ c (sliding), (18)
with µ = µ(g s ) from Eq. [START_REF] Chong | Analysis of bone-prosthesis interface micromotion for cementless tibial prosthesis fixation and the influence of loading conditions[END_REF]. Fig. 3 illustrates the variation of the tangential contact traction [START_REF] Duong | A concise frictional contact formulation based on surface potentials and isogeometric discretization[END_REF] for θ = θ lin , θ = θ max , and some θ larger than θ max .

While the two analytical models presented in Sec. 3.2 and Sec. 3.3 both assume circular crack propagation from outside to inside and determine the stick/slip transition by a critical (or crack) radius, the analytical model from Ref. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] imposes the crack radius and computes the corresponding rotation angle, while the new model directly imposes the rotation angle. 18)) at different rotation angles.

Numerical setup

As in the experiments, we consider a coin-shaped cylindrical implant with dimensions R = 2.5 mm and H = 3 mm. The bone sample is modeled as a rectangular cuboid with dimensions 12.5 × 12.5 × 5 mm. The implant is positioned at the center of the upper bone surface. The bodies are meshed according to the parameters given in Tab. 1, where n e denotes the number of elements of the body/surface and n gp denotes the number of Gauss-points per element. While the bulk is discretized with linear Lagrangian shape functions, the contact surfaces are discretized with non-uniform rational B-Splines (NURBS) (see Appendix A.3). The finite element mesh is pictured in Fig. 2b. To justify this coarse discretization, a refinement analysis of the mesh and the load step size is presented in Appendix C.

Contact is computed with a penalty regularization, and the corresponding penalty parameter is set according to the Young's modulus of bone to = E b /L 0 , with L 0 = 0.01 m.

The lower surface of the bone block is fixed in all directions, while the sides of the bone block are fixed in their corresponding normal direction. Our modified Coulomb's friction model requires a contact pressure during sliding, see Eq. ( 8) and [START_REF] Duong | A concise frictional contact formulation based on surface potentials and isogeometric discretization[END_REF]. We generate this pressure by applying a uniform vertical displacement d at the upper surface of the implant. Then, the implant is rotated by 10 • around its central axis with a constant load step size of 0.01 For both bodies, the Neo-Hookean material model of Eq. ( 1) is used. The material properties for the implant are those of titanium alloy (Ti-6Al-4V: E i = 113 GPa, ν i = 0.3). The material properties of the bone and the friction parameters µ ub , µ b , a s , b s , and the contact pressure p have to be determined by a parameter study and are presented in Sec. 4.1. All simulations were performed with an in-house, MATLAB-based solver3 .

Results

In the following, the results obtained with the new analytical model and numerical study are presented and compared to the experimental and analytical results corresponding to the reference study. First, the parameter estimation and the subsequent error estimation for homogeneous osseointegration are presented. Second, the torque-per-angle curves of the different models are compared and the debonding behavior of the implant is discussed. Third, the work of adhesion and frictional energy loss of the models are compared. Last, several cases of partial osseointegration are presented and compared with the homogeneous case.

Parameter calibration

During the parameter estimation stage, the Poisson ratio of bone is fixed to ν b = 0.3. The remaining parameters G b , d, µ ub , a s , b s are determined by minimizing the mean relative error

e mp T = mean θ∈[0 • ,10 • ] T exp (θ) -T (θ) T exp (θ) , (19)
where T exp is defined as the torque over the rotation angle θ obtained from the corresponding experiment. The shear modulus G b is calibrated using the initial slope of the linear part of the torque-per-angle curve (i.e. T (θ ≤ θ lin )) (e.g. see Fig. 4a). The other parameters depend on the friction coefficient µ b , which is investigated at the fixed values µ b = [0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5], which corresponds to the typical range found in the literature [START_REF] Biemond | Frictional and bone ingrowth properties of engineered surface topographies produced by electron beam technology[END_REF][START_REF] Grant | Artificial composite bone as a model of human trabecular bone: the implant-bone interface[END_REF][START_REF] Rancourt | Friction properties of the interface between porous-surfaced metals and tibial cancellous bone[END_REF][START_REF] Shirazi-Adl | Experimental determination of friction characteristics at the trabecular bone and porous-coated metal interface in cementless implants[END_REF]. The vertical displacement d is calibrated from steady state sliding, which is considered here to occur for θ > 9 • . We thus minimize Eq. ( 19) for θ ∈ [9 • , 10 • ] to calibrate d. There, the implant is assumed to be fully debonded and thus, tangential sliding contact is governed by Coulomb's law [START_REF] Brånemark | Biomechanical characterization of osseointegration: an experimental in vivo investigation in the beagle dog[END_REF] with µ = µ b . Once d is determined, it is considered constant for all rotation angles θ.

The initial friction coefficient µ ub is calibrated using T max (as t max t = µ ub p) (see Tab. 2). Finally, the parameters a s and b s are then determined by minimizing Eq. ( 19) for the whole torque-per-angle curve.

The bone shear modulus can be estimated from the new analytical model. Due to the nonhomogeneous pressure distribution, the displacement d and subsequently all other parameters are determined using the numerical model. However, adjusting d for the new analytical model of Sec. 3.3 leads to the same estimated values for µ ub , a s , and b s .

Tab. 3 shows the estimated shear moduli and the corresponding Young's moduli obtained from the numerical parameter estimation for the two data sets, compared to the computed shear moduli obtained in Ref. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF]. The estimated shear moduli of 7 and 8 GPa are higher than the reported values of 2-6 GPa [START_REF] Sharma | Comparative study of locational variation in shear and transverse elastic modulus of buffalo cortical bone[END_REF][START_REF] Tang | Shear deformation and fracture of human cortical bone[END_REF]), while the corresponding Young's moduli of 18 and 4a and 4b of Mathieu et al. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF], respectively. θ lin denotes the limit of the elastic part of the deformation, while θ max is the location of the maximum torque T max . T ∞ denotes the torque for a fully debonded implant and is taken at 10 • .

data set T max [N m] T ∞ [N m] θ lin θ max 1 0
21 GPa are in good agreement with experimental data from the literature [START_REF] Novitskaya | Recent advances on the measurement and calculation of the elastic moduli of cortical and trabecular bone: a review[END_REF][START_REF] Rho | Young's modulus of trabecular and cortical bone material: ultrasonic and microtensile measurements[END_REF] and previous studies [START_REF] Vayron | Variation of biomechanical properties of newly formed bone tissue determined by nanoindentation as a function of healing time[END_REF][START_REF] Vayron | Nanoindentation measurements of biomechanical properties in mature and newly formed bone tissue surrounding an implant[END_REF]. In Ref. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] the model was fit to match the peak and the decrease in torsion, which results in a considerable error for the shear modulus, as shown in Tab. 3 and Figs. [START_REF] Bishop | The influence of bone damage on press-fit mechanics[END_REF].

As our proposed model allows for more control, it is possible to match more characteristics of the torque curve, such as the initial slope, the peak, the softening and the torque for complete debonding.

Tab. 4 shows the results of the numerical parameter estimation of the different data sets and the corresponding mean percentage error, for the chosen levels of friction coefficient µ b . For the analyzed µ b , the corresponding initial friction coefficient µ ub lies between 0.29 and 0.73 and agrees well with the values of 0.28-1.1 reported in the literature [START_REF] Biemond | Frictional and bone ingrowth properties of engineered surface topographies produced by electron beam technology[END_REF][START_REF] Damm | Friction coefficient and effective interference at the implant-bone interface[END_REF][START_REF] Grant | Artificial composite bone as a model of human trabecular bone: the implant-bone interface[END_REF][START_REF] Rancourt | Friction properties of the interface between porous-surfaced metals and tibial cancellous bone[END_REF][START_REF] Shirazi-Adl | Experimental determination of friction characteristics at the trabecular bone and porous-coated metal interface in cementless implants[END_REF][START_REF] Zhang | Interfacial frictional behavior: cancellous bone, cortical bone, and a novel porous tantalum biomaterial[END_REF]. [START_REF] Novitskaya | Recent advances on the measurement and calculation of the elastic moduli of cortical and trabecular bone: a review[END_REF][START_REF] Rho | Young's modulus of trabecular and cortical bone material: ultrasonic and microtensile measurements[END_REF][START_REF] Sharma | Comparative study of locational variation in shear and transverse elastic modulus of buffalo cortical bone[END_REF] Table 3: Results for the parameter estimation of the two data sets for values independent of the friction coefficient µ b . Shown are results from the analytical model presented in [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] (see Eqs. ( 15),( 16)) and the present numerical solution. The estimated parameter is the bone shear modulus G b . The Young's modulus E b then follows from Eq. ( 2).

Increasing data set model G b [GPa] E b [GPa] literature
µ b results in a smaller normal force and a higher µ ub . Calibrating the parameters using T max and T ∞ results in almost identical curves for all tested values of µ b . Therefore, for the second data set only µ b = 0.3 [START_REF] Rancourt | Friction properties of the interface between porous-surfaced metals and tibial cancellous bone[END_REF] and µ b = 0.4 [START_REF] Grant | Artificial composite bone as a model of human trabecular bone: the implant-bone interface[END_REF][START_REF] Shirazi-Adl | Experimental determination of friction characteristics at the trabecular bone and porous-coated metal interface in cementless implants[END_REF] are investigated, which are the values most commonly reported for the interface between cortical bone and polished metal implants.

As shown in Tab. 4, the parameters a s and b s only depend on the shear modulus and on the overall shape of the torque curve, i.e., the width of the peak and the slope of T (θ). The computed values for a s indicate that adhesion takes effect for a sliding distance up to 22 and 26 µm, respectively. When these values are exceeded, the implant starts to debond, which is indicated by a decreasing friction coefficient (see Fig. 6b). This observation is in accordance with the reported threshold for micro-motion of the BII, where no deformation occurs. In most studies, a value of up to 50 µm is reported [START_REF] Bragdon | Differences in stiffness of the interface between a cementless porous implant and cancellous bone in vivo in dogs due to varying amounts of implant motion[END_REF][START_REF] Fitzpatrick | Computationally efficient prediction of boneimplant interface micromotion of a cementless tibial tray during gait[END_REF], while values exceeding 150 µm have shown to inhibit bone growth and promote bone loss [START_REF] Jasty | In vivo skeletal responses to porous-surfaced implants subjected to small induced motions[END_REF][START_REF] Pilliar | Observations on the effect of movement on bone ingrowth into porous-surfaced implants[END_REF]. However, these values were only reported for normal displacement and may vary for tangential displacement.

data set 15),( 16)) and the present numerical solution. The estimated parameters are the enforced normal displacement d, corresponding average contact pressure p, friction coefficient for the unbroken state µ ub , sliding threshold a s , and transition factor b s .

µ b d [µm] p [MPa] µ ub a s [µm] b s e mp T [

Torque curves and debonding behavior

The corresponding curves representing the variation of the torque as a function of the angle of rotation are shown in Figs. 4a and4b. The resulting torque obtained with the new analytical model described in Sec. 3.3 and the numerical solution of the proposed friction model of Sec. 2.2.3 (using µ b = 0.4) are very close. While the analytical solution shows a slightly closer fit to the experimental data before the peak in torque (resulting in an error of 2.18 % for the first data set and 2.83 % for the second data set), the numerical solution provides a better estimation of the behavior after debonding.

The errors obtained with the numerical solution are given in Tab. 4 and are compared to the analytical model from the reference study. Overall, the numerical solutions yield the best agreement with the experimental data, especially concerning the initial slope (i.e. stiffness) of the torque and the decrease after its peak. The torque curves show a flat plateau at the peak, which comes from the behavior induced by the sliding threshold a s . Increasing a s induces an offset of the debonding process and thus, results in an elongated peak. Another difference is shown in Fig. 4b for θ > 2.5 • , where the decrease of the torque is not exactly reproduced. A different transition function φ may allow for a closer fit there.

The top row in Fig. 5 shows the distribution of the friction coefficient µ in the contact area for different angles of rotation. A transition zone (characterized by µ b < µ < µ ub ), which may be understood as a crack front, cf. Ref. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] and Fig. 2a, appears at 1.13 • . The transition zone propagates inward in the radial direction from the external radius R into the center, which corresponds to the crack mode assumed by the analytical models. The transition zone can also be observed in Fig. 6a, which shows the value of the friction coefficient as a function of the implant radius for different angles of rotation. When the rotation angle increases, the width of the transition zone decreases. It becomes apparent, that for the numerical model no full debonding is achieved after a rotation of 10 • , since for the nodes close to the center of the implant's contact surface the appropriate sliding distance to start the transitioning of µ has not been reached yet. In fact, under perfect twisting conditions, g s remains zero at the center of the implant. Thus the center of the implant will never debond for (perfect) twisting.

The bottom row of Fig. 5 shows the distribution of the sliding distance over the contact surface for different angles of rotation. Although the body starts sliding before a twisting angle of 1 • , the factor a s prevents a change of the friction coefficient until a sliding distance of 22 µm is reached. This is also shown in Fig. 6b, which illustrates the variation of the friction coefficient as a function of the sliding distance g s for r = R and data set 1. The friction coefficient stays constant as µ = µ ub for g s ≤ a s and then decreases, until it reaches µ = µ b at g s = 66 µm. experiment from [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] analytical solution from [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF]; see Eqs. ( 15),( 16) analytical solution of Eq. ( 18) present numerical solution experiment from [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] analytical solution from [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF]; see Eqs. ( 15),( 16) analytical solution of Eq. ( 18) present numerical solution

• 2 • 3 • 5 • 8 • 10 • µ 0.58 0.5 0.4 1 • 2 • 3 • 5 • 8 • 10 • 0 0.1 0.22 g s [mm]

Work of adhesion

Due to the poor agreement between the experimental results and the analytical model developed in Mathieu et al. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF], an energetic approach was proposed to determine the dissipated frictional energy E fric , the work of adhesion W adh , and the total energy necessary for debonding E deb , 4which are given by

E deb = θ=10 • θ=0 • T (θ) dθ, E fric = θ=10 • θ=θmax T ∞ dθ, W adh = E deb -E fric . (20)
Based on the experimental results, θ = 10 • was chosen to be the angle of rotation where the implant was completely debonded from the bone, indicated by a constant torque T ∞ . The area-specific average work of adhesion wadh is then given by

wadh = W adh πR 2 φ0 , ( 21
)
where φ0 is the average of the initially osseointegrated area, see also Sec. 4.4. A visual analysis of the implant surfaces after debonding yielded φ0 = 0.73 for the first data set (see Fig. 7a) and φ0 = 0.72 for the second.

Tab. 5 gives the results for the different energies with respect to the values of φ0 determined by the surface analysis. In the cases presented here, the analytical and the numerical models use φ 0 (x) = 1 ∀ x. Therefore, the true area-specific work of adhesion for these models w * adh was computed by using φ0 = 1 in Eq. ( 21). Results for wadh where φ0 = 1 is used during the simulation are presented in Sec. 4.4.

For both data sets, the analytical by Mathieu et al. underestimates the total debonding work and the area-specific adhesion work, while the numerical results with the modified Coulomb's law lead to very good agreement with the experimental data. The analytical solution with the modified Coulomb's law yields less accurate results for the second data set.

For simplicity, the results in the remaining part of this work were obtained based only on the estimated parameters for the first data set and µ b = 0.4 (see Tab. 4).

data model E deb [N m] E fric [N m] W adh [N m] wadh [N m -1 ] w * adh [N m -1 ]

Partial osseointegration

Mathieu et al. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] showed that part of the limitations of their model lies in the assumption of a fully bonded surface at the beginning of the experiment. The analysis of the implants' surfaces after debonding indicated that most likely, full osseointegration was not achieved. This resulted in regions where no bone tissue was initially attached to the implant surface, as seen in Fig. 7a, and thus no adhesive or frictional effects can take place. This state of partial osseointegration can be readily investigated with the proposed numerical model. Therefore, an analysis assuming inhomogeneous initial bonding was performed, with various distributions φ 0 (x) in model [START_REF] Corbett | Three-dimensional isogeometrically enriched finite elements for frictional contact and mixed-mode debonding[END_REF].

In order to analyze the debonding behavior of the BII and to determine the influence of the percentage and the distribution of osseointegration, different cases for a fine mesh with 400 contact elements on the implant surface were constructed: First, a bonding pattern based on Fig. 7a, assuming that light gray areas indicate osseointegration, was reconstructed. From a visual inspection of the photograph, the sample in Fig. 7a

Partial osseointegration for different material parameters

When the model parameters are fitted for every pattern (see Tab. 6), the torque-angle curve only weakly depends on the various distributions of osseointegration, as seen in Fig. 9, resulting in a minimum error of 1.95 % for the reconstructed pattern, 1.99 % for the circular pattern and 2.15 % for the random distribution, respectively. While the error is slightly larger compared to the reconstructed pattern, the circular pattern has a closer fit to the experimental data for θ < θ max , which leads to a small improvement of the prediction of wadh . Since the model parameters are unknown, fitting the parameters to the experimental data leads to similar curves for every presented case. It can be seen, that out of the different osseointegration patterns, the random pattern is the closest to the results obtained with full initial bonding. The random pattern still has a balanced degree of osseointegration over the implant radius, while pattern (b) and (c) are only weakly (or not) bonded on the outer part of the implant and thus, a bigger difference is seen in the beginning of the torque curve. Overall, the aforementioned results lead to the conclusion that without knowing the friction coefficients, the torque-per-angle curve does not provide sufficient information on the degree and distribution of osseointegration of the BII. 7).

bonding φ0 d [µm] µ ub a s [µm] b s e mp T [%] W adh [N m] wadh [N m -
As the state of an element does not depend on the states of neighboring elements, the total sliding distance is not affected by inhomogeneous initial bonding. Therefore, a similar propa-gation of the transition zone as in the homogeneous case can be seen in Fig. 8. In contrast to the torque, partial osseointegration has a larger effect on the model parameters, as shown in Tab. 6. In general, partial osseointegration leads to the estimation of higher µ ub , that are still well within the range of reported values in the literature. In addition, the transition time is notably lower and for the distinct patterns (b) and (c), the sliding threshold is lower, as well.

1 • 2 • 3 • 4 • 5 • µ 0.8 0.6 0.4 1 • 2 • 3 • 4 • 5 • µ 0.97 0.69 0.4

Partial osseointegration for equal material parameters

The effect of the different osseointegration patterns becomes much more distinct when studying the results obtained for equal material parameters for bone tissue and the implant, as presented in Fig. 10. Here, the parameters for the homogeneous case of the first data set (µ b = 0.4) were used (see Tabs. 3 and4) with φ0 = 0.55 for all four patterns. Using partial osseointegration only affects the part of the deformation where (tangential) adhesive forces are expected to play a mayor role, represented by the peak in the torque curve. In general, partial osseointegration patterns lead to a lower T max and differences in the softening of the curve. The initial slope of the torque curve and T ∞ are not affected by the different osseointegration patterns. Applying φ 0 = 0.55 to all contact elements in the homogeneous bonding case produces the same result as using a randomly distributed pattern with φ0 = 0.55, where φ 0 = 1 for the osseointegrated elements. The distinct patterns (b) and (c) with weak or no boding on the outer part of the implant produce a lower peak and longer softening periods.

Discussion and Limitations

This work aims to provide a circular friction law to model the debonding of the BII. The mechanical model is incorporated into an analytical model and is implemented into an accurate and efficient contact algorithm with NURBS-enriched finite elements and thus, allows to predict the debonding of osseointegrated implants numerically. The model only depends on four physiological parameters (µ ub , p, a s , and b s ), that can be determined by a parameter study based on experimental results. Moreover, it allows for reasonable estimations of other parameters, such as the bone's Young's modulus and adhesion energy. However, due to the limited experimental data analyzed here, no a-priori determination of the input parameters can be made yet. While the analytical solution of the proposed friction law already leads to good estimates, the results obtained by the finite element simulation are more accurate as they can also account for non-linear material behavior, large deformations, and partial osseointegration.

One difficulty associated with the present study arises from the contact pressure, which is required to establish frictional contact. In Mathieu et al. [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF], it was reported that although the experimental pressure was minimized, it could not be completely excluded nor measured. Since the experimental torque does not go to zero for fully debonded implants, a normal pressure is likely to remain applied to the implants. As friction coefficients and normal force are unknown, no statement can be made about the accuracy of the estimated contact pressure. Furthermore, in the beginning, p accounts for applied normal forces as well as adhesive forces due to chemical and mechanical bonding. The measurement or elimination of applied normal forces would clearly determine the friction coefficient for the broken state and thus also for the unbroken state of a certain sample. Therefore, an improvement of the testing machine used for the experimental measurements is needed.

While the estimated parameters are within a reasonable range compared to the literature (see e.g. Refs. [START_REF] Damm | Friction coefficient and effective interference at the implant-bone interface[END_REF][START_REF] Novitskaya | Recent advances on the measurement and calculation of the elastic moduli of cortical and trabecular bone: a review[END_REF][START_REF] Shirazi-Adl | Experimental determination of friction characteristics at the trabecular bone and porous-coated metal interface in cementless implants[END_REF]), and result in a good qualitative and quantitative representation of the torque-angle curve, there are still visible differences, in particular at the beginning of the peak. Following the argumentation of the reference study, these differences might be explained by the assumption of a full initial bonding between bone and implant, while an initial bonding of 30-70 % is reported in the literature [START_REF] Brånemark | Biomechanical characterization of osseointegrazion during healing: an experimental in vitro study in the rat[END_REF][START_REF] Marin | Biomechanical and histomorphometric analysis of etched and non-etched resorbable blasting media processed implant surfaces: an experimental study in dogs[END_REF]. Accounting for inhomogeneous osseointegration in our model has shown an influence on the torque curves and the model parameters, such as the friction coefficient for the unbroken state. As these values were not or cannot be measured experimentally yet, it is assumed that the parameters obtained by incorporating imperfect osseointegration are more accurate. Furthermore, a partial bonding will most likely lead to a more complex crack front and propagation than what is assumed by the analytical and numerical models.

Other factors that were not taken into account in the present work are the roughness of bone and implant surfaces, as well as wear and debris. Furthermore, only a healing time of 7 weeks was considered, while a comparison of different healing times in terms material parameters and of osseointegration would be interesting. Thus, a study on the influence of the surface roughness and healing time is planned for future work.

Another interesting aspect for future work is the application of the model to actual implant and bone geometries, for example in artificial hip joints and dental implants. The latter have been recently examined by Rittel et al. [START_REF] Rittel | Modeling the effect of osseointegration on dental implant pullout and torque removal tests[END_REF] to study the influence of partial osseointegration on implant stability and cohesive failure. In addition, only torsional debonding was tested in this work while other loading conditions, such as push-in and pull-out of the implant are more commonly analyzed (see e.g. Refs. [START_REF] Berahmani | An experimental study to investigate biomechanical aspects of the initial stability of press-fit implants[END_REF][START_REF] Bishop | The influence of bone damage on press-fit mechanics[END_REF][START_REF] Damm | Friction coefficient and effective interference at the implant-bone interface[END_REF][START_REF] Wennerberg | Nanostructures and hydrophilicity influence osseointegration: a biomechanical study in the rabbit tibia[END_REF]).

Although, especially trabecular bone is highly non-linear and anisotropic, the choice of an isotropic, non-linear elastic Neo-Hookean material model here leads to reasonable results for the planar mode III debonding of titanium and cortical bone, due to the small deformation. It has to be investigated if this still holds for, e.g. pull-out tests and contact with trabecular bone. Especially for the modeling of cohesive failure of bone a fracture model would be needed.

Conclusion

This paper presents an extension of the classical Coulomb's law to a non-constant, statedependent friction coefficient and its application to the debonding of homogeneous and partially osseointegrated implants. An analytical model and numerical solutions with the new friction model are compared to experimental and analytical results of the reference study [START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF]. Overall it is shown, that assuming a smooth transition from an unbroken to a broken state, characterized by a decreasing friction coefficient during the debonding process leads to good agreement with experimental data with both, an analytical and a numerical approach. While the analytical model is simple, it is an efficient way to provide initial estimates for this kind of experiment. The numerical results on the other hand are more accurate and allow for more complex material behavior, stress distribution, and (partial) osseointegration. Both approaches enable the estimation of several parameters of the BII. The proposed friction model is expected to help in understanding the debonding phenomena at the BII. An extension to adhesive friction [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF] and rough surfaces, as well as hip implants, will be considered in future works.

A.1 Contact surface description

Considering a 3D body B, its boundary ∂B, such as its contact surface ∂ c B, can be described by the mapping x = x(ξ), ξ ∈ P,

which maps a point ξ = ξ 1 , ξ 2 lying in the 2D parameter space P to the surface point

x ∈ ∂ c B.
Then, a set of covariant tangent vectors on ∂ c B can be defined as

a α := ∂x ∂ξ α , (α = 1, 2) (23)
with the normal vector

n := a 1 × a 2 ||a 1 × a 2 || . ( 24
)
With Eqs. ( 23) and ( 24), the contact surface can then be characterized by the surface metric components

a αβ = a α • a β , (25)
The dual tangent vectors a α defined by the relation a α • a β = δ α β can be related to the covariant tangent vectors in Eq. ( 23) by a α = a αβ a β . To track changes of ∂ c B during deformation, one chooses a reference configuration ∂ c B 0 , where the tangent vectors A α and the surface metric tensor components A αβ can be defined as in Equations ( 23) and [START_REF] Jasty | In vivo skeletal responses to porous-surfaced implants subjected to small induced motions[END_REF], respectively.

A.2 Contact kinematics and contact traction

To formulate frictional contact between two bodies B 1 and B 2 , we consider interactions between a so called slave point x k ∈ ∂ c B k (k = 1 or 2) and a neighboring master surface ∂ c B (l = 2 or 1). For the classical full-pass contact algorithm (see Laursen and Simo [START_REF] Laursen | A continuum-based finite element formulation for the implicit solution of multibody, large deformation frictional contact problems[END_REF]) k is either set to 1 or 2. Further, we assume point interaction, i.e., x k can only interact with at most one point x ∈ ∂ c B at a given time.

∂ c B x k g e x ( ξ) n a 2 ∂ c B k g en g et g n e x ( ξn ) C a 1
Figure 11: Frictional contact kinematics at current time t n+1 : current slave point x k , current position of the previous interacting point x ( ξn ), current interacting point x ( ξ), the current elastic gap vector g e and its components, previous elastic gap g n e , and the sliding path C. Adopted and modified from Ref. [START_REF] Duong | A concise frictional contact formulation based on surface potentials and isogeometric discretization[END_REF].

In the following, for the sake of conciseness, all variables without superscript n are evaluated at the current time t n+1 , if not stated otherwise. To characterize the interaction between the two bodies we consider a penalty formulation, with the penalty parameter (considered equal in normal and tangential direction). In order to determine the contact traction at the current time step, according to the contact formulation of Ref. [START_REF] Duong | A concise frictional contact formulation based on surface potentials and isogeometric discretization[END_REF], a so-called interacting (elastic) gap vector g e ( ξ) is introduced (see Fig. 11). This gap vector is defined between the current slave point x k and the so-called current interacting point x ( ξ) on the master surface ∂ c B (defined below), i.e.

g e ( ξ) := x k -x ( ξ).

The current gap vector can be further into a tangential and normal component g e ( ξ) = g en + g et , with

g en ( ξ) := (n ⊗ n) g e , g et ( ξ) := (a α ⊗ a α ) g e . (27)
According to the penalty formulation, the total (normal and tangential) frictional contact traction is proportional to the interacting gap vector g e ( ξ), according to

t c = g e , (28)
which follows from using the contact potential W := 1 2 g e •g e . At initial contact, the interacting point x ( ξ) is equal to the closest projection point of x k . During sticking, the current interacting point is equal to the previous interacting point ξn . Therefore, for sticking, the current contact gap vector g e is determined from Eq. ( 26) with ξ = ξn . During sliding, the current interacting point ξ is the solution of the kinematic constraint equation

f g ( ξ) := g et -g max et = 0, (29)
in the current configuration (see Ref. [START_REF] Duong | A concise frictional contact formulation based on surface potentials and isogeometric discretization[END_REF]). g e then follows from Eq. ( 26). The critical value during sliding g max et is defined by the chosen friction law. For example, for Coulomb's friction, it is defined as

g max et = µ g en τ , (30)
where τ denotes the sliding direction and can be computed by projecting the previous interacting gap g n e onto the tangent plane at the current interaction point x ( ξ), i.e.

τ = (a α ⊗ a α ) g n e (a β ⊗ a β ) g n e , (31)
where a α , a α , are evaluated at x ( ξ).

The computation of the friction coefficient µ in Eq. ( 13) requires the knowledge of the accumulated sliding distance g s . Here, we approximate Eq. ( 12) by accumulating the distance from the initial interacting point to the current interacting point, i.e.

g n+1 s ≈ n+1 i=1 x ξi -x ξi-1 . (32)

A.3 Finite element discretization

For the finite element formulation, the contact bodies are discretized into n el elements. The geometry of an element e in the current configuration can be interpolated from the positions of the elemental nodes, or control points x e , as x = N e (ξ) x e , ξ ∈ P,

where N e := [N 1 I, N 2 I, ..., N n el I] denotes the element shape function array. Following the 3D contact enrichment approach of Corbett and Sauer [START_REF] Corbett | NURBS-enriched contact finite elements[END_REF][START_REF] Corbett | Three-dimensional isogeometrically enriched finite elements for frictional contact and mixed-mode debonding[END_REF], the bulk of the bodies is discretized by linear elements with standard Lagrange interpolation following Sauer [START_REF] Sauer | Enriched contact finite elements for stable peeling computations[END_REF], while the contact surfaces are discretized by non-uniform rational B-Splines (NURBS) interpolation (see, e.g., Hughes et al. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]). The NURBS basis functions can be computed element-wise by employing the Bézier extraction operator C e of Borden et al. [START_REF] Borden | Isogeometric finite element data structures based on Bezier extraction of NURBS[END_REF] and the Bernstein polynomials B. Then, the shape function of control point A can be computed by

N A (ξ 1 , ξ 2 ) = w A N e A (ξ 1 , ξ 2 ) n A=1 w A N e A (ξ 1 , ξ 2 ) , (34)
with the corresponding weight w A , and the B-spline basis functions

Ne = { N e A } ne A=1 , Ne (ξ 1 , ξ 2 ) = C e ξ 1 B(ξ 1 ) ⊗ C e ξ 2 B(ξ 2 ). (35)
The finite element forces acting on the slave surface element e and the master surface element ê are then given by

f e c := Γ e
N T e t c da (slave),

f ê c := - Γ ê N T ê ( ξ) t c da (master), (36)
where Γ denotes the surface element domains on the corresponding surfaces. As NURBS elements are at least C 1 continuous over the element boundaries, this gives us an accurate representation of the geometry, as well as an efficient and stable framework for contact computations.

B Derivation of the new analytical model with modified Coulomb's friction

Given the modified Coulomb's law in Sec. 2.2.3, a new analytical model can be derived. In general, the torque T , as a function of the rotation angle θ, can be computed analytically, as

T (θ) = 2π
Assuming bone and implant to be linear elastic bodies, we know that for sticking, the tangential traction will be proportional to the applied rotation angle θ and the current radius r, i.e.

σ stick θz (θ, r) = λθr, λ = const, [START_REF] Rho | Young's modulus of trabecular and cortical bone material: ultrasonic and microtensile measurements[END_REF] where λ is a constant stress per length. Thus, the torque for the sticking region becomes We also know that for linear elasticity dT (θ) dθ = C, θ ≤ θ lin ,

T (θ) = 2π
where C is the effective shear stiffness of the system, that is approximately proportional to the shear modulus of bone. If we equate Eqs. ( 41) and ( 42), we obtain

λ = 2C πR 4 , (43)
and inserting λ into Eq. ( 39) yields σ stick θz (θ, r) = 2C πR 4 θr.

Applying the definition in Eq. [START_REF] Rancourt | Friction properties of the interface between porous-surfaced metals and tibial cancellous bone[END_REF] to Eq. ( 44) yields

σ stick θz (θ lin , R) = 2C πR 4 θ lin R = t max t , (45)
i.e. 2C πR 3 =

t max t θ lin . (46)
Together with Eq. ( 17) we finally obtain the tangential traction component for sticking

σ stick θz (θ, r) = θr t max t θ lin R = r c(θ) t max t . (47)
In the sliding region r ≥ c, the tangential traction follows the modified Coulomb's law

σ slide θz = µp, (48)
with µ = µ(g s ) from Eq. [START_REF] Chong | Analysis of bone-prosthesis interface micromotion for cementless tibial prosthesis fixation and the influence of loading conditions[END_REF]. Here, we assume p to be homogeneously distributed along the contact surface. As long as a part of the interface remains sticking g s = 0, otherwise g s = r (θ -θ max + θ lin ) .

The tangential contact traction as a function of the rotation angle can then be summarized by combining Eq. ( 47) and Eq. ( 48) to Eq. ( 18).

C Mesh and load step sensitivity

To analyze the convergence behavior, three different finite element meshes were constructed, denoted coarse, medium, and fine with 3,390, 12,354, and 47,130 degrees of freedom, respectively. In addition, different load step sizes [0.1 T exp (θ) -T (θ)

T max , (50)
where here T max is the maximum torque obtained by the numerical solution.

Fig. 13 shows the convergence behavior of the different meshes. It can be seen that e mp T reaches its limit for all meshes after 1,000 load steps to 0.0175, 0.0171, and 0.017, respectively. This is also the case for the mean percentage error shown in Tab. 7, with its lowest value of 2.176 % for the fine mesh and the highest value of 2.241 % for the coarse mesh. In addition, the error is increasing with the number of load steps for the coarse mesh. This stems from the coarse resolution of the peak for larger load steps and thus leading to the torque values to be closer to the experimental data. It should be noted, that the mesh size has a small effect on the outcome of the parameter study and thus, both e mp

T and e rel T can be further minimized by performing a separate parameter estimation for each mesh.

( a )

 a Variation of the sliding threshold a s for b s = 1. Variation of the transition factor b s for a s = 0.22.

Figure 1 :

 1 Figure 1: Modified Coulomb's law: Tangential contact traction t t /t max t as a function of the sliding distance g s for varying sliding threshold a s and transition factor b s for φ 0 = 1.

  Schematic representation of the setup and the crack propagation in the uniformly bonded case. (b) NURBS-enriched finite element mesh showing an enlargement of the contact surface.

Figure 2 :

 2 Figure 2: Illustrations for the analytical and numerical setup.

Figure 3 :

 3 Figure 3: New analytical model: Normalized tangential traction σ θz /t max t as a function of the radius r/R of the contact surface (see Eq. (18)) at different rotation angles.

( b )

 b Second data set.

Figure 4 :

 4 Figure 4: Variation of the Torque T as a function of the imposed rotation angle θ.

1

 1

Figure 5 :θ

 5 Figure 5: Evolution of the debonding of the BII: Value of the friction coefficient µ (top) and the sliding distance g s (bottom) on the contact area of the implant for angles of rotation.

Figure 6 :

 6 Figure 6: Evolution of the debonding of the BII: Behavior of the friction coefficient and the transition zone for the first data set.

  is considered to have average osseointegration φ0 = 0.55.Second, to compare the influence of osseointegration patterns, two other patterns with φ0 = 0.55 were constructed. A simple circular pattern, where only the center part of the interface is integrated and last, a random distribution. For patterns (b)-(d), φ 0 = 1 within the light gray surface elements. The corresponding osseointegration patterns are shown in Fig. 7. (a) Original pattern [32]. (b) Reconstructed pattern. (c) Circular pattern. (d) Random pattern.

Figure 7 :

 7 Figure 7: Original and constructed osseointegration patterns with average osseointegration φ0 = 0.55. Light areas represent full osseointegration (φ 0 = 1), while dark areas represent no osseointegration (φ 0 = 0). For the corresponding elements in the reconstructed and artificial patterns φ e 0 is set to 1.

Figure 8 :Figure 9 :

 89 Figure 8: Variation of the friction coefficient µ at the BII for different angles of rotation for different patterns of partial osseointegration.

Figure 10 :

 10 Figure 10: Partial osseointegration: Variation of the Torque T as a function of the imposed rotation angle θ for different patterns of the initial bonding state of the interface with φ0 = 0.55 and the same material parameters for all patterns. T ∞ = 0.0368 N m according to data set 1 (see Tab. 2).

r r 2 •

 2 σ θz (r, θ, g s ) dr.[START_REF] Pilliar | Observations on the effect of movement on bone ingrowth into porous-surfaced implants[END_REF] Due to symmetry, the tangential traction component σ θz is distributed radially symmetric along the BII, while the radial traction component σ zr is zero. With the definition of the critical radius c from Eq. (17) we have c(θ lin ) = R and σ θz (θ lin , c(θ lin )) := t max t .

R 0 r 2 •Figure 12 :

 0212 Figure 12: New analytical model: Critical radius for stick/ slip transition c as a function of the imposed rotation angle θ from Eq. (17).

  load

Figure 13 :

 13 Figure 13: Mesh sensitivity: Mean relative error e rel T according to Eq. (50) for different configurations of data set 1 (µ b = 0.4).

Table 1 :

 1 All computations in Secs. 4.1-4.3 use homogeneous initial bonding (φ 0 (x) = 1, ∀x ∈ ∂ c B), while Sec. 4.4 presents cases with inhomogeneous initial bonding. Parameters

	body	n e type of shape fcts.	n gp
	implant bulk	18	linear Lagrange	2 × 2 × 2
	bone bulk	450	linear Lagrange	2 × 2 × 2
	lower implant surface	9	quadratic NURBS	5 × 5
	upper bone surface	225 quadratic NURBS	5 × 5

• . of the finite element mesh: Number of elements n e , type of shape functions and number of Gauss-points per element n gp for the two bodies and their contact surfaces.

Table 2 :

 2 Data sets used for the parameter estimation. Data set 1 and 2 correspond to the data shown in Figure

		.0538	0.0368	0.13 • 1.13 •
	2	0.0595	0.0444	0.11 • 2.02 •

Table 4 :

 4 Results for the parameter estimation of the two data sets. Shown are results from the analytical model presented in Ref.[START_REF] Mathieu | Mode III cleavage of a coin-shaped titanium implant in bone: Effect of friction and crack propagation[END_REF] (see Eqs. (

	%]

Table 5 :

 5 Total debonding energy E deb , frictional energy E fric , work of adhesion W adh , and corresponding area-specific works of adhesion wadh and w * adh for the different models and data sets.

	set
		exp.	0.0070	0.0056	0.0014	98	71
	1	ana. Ref. [32] ana. Eq. (18)	0.0066 0.0071	0.0056 0.0057	0.0010 0.0014	70 98	51 72
		sim.	0.0071	0.0057	0.0014	98	72
		exp.	0.0088	0.0060	0.0028	198	143
	2	ana. Ref. [32] ana. Eq. (18)	0.0080 0.0087	0.0064 0.0067	0.0016 0.0020	120 141	81 102
		sim.	0.0089	0.0062	0.0027	191	138

Table 6 :

 6 Change in model parameters and results for implants with partial initial bonding compared to homogeneous bonding (see Fig.

	1 ]

  • , 0.05 • , 0.02 • , 0.01 • , 0.005 • , 0.004 • ] were investigated, corresponding to a number of load steps of [100, 200, 500, 1,000, 2,000, 2,500], respectively. For the parameters, data set 1 with µ b = 0.4 was chosen. To compare the different setups, we define the mean relative torque error

	e rel T = mean θ∈[0,10 • ]

Table 7 :

 7 Mesh sensitivity: Mean percentage error e mpT according to Eq. (19) for different configurations of data set 1 (µ b = 0.4).

				0.0178	coarse
	step 0.10 •	coarse medium 0.02233 0.02195 0.02186 fine		0.0176	medium fine
	0.05 • 0.02 •	0.02242 0.02197 0.02185 0.02239 0.02191 0.02177	e rel T	0.0174
	0.01 • 0.005 •	0.02240 0.02190 0.02176 0.02241 0.02190 0.02176		0.0172
	0.004 •	0.02241 0.02190 0.02176		0.0170
				0	500 1,000 1,500 2,000 2,500
					load steps

In principle Eq. (13) can be brought into the form of Eq. (10) if a dot is applied to Eq. (13) and then gs is eliminated by the inverse function of[START_REF] Corbett | Three-dimensional isogeometrically enriched finite elements for frictional contact and mixed-mode debonding[END_REF].

In principle, φ0 follows from an evolution law that describes the healing/osseointegration process (φ0 increasing from 0 to 1). However, this is not considered here. Here, only the (further) evolution of φ during debonding (φ decreasing from 1 to 0) is studied.

(R2018b, The MathWorks, Natick, MA, USA)

The strain energy inside the bodies, which generally should also be accounted for, is negligible in this case.

Acknowledgments

This work has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No 682001, project ERC Consolidator Grant 2015 BoneImplant).

Conflict of interest

The authors declare that they have no conflict of interest.

Appendix

A Contact formulation

In this section the basic contact mechanics and the numerical discretization are briefly summarized. For the formulation of sliding contact, one can use elastoplasticity theory in order to incorporate friction response in the tangential direction (see e.g. Refs. [START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF][START_REF] Wriggers | Computational contact mechanics[END_REF][START_REF] Sauer | An unbiased computational contact formulation for 3d friction[END_REF]). In the context of this work, we use another approach that is based on the surface potential-based contact formulations of Duong and Sauer [START_REF] Duong | A concise frictional contact formulation based on surface potentials and isogeometric discretization[END_REF] and Sauer and de Lorenzis [START_REF] Sauer | A computational contact formulation based on surface potentials[END_REF]. The reader is referred to Corbett and Sauer [START_REF] Corbett | NURBS-enriched contact finite elements[END_REF][START_REF] Corbett | Three-dimensional isogeometrically enriched finite elements for frictional contact and mixed-mode debonding[END_REF] for the NURBS-enriched contact finite element discretization.


X Iltis

H Palancher

F Vanni

J Allenou

B Stepnik

A Leenaers

S Van Den Berghe

I Glagolenko

D D Keiser

Characterization of fresh EMPIrE and SEMPER FIDELIS plates with PVD-coated U(Mo) particles

The HERACLES group and the US-DOE work jointly to develop dispersed U(Mo)/Al as LEU fuel for conversion of high performance research reactors. Within this frame, two irradiation programs are in progress. In the first, EMPIrE, mini-plates are tested in the ATR reactor (USA) and in the second, SEMPER FIDELIS, full-size plates are irradiated in BR2 (Belgium). In both experiments, U(Mo)/Al plates with optimized microstructure are tested under aggressive conditions. This paper focuses on analyses made at CEA Cadarache on seven fresh plates made of atomized particles, with or without Mo homogenization, and with ZrN coating. Five EMPIrE mini-plates and two SEMPER FIDELIS full-size plates were examined by Optical Microscopy (OM), Scanning Electron Microscopy (SEM) and Energy Dispersive X-ray Spectroscopy (EDS). A particular attention is paid to the integrity of the ZrN coating (thickness, cracks…) and to the U(Mo) particles microstructure.

Introduction

The use of low enriched uranium (LEU) nuclear fuels for research reactors (material testing reactors and neutron sources) is worldwide encouraged. A fuel core made of a dispersion of U(Mo) alloy (with 7 to 10 wt% Mo) particles within an Al matrix, colaminated between two aluminium alloy plates, constitutes one of the candidates for this type of fuel. Under aggressive irradiation conditions, the performances of U(Mo)/Al dispersed fuels are limited by an interaction process between U(Mo) particles and the Al matrix, which induces a significant swelling of the fuel plates [START_REF] Van Den | Review of 15 years of high-density low-enriched UMo dispersion fuel development for research reactors in Europe[END_REF]. A 1 µm thick ZrN layer can strongly decrease this interaction but, at high burn-up, a swelling acceleration is still observed and is attributed to a recrystallization mechanism of the fuel which induces both U(Mo) grain refinement and precipitation of large (micrometre sized) fission gas bubbles [START_REF] Leenaers | Fuel swelling and interaction layer formation in the SELENIUM Si and ZrN coated U(Mo) dispersion fuel plates irradiated at high power in BR2[END_REF]. A way to delay this recrystallization and, thus, the associated swelling under irradiation, could consist in modifying the U(Mo) particle microstructure (which is an "as-solidified" microstructure, in powders produced by an atomization process), using homogenization thermal treatments, as pointed out in the frame of the KOMO-5 experiment [START_REF] Kim | In-pile test results of silicide or U-nitride coated U-7Mo particle dispersion fuel in Al[END_REF].

Verifying the benefits of a homogenization heat treatment, followed by a ZrN coating of the particles, on the swelling behaviour of U(Mo)/Al dispersed fuel, is one of the main goals of the EMPIrE and SEMPER FIDELIS irradiation experiments [START_REF] Glagolenko | EMPIrE -European mini-plate irradiation experiment in the Advanced Test Reactor (ATR)[END_REF] [START_REF] Stepnik | Manufacturing of the Semper Fidelis UMo irradiation experiment[END_REF]. Both in-pile tests are currently ongoing in the ATR (Idaho, USA) and BR-2 (Mol, Belgium) respectively. This work is aimed to characterize fresh fuel plates from both experiments. We have chosen to focus on PVD-coated plates, made with atomized U(Mo) particles produced by KAERI. If only some powder batches were homogenized at 1000°C for one hour, all of them were coated by PVD with ZrN in conditions similar to those used for SELENIUM experiment [START_REF] Leenaers | Surface Engineering of low enriched uranium-molybdenum[END_REF]. Both powder treatments (high temperature annealing and powder coating) were performed at SCK-CEN. Finally, plates were manufactured by CERCA using the conventional pictureframe technique.

Experimental details

Materials

Table 1 recaps the main characteristics of the seven plates examined here. Four of them are made with heat treated particles (three mini-plates and one full-size plate). The two types of plates considered here (EMPIrE and SEMPER-FIDELIS) are both cladded with AG3NE aluminium sheets. Mini and full-size plates were manufactured using close conditions. A 9×28 mm 2 piece was taken from each plate at CERCA and sent to CEA Cadarache. There, two samples per piece were cut and then mechanically polished: -a 9×9 mm 2 square sample, polished parallel to the fuel-cladding interface, down to the middle part of the fuel zone, for XRD (X-Ray Diffraction) and SEM (Scanning Electron Microscopy) coupled with EDS (Energy Dispersive Spectroscopy) examinations, -a 9 mm long cross-section, for OM (Optical Microscopy) and complementary SEM+EDS examinations, especially close to the fuel-cladding interface.

Characterization methods

OM examinations were performed with an Olympus DSX500 opto-digital microscope. SEM examinations were carried out mainly in backscattered electrons mode (BSE), with a FEI Nova Nano SEM 450, equipped with an Oxford Instruments EDS system. The mean Mo content in U(Mo) particles was determined by EDS, on 15 particles per plate [START_REF] Iltis | Microstructure of as-atomized and annealed U-Mo7 particles: A SEM/EBSD study of grain growth[END_REF].

The STREAM Olympus software was used to analyse BSE images in order to measure:

-the U(Mo) particles sphericity (defined as the squared ratio between the width and the length of a particle, excluding ZrN coating), on about 100 particles per plate, with a minimal diameter of 20 µm, -the thickness of ZrN coating, on about 30 particles per plate.

In order to minimize the potential measurement error due to bias cutting effects, only particles with diameters larger than 70 µm were considered for ZrN thickness measurements.

Those with significantly damaged coatings (partly delaminated and/or fragmented) were also disregarded. On each particle, ZrN layer boundaries were detected thanks to the grey level variations at the Al/ZrN and ZrN/U(Mo) interfaces (Figure 1a). After that, the coating thickness was measured at 100 different locations regularly distributed around the particle (Figure 1b). Statistical results presented in section 3.3 are based on mean values obtained from each particle.

Results

Macrographs

Figure 2 gathers four representative macrographs, obtained by optical microscopy on crosssections from two EMPIrE and two SEMPER FIDELIS plates, each couple of plates comprising one plate made with heat treated particles. At this scale, the only noticeable difference is that EMPIrE mini-plates seem to present a more heterogeneous distribution of U(Mo) particles within the fuel core (in accordance with less demanding manufacturing specifications concerning its homogeneity).

U(Mo) particles characteristics

Figure 3 shows SEM micrographs taken on each plate, in imaging conditions allowing grains visualization within U(Mo) particles. As expected, the three plates made with as-atomized particles exhibit typical columnar and/or cellular solidification microstructures (Figures 3a, 3b and 3g). The four other ones are characterized by large grains which reach a size of the order of several tens of microns, as expected after a 1000°C-1h annealing [7][8]. In some cases, single crystalline particles are encountered (see for example Figure 3d: particle on the left hand side). Grains are more or less difficult to distinguish, since they almost systematically contain lamellas with darker contrasts. Such features were already observed in certain particles from heat treated KOMO-5 powder batches and were interpreted as Widmanstätten microstructures which could appear depending on the cooling rate conditions [START_REF] Iltis | Microstructure of as-atomized and annealed U-Mo7 particles: A SEM/EBSD study of grain growth[END_REF]. They were also present in certain particles from heat treated EMPIrE powder batches (even if this point was not illustrated in [START_REF] Iltis | Microstructure of as-atomized and annealed U-Mo7 particles: A SEM/EBSD study of grain growth[END_REF]), but to a lower extend than in the EMPIrE plates. Batch to batch variabilities linked to cooling conditions after homogenization and/or an incidence of the thermal history during plates manufacturing could explain why these lamellas are numerous in plates made with heat treated U(Mo) particles. XRD analyzes are in progress and should give clues about the origin of these features. Table 2 recaps the results of sphericity measurements on U(Mo) particles, by image analysis, and of EDS semi-quantitative measurements of Mo content in these particles, for the seven studied plates. The grey colored lines, in this table as in Table 1, correspond to plates made with non-homogenized particles.

Circularity measurements were already performed on several PVD-coated EMPIrE powder batches [START_REF] Vanni | Characterizations of the LEU powders used in EMPIrE irradiation experiment[END_REF]. They tended to evidence a deformation of originally round atomized particles, which was probably linked to their swirling in the coating reactor, and didn't seem to vary significantly with their metallurgical state. However, in EMPIrE plates, except for EMP-819 (probably because of sampling effects), the sphericity of particles seems to be related to this state: it is close to 0.9, for as-atomized particles, while this value is only about 0.75 in heat treated ones. These results tend to show that homogenization leads to a softening of U(Mo) particles, linked to grain growth, which favors their deformation during rolling. This difference of shape between annealed and asatomized particles is not observed in the case of SEMPER FIDELIS plates, which are both characterized by a 0.71 mean sphericity value. This value suggests that the full-size plates rolling conditions led to a slightly greater deformation of U(Mo) particles, whatever their mechanical properties, compared to mini-plates.

Mean particle Mo content is almost the same in all plates, i.e. about 6.3 wt% (standard deviation: 0.2 wt%). This amount is fully consistent with those measured, with the same EDS equipment and the same analysis procedure, in several EMPIrE powder batches [START_REF] Iltis | Microstructure of as-atomized and annealed U-Mo7 particles: A SEM/EBSD study of grain growth[END_REF]. It is underestimated, compared to the content given by KAERI (7 wt%), very probably because of a systematic error due to virtual standards implemented in the EDS software.

ZrN coating characteristics

General features Damaged coatings:

In all studied plates, the ZrN coating around U(Mo) particles appears to be damaged, whereas few cracks could be evidenced on corresponding powders. It is difficult to determine precisely a representative amount of defective coatings, since it can largely vary from place to place, in a same sample, depending in particular on the local density of particles. In any case, there are never more than 50% of the particles exhibiting a perfectly continuous and adherent coating. This ratio can even fall below 30%, in certain cases, especially close to the fuel-cladding interface.

Figures 4 and5 illustrate common degradation modes for ZrN coatings. The first one consists in an obvious delamination of the coating, around certain particles, in a preferential direction which is parallel to the rolling one. Figure 4a corresponds to an area where this phenomenon is particularly clear (the rolling direction is indicated by the white arrow). It is accompanied by the formation of a void between the ZrN layer (which remains unexpectedly "stuck" on the matrix) and the particle, this void reaching frequently several microns in thickness. This type of damage was mainly observed in EMPIrE mini-plates, manufactured with non-annealed U(Mo) particles (i.e. EMP-711 and EMP-717 plates). This suggests that rolling conditions of mini-plates, combined with relatively hard (i.e. poorly deformable) asatomized particles could contribute to such coating delamination. The second degradation mode is directly linked to the ZrN thickness. As it will be quantified in the next section, this thickness can vary by a factor two, from a particle to another one, and its mean value can also vary by the same order of magnitude, in the set of studied plates. Figures 4b and4c is very probably due to deformation incompatibilities between the coating and its substrate, during plate manufacturing and in particular during temperature transients as proposed in literature [START_REF] Keiser | Microstructural characterization of a thin film ZrN diffusion barrier in an As-fabricated U-7Mo/Al matrix dispersion fuel plate[END_REF].

When particles are significantly deformed and/or in direct contact, as illustrated by Figure 5a, the ZrN coating frequently peels off, or breaks down into small pieces, which can be dispersed in the aluminium matrix. U(Mo) particle distribution heterogeneity, in the fuel core, thus strongly influences the ZrN layers physical state.

Oxidized particles:

Another factor has a large influence on the mechanical stability of ZrN coatings: the presence of a brittle, highly circumferentially cracked oxidized layer, on the surface of some U(Mo) particles, below the coating, considerably reduces the adhesion of the ZrN layer. Such oxidized particles are preferentially observed close to the interface with the cladding, in all plates, as illustrated by Figure 5b and Figure 6, but some can also be encountered in the middle of the fuel core, in certain plates, as it will be discussed later.

EDS profiles presented in Figure 6d show that the two analyzed ZrN coatings are partly oxidized: their oxidation probably occurred during the storage of the coated powder before plate manufacturing, as proposed in [START_REF] Keiser | Microstructural characterization of a thin film ZrN diffusion barrier in an As-fabricated U-7Mo/Al matrix dispersion fuel plate[END_REF]. In Figure 6d, two ZrN layers which cover two different particles can be seen: in the following, we will focus on the U(Mo) particle located on the left hand side. Two sub-layers are found between the U(Mo) core and its coating: the outer sub-layer contains both nitrogen and oxygen, whereas the inner, which is thicker, mainly corresponds to an oxide (with a composition close to UO 2 ).

Severely oxidized particles, such as that presented in Figure 6, were already encountered in annealed EMPIrE powders. This allows concluding that their oxidation took place during their heat treatment. That is the reason why such particles can be found in the middle part of the fuel core, in plates made with homogenized U(Mo) particles. Their density varies from plate to plate and is particularly high in the EMPIrE 803 one, in accordance with observations previously carried out on the powder batch used to manufacture it. As already stated, oxidized particles are also present in all plates, close to the fuel-cladding interface, one side of the plate being sometimes more affected than the other side. This localization indicates that they probably oxidized during the plate hot rolling. The fraction of oxidized particles in the plates evolves basically as follows (from plates with numerous oxidized particles to ones with very few ones): EMPIrE mini-plates made with homogenized U(Mo) particles (EMP-803, EMP-819, EMP-828) > EMPIrE mini-plates made with as-atomized particles (EMP-711, EMP-717) > SEMPER FIDELIS plate made with homogenized particles (SF-202) > SEMPER FIDELIS plate made with as-atomized ones (SF-402). In the most oxidized plates, small UO 2 fragments, sometimes still partially coated with ZrN, can be locally dispersed in the matrix, especially close to the cladding. The white arrow, in Figure 7c, points such a fragment. Finally, it is interesting to remind that, before plate fabrication, PVD-coated U(Mo) batches also comprise large ZrN flakes and much smaller spherical ZrN particles (Figure 7a). The fraction of flakes and spherules can vary, from a powder batch to the other and, thus, from a plate to the other too. In the seven studied plates, it remains relatively low and comparable. ZrN spherules often agglomerate onto U(Mo) particles (Figure 7b), whereas ZrN flakes are found broken into pieces in the Al matrix (Figure 7c). ZrN fragments and spherules tend to form layers, roughly parallel to the rolling direction, especially close to the cladding interface (Figure 7d).

Thickness measurements

Table 3 summarizes the ZrN coating thicknesses measured on the seven studied plates, with the method previously illustrated in Figure 1. In all cases, mean values range from 1 to 2 µm, the thinnest layer being found in the SF-202 plate and the thickest in the EMP-819 mini-plate. Standard deviations are of the same order of magnitude for each sample.

Conclusion

In this study, five mini-plates, from EMPIrE test matrix, and two full-size plates, from SEMPER FIDELIS one, were examined by Optical Microscopy (OM), Scanning Electron Microscopy (SEM) and Energy Dispersive X-ray Spectroscopy (EDS). All of them were made with atomized U(Mo) particles (from KAERI), homogenized or not at 1000°C, and PVDcoated with a ZrN layer.

OM macrographs showed that the fuel core thickness and the distribution of U(Mo) particles within it were more irregular in EMPIrE mini-plates. U(Mo) particles microstructure was similar in both types of plates and corresponded either to a fine solidification microstructure, for plates made with as-atomized particles, or to relatively large equiaxed grains containing lamellas, for plates made with annealed powder.

In EMPIrE plates, the sphericity of particles seemed to be related to their metallurgical state (as-atomized or annealed). This difference was not observed in the case of SEMPER FIDELIS plates, which were both characterized by a lower particle mean sphericity value, suggesting that the full-size plates rolling conditions led to a slightly greater deformation of U(Mo) particles, whatever their mechanical properties, compared to mini-plates. Mean particle Mo content, measured by EDS, was almost the same in all plates.

In all studied plates, coated particles often exhibited damaged ZrN layers. The following main types of defects were identified: (i) a delamination of the coating, around certain particles, in a preferential direction parallel to the rolling one, (ii) the development of radial cracks in ZrN especially in the thicker layers, (iii) multiple cracks in the coating on significantly deformed particles or on those in direct contact, (iv) a considerable reduction of the adhesion of the ZrN layer when an underlying oxidized layer was present at the surface of U(Mo) particles, as often observed close to the fuel-cladding interface.

All of these damages were difficult to quantify. Delaminations occurred preferentially in plates made with as-atomized particles, and seemed to be more marked in EMPIrE ones. Radial cracks were closely related to ZrN thickness, which varied between nearly 1 and 2 µm in both sets of plates. Damages linked to contacts were high in all plates. Finally, those linked to the powder oxidation were more numerous in plates made with homogenized particles (as heat treatment at 1000°C induced the oxidation of some particles), and especially in EMPIrE ones.

Figure 1 :

 1 Figure 1: Principle of ZrN coating thickness measurement, on BSE images: (a) detection of ZrN layer boundaries, (b) measurement of ZrN layer thickness at 100 locations around the particle.

Figure 2 :

 2 Figure 2: OM macrographs of (a) EMP-717, (b) EMP-828, (c) SF-202, (d) SF-402 plates.

  717 mini-plate, made with as-atomized particles (b) EMPIrE 828 mini-plate, made with annealed particles (c) SEMPER FIDELIS 202 full-size plate, made with annealed particles (d) SEMPER FIDELIS 402 full-size plate, made with as-atomized particles

Figure 3 :

 3 Figure 3: U(Mo) particles microstructure, observed by SEM (BSE mode) in (a) EMP-711, (b) EMP-717, (c) EMP-803, (d) EMP-819, (e) EMP-828, (f) SF-202, (g) SF-402 plates.

Figure 4 :

 4 Figure 4: SEM micrographs taken on (a) EMP-717 (cross-section), (b) SF-202, (c) EMP-819 plates.

  compare two particles coated by respectively a one and a two microns thick layer. The second one undoubtedly presents more radial cracks. This cracking

Figure 5 :

 5 Figure 5: SEM micrographs taken on cross-sections of (a) EMP-828, (b) EMP-819 plates (white arrows: rolling direction -dashed line: approximate location of fuel-cladding interface).

Figure 6 :

 6 Figure 6: SEM/EDS analysis of an oxidized U(Mo) particle in SF-202 plate (a) general view (white arrow: rolling directiondashed line: approximate position of fuel-cladding interface), (b) and (c) location of the analysis line, (d) EDS profiles along the yellow line drawn in (c).

Figure 7 :

 7 Figure 7: SEM micrographs taken on: (a) a PVD-coated powder from one EMPIrE batch, (b) SF-202 plate, (c) EMP-819 plate, (d) EMP-828 plate (white arrow: rolling directiondashed line: approximate position of fuel-cladding interface).

Table 1 :

 1 Main characteristics of the plates.

	Irradiation program	Plate designation	Mo homogenization
	EMPIrE	EMP-711	No
	(mini-plates)	EMP-717	No
		EMP-803	Yes
		EMP-819	Yes
		EMP-828	Yes
	SEMPER FIDELIS	SF-202	Yes
	(full-size plates)	(full designation: FIDJ0202)
		SF-402	No
		(full designation: FIDJ0402)

Table 2 :

 2 Sphericity and Mo content of U(Mo) particles.

	Plate designation	Mean sphericity	Mean Mo content (wt%)
	EMP-711	0.90	6.2
	EMP-717	0.85	6.2
	EMP-803	0.76	6.5
	EMP-819	0.85	6.3
	EMP-828	0.73	6.3
	SF-202	0.71	6.3
	SF-402	0.71	6.3

Table 3 :

 3 ZrN coating thickness measurements.

	ZrN coating thickness (µm)

Acknowledgements:

A. Sanchez and D. Drouan, from CEA Cadarache, are warmly thanked for their careful preparation of the samples and for their examinations by optical microscopy. D. Drouan contributed also to SEM and EDS characterizations, with N. Tarisien. They both produced a large number of results, in a record time! Finally, many thanks to N. Tarisien for his very good and huge work in image analysis. This project has received funding from the EURATOM research and training program 2014-2018 under grant agreement n°661935.