Medical Engineering and Physics 37 (2015) 808–812

Contents lists available at ScienceDirect

Medical Engineering and Physics journal homepage: www.elsevier.com/locate/medengphy

Technical note

Comparisons of node-based and element-based approaches of assigning bone material properties onto subject-specific finite element models G. Chen a,∗, F.Y. Wu a, Z.C. Liu a, K. Yang a, F. Cui b a b

School of Civil and Transportation Engineering, Guangdong University of Technology, Guangzhou 510006, China Institute of High Performance Computing, A∗ STAR, Singapore

a r t i c l e

i n f o

Article history: Received 21 January 2015 Revised 8 April 2015 Accepted 2 May 2015

Keywords: Bone Computed tomography Finite element models Material assignment Node-based approaches Element-based approaches

a b s t r a c t Subject-specific finite element (FE) models can be generated from computed tomography (CT) datasets of a bone. A key step is assigning material properties automatically onto finite element models, which remains a great challenge. This paper proposes a node-based assignment approach and also compares it with the element-based approach in the literature. Both approaches were implemented using ABAQUS. The assignment procedure is divided into two steps: generating the data file of the image intensity of a bone in a MATLAB program and reading the data file into ABAQUS via user subroutines. The node-based approach assigns the material properties to each node of the finite element mesh, while the element-based approach assigns the material properties directly to each integration point of an element. Both approaches are independent from the type of elements. A number of FE meshes are tested and both give accurate solutions; comparatively the node-based approach involves less programming effort. The node-based approach is also independent from the type of analyses; it has been tested on the nonlinear analysis of a Sawbone femur. The nodebased approach substantially improves the level of automation of the assignment procedure of bone material properties. It is the simplest and most powerful approach that is applicable to many types of analyses and elements. © 2015 IPEM. Published by Elsevier Ltd. All rights reserved.

1. Introduction Subject-specific finite element (FE) models have been commonly employed to study the mechanical behaviours of bone structures. Both model geometries and material properties of subject-specific FE models can be automatically generated from computed tomography (CT) datasets of a bone. Subject-specific FE models can be created by two strategies: (1) voxel-based and (2) geometry-based [1,2]. For voxel-based models, FE datasets, including bone material properties, are directly generated from CT images using a commercial code, e.g., Simpleware (Simpleware Ltd., Exeter, UK). The apparent disadvantages of this strategy include low accuracies in FE model geometries and bone material properties [2]. The geometry-based strategy is more commonly used in bone simulations [1], and the standard procedures for this strategy were reported, for example, in Taddei et al. [3] and Kluess et al. [4]. Generation of FE model geometries can be realised using a couple of commercial software packages unitedly, as reported in Messmer et al. [5]; however, assigning material properties onto FE models usually needs an in-house program code (e.g., BONEMAT) to connect the CT dataset and FE datasets [6,7]. The automation



Corresponding author. Tel.: +86 136 6248 3527 E-mail address: [email protected] (G. Chen).

http://dx.doi.org/10.1016/j.medengphy.2015.05.006 1350-4533/© 2015 IPEM. Published by Elsevier Ltd. All rights reserved.

of assigning material properties is influenced by the implementation of assignment approaches [8] and great efforts have been made to improve the level of assignment automation and accuracy of material properties [3,4,6,7,9–12]. Several assignment approaches have been reported for mapping bone material properties from CT images onto finite element models for stress analyses [4,6–9,13,14]. In the early stage, the researchers developed in-house assignment programs [7,14], which usually connect the CT dataset and FE dataset to assign the Young’s modulus to each element (i.e., element-based). As each element type has a different data structure in the FE models, hence the element-based assignment approach is dependent on the type of elements used; consequently, the implementation of this approach needs complex programming. Obviously the implementation of the element-based approach is also dependent on the type of analyses, because a different analysis owns a different data structure in the FE models. To make the assignment process independent from the type of elements, Helgason et al. [8] and Kluess et al. [4] assigned the intensity value of CT images as the “temperature” at the FE nodes (i.e., node-based) and made the bone material properties (i.e., Young’s modulus) dependent on the “temperature”. This approach is independent from the type of elements but dependent on the type of analyses. The obvious disadvantage of all the above-mentioned approaches is that, assignment of material properties from the CT dataset to a FE model of a bone involves

