Meccanica (2007) 42:85–93 DOI 10.1007/s11012-006-9024-8
Numerical and experimental structural analysis of trabecular architectures F. Cosmi · D. Dreossi
Received: 25 March 2006 / Accepted: 19 July 2006 / Published online: 6 December 2006 © Springer Science+Business Media B.V. 2006
Abstract This paper presents the application of the Cell Method (CM) to the static analysis of 3D structures obtained from micro-computed tomography reconstructions of trabecular bone. The CM is a recently introduced numerical method, based on a direct discrete formulation of equilibrium equations, which is particularly promising for the analysis of complex structures. In fact, due to the direct discrete approach employed, no restriction is imposed by differentiability conditions and the characteristic length of the elementary cell of the discretization can be of the same order of magnitude as the heterogeneities of the structure. The same 3D microstructures used for the numerical simulations were reproduced by means of a rapid prototyping process, by selective laser sintering of polyamide powder. The compression elastic modulus of the replicas was experimentally determined and used as parameter for comparison with the simulations results. The experimental values are in good agreement with the numerical ones, thus validating the methodology employed. Keywords Cell method · Micro-CT · Rapid prototyping · Structural analysis · Trabecular bone F. Cosmi (B) · D. Dreossi Università di Trieste, Dipartimento di Ingegneria Meccanica, via A. Valerio 10, 34127 Trieste, Italy e-mail:
[email protected]
1 Introduction The macro-scale mechanical properties of several materials encountered in both biological tissues and industrial applications are a consequence of their complex micro-scale structure. A typical and important example is the cancellous bone which consists of an irregular lattice work of thin plates called trabeculae; within the trabeculae there are small spaces called lacunae that contain red bone marrow. For example, a digital radiography of the trabecular structure at the human metacarpus base is visible in Fig. 1. It is widely accepted that the mechanical strength of cancellous bone is related not only to the mineral content, but also to the trabecular micro-architecture arrangement [1–3]. There is a clinical relevance of the factors determining the mechanical properties of cancellous bone structure on osteoporosis, osteoarthritis diagnosis and orthopaedic implants. Several authors have addressed the problem of predicting the elastic properties of the bone mass architecture, since apparent density or mineral volume fraction are not sufficient for the complete characterization of a porous material of this structural complexity. A second rank symmetric tensor, the fabric tensor, was introduced by Cowin [4] as a measure of microstructure arrangement in a porous solid, assuming that any anisotropy exhibited is due only
86
Fig. 1 Digital radiography of the cancellous bone at the human metacarpus base. The complex trabecular structure is visible
to the microstructure geometry while the matrix material is isotropic. The consequence of the anisotropy exhibited by trabecular bone is that even the elastic behaviour description becomes a prohibitive task, involving up to 21 elastic constants. The experimental work of Brear et al. [5] showed that there is a linear relationship between the compressive strength of cancellous bone and its Young modulus, suggesting that inferences can be drawn between the failure properties and the elastic modulus measures of trabecular bone. Hodgskinson and Currey [6] could explain approximately 90% in the variance in Young’s modulus in cubes of human and non-human cancellous bone, using four variables: apparent density, mineral volume fraction of the bone tissue and two parameters describing trabecular architecture, fabric and connectivity, which were defined and computed on the bone cube sides. More recently, with the recent increasing performance of in vitro medical imaging techniques, the acquisition of bone volumes with accuracy that enables the definition of trabecular architecture has become possible [7]. It is therefore possible to employ micro-scale numerical models for the so-called computational material testing [8]. In literature, the 3D numerical analysis of bone structures has been based on a Finite Element Method (FEM) approach. Starting from the first models consisting of idealized rods arranged in an irregular lattice, more accurate representations have been developed using micro-FEM analysis
Meccanica (2007) 42:85–93
[9–14]. Although FEM represents a powerful and widely accepted method, there are some drawbacks involved in its application. In fact, the use of FEM implies that the size of the mesh must be smaller than any characteristic length involved in the problem [15]. This leads to a large number of elements when modelling a complex structure. Even though appropriate techniques have been developed [16], memory allocation and computing time requirements still represent a heavy burden in FEM modelling of trabecular bone. This paper presents the results of numerical simulations performed on 3D structures obtained from micro-computed tomographies of trabecular bone. The same 3D microstructures used in the numerical simulations were reproduced by means of a rapid prototyping process by selective laser sintering of polyamide powder. The apparent Young modulus of the examined microstructure was computed simulating a compression test and assuming elastic, homogeneous and isotropic properties for the base material. The numerical method used in this work is the Cell Method (CM) [17]. A brief description of the method and details of the simulations are given in the following section. The replicas, in which the side length of the original structure was magnified several times, were used in mechanical tests as a validation tool for the numerical simulations.
2 Materials 2.1 Samples Three cylindrical samples of trabecular bone extracted from different sites on the same pig shoulder are the starting point of this investigation (Fig. 2). Wolff’s law states that bone is deposited and reabsorbed in accordance with the stresses placed upon it [18]. The differences in both intensity and direction of the load are reflected in different organizations of the bone microstructure of the samples. Samples (1) and (2) were extracted from in vivo weight bearing regions, while sample (3) came from an in vivo weakly loaded region (Fig. 2). The samples measure 8–10 mm in diameter and
Meccanica (2007) 42:85–93
87
Fig. 2 Bone samples: (a) and (b) from in vivo weight bearing regions, (c) from in vivo weakly loaded region
5–10 mm in height. The mineral matrix was extracted from the biological tissue removing all the lipid and protein part by means of a chemical treatment (formaldehyde and hypochlorite).
A subset of these volumes (2003 voxel) was processed and used as input for the numerical simulations with the CM and elaborated for the rapid prototyping of the replicas.
2.2 Data acquisition
2.3 SLS replicas
Accuracy of structural computation strongly depends on the spatial resolution of the acquisition system. The samples were subjected to µCT (Computed Tomography at sub-millimetric scale) with synchrotron radiation at SYRMEP (ELETTRA, Trieste, Italy). SYRMEP is a beamline devoted to biomedical applications and employs a laminar X-ray beam coming out of one of ELETTRA bending magnets. A CCD camera equipped with a 20 µm Gadolinium Oxy-sulfide screen was used as 2D detector, providing wide dynamic range, linearity over many orders of magnitude, low-noise levels and small pixel size. The experimental set-up is shown in Fig. 3. By rotating the object in front of the detector, a total number of 720 planar projections were collected as the sample was rotated from 0 to 180o . These projections were used to obtain the cross-sectional images (tomographies) of the specimens by using a filtered back-projection algorithm [19]. By stacking the 2D images, a 3D tomographic reconstruction with a virtual spatial resolution of approximately 14 µm was obtained (Fig. 4). The trabecular structure can therefore be visualized in 3D: as already mentioned, different architectures, consequent to different loading conditions, are clearly visible in the reconstructed images of the samples. In Fig. 5, the sample extracted from the in vivo weakly loaded region exhibits a structure which is thinner and more porous than the other samples.
The 3D numerical analysis of bone structures can lead to results that differ from experimental ones because the tests performed on bone samples are affected by intrinsic difficulties such as specimen geometry, friction at endplates, structural end phenomena, storage, continuum assumption, viscoelasticity and temperature effects [20]. All these effects are amplified by the uniqueness of the sample and the impossibility to repeat the test. The consequences are errors in stiffness and strength determination, which can both mask real physical changes in the microstructure and give mismatches with the numerical analyses results. With the aim of overcoming these difficulties, the 3D reconstructions of the trabecular samples
Fig. 3 The experimental set-up at SYRMEP beamline (ELETTRA, Trieste, Italy)
88
Meccanica (2007) 42:85–93
Fig. 4 Reconstruction of the trabecular bone (2.1 × 2.1 × 2.1 mm3 )
have been reproduced by means of a Selective Laser Sintering (SLS) rapid prototyping process and the experimental tests performed on the replicas were used as a validation tool for the numerical simulations. Rapid prototyping is known to be a technology that can be used to manufacture complex 3D structures in a relatively homogeneous and wellcharacterized material, by progressively building up wafer-thin horizontal cross-sections (slices), one on top of another. Out of the several processes and materials that can be used for rapid prototyping, SLS appeared to eligible for trabecular bone reconstruction. The SLS process employs a powder, which can be fused by the heat produced by a laser. Figure 6 shows the process. The model is built on a platform which can be lowered a small amount so that the powder can be spread in the recess. The entire fabrication chamber is sealed and maintained at a temperature just below the melting point of the plastic powder. A laser directed across by a scanning mirror elevates the temperature slightly in order to sinter the powder while drawing the shape of the
cross-section. The process is then repeated. Complex structures with internal open cavities, such as the trabecular architecture, can be built [21]. In this process, the conversion of the reconstructed 3D image into a CAD object represents a critical step, as discussed in [22] and a dedicated procedure [23] was developed by the authors for this purpose. Polyamide powder and a CO2 laser were employed. The laser spot dimension was 0.6 mm and the layer thickness 0.15 mm. However, when hitting the powder, the area influenced and sintered by the heat was 0.75 mm in diameter. Moreover, when during the creation process the powder is heated and then cools down to room temperature, the parts shrink but the magnitude of the contraction is not uniform in all directions. The problems due to the solid contraction effect are automatically solved by the software, which scales and adjusts the laser trajectory so that the resulting object, at room temperature, will have the desired dimensions. Creation rate can be estimated as 15 mm/h. The replicas magnified the original structures side from 2.8 to 50 mm. The magnification factor was chosen
Fig. 5 3D reconstructions of trabecular bone from: (a) and (b) load bearing sites, (c) weakly loaded site (2.8 × 2.8 × 2.8 mm3 )
Meccanica (2007) 42:85–93
89
Fig. 6 The SLS process: (a) the platform is lowered; (b) the powder is applied; (c) the powder is sintered by exposure to CO2 laser
as a compromise between precision and costs. One of the replicas obtained from the micro-CT reconstructions previously discussed is shown in Fig. 7. 3 The Cell Method The CM is a recently introduced numerical method, which can be regarded as an alternative to more widely used numerical methods. A brief description follows, while details on the formulation can be found in [24, 25]. The CM comes from two considerations: (1)
(2)
A classification of the variables involved in a field problem into configuration, source and energy variables. This classification can be applied to any physical problem, and a collection of diagrams can be found at [26]. For structural problems, configuration variables include displacements, velocities, strain tensor, while forces, momenta and stress tensor are source variables. The recognition that variables of different classes are associated to different geometrical elements of the mesh.
More precisely, the mesh obtained by connecting the nodes is called the primal complex of cells,
and is analogous to the FEM one. But besides this mesh, the CM states that each node is surrounded by a second cell, belonging to a geometrical complex in which each cell constitutes an influence region of the node (Fig. 8). Different choices are possible for the geometrical complexes of cells, the only condition being that one and only one node falls inside a dual cell. In this work tetrahedral primal cells have been adopted, while the dual complexes have been always obtained by connecting the barycentres of the primal cells and of their sides and edges, as shown in Fig. 8. Thus, configuration variables are linked to geometrical elements belonging to a primal complex of cells and source variables are linked to the geometrical elements of the dual complex of cells. Tetrahedral cells have been used as primal cells in this work, because they allow the best curve shapes modelling. In this paper, an affine (linear function) interpolation is assumed, so that: {ε}c = [B]c {u}c ,
(1)
where {u} collects the nodal displacements. It can be easily verified that with this choice of interpolation function and cells, the matrix [B]c is the same as in FEM.
Fig. 7 From the micro-CT reconstruction (left) to the sintered replica (right)
90
Meccanica (2007) 42:85–93 Two dimensions
Three dimensions
Primal cells
Dual cells
Fig. 8 Examples of possible primal and dual cells. Tetrahedral cells, like the 3D cells shown here, have been used in this work
The stress tensor components can then be expressed by introducing the constitutive equation. Assuming a linear-elastic isotropic behaviour of the material: {σ }c = [D]c {ε}c = [D]c [B]c {u}c ,
(2)
where {σ }c collects the stress components and [D]c represents the Hooke’s law for the primal cell. It is then possible to express the forces acting on each side of the dual cell and obtain: {T}c = −Vc [B]T c [D]c [B]c {u}c .
(3)
Equilibrium condition is then written by adding up all the contributions and the external forces acting on the dual cell, so that balance equations are written for each influence region without applying a differential formulation. Since each node is surrounded by a dual cell, the equilibrium condition for each dual cell induces the same number of equations as in FEM. These equations can be rearranged in a linear system in the usual form [K] {u} = {F},
(4)
that can be solved with the usual methods. In general, accuracy and convergence rate of the numerical results obtained with the CM in structural problems are comparable with, and in some cases better than, the ones achieved with Finite
Element [24, 25], but the approaches on which the two methods are developed are different. In fact, the FEM writes balance laws for a point, the node, introducing derivatives of the variables of the approximated field. Thus, a two-step process is required: the problem is first formulated as a differential problem, and discretization is applied as a second step, for the numerical solution. Since the differentiation operation is not possible in presence of discontinuities, special elements must be developed in FEM to overcome this problem: as a general statement, the use of FEM implies the size of the mesh to be smaller than any characteristic length involved in the problem [15]. On the contrary, the CM writes balance equations for a region, the dual cell, not for a point: a direct discrete formulation is used. Models where heterogeneities are the same size of the cell and the constitutive matrix varies freely from one cell to the neighbour can be solved. Cell Method models of industrial materials such as sintered alloys and composites have already produced results in good agreement with the experimental ones [27–29]. In this work, the volume occupied by the examined structure was meshed with 812,905 tetrahedral primal cells (141,982 nodes). The same mesh was used in all the simulations. The mechanical behaviour of each cell was assumed elastic, homogeneous
Meccanica (2007) 42:85–93
91
Fig. 9 From the micro-CT reconstruction (left) to the numerical model (right)
and isotropic. The elastic modulus of each cell was determined by a scaling procedure based on the grey level in the tetrahedra vertexes and barycentres, as shown in Fig. 9. The assumed elastic modulus of the base material was 700 MPa, with Poisson ratio ν = 0.3, in agreement with the mechanical properties of the base material used to manufacture the replicas. Simulated mechanical compression tests were performed along the axes of the model, and the apparent elastic modulus of the examined bone structure was computed along the three coordinate axes. A factor that must be pointed out is the influence of data processing. The raw data matrix must be processed in two different ways in order to obtain the tetrahedral cells of the numerical model and the 3D structure required for the sintering process. To ensure the better possible match between experimental and simulation results, the characteristic length of the discretization of the numerical model (around 56 µm) was chosen to be of the same order of magnitude of the equivalent resolution due to the manufacturing process (approximately 66 µm).
The numerical results, shown in Table 1, were obtained by simulating the compression tests. The computation time was 2.5 h for each axis on a ordinary PC. The experimental results were obtained as average of two compression tests performed at University of Trieste with a Shimadzu machine on the sintered replicas. As an example, the stress/strain curves obtained in the compression tests along z-axis for sample (a) are shown in Fig. 10. The experimental and numerical results are compared in Fig. 11. Except in a couple of cases, the discrepancy was less than 5% in average. This discrepancy can be explained by the intrinsic limitations in the SLS procedure used, particularly with regard to the variability of manufacturing parameters. In fact, the prototypes were manufactured simultaneously and were therefore not expected to differ in powder composition, but a posteriori it was recognized that differences among the replicas occurred during the sintering process, as the temperatures achieved in the central part of the chamber are different from the temperatures in regions close to the lateral
4 Results and discussion The same reconstructed trabecular structure volumes were used for simulations and experimental tests. The parameter of the comparison was the compressive elastic modulus of the trabecular structure, computed along each axis of the reference frame. The triad was parallel to the reconstructed volume sides, with z-axis in the core sampling direction.
Fig. 10 Sample (1): stress/strain curves obtained in the compression tests along z-axis
92
Meccanica (2007) 42:85–93
Table 1 Numerical and experimental results(MPa) Simulations
Sample(a) Sample(b) Sample(c)
Tests
Ex
Ey
Ez
Ex
Ey
Ez
166 144 81
146 166 97
182 180 64
173 120 79
134 167 99
171 192 80
Fig. 11 Trabecular structure elastic modulus. Solid experimental, striped simulation
walls of the chamber. Consequently, prototypes created in different positions within the chamber exhibit different sintering degrees, so that the percentage of fused (sintered) powder may be slightly different among samples. The mechanical properties of the sintered replicas were recognized to depend on the temperatures achieved during manufacturing and hence on the position within the chamber. This can be considered the principal factor explaining the discrepancies between numerical and experimental results.
5 Conclusions The simulations results were in good agreement with the experimental tests, thus validating the application for the analysis of trabecular microstructures of the numerical model developed with the CM. The robustness exhibited by the CM in consequence of its direct discrete approach, allowed the characteristic length of the model to be of the same order of magnitude of the discontinuities of the structure. On one side, this consented the use of a relatively limited number of nodes, on the other side it was possible to create a single mesh for the three samples, which filled the entire examined volume. The mechanical properties of each
cell were determined by a scaling procedure based on the grey level in the cell itself, thus speeding up the numerical model creation process. The analysis requires a reasonably short computation time and does not need to be performed on a dedicated workstation. An important consequence will be the possibility to investigate bigger portions of bone volumes at trabecular level. Acknowledgements The Area Science Park project SISTER, founded by FVG Region, financially supported the work presented in this paper. The authors are also thankful to prof. Franco Vittur who provided the samples examined and to SYRMEP research group for support during images acquisition at Elettra (Trieste). Patent: PCT WO03/082118 2002, “Method to identify the mechanical properties of a material” deposited by University of Trieste.
References 1. Consensus Development Conference (1993) Diagnosis, prophylaxis and treatment of osteoporosis. AM J Med 94:646–650 2. Kleerekoper M, Villaneva AR, Stanciu J, Sudhaker Rao D, Parfitt AM (1985) The role of three dimensional trabecular microstructure in the pathogenesis of vertebral compression fractures. Calcif Tissue Int 37(Suppl):S594–S597 3. Uchiyama T, Tanizawa T, Muramatsu H, Endo N, Takahashi HE, Hara T (1999) Three-dimensional microstructural analysis of human trabecular bone
Meccanica (2007) 42:85–93
4. 5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
in relation to its mechanical properties. Bone 25(4): 487–491 Cowin SC (1985) The relationship between the elasticity tensor and the fabric tensor. Mech Mater 4:137-147 Brear K, Currey JD, Raines S, Smith KJ (1988) Density and temperature effects on some mechanical properties of cancellous bone. Eng Med 17: 163–167 Hodgskinson R, Currey JD (1990) The effect of variation in structure on the Young’s modulus of cancellous bone: a comparison of human and non-human material. J Eng Med PIME part H 204:115–121 ¨ ¨ Muller R, Hildebrand T, Ruegsegger P (1994) Noninvasive bone biopsy: a new method to analyse and display the three-dimensional structure of trabecular bone. Phys Med Biol 39:145–164 Zohdi TD, Wriggers P (2001) Computational micromacro material testing. Arch Comput Methods Eng 8:131–228 Yeh OC, Keaveny TM (1999) Biomechanical effects of intraspecimen variations in trabecular architecture: a three-dimensional finite element study. Bone 25:223–228 Van Rietbergen B, Weinans H, Huiskes R, Odgaar A (1995) A new method to determine trabecular bone elastic properties and loading using micromechanical finite-element models. J Biomech 28:69–81 Van Rietbergen B, Odgaar A, Kabel J, Huiskes R (1996) Direct mechanics assessment of elastic symmetries and properties of trabecular bone architecture. J Biomech 29:1653–1657 ¨ Ulrich D, van Rietbergen B, Weinans H, Ruegsegger P (1998) Finite element analysis of trabecular bone structures: a comparison of image-based meshing techniques. J Biomech 31:1187–1192 Niebur GL, Feldstein MJ, Yuen JC, Chen TJ, Keaveny TM (2000) High resolution finite element models with tyssue strength asymmetry accurately predict failure of trabecular bone. J Biomech 33:1575–1583 ¨ Hommiga J, Huiskes R, Van Rietbergen B, Ruegsegger P, Weinans H (2001) Introduction and evaluation of a gray-value voxel conversion technique. J Biomech 34:513–517 Roux S (1990) Continuum and discrete description of elasticity. In:Herrmann HJ, Roux S (eds) Statistical models for the fracture of disordered media. Elsevier, North Holland, 109–113
93 16. Van Rietbergen B, Weinans H, Huiskes R (1996) Computational strategies for iterative solutions of large FEM applications employing voxel data. Int J Num Methods Eng 39:2743–2767 17. Tonti E (2001) A direct discrete formulation of field laws: the cell method. CMES, Comput Model Eng Sci 2:237–258 18. Prendergast PJ, Huiskes R (1995) The biomechanics of Wolff’s law: recent advances. J Med Sci 164(2):152–154 19. Kak AC, Slaney M (1988) Principles of computerized tomographic imaging. IEEE Press, New york 20. Odgaard A (1997) Three-dimensional methods for quantification of cancellous bone architecture. Bone 4:315–328 21. Smith T (1999) Layer manufacture. A new form of 3-D visualization for artists. http://www.agocg.ac.uk/train/hitch/hitch2.pdf 22. Borah B, Gross GJ, Dufresne TE, Smith TS, Cockman MD, Chmielewski PA, Lundy MW, Hartke JR, Sod EW (2001) Three-dimensional microimaging (MRµI and µCT), finite element modeling, and rapid prototyping provide unique insights into bone architecture in osteoporosis. Anat Rec (New Anat) 206: 101–110 23. Cosmi F, Dreossi D, Scalici M (2004) Bone sintering for strength assessment: first results in BOSSA project. Trans FAMENA 28:35–42 24. Cosmi F (2001) Numerical solution of Plane elasticity problems with the cell method. CMES, Comput Model Eng Sc 2:365–372 25. Cosmi F (2005) Elastodynamics wih the cell method. CMES, Comput Model Eng Sci 8:191–200 26. Algebraic Formulation of Physical Fields (2006) http://www.dic.units.it/perspage/discretephysics/ 27. Cosmi F, Di Marino F (2001) Modelling of the mechanical behaviour of porous materials: a new approach. Acta Bioeng Biomech 3:55–65 28. Cosmi F (2004) Two-dimension estimate of effective properties of solid with random voids. Theor Appl Fract Mech, Elsevier Sci 42:183–186 29. Cosmi F (2003) Numerical modeling of porous materials mechnical behaviour with the cell method. In: Proceedings of Second MIT conference on computational fluid Massachusetts Institute of Technology, and solid mechanics, 17–20, June 2003, Cambridge, MA, USA
