FE model of intact lumbar spine (Intact model)

To create a three-dimensional FE model, computed tomography scan DICOM files of the L1 to L5 lumbar spine of a middle-aged male were obtained at 1-mm intervals. The commercially available visualization software Amira 3.1.1 (Mercury Computer Systems, Inc., Berlin, Germany) was used to describe cross-section contours of each spinal component in accordance with gray scale value (Figure 3.1). Then, the three-dimensional surface geometries were constructed through sequential processed cross-section contours as shown in Figure 3.2 A. Each spinal component was exported as a Drawing eXchange Format (DXF) file and converted to the Initial Graphics Exchange Specification (IGES) file as shown in Figure 3.2 B. The FE analysis software ANSYS 9.0 (ANSYS Inc., Canonsburg, PA) was used to reconstruct the FE model by converting the IGES file to ANSYS Parametric Design Language (APDL) code in Figure 3.2 C. The INT model was an osseo-ligamentous lumbar spine, which included the vertebrae, intervertebral discs, endplates, posterior bony elements, and all seven ligaments (Figure 3.3 A).

An eight-node solid element (SOLID185) was used for modeling the cortical bone, cancellous bone, posterior bony element, cartilage endplate, and annulus ground substance.

The cortical bone and cancellous bone were assumed to be homogeneous and transversely

isotropic [81]. The posterior bony element and cartilage endplate were assumed to be homogeneous and isotropic [82]. The intervertebral disc consisted of annulus ground substance, nucleus pulposus and collagen fibers embedded in the ground substance. The nonlinear annulus ground substance was simulated by using a hyper-elastic Mooney-Rivlin formulation [83][84]. The collagen fibers simply connected between nodes on adjacent endplates to create an irregular criss-cross configuration. These irregular angles of collagen fibers were oriented within the range of the Marchand’s [85] study. In the radial direction, twelve double cross-linked fiber layers were defined to decrease elastic strength proportionally from the outermost layer to the innermost. Therefore, the collagen fibers in different annulus layers were weighted (elastic modulus at the outermost layers 1-3: 1.0, layers 4-6: 0.9, layers 7-9: 0.75, and at the innermost layers 10-12: 0.65; cross sectional areas at the outermost layers 1-3: 1.0, layers 4-6: 0.78, layers 7-9: 0.62, and at the innermost layers 10-12: 0.47) based on previous studies [86][87]. The nucleus pulposus was modeled as an incompressible fluid with a bulk modulus of 1666.7 MPa by eight-node fluid elements (FLUID80) [81]. The 43 % of the cross-sectional area in the disc was defined as the nucleus, which was within the range of the study by Panagiotacopulos (30-50 %) [88] Therefore, approximately 47 % to 49 % disc volume was assigned to nucleus pulposus. All seven ligaments and collagen fibers were simulated by two-node bilinear link elements (LINK10) with uniaxial tension resistance only, which were arranged in an anatomically correct direction [89]. The cross-sectional area of each ligament was obtained from previous studies [82][87][90][91], and material properties of the spine are listed in Table 3.1. The facet joint was treated as having sliding contact behavior using three-dimensional eight-node surface-to-surface contact elements (CONTA174), which may slide between three-dimensional target elements (TARGE170). The coefficient of friction was set at 0.1[92].

The initial gap between a pair of facet surfaces was kept within 0.5 mm as shown in Figure 3.3 (b) [81]. The stiffness of the spinal structure changes depending on the contact status, so

the s NT model co

ach spinal c

on in ANS study. In ad changing th tion in AN onsisted of

component w file to crea

32

adopted to element’s ual element

chosen to ments and 9

d from com l-related con

account fo shape will t stiffness.

solve this 4,162 node

mputed tomo ntours.

or the chan change afte

Fi re

igure 3.2: M econstructed surface g

Modeling pr d through se geometry wa

rocess of the equential pr as exported

e L3 vertebr rocessed co d to the DXF

ra: (A) surf mputed tom F file; (C) F

(A)

(B)

(C) face geometr mography sc FE model of

tries of verte can DICOM f the L3 vert

ebra were M file; (B)

tebra.

Figgure 3.3: Finnite elemen trans

nt model of t sverse view

34

the L1 to L ws of facet jo

5 segments oint curvatu

(B) is shown: ( ure and gap.

Facet jo

(A) intact m oint 

(A)

model; (B)

Table 3.1: Material properties used in the FE model

Material Element type Young’s

modulus

Cortical 8node-Solid 185 Ex=11300 Ey=11300 Cancellous 8node-Solid 185 Ex=140

Ey=140

Nucleus pulposus Annulus Ground substance

8node-Fluid 80 8node-Solid 185

1666.7

Ligaments^* ALL

2node-Link 10

7.8 - 24

*ALL, anterior longitudinal ligament; PLL, posterior longitudinal ligament; TL, transverse ligament; LF, ligamentum flavum; ISL, interspinous ligament; SSL, supraspinous ligament;

CL, capsular ligament.

36 

In order to get reliable data, convergence test were conducted. Three mesh densities (coarse model: 4,750 elements / 4,960 nodes; normal model: 27,244 elements / 30,630 nodes;

finest model: 112,174 elements / 94,162 nodes) were selected to test ROM in the intact model (Figure 3.4). The boundary and loading conditions of the test were that the inferior surface of L5 vertebra was fixed, and 10 N-m moment and a 150 N preload were applied to the superior surface of L1 vertebra.

Compared with normal model and finest model, the variation of ROM was within 1.03%

in flexion (less than 0.2^o), 4.39% in extension (less than 0.5^o), 0.01% in axial rotation (less than 0.2^o), and 0.001% in lateral bending (less than 0.1^o). From the simulation results, the normal model only required fewer computational times to complete. However, several contact surfaces in facet joint have stress concentration owed to the lower smooth geometry for fewer elements and nodes. Therefore, the finest mesh density was selected in this study.

Figu res re

ure 3.4: Con sult of motio esult of mot

nvergence te on changes tion change

est of the in under flexi s under axia

(B)

(D) ntact model:

ion; (C) resu al rotation;

bendin

: (A) three m ult of motio

(E) result o ng.

mesh densit on changes u

f motion ch

(

ties were se under exten hanges unde

(A)

(C)

(E) lected; (B) nsion; (D)

er lateral

38