G. Chen et al. / Medical Engineering and Physics 37 (2015) 808–812

complex programming, as these methods directly modify the FE dataset. Recently, the assignment procedure has been substantially simplified by Chen et al. [9]. Instead of modifying the dataset of FE models, their assignment approach resorts the special functionality of commercial codes and assigns the inhomogeneous material properties via a user-defined material. Their assignment procedure is divided into two steps: generating the data file of the image intensity of a bone in a MATLAB (MathWorks, Inc., Natick, MA, USA) program and reading the data file into ABAQUS (SIMULIA, Inc., Providence, USA) via user subroutines. The material properties are directly assigned to each integration point of an element (element-based). This assignment approach is independent from the type of elements [9]; however, it is dependent on the type of analyses. For a normal linear elastic analysis, the approach involves very little programming; nevertheless, for complex analyses, e.g., plastic analyses, the users may have to develop complicated subroutines to realise their analyses. To avoid complicated programming, an assignment approach, which is independent from the type of elements and the type of analyses, is desired. This paper proposes a node-based approach to assign bone material properties, in which material properties are defined on the FE nodes by using a field variable and the field variable is calculated from the Hounsfield units (HUs) of bone CT images. It also incorporates the idea of splitting the assignment procedure into two steps [9], i.e., the image intensity of a bone is generated in a MATLAB program and saved in a data file which is read into ABAQUS via user subroutines. The node-based approach will be tested on different mesh types and different analyses (linear and nonlinear). The assignment accuracy of this approach will be verified by comparisons with the elementbased approach proposed by Chen et al. [9]. 2. Methods A Sawbone femur model (Pacific Research Laboratories, Inc., Washington, United Sates) was scanned using a Brilliance 64 CT scanner (Philips, Best, the Netherlands), operating at 120 KVP with a 0.23 × 0.23 mm in-plane pixel size and with a slice thickness of 0.67 mm. The images were saved in 526 DICOM files. Both the geometry and material properties of a FE model were derived from the CT dataset of the Sawbone femur. 2.1. Acquisition of bone density from CT images All CT images of the Sawbone femur were first processed by an in-house MATLAB program [9] to obtain the intensity of the CT images and the values were saved into a data file. In the MATLAB program, three tasks were performed: (1) 2D segmentation of bone for each slice; (2) correction of partial volume effect (PVE) on the boundary pixels; (3) output of the pre-processed images and their HUs. To extract an accurate FE model geometry from CT images, a multiplethreshold segmentation method was used. The FE model geometry generated in this segmentation method is more accurate than that created by single-threshold segmentation. The multiple-threshold segmentation techniques and correction of PVE used in this paper had been reported in Chen et al. [9]. The pre-processed CT images were used for creation of Sawbone geometry in Amira (Mercury Computer Systems, Inc., Chelmsford, USA) and the HUs of CT images were saved into a TEXT file, which will be read by ABAQUS. 2.2. Creation of FE model geometry from CT images A standard procedure to create virtual 3D models of the bones, using commercial software packages, has been reported in Messmer et al. [5]. First, the visualisation software Amira was used to extract the external contour of the Sawbone femur from the pre-processed

(a)

809

(b) A pixel in a CT slice

A node in an element

A FE element

A Gauss integration point

Fig. 1. Assigning material properties from CT dataset to a FE model. (a) The FE geometry imported into ABAQUS. The arrow indicates the loading direction. (b) Assigning material properties from CT dataset to a Gauss integration point or a node.

CT images; the results were saved in a STL file. This surface model was then converted into a NURBS model using the software Rapidform (INUS Technology, Seoul, Korea). A solid was then created from the NURBS model in Rapidform and saved as an IGES file. This IGES file was imported by ABAQUS for creation of the geometry of the FE model (Fig. 1a). 2.3. Finite element model The femur model was fixed at the distal end and subjected to lateral forces of 400 N at the proximal end along the X direction (Fig. 1a). To avoid apparent stress concentration, the loads were applied at the four nodes which were located at (−105.5, 158.4, 508.6), (−96.7, 168.4, 536.1), (−93.2, 140.5, 546.8), and (−102.1, 130.1, 519.2). The four-node tetrahedral elements were used for the major part of the bone and the average aspect ratio was 1.6. Both the fournode tetrahedral elements and eight-node hexahedral elements were tested in the middle part of the bone where the comparisons were made. Two different element sizes, rough (2.0 mm) and fine (1.0 mm), were used for stress analysis. For the nonlinear stress analyses, the yield stress of the cortical part was σy = 198.0 MPa, the tangential modulus was E p = 0.14E. The von Mises yield criterion and isotropic hardening law were used [15]. The loading position and direction were kept as above and the value of the load was doubled. The material properties were assigned by either the elementbased or node-based assignment approaches onto the FE models during analyses, which were explained in the following. 2.4. Assignment of material properties to finite element model ABAQUS read in the text file which saves the HUs via a subroutine “UEXTERNALDB” at the beginning of the analysis to calculate the material properties of the Sawbone femur. •

Element-based assignment approach [9]. A user-defined material was assigned for the bone. For each Gauss integration point, the pixels enclosing this point were searched (Fig. 1b) during the analysis and the average HU of the surrounding (3 × 3 × 3) pixels was calculated. This average HU was assigned to a solution-dependent state variable as the initial value at this integration point through a user-subroutine “SDVINI” in ABAQUS. The values of the solutiondependent state variable at integration points were used to evaluate the Young’s modulus in the user-subroutine “UMAT”, which calculates the density, stiffness, strain, and stress at each integration point [16].

810 •

G. Chen et al. / Medical Engineering and Physics 37 (2015) 808–812

Node-based assignment approach. The material properties of the bone were defined to be dependent on a pre-defined field variable. For each node, the pixels enclosing this node were searched (Fig. 1b) and the average HU of the surrounding (3 × 3 × 3) pixels was calculated. This average HU was assigned to the pre-defined field variable through a user-subroutine “UFIELD” in ABAQUS. ABAQUS will automatically interpolate its values at Gauss integration points from the nodal values. Users can define the relations between the Young’s modulus and the pre-defined field variable. The Young’s modulus of the bone was calculated from this predefined field variable.

2.5. Acquisition of bone density and Young’s modulus from HUs of CT images The bone density (ρ ) was calculated from the Hounsfield unit values (H) as ρ = a · H + b, where a and b were the calibration coefficients and ρ was in g/cm3 . For the Sawbone femur, the bone density at the cortical part was 1.7 g/cm3 [17]. Ten pixels at the cortical part were selected and the average H was calculated as 1110. Hence, a was set to 1.532 × 10−3 and b zero. The Young’s modulus (E) was related to the bone density (ρ ) by E = Emax · (ρ /ρmax )2 , where ρmax = 1.7 g/cm3 and Emax = 16.9 GPa which represented the maximum values of the apparent density and Young’s modulus of the cortical part [17]. 2.6. Tests of node-based and element-based assignment approaches

Fig. 2. Mapping results of bone density obtained by node-based and element-based assignment approaches. A contour plot in a XY plane of diaphyseal for rough and fine meshes. (A) Element-based approach. (B) Node-based approach. (For interpretation of the references to colour in the text, the reader is referred to the web version of this article.)

First, both the element-based and node-based assignment approaches were tested on linear stress analyses. Two different elements, tetrahedral and hexahedral, were used respectively, and the results were obtained and compared for the assignment approaches. Second, the node-based assignment approach was examined on a nonlinear stress analysis to demonstrate its applicability to various types of analyses. 3. Results First, both the node-based and element-based assignment approaches were applied to a linear analysis. The bone density and von Mises stress were obtained and compared for the two approaches. Furthermore, the node-based assignment approach was applied to a nonlinear analysis, where the element-based approach was difficult to be implemented. 3.1. Bone density obtained by node-based and element-based assignment approaches The mapping results of bone density for the two meshes were shown in Fig. 2, which was a cut-view along the Z axis in the diaphyseal part. The node-based and element-based assignment approaches give very similar results for both the rough and fine meshes. The rough mesh gives a slightly less smooth distribution as the elements were large. The average bone density for the diaphyseal part was calculated from 10 typical elements, which were randomly selected on the boundary with red colour. The results of the average bone density (in g/cm3 ) for the four cases were 1.61 (element-based approach, rough), 1.60 (element-based approach, fine), 1.62 (node-based approach, rough), 1.61 (node-based approach, fine) respectively. The errors in the solutions of all meshes were in the range of 1–2%. The average error between the element-based approach and nodebased approach was less than 1%. For the comparisons between different meshes, the locations of the 10 selected elements were approximately matching to each other.

Fig. 3. von Mises stress obtained with two assignment approaches using tetrahedral elements. A contour plot in a XY plane of diaphyseal for both rough and fine meshes. (A) Element-based approach. (B) Node-based approach. (For interpretation of the references to colour in the text, the reader is referred to the web version of this article.)

3.2. Stress obtained by node-based and element-based assignment approaches The contour plots of von Mises stress for the four meshes were shown in Fig. 3. The distributions of von Mises stress obtained with the node-based and element-based assignment approaches were

G. Chen et al. / Medical Engineering and Physics 37 (2015) 808–812

811

4. Discussion

Fig. 4. von Mises stress obtained by two assignment approaches using hexahedral elements. A contour plot in a XY plane of diaphyseal. (A) Element-based approach. (B) Node-based approach.

Fig. 5. Node-based assignment approach applied to non-linear stress analyses. The plastic strain in a cross section in diaphyseal.

very similar. The rough meshes give a slightly less smooth distribution than the fine meshes. The errors in the solutions of all meshes were in the range of 1–2% (for 10 typical elements compared). The average error between the element-based approach and node-based approach was approximately 1%. The relative errors were calculated for 10 elements randomly selected on the boundaries in the red regions. 3.3. Application to different types of elements The contour plots of von Mises stress, shown in Fig. 4, were obtained with hexahedral elements. The distributions of von Mises stress obtained by the two assignment approaches were very similar. They were also very similar to those obtained with tetrahedral elements (Fig. 3). The average error for 10 typical elements between the elementbased approach and node-based approach was approximately 1%. The errors in the solutions of two meshes (hexahedral and tetrahedral) were in the range of 1–2%. The relative errors were calculated for 10 elements randomly selected on the boundaries in the red regions. 3.4. Application of node-based approach to nonlinear analysis Yielding occurred in the cortical part at approximately 1/4 length of the femur from the distal ends. Fig. 5 showed the contour plot of the plastic strain in a cross section, on which the plastic strain had the maximum value.

Two strategies had been developed to introduce the spatial variations of bone material properties in FE models in the simulations of bone remodelling. One was the element-based, i.e., the material properties were defined on each Gauss integration point [18–21]; the other was the node-based, i.e., the material properties were defined on each node [20–23]. The element-based strategy has been applied by Chen et al. [9] to assigning material properties onto a FE model of a bone using its CT data, which substantially simplified the assignment procedure of material properties. As this assignment approach is dependent on the type of analyses, for nonlinear analyses, complicated programming will be involved. Hence, the question is whether the node-based strategy is applicable to assigning material properties onto a FE model from the CT dataset of a bone and whether it would behave better for complex analyses. This paper has tested the node-based strategy to assign the material properties of a bone to a FE model in ABAQUS. The proposed node-based approach successfully integrates the advantages of the assignment approaches previously proposed in literature. First, following Chen et al. [9], it divides the assignment procedure into two steps and utilises the special functionality of ABAQUS to read in the bone material properties, which avoids advanced programming for the modifications of the FE dataset. Second it inherits the idea of Helgason et al. [8] to make bone material properties dependant on a pre-defined field variable, which completely avoids the programming effort in calculations of stiffness matrix as in Chen et al. [9]. Successful combination of these two ideas makes the proposed approach independant from both the type of the elements and the type of the analyses. To authors’ knowledge, it is the simplest and most powerful approach that is applicable to various types of analyses and elements. It improves the level of automation of the assignment procedure significantly. The user subroutines, provided with this paper, are suitable to various types of elements and analyses; the users can apply them to their applications with minor changes in redefining the relations presented in Section 2.5. The proposed strategy is illustrative. The distributions of bone density and Young’s modulus can be visualised in ABAQUS CAE [9]. Comparisons have been made for the node-based and elementbased approaches. Both work for the linear elastic analyses flawlessly and give very close numerical results. The accuracy of the element-based approach proposed by Chen et al. [9] has been verified with a bone phantom. Hence, the comparisons of of the numerical results of the node-based and element-based approaches demonstrate that both approaches can accurately assign the inhomogeous material properties of bones onto FE models. Both approaches are independent from the type elements used in the FE models. The element-based approach assigns the HUs directly onto Gauss integration points and the node-based approach onto nodes. The elementbased approach relies on user-defined materials which involves the calculations of the stiffness matrix at Gauss integration points and requests advanced programming for complex analyses. While the nodebased approach assigns HUs onto each node directly as a field variable and makes the bone material properties dependant on this field variable, hence, it is independent from both the type of elements and the type of analyses. Creation of subject-specific models from image data is still expensive in terms of time, and usually requires some degree of user intervention [24]. It is significant to simplify and improve the level of automation of the assignment procedure. The examples presented in this study demonstrated that the proposed node-based approach can automatically and accurately assign bone material properties into FE models of various analyses without manual manipulations. This confirms the observations of Helgason et al. [8] that, the node-based method performs better than the element-based method in extraction of bone material properties from the CT data. The programming

812

G. Chen et al. / Medical Engineering and Physics 37 (2015) 808–812

References

Fig. 6. Mapping results of bone density obtained by node-based assignment approach. (A) Distal end; (B) proximal end.

effort is trivial and the codes are applicable to various types of element and types of analyses. The proposed node-based approach can be easily modified to fit users’ situations. It is based on ABAQUS. Obviously the ideas are applicable to other commercial FE codes. For other FE codes, new author-defined subroutines or programs need to be developed. However, once the programs have been developed, the users do not need to make major changes. FE models of long bones are generally limited by the ability to truly capture the cortical thickness in the metaphyseal and epiphyseal regions [25,26] and different formulas should be resorted to describe the relation of the Young’s modulus and bone density in cortical and trabecular regions. Since the aim of this paper is rather to propose a node-based assignment approach than to provide a novel CT-based FE mothod, it did not separate the cortical and trabecular parts in image segmentation. However, the mapping results of bone density obtained by the node-based assignment approach illustrate that the cortical and trabecular parts can be noticeably separated (Fig. 6). The bone density for the cortical part is in the range 1.2–2.0 g/cm3 and for the trabecular part below 0.4 g/cm3 . In summary, the proposed node-based approach can automatically and accurately assign the material properties onto FE models of a bone without manual manipulations. It significantly improves the level of automation and substantially simplifies the assignment procedure. It can be easily implemented and is suitable to various types of analyses and types of elements. Conflict of interest None declared. Ethical approval Not required. Acknowledgement This research was supported by the project (no. 31470908) of the National Natural Science Foundation of China. Supplementary Materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.medengphy.2015.05.006.

[1] Lengsfeld M, Schmitt J, Alter P, Kaminsky J, Leppek R. Comparison of geometrybased and CT voxel-based finite element modelling and experimental validation. Med Eng Phys 1998;20:515–22. [2] Kourtis LC, Carter DR, Kesari H, Beaupre GS. A new software tool (VA-BATTS) to calculate bending, axial, torsional and transverse shear stresses within bone cross sections having inhomogeneous material properties. Comput Methods Biomech Biomed Eng 2008;11:463–76. [3] Taddei F, Cristofolini L, Martelli S, Gill HS, Viceconti M. Subject-specific finite element models of long bones: an in vitro evaluation of the overall accuracy. J Biomech 2006;39:2457–67. [4] Kluess D, Souffrant R, Mittelmeier W, Wree A, Schmitz K-P, Bader R. A convenient approach for finite-element-analyses of orthopaedic implants in bone contact: modeling and experimental validation. Comput Methods Programs Biomed 2009;95:23–30. [5] Messmer P, Matthews F, Jacob AL, Kikinis R, Regazzoni P, Noser H. A CT database for research, development and education: concept and potential. J Digit Imaging 2007;20:17–22. [6] Zannoni C, Mantovani R, Viceconti M. Material properties assignment to finite element models of bone structures: a new method. Med Eng Phys 1998;20:735– 40. [7] Taddei F, Pancanti A, Viceconti M. An improved method for the automatic mapping of computed tomography numbers onto finite element models. Med Eng Phys 2004;26:61–9. [8] Helgason B, Taddei F, Pálsson H, Schileo E, Cristofolini L, Viceconti M, et al. A modified method for assigning material properties to FE models of bones. Med Eng Phys 2008;30:444–53. [9] Chen G, Schmutz B, Epari D, Rathnayaka K, Ibrahim S, Schuetz MA, et al. A new approach for assigning bone material properties from CT images onto finite element models. J Biomech 2010;43:1011–15. [10] Unnikrishnan GU, Morgan EF. A new material mapping procedure for quantitative computed tomography-based, continuum finite element analyses of the vertebra. J Biomech Eng 2011;133 071001. [11] Shim VB, Battley M, Anderson IA, Munro JT. Validation of an efficient method of assigning material properties in finite element analysis of pelvic bone. Comput Methods Biomech Biomed Eng 2014;18:1495–9. [12] Kim YH, Kim J-E, Eberhardt AW. A new cortical thickness mapping method with application to an in vivo finite element model. Comput Methods Biomech Biomed Eng 2014;17(9):997–1001. [13] Cattaneo P, Dalstra M, Frich L. A three-dimensional finite element model from computed tomography data: a semi-automated method. Proc Inst Mech Eng Part H: J Eng Med 2001;215:203–12. [14] Wijayathunga VN, Jones AC, Oakland RJ, Furtado NR, Hall RM, Wilcox RK. Development of specimen-specific finite element models of human vertebrae for the analysis of vertebroplasty. Proc Inst Mech Eng Part H: J Eng Med 2008;222:221–8. [15] Gong H, Zhang M, Fan Y, Kwok W, Leung P. Relationships between femoral strength evaluated by nonlinear finite element analysis and BMD, material distribution and geometric morphology. Ann Biomed Eng 2012;40:1575–85. [16] Zienkiewicz OC, Taylor RL. The finite element method. 5th ed. ButterworthHeinemann; 2000. [17] Cristofolini L, Viceconti M, Cappello A, Toni A. Mechanical validation of whole bone composite femur models. J Biomech 1996;29:525–35. [18] Doblare M, Garcia JM. Application of an anisotropic bone-remodelling model based on a damage–repair theory to the analysis of the proximal femur before and after total hip replacement. J Biomech 2001;34:1157–70. [19] Beaupre GS, Orr TE, Carter DR. Approach for time-dependent bone modeling and remodeling—application. A preliminary remodeling simulation. J Orthop Res 1990;8:662–70. [20] Chen G, Pettet G, Pearcy M, McElwain DLS. Comparison of two numerical approaches for bone remodelling. Med Eng Phys 2007;29:134–9. [21] Sharma GB, Debski RE, McMahon PJ, Robertson DD. Adaptive glenoid bone remodeling simulation. J Biomech 2009;42:1460–8. [22] Jacobs CR, Levenston ME, Beaupre GS, Simo JC, Carter DR. Numerical instabilities in bone remodeling simulations: the advantages of a node-based finite element approach. J Biomech 1995;28:449–59. [23] Levenston ME. Temporal stability of node-based internal bone adaptation simulations. J Biomech 1997;30:403–7. [24] Wilcox RK. The influence of material property and morphological parameters on specimen-specific finite element models of porcine vertebral bodies. J Biomech 2007;40:669–73. [25] Yosibash Z, Padan R, Joskowicz L, Milgrom C. A CT-based high-order finite element analysis of the human proximal femur compared to in-vitro experiments. J Biomech Eng 2007;129:297–309. [26] Trabelsi N, Yosibash Z, Milgrom C. Validation of subject-specific automated p-FE analysis of the proximal femur. J Biomech 2009;42:234–41.

Comparisons of node-based and element-based approaches of assigning bone material properties onto subject-specific finite element models.

Subject-specific finite element (FE) models can be generated from computed tomography (CT) datasets of a bone. A key step is assigning material proper...
2MB Sizes 1 Downloads 9 Views