Journal of Biomedical Science and Engineering
Vol. 5  No. 5 (2012) , Article ID: 19394 , 16 pages DOI:10.4236/jbise.2012.55030

A model of bone adaptation as a topology optimization process with contact

António Andrade-Campos, António Ramos, José A. Simões

Department of Mechanical Engineering, Centre for Mechanical Technology & Automation, University of Aveiro, Campus Universitário de Santiago, Aveiro, Portugal


Received 7 January 2012; revised 22 February 2012; accepted 2 March 2012

Keywords: Bone remodelling; topology optimization; finite element; contact; osteointegration


Topology optimization is presently used in most diverse scientific, technologic and industrial areas, including biomechanics. Bone remodelling models and structural optimization has mutually provided inspiration for new developments in biomechanics and biomedicine. Considering that bone has the ability to adapt its internal structure to mechanical loading (Wolff’s law and Roux’s paradigm), it is possible to model the behaviour of the bone structure by the use of a topology optimization methodology whose optimization variables can be the relative densities and the orthotropic directions. In this work, the internal bone adaptation of a proximal femur is considered. The bone-remodelling scheme is numerically described by a time-dependent evolutionary procedure with anisotropic material parameters. The remodelling rate equation is obtained from the structural optimization task of maximizing the stiffness subject to a biological cost associated with metabolic maintenance of bone tissue in time. The situation of multiple load conditions is considered for a three-dimensional finite element model of the proximal femur. The bone density distribution of a real femur is used as the initial design for the onset of the remodelling mechanism. Examples of bone adaptation resulting from load changes are presented. The three-dimensional finite element model of the proximal femur with initial bone density distribution was adapted to implant a cementless stem. A remeshing technique is used to assign the bone relative density distribution to the new geometry and mesh. The time adaptation of the bone is assessed considering contact with friction at the bone-stem interface. Results of bone density evolution and osteointegration distribution are obtained.


Since the works of Michell and Bendsøe and Kikuchi [1], topology optimization has become an effective design methodology to obtain lighter and efficient structures. Until today, diverse topology optimization methods have been developed and presented. The methods based on the SIMP (Solid Isotropic Microstructure with Penalization) technique, originally introduced by Bendsøe [1], have had large acceptance due to its efficiency and simplicity.

Currently, topology optimization methods are being used in biomechanics and biomedical engineering, namely in bone-remodelling schemes and bone-adaptation models. Considering that bone has the ability to adapt its internal structure to mechanical loads and to changes of the load environment (paradigm of Roux and/or Wolff law [2]), it is possible to mathematically model the behaviour of the bone structure through a topological optimization methodology whose optimization variable is the relative or apparent density.

Many different mathematical and phenomenological models have described the remodelling process of the internal structure of bone. These models have in common the definition of an equilibrium state based on energy levels, on stress levels or on a reference state density. Generally, these models present an equation of the evolution of the non-equilibrium state for the equilibrium state through the change of state variables as, for instance, the local densities [3]. The remodelling process occurs when the bone senses a stimulus originated from a change of external loads [2,3].

The remodelling process also occurs when an orthopaedic implant is present. The insertion of a stem changes the equilibrium state and, consequently, compels the bone to achieve a new equilibrium state. Generally, it is observed a loss of bone when an orthopaedic implant is inserted in the bone structure. This loss of bone can lead to bone fracture and to implant loosening. The arthroplasty postoperative process is an evolutionary process. To analyse numerically this process it is necessary to consider the bone remodelling time process together with the changes of the implant-bone interface conditions.

In 1993, Sadegh et al. presented a bone remodelling theory to predict the ingrowth of bone into implant cavities. Considering a very small dimensional scale, the model predicted the bone ingrowth in screwed implants. More recently, Fernandes [4,5] and Bagge [3] presented models of bone adaptation of the proximal femur. Both models consider a three-dimensional finite element model that evolves as a topology optimization scheme. Although the presented models could achieve the bone density distribution, none of them consider the inclusion of a prosthetic device. Both analyses use a single hexahedral finite element geometric model.

The osteointegration1 process (also called osseo-integration) over an implant surface was also investigated using a finite element model by Spears et al. in 2000. The contact interface considers contact elements that account large sliding contact movement and total bonding was considered possible for elements with a micromotion of less than 40 μm.

The use of elastic springs to model the bone implant interface was introduced by Egan and Marsden in 2001. The stiffness of an interface bonding spring is progressively increased at the point where the deformation (proportional to the displacement) induced by the external loads remained below a given threshold.

A similar progressive interface bonding has been presented by Fernandes et al. [6,7] together with the bone remodelling model of the whole femur after total hip replacement. Changes in the interface conditions in conjunction with a global optimization criterion based in Microstructural homogenisation and a cost minimization was used. Although the previous model considers not only the bonding process but also de-bonding process, the time was not correlated and the loading conditions applied were the same as the healthy bone that, for an implanted femur, does not correspond to the reality [7-9]. Additionally, the numerical solutions were obtained starting with a homogeneous solution of the trabecular bone with a relative density equal to 30%.

A purely biomechanical model, intended to predict the long-term secondary stability of the implant starting from the biomechanical stability immediately after the operation, was also presented by Viceconti et al. [10-13]. This model, also compared with a simpler one [13], uses a finite element model discretized in three-dimension hexahedral elements and a continuous rule-based adaptation scheme as a dynamic system. Although this model describes the whole bone-implant model, it does not consider the change in the microstructural apparent density of the bone that can affect the interface and overall stiffness of the bone.

A recently model for living interfaces with bone implants was introduced by Moreo et al. [14]. This model reproduces the evolutive behaviour of bony interfaces accounting with bone ingrowth and deterioration by means of a finite element continuum damage mechanics formulation. It was observed that stiffer stems (such as stainless steel and Titanium stems) improve the primary fixation comparing with softer ones (Polyacetal stems). Even though the results of the simulation agree reasonably with clinical observations, the model presented exhibits some limitations associated to the lack and uncertainly of the different model parameters.

