Micromechanical modeling of the contact stiffness of an osseointegrated bone–implant interface

Background The surgical success of cementless implants is determined by the evolution of the biomechanical properties of the bone–implant interface (BII). One difficulty to model the biomechanical behavior of the BII comes from the implant surface roughness and from the partial contact between bone tissue and the implant. The determination of the constitutive law of the BII would be of interest in the context of implant finite element (FE) modeling to take into account the imperfect characteristics of the BII. The aim of the present study is to determine an effective contact stiffness \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( {K_{c}^{\text{FEM}} } \right)$$\end{document}KcFEM of an osseointegrated BII accounting for its micromechanical features such as surface roughness, bone–implant contact ratio (BIC) and periprosthetic bone properties. To do so, a 2D FE model of the BII under normal contact conditions was developed and was used to determine the behavior of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_{c}^{\text{FEM}}$$\end{document}KcFEM. Results The model is validated by comparison with three analytical schemes based on micromechanical homogenization including two Lekesiz’s models (considering interacting and non-interacting micro-cracks) and a Kachanov’s model. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_{c}^{\text{FEM}}$$\end{document}KcFEM is found to be comprised between 1013 and 1015 N/m3 according to the properties of the BII. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_{c}^{\text{FEM}}$$\end{document}KcFEM is shown to increase nonlinearly as a function of the BIC and to decrease as a function of the roughness amplitude for high BIC values (above around 20%). Moreover, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_{c}^{\text{FEM}}$$\end{document}KcFEM decreases as a function of the roughness wavelength and increases linearly as a function of the Young’s modulus of periprosthetic bone tissue. Conclusions These results open new paths in implant biomechanical modeling since this model may be used in future macroscopic finite element models modeling the bone–implant system to replace perfectly rigid BII conditions.


Background
Endosseous cementless implants have been used clinically for more than 40 years and have allowed considerable progresses in dental, maxillofacial and orthopedic surgery, to replace missing organs or to restore joints functionality. Despite their routine clinical use, implant failures still occur and remain difficult to predict. One of the main determinants of the surgical success lies in the determination of the implant stability [1]. The primary stability occurs during implant surgery. It is a phenomenon of biomechanical nature related to bone quality at the implant site and to the implant properties. A good implant primary stability is a necessary condition to obtain implant osseointegration. Secondary stability is obtained after a certain healing period and corresponds to the initial stability reinforced by newly formed bone production and maturation at the bone-implant interface [2,3]. Osseointegration phenomena [4] are stimulated by the application of mechanical stimuli to the bone-implant interface (BII). Moreover, the implant surface roughness, which is obtained using different processes such as, for example, sand blasting [5], plasma spraying [6], or laser blasting [7], is known to strongly influence the quality of osseointegration phenomena [8]. When the primary stability is not sufficient, micro-movements may appear, preventing good healing conditions and leading to the formation of fibrous tissue and eventually to surgical failure [9,10]. Relative micromotions between the implant and bone tissue should not exceed around 150 µm, because it leads to fibrous tissue formation rather than bone ongrowth [11,12]. Fibrous tissue may develop instead of an osseointegrated interface when there are excessive interfacial micromotions early after surgery [13,14]. However, micromotions at a relatively low level may be responsible for biomechanical stimulation of bone remodeling.
Consequently, it is important to understand the biomechanical behavior of the BII at the microscopic scale, which depends on the bone geometrical and material properties as well as on the implant surface roughness. Experimental approaches remain of limited interest to retrieve the main determinant of the micromechanics of the BII because many parameters such as the bone-implant contact (BIC) ratio as well as bone material properties are difficult to assess and are likely to vary in parallel. Despite the development of acoustical methods [15,16] to retrieve information on the BII properties, it remains difficult to employ noninvasive techniques. Classical imaging techniques such as magnetic resonance imaging or X-ray microcomputed tomography cannot be used in vivo due to diffraction effects related to the presence of the titanium [17,18].
Finite element (FE) models have been widely used to model the implant biomechanical behavior at the organ scale in the context of primary and secondary stability, but often the BII is modeled as a perfect interface, i.e., continuity of stresses and displacements at the interface [19]. In particular in the context of dental implantology, microfinite element analyses were applied to images obtained using X-ray microcomputed tomography [20,21], which allowed to assess the strain and stress field around the implant. However, the BII was often considered as fully bonded. Some groups modeled the BII during osseointegration as an interphase, considering a thin layer described by the Drucker-Prager plasticity model [22]. The use of springs to model the BII was introduced by Egan and Marsden [23], who described load transfer at the BII by a network of linear springs whose stiffness can vary in time.
However, it remains difficult to relate the macroscopic behavior of the BII in terms of contact stiffness [23] or continuous mechanical properties [22] to its micromechanical properties such as the BIC and bone properties. Establishing such a relationship would be interesting in the context of large-scale finite element (FE) simulations because it could allow to replace the BII by an interface constitutive law (soft imperfect interface [24]) instead of considering a rigid interface behavior (perfect interface). Such an approach will allow to scale the BII microscopic properties up to the macroscopic behavior in a multiscale FE-simulation framework.
The aim of this study is to develop a multiscale model of the biomechanical behavior of an osseointegrated BII, which is modeled as an imperfect interface, taking into account its microscopic properties: the BIC ratio, the implant surface roughness and the bone properties. An effective incremental contact stiffness in normal direction is numerically derived and compared to analytical predictions.

