﻿
 UFDC Home myUFDC Home  |   Help
<%BANNER%>

A Two-dimensional human spine simulation and three-dimensional spine model construction

PAGE 1

A TWO-DIMENSIONAL HUMAN SPINE SIMULATION AND THREEDIMENSIONAL SPINE MODEL CONSTRUCTION By JIANZHI LIU A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2001

PAGE 2

ii ACKNOWLEDGMENTS Many thanks are in order for Dr. Andrew J. Kurdila for his continual assistance, expertise, and leadership. I would like to thank Dr. Andrew J. Rapoff and Dr. B. J. Fregly for their assistance and expertise in helping me on this project. I would also like to thank all the faculty and students who helped at different times during my research.

PAGE 3

iii TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ ... ii LIST OF TABLES ................................ ................................ ................................ ............... v LIST OF FIGURES ................................ ................................ ................................ ........... vii ABSTRACT ................................ ................................ ................................ ......................... x CHAPTERS 1 INTRODUCTION AND LITERATURE REVIEW ................................ ......................... 1 1.1 Introduction ................................ ................................ ................................ ............... 1 1.2 Literature Review ................................ ................................ ................................ ...... 3 1.3 Goal and Outline ................................ ................................ ................................ ....... 6 2 MULTIBODY DYNAMICS THEORY ................................ ................................ ........... 8 2.1 Two-Dimensional Formulation Example ................................ ................................ .. 8 2.2 Identification Problem in the Two-Dimensional Example ................................ .... 15 3 NUMERICAL EXPERIMENT ................................ ................................ ....................... 22 3.1 Multibody Dynamics Simulation ToolsPromechanica ................................ ....... 22 3.2 Forward Simulation with Custom Load ................................ ................................ .. 24 3.3 Identification of Pin Joint Position Using Inverse Simulation ................................ 36 4 SCOLIOSIS MODEL DATABASE CONSTRUCTION ................................ ............... 48 4.1 Geometry ................................ ................................ ................................ ................. 48 4.2 Stiffness ................................ ................................ ................................ ................... 58 4.3 A Simple Demonstration of the Stiffness Database ................................ ................ 75 5 CONCLUSIONS AND RECOMMENDATIONS ................................ ......................... 77

PAGE 4

iv APPENDICES A GEOMETRY OF MODELING NODES ................................ ................................ ....... 78 B ATTACHING POSITIONS OF SCREW & HOOKS ................................ .................... 95 LIST OF REFERENCES ................................ ................................ ................................ 112 BIOGRAPHICAL SKETCH ................................ ................................ ........................... 114

PAGE 5

v LIST OF TABLES Table Page 3-1 8 Springs identification ................................ ................................ ................................ .... 34 3-2 16 Springs identification ................................ ................................ ................................ .. 34 4-1 Scoliosis curve data ................................ ................................ ................................ ......... 50 4-2 Coordinates of origin nodes ................................ ................................ ............................. 51 A1 Lumbar 5 ................................ ................................ ................................ .......................... 78 A2 Lumbar 4 ................................ ................................ ................................ .......................... 79 A3 Lumbar 3 ................................ ................................ ................................ .......................... 80 A4 Lumbar 2 ................................ ................................ ................................ .......................... 81 A5 Lumbar 1 ................................ ................................ ................................ .......................... 82 A6 Thoracic 12 ................................ ................................ ................................ ...................... 83 A7 Thoracic 11 ................................ ................................ ................................ ...................... 84 A8 Thoracic 10 ................................ ................................ ................................ ...................... 85 A9 Thoracic 9 ................................ ................................ ................................ ........................ 86 A10 Thoracic 8 ................................ ................................ ................................ ...................... 87 A11 Thoracic 7 ................................ ................................ ................................ ...................... 88 A12 Thoracic 6 ................................ ................................ ................................ ...................... 89 A13 Thoracic 5 ................................ ................................ ................................ ...................... 90 A14 Thoracic 4 ................................ ................................ ................................ ...................... 91

PAGE 6

vi A15 Thoracic 3 ................................ ................................ ................................ ...................... 92 A16 Thoracic 2 ................................ ................................ ................................ ...................... 93 A17 Thoracic 1 ................................ ................................ ................................ ...................... 94 B1 Screw & Hooks on L5 (mm) ................................ ................................ ............................ 95 B2 Screw & Hooks on L4 (mm) ................................ ................................ ............................ 96 B3 Screw & Hooks on L3 (mm) ................................ ................................ ............................ 97 B4 Screw & Hooks on L2 (mm) ................................ ................................ ............................ 98 B5 Screw & Hooks on L1 (mm) ................................ ................................ ............................ 99 B6 Screw & Hooks on T12 (mm) ................................ ................................ .......................... 100 B7 Screw & Hooks on T11 (mm) ................................ ................................ .......................... 101 B8 Screw & Hooks on T10 (mm) ................................ ................................ .......................... 102 B9 Screw & Hooks on T9 (mm) ................................ ................................ ............................ 103 B10 Screw & Hooks on T8 (mm) ................................ ................................ .......................... 104 B11 Screw & Hooks on T7 (mm) ................................ ................................ .......................... 105 B12 Screw & Hooks on T6 (mm) ................................ ................................ .......................... 106 B13 Screw & Hooks on T5 (mm) ................................ ................................ .......................... 107 B14 Screw & Hooks on T4 (mm) ................................ ................................ .......................... 108 B15 Screw & Hooks on T3 (mm) ................................ ................................ .......................... 109 B16 Screw & Hooks on T2 (mm) ................................ ................................ .......................... 110 B17 Screw & Hooks on T1 (mm) ................................ ................................ .......................... 111

PAGE 7

vii LIST OF FIGURES Figure Page 1-1 Human Spine ................................ ................................ ................................ .................... 1 2-1 Two vertebral body system ................................ ................................ .............................. 9 2-2 Two-body example illustration ................................ ................................ ........................ 11 2-3 Identification procedure illustration ................................ ................................ ................. 17 2-4 Numerical simulation procedure illustration ................................ ................................ ... 21 3-1 Human spine modeling procedure ................................ ................................ ................... 24 3-2 A two vertebral body model built in Promechanica ................................ ........................ 25 3-3 Generalized forces in two-dimensional plane ................................ ................................ .. 26 3-4 Y-translational trajectory numerical results from Promechanica ................................ .... 28 3-5 Y-translational trajectory numerical results from Matlab model ................................ ..... 28 3-6 Z-translational trajectory numerical results from Promechanica ................................ ..... 29 3-7 Z-translational trajectory numerical results from Matlab model ................................ ..... 29 3-8 X-rotational trajectory numerical results from Promechanica ................................ ......... 30 3-9 X-rotational trajectory numerical results from Matlab model ................................ ......... 30 3-10 Y-translational trajectory from Promechanica model ................................ .................... 31 3-11Y-translational trajectory from Matlab model ................................ ................................ 31 3-12 Z-translational trajectory from Promechanica model ................................ .................... 32 3-13 Z-translational trajectory from Matlab model ................................ ............................... 32 3-14 X-rotational trajectory from Promechanica model ................................ ........................ 33

PAGE 8

viii 3-15 X-rotational trajectory from Matlab model ................................ ................................ ... 33 3-16 Translation movement simulation ................................ ................................ .................. 35 3-17 Rotation movement simulation ................................ ................................ ...................... 35 3-18 vertebral body connected by pin ................................ ................................ .................... 36 3-19 Probability density distribution of spring stiffness ................................ ........................ 38 3-20 Cumulative distribution of spring stiffness ................................ ................................ .... 38 3-21 Probability density distribution of spring stiffness ................................ ........................ 39 3-22 Cumulative distribution of spring stiffness ................................ ................................ .... 39 3-23 Probability density distribution of spring stiffness ................................ ........................ 40 3-24 Cumulative distribution of spring stiffness ................................ ................................ .... 40 3-25 Probability density distribution of spring stiffness ................................ ........................ 41 3-26 Cumulative distribution of spring stiffness ................................ ................................ .... 41 3-27 Probability density distribution of spring stiffness ................................ ........................ 42 3-28 Cumulative distribution of spring stiffness ................................ ................................ .... 42 3-29 Probability density distribution of spring stiffness ................................ ........................ 43 3-30 Cumulative distribution of spring stiffness ................................ ................................ .... 43 3-31 Probability density distribution of spring stiffness ................................ ........................ 44 3-32 Cumulative distribution of spring stiffness ................................ ................................ .... 44 3-33 Probability density distribution of spring stiffness ................................ ........................ 45 3-34 Cumulative distribution of spring stiffness ................................ ................................ .... 45 3-35 Probabilty density distribution of different spring numbers ................................ .......... 46 3-36 The cumulative distribution of different spring numbers ................................ .............. 46 3-37 Error vs. different number of springs ................................ ................................ ............ 47 4-1 Spine Model in Promechanica Motion ................................ ................................ ............ 52 4-2 A Vertebral Model with Complex Surface ................................ ................................ ...... 53

PAGE 9

ix 4-3 A Simple Geometry Vertebral Model ................................ ................................ .............. 53 4-4 Right Lateral View of the Primary Nodes ................................ ................................ ...... 54 4-5 Top View of the Primary Nodes ................................ ................................ ...................... 55 4-6 Bottom View of the Primary Nodes ................................ ................................ ................ 56 4-7 Simplified Geometry Spine Model ................................ ................................ .................. 57 4-8 Illustration Stiffness Property in Three-dimensional Coordinates System ...................... 58 4-9 Flextion-extension between L1 and L2 ................................ ................................ ........... 60 4-10 Flextion-extension between L2 and L3 ................................ ................................ ......... 61 4-11 Flextion-extension between L3 and L4 ................................ ................................ ......... 62 4-12 Flextion-extension between L4 and L5 ................................ ................................ ......... 63 4-13 Flextion-extension between L5 and S1 ................................ ................................ ......... 64 4-14 Right and Left Axial Torque between L1 and L2 ................................ ......................... 65 4-15 Right and Left Axial Torque between L2 and L3 ................................ ......................... 66 4-16 Right and Left Axial Torque between L3 and L4 ................................ ......................... 67 4-17 Right and Left Axial Torque between L4 and L5 ................................ ......................... 68 4-18 Right and Left Axial Torque between L5 and S1 ................................ ......................... 69 4-19 Left and Right Bending Moment between L1 and L2 ................................ ................. 70 4-20 Left and Right Bending Moment between L2 and L3 ................................ ................. 71 4-21 Left and Right Bending Moment between L3 and L4 ................................ ................. 72 4-22 Left and Right Bending Moment between L4 and L5 ................................ ................. 73 4-23 Left and Right Bending Moment between L5 and S1 ................................ ................. 74 4-24 Load-Displacement Compare between Database and Reconstruction ......................... 76

PAGE 10

x Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science A TWO-DIMENSIONAL HUMAN SPINE SIMULATION AND THREEDIMENSIONAL SPINE MODEL CONSTRUCTION By Jianzhi Liu December 2001 Chairman: Dr. Andrew Kurdila Major Department: Aerospace Engineering, Mechanics, and Engineering Science This thesis describes two different, but complementary modeling tasks for the human spine. Methodologies for modeling and simulation of the motion of the human spine vary dramatically in complexity. While providing detailed stress and strain field representations, full finite element modeling of the human spine can be computationally expensive. Alternatively, nonlinear multibody dynamics representations are often used because of their simplicity. These formulations employ rigid models of vertebrae interconnected via conventional mechanical joints. However, it is well documented that inter-vertebral motion can depart significantly from conventional mechanical joint constraints. We present an identification methodology that employs a relaxation technique in which joint mechanical properties are represented via a probability measure. In addition we describe an effort to obtain accurate three-dimensional models of the human spine. These efforts are motivated by considering scoliosis. Scoliosis is

PAGE 11

xi defined as abnormal lateral curvature of the spine. It is usually considered as a threedimensional deformity, because axial rotation will always accompany the lateral curvature. The correction of the deformity is required when the patient risks severe deformity. A three-dimensional spine model is constructed and a computer simulation tool is provided.

PAGE 12

1 CHAPTER1 INTRODUCTION AND LITERATURE REVIEW Chapter 1 first talks about some fundamental concepts and issues related to the human spine modeling, then discusses some basic information about scoliosis. An outline of the following chapter is given last. 1.1 Introduction Because each spine is unique and reacts in a slightly different way, the dynamics of the human spine are difficult to quantify. Therefore, in order to study the dynamics of the spine, it is necessary to develop a mathematical model of the human spine that is general enough to encompass the possible variations in form and function, but specific enough to react realistically under loading and stress. Human vertebral column (Figure 1.1) is a complex structure and mechanism. It is made up of bony vertebrae interconnected by discs and ligaments. The normal spine consists of seven cervical, twelve thoracic, five lumbar, five sacral and five coccygeal vertebrae. The sacral and coccygeal vertebrae are usually fused. In the lateral plane, the normal spine exhibits physiologic cervical lordosis, Figure 1-1 Human Spine cervical vertebrae thoracic vertebrae sacrum coccyx cervical lordosis lumbar lordosis thoracic kyphosis

PAGE 13

2 thoracic kyphosis, and lumbar lordosis. They increase the flexibility and shockabsorbing capacity of the vertebral column, while maintaining the adequate stiffness and stability. Although a vertebra is not actually a rigid body, usually it can be considered to behave as such, because deformations within the vertebra are insignificant compared to the deformations of the connective tissue. Two adjacent vertebrae, the ligaments, intervertebral disc and the facets form a functional spinal unit, or a motion segment. It exhibits biomechanical characteristics similar to those of the entire spine. Each of the motion segments has six degree-of-freedom (DOF), and the whole spine can produce three movements: flexion extension, lateral bending, and axial rotation. We simply summarize them as the combinations of rotational and translational movements. The motions in human spine are coupled, where translation or rotation along or about one single axis is consistently associated with translation or rotation along or about a second axis. Some literature has shown these coupled motions clearly [Oxa92], [Pan94]. The coupled motion can lead to computation difficulties while modeling a complex structure like the human spine. Another fundamental issue related to our research is scoliosis. Scoliosis is defined as abnormal lateral curvature of the spine. It is usually considered as a threedimensional deformity, because axial rotation will always accompany the lateral curvature. Most scoliotic curves will stabilize if not treated, but a few will progress. Mild scoliotic curves result in cosmetic deformities of the trunk. These can create psychological problems for the patients who have them. Severe scoliotic curves may alter balance and coordination, interfere with the function of internal organs, allow

PAGE 14

3 premature degeneration of the spinal column, and cause deterioration of neurological function. Thus, the correction of the deformity is required when the patient risks severe deformity. The purpose of this correction is to stabilize the spine until a solid fusion occurs in the treated area. Once the spine has been fused, it can once again perform its function without the aid of the instruments. Hence, it is necessary to generate a general model that can represent these problems related to the human spine. Many efforts haven been made in the pass. They can be categorized into different sub areas according to the researchers' interests. 1.2 Literature Review Multiple biomechanical models have been employed for the spinal modeling. Most can be grouped into two classes. These classes are computational models and experimental models. The computational models can be used to provide insight into the mechanisms by which a particular motion occurs or can be used to investigate the complex dynamics of multi-linked segments during motion. Also, these models are used to undertake analyses of activities, which would otherwise be impossible to carry out experimentally. Computational modeling can be classified into two categories: one includes the numerical models, which usually use finite element methods. The other includes the analytical models, where the geometry is simplified and requires few boundary conditions. Both categories of models have their strengths and weaknesses, but the inherent simplicity of analytical models often leads to an intuitive understanding of human spine mechanisms.

PAGE 15

4 Not surprisingly, the computational cost and accuracy of these methods vary greatly. Detailed finite element discretisations of intervertebral discs can provide finescale stress and strain field approximations. Such highly resolved fields are obtained, of course, at significant computational cost [Bel73] [Bel74] [Kul75] [Kul76] [Sch73]. A recent study reports on the development and validation of a nonlinear viscoelastic FE model that qualifies the mechanical responses of the L2/L3 motion segment subset to time varying external loads. The vertebral bodies are represented as rigid bodies, and intervertebral discs are modeled via nonlinear finite elements. An alternative choice makes use of nonlinear multibody dynamics formulations [Buc96] [Hur94]. In these formulations, the vertebral bodies are assumed to be rigid and interconnected by a collection of conventional mechanical joints. The advantages of this formulation include the fact that these models are relatively low dimensional. It can often provide accurate predictions of large amplitude articulated motion. Optimization techniques are usually used in the nonlinear dynamics formulation, where objective functions are chosen, and a set of parameters is found that minimizes that function. One earlier example is the model developed by Schultz [Sch73]. In his study, Schultz described his two objective functions, which minimized the compression on the motion segment. In addition to computational modeling that determines the physics of spinal motion, a significant number of clinical experimental studies of spinal motion have appeared. Some of them are relatively old and represent early efforts to validate specific spine models [Kul75]. Other earlier experimental studies have estimated effective

PAGE 16

5 compliance parameters for motion segments. These are subsequently employed to provide more accurate lumped simulation models [Wil83]. There are some other experimental studies of spinal motion that have been carried out independent of numerical simulation, solely to examine and characterize the physics of spinal motion. Panjabi et al [Oxa91] [Pan94] investigated a three-dimensional human spine structure that exhibits non-linear viscoelastic mechanical behavior. As advanced measurement technologies have become more and more common, the experimental techniques and protocols have likewise matured. A hybrid approach called EMG assisted optimization (EMGAO) was developed by using both optimization and an EMG assisted approach. In other studies [Har98] [Ros98], lateral radiographs are employed to characterize spinal kinematics and frequency domain characterization of the spinal motion is studied. Despite the progress made in computational modeling and experimental methodologies to refine models, significant technical issue remains unresolved. As noted by M. Klenberger in [Kle98] Biological material data remains one of the limiting factors for biomechanical modeling. Existing data are comprised mostly of quasi-static force-displacement measurements with calculations of linear stiffness. Normal variance of biological tissues typically yields a standard deviation for stiffness values of at least 20 percent. No matter how refined an element mesh might be, the accuracy of the model can never be better than 20 percent. Thus, the uncertainty in data derived from biological specimens has lead to increased interest in identification and inverse methods, because identification and inverse methods by their nature are developed to accommodate some level of noise

PAGE 17

6 uncertainty and poorly conditioned data. Research [Wol94] [Ris97] [Kle98] [Xu99] provides recent examples in which identification and inverse methods have been developed in biomechanics. 1.3 Goal and Outline The objective of our research is to derive an identification method for human spinal motion that is based on a relaxation of the original governing equations. This approach is based on a simple observation. While some human joints are well represented by simple conventional mechanical joints, the relative motion of human vertebrae is not amenable to such a straightforward representation. In the following chapters, we present an identification methodology that employs a relaxation technique in which joint mechanical properties are represented via the probability measure. Numerical studies are used to study the performance of this approximation framework. Chapter 2 discusses the fundamentals of the multibody dynamics formulation. The study of the dynamics of the multi-body system includes a development of the governing equations of the motion. In the model, each spine motion segment has six degree of freedom (DOF), and the whole human spine can undergo the translational and rotational movements. Based on this observation, the governing equations of motion (EOM) of the human spine are developed by using Lagranges Equation. Based on the identification theory and algorithm introduced in chapter 2, the numerical experiments are implemented in chapter 3. In this chapter, several cases of identification are studied. A forward numerical simulation and an inverse identification numerical experiment are discussed. Most of the simulation work has been achieved using Promechanica, which is a multibody dynamics simulation software.

PAGE 18

7 In chapter 4, some basic tools for scoliosis research issues have been developed. These tools include a spine (L5-T1) model constructed using the simulation software Promechanica in chapter 4. To facilitate the future development of spinal modes incorporating rod, hook and screw assembles, a spinal simulation database was developed that includes both geometry and attachment parts. The attachment points of pedicle screw and hooks are also included in this geometry database. In addition, a stiffness database has been constructed from the published literature. At the end of this chapter, a numerical simulation experiment is studied to verify the stiffness properties. The detailed geometry database and stiffness database can be found in the Appendix. Chapter 5 provides concluding remarks and suggestions for future work.

PAGE 19

8 CHAPTER 2 MULTIBODY DYNAMICS THEORY In the first part of chapter 2, a quick review of fundamental multibody dynamics theory is given. The governing equations for the model are derived based on Lagranges Equations. In the second part, we discuss the basic procedures of numerical simulation. 2.1 Two-Dimensional Formulation Example Before talking about several aspects of the theoretical background for the relaxed formulation used in the multibody identification problem, we derive the equations of motion. In this section we first show that classical Lagrangian methods of dynamics are reasonable to model the dynamics of the human spine. Generally, we can describe the configuration of a rigid body in space by using six coordinates. Three of them describe the body translation; three other coordinates define the orientation of the body. In two-dimensional analysis, the configuration of a rigid body system can be identified by using three coordinates; two of them define the translation of the points on the body, and one coordinate defines the orientation of the body with respect to a selected inertial frame of reference. Consider one rigid body in figure 2-1. From the viewpoint of rigid body kinematics, the global position of arbitrary point P i on the rigid body I can be defined as i ci i rR g =+ (2-1) where i r is the global position of point P i c i R is the global position vector of the origin O i

PAGE 20

9 Figure 2-1 Two vertebral body system of the body I in the fixed reference frame B i, and i g is the position vector of point P i with respect to O i. So the whole body configuration can be completely defined, provide the vector components of c i R i g in the right side of the equation are known. The vector i r and c i R define in the global coordinates, and i g is defined in the body fixed local coordinates. We need to express the vector i g in terms of the global systems. This is done by defining the orientation transformation of the body fixed coordinates with respect to the global frame of reference. For example, the model depicted in figure 2-1 is comprised of two vertebral bodies, denoted as I and J. Each body has an associated body fixed reference frame B i and B j whose orientation is defined via a change of basis expressed in terms of an orientation matrix ( ) [ ] i L q as follows. ( ) [ ] ( ) [ ] e L b e L b j j i i q q = = ( ) [ ] = 3 2 1 3 2 1 e e e L b b b i i i i q (2-2) I J c i R