Almost all the bone ingrowth models that use finite element analysis consider a constant finite element discretization and the bone ingrowth is characterized by an increase of stiffness. Recently, Reina et al. [15,16] tried to model bone ingrowth in distraction osteogenesis considering mechanical and biological factors. Although this model considers a remeshing technique that allow the creation of new elements, the analysis was limited to a two-dimension rectangular geometric model and some biological issues are not clear.

Also recently, the bone remodelling model based on topology optimization was extended to a hierarchical model by Coelho et al. [17,18]. This hierarchical model characterizes both the macro and microstructure of the bone by means of a multiscale asymptotic homogenisation technique. However, this study does not consider an insertion of any prosthetic device and the implant-bone interface is not studied.

The development of biomechanical models is a hard task. The major difficulty in the development of such models is the lack of experimental quantitative results, specifically in vivo human data results. As little or no experimental methods are currently available to investigate the relationship between primary stability and osteointegration in vivo, all published studies on this topic are based on numerical simulations. However, some in vivo and ex vivo observations are available for the portion of implant surface that has some osteo-bond but not necessary fully bonded (or fully osteointegrated) [19,20].

The objective of this work is to present a model of bone adaptation that accounts the time and the change of stimulus. These stimuli are represented by mechanical loads. In order to obtain a realistic bone remodelling simulation, the initial structure needs to be approximately equivalent to a femur. Thus, a multiload topology optimization process is applied with the purpose of finding the optimum initial distribution of densities corresponding to the proximal region of a healthy femur. Different finite element meshes are used.

The bone adaptation model is also applied in the case of an implanted femur considering the bone and implant stem surfaces in contact. The bone ingrowth process is simulated through a numerical model able to modify the contact conditions during the remodelling process. Therefore, a global optimization criterion with contact formulation and bone ingrowth is used. The contact conditions are modified based on the stem-bone relative displacement and stresses. The numerical initial conditions were retrieved from the femur bone just before the medical surgery. A remeshing technique is employed in order to transfer the properties from the intact femur to the implanted-femur.


2.1. Material Description

Due to its porous and adaptive structure, bone presents an orthotropic behaviour whose orthotropic coefficients change with the bone region and evolve with the remodelling process. These aspects imply that modelling the bone microstructure is an extreme hard task. Many boneremodelling models use a power-law relation between the apparent density and the isotropic elasticity modulus of cancellous bone [2]:


where E is the isotropic elasticity matrix, ρ is the density and Ebone is the stiffness matrix of the bone structure. The previous relation has the advantage of being simple and computationally very inexpensive, but anisotropy is not considered. The power-law relation is similar to the use of the SIMP method in the bone adaptation process.

Other bone-remodelling models use polynomial functions that characterize the stiffness matrix coefficients as a function of the density and can be written as [3]:


This model, which already accounts for anisotropy, considers a constant microstructure whose material properties are symmetrically cubic. This approach is also computationally inexpensive.

An orthotropic material model can be obtained using homogenization methods such as, for example, the asymptotic homogenization method. In this method, the cancellous bone structure is idealized by a microstructural cell assuming that the material can be reproduced by a periodic repetition of the cell. Taking in account the shape of the bone microstructure, the material properties are a function of the mechanical characteristics of the solid material and the volume content in the cell [21-24]:


where x and y are the properties at the macro (bone) and microscale (cell), respectively, and ε is the asymptotic homogenization constant that corresponds to the relation of the characteristic dimensions between the micro and macroscale. The use of two scales related with a small dimension parameter corresponds to a particular application of the multiple-scale method. This method, that assumes a total periodicity of the cells, requires the calculation of the cell microstructure properties (through a finite element analysis) for each integration point of the macrostructure, leading to a very computationally expensive process. Although this method considers full anisotropy, generally, the used unit cells take parallelepiped geometries that are quite different from the real trabecular bone microstructure, and can be considered non-realistic.

The difficulty of description of the orthotropic behaviour of the bone material is well-known and evidenced by the numerous mentioned models. In this work, despite the disadvantages presented, the asymptotic homogenization method is used. To achieve the material properties as continuous functions of the relative volume fraction, the cell structure, shown in figure 1(a), has been used. However, the material properties at tissue levels are assumed isotropic [3]. In order to decrease the CPU time required by this method, the homogenized material properties are fitted to polynomials (see figure 1(b)). For more details concerning the implementation of this method, see [22-24].


Figure 1. The microstructure of the bone is represented by a periodic unit cell. The material properties are calculated through the asymptotic homogenization method and fitted to polynomials functions of the relative volume fraction.

2.2. Bone Remodelling as an Optimization Process

The initial distribution of the apparent density of bone needs to be realistic. Therefore, it can be obtained by a CT scan (computerized tomography scan) or by stiffness maximization of the structure [3] , if Wolff’s law is considered. The maximization of the stiffness of a structure corresponds to the minimization of the strain energy. Assuming that bone adapts to its mechanical environment, the bone remodelling process consists of the computation of bone relative density at each point by the solution of an optimization problem formulated in the continuum universe. Numerically, the continuum can be approximated by finite quantities. Therefore, the process of finding a reasonable structure can be an optimization task that considers the three-dimensional domain discretized in finite elements and whose optimization variables are the relative volume fractions x of each element. These variables are constants in each element and take the values of 0 (void) to 1 (total density, cortical bone). The optimization process is subject to a volume constraint and to the finite element equations for each load. The optimization problem can be defined as


with the following constraints and conditions


and, in the case of multiple solids, contact conditions for the bone/stem interface


Πj is the strain energy for the jth load case, w represents the weight factor and Ne and Nl are the number of elements and load cases respectively. The maximum admissible volume and the volume in each finite element are represented respectively as and. In each iteration, a finite element analysis (FEA) is performed. The global stiffness matrix K and the load vector F are known for each load case and the global and element displacement vector, U and ue, respectively, are found.

The contact problem with Coulomb friction in eq.6 considers a penalization approach and the inclusion of an elastic reversible term for the tangential displacements (or velocities). In this equation, the subscripts N and T denote normal and tangential properties, respectively. Therefore, tN and tT are the normal and tangential contact forces applied on the unilateral contact interface. The first expressions in 6 define the normal and tangential displacements of the contact boundary through projections of the displacement u on the boundary normal n. The function g is the actual gap between the two bodies in contact while s is the initial gap. The second line defines the penalisation-type constitutive equation for normal traction, where the operator represents the positive value of the argument. The tangential velocity is divided into the elastic and sliding component (as in perfect elastoplasticity). The sliding velocity is proportional to the derivative of the sliding potential f defined in the last line of 6 (coulomb model).

