- Email: [email protected]
Investigation of cancellous bone architecture using structural optimisation
ARTICLE IN PRESS
Journal of Biomechanics 41 (2008) 629–635 www.elsevier.com/locate/jbiomech www.JBiomech.com
Investigation of cancellous bone architecture using structural optimisation Hyunsun A. Kim, Paul J. Clement, James L. Cunningham Department of Mechanical Engineering, University of Bath, Bath BA2 7AY, UK Accepted 28 September 2007
Abstract The optimality of bone’s internal architecture and its remodelling is investigated by applying topological optimisation. The os-calcis (heel bone) which is subject to a relatively simple loading environment was used as the test problem. The solution is compared with X-ray and CT images of the os-calcis and is demonstrated to agree favourably with the bone’s typical trabecular orientation in the coronal plane. The optimality of bone remodelling is further investigated by applying a perturbation during optimisation which leads to a different optimum topology. This is compared with in-vivo experimental results from a model in which the same perturbation was applied. The findings indicate that there may exist more than one optimum state for the same mechanical condition and that structural optimisation is able to find these multiple states. It is therefore suggested that structural optimisation can be used to investigate the mechanisms of bone remodelling. r 2007 Elsevier Ltd. All rights reserved. Keywords: Bone remodelling; Bone adaptation; Topology optimisation
1. Introduction The form of adult bone is regulated by two main factors. Firstly, the predetermined genetic template that explains individual differences in both the form of individual bones and the general size of bones between individuals. Secondly, the fact that bone as a tissue is able to react to changes in the level of mechanical loading and will respond to both the magnitude and pattern of loading. Bone mass and architecture are regulated by physical activity, and significant changes in bone mass have been observed when physical activity has been reduced or increased (Lanyon, 1987). In cortical bone, the remodelling response is seen as a thickening or thinning of the area of the cortex in which the strains are above or below the optimal level (Lanyon et al., 1982). In cancellous bone, the resorptive remodelling response observed in response to disuse is a thinning of existing trabeculae and a loss of trabecular bone volume
Corresponding author. Tel.: +44 1225 383374; fax: +44 1225 386928.
E-mail address: [email protected] (H.A. Kim). 0021-9290/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2007.09.036
with eventual loss of connectivity between trabeculae (Thomas et al., 1996). The well-defined internal architecture of trabecular bone has been thought to represent an optimum biological response to its mechanical environment. What is not clearly established, however, is to what problem it is an optimal response. In this context, the understanding of the internal architecture of bone is an inverse problem of optimisation, and therefore structural optimisation methods can be utilised to investigate the complicated bone structures and the remodelling mechanism. The structural optimum is usually defined as the minimum compliance design hence the stiffest design, and therefore the local sensitivity becomes the strain energy density. Indeed, many current bone modelling hypotheses often employ strain energy density as the stimulus (Huiskes et al., 2000). Topological optimisation is the most general form of structural optimisation where the design variable is the existence of material, which is binary and discrete. A popular approach to this discrete optimisation problem is to relax the discrete variable to a continuous one analogous to material density. A finite element (FE)
ARTICLE IN PRESS 630
H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635
method is commonly employed to compute design sensitivity and the piecewise constant elemental density is updated according to this local sensitivity value. In order to gain a discrete solution without the intermediate densities, elemental densities are usually penalised by a power-law approach. This is often referred to as Solid Isotropic Materials with Penalisation or SIMP. Whilst this power-law-based penalisation may not give a completely discrete solution, it is understood that an appropriate penalisation can indicate a general topology of a good design for the given mesh (Bensøe, 1995). SIMP has been investigated for the study of bone remodelling. Huiskes (2000) summarised such efforts as limited, as they tended to report contours of the bone density patterns which could not offer much understanding of the internal architecture. Prendergast and Weinans (1998) applied SIMP to the well-known MBB-like beam problem of a simply supported beam of aspect ratio 6 with a central downward load. Their numerical experiments showed that applying a perturbation to a near-optimum solution leads to a sub-optimal solution different from the global optimum. This is consistent with the general understanding that the convergence of optimisation is dependent on the initial design and the fully populated uniform density domain is recommended to obtain the global optimum solution. They concluded, therefore, that structural optimisation would not offer an insight into the mechanism of bone remodelling in such cases. The aim of this paper is to further investigate the relevance of structural optimisation to the internal architecture of bone and its remodelling mechanism. We employ an alternative topological optimisation technique based on fully stressed design, instead of the minimum compliance design. The formulation and algorithm of the optimisation strategy, namely evolutionary structural optimisation (ESO), is outlined in the following section. We then apply it to a simple two-dimensional finite element model of an ovine (sheep) os-calcis and compare with the in-vivo architecture of the os-calcis in the sagittal plane, firstly with the normal trabecular architecture, secondly with the adapted bone response resulting from a perturbation to the normal trabecular structure. 2. Topology optimisation 2.1. Topology optimisation problem Topology optimisation is the most general form of structural optimisation in engineering, where the shape of the boundaries that defines a design is optimised as well as the number of boundaries, i.e. creating holes. It typically considers design problems in a continuum environment, such as Fig. 1, where the continuum domain is represented by O enclosed by a boundary, G, with body forces, f, the boundary traction, t on GtCG, and support on GuCG. The optimisation problem is formulated as a material distribution problem where the design variable, xi is either present
Γ
Ω
t xi Γt
f
Γu
Fig. 1. Generalised topology design problem in a two-dimensional domain.
or absent, xiA{1, 0}. The implementation typically employs the presence or absence of finite elements as a design variable. The most popular formulation of a structural optimisation problem is to minimise compliance, which in turn maximises the global stiffness, as written in a discrete form in Eq. (1). The constraints represent the finite element elastic equilibrium and the volume limit: min Cðxi ; ui Þ ¼ f T u subject to Kðxi Þu ¼ f
(1)
V pV max where C denotes total compliance, f and u represent the load and displacement vectors, respectively. The stiffness matrix K is dependent on the existence of material xi. V is the total volume of the design which is to be less than or equal to the maximum specified volume, Vmax (Bensøe, 1995). As the most generalised form of structural optimisation, topology optimisation is recognised to be the most versatile and least dependent on the initial design relative to sizing or shape optimisation. However, it is also considered to be the most difficult, as the problem is inherently discrete where the design variable xi is binary in nature. Discrete problems are usually non-convex and most of the existing methods are not applicable. Many efforts were seen in the last few decades to address this challenge and a number of methods have emerged. These methods can be categorised broadly into two categories. One is to relax the discrete design variables into continuous variables such that the existing methods are applicable to find the solution, then the continuous solution is discretised to give a binary solution with a uniform density or stiffness distribution which is more meaningful. These methods are considered to be more mathematically rigorous; however, a recent study showed that it is more sensitive to numerical techniques and parameters (Edwards et al., 2007). The other approach is to directly treat the discrete design variables based on engineering design practise or intuitive rules. These heuristic methods do not employ the more traditional mathematical programming techniques; however, they are often observed to be mathematically equivalent to the mathematical programming approaches and appear to be
ARTICLE IN PRESS H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635
numerically more stable relative to the continuous approaches (Edwards et al., 2007). 2.2. Optimisation method The method considered in this paper is ESO (Xie and Steven, 1997). This heuristic method has been shown to minimise a global parameter of the total compliancevolume product and is considered equivalent to sequential linear programming (Tanskanen, 2002). This study thus indicates that ESO may find a local optimum solution, similar to the mathematical programming approaches. It employs the finite element method to compute the efficiency of material, which is defined by elemental stress, and slowly removes inefficient material and adds material around high-sensitivity areas. A filtering method based on image processing techniques is applied to the sensitivities to prevent a chequerboard formation (Diaz and Sigmund, 1995). The ESO algorithm can thus be summarised as follows: (1) A fully populated design domain is specified as the initial design with loads and boundary conditions in the finite element environment. (2) Finite element analysis is carried out to determine the elemental von Mises stress of the current design. (3) Filtering is applied to von Mises stress to compute elemental sensitivities: P
wi si X i s¯ e ¼ P wi si , wi i
(2)
(6) The modified design is checked for convergence, i.e. the objective function value remains within a given small tolerance over a specified number of iterations, such as 5. If a convergence is reached, the optimisation algorithm is terminated. If not, it returns to step 2 with an increased MR (given by Eq. (6)) and continues to iterate: MR ¼ MR þ ER;
where element e represents the element under consideration and elements i represent all elements within a given radius, r. The weighting factor wi is based on the distance between element, e and the surrounding elements, i, i.e. dist(e, i):
where ER is the evolutionary rate, typically a small value such as 0.01, which controls the rate of modification. A drawback common to heuristic approaches is the typically high computational cost relative to the more mathematically rigorous methods. This is attributed to its zeroth order nature where the magnitudes of the sensitivities do not govern the convergence rate unlike the mathematical programming approaches, which in turn require higher order sensitivity computations. Due to this lack of dependency on high-order sensitivities, ESO is often perceived to be intuitive and easy to implement. Furthermore, a recent study has shown that ESO tends to be less sensitive to the numerical implementation and more independent of parameter settings and initial design domain, provided that a reasonable size mesh density is selected (Edwards et al., 2007). For this study, the algorithm was implemented using ANSYS Parametric Design Language with the ANSYS finite element linear static solver (ANSYS Inc., 2003).
3.1. Finite element model of ovine os-calcis The ovine os-calcis was selected as a test problem due to its relatively simple shape and loading (Su et al., 1999) and
(3)
(4) Elements that satisfy Eq. (4) are removed: se pMR RR smaximum ,
(4)
where se is the elemental von Mises stress; MR is the modification ratio, typically a small value such as 0.01, this ensures that only small modifications are made to the design at each iteration; RR is the rejection ratio, 0o RRp1.0; and smaximum is the maximum elemental von Mises stress. (5) Elements are added around element e when Eq. (5) is satisfied: se Xð1 MR ARÞ smaximum , where AR is the addition ratio, 0o ARp1.0.
(6)
3. Investigation of ovine os-calcis
i
wi ¼ r distðe; iÞ.
631
(5) Fig. 2. Ovine os-calcis.
ARTICLE IN PRESS 632
H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635
because detailed data was available from a previous study on its internal architecture both before and after a surgical insult. A typical 2D geometry of an ovine oscalcis in Fig. 2 was used. Based on a mesh convergence study, the 2D geometry was modelled with a structured mesh of 108 51 linear quadrilateral elements. Fig. 3 shows the finite element model where the dark area represents the elements outside the design domain and the white area represents the fully populated initial design. The structural configuration was assumed to be similar to a cantilever beam, where the boundary along the left-hand side was fixed both in both translation and rotation. The loading condition was simplified to a single load of bending induced by a vertical load applied through the Achilles tendon. The magnitude of the applied load (700 N) was estimated as being between two and three times the measured vertical ground reaction force during a normal gait cycle (Tapper et al., 2006) for a typical body weight of 70 kg. Isotropic material properties were specified, with Young’s modulus of 20 GPa and Poisson’s ratio of 0.3 (Guo, 2001).
3.2. Investigation of structural optimality Optimisation was applied by allowing only element removals with ER ¼ 0.001 and RR ¼ 0.01. At iteration 100, an optimum solution of Fig. 4 was obtained with the structural volume reduction of 30%. A typical cantilevered beam solution with an internal truss-like configuration was apparent, and a relatively hollow section near the fixed boundary was also observed. These features were similar to the internal architecture of the ovine os-calcis, an example of which can be seen in Fig. 2. The internal trabecular orientation of Fig. 4 was observed to be approximately 401, which compared closely with the inclination of an internal branch at a similar location in the optimum solution, Fig. 4, of 451. The single vertical applied load is a grossly simplified representation of the realistic loading applied to the oscalcis, and indeed the fine details were lost in the ESO solution as a result. In order to investigate the effects of the loading, the angle at which the load was applied was varied by 7201 from the vertical. The numerical tests showed that
Fig. 3. Finite element model of ovine os-calcis domain.
Fig. 4. Optimum solution for a vertical load.
ARTICLE IN PRESS H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635
633
this primarily influenced the size of the upper and the lower cortical surfaces but the general orientation of the internal grillage struts remained largely unaffected (Fig. 5). Therefore, it is evident that topological optimisation based on the localised stress response can offer an insight into the internal bone architecture. 3.3. Investigation of os-calcis with a drilled section
Fig. 5. Optimum solutions with load variations: (a) load at +201 and (b) load at 201.
Following on from obtaining the internal trabecular structure from the fully populated domain, a perturbation was applied by simulating the drilled-out centre section of the os-calcis at iteration 150, as shown in Fig. 6. Optimisation was then applied to this new structure with a hollowed centre section, where elements were allowed to be both added and removed, with MR ¼ 0.1, AR ¼ 0.6, RR ¼ 0.4 and ER ¼ 0.01. This led to a new optimum solution of Fig. 7. The internal grillage structure was not re-formed, instead the upper and lower beams were increased in their sizes. The average stress level in the solution with a drilled section, Fig. 7, was higher, which accounted for its lower volume. However, the stress distribution of the two
Fig. 6. Ovine os-calcis with drilled section.
Fig. 7. ESO optimum solution for os-calcis with drilled section.
ARTICLE IN PRESS 634
H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635
solutions, Figs. 4 and 7 exhibit a similar level of even distribution of stress which is the objective. As the optimum solutions are dependent on the initial design, the optimisation problems are not convex but have at least two local optima. As the optimisation method cannot guarantee to find the global solution for non-convex problems, it is seen that structural optimisation leads to a local optimum. However, this finding is, in fact, consistent with in-vivo experimental results of a study aimed at determining the remodelling response of the os-calcis to a significant loss of trabecular bone (Cunningham et al., 2005). The experiment involved the removal of a significant amount of the internal trabeculae of the right os-calcis. The animal was returned to normal husbandry, with unrestricted weight-bearing and the bone remodelling process was observed over a 15-month period. Fig. 8 shows the sagittal sections obtained from a micro-CT scan at 15 months postoperation for the both unoperated and operated os-calcis, parts (a) and (b), respectively. It was observed, as shown in Fig. 8(b), that the trabeculae had not regenerated, rather the os-calcis remained hollow with an increased density and thickness of the cortical bone. Such losses of trabeculae have also been observed in other cases, such as Li et al. (2003), where trabeculae were not restored in
older adult dogs, instead the thickness of the cortical bone was increased. Both these features are prominent in the structural optimum solution of Fig. 7, where the internal trabeculae did not reform but the thickness of the upper and lower beams have increased to redistribute the applied load, this being evident in the increased density and thickness of the cortical bone in Fig. 8(b). It appears that nature’s bone remodelling process may also be of local optimisation, rather than global optimisation. Whilst this needs further research to draw a meaningful conclusion, it can be seen that structural optimisation may be considered a useful investigative tool for bone remodelling. 4. Conclusion In this paper, a heuristic topological optimisation method, ESO, was applied to investigate the internal architecture of bone. The initial findings are consistent with the previous research results that the mechanical environment is a dominant factor in the trabecular configuration which tends to resemble the structurally optimum grillage structure. Also under the assumption of isotropy, stressdriven optimisation offers a good insight into the biological response to the loading, as an alternative to compliancebased formulation. When a perturbation was applied to an optimum solution by drilling out the centre section of an ovine oscalcis, optimisation found a different solution where the internal trabeculae were not reformed, but rather the density and thickness of the cortical bone were increased. This is because the optimisation algorithm is based on the local response and the solution is dependent on the initial design, and hence it appears that the biological response is also based on local response. The comparison with in-vivo experimental results revealed that this local optimum is indeed what was observed at 15 months post-operation. This suggests that there may exist more than one equilibrium state and remodelling may therefore not always aim for the global optimum but perhaps a local optimum. To this extent, structural optimisation may be considered a valid investigative tool for bone remodelling. Conflict of interest We declare that there are no conflicts of interest of any of the authors of this paper in the work contained within. References
Fig. 8. Os-calcis micro-CT images 15 months post-operation: (a) unoperated left leg: intact showing normal trabeculae and (b) operated right leg: internal trabeculae has not reformed.
ANSYS, Inc., 2003. ANSYS 7.1 On-line Manual. Bensøe, M.P., 1995. Optimization of Structural Topology, Shape and Material. Springer, Berlin. Cunningham, J.L., Smith, T.J., Blunn, G.W., Goodship, A.E., 2005. Following damage, trabecular architecture is not restored by functional loading. In: Transactions of the 51st Annual Meeting of the Orthopaedic Research Society.
ARTICLE IN PRESS H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635 Diaz, A.R., Sigmund, O., 1995. Checkerboard patterns in layout optimization. Structural Optimization 10, 40–45. Edwards, C.S., Kim, H.A., Budd, C.J., 2007. An evaluative study of ESO and SIMP for optimising a cantilever tie-beam. Structural and Multidisciplinary Optimization, On-Line First. /http://dx.doi.org/ 10.1007/s00158-007-0102-xS. Guo, X.E., 2001. Mechanical properties of cortical bone and cancellous bone tissue. In: Cowin, S.C. (Ed.), Bone Mechanics Handbook. CRC Press, Boca Raton, FL, pp. 10–23. Huiskes, R., 2000. If bone is the answer, then what is the question? Journal of Anatomy 197, 145–156. Huiskes, R., Rulmerman, R., van Lenthe, G.H., Janssen, J.D., 2000. Effects of mechanical forces on maintenance and adaptation of form in trabecular bone. Nature 405, 706–740. Lanyon, L.E., 1987. Functional strain in bone tissue as an objective and controlling stimulus for adaptive bone remodelling. Journal of Biomechanics 20, 1083–1093. Lanyon, L.E., Goodship, A.E., Pye, C.J., MacFie, J.H., 1982. Mechanically adaptive bone remodelling. Journal of Biomechanics 15, 141–154. Li, C.Y., Laudier, D., Schaffler, M.B., 2003. Remobilization restores cancellous bone mass but not microarchitecture after long term disuse
635
in older adult dogs. In: Transactions of the 48th Annual meeting of the Orthopaedic Research Society, New Orleans. Prendergast, P.J., Weinans, H, 1998. Tissue adaptation as a discretedynamical process in time and space. In: Pedersen, P., Bensøe, M.P., (Eds.), Proceedings of the IUTAM Symposium on Synthesis in Bio Solid Mechanics. Springer, Berlin, pp. 321–332. Su, S.C., Skedros, J.G., Bachus, K.N., Bloebaum, R.D., 1999. Loading conditions and cortical bone construction of an artiodactyl calcaneus. Journal of Experimental Biology 202, 3239–3254. Tanskanen, P., 2002. The evolutionary structural optimization method: theoretical aspects. Computer Methods in Applied Mechanics and Engineering 191, 5485–5498. Tapper, J.E., Fukushima, S., Azuma, H., Thornton, G.M., Ronsky, J.L., Shrive, N.G., Frank, C.B., 2006. Dynamic in vivo kinematics of the intact ovine stifle joint. Journal of Orthopaedic Research 24, 782–792. Thomas, T., Vico, L., Skerry, T.M., Caulin, F., Lanyon, L.E., Alexandre, C., Lafage, M.-H., 1996. Architectural modifications and cellular response during disuse-related bone loss in calcaneus of the sheep. American Journal of Physiology 80, 198–202. Xie, Y.M., Steven, G.P., 1997. Evolutionary Structural Optimisation. Springer, Berlin.
Journal of Biomechanics 41 (2008) 629–635 www.elsevier.com/locate/jbiomech www.JBiomech.com
Investigation of cancellous bone architecture using structural optimisation Hyunsun A. Kim, Paul J. Clement, James L. Cunningham Department of Mechanical Engineering, University of Bath, Bath BA2 7AY, UK Accepted 28 September 2007
Abstract The optimality of bone’s internal architecture and its remodelling is investigated by applying topological optimisation. The os-calcis (heel bone) which is subject to a relatively simple loading environment was used as the test problem. The solution is compared with X-ray and CT images of the os-calcis and is demonstrated to agree favourably with the bone’s typical trabecular orientation in the coronal plane. The optimality of bone remodelling is further investigated by applying a perturbation during optimisation which leads to a different optimum topology. This is compared with in-vivo experimental results from a model in which the same perturbation was applied. The findings indicate that there may exist more than one optimum state for the same mechanical condition and that structural optimisation is able to find these multiple states. It is therefore suggested that structural optimisation can be used to investigate the mechanisms of bone remodelling. r 2007 Elsevier Ltd. All rights reserved. Keywords: Bone remodelling; Bone adaptation; Topology optimisation
1. Introduction The form of adult bone is regulated by two main factors. Firstly, the predetermined genetic template that explains individual differences in both the form of individual bones and the general size of bones between individuals. Secondly, the fact that bone as a tissue is able to react to changes in the level of mechanical loading and will respond to both the magnitude and pattern of loading. Bone mass and architecture are regulated by physical activity, and significant changes in bone mass have been observed when physical activity has been reduced or increased (Lanyon, 1987). In cortical bone, the remodelling response is seen as a thickening or thinning of the area of the cortex in which the strains are above or below the optimal level (Lanyon et al., 1982). In cancellous bone, the resorptive remodelling response observed in response to disuse is a thinning of existing trabeculae and a loss of trabecular bone volume
Corresponding author. Tel.: +44 1225 383374; fax: +44 1225 386928.
E-mail address: [email protected] (H.A. Kim). 0021-9290/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2007.09.036
with eventual loss of connectivity between trabeculae (Thomas et al., 1996). The well-defined internal architecture of trabecular bone has been thought to represent an optimum biological response to its mechanical environment. What is not clearly established, however, is to what problem it is an optimal response. In this context, the understanding of the internal architecture of bone is an inverse problem of optimisation, and therefore structural optimisation methods can be utilised to investigate the complicated bone structures and the remodelling mechanism. The structural optimum is usually defined as the minimum compliance design hence the stiffest design, and therefore the local sensitivity becomes the strain energy density. Indeed, many current bone modelling hypotheses often employ strain energy density as the stimulus (Huiskes et al., 2000). Topological optimisation is the most general form of structural optimisation where the design variable is the existence of material, which is binary and discrete. A popular approach to this discrete optimisation problem is to relax the discrete variable to a continuous one analogous to material density. A finite element (FE)
ARTICLE IN PRESS 630
H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635
method is commonly employed to compute design sensitivity and the piecewise constant elemental density is updated according to this local sensitivity value. In order to gain a discrete solution without the intermediate densities, elemental densities are usually penalised by a power-law approach. This is often referred to as Solid Isotropic Materials with Penalisation or SIMP. Whilst this power-law-based penalisation may not give a completely discrete solution, it is understood that an appropriate penalisation can indicate a general topology of a good design for the given mesh (Bensøe, 1995). SIMP has been investigated for the study of bone remodelling. Huiskes (2000) summarised such efforts as limited, as they tended to report contours of the bone density patterns which could not offer much understanding of the internal architecture. Prendergast and Weinans (1998) applied SIMP to the well-known MBB-like beam problem of a simply supported beam of aspect ratio 6 with a central downward load. Their numerical experiments showed that applying a perturbation to a near-optimum solution leads to a sub-optimal solution different from the global optimum. This is consistent with the general understanding that the convergence of optimisation is dependent on the initial design and the fully populated uniform density domain is recommended to obtain the global optimum solution. They concluded, therefore, that structural optimisation would not offer an insight into the mechanism of bone remodelling in such cases. The aim of this paper is to further investigate the relevance of structural optimisation to the internal architecture of bone and its remodelling mechanism. We employ an alternative topological optimisation technique based on fully stressed design, instead of the minimum compliance design. The formulation and algorithm of the optimisation strategy, namely evolutionary structural optimisation (ESO), is outlined in the following section. We then apply it to a simple two-dimensional finite element model of an ovine (sheep) os-calcis and compare with the in-vivo architecture of the os-calcis in the sagittal plane, firstly with the normal trabecular architecture, secondly with the adapted bone response resulting from a perturbation to the normal trabecular structure. 2. Topology optimisation 2.1. Topology optimisation problem Topology optimisation is the most general form of structural optimisation in engineering, where the shape of the boundaries that defines a design is optimised as well as the number of boundaries, i.e. creating holes. It typically considers design problems in a continuum environment, such as Fig. 1, where the continuum domain is represented by O enclosed by a boundary, G, with body forces, f, the boundary traction, t on GtCG, and support on GuCG. The optimisation problem is formulated as a material distribution problem where the design variable, xi is either present
Γ
Ω
t xi Γt
f
Γu
Fig. 1. Generalised topology design problem in a two-dimensional domain.
or absent, xiA{1, 0}. The implementation typically employs the presence or absence of finite elements as a design variable. The most popular formulation of a structural optimisation problem is to minimise compliance, which in turn maximises the global stiffness, as written in a discrete form in Eq. (1). The constraints represent the finite element elastic equilibrium and the volume limit: min Cðxi ; ui Þ ¼ f T u subject to Kðxi Þu ¼ f
(1)
V pV max where C denotes total compliance, f and u represent the load and displacement vectors, respectively. The stiffness matrix K is dependent on the existence of material xi. V is the total volume of the design which is to be less than or equal to the maximum specified volume, Vmax (Bensøe, 1995). As the most generalised form of structural optimisation, topology optimisation is recognised to be the most versatile and least dependent on the initial design relative to sizing or shape optimisation. However, it is also considered to be the most difficult, as the problem is inherently discrete where the design variable xi is binary in nature. Discrete problems are usually non-convex and most of the existing methods are not applicable. Many efforts were seen in the last few decades to address this challenge and a number of methods have emerged. These methods can be categorised broadly into two categories. One is to relax the discrete design variables into continuous variables such that the existing methods are applicable to find the solution, then the continuous solution is discretised to give a binary solution with a uniform density or stiffness distribution which is more meaningful. These methods are considered to be more mathematically rigorous; however, a recent study showed that it is more sensitive to numerical techniques and parameters (Edwards et al., 2007). The other approach is to directly treat the discrete design variables based on engineering design practise or intuitive rules. These heuristic methods do not employ the more traditional mathematical programming techniques; however, they are often observed to be mathematically equivalent to the mathematical programming approaches and appear to be
ARTICLE IN PRESS H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635
numerically more stable relative to the continuous approaches (Edwards et al., 2007). 2.2. Optimisation method The method considered in this paper is ESO (Xie and Steven, 1997). This heuristic method has been shown to minimise a global parameter of the total compliancevolume product and is considered equivalent to sequential linear programming (Tanskanen, 2002). This study thus indicates that ESO may find a local optimum solution, similar to the mathematical programming approaches. It employs the finite element method to compute the efficiency of material, which is defined by elemental stress, and slowly removes inefficient material and adds material around high-sensitivity areas. A filtering method based on image processing techniques is applied to the sensitivities to prevent a chequerboard formation (Diaz and Sigmund, 1995). The ESO algorithm can thus be summarised as follows: (1) A fully populated design domain is specified as the initial design with loads and boundary conditions in the finite element environment. (2) Finite element analysis is carried out to determine the elemental von Mises stress of the current design. (3) Filtering is applied to von Mises stress to compute elemental sensitivities: P
wi si X i s¯ e ¼ P wi si , wi i
(2)
(6) The modified design is checked for convergence, i.e. the objective function value remains within a given small tolerance over a specified number of iterations, such as 5. If a convergence is reached, the optimisation algorithm is terminated. If not, it returns to step 2 with an increased MR (given by Eq. (6)) and continues to iterate: MR ¼ MR þ ER;
where element e represents the element under consideration and elements i represent all elements within a given radius, r. The weighting factor wi is based on the distance between element, e and the surrounding elements, i, i.e. dist(e, i):
where ER is the evolutionary rate, typically a small value such as 0.01, which controls the rate of modification. A drawback common to heuristic approaches is the typically high computational cost relative to the more mathematically rigorous methods. This is attributed to its zeroth order nature where the magnitudes of the sensitivities do not govern the convergence rate unlike the mathematical programming approaches, which in turn require higher order sensitivity computations. Due to this lack of dependency on high-order sensitivities, ESO is often perceived to be intuitive and easy to implement. Furthermore, a recent study has shown that ESO tends to be less sensitive to the numerical implementation and more independent of parameter settings and initial design domain, provided that a reasonable size mesh density is selected (Edwards et al., 2007). For this study, the algorithm was implemented using ANSYS Parametric Design Language with the ANSYS finite element linear static solver (ANSYS Inc., 2003).
3.1. Finite element model of ovine os-calcis The ovine os-calcis was selected as a test problem due to its relatively simple shape and loading (Su et al., 1999) and
(3)
(4) Elements that satisfy Eq. (4) are removed: se pMR RR smaximum ,
(4)
where se is the elemental von Mises stress; MR is the modification ratio, typically a small value such as 0.01, this ensures that only small modifications are made to the design at each iteration; RR is the rejection ratio, 0o RRp1.0; and smaximum is the maximum elemental von Mises stress. (5) Elements are added around element e when Eq. (5) is satisfied: se Xð1 MR ARÞ smaximum , where AR is the addition ratio, 0o ARp1.0.
(6)
3. Investigation of ovine os-calcis
i
wi ¼ r distðe; iÞ.
631
(5) Fig. 2. Ovine os-calcis.
ARTICLE IN PRESS 632
H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635
because detailed data was available from a previous study on its internal architecture both before and after a surgical insult. A typical 2D geometry of an ovine oscalcis in Fig. 2 was used. Based on a mesh convergence study, the 2D geometry was modelled with a structured mesh of 108 51 linear quadrilateral elements. Fig. 3 shows the finite element model where the dark area represents the elements outside the design domain and the white area represents the fully populated initial design. The structural configuration was assumed to be similar to a cantilever beam, where the boundary along the left-hand side was fixed both in both translation and rotation. The loading condition was simplified to a single load of bending induced by a vertical load applied through the Achilles tendon. The magnitude of the applied load (700 N) was estimated as being between two and three times the measured vertical ground reaction force during a normal gait cycle (Tapper et al., 2006) for a typical body weight of 70 kg. Isotropic material properties were specified, with Young’s modulus of 20 GPa and Poisson’s ratio of 0.3 (Guo, 2001).
3.2. Investigation of structural optimality Optimisation was applied by allowing only element removals with ER ¼ 0.001 and RR ¼ 0.01. At iteration 100, an optimum solution of Fig. 4 was obtained with the structural volume reduction of 30%. A typical cantilevered beam solution with an internal truss-like configuration was apparent, and a relatively hollow section near the fixed boundary was also observed. These features were similar to the internal architecture of the ovine os-calcis, an example of which can be seen in Fig. 2. The internal trabecular orientation of Fig. 4 was observed to be approximately 401, which compared closely with the inclination of an internal branch at a similar location in the optimum solution, Fig. 4, of 451. The single vertical applied load is a grossly simplified representation of the realistic loading applied to the oscalcis, and indeed the fine details were lost in the ESO solution as a result. In order to investigate the effects of the loading, the angle at which the load was applied was varied by 7201 from the vertical. The numerical tests showed that
Fig. 3. Finite element model of ovine os-calcis domain.
Fig. 4. Optimum solution for a vertical load.
ARTICLE IN PRESS H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635
633
this primarily influenced the size of the upper and the lower cortical surfaces but the general orientation of the internal grillage struts remained largely unaffected (Fig. 5). Therefore, it is evident that topological optimisation based on the localised stress response can offer an insight into the internal bone architecture. 3.3. Investigation of os-calcis with a drilled section
Fig. 5. Optimum solutions with load variations: (a) load at +201 and (b) load at 201.
Following on from obtaining the internal trabecular structure from the fully populated domain, a perturbation was applied by simulating the drilled-out centre section of the os-calcis at iteration 150, as shown in Fig. 6. Optimisation was then applied to this new structure with a hollowed centre section, where elements were allowed to be both added and removed, with MR ¼ 0.1, AR ¼ 0.6, RR ¼ 0.4 and ER ¼ 0.01. This led to a new optimum solution of Fig. 7. The internal grillage structure was not re-formed, instead the upper and lower beams were increased in their sizes. The average stress level in the solution with a drilled section, Fig. 7, was higher, which accounted for its lower volume. However, the stress distribution of the two
Fig. 6. Ovine os-calcis with drilled section.
Fig. 7. ESO optimum solution for os-calcis with drilled section.
ARTICLE IN PRESS 634
H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635
solutions, Figs. 4 and 7 exhibit a similar level of even distribution of stress which is the objective. As the optimum solutions are dependent on the initial design, the optimisation problems are not convex but have at least two local optima. As the optimisation method cannot guarantee to find the global solution for non-convex problems, it is seen that structural optimisation leads to a local optimum. However, this finding is, in fact, consistent with in-vivo experimental results of a study aimed at determining the remodelling response of the os-calcis to a significant loss of trabecular bone (Cunningham et al., 2005). The experiment involved the removal of a significant amount of the internal trabeculae of the right os-calcis. The animal was returned to normal husbandry, with unrestricted weight-bearing and the bone remodelling process was observed over a 15-month period. Fig. 8 shows the sagittal sections obtained from a micro-CT scan at 15 months postoperation for the both unoperated and operated os-calcis, parts (a) and (b), respectively. It was observed, as shown in Fig. 8(b), that the trabeculae had not regenerated, rather the os-calcis remained hollow with an increased density and thickness of the cortical bone. Such losses of trabeculae have also been observed in other cases, such as Li et al. (2003), where trabeculae were not restored in
older adult dogs, instead the thickness of the cortical bone was increased. Both these features are prominent in the structural optimum solution of Fig. 7, where the internal trabeculae did not reform but the thickness of the upper and lower beams have increased to redistribute the applied load, this being evident in the increased density and thickness of the cortical bone in Fig. 8(b). It appears that nature’s bone remodelling process may also be of local optimisation, rather than global optimisation. Whilst this needs further research to draw a meaningful conclusion, it can be seen that structural optimisation may be considered a useful investigative tool for bone remodelling. 4. Conclusion In this paper, a heuristic topological optimisation method, ESO, was applied to investigate the internal architecture of bone. The initial findings are consistent with the previous research results that the mechanical environment is a dominant factor in the trabecular configuration which tends to resemble the structurally optimum grillage structure. Also under the assumption of isotropy, stressdriven optimisation offers a good insight into the biological response to the loading, as an alternative to compliancebased formulation. When a perturbation was applied to an optimum solution by drilling out the centre section of an ovine oscalcis, optimisation found a different solution where the internal trabeculae were not reformed, but rather the density and thickness of the cortical bone were increased. This is because the optimisation algorithm is based on the local response and the solution is dependent on the initial design, and hence it appears that the biological response is also based on local response. The comparison with in-vivo experimental results revealed that this local optimum is indeed what was observed at 15 months post-operation. This suggests that there may exist more than one equilibrium state and remodelling may therefore not always aim for the global optimum but perhaps a local optimum. To this extent, structural optimisation may be considered a valid investigative tool for bone remodelling. Conflict of interest We declare that there are no conflicts of interest of any of the authors of this paper in the work contained within. References
Fig. 8. Os-calcis micro-CT images 15 months post-operation: (a) unoperated left leg: intact showing normal trabeculae and (b) operated right leg: internal trabeculae has not reformed.
ANSYS, Inc., 2003. ANSYS 7.1 On-line Manual. Bensøe, M.P., 1995. Optimization of Structural Topology, Shape and Material. Springer, Berlin. Cunningham, J.L., Smith, T.J., Blunn, G.W., Goodship, A.E., 2005. Following damage, trabecular architecture is not restored by functional loading. In: Transactions of the 51st Annual Meeting of the Orthopaedic Research Society.
ARTICLE IN PRESS H.A. Kim et al. / Journal of Biomechanics 41 (2008) 629–635 Diaz, A.R., Sigmund, O., 1995. Checkerboard patterns in layout optimization. Structural Optimization 10, 40–45. Edwards, C.S., Kim, H.A., Budd, C.J., 2007. An evaluative study of ESO and SIMP for optimising a cantilever tie-beam. Structural and Multidisciplinary Optimization, On-Line First. /http://dx.doi.org/ 10.1007/s00158-007-0102-xS. Guo, X.E., 2001. Mechanical properties of cortical bone and cancellous bone tissue. In: Cowin, S.C. (Ed.), Bone Mechanics Handbook. CRC Press, Boca Raton, FL, pp. 10–23. Huiskes, R., 2000. If bone is the answer, then what is the question? Journal of Anatomy 197, 145–156. Huiskes, R., Rulmerman, R., van Lenthe, G.H., Janssen, J.D., 2000. Effects of mechanical forces on maintenance and adaptation of form in trabecular bone. Nature 405, 706–740. Lanyon, L.E., 1987. Functional strain in bone tissue as an objective and controlling stimulus for adaptive bone remodelling. Journal of Biomechanics 20, 1083–1093. Lanyon, L.E., Goodship, A.E., Pye, C.J., MacFie, J.H., 1982. Mechanically adaptive bone remodelling. Journal of Biomechanics 15, 141–154. Li, C.Y., Laudier, D., Schaffler, M.B., 2003. Remobilization restores cancellous bone mass but not microarchitecture after long term disuse
635
in older adult dogs. In: Transactions of the 48th Annual meeting of the Orthopaedic Research Society, New Orleans. Prendergast, P.J., Weinans, H, 1998. Tissue adaptation as a discretedynamical process in time and space. In: Pedersen, P., Bensøe, M.P., (Eds.), Proceedings of the IUTAM Symposium on Synthesis in Bio Solid Mechanics. Springer, Berlin, pp. 321–332. Su, S.C., Skedros, J.G., Bachus, K.N., Bloebaum, R.D., 1999. Loading conditions and cortical bone construction of an artiodactyl calcaneus. Journal of Experimental Biology 202, 3239–3254. Tanskanen, P., 2002. The evolutionary structural optimization method: theoretical aspects. Computer Methods in Applied Mechanics and Engineering 191, 5485–5498. Tapper, J.E., Fukushima, S., Azuma, H., Thornton, G.M., Ronsky, J.L., Shrive, N.G., Frank, C.B., 2006. Dynamic in vivo kinematics of the intact ovine stifle joint. Journal of Orthopaedic Research 24, 782–792. Thomas, T., Vico, L., Skerry, T.M., Caulin, F., Lanyon, L.E., Alexandre, C., Lafage, M.-H., 1996. Architectural modifications and cellular response during disuse-related bone loss in calcaneus of the sheep. American Journal of Physiology 80, 198–202. Xie, Y.M., Steven, G.P., 1997. Evolutionary Structural Optimisation. Springer, Berlin.