 Research
 Open Access
 Published:
Micromechanical modeling of the contact stiffness of an osseointegrated bone–implant interface
BioMedical Engineering OnLine volume 18, Article number: 114 (2019)
Abstract
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 \(\left( {K_{c}^{\text{FEM}} } \right)\) 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 \(K_{c}^{\text{FEM}}\).
Results
The model is validated by comparison with three analytical schemes based on micromechanical homogenization including two Lekesiz’s models (considering interacting and noninteracting microcracks) and a Kachanov’s model. \(K_{c}^{\text{FEM}}\) is found to be comprised between 10^{13} and 10^{15} N/m^{3} according to the properties of the BII. \(K_{c}^{\text{FEM}}\) 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, \(K_{c}^{\text{FEM}}\) 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, micromovements 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 Xray 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 Xray 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 largescale 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 FEsimulation 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_{c}^{\text{FEM}}\) 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.
Voigt–Reuss bounds
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.
Comparisons with the analytical models
Figure 2 shows the variation of the effective contact stiffness \(K_{c}\) obtained with the finite element model (gray solid line), the Kachanov’s model (black solid line), the Lekesiz’s model with crack interaction (black dashed line) and the Lekesiz’s model without crack interaction (gray dashed line) as a function of the roughness amplitude ∆. All other parameters are taken equal to their reference value. The effective contact stiffness is shown to decrease as a function of ∆ for the FEM and for the Kachanov’s model. Moreover, as expected, the results obtained with both Lekesiz’s models do not depend on ∆.
Figure 3 shows the variation of the effective contact stiffness \(K_{c}\) obtained with the finite element model (gray solid line), the Kachanov’s model (black solid line), the Lekesiz’s model with crack interaction (black dashed line) and the Lekesiz’s model without crack interaction (gray dashed line) as a function of the roughness wavelength λ. All other parameters are taken equal to their reference values. The effective contact stiffness is shown to decrease as a function of λ for all four models.
Figure 4 shows the variation of the effective contact stiffness \(K_{c}\) obtained with the finite element model (gray solid line), the Lekesiz’s model with crack interaction (black dashed line) and the Lekesiz’s model without crack interaction (gray dashed line) as a 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 cracklike noncontact 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.
Parametrical analyses
Figure 5 shows the variation of the numerical effective contact stiffness \(K_{c}^{\text{FEM}}\) as a function of the BIC for different values of ∆. The wavelength λ and the bone Young’s modulus E_{b} are equal to their reference values 80 µm and 2 GPa, respectively. \(K_{c}^{\text{FEM}}\) is shown to increase as a function of the BIC for all values of ∆ considered. The sensitivity of \(K_{c}^{\text{FEM}}\) to BIC variations is more important for lower values of ∆ and high values of the BIC. For low BIC values, \(K_{c}^{\text{FEM}}\) weakly depends on ∆; while for high BIC values, \(K_{c}^{\text{FEM}}\) decreases significantly as a function of ∆.
Figure 6 shows the variation of the numerical effective contact stiffness \(K_{c}^{\text{FEM}}\) as a function of the BIC for different values of λ. The roughness amplitude ∆ and the bone Young’s modulus E_{b} are equal to their reference values 5 µm and 2 GPa, respectively. \(K_{c}^{\text{FEM}}\) is shown to increase as a function of the BIC for all values of λ considered. Moreover, \(K_{c}^{\text{FEM}}\) decreases significantly as a function of λ for all BIC values.
Figure 7 shows the variation of the numerical effective contact stiffness \(K_{c}^{\text{FEM}}\) as a function of the bone Young’s modulus E_{b} for BIC = 50% and BIC = 100%. The roughness amplitude ∆ and the roughness wavelength λ are equal to their reference values 5 µm and 80 µm, respectively. \(K_{c}^{\text{FEM}}\) is shown to increase linearly as a function of E_{b} for all BIC values considered. The slope of the variation of \(K_{c}^{\text{FEM}}\) as a function of E_{b} is higher for BIC = 100% compared to the case where BIC = 50%.
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 microroughness. 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 FEbased interface model when the aspect ratio of the noncontact zone \({b \mathord{\left/ {\vphantom {b a}} \right. \kern0pt} a}\) tends to zero, which corresponds to the case where the noncontact zones can be modeled as internal microcracks. The results obtained with Lekesiz’s models were also employed to highlight the effect of the interaction between noncontact zones. Kachanov’s scheme was used to describe the asymptotic limit behavior when the aspect ratio of the noncontact zone \({b \mathord{\left/ {\vphantom {b a}} \right. \kern0pt} a}\) tends to 1 which corresponds to the case where the noncontact 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, \(\nu\)) [25]. In the Hertzian theory, for a circular contact area S, the incremental normal contact stiffness \(K_{N}^{c}\) per contact area, reads as follows [26]:
Equation (1) shows that the contact stiffness \(K_{N}^{c}\) is proportional to the material Young’s modulus E, which explains the linearity of \(K_{c}^{\text{FEM}}\) 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_{c}^{\text{FEM}}\) 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 fibercomposite 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_{c}^{\text{FEM}}\) 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_{c}^{\text{FEM}}\) decreases as a function of ∆, for BIC > 20%. This result may be explained by the increase of the size of the cavity in the ydirection when ∆ increases, leading to lower rigidity of the system compared to a cavity with a lower aspect ratio (defined by \({b \mathord{\left/ {\vphantom {b a}} \right. \kern0pt} 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_{c}^{\text{FEM}}\) 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 bimaterial 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_{c}^{\text{FEM}}\) 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 xdirection) when the wavelength λ increases. This behavior can be analytically explained by the expression of the contact length per unit length \(\frac{{L_{P} }}{\lambda } = \frac{1}{\pi }E\left( {  \frac{{4\pi^{2} \Delta^{2} }}{{\lambda^{2} }}} \right)\), 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 noncontact zone \({b \mathord{\left/ {\vphantom {b a}} \right. \kern0pt} a}\) tending towards 0. Such situation corresponds in turn to the microcrack 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 noncontact zone \({b \mathord{\left/ {\vphantom {b a}} \right. \kern0pt} 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 noncontact zones (i.e., \(\lambda  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} \left( {\frac{a}{\lambda }} \right)\) expressed by Eq. (16) is equal to 1. Physically, this finding is in agreement with the theory of the microcracked 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} \left( {\frac{a}{\lambda }} \right)\).
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 \(\sigma_{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 \(\sigma_{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 FEbased 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 bone–implant 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 microroughness.
Conclusions
The present study proposes a nonlinear spring model for an osseointegrated bone–implant interface in normal contact conditions based on micromechanical modeling. A 2D FE model accounting for implant microroughness and bone ongrowth is used to obtain an effective incremental contact stiffness \(K_{c}^{\text{FEM}}\) in normal direction. Several parametrical analyses have been carried out to investigate the effects of the microroughness parameters Δ and \(\lambda\), the BIC ratio and the bone Young’s modulus E_{b}, on \(K_{c}^{\text{FEM}}\). 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 closedform 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 largescale 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 nosliding 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 \subset {\mathbb{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.
Finite element modeling
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 wavylike surface, of amplitude 2Δ and wavelength λ, is adopted:
Since the outofplane 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. Noncontact 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 \(\Delta ,\lambda\) 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\left( z \right) = E\left( {\frac{\pi }{2}z} \right)\) and \(E = \left( {zm} \right) = \int_{0}^{z} {\sqrt {1  m\sin^{2} \left( t \right)} {\text{d}}t}\) 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 (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 \({\mathbf{u}} \cdot {\mathbf{e}}_{{\mathbf{x}}} = u_{x} = 0\) holds on the boundaries parallel to e_{y}, where \({\mathbf{u}}\) is the displacement vector. The lower boundary (y = − H) of the bone domain is fixed \(({\mathbf{u}} = 0)\). A uniform normal tension \(\sigma_{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_{c}^{\text{FEM}}\). 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 onedimensional 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_{c}^{\text{FEM}}\)) was considered. The stiffness of the onedimensional bone–implant system K_{eq} described above is determined numerically following:
where \(\sigma_{yy} = \left( {{\varvec{\upsigma}} {\mathbf{e}}_{{\mathbf{y}}} } \right) \cdot {\mathbf{e}}_{{\mathbf{y}}}\) and \({\text{u}}_{y}\) are, respectively, the stress and displacement along e_{y} and \(\left\langle \bullet \right\rangle = \frac{2}{\lambda }\int { \bullet \,{\text{d}}x}\) represents the average on the loaded boundary (at y = \(2\Delta\) + H, see Fig. 1). As a result, the numerical contact stiffness of the BII is given by:
where \(K_{j} = \frac{{\lambda_{j} + 2\mu_{j} }}{\delta }; \; j = i,b\) (the subscript i and b denotes the implant and bone tissue, respectively). In Eq. (6), \(\lambda_{j}\) and \(\mu_{j}\) are the Lamé parameters of the implant and bone materials and \(\delta = H  \Delta\) 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 \({\text{R}}_{\text{a}} \cong 2\Delta /\pi\). 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 secondorder triangular Lagrange elements depending on the geometrical configuration of the microroughness 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 secondorder 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 noncontact 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 \({\text{u}}_{y}^{a}\) calculated in terms of the displacement \({\text{u}}_{y}\) is \(\frac{{\left\ {{\text{u}}_{y}^{e}  {\text{u}}_{y}^{a} } \right\_{{L^{2} }} }}{{\left\ {{\text{u}}_{y}^{e} } \right\_{{L^{2} }} }} = 0. 0 0 6 {\text{\%,}}\) where \(\left\ \bullet \right\_{{L^{2} }}\) indicates the L^{2} norm.
Analytical approaches
The noncontact zone of the BII (Fig. 9a) may be approximated by an internal microcrack [50] when \({b \mathord{\left/ {\vphantom {b a}} \right. \kern0pt} a} \to 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) noninteracting microcracks.
Conversely, the noncontact zone may be better approximated by an elliptical cavity [52] when its aspect ratio \({b \mathord{\left/ {\vphantom {b a}} \right. \kern0pt} a} \to 1\), as shown in Fig. 2c. Kachanov’s scheme considers elliptical holes randomly oriented embedded into an isotropic homogenized material [53] in the noninteracting 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 multiphase 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 noncontact zone (for BIC < 100%). For BIC = 100%, a twophase 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_{\text{void}}\), respectively, and were related to the BIC, ∆ and λ. The analytical contact stiffness given by the Voigt bound writes [24]:
where \(\kappa_{V}\) and \(\mu_{V}\) are the equivalent bulk and shear moduli, respectively, given as a function of the bulk moduli and shear moduli of bone tissue \(\kappa_{b} , \mu_{b}\); implant \(\kappa_{i} , \mu_{i}\) and void \(\kappa_{\text{void}} , \mu_{\text{void}}\), as follows:
Similarly, the analytical contact stiffness given by the Reuss bound is given by:
where \(\kappa_{R}\) and \(\mu_{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 \(\kappa_{\text{void}} = \varphi \, \kappa_{b}\) and \(\mu_{\text{void}} = \varphi \,\mu_{b}\), with \(\varphi\) = 10^{−9}, following what has been done in [55]. Note that hanging the value of \(\varphi\) 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 noncontact zone modeled as an internal microcrack of length 2a. Note that 2a = λ/2 for BIC = 50% (see Fig. 2b). Lekesiz provided a closedform 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 account the effects of the interaction between cracks. Two models taken from [51] were considered in this study.
Interacting microcracks When the distance between noncontact zones is small compared to λ, these noncontact zones can be assumed to interact with each other. In this case, the effective contact stiffness \(K_{{c, {\text{int}}}}^{\text{LEK}}\) is:
where the elastic dissimilarity function \(M_{k} \left( {\alpha , \beta } \right)\) is expressed in terms of the Dundurs’ parameters \(\alpha\) and \(\beta\) within the assumption of plane strain [57] following:
where \(\in\) is the oscillation index and \(E_{j}^{*} = \frac{{E_{j} }}{{1  \nu^{2} }}; \quad j = i,b\). The interaction parameter \(I_{k} \left( {\frac{a}{\lambda }, \in } \right)\) can be approximated by the homogeneous case \(I_{k} \left( {\frac{a}{\lambda }} \right)\) when \(\left \in \right < 0.05\). In the proposed model, the oscillation index is equal to \(\left \in \right = 0.039\). Therefore, the interaction parameter is:
where \(\frac{a}{\lambda }\) represents the linear crack density, as sketched in Fig. 9.
The parameter \(K_{{c ,\,{\text{nonint}}}}^{\text{homog}}\) in Eq. (11) corresponds to the contact stiffness in the case of noninteracting cracks embedded in a homogeneous domain, assumed to be the implant material, and reads:
where \(\gamma_{i} = 3  4\nu\) is the Kolosov’s constant for the implant material.
Noninteracting microcracks When the distance between noncontact zones is relatively large compared to λ, the interaction parameter \(I_{k} \left( {\frac{a}{\lambda }} \right)\) may be taken equal to 1, which leads to the effective contact stiffness \(K_{{c , \,{\text{nonint}}}}^{\text{LEK}}\) in the noninteraction case:
Kachanov’s model
As shown in Fig. 9c, in the Kachanov’s model, the noncontact zone is modeled by a circumscribed elliptical cavity with \(a\) and \(b\) being the semimajor and semiminor axes, respectively. This model neglects interactions between voids. The aforementioned elliptical cavity is embedded in a homogenized medium whose mechanical properties are given by the homogenization of the isotropic mechanical properties of implant and bone materials [58]. The closedform analytical expression for the effective Young’s modulus \(E_{\text{eff}}\) of a twodimensional isotropic matrix with an embedded elliptical hole randomly oriented is [53]:
where the microstructural parameters p and q are the porosity and the average eccentricity over the bone–implant area \(A = 2\Delta \lambda\) and are expressed by:
In Eq. (19), E_{0} corresponds to the equivalent Young’s modulus of the homogenized matrix in which the elliptical hole is embedded and based on plane strain assumption, it is given by:
where \(E_{i}^{*}\) and \(E_{b}^{*}\) are defined above. The effective contact stiffness of the BII obtained with the Kachanov’s model is given by:
Availability of data and materials
The datasets used, generated and/or analyzed during the current study are available from the corresponding author on reasonable request.
Abbreviations
 BII:

bone–implant interface
 FE(M):

finite elements (method)
 BIC:

bone–implant contact ratio
 ROI:

region of interest
 I:

interacting microcracks (Lekesiz’s scheme)
 NI:

noninteracting microcracks (Lekesiz’s scheme)
References
 1.
Haiat G, Wang HL, Brunski J. Effects of biomechanical properties of the boneimplant interface on dental implant stability: from in silico approaches to the patient’s mouth. Annu Rev Biomed Eng. 2014;16:187–213. https://doi.org/10.1146/annurevbioeng071813104854.
 2.
Abrahamsson I, Berglundh T, Linder E, Lang NP, Lindhe J. Early bone formation adjacent to rough and turned endosseous implant surfaces. An experimental study in the dog. Clin Oral Implant Res. 2004;15(4):381–92.
 3.
Berglundh T, Abrahamsson I, Albouy JP, Lindhe J. Bone healing at implants with a fluoridemodified surface: an experimental study in dogs. Clin Oral Implant Res. 2007;18(2):147–52.
 4.
Albrektsson T, Johansson C. Osteoinduction, osteoconduction and osseointegration. Eur Spine J. 2001;10(2):S96–101.
 5.
Wennerberg A, Jimbo R, Stübinger S, Obrecht M, Dard M, Berner S. Nanostructures and hydrophilicity influence osseointegration: a biomechanical study in the rabbit tibia. Clin Oral Implant Res. 2014;25(9):1041–50.
 6.
Karachalios T. Boneimplant interface in orthopedic surgery: basic science to clinical applications. Berlin: Springer; 2013.
 7.
Götz HE, Müller M, Emmel A, Holzwarth U, Erben RG, Stangl R. Effect of surface finish on the osseointegration of lasertreated titanium alloy implants. Biomaterials. 2004;25(18):4057–64.
 8.
Rønold HJ, Ellingsen JE. Effect of microroughness produced by TiO2 blasting—tensile testing of bone attachment by using coinshaped implants. Biomaterials. 2002;23(21):4211–9.
 9.
Heller AL, Heller RL. Clinical evaluations of a poroussurfaced endosseous implant system. J Oral Implantol. 1996;22(3–4):240–6.
 10.
Rangert B, Sennerby L, Meredith N, Brunski J. Design, maintenance and biomechanical considerations in implant placement. Dental Update. 1997;24(10):416–20.
 11.
Pilliar RM, Lee JM, Maniatopoulos C. Observations on the effect of movement on bone ingrowth into poroussurfaced implants. Clin Orthop Relat Res. 1986;208:108–13.
 12.
Søballe K, Hansen ES, Rasmussen H, Jorgensen PH, Bunger C. Tissue ingrowth into titanium and hydroxyapatitecoated implants during stable and unstable mechanical conditions. J Orthop Res. 1992;10(2):285–99.
 13.
Duyck J, Vandamme K, Geris L, Van Oosterwyck H, De Cooman M, Vandersloten J, Puers R, Naert I. The influence of micromotion on the tissue differentiation around immediately loaded cylindrical turned titanium implants. Arch Oral Biol. 2006;51(1):1–9.
 14.
Orlik J, Zhurov A, Middleton J. On the secondary stability of coated cementless hip replacement: parameters that affected interface strength. Med Eng Phys. 2003;25(10):825–31.
 15.
Michel A, Bosc R, Vayron R, Haiat G. Assessing the acetabular cup implant primary stability by impact analysis: a cadaveric study. PLoS ONE. 2016;11(11):e0166778.
 16.
Vayron R, Soffer E, Anagnostou F, Haïat G. Ultrasonic evaluation of dental implant osseointegration. J Biomech. 2014;47(14):3562–8.
 17.
Albrektsson T, Dahl E, Enbom L, Engevall S, Engquist B, Eriksson AR, Feldmann G, Freiberg N, Glantz PO, Kjellman O. Osseointegrated oral implants: a Swedish multicenter study of 8139 consecutively inserted Nobelpharma implants. J Periodontol. 1988;59(5):287–96.
 18.
Hecht S, Adams WH, Narak J, Thomas WB. Magnetic resonance imaging susceptibility artifacts due to metallic foreign bodies. Vet Radiol Ultrasound. 2011;52(4):409–14.
 19.
Nguyen VH, Rosi G, Naili S, Michel A, Raffa ML, et al. Influence of anisotropic bone properties on the biomechanical behavior of the acetabular cup implant: a multiscale finite element study. Comput Methods Biomech Biomed Eng. 2017;20(12):1312–25.
 20.
Du J, Lee JH, Jang AT, Gu A, HossainiZadeh M, Prevost R, Curtis DA, Ho SP. Biomechanics and strain mapping in bone as related to immediatelyloaded dental implants. J Biomech. 2015;48(12):3486–94.
 21.
Mao Q, Su K, Zhou Y, HossainiZadeh M, Lewis GS, Du J. Voxelbased microfinite element analysis of dental implants in a human cadaveric mandible: tissue modulus assignment and sensitivity analyses. J Mech Behav Biomed Mater. 2019;94:229–37.
 22.
Lutz A, Nackenhorst U. Numerical investigations on the osseointegration of uncemented endoprostheses based on bioactive interface theory. Comput Mech. 2012;50(3):367–81.
 23.
Egan J, Marsden D. A spring network model for the analysis of load transfer and tissue reactions in intramedullary fixation. Clin Biomech. 2001;16(1):71–9.
 24.
Benveniste Y, Miloh T. Imperfect soft and stiff interfaces in twodimensional elasticity. Mech Mater. 2001;33(6):309–23.
 25.
Mindlin R. Compliance of elastic bodies in contact. J Appl Mech. 1949;16:259–68.
 26.
Sevostianov I, Kachanov M. Contact of rough surfaces: a simple model for elasticity, conductivity and crossproperty connections. J Mech Phys Solids. 2008;56(4):1380–400.
 27.
Abbakumov K. Scattering of plane elastic waves on a microrough interface between solid media. Russ J Nondestr Test. 2017;53(7):475–84.
 28.
Shi X, Polycarpou AA. Measurement and modeling of normal contact stiffness and contact damping at the meso scale. J Vib Acoust. 2005;127(1):52–60.
 29.
Sadegh A, Luo G, Cowin S. Bone ingrowth: an application of the boundary element method to bone remodeling at the implant interface. J Biomech. 1993;26(2):167–82.
 30.
Pensée V, Kondo D, Dormieux L. Micromechanical analysis of anisotropic damage in brittle materials. J Eng Mech. 2002;128(8):889–97.
 31.
Deude V, Dormieux L, Kondo D, Maghous S. Micromechanical approach to nonlinear poroelasticity: application to cracked rocks. J Eng Mech. 2002;128(8):848–55.
 32.
Pichler B, Hellmich C, Mang AH. A combined fracturemicromechanics model for tensile strainsoftening in brittle materials, based on propagation of interacting microcracks. Int J Numer Anal Methods Geomech. 2007;31(2):111–32.
 33.
Ciavarella M. Adhesive rough contacts near complete contact. Int J Mech Sci. 2015;104:104–11.
 34.
Baggi L, Cappelloni I, Maceri F, Vairo G. Stressbased performance evaluation of osseointegrated dental implants by finiteelement simulation. Simul Model Pract Theory. 2008;16(8):971–87.
 35.
Ghosh R, Gupta S. Bone remodelling around cementless composite acetabular components: the effects of implant geometry and implant–bone interfacial conditions. J Mech Behav Biomed Mater. 2014;32:257–69.
 36.
Haïat G, Naili S, Grimal Q, Talmant M, Desceliers C, Soize C. Influence of a gradient of material properties on ultrasonic wave propagation in cortical bone: application to axial transmission. J Acousti Soc Am. 2009;125(6):4043–52.
 37.
Rourke DO, AlDirini R, Taylor M. Primary stability of a cementless acetabular cup in a cohort of patientspecific finite element models. J Orthop Res. 2017;36(3):1012–23.
 38.
Michel A, Nguyen VH, Bosc R, Vayron R, Hernigou P, Naili S, Haiat G. Finite element model of the impaction of a pressfitted acetabular cup. Med Biol Eng Compu. 2017;55(5):781–91.
 39.
Vayron R, Matsukawa M, Tsubota R, Mathieu V, Barthel E, Haiat G. Evolution of bone biomechanical properties at the micrometer scale around titanium implant as a function of healing time. Phys Med Biol. 2014;59(6):1389.
 40.
Cowin SC. Bone mechanics handbook. Boca Raton: CRC Press; 2001.
 41.
Fritsch A, Hellmich C. Universal’ microstructural patterns in cortical and trabecular, extracellular and extravascular bone materials: micromechanicsbased prediction of anisotropic elasticity. J Theor Biol. 2007;244(4):597–620.
 42.
Mathieu V, Vayron R, Soffer E, Anagnostou F, Haiat G. Influence of healing time on the ultrasonic response of the boneimplant interface. Ultrasound Med Biol. 2012;38(4):611–8.
 43.
Pontes AEF, Ribeiro FS, Iezzi G, Pires JR, Zuza EP, Piattelli A, Marcantonio E. Boneimplant contact around crestal and subcrestal dental implants submitted to immediate and conventional loading. Sci World J. 2014;2014:606947.
 44.
Scarano A, Degidi M, Iezzi G, Petrone G, Piattelli A. Correlation between implant stability quotient and boneimplant contact: a retrospective histological and histomorphometrical study of seven titanium implants retrieved from humans. Clin Implant Dentist Relat Res. 2006;8(4):218–22.
 45.
Rojek J, Telega JJ. Numerical simulation of boneimplant systems using a more realistic model of the contact interfaces with adhesion. J Theor Appl Mech. 1999;1999:659–86.
 46.
Heriveaux Y, Nguyen VH, Haiat G. Reflection of an ultrasonic wave on the boneimplant interface: a numerical study of the effect of the multiscale roughness. J Acoust Soc Am. 2018;144(1):488–99.
 47.
Byrd PF, Friedman MD. Handbook of elliptic integrals for engineers and physicists, vol. 67. Berlin: Springer; 2013.
 48.
Mishra AK, Davidson JA, Poggie RA, Kovacs P, FitzGerald TJ. Mechanical and tribological properties and biocompatibility of diffusion hardened Ti13Nb13Zr—a new titanium alloy for surgical implants., Medical applications of titanium and its alloys: the material and biological issuesWest Conshohocken: ASTM International; 1996.
 49.
Wazen RM, Currey JA, Guo H, Brunski JB, Helms JA, Nanci A. Micromotioninduced strain fields influence early stages of repair at boneimplant interfaces. Acta Biomater. 2013;9(5):6663–74.
 50.
Raffa ML, Lebon F, Vairo G. Normal and tangential stiffnesses of rough surfaces in contact via an imperfect interface model. Int J Solids Struct. 2016;87:245–53.
 51.
Lekesiz H, Katsube N, Rokhlin SI, Seghi RR. Effective spring stiffness for a planar periodic array of collinear cracks at an interface between two dissimilar isotropic materials. Mech Mater. 2011;43(2):87–98.
 52.
Kachanov M. Elastic solids with many cracks and related problems. Adv Appl Mech. 1993;30:259–445.
 53.
Tsukrov I, Kachanov M. Effective moduli of an anisotropic material with elliptical holes of arbitrary orientational distribution. Int J Solids Struct. 2000;37(41):5919–41.
 54.
Sevostianov I, Kachanov M. Noninteraction approximation in the problem of effective properties. In: Kachanov M, Sevostianov I, editors. Effective properties of heterogeneous materials. Netherlands: Springer; 2013. p. 1–95.
 55.
Andreassen E, Andreasen CS. How to determine composite material properties using numerical homogenization. Comput Mater Sci. 2014;83:488–95.
 56.
England A. A crack between dissimilar media. J Appl Mech. 1965;32(2):400–2.
 57.
Dundurs J. Effect of elastic constants on stress in a composite under plane deformation. J Compos Mater. 1967;1(3):310–22.
 58.
Eshelby JD. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc R Soc Lond A. 1957;241(1226):376–96.
Acknowledgements
Not applicable.
Funding
This work was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No 682001, project ERC Consolidator Grant 2015 BoneImplant).
Author information
Affiliations
Contributions
MLR performed the computation and wrote the first draft. VHN worked on the computation and the study design. GH supervised the work. All authors read and approved the final manuscript.
Corresponding author
Correspondence to Guillaume Haiat.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Raffa, M.L., Nguyen, V. & Haiat, G. Micromechanical modeling of the contact stiffness of an osseointegrated bone–implant interface. BioMed Eng OnLine 18, 114 (2019). https://doi.org/10.1186/s1293801907333
Received:
Accepted:
Published:
Keywords
 Bone–implant interface
 Contact
 Roughness
 Homogenization
 Finite element modeling