The FEA can be executed with the aid of an exterior FE program. FE programs ABAQUS® [25] or MSC. MARC® [26] were used in this work. In order to change the properties in each element, a user routine UMAT (ABAQUS 2006) was developed for the commercial program ABAQUS®. For the MSC.MARC® a user routine ANELAS [26] was used. The contact problem is solved using a surface-based contact model being the stem the master surface and the bone the slave surface. As mentioned, the penalty method is used for the tangential contact condition.

It is assumed that the design variables x (relative density) are constants within each finite element. The optimization problem, achieved through an iterative procedure, is solved using the MMA method (Method of Moving Asymptotes) [1] or the optimality criteria (OC) [1,3,27]. The OC algorithm is easy to implement and very efficient for a single constraint. For multiple constraints, the OC cannot be used and, therefore, this is the main motivation for switching to MMA. On the standard stiffness optimization problem, the MMA can be carefully tuned to be more efficient than the OC. However, in the standard settings they are very similar. The OC method, based in the optimum Lagrange function, uses a fixed-point updating scheme. It can be written as


where η is a numerical damping parameter and k represents the iteration number. It is expected regions with high energy with low stiffness. Therefore, in order to avoid numerical problems, the previous fix-point type updating scheme for the densities should have boundaries. Considering that xmin = 0 to avoid singularities in the finite element stiffness matrix, it was defined xmin = 0.001. Eq.7 is hence computationally defined as:


The Be parameter is defined as


where λ is the Lagrange multiplier parameter that results from the volume constraint. A local optimum is reached if Be = 1 for densities xmin < x < 1. This update scheme adds material to areas with specific strain energy higher than λ (Be > 1) and removes it if the energy is below this value [1]. The variable η is a tuning parameter and ξ a move limit. Both η and ξ control the changes that can happen in each iteration. They can evolve with the procedure in order to enhance the efficiency of the method.

The update of xe depends on the present value of the Lagrange multiplier λ and, as a consequence, should be adjusted in order to satisfy the volume constraint. The value of λ can be found by a bi-sectioning algorithm considering that the total volume, with the actual relative densities, is a continuous and decreasing function of the multiplier λ.

If a direct implementation of a topology optimization method is made, patches of checkerboard pattern appear often in solutions. Restriction methods reduce or remove the checkerboarding in some cases, but not in all of them. With the aim of preventing the numerical problem of checkerboard, the filter method suggested by Sigmund [1] is used. In this method, the sensitivity of each element Be is affected by the sensibilities of the neighbouring elements. The change in the sensitivity is proportional to the distance between the elements. The scheme works by modifying the element sensitivities of the compliance as follows [1]:


The convolution operator (weight factor) is calculated as


where the operator is defined as the distance between the centre of element e and the centre of element f. Note that the convolution operator is zero outside the filter area and decays linearly with the distance from element f. The modified sensitivities of Eq.10 are used in the optimization problem using the OC or MMA update method.

In the case of the bone adaptation process, the material behaviour must be formulated progressively with time. Therefore, eq.7 is used as the remodelling rate equation but it was considered that each iteration is a time step. Hence,


where t represents the time. In the previous equation, if the Lagrange multiplier is kept constant during the remodelling process, bone can form and resort material. Therefore, in the remodelling process, λ is denoted osteometabolic cost and is related with the maintenance of bone mass. This parameter, 0 < λ < 5, should be dependent of the age, degree of bone osteoporoses, etc., i.e.,


when the strain energy sensitivity is not in equilibrium with the element constant there is a change in the bone relative volume. In each element, a response Be greater than the one hypertrophy is assigned and Be less than one atrophy is assigned [3]. The equilibrium level, where no resorption or apposition takes place, is achieved when the elements with intermediate relative volume fractions have Be equal to the unity.

In this work, the limit of the change in volume fraction in each time increment was based on the work of Jee [28]. This author verified that it can take place a maximum bone turnover rate of 7.6%/year (≈ 0.021%/day), considering a total bone loss in 13 years. As a result, if the stimuli in an element leads to a local bone changing greater than the limit, the value of Be is numerically damped with aid of η [3]. The parameter η (0 < η < 1), that influences the sensitivity of all the elements in a proportional way in each time increment, can be found by a bisection algorithm.

2.3. Osteointegration Model

The osteointegration model intends to reproduce the adaptation of the bone considering the insertion of a stem. This model simulates the changes of the bone/stem interface conditions.

During surgery, a part of the femoral bone is removed and stem is inserted into the femoral canal [29,30]. During the first days, the phenomenon of integration of a biologic bone material in the non-biologic metallic part is observed. If a coated uncemented prosthesis is inserted, it is expected bone ingrowth in the coated zones after some postoperative time (see figure 2).

In this work, a contact formulation is used together with a methodology to detect where bone ingrowth exists and evolutes. At each time increment, the contact conditions are updated considering a relative displacement criteria. These criteria are based on the studies performed by Fernandes et al. [6,7].

The osteointegration model aims the simulation of the bone/stem interaction. After implantation of the prosthesis, bone begins to attach to the coated stem wall and the prosthesis retains stability. However, bone ingrowth can be inhibited in certain zones due to high relative displacements.

Initially, the osteointegration is considered null (see figure 3). In this post-operative condition, contact with friction at the coated surfaces and frictionless contact at the non-coated surfaces is considered. For each time step, the relative displacements between the bone and the stem are computed. If an effective contact between bone and implant is verified (therefore, the normal contact force is greater than one) and the tangential relative displacement is inferior to a threshold value, osteointegration increases. Numerically, the osteointegration is progressively augmented with the use of connector elements of the Cartesian type [25], increasing the connector element stiffness. The conditions for progressive osteointegration can be resumed by the following conditions:

I. Contact within the stem and bone material (no separation of the surfaces): Fcontact > 0 → COPEN = 0;