PAGE 21

10 In this equation, k i b is the k th basis vector defining body i k e is a basis of the inertial frame E, and i q is the angle measures the rotation of the body frame relative to the inertial frame. The location of the center of mass of body B i and B j is given by 12 i c iii i e X R X e Ye Y = += (2-3) 12 j c jjj j e X R X e Ye Y =+= (2-4) respectively. Similar expression for velocity and acceleration of the arbitrary point on the rigid body can be obtained. Thus the position, velocity, and acceleration of an arbitrary point on the rigid body can be written in terms of these coordinates. Equation (2-1) implies that the general motion of a rigid body is equivalent to a translation of one point, usually the original point of the body fixed frame, plus a rotation. From the viewpoint of kinetics, we can describe the kinetic energy of a rigid body as the sum of translational kinetic energy and the rotational energy. The total kinetic energy of a system is the sum of its components kinetic energy. In the absence of any explicit constraint between the two bodies, the kinetic energy of the system is easily shown to be 2 2 2 2 22 22 2 222 1 1 11 2 2 22 1 2 i i i i i jj j jj i i i i ii ii i jjj j j j j j j T m X Y J m X YJ X m Y m J X Y XY m X m Y J qq q qq q = + + + ++ = (2-5)

PAGE 22

11 where i m is the mass of body I and i J is the inertia of body I about the 3-axis passing through its center of mass. and XY q & && are the translational and rotational velocity of the body I They are all defined in the global coordinates system. The same description can be applied to the other body J Figure 2-2 Two-body example illustration If we sought to enforce a conventional revolute joint at the point denoted as s in figure (2-2), and assume it is frictionless, so no work is done by this constraint, the kinematics constraint could be expressed as [ ] [ ] + = = F ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 2 1 2 1 s s l Y X s s l Y X s R s R s j j T j j j i i T i i i P j P i g g q g g q (2-6) In this 2D and two-body system, ,1 i g and ,1 j g are defined as the translational position of the arbitrary points P i and P j in X direction, and ,2 i g and ,2 j g as the translational position of the arbitrary points P i and P j in the Y direction. Both of them are the functions depending on the position of the revolute joint s. 1 e 2 e 1 b 2 b 2 b 1 b p i p j ( s) R i R j S

PAGE 23

12 One method to solve this problem is to use Lagrange multiplier equation. A classical Lagrange multiplier formulation of the equations would be k M i i k i k k k Q q q V q T q T dt d = F + + = 1 l & N k L 1 = (2-7) where k Q are the generalized forces, k q are the generalized coordinates. By defining potential energy ) ( ) ( ) ( 2 1 ds s s V T m m G F F = (2-8) instead of using Lagrange multiplier, we could use the standard Lagranges equation to solve an unconstrained problem. k k kkk V d T TV Q d t q qqq m -+ += & 1 kN = L (2-9) where m is the probability measure. In this case, the kinetic energy is the same as that in the classical Lagranges multiplier formulation. We just need to calculate the potential energy part ( ) 1 ( ) () () 2 1 ( ( ) () ( ) ( ) () ( )) 2 T TT V ss ds qq s s d s s s ds qq m m mm G GG = FF F = F F +F (2-10) Alternatively, this equation can be written in terms of the summation convention as F F + F F = F F = G G G ) ( ) ( ) ( ) ( ) ( ) ( 2 1 ) ( ) ( ) ( 2 1 ds q s s ds s q s ds s s q q V j j j i j j j i i m m m m

PAGE 24

13 ) ( ) ( ) ( ds s q s q V j i j i m m G F F = (2-11) In this form, it would be very difficult, if not impossible, to solve the above relaxed formulation. The probability measure m in equation (2-10), (2-11) could be an element of an infinite dimensional space. To make the problem tractable, we employ a simple approximation of the probability measure m ( ) ( ) = = N l S l l K 1 d m (2-12) In this equation, sl d is the dirac measure concentrated at l s That is ( ) 1 0 l sl sA SA otherwise = (2-13) From a mathematical standpoint, we approximate m by k m which is the weighted sum of dirac functions concentrated at m k s k L 1 = From a practical standpoint, this approximation can be thought as seeking a distribution of spring elements along the curve between the two rigid bodies. The expression in equation (2-11) can be simplified to ( ) ( ) ( ) ( ) 1 1 N jl l jl l ii T N l ll l VS KS qq S KS q m = = F =F F =F (2-14) Recall the constraint function F is given as [ ] [ ] + = D F ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 2 1 2 1 s s l Y X s s l Y X s R s R s j j T j j j i i T i i i P j P i g g q g g q (2-15)

PAGE 25

14 Recall orientation transformation of the basis, the transform matrix [( )] L q can be expressed in the two dimensional plane as follows, = 3 2 1 3 2 1 1 0 0 0 cos sin 0 sin cos e e e b b b i i i i i i i q q q q (2-16) = 3 2 1 3 2 1 1 0 0 0 cos sin 0 sin cos i i i i i i i b b b e e e q q q q (2-17) Plug the basis transform equations (2-16) (2-17) into the constraint equation, the constraint vector F becomes + = F ) ( ) ( cos sin sin cos ) ( ) ( cos sin sin cos ) ( 2 1 2 1 s s Y X s s Y X s j j j j j j j j i i i i i i i i g g q q q q g g q q q q (2-18) So, we can describe the constraint function in the inertial frame E, and then calculate the Jacobian j i q F in equation (2-14) by definition, F F F F F F F F F F F F = F j j j i i i j j j i i i j i Y X Y X Y X Y X q s q q q q 2 2 2 2 2 2 1 1 1 1 1 1 ) ( (2-19) where the generalized coordinates j q are defined as { } j j j i i i Y X Y X q q q = (2-20) They are the translational and rotational basis in a two-dimensional coordinates system as we mentioned above. If we plug the constraint function into equation (2-19), differentiate the vector components separately, we get the Jacobian matrix

PAGE 26

15 + + = F ) ) sin( ) cos( ( 1 0 ) ) sin( ) (cos( 1 0 ) ) cos( ) (sin( 0 1 ) ) cos( ) sin( ( 0 1 2 1 2 1 2 1 2 1 j j j j i i i i j j j j i i i i j i q g q g q g q g q g q g q g q g q (2-21) At this point, we have derived both the kinetic energy and the potential energy for this two-body system using Lagranges Equations k k kk V d TT Q d t q qq m -+= & (2-22) The final set of equations for the relaxed formulation for the two vertebral bodies depicted in figure (2-1) is 1 2 1 ,2 1 2 1 ,2 10 01 ( sin ( ) cos ( )) (cos ( ) sin ( )) 10 01 (sin ( ) cos ( ) ) ( cos ( ) sin ( )) i i i i i i i i i i i ii i l j j j j j j j j j j j jj j m X m Y J K m X m Y J q g q g q g qg q q g q g q g qg q -+ + -+ && && && && && && ( ) ( ) ( ) ( ) 1 1 2 ,1 ,1 3 ,2 ,2 4 5 6 cos ( ) sin () cos ( ) sin () sin ( ) cos () sin ( ) cos () N l jjj jl il iii j jj jl il i ii Q Q X s s XQ Y s s YQ Q Q qq g g qq qq g g qq = + -= (2-23) 2.2 Identification Problem in the Two-Dimensional Example The system identification issue in our research is mainly to derive a mathmatical model for the human spine based on clinically measured data. By essentially adjusting model parameters within this given model, we try to make the output of the simulation system coincide with the experimental data as well as possible. In the above section, we have derived the governing equations of the two-body system. But we wish to know how

PAGE 27

16 well this equation can describe the two-body example. In other words, we want to find the a set of the suitable probability measures that can fit our model with the experimental data. As we know, the mathematical model describes the relationship between the input signals and output signals. Considering the example we have, the inputs would be the different forces and moments applied on the vetebral bodies, while the output would be the vetebral bodys orientation, position and movement. The outputs are affected not only by the system characterization, but also by the input signals. We denote the inputs by u(t) and outputs by y(t) respectively. The relationship is depicted in figure 2-3. All these signals are functions of time, and the value of the input at time t will be denoted as u(t) Often, in the identification context, only discrete-time points are considered, since in the clinical experiment, the instruments used in the experiments typically record the signal just at discrete time instants. The identification of spine model is based on these experimental input and output data. Our human spine system identification is comprised of the following steps: Design the clinical experiment and collect input-output data from the experimental procedures. This is the first stepping stone in going from an experiment to the theratical model. However, these experimental data often need to be polished by removing trends and outliers, and selecting useful information from the original data. Many research has been done in this field, as noted in the literature we reviewed in the introduction.

PAGE 28

17 Select and define a spine model mathematically to respresent the experiment. There are a set of candidate system descriptions, within which a model is to be found. Compute the best spine model characteristic parameters according to the input-output clinical experimental data and a given criterion of fit. Then, examine the obtained models properties. In figure 2-3, the output y(t) is the system response to the input signal u(t), while d(t) is the desired output. In spine modeling, d(t) can be considered as the experiment data we obtained from the clinical experiment. By comparing the output of the given model with the desired output, we can adjust the characteristic parameter within the mathematical model, then minimize the error e(t) between the desired output d(t) and the output of identification model. Figure 2-3 Identification procedure illustration In our simple example, the identification problem is a minimization problem over the probability measure space. In this two-body vertebral segment system, we need to find a stiffness distribution that could reconstruct the clinical experimental data. Suppose we have measurements (,) yy & (which stand for the displacement and velocity measurements in our case). Then, we define the performance function as MODEL COMPARE u(t) y(t) d(t) e(t)

PAGE 29

18 1 2 2 2 12 ( ) (( ) ,( ) ) (( ) ,( )) Z Z J y H q q yH qq m m m mm = +& && (2-24) The functions 1 H and 2 H estimate the experimental data from the model for a given m So the the performance function defines the error E as described in figure 2-3. In this work we try to identify the stiffness coefficients of the system that minimize the performance equation. Recall the governing function we have derived in the first session. Its difficult to get a analytical solution in the form of 1 H and 2 H that we can substitute it into equation (2-24) directly. Numerical analysis is widely used today to solve this problem. It can provide the numerical answer while the analytical answer is difficult to obtain. There are two major computing procedures employed in solving this problem. One is to find a numerical solution to the ordinary differential equations arising from solving the Lagranges equations. Another is to use a least-squares method to calculate the spring stiffness coefficients. We obtained the governing equation of motion in the form of equation (2-23). This is a set of second-order ordinary differential equations coupled to each other. However, these second-order equations can be reduced to a system of simutaneous firstorder equations. Lets express one of the second-order differential equation as a example, 2 2 ( ,) d X dX ftX d t dt = (2-25) The initial values can be described as 00 () X tX = and 00 () X tX = We can convert this to a pair of first-order equations by the simple expedient of defining the derivative as a second function.

PAGE 30

19 1 dX Y dt = 00 () X tX = (2-26) 1 1 ( ,,) dY ftXY dt = 1 00 () Y tX = (2-27) This pair of first-oder equation is equivalent to the original equations. Similar equivalent first-order equations can be obtained for the other second order differential equations. So, equation (2-26) can be coverted into { } 1 1 2 2 3 3 1 4 4 1 5 5 2 6 6 3 7 4 8 5 9 6 10 11 12 1 1 1 (,) 1 1 0 1 0 0 0 1 0 1 1 1 1 1 i i N i l l l j j j m m Y Q Y Q J Y Q K f qt Y Q m Y Y Q Y Y Q m Y Y Y Y m Y Y Y Y Y Y = =+ rr 0 (2-28) where the vector i Y contains six generalized coordinates q r namely ( ) ,, ,, ,, ,, ,, ii i jj j ii i jjj X Y X Y X Y XY qqqq && & & && And, { } 1 (,) N l l l K fqt = rr stands for the Jacobian matrix j i q F multipied by the constraint vector F as follows.

PAGE 31

20 ( ) ( ) 1 2 1 ,2 1 2 1 ,2 ,1 ,2 10 01 ( sin ( ) cos ( ) ) (cos ( ) sin ( )) 10 01 (sin ( ) cos ( ) ) ( cos ( ) sin ( )) cos ( ) sin () sin ( ) cos () i i i i i i ii j j j j j j jj il iii il i ii s X s Y q g q g q g qg q g q g q g qg g qq g qq -+ -+ + ( ) ( ) ,1 ,2 cos ( ) sin () sin ( ) cos () jjj jl j jj jl X s Y s qq g qq g -(2-29) At this point, we can apply numerical methods to solve these pairs of first-order equations. Another numerical issue is related to the identification procedure, where we try to fit the clinical experiment data into the given mathematical model. Obviously, the problem involved is one of fitting nonlinear data. The overall algorithm describing the numerical solution procedure is described in figure 2-4. First, an initial guess of the probability measure vector 0 m is made. That is, the initial value of the spring stiffness vector 0 k is given. By using the ODE solver in Matlab, we can obtain the numerical approximation at the output data. We then compare these data with the clinical experiment data. By checking the performance function, or error measure to a certain threshold, we decide to stop or go on to the next iteration step. The last results in this iterative process are estimates of the identified characterization parameters. That is, we obtain estimates of the spring stiffness distribution between the two vertebral bodies. Detailed numerical simulation procedures and simulation results will be explained in the following chapters.

PAGE 32

21 Figure 2-4 Numerical simulation procedure illustration Experiment Data Least Square Controller Minimum? ODE Solver Approximation Data Stop

PAGE 33

22 CHAPTER 3 NUMERICAL EXPERIMENT In this chapter, we first give an introduction to common terms for the simulation tool Promechanica. Then, two numerical experiments are discussed in detail. The first one is a forward numerical simulation, and the second is an inverse dynamics identification experiment. 3.1 Multibody Dynamics Simulation ToolsPromechanica In this chapter, we will describe the development of a computer model of the human spine using the commercial dynamics program Promechanica. It has a particular high level language to allow modeling of rigid bodies, joints and constraints. It is possible using this language to simulate the movement of mechanisms. Additionally it is possible to measure forces and moments at required points. Since these quantities will be used in the remainder of the text, we introduce the basic terms within the programming language. Terminology: Body: Representation of a rigid body within the computer model. Connection: Joint used to connect two or more bodies together (pins, cylinder, ball, bearing, universal joint, sliders, free, weld, gimbals, planar and six degree of freedom). Load: Applies loads to points, joint axis or between two points on different parts (gravity, fixed loads, moments, static and dynamic friction).

PAGE 34

23 Driver: Used to specify the acceleration, velocity or position of a part or part joint axis (ramps or constant functions, cosine, table or combinations of these. Measure: Used to monitor changes in model parameters during the analysis. Mass primitive: A geometric shape (brick, cone, cylinder, plate, and sphere) which can be added to a part, changing its mass, center of mass and moments/products of inertia. The geometry of bodies are created by connecting points specified in the body fixed local coordinate system (LCS). Mass and inertial properties can then be assigned to the bodies either explicitly by the user or by using mass primitives and material defined within the program. These bodies can be connected to each other and to a ground point using a variety of joints to form a mechanical system. Drivers and measures can be assigned to each joint axis if so required. Surface geometry allows a more realistic visualization of the model. Surface is attached automatically to mass primitives and can be attached manually to bodies. The mass and inertial properties of bodies are entered explicitly by the analyst. Also there are some other capabilities to tailor the kinematic and dynamic aspects of the mechanism design. One important capability is the custom load, which is created outside of Promechanica Motion by the analyst. This is convenient for a user to model a more complex load. To analyze the stiffness characteristics of the human spine and simulate its movements, we would like to describe the joint characteristics between two vertebral bodies as a distribution of springs. There are no such joints can be used directly inside the Promechanica from it lists of build-in entities. The custom load is good way that we can define our specific model by using analyst-developed loads.

PAGE 35

24 3.2 Forward Simulation with Custom Load Recall the main procedure in the human spine modeling as described in Fig 3-1. First, we obtained the data from the clinical experiments. By using an approximation method, we can identify the mechanical characteristic as parameters in the model we construct. Then, we input these identified parameters into Promechanica, to simulate the motion of the experiment. Figure 3-1 Human spine modeling procedure If we consider the approximation procedure as the inverse dynamics problem, the Promechanica simulation provides the direct dynamics problem. The direct dynamics problem, being the inverse of inverse dynamics solutions, should provide a way to check the inverse dynamics problem.[DLC1997] That is, if the approximate solution, calculated using our mathematical model in chapter 2, is used as input to Promechanica simulation with the same model, the outputs of the simulation should match the data as Experiment Data Approximation Characteristic Parameter Promechanica Simulation

PAGE 36

25 the input to the identification procedure. If instead, the output of simulation and the input to the inverse dynamics problem are different, then an error exists either in the fundamental formulation or the implementation of either the inverse dynamics or direct dynamics solution techniques. A forward numerical experiment was designed for this consistency check. Figure 3-2 A two vertebral body model built in Promechanica Consider the two-body system example we used in driving the governing equations. We construct the same two-body system example in Promechanica Motion by using published data together with information on the spine geometry imported from a commercially available source. Figure 3-2 shows the two vertebral bodies with geometry attached. In this forward simulation, we simply compare two mathematical models. One is the Matlab model based on the governing equations derived in chapter2. Another is the one constructed in Promechanica. The same external forces are applied to identified model and Promechanica model. Also, we input the same mechanical parameters to these two models, like the same number of springs, the same attachment position, and the

PAGE 37

26 same stiffness probability measurements. The outputs of both systems are time responses to the same external forces. From the clinic viewpoint, the Promechanica simulation output can be considered as the clinical experimental data. By comparing the outputs of these two models, we can find if the Matlab model matches the Promechanica model. Figure 3-3 Generalized forces in two-dimensional plane The external forces were added along the three generalized coordinates in this two-dimensional example. They are generalized forces Q1 and Q2 along the Z and Y translational direction, and one generalized moment Q3 about the X-axis, as shown in the above figure 3-3. The generalized forces could have different expressions. For example, in a simulation using 8 discrete springs, we use where f l f h stand for the low frequency part and high frequency part in the external forces. In a 16-spring rotational movement simulation, we use y x z Q 1 Q 2 Q 3 ( ) ( ) q p q p + + + = h l f A f A Q 2 cos 2 cos 2 1 1 ( ) ( ) q p q p + + + = h l f B f B Q 2 cos 2 cos 2 1 2 ( ) ) 2 cos( 2 cos 2 1 3 q p q p + + + = h l f C f C Q (3-1)

PAGE 38

27 With different coefficients A i B i C i D i E i we can give different force inputs. After the simulations have been run, we compare the numerical results from both models as shown figures 3-4 through 3-15. 2 1 1 1 1 1 2 1 1 + + + + = x E x D C x B x A Q 2 2 1 2 2 2 2 2 2 + + + + = x E x D C x B x A Q 2 3 1 3 3 3 2 3 3 + + + + = x E x D C x B x A Q (3-2)

PAGE 39

28 8-spring numerical simulation comparison 0 50 100 150 200 250 -1.5 -1 -0.5 0 0.5 1 1.5 Figure 3-4 Y-translational trajectory numerical results from Promechanica 0 50 100 150 200 250 -1.5 -1 -0.5 0 0.5 1 1.5 Figure 3-5 Y-translational trajectory numerical results from Matlab model

PAGE 40

29 0 50 100 150 200 250 37 37.5 38 38.5 39 39.5 40 Figure 3-6 Z-translational trajectory numerical results from Promechanica 0 50 100 150 200 250 37 37.5 38 38.5 39 39.5 40 Figure 3-7 Z-translational trajectory numerical results from Matlab model

PAGE 41

30 0 50 100 150 200 250 -1.5 -1 -0.5 0 0.5 1 1.5 Figure 3-8 X-rotational trajectory numerical results from Promechanica 0 50 100 150 200 250 -1.5 -1 -0.5 0 0.5 1 1.5 Figure 3-9 X-rotational trajectory numerical results from Matlab model

PAGE 42

31 16-spring numerical simulation comparison 0 50 100 150 200 250 -6 -4 -2 0 2 4 6 8 10 12 x 10 -5 Figure 3-10 Y-translational trajectory from Promechanica model 0 50 100 150 200 250 -6 -4 -2 0 2 4 6 8 10 12 14 x 10 -5 Figure 3-11Y-translational trajectory from Matlab model

PAGE 43

32 0 50 100 150 200 250 38.299 38.3 38.301 38.302 38.303 38.304 38.305 38.306 38.307 38.308 Figure 3-12 Z-translational trajectory from Promechanica model 0 50 100 150 200 250 38.299 38.3 38.301 38.302 38.303 38.304 38.305 38.306 38.307 38.308 38.309 Figure 3-13 Z-translational trajectory from Matlab model

PAGE 44

33 0 50 100 150 200 250 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 x 10 -4 Figure 3-14 X-rotational trajectory from Promechanica model 0 50 100 150 200 250 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 Figure 3-15 X-rotational trajectory from Matlab model

PAGE 45

34 In many different simulation cases, the numerical trajectory from the Matlab model matches the results from the Promechanica Motion simulation, as shown in the above figures. At last, we used the forward simulation capability to solve an identification problem. That is, we use Promechanica Motion to get an output trajectory. This data can be considered as the clinical, experimental data. Then, the experimental trajectory data is supplied to the identification procedure in Matlab. We can identify the spring stiffness distribution that best represents the trajectory obtained from Promechanica. This is the same as the identification procedure we discussed in the chaper 3. We compare the orginal input spring stiffness coefficients use to generate the experimental trajectory and the identified stiffness coefficients in the following tables. Table 3-1 8 Springs identification Spring 1 2 3 4 5 6 7 8 Input 0.2 0.2 0.1 0.3 0.15 0.12 0.08 0.22 Identify 0.1962 0.1989 0.1045 0.29 0.152 0.1227 0.0818 0.2239 Table 3-2 16 Springs identification Spring 1 2 3 4 5 6 7 8 Input 0.3 0.15 0.12 0.08 0.22 0.25 0.16 0.3 Identify 0.2913 0.1413 0.11 0.0743 0.2151 0.2506 0.1608 0.3048 Spring 9 10 11 12 13 14 15 16 Input 0.1 0.07 0.09 0.18 0.2 0.2 0.2 0.1 Identify 0.1052 0.0776 0.0986 0.19 0.2099 0.1915 0.1919 0.1071

PAGE 46

35 When we input the identified spring stiffness into the simualtion software Promechanica, we can visualize of the verte bral body movent. Figure 3-16 Translation movement simulation Figure 3-17 Rotation movement simulation

PAGE 47

36 3.3 Identification of Pin Joint Position Using Inverse Simulation The last examples considered trajectory matching problems. That is, we sought the distributions to minimize the difference between observed data and calculated outputs. One of the critical properties of a relaxed formulation is that it should be able to represent the classical results as a special case. A numerical experiment is used to show this is the case. Suppose we use a pin joint to connect two vertebral bodies. We want to know if we can use the relaxed formulation to approximate this kind of the classical mechanical joint. Obviously, this is an identification problem. Figure 3-18 vertebral body connected by pin Similar to the previous numerical experiment, we use the simulation tool Promechanica to construct a two-dimensional, two-body system. By importing a two vertebral segment system, we use a pin joint to connect the vertebrae together as shown in figure 3-18. The center of the pin is at a certain point between the two bodies in the sagittal plane. Since our spine model is a three-dimensional geometry model, we use 2 e 3 e 2 i b 2 i b 2 j b 2 j b A B

PAGE 48

37 drivers to apply constraints on three coordinates (one translational and two rotational), and thus simplify the three dimensional system to a two dimensional system. After the dynamical simulation, we get the trajectories of the movement for the system connected via the classical pin joint. Then we want to find a distribution of springs that can replace the pin joint and can obtain the same system outputs. The attachment points for each spring are fixed on the surface of each vertebral body. The interval distance between these attachment points is equal. We can make the distance finer by adding more locations to connect the vertebral bodies. By using the identification method as we discuss in the chapter three, we can identify the stiffness coefficient value of each spring.

PAGE 49

38 The 2 spring forces case 10 15 20 25 30 35 40 0 0.5 1 1.5 2 2.5 x 10 -14 Discrete Probability Density Function Discrete Probability Density Function 2 Springs Case Position of Spring Attaching Point marked by '*' Probability Density Distribution of Springs data1 data2 data3 data4 A B Figure 3-19 Probability density distribution of spring stiffness 10 15 20 25 30 35 40 0 1 2 3 4 5 x 10 -14 Cumulative Distribution of Spring Stiffness Cumulative Distribution of Spring Stiffness 2 Springs Case Position of Springs Attaching Points Cumulative Distribution of Spring Stiffness A B Figure 3-20 Cumulative distribution of spring stiffness

PAGE 50

39 The 4 spring forces case 10 15 20 25 30 35 40 0 0.2 0.4 0.6 0.8 1 1.2 1.4 x 10 -3 Discrete Probability Density Function Discrete Probability Density Function 4 Springs Case Position of Spring Attaching Point Marked By '*' Probability Density Distribution of Springs A B Figure 3-21 Probability density distribution of spring stiffness 10 15 20 25 30 35 40 0 0.5 1 1.5 2 2.5 3 x 10 -3 Cumulative Distribution of Spring Stiffness Cumulative Distribution of Spring Stiffness 4 Springs Case Position of Springs Attaching Points Cumulative Distribution of Spring Stiffness A B Figure 3-22 Cumulative distribution of spring stiffness

PAGE 51

40 The 8 spring forces case 10 15 20 25 30 35 40 0 0.5 1 1.5 2 2.5 x 10 -3 Discrete Probability Density Function Discrete Probability Density Function 8 Springs Case Position of Spring Attaching Point Marked By '*' Probability Density Distribution of Springs A B Figure 3-23 Probability density distribution of spring stiffness 10 15 20 25 30 35 40 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 -3 Cumulative Distribution of Spring Stiffness Cumulative Distribution of Spring Stiffness 8 Springs Position of Springs Attaching Points Cumulative Distribution of Spring Stiffness A B Figure 3-24 Cumulative distribution of spring stiffness

PAGE 52

41 The 16 spring forces case 10 15 20 25 30 35 40 0 0.5 1 1.5 2 2.5 x 10 -3 Discrete Pobability Density Function Discrete Pobability Density Function 16 Springs Case Position of Spring Attaching Point Marked By '*' Probability Density Distribution of Springs A B Figure 3-25 Probability density distribution of spring stiffness 10 15 20 25 30 35 40 0 1 2 3 4 5 6 7 x 10 -3 Cumulative Distribution of Spring Stiffness Cumulative Distribution of Spring Stiffness 16 Springs Case Position of Springs Attaching Points Cumulative Distribution of Spring Stiffness A B Figure 3-26 Cumulative distribution of spring stiffness

PAGE 53

42 The 32 spring forces case 10 15 20 25 30 35 40 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 -3 Discrete Probability Density Function Discrete Probability Density Function 32 Springs Case Probability Density Distribution of Springs Position of Spring Attaching Point Marked By '*' A B Figure 3-27 Probability density distribution of spring stiffness 10 15 20 25 30 35 40 0 0.002 0.004 0.006 0.008 0.01 0.012 Cumulative Distribution of Springs Stiffness Cumulative Distribution of Springs Stiffness 32 Springs Case Position of Springs Attaching Points Cumulative Distribution of Springs Stiffness A B Figure 3-28 Cumulative distribution of spring stiffness

PAGE 54

43 The 64 spring forces case 10 15 20 25 30 35 40 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 -3 Discrete Probability Density Function Discrete Probability Density Function 64 Springs Cases Position of Spring Attaching Point Probability Density DiDistribution of Springs A B Figure 3-29 Probability density distribution of spring stiffness 10 15 20 25 30 35 40 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 Cumulative Distribution of Spring Stiffness Cumulative Distribution of Spring Stiffness 64 Springs Case Position of Springs Attaching Points Cumulative Distribution of Spring Stiffness A B Figure 3-30 Cumulative distribution of spring stiffness

PAGE 55

44 The 128 spring forces case 10 15 20 25 30 35 40 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 -3 Discrete Probability Density Function Discrete Probability Density Function 128 Spring Cases Position of Spring Attaching Point Probability Density Distribution of Springs A B Figure 3-31 Probability density distribution of spring stiffness 10 15 20 25 30 35 40 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Cumulative Distribution of Springs Stiffness Cumulative Distribution of Springs Stiffness 128 Springs Case Position of Springs Attaching Points Cumulative Distribution of Springs Stiffness A B Figure 3-32 Cumulative distribution of spring stiffness

PAGE 56

45 The 256 spring forces case 10 15 20 25 30 35 40 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Discrete Probability Density Function Discrete Probability Density Function 256 Springs Case Position of Spring Attaching Point probability Density Distribution of Springs A B Figure 3-33 Probability density distribution of spring stiffness 10 15 20 25 30 35 40 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Cumulative Distribution Function of Spring Stiffness Cumulative Distribution Function of Spring Stiffness 256 Springs Case Position of Springs Attaching Points Cumulative Distribution of Spring Stiffness A B Figure 3-34 Cumulative distribution of spring stiffness

PAGE 57

46 If we consider the spring stiffness coefficients as the probability function, the value of stiffness can be normalized between 0 and 1. By plotting all the spring coefficients probability distribution and cumulative distribution are obtained. We can see the identified numerical results clearly. With a number of springs connecting the two bodies, the lager stiffness coefficients are more likely move to the position where the pin joint exactly is. 8 20 32 44 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Spring Stiffness Distribution Spring Stiffness Distribution 2,4,8,16,32,64,128,256 springs case Position of Springs Attaching Points Spring Stiffness Distribution Figure 3-35 Probabilty density distribution of different spring numbers 10 15 20 25 30 35 40 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Cumulative Distribution of Spring Stiffness Comparation Cumulative Distribution of Spring Stiffness Comparation 2,4,8,16,32,64,128,256 springs cases position of springs attaching points cumulative distribution of spring stiffness A B Figure 3-36 The cumulative distribution of different spring numbers

PAGE 58

47 Instead of the pin joint between the two vertebral boides, we input these different spring stiffness distribution into a Promechanica Motion model with the same external forces, and simulated the system to get the resulting trajectories of the system. By comparing the error between the original system and the identified system, we found that better correlation is achieved using a finer distribution vertebral segment. The following figures show the relationship between the approximation error and different number of the springs. 0 10 20 30 40 50 60 70 23.5 24 24.5 25 25.5 26 26.5 27 27.5 28 28.5 Figure 3-37 Error vs. different number of springs The identification above shows that using a finer spring distribution yields a smaller approximation error between identified model and the ideal model. Also, the position of the pin joint can be identified more clearly. This simulation illustrates an important fact: the joint relaxation technique includes, as a special case, the classical pin joint constraint.

PAGE 59

48 CHAPTER 4 SCOLIOSIS MODEL DATABASE CONSTRUCTION In chapter 4, we summarize the development of a structural and geometric database that will support the future study of scoliosis. We begin with a discussion of the spine stiffness database construction. Finally, we show an example of a two-dimensional vertebral segment whose stiffness characteristics are modeled using the database. The tables in that constitute the database can be found in the Appendix. 4.1 Geometry In order to construct a mathematical model of the human spine, we need to construct a three-dimensional scoliosis model that can be used in Promechanica geometry. From the literature, we have found many types of human spine models that have been published. Based on these published models, we have built an electronic database. A three-dimensional spine model and coordinate system is constructed for use with the dynamics simulation tool Promechanica. This simulation package is selected because scoliosis is a three-dimensional deformity. It makes the spine deform in a way that requires a six DOF model. However, by adding the corresponding constraints and drivers, we can still simplify it to a two-dimensional model in some specific cases. The spine can be idealized as a collection of rigid bodies interconnected by some joints or elements. The vertebrae are modeled as rigid bodies, since their deformation is negligible compared to the intervertebral tissue. The geometry of an individual vertebral body can be constructed from the published data by using the Promechanica geometry

PAGE 60

49 attachment function. In the global coordinates system, the sacrum will be considered as the fixed ground of the whole model. Each vertebral body above it, from L5 to T1, can be imported individually by using the subassembly function. In each subassembly system a local coordinate system is defined. Every vertebral body is imported into its corresponding subassembly system. We can describe the geometry characteristics in this local subassembly system. Then, we import subassemblies into the global system one by one, according to the published data defining the shape of scoliosis. The global coordinates of each primary node in the subassembly system can be computed by using the transform matrix we derived in chapter 3. The orientation values in the sagittal plane are derived from the published data. By inputting the mechanical properties, like the mass, and other connecting joint constraints, a spine model is constructed in Promechanica. In table 4-1, some geometry data are given that characterize a spine that is scoliotic.

PAGE 61

50 Table 4-1 Scoliosis curve data Vertebra Sagittal Orientation (deg) Vertebra Height (mm) Inferior Disc Height (mm) L5 -8 22 15.7 L4 -2 28 14 L3 -11 26 12.2 L2 -21 26 11.4 L1 -21 24 10 T12 -18 24 8.4 T11 -17 20 6.8 T10 -17 20 5.1 T9 -15 20 4.7 T8 -8 20 4.5 T7 0 20 4.0 T6 3 20 3.2 T5 4 20 2.6 T4 9 20 2.2 T3 15 18 2.7 T2 20 18 3.1 T1 24 16 4.5

PAGE 62

51 Table 4-2 Coordinates of origin nodes Vertebra globalX globalY globalZ L5 0 0 0 L4 0 -31.19 37.46 L3 0 -32.05 77.73 L2 0 -41.45 112.73 L1 0 -51.50 145.50 T12 0 -64.00 175.40 T11 0 -68.02 202.95 T10 0 -76.32 226.65 T9 0 -83.62 250.35 T8 0 -88.72 274.15 T7 0 -90.72 298.05 T6 0 -90.72 321.25 T5 0 -89.52 343.75 T4 0 -87.92 365.95 T3 0 -81.28 387.22 T2 0 -75.78 407.63 T1 0 -66.04 427.90

PAGE 63

52 Table 4-2 shows the three-dimensional coordinates of the original nodes of each local subassembly systems. The subassembly systems are imported into the global system by pointing their original nodes and orientations. Using this information the global coordinates of all the individual points can be calculated. The relative position of any arbitrary points on the superior and inferior vertebral bodies is then found in this global system via the transformation matrix between the local coordinate systems and global coordinate system. A spine model form T1 to L5 is generated in Promechanica, as shown in Fig 4-1 Fig 4-1 Spine Model in Promechanica Motion To make a simpler geometry model, we pick up 17 primary nodes, which could describe the basic geometric outline of vertebra. And the body of vertebra is simplified to a cylinder. The other primary nodes can be picked up on the pedicle, superior articular

PAGE 64

53 process, transverse process, and inferior articular process. Using spline curves to link these nodes, we can get a basic shape of a vertebral body. Fig 4-2 A Vertebral Model with Complex Surface Fig 4-3 A Simple Geometry Vertebral Model

PAGE 65

54 Fig 4-4 Right Lateral View of the Primary Nodes 12 6 3 4 2 1 8 16 15 14 10 17

PAGE 66

55 Fig 4-5 Top View of the Primary Nodes 1 2 6 7 8 11 12 15 5 3 16

PAGE 67

56 Fig 4-6 Bottom View of the Primary Nodes 4 3 6 5 14 10 9 13 17 16 1

PAGE 68

57 Fig 4-7 Simplified Geometry Spine Model A simplified geometry spine model is constructed as shown in Fig 4-7. One objective of this biomechanical simulation is to study the implants that attach to the bone. We have derived the position of the pedicle screws, transverse process hooks, and lumina hooks based on these primary nodes. Because the internal fixation devices used in spinal surgery are not strictly sensitive to their position, we can model the geometry and the position of these implant devices by some simple computations using our biomedical experiment database. The geometric database describing the orientation of the screws and hooks are described in the Appendix.

PAGE 69

58 4.2 Stiffness Many researchers have studied on the stiffness properties of the connective tissues between the vertebral bodies. The literature review showed that different stiffness values have been reported in different clinical experiments. The diversity in results is partly because of the differing experimental methods and instruments used in these experiments. Also, the experimental data will be different when testing different specimens. Among these experiments, Panjabi et al made extensive studies of the nonlinear viscoelastic mechanical behavior of the human spine. The three-dimensional behavior of the spine is presented in the form of average load-displacement curves for different vertebral bodies. After examining these discrete experiment data, we have identified and tabulated a nonlinear load-displacement function to describe the stiffness characteristics of each Fig 4-8 Illustration Stiffness Property in Three-dimensional Coordinates System

PAGE 70

59 vertebral level. By inputting these simplified stiffness curves into the model we constructed in Promechanica, we can simulate the clinical experiment and reconstruct the relationship between the load and vertebrae movements using numerical simulation tools. This simple test illustrates the utility of the database. Fig 4-8 illustrates the three-dimensional coordinates system that is necessary to describe the three-dimensional physical properties of the spine. The origin of the coordinate system is located at the infer posterior corner in the mid-sagittal plane of the vertebral body. The X-axis is directed to the left and is perpendicular to the sagittal plane. The Y-axis is directed superiorly, and the Z-axis is directed anteriorly. The broad arrows show the pure moments and thin, straight arrows show the translations on the three rotational axes and three translational axes respectively. In the following figures, we show the load-displacement curves between each level of the lumbar spine in three-dimensional system. They are from the location between first and second lumbar vertebrae through location between the fifth lumbar and first sacral vertebrae. For each location, the load-displacement relationships are described according to the different type of loads applied on the vertebrae, pure moments of flexion-extension, lateral axial torque, and bilateral bending, respectively.

PAGE 71

60 Load-displacement curves are shown at each of the five locations on the vertebral lumbar due to the application of flexion and extension moments. The main motion, defined as the rotation in the direction of the applied moment, is indicted by the heavier line. The figures also show the coupled motions, defined as all rotations and translations other than the main rotation. -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 Extension Flexion MOMENT (Nm) MOTION (degree or mm) L1-L2 Fig 4-9 Flextion-extension between L1 and L2

PAGE 72

61 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 Extension Flexion MOMENT (Nm) MOTION (degree or mm) L2-L3 Fig 4-10 Flextion-extension between L2 and L3

PAGE 73

62 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 Extension Flexion MOMENT (Nm) MOTION (degree or mm) L3-L4 Fig 4-11 Flextion-extension between L3 and L4

PAGE 74

63 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 Extension Flexion MOMENT (Nm) MOTION (degree or mm) L4-L5 Fig 4-12 Flextion-extension between L4 and L5

PAGE 75

64 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 Extension Flexion MOMENT (Nm) MOTION (degree or mm) L5-S1 Fig 4-13 Flextion-extension between L5 and S1

PAGE 76

65 Load-displacement curves are next shown at each of the five locations on the lumbar spine due to the application of right and left axial torques. The main motion, defined as the rotation in the direction of the applied moment, is indicted by the heavier line. On the following figures, also shown are the coupled motions, defined as all rotations and translations other than the main rotation. -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -4 -3 -2 -1 0 1 2 3 4 Right Torque Left Torque MOMENT (Nm) MOTION (degree or mm) L1-L2 Fig 4-14 Right and Left Axial Torque between L1 and L2

PAGE 77

66 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -4 -3 -2 -1 0 1 2 3 4 Right Torque Left Torque MOMENT (Nm) MOTION (degree or mm) L2-L3 Fig 4-15 Right and Left Axial Torque between L2 and L3

PAGE 78

67 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -4 -3 -2 -1 0 1 2 3 4 Right Torque Left Torque MOMENT (Nm) MOTION (degree or mm) L3-L4 Fig 4-16 Right and Left Axial Torque between L3 and L4

PAGE 79

68 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -4 -3 -2 -1 0 1 2 3 4 Right Torque Left Torque MOMENT (Nm) MOTION (degree or mm) L4-L5 Fig 4-17 Right and Left Axial Torque between L4 and L5

PAGE 80

69 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -4 -3 -2 -1 0 1 2 3 4 Right Torque Left Torque MOMENT (Nm) MOTION (degree or mm) L5-S1 Fig 4-18 Right and Left Axial Torque between L5 and S1

PAGE 81

70 Load-displacement curves below correspond to each of the five locations on the lumbar spine due to the application of left and right lateral bending moments. The main motion, defined as the rotation in the direction of the applied moment, is indicted by the heavier line. The figures also show the coupled motions, defined as all rotations and translations other than the main rotation. -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -8 -6 -4 -2 0 2 4 6 8 Left Bending Right Bending MOMENT (Nm) MOTION (degree or mm) L1-L2 Fig 4-19 Left and Right Bending Moment between L1 and L2

PAGE 82

71 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -8 -6 -4 -2 0 2 4 6 8 Left Bending Right Bending MOMENT (Nm) MOTION (degree or mm) L2-L3 Fig 4-20 Left and Right Bending Moment between L2 and L3

PAGE 83

72 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -8 -6 -4 -2 0 2 4 6 8 Left Bending Right Bending MOMENT (Nm) MOTION (degree or mm) L3-L4 Fig 4-21 Left and Right Bending Moment between L3 and L4

PAGE 84

73 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -8 -6 -4 -2 0 2 4 6 8 Left Bending Right Bending MOMENT (Nm) MOTION (degree or mm) L4-L5 Fig 4-22 Left and Right Bending Moment between L4 and L5

PAGE 85

74 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -8 -6 -4 -2 0 2 4 6 8 Left Bending Right Bending MOMENT (Nm) MOTION (degree or mm) L5-S1 Fig 4-23 Left and Right Bending Moment between L5 and S1

PAGE 86

75 4.3 A Simple Demonstration of the Stiffness Database After the spring stiffness database for the spine (Lumbar and thoracic) has been built, we simply approximate the stiffness coefficients as the nonlinear function for each vertebral level. When we input these nonlinear functions into a model in Promechanica, we can simulate the clinical experiment with the same external moments. Take an example consisting of an L5-S1 two-body system. Define X as the relative angular rotation between the bodies. The unit of angle is degree; the unit of displacement is Nm. We input the conditional spring stiffness coefficients as followings: Phase1: -5.8 < X < -5 (Degree) Initial condition = (Degree) Spring stiffness = 5/0.8 (Nm/mm) Phase2: -5 < X < -1.7 (Degree) Initial condition = .7 (Degree) Spring stiffness = 5/3.3 (Nm/mm) Phase3: -1.7 < X < 1.7 (Degree) Initial condition (Degree) Spring stiffness = 0 (Nm/mm) Phase4: 1.7 < X < 6.8 (Degree) Initial condition = 1.7 (Degree) Spring stiffness = 2.5/5.1 (Nm/mm) Phase5: 6.8 < X < 8.5 (Degree) Initial condition = 6.8 (Degree) Spring stiffness = 7.5/1.7 (Nm/mm)

PAGE 87

76 Similar to the two-body system we built in the chapter 4, we construct the twobody S1-L5 system in Promechanica, fixing S1 on the ground, using a 6 DOF joint to connect S1 and L5. We add a spring on the rotational joint axis between the two bodies. By adding the appropriate constraints to the system, we can keep the movement of L5 with admissible limits defined in the database. It should be noted that the nonlinear stiffness profiles may bring very large forces into the system. Also, at the point where we change the stiffness coefficient suddenly, we may also bring large forces or high frequency forces to bear on system. So, we add a damper into the system, (c = 10) to make the system more stable. The reconstruction result is depicted in the Figure 4-24. Other stiffness database input code can be found at appendix. -10 -8 -6 -4 -2 0 2 4 6 8 10 -10 -8 -6 -4 -2 0 2 4 6 8 10 database simulation MOMENT Angle Extension/Flexion Fig 4-24 Load-Displacement Compare between Database and Reconstruction

PAGE 88

77 CHAPTER 5 CONCLUSIONS AND RECOMMENDATIONS This thesis presented a method for simulation of the human spine and described the creation of a database for human spine simulation. By using fundamental multibody dynamics theory and computer simulation software, a forward numerical simulation and an inverse dynamics identification experiment are described. The identification method for human spinal motion is based on a relaxation of the original governing equations derived using Lagranges method. The results of the numerical experiment show that the position of the pin joint can be identified as a special case of the identified probability measure. The reconstructed simulation results match the original simulation successfully. Using computer simulation program Promechanica, we also constructed a structural and geometric human spine database based on some published human spine data. The stiffness characteristics of intervertebral bodies were incorporated into the database. Experimental results from the literature were reconstructed from the database by using Matlab. The database should provide the foundation for future scoliosis model simulation and analysis.

PAGE 89

78 APPENDIX A GEOMETRY OF MODELING NODES Table A1 Lumbar 5 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 0 0 2 0 16 11 0 14.3 13.1 3 0 32 0 0 31.7 4.5 4 0 16 -11 0 17.4 -8.7 5 46 43 3 46 42.2 9 6 -46 43 3 -46 42.2 9 7 3 51 0 3 50.5 7.1 8 -3 51 0 -3 50.5 7.1 9 3 52 -11 3 53 -3.6 10 -3 52 -11 -3 53 -3.6 11 21 47 9 21 45.3 15.5 12 -21 47 9 -21 45.3 15.5 13 21 47 -12 21 48.2 -5.3 14 -21 47 -12 -21 48.2 -5.3 15 0 61 -5 0 61.1 3.6 16 0 72 -15 0 73.4 -4.8 17 0 61 -20 0 63.2 -11.3

PAGE 90

79 Table A2 Lumbar 4 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 -3 38.5 2 0 20 14 0 17.5 52.1 3 0 40 0 0 37 37.4 4 0 20 -14 0 16.5 24.1 5 45 43 8 45 52.2 44.9 6 -45 43 8 -45 52.2 44.9 7 3 55 0 3 55 36.8 8 -3 55 0 -3 55 36.8 9 3 58 -14 3 58.5 22.6 10 -3 58 -14 -3 58.5 22.6 11 21 61 14 21 58.5 50.7 12 -21 61 14 -21 58.5 50.7 13 21 61 -16 18 45.3 15.5 14 -21 61 -16 -18 45.3 15.5 15 0 73 -1 0 69.9 35.3 16 0 92 -17 0 88.3 18.6 17 0 73 -21 0 69.2 15.3

PAGE 91

80 Table A3 Lumbar 3 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 0.3 77.7 2 0 19 13 0 21.4 84.8 3 0 37 0 0 36.6 70.6 4 0 19 -13 0 16.5 61.3 5 52 62 5 52 62.1 70.8 6 -52 62 5 -52 62.1 70.8 7 3 52 0 3 51.3 67.8 8 -3 52 0 -3 51.3 67.8 9 3 55 -10 3 52.4 57.4 10 -3 55 -10 -3 52.4 57.4 11 16 58 13 16 59.7 79.4 12 -16 58 13 -16 59.7 79.4 13 16 58 -23 18 58.4 50.7 14 -16 58 -23 -18 58.4 50.7 15 0 76 0 0 74.3 63.2 16 0 97 -9 0 93.3 50.4 17 0 76 -18 0 71.5 45.5

PAGE 92

81 Table A4 Lumbar 2 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 9.7 113.7 2 0 18 13 0 31.2 119.4 3 0 36 0 0 43.3 100.8 4 0 18 -13 0 21.9 95.1 5 48 60 7 48 68.2 98.8 6 -48 60 7 -48 68.2 98.8 7 3 51 0 3 57.3 95.5 8 -3 51 0 -3 57.3 95.5 9 3 53 -9 3 56 86.3 10 -3 53 -9 -3 56 86.3 11 16 58 16 16 69.6 107.9 12 -16 58 16 -16 69.6 107.9 13 15 58 -21 13 59.7 79.4 14 -15 58 -21 -13 59.7 79.4 15 0 71 4 0 77.4 92 16 0 91 -9 0 91.4 72.7 17 0 71 -14 0 71 75.2

PAGE 93

82 Table A5 Lumbar 1 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 22.2 146.4 2 0 18 12 0 43.3 151.2 3 0 36 0 0 55.8 133.5 4 0 18 -12 0 34.7 110.5 5 42 54 6 42 74.8 132.7 6 -42 54 6 -42 74.8 132.7 7 3 50 0 3 68.9 128.5 8 -3 50 0 -3 68.9 128.5 9 3 51 -7 3 67.3 121.6 10 -3 51 -7 -3 67.3 121.6 11 16 56 16 16 80.2 141.3 12 -16 56 16 -16 80.2 141.3 13 14 56 -20 13 69.6 107.9 14 -14 56 -20 -13 69.6 107.9 15 0 69 4 0 88.1 125.4 16 0 89 -3 0 104.2 111.8 17 0 69 -12 0 82.3 110.5

PAGE 94

83 Table A6 Thoracic 12 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 34.7 176.3 2 0 16 12 0 53.6 182.8 3 0 32 0 0 65.1 166.4 4 0 16 -12 0 46.2 160 5 25 55 5 25 88.5 164.1 6 -25 55 5 -25 88.5 164.1 7 3 50 4 3 83.5 164.7 8 -3 50 4 -3 83.5 164.7 9 3 48 -4 3 79.1 157.7 10 -3 48 -4 -3 79.1 157.7 11 15 46 15 15 83.1 176.4 12 -15 46 15 -15 83.1 176.4 13 12 50 -20 13 80.2 141.3 14 -12 50 -20 -13 80.2 141.3 15 0 62 0 0 93.7 157.2 16 0 76 -7 0 104.8 146.2 17 0 62 -15 0 89 142.9

PAGE 95

84 Table A7 Thoracic 11 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 43.6 203.7 2 0 17 10 0 62.8 208.3 3 0 34 0 0 76.1 193.8 4 0 17 -10 0 56.9 189.2 5 30 60 10 30 103.9 195.7 6 -30 60 10 -30 103.9 195.7 7 3 47 8 3 90.9 197.6 8 -3 47 8 -3 90.9 197.6 9 3 48 -1 3 89.2 188.7 10 -3 48 -1 -3 89.2 188.7 11 15 44 14 15 89.8 204.2 12 -15 44 14 -15 89.8 204.2 13 15 48 -13 15 86.1 176.4 14 -15 48 -13 -15 86.1 176.4 15 0 62 -1 0 102.6 184.6 16 0 77 -10 0 114.3 171.6 17 0 62 -14 0 98.8 172.2

PAGE 96

85 Table A8 Thoracic 10 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 51.9 227.4 2 0 16 10 0 70.1 232.3 3 0 32 0 0 82.5 218 4 0 16 -10 0 64.3 213.2 5 30 62 4 30 112.4 213.1 6 -30 62 4 -30 112.4 213.1 7 3 45 5 3 96.4 219 8 -3 45 5 -3 96.4 219 9 3 46 -2 3 95.3 212 10 -3 46 -2 -3 95.3 212 11 13 44 15 13 98.4 228.9 12 -13 44 15 13 98.4 228.9 13 15 46 -10 15 92.8 204.2 14 -15 46 -10 -15 92.8 204.2 15 0 63 -8 0 109.8 201.3 16 0 77 -21 0 119.4 184.8 17 0 63 -21 0 106 188.9

PAGE 97

86 Table A9 Thoracic 9 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 59.2 251 2 0 16 10 0 77.2 256.5 3 0 32 0 0 90.1 242.7 4 0 16 -10 0 72.1 237.2 5 33 66 0 33 122.9 233.9 6 -33 66 0 -33 122.9 233.9 7 3 50 5 3 108.8 242.9 8 -3 50 5 -3 108.8 242.9 9 3 52 -2 3 108.9 235.6 10 -3 52 -2 -3 108.9 235.6 11 12 45 14 12 106.3 252.9 12 -12 45 14 -12 106.3 252.9 13 13 48 -10 13 101.4 228.9 14 -13 48 -10 -13 101.4 228.9 15 0 68 -12 0 121.8 221.8 16 0 83 -24 0 133.1 206.3 17 0 68 -24 0 118.7 210.2

PAGE 98

87 Table A10 Thoracic 8 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 64.3 274.9 2 0 16 10 0 81.5 282.6 3 0 32 0 0 96 270.4 4 0 16 -10 0 78.7 262.8 5 33 64 3 33 128.1 268.9 6 -33 64 3 -33 128.1 268.9 7 3 50 5 3 114.5 272.9 8 -3 50 5 -3 114.5 272.9 9 3 52 -2 3 115.5 265.7 10 -3 52 -2 -3 115.5 265.7 11 12 46 13 12 111.7 281.3 12 -12 46 13 -12 111.7 281.3 13 12 48 -12 12 109.3 252.9 14 -12 48 -12 -12 109.3 252.9 15 0 68 -16 0 129.4 249.6 16 0 82 -35 0 140.6 228.8 17 0 68 -26 0 128 239.7

PAGE 99

88 Table A11 Thoracic 7 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 66.3 298.8 2 0 15 10 0 81.3 308.8 3 0 30 0 0 96.3 298.8 4 0 15 -10 0 81.3 288.8 5 34 58 5 34 124.3 303.8 6 -34 58 5 -34 124.3 303.8 7 3 45 5 3 111.3 303.8 8 -3 45 5 -3 111.3 303.8 9 3 47 -2 3 113.1 296.8 10 -3 47 -2 -3 113.1 296.8 11 11 43 13 11 109.3 311.8 12 -11 43 13 -11 109.3 311.8 13 12 45 -9 12 114.7 281.3 14 -12 45 -9 -12 114.7 281.3 15 0 63 -17 0 129.3 281.8 16 0 76 -32 0 142.3 266.8 17 0 63 -28 0 129.3 270.8

PAGE 100

89 Table A12 Thoracic 6 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 66.3 322 2 0 15 10 0 80.8 322.8 3 0 30 0 0 96.3 323.6 4 0 15 -10 0 81.8 312.8 5 34 58 6 34 123.9 331 6 -34 58 6 -34 123.9 331 7 3 46 6 3 111.9 330.4 8 -3 46 6 -3 111.9 330.4 9 3 47 -2 3 113.3 322.4 10 -3 47 -2 -3 113.3 322.4 11 10 43 13 10 108.6 337.4 12 -10 43 13 -10 108.6 337.4 13 11 45 -8 11 112.3 311.8 14 -11 45 -8 -11 112.3 311.8 15 0 63 -13 0 129.9 312.3 16 0 77 -29 0 144.7 297 17 0 63 -24 0 130.5 301.3

PAGE 101

90 Table A13 Thoracic 5 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 65.1 344.5 2 0 14 10 0 78.4 355.4 3 0 28 0 0 93 346.5 4 0 14 -10 0 79.8 335.5 5 33 54 7 33 118.5 355.3 6 -33 54 7 -33 118.5 355.3 7 3 45 5 3 109.6 352.6 8 -3 45 5 -3 109.6 352.6 9 3 46 -2 3 111.1 345.7 10 -3 46 -2 -3 111.1 345.7 11 10 41 13 10 105.1 360.3 12 -10 41 13 -10 105.1 360.3 13 10 44 -9 10 111.6 337.3 14 -10 44 -9 -10 111.6 337.3 15 0 62 -15 0 128 333.9 16 0 68 -38 0 135.6 311.4 17 0 62 -32 0 129.2 316.9

PAGE 102

91 Table A14 Thoracic 4 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 63.5 366.7 2 0 12 10 0 73.8 378.4 3 0 24 0 0 87.2 370.4 4 0 12 -10 0 76.9 358.7 5 36 48 12 36 109 386.1 6 -36 48 12 -36 109 386.1 7 3 41 6 3 103.1 379 8 -3 41 6 -3 103.1 379 9 3 43 -2 3 106.3 371.4 10 -3 43 -2 -3 106.3 371.4 11 12 36 17 12 96.4 389.1 12 -12 36 17 -12 96.4 389.1 13 10 42 -8 10 108.1 360.3 14 -10 42 -8 -10 108.1 360.3 15 0 61 -6 0 124.7 370.3 16 0 74 -23 0 140.2 355.6 17 0 61 -21 0 127 355.5

PAGE 103

92 Table A15 Thoracic 3 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 59.3 387.9 2 0 11 9 0 67.6 399.4 3 0 22 0 0 80.5 393.6 4 0 11 -9 0 72.3 382.1 5 37 44 14 37 98.2 412.8 6 -37 44 14 -37 98.2 412.8 7 3 38 6 3 94.4 403.5 8 -3 38 6 -3 94.4 403.5 9 3 39 -2 3 97.5 396.1 10 -3 39 -2 -3 97.5 396.1 11 14 31 18 14 84.6 413.3 12 -14 31 18 -14 84.6 413.3 13 12 37 -8 12 99.4 389.1 14 -12 37 -8 -12 99.4 389.1 15 0 56 -6 0 114.9 396.6 16 0 71 -20 0 133.1 387 17 0 56 -19 0 118.3 384.1

PAGE 104

93 Table A16 Thoracic 2 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 53.8 408.3 2 0 10 9 0 60.1 420.2 3 0 20 0 0 72.6 415.1 4 0 10 -9 0 66.3 403.3 5 40 38 17 40 83.7 437.3 6 -40 38 17 -40 83.7 437.3 7 3 33 6 3 82.8 425.2 8 -3 33 6 -3 82.8 425.2 9 3 36 -2 3 88.3 418.7 10 -3 36 -2 -3 88.3 418.7 11 17 27 18 17 73 434.4 12 -17 27 18 -17 73 434.4 13 14 32 -8 14 87.6 413.3 14 -14 32 -8 -14 87.6 413.3 15 0 50 -6 0 102.8 419.8 16 0 63 -19 0 119.5 412 17 0 50 -19 0 107.3 407.5

PAGE 105

94 Table A17 Thoracic 1 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 46.5 428.5 2 0 9 8 0 51.5 439.5 3 0 18 0 0 62.9 435.8 4 0 9 -8 0 58 424.9 5 43 30 14 43 68.2 453.5 6 -43 30 14 -43 68.2 453.5 7 3 34 6 3 75.1 447.8 8 -3 34 6 -3 75.1 447.8 9 3 34 -2 3 78.4 440.5 10 -3 34 -2 -3 78.4 440.5 11 20 26 16 20 63.7 453.7 12 -20 26 16 -20 63.7 453.7 13 17 30 -8 17 76 434.4 14 -17 30 -8 -17 76 434.4 15 0 51 -5 0 95.1 444.7 16 0 66 -17 0 113.7 439.8 17 0 51 -17 0 100 433.7

PAGE 106

95 APPENDIX B ATTACHING POSITIONS OF SCREW & HOOKS Table B1 Screw & Hooks on L5 (mm) Node Local_X Local_Y Local_Z 1 18.9 43 7.2 2 -18.9 43 7.2 3 33.5 43 11 4 -33.5 43 11 5 33.5 43 0 6 -33.5 43 0 7 3 51 0 8 -3 51 0 9 3 52 -11 10 -3 52 11

PAGE 107

96 Table B2 Screw & Hooks on L4 (mm) Node Local_X Local_Y Local_Z 1 13.8 43 7.3 2 -13.8 43 7.3 3 33 43 14 4 -33 43 14 5 33 43 0 6 -33 43 0 7 3 55 0 8 -3 55 0 9 3 58 -14 10 -3 58 -14

PAGE 108

97 Table B3 Screw & Hooks on L3 (mm) Node Local_X Local_Y Local_Z 1 20.7 62 10.6 2 -20.7 62 10.6 3 34 62 13 4 -34 62 13 5 34 62 0 6 -34 62 0 7 3 52 0 8 -3 52 0 9 3 55 -10 10 -3 55 -10

PAGE 109

98 Table B4 Screw & Hooks on L2 (mm) Node Local_X Local_Y Local_Z 1 20.2 60 10.5 2 -20.2 60 10.5 3 32 60 13 4 -32 60 13 5 32 60 0 6 -32 60 0 7 3 51 0 8 -3 51 0 9 3 53 -9 10 -3 53 -9

PAGE 110

99 Table B5 Screw & Hooks on L1 (mm) Node Local_X Local_Y Local_Z 1 17.7 54 8.8 2 -17.7 54 8.8 3 29 54 12 4 -29 54 12 5 29 54 0 6 -29 54 0 7 3 50 0 8 -3 50 0 9 3 51 -7 10 -3 51 -7

PAGE 111

100 Table B6 Screw & Hooks on T12 (mm) Node Local_X Local_Y Local_Z 1 19.1 55 10 2 -19.1 55 10 3 20 55 12 4 -20 55 12 5 20 55 0 6 -20 55 0 7 3 50 4 8 -3 50 4 9 3 48 -4 10 -3 48 -4

PAGE 112

101 Table B7 Screw & Hooks on T11 (mm) Node Local_X Local_Y Local_Z 1 17.1 60 8.6 2 -17.1 60 8.6 3 22.5 60 10 4 -22.5 60 10 5 22.5 60 0 6 -22.5 60 0 7 3 47 8 8 -3 47 8 9 3 48 -1 10 -3 48 -1

PAGE 113

102 Table B8 Screw & Hooks on T10 (mm) Node Local_X Local_Y Local_Z 1 17.3 62 9.4 2 -17.3 62 9.4 3 21.5 62 10 4 -21.5 62 10 5 21.5 62 0 6 -21.5 62 0 7 3 45 5 8 -3 45 5 9 3 46 -2 10 -3 46 -2

PAGE 114

103 Table B9 Screw & Hooks on T9 (mm) Node Local_X Local_Y Local_Z 1 17.9 66 10 2 -17.9 66 10 3 22.5 66 10 4 -22.5 66 10 5 22.5 66 0 6 -22.5 66 0 7 3 50 5 8 -3 50 5 9 3 52 -2 10 -3 52 -2

PAGE 115

104 Table B10 Screw & Hooks on T8 (mm) Node Local_X Local_Y Local_Z 1 16.9 64 9.7 2 -16.9 64 9.7 3 22.5 64 10 4 -22.5 64 10 5 22.5 64 0 6 -22.5 64 0 7 3 50 5 8 -3 50 5 9 3 52 -2 10 -3 52 -2

PAGE 116

105 Table B11 Screw & Hooks on T7 (mm) Node Local_X Local_Y Local_Z 1 15.9 58 9.4 2 -15.9 58 9.4 3 22.5 58 10 4 -22.5 58 10 5 22.5 58 0 6 -22.5 58 0 7 3 45 5 8 -3 45 5 9 3 47 -2 10 -3 47 -2

PAGE 117

106 Table B12 Screw & Hooks on T6 (mm) Node Local_X Local_Y Local_Z 1 16 58 9.4 2 -16 58 9.4 3 22 58 10 4 -22 58 10 5 22 58 0 6 -22 58 0 7 3 46 6 8 -3 46 6 9 3 47 -2 10 -3 47 -2

PAGE 118

107 Table B13 Screw & Hooks on T5 (mm) Node Local_X Local_Y Local_Z 1 16.2 54 9.3 2 -16.2 54 9.3 3 21.5 54 10 4 -21.5 54 10 5 21.5 54 0 6 -21.5 54 0 7 3 45 5 8 -3 45 5 9 3 46 -2 10 -3 46 -2

PAGE 119

108 Table B14 Screw & Hooks on T4 (mm) Node Local_X Local_Y Local_Z 1 16.7 48 9.6 2 -16.7 48 9.6 3 24 48 10 4 -24 48 10 5 24 48 0 6 -24 48 0 7 3 41 6 8 -3 41 6 9 3 43 -2 10 -3 43 -2

PAGE 120

109 Table B15 Screw & Hooks on T3 (mm) Node Local_X Local_Y Local_Z 1 17.4 44 8.6 2 -17.4 44 8.6 3 25.5 44 9 4 -25.5 44 9 5 25.5 44 0 6 -25.5 44 0 7 3 38 6 8 -3 38 6 9 3 39 -2 10 -3 39 -2

PAGE 121

110 Table B16 Screw & Hooks on T2 (mm) Node Local_X Local_Y Local_Z 1 17.9 38 8.1 2 -17.9 38 8.1 3 28.5 38 9 4 -28.5 38 9 5 28.5 38 0 6 -28.5 38 0 7 3 33 6 8 -3 33 6 9 3 36 -2 10 -3 36 -2

PAGE 122

111 Table B17 Screw & Hooks on T1 (mm) Node Local_X Local_Y Local_Z 1 18.2 30 6.3 2 -18.2 30 6.3 3 31.5 30 8 4 -31.5 30 8 5 31.5 30 0 6 -31.5 30 0 7 3 34 6 8 -3 34 6 9 3 34 -2 10 -3 34 -2

PAGE 123

112 LIST OF REFERENCES [Bel73] Belytschko, T. B., Andriacchi, T.P., Schultz, A. B. and Galante, J. O., Analog studies of Forces in The Human Spine: Computational Techniques, J. Biomechanics, 1973, Vol. 6, pp. 361-371. [Bel74] Belytschko, T. B., Kulak, R.F., Schultz, A. B., and Galante, J. O., Finite Element Stress Analysis of An Intervertebral Disc, J. Biomechanics, 1974, Vol. 7, pp. 277-285. [Buc96] Buckley, M.A., and Johnson, G.R., Computer Simulation of the Dynamics of a Human Arm and Orthosis Linkage Mechanism, Proceedings Institution of Mechanical Engineers, 1996, Vol. 211, Port H, pp. 349-357. [Ebe99] Eberhard, P., Spagele, T., and Gollhoter, Investigations for the Dynamical Analysis of Human Motion, Multibody System Dynamics, 1999, Vol. 3, pp.1-20. [Har98] Harvey, S. B. and Hukins, D. W. L., Measurement of Lumbar Spinal FlexionExtension Kinematics from Lateral Radiograph: Simulation of the Effects of Out-ofPlane Movement and Errors in Reference Point Placement, Medical Engineering and Physics, 1998, Vol. 20, pp. 403-409. [Hur94] Hurmuzlu, Y., Basdogam, C. and Carollo, J.J., Presenting Joint Kinematics of Human Location Using Phase Portraits and Poincare Maps, J. Biomechanics, 1994, Vol. 27, No. 12, pp. 1495-1499. [Kle98] Kleinberger, M. K., Computational Modeling of Central Spine Biomechanics Frontiers in Head and Neck Treatment, N. Yoganandan et al. (eds.) 1998, pp. 398-408. [Kul75] Kulak, R.F., Schulz, A. B., Belytschko, T. B., Galante, J. O., Biomechanical Characteristics of Vertebral Motion Segments and Intervertebral Discs, Symposium on the Lumbar Spine, Orthopedic Clinics of North America, 1975, Vol. 6, No. 1, January, pp. 121-133. [Kul76] Kulak, R.F., Belytschko, T. B., Schulz, A.B., and Galante, J. O., Nonlinear Behavior of the Human, Intervertebral Disc Under Axial Load, J. Biomechanics, 1976, Vol. 9, pp. 377-386. [Kuo98] Kuo, A. D., A Least Squares Estimation Approach to Improving the Precision of Inverse Dynamics Computations, Transactions of The ASME, 1998, Vol. 120, pp. 148-159.

PAGE 124

113 [Oxa92] Oxaland, T. R., Lin, R-M, and Panjabi, M. M., Three-dimensional Mechanical Properties of the Thoracolumbar Junction, Journal of Orthopaedic Research 10, 1992, pp. 573-580. [Pan94] Panjabi, M. M., Oxland, T. R., Yamamoto I., and Crisco, J. J., Mechanical Behavior of the Human Lumbar and Lumbosacral Spine as Shown by Three-Dimensional Load-Displacement Curves, Journal of Bone and Joint Surgery, March, 1994, pp. 413424. [Pop93] Pope, M.H. and Novotny, J.E., Spinal Biomechanics, Journal of Biomechanical Engineering, 1993, Vol. 115, pp. 569-574. [Ris97] Risher, D. W., Schutte, L. M., and Runge, C. F., The Use of Inverse Dynamics Solutions in Direct Dynamics Simulations, Journal of Biomechanical Engineering, Nov. 1997, Vol. 119, pp. 417-422. [Ros98] Rosfedt, M., Ekstom, L., Broman, H. and Hansson, T. Axial Stiffness of Human Lumbar Motion Segments, Force Dependence, Journal of Mechanics, 1998, Vol. 31, pp. 503-509. [Sch73] Schulz, A. B., Belytschko, T .B., Andriacchi, T.P., and Galante, J. O., Analog Studies of Forces in The Human Spine: Mechanical Properties and Motion Segment Behavior, J. Biomechanics, 1973, Vol. 6, pp. 373-383. [Wil83] Williams, J.L., and Belytschko, T. B., A Three-Dimensional Model of The Human Cervical Spine for Impact Simulation, J. of Biomechanical Engineering, 1983, Vol. 105, pp. 321-331. [Wol94] Woltering, H. J., Long, K. L., Osterbauer, P. J., and Fuhr, A. W., Instantaneous Helical Axis Estimation From 3-D Video Data in Neck Kinematics for Whiplash Diagnostics, J. Biomechanics, 1994, Vol. 27, No. 12, pp. 1415-1432. [Xu99] Xu, Y. and Hollerbach, J. M., A Robust Ensemble Data Method for Identification of Human Joint Mechanical Properties During Movement, IEEE Transactions on Biomedical Engineering, 1999, Vol. 46, No. 4, pp. 409-419. [Zha98] Zhang Y. M., Vool, M. J., Wang, M. and Johnson, J. R., Intervertebral Measurement of Lumbar Segmenral Motion With a New Measuring Device, Medical Engineering and Physics, 1998, Volume 20, pp. 139-148.

PAGE 125

114 BIOGRAPHICAL SKETCH Jianzhi Liu was born in Qingdao, Shan Dong province, P.R. China, in 1971. In 1995, he graduated from Mechanical Engineering Department of Shan Dong University of Technology, Jinan, Shan Dong province, P.R. China, with a Bachelor of Science degree in mechanical engineering. In 1998, he graduated from Beijing University of Posts & Telecommunications, Beijing, P. R. China, with a Master of Engineering degree. In the fall of 1998, he joined the Dynamics and Control Group of Professor Kurdila at the Department of Aerospace Engineering, Mechanics, and Engineering Science of the University of Florida. He earned his Master of Science degree in engineering mechanics in August 2001.

Material Information

Title: A Two-dimensional human spine simulation and three-dimensional spine model construction
Physical Description: xi, 114
Language: English
Creator: Liu, Jianzhi ( Dissertant )
Kurdila, Andrew J. ( Thesis advisor )
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2001

Subjects

Subjects / Keywords: Aerospace Engineering, Mechanics, and Engineering Science thesis, M.S   ( local )
Dissertations, Academic -- UF -- Aerospace Engineering, Mechanics, and Engineering Science   ( local )

Notes

Abstract: This thesis describes two different, but complementary modeling tasks for the human spine. Methodologies for modeling and simulation of the motion of the human spine vary dramatically in complexity. While providing detailed stress and strain field representations, full finite element modeling of the human spine can be computationally expensive. Alternatively, nonlinear multibody dynamics representations are often used because of their simplicity. These formulations employ rigid models of vertebrae interconnected via conventional mechanical joints. However, it is well documented that inter-vertebral motion can depart significantly from conventional mechanical joint constraints. We present an identification methodology that employs a relaxation technique in which joint mechanical properties are represented via a probability measure. In addition we describe an effort to obtain accurate three-dimensional models of the human spine. These efforts are motivated by considering scoliosis. Scoliosis is defined as abnormal lateral curvature of the spine. It is usually considered as a three-dimensional deformity, because axial rotation will always accompany the lateral curvature. The correction of the deformity is required when the patient risks severe deformity. A three-dimensional spine model is constructed and a computer simulation tool is provided.
Subject: spine -- multibody dynamics -- relaxation -- identification -- biomechanics -- Pro/Mechanica
General Note: Title from title page of source document.
General Note: Document formatted into pages; contains xi, 114 p.; also contains graphics.
General Note: Includes vita.
Thesis: Thesis (M.S.)--University of Florida, 2001.
Bibliography: Includes bibliographical references.
General Note: Text (Electronic thesis) in PDF format.

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
System ID: UFE0000329:00001

Material Information

Title: A Two-dimensional human spine simulation and three-dimensional spine model construction
Physical Description: xi, 114
Language: English
Creator: Liu, Jianzhi ( Dissertant )
Kurdila, Andrew J. ( Thesis advisor )
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2001

Subjects

Subjects / Keywords: Aerospace Engineering, Mechanics, and Engineering Science thesis, M.S   ( local )
Dissertations, Academic -- UF -- Aerospace Engineering, Mechanics, and Engineering Science   ( local )

Notes

Abstract: This thesis describes two different, but complementary modeling tasks for the human spine. Methodologies for modeling and simulation of the motion of the human spine vary dramatically in complexity. While providing detailed stress and strain field representations, full finite element modeling of the human spine can be computationally expensive. Alternatively, nonlinear multibody dynamics representations are often used because of their simplicity. These formulations employ rigid models of vertebrae interconnected via conventional mechanical joints. However, it is well documented that inter-vertebral motion can depart significantly from conventional mechanical joint constraints. We present an identification methodology that employs a relaxation technique in which joint mechanical properties are represented via a probability measure. In addition we describe an effort to obtain accurate three-dimensional models of the human spine. These efforts are motivated by considering scoliosis. Scoliosis is defined as abnormal lateral curvature of the spine. It is usually considered as a three-dimensional deformity, because axial rotation will always accompany the lateral curvature. The correction of the deformity is required when the patient risks severe deformity. A three-dimensional spine model is constructed and a computer simulation tool is provided.
Subject: spine -- multibody dynamics -- relaxation -- identification -- biomechanics -- Pro/Mechanica
General Note: Title from title page of source document.
General Note: Document formatted into pages; contains xi, 114 p.; also contains graphics.
General Note: Includes vita.
Thesis: Thesis (M.S.)--University of Florida, 2001.
Bibliography: Includes bibliographical references.
General Note: Text (Electronic thesis) in PDF format.

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
System ID: UFE0000329:00001

Full Text

A TWO-DIMENSIONAL HUMAN SPINE SIMULATION AND THREE-
DIMENSIONAL SPINE MODEL CONSTRUCTION

By

JIANZHI LIU

A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE

UNIVERSITY OF FLORIDA

2001

ACKNOWLEDGMENTS

Many thanks are in order for Dr. Andrew J. Kurdila for his continual assistance,

expertise, and leadership. I would like to thank Dr. Andrew J. Rapoff and Dr. B. J.

Fregly for their assistance and expertise in helping me on this project. I would also like

to thank all the faculty and students who helped at different times during my research.

A C K N O W L E D G M E N T S ................................................................................................... ii

LIST OF TABLES ............................... ..................................

LIST OF FIGURES ..................................... ................ ............... ... ..... .. vii

A B STR A C T .................................................... ........ ........... ................. .. x

CHAPTERS

1 INTRODUCTION AND LITERATURE REVIEW....................................................1

1.1 In tro du ctio n ...................................... ................................. ............... 1
1.2 Literature Review .......................... ........... ..... 3
1.3 G o al an d O utlin e ............................................ ........................... 6

2 MULTIBODY DYNAMICS THEORY....... ...... ..............................................8

2.1 Two-Dimensional Formulation Example................ ..... ............. 8
2.2 Identification Problem in the Two-Dimensional Example .................................... 15

3 NUM ERICAL EXPERIM ENT ............................................................ ............. .22

3.1 Multibody Dynamics Simulation Tools-Promechanica..................................... 22
3.2 Forward Simulation with Custom Load.......................................................... 24
3.3 Identification of Pin Joint Position Using Inverse Simulation............................ 36

4 SCOLIOSIS MODEL DATABASE CONSTRUCTION ..............................................48

4 .1 G eom etry ..................................................... 4 8
4.2 Stiffness............................. ....... ............................... 58
4.3 A Simple Demonstration of the Stiffness Database...................... .............. 75

5 CONCLUSIONS AND RECOMMENDATIONS ............................... ................77

APPENDICES

A GEOMETRY OF MODELING NODES .....................................................................78

B ATTACHING POSITIONS OF SCREW & HOOKS.................................................95

L IST O F R E FE R E N C E S ...................................................................... ..................... 112

BIOGRAPHICAL SKETCH ........................................................................114

iv

LIST OF TABLES

Table Page

3-1 8 Springs identification........................................................................... ....................34

3-2 16 Springs identification.......................................................................... ...................34

4-1 Scoliosis curve data ................... .................................... .. 50

4-2 Coordinates of origin nodes............................................... .........................................51

A l Lum bar 5 ................................................................78

A 2 Lum bar 4 .......... ...................................................................79

A 3 Lum bar 3 ................................................................80

A 4 Lum bar 2 ........................................................................ 81

A 5 Lum bar 1 ..................................................................................82

A 6 Thoracic 12 ................................................................83

A 7 Thoracic 11 ................................................................84

A 8 Thoracic 10 ................................................................85

A 9 Thoracic 9 ................................................................86

A 10 Thoracic 8 ................................................................87

A l Thoracic 7 ................................................................88

A 12 Thoracic 6 ................................................................89

A 13 Thoracic 5 ................................................................90

A 14 Thoracic 4 ..............................................................................91

V

A 15 T h oracic 3 ................................................................92

A 16 T h oracic 2 ................................................................9 3

A 17 T h oracic 1 ................................................................94

B l Screw & H ooks on L 5 (m m ) ................................................. ............................... 95

B2 Screw & H ooks on L4 (m m ) ................................................. ............................... 96

B 3 Screw & H ooks on L3 (m m ) ................................................. ............................... 97

B4 Screw & H ooks on L2 (m m ) ................................................. ............................... 98

B 5 Screw & H ooks on L (m m ) ................................................. ............................... 99

B 6 Screw & H ooks on T 12 (m m ) ................................................ .............................. 100

B 7 Screw & H ooks on T (m m ) .............................. ... ........................................101

B 8 Screw & H ooks on T 10 (m m ) ................................................ .............................. 102

B 9 Screw & H ooks on T9 (m m ) ................................................. ............................... 103

B 10 Screw & H ooks on T8 (m m ) ................................................ .............................. 104

B 11 Screw & H ooks on T7 (m m ) ........................................ ...................................... 105

B 12 Screw & H ooks on T6 (m m ) ................................................ .............................. 106

B 13 Screw & H ooks on T5 (m m ) ................................................ .............................. 107

B 14 Screw & H ooks on T4 (m m ) ................................................ .............................. 108

B 15 Screw & H ooks on T3 (m m ) ................................................ .............................. 109

B 16 Screw & H ooks on T2 (m m ) ............................................... .................................... 110

B 17 Screw & H ooks on T (m m ) .................................... ............................................

LIST OF FIGURES

Figure Page

1-1 H um an Spine .................................................. 1

2-1 Tw o vertebral body system ......... .................................... .................... ............... 9

2-2 Tw o-body exam ple illustration.......................... ................................... ............... 11

2-3 Identification procedure illustration.................................................... .. ................. 17

2-4 N um erical sim ulation procedure illustration ........................................ .....................21

3-1 Human spine modeling procedure .............. ..................................... 24

3-2 A two vertebral body model built in Promechanica.............. .......................... 25

3-3 Generalized forces in two-dimensional plane.................................. .......... ......... 26

3-4 Y-translational trajectory numerical results from Promechanica.............. ...............28

3-5 Y-translational trajectory numerical results from Matlab model............... .................28

3-6 Z-translational trajectory numerical results from Promechanica...............................29

3-7 Z-translational trajectory numerical results from Matlab model................. ............29

3-8 X-rotational trajectory numerical results from Promechanica................ .. ...............30

3-9 X-rotational trajectory numerical results from Matlab model............... ................ 30

3-10 Y-translational trajectory from Promechanica model ...............................................31

3-11 Y-translational trajectory from Matlab model ....................................... ...........31

3-12 Z-translational trajectory from Promechanica model............................... ...............32

3-13 Z-translational trajectory from M atlab model .................................... ............... 32

3-14 X-rotational trajectory from Promechanica model........................................................33

3-15 X-rotational trajectory from Matlab model ....................................... ............... 33

3-16 Translation m ovem ent sim ulation....................................................... ............... 35

3-17 R otation m ovem ent sim ulation............................................. ................................... 35

3-18 vertebral body connected by pin......................................................... ............... 36

3-19 Probability density distribution of spring stiffness.....................................................38

3-20 Cumulative distribution of spring stiffness......................................... ............... 38

3-21 Probability density distribution of spring stiffness.....................................................39

3-22 Cumulative distribution of spring stiffness......................................... ............... 39

3-23 Probability density distribution of spring stiffness...................................................40

3-24 Cumulative distribution of spring stiffness......................................... ............... 40

3-25 Probability density distribution of spring stiffness.....................................................41

3-26 Cumulative distribution of spring stiffness......................................... ............... 41

3-27 Probability density distribution of spring stiffness...................................................42

3-28 Cumulative distribution of spring stiffness......................................... ............... 42

3-29 Probability density distribution of spring stiffness...................................................43

3-30 Cumulative distribution of spring stiffness......................................... ............... 43

3-31 Probability density distribution of spring stiffness...................................................44

3-32 Cumulative distribution of spring stiffness......................................... ............... 44

3-33 Probability density distribution of spring stiffness...................................................45

3-34 Cumulative distribution of spring stiffness......................................... ............... 45

3-35 Probabilty density distribution of different spring numbers........................................46

3-36 The cumulative distribution of different spring numbers............................................46

3-37 Error vs. different num ber of springs.................................................. ..... .......... 47

4-1 Spine M odel in Prom echanica M otion........................................................ ... ........... 52

4-2 A Vertebral Model with Complex Surface...................................................................53

4-3 A Sim ple Geom etry V ertebral M odel........................................ .......................... 53

4-4 Right Lateral View of the Primary Nodes ........................................... ............... 54

4-5 Top View of the Prim ary N odes ...................................................................... 55

4-6 Bottom V iew of the Prim ary N odes ................................... ........................................... 56

4-7 Sim plified G eom etry Spine M odel ...................................................................... ...... 57

4-8 Illustration Stiffness Property in Three-dimensional Coordinates System....................58

4-9 Flextion-extension between L1 and L2...................................... ......................... 60

4-10 Flextion-extension between L2 and L3................................ ................................. 61

4-11 Flextion-extension betw een L3 and L4...................................... ....................... 62

4-12 Flextion-extension between L4 and L5................................................... ............... 63

4-13 Flextion-extension betw een L5 and S1...................................... ........ ............... 64

4-14 Right and Left Axial Torque between L1 and L2 ....................................................65

4-15 Right and Left Axial Torque between L2 and L3 ............... ............................... 66

4-16 Right and Left Axial Torque between L3 and L4 ....................................................67

4-17 Right and Left Axial Torque between L4 and L5 ............... ............................... 68

4-18 Right and Left Axial Torque between L5 and S ................................................69

4-19 Left and Right Bending Moment between L1 and L2...............................................70

4-20 Left and Right Bending Moment between L2 and L3 ...........................................71

4-21 Left and Right Bending Moment between L3 and L4...............................................72

4-22 Left and Right Bending Moment between L4 and L5 ...........................................73

4-23 Left and Right Bending Moment between L5 and S .............................................74

4-24 Load-Displacement Compare between Database and Reconstruction .......................76

Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science

A TWO-DIMENSIONAL HUMAN SPINE SIMULATION AND THREE-
DIMENSIONAL SPINE MODEL CONSTRUCTION

By

Jianzhi Liu

December 2001

Chairman: Dr. Andrew Kurdila
Major Department: Aerospace Engineering, Mechanics, and Engineering Science

This thesis describes two different, but complementary modeling tasks for the

human spine. Methodologies for modeling and simulation of the motion of the human

spine vary dramatically in complexity. While providing detailed stress and strain field

representations, full finite element modeling of the human spine can be computationally

expensive. Alternatively, nonlinear multibody dynamics representations are often used

because of their simplicity. These formulations employ rigid models of vertebrae

interconnected via conventional mechanical joints. However, it is well documented that

inter-vertebral motion can depart significantly from conventional mechanical joint

constraints. We present an identification methodology that employs a relaxation

technique in which joint mechanical properties are represented via a probability measure.

In addition we describe an effort to obtain accurate three-dimensional models of

the human spine. These efforts are motivated by considering scoliosis. Scoliosis is

defined as abnormal lateral curvature of the spine. It is usually considered as a three-

dimensional deformity, because axial rotation will always accompany the lateral

curvature. The correction of the deformity is required when the patient risks severe

deformity. A three-dimensional spine model is constructed and a computer simulation

tool is provided.

CHAPTER
INTRODUCTION AND LITERATURE REVIEW

Chapter 1 first talks about some fundamental concepts and issues related to the

human spine modeling, then discusses some basic information about scoliosis. An

outline of the following chapter is given last.

1.1 Introduction

Because each spine is unique and reacts in a slightly different way, the dynamics

of the human spine are difficult to quantify. Therefore, in
cervical -
order to study the dynamics of the spine, it is necessary to vertebrae .- cervical
lordosis
develop a mathematical model of the human spine that is Y

general enough to encompass the possible variations in ,
thoracic 6
form and function, but specific enough to react realistically vertebrae
thoracic
kviphosis

Human vertebral column (Figure 1.1) is a complex

structure and mechanism. It is made up of bony vertebrae lumbar
lordosis
interconnected by discs and ligaments. The normal spine

consists of seven cervical, twelve thoracic, five lumbar, sacrum

five sacral and five coccygeal vertebrae. The sacral and coccyx
coccyx

coccygeal vertebrae are usually fused. In the lateral plane,

the normal spine exhibits physiologic cervical lordosis,

Figure 1-1 Human Spine

thoracic kyphosis, and lumbar lordosis. They increase the flexibility and shock-

absorbing capacity of the vertebral column, while maintaining the adequate stiffness and

stability.

Although a vertebra is not actually a rigid body, usually it can be considered to

behave as such, because deformations within the vertebra are insignificant compared to

the deformations of the connective tissue. Two adjacent vertebrae, the ligaments, inter-

vertebral disc and the facets form a functional spinal unit, or a motion segment. It

exhibits biomechanical characteristics similar to those of the entire spine. Each of the

motion segments has six degree-of-freedom (DOF), and the whole spine can produce

three movements: flexion extension, lateral bending, and axial rotation. We simply

summarize them as the combinations of rotational and translational movements.

The motions in human spine are coupled, where translation or rotation along or

about one single axis is consistently associated with translation or rotation along or about

a second axis. Some literature has shown these coupled motions clearly [Oxa92],

[Pan94]. The coupled motion can lead to computation difficulties while modeling a

complex structure like the human spine.

Another fundamental issue related to our research is scoliosis. Scoliosis is

defined as abnormal lateral curvature of the spine. It is usually considered as a three-

dimensional deformity, because axial rotation will always accompany the lateral

curvature. Most scoliotic curves will stabilize if not treated, but a few will progress.

Mild scoliotic curves result in cosmetic deformities of the trunk. These can create

psychological problems for the patients who have them. Severe scoliotic curves may

alter balance and coordination, interfere with the function of internal organs, allow

premature degeneration of the spinal column, and cause deterioration of neurological

function. Thus, the correction of the deformity is required when the patient risks severe

deformity. The purpose of this correction is to stabilize the spine until a solid fusion

occurs in the treated area. Once the spine has been fused, it can once again perform its

function without the aid of the instruments.

Hence, it is necessary to generate a general model that can represent these

problems related to the human spine. Many efforts haven been made in the pass. They

can be categorized into different sub areas according to the researchers' interests.

1.2 Literature Review

Multiple biomechanical models have been employed for the spinal modeling.

Most can be grouped into two classes. These classes are computational models and

experimental models.

The computational models can be used to provide insight into the mechanisms by

which a particular motion occurs or can be used to investigate the complex dynamics of

multi-linked segments during motion. Also, these models are used to undertake analyses

of activities, which would otherwise be impossible to carry out experimentally.

Computational modeling can be classified into two categories: one includes the numerical

models, which usually use finite element methods. The other includes the analytical

models, where the geometry is simplified and requires few boundary conditions. Both

categories of models have their strengths and weaknesses, but the inherent simplicity of

analytical models often leads to an intuitive understanding of human spine mechanisms.

Not surprisingly, the computational cost and accuracy of these methods vary

greatly. Detailed finite element discretisations of intervertebral discs can provide fine-

scale stress and strain field approximations. Such highly resolved fields are obtained, of

course, at significant computational cost [Bel73] [Bel74] [Kul75] [Kul76] [Sch73]. A

recent study reports on the development and validation of a nonlinear viscoelastic FE

model that qualifies the mechanical responses of the L2/L3 motion segment subset to

time varying external loads. The vertebral bodies are represented as rigid bodies, and

intervertebral discs are modeled via nonlinear finite elements.

An alternative choice makes use of nonlinear multibody dynamics formulations

[Buc96] [Hur94]. In these formulations, the vertebral bodies are assumed to be rigid and

interconnected by a collection of conventional mechanical joints. The advantages of this

formulation include the fact that these models are relatively low dimensional. It can

often provide accurate predictions of large amplitude articulated motion. Optimization

techniques are usually used in the nonlinear dynamics formulation, where objective

functions are chosen, and a set of parameters is found that minimizes that function. One

earlier example is the model developed by Schultz [Sch73]. In his study, Schultz

described his two objective functions, which minimized the compression on the motion

segment.

In addition to computational modeling that determines the physics of spinal

motion, a significant number of clinical experimental studies of spinal motion have

appeared. Some of them are relatively old and represent early efforts to validate specific

spine models [Kul75]. Other earlier experimental studies have estimated effective

compliance parameters for motion segments. These are subsequently employed to

provide more accurate lumped simulation models [Wil83].

There are some other experimental studies of spinal motion that have been carried

out independent of numerical simulation, solely to examine and characterize the physics

of spinal motion. Panjabi et al [Oxa91] [Pan94] investigated a three-dimensional human

spine structure that exhibits non-linear viscoelastic mechanical behavior. As advanced

measurement technologies have become more and more common, the experimental

techniques and protocols have likewise matured. A hybrid approach called EMG assisted

optimization (EMGAO) was developed by using both optimization and an EMG assisted

approach. In other studies [Har98] [Ros98], lateral radiographs are employed to

characterize spinal kinematics and frequency domain characterization of the spinal

motion is studied.

Despite the progress made in computational modeling and experimental

methodologies to refine models, significant technical issue remains unresolved. As noted

by M. Klenberger in [Kle98]

"... Biological material data remains one of the limiting factors for biomechanical

modeling. Existing data are comprised mostly of quasi-static force-displacement

measurements with calculations of linear stiffness. Normal variance of biological

tissues typically yields a standard deviation for stiffness values of at least 20

percent. No matter how refined an element mesh might be, the accuracy of the model

can never be better than 20 percent."

Thus, the uncertainty in data derived from biological specimens has lead to

increased interest in identification and inverse methods, because identification and

inverse methods by their nature are developed to accommodate some level of noise

uncertainty and "poorly conditioned" data. Research [Wol94] [Ris97] [Kle98] [Xu99]

provides recent examples in which identification and inverse methods have been

developed in biomechanics.

1.3 Goal and Outline

The objective of our research is to derive an identification method for human

spinal motion that is based on a relaxation of the original governing equations. This

approach is based on a simple observation. While some human joints are well

represented by simple conventional mechanical joints, the relative motion of human

vertebrae is not amenable to such a straightforward representation. In the following

chapters, we present an identification methodology that employs a relaxation technique in

which joint mechanical properties are represented via the probability measure.

Numerical studies are used to study the performance of this approximation framework.

Chapter 2 discusses the fundamentals of the multibody dynamics formulation.

The study of the dynamics of the multi-body system includes a development of the

governing equations of the motion. In the model, each spine motion segment has six

degree of freedom (DOF), and the whole human spine can undergo the translational and

rotational movements. Based on this observation, the governing equations of motion

(EOM) of the human spine are developed by using Lagrange's Equation.

Based on the identification theory and algorithm introduced in chapter 2, the

numerical experiments are implemented in chapter 3. In this chapter, several cases of

identification are studied. A forward numerical simulation and an inverse identification

numerical experiment are discussed. Most of the simulation work has been achieved

using Promechanica, which is a multibody dynamics simulation software.

In chapter 4, some basic tools for scoliosis research issues have been developed.

These tools include a spine (L5-T1) model constructed using the simulation software

Promechanica in chapter 4. To facilitate the future development of spinal modes

incorporating rod, hook and screw assembles, a spinal simulation database was developed

that includes both geometry and attachment parts. The attachment points of pedicle

screw and hooks are also included in this geometry database. In addition, a stiffness

database has been constructed from the published literature. At the end of this chapter, a

numerical simulation experiment is studied to verify the stiffness properties. The detailed

geometry database and stiffness database can be found in the Appendix.

Chapter 5 provides concluding remarks and suggestions for future work.

CHAPTER 2
MULTIBODY DYNAMICS THEORY

In the first part of chapter 2, a quick review of fundamental multibody dynamics

theory is given. The governing equations for the model are derived based on Lagrange's

Equations. In the second part, we discuss the basic procedures of numerical simulation.

2.1 Two-Dimensional Formulation Example

Before talking about several aspects of the theoretical background for the relaxed

formulation used in the multibody identification problem, we derive the equations of

motion. In this section we first show that classical Lagrangian methods of dynamics are

reasonable to model the dynamics of the human spine.

Generally, we can describe the configuration of a rigid body in space by using six

coordinates. Three of them describe the body translation; three other coordinates define

the orientation of the body. In two-dimensional analysis, the configuration of a rigid

body system can be identified by using three coordinates; two of them define the

translation of the points on the body, and one coordinate defines the orientation of the

body with respect to a selected inertial frame of reference.

Consider one rigid body in figure 2-1. From the viewpoint of rigid body

kinematics, the global position of arbitrary point Pi on the rigid body I can be defined as

r' = 1R+y (2-1)

where r' is the global position of point Pi, Rc is the global position vector of the origin Oi

Figure 2-1 Two vertebral body system

of the body I in the fixed reference frame Bi, and / is the position vector of point Piwith

respect to Oi. So the whole body configuration can be completely defined, provide the

vector components of R,, y' in the right side of the equation are known. The vector

r' and R, define in the global coordinates, and / is defined in the body fixed local

coordinates. We need to express the vector y in terms of the global systems. This is

done by defining the orientation transformation of the body fixed coordinates with

respect to the global frame of reference. For example, the model depicted in figure 2-1 is

comprised of two vertebral bodies, denoted as I and J. Each body has an associated body

fixed reference frame Bi and Bj, whose orientation is defined via a change of basis

expressed in terms of an orientation matrix [L(,0)] as follows.

b, I = [L(0,) b*, I
S b ,,2 [L(O j2 (2-2)
b = [L 0i k Le3

In this equation, bk is the kth basis vector defining body i, e~ is a basis of the

inertial frame E, and O, is the angle measures the rotation of the body frame relative to

the inertial frame. The location of the center of mass of body Bi and Bj is given by

R 1e + (2-3)

R = Xe+ Y,2 (2-4)
R c=X +Y = (2-4)

respectively. Similar expression for velocity and acceleration of the arbitrary point on the

rigid body can be obtained. Thus the position, velocity, and acceleration of an arbitrary

point on the rigid body can be written in terms of these coordinates. Equation (2-1)

implies that the general motion of a rigid body is equivalent to a translation of one point,

usually the original point of the body fixed frame, plus a rotation. From the viewpoint of

kinetics, we can describe the kinetic energy of a rigid body as the sum of translational

kinetic energy and the rotational energy. The total kinetic energy of a system is the sum

of its components' kinetic energy. In the absence of any explicit constraint between the

two bodies, the kinetic energy of the system is easily shown to be

1 2 .2 .2 .2 .2 .2
T=-m X +YY +-Jr6 +-m X +Y +-J 6
2 2 2 j J J+ 2 J J

X
m *

m, (2-5)
1 2 2 2 2 2 J 0
= X Y,0 X0 YJ

Yi

where m, is the mass of body I and J, is the inertia of body I about the 3-axis passing

through its center of mass. X, Y and 0 are the translational and rotational velocity of the

body I. They are all defined in the global coordinates system. The same description can

be applied to the other body J.

e2

f \P '--
."Ri

Figure 2-2 Two-body example illustration

If we sought to enforce a conventional revolute joint at the point denoted as "s" in

figure (2-2), and assume it is frictionless, so no work is done by this constraint, the

kinematics constraint could be expressed as

D(s) = P(s)- RP (s)

X, + [1(0, )TY 1(S) X [(O y,,, (s) (2-6)

In this 2D and two-body system, y,,1 and yY, are defined as the translational

position of the arbitrary points Pi and Pj in X direction, and Y,2 and y,,2 as the

translational position of the arbitrary points Pi and Pj in the Y direction. Both of them are

the functions depending on the position of the revolute joint "s".

One method to solve this problem is to use Lagrange multiplier equation. A

classical Lagrange multiplier formulation of the equations would be

d IT IT T V m
-- +--y '2 = Qk k = ...N (2-7)
dt 8 qk aqk I=1 ak

where Qk are the generalized forces, q, are the generalized coordinates. By defining

potential energy

V, = I (s)c(s)s(ds) (2-8)
2

instead of using Lagrange multiplier, we could use the standard Lagrange's equation to

solve an unconstrained problem.

d 8T 8T 8V 8V,
T -- + + =Q Qk k=l1.-N (2-9)
dt qk) c/ c O qk

where yu is the probability measure. In this case, the kinetic energy is the same as that in

the classical Lagrange's multiplier formulation. We just need to calculate the potential

energy part

SY 1 T d (s)(s )j(ds)2
q 2 (2-10)

= (f (T (s)) (s)u(ds) + (s) (S(ds))
2 Oq 8q

Alternatively, this equation can be written in terms of the summation convention as

IVY I a D, (s)2l (s)t(ds)
8q, 2 aq, (
1 Lj0 (s) BO(s)
f (s)p(ds) + (s) -(ds)
2 F8q, O8q

V ( (0, (s)
-=-- (s)u(ds) (2-11)
Oq, 8q,

In this form, it would be very difficult, if not impossible, to solve the above relaxed

formulation. The probability measure / in equation (2-10), (2-11) could be an element

of an infinite dimensional space. To make the problem tractable, we employ a simple

approximation of the probability measure yu.

N
/(.)= ZK1 () (2-12)
l=1

In this equation, s is the dirac measure concentrated at s,. That is

S1 s, EA
s, (A) = (2-13)

From a mathematical standpoint, we approximate p by /pk which is the weighted sum of

dirac functions concentrated at sk,k = 1.. m. From a practical standpoint, this

approximation can be thought as seeking a distribution of spring elements along the curve

between the two rigid bodies. The expression in equation (2-11) can be simplified to

aV 5K (SI)
8q, Z=1 aq,
rT (2-14)

Recall the constraint function is given as

x(s)R_ (s) R_ (s)
+X, y[(O ,{ (s){ X [ ]T (s) (2-15)
Y, Y,,2 (S)j [ ) Yj, 2 (S)

Recall orientation transformation of the basis, the transform matrix [L(O)] can be

expressed in the two dimensional plane as follows,

&b,i cosO, sin6, 0 ^1
Lbo,2 s-in 0 cosO, O 02J (2-16)
b,,3 0 0 1 le3

rel\ cos0, -sin, 0 b,,,
e = sin 0, cos 0H bb (2-17)
e3 0 0 1 b],3

Plug the basis transform equations (2-16) (2-17) into the constraint equation, the

constraint vector D becomes

X, +cosO, -sinO, l ,(s) X ~rX _cosO -sin0 Jy,(s)
D(S)= Lsine, cosO, [Z,2(S)J l Ls in, cos, jy,2(S)

(2-18)

So, we can describe the constraint function in the inertial frame E, and then calculate the

Jacobian in equation (2-14) by definition,
Oqj

a0, (s) ax, OY, a0, aQ x ayj a0
q 02 002 002 02 02 2(2-19)
ax, aY, ao, aOx aY a0o

where the generalized coordinates q are defined as

q = {X, Yo, X, YJ 8, (2-20)

They are the translational and rotational basis in a two-dimensional coordinates system as

we mentioned above. If we plug the constraint function into equation (2-19),

differentiate the vector components separately, we get the Jacobian matrix

S 0 (-sin(60)y1-cos(6)y2) -1 0 (sin(8 )y, +cos(8 )y2 )
qJ LO 1 (cos(O,)y, sin(O,)y,2) 0 -1 (-cos(O,)y,, +sin(8 )Y,2)y
(2-21)

At this point, we have derived both the kinetic energy and the potential energy for

this two-body system using Lagrange's Equations

d cT a T +V
d O- + OQk (2-22)
dt 0q ) qk dqk

The final set of equations for the relaxed formulation for the two vertebral bodies

depicted in figure (2-1) is

X 1 0
m, 0 1
J, 0 (-sin()y cos()y2 ) (cos(O)y -sin(O)y,)
mI 1 -1 0
m Y 0 -1
J, ; (sin(O)y, +cos(O)Y,2) (-cos(O)Y, +sin(O),2))
Q1
Q2
jp1-cos(o) -sin(OfQ f1"l 3 cos(O) -sin(O.) Yi), Q
LJsin(O) cos(O) y, (SI) YJ sin(QO) cos(O]) Y%,2(s)Jj K

(2-23)

2.2 Identification Problem in the Two-Dimensional Example
The system identification issue in our research is mainly to derive a mathematical

model for the human spine based on clinically measured data. By essentially adjusting

model parameters within this given model, we try to make the output of the simulation

system coincide with the experimental data as well as possible. In the above section, we

have derived the governing equations of the two-body system. But we wish to know how

well this equation can describe the two-body example. In other words, we want to find

the a set of the suitable probability measures that can fit our model with the experimental

data.

As we know, the mathematical model describes the relationship between the input

signals and output signals. Considering the example we have, the inputs would be the

different forces and moments applied on the vetebral bodies, while the output would be

the vetebral body's orientation, position and movement. The outputs are affected not

only by the system characterization, but also by the input signals. We denote the inputs

by u(t) and outputs by y(t), respectively. The relationship is depicted in figure 2-3. All

these signals are functions of time, and the value of the input at time t will be denoted as

u(t). Often, in the identification context, only discrete-time points are considered, since

in the clinical experiment, the instruments used in the experiments typically record the

signal just at discrete time instants. The identification of spine model is based on these

experimental input and output data. Our human spine system identification is comprised

of the following steps:

Design the clinical experiment and collect input-output data from the

experimental procedures. This is the first stepping stone in going from an

experiment to the theratical model. However, these experimental data often

need to be polished by removing trends and outliers, and selecting useful

information from the original data. Many research has been done in this field,

as noted in the literature we reviewed in the introduction.

Select and define a spine model mathematically to represent the experiment.

There are a set of candidate system descriptions, within which a model is to be

found.

Compute the best spine model characteristic parameters according to the

input-output clinical experimental data and a given criterion of fit. Then,

examine the obtained model's properties.

In figure 2-3, the output y(t) is the system response to the input signal u(t), while

d(t) is the desired output. In spine modeling, d(t) can be considered as the experiment

data we obtained from the clinical experiment. By comparing the output of the given

model with the desired output, we can adjust the characteristic parameter within the

mathematical model, then minimize the error e(t) between the desired output d(t) and the

output of identification model.

Sd(t)

u(t) ,,-;'-7 y(t)
S olEL >- COMPARE

e(t)

Figure 2-3 Identification procedure illustration

In our simple example, the identification problem is a minimization problem over

the probability measure space. In this two-body vertebral segment system, we need to

find a stiffness distribution that could reconstruct the clinical experimental data. Suppose

we have measurements (y,y) (which stand for the displacement and velocity

measurements in our case). Then, we define the performance function as

J(i)= y- H,(q(u),q()) + j1 -H2(q (), q())2 (2-24)

The functions H1 and H2 estimate the experimental data from the model for a

given Pu. So the the performance function defines the error E as described in figure 2-3.

In this work we try to identify the stiffness coefficients of the system that minimize the

performance equation.

Recall the governing function we have derived in the first session. It's difficult to

get a analytical solution in the form of H1 and H2 that we can substitute it into equation

(2-24) directly. Numerical analysis is widely used today to solve this problem. It can

provide the numerical answer while the analytical answer is difficult to obtain. There are

two major computing procedures employed in solving this problem. One is to find a

numerical solution to the ordinary differential equations arising from solving the

Lagrange's equations. Another is to use a least-squares method to calculate the spring

stiffness coefficients.

We obtained the governing equation of motion in the form of equation (2-23).

This is a set of second-order ordinary differential equations coupled to each other.

However, these second-order equations can be reduced to a system of simutaneous first-

order equations. Let's express one of the second-order differential equation as a

example,

d2X dX
-= f(t, ) (2-25)
dt dt

The initial values can be described as X(to) = X0 and X'(to)= Xo. We can convert this

to a pair of first-order equations by the simple expedient of defining the derivative as a

second function.

19

dX
t= Y, X(t) = X,, (2-26)
dt

dY
= f(tX, YO), Y,(to)= Xo (2-27)
dt

This pair of first-oder equation is equivalent to the original equations. Similar equivalent

first-order equations can be obtained for the other second order differential equations.

So, equation (2-26) can be covered into

1
m,
1
Y,
Y2 Q2

Y6 12 Q6
Y m, YO 0

Y, 1 Y, 0

1Y 1 0
Y12 1 0
1
1
1 (2-28)

where the vector Y, contains six generalized coordinates q, namely

(x, ,,,X,Y,,,0,,X,, Y,)

And, ZK,( Q,t), stands for the Jacobian matrix multipied by the constraint vector
1=1 Sqj

D as follows.

1 0
0 1
(-sin(6),1 cos(OQ)y,) (cos(,)y, sin(,)y,2)
x
-1 0
0 -1
(sin(0,)y,, + cos(8,)y,2) (- cos(,)y,, + sin(O)y,,2)

{i ] cos(qs ) -sin(6O ) (s,) XJ cos(8,) -sin(O ) g (si)
IYJ sin(8O) cos(,) 'Y(s)J sin(O,) cos(O ) r,,2(sl)j

(2-29)

At this point, we can apply numerical methods to solve these pairs of first-order

equations.

Another numerical issue is related to the identification procedure, where we try to

fit the clinical experiment data into the given mathematical model. Obviously, the

problem involved is one of fitting nonlinear data. The overall algorithm describing the

numerical solution procedure is described in figure 2-4.

First, an initial guess of the probability measure vector/ pis made. That is, the

initial value of the spring stiffness vector ko is given. By using the ODE solver in

Matlab, we can obtain the numerical approximation at the output data. We then compare

these data with the clinical experiment data. By checking the performance function, or

error measure to a certain threshold, we decide to stop or go on to the next iteration step.

The last results in this iterative process are estimates of the identified characterization

parameters. That is, we obtain estimates of the spring stiffness distribution between the

two vertebral bodies. Detailed numerical simulation procedures and simulation results

will be explained in the following chapters.

Figure 2-4 Numerical simulation procedure illustration

4::::::::::

CHAPTER 3
NUMERICAL EXPERIMENT

In this chapter, we first give an introduction to common terms for the simulation

tool Promechanica. Then, two numerical experiments are discussed in detail. The first

one is a forward numerical simulation, and the second is an inverse dynamics

identification experiment.

3.1 Multibody Dynamics Simulation Tools-Promechanica

In this chapter, we will describe the development of a computer model of the

human spine using the commercial dynamics program Promechanica. It has a particular

high level language to allow modeling of rigid bodies, joints and constraints. It is

possible using this language to simulate the movement of mechanisms. Additionally it is

possible to measure forces and moments at required points. Since these quantities will be

used in the remainder of the text, we introduce the basic terms within the programming

language.

Terminology:

Body: Representation of a rigid body within the computer model.

Connection: Joint used to connect two or more bodies together (pins, cylinder,

ball, bearing, universal joint, sliders, free, weld, gimbals, planar and six degree of

freedom).

Load: Applies loads to points, joint axis or between two points on different

parts (gravity, fixed loads, moments, static and dynamic friction).

Driver: Used to specify the acceleration, velocity or position of a part or part

joint axis (ramps or constant functions, cosine, table or combinations of these.

Measure: Used to monitor changes in model parameters during the analysis.

Mass primitive: A geometric shape (brick, cone, cylinder, plate, and sphere)

which can be added to a part, changing its mass, center of mass and

moments/products of inertia.

The geometry of bodies are created by connecting points specified in the body

fixed local coordinate system (LCS). Mass and inertial properties can then be assigned to

the bodies either explicitly by the user or by using mass primitives and material defined

within the program. These bodies can be connected to each other and to a ground point

using a variety of joints to form a mechanical system. Drivers and measures can be

assigned to each joint axis if so required. Surface geometry allows a more realistic

visualization of the model. Surface is attached automatically to mass primitives and can

be attached manually to bodies. The mass and inertial properties of bodies are entered

explicitly by the analyst.

Also there are some other capabilities to tailor the kinematic and dynamic aspects

of the mechanism design. One important capability is the custom load, which is created

outside of Promechanica Motion by the analyst. This is convenient for a user to model a

more complex load. To analyze the stiffness characteristics of the human spine and

simulate its movements, we would like to describe the joint characteristics between two

vertebral bodies as a distribution of springs. There are no such joints can be used directly

inside the Promechanica from it lists of build-in entities. The custom load is good way

that we can define our specific model by using analyst-developed loads.

3.2 Forward Simulation with Custom Load

Recall the main procedure in the human spine modeling as described in Fig 3-1.

First, we obtained the data from the clinical experiments. By using an approximation

method, we can identify the mechanical characteristic as parameters in the model we

construct. Then, we input these identified parameters into Promechanica, to simulate the

motion of the experiment.

Figure 3-1 Human spine modeling procedure

If we consider the approximation procedure as the inverse dynamics problem, the

Promechanica simulation provides the direct dynamics problem. "The direct dynamics

problem, being the inverse of inverse dynamics solutions, should provide a way to check

the inverse dynamics problem."[DLC1997] That is, if the approximate solution,

calculated using our mathematical model in chapter 2, is used as input to Promechanica

simulation with the same model, the outputs of the simulation should match the data as

the input to the identification procedure. If instead, the output of simulation and the input

to the inverse dynamics problem are different, then an error exists either in the

fundamental formulation or the implementation of either the inverse dynamics or direct

dynamics solution techniques. A forward numerical experiment was designed for this

consistency check.

Figure 3-2 A two vertebral body model built in Promechanica

Consider the two-body system example we used in driving the governing

equations. We construct the same two-body system example in Promechanica Motion by

using published data together with information on the spine geometry imported from a

commercially available source. Figure 3-2 shows the two vertebral bodies with geometry

attached.

In this forward simulation, we simply compare two mathematical models. One is

the Matlab model based on the governing equations derived in chapter. Another is the

one constructed in Promechanica. The same external forces are applied to identified

model and Promechanica model. Also, we input the same mechanical parameters to

these two models, like the same number of springs, the same attachment position, and the

same stiffness probability measurements. The outputs of both systems are time responses

to the same external forces. From the clinic viewpoint, the Promechanica simulation

output can be considered as the clinical experimental data. By comparing the outputs of

these two models, we can find if the Matlab model matches the Promechanica model.

z
Q1

Q2

X
Figure 3-3 Generalized forces in two-dimensional plane

The external forces were added along the three generalized coordinates in this

two-dimensional example. They are generalized forces Q1 and Q2 along the Z and Y

translational direction, and one generalized moment Q3 about the X-axis, as shown in the

above figure 3-3. The generalized forces could have different expressions. For example,

in a simulation using 8 discrete springs, we use

Q1 = A1 cos(27rft + )+ A2 cos(2rfh + )
Q2 = Bi COS(27rft + 0)+ B2 CO (2fh +0)
Q3 = C1cos(27rf, +O)+C2 COs(27rfh +0) (3-1)

where fl, fh stand for the low frequency part and high frequency part in the external

forces. In a 16-spring rotational movement simulation, we use

27

Q1 = A1x2 + Bx + C1 + Dx1 + E1x-2
Q2 = A2x2 + B2x+ C2 + D2x 1 + E2x2
Q3 = A3x2 + B3 + C3 + D3x 1 + E3 2 (3-2)

With different coefficients Ai, Bi, Ci, Di, Ei, we can give different force inputs. After the

simulations have been run, we compare the numerical results from both models as shown

figures 3-4 through 3-15.

8-sDrina numerical

simulation comDarison

50 100 150 200

250

Figure 3-4 Y-translational trajectory numerical results from Promechanica

50 100 150 200

250

Figure 3-5 Y-translational trajectory numerical results from Matlab model

40

39 5

39

38 5

38

37 5

37
0 50 100 150 200 250

Figure 3-6 Z-translational trajectory numerical results from Promechanica

40

39 5

39

38 5 -

38

37 5

37
0 50 100 150 200 250

Figure 3-7 Z-translational trajectory numerical results from Matlab model

15

05

0

-05

-1 5
0 50 100 150 200 250

Figure 3-8 X-rotational trajectory numerical results from Promechanica

15

05

0

-0 5 -

-1

-1 5
0 50 100 150 200 250

Figure 3-9 X-rotational trajectory numerical results from Matlab model

16-sDrina numerical simulation comDarison

x 10

50 100 150 200

250

Figure 3-10 Y-translational trajectory from Promechanica model

x 10

50 100 150 200

250

Figure 3-11Y-translational trajectory from Matlab model

38 308

38 307 -

38 306

38 305

38 304

38 303

38 302

38 301

38 3

38 299
0 50 100 150 200 250

Figure 3-12 Z-translational trajectory from Promechanica model

38 309

38 308 -

38 307

38 306

38 305

38 304

38 303

38 302

38 301

38 3

38 299
0 50 100 150 200 250

Figure 3-13 Z-translational trajectory from Matlab model

x 10
1 6

14

12

08

06

04

02

0 50 100 150 200 250

Figure 3-14 X-rotational trajectory from Promechanica model

0 016

0 014

0 012

0 01

0 008

0 006 -

0 004 -

0 002

0
0 50 100 150 200 250

Figure 3-15 X-rotational trajectory from Matlab model

In many different simulation cases, the numerical trajectory from the Matlab

model matches the results from the Promechanica Motion simulation, as shown in the

above figures. At last, we used the forward simulation capability to solve an

identification problem. That is, we use Promechanica Motion to get an output trajectory.

This data can be considered as the clinical, experimental data. Then, the experimental

trajectory data is supplied to the identification procedure in Matlab. We can identify the

spring stiffness distribution that best represents the trajectory obtained from

Promechanica. This is the same as the identification procedure we discussed in the

chapter 3. We compare the original input spring stiffness coefficients use to generate the

experimental trajectory and the identified stiffness coefficients in the following tables.

Table 3-1 8 Springs identification

Spring 1 2 3 4 5 6 7 8

Input 0.2 0.2 0.1 0.3 0.15 0.12 0.08 0.22

Identify 0.1962 0.1989 0.1045 0.29 0.152 0.1227 0.0818 0.2239

Table 3-2 16 Springs identification

Spring 1 2 3 4 5 6 7 8

Input 0.3 0.15 0.12 0.08 0.22 0.25 0.16 0.3

Identify 0.2913 0.1413 0.11 0.0743 0.2151 0.2506 0.1608 0.3048

Spring 9 10 11 12 13 14 15 16

Input 0.1 0.07 0.09 0.18 0.2 0.2 0.2 0.1

Identify 0.1052 0.0776 0.0986 0.19 0.2099 0.1915 0.1919 0.1071

When we input the identified spring stiffness into the simualtion software

Promechanica, we can visualize of the verte bral body movent.

Figure 3-16 Translation movement simulation

Figure 3-17 Rotation movement simulation

3.3 Identification of Pin Joint Position Using Inverse Simulation

The last examples considered trajectory matching problems. That is, we sought

the distributions to minimize the difference between observed data and calculated

outputs.

One of the critical properties of a relaxed formulation is that it should be able to

represent the classical results as a special case. A numerical experiment is used to show

this is the case. Suppose we use a pin joint to connect two vertebral bodies. We want to

know if we can use the relaxed formulation to approximate this kind of the classical

mechanical joint. Obviously, this is an identification problem.

b1,2

e,

e2

Figure 3-18 vertebral body connected by pin

Similar to the previous numerical experiment, we use the simulation tool

Promechanica to construct a two-dimensional, two-body system. By importing a two

vertebral segment system, we use a pin joint to connect the vertebrae together as shown

in figure 3-18. The center of the pin is at a certain point between the two bodies in the

sagittal plane. Since our spine model is a three-dimensional geometry model, we use

drivers to apply constraints on three coordinates (one translational and two rotational),

and thus simplify the three dimensional system to a two dimensional system. After the

dynamical simulation, we get the trajectories of the movement for the system connected

via the classical pin joint.

Then we want to find a distribution of springs that can replace the pin joint and

can obtain the same system outputs. The attachment points for each spring are fixed on

the surface of each vertebral body. The interval distance between these attachment points

is equal. We can make the distance finer by adding more locations to connect the

vertebral bodies. By using the identification method as we discuss in the chapter three,

we can identify the stiffness coefficient value of each spring.

38

The 2 spring forces case

Discrete Probability Density Function
2 Springs Case

data
-' data2
data3
[ data4

x10-14

A 10

15 20 25 30 35
Position of Spring Attaching Point marked by

40 B

Figure 3-19 Probability density distribution of spring stiffness

Cumulative Distribution of Spring Stiffness
x10-14 2 Springs Case

A 10 15

20 25 30 35 40 B
Position of Springs Attaching Points

Figure 3-20 Cumulative distribution of spring stiffness

I I I I

x 10
x10

39

The 4 spring forces case

Discrete Probability Density Function
4 Springs Case

11 20 25 '30 35
Position of Spring Attaching Point Marked By '*

Figure 3-21 Probability density distribution of spring stiffness

Cumulative Distribution of Spring Stiffness
S10-3 4 Springs Case

20 25 30 3
Position of Springs Attaching Points

Figure 3-22 Cumulative distribution of spring stiffness

x 10
x10

40

The 8 spring forces case

Discrete Probability Density Function
8 Springs Case

05
0

A 10 1 2 B
Position of Spring Attaching Point Marked By '*

Figure 3-23 Probability density distribution of spring stiffness

Cumulative Distribution of Spring Stiffness
x 103 8 Springs
4 5 -- ----- | ---- ---- | --------------- ----

4

3 35

a 3
Q- -
Co
U)
S 25

2

S15

O
05

0 --1-------------------------------------------------------- ________________
A 10 15 20 25 30 35 40 B
Position of Springs Attaching Points

Figure 3-24 Cumulative distribution of spring stiffness

41

The 16 scoring forces case

Discrete Pobability Density Function
16 Springs Case

x 1 0
x10

t bu

A 10 15 20 25 30 35 40 B
Position of Spring Attaching Point Marked By '*

Figure 3-25 Probability density distribution of spring stiffness

Cumulative Distribution of Spring Stiffness
x 10-3 16 Springs Case

6 -

Co

0
C

E
S1.

A 10 15 20 25 30 35 40 B
Position of Springs Attaching Points

Figure 3-26 Cumulative distribution of spring stiffness
x 3

Figure 3-26 Cumulative distribution of spring stiffness

42

The 32 spring forces case

Discrete Probability Density Function
32 Springs Case

x 10
x10

35

-o
S 3

25

S 2

015

o

053

O
10 15 20 25 30 35 40
Probability Density Distribution of Springs

Figure 3-27 Probability density distribution of spring stiffness
n-

Pr b b lt e st iti u i n o p i g

Cumulative Distribution of Springs Stiffness
32 Springs Case

0 012

0 01
fc

U)

S0 008

0

S0 006

O 0 004

S0 002

0

0

A 10 15 20 25 30 35 40
Position of Springs Attaching Points

Figure 3-28 Cumulative distribution of spring stiffness

. . . . . . ..I I I I I

43

The 64 sDrinu forces case

Discrete Probability Density Function
64 Springs Cases

x 10
x10

A 10 15 20 25 30 35 40 B
Position of Spring Attaching Point

Figure 3-29 Probability density distribution of spring stiffness

Cumulative Distribution of Spring Stiffness
64 Springs Case
0 02

0 0 18 - - -

0 016
o 0 016 .

0 014

001
S0 012

0 008

S0 006

E0 004

0 002

A 10 15 20 25 30 35 40 B
Position of Springs Attaching Points

Figure 3-30 Cumulative distribution of spring stiffness

x 10

44

The 128 spring forces case

Discrete Probability Density Function
128 Spring Cases

A 10 15 20 25 30 35 40 B
Position of Spring Attaching Point

Figure 3-31 Probability density distribution of spring stiffness

Cumulative Distribution of Springs Stiffness
128 Springs Case
0 04

0 035

c
0 03 .
U)

0 025

0 02
VQ

0 0 015 .
o o

FO 001 .

0 005

A 10 15 20 25 30 35 40 B
Position of Springs Attaching Points

Figure 3-32 Cumulative distribution of spring stiffness

The 256 spring forces case

Discrete Probability Density Function
256 Springs Case

0 005

0
- 04

C

0 03

0 02

S0 01

0
A 10

15 20 25 30
Position of Spring Attaching Point

35 40

Figure 3-33 Probability density distribution of spring stiffness

Cumulative Distribution Function of Spring Stiffness
256 Springs Case

0 35

uE

S0 25 .

0
02

- 015
C

> 01

E
0 0 05

15 20 25 30 35 40 B
Position of Springs Attaching Points

Figure 3-34 Cumulative distribution of spring stiffness

46

If we consider the spring stiffness coefficients as the probability function, the

value of stiffness can be normalized between 0 and 1. By plotting all the spring

coefficients probability distribution and cumulative distribution are obtained. We can see

the identified numerical results clearly. With a number of springs connecting the two

bodies, the lager stiffness coefficients are more likely move to the position where the pin

joint exactly is.

Spring Stiffness Distribution
2,4,8,16,32,64,128,256 springs case
.

0 07

006

a 005

S004

0 03

S002

001

8

Figure 3-35 Probabilty

0 35

03 -

S0 25

02

S0 15

01

E
C- 0 05

20 32 44
Position of Springs Attaching Points

density distribution of different spring numbers

Cumulative Distribution of Spring Stiffness Comparation
2,4,8,16,32,64,128,256 springs cases

20 25 30
position of springs attaching points

Figure 3-36 The cumulative distribution of different spring numbers

Instead of the pin joint between the two vertebral boides, we input these different

spring stiffness distribution into a Promechanica Motion model with the same external

forces, and simulated the system to get the resulting trajectories of the system. By

comparing the error between the original system and the identified system, we found that

better correlation is achieved using a finer distribution vertebral segment. The following

figures show the relationship between the approximation error and different number of

the springs.

28 5

28

27 5

27

26 5

26

25 5

25

24 5

24

23 5
0 10 20 30 40 50 60 70

Figure 3-37 Error vs. different number of springs

The identification above shows that using a finer spring distribution yields a

smaller approximation error between identified model and the ideal model. Also, the

position of the pin joint can be identified more clearly. This simulation illustrates an

important fact: the joint relaxation technique includes, as a special case, the classical pin

joint constraint.

CHAPTER 4
SCOLIOSIS MODEL DATABASE CONSTRUCTION

In chapter 4, we summarize the development of a structural and geometric

database that will support the future study of scoliosis. We begin with a discussion of the

spine stiffness database construction. Finally, we show an example of a two-dimensional

vertebral segment whose stiffness characteristics are modeled using the database. The

tables in that constitute the database can be found in the Appendix.

4.1 Geometry

In order to construct a mathematical model of the human spine, we need to

construct a three-dimensional scoliosis model that can be used in Promechanica

geometry. From the literature, we have found many types of human spine models that

have been published. Based on these published models, we have built an electronic

database. A three-dimensional spine model and coordinate system is constructed for use

with the dynamics simulation tool Promechanica. This simulation package is selected

because scoliosis is a three-dimensional deformity. It makes the spine deform in a way

that requires a six DOF model. However, by adding the corresponding constraints and

drivers, we can still simplify it to a two-dimensional model in some specific cases.

The spine can be idealized as a collection of rigid bodies interconnected by some

joints or elements. The vertebrae are modeled as rigid bodies, since their deformation is

negligible compared to the intervertebral tissue. The geometry of an individual vertebral

body can be constructed from the published data by using the Promechanica geometry

attachment function. In the global coordinates system, the sacrum will be considered as

the fixed ground of the whole model. Each vertebral body above it, from L5 to T1, can

be imported individually by using the subassembly function. In each subassembly

system a local coordinate system is defined. Every vertebral body is imported into its

corresponding subassembly system. We can describe the geometry characteristics in this

local subassembly system. Then, we import subassemblies into the global system one by

one, according to the published data defining the shape of scoliosis. The global

coordinates of each primary node in the subassembly system can be computed by using

the transform matrix we derived in chapter 3. The orientation values in the sagittal plane

are derived from the published data. By inputting the mechanical properties, like the

mass, and other connecting joint constraints, a spine model is constructed in

Promechanica. In table 4-1, some geometry data are given that characterize a spine that

is scoliotic.

Table 4-1 Scoliosis curve data
Vertebra Sagittal Vertebra Inferior Disc
Orientation (deg) Height (mm) Height (mm)

L5 -8 22 15.7

L4 -2 28 14

L3 -11 26 12.2

L2 -21 26 11.4

L1 -21 24 10

T12 -18 24 8.4

T11 -17 20 6.8

T10 -17 20 5.1

T9 -15 20 4.7

T8 -8 20 4.5

T7 0 20 4.0

T6 3 20 3.2

T5 4 20 2.6

T4 9 20 2.2

T3 15 18 2.7

T2 20 18 3.1

T1 24 16 4.5

Table 4-2 Coordinates of origin nodes
Vertebra globalX globalY globalZ

L5 0 0 0

L4 0 -31.19 37.46

L3 0 -32.05 77.73

L2 0 -41.45 112.73

L1 0 -51.50 145.50

T12 0 -64.00 175.40

T11 0 -68.02 202.95

T10 0 -76.32 226.65

T9 0 -83.62 250.35

T8 0 -88.72 274.15

T7 0 -90.72 298.05

T6 0 -90.72 321.25

T5 0 -89.52 343.75

T4 0 -87.92 365.95

T3 0 -81.28 387.22

T2 0 -75.78 407.63

T1 0 -66.04 427.90

Table 4-2 shows the three-dimensional coordinates of the original nodes of each

local subassembly systems. The subassembly systems are imported into the global

system by pointing their original nodes and orientations. Using this information the

global coordinates of all the individual points can be calculated. The relative position of

any arbitrary points on the superior and inferior vertebral bodies is then found in this

global system via the transformation matrix between the local coordinate systems and

global coordinate system.

A spine model form T1 to L5 is generated in Promechanica, as shown in Fig 4-1

Fig 4-1 Spine Model in Promechanica Motion

To make a simpler geometry model, we pick up 17 primary nodes, which could

describe the basic geometric outline of vertebra. And the body of vertebra is simplified

to a cylinder. The other primary nodes can be picked up on the pedicle, superior articular

process, transverse process, and inferior articular process. Using spline curves to link

these nodes, we can get a basic shape of a vertebral body.

Fig 4-2 A Vertebral Model with Complex Surface

Fig 4-3 A Simple Geometry Vertebral Model

8
~ 15 ---

16

< ^__^'

-0-------
3

14

------------ 1

Fig 4-4 Right Lateral View of the Primary Nodes

I W-1 %,i-,f-l-, i~h.-O nn-v Civi-I Maim r-rF- FqrTEl

16
7I--,--~

Fig 4-5 Top View of the Primary Nodes

I W-11%,i-...., n........O ,n=,, l.....n. m r)r- FqrEi 3

16

-C1

.4

I
- - -I

Fig 4-6 Bottom View of the Primary Nodes

I W-1 %,i-,f-l-, i~h.-O nn-v flivim.nil Mam r-F-I qrTE

Fig 4-7 Simplified Geometry Spine Model

A simplified geometry spine model is constructed as shown in Fig 4-7. One

objective of this biomechanical simulation is to study the implants that attach to the bone.

We have derived the position of the pedicle screws, transverse process hooks, and lumina

hooks based on these primary nodes. Because the internal fixation devices used in spinal

surgery are not strictly sensitive to their position, we can model the geometry and the

position of these implant devices by some simple computations using our biomedical

experiment database. The geometric database describing the orientation of the screws

and hooks are described in the Appendix.

4.2 Stiffness
Many researchers have studied on the stiffness properties of the connective tissues

between the vertebral bodies. The literature review showed that different stiffness values

have been reported in different clinical experiments. The diversity in results is partly

because of the differing experimental methods and instruments used in these experiments.

Also, the experimental data will be different when testing different specimens. Among

these experiments, Panjabi et al made extensive studies of the nonlinear viscoelastic

mechanical behavior of the human spine. The three-dimensional behavior of the spine is

presented in the form of average load-displacement curves for different vertebral bodies.

After examining these discrete experiment data, we have identified and tabulated a

nonlinear load-displacement function to describe the stiffness characteristics of each

Y

t

TMY

RY

SRZ
Txl C

Fig 4-8 Illustration Stiffness Property in Three-dimensional Coordinates System

vertebral level. By inputting these simplified stiffness curves into the model we

constructed in Promechanica, we can simulate the clinical experiment and reconstruct the

relationship between the load and vertebrae movements using numerical simulation tools.

This simple test illustrates the utility of the database.

Fig 4-8 illustrates the three-dimensional coordinates system that is necessary to

describe the three-dimensional physical properties of the spine. The origin of the

coordinate system is located at the infer posterior corner in the mid-sagittal plane of the

vertebral body. The X-axis is directed to the left and is perpendicular to the sagittal

plane. The Y-axis is directed superiorly, and the Z-axis is directed anteriorly. The broad

arrows show the pure moments and thin, straight arrows show the translations on the

three rotational axes and three translational axes respectively.

In the following figures, we show the load-displacement curves between each

level of the lumbar spine in three-dimensional system. They are from the location

between first and second lumbar vertebrae through location between the fifth lumbar and

first sacral vertebrae. For each location, the load-displacement relationships are

described according to the different type of loads applied on the vertebrae, pure moments

of flexion-extension, lateral axial torque, and bilateral bending, respectively.

60

Load-displacement curves are shown at each of the five locations on the vertebral

lumbar due to the application of flexion and extension moments. The main motion,

defined as the rotation in the direction of the applied moment, is indicted by the heavier

line. The figures also show the coupled motions, defined as all rotations and translations

other than the main rotation.

L1-L2

12

10

8

6

4

V

-- -- - -- - -- -
------ I I I I I I

0
-2 ------- -----
-2 ...........
0 I i- i

-..-- .... .....
............. .............. .............. .............. .............. I .............. .............. .............. .............. .............. i .............. .............

- ----- -
--------------- ------- ------- -------------- ------- ------- ---------------------- ------

1--------------------------------------------------------------------------------- 1------

-6

-8

-10

-19

-12 -10 -8 -6
Extension

-4 -2 0 2 4 6 8 10 12
Flexion MOMENT (Nm)

] !A

Fig 4-9 Flextion-extension between L1 and L2

J*-

61

L2-L3
12

1 0------------------------ --------------- --------

S -.......- ..................................................................................--

0

-4 ............. ..... ..... .............. ----. ... ........... ........... .............. .............. .............. .............. .............
S-2 --------.... ....... ............. ............ ............ ...
O. -------0---e--e ----- ------
-12 I I I- I

-6 -. . I ........ .... .... ... ....-. ... ... ... -- -- -- -- -- -- -- -- I-- -- -- -- -
o A

o -4

-6

-8

-10

-12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12
Extension Flexion MOMENT (Nm)

Fig 4-10 Flextion-extension between L2 and L3

62

L3-L4

............. ........... ......... .............. ....... T .............. r ... ........... ... ......... ...... .............. .. ................ .....

i i i

I --- .....i ..............................I..................................

.- . . ...--- - .. . . ..-.. . . .. - - -- -I- - - -- - -- - -

-- - - - - - - - - - - - - - - - -- - - -

- - - - - - - - - - - -- - - - - - - - - - -

-10 -8 -6 -4 -2 0 2 4 6 8 10 12

Extension Flexion MOMENT (Nm)

Fig 4-11 Flextion-extension between L3 and L4

-6

-8

-10

-12
-1

-

63

L4 L5

14
T. T .r

A
*

-10 -8 -6

-4 -2 0 2 4 6 8 10 12

Extension Flexion MOMENT (Nm)

Fig 4-12 Flextion-extension between L4 and L5

-12L
-12

64

L5-S1

-............. -.............. -.............. -.............. -....... .............. .............. ............ ..............i .............. .............. ............

........... ............. ... ........ ---- ....... ..... .............. .............. .............. .............. .............

- - ---- -& -- - - - - - -
T r 1

------

- - - - - - - - i - - - - - - - - --i- - - -

-10 -8 -6

-4 -2 0 2 4 6 8 10 12

Extension Flexion MOMENT (Nm)

Fig 4-13 Flextion-extension between L5 and S1

-12 L
-12

Load-displacement curves are next shown at each of the five locations on the

lumbar spine due to the application of right and left axial torques. The main motion,

defined as the rotation in the direction of the applied moment, is indicted by the heavier

line. On the following figures, also shown are the coupled motions, defined as all

rotations and translations other than the main rotation.

E
E 1

0)
0)
0L

7a

I-
O
-1
O
2

.A

L1-L2

_ .......... .............. ............ .............. .............. .. --............... --------_

----------------------........-------------------- ---------- ----^^---+------- --------J---- -

I II

-12 -10 -8 -6
Right Torque

-4 -2 0 2 4 6 8 10 12
Left Torque MOMENT (Nm)

Fig 4-14 Right and Left Axial Torque between L1 and L2

L2-L3

........ ............- .............. -.............. .............. ............................ .............. .............. .............. .............. i.......... _

3D ---- -- -- -- -- --"S -- -- -- -- ----- -^
E

o 0

z

-12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12

Right Torque Left Torque MOMENT (Nm)

Fig 4-15 Right and Left Axial Torque between L2 and L3
I--

-2 .--

Fig 4-15 Right and Left Axial Torque between L2 and L3

67

L3-L4

i--------------

............. ............. .............. .............. .... .............. ............................................... .............. ............. ...
1 T F F1 1 1

V4

-10 -8 -6
Right Torque

-4 -2 0 2 4 6 8 10 12
Left Torque MOMENT (Nm)

Fig 4-16 Right and Left Axial Torque between L3 and L4

68

L4-L5

E

-2 -------- --- -------------------

4 ---- ---- ---------------------- -----i-----------------------------
-2 ......................................... ----- ----------- I.............. .............. ..............I .............. -- .............. .............. .
i- i

I II I I
-o

-4

-12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12
Right Torque Left Torque MOMENT (Nm)

Fig 4-17 Right and Left Axial Torque between L4 and L5

69

L5-S1

-------------- -------------- ---------------------------- --- -------- ------------- -------- -
............ ............ ......... ..... ....... ....... .............. .............. .. ............ .............. .. .... .......-

I I I I I I_ I I I I

2 -10 -8 -6 -4 -2 0 2 4 6 8 10 12
Right Torque Left Torque MOMENT (Nm)

Fig 4-18 Right and Left Axial Torque between L5 and S1

70

Load-displacement curves below correspond to each of the five locations on the

lumbar spine due to the application of left and right lateral bending moments. The main

motion, defined as the rotation in the direction of the applied moment, is indicted by the

heavier line. The figures also show the coupled motions, defined as all rotations and

translations other than the main rotation.

8

6

4

E

0

1--
0)
01

-2
I-

-4

-6

L1-L2

............. .............. --------------............... -------------- .............. .............. -------------- .............. -------------- .............. -............

-

-

-
i i i i 0
------------ --------------_------------- ............... I . -- ..................... ---- --- --- --- --- --- I--- .............. I -- I ......... ....
- ------------ -- -- -------------- -- -- -------------- -- - - - -------------- -- -- -------------- -- -- ------------... .. .. .. .. .. .. .. .. .

-10 -8 -6
Left Bending

-4 -2 0 2 4 6
Right Bending

8 10 12
MOMENT (Nm)

Fig 4-19 Left and Right Bending Moment between L1 and L2

71

L2-L3

............. I.............. .............. .............. .....i. ...... .. .............. T--T--.. ............. ...........

-10 -8 -6
Left Bending

-4 -2 0 2 4 6 8 10 12

Right Bending

MOMENT (Nm)

Fig 4-20 Left and Right Bending Moment between L2 and L3

72

L3-L4

6 ------- ---- -------------------- ----------- -----------------

4-

E

---- .---. .. .. .. .. .... ....... .. .... --.-
E

0

-4 .. .. .. i .. ------- - ............
-- -

-2-

I- I
-4-

-6

-8
-12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12
Left Bending Right Bending MOMENT (Nm)

Fig 4-21 Left and Right Bending Moment between L3 and L4

73

L4-L5

6 ------------------------------- -------------------------------

E
I i i i i i i i i .

i

S 2 ---------.......... ....... ---.-------------7-.------------.- --------------------------------------------------.....---
0

I-

-6-
-4 --------- ---------- ------------ -- -

-8
-12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12
Left Bending Right Bending MOMENT (Nm)

Fig 4-22 Left and Right Bending Moment between L4 and L5

74

L5-S1

S. ........... ............. .............- ............- .............. ............. ............. ...... ...... 4 ......... ............. I ..... ............
6 ----- ------ ------ --------------------------------- ------------ ------------- -------------- ---- -----------

o -I i-------- ----------- ---

-12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12.............
2-8

Left Bending Right Bending MOMENT (Nm)

Fig 4-23 Left and Right Bending Moment between L5 and S1

4.3 A Simple Demonstration of the Stiffness Database

After the spring stiffness database for the spine (Lumbar and thoracic) has been

built, we simply approximate the stiffness coefficients as the nonlinear function for each

vertebral level. When we input these nonlinear functions into a model in Promechanica,

we can simulate the clinical experiment with the same external moments.

Take an example consisting of an L5-S1 two-body system. Define X as the

relative angular rotation between the bodies. The unit of angle is degree; the unit of

displacement is Nm. We input the conditional spring stiffness coefficients as followings:

Phasel: -5.8 < X < -5 (Degree)

Initial condition = -5 (Degree)

Spring stiffness = 5/0.8 (Nm/mm)

Phase2: -5 < X < -1.7 (Degree)

Initial condition = -1.7 (Degree)

Spring stiffness = 5/3.3 (Nm/mm)

Phase3: -1.7 < X < 1.7 (Degree)

Initial condition (Degree)

Spring stiffness = 0 (Nm/mm)

Phase4: 1.7 < X < 6.8 (Degree)

Initial condition = 1.7 (Degree)

Spring stiffness = 2.5/5.1 (Nm/mm)

PhaseS: 6.8 < X < 8.5 (Degree)

Initial condition =6.8 (Degree)

Spring stiffness = 7.5/1.7 (Nm/mm)

Similar to the two-body system we built in the chapter 4, we construct the two-

body S1-L5 system in Promechanica, fixing S1 on the ground, using a 6 DOF joint to

connect S1 and L5. We add a spring on the rotational joint axis between the two bodies.

By adding the appropriate constraints to the system, we can keep the movement of L5

with admissible limits defined in the database. It should be noted that the nonlinear

stiffness profiles may bring very large forces into the system. Also, at the point where we

change the stiffness coefficient suddenly, we may also bring large forces or high

frequency forces to bear on system. So, we add a damper into the system, (c = 10) to

make the system more stable. The reconstruction result is depicted in the Figure 4-24.

Other stiffness database input code can be found at appendix.

1 0 T--- --- --- i -- i -- ---- ---- ---- --1- ---

6.--

Ssimulition

Angle ---------------

-10 -8 -6 -4 -2 0 2 4 6 8 10
Extension/Flexion MOMENT

Fig 4-24 Load-Displacement Compare between Database and Reconstruction

CHAPTER 5
CONCLUSIONS AND RECOMMENDATIONS

This thesis presented a method for simulation of the human spine and described

the creation of a database for human spine simulation. By using fundamental multibody

dynamics theory and computer simulation software, a forward numerical simulation and

an inverse dynamics identification experiment are described. The identification method

for human spinal motion is based on a relaxation of the original governing equations

derived using Lagrange's method. The results of the numerical experiment show that the

position of the pin joint can be identified as a special case of the identified probability

measure. The reconstructed simulation results match the original simulation successfully.

Using computer simulation program Promechanica, we also constructed a

structural and geometric human spine database based on some published human spine

data. The stiffness characteristics of intervertebral bodies were incorporated into the

database. Experimental results from the literature were reconstructed from the database

by using Matlab. The database should provide the foundation for future scoliosis model

simulation and analysis.

APPENDIX A
GEOMETRY OF MODELING NODES

Table Al Lumbar 5
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 0 0
2 0 16 11 0 14.3 13.1
3 0 32 0 0 31.7 4.5
4 0 16 -11 0 17.4 -8.7
5 46 43 3 46 42.2 9
6 -46 43 3 -46 42.2 9
7 3 51 0 3 50.5 7.1
8 -3 51 0 -3 50.5 7.1
9 3 52 -11 3 53 -3.6
10 -3 52 -11 -3 53 -3.6
11 21 47 9 21 45.3 15.5
12 -21 47 9 -21 45.3 15.5
13 21 47 -12 21 48.2 -5.3
14 -21 47 -12 -21 48.2 -5.3
15 0 61 -5 0 61.1 3.6
16 0 72 -15 0 73.4 -4.8
17 0 61 -20 0 63.2 -11.3

Table A2 Lumbar 4
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 -3 38.5
2 0 20 14 0 17.5 52.1
3 0 40 0 0 37 37.4
4 0 20 -14 0 16.5 24.1
5 45 43 8 45 52.2 44.9
6 -45 43 8 -45 52.2 44.9
7 3 55 0 3 55 36.8
8 -3 55 0 -3 55 36.8
9 3 58 -14 3 58.5 22.6
10 -3 58 -14 -3 58.5 22.6
11 21 61 14 21 58.5 50.7
12 -21 61 14 -21 58.5 50.7
13 21 61 -16 18 45.3 15.5
14 -21 61 -16 -18 45.3 15.5
15 0 73 -1 0 69.9 35.3
16 0 92 -17 0 88.3 18.6
17 0 73 -21 0 69.2 15.3

Table A3 Lumbar 3
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 0.3 77.7
2 0 19 13 0 21.4 84.8
3 0 37 0 0 36.6 70.6
4 0 19 -13 0 16.5 61.3
5 52 62 5 52 62.1 70.8
6 -52 62 5 -52 62.1 70.8
7 3 52 0 3 51.3 67.8
8 -3 52 0 -3 51.3 67.8
9 3 55 -10 3 52.4 57.4
10 -3 55 -10 -3 52.4 57.4
11 16 58 13 16 59.7 79.4
12 -16 58 13 -16 59.7 79.4
13 16 58 -23 18 58.4 50.7
14 -16 58 -23 -18 58.4 50.7
15 0 76 0 0 74.3 63.2
16 0 97 -9 0 93.3 50.4
17 0 76 -18 0 71.5 45.5

Table A4 Lumbar 2
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 9.7 113.7
2 0 18 13 0 31.2 119.4
3 0 36 0 0 43.3 100.8
4 0 18 -13 0 21.9 95.1
5 48 60 7 48 68.2 98.8
6 -48 60 7 -48 68.2 98.8
7 3 51 0 3 57.3 95.5
8 -3 51 0 -3 57.3 95.5
9 3 53 -9 3 56 86.3
10 -3 53 -9 -3 56 86.3
11 16 58 16 16 69.6 107.9
12 -16 58 16 -16 69.6 107.9
13 15 58 -21 13 59.7 79.4
14 -15 58 -21 -13 59.7 79.4
15 0 71 4 0 77.4 92
16 0 91 -9 0 91.4 72.7
17 0 71 -14 0 71 75.2

Table A5 Lumbar 1
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 22.2 146.4
2 0 18 12 0 43.3 151.2
3 0 36 0 0 55.8 133.5
4 0 18 -12 0 34.7 110.5
5 42 54 6 42 74.8 132.7
6 -42 54 6 -42 74.8 132.7
7 3 50 0 3 68.9 128.5
8 -3 50 0 -3 68.9 128.5
9 3 51 -7 3 67.3 121.6
10 -3 51 -7 -3 67.3 121.6
11 16 56 16 16 80.2 141.3
12 -16 56 16 -16 80.2 141.3
13 14 56 -20 13 69.6 107.9
14 -14 56 -20 -13 69.6 107.9
15 0 69 4 0 88.1 125.4
16 0 89 -3 0 104.2 111.8
17 0 69 -12 0 82.3 110.5

Table A6 Thoracic 12
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 34.7 176.3
2 0 16 12 0 53.6 182.8
3 0 32 0 0 65.1 166.4
4 0 16 -12 0 46.2 160
5 25 55 5 25 88.5 164.1
6 -25 55 5 -25 88.5 164.1
7 3 50 4 3 83.5 164.7
8 -3 50 4 -3 83.5 164.7
9 3 48 -4 3 79.1 157.7
10 -3 48 -4 -3 79.1 157.7
11 15 46 15 15 83.1 176.4
12 -15 46 15 -15 83.1 176.4
13 12 50 -20 13 80.2 141.3
14 -12 50 -20 -13 80.2 141.3
15 0 62 0 0 93.7 157.2
16 0 76 -7 0 104.8 146.2
17 0 62 -15 0 89 142.9

Table A7 Thoracic 11
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 43.6 203.7
2 0 17 10 0 62.8 208.3
3 0 34 0 0 76.1 193.8
4 0 17 -10 0 56.9 189.2
5 30 60 10 30 103.9 195.7
6 -30 60 10 -30 103.9 195.7
7 3 47 8 3 90.9 197.6
8 -3 47 8 -3 90.9 197.6
9 3 48 -1 3 89.2 188.7
10 -3 48 -1 -3 89.2 188.7
11 15 44 14 15 89.8 204.2
12 -15 44 14 -15 89.8 204.2
13 15 48 -13 15 86.1 176.4
14 -15 48 -13 -15 86.1 176.4
15 0 62 -1 0 102.6 184.6
16 0 77 -10 0 114.3 171.6
17 0 62 -14 0 98.8 172.2

Table A8 Thoracic 10
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 51.9 227.4
2 0 16 10 0 70.1 232.3
3 0 32 0 0 82.5 218
4 0 16 -10 0 64.3 213.2
5 30 62 4 30 112.4 213.1
6 -30 62 4 -30 112.4 213.1
7 3 45 5 3 96.4 219
8 -3 45 5 -3 96.4 219
9 3 46 -2 3 95.3 212
10 -3 46 -2 -3 95.3 212
11 13 44 15 13 98.4 228.9
12 -13 44 15 13 98.4 228.9
13 15 46 -10 15 92.8 204.2
14 -15 46 -10 -15 92.8 204.2
15 0 63 -8 0 109.8 201.3
16 0 77 -21 0 119.4 184.8
17 0 63 -21 0 106 188.9

Table A9 Thoracic 9
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 59.2 251
2 0 16 10 0 77.2 256.5
3 0 32 0 0 90.1 242.7
4 0 16 -10 0 72.1 237.2
5 33 66 0 33 122.9 233.9
6 -33 66 0 -33 122.9 233.9
7 3 50 5 3 108.8 242.9
8 -3 50 5 -3 108.8 242.9
9 3 52 -2 3 108.9 235.6
10 -3 52 -2 -3 108.9 235.6
11 12 45 14 12 106.3 252.9
12 -12 45 14 -12 106.3 252.9
13 13 48 -10 13 101.4 228.9
14 -13 48 -10 -13 101.4 228.9
15 0 68 -12 0 121.8 221.8
16 0 83 -24 0 133.1 206.3
17 0 68 -24 0 118.7 210.2

Table A10 Thoracic 8
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 64.3 274.9
2 0 16 10 0 81.5 282.6
3 0 32 0 0 96 270.4
4 0 16 -10 0 78.7 262.8
5 33 64 3 33 128.1 268.9
6 -33 64 3 -33 128.1 268.9
7 3 50 5 3 114.5 272.9
8 -3 50 5 -3 114.5 272.9
9 3 52 -2 3 115.5 265.7
10 -3 52 -2 -3 115.5 265.7
11 12 46 13 12 111.7 281.3
12 -12 46 13 -12 111.7 281.3
13 12 48 -12 12 109.3 252.9
14 -12 48 -12 -12 109.3 252.9
15 0 68 -16 0 129.4 249.6
16 0 82 -35 0 140.6 228.8
17 0 68 -26 0 128 239.7

Table Al Thoracic 7
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 66.3 298.8
2 0 15 10 0 81.3 308.8
3 0 30 0 0 96.3 298.8
4 0 15 -10 0 81.3 288.8
5 34 58 5 34 124.3 303.8
6 -34 58 5 -34 124.3 303.8
7 3 45 5 3 111.3 303.8
8 -3 45 5 -3 111.3 303.8
9 3 47 -2 3 113.1 296.8
10 -3 47 -2 -3 113.1 296.8
11 11 43 13 11 109.3 311.8
12 -11 43 13 -11 109.3 311.8
13 12 45 -9 12 114.7 281.3
14 -12 45 -9 -12 114.7 281.3
15 0 63 -17 0 129.3 281.8
16 0 76 -32 0 142.3 266.8
17 0 63 -28 0 129.3 270.8

Table A12 Thoracic 6
node localX localY localZ globleX globleY globleZ

1 0 0 0 0 66.3 322
2 0 15 10 0 80.8 322.8
3 0 30 0 0 96.3 323.6
4 0 15 -10 0 81.8 312.8
5 34 58 6 34 123.9 331
6 -34 58 6 -34 123.9 331
7 3 46 6 3 111.9 330.4
8 -3 46 6 -3 111.9 330.4
9 3 47 -2 3 113.3 322.4
10 -3 47 -2 -3 113.3 322.4
11 10 43 13 10 108.6 337.4
12 -10 43 13 -10 108.6 337.4
13 11 45 -8 11 112.3 311.8
14 -11 45 -8 -11 112.3 311.8
15 0 63 -13 0 129.9 312.3
16 0 77 -29 0 144.7 297
17 0 63 -24 0 130.5 301.3