Physical modeling with orthotropic material based on harmonic fields

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

journal homepage: www.intl.elsevierhealth.com/journals/cmpb

Physical modeling with orthotropic material based on harmonic fields Sheng-Hui Liao a,∗ , Bei-Ji Zou a , Jian-Ping Geng b , Jin-Xiao Wang c , Xi Ding c a b c

School of Information Science and Engineering, Central South University, Changsha 410083, China Institute for Biomedical Engineering, Nanjing University of Technology, Nanjing, China First Affiliated Hospital of Wenzhou Medical College, Wenzhou, China

a r t i c l e

i n f o

a b s t r a c t

Article history:

Although it is well known that human bone tissues have obvious orthotropic material

Received 21 September 2010

properties, most works in the physical modeling field adopted oversimplified isotropic or

Received in revised form

approximated transversely isotropic elasticity due to the simplicity. This paper presents a

7 April 2011

convenient methodology based on harmonic fields, to construct volumetric finite element

Accepted 24 April 2011

mesh integrated with complete orthotropic material. The basic idea is taking advantage

Keywords:

most bone tissues is compatible with the trajectory of the maximum material stiffness.

Finite element modeling

First, surface harmonic fields of the longitudinal axis direction for individual bone mod-

of the fact that the longitudinal axis direction indicated by the shape configuration of

Biomaterials

els were generated, whose scalar distribution pattern tends to conform very well to the

Anisotropy

object shape. The scalar iso-contours were extracted and sampled adaptively to construct

Orthotropic material

volumetric meshes of high quality. Following, the surface harmonic fields were expanded

Harmonic field

over the whole volumetric domain to create longitudinal and radial volumetric harmonic fields, from which the gradient vector fields were calculated and employed as the orthotropic principal axes vector fields. Contrastive finite element analyses demonstrated that elastic orthotropy has significant effect on simulating stresses and strains, including the value as well as distribution pattern, which underlines the relevance of our orthotropic modeling scheme. © 2011 Elsevier Ireland Ltd. All rights reserved.

1.

Introduction

The validity of the simulation result of physical finite element (FE) model is mainly dependent upon the geometrical similarity, as well as the material similarity of the model to the real structure of the object to be analyzed. Recent advances in the field of biomechanical modeling have made some interesting tools available [1]. One fully automatic approach allows direct generation of voxel meshes from CT data [2], but has shortcomings of inaccuracy in representing the geometry, and the large number of elements required to build the

∗

model [3]. Due to its rigorous adaptability to structures with geometric complexities, tetrahedral meshes are always quite popular and generated from medical images by various methods [4]. Besides, because a good hexahedral mesh can vastly reduce the number of elements and often perform better than tetrahedral mesh in the analysis stage, several authors tried to develop unstructured hexahedral models from medical images [5]. But almost all of those modeling techniques did not pay enough attention to the assignment of proper material properties to the FE model, which is also a fundamental step to ensure predictive accuracy. Accurate modeling of biological

Corresponding author. Tel.: +86 0731 85714036; fax: +86 0731 85714036. E-mail addresses: [email protected], [email protected] (S.-H. Liao). 0169-2607/$ – see front matter © 2011 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2011.04.005

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

tissues, such as bone related organs, is a very difficult task because of their inherent inhomogeneous and anisotropic characteristics [6]. The inhomogeneous property is in some sense “directly” accessible via the CT data based on the relationship between CT numbers, Hounsfield Unit (HU), and bone material properties [7,8]. Unfortunately, this is not the case for the trajectories of anisotropy.

2.

Background

Anisotropic material properties of bone tissues are important for related mechanobiology studies. For example, local anisotropy and regional variations in bone elastic properties can have pronounced effects on the relationship between stress and strain patterns, thus is crucial for both accurate mechanical models and to understand potential structural adaptation in the bone tissue [9]. Many researchers investigated the anisotropic material property of bone tissues by various experimental methods [9–14]. It is commonly accepted that the kind of bone anisotropy is the complex orthotropy, although can be approximated by the simplified transversely isotropic elasticity in some cases. In addition, the direction of maximum stiffness is primarily oriented in the longitudinal axis direction indicated by the shape configuration of most bone tissues, and the elastic modulus is greatest in this longitudinal orientation compared to radial and tangential (circumferential) orientations [9,10]. On the other hand, without an automatic and well defined assignment scheme, it is hard to generate the orthotropic principal axes vector fields (orientation of the principal symmetry axes of orthotropy), which change from point to point inside of bone tissues. And most work done in the physical modeling field adopted oversimplified isotropic material law due to its simplicity, which in fact should use orthotropic elasticity for bone tissues. Several authors incorporated simplified transversely isotropic or orthotropic elasticity into their FE models, but limited to only a segment of bone diaphysis, whose longitudinal axis is almost on a same straight line direction [15–18], which is convenient to define one global orthotropy orientation coordinate system. Some other authors tried to respect more real anatomical situation with extensive manual work. Based on the orthotropic properties of cortical bone measured by SchwartzDabney and Dechow [13], Natali et al. build one segment of mandible FE model used in dental biomechanics [19]. The material axis 1 of cortical bone is aligned along the distal–mesial direction, while axes 2 and 3 are defined by mutually orthogonal vectors and depend on the assumed location as reported in [13]. Apicella et al. extended this work to the entire mandible structure [20]. The cortical bone is considered as piecewise heterogeneous, and divided into 62 sites. Each site has its own thickness, orthotropic constant values and its own orientation of the maximum stiffness direction according to the proposed experimental mechanical characterization [13]. Some authors implemented invasive methods to manually assign the principal axes of orthotropy to FE models. Wirtz et al. [21] cut the cadaveric femurs in 2 mm slices, which were scanned into an image processing system. Then the principal