II. The relative tangential displacement of the surfaces in all directions should be null (or inferior to a threshold value): CSLIPsurf < CSLIPlimit.

Figure 2. Artificial hip replacement of an uncemented prosthesis. A bone ingrowth is observed in the coated zones of the stem increasing the stiffness of the bone/ stem interface.

Figure 3. Artificial hip replacement of an uncemented prosthesis. A bone ingrowth is observed in the coated zones of the stem increasing the stiffness of the bone/stem interface.

In the following time steps, if the osteointegration conditions are not verified, the bone/stem connection is ended and the initial condition is replaced. The full integration condition is achieved when the osteointegration is augmented for a chosen number of consecutively time steps. The number of time steps (iterations) necessary to reach the full osteointegration connection characterizes the bone ingrowth rate. In this case, the contact condition is transformed into a type of glue connection, where the relative displacements are constrained. However, the full osteointegration condition is verified for every step and can be destroyed. That is the case of high levels of stress. Numerically, to maintain the full osteointegration condition, the following condition must be also verified:

I. The contact stresses within the stem and bone material are inferior to some limit values (otherwise it will break the connection): σcontact < σlimit → CPress < CPressLimit and CShear < CShearLimit.

If the previous condition is not verified, the interface and contact conditions must be updated to the condition of non-osteointegration (initial condition).

3. Validation—Structural Problems with Contact Conditions

The formulation presented from eq.4 to Eq.10 is validated through structural problems introduced by other authors. Four different structural problems with contact conditions were addressed and herby presented. The microstructure is considered isotropic for these problems.

Both FE programs Msc.Marc® and Abaqus® were used in the following problems. The final results obtained were similar. However, it was observed that the analyses with the aid of Msc.Marc® are faster. The main advantage of using the latter software is the possibility of creation of an autonomous single execution file that can be directly runned in other computers. As these problems are single constraint topology optimization problems (single volume constraint), both OC and MMA methods can be used with similar results and efficiency.

3.1. Rectangular Domain with Two Contact Surfaces

The problem of a rectangular domain was addressed by Desmorat [31]. A 2D case of a 2 × 1 rectangle with two internal contact zones, fixed at one end and loaded vertically at the centre of the other end is considered (see figure 4(a)). Figure 4(b) presents the optimization results on the coarse mesh structure with the two contact surfaces taking into account the effective contact or noncontact cases. By opposition with the result if no contact surfaces were assumed, this result shows asymmetric material distribution, an open contact zone on the upper contact surface and an effective contact zone on the lower contact surface. This observation fully agrees with the results presented by Desmorat [31].

3.2. Linked Bars

The linked bars problem, described in figure 5(a), considers the case of two 3 × 0.5 linked bars with two internal contact zones, fixed at one end and loaded vertically on the other end (see [31] for more details). Figure 5(b) presents the optimization result on the deformed structure for a volume constraint of 20% and taking into account the effective contact cases. Again, the results obtained are in agreement with the results of Desmorat [31] and, comparing with the solution without contact conditions, a concentration of material on the left of the link line and a different material distribution in the central part of the structure can be observed.

3.3. L-Shaped Test

The purpose of this benchmark is the verification of the behaviour in the interior corner of a L-shaped domain [32]. The geometry and boundary conditions are displayed in figure 6(a). A single load (upper traction of 100 MPa) pushes the structure upward while the contact condition is defined on the lateral support. The optimization result corresponding to the minimum strain energy is presented in figure 6(b). This result is in agreement with the results presented by Fancello [32]. It shows a concentration of material at the interior corner of the L-shaped domain and in the load surface.

(a) (b)

Figure 4. Artificial hip replacement of an uncemented prosthesis. A bone ingrowth is observed in the coated zones of the stem increasing the stiffness of the bone/stem interface.

(a) (b)

Figure 5. Linked bars with two contact surfaces: (a) Geometry and boundary conditions [31]; (b) optimization result for a coarse mesh and a volume constraint of 20%.


Figure 6. L-shaped test with contact surface: (a) Geometry and boundary conditions [32]; (b) optimization results on the deformed structure for a volume constraint of 40%.

3.4. Flexible Pinned Union

The union of three blocks with a flexible pin is modelled (see figure 7(a)). The 3D lateral blocks are constraint in the left edge and a horizontal load of 1000 N is applied to the central block. The lateral block dimensions are 75 × 60 × 6 mm, the central block is 90 × 60 × 12 mm and the dimensions of the cylindrical pinned are 20 × 25 mm [33] . Between the blocks, in the lateral surfaces, there are 30 mm diameter rings with 0.5 mm of thickness. All the solids are deformable but only the blocks belong to the design domain. Due to symmetry conditions, only a quarter of the geometry is modelled. The material employed is steel except for the rings that are of polymeric material (elasticity modulus of 2 GPa). Frictionless contact conditions were considered between the pin, the rings and the blocks.

Figures 7(b) and (c) show the results obtained for a volume constraint of 30%. As expected, the material is concentrated in areas where strain energies are higher, i.e., the material that wraps the pin and the distribution is made into the constraints (the lateral block case) or into the applied load (the central block case). These results are similar to the results presented by Folgado [33] validating the presented formulation.


The formulation presented in Section 2.3 is validated through structural problems introduced by the authors.

Progressive Contact Integration of a Compression Specimen

This problem characterizes the behaviour of a specimen subject to uniaxial compression loading. The contact surface of the specimen has a coated area, assembled by the area of four elements that allows material integration (see figure 8). Initially, while the material integration is null, the behaviour of the compressed material is equal to an ordinary compressed specimen. The contour of von Mises equivalent stress in the compressed specimen and in the exterior material can be seen in figure 8 for the initial tine step. For this step, the tangential relative displacement for all nodes that belong to the material integration surface is 1 × 107 < CSLIP < 1.2 × 103 and the contact pressure is positive, evidencing that there is an effective contact between the two surfaces. Considering the threshold value for the relative tangential displacement of CSLIPLimit = 5 × 105 m, only the nodes that present relative tangential displacement lower than this value will start material integration and their connector stiffness will increase. The values of CPRESSLimit and CSHEARLimit are 0.8 and 35 MPa, respectively.


