- Email: [email protected]
Clarification of the mechanical behaviour of spinal motion segments through a three-dimensional poroelastic mixed finite element model
1350-4533(95)00027-5
Clarification of the mechanical behaviour of spinal motion segments through a threedimensional poroelastic mixed finite element model J.S.S. Wu and J.H.
Chen*
Department of Mechanical Engineering, National Chung-Hsing University, Taichung, Taiwan 402, ROC; *Department of Applied Physics, Chung-Yuen Christian University, Chung-Li, Taiwan, ROC
ABSTRACT The purpose of this study is to c1nr-i~ the mechanical hehn-iliour (d ,spinul motion .tegments through a proper num&cal model. The model constructed ran give correct irtformntion and provide medical fields with valuble Lguidnncr irr .solving clinical problems orcurrkg in the sfiine. A three-dimensional pornelastic finite element model of sfiinnl motion segments is r.onstructed and a mixed formulation is introdured. The
Keywords:
Spinal
segment,
Med.
Phys.,
1996,
Eng.
poroelastic Vol.
18, 215-221,
models,
mixed
elements.
JD FEM
April
INTRODUCTION The spine is a biomechanical structure composed of 24 articulated columns of vertebrae; each of which is able to move through a joint of intervertebral disc (MI). Each IVD is a symphysis between two vertebral bodies (VBs) and formed with three parts; the nucleus pulposus, the annulus fibrosis and the hyaline cartilage end plates (EPs) on the opposing surfaces of each VB. The posterior facet joint (FJ) is of the diarthrodial type; it serves to guide and steady the movements of the VBs on one another. Each pair of FJs contacts with the upper and the lower FJs accordingly. Ligaments surround both the VB and the IVTI. The mechanism of a spinal motion segment (SMS) is defined as the composition of the vertebral column, the intervertebral joint, the FJs and ligaments. Changes from normal mechanical behaviour in the SMS may cause the pathogenesis Correspondence ical Engineering, wan 402, ROC.
tinite
to: Professor JSS Wu, PhD, Department of MechanNational Chung-Hsing University, Taichung, Tai-
of herniation and degeneration of the IVD, EP fracture, spinal deformity, low back pain syndromes, etc. Thus, understanding of mechanical behaviour in SMSs is quite important in relation to clinical problems occurring in the spine. Schultz et aZ.’ discussed mechanical properties of mathematical models of three-dimensional SMSs. A review of the mechanics of lumbar discs, towards a better understanding of low back pain, was made by Nachemson in 1975”. The fundamentals of the biomechanics of lumbar spine and its clinical significance were described by Koreska et ab” A numerical model of an intervertebral disc was developed by Belytschko et al. in 19744 for a three-dimensional stress analysis where the vertebrae were idealized as rigid bodies, while the discs, ligaments, and connective tissues were represented by deformable elements. Berkson et uZ.~ in 1979 studied the mechanical behaviour of 42 fresh human cadaver lumbar segments in compression, and in anterior, posterior, and lateral shear. Hakim and King’ used three-dimensional finite element models (FEMs) to analyse the
Spinal
segment
proelastic
mixed model: J.S.S. Wu and J H. Chen
dynamic response of vertebra, with experimental verification. Many investigations have been presented using poroelastic models to study the mechanical behaviour of tissues recently’*8. All components in the SMS contain different amounts of fluid in the pores of their solid tissues and bones. Solid models do not provide characteristics such as compressibility of the disc, long term creep response, etc. in the study of the mechanical behaviour of the porous structures of SMSs. Assuming each component of the SMS is porous, can provide a physical basis for determining the overall mechanical behaviour of the SMS under various loadings. Mathematical simulation of the SMSs using poroelastic axi-symmetric models can provide further insight ’ - l4 . However, since axi-symmetric models are difficult to apply to the complete geometry of SMSs, e.g., the non-circular transverse section of the SMS and the FJ part, correct mechanical behaviour of the objects is always misrepresented. Generally, the FJ exhibits strong support for the SMS that may affect the overall structural behaviour of the SMS significantly. For these reasons, a three-dimensional poroelastic FEM of a representative SMS, say, lumbar 5, is made to study the long term creep behaviour of the SMS here. The study can give correct information of mechanical behaviour in the SMS and provide medical areas valuable guidance in solving clinical problems occurring in the spine. The derivation of the constitutive relations based on energy laws is first briefly discussed and physical parameters in the equations of motion of porous media are interpreted. Then, the spatial discretization of the U-WV mixed finite element procedure is introduced. (Notations are defined in the Appendix.) Iteration in the computing process is performed in order to meet the specific boundary requirements on the contact surface of FJs and ligaments. In the final part, results such as the deformed shape, the fluid flow field and the failure stress distribution, are given and discussed. DEFINITION
OF DISPLACEMENTS
A general theory of threedimensional consolidation problems was discussed by Biot in 19’7’1’5 where the constitutive relations of porous media were discussed. An alternate form of the constitutive relations of porous media was also discussed by Zienkiewicz and Shiomi”. Here, the constitutive relations of porous media are established based on the energy laws. The poroelastic structures are composed of a ‘solid skeleton’ and an ‘interstitial fluid’ where the voids of the solid skeleton are saturated with fluid. The displacement field of the solid skeleton is defined by a column vector u, u = {ur u, ~3) where subscripts 1, 2 and 3 indicate parameters on the Cartesian coordinates in X-, p and zdirections, respectively. The displacement of fluid is defined as U and U = {U, U, U3}. Thus, the fluid volume displaced through pores perpendicular to the x,-axis is equal to nU, where n is called the porosity, defined as the ratio of the fluid volume
216
to the volume of total medium (n = V,/ V,) . A more significant vector, w and w = {w, w, w3} is defined as the displacement field of the pore fluid relative to the solid skeleton; wi( i = 1,2,3) is interpreted as the amount of fluid displaced through pores of the total area perpendicular to the i-th direction. So, the actual mean displacement of fluid through pores is equal to U-u and U-u = (l/n)w. STRESS-STRAIN
RJZATIONS
In the bi-phasic structure, strains are generated in both the solid skeleton and the fluid portions. The strain of the solid skeleton is expressed by e with a linear strain-displacement form of e = L u. Here, e is a column vector specified by {e,, q, es, %, e13 e,,}; L is the general strain-displacement operator in solid mechanics”. The fluid volumetric strain is given by 5, defined by the change of the fluid volume in the bulk volume of the medium, and 5 = VTw where V = LTm = @/ax, a/a% a/ax,1 and m = (1 1 10 0 O}. Stresses at a point in porous media are generated in both the solid and the fluid phases. Such stresses are represented by the total stress u and the hydrostatic pore fluid pressure 7~. The total stress u, is defined as forces acting on the bulk area of the medium, and the pore fluid pressure, n, is defined as forces acting on the fluid area. An extra parameter, cy, called compressibility of the ‘solid matrix’ (a solid including the effect from fluid), is equal to l-(&/K,), where & is the bulk modulus of the solid skeleton and K, is the bulk modulus of the solid matrix. The total stress u, accounting for total forces acting on the overall medium, is the sum of stresses in the solid skeleton (a’) and in the fluid (n). That is u=u’+amn,
(1) stress-strain
where a’ = De and D is the general matrix in solid mechanics”. The fluid strain is defined as n/Q where Q is the bulk modulus of fluid. For the hydrostatic case considered here, Q can be a scalar quantity onlvg*“. The fluid strain comes from two sources: (1)’ the fluid volumetric strain, 5, and (2) the excessive fluid strain, amTe, and is given by rr=a Substitution yields
Qm’e+
of equation
Qt (2) into
u=(D+a’@nmT)e+QJ.
(2)
equation
(1) (3)
Equations (2) and (3) are the general stress-strain relations of porous media which may be used to formulate the boundary value problem of such structures accordingly. ENERGY
THEORY
The total otential energy, II, of a structure is defined as R = E+ W, where E is the strain energy and W, is the potential energy of the applied loads. The external work, W, equals -W,. The principle of stationary potential energy gives
.Spnnl
6n=6E+6Wp=i3E-6W=0.
(4)
Generally, the study of mechanical behaviour of the SMSs involves consideration of the motion of stress waves that is time-dependent. Such timecreep and transient dependent cases, including cases depending on the applied loading rate, may raise the responding stress levels significantly. However, for the case of quite slow loading rate, the so-called long term creep case, or the static case for elastic structures over a long time interval, deformed rates may be neglected”. Here, the long term creep case only is introduced. Therefore, the derivation of the stiffness matrix [r;rl is provided. The variation of strain energy in the overall system, 6E:,, is derived as
vgmmt
proelnctic
tions (2), (3)) (5) and we obtain approach, manipulation:
Vn=
(l/k)*+
p&J;
/-fiindV=
jS(VTw)TmlI'.
(7) and the
using matrix
1%‘~ and ,I. H. Chev
the Galerkin form, after
F: 0 FZ
i I
(8)
where:
Y,,,=I B! DB,,dV &., = B,‘;cYN,dV I = N;B,,,dV
(61
(7)
Substituting equations (5) and (7) into equation (4)) two sets of equilibrium equations for the bi-phasic structures can be derived. Using the expression of stresses in equations (2) and (3) for u and 7~, and discretizing them into the proper form based on pertinent variables, the problem can then be solved. Details concerning the derivation of dynamic and kinetic relations of orous media can be found in the report by Wu’ R).
u-WV MIXED
[.,\.S.
con-
where pf is the fluid density and the dot specifies the time rate. The strain energy in fluid can be written in the form
SE,=
motirl:
.--
&, Darcy’s law when accounting for the general stant ‘permeability’, k, in fluid is”
mixrd
and B,, (= L) and B,v (= L m” = V) are the straindisplacement matrices of ZLand ZLI. Equation (8) can be solved and must be used when the solid skeleton and fluid are both incompressible, i.e., ~~41, Q+m and I<71*= 0. When Q is infinite, equation (8) cannot be solved by the Gaussian elimination procedure due to the singularity. In such case, reorder of parameters may be needed resulting in a non-symmetric form. Such non-symmetric fbrm is not encouraged in FEM analysis. In addition, when Qis finite for the compressible case, vector n* may be eliminated. Thus, we must have the equation. from the second row of equation (8), as m- = K$ (K&q* + K,,w* ) Substituting equations the or’ term, we have:
FORMULATION
Now, the Galerkin approachI is applied to derive proper formulas. Many alternative forms and solutions based on different parameters of porous media, have been studied and demonstrated”-” for axi-symmetric poroelastic problems. In a general finite element process, stresses are postinterpolated on the Gaussian points from the nodal displacements. Such post-interpolation procedures fit the resulting stresses poorly”. In addition, the surface condition of the pore fluid pressure in SMSs is quite important for these interactive, bi-phasic structures. Thus, a threedimensional mixed form, say u-WTT, is recommended here. This strategy may not only allow imposition and derivation of the pore fluid pressure on finite element nodes directly, but also include the effect of interaction between the pressure and displacements completely’H. In order to derive the finite element mixed form, we use the spatial semi-discretization with interpolation functions such as u = N,,u*, w = N,w*, rr = N,rr*, etc. where superscript * indicates finite element nodal values. Applying equa-
(9)
(9) and (8) and eliminating
(10) where:
Equation (10) is an alternative form of u-7~7~formulation where the benefit of the direct imposition and derivation of the pore fluid pressure can be obtained on finite element nodes. Equations (9) and (10) are both symmetric; thus, a partition frontal algorithm is implemented to carry out solutions for the large porous structure of SMSs. Once the displacements and fluid pressures are derived, total stresses may then be obtained by equation (3). (The total stress is interpolated on the Gaussian points first and then extrapolated to the elemental nodes by a specific designed algorithmlH.) The stress distributed in the solid skeleton is mostl) concerned in consideration of the failure of the poroelastic structure. The solid skeleton stress, u’, called the effective stress, is dcfincd as”’
217
Spinal
segment poroelastic
mixed model: J.S.S. Wu and J. H. Chen
u’=u-am7r
(13)
where
In this study, the von Mises criterion as a measure of failure of materials is used to quantify stresscomponents to arrive at one stress level on each node. The von Mises criterion applied to the solid skeleton may be expressed asI’: u;u
=
J
$(a:
- cq
+ $(a;, - 0:)’
+ t(a:
- u;,*
+ 375 + 3T; + 3TZI (14)
where a:,is called the effective von Mises stress, used to estimate the stresslevels generated in the SMSs in the following demonstrations. FINITE
ELEMENT
MODEL
Usually, each SMS contains three major parts: the VB, the IVD and the FJ. The VB is mainly constructed by an articulated column of the cancellous bone surrounded by a relatively stiff cortical shell. In the IVD, the pulposus, hydrated nucleus is centrally located and is relatively soft compared to the surrounding part of the fibrosus annulus. A porous, thin but rather stiff cartilage part, the EP, is located between the VB and the IVD. A pair of FJs is located at the posterior side of each VB. The basic function of the FJs is to guide and to steady the movements of each segment. In addition, there are ligaments surrounding the VB and the IVD. Here, the geometry of the SMS model is constructed from CT scanning images sectioned with the uniform space of 1 mm on a lumbar 5 SMS. This L5 SMS specimen was taken from a healthy man aged 46 with a height of 1.75 m and weight 160 pounds. The FEM of the L5 SMS is generated with brick elements and its degenerated elements automatically. Only one half of the height of IVD is included in the model for consideration of the opposing motion of the model at the mid-section of the IVD’,‘. Figure 1 shows the FEM and the reconstructed solid model. The definition of the height of each trans-
Figure
218
1
The finite
element
model
and the solid model
of the L5 spinal
verse section is shown in Fz@re 2. The above two figures will be used to calculate the deformed shapes and to show the stress distribution in the model. The L5 SMS model studied here is porous and elastic. Thus, material properties for each component include those in the solid skeleton, fluid and the solid matrix, see equation (8). The percentages of fluid contents (porosity n) in the cancellous bone, cortical shell, the EP and the FJs are relatively low, and may be ignored when compared to that in the MS~‘“,‘4. The porosity in the annulus and nucleus are both taken to be 0.70. All components including cancellous bone, cortical shell, EP, nucleus, annulus and FJs, have the same compressibility (Y taken as 0.9. The rest of the poroelastic material properties, such as the bulk modulus of the fluid Q are calculated from the fluid content and compressibility of the solid matrix”. Young’s modulus for each low fluid content is taken from the report by Belytschko in 19784. All other Young’s moduli are derived by a modified nonlinear optimization technique on experimental results”,18. Ligaments are simulated by an axial linear spring element with a constant of 15 N/mm correctly applied to nodes surrounding the VB and the IVD”~‘“. FJs take the same as for the cortical shell. material properties Material properties for each component are listed in Table 1. It is generally known that the human spine may transfer loads from’ the head to the trunk and to
Figure
motion
2
Definition
segment
of the section
height
-
30 26 22 18
z
;;t
50 46 42
38 34
Table 1 Material properties of the poroelastic of the 15 spinal motion segment Constituent
a
Q N/mm”
finite
I:’ N/mm”
element
n
model
v
flow and restricted to move horizontally. Each node on the contact surfaces of FJs is properly selected and its slope is estimated by interpolating the position of its neighbouring nodes in space. In the current FEM model, 1362 elements and 2144 nodes are built up and total degree of freedom is over 15,000. So, a specific frontal width optimization technique has been studied and implemented before solving the problem by the partition frontal solver here”. ITERATIONS
the pelvis which allow physiological motion among these three parts. So, the model set up here, does not relate to any postures of the human body but only considers forces transmitted fromthe earlier to the later segments. The upper surface of the VB is made to be the main region to receive loads. In this study, four loading cases are simulated: full compression, anterior bend, posterior bend and lateral bend. These four loading cases may be explained as the possible actions on the lumbar SMS caused by body motions. A uniform pressure of 0.4309 N/mm2* is applied to the upper surface of the VB’.“‘, see F@~re?. (This loading pressure has been used for the study of the mechanical behaviour of axi-symmetric models of SMS by the first author before. This report will compare some results with the axi-symmetric models.) Boundary conditions imposed on the structure are: (1) Nodes on the contact surface of facet joints are assumed to be inclined, (2) the bottom of the IVD in the FEM is impervious to the fluid
(a) Full
(c) Anterior
Compression
Bend
FOR
SPECIFIC
STRUCTURES
In computation, inclined boundary conditions of any compressive nodes on the contact surface of FJs are imposed and the element of ligaments are put on any tensile nodes located on the outer surface of the VB and the IVD. The strategy of accounting such restrictions in the computational procedure is explained as follows. In the initial step of computation, all nodes on the contact surface of I?Js are assumed to be compressive (inclined boundaries are imposed), and ligaments are assumed not appear. After the computation, if any nodal stresses on the FJ contact surface and the outer surface of the VB and the IVD are examined with different phases (compression or tension) from the latest assumption, a seconda? computation is thus necessary. In the second computation, adjustment of the boundary conditions is first considered. That is, (1) inclined boundary conditions of any compressive nodes on the contact surface of the FJs are imposed or otherwise released, and (2) the ligament elements are put on anv tensile nodes on the outer surface of the
(b) Posterior
Bend
(d) Lateral Bend
219
Spinal
segment poroelastic
mixed model: J.S.S. Wu andJH.
Chen
VB and the IVD or otherwise are taken away. Nodal stresses on the contact surface and the outer surface of the VB and the IVD are checked after the second computation. Such iterative computation continues until the nodal stress phase (compression or tension) on the contact surface of FJ and the outer surface of the VB and the IVD are maintained. RESULTS
AND
DISCUSSIONS
The mechanical behaviour of an L5 SMS is studied through a three-dimensional poroelastic mixed FEM. Structural analysis of the SMS considers certain displacements and stress levels as a measure of the failure of materials under specific loading and boundary conditions. Here, only the long time creep response is demonstrated. Figure 3 shows the deformed solid shapes of the model for cases of (a) full compression, (b) posterior bend, (c) anterior bend and (d) lateral bend. It clearly shows that the deformed shape of each component in the SMS is never symmetric with respect to the central line of either the VB or the IVD. Thus, the mechanical behaviour of the SMS is difficult to interpret by an axi-symmetric structural model. There is no doubt that a three-dimensional model of the SMSs is valuable in providing clinical information on the spine. Table 2 lists the deformed sizes of the model. The disc bulge and vertebral vertical deflection calculated by the axi-symmetric model under full compressive load are 0.254 mm. and 0.385 mm., respectively’,‘; so, the results shown here are physically significant. The largest disc bulge and vertebral vertical deflection appear in the fully compressive case. The smallest posterior and anterior bulges and vertebral vertical deflection appear in the lateral bending case, no obvious bulging is found for this case. The above phenomena are caused by the strong support from the FJs; so, the disc bulge, the EP and the vertebral deflections are obviously not able to be read from axi-symmetric models. Since the fluid content in the IVD is relatively high compared to others, the motion of fluid in the IVD is discussed here. Also, this may be considered as a representation in the in vitro open top case ‘s9. The fluid flow field in the IVD is presented in Figure 4 for various loading cases. Top and side views for each loading case are shown by Table Loading
2
Deflections
in the
Disc
Bulge
computing Vertebral vertical deflection
CXX
a b :
1.49 0.99 -0.99 0.72
(P) (P) (P) (L)
1.36 0.49 -0.72 -0.49
results
(A) (A) (A) (R)
2.63 1.66 -1.66 0.63
in units
of mm
End plate deflection
(A) (A) (A) (L)
2.30 1.41 -0.63 -1.41
(A) (A) (A) (R)
“+” in disc bulge specifies bulge out and “-” is the bulge vertical deflections specifies downward and *-” is upward. anterior, (P) is the posterior, (L) is the left, and (R) is the B, A, D and C in F@UP I, respectively)
220
vertical
in. “+” in (A) is the right (see
separate plots. Note that the convex of the top view specifies the posterior side. General comment from the study of the axi-symmetric models has shown that fluid will flow out through the EP under external compression; however, the fluid does not flow out of the IVD through the EP for the fully compressive case, as shown in Figure 4. From current demonstration and in view of structural mechanics, fluid may flow in or out through the IVD depending on the external loading cases. The mechanical behaviour of the SMS involves solid-fluid interaction; so, the fluid motion is deeply affected by the deformed shapes of the SMS. The fluid flow fields shown here may suggest reconsideration of many related works. Figure5 shows the von Mises stress contour at selected, significant sections A-B and C-D. Values labeled on each contour line have been normalized with respect to the surface pressure. The von Mises stress levels for the anterior bending case are the same as for the posterior bending case because it is just a magnitude of failure criterion. In the figure, it is obvious that the cortical shell at the anterior side experiences high stress levels. Especially, high level stress (above 7.0) is found at the mid-height, anterior of the VB. The stress level in the cancellous bone and the IVD is quite low (beween 0.0 and 1.0). Results for the lateral bending case also show that the stress level is high in the cortical shell and quite low within the cancellous bone and the intervertebral disc. However, the lateral bending case may produce smaller stress than the other cases on the cortical shell. Stress levels in the transverse direction of the SMS are chosen from 7 mm, 11 mm, and 26 mm sections. These selected sections are used to study the von Mises stress levels in the IVD, EP and VB, respectively. Figure 6 shows that the stress level in the nucleus is quite low and is higher in the annulus than in the nucleus for all cases. The fully compressive case gives higher stress in the annulus than the other three cases. However, a comparable level of stress is generated in the FJs for both the fully compressive and the lateral bending cases. No significant levels may be found for other two loading cases. In the FJs, stress is higher close to the anterior side than the posterior side. The stress level in. the EP is high, about 4-8 times to that in the annulus of the IVD, as shown in Figure 7. The stress level increases in the FJs as the height of the FJs increases. The fully compressive case gives higher stress than the other three cases in the EP. However, comparable level of stress is generated in the FJs for both the fully compressive and the lateral bending cases. No significant levels in the FJs can be found for the other two loading cases. Again, stress is higher closer to the anterior side than the posterior side in the FJs. F&UW 8 shows that the stress level in the cancellous bony part of the VB is quite low but about 8-10 times higher in both the cortical shell and the FJs. The anterior end of the VB exhibits high stress levels for all loading cases except lateral bend. Pore fluid pressures obtained in the disc
cr-rln I 1
L.C. a
L.C. b
-.
\ II’
j
~\--I
.-
L.C. c
-! 1, -i:i.l' , 1
8.
)
-’
L.C. a
L.Cs. b and c i
Di
L.C. d
221
Spinal
segment poroelastic
mixed mode1:JS.S.
0.5‘1,
\
Wu andJH.
/OS )
L.C. a Figure
6
Chen
Stress contours
on section
L.Cs. b and c 7 mm
(in N disc) for various
L.C. a Figure
7
Stress contours
are quite low shown here.
on section
(below
cases
L.Cs. b and c 11 mm
0.001)
(in endplate)
which
for various
are
not
CONCLUSION In the above, we have shown that the mechanical behaviour of the SMSs for various loading conditions can be clarified through current study by a threedimensional poroelastic FEM. Only the long time creep response is demonstrated here. From the mechanical point of view the facet joint plays an important role in the overall mechanical behaviour of the SMSs; therefore, the axi-symmetric model of the SMS is difficult to work cor-
222
loading
L.C. d
loading
L.C. d
cases
rectly because the FJs are- hard to include. general results may be concluded as:
Some
1. The disc bulge, the EP and the vertebral deflections must be read from the threedimensional model instead of axi-symmetric models. the IVD 2. Fluid may flow in or out through depending on the external loading cases. 3. The fluid motion is deeply affected by the deformed shapes of the SMS. 4. Stresses induced in the VB are always larger at the anterior side than the posterior side. 5. Stresses induced in the cortical shell are much higher than in the cancellous bone.
,11 1/-.. ‘\._A[ I .5
L.&T
6.
Stresses concentrate on the EP, cortical shell and the FJs. However, the highest stress occurs in the FJs. 7. The stress level in the MI is quite low but higher in the annulus than in the nucleus. a. Stress levels generated for the lateral bending case are lowest in all loading cases. 9. . The pore fluid pressure in the IVD is quite low. 10. FJs play an important role in the mechanical behaviour of the SMSs; so must be included in the analysis. Since the model set up here has been found valuable from the presentation, it can give correct information of mechanical behaviour in the SMS and provide medical areas valuable guidance in solving clinical problems occurring in the spine. Our goal is to extend current work using the substructure technique to simulate each of the SMSs; this will form a complete and widely usable model of the human spine to solve more sophisticated clinical problems.
and c
ACKNOWLEDGEMENTS The authors would like to express thanks for the project funding from National Science Council of the ROC under grant (NSC81-0401-EOOFi-07, 1992-1993). Valuable discussions related to pathological problems in lumbar segments with Dr Thomas MD Poon are highly appreciated.
REFERENCES I. Schultz AB, Belytschko T. Andriacchi TP. Analog studies of forces in the human spine: mechanical properties and motion segment behavior. J B&EC/L 1973; 6: 373-83. AI., Schultz AB, Berkson MH. Mechanical 2. Nachemson properties of human lumbar spine motion segments, influences of ages,’ sex, disc level, and degeneration. Spine (European edition) 1979; 4: l-5. 3. Koreska J, Roberson D, Mills RH, Gibson DA, Allisser AM. Biomechanics of the lumbar spine and its clinical significance. Orth Clinics ofN America: 1977; S(1): 121-33. 4. Belytschko T, Kulak RF, Schultz AB. Finite element analysis on an intervertebral disc. J Biomech 1974; 7: 277-85. 5. Berkson MH, Schultz AR, Warwick DN, Nachemson AL. Mechanical properties of human lumbar spine motion
223
Spinal
6. 7.
8.
9.
10.
11.
12.
13.
14. 15. 16. 17. 18.
224
segment poroel~tic
mixed model: J.S.S. Wu and J.H.
Chen
segment-Part III response in compression and shear; influence of gross morphology. J Biomech Engng, ASME Truns 1979; 101: 53-7. Hakim NS, Ring AI. A three dimensional finite element dynamic response analysis of a vertebra with experimental verification. J Biowch 1979; 12: 277-92. Laible JP, Pflaster D, Simon BR, Krag MH, Pope M, Haugh LD. A dynamic material parameter estimation procedure for soft tissue using a poroelastic finite element model, J Biomech Eng, 1994; 116: 19-29. Laible JP, Pflaster DS, Krag MH, Simon BR, Haugh LD. A poroelastic-swelling finite element model with application to the intervertebral disc, Spine, 1993; 18(5): 659-70. Simon BR, Wu JSS, Carlton MW, Kazarian LE, France EP, Evans JH, Zienkiewicz OC. Poroelastic dynamic structural models of rhesus spinal motion segments. Spine (European edition) 1985; lO(6): 494-507. Simon BR, Wu JSS, Carlton MW, Evans JH, Kazarian LE. Structural models for human spinal motion segments based on poroelastic view of the intervertebral disc. J Biomech lkgng 1985; 107: 327-35. Wu JSS, Huang JC, Lee CM, Simon BR: Derivation of the physical parameters of spinal motion segment based on porous medium theory. J Biompd Engng Sot ROC 1986; 6( 1): 55-66. Wu JSS, Huang JC, Lee CM, Simon BR. Dynamic analysis of axi-symmetric finite element models of spinal motion segments under anti-symmetric loads by a mixed procedure. J Biomed &gng Sot ROC 1986; 6(3): 157-81. Wu JSS. Dynamic analysis of poroelastic finite element models of rhesus spinal motion segments by the mixed procedure. Proc National Science Council ROC 1988; 12(6): 385-99. Wu JSS. Simulation of the structural model of spinal motion segments using fully mixed procedure. J Biomed Engng Sot ROC 1989; 8(2): 77-94. Biot MA. General theory of three-dimensional consolidation. JApp Physics 1971; 2: 155-64. Zienkiewicz OC, Shiomi T. Dynamic behavior of saturated porous media. Dept. Report CR/431/82, Dept. Chem Engng, Univ College, Swansea, October, 1982. Zienkiewicz OC, Tayler RL. The Finite EZement Method, 4th Ed., Volume 1: Basic formulations and linear problems, McGraw-Hill, 1989. Wu JSS. Generation of a three-dimensional poroelastic finite element model of the spinal motion segments by
energy theory. Project Report (NSC81-0401-E005-07, 19921993), Taipei, Taiwan, ROC. August 1993. 19. Wu JSS, Wang D. An improved three-dimensional poroelastic finite element model of the human spine. ChineseJ Med Biol Engng 1993; 13(2): 103-24.
APPENDIX:
NOMENCLATUREl
Strain-displacement matrix for u Strain-displacement matrix for w Solid strains : Strain energy 1, I(2, K3 Alternate stiffness matrices Bulk modulus of solid skeletion Bulk modulus of solid matrix sU,,,.,L Stiffness matrices Strain-displacement operator m Equivalent Kronecker delta Porosity of the medium Interpolation function for u iJ Interpolation function for w N:: Bulk modulus of fluid Q Solid displacements b Fluid displacements V Volume W Relative fluid displacements Applied load potential energy w, W Work by external load Compressibility of solid matrix Fluid strain Total potential energy 7T Pore fluid stress von Mises solid skeleton stress ck!VM u Solid stresses Solid skeleton stresses Differential operator Bu
BW
superscripts T
*
Transpose matrix Finite element nodal
values
Subscripts
1,2,3 f S
t
Parameters in x-, y, r-directions Fluid Solid skeleton Total medium
Clarification of the mechanical behaviour of spinal motion segments through a threedimensional poroelastic mixed finite element model J.S.S. Wu and J.H.
Chen*
Department of Mechanical Engineering, National Chung-Hsing University, Taichung, Taiwan 402, ROC; *Department of Applied Physics, Chung-Yuen Christian University, Chung-Li, Taiwan, ROC
ABSTRACT The purpose of this study is to c1nr-i~ the mechanical hehn-iliour (d ,spinul motion .tegments through a proper num&cal model. The model constructed ran give correct irtformntion and provide medical fields with valuble Lguidnncr irr .solving clinical problems orcurrkg in the sfiine. A three-dimensional pornelastic finite element model of sfiinnl motion segments is r.onstructed and a mixed formulation is introdured. The
Keywords:
Spinal
segment,
Med.
Phys.,
1996,
Eng.
poroelastic Vol.
18, 215-221,
models,
mixed
elements.
JD FEM
April
INTRODUCTION The spine is a biomechanical structure composed of 24 articulated columns of vertebrae; each of which is able to move through a joint of intervertebral disc (MI). Each IVD is a symphysis between two vertebral bodies (VBs) and formed with three parts; the nucleus pulposus, the annulus fibrosis and the hyaline cartilage end plates (EPs) on the opposing surfaces of each VB. The posterior facet joint (FJ) is of the diarthrodial type; it serves to guide and steady the movements of the VBs on one another. Each pair of FJs contacts with the upper and the lower FJs accordingly. Ligaments surround both the VB and the IVTI. The mechanism of a spinal motion segment (SMS) is defined as the composition of the vertebral column, the intervertebral joint, the FJs and ligaments. Changes from normal mechanical behaviour in the SMS may cause the pathogenesis Correspondence ical Engineering, wan 402, ROC.
tinite
to: Professor JSS Wu, PhD, Department of MechanNational Chung-Hsing University, Taichung, Tai-
of herniation and degeneration of the IVD, EP fracture, spinal deformity, low back pain syndromes, etc. Thus, understanding of mechanical behaviour in SMSs is quite important in relation to clinical problems occurring in the spine. Schultz et aZ.’ discussed mechanical properties of mathematical models of three-dimensional SMSs. A review of the mechanics of lumbar discs, towards a better understanding of low back pain, was made by Nachemson in 1975”. The fundamentals of the biomechanics of lumbar spine and its clinical significance were described by Koreska et ab” A numerical model of an intervertebral disc was developed by Belytschko et al. in 19744 for a three-dimensional stress analysis where the vertebrae were idealized as rigid bodies, while the discs, ligaments, and connective tissues were represented by deformable elements. Berkson et uZ.~ in 1979 studied the mechanical behaviour of 42 fresh human cadaver lumbar segments in compression, and in anterior, posterior, and lateral shear. Hakim and King’ used three-dimensional finite element models (FEMs) to analyse the
Spinal
segment
proelastic
mixed model: J.S.S. Wu and J H. Chen
dynamic response of vertebra, with experimental verification. Many investigations have been presented using poroelastic models to study the mechanical behaviour of tissues recently’*8. All components in the SMS contain different amounts of fluid in the pores of their solid tissues and bones. Solid models do not provide characteristics such as compressibility of the disc, long term creep response, etc. in the study of the mechanical behaviour of the porous structures of SMSs. Assuming each component of the SMS is porous, can provide a physical basis for determining the overall mechanical behaviour of the SMS under various loadings. Mathematical simulation of the SMSs using poroelastic axi-symmetric models can provide further insight ’ - l4 . However, since axi-symmetric models are difficult to apply to the complete geometry of SMSs, e.g., the non-circular transverse section of the SMS and the FJ part, correct mechanical behaviour of the objects is always misrepresented. Generally, the FJ exhibits strong support for the SMS that may affect the overall structural behaviour of the SMS significantly. For these reasons, a three-dimensional poroelastic FEM of a representative SMS, say, lumbar 5, is made to study the long term creep behaviour of the SMS here. The study can give correct information of mechanical behaviour in the SMS and provide medical areas valuable guidance in solving clinical problems occurring in the spine. The derivation of the constitutive relations based on energy laws is first briefly discussed and physical parameters in the equations of motion of porous media are interpreted. Then, the spatial discretization of the U-WV mixed finite element procedure is introduced. (Notations are defined in the Appendix.) Iteration in the computing process is performed in order to meet the specific boundary requirements on the contact surface of FJs and ligaments. In the final part, results such as the deformed shape, the fluid flow field and the failure stress distribution, are given and discussed. DEFINITION
OF DISPLACEMENTS
A general theory of threedimensional consolidation problems was discussed by Biot in 19’7’1’5 where the constitutive relations of porous media were discussed. An alternate form of the constitutive relations of porous media was also discussed by Zienkiewicz and Shiomi”. Here, the constitutive relations of porous media are established based on the energy laws. The poroelastic structures are composed of a ‘solid skeleton’ and an ‘interstitial fluid’ where the voids of the solid skeleton are saturated with fluid. The displacement field of the solid skeleton is defined by a column vector u, u = {ur u, ~3) where subscripts 1, 2 and 3 indicate parameters on the Cartesian coordinates in X-, p and zdirections, respectively. The displacement of fluid is defined as U and U = {U, U, U3}. Thus, the fluid volume displaced through pores perpendicular to the x,-axis is equal to nU, where n is called the porosity, defined as the ratio of the fluid volume
216
to the volume of total medium (n = V,/ V,) . A more significant vector, w and w = {w, w, w3} is defined as the displacement field of the pore fluid relative to the solid skeleton; wi( i = 1,2,3) is interpreted as the amount of fluid displaced through pores of the total area perpendicular to the i-th direction. So, the actual mean displacement of fluid through pores is equal to U-u and U-u = (l/n)w. STRESS-STRAIN
RJZATIONS
In the bi-phasic structure, strains are generated in both the solid skeleton and the fluid portions. The strain of the solid skeleton is expressed by e with a linear strain-displacement form of e = L u. Here, e is a column vector specified by {e,, q, es, %, e13 e,,}; L is the general strain-displacement operator in solid mechanics”. The fluid volumetric strain is given by 5, defined by the change of the fluid volume in the bulk volume of the medium, and 5 = VTw where V = LTm = @/ax, a/a% a/ax,1 and m = (1 1 10 0 O}. Stresses at a point in porous media are generated in both the solid and the fluid phases. Such stresses are represented by the total stress u and the hydrostatic pore fluid pressure 7~. The total stress u, is defined as forces acting on the bulk area of the medium, and the pore fluid pressure, n, is defined as forces acting on the fluid area. An extra parameter, cy, called compressibility of the ‘solid matrix’ (a solid including the effect from fluid), is equal to l-(&/K,), where & is the bulk modulus of the solid skeleton and K, is the bulk modulus of the solid matrix. The total stress u, accounting for total forces acting on the overall medium, is the sum of stresses in the solid skeleton (a’) and in the fluid (n). That is u=u’+amn,
(1) stress-strain
where a’ = De and D is the general matrix in solid mechanics”. The fluid strain is defined as n/Q where Q is the bulk modulus of fluid. For the hydrostatic case considered here, Q can be a scalar quantity onlvg*“. The fluid strain comes from two sources: (1)’ the fluid volumetric strain, 5, and (2) the excessive fluid strain, amTe, and is given by rr=a Substitution yields
Qm’e+
of equation
Qt (2) into
u=(D+a’@nmT)e+QJ.
(2)
equation
(1) (3)
Equations (2) and (3) are the general stress-strain relations of porous media which may be used to formulate the boundary value problem of such structures accordingly. ENERGY
THEORY
The total otential energy, II, of a structure is defined as R = E+ W, where E is the strain energy and W, is the potential energy of the applied loads. The external work, W, equals -W,. The principle of stationary potential energy gives
.Spnnl
6n=6E+6Wp=i3E-6W=0.
(4)
Generally, the study of mechanical behaviour of the SMSs involves consideration of the motion of stress waves that is time-dependent. Such timecreep and transient dependent cases, including cases depending on the applied loading rate, may raise the responding stress levels significantly. However, for the case of quite slow loading rate, the so-called long term creep case, or the static case for elastic structures over a long time interval, deformed rates may be neglected”. Here, the long term creep case only is introduced. Therefore, the derivation of the stiffness matrix [r;rl is provided. The variation of strain energy in the overall system, 6E:,, is derived as
vgmmt
proelnctic
tions (2), (3)) (5) and we obtain approach, manipulation:
Vn=
(l/k)*+
p&J;
/-fiindV=
jS(VTw)TmlI'.
(7) and the
using matrix
1%‘~ and ,I. H. Chev
the Galerkin form, after
F: 0 FZ
i I
(8)
where:
Y,,,=I B! DB,,dV &., = B,‘;cYN,dV I = N;B,,,dV
(61
(7)
Substituting equations (5) and (7) into equation (4)) two sets of equilibrium equations for the bi-phasic structures can be derived. Using the expression of stresses in equations (2) and (3) for u and 7~, and discretizing them into the proper form based on pertinent variables, the problem can then be solved. Details concerning the derivation of dynamic and kinetic relations of orous media can be found in the report by Wu’ R).
u-WV MIXED
[.,\.S.
con-
where pf is the fluid density and the dot specifies the time rate. The strain energy in fluid can be written in the form
SE,=
motirl:
.--
&, Darcy’s law when accounting for the general stant ‘permeability’, k, in fluid is”
mixrd
and B,, (= L) and B,v (= L m” = V) are the straindisplacement matrices of ZLand ZLI. Equation (8) can be solved and must be used when the solid skeleton and fluid are both incompressible, i.e., ~~41, Q+m and I<71*= 0. When Q is infinite, equation (8) cannot be solved by the Gaussian elimination procedure due to the singularity. In such case, reorder of parameters may be needed resulting in a non-symmetric form. Such non-symmetric fbrm is not encouraged in FEM analysis. In addition, when Qis finite for the compressible case, vector n* may be eliminated. Thus, we must have the equation. from the second row of equation (8), as m- = K$ (K&q* + K,,w* ) Substituting equations the or’ term, we have:
FORMULATION
Now, the Galerkin approachI is applied to derive proper formulas. Many alternative forms and solutions based on different parameters of porous media, have been studied and demonstrated”-” for axi-symmetric poroelastic problems. In a general finite element process, stresses are postinterpolated on the Gaussian points from the nodal displacements. Such post-interpolation procedures fit the resulting stresses poorly”. In addition, the surface condition of the pore fluid pressure in SMSs is quite important for these interactive, bi-phasic structures. Thus, a threedimensional mixed form, say u-WTT, is recommended here. This strategy may not only allow imposition and derivation of the pore fluid pressure on finite element nodes directly, but also include the effect of interaction between the pressure and displacements completely’H. In order to derive the finite element mixed form, we use the spatial semi-discretization with interpolation functions such as u = N,,u*, w = N,w*, rr = N,rr*, etc. where superscript * indicates finite element nodal values. Applying equa-
(9)
(9) and (8) and eliminating
(10) where:
Equation (10) is an alternative form of u-7~7~formulation where the benefit of the direct imposition and derivation of the pore fluid pressure can be obtained on finite element nodes. Equations (9) and (10) are both symmetric; thus, a partition frontal algorithm is implemented to carry out solutions for the large porous structure of SMSs. Once the displacements and fluid pressures are derived, total stresses may then be obtained by equation (3). (The total stress is interpolated on the Gaussian points first and then extrapolated to the elemental nodes by a specific designed algorithmlH.) The stress distributed in the solid skeleton is mostl) concerned in consideration of the failure of the poroelastic structure. The solid skeleton stress, u’, called the effective stress, is dcfincd as”’
217
Spinal
segment poroelastic
mixed model: J.S.S. Wu and J. H. Chen
u’=u-am7r
(13)
where
In this study, the von Mises criterion as a measure of failure of materials is used to quantify stresscomponents to arrive at one stress level on each node. The von Mises criterion applied to the solid skeleton may be expressed asI’: u;u
=
J
$(a:
- cq
+ $(a;, - 0:)’
+ t(a:
- u;,*
+ 375 + 3T; + 3TZI (14)
where a:,is called the effective von Mises stress, used to estimate the stresslevels generated in the SMSs in the following demonstrations. FINITE
ELEMENT
MODEL
Usually, each SMS contains three major parts: the VB, the IVD and the FJ. The VB is mainly constructed by an articulated column of the cancellous bone surrounded by a relatively stiff cortical shell. In the IVD, the pulposus, hydrated nucleus is centrally located and is relatively soft compared to the surrounding part of the fibrosus annulus. A porous, thin but rather stiff cartilage part, the EP, is located between the VB and the IVD. A pair of FJs is located at the posterior side of each VB. The basic function of the FJs is to guide and to steady the movements of each segment. In addition, there are ligaments surrounding the VB and the IVD. Here, the geometry of the SMS model is constructed from CT scanning images sectioned with the uniform space of 1 mm on a lumbar 5 SMS. This L5 SMS specimen was taken from a healthy man aged 46 with a height of 1.75 m and weight 160 pounds. The FEM of the L5 SMS is generated with brick elements and its degenerated elements automatically. Only one half of the height of IVD is included in the model for consideration of the opposing motion of the model at the mid-section of the IVD’,‘. Figure 1 shows the FEM and the reconstructed solid model. The definition of the height of each trans-
Figure
218
1
The finite
element
model
and the solid model
of the L5 spinal
verse section is shown in Fz@re 2. The above two figures will be used to calculate the deformed shapes and to show the stress distribution in the model. The L5 SMS model studied here is porous and elastic. Thus, material properties for each component include those in the solid skeleton, fluid and the solid matrix, see equation (8). The percentages of fluid contents (porosity n) in the cancellous bone, cortical shell, the EP and the FJs are relatively low, and may be ignored when compared to that in the MS~‘“,‘4. The porosity in the annulus and nucleus are both taken to be 0.70. All components including cancellous bone, cortical shell, EP, nucleus, annulus and FJs, have the same compressibility (Y taken as 0.9. The rest of the poroelastic material properties, such as the bulk modulus of the fluid Q are calculated from the fluid content and compressibility of the solid matrix”. Young’s modulus for each low fluid content is taken from the report by Belytschko in 19784. All other Young’s moduli are derived by a modified nonlinear optimization technique on experimental results”,18. Ligaments are simulated by an axial linear spring element with a constant of 15 N/mm correctly applied to nodes surrounding the VB and the IVD”~‘“. FJs take the same as for the cortical shell. material properties Material properties for each component are listed in Table 1. It is generally known that the human spine may transfer loads from’ the head to the trunk and to
Figure
motion
2
Definition
segment
of the section
height
-
30 26 22 18
z
;;t
50 46 42
38 34
Table 1 Material properties of the poroelastic of the 15 spinal motion segment Constituent
a
Q N/mm”
finite
I:’ N/mm”
element
n
model
v
flow and restricted to move horizontally. Each node on the contact surfaces of FJs is properly selected and its slope is estimated by interpolating the position of its neighbouring nodes in space. In the current FEM model, 1362 elements and 2144 nodes are built up and total degree of freedom is over 15,000. So, a specific frontal width optimization technique has been studied and implemented before solving the problem by the partition frontal solver here”. ITERATIONS
the pelvis which allow physiological motion among these three parts. So, the model set up here, does not relate to any postures of the human body but only considers forces transmitted fromthe earlier to the later segments. The upper surface of the VB is made to be the main region to receive loads. In this study, four loading cases are simulated: full compression, anterior bend, posterior bend and lateral bend. These four loading cases may be explained as the possible actions on the lumbar SMS caused by body motions. A uniform pressure of 0.4309 N/mm2* is applied to the upper surface of the VB’.“‘, see F@~re?. (This loading pressure has been used for the study of the mechanical behaviour of axi-symmetric models of SMS by the first author before. This report will compare some results with the axi-symmetric models.) Boundary conditions imposed on the structure are: (1) Nodes on the contact surface of facet joints are assumed to be inclined, (2) the bottom of the IVD in the FEM is impervious to the fluid
(a) Full
(c) Anterior
Compression
Bend
FOR
SPECIFIC
STRUCTURES
In computation, inclined boundary conditions of any compressive nodes on the contact surface of FJs are imposed and the element of ligaments are put on any tensile nodes located on the outer surface of the VB and the IVD. The strategy of accounting such restrictions in the computational procedure is explained as follows. In the initial step of computation, all nodes on the contact surface of I?Js are assumed to be compressive (inclined boundaries are imposed), and ligaments are assumed not appear. After the computation, if any nodal stresses on the FJ contact surface and the outer surface of the VB and the IVD are examined with different phases (compression or tension) from the latest assumption, a seconda? computation is thus necessary. In the second computation, adjustment of the boundary conditions is first considered. That is, (1) inclined boundary conditions of any compressive nodes on the contact surface of the FJs are imposed or otherwise released, and (2) the ligament elements are put on anv tensile nodes on the outer surface of the
(b) Posterior
Bend
(d) Lateral Bend
219
Spinal
segment poroelastic
mixed model: J.S.S. Wu andJH.
Chen
VB and the IVD or otherwise are taken away. Nodal stresses on the contact surface and the outer surface of the VB and the IVD are checked after the second computation. Such iterative computation continues until the nodal stress phase (compression or tension) on the contact surface of FJ and the outer surface of the VB and the IVD are maintained. RESULTS
AND
DISCUSSIONS
The mechanical behaviour of an L5 SMS is studied through a three-dimensional poroelastic mixed FEM. Structural analysis of the SMS considers certain displacements and stress levels as a measure of the failure of materials under specific loading and boundary conditions. Here, only the long time creep response is demonstrated. Figure 3 shows the deformed solid shapes of the model for cases of (a) full compression, (b) posterior bend, (c) anterior bend and (d) lateral bend. It clearly shows that the deformed shape of each component in the SMS is never symmetric with respect to the central line of either the VB or the IVD. Thus, the mechanical behaviour of the SMS is difficult to interpret by an axi-symmetric structural model. There is no doubt that a three-dimensional model of the SMSs is valuable in providing clinical information on the spine. Table 2 lists the deformed sizes of the model. The disc bulge and vertebral vertical deflection calculated by the axi-symmetric model under full compressive load are 0.254 mm. and 0.385 mm., respectively’,‘; so, the results shown here are physically significant. The largest disc bulge and vertebral vertical deflection appear in the fully compressive case. The smallest posterior and anterior bulges and vertebral vertical deflection appear in the lateral bending case, no obvious bulging is found for this case. The above phenomena are caused by the strong support from the FJs; so, the disc bulge, the EP and the vertebral deflections are obviously not able to be read from axi-symmetric models. Since the fluid content in the IVD is relatively high compared to others, the motion of fluid in the IVD is discussed here. Also, this may be considered as a representation in the in vitro open top case ‘s9. The fluid flow field in the IVD is presented in Figure 4 for various loading cases. Top and side views for each loading case are shown by Table Loading
2
Deflections
in the
Disc
Bulge
computing Vertebral vertical deflection
CXX
a b :
1.49 0.99 -0.99 0.72
(P) (P) (P) (L)
1.36 0.49 -0.72 -0.49
results
(A) (A) (A) (R)
2.63 1.66 -1.66 0.63
in units
of mm
End plate deflection
(A) (A) (A) (L)
2.30 1.41 -0.63 -1.41
(A) (A) (A) (R)
“+” in disc bulge specifies bulge out and “-” is the bulge vertical deflections specifies downward and *-” is upward. anterior, (P) is the posterior, (L) is the left, and (R) is the B, A, D and C in F@UP I, respectively)
220
vertical
in. “+” in (A) is the right (see
separate plots. Note that the convex of the top view specifies the posterior side. General comment from the study of the axi-symmetric models has shown that fluid will flow out through the EP under external compression; however, the fluid does not flow out of the IVD through the EP for the fully compressive case, as shown in Figure 4. From current demonstration and in view of structural mechanics, fluid may flow in or out through the IVD depending on the external loading cases. The mechanical behaviour of the SMS involves solid-fluid interaction; so, the fluid motion is deeply affected by the deformed shapes of the SMS. The fluid flow fields shown here may suggest reconsideration of many related works. Figure5 shows the von Mises stress contour at selected, significant sections A-B and C-D. Values labeled on each contour line have been normalized with respect to the surface pressure. The von Mises stress levels for the anterior bending case are the same as for the posterior bending case because it is just a magnitude of failure criterion. In the figure, it is obvious that the cortical shell at the anterior side experiences high stress levels. Especially, high level stress (above 7.0) is found at the mid-height, anterior of the VB. The stress level in the cancellous bone and the IVD is quite low (beween 0.0 and 1.0). Results for the lateral bending case also show that the stress level is high in the cortical shell and quite low within the cancellous bone and the intervertebral disc. However, the lateral bending case may produce smaller stress than the other cases on the cortical shell. Stress levels in the transverse direction of the SMS are chosen from 7 mm, 11 mm, and 26 mm sections. These selected sections are used to study the von Mises stress levels in the IVD, EP and VB, respectively. Figure 6 shows that the stress level in the nucleus is quite low and is higher in the annulus than in the nucleus for all cases. The fully compressive case gives higher stress in the annulus than the other three cases. However, a comparable level of stress is generated in the FJs for both the fully compressive and the lateral bending cases. No significant levels may be found for other two loading cases. In the FJs, stress is higher close to the anterior side than the posterior side. The stress level in. the EP is high, about 4-8 times to that in the annulus of the IVD, as shown in Figure 7. The stress level increases in the FJs as the height of the FJs increases. The fully compressive case gives higher stress than the other three cases in the EP. However, comparable level of stress is generated in the FJs for both the fully compressive and the lateral bending cases. No significant levels in the FJs can be found for the other two loading cases. Again, stress is higher closer to the anterior side than the posterior side in the FJs. F&UW 8 shows that the stress level in the cancellous bony part of the VB is quite low but about 8-10 times higher in both the cortical shell and the FJs. The anterior end of the VB exhibits high stress levels for all loading cases except lateral bend. Pore fluid pressures obtained in the disc
cr-rln I 1
L.C. a
L.C. b
-.
\ II’
j
~\--I
.-
L.C. c
-! 1, -i:i.l' , 1
8.
)
-’
L.C. a
L.Cs. b and c i
Di
L.C. d
221
Spinal
segment poroelastic
mixed mode1:JS.S.
0.5‘1,
\
Wu andJH.
/OS )
L.C. a Figure
6
Chen
Stress contours
on section
L.Cs. b and c 7 mm
(in N disc) for various
L.C. a Figure
7
Stress contours
are quite low shown here.
on section
(below
cases
L.Cs. b and c 11 mm
0.001)
(in endplate)
which
for various
are
not
CONCLUSION In the above, we have shown that the mechanical behaviour of the SMSs for various loading conditions can be clarified through current study by a threedimensional poroelastic FEM. Only the long time creep response is demonstrated here. From the mechanical point of view the facet joint plays an important role in the overall mechanical behaviour of the SMSs; therefore, the axi-symmetric model of the SMS is difficult to work cor-
222
loading
L.C. d
loading
L.C. d
cases
rectly because the FJs are- hard to include. general results may be concluded as:
Some
1. The disc bulge, the EP and the vertebral deflections must be read from the threedimensional model instead of axi-symmetric models. the IVD 2. Fluid may flow in or out through depending on the external loading cases. 3. The fluid motion is deeply affected by the deformed shapes of the SMS. 4. Stresses induced in the VB are always larger at the anterior side than the posterior side. 5. Stresses induced in the cortical shell are much higher than in the cancellous bone.
,11 1/-.. ‘\._A[ I .5
L.&T
6.
Stresses concentrate on the EP, cortical shell and the FJs. However, the highest stress occurs in the FJs. 7. The stress level in the MI is quite low but higher in the annulus than in the nucleus. a. Stress levels generated for the lateral bending case are lowest in all loading cases. 9. . The pore fluid pressure in the IVD is quite low. 10. FJs play an important role in the mechanical behaviour of the SMSs; so must be included in the analysis. Since the model set up here has been found valuable from the presentation, it can give correct information of mechanical behaviour in the SMS and provide medical areas valuable guidance in solving clinical problems occurring in the spine. Our goal is to extend current work using the substructure technique to simulate each of the SMSs; this will form a complete and widely usable model of the human spine to solve more sophisticated clinical problems.
and c
ACKNOWLEDGEMENTS The authors would like to express thanks for the project funding from National Science Council of the ROC under grant (NSC81-0401-EOOFi-07, 1992-1993). Valuable discussions related to pathological problems in lumbar segments with Dr Thomas MD Poon are highly appreciated.
REFERENCES I. Schultz AB, Belytschko T. Andriacchi TP. Analog studies of forces in the human spine: mechanical properties and motion segment behavior. J B&EC/L 1973; 6: 373-83. AI., Schultz AB, Berkson MH. Mechanical 2. Nachemson properties of human lumbar spine motion segments, influences of ages,’ sex, disc level, and degeneration. Spine (European edition) 1979; 4: l-5. 3. Koreska J, Roberson D, Mills RH, Gibson DA, Allisser AM. Biomechanics of the lumbar spine and its clinical significance. Orth Clinics ofN America: 1977; S(1): 121-33. 4. Belytschko T, Kulak RF, Schultz AB. Finite element analysis on an intervertebral disc. J Biomech 1974; 7: 277-85. 5. Berkson MH, Schultz AR, Warwick DN, Nachemson AL. Mechanical properties of human lumbar spine motion
223
Spinal
6. 7.
8.
9.
10.
11.
12.
13.
14. 15. 16. 17. 18.
224
segment poroel~tic
mixed model: J.S.S. Wu and J.H.
Chen
segment-Part III response in compression and shear; influence of gross morphology. J Biomech Engng, ASME Truns 1979; 101: 53-7. Hakim NS, Ring AI. A three dimensional finite element dynamic response analysis of a vertebra with experimental verification. J Biowch 1979; 12: 277-92. Laible JP, Pflaster D, Simon BR, Krag MH, Pope M, Haugh LD. A dynamic material parameter estimation procedure for soft tissue using a poroelastic finite element model, J Biomech Eng, 1994; 116: 19-29. Laible JP, Pflaster DS, Krag MH, Simon BR, Haugh LD. A poroelastic-swelling finite element model with application to the intervertebral disc, Spine, 1993; 18(5): 659-70. Simon BR, Wu JSS, Carlton MW, Kazarian LE, France EP, Evans JH, Zienkiewicz OC. Poroelastic dynamic structural models of rhesus spinal motion segments. Spine (European edition) 1985; lO(6): 494-507. Simon BR, Wu JSS, Carlton MW, Evans JH, Kazarian LE. Structural models for human spinal motion segments based on poroelastic view of the intervertebral disc. J Biomech lkgng 1985; 107: 327-35. Wu JSS, Huang JC, Lee CM, Simon BR: Derivation of the physical parameters of spinal motion segment based on porous medium theory. J Biompd Engng Sot ROC 1986; 6( 1): 55-66. Wu JSS, Huang JC, Lee CM, Simon BR. Dynamic analysis of axi-symmetric finite element models of spinal motion segments under anti-symmetric loads by a mixed procedure. J Biomed &gng Sot ROC 1986; 6(3): 157-81. Wu JSS. Dynamic analysis of poroelastic finite element models of rhesus spinal motion segments by the mixed procedure. Proc National Science Council ROC 1988; 12(6): 385-99. Wu JSS. Simulation of the structural model of spinal motion segments using fully mixed procedure. J Biomed Engng Sot ROC 1989; 8(2): 77-94. Biot MA. General theory of three-dimensional consolidation. JApp Physics 1971; 2: 155-64. Zienkiewicz OC, Shiomi T. Dynamic behavior of saturated porous media. Dept. Report CR/431/82, Dept. Chem Engng, Univ College, Swansea, October, 1982. Zienkiewicz OC, Tayler RL. The Finite EZement Method, 4th Ed., Volume 1: Basic formulations and linear problems, McGraw-Hill, 1989. Wu JSS. Generation of a three-dimensional poroelastic finite element model of the spinal motion segments by
energy theory. Project Report (NSC81-0401-E005-07, 19921993), Taipei, Taiwan, ROC. August 1993. 19. Wu JSS, Wang D. An improved three-dimensional poroelastic finite element model of the human spine. ChineseJ Med Biol Engng 1993; 13(2): 103-24.
APPENDIX:
NOMENCLATUREl
Strain-displacement matrix for u Strain-displacement matrix for w Solid strains : Strain energy 1, I(2, K3 Alternate stiffness matrices Bulk modulus of solid skeletion Bulk modulus of solid matrix sU,,,.,L Stiffness matrices Strain-displacement operator m Equivalent Kronecker delta Porosity of the medium Interpolation function for u iJ Interpolation function for w N:: Bulk modulus of fluid Q Solid displacements b Fluid displacements V Volume W Relative fluid displacements Applied load potential energy w, W Work by external load Compressibility of solid matrix Fluid strain Total potential energy 7T Pore fluid stress von Mises solid skeleton stress ck!VM u Solid stresses Solid skeleton stresses Differential operator Bu
BW
superscripts T
*
Transpose matrix Finite element nodal
values
Subscripts
1,2,3 f S
t
Parameters in x-, y, r-directions Fluid Solid skeleton Total medium