537

directions of stiffness of cancellous and compact bone were determined manually using the orientation of the trabecular structures and the Haversian system. Baca et al. [22] proposed to stain the femoral bone with red and green Indian ink depicting the direction of Haversian lamellar systems. Grinding was repeated several times with abrasive papers of increasing fineness. Repeated grinding exposed deeper layers and revealed canal networks which provided the picture of osteons spatial organization. These methods require a lot of time and manual effort, and obviously not suitable for clinical biomechanical analysis for patient-specific models. Recently, Liao et al. presented an ad hoc meshing tool to construct anisotropic model for patient-specific mandible, and employed projection and interpolation scheme to generate the trajectories of anisotropic material [23,24]. The approach is basically an interaction tool and the quality of the meshing resulting is largely determined by the experience of the user. In addition, it is only suitable for the out layer of the mandible, e.g., the compact bone, to build the orthotropic principal axes vector fields, while not suitable for the internal layer because of the limitation of projection and interpolation scheme. And simplified transversely isotropic elasticity is used for the internal cancellous bone. With intensive manual efforts, Korioth & Hannam built a nice mandible FE model, and the orthotropic material properties for cortical bone were assigned by use of a curve running along the entire lower jaw border and representing the long axis of the jaw [25]. Compared to the method by Liao et al. [23], using one curve seems too coarse to capture the mandibular geometry, especially for the posterior regions, like condylar process and ramus. The purpose of this study is to investigate a more convenient methodology for physical FE modeling of bone tissues, integrated with complete and continuous orthotropic material. The new approach is reliable and effective, and is almost automatic. The long term goal of our research is to provide a practical solution for the purpose of biomechanical simulations in the clinic.

3.

Design considerations

Our methodology is based on two basic observations: the first is that the longitudinal axis direction indicated by the shape configuration of most bone tissues is compatible with the trajectory of maximum material stiffness [9,10]; the second is that the scalar distribution pattern of harmonic field, which has scalar iso-contours and gradient vectors of high quality, tends to conform very well to the shape of object [26]. That is, although not by physical observations in a strict sense, the orthotropic principal axes can be determined utilizing the harmonic fields. Then, both meshing and orthotropic material assignment could be considered together to design a novel automatic mesh generator (AMG), if we can design appropriate harmonic fields to drive the modeling process. At a high level, our orthotropic FE modeling framework consists of the following basic steps: (1) reconstruct individual bone surface models from CT data; (2) construct surface harmonic fields of the longitudinal axis direction for bone models; (3) distribute iso-contours adaptively and construct volumetric meshes; (4) generate a longitudinal volumetric harmonic

538

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

field by extending the surface harmonic field to the volumetric mesh, from which the gradient vector field is calculated and employed as the principal axis field of maximum material stiffness (longitudinal orientation); subsequently, trace center flow curves of the longitudinal volumetric harmonic field to create a radial volumetric harmonic field, from which the radial and tangential (circumferential) principal axis fields are generated; (5) establish the final orthotropic physical FE model, set the inhomogeneous elastic parameters via the CT data if required. The detailed process of the method is discussed in following subsections. Then, two groups of contrastive FE analysis between the complete orthotropic model and comparing isotropic model are conducted to demonstrate the influence of the elastic orthotropy.

4.

Description of method

4.1.

Reconstruct surface model from CT