Figure 7. Flexible pinned union test: (a) Geometry and boundary conditions [33] ; (b) optimization results on the deformed structure for a volume constraint of 30%. Elements with relative densities inferiors of 0.1 were removed; (c) a of the solution of the optimization domain.

Figure 8. Compression specimen with an integration surface: (a) Geometry and boundary conditions; (b) Equivalent stress contour on the deformed geometry in the initial step, without any material integration.

By the charts of figure 9, it can be seen that material integration (osteointegration if the specimen is considered as a bone) occurs at nodes number 33, 34 and 35 (nodes in the middle of the width of the coated surface, as represented in figure 10). These nodes achieve full integration at the 20th time step and maintain it to the end of the test.

In the first 6 time steps, the integration of nodes 20 and 46 (the extreme bottom nodes of the coated surface) oscillates due to the checking of the integration conditions. During that time, the relative tangential displacement decreases towards step 5 but oscillates between steps 5 and 8 (see details in CSLIP chart). The material integration conditions are successively fulfilled after step 6 and the full integration is achieved in time step 26. To note that the increasing of stiffness of the connector of nodes 33, 34 and 35 aid the decreasing of relative tangential displacement of theses nodes.

Nodes 21, 22 and their symmetric never attain integration. The relative tangential displacements (CSLIP) chart in figure 9 shows that, although with some oscillations, the relative tangential displacements never decrease to a value inferior to the threshold, therefore never checking the integration conditions.

This problem characterizes well the possibility of three types of integration nodes: i) nodes that immediately start to integrate; ii) nodes whose integration depends of the integration of other nodes and iii) nodes that integration is never achieved. These occurrences take place in the osteointegration process, particularly in the bone/stem interface.

Figure 9. Time evolution of osteointegration and relative tangential displacement (CSLIP) of the nodes belonging to the coated area of figure 10.


Figure 10. Compression specimen with an integration surface: (a) Contour of the area with osteointegration. The value of one correspond to full integration; (b) displacement contour on the deformed geometry in the final step with material integration.


5.1. Geometric Model and Domain Discretization

The geometry of the proximal femur and the boundary conditions used in this work are represented in figure 11a. The multiple loads applied correspond to 2.5 times of the medium weight of a human body (approximately 700 N) and can be seen in table 1 [21]. The model is constrained in all three directions at the nodes of the outer rim of the bottom element layer. Figure 11(b) to f show the different domain discretization used in this work. Figures 11(b), (c) and (d) show the original geometry discretized in 1158, 5224 and 9408 hexahedral finite elements, respectively. Figures 11(e) and (f) present the discretization of 18001 and 143904 tetrahedral finite elements, respectively. Both elements are solid element with linear interpolation. It was assumed that bone tissue has the mechanical properties of compact bone. Therefore, the dense compact bone corresponds to a cellular material with relative density equal to 1 and trabecular bone has values less than 1. The elasticity modulus for dense compact bone is 20 GPa. It was also assumed a Poisson ratio of 0.3.

In the optimization process, the three load cases have the same weight and, to ensure consistency between the energy of one load and the weighted energy, it must hold that [3] . The starting point for the optimization process is uniform relative volume distribution.

Table 1. Load cases for the intact and implanted femurs.


Figure 11. (a) Geometric model and boundary conditions of the proximal part of a femur. Discretization in (b)-(d) hexahedral and (e)-(f) tetrahedral finite elements.

5.2. Initial Distribution of Material

Results of the optimization process using three values of volume constraint are given in this section. It was considered different volume fractions (40%, 55% and 70%) and it was taken in account the domain discretizations represented in figures 11(b) to (f). The optimization process that leads to a similar result to the apparent densities distribution of a human femur was the process using 55% volume fraction. The use of different meshes results in analogous general results, even though particular differences can be observed. Figure 12 shows the material distribution after 100 iterations with the total volume fraction of 55%. The distribution of the relative density is interpolated from the elements to the nodes.

The results of figure 12(a) to f show agreement between the numerical density distribution and the apparent density of the real bone. Although the shaft cavity is well represented, the densities distribution in the femoral head should be more homogeneous for the case of tetrahedral elements. This fact can be attributed to the applied load cases and the type of elements used and their importance in the numerical process.


Figure 12. Relative density distribution of the trabecular bone using a volume fraction of 55% (E ≈ 1.185 × 104 N mm and λ = 1.284 N·mm2). (c)-(d) Longitudinal and (e)-(f) transversals cuts of the femur’s model for hexahedral and tetrahedral finite elements, respectively.

5.3. Bone Adaptation by Stimuli Changes - Results

After the achievement of the initial relative densities distribution, the bone adaptation to the stimuli variations is analyzed. The osteometabolic value (λ) is considered constant and equal to the value found in the initial analysis (λ = 1.284 N·mm2). It is used the initial densities configuration as well. In the first case, and to illustrate the time dependence of the bone adaptation process, the load cases (stimulations) are multiplied by 1.5 (corresponding to approximately 4 times the weight of the human body) during 9000 days. In the second case, the stimuli are reduced to 2/3 during 3000 days and, afterwards, increased 1.5 times. This case demonstrates bone adaptation to stimuli inversion. The third case should show the reversibility of the process. From the first case solution (figure 13(a)), the stimulations are then kept as the initials during more 24,000 days.

In the three analyses, a maximum bone turnover rate of 0.021%/day is defined [28]. A sequence of multiload cases is made in order to consider a regular day. The sensitivities are therefore averaged [3]. The number of times the bone is loaded has no influence in the numeric result when the limit of the bone changing is reached. In that case, a h value is found to bound the maximum bone turnover rate and the sensitivities are numerically damped (see eq.12). Each iteration is considered a day.

Figures 13(a), (b) and (c) show the time evolution of the bone volume fraction for the three cases, respectively. For case 1, figure 13(a) shows an increasing of ≈9% of the total volume within the 9000 days. It should be notice that there are small oscillations around the volume-days curve. This case shows the influence of time and stimuli in the bone adaptation process. Figure 13(b) shows an inversion of the bone adaptation process with the inversion of the load history. At first, the volume has decreased 5% and, when the body weight has raised to four times, the volume has increased to 63%. The reversibility of the process is studied in case 3. It is known that bone tissue adaptation is an irreversible process. However, the numerical process used in this work shows a reversible process.


