Modeling the debonding process of osseointegrated implants due to coupled adhesion and friction

Abstract

Cementless implants have become widely used for total hip replacement surgery. The long-term stability of these implants is achieved by bone growing around and into the rough surface of the implant, a process called osseointegration. However, debonding of the bone–implant interface can still occur due to aseptic implant loosening and insufficient osseointegration, which may have dramatic consequences. The aim of this work is to describe a new 3D finite element frictional contact formulation for the debonding of partially osseointegrated implants. The contact model is based on a modified Coulomb friction law by Immel et al. (2020), that takes into account the tangential debonding of the bone-implant interface. This model is extended in the direction normal to the bone-implant interface by considering a cohesive zone model, to account for adhesion phenomena in the normal direction and for adhesive friction of partially bonded interfaces. The model is applied to simulate the debonding of an acetabular cup implant. The influence of partial osseointegration and adhesive effects on the long-term stability of the implant is assessed. The influence of different patient- and implant-specific parameters such as the friction coefficient $\mu _\text {b}$, the trabecular Young’s modulus $E_\text {b}$, and the interference fit $I\!F$ is also analyzed, in order to determine the optimal stability for different configurations. Furthermore, this work provides guidelines for future experimental and computational studies that are necessary for further parameter calibration.

1 Introduction

Uncemented acetabular cup implants have become the gold standard for acetabular replacement (Small et al. 2013) in the context of total hip arthroplasty (Philippot et al. 2009). However, aseptic implant loosening is a major complication of total hip arthroplasty and is mostly due to issues related to implant stability. The primary (or initial) stability of acetabular cup implants is established during the surgery and is governed by mechanical factors, such as surface interlocking and bone quality. Secondary (or long-term) stability is achieved several weeks or months after surgery, through the formation and maturation of newly formed bone tissue at the bone-implant interface, a process called osseointegration (Albrektsson et al. 1981). While the evolution of secondary implant stability is governed by complex biochemical processes, the mechanical behavior of the bone-implant interface remains crucial for the surgical outcome (Gao and Haïat 2019).

Pull-out, lever-out, and torsional debonding tests have been widely used to evaluate the implant fixation by recording the force–displacement curve, the maximum removal force or the shear strength of the bone-implant interface (Søballe 1993; Brunski et al. 2000; Chang et al. 2010; Trisi et al. 2011; Mathieu et al. 2012a), which have been correlated with histological assessments in animal studies (Johansson and Albrektsson 1991; Haïat et al. 2014). A simple animal model consisting in considering coin-shaped implants has been used by different groups (Rønold and Ellingsen 2002; Mathieu et al. 2012a; Gao and Haïat 2019) in order to investigate the evolution of the secondary stability of osseointegrated implants. The interest of such an approach is to work in a standardized configuration with a planar bone-implant interface, allowing to carry out efficient and precise experimental testing. However, the influence of biological as well as mechanical factors on the long-term stability and the size and shape of joint implants, makes experimental testing of this kind of implants difficult and at present, such studies are lacking in the literature (Viceconti et al. 2004; Helgason et al. 2008). Therefore, there is a high demand in reliable numerical models that can predict long-term stability of implants to aid in implant conception but also to choose the best implantation technique in a patient specific manner. The advantage of numerical models is that all parameters can be precisely controlled, which is not the case when working with animal models. The difficulty of predicting implant stability arises from the complex nature of the bone-implant interface, which is related to (1) the implant surface roughness, (2) the non-homogeneous contact between bone and implant, (3) local adhesion and friction, (4) the time dependence of peri-implant bone properties (Mathieu et al. 2012a), and (5) the widely varying loading scenarios 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 (Shirazi-Adl et al. 1993).

Most finite element (FE) models that study implant behavior assume the bone-implant interface to be either perfectly bonded or fully sliding (Gupta et al. 2010; Tomaszewski et al. 2010; Chou et al. 2014; Demenko et al. 2016; Rittel et al. 2017; Mondal and Ghosh 2019). While perfectly bonded contact conditions do not allow for debonding to occur, interface behavior that is modeled as frictionless or by classical Coulomb’s friction cannot fully represent the nonlinear interface behavior of the bone-implant interface before nor after osseointegration occurs (Dammak et al. 1997; Viceconti et al. 2004). Furthermore, it was shown that the implants are never fully osseointegrated, rather the bone-to-implant contact ratio is typically only between 30 and 70 % after healing (Brånemark et al. 1997; Marin et al. 2010). Therefore, imperfect osseointegration and its influence on the evolution of implant stability must be considered. A common approach is to model imperfect osseointegration by setting osseointegrated contact regions to be perfectly bonded, while not integrated contact regions follow frictionless or Coulomb’s friction contact (Spears et al. 2001; Viceconti et al. 2004; Helgason et al. 2009; Galloway et al. 2013). Another approach is to consider a varying degree of osseointegration and to adjust the material properties of the bone-implant interface, while keeping the interface fully bonded (Kurniawan et al. 2012) or by varying the friction coefficient of the bone-implant interface from $\mu =0$ for fully unbonded to $\mu =\infty$ for perfectly osseointegrated surfaces (Korabi et al. 2017). However, although interesting, these models cannot represent the adaptive changes of the bone-implant interface, and debonding is usually characterized by excessive stress or strain at the bone-implant interface or in the bone, without modeling the actual separation between bone and implant and local changes of contact conditions.

To the best of the authors’ knowledge, there is currently no other continuum mechanics model that i) takes into account friction and adhesion in normal and tangential direction without using fully bonded interfaces, ii) takes into account imperfect/partial osseointegration, iii) models adhesive debonding in normal and tangential direction, and iv) was used to model implant debonding on the macroscopic scale. So far, there is only one similar model we know of, which is the cohesive debonding model of Rittel et al. (2018) that was used to model the debonding of partially integrated dental implants. There, a tie constraint was applied to the bone-implant interface, so that cohesive failure occurs in the bone tissue around the bone-implant interface. Partial osseointegration was modeled by defining a relative osseointegrated area with random distribution and considering non-integrated areas in frictional contact. However, comparison is difficult as there are considerable differences between adhesive and cohesive debonding, the choice of geometries, and contact conditions.

The aim of this work is to propose a phenomenological model for the frictional contact behavior of debonding osseointegrated implants. In the previous work by our group (Immel et al. 2020), the adhesive failure and tangential debonding of partially osseointegrated coin-shaped implants was modeled by proposing a modified Coulomb’s law calibrated from experimental data. In this work, this modified Coulomb’s friction law for tangential debonding is first extended by a cohesive zone model in order to account for debonding in the normal direction as well as adhesive friction. Second, this updated contact model is demonstrated on an osseointegrated, coin-shaped implant using mode I, III, and mixed mode debonding. Third, the modified Coulomb’s friction law is then applied to simulate the debonding of a 3D, osseointegrated acetabular cup implant in different removal tests. The implant stability is quantified by assessing the removal force/torque and the biomechanical determinants of the long-term stability, such as primary stability and degree of osseointegration, are assessed.

The remainder of this work is structured as follows: Section 2 describes the governing equations and the contact formulation for normal and tangential debonding. In Section 3, the demonstration with a coin-shaped implant is presented. The extended contact model is then applied to simulate removal experiments of osseointegrated acetabular cup implants in Section 4. Section 5 discusses the results obtained with the modified Coulomb’s law and its new extension and provides perspectives and guidelines for future experimental and numerical studies. Last, Section 6 gives a conclusion and an outlook to possible extensions.

2 Models and Methods

This section discusses the governing equations describing the contact behavior, including a summary of the modified Coulomb’s friction law (Immel et al. 2020) and its extension in normal direction. The resulting equations are discretized within a finite element framework to obtain a numerical solution. The readers are referred to Sauer and De Lorenzis (2013, 2015), Duong and Sauer (2019) and Immel et al. (2020) for a more detailed derivation of the considered contact formulation and its finite element implementation.

2.1 Material law

Throughout this work, a hyper-elastic Neo-Hookean material model is used for all bodies. The stress–strain relation for the Cauchy stress $\varvec{\sigma }$ of this model is defined by (Zienkiewicz and Taylor 2005)

$$\begin{aligned} \varvec{\sigma } = \frac{\Lambda }{J} \left( \mathrm {ln} J \right) {\varvec{I}} + \frac{G}{J} \left( {\varvec{b}} - {\varvec{I}} \right) , \end{aligned}$$

(1)

with the volume change J, the identity tensor ${\varvec{I}}$, and the left Cauchy–Green tensor ${\varvec{b}}$. The Lamé parameters G (shear modulus) and $\Lambda$ can be expressed in terms of Young’s modulus E and Poisson's ratio $\nu$, by

$$\begin{aligned} G = \frac{E}{2(1+\nu )} \quad \mathrm {and} \quad \Lambda = \frac{2G \nu }{1 - 2\nu }. \end{aligned}$$

(2)

2.2 Contact kinematics and contact traction

The contact traction ${\varvec{t}}_\text {c}$ can be decomposed into a normal and tangential component, i.e.,

$$\begin{aligned} {\varvec{t}}_\text {c} = {\varvec{t}}_\text {n} + {\varvec{t}}_\text {t}. \end{aligned}$$

(3)

The normal traction is proportional to the negative contact pressure p, i.e.,

$$\begin{aligned} {\varvec{t}}_\text {n} = -p \, {\varvec{n}} \end{aligned}$$

(4)

where ${\varvec{n}}$ is the outward surface normal. The contact pressure p can be modeled by the penalty method, i.e.,

