Journal of Biomechanics 33 (2000) 255}259
Technical note
Improving the local solution accuracy of large-scale digital image-based "nite element analyses Guillaume T. Charras, Robert E. Guldberg* School of Mechanical Engineering, Georgia Institute of Technology, 281 Ferst Drive NW, Atlanta, GA 30332-0405, USA Received 4 June 1998; accepted 28 June 1999
Abstract Digital image-based "nite element modeling (DIBFEM) has become a widely utilized approach for e$ciently meshing complex biological structures such as trabecular bone. While DIBFEM can provide accurate predictions of apparent mechanical properties, its application to simulate local phenomena such as tissue failure or adaptation has been limited by high local solution errors at digital model boundaries. Furthermore, re"nement of digital meshes does not necessarily reduce local maximum errors. The purpose of this study was to evaluate the potential to reduce local mean and maximum solution errors in digital meshes using a post-processing "ltration method. The e!ectiveness of a three-dimensional, boundary-speci"c "ltering algorithm was found to be mesh size dependent. Mean absolute and maximum errors were reduced for meshes with more than "ve elements through the diameter of a cantilever beam considered representative of a single trabecula. Furthermore, mesh re"nement consistently decreased errors for "ltered solutions but not necessarily for non-"ltered solutions. Models with more than "ve elements through the beam diameter yielded absolute mean errors of less than 15% for both Von Mises stress and maximum principal strain. When applied to a high-resolution model of trabecular bone microstructure, boundary "ltering produced a more continuous solution distribution and reduced the predicted maximum stress by 30%. Boundary-speci"c "ltering provides a simple means of improving local solution accuracy while retaining the model generation and numerical storage e$ciency of the DIBFEM technique. ( 2000 Elsevier Science Ltd. All rights reserved. Keywords: Finite elements; Imaging; Trabecular bone
1. Introduction The inordinate time and e!ort required to generate three-dimensional (3D) "nite element models (FE) of complex biological structures has motivated the development of e$cient mesh generation techniques based upon reconstructed image data. One image-based modeling method directly converts individual image voxels into 3D hexahedral "nite elements by assigning element connectivity and material properties (Fyhrie et al., 1992; Hollister and Kikuchi, 1992; Keyak et al., 1990; Van Rietbergen et al., 1995). Digital image-based "nite element modeling (DIBFEM) allows fully automated generation of wellconditioned FE meshes regardless of structural complexity. Furthermore, this method can be combined with e$cient numerical storage strategies and iterative equa-
* Corresponding author. Tel.: #1-404-894-6589; fax: #1-404-8942291. E-mail address:
[email protected] (R.E. Guldberg)
tion solvers to analyze stresses, strains, and tissue properties in large-scale models containing several million "nite elements. DIBFEM has been applied primarily to model various levels of bone structure based on reconstructed computed tomography data (Guldberg et al., 1997; Hollister et al., 1994; Huiskes and Hollister, 1993; Keyak et al., 1993; Mullender et al., 1998; Van Rietbergen et al., 1995). The cost of such e$cient and robust mesh generation of complex structures is a digitized approximation of curved boundaries. The inexact boundary representation characteristic of digital FE models produces higher local solution errors relative to traditional smooth boundary models (Hollister and Riemer, 1993; Jacobs et al., 1993). The e!ects of digital element size, position, and orientation on local solution accuracy were recently examined (Guldberg et al., 1998). Analysis of error distributions revealed that local stress predictions within individual boundary elements oscillated approximately about a known theoretical solution. These results explain the observation that digital models provide accurate
0021-9290/00/$ - see front matter ( 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 1 - 9 2 9 0 ( 9 9 ) 0 0 1 4 1 - 4
256
G.T. Charras, R.E. Guldberg / Journal of Biomechanics 33 (2000) 255}259
averaged solutions (e.g. apparent modulus) while simultaneously producing high local solution errors. Recently, an error reduction strategy was introduced that smoothes digital meshes as a pre-processing step (Camacho et al., 1997). Although this technique improves solution accuracy, the large-scale analyses typically performed with fully digital models may not be possible with the smoothed meshes due to the altered geometry of the surface elements. The application of DIBFEM to simulate local phenomena such as tissue failure or adaptation is currently limited by high local solution errors. The purpose of this study was to investigate a post-processing strategy to improve the local accuracy of digital models using a 3D boundary-speci"c "ltering method.
2. Methods The theoretical stress and strain solutions for a 3D circular cantilever beam with linear elastic, isotropic material properties were used to quantify the local accuracy of digital model solutions (Timoshenko and Goodier, 1951). Digital FE models were generated at several mesh resolutions and orientations as described previously (Guldberg et al., 1998). Stress and strain predictions were obtained using VOX3D (Voxel Computing, Inc., Ann Arbor, MI), an iterative element by element preconditioned conjugate gradient FE solver. Predicted and theoretical solutions were determined at the centroid of each element. Solution error was de"ned as the di!erence between the theoretical (p ) and preTH dicted (p ) solutions normalized by the mean theoretical FE solution (1). This normalized de"nition re#ects the relative importance of the local errors and avoids overemphasizing small solution di!erences in regions of low theoretical stress or strain. The mean absolute and maximum errors were determined within a consistent midlength section for each beam model, p !p FE 100 %Error" TH p TH
(1)
A 3D boundary-speci"c "ltering algorithm was developed that allowed several parameters to be varied. In a three-dimensional cubic construct of elements, each element node may be connected to a maximum of eight elements. Nodes connected to less than eight elements were considered to be boundary nodes. Boundary elements could be de"ned as containing at least one, two, or four boundary nodes. Neighboring elements could be de"ned as possessing one, two, or four nodes in common with the boundary element being processed. The "lter de"nition in terms of neighborhood size and weighting was also varied.
3. Results Although several boundary-speci"c "lter de"nitions investigated improved local solution accuracy, results are shown for a single "lter de"nition only. An arithmetic mean "lter that de"ned a boundary element as containing at least two boundary nodes and neighboring elements as having at least two common nodes with the boundary element was found to be e!ective. Filtered values for individual elements with n neighbors were therefore calculated according to Eq. (2) in which p and p represent the un"ltered stress predicU/&*-5%3%$ * tions for the element and its n neighbors, respectively, p #+n p i/1 i . p " U/&*-5%3%$ F*-5%3%$ n#1
(2)
Oscillatory deviations from the theoretical solutions were evident on the boundary of the digital FE models prior to "ltering. Qualitatively, boundary-speci"c "ltering markedly improved the predicted distribution of Von Mises stress and principal strain. Quantitatively, "ltering reduced both mean absolute and maximum errors, especially at higher mesh resolutions Fig. 1. The bene"cial e!ects of "ltering were most evident for mesh resolutions greater than "ve elements through the beam thickness. The boundary-speci"c "lter was also applied to a 3D model of trabecular bone microstructure generated from microcomputed tomography image data. Fig. 2 shows the predicted Von Mises stress distribution due to a uniaxial compressive load before and after "ltering. Filtering produced a more continuous predicted stress distribution on trabecular surfaces and reduced the predicted maximum stress by approximately 30%.
4. Discussion Digital image-based FE methods o!er the potential to e$ciently generate complex 3D meshes of biological structures. Predictions of trabecular bone apparent and average tissue level properties using DIBFEM agree well with experimental data (Hollister et al., 1994, Van Rietbergen et al., 1995). However, the jagged mesh boundaries characteristic of digital models limit the accuracy of local stress and strain predictions at tissue surfaces. Oscillatory boundary solution patterns have suggested that numerical smoothing techniques might reduce solution artifacts (Guldberg et al., 1998; Hollister and Riemer, 1993). The e!ectiveness of error reduction was evaluated based on calculations of maximum and mean absolute errors within a mid-length section of a cantilever beam. The mean absolute error allowed the e!ects of boundary-speci"c "lters to be isolated from random averaging of positive and negative errors. In practical use, investiga-
G.T. Charras, R.E. Guldberg / Journal of Biomechanics 33 (2000) 255}259
257
Fig. 1. Fitted functions showing mean absolute and maximum Von Mises Stress errors (1A and 1B) and principal strain errors (1C and 1D) as a function of digital mesh resolution. Boundary-speci"c "ltering decreased both mean absolute and maximum errors for beams with more than "ve elements across the beam diameter.
tors may be most interested in the solution accuracy within a localized region such as through the thickness of single trabecula. Fig. 3 therefore plots absolute mean stress and strain errors at several mesh resolutions and orientations. As expected, the absolute mean errors are about 10% lower than the mean absolute errors. A potential limitation of this study is that the reported error magnitudes were determined using a single circular cross-section cantilever beam structure. Unlike trabecular bone models, this problem has a known theoretical solution that provides a basis for quantifying local FE solution errors. The use of DIBFEM to model other structures under other loading conditions yields di!erent error magnitudes. Axial rather than #exural loading conditions, for example, produce less extreme boundary errors. However, the oscillatory pattern of boundary errors is consistently observed and therefore error reduction strategies shown to be e!ective on the cantilever beam should be generally applicable. The errors that a!ect digital image-based "nite element predictions stem from several di!erent origins. At higher resolutions, oscillatory boundary errors resulting from the digital approximation of curved surfaces are dominant. At lower resolutions, additional errors can be introduced due to volume errors or numerical sti!ening
errors for structures loaded in bending. Numerical sti!ening errors in#uence all types of insu$ciently discretized "nite element meshes and have motivated the development of reduced integration element formulations (Bathe, 1982). Although DIBFEM may be utilized to solve models containing several million elements, microstructural features such as single trabecula often remain insu$ciently discretized. Reduced integration may be considered as a potential means to compensate for numerical sti!ening artifacts in these models. However, the bene"cial e!ects of reduced integration element formulations are directionally dependent while, for many biomechanical problems, the orientation between the structure and digital "nite elements cannot be speci"ed. This study con"rmed (data not shown) that reduced integration is not an e!ective error reduction strategy for cases in which the digital "nite elements are rotated with respect to the modeled structure. For aligned structures, however, reduced integration decreased local errors by about 10%, supporting further investigation of non-directionally dependent reduced integration schemes. The present study demonstrates that boundary-speci"c "ltering can be utilized to markedly improve the distribution of digital model solutions and reduce maximum errors, thereby addressing a signi"cant shortcoming
258
G.T. Charras, R.E. Guldberg / Journal of Biomechanics 33 (2000) 255}259
Fig. 2. Application of boundary-speci"c "ltering to the Von Mises stress solution for a 3D digital model of trabecular bone microstructure loaded in compression. The contour range shown was selected to illustrate that "ltering produced a more continuous stress on trabecular surfaces. Filtering also reduced high-stress concentrations by approximately 30%.
of DIBFEM. Importantly, mesh re"nement consistently decreased errors for "ltered solutions but not necessarily for non-"ltered solutions. Hollister and Riemer (1993) reported Gaussian "lters were not e!ective at reducing DIBFEM solution errors, however, these "lters were applied to the entire digital model. The speci"city of the "ltering algorithm to the digital model boundary may therefore be critical. Using this approach, relatively simple post-processing "lters may be applied to DIBFEM solutions for applications that require more accurate estimates of local stress or strain distributions. Fig. 3. Absolute mean Von Mises stress and principal strain errors as a function of digital mesh resolution. The "tted functions suggest that models with more than "ve elements through the beam diameter yielded absolute mean errors of less than 15% for both Von Mises stress and maximum principal strain.
Acknowledgements This work was supported by The Whitaker Foundation. Finite element modeling and analysis software
G.T. Charras, R.E. Guldberg / Journal of Biomechanics 33 (2000) 255}259
provided by Voxel Computing, Inc., Ann Arbor, Michigan. The authors would like to thank Scott Hollister for his advise and review of the manuscript. References Bathe, K.J., 1982. In Finite Element Procedures in Engineering Analysis, 1st Edition. Prentice Hall, Eaglewood Cli!s, NJ. Camacho, D.L.A., Hopper, R.H., Lin, G.M., Myers, B.S., 1997. An improved method for "nite element mesh generation of geometrically complex structures with application to the skullbase. Journal of Biomechanics 30 (10), 1067}1070. Fyhrie, D.P., Hamid, M.S., Kuo, R.F., Lang, S.M., 1992. Direct three-dimensional "nite element analysis of human vertebral cancellous bone. Orthopaedic Research Society Transactions 12, 551. Guldberg, R.E., Caldwell, N.J., Guo, X.E., Goulet, R.W., Hollister, S.J., Goldstein, S.A., 1997. Mechanical stimulation of tissue repair in the hydraulic bone chamber. Journal of Bone and Mineral Research 12 (8), 1295}1302. Guldberg, R.E., Hollister, S.J., Charras, G.T., 1998. The accuracy of digital image-based "nite element models. Journal of Biomechanical Engineering 120, 289}295. Hollister, S.J., Kikuchi, N., 1992. Direct analysis of trabecular bone sti!ness and tissue level mechanics using an element-by-element homogenization method. Orthopaedic Research Society Transactions 12, 559.
259
Hollister, S.J., Riemer, B.A., 1993. Digital image based "nite element analysis for bone microstructure using conjugate gradient and Gaussian "lter techniques. Mathematical Methods in Medical Imaging II}SPIE Proceedings 2035, 95}106. Hollister, S.J., Brennan, J.M., Kikuchi, N., 1994. A homogenization procedure for calculating trabecular bone sti!ness and tissue level stress. Journal of Biomechanics 27, 433}444. Huiskes, R., Hollister, S.J., 1993. From structure to process, from organ to cell: recent developments of FE-analysis in orthopaedic biomechanics. Journal of Biomedical Engineering 115, 520}527. Jacobs, C.R., Mandell, J.A., Beaupre, G.S., 1993. A comparative study of automatic "nite element mesh generation techniques in orthopaedic biomechanics. Bioengineering Conference, ASME BED 24, 512}514. Keyak, J.H., Meagher, J.M., Skinner, H.B., Mote, C.D., 1990. Automated three-dimensional "nite element modelling of bone: a new method. Journal of Biomedical Engineering 12 (5), 389}397. Keyak, J.H., Fourkas, M.J., Meagher, J.M., Skinner, H.B., 1993. Validation of an automated method of three-dimensional "nite element modelling of bone. Journal of Biomedical Engineering 15 (6), 505}509. Mullender, M., van Rietbergen, B., RuK egsegger, P., Huiskes, R., 1998. E!ect of mechanical set point of bone cells on mechanical control of trabecular bone architecture. Bone 22 (2), 125}131. Timoshenko, S., Goodier, J.N., 1951. Theory of Elasticity, 2nd Edition. McGraw-Hill, New York. Van Rietbergen, B., Weinans, H., Huiskes, R., Odgaard, A., 1995. A new method to determine trabecular bone elastic properties and loading using micro-mechanical "nite-element models. Journal of Biomechanics 28 (1), 69}81.