Figure 13. Time evolution of the volume fraction. (a) Case 1: the stimuli (loads) are multiplied by 1.5 during 9000 days. (b) Case 2: the stimuli (loads) are reduced by 2/3 during 3000 days and, subsequently, increased 1.5 times. (c) Case 3: from case 1, the stimuli (loads) are held as the initials during 24,000 days.

5.4. Numerical Surgery and Remeshing

For the continuation of this work, a numerical surgery was performed in order to obtain an implanted femur. During real surgery, a part of the femoral bone is removed. First, the femoral head is removed by cutting through the femoral neck with a power saw. Then, special rasps are used to shape the hollow femur to the exact shape of the metal stem of the femoral component. Once the size and shape are satisfactory, the stem is inserted into the femoral canal. In the uncemented variety of femoral component the stem is held in place by the tightness of the fit interference into the bone (MMG 1999). These operations were numerically conducted with the virtual femur used in the previous sections. A titanium stem similar to the Trilock©/Dual-Lock© of DePuy Ortopedics® was used (see figure 14).

The pos-operative geometry of the femur was discretized in 7224 hexahedral elements whose mesh is different of the original mesh (figure 15(a)). Therefore, a technique for transference of properties (or remeshing technique) should be employed. The technique employed is based in the distance between the centre of the new element and the centre of the old elements, and can be written as (see figure 15(b))


where V is the element volume and dist define the distance between the centres of elements. This technique is simple and effective. Figure 15(c) shows the transference of the relative density between the different meshes.

(a) (b) (c) (d)

Figure 14. A numerical surgery was performed in order to prepare the femur to receive a prosthesis. (a) Original femur; (b) femur after the removal of the femoral head and (b) the geometry of the titanium stem; (d) The coated and non-coated area of the titanium stem (E ≈ 115 GPa).


Figure 15. The use of different meshes requires the use of a remeshing technique. (a) Domain discretization of the femur before and after the numerical surgery. The original 9408 hexahedral elements properties must be converted to the new 7224 elements; (b) Relative density contours for the old and new domain discretization; (c) The remeshing technique scheme.

5.5. New Boundary Conditions

The boundary conditions affect significantly the results. A modification of the applied loads direction and intensity in the implanted femur must be considered. The hip reaction joint (HRJ) and muscles forces are modified due to geometric alterations created by each type of arthroplasty. These alterations are due to the change of prosthesis head position and to the changes in direction and intensity of the joint contact forces [34,35]. The different position of the stems depends not only on the geometry of stem but also on the technical procedure of the total hip arthroplasty. The prosthesis does not restore the original position of the femur head, which makes the bending moments different from the physiological original ones. Therefore, the loads and momentum equilibrium system must be re-established. From the work of Ramos [35], considering the muscle forces unchanged, it was showed that the conditions that must be respected are the force equilibrium in the z-axis and the moments equilibrium in the x and y rotation axis [35]. The solution of the force and moment equilibrium system of equations is listed in table 1.

Other authors [7,9,10] also performed the correction of the muscles force when carrying out in vitro studies of different behaviour stems.

5.6. Osteointegration and Bone Remodelling Results

For the remodelling and bone ingrowth test, a tangential relative displacement threshold of CSLIPLimit = 50 μm is used. The conditions to maintain full osteointegration are CPressLimit = 0.8 MPa and CShearLimit = 35 MPa. The initial condition for the coated interface is set to contact with friction (coefficient of 0.3). The non-coated surfaces are considered as frictionless contact surfaces.

The osteointegration results (bone ingrowth), after 100 time increments, can be seen in figure 16. The blue regions indicate where the bone ingrowth does not occur - there is a separation or a relative tangential displacement greater than 50 μm. The bone ingrowth and full osteointegration regions are in red. The green and orange regions indicate regions of bone ingrowth but with no full integration.

Table 2 summarizes the full osteointegration results. It should be noticed that there are zones with and without bone ingrowth on the coated zones. The results show higher bone osteointegration at regions 1 and 2. This can be explained by the direction of the applied load cases. However, for the time considered, bone integration was not achieved in half of the coated zone. These considerations corroborate clinical observations [7].

(a) (b) (c) (d) (e)

Figure 16. Osteointegration results after 100 time increments: (a) Full model and scale; Details of (b) lateral left and bottom; (c) Bottom and lateral right; (d) Top and lateral right and (e) Top and lateral left coated surfaces of osteointegration.

Table 2. Surface areas of osteointegration and bone ingrowth.

Figure 17 shows the bone density distribution for the implanted proximal femur after 100 time increments. Although the overall results show the loss of bone on a proximal region of the femur, a final relative volume fraction of approximately 70% was obtained. This high volume fraction must be evaluated considering that a fraction of the original bone was removed in the numerical surgery. The fraction removed (femur head) has low mean values of apparent density, therefore the remained part have higher mean volume fraction than the original bone. The region that decreased the initial global mean volume fraction, such as the intramedular cavity, is now occupied by the stem.

Comparing figures 17 and 12, it is possible to observe the influence of the prosthesis in the density distribution. First, it is observed a global redistribution of density. In figure 12, the zone with greater density is the diaphyseal zone, decreasing to the metaphyseal zone. Although the thickness of the cortical is thinner, in the implanted femur, there is no obvious decreasing of the volume fraction from the diaphyseal to the metaphyseal region near the region 4 of the implant. Relatively high bone density is seen in the region of front contact of the stem. The region where the head of the femur was ressected also shows high bone density. The density evolution in this region is due to the proximity to the muscle forces (the major trocanter) and to the stem. Therefore, this region is subjected to high stresses. Some of the previous observations and the results of figure 17 are similar to the results presented in [6,7] thought the shape and size of the stem are different.

Although a lack of clinical observations subsists, it can be stated that the prediction of the model is consistent with some clinical observations of the bone relative density distribution. The model presented in this work shows that osteointegration does not occur in all coated surface. This result is consistent with clinical observations [36-40] of 20% of fully osteointegration from partial coated recovered stems. The numerical results obtain by Viceconti et al. [12] to predict the secondary stability show higher


Figure 17. (a) Apparent density distribution of the implanted trabecular bone; (b) Transversals and (c) saggital longitudinal cuts of the femur’s model. The final volume fraction obtained was of approximately 70%.

values for the percentage of bonded surface, including partial or fully-bonded points. This difference can be attributed to the consideration of the bone remodelling scheme by this work, to the difference in the shape and size of the prosthesis and the coated zone. The numerical test conditions, such as the values for micromotions thresholds, were also considered different.


The remodelling scheme presented in this work is time dependent and reversible. The bone-remodelling scheme is numerically described by a time-dependent evolutionary procedure with anisotropic material parameters. The remodelling rate equation is obtained from the structural optimization task of maximizing the stiffness subject to a biological cost associated with metabolic maintenance of bone tissue in time. Results of bone adaptation process were presented. In this process, when the stimuli are increased, there is a production of bone mass. When they are reduced, the bone mass is absorbed. The model presented is in agreement with the generality of osteobiologic phenomenology. The time dependence is defined by the maximum density variation of bone. The maximum value of volume variation found by [28] was used in this work.

In this work, a computational model for bone remodelling around cementless (press-fit) stems with in growth control was also presented. The bone-prosthesis interface is modelled using a penalty contact formulation and the bone ingrowth regions are predicted using relative displacements and stresses criteria. During the process all the bone density, contact conditions and criteria verifications are updated. Therefore, the interface stiffness can evolve to reflect the bone ingrowth phenomenon. The interface stiffness increases with the osteointegration. When reaching full osteointegration, a complete glue interface connection is achieved.

The model presented was applied to a three-dimension finite element model of an implanted femur whose initial conditions, such as the bone density, were retrieved from the previous intact bone remodelling process by means of a remeshing technique. The information about the distribution of bone ingrowth and bone density can be given by the results. The bone-remodelling scheme shows that the introduction of a femoral prosthesis compels to a new bone density distribution. Although the distribution of bone densities obtained shows small loss of bone in proximal femur, it was observed an increasing of the bone density in the neighbourhood of the prosthesis.

Although a lack of clinical observations subsists, the predictions of the model agree with some qualitative clinical observations of the bone relative density distribution and agree with computational results presented by other authors.


The authors acknowledge the financial support by FCT—Fundação para a Ciência e Tecnologia—by the project SFRH/PTDC/EME-PME/ 68975/2006. The authors also thank Professor J. Folgado from IDMEC-IST for the help in some numerical details. The article was partially presented at the International Conference on Engineering Optimization, Rio de Janeiro, Brazil, 2008.


  1. Bensøe, M.P. and Sigmund, O. (2003) Topology optimization—Theory, methods and applications. SpringerVerlag, BH.
  2. Huiskes, R. (2000) If bone is the answer, then what is the question? Journal of Anatomy, 197, 145-156. doi:10.1046/j.1469-7580.2000.19720145.x
  3. Bagge, M. (2000) A model of bone adaptation as an optimization process. Journal of Biomechanics, 33, 1349- 1357. doi:10.1016/S0021-9290(00)00124-X
  4. Fernandes, P., Rodrigues, H. and Jacobs, C. (1999) A model of bone adaptation using a global optimisation criterion based on the trajectorial theory of wolff. Computer Methods in Biomechanics and Biomedical Engineering, 2, 125-138. doi:10.1080/10255849908907982
  5. Fernandes, P., Guedes, J.M. and Rodrigues, H. (1999) Topology optimization of 3D linear elastic structures with a constraint on perimeter. Computers & Structures, 73, 583-594. doi:10.1016/S0045-7949(98)00312-5
  6. Fernandes, P.R., Folgado, J., Jacobs, C. and Pellegrini, V. (2002) A contact model with ingrowth control for bone remodelling around cementless stems. Journal of Biomechanics, 35, 167-176. doi:10.1016/S0021-9290(01)00204-4
  7. Folgado, J., Fernandes, P.R., Guedes, J.M. and Rodrigues, H. (2004) Evaluation of osteoporotic bone quality by a computational model for bone remodelling. Computers & Structures, 82, 1381-1388. doi:10.1016/j.compstruc.2004.03.033
  8. Stolk, J., Verdonschot, N. and Huiskes, R. (2001) Hip joint and abductor-muscle forces adequately represent in vivo loading of a cemented total hip reconstruction. Journal of Biomechanics, 34, 917-926. doi:10.1016/S0021-9290(00)00225-6
  9. Waide, V., Cristofolini, L., Stolk, J., Verdonschot, N. and Toni, A. (2003) Experimental investigation of bone remodeling using composite femurs. Clinical Biomechanics, 18, 523-536. doi:10.1016/S0268-0033(03)00072-X
  10. Viceconti, M., Cavalloti, G., Andrisano, A.O. and Toni, A. (1995) Discussion on the design of a hip joint simulator. Medical Engineering & Physics, 18, 234-240. doi:10.1016/1350-4533(95)00026-7
  11. Viceconti, M., Muccini, R., Bernakiewicz, M., Baleani, M. and Cristofolini, L. (2000) Large-sliding contact elements accurately predict levels of bone-implant micromotion relevant to osseointegration. Journal of Biomechanics, 33, 1611-1618. doi:10.1016/S0021-9290(00)00140-8
  12. Viceconti, M., Pancanti, A., Dotti, M., Traina, F., Cristofolini, L. (2004) Effect of the initial implant fitting on the predicted secondary stability of a cementless stem. Medical and Biological Engineering and Computing, 42, 222-229. doi:10.1007/BF02344635
  13. Viceconti, M., Ricci, S., Pancanti, A. and Cappello, A. (2004b) Numerical model to predict the long-term mechanical stability of cementless orthopaedic implants. Medical and Biological Engineering and Computing, 42, 747-753. doi:10.1007/BF02345207
  14. Moreo, P., Péreza, M.A., García-Aznar, J.M. and Doblaré, M. (2007) Modelling the mechanical behaviour of living bony interfaces. Computer Methods in Applied Mechanics and Engineering, 196, 3300-3314. doi:10.1016/j.cma.2007.03.020
  15. Reina, E., Gómez-Benito, M.J., Garcia-Aznar, J.M., Dominguez, J. and Doblaré, M. (2007) Modelo computacional de distracción ósea. Proceeding of CEMNI/ CILAMCE Conference, Sociedad Española de Métodos Numéricos en Ingeniería (SEMNI), CD-ROM paper 132, 1-10.
  16. Reina, E., Gómez-Benito, M.J., Garcia-Aznar, J.M., Dominguez, J. and Doblaré, M. (2008) Modeling distraction osteogenesis: Analysis of the distraction rate. Biomechanics and Modeling in Mechanobiology, 8, 323-335.
  17. Coelho, P., Fernandes, P.R., Guedes, J.M. and Rodrigues, H.C. (2008) A hierarchical model for concurrent material and topology optimization of three-dimensional structures. Structural and Multidisciplinary Optimization, 35, 107- 115. doi:10.1007/s00158-007-0141-3
  18. Coelho, P.G., Fernandes, P.R., Rodrigues, H.C., Cardoso, J.B. and Guedes, J.M. (2009) Numerical modeling of bone tissue adaptation—A hierarchical approach for bone apparent density and trabecular structure. Journal of Biomechanics, 42, 830-837. doi:10.1016/j.jbiomech.2009.01.020
  19. Hardy, D.C., Frayssinet, P., Guilhem, A., Lafontaine, M.A. and Delince, P.E. (1991) Bonding of hydroxyapatite-coated femoral prostheses. Histopathology of specimens from four cases. Journal of Bone & Joint Surgery, 73, 732-740.
  20. Tonino, A.J., Therin, M. and Doyle, C. (1999) Hydroxyapatitecoated femoral stems. Histology and histomorphometry around five components retrieved at post mortem. Journal of Bone & Joint Surgery (British Volume), 81, 148-154. doi:10.1302/0301-620X.81B1.8948
  21. Fernandes, P. (1998) Optimização da Topologia de Estruturas Tridimensionais. PhD Thesis, Instituto Superior Técnico, UTL.
  22. Oliveira, J.A., Pinho-da-Cruz, J.A., Andrade-Campos, A. and Teixeira-Dias, F. (2006) Modelling the elastic behaviour of composites materials with asymptotic expansion homogenisation. Proceeding of COMPTEST, 151- 152.
  23. Pinho-da-Cruz, J., Oliveira, J.A., Teixeira-Dias, F. (2009) Asymptotic homogenisation in linear elasticity. Part I: Mathematical formulation and finite element modelling. Computational Materials Science, 45, 1073-1080. doi:10.1016/j.commatsci.2009.02.025
  24. Oliveira, J.A., Pinho-da-Cruz, J. and Teixeira-Dias, F. (2009) Asymptotic homogenisation in linear elasticity. Part II: Finite element procedures and multiscale applications. Computational Materials Science, 45, 1081-1096. doi:10.1016/j.commatsci.2009.01.027
  25. ABAQUS version 6.7 user manual (2006) Hibbitt. Karlsson & Sorenson. Inc.
  26. MSC.MARC (2005) Theory and User Information and User subroutines and special Routines. Msc. Software Corporation.
  27. Svanberg, K. (1987) The method of moving asymptotes – A new method for structural optimization. International Journal for Numerical Methods in Engineering, 24, 358-373. doi:10.1002/nme.1620240207
  28. Jee, W.S.S (1953) The skeletal tissues, histology: Cell & tissue biology. 5th Edition, Elsevier, Amsterdam.
  29. Medical Multimedia group (1999) A Patient’s Guide to Artificial Hip Replacement. Accessed 15 Jan 2008
  30. Encyclo, online encyclopedia (2009), encyclo MMIX
  31. Desmorat, B. (2007) Structural rigidity optimization with frictionless unilateral contact. International Journal of Solids and Structures, 44, 1132-1144. doi:10.1016/j.ijsolstr.2006.06.010
  32. Fancello, E.A. (2006) Topology optimization for minimum mass design considering local failure constraints and contact boundary conditions. Structural and Multidisciplinary Optimization, 32, 229-240. doi:10.1007/s00158-006-0019-9
  33. Folgado, J. (2004) Modelos computacionais para análise e projecto e próteses ortopédicas. PhD Thesis, IST-UTL.
  34. Ramos, A. (2006) Estudo numérico e experimental de uma nova componente Femural da prótese de Anca Cimentada. PhD Thesis, Universidade de Aveiro.
  35. Ramos, A., Fonseca, F. and Simões, J.A. (2006) Simulation of physiological loading in total hip replacements. Journal of Biomechanical Engineering, 128, 579-588. doi:10.1115/1.2205864
  36. Ramos, A. and Simões, J.A. (2006) Tetrahedral versus hexahedral finite elements in numerical modelling of the proximal femur. Medical Engineering & Physics, 28, 916-924. doi:10.1016/j.medengphy.2005.12.006
  37. Egan, J.M. and Marsden, D.C. (2001) A spring network model for the analysis of load transfer and tissue reactions in intramedullary fixation. Clinical Biomechanics, 16, 71-79. doi:10.1016/S0268-0033(00)00070-X
  38. Sadegh, A.M., Luo, G.M. and Cowin, S.C. (1993) Bone ingrowth: An application of the boundary element method to bone remodeling at the implant interface. Journal of Biomechanics, 26, 167-182. doi:10.1016/0021-9290(93)90047-I
  39. Spears, I.R., Pfleiderer, M., Schneider, E., Hille, E., Bergmannn, G. and Morlock, M.M. (2000) Interfacial conditions between a press-fit acetabulax cup and bone during daily activities: Implications for achieving bone in-growth. Journal of Biomechanics, 33, 1471-1477. doi:10.1016/S0021-9290(00)00096-8
  40. Keaveny, T. and Bartel, D. (1993) Mechanical consequences of bone ingrowth in a hip prosthesis inserted without cement. Journal of Bone & Joint Surgery, 77, 911-923.


1The prefix osteo, meaning bone, came from the Greek osteo. This term, used as an English prefix, is more accurate than the Latin term osseo.