The starting point of our investigation is well reproduced individual geometry surfaces of various bone tissues, including the separation between compact and cancellous bone if necessary, such as the human mandible model illustrated in Fig. 1. These surfaces are usually reconstructed from CT data by segmentation and reconstruction algorithms such as the Discretizied Marching Cube algorithm. Each surface is represented as a closed manifold triangular mesh MS = (V,E,F), where V = {vi|1 ≤ i ≤ N} denotes the set of vertices, E = {(i, j)|1 ≤ i, j ≤ N, i = / j} denotes the set of edges, / j = / k} – the set of triangle and F = {(i, j, k)|1 ≤ i, j, k ≤ N, i = faces. Each vertex is assigned a position in a 3D Euclidean space. Because original surfaces contain numerous poorly shaped triangles, and are too large by several orders of magnitude to create FE model directly, a convenient post-processing is investigated to address the problems as described subsequently.

4.2.

Generate longitudinal surface harmonic field

To take advantage of the fact that the longitudinal axis direction indicated by the shape configuration of most bone tissues is compatible with the trajectory of maximum material stiffness, we need to design a longitudinal scalar field whose scalar distribution pattern tends to conform very well to the shape of surface. Then, from which the iso-contours can be extracted to re-mesh the surface, and the gradient vectors can be calculated as the principal orientation of maximum material stiffness. Specifically, we are interested in constructing a harmonic field f such that f = 0, where  is the Laplacian operator, subject to the Dirichlet boundary conditions that vertices in the set VC ∈ V of k ≥ 2 constrained vertices taking on prescribed values. The standard definition of the Laplacian operator on a piecewise linear mesh M is the umbrella operators fi =



wij (fi − fj ),

j ∈ N(i)

(1)

where j ∈ N(i) is the set of vertices adjacent to vertex i and wij is a scalar weight assigned to the edge (i,j). The standard choice for the weights wij of surface mesh are the discrete harmonic weights wij = 1/2(cot ˛ij + cot ˇij ), where ˛ij and ˇij are the angles opposite the edge [27]. Assembling the vertices’ function values fi into an n-vector f, the Laplacian operator can be written in matrix form LF = 0, where the elements of the n × n matrix L are given by

Lij =

⎧ ⎨ k ∈ N1 (i) wik , ⎩

if i = j

−wij , if j ∈ N(i) 0, otherwise

(2)

Eliminating rows and columns corresponding to the constraint vertices by bringing them to the right-hand side yields a linear system of the form Ax = b, consisting of positive definite sparse matrix A and a right-hand-side vector b. Iterative methods such as the preconditioned conjugate gradient method are very efficient to obtain the solution of the linear system. One of the primary properties of harmonic field is having no local extrema other than at constrained vertices. Furthermore, if all constrained minima are assigned the same global minimum value and all constrained maxima are assigned same maximum value, then all the constraints will be guaranteed to be extrema in the resulting field, and whose gradient vector field will converge precisely at the specified constraint points and flow smoothly everywhere else [21]. In addition, x is a discretization of the Laplace–Beltrami operator for triangular mesh, and approximates the mean curvature normal, so that the scalar distribution pattern of harmonic field tends to conform very well to the shape of surface. These characteristics are the very properties we are looking for to our application. In other words, we just need to place the minima and maxima constraints on the end tips of longitudinal axis direction on bone model to generate the longitudinal harmonic field to drive the meshing scheme. For the mandible, the direction of maximum stiffness is oriented along the long axis of the bone [9], although it is an irregular ‘curve’ long bone [28]. So the tips of condyle process and coronoid process of left and right sides are selected as the minima and maxima constraints respectively (the same process is applied for internal surface), such as the little red and blue balls on mandible ends in Fig. 2. The resulting harmonic field is illustrated in Fig. 2(a), mapping the scalar values to colors (red represents small value and blue represents large value), and the iso-contours are sampled by even scalar values. Part of the gradient vector field is illustrated in Fig. 2(b). We can see that: on one hand, although the underling surface has irregular mesh sampling (white mesh edges in Fig. 2(b)), the iso-contours and the gradient flow lines still follow the shape of the surface very well; on the other hand, the iso-contours sampled by even scalar values can not be used to construct new mesh as they distribute too densely on some regions (e.g., near the constraints points and the mandibular oblique line), whereas too sparsely on some other regions (e.g., near the mandibular notch and mandibular angle). Furthermore, current gradient vector field is only on surface, and need to be expanded to the volumetric domain to represent the principal axis field of maximum material stiffness.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

539

Fig. 1 – Individual geometry of human mandible.

4.3.

Create volumetric mesh

As discussed before, it is necessary to improve the distribution pattern of the iso-contours. Our solution is to first compute the shortest path between the minima and maxima constraints using the multi-source Dijkstra’s Algorithm, e.g., from the left coronoid process to the right coronoid process for mandible, then trace a surface gradient flow curve (use method by Dong et al. [26]) by using one path vertex as a seed point, such as the red flow curve in Fig. 3(a). Obviously, this gradient flow curve is also almost the shortest, and the scalar values on the curve increase rapidly. Thus, the algorithm selects seed points by even distance or any user define function, and gets the corresponding scalar values to generate iso-contours, so that there will not occur any too dense contour distribution. Subsequently, the local distance between adjacent iso-contours is checked, if the distance exceeds some threshold value, a middle seed point is generated to trace a new iso-

contour, which is not closed of course. This process amends those too sparse iso-contour distributions, such as shown in Fig. 3(a). Following, these iso-contours are re-sampled with new mesh key points, which are distributed adaptively considering factors of curvature, biomechanics significant region, and so on. Note that local surface harmonic scalar values of the iso-contours are assigned to corresponding key points. Then, these key points are added temporarily to the original surface as site centers to build the dual Voronoi cell graph, such as shown in Fig. 3(b), which is used to construct the new remeshing surface. Afterwards, the volumetric mesh is constructed by a two step workflow (only the first step if there is no internal separation surface). The first step is to build the volumetric mesh for the compact bone (the region enclosed by both the out and internal surface) from the previous re-meshing surfaces, by a Boundary Constrained Delaunay Triangulation proce-

Fig. 2 – (a) Longitudinal surface harmonic field and iso-contours sampled by even scalar values. (b) Local view of gradient vector field.

540

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

Fig. 3 – (a) Adaptively distribute iso-contours. (b) Generate new mesh key points and the dual Voronoi cell graph.

dure. And the algorithm employs a simple centroid point insertion strategy, if necessary, to ensure the quality of the produced tetrahedral elements. The second step is to construct the interior mesh for cancellous bone by an Advancing

front algorithm, with the internal boundary surface of the produced compact bone mesh as the initial active front. An example of resulting volumetric mesh for mandible is given in Fig. 4.

Fig. 4 – Resulting volumetric meshes for mandible compact bone (a and b) and cancellous bone (c and d). (a) and (c) Demonstrate boundary surface meshes, (b) and (d) depict shrinking tetrahedral elements.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

4.4.

Construct orthotropic principal axes vector fields

For now, we have tetrahedral mesh MV = (V, E, F, T), where V, E, F denote the set of vertices, edges, and triangles, and T = {(i, j, p, q)|1 ≤ i, j, p, q ≤ N, i = / j = / p = / q} denotes the set of tetrahedrons. As discussed before, we want to create volumetric harmonic fields to generate the orientation vector fields of the orthotropic material, which needs the discrete volumetric Laplacian operator. However, unlike the surface Laplacian operator, there are only a few published studies that describe the volumetric Laplacian operator [29]. The authors [30] recently deduced the volumetric Laplacian operator directly based on its definition f ≡ div(f), where  is the volumetric gradient operator and div(·) is the divergence operator, and set n the edge weight as wij = 1/6 k=1 lk cot(k ). It should be pointed out that when solving a Possion problem, e.g., volumetric mesh deformation [30], the exact coefficients of the Laplacian operator is important, or else leads to some unnatural shrinking/expanding effect. While, for solving a harmonic problem, the constant coefficients can be discarded. Using the form of formula (1), we can set the volumetric mesh edge weight as n wij = l cot(k ), and substitute into formula (2) to genk=1 k erate volumetric harmonic field, with appropriate boundary conditions. The first volumetric harmonic field is generated directly by extending the outside surface harmonic field to the volumetric mesh. Remember that those new mesh key points on outside surface record the local values of surface harmonic field, then all the corresponding boundary vertices of volumetric mesh are selected as constrained vertices and assigned these scalar values. In other words, the outside surface harmonic field is utilized as the Dirichlet boundary conditions to construct a longitudinal volumetric harmonic field. Fig. 5(a) illustrates the longitudinal volumetric harmonic field for the mandible. Note that its scalar distribution pattern still conforms very well to the model shape. In addition, the volumetric gradient vectors coincide with the trajectories of maximum stiffness, and is employed directly as the orthotropic longitudinal principal axis field, as shown in Fig. 5(b). Following, we need to create a radial volumetric harmonic field, which requires center axis curves of the volumetric domain, as well as the boundary surfaces to be used as Dirichlet boundary conditions. Instead of extracting the medial axis of the object, which is difficult and slow, the algorithm traces the center flow curves of the longitudinal volumetric harmonic field, which is simple and fast. And the vector flow tracing method on surface mesh by Dong et al. [26] is generalized to volumetric mesh in our case, such as shown in Fig. 5(c). The constrained minima or maxima vertices on internal separation surface are used as the flow seeds, which need no more interaction inputs. As the path points of flow curves usually locate on faces or edges of internal tetrahedrons other than vertices, a simple refinement scheme is employed to temporarily add them as mesh vertices. Then, with these flow vertices constrained with value 0, and outside mesh boundary vertices constrained with value 1, a radial volumetric harmonic field can be readily constructed. However, it should be noted that the radial volumetric gradient field does not guarantee to be fully orthogonal to the longitudinal gradient

541

field. So for every tetrahedron, the algorithm project its radial gradient vector onto a local plane with the longitudinal gradient vector as normal, to generate the final orthotropic radial principal axis field, such as these green vectors depicted in Fig. 5(d). At last, the tangential (circumferential) material orientation could be derived from the cross product of the longitudinal and radial orientations, such as these blue vectors depicted in Fig. 5(d).

4.5.

Programming environment and hardware

All these techniques have been implemented by Visual C++ 2008 (Microsoft Inc., Seattle, USA) as a modular software tool in the USIS system, running on a common Pentium IV personal computer.

5.

Status report

5.1.

Contrastive finite element analysis

To reveal the influence of elastic orthotropy for physical FE simulation, we give two groups of contrastive FE analysis for the mandible, embedding 2 ITI threaded implant CAD models of standard size in the posterior zone such as shown in Fig. 6(a), with assumed perfect bonding at all interface between the bone and implants. The implant model is made of titanium material with E = 110,000 MPa,  = 0.35. The first experiment is to investigate the ‘pure’ influence of orthotropy, using homogeneous material property. First, a fully orthotropic FE model based on above methodology, whose orthotropic elastic coefficients are listed in Table 1 [9–12], and a comparing FE model with effective isotropic elastic material, E = 16.5 GPa and  = 0.34 for compact bone, E = 0.482 GPa and  = 0.26 for cancellous bone, are constructed. For each model, the loading is simulated by applying an oblique load (vertical load of 200 N and horizontal load of 40 N) from buccal to lingual to each implant. And boundary constraint conditions include the joint surface of the condyles and the attachment regions of masticatory muscles, such as illustrated in Fig. 6(b). As Baca et al. recently reported that that realistic material property assignment (orthotropy) is very important for the FE analyses of small bone specimens, whereas in global FE analyses, this assignment can be omitted, if the inhomogeneous material model was used [22]. We thus design a second experiment to investigate the influence of orthotropy with inhomogeneous material property. For each mesh element, the effective density () is derived by a linear relationship from the average CT HU value [17]. Then, following functions are used to set the inhomogeneous elastic coefficients as Liao et al. [23]. For the isotropic model:

EC = 20653.09 , C = 0.3, ET = 19041.64 , T = 0.3,

(3)

542

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

Fig. 5 – (a) Longitudinal volumetric harmonic field. (b) Orthotropic longitudinal principal axis vector field. (c) Center flow curves of longitudinal field. (d) Orthotropic principal orientations of mandible: the longitudinal, radial, and tangential (circumferential) directions are depicted by red, green and blue vectors. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

Fig. 6 – (a) Simulation of implants embedding in mandible posterior zone. (b) Oblique loading and boundary constraint conditions.

543

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

Table 1 – Orthotropic elastic coefficients for compact (Comp.) and cancellous (Canc.) bone of human mandiblea .

Comp. Canc. a

E1

E2

E3

G12

G13

G23

12

13

12.5 0.511

17.9 0.114

26.6 0.907

4.5 0.078

5.3 0.434

7.1 0.081

0.18 0.22

0.31 0.31

Ei represents Young’s modulus (GPa); Gij represents shear modulus (GPa); ij represents Poisson’s ratio. The 1-direction is radial, the 2-direction is tangential (circumferential), and the 3-direction is axial (longitudinal).

where C and T represent compact and cancellous (trabecular) bone. The functions for orthotropic model are: EC 1

23 0.28 0.30

=

23141.57 ,

EC 2

=

23141.57 ,

EC 3

=

20653.09 ,

1.78 , ET = 11571.78 , ET = 19041.64 , ET 1 = 1157 2 3

G12 =

G12max 2 2 max

, G23 =

G23max 2 2 max

, G31 =

G31max 2 2 max

,

(4)

12 = 0.4, 23 = 0.25, 31 = 0.25, where G12max = 5.71 MPa, G23max = 7.11 MPa, G31max = 6.58 MPa, and max is the maximum density of all mesh elements, whose value depends on individual CT data. Then, same loading conditions are applied as the first experiment. It should be pointed out that these mapping functions are mainly for femur bone, but the objective here is only to check out whether in global FE analyses the assignment of orthotropy can be omitted or not if the inhomogeneous material model is used, thus the configuration here is appropriate for comparative studies. The results of homogeneous and inhomogeneous groups were processed separately. For the homogeneous configuration, the principal stresses and principal strains surrounding the implants, cortical and cancellous bone separately, of the orthotropic and isotropic models were compared. The result demonstrated that the orthotropy for cortical bone increased almost all of the principal stresses and principal strains compared to the isotropic case, the maximum increasing magnitude reached up to 68% for principal stresses, and 53% for principal strains. On the other hand, orthotropy decreased 1st principal tensile stress in cortical bone by 10% to 17%. For the cancellous bone, orthotropy also increased almost all of the principal stresses and principal strains compared to isotropic case, the maximum increasing magnitude reached to 160% for principal stresses and 60% for principal strains. Compared to the simplified transversely isotropic elasticity used in [19], the orthotropy results much higher stresses for the cancellous bone. In addition, there were many differences among the stress and strain distribution patterns between the orthotropic and isotropic cases, such as shown in Fig. 7. For the inhomogeneous configuration, the maximum values for displacement, principal stresses, principal strains and shear strains surrounding the implants of the two FE models were compared. The result demonstrated that the orthotropy increased the maximum value of displacement up to 2.5 times, and increased three principal stresses and principal strains up to 1.2 to 5.8 times, compared to the isotropic case, which absolutely could not be ignored. The maximum values for shear strains revealed less differences between the orthotropic and isotropic models, but there were relatively large differences among the distribution patterns like the principal stresses and strains, such as shown in Fig. 8.

5.2.

Methodology strengths and limitations

To the best of the authors’ knowledge, there is almost no published study that develops an AMG technique integrated with complete and continuous orthotropic material characteristics assignment function in the biomechanics and physical modeling fields. The results demonstrated that the volumetric meshes produced by our method provide good fit to the original complex bone geometries, with a geometric root mean square (RMS) error of less than 0.4 mm and peak errors less than 2.5 mm, and have nice mesh structures, such as shown in Fig. 4, which properly account for the preferential direction dictated by the shape configuration as well as the orthotropic material characteristics. The element shape quality measure, using the mean-ratio metric [31] (the value range varies from 0 to 1, attaining 1 for the ideal element shape), demonstrates that the resulting mesh usually has an overall average value of more than 0.85. A few poorly-shaped elements might exist in the interior (center region) of the mesh, with the minimum value closes to 0.45, because the Advancing front algorithm is employed for the cancellous bone region. But this has minimal influence on the analysis results in our numerical experiments. To assess the necessary level of mesh density, a group of meshes with different average element size from 0.5 mm to 5 mm were generated with local mesh density control (by adaptive iso-contour and key point distribution). It turns out that an average element size of 2 mm ensures that the model is already in the asymptotic region of the convergence and represents a good trade off between numerical accuracy and computational efficiency. Most importantly, the volumetric harmonic fields guarantee a high degree of smoothness of the resulting principal axes vector fields of orthotropic material. The computation is very stable and robust through out the whole object. In contrast, the previous projection and interpolation scheme is limited to outside compact bone region [23–25]. In addition, the whole process is automatic except only a few extreme constraint points input for the construction of surface harmonic field at the beginning. The orthotropic material property assignment to femur has been adopted in recent studies [21,22], using some invasive methods of manual cutting or grinding schemes. However, these methods not only require a lot of time and manual effort but also include error prone procedures. Fig. 9 shows an orthotropic femur FE model constructed by our methodology, which needs only about 1 min on a standard desktop computer with Intel P4 3.0 GHz processor and is almost automatic. We also acknowledge the limitations of our study. As a global and automatic method, it is very convenient to use and generate the whole orthotropic physical model. On the

544

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

Fig. 7 – Three principal strains in cancellous bone, using homogeneous material. Left column is for isotropic and right for orthotropic.

other hand, we notice that in some area, the orthotropic principal axes orientations may be not completely identical to those measure results by experimental methods. Of course, all those previous anisotropic modeling methods have similar problems. Yang et al. [18] pointed out that bone structure is customarily recognized as confirming to Wolff’s law, which is essentially the observation that bone changes its external shape and internal architecture in response to strains and stresses acting on it. The trabecular structure and haversian system strongly coincide with the principal stress track in femur mainly in order to bear tensile or compressive stresses and avoid being sheared. The principal stress track is consistent with the load transfer direction, from the femur head, along the direction of the neck and the stem, to the distal femur. Therefore, the directions of the femoral neck and stem can be approximately regarded as the principal material orientations in femur. In the femoral neck, the principal axis

has an angle 120◦ to the global Z-axis; in the femoral stem, the principal axis has an approximate angle 12◦ to Z-axis. For the entire mandible model of Apicella et al. [20], the cortical bone is divided into 62 sites. Each site has its own thickness, orthotropic constant values and its own orientation of the maximum stiffness direction. The discontinuous configuration between piecewise regions will inevitably lead to some numerical deviation. Another limitation is that, in our current FE model, some underlying kinematical conditions are simplified, e.g., not including the nonlinear visco-elatic behavior of the periodontal ligament for the mandible [32]. We intend to add these material assignment functions into the AMG software in the near future, which is a straightforward task. About the validation of the orthotropic models, Korioth & Hannam did validation tests for a freshly dissected mandible by comparing the predicted principal strains from orthotropic

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

545

Fig. 8 – Three shear strains (left three columns) and three principal strains (right three columns) in cancellous bone, using inhomogeneous material. First row (a–f) is for isotropic and second row (g–l) for orthotropic.

FE analysis with the strain readings measured experimentally, and the validation revealed small differences between observed and predicted principal strains [25]. While, our models are constructed from in vivo human CT data, a direct check of stress and strain at the simulated model itself is not available physically. Therefore, related activities are beyond the scope of current paper. Anyway, we have resorted to some clinical experts for qualitative evaluation, conformed that the longitudinal axis direction generated by our method is very close to the bone Haversian system, and believe that the orthotropic FE modeling methodology in this study brings us one important step closer toward realizing realistic biomechanics and physical simulation for the mechanobiology field.

6.

Lessons learned

Information on the orientation of material axes is especially important for studies of regional differences and questions relating bone function of bone adaptation and growth, as local anisotropic behavior and regional variations in bone tissues can have pronounced effects on the relationship between stress and strain patterns. Apicella et al. [20] tested two cubes with orthotropic materials, one cube has its maximum stiffness parallel to the applied load, whereas the other has maximum stiffness normal to the applied load. The analysis showed that even if the structures share similar stress values and distribution, the strain state is significantly different and

is governed by the orientation of the maximum stiffness. As the bone adaptation and remodeling based on the trajectorial theory of Wolff is strongly sensible to the computed strains component values [33], an error occurring in computing the strain state through FEA would also have errors in determining bone density or material properties, which would alter the algorithms reliability. To address this problem, this paper has presented an automatic methodology based on harmonic fields, to construct high quality volumetric FE mesh integrated with three inherent principal axes vector fields of orthotropic material. Contrastive FE analyses, including homogeneous and inhomogeneous configurations, demonstrate that elastic orthotropy has significant effect on simulating stresses and strains, including the value as well as distribution pattern, which underlines the relevance of our orthotropic modeling scheme. Note that in the work of Baca et al. [22], the single force loading was applied almost along the axis direction of the femur, which is also the direction of maximum stiffness for orthotropy. And the function for the Young’s modulus for the maximum stiffness direction is just the same as the function for the Young’s modulus of the isotropic model. Thus, it is understandable that the maximum displacement values compared did not demonstrate large difference. In addition, they only compared the maximum displacement values, while, examinations of strains and stresses, especially the distribution patterns, are often considered to be more important, as they are directly related to the functional adaptation of bone

546

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

Fig. 9 – (a) Femur volumetric mesh. (b) Longitudinal volumetric harmonic field and center flow curves. (c and d) Three orthotropic principal axes vector fields.

tissue. More recently, Yang et al. [18] investigated differences between isotropic and orthotropic inhomogeneous material models of femur. The results showed that great differences exist in many femur regions. Apicella et al. [20] also pointed out the mandible structure is significantly sensible to compact bone orthotropy and thickness in the condylar neck, anterior margin of ramus, inferior portion of anterior margin of the ramus, retro molar area on facial side, middle portion of the corpus in molars area, and anterior margin of the ramus on lingual side.

7.

Future plans

To solve the problem that in some area, the orthotropic principal axes orientations may be not completely identical to those measure results by experimental methods, we are planning to add a control function that local axe orientation constraints could be set in some area by the user if necessary, e.g., according to those experimental measurements. These orientation vectors are to be used as some constraint conditions when calculating the harmonic fields, thus can provide much better consistency with the experimental measurements. In addi-

tion, currently limited orthotropic FE analysis was conducted in this paper. Further detailed researches are surely needed to study the influent of orthotropic material on the mechanical function and biological adaptation of bone tissue.

Acknowledgements The research was supported by the National natural science foundations (No.60903136, No.60970098), Special postdoctoral science foundation of China (No.201003723), Foreign cooperation project of Wenzhou technology division (H20080026) and Wenzhou medical college (YXYZD5).

references

[1] J.W. Fernandez, P.J. Hunter, An anatomically based patient-specific finite element model of patella articulation: towards a diagnostic tool, Biomechanics and Modeling in Mechanobiology 4 (2005) 20–38. [2] M. Keyak, S. Rossi, K. Jones, H. Skinner, Prediction of femoral fracture load using automated finite element modeling, Journal of Biomechanics 31 (1998) 125–133.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 536–547

[3] M. Viceconti, L. Bellingeri, L. Cristofolini, A comparative study on different methods of automatic mesh generation of human femurs, Medical Engineering & Physics 20 (1998) 1–10. [4] J. Gao, W. Xu, Z. Ding, 3D finite element mesh generation of complicated tooth model based on CT slices, Computer Methods and Programs in Biomedicine 82 (2006) 97–105. [5] Y. Zhang, C. Bajaj, Adaptive and quality quadrilateral/hexahedral meshing from volumetric data, Computer Methods in Applied Mechanics and Engineering 95 (2006) 942–960. [6] S.C. Cowin, Bone Mechanics Handbook, Second ed., CRC Press, Boca Raton, FL, 2003. [7] D.C. Wirtz, N. Schiffers, T. Pandorf, K. Radermacher, Critical evaluation of known bone material properties to realize anisotropic FE-simulation of the proximal femur, Journal of Biomechanics 33 (2000) 1230–1325. [8] F. Taddei, A. Pancanti, M. Viceconti, An improved method for the automatic mapping of computed tomography numbers onto finite element models, Medical Engineering & Physics 26 (2004) 61–69. [9] P.C. Dechow, D.H. Chung, M. Bolouri, Relationship between three-dimensional microstructure and elastic properties of cortical bone in the human mandible and femur, in: Primate Craniofacial Function and Biology, New York: Springer US Press, 2008, pp. 265–292. [10] P.C. Dechow, W.L. Hylander, Elastic properties and masticatory bone stress in the macaque mandible, American Journal of Physical Anthropology 112 (2000) 553–574. [11] A.M. O’mahony, J.L. Williams, J.O. Katz, P. Spencer, Anisotropic elastic properties of cancellous bone from a human edentulous mandible, Clinical Oral Implants Research 11 (2000) 415–421. [12] C.L. Schwartz-Dabney, P.C. Dechow, Edentulation alters material properties of cortical bone in the human mandible, Journal of Dental Research 81 (2002) 613–617. [13] C.L. Schwartz-dabney, P.C. Dechow, Variations in cortical material properties throughout the human dentate mandible, American Journal of Physical Anthropology 120 (2003) 252–277. [14] Z. Fan, J.G. Swadener, J.Y. Rho, M.E. Roy, G.M. Pharr, Anisotropic properties of human tibial cortical bone as measured by nanoindentation, Journal of Orthopaedic Research 20 (2002) 806–810. [15] A.M. O’mahony, J.L. Williams, P. Spencer, Anisotropic elasticity of cortical and cancellous bone in the posterior mandible increases peri-implant stress and strain under oblique loading, Clinical Oral Implants Research 12 (2001) 648–657. [16] H.L. Huang, C.L. Lin, C.C. Ko, C.H. Chang, Stress analysis of implant-supported partial prostheses in anisotropic mandibular bone: in-line versus offset placements of implants, Journal of Oral Rehabilitation 33 (2006) 501–508. [17] L. Peng, J. Bai, X. Zeng, Y. Zhou, Comparison of isotropic and orthotropic material property assignments on femoral finite element models under two loading conditions, Medical Engineering & Physics 28 (2006) 227–233. [18] H. Yang, X. Ma, T. Guo, Some factors that affect the comparison between isotropic and orthotropic

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

547

inhomogeneous finite element material models of femur, Medical Engineering & Physics 32 (2010) 553–560. A.N. Natali, E.L. Carniel, P.G. Pavan, Investigation of bone inelastic response in interaction phenomena with dental implants, Dental Materials 24 (2008) 561–569. D. Apicella, R. Aversa, F. Ferro, D. Ianniello, A. Apicella, The importance of cortical bone orthotropicity, maximum stiffness direction and thickness on the reliability of mandible numerical models, Journal of Biomedical Materials Research Part B: Applied Biomaterials 93B (2010) 150–163. D.C Wirtz, T. Pandorf, F. Portheine, K. Radermacher, N. Schiffers, Concept and development of an orthotropic FE model of the proximal femur, Journal of Biomechanics 36 (2003) 289–293. V. Baca, Z. Horak, P. Mikulenka, V. Dzupa, Comparison of an inhomogeneous orthotropic and isotropic material models used for FE analyses, Medical Engineering & Physics 30 (2008) 924–930. S.-H. Liao, R.-F. Tong, J.-X. Dong, Anisotropic finite element modeling for patient specific mandible, Computer Methods and Programs in Biomedicine 88 (2007) 197–209. S.-H. Liao, R.-F. Tong, J.-X. Dong, Influence of anisotropy on peri-implant stress and strain in complete mandible model from CT, Computerized Medical Imaging and Graphics 32 (2008) 53–60. T.W. Korioth, A.G. Hannam, Deformation of the human mandible during simulated tooth clenching, Journal of Dental Research 73 (1994) 56–66. S. Dong, S. Kircher, M. Garland, Harmonic functions for quadrilateral remeshing of arbitrary manifolds, Computer Aided Geometric Design 22 (2005) 392–423. U. Pinkall, K. Polthier, Computing discrete minimal surfaces and their conjugates, Experimental Mathematics 2 (1993) 15–36. R.B. Ashman, W.C. Van Buskirk, The elastic properties of a human mandible, Advances in Dental Research 1 (1987) 64–67. M. Meyer, M. Desbrun, P. Schroder, A.H. Barr, Discrete differential-geometry operators for triangulated 2-manifolds, in: Proceedings of VisMath’02, Berlin, 2002, pp. 35–57. S.-H. Liao, R.-F. Tong, J.-X. Dong, Gradient field based inhomogeneous volumetric mesh deformation for maxillofacial surgery simulation, Computers & Graphics 33 (2009) 424–432. J. Dompierre, P. Labbe, F. Guibault, R. Camarero, Proposal of benchmarks for 3d unstructured tetrahedral mesh optimization, in: Proceedings of 7th International Meshing Roundtable, Dearborn, 1998, pp. 459–478. R. Sorrentino, D. Apicella, C. Riccio, E. Gherlone, et al., Nonlinear visco-elastic finite element analysis of different porcelain veneers configuration, Journal of Biomedical Materials Research Part B: Applied Biomaterials 91B (2009) 727–736. R. Aversa, D. Apicella, L. Perillo, R. Sorrentino, et al., Non-linear elastic three-dimensional finite element analysis on the effect of endocrown material rigidity on alveolar bone remodeling process, Dental Materials 25 (2009) 678–690.