$$\begin{aligned} p = -\epsilon _\text {n} g_\text {n}, \quad \end{aligned}$$

(5)

which is active when the bodies penetrate, i.e., the normal gap $g_\text {n} = {\varvec{g}}_\text {n}\cdot {\varvec{n}}$ becomes negative ($g_\text {n} < 0$). Here, ${\varvec{g}}_\text {n}$ is the normal contact gap vector between surface points on the two bodies, and $\epsilon _\text {n}$ is the penalty parameter.

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

$$\begin{aligned} f_\text {s} {\left\{ \begin{array}{ll} <0, &{} \text {for sticking},\\ =0, &{} \text {for sliding}. \end{array}\right. } \end{aligned}$$

(6)

It can be formulated as

$$\begin{aligned} f_\text {s} = \Vert {\varvec{t}}_\text {t} \Vert - t^\text {slide}_\text {t}, \end{aligned}$$

(7)

where $t^\text {slide}_\text {t}>0$ is a limit value for the tangential traction that corresponds to the magnitude of the tangential traction during sliding. 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. One of the simplest but most commonly used sliding laws is Coulomb’s friction law, which states

$$\begin{aligned} t^\text {slide}_\text {t} = \mu \, p, \end{aligned}$$

(8)

where the friction coefficient $\mu$ is a material parameter that can change locally ($\mu = \mu ({\varvec{x}})$), but does not depend on p or $g_\text {n}$.

2.2.1 State function

To model the current bonding state at each point ${\varvec{x}}$ on the interface and its transition from fully bonded to fully broken, we use a smooth state function $\phi$ (Immel et al. 2020). This state function depends on the initial bonding state $\phi _0$, the model parameters $a_\text {s}, b_\text {s}$ and a damage parameter $g_\text {d}({\varvec{x}})$, i.e.,

$$\begin{aligned} \displaystyle \begin{aligned}&\phi ({\varvec{x}}, g_\text {d}({\varvec{x}})) = \phi _0({\varvec{x}}) \cdot \\&{\left\{ \begin{array}{ll} 1, &{} c_1 < 1, \\ \frac{1}{2} - \frac{1}{2} \mathrm {sin} \left( \frac{\pi }{2b_\text {s}} \left( c_1 - b_\text {s} - 1 \right) \right) , &{} 1 \le c_1 \le c_2,\\ 0, &{} c_1 > c_2, \end{array}\right. } \end{aligned} \end{aligned}$$

(9)

where $c_1=\frac{g_\text {d}({\varvec{x}})}{a_\text {s}}$ and $c_2=1+2b_\text {s}$. Here, $a_\text {s}$ denotes a displacement threshold, while $b_s$ defines the size of the transition zone from a fully bonded to a fully debonded state. The initial bonding state $\phi _0$ can range between 0 (no initial bonding) and 1 (full initial bonding). A value in between indicates imperfect osseointegration, i.e., $\phi _0=0.5$ denotes that the level of osseointegration has reached 50%.

The damage parameter $g_\text {d}({\varvec{x}})$ is chosen to be additively composed of an accumulated irreversible tangential slip $g_\text {s}$ and the accumulated normal gap $g_\text {sn}$, i.e.,

$$\begin{aligned} g_\text {d}= g_\text {s} + g_\text {sn}. \end{aligned}$$

(10)

Numerically, $g_\text {s}$ and $g_\text {sn}$ can be approximated as (see Fig. 1 and Immel et al. (2020), Appendix A.2)

$$\begin{aligned} \begin{aligned} g_\text {s}^{n+1} \approx&\sum _{i=1}^{n+1} \Vert {\varvec{x}}_\ell \left( \hat{\varvec{\xi }}^i \right) - {\varvec{x}}_\ell \left( \hat{\varvec{\xi }}^{i-1} \right) \Vert ,\\ g_\text {sn}^{n+1} \approx&\sum _{i=1}^{n+1} \Vert g^{i}_\text {en} - g^{i-1}_\text {en} \Vert , \end{aligned} \end{aligned}$$

(11)

where

$$\begin{aligned} \begin{aligned} g^i_\text {en}&:= {\varvec{g}}^i_\text {en} \cdot {\varvec{n}}, \\ {\varvec{g}}^i_\text {en}&:= ({\varvec{n}} \otimes {\varvec{n}}) \, {\varvec{g}}^i_\text {e},\\ {\varvec{g}}^i_\text {e}&:= {\varvec{x}}^i_k - {\varvec{x}}^i_\ell (\hat{\varvec{\xi }}). \end{aligned} \end{aligned}$$

(12)

Here, subscripts i and $i-1$ denote the current and previous quantities at time step $t_i$, analogously to the quantities at $t_n+1$ shown in Figure 1. By definition, $g_\text {s}$ and $g_\text {sn}$ are always positive and monotonically increasing, even when the loading reverses. The modified Coulomb’s law only depends on $g_\text {s}$, while the extension also takes into account $g_\text {sn}$. This implies that during sticking, $g_\text {d}$ only increases if $g_\text {sn}$ increases. Otherwise, there is no change in the debonding state.

As this is a local model, where $\phi$ 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. (9), the state variable $\phi$ determines whether a point is in an unbroken ($\phi =1$), partially broken ($0<\phi <1$), or fully broken state ($\phi =0$). The current bonding state determines the contact traction components, which is explained in the following.

Fig. 1

Frictional contact kinematics at current time $t_{n+1}$: current slave point ${\varvec{x}}_k$, current position of the previous interacting point $\displaystyle {\varvec{x}}_\ell (\hat{\varvec{\xi }}^n)$, current interacting point ${\varvec{x}}_\ell (\hat{\varvec{\xi }})$, the current elastic gap vector ${\varvec{g}}_\text {e}$ and its normal and tangential components, previous elastic gap ${\varvec{g}}^n_\text {e}$, and the sliding path ${\mathcal {C}}$. Adopted from Immel et al. (2020)

2.2.2 Adhesive friction and debonding

To account for normal adhesion and debonding in the extension of the modified Coulomb’s law, the normal traction (4) is extended by an exponential cohesive zone model (CZM) (Xu and Needleman 1992) (see Sauer (2016)), i.e.,

$$\begin{aligned} \displaystyle \begin{aligned}&{\varvec{t}}_\text {n} = \\&{\left\{ \begin{array}{ll} {\varvec{0}}, &{} g_\text {b} \ge g_\text {n} \text { or } \phi =0,\\ \frac{\phi _0 t_0 g_\text {n}}{g_0} \,\text {exp}\left( 1-\frac{g_\text {n}}{g_0}\right) {\varvec{n}}, &{} 0 \le g_\text {n}<g_\text {b} \text { and } \phi > 0, \\ -\epsilon _\text {n} \, {\varvec{g}}_\text {n}, &{} g_\text {n} < 0, \end{array}\right. } \end{aligned} \end{aligned}$$

(13)

where $t_0$ is the maximum positive normal traction, $g_0$ is the contact distance, where the maximum traction $t_0$ occurs, and $g_\text {b}$ is a cutoff distance, where contact is lost. The parameters $t_0, g_0, g_\text {b}$ depend on the interface. The normal traction model (13) is illustrated in Figure 2a.

Equation (13) implies that, when pulling the contact surfaces apart in normal direction, as long as the point remains fully or partially bonded ($\phi > 0$) the normal traction keeps increasing until $g_\text {n} = g_\text {b}$. As soon as the point is fully debonded ($\phi =0$) or the normal gap is $g_\text {n}>g_\text {b}$, the contact is lost and the normal traction component becomes ${\varvec{t}}_\text {n}={\varvec{0}}$. The sharp drop in the normal traction at $g_\text {b}$ is motivated by observations from experimental pull-out tests of osseointegrated, coin-shaped implants (Rønold and Ellingsen 2002; Nonhoff et al. 2015).

In this work, the modified Coulomb’s law introduced by Immel et al. (2020) is used to model the tangential debonding of osseointegrated interfaces. There, the friction coefficient $\mu$ is modeled as a function of the scalar state variable $\phi$, as

$$\begin{aligned} \mu := \mu (\phi ) = \phi \, \mu _\text {ub} + \left( 1 - \phi \right) \mu _\text {b}, \end{aligned}$$

(14)

where $\mu _\text {ub}$ and $\mu _\text {b}$ are the friction coefficient for the unbroken (initial) and broken state, respectively, that are weighted according to the state variable. The friction function (14) is illustrated in Figure 3. To enable sliding for tensile normal traction in the present extension of the contact model, the tangential sliding limit is shifted by

$$\begin{aligned} t^\text {slide}_\text {t} = \mu (\phi ) \left( t_0 - t_\text {n}\right) , \end{aligned}$$

(15)

according to Mergel et al. (2019, 2021).

The slope of the function $t_\text {n}(g_\text {n})$ and $t^\text {slide}_\text {t}(\mu )$ at $g_\text {n}=0$ depends solely on the choice of the parameters. It is smooth when

$$\begin{aligned} \epsilon _\text {n} = \phi _0 \, e \frac{t_0}{g_0}, \end{aligned}$$

(16)

(where e is Euler’s number), otherwise it is discontinuous. A comparison of the standard Coulomb’s law and the proposed extended modified Coulomb’s law based on adhesive friction is shown in Figure 2. A list of all constant parameters of the modified Coulomb friction law with their respective mode of determination is given in Table 1.

Fig. 2

Illustration of cohesive zone model (a) and the extended modified Coulomb’s law with adhesive friction (b) for $\phi _0 > 0$

Fig. 3

Illustration of the friction function $\mu (\phi )$ based on Eqs. (9) and (14) with respect to the accumulated deformation $g_\text {d}$. Parameter $a_\text {s}$ denotes a displacement threshold

Table 1 List of the constant parameters of the modified Coulomb friction law with their mode of determination

3 Application to coin-shaped implants

To demonstrate the new contact formulation (13) and (15) a simple implant model of an osseointegrated coin-shaped implant (CSI) (Mathieu et al. 2011; Vayron et al. 2012; Mathieu et al. 2012a, b; Vayron et al. 2014; Fraulob et al. 2020b) is used to simulate different debonding modes.

3.1 Setup

We consider the same basic setup as in Immel et al. (2020), which is briefly summarized below. The CSI is a short cylinder with radius $R_\text {i}=2.5$ mm and height $H_\text {i}=3$ mm. The bone sample is modeled as a rectangular cuboid with dimensions 12.5 $\times$ 12.5 $\times$ 5 mm. The implant is positioned at the center of the upper bone surface. 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_\text {i}=113$ GPa, $\nu _\text {i}=0.3$). The elastic properties of the bone block are $E_\text {b}=18$ GPa, $\nu _\text {b}=0.3$. Furthermore, all materials are assumed to be homogeneous and isotropic and both contact surfaces are assumed to be perfectly flat.

The bodies are meshed according to the parameters given in Table 2, where $n_\text {el}$ denotes the number of elements of the body/surface and $n_\text {gp}$ denotes the number of Gauss points per element. While the bulk is discretized with linear Lagrange shape functions, the contact surfaces are discretized with enriched contact finite elements based on quadratic non-uniform, rational B-splines (NURBS) (Corbett and Sauer 2014, 2015).

Table 2 CSI debonding: Number of finite elements $n_\text {el}$, type of shape functions and number of Gauss points per element $n_\text {gp}$ for the two bodies and their contact surfaces. Adopted from Immel et al. (2020)

The parameters of the state function (9) are chosen to be $a_\text {s}=22$ ${\upmu }$m, $b_\text {s}=0.74$, $\mu _\text {ub}=0.44$, and $\mu _\text {b}=0.3$, based on Immel et al. (2020). The initial osseointegration is constant across the bone-implant interface and is set to be $\phi _0=1$ (perfectly integrated). Due to the lack of experimental data, the cohesive zone parameters $g_0, g_\text {b}$ are set to $g_\text {b}=g_0=3a_\text {s}$, for simplicity.

The maximum traction of the cohesive zone model, $t^*_0=1.8$ MPa, is calibrated based on the results of Rønold and Ellingsen (2002) for polished, titanium CSIs. In that experimental study, CSIs with different surface roughness were implanted into rabbit tibiae and allowed to osseointegrate for 10 weeks. Then, the implants were removed together the surrounding bone. The bone and implant parts were fixed into a tensile test machine and the implant was pulled constantly in the normal direction until it was completely debonded from the bone. For the polished CSI an average degree of osseointegration of $\phi _0=0.26$ and an average maximum pull-out force of 9 N were determined, which results in an approximate 35 N for $\phi =1$. All parameters of the contact model and their values are listed in Table 3.

Table 3 List of the constant parameters of the modified Coulomb friction law and the chosen values

The boundary conditions and considered test configurations are shown in Figure 4. The lower surface of the bone block is fixed in all directions. In this work, only quasi-static conditions are considered. First, the implant is pressed into the bone block until a normal reaction force of -70 N is reached, as is done in Immel et al. (2020) for the corresponding parameter set $a_\text {s}=22$ ${\upmu }$m, $b_\text {s}=0.74$, $\mu _\text {ub}=0.44$, and $\mu _\text {b}=0.3$. Then, for the first three test cases, full and homogeneous initial osseointegration ($\phi _0=1$) of the bone-implant interface is applied (Fig. 4(a) and (c)). For test cases with tension, the implant is then pulled in normal direction until an average normal reaction force of 20 N is reached (Fig. 4(b)). Last, debonding with no initial pressure or tension is considered (Fig. 4(d)).

Fig. 4

CSI debonding: Illustrations of the boundary conditions for different debonding cases. (a) Debonding under an initial compression of -70 N, either in mode II or III. (b) Debonding under an initial tension of 20 N, either in mode II or III. Mixed mode debonding with initial compression of -70 N (ci) and without initial contact force (di), and under various loading angles $\alpha$. Mixed mode debonding with initial compression of -70 N (cii) and without initial contact force (dii), under loading angle $\alpha =45^\circ$ and various CZM parameter $t_0$ in Eq. (13). Mixed mode debonding with initial compression of -70 N (ciii) and without initial contact force (diii), under loading angle $\alpha =45^\circ$ and various initial degrees of osseointegration $\phi _0$. Altogether five test cases are investigated as discussed in the text

Then, the new contact model is examined for five different debonding test cases:

• 1. mode II: the upper implant surface is moved in x-direction under constant compression (mode IIa) or tension (mode IIb).

• 2. mode III: the upper surface of the implant is rotated around its z-axis under constant compression (mode IIIa) or tension (mode IIIb).

• 3. mode I+II: the upper implant surface simultaneously pulled along the z-axis and in x-direction, corresponding to an angle $\alpha =30, 45,$ or 60$^\circ$. This is performed with initial compression (mode I+IIci) and without initial contact force (mode I+IIdi).

• 4. mode I+II: the upper implant surface is simultaneously pulled along the z-axis and in x-direction, corresponding to an angle of $\alpha =45^\circ$ for different choices of $t_0\in [t^*_0/2, t^*_0, 2t^*_0]$. This is performed with initial compression (mode I+IIcii) and without initial contact force (mode I+IIdii).

• 5. mode I+II: the upper implant surface is simultaneously pulled along the z-axis and in x-direction, corresponding to an angle of $\alpha =45^\circ$ for increasing degrees of initial osseointegration $\phi _0\in [0, 0.25, 0.5, 0.75, 1]$. This is performed with initial compression (mode I+IIciii) and without initial contact force (mode I+IIdiii).

All simulations are performed with an in-house, MATLAB-based solver (R2019b, The MathWorks, Natick, MA, USA). Contact is computed with a penalty regularization, and the corresponding penalty parameter is set to $\epsilon _\text {n}=\epsilon _\text {t}= E_\text {b}/L_0$, with $L_0=0.01$ m. The step size for all load cases is $\Delta u=0.65$ ${\upmu }$m (for applied displacement loads and $\Delta \theta =0.1^\circ$ for applied rotations).

3.2 Results

In the following, the results of the debonding tests for the CSI, in terms of load–displacement curves, obtained with the modified Coulomb’s law (MC) and its new extension to adhesive friction (EMC) are presented and compared with each other.

3.2.1 Test 1: Mode II debonding

Figure 5 shows the normal and tangential reaction forces $F_z$ and $F_x$ for debonding and possible subsequent sliding in (tangential) x-direction under prescribed constant compression or tension. For a constant compression of −70 N, the slope of the curve of the tangential reaction force is identical for the MC and the EMC. The maximum tangential reaction force increases from 30 N to 45 N for the EMC.

Fig. 5

CSI debonding: Variation of the normal reaction $F_z$ (top) and the tangential reaction force $F_x$ (bottom) as a function of the tangential displacement $u_x$ for mode II debonding under constant compression (mode IIa) or constant tension (mode IIb). The results illustrate the difference between the modified Coulomb’s law (MC) and the new extension (EMC)

For a constant tension of 20 N, the tangential reaction force reaches up to 6 N before decreasing and dropping to 0 because of the absence of contact. The maximum tangential reaction force under tension is smaller than under compression, due to the decrease in $\phi$ stemming from the accumulated deformation in normal direction before the debonding started (due to pulling the implant back up before sliding). The contact is lost abruptly after the limit for the accumulated deformation $g_\text {d}$ is reached, due to the positive contact gap at the bone-implant interface.

3.2.2 Test 2: Mode III debonding

Figure 6 shows the normal reaction force $F_z$ and the debonding torque $M_z$ for mode III debonding due to rotation around the implant’s (normal) $z-$axis under prescribed constant compression or tension for the considered contact laws. For a constant compression of -70 N, the slope of the torque curve is identical for both contact laws. The maximal torque increases by 0.027 Nm (about 50%) when including normal adhesion (EMC). For a constant tension of 20 N the torque reaches 0.011 Nm and then decreases down to zero due to loss of contact. This loss is gradual, starting in the external region of the cylinder and propagating inward to its center. These results emphasize the fact that torque tests yield a stable crack propagation, which is particularly interesting when it comes to assessing the effective adhesion energy of the bone-implant interface (Mathieu et al. 2012a).

Fig. 6

CSI debonding: Variation of the normal reaction $F_z$ (top) and the torque $M_z$ (bottom) as a function of the rotation angle $\theta$ for mode III debonding under constant compression (IIIa) or constant tension (IIIb). The results illustrate the difference between the modified Coulomb’s law (MC) and its new extension (EMC)

3.2.3 Test 3: Mode I+II debonding for varying angles

Figure 7 shows the normal and tangential reaction forces $F_z$ and $F_x$ for mixed mode debonding under different angles $\alpha$ (mode I+IIci) starting from an initial contact pressure, based on the MC and the EMC. The normal reaction force $F_z$ increases linearly until it reaches zero. For each angle, the slope of the reaction force curve is identical for both considered contact laws, respectively. In case of the EMC, the reaction force becomes positive at some point and follows the debonding curve of cohesive zone model (13) seen in Fig. 2a. In all presented cases, the debonding occurs because the maximal normal gap exceeds $g_\text {b}$, due to the prescribed upward movement. Therefore, increasing the debonding angle $\alpha$ decreases the amount of tangential deformation necessary for debonding, i.e., where contact is lost and the reaction force becomes zero.

The tangential reaction force $F_x$ increases linearly until the respective sliding limit is reached. Then the implant starts sliding and the tangential reaction force decreases linearly, as long as the normal force $F_z$ is still negative. When the normal reaction force reaches zero, the cases with the MC show zero tangential reaction force, as there is no contact anymore. For the cases with the EMC, the bone-implant interface has not fully debonded yet and thus, there is still a normal (adhesive) contact force building up. As a result, the tangential reaction force decreases nonlinearly until it reaches zero.

Fig. 7

CSI debonding: Variation of the normal reaction force $F_z$ (top) and the tangential reaction force $F_x$ (bottom) as a function of the tangential displacement $u_x$ for mixed mode debonding, starting from an initial contact pressure (mode I+IIci). The results illustrate the difference between the modified Coulomb’s law (MC) and its new extension (EMC)

Fig. 8

CSI debonding: Variation of the normal reaction force $F_z$ (top) and the tangential reaction force $F_x$ (bottom) as a function of the tangential displacement $u_x$ for mixed mode debonding, starting from zero contact force (mode I+IIdi). The results show the extended modified Coulomb’s law (EMC)

While the maximum normal reaction force is the same for all three tested angles, the maximum tangential reaction force decreases for an increasing debonding angle. The respective maximal tangential reaction force for each case is around 1.5 times higher for the EMC compared to the MC, due to the shift in tangential contact traction (15).

Figure 8 shows the normal and tangential reaction forces $F_z$ and $F_x$ for mixed mode debonding under different angles $\alpha$ (mode I+IIdi), with no initial contact force, based on the EMC. The different curves of the normal reaction force $F_z$ are identical to the positive part of the curves in Figure 7. As there is no initial compression or tension in the beginning of this loading case, all curves begin at the origin. Similarly, the different curves for the tangential reaction force $F_x$ are identical to the exponential part of the curves in Figure 7 and are shifted by the same displacement toward the origin, respectively, as the normal reaction force.

3.2.4 Test 4: Mode I+II debonding for varying CZM parameter $t_0$

Figure 9 shows the normal and tangential reaction force $F_z$ and $F_x$ for mixed mode debonding under $\alpha =45^\circ$ as a function of the tangential displacement with different values of the maximal CZM traction $t_0$ when considering adhesive friction (mode I+IIcii). As expected, the maximum normal and tangential reaction forces are proportional to $t_0$ (see Eqs. (13) and (15)). Furthermore, the slopes of $F_z$ and $F_x$ at the transition from compression to tension ($F_x=0$, $u_x=6.5$ mm) are smooth for about $t_0=2t^*_0$ (as $2t^*_0 \approx \epsilon _\text {n}g_0/\phi _0 e$ (see Eq. 16)).

Figure 10 shows the normal and tangential reaction forces $F_z$ and $F_x$ for mixed mode debonding under $\alpha =45^\circ$, with no initial contact force, as a function of the tangential displacement for different values of the maximal CZM traction $t_0$ when considering adhesive friction (mode I+IIdii). The different curves of the normal reaction force $F_z$ are identical to the positive part of the curves in Figure 9. As there is no initial contact force in the beginning of this loading case, all curves begin at the origin. Similarly, the different curves of the tangential reaction force $F_x$ are identical to the exponential part of the curves in Figure 9 and are shifted by the same displacement toward the origin, respectively, as the normal reaction force.

Fig. 9

CSI debonding: Variation of the normal reaction force $F_z$ (top) and the tangential reaction force $F_x$ (bottom) as a function of the tangential displacement for mixed mode debonding under $\alpha =45^\circ$ and varying CZM parameter $t_0$, starting from an initial contact pressure (mode I+IIcii). The results show the extended modified Coulomb’s law (EMC)

Fig. 10

CSI debonding: Variation of the normal reaction force $F_z$ (top) and the tangential reaction force $F_x$ (bottom) as a function of the tangential displacement for mixed mode debonding under $\alpha =45^\circ$ and varying CZM parameter $t_0$, starting from zero contact force (mode I+IIdii). The results show the extended modified Coulomb’s law (EMC)

3.2.5 Test 5: Mode I+II debonding for varying degree of osseointegration

Figure 11 shows the normal and tangential reaction force $F_z$ and $F_x$ for mixed mode debonding under 45$^\circ$ as a function of the tangential displacement with increasing degree of osseointegration $\phi _0$ when considering adhesive friction (mode I+IIciii). Increasing the degree of initial osseointegration increases the peak magnitude in the normal and tangential reaction force. This is due to the fact that in this test case, debonding occurs first due to $g_\text {n} > g_\text {b}$ and not due to exceeding the limit of the deformation of the interface $g_\text {d} = a_\text {s}(1+2b_\text {s})$ (see Eq. (9)). The maximal normal and tangential reaction force increases proportionally with increasing $\phi _0$, while the tangential displacement necessary for debonding remains the same.

Fig. 11

CSI debonding: Variation of the normal reaction force $F_z$ (top) and the tangential reaction force $F_x$ (bottom) as a function of the tangential displacement for mixed mode debonding under $\alpha =45^\circ$ and initial degree of osseointegration $\phi _0$, starting from an initial contact pressure (mode I+IIciii). The results show the extended modified Coulomb’s law (EMC)

Fig. 12

CSI debonding: Variation of the normal reaction force $F_z$ (top) and the tangential reaction force $F_x$ (bottom) as a function of the tangential displacement for mixed mode debonding under $\alpha =45^\circ$ and initial degree of osseointegration $\phi _0$, starting from zero contact force (mode I+IIdiii). The results show the extended modified Coulomb’s law (EMC)

Figure 12 shows the normal and tangential reaction forces $F_z$ and $F_x$ for mixed mode debonding under $\alpha =45^\circ$, with no initial contact force, as a function of the tangential displacement with increasing degree of osseointegration $\phi _0$ when considering adhesive friction (mode I+IIdiii). The different curves of the normal reaction force $F_z$ are identical to the positive part of the curves in Figure 11. As there is no initial contact force in the beginning of this loading case, all curves begin at the origin. Similarly, the different curves of the tangential reaction force $F_x$ are identical to the exponential part of the curves in Figure 11 and are shifted by the same displacement toward the origin, respectively, as the normal reaction force.

4 Application to acetabular cup implants

Both, the MC and the EMC have been examined for a simple implant model in Section 3. Now, both models are applied to simulate the debonding of a 3D, osseointegrated acetabular cup implant (ACI), under different removal conditions, similar to the simulations in Raffa et al. (2019) and Immel et al. (2021a). However, the aforementioned works only considered primary stability. Here, the implant’s secondary stability is considered and is quantified by assessing the removal force/torque. Furthermore, the biomechanical determinants of the long-term stability, such as primary stability and initial degree of osseointegration, are assessed. The results of the modified Coulomb’s law (MC) and its extension to adhesive friction (EMC) are compared to assess the importance of adhesive effects for long-term stability because it allows to distinguish the influence of primary stability and osseointegration phenomena on the secondary stability.

4.1 Setup

A simple cylindrical block is considered, as it is a suitable simplification of the pelvis geometry that qualitatively captures the relevant contact conditions. The same geometry of the ACI including the ancillary used in Raffa et al. (2019) and Immel et al. (2021a) is considered herein and is briefly summarized in the following. An idealized bone block with the same properties as in Raffa et al. (2019) is used in order to calibrate the model and compare results. The bone block is modeled as a cylinder with a radius of 50 mm and a height of 40 mm. A hemi-spherical cavity is cut into the cylinder with a radius $R_\text {b}$ based on the fixed radius of the implant $R_\text {i}$ and the chosen interference fit $I\!F$, i.e., $R_\text {b}=R_\text {i}-I\!F/2$. The edge of the cavity has a fillet radius of 2 mm.

As in Section 3, the bodies are meshed with surface-enriched hexahedral elements according to the parameters given in Table 4. The finite element mesh is shown in Figure 13. A refinement analysis of the mesh and the load step size is performed to ascertain mesh convergence (see Appendix A) for the reference case (see Section 4.1.1).

Fig. 13

(a) Finite element mesh of the ACI, bone block, and ancillary (Immel et al. 2021b). (b) Bottom view of the ACI. (c) Top view of the bone block. A very fine mesh is used around the rim of the cavity, as the contact forces are expected to vary most strongly there

Table 4 ACI debonding: Number of finite elements $n_\text {el}$, type of shape functions and number of Gauss points per element $n_\text {gp}$ for the two bodies and their contact surfaces

4.1.1 Model parameters

The ancillary and the ACI are assumed to be made of stainless steel ($E_\text {a}=211$ GPa, $\nu _\text {a}=0.3$) and titanium alloy (Ti-Al6-V4; $E_\text {i}=113$ GPa, $\nu _\text {i}=0.3$), respectively. The bone block is assumed to consist only of trabecular bone tissue ($\nu _\text {b}=0.3$ (Yew et al. 2006)), without an outer cortical layer. This is in accordance with a previous study (Immel et al. 2021a) and findings in the literature (Anderson et al. 2005; Phillips et al. 2007; Watson et al. 2017) that indicate that the reaming of the hip performed during surgery may completely remove cortical bone tissue from the contact area. For all bodies, the Neo-Hookean material model of Eq. (1) is used. Furthermore, all materials are assumed to be homogeneous and isotropic and both contact surfaces are assumed to be perfectly smooth.

Table 5 List of the constant parameters of the modified Coulomb friction law and the chosen values, as well as the range for the varied parameters

In this work, the effect of various biomechanical properties of the bone-implant system on the ACI long-term stability is assessed. Therefore, different degrees of osseointegration from $0-100$% with homogeneous distribution over the contact surfaces are considered. Furthermore, the influence of varying bone stiffness $E_\text {b}=0.1-0.6$ GPa (Phillips et al. (2007); Janssen et al. (2010)), interference fit $I\!F=0-2.0$ mm ( Kwong et al. (1994); Macdonald et al. (1999); Spears et al. (1999)), and friction coefficient $\mu _\text {b}=0-0.7$ (Dammak et al. (1997); Spears et al. (1999); Novitskaya et al. (2011)) on the long-term stability are analyzed. The corresponding friction coefficient $\mu _\text {ub}=0.15-1$ is taken from Table 4 from Immel et al. (2020) and is roughly 1.5 times higher than $\mu _\text {b}$. Based on previous studies (Raffa et al. 2019; Immel et al. 2021a) the parameter set of $E^*_\text {b}= 0.2$ GPa (Phillips et al. (2007)), $I\!F^*=1$ mm (Kwong et al. (1994)), and $\mu ^*_\text {b}=0.3$ (Dammak et al. (1997)) is denoted as the reference case and marked with $*$. The parameters of the state function (9) are chosen to be $a_\text {s}=128$ ${\upmu }$m and $b_\text {s}=1.84$, which does not affect the maximum of the removal force/torque. Due to the lack of experimental data, the values of $a_\text {s}$ and $b_\text {s}$ are chosen large enough so that the debonding process is visible and a removal force/torque can be identified (Fig. 16). The coefficients of the cohesive zone model $t_0=t^*_0=1.8$ MPa and $g_\text {b}=g_0=3a_\text {s}$ are calibrated based on the results of Rønold and Ellingsen (2002) for polished CSI, as was done in Section 3.1. All parameters of the contact model and their values, as well as the studied parameters are listed in Table 5.

4.1.2 Boundary and loading conditions and solver settings

The bone block is fixed in all directions at the bottom surface. As before, only quasi-static conditions are considered. The simulations of implant insertion and subsequent removal are comprised of three stages (see Fig. 14):

• 1. insertion: the implant is inserted vertically into the cavity, by pushing the upper surface of the ancillary in negative z-direction, until the reaction force at the top of the ancillary reaches $F_0=-2500$ N, similar to values found in the literature (Sotto-Maior et al. 2010; Souffrant et al. 2012; Le Cann et al. 2014) and to what was done in previous studies (Raffa et al. 2019; Immel et al. 2021a). The downward displacement attained at the top of the ancillary for $F_0=-2500$ N is denoted $d_0$. It depends on the considered parameters $E_\text {b}, I\!F$, and $\mu _\text {b}$ and thus changes for each case, i.e., $d_0=d_0$($\mu _\text {b}$, $I\!F$, $E_\text {b}$, $F_0$).

• 2. osseointegration: the contact surfaces are assumed to be homogeneously osseointegrated with an initial degree of osseointegration varying from $\phi _0\in [0, 0.25, 0.5, 0.75, 1].$

• 3. removal: the implant is removed either

  • global¹ mode I: by displacing the upper surface of the ancillary in positive z-direction by -$d_0$,

  • global mode II: by displacing the center of the upper surface of the ancillary in positive x-direction by $d_0$,

  • mode III: by rotating the upper surface of the ancillary around its z-axis by $\theta =10^\circ$.

The three simulation stages are shown in Figure 14(a) and the different removal cases are illustrated in Figure 14(b). The example of mode I debonding is shown, with the final output of the load–displacement curve inside the red square (Fig. 14(a)).

The stability of the configuration is then assessed by determining the maximum pull-out force in normal direction, $F^\text {max}_z$, the maximum pull-out force in tangential direction, $F^\text {max}_x$, or the maximum debonding torque $M^\text {max}_\text {z}$.

Contact is computed with a penalty regularization, and the corresponding penalty parameter is chosen based on the Young’s modulus of trabecular bone as $\epsilon _\text {n} = \epsilon _\text {t} = E_\text {b}/L_0$, with $L_0=R_\text {I}=0.0255$ m corresponding to the radius of the implant. The number of load steps for the different simulation stages are: $l_1=l_{3\text {.modeI}}=l_{3.\text {modeIII}}=100$ and $l_{3.\text {modeII}}=1000.$ All simulations were performed with an in-house, MATLAB-based solver (R2019b, The MathWorks, Natick, MA, USA) with MATLAB’s own parallelization. Computations were performed on the RWTH Compute Cluster (Intel HNS2600BPB, Platinum 8160) with 20 cores. The average computing time for the different contact laws and loading cases is listed in Table 8. The computing time is sensitive to the parameter combination. Parameter combinations that produce high pull-out forces/debonding torque have a longer computing time. The difference in computing time between the debonding tests and the contact models is discussed in Section 5.3.

Fig. 14

ACI debonding: (a) Illustration of the three simulation stages and (b) the three removal tests

4.2 Debonding without adhesion in normal direction

First, the results of the removal tests, in terms of load–displacement curves and pull-out force/ debonding torque, obtained with the MC are presented. The results with the EMC follow in Section 4.3, and a comparison is given in Section 5.1.

4.2.1 Global mode I: Normal pull-out test

Figure 15 (a) shows the normal reaction force $F^*_z$ for the reference case, which increases and reaches a peak at a displacement of 0.25–0.32 mm and then slowly decreases to zero. This maximum coincides with the start of the decrease of the average degree of osseointegration of the bone-implant interface ${\bar{\phi }}$ (Fig. 15, (b)). At a displacement of 1.07 mm, the reaction force becomes independent from $\phi _0$. At this point, the bone-implant interface is completely debonded (${\bar{\phi }}=0$) and only pure Coulomb’s friction is taking place until the contact at the bone-implant interface is lost completely at a displacement of about 4.25 mm. In this test, osseointegration only affects the magnitude of the peak, while the overall slope of the load–displacement curve remains unchanged when increasing the initial degree of osseointegration. The location of the peak does not change significantly with increasing $\phi _0$.

Due to the lack of experimental data for this work, the values of $a^*_\text {s}=128$ ${\upmu }$m and $b^*_\text {s}=1.84$ are chosen large enough so that the debonding process is visible and a removal force/torque can be identified. The effect of changing the value of $a_\text {s}$ and $b_\text {s}$ on $F^*_z(\phi _0=1)$ is shown in Figure 16. Naturally, both parameters have no effect on the mechanical behavior before debonding and on the maximum pull-out force. Decreasing or increasing $a_\text {s}$ and $b_\text {s}$ decreases or increases the amount of deformation that is necessary for the interface to fully debond (about 0.7, 0.75, 1.1, 1.75, 1.9 mm, respectively, see Fig. 16). After debonding, only pure Coulomb’s friction takes place until the contact between bone and implant is lost completely (after a displacement of about 4.25 mm).

Figure 17(a)–(c) shows the maximum normal pull-out force $F^\text {max}_z$ as a function of the interference fit $I\!F$, trabecular bone stiffness $E_\text {b}$, friction coefficient $\mu _\text {b}$, for different values of the initial degree of osseointegration $\phi _0$. The results obtained for $F^\text {max}_z$ with $\phi _0=0$ are identical to the results from Raffa et al. (2019). First, the pull-out force increases as a function of $I\!F, E_\text {b}, \mu _\text {b}$, then reaches a peak, and eventually decreases. The maximum value of the pull-out force is obtained for around $I\!F=1.4$ mm, $E_\text {b}=0.4$ GPa, and $\mu _\text {b}=0.6$. For $\mu _\text {b} \le 0.15$ the pull-out force is zero, for all degrees of osseointegration.

Fig. 15

Normal debonding without adhesion for the reference case: (a) Variation of the normal force $F^*_z$ a function of the initial degree of osseointegration $\phi _0$. The maximum pull-out force $F^{\text {max}*}_z$ is marked with $*$. (b) Average degree of osseointegration of the bone-implant interface ${\bar{\phi }}$ for an initial degree of osseointegration $\phi _0=1$

Fig. 16

Normal debonding without adhesion: Variation of the normal force $F^*_z$ for $\phi _0=1$ as a function of the MC parameters $a_\text {s}$ and $b_\text {s}$,

Fig. 17

Normal debonding without adhesion: Variation of the maximum normal pull-out force $F^\text {max}_z$ as a function of the initial degree of osseointegration $\phi _0$ and (b) the interference fit $I\!F$, (c) the trabecular Young’s modulus $E_\text {b}$, and (d) the friction coefficient $\mu _\text {b}$. The reference case is marked with $*$

4.2.2 Global mode II: Tangential pull-out test

The tangential reaction force $F^*_x$ for the reference case increases and reaches a peak at a displacement of about 75 ${\upmu }$m and then slowly decreases to zero (Fig. 18(a)). The average degree of osseointegration of the bone-implant interface ${\bar{\phi }}$ starts to decrease already beyond 34 ${\upmu }$m (Fig. 18 (b)).

Fig. 18

Tangential debonding without adhesion for the reference case: (a) Variation of the tangential force $F^*_x$ as s function of the initial degree of osseointegration $\phi_0$. The maximum pull-out force $F^{\text {max}*}_x$ is marked with $*$. (b) Average degree of osseointegration of the bone-implant interface ${\bar{\phi }}$ for an initial degree of osseointegration $\phi _0=1$

At a displacement of about 0.3 mm, the reaction force becomes independent from $\phi _0$. Similarly to the normal pull-out test, increased osseointegration only affects the magnitude of the tangential pull-out force, while the location of the peaks and the initial slope of the curves for different degrees of osseointegration are very similar.

Figure 19(a)–(c) shows the maximum tangential pull-out force $F^\text {max}_x$ as a function of the interference fit $I\!F$, trabecular bone stiffness $E_\text {b}$, friction coefficient $\mu _\text {b}$, for different values of the initial degree of osseointegration $\phi _0$. First, the pull-out force increases as a function of $I\!F, E_\text {b}$ and $\mu _\text {b}$, then reaches a peak, and eventually decreases. The maximum value of the pull-out force is obtained for around $I\!F=1.4$ mm, $E_\text {b}=0.4$ GPa, and $\mu _\text {b}=0.6$ – the same values as for the normal pull-out test. For $\mu _\text {b} \le 0.15$ the pull-out force is zero for all degrees of osseointegration. Tangential pull-out forces are roughly one magnitude lower than the corresponding normal pull-out force, which agrees with observations from clinical practice. During surgery, after the insertion of the ACI, surgeons often attempt to lever out an acetabular cup to test the seating of the ACI. That is, the surgeon applies a tangential force, such as is considered here, instead of a normal force since normal pull-out would require too much force.

Fig. 19

Tangential debonding without adhesion: Variation of the maximum tangential pull-out force $F^\text {max}_x$ as a function of the initial degree of osseointegration $\phi _0$ and (a) the interference fit $I\!F$, (b) the trabecular Young’s modulus $E_\text {b}$, and (c) the friction coefficient $\mu _\text {b}$. The reference case is marked with $*$

4.2.3 Mode III: Torsional debonding test

Figure 20 (a) shows the debonding torque $M^*_z$ as a function of the rotation angle for different values of $\phi _0$ and the reference case. The torque increases, reaches a peak at an angle of about 3$^\circ$, and then decreases to reach a constant torque of 47 Nm at about 4.5$^\circ$ due to the present compressive normal force. The degree of osseointegration starts to decrease at an angle of about 2.6$^\circ$ and becomes zero at about 4.5$^\circ$ (Fig. 20 (b)). As for the normal and the tangential pull-out cases, only the magnitude of the peak of the load–displacement curve is affected when increasing the degree of osseointegration $\phi _0$.

Figure 21(a)–(c) shows the variation of the maximum debonding torque $M^\text {max}_z$ as a function of the parameters $I\!F, E_\text {b}, \mu _\text {b}$, and $\phi _0$. First, the torque increases with increasing parameter $I\!F, E_\text {b}, \mu _\text {b}$, reaches a peak, and then decreases. The maximum values of the torque are obtained around $I\!F=1.4$ mm, $E_\text {b}=0.4$ GPa, and $\mu _\text {b}=0.6$, which correspond to the same parameters as for the pull-out tests. For the interference fit $I\!F$, a larger plateau for $I\!F=1.0-1.5$ mm as compared to the pull-out tests is obtained. The maximum torque obtained for $\phi _0=1$ is 55 Nm, 68 Nm, and 65 Nm, respectively.

Fig. 20

Torsional debonding without adhesion for the reference case: (a) Variation of the debonding torque $M^*_z$ a function of the initial degree of osseointegration $\phi _0$. The maximum debonding torque $M^{\text {max}*}_z$ is marked with $*$. (b) Average degree of osseointegration of the bone-implant interface ${\bar{\phi }}$ for an initial degree of osseointegration $\phi _0=1$

Fig. 21

Torsional debonding without adhesion: Variation of the maximum debonding torque $M^\text {max}_z$ as a function of the initial degree of osseointegration $\phi _0$ and (a) the interference fit $I\!F$, (b) the trabecular Young’s modulus $E_\text {b}$, and (c) the friction coefficient $\mu _\text {b}$. The reference case is marked with $*$

4.3 Debonding with adhesion in normal direction and adhesive friction

The results corresponding to the load–displacement curves and pull-out force/ debonding torque obtained with the three removal tests and with the EMC are presented below. In addition to the modified Coulomb’s friction law (14), the EMC includes a CZM in normal direction (13) and adhesive friction (15).

4.3.1 Global mode I: Normal pull-out test

Figure 22(a) shows the variation of the normal reaction force $F^*_z$ as a function of the tangential displacement and the initial degree of osseointegration $\phi _0$ for the reference case. The normal reaction force increases, reaches a peak at a displacement of about 0.25 mm and then decreases. The effect of osseointegration and adhesive friction on the load–displacement curve is more pronounced than for the MC. This can be seen as the magnitude increase of the pull-out forces is higher and the peaks are wider (cf. Figs. 15(a) and 22). In contrast to the MC, here, $F^*_z$ depends on $\phi _0$ throughout the whole debonding process, which is due to the adhesion in normal direction.

Fig. 22

Normal debonding with adhesive friction: Variation of the normal force $F^*_z$ as a function of the initial degree of osseointegration $\phi _0$ for the reference case. The maximum pull-out force $F^{\text {max}*}_z$ is marked with $*$

Fig. 23

Normal debonding with adhesive friction: Variation of the maximum normal pull-out force $F^\text {max}_z$ as a function of the initial degree of osseointegration $\phi _0$ and (a) the interference fit $I\!F$, (b) the trabecular Young’s modulus $E_\text {b}$, and (c) the friction coefficient $\mu _\text {b}$. The reference case is marked with $*$

However, the initial slope of the normal reaction force curves does not change significantly when increasing $\phi _0$. Compared to the results obtained when considering only tangential debonding (see Fig. 15 (a)), some small oscillations after the peak are visible.

Figure 23(a)–(c) shows the maximum normal pull-out force $F^\text {max}_z$ as a function of the parameters $I\!F, E_\text {b}, \mu _\text {b}$, and $\phi _0$. The slopes of the different curves of pull-out forces are similar to the ones obtained with the MC (cf. Section 4.2.1), with the peak values obtained for the same values of $I\!F, E_\text {b},$ and $\mu$. For $\mu _\text {b} \le 0.15$ the pull-out force remains equal to zero, regardless of the degree of osseointegration.

4.3.2 Global mode II: Tangential pull-out test

Figure 24(a) shows the tangential reaction force $F^*_x$ as a function of the tangential displacement for different values of $\phi _0$. The effect of osseointegration and adhesive friction on the load–displacement curve is more pronounced with the EMC than with the MC. As for the normal pull-out test, $F^*_x$ remains dependent on $\phi _0$ throughout the whole debonding process.

Fig. 24

Tangential debonding with adhesive friction: Variation of the tangential force $F^*_x$ a function of the initial degree of osseointegration $\phi _0$ for the reference case. The maximum pull-out force $F^{\text {max}*}_x$ is marked with $*$

Fig. 25

Tangential debonding with adhesive friction: Variation of the maximum tangential pull-out force $F^\text {max}_x$ as a function of the initial degree of osseointegration $\phi _0$ and (a) the interference fit $I\!F$, (b) the trabecular Young’s modulus $E_\text {b}$, and (c) the friction coefficient $\mu _\text {b}$. The reference case is marked with $*$

The peak in tangential reaction force is reached at a displacement of about 0.08 mm. Furthermore, the increase in magnitude is considerably larger than for the MC, while remaining roughly one magnitude lower than the results for the normal pull-out case with adhesive friction. Here, larger oscillations are visible, which are discussed in Section 5.3.

Figure 25(a)–(c) shows the variation of the maximum tangential pull-out force $F^\text {max}_x$ as a function of the parameters $I\!F, E_\text {b}, \mu _\text {b}$, and $\phi _0$. While the peaks in tangential pull-out force are obtained for the same parameters as before, the slope of $F^\text {max}_x$ as a function of all parameters ($I\!F, E_\text {b}, \mu _\text {b}$) depends on the initial degree of osseointegration. As before, for $\mu _\text {b} \le 0.15$ the tangential pull-out force remains zero independent of the degree of osseointegration.

4.3.3 Mode III: Torsional debonding test

Figure 26(a) shows the variation of the debonding torque $M^*_z$ as a function of the rotation angle and the initial degree of osseointegration $\phi _0$ for the reference case. In contrast to the results obtained with the MC, the peak in torque is obtained at a rotation angle of approximately $\theta =3.5^\circ$. Then the torque decreases to a constant value due to the present compressive normal force.

Fig. 26

Torsional debonding with adhesive friction: Variation of the debonding torque $M^*_z$ a function of the initial degree of osseointegration $\phi _0$ for the reference case. The maximal torque $M^{\text {max}*}_z$ marked with $*$

Fig. 27

Torsional debonding with adhesive friction: Variation of the maximal debonding torque $M^\text {max}_z$ as a function of the initial degree of osseointegration $\phi _0$ and (a) the interference fit $I\!F$, (b) the trabecular Young’s modulus $E_\text {b}$, and (c) the friction coefficient $\mu _\text {b}$. The reference case is marked with $*$

When considering adhesive friction, the torque after full debonding does not reach the same constant values for each $\phi _0$ due the shift in the tangential sliding threshold (15), that depends on $\phi _0$.

Figure 27(a)–(c) shows the variation of the maximum debonding torque $M^\text {max}_z$ as a function of the parameters $I\!F, E_\text {b}, \mu _\text {b}$, for different values of $\phi _0$. In contrast to the pull-out test, the removal torque curves are very similar to the corresponding results obtained with the modified Coulomb’s law. Peaks in torque are obtained for the same values of $I\!F, E_\text {b}, \mu _\text {b}$ as for the pull-out tests and the modified Coulomb’s law.

5 Discussion

This work studies the contact and debonding behavior between implant and bone using a new adhesive friction model that accounts for partial osseointegration. The new extension to adhesive friction is first demonstrated on a simple model of an osseointegrated implant, following previous studies (Rønold and Ellingsen 2002; Rønold et al. 2003; Fraulob et al. 2020a, b, c; Immel et al. 2020). Then, both the original and the extended debonding models are applied to the debonding of a partially osseointegrated acetabular cup implant, which corresponds to a situation of clinical interest. The effect of increasing the osseointegration level on implant stability is examined by analyzing the behavior of the maximum removal force/torque, for three patient- and implant-dependent parameters: $I\!F, E_\text {b},$ and $\mu _\text {b}$. Overall, both debonding models provide reasonable qualitative estimates of long-term stability with higher estimates of implant stability for the extension to adhesive friction.

5.1 Comparison of the modified and extended Coulomb’s law with respect to their biomechanical relevance

When keeping parameters $E_\text {b}, I\!F$ and $\mu _\text {b}$ fixed, increasing only $\phi _0$ results in an almost linear increase of the removal force/torque. Increasing $\phi _0$ has the largest effect on the tangential pull-out force and the least on the removal torque (see Fig. 28). In general, $\mu _\text {b}$ has the highest influence on the removal force/torque and $I\!F$ the least (see Figs. 17, 19, 21, 23, 25, 27).

Figure 28 shows the ratio between the maximum removal forces/torque obtained for perfect initial osseointegration ($\phi _0=1$) and no initial osseointegration ($\phi _0=0$) ($F^\text {max}(\phi _0=1)/F^\text {max}(\phi _0=0)$) for the studied parameters and removal tests when considering both proposed models. The relative variation of the pull-out force/debonding torque obtained by considering the modified Coulomb’s law is qualitatively similar when varying $I\!F$ and $E_\text {b}$, with values ranging between 38 and 62%, with a slightly higher increase of the reaction force for lower values of $I\!F$ and $E_\text {b}$. Concerning the friction coefficient $\mu _\text {b}$, the modified Coulomb’s law shows the largest effect on the pull-out force/torque for a value of $\mu _\text {b}=0.2$. This effect then decreases when increasing the friction coefficient. The increase of the maximum pull-out force is much higher for the EMC compared to the MC with values ranging between 46 and 172%. In addition, osseointegration modeled with the EMC leads to a larger increase of the maximum removal force/torque for low parameter values, which corresponds to low initial stability but high contact area.

While the relative variation of $F^\text {max}_x, F^\text {max}_z, M^\text {max}_z$ produced by the two debonding models due to changes of $\mu _\text {b}$ are very similar, the slopes of the curves in Figure 28 for $I\!F$ and $E_\text {b}$ show considerable differences between the two contact models. The MC only has a small effect on the maximum torque for all observed parameters with a total increase in torque of 7–15%. The present extension produces a higher increase in the maximum torque of 21–35%, due to the shift in the tangential traction Eq. (15). Overall, the effect on the maximum torque remains considerably lower compared to the pull-out tests, as no contact is lost during the torsion test. Table 3 lists the average percentage increase in the maximum pull-out forces $F^\text {max}_z, F^\text {max}_x$ and debonding torque $M^\text {max}_z$ from 0 to 100% osseointegration for interference fit $I\!F$, Young’s modulus $E_\text {b}$, and friction coefficient $\mu _\text {b}$ for both contact laws.

Both presented contact models produce reasonable estimates for the long-term stability of the ACI, when compared to existing numerical results for the initial stability (Raffa et al. 2019) (see Figures 15 and 22, $\phi _0=0$). Overall, the maximum pull-out forces $F^\text {max}_x, F^\text {max}_z$ and the debonding torque $M^\text {max}_z$ all increase nearly linearly with increasing degree of osseointegration $\phi _0$ for every chosen parameter $I\!F, \mu _\text {b}, E_\text {b}$. In this work, osseointegration is shown to significantly increase implant stability (see Fig. 28 and Table 6). However, the dependence of the maximum pull-out force/debonding torque on the different parameter sets remains essentially the same as for primary stability.

Fig. 28

Ratio between the maximum removal forces/torque obtained for perfect initial osseointegration ($\phi _0=1$) and no initial osseointegration ($\phi _0=0$) as a function of (a) the interference fit $I\!F$, (b) the trabecular Young’s modulus $E_\text {b}$, and (c) the friction coefficient $\mu _\text {b}$ for the different removal tests. Shown are results for the modified Coulomb’s law (MC) and the new extension (EMC). The reference case is marked with $*$. Some results for $I\!F=0$ mm, $E_\text {b}=0$, and $\mu _\text {b}=0-0.15$ are omitted, as there is no measurable increase in removal force/debonding torque

Table 6 Average percentage increase in the maximum pull-out forces $F^\text {max}_z, F^\text {max}_x$ and debonding torque $M^\text {max}_z$ from 0 to 100% osseointegration for interference fit $I\!F$, Young’s modulus $E_\text {b}$, and friction coefficient $\mu _\text {b}$ for the modified Coulomb’s law (MC) and its new extension (EMC)

The two presented contact models indicate that poor initial stability will lead to poor or suboptimal long-term stability, which emphasizes the crucial role of primary stability for the implant outcome. This finding is in agreement with the literature, where initial stability is determined as the governing factor of long-term stability (Pilliar et al. 1986; Engh et al. 1992, 2004; Rittel et al. 2018), as the mechanical conditions at the bone-implant interface have a significant effect on bone growth and remodeling. Furthermore, the present extension has a higher effect on poor initial stability, stressing the importance of adhesion for low initial stability.

Both presented debonding models also allow the assessment of how loading that does not result in complete debonding affects the remaining osseointegration state $\phi$ of the bone-implant interface (e.g. Fig. 15 (b)). Future studies that couple the EMC with cyclic loading and bone growth and remodeling could, for example, provide answers on how daily loading affects the bonding state of the interface during and after healing.

5.2 Comparison with similar studies

Since most numerical studies that model osseointegrated interfaces assume perfectly bonded surfaces and thus, do not simulate the actual debonding of the interface, only few comparisons with existing work can be made. One comparable work is the study of Rittel et al. (2018), where the influence of partial osseointegration on dental implant stability and cohesive failure was studied. There, a tie constraint was applied to parts of the bone-implant interface throughout the simulation, such that bone-implant debonding occurred as cohesive failure in the bone around the bone-implant interface. Partial osseointegration was modeled by defining a relative osseointegrated area with a random distribution and restricting non-integrated areas to frictional contact. One key finding of the study of Rittel et al. (2018) was that none of their removal tests was able to distinguish osseointegration above 20% and that the torque test was more accurate than a pull-out test in determining the degree of osseointegration. Based on these findings, it was concluded that osseointegration of only 20% of the bone-implant interface provides sufficient long-term stability. In the present study, opposite findings are obtained. Here, all considered debonding tests show consistent increase in stability for increasing initial degree of osseointegration. Furthermore, osseointegration showed the least effect on the debonding torque and the highest for mode II debonding. The difference between the two studies might stem from the difference between the cohesive failure model of Rittel et al. (2018) and the adhesive failure models presented here, and/or the difference in geometry and contact conditions. Further studies and especially experimental testing, as proposed in Section 3, are necessary in order to calibrate and validate the proposed contact models.

5.3 Numerical stability

Mesh convergence was investigated for the reference case and the modified Coulomb’s friction law (see Appendix A). The load–displacement curves obtained when considering adhesive friction (Fig. 22 and 24(a)) show oscillations in the reaction force after the peak and require an increased number of Newton–Raphson iterations and thus, increased computing (see Appendix A). In the cases of normal and tangential debonding, the added adhesion in normal direction results in alternating sticking and sliding (so called stick-slip motion), producing oscillations in the results. The quasi-static assumption used in this work is not suitable in those cases and a dynamic simulation should be performed instead to account for the inertia in the system.

Due to the lack of experimental data and comparable numerical results, the a priori assessment of the choice of mesh, boundary conditions and relevance of inertia, remains difficult and thus the results can only provide a qualitative statement of the relevance of the analyzed parameters on implant stability.

5.4 Perspectives and guidelines for future work

In the following, perspectives for future extensions and applications of the proposed bone and contact models are discussed. Furthermore, we state guidelines for future experimental testing, in order to obtain relevant data to calibrate and validate the proposed models.

5.4.1 Bone modeling perspectives

This work uses idealized bone geometries. This was done in order to use results from Raffa et al. (2019) as calibration for cases with $\phi _0=0$. Further, our work focuses on the contact behavior of the osseointegrated bone-implant interface. The contact geometry and contact conditions of the hemispherical cavity are very similar to a generic pelvis. While the simplified bone geometry is a justified simplification in this work, an analysis of, for example, different pelvis shapes and defects on the contact behavior of the bone-implant interface would be of clinical relevance.

The bone block was modeled with trabecular bone without a cortical layer and the bone was rigidly fixed at the entire bottom surface. The absence of cortical bone in the contact area is in accordance with a previous study (Immel et al. 2021a) and findings in the literature (Anderson et al. 2005; Phillips et al. 2007; Watson et al. 2017) that indicate that the reaming performed during surgery may completely remove cortical bone tissue from the contact area. Future studies should include a cortical layer and study the influence and effects of cortical layer thickness and lack thereof on the acetabular cup stability. Due to the simplified setup, the influence of muscle tissue and ligaments on the deformation behavior and load response was neglected as well, which is in agreement with what is commonly done in the literature (Hao et al. 2011; Clarke et al. 2013). However, it has been shown that muscles and ligaments have to be taken into account when analyzing the stress distribution inside the acetabulum (Shirazi-Adl et al. 1993), which is beyond the scope of the present study. Future studies should consider more realistic and physiological geometries and boundary conditions to improve the accuracy of the numerical results and provide more reliable estimations of implant stability.

No actual bone ingrowth or bone remodeling was modeled and homogeneous osseointegration over the whole bone-implant interface was assumed. In reality, only certain parts of the bone-implant interface are osseointegrated depending on the contact conditions, such as contact stress, micromotion, and initial gap. In addition, initial gaps after surgery might be filled with bone tissue during the healing phase and thus increase the contact area and bonding strength over time. In future works, the presented debonding models should be coupled with suitable osseointegration models and bone remodeling algorithms (Caouette et al. 2013; Mukherjee and Gupta 2016; Chanda et al. 2020; Martin et al. 2020), to achieve a more reliable assessment of implant long-term stability. These models should account for pressure- and micromotion-depended bone apposition and resorption, as well as changes in the contact gap and the maturation of new bone tissue, for example, by changing the bone’s elastic properties with respect to healing time. Furthermore, due to bone growth and the change in elastic properties of the bone during osseointegration and remodeling, the stress inside the bone changes during the healing process and might be significantly different after healing compared to the state directly after surgery. As the change in stress can significantly affect secondary stability, remodeling-related effects should be considered in future works.

5.4.2 Contact modeling perspectives

This work neglects the roughness of the implant surface and of the reamed bone cavity. While a simple modeling of rough surfaces by adjusting $\mu _\text {b}, \mu _\text {ub}, t_0$ is possible, the explicit modeling of rough surfaces should be considered in future works, as surface roughness affects initial stability and osseointegration and thus also long-term stability. Furthermore, due to the rise of additive manufacturing in implantology, complex implant surface topologies become more and more relevant and should be studied.

The CZM in Eq. (13) is modeled with a sharp drop in $t_\text {n}$ at $g_\text {n}=g_\text {b}$. Future studies should explore CZM models that depend on $\phi$ instead of $\phi _0$ and have a smooth decline in $t_\text {n}$ for $g_\text {n} > g_\text {b}$.

The removal force/debonding torque was chosen as determinants of long-term stability. The stress distribution could be used as another determinant, as is done in other works (Janssen et al. 2010; Rourke and Taylor 2020). However, the stress distribution inside the bone changes during healing and osseointegration, as the mechanical properties of the bone tissue change when the new bone tissue mineralizes. This makes comparisons of stress fields of initial stability and secondary stability scenarios difficult, when this temporal change is not accounted for.

As in previous studies by our group (Raffa et al. 2019; Immel et al. 2020, 2021a), a quasi-static configuration was considered, and all dynamic aspects were neglected, similarly to what was done in comparable works (Spears et al. 2001; Le Cann et al. 2014; Raffa et al. 2019). Note that a previous study focuses on the insertion process of an acetabular cup implant by considering dynamic modeling (Michel et al. 2017), which is important when modeling the insertion by hammer impacts. Furthermore, the stick-slip results with the present extended contact model (see Figs. 22 and 24) indicate that dynamic simulations become necessary when considering high frictional and adhesive forces.

5.4.3 Experimental perspectives

Model EMC depends on two additional physiological parameters $t_0, g_0$ that can be determined based on some of the few experimental results available in the literature (Rønold and Ellingsen 2002). However, to the best of our knowledge, no suitable measurements have been obtained for osseointegrated acetabular cup implants yet, which is why we calibrated our models with measurements for coin-shaped implants instead. Future experimental tests of osseointegrated implants under mixed mode or mode III debonding under constant tension (as presented in Section 3.1) can provide important insight on the adhesive behavior of the osseointegrated interface to calibrate and validate the proposed debonding models.

The strong influence of biological as well as mechanical factors and the bone geometry on the long-term stability make validation of the presented numerical models difficult. At present, experimental studies that provide sufficient information on the behavior and stability of the partially osseointegrated bone-implant interface are lacking in the literature (Helgason et al. 2008). We suggest to perform mixed mode debonding and mode III debonding under constant tension, as demonstrated in Section 3.1. These results would provide important information on the debonding behavior of osseointegrated interfaces and allow to further calibrate and validate the extension of the modified Coulomb’s law. Further computational studies cannot reliably provide more insight on the in vivo behavior, as the level of sophistication of the models is beyond the point of verification with current in vivo, ex vivo, and even some in vitro measurement techniques (Taylor and Prendergast 2015). Therefore, it becomes more and more difficult to reliably assess the performance of numerical models for the bone-implant interface. If FE models are to be trusted and accepted by clinicians, they need to demonstrate that they are capable of predicting realistic in vivo behavior. Thus, further development of experimental measurement techniques and quantification of relevant biomechanical metrics (e.g., stress-strain behavior, micromotion, friction, adhesion, and debonding under tension) is essential to provide the data necessary to develop and improve numerical models. However, the development of new and more accurate experimental machinery and techniques that are able to provide the necessary data is difficult and time consuming and provides a constant challenge. While experimental and numerical methods keep improving, a certain acceptance that FE studies may not be representative of the in vivo conditions yet but are an approximate model, needs to be established.

6 Conclusion

This work presents a new extended debonding model for the bone-implant interface, which can describe the debonding behavior of osseointegrated acetabular cup implants and thus assess their stability. In addition to the modified Coulomb’s law of Immel et al. (2020), it includes a cohesive zone model in normal direction and adhesive friction in tangential direction.

The modified Coulomb’s law and its extension show that friction and adhesion increase the pull-out force/debonding torque of osseointegrated implants and thus are relevant for long-term stability. Furthermore it is shown that, while osseointegration increases implant secondary stability, a sufficient primary stability remains crucial for long-term stability, which is in agreement with the literature. These findings underline the importance of the development of surgical decision support systems such as the surgical hammer instrumented with a force sensor to measure the displacement of an osteotome or implant and determine when full insertion has taken place (Michel et al. 2016a, b; Dubory et al. 2020; Lomami et al. 2020) or contactless vibro-acoustic measurement devices that can monitor implant seating (Goossens et al. 2021). Coupling simulations of initial stability, subsequent osseointegration and bone remodeling, and long-term stability and debonding can provide more reliable assessments of implant stability and aid in implant conception and individual patient treatment. Furthermore, a future detailed study would be able to answer how cyclic loading affects the bonding state of the interface during and after healing. Last, this work provides directions for important experimental testing of osseointegrated coin-shaped implants. Mixed mode debonding and mode III debonding under constant tension could provide important information on the debonding behavior of osseointegrated interfaces and allow for further calibration and validation of the proposed contact models.

Appendix A Convergence study

To analyze the convergence behavior of the modified Coulomb’s law applied to the acetabular cup implant, five different finite element meshes for the bone block were constructed with increasing refinement of the elements in x and $y-$direction. The number of elements for each mesh are shown in Table 7. The reference case ($E^*_\text {b}=0.2$ GPa, $\mu ^*_\text {b}=0.3$, $I\!F^*=1.0$ mm) was chosen as the parameter set, and the maximum normal pull-out force $F^{\text {max}*}_z(\phi _0=1)$ was chosen as the target value.

Table 7 Number of elements of the finite element meshes

Figure 29 shows the convergence behavior of the maximum pull-out force. It decreases with increasing number of elements, while the computing time increases exponentially. The estimated exact value of $F^{\text {max}*}_z(\phi _0=1)=659$ N. The corresponding relative error of the maximum pull-out force of the meshes 1 to 4 is 10.5, 6.1, 2.1, and 0.4%, respectively. Due to the focus on contact problems, the mesh was only refined in $x-$ and $y-$direction. The number of elements in $z-$direction remains 5, as a previous convergence study showed no measurable improvement with further refinement in $z-$direction. Due to the lack of comparable experimental data and the amount of computations (58 with standard Coulomb’s law, 435 with each of the two contact models) necessary for Section 4, mesh 3 was deemed to have a reasonable ratio between accuracy and computing time and was thus chosen for all computations in Section 4.

Fig. 29

Maximum normal pull-out force $F^{\text {max}*}_z$ and corresponding computing time with respect to number of contact elements on the bone block. Mesh 3, which is used for the computations in this work, is marked with $\star$

Table 8 Average computing time, average number of Newton–Raphson steps, and corresponding number of load steps for the simulations with the modified Coulomb’s law (MC) and its extension to adhesive friction (EMC) for the different loading cases and mesh 3

The average computing time with both contact models for the different loading cases performed with mesh 3 is listed in Table 8. It should be noted that the chosen parameter combination has an influence on the computing time. Parameter combinations that produce high pull-out forces also require more computing time.