Results
In "Voigt-Reuss bounds" section, the numerical contact stiffness K FEM c is compared with the analytical stiffness obtained from the Voigt-Reuss bounds. In "Comparisons with the analytical models" section, the results obtained with the FEM (described in "Finite element modeling" section) and the three analytical models (described in "Analytical approaches" section) are compared. Finally, the parametrical analyses investigating the effect of the BIC, of the roughness amplitude ∆, of the roughness wavelength λ and of the bone Young's modulus E b are presented in "Parametrical analyses" section. Figure 1 shows the variation of the effective contact stiffness obtained using the numerical approach described in "Numerical resolution" section as a function of the BIC in the reference configuration (∆ = 5 µm, λ = 80 µm, E b = 2 GPa), which is compared with the Voigt-Reuss bounds. As expected, the numerical results are comprised between Voigt (upper) and Reuss (lower) bounds for different BIC scenarios, which constitutes a first validation of the numerical model.   function of the BIC. In this case, the numerical stiffness has been calculated assuming ∆ = 1 µm to compare the FE solution with the analytical solutions by Lekesiz, which are valid in the case of crack-like non-contact zones (see "Lekesiz's models" section). All other parameters in the FE model are taken equal to their reference values. The effective contact stiffness is shown to increase as a function of the BIC for all three models.

Discussion
The originality of this work is to provide a homogenized interface model able to obtain an effective normal contact stiffness of the BII accounting for bone ingrowth and implant micro-roughness. The proposed 2D numerical model was validated by comparison with five analytical models. The Lekesiz's schemes were employed to describe the asymptotic limit behavior of the proposed FE-based interface model when the aspect ratio of the non-contact zone b a tends to zero, which corresponds to the case where the non-contact zones can be modeled as internal micro-cracks. The results obtained with Lekesiz's models were also employed to highlight the effect of the interaction between non-contact zones. Kachanov's scheme was used to describe the asymptotic limit behavior when the aspect ratio of the non-contact zone b a tends to 1 which corresponds to the case where the non-contact zones can be modeled as elliptical holes. The results presented herein show that the roughness parameters (∆ and λ), the bone Young's modulus Eb as well as the BIC have a significant effect on the effective contact stiffness of the BII. Figure 7 shows that the effective normal contact stiffness increases linearly as a function of the bone Young's modulus, which may be explained by considering the Hertzian contact problem between two elastic spheres made of the same isotropic material (E, ν ) [25]. In the Hertzian theory, for a circular contact area S, the incremental normal contact stiffness K c N per contact area, reads as follows [26]: Equation (1) shows that the contact stiffness K c N is proportional to the material Young's modulus E, which explains the linearity of K FEM c as a function of the bone Young's modulus E b obtained in Fig. 7. Moreover, Eq. (1) shows that the normal contact stiffness increases as a function of S, which is directly related to the BIC, thus explaining that the slope of the variation of K FEM c as a function of E b is higher for BIC = 100% compared to the case where BIC = 50%.
Voigt-Reuss bounds are obtained by a rule of mixtures as a weighted average, and they provide the theoretical widest upper-and lower bounds used to predict various properties (mechanical, thermal, etc.) of fiber-composite materials and porous materials. The effective elastic properties obtained by a homogenization technique must fall into these bounds. Accordingly, results shown in Fig. 1 constitute a first validation for the proposed homogenized numerical model. In addition, the analytical stiffness obtained by both Lekesiz's and Kachanov's schemes are comprised within Voigt-Reuss bounds as can be highlighted by comparing Figs. 1, 2, 3 and 4.
The results shown in Figs. 1, 4, 5 and 6 indicate that the normal contact stiffness obtained by FE analyses K FEM c always increases as a function of the BIC. Note that a similar behavior is obtained in the Reuss model (see Fig. 1), as well as with the Lekesiz's models (see Fig. 4). These results may be explained by the fact that increasing of the BIC leads to a more important contact area, which is known to increase contact rigidity nonlinearly (refer to the Hertzian contact model Eq. (1), for instance). Moreover, increasing the BIC leads to a decrease of the void fraction (see Fig. 8), which also leads to a more rigid interface behavior. The results shown in Figs. 2 and 5 indicate that the contact stiffness K FEM c decreases as a function of ∆, for BIC > 20%. This result may be explained by the increase of the size of the cavity in the y-direction when ∆ increases, leading to lower rigidity of the system compared to a cavity with a lower aspect ratio (defined by b a ). Note that similar results were obtained with the Kachanov's model (see Fig. 2), which may be explained by the same analysis. Other authors have also found a higher normal contact stiffness for lower values of the roughness amplitude with analytical models [27] and experimental approaches [28].
For BIC values lower than around 20%, the contact stiffness K FEM c weakly depends on the amplitude ∆, as shown in Fig. 5. In this case, the overall contact is concentrated near the peaks of the implant surface (y = 0 and y = λ/2 in Fig. 8), which present an initially horizontal slope due to the sinusoidal geometry. Therefore, for low values of BIC, such situation may be approximated by a "flat punch" configuration, which does not depend on variations of ∆ because the effective contact area does not vary significantly.
When the BIC is equal to 100%, the porosity of the BII vanishes (see Fig. 8); so, the numerical interface model acts like a bi-material layer. In this case, as shown in Fig. 5, the overall stiffness is shown to increase when ∆ decreases, as expected for a thin elastic layer [24].
The results shown in Fig. 6 indicate that the contact stiffness K FEM c decreases as a function of λ, which may be explained by the decrease of the contact area (or contact length in our 2D model) per unit length (along the x-direction) when the wavelength λ increases. This behavior can be analytically explained by the expression of the contact length per unit length L P = 1 π E − 4π 2 � 2 2 , obtained by Eq. (4) after some algebra.
The numerical and analytical results may be compared in specific configurations. First, the cases where (i) ∆ tends towards 0 and (ii) λ increases significantly both correspond to an aspect ratio of the non-contact zone b a tending towards 0. Such situation corresponds in turn to the micro-crack of the Lekesiz's model, which explains why the effective contact stiffness obtained with the FEM tends towards that obtained with the Lekesiz's model with crack interaction when ∆ tends towards 0 (see Fig. 2) as well as when λ > 200 µm (see Fig. 3). Second, the cases where (i) ∆ increases significantly and (ii) λ tends towards 0 correspond to aspect ratio of the non-contact zone b a tending towards 1. This situation corresponds in turn to Kachanov model, which explains why the effective contact stiffness obtained with the FEM tends towards that obtained with the Kachanov model for λ < 30 µm (see Fig. 3) as well as when ∆ > 15 µm (see Fig. 2). Figure 4 shows that for relatively low BIC values (below around 50%), the results obtained with the Lekesiz's model with crack interaction (I) and with the FE model are in good agreement, while the results obtained with the Lekesiz's model without crack interaction (NI) overestimates the contact stiffness obtained with the FE model. When considering high BIC values within the physiological range (50% < BIC < 75%), the results obtained with the Lekesiz's model without crack interaction (NI) and with the FE model are in good agreement. This result may be explained by the fact that the distance between two contiguous non-contact zones (i.e., − 2a ) increases as a function of the BIC (see Fig. 8). Note that for BIC > 50% the analytical solutions obtained by models with (I) and without (NI) interactions overlap because the interaction parameter I k a expressed by Eq. (16) is equal to 1. Physically, this finding is in agreement with the theory of the micro-cracked media, which establishes that the interaction between cracks could be neglected when the distance between two contiguous cracks is greater or equal to their length [52]. Therefore, when BIC values are relatively low, a suitable choice is to consider interacting cracks, while cracks interactions may be neglected for relatively high BIC values. Note that when the BIC tends towards 100% (which corresponds to a tending to 0), both Lekesiz's models becomes invalid giving an infinite contact stiffness because of the interaction parameter I k a .
Several parameters were chosen empirically. First, the choice of the thickness H = 500 µm of the layers representing bone tissue and the implant was made to find a compromise between a sufficiently low value to obtain reasonable computational time and a sufficiently high value so that the results do not depend on H for all configurations. Accordingly, a parametric study on H was carried out in the most unfavorable case that is for λ = 350 µm. The numerical contact stiffness was found not to depend on H for H > 300 µm.
Second, all numerical simulations have been carried out considering a uniform normal tensile stress σ 0 = 25 MPa (see Fig. 8), which was chosen because such amount of stress has been shown to be obtained in clinical configuration [29]. However, it has been verified that the results do not depend on the choice of the value of σ 0 , which is explained by the linear assumptions of the model. This study has several limitations. First, the FE model was developed in 2D and the BII description should be realized considering the three dimensions, which could be done using approaches such as, for example, 3D analytical approaches based on the Eshelby theory [30][31][32]. Second, a sinusoidal function is used to describe the implant surface roughness, similarly as what was done by [46]. This sinusoidal description of the implant surface constitutes a strong approximation and considering the real surface texture is likely to lead to different results. Note that comparable approaches using wavy surfaces have already been developed in contact mechanics in the past  [33]. Third, only normal stresses were considered and assessing the influence of shear stresses would be of interest since it corresponds to a situation of clinical interest. Fourth, bone material properties were assumed to be linearly elastic, homogeneous and isotropic, and fluid structure interactions were neglected, similarly as what was done in various previous FE-based numerical studies [34][35][36][37][38]. However, newly formed bone properties are heterogeneous [39] and viscoelastic, which was neglected herein. The assumption of bone homogeneity implies a simplified modeling of the bone microstructure, which may include cavities at various scales [40,41]. Note that adopted analytical approaches can be extended to anisotropy of the bone by considering microcracks and pores embedded into an anisotropic matrix. Moreover, a drymaterial model was employed and fluid-structure interaction effects were neglected; these assumptions were necessary to understand the main contact mechanisms occurring at the BII, by simplifying as coherently as possible the biomechanical problem. Note that in vivo, fully bonded interfaces are not likely to occur since values of the BIC are typically comprised between around 30 and 80% [42][43][44]. Fifth, the proposed model should be validated experimentally. However, accurate data on boneimplant contact stiffness are scarce because of the difficulty of controlling the implant surface as well as the bone distribution around the implant. However, a validation of the proposed numerical model was carried out by comparison with three analytical models (and Voigt-Reuss bounds). Sixth, adhesion phenomena at the BII [45] that may be important in particular at early stage of osseointegration were not taken into account, which should be considered in the future. Despite the aforementioned limitations, this work constitutes the first attempt to determine an effective normal contact stiffness for an osseointegrated BII, which is able to account for the microstructural effects of periprosthetic bone properties and implant micro-roughness.

Conclusions
The present study proposes a nonlinear spring model for an osseointegrated boneimplant interface in normal contact conditions based on micromechanical modeling. A 2D FE model accounting for implant micro-roughness and bone ongrowth is used to obtain an effective incremental contact stiffness K FEM c in normal direction. Several parametrical analyses have been carried out to investigate the effects of the micro-roughness parameters Δ and , the BIC ratio and the bone Young's modulus E b , on K FEM c . Comparisons with analytical schemes based on micromechanical homogenization (Eshelby's problem [58]) have been employed to validate the FE model, as well as to provide an insight on the advantages and limitations of closed-form analytical solutions for the estimation of the effective contact stiffness of BII. The proposed nonlinear spring modeling strategy for BII can allow to overcome computational difficulties to account for BII evolving microstructure within large-scale finite element (FE) simulations.
Future work should focus on 3D modeling (considering actual surface profiles), on the tangential movement and on studying the biomechanical reaction terms which mimic osseointegration phenomena.

Methods and models
Let the general framework of the contact problem be introduced by considering two continuous bodies (representing in the following implant and bone tissue) comprising linearly elastic isotropic materials, in no-sliding contact via a rough surface under a tensile loading condition (Fig. 8). Let a Cartesian frame (O, e x , e y , e z ) be introduced, with x, y and z the corresponding coordinates. Let S ⊂ R 2 be the nominal contact area between the two bodies. Then, the normal incremental contact stiffness per unit nominal contact area in S is defined as: where dw is the increment of the relative displacement at the contacting interface region in normal (i.e., along e y ) direction, and dF N is the increment of the normal force transmitted through the unit contact area.

Bone-implant interface modeling
The contact region of interest (ROI) comprises two subdomains corresponding to the bone tissue and to the implant with the same thickness H = 500 µm (along e y ). Similarly as what was done by [46], a simplistic idealization of the contacting rough implant surface via a sinusoidal wavy-like surface, of amplitude 2Δ and wavelength λ, is adopted: Since the out-of-plane dimension of the ROI (i.e., along e z ) is assumed to be infinite and applied load is homogeneous, an assumption of plane strain is considered herein, in the plane (e x , e y ).
The BIC ratio (i.e., the ratio of the bone in direct contact with the implant surface) which is assumed to be comprised between 5 and 100% depends on the bone ongrowth onto the implant surface and is an input parameter of the model. Non-contact zones between bone and the implant determine the microstructural voids of the BII. The dimension of these voids along e y (respectively, along e x ) is noted 2b (respectively, 2a) and depends on , and BIC. Note that for the proposed model, the BIC is geometrically given by: where L P (respectively, L T ) is the arc length of the implant surface in contact with bone tissue (respectively, the total arc length of the implant boundary), and E(z) = E π 2 |z and E = (z|m) = z 0 1 − m sin 2 (t)dt represent the complete and incomplete elliptic integral of the second kind [47], respectively.
The material properties of the implant and bone tissue are assumed to be linear elastic, isotropic and homogeneous. The implant is assumed to be made of titanium alloy (2) 18:114 (TiAl6V4) with a Young's modulus E i = 113 GPa [48]. The bone Young's modulus E b is varied between 1 and 20 GPa, to simulate the increase of the stiffness of the healing periprosthetic bone tissue due to osseointegration phenomena [49]. All materials are assumed to have a Poisson ratio ν of 0.3. Figure 8 shows the boundaries conditions for the proposed FE model. A symmetric boundary condition u · e x = u x = 0 holds on the boundaries parallel to e y , where u is the displacement vector. The lower boundary (y = − H) of the bone domain is fixed (u = 0) . A uniform normal tension σ 0 = 25 MPa is applied on the upper boundary of the implant domain (y = 2∆ + H) for all FE analyses. The choice of these model parameters is discussed in "Discussion" section. Continuity conditions of the displacement and of the traction (bounded interface condition) hold at the contacting surfaces between bone and implant.

Numerical resolution
The goal of the proposed FE simulations is to derive the incremental contact stiffness in the normal direction of the BII, denoted K FEM c . To do so, a homogenization model of the interphase representing the bone-implant contacting zone (see Fig. 8) was developed to derive the normal incremental contact stiffness. An equivalent one-dimensional system comprising three springs in series representing the implant domain (spring stiffness K i ), the bone domain (spring stiffness K b ) and the BII (spring stiffness K FEM c ) was considered. The stiffness of the one-dimensional bone-implant system K eq described above is determined numerically following: where σ yy = σe y · e y and u y are, respectively, the stress and displacement along e y and �•� = 2 • dx represents the average on the loaded boundary (at y = 2 + H, see Fig. 1). As a result, the numerical contact stiffness of the BII is given by: where K j = j +2µ j δ ; j = i, b (the subscript i and b denotes the implant and bone tissue, respectively). In Eq. (6), j and µ j are the Lamé parameters of the implant and bone materials and δ = H − � corresponds to the size of the implant and bone domain. Equation 6 is used to obtain the numerical contact stiffness of the BII for each FEbased simulation. Several parametrical analyses have been carried out to investigate the influence of the BIC (between 5 and 100%), of the roughness amplitude ∆ (between 0.1 and 20 µm), of the roughness wavelength λ (between 20 and 350 µm) and of the bone Young's modulus E b (between 1 and 20 GPa). The ranges of variation of ∆ and λ were chosen based on the values of the arithmetic mean roughness (R a ) and the mean spacing (S m ) measured on real titanium implants [10], noting that R a ∼ = 2�/π . The reference configuration corresponds to the following parameters: BIC = 50%, ∆ = 5 µm, λ = 80 µm and E b = 2 GPa. The numerical analyses have been carried out using the Comsol Multiphysics ® simulation software (Stockholm, Sweden). The total mesh comprises approximately 9300 second-order triangular Lagrange elements depending on the geometrical configuration of the micro-roughness parameters ∆ and λ (e.g., 9274 in the reference configuration). The global system comprises about 38,000° of freedom (e.g., 37,742 in the reference configuration). For the sake of regularity, a mapped second-order quadrangular mesh has been chosen only for the case BIC = 100%, as shown in Fig. 8. The interpolation functions for the displacement field are quadratic. The mesh size has to be finer around the tip of the non-contact zone, as shown in Fig. 8, to describe the stress localization, especially for a smaller interface thickness 2∆. To this aim, a convergence study was performed to choose the mesh size at the tip in the most unfavorable configuration that is for ∆ = 0.1 µm. The chosen mesh has a smallest element size of 3 × 10 −5 µm. The local error on the approximated solution u a y calculated in terms of the displacement u y is

Analytical approaches
The non-contact zone of the BII (Fig. 9a) may be approximated by an internal microcrack [50] when b a → 0 , as shown in Fig. 2b. Lekesiz's schemes consider an array of planar cracks at the interface between two dissimilar isotropic materials [51] under two hypotheses assuming i) interacting and ii) non-interacting micro-cracks. Conversely, the non-contact zone may be better approximated by an elliptical cavity [52] when its aspect ratio b a → 1 , as shown in Fig. 2c. Kachanov's scheme considers elliptical holes randomly oriented embedded into an isotropic homogenized material [53] in the non-interacting approximation [54].

Voigt and Reuss bounds
The Voigt and Reuss bounds represent the upper and the lower limits for the homogenized mechanical properties of a multi-phase composite material, respectively. These analytical bounds have been calculated considering the region of interest shown in Fig. 9a and constituted by the implant, bone tissue and the void corresponding to the non-contact zone (for BIC < 100%). For BIC = 100%, a two-phase medium made of bone and implant was considered. The area fractions of bone tissue, implant and void were noted f b , f i and f void , respectively, and were related to the BIC, ∆ and λ. The analytical contact stiffness given by the Voigt bound writes [24]: where κ V and µ V are the equivalent bulk and shear moduli, respectively, given as a function of the bulk moduli and shear moduli of bone tissue κ b , µ b ; implant κ i , µ i and void κ void , µ void , as follows: Similarly, the analytical contact stiffness given by the Reuss bound is given by: where κ R and µ R are the equivalent bulk modulus and shear modulus, respectively, given by: In Eqs. (8) and (10), the bulk modulus and shear modulus of the void are assumed to be equal to those of a very soft material of which the elastic moduli are taken as κ void = ϕ κ b and µ void = ϕ µ b , with ϕ = 10 −9 , following what has been done in [55]. Note that hanging the value of ϕ between 10 −12 and 10 −6 does not affect the results.

Lekesiz's models
Lekesiz's scheme [51] considers two different media separated by a partial non-contact zone modeled as an internal micro-crack of length 2a. Note that 2a = λ/2 for BIC = 50% (see Fig. 2b). Lekesiz provided a closed-form analytical expression for the effective spring stiffness of a planar periodic array of collinear cracks at the interface between two dissimilar isotropic materials. Lekesiz's scheme is based on an elastic analysis under normal loading in the framework of the open crack model [56]. Lekesiz also provided an analytical expression of the effective contact stiffness taking into