- Email: [email protected]
Analysis of large head-neck motions
0021 9290.79/021 I 0221 802 00/O
J Btmnechunics Vol 12. pp. 21 I-222 Q PergamonPress Ltd 1979 Printedm Great Britam
ANALYSIS OF LARGE HEAD-NECK
MOTIONS*
JAN G. RErn:Kt and WERNERGOLDSMITH College
of Engineering,
Berkeley,
University of California, California, U.S.A.
Abstract An eleven-element lumped-parameter two-dimensional model of the human head -neck system was developed to predict the motion of the head and neck under impact loadings and upper torso accelerations that represent “whiplash” conditions. This model incorporales mechanical analogies for the intervertebral disks. muscles, ligaments and articular facets. These simulations are composed ofsprings and dashpots with both linear and nonlinear moduli. The forces developed in these substructures are calculated as functions of the deformation of the system, and employed in the equations of motion. The model was tested by comparing the experimental and analytical results for the head motion determined from force inputs and data reported by other researchers for three different experiments. In general, highly satisfactory agreement was found to exist between predictions and data despite parameter matching problems. Theforces and stresses extant in the individual subsystems were also examined for one of the test runs. Force histories were plotted for elements which were representative of the maximum force levels calculated. A procedure was also developed to approximate the strain in the spinal nerve. Strains greater than the current biological material test ranges were predicted, indicating the need for extended testing of these materials and the acquisition of more accurate element strength parameters. Methods to improve and expand the model are discussed. with particular reference to the extension of the model to three dimensions.
and Roberts (1968), Rybicki and Hopper (l971), Landkof et al. (1976) and Li er al. (1971). A number of related discrete systems subjected primarily to +G, accelerations have been evaluated by Toth (1966). Hopkins (1971), Kaleps et al. (l971), Orne and Liu (1971), Prasad and King (1974), Huston et al. (1976) and Huston et al. (1978). Summaries of investigations involving loads and motions for the cervical region are given in King and Chou (1976), McElhaney CI al. (1976) and Goldsmith (1979). Only Landkof er u[. (1976), Prasad er ul. (1974), Huston er ul. (1976) and Huston et al. (1978) compared their predicted numerical results to experimental data. The current approach differs from previous analyses in the nature of its mathematical considerations which allow an individual evaluation of the loads acting on all elements incorporated. The present technique provides good agreement with experimental results when compared to other predictions and, further, has the advantage of permitting relatively simple extensions to incorporate both a greater number and more complicated representation of the elements. Also. the capability to develop a three-dimensional model using the present techniques can be effectively implemented.
This paper is concerned with the development of a two-dimensional lumped-parameter model of the motion of a head-neck system under impact or impulsive loading and comparison of its predictions with measured results. The basic elements employed in this analysis included a rigid-body representation of the head and the first 10 vertebrae, a viscoelastic model of the muscles, and an elastic scheme for the ligaments. Disks are represented by combinations of linear springs and dampers; a similar model is used for the articular facets together with an appropriate contact mechanism between their respective surfaces. Impact loads to the head as well as sudden accelerations of the torso were utilized as inputs to the system. Previous work in this area has consisted of experimental as well as analytical and/or numerical work. In the first category, measurements were made on inanimate models (Landkof et uI., 1976; Mertz and Patrick, 1971),animals(Hodgson, 1970: Wickstromer u/., 1965), cadavers (Prasad er al.. 1974) and human volunteers (Ewing and Thomas, 1972; Muzzy and Lustick, 1976). Measurements have included kinematic variables as well as forcing inputs including acceleration, impact and impulsive loadings. Models have consisted necessarily of relatively crude continuum representations with closed-form solutions obtained by Liu and Murray (1966) and Shirazi (1971). and finite difference solutions by Terry
* Received f
15 June
SYSTEM
The present model employs vertebra T3 as its reference for the head-neck system. This vertebra and the torso at and below the T3 level are considered to be a rigid base and the combination is denoted as mass segment 1. The remaining IO rigid mass segments are vertebrae Tl, T2, Cl-C7 and the head. The mass and inertia of these segments were derived from the data of
1978.
Current address: P.O. Box IlZA, RD 1, Mohrsville, PA
19541, U.S.A. IS\,
I2
i
/,
211
212
JAN G. REBER and WERNERGOLDSMITH
HEAD
/
DISK STIFFNESS STRUCTURE DISK DAMPING STRUCTURE
MUSCLE ELEMENTS (TYP),
AL 'URE
/
I
I
TO&O
Fig. 1. Vertebral composite with typical soft tissue element models.
Prasad and King (1974) and included the total mass of the neck segment at the appropriate level. The soft tissues which connect these mass segments, consisting of the disks, ligaments, muscles and the synovial joints between the vertebral articular facets were modeled by combinations of springs and dashpots attached to the vertebral bodies at the appropriate points. Figure 1 shows the vertebral bodies employed, together with one representation of each type of soft tissue attachment.
Each type of element is represented by the same model configuration consisting of a combination of springs and dashpots. The individual elements within one element type differ from each other only in their geometric and resistive parameters. Table 1 illustrates the various element type configurations together with their constitutive equations. Figure 2 presents an anatomical section of the skull and vertebrae Cl-C3 together with the equivalent model showing typical ligament, disk and muscle replacements and their attachments to the hard tissue members. A total of 76 muscular elements were utilized; attachment points for the soft tissues onto the head and vertebrae C4 and T3, considered typical, are depicted in Fig. 3. Table 2 lists the number and type of attachment points on each mass segment. The masses and inertial moments for each mass segment are given in Table 3. The stiffnesses of the linear springs employed, which represent the synovial joint and the disk compliances, ranged from 4.5 to 445.0KN/m. The constants of the damping elements that modeled the viscous effects of the disks and the synovial joints ranged from 0.015 to 3 N/m; the ligament areas varied from 5 to 6400 mm’ and the muscle areas ranged from 10 to 640 mm’. A complete listing of the parameters used for the constitutive equations is given in Reber (1978). The above data were obtained from a variety of sources, with particular emphasis on values determined from a physical replica of the neck currently under construction at the University of California, Berkeley. Specifically, this physical model provided the dimensions of the vertebrae, their centers of gravity, and the distances from the center of gravity to the various soft tissue attachment points.
Table 1. Substructure representation with constitutive equations for axial motion Structure Type and Element Description
Modeling Representation
F,d Ligament Nonlinear
NUSCle
Kelvin-Voigt Solid with Nonlinear Spring
t
Constitutive Equations (F in Newtons. Area in mn') F = Area[3x1051n(3.45~+l) + 2.77~lO~(e'.~~'-l)]
F = Area[4477e(e3"'
2
-l)]+VC
F = (Kl+K2+K3)d+(K4+K5)' Disk Stiffness Structure with Linear Elanents
Disk Damping Structure with Linear Elements
Articular Facet Structure with Linear Elements
.~(2L,)2(L2+d)2]'+L;l~Isi"(e)
F = ClV+(C2+C3)* L2+d)sin(e)V [(2Ll)2+(L2 +d214
F = (K, +K2)d +(Cl+C2)V(L4+d)sins [L32+
(L,+d121’
[(2L,)*
Large head-neck
Fig. 2. Anatomical
and equivalent
model of lower head and upper cervical region showing replacement.
1%
--_--
i:s”i
213
motions
HEAD __-
---
The ligaments of the physical model were simulated by strips of natural rubber. A general nonlinear stress-strain equation was developed for this material by experimentation and data curve fitting. The crosssectional areas of the natural rubber strips actually employed in the physical model were then measured and tabulated. The parameters of the disk model were approximated by means ofdata tabulated by Prasad and King (1974), who presented overall linear axial and rotational stiffness factors for all intervertebral disks. This reference also provided the total articular facet stiffness parameters used for the synovial joint model. The nonlinear stress-strain relation employed for the spring element in the Kelvin-Voigt model used for the musculature was derived by curve-fitting a stress-elongation relation reported by Yamada (1970), who tested human abdominal and neck muscles. The dissipative constants used for the dash-
Table 2. Soft tissue attachments per mass segment:
+Yz
8 TORSO
+ M2T
Fig. 3. Attachment points for the head. C4. T3 and torso. I I Multiae attachments at one kndpoint A Articular facet stiffness D Disk stiffness L Ligament element M Muscle element R Pure damping element.
number
and type
+M31
T3
soft tissue element
~____._ Mass no.
Pure stiffness
1 2 3 4 5 6 7 8 9 IO 11
12 12 12 12 12 12 12 I4 I4 6
6
Pure damping 3 6 6 6 6 6 6 6 5 6 4
Muscle element 37 9 9 9 10 11 IO 10 9 11 21
214
JAN G. REBER and WERNER GOLDSMITH
Table 3. Effective
masses and inertial
moments
Mass no.
Mass (kg)
Moment of inertia kg - m’ x IO’
1 2 3 4 5 6 7 8 9 10 II
1.860 1.860 0.772 0.772 0.772 0.772 0.172 0.772 2.189 2.189 3.454
28.81 20.79 7.457 4.519 3.390 3.390 2.825 1.808 1.469 1.469 318.2
pot elements in the muscle model were obtained by deriving the stress equation of the model as a function of strain (using the stress-strain relation of the nonlinear spring) and strain rate. The resulting equation was solved for the strain rate, ri, and the damping constant was varied until the solution approximated the results obtained by Ohnishi (1963) for creep in muscle fibers.
requires that the composite structure can be separated into elements in which the properties of interest can be assumed to be homogeneous. (ii) The treatment of the vertebrae as rigid bodies as implied by (i). This requires that the actual strains induced in the other elements of the system be much greater than those produced in the vertebrae. (iii) The dynamics of the actual system can be approximated by a two-dimensional analysis. This requires mid-sagital plane symmetry and specification of the orientation of the body elements parallel to this plane. (iv) The anisotropic physical properties of complicated body substructures can be approximated by combinations of uniaxial elements. This is accomplished by inclining the individual modeling elements of the substructures at different angles so as to allow mutually perpendicular motions of the total substructure to induce different defqrmations in the individual elements.
STRUCTURAL ANALYSIS AND EQUATIONS OF MOTION
Reference frames ASSUMPTIONS
Certain assumptions have been employed to reduce the complexity of the structure. These are: (i) Lumped-parameter-analysis applicability. This
The formulation of this model requires the use of three sets of two-dimensional reference frames: two sets of inertial frames and one moving frame set (Fig. 4). The first inertial frame is global and locates the origins of the second inertial frame set (referred to as
Y
Sm MOVING FRAME
m
\
I
MASS m ,
-
, GLOBAL
ORIGIN
Fig. 4. Reference frame relationships.
-.
715
Large head -neck motions
the local inertial frames) with respect to a common This point is arbitrarily positioned 52.5 mm anteriorly and 7.5 mm inferiorly to the center of gravity of vertebra T3. The local inertial frames (one for each mass segment) are located with their origins at the centers of gravity of the mass segments when the system is initially at its no-load equilibrium position. These frames serve as references for the motion of the individual mass bodies; the inclinations of their axes are identical to that of the axes of the global frame. This implies that only x and J coordinates are necessary to locate each local frame within the global plane. The third set of frames is attached to the mass segments, and translates and rotates with the segments. The geometry of the vertebral bodies and the attachment points of the spring-dashpot systems are defined within these moving frames delineated by two translational and one rotational coordinate relative to the local inertial frames. This arrangement determines the position of any point in the model by the addition of the vectors which locate the frames described above with respect to the global origin.
where
point.
Equations
of motion
The motion of the system is determined by the application of Newton’s second law. To evaluate the forces developed in the individual elements, the total deflection of each spring and the velocity difference between the endpoints of the damping units must be computed. The components of the displacement, dij, between any two arbitrary points,,j and k, on masses i and m, respectively, are given by
vi; = ctj = rij fli sin(Oi + iJij) - ci - r,& 0, sin(8,
+ $,,&) + &
(4)
and “Tj f $ = rmk 0, cos(U,
+ IJJ,~)
+ C, - [rij Bi cos(0, + tiij)
+ Pi].
The velocity, Vc. which is directed along dij (the line of action of the two damping elements), is the projection of Vii on dij and is given by
dij vt=Vij’m,
or
(5)
The corresponding
force, .Fij, for linear damping
is
.P,j = CijV$
(6)
where Cij is a damping constant. The torques generated by these forces about the centers ofgravity of the individual vertebral bodies are expressed as the vector cross product of the forces and the distance to the points of application. for example Ti, = rij x f,, = rijfTj cos(6, + tiij) - rijffj
sin (Oi + I/I,,).
(7)
The equations of motion are developed by summation of the forces caused by all the elements attached to each mass segment or
tr:, = rrnkcos(Q, + tirnk) + u, + p: - [ri,Cos(O, +
$ij)
+
Ui +
pf].
and
(1) d:, = rmt sin (0, + $I& + L’, + -
[rij
Sin(Bi
+
$ij)
+
I‘, +
p’,
pf],
where r and j/ are the polar coordinates of the pointsj and k in their respective moving reference frames for masses i and m. The coordinates u, L’and 0 locate the origins of the moving frames with respect to those of the local inertial frames and p” and pY are the components of the displacement from the global origin to the respective local inertial origins (Fig. 4). The force fij developed in the stiffness element attached to point ij is then expressed in terms of this displacement as fij = Jj(Kij, dij).
to point j on mass i is
4, = 0;.
(8)
where Mi and Ii are the mass and inertial moment of the ith mass segment, a and a are linear and rotational accelerations, and ‘/,, Y and 0 are the two linear and one rotational state velocities, respectively. There are 11 sets of such relations, one for each mass segment (i = 1-l 1) in the system, providing a set of 66 first-order nonlinear differential equations describing the system motion.
(2)
where xj represents a functional relationship with respect to point j on mass i and Kij represents the stiffness parameters of the element. The magnitude of the velocity, Vij, of point k on mass m with respect
I& = ‘N,, Ci = Y ‘,, and
Substructure
delineation
The individual substructures used in the model are listed in Table 1 along with the constitutive equations employed. The ligaments are represented as nonlinear springs and the muscles as Kelvin-Voigt solids. The parameter values of these elements were determined as described previously. The nonlinear anisotropic be-
216
JAN G. REBERand WERNERGOLDSMITH
havior under strain exhibited by the intervertebral disks was modeled by the structures indicated in Table 1. The two crosslinked members allow sufficient freedom in the choice of stiffness and damping functions to permit good matching of the overall response of the unit to the actual disk force-deflection curves. As a first approximation, all the elements in these models were assumed to be linear. Because of geometric considerations, the complete response of the structure is nonlinear and anisotropic, a necessity for proper disk representation. The synovial joints at the interfaces of the articular facets are represented in the model by two springs and two dampers (Table 1). One disadvantage of this model is the occurrence of a phenomenon termed “snap through.” Because the forces are calculated strictly from the deflection and velocities of the individual elements without regard to the actual mass boundaries, there is a possibility that, in the current analytical system, the articulating surfaces of the facets could “snap through” each other. To circumvent this problem, a contact criterion was developed for the articular facet surfaces. Contact between the two line segments representing the two surfaces of the articular facets was considered to occur only if these intersected within the segment lengths (Fig. 5). This criterion can be expressed mathematically by IG
where R indicates a displacement with respect to the local inertial frame referring to mass segment i, subscripts 1 and 2 denote the two line segment endpoints on lines i and m, and subscript “con” designates the point of intersection. If both these inequalities are satisfied, contact has occurred. For this analysis, the coordinates of R,,, must be determined from the equations for the two lines containing the relevant segments, yielding in terms of slope m,
R&,6,= $R:..
- R:,) + R:,.
(10)
‘I (F sin a)’ dy
+
1
2
A2E2 L2 [F(L2
s0
-
y)cosa12dy E212
Ln(F~~sa)2 .
,
mi
These equations are evaluated in the computer program at every time step for each facet junction and the contact criterion is checked. To calculate the force due to contact, the pedicles supporting the facets are assumed to be elasticallydeformable objects, rigidly attached to the body of the vertebra. The pedicles are represented as two round bars as shown in cross-section in Fig. 6. The principle of strain energy is applied to this representation to deduce an equivalent stiffness for the deflection at the facets. The total strain energy, V, in the upper pedicle is evaluated as
I0
“i
mj+-
and
fJ=)
I&-R:,
1
>( >
mi
-K,
and
Rml + R&m,,, R&, - R& + -
R&. =
dx
+t s0 L8 [(F
AI&
sin a)&
- x) - F(cos a)L,]’
dx
+* s 0
~511 (11)
= LOCAL INERTIAL ORIGIN i
“CONTACT”
where A is cross-sectional area, E is Young’s Modulus and F is a force (dimensions L are found in Fig. 6). The value of V for the lower pedicle is given by changing the limits in the first two integrals to read from -L/2 to zero and inverting the sign of the second term in the last integral. From strain energy theory
(12) and after integration and partial differentiation of equation (1 l), the deflection, 6, is found to be equal to the product of the applied force, F, and a constant. The reciprocal of this constant is interpreted as an effective stiffness, K,,. The contact force generated is then evaluated as F = 6K,,,. “NO
CONTACT”
Fig. 5. Contact
criterion.
(13)
The deflection 6 is interpreted as being half the perpendicular distance from the facet endpoint closest
217
Large head-neck motions
neck as a three-parameter viscoelastic solid calculate the this model to and used head-neck-junction displacement that was then compared to the experimental data. Figure 7 shows the head accelerations predicted by the present model for which no corresponding test data is available. Figure 8 shows the predicted motion of the head-neck junction and the results obtained by Landkof (1976). It may be noted that the present model quite faithfully reproduces the observed motion of the junction. Huston et al. (1976) presented data from experiments with seated cadavers subjected to pendulum
artificial
/
E, .I, . A,
/, / 1
Ll LOWER ARTICULAR
/_
J
LI
/”
PROCESS
rl --
/
---1 E, 91, ,A, _i
I
r2
/
L2
UPPERARTICULARPROCESS
Fig. 6. Model approximating pedicle stiffnesses.
to the intersection point to the other intersecting line segment (Fig. 5).
RESULTS AND DISCUSSION t, TIME, ms
A computer program was written incorporating the above equations in Fortran IV that was solved on a CDC 6400 system employing a fourth-order Runge-Kutta scheme for the integration. The inputs from three sets of experimental runs were used in the computational procedure to predict the kinematic histories of these tests. A comparison of the predicted and experimental results allows an assessment of the validity of the model. The first comparison was based on the data of Landkof et al. (1976) who conducted non-destructive radial impact experiments on a thin water-filled spherical lucite shell attached to a viscoelastic articular artificial neck mounted on a rigid base by a suspended thin shell of aluminum. The input force to the system was measured by a sandwiched quartz crystal yielding a contact variation given by F(t), N =
nt
1272 sin .s ( >
o
= 0
t 2 7,
Acceleration components of the head for Run 1.
(14)
where r is the contact duration of 2.5 ms. Among other measurements, the displacement history of the junction between the lucite sphere and the artificial neck was determined. Landkof (1976) approximated the
1, TIME,ms
Fig. 8. Displacement components of Cl for Run 1.
JANG. REBER and WERNER GOLDSMITH
218
head impact and a corresponding open-chain type system analysis. The representation used for the muscles and ligaments is similar to that of the present work. The reference does not indicate the methodology of representing the attachments and uses a different scheme for denoting the disks and articular facets than employed here. Figure 9 portrays the reported pulse input and Fig. 10 displays the experimental displacement, velocity and acceleration of the head in the x-direction as well as that predicted both by Huston ef al. (1976) and by the present analysis. The
-
PRESENT
---
EXPERIMENT,
I
0
MODEL HUSTON
2
3
t , TIME,
Fig. 9. Input
IO a
et 01
4
5
rns
pulse for Run 2.
-
-
PRESENT MODEL - MODEL, HUSTON et 01
-
-
-
EXPERIMENT.HUSTON
agreement of the present model can only be regarded as fair with respect to magnitude when compared with the prediction reported by Huston er al. (1976), although the shapes of the responses are generally similar. This could be a result of parameter inaccuracies in the present model with respect to Huston’s experiment such as that involving the mass and inertia of the head. However, the experimental acceleration curve reported by Huston seems inconsistent with the reported input force. The reported acceleration curve implies that the maximum positive x component of force occurs at 6 ms, after the end of the 5 ms duration impact. A second input force cpnsisting of a linear rise to 8 KN in 5 ms followed by a sudden drop to zero force level was used in the present analysis to attempt to duplicate this effect. The present model again predicts acceleration reversal after the end of impact as seen in Fig. II. However, it can be noted that the x components of velocity and displacement are now in much better agreement with the experimental data due to the large impulse imparted by the ramp input function. The last comparison will be made with the data of Ewing and Thomas (1972) where human volunteers strapped in a sled on a high-speed track were subjected to accelerations at levels comparable to those causing hyperextensions. The data, obtained by means of accelerometers and high-speed framing camera photographs and reported in terms of the motion of the head relative to that at the top of the torso, are considered to be very reliable. The input is given as a linear rise to an acceleration of 7.4 g at 14.2 ms followed by a linear decay to zero at 340ms. Figures 12, 13 and 14 portray experimental accelerations, velocities and displacements of the head from
et d
q-/
,
,
,
,
ar
‘r_t,-..------‘--;: 0 I
IO
5 I
I
t.
15 I
20 I
25 ,
TIME.ms
Fig. 10. Comparison of predicted and measured histories of the horizontal (x) components of displacement, velocity and acceleration
for Run. 2. Original
input.
t. Fig.
IL Results corresponding
TIME,
fns
to Fig. 10 for
revised
input.
219
Large head- neck motions
__
MuGEL
---
FXPERIMEhl
I
0 0 ,
//
100 / t,
Fig. 12. Comparison of predicted and measured acceleration components of the head for Run 3.
3r
-
MODEL
---
EXPERIMENT
1
,
50
TIME,
150
I-zoo
250
rns
Fig. 14. Comparison of predicted and measured displacement components of the head relative to T3 (base) for Run 3.
T!L pq#-3 L
;_;i_i:‘!-:;,,;:+
,
..__,’ 3 =
-2 1 0 I
50
100 (
150 ,
200
250
I, TIME,ms Fig. 13. Comparison of predicted and measured velocity components of the head relative to T3 (base) for Run 3.
Fig. 15. Deformationofthestructurefor
Ewing to Thomas (1972) compared to the predictions of the present analysis. Agreement is considered to be generally good except for high frequency spikes in the curves which might be due to the contact criterion built into the articular facet analysis. Figure 15 shows the predicted maximum excursion of the spinal column for this case where numerous contacts between the articular facets are occurring. A force analysis was executed for the run corresponding to the data of Ewing and Thomas (1972) for three instants during the
motion, including the time of maximum excursion of the spine. The force histories for the two ligament and muscle elements, that exhibit the largest magnitudes for all such elements at these three instances, respectively, have been plotted in their entirety in Figs. 16 and 17. The stress peaks of 10 and 4 MPa correspond to strains of 3.5 and 1.1, respectively. The strain for the first case is believed to be too high by at least a factor of 2 due to an inaccurate representation of the ligament area in the physical model.
Run 3 at
f = 140ms.
220
JAN G.
REBERand WERNERGOLDSMITH
limit beyond which the soft tissues should not be loaded without serious potential for injury. It was found from the present analysis that the supraspinous ligament as well as the three ligaments attaching Cl to the head were loaded beyond this bound by as much as a factor of 10 for many ms, as indicated in Fig. 16. The region between the head and Cl also exceed the stated limit in the case of muscles, as well as the longus colli muscle. Typical computed load histories for the disks are presented in Figs. 18, 19 and 20; the maximum loads are within the test ranges reported in the literature (Kulak er al., 1976). The sharp shear reversal noted in Fig. 18 may very well represent a potential source of damage.
LIGAMENT ELEMENT IO -
5-
0
I
-I
a” “c
6
::
LIGAMENT NO. 109
ELEMENT
50 I
100 1
nw kl? d
3
0
I 0 I
t,
150 I
200 I
IO ---AXIAL
250 I
I-’ ,’
---SHEAR [“’
TIME.ms
Fig. 16. Representative histories involving maximum stresses
,J
5
/ ‘I \,
‘1
\ ‘_. ‘1 ‘J\\
of two ligament elements.
MUSCLE ELEMENl NO. 104
n
Fig. 18. Axial and shear force histories for the intervertebral disk at level T2-T3.
-
2r MUSCLE ELEMENT NO. 124
z 2 &I d
I-
I
0.
-I-
O I
50 I
100 I t,
150 1
200 1
250
TIME,ms
Fig. 19. Axial and shear force histories for the intervertebral disk at level C7-Tl.
Fig. 17. Representative histories involving maximum stresses of two muscle elements. ---
Replications for the ligaments in the physical model currently under construction consist of natural rubber, and this material as well as excised muscles (King and Chou, 1976) were tested in quasi-static tension to an upper limit of about 0.5 MPa, well into the nonlinear stress-strain regions of these materials. Although not definitive from a physiological damage point of view, this magnitude was felt to be a representative upper
AXIAL SHEAR
Fig. 20. Axial and shear force histories for the intervertebral disk at level C4-C5.
221
Large head-neck motions One of the most frequent complaints of injury in the cervical region involves trauma to nerves, particularly to those within the spinal cord. As an estimate of mechanical damage due to excessive elongation, a criterion might be placed on the permissible extension of the spinal cord. In the present model, this extension is hypothesized to be the same as that of the posterior longitudinal ligament which is also contained within the vertebral foramen and is explicitly represented in the present model. An examination of the angular deformation and strain experienced by the spinal cord on the basis of the above hypothesis points to three regions whose loadings indicate that they are most likely damaged before others. The first, between C6 and C7, exhibits the largest positive strain of I.52 together with a large angular dislocation of 52.4” at a time 100 ms beyond the initiation of loading. The other two regions are between C2 and C3, and the head and C 1, respectively, where the strains are of the order of 0.4 coupled with angular deformations of the order of 30”. The present system can be modified with respect to the number of elements and element parameters by adjustment of the input data to the system and may also be extended to three dimensions. The latter would require the derivation of general displacement and velocity equations in three dimensions and the inclusion of the third linear acceleration and the two rotational acceleration components to the equations of motion. It is estimated that the modifications mentioned above would increase computation time by a factor of eight.
CONCLUSIONS
A lumped parameter model of the head and the cervical spine region encompassing the first 9 vertebrae fixed on a rigid base has been developed and solved by a fourth-order Runge-Kutta integration scheme. The elements were modeled by combinations of rigid masses, nonlinear springs and linear dashpots and comprise the head, the vertebral bodies, the disks, ligaments, muscles and articular facets. The program was evaluated for three situations where experimental results had been obtained by other investigators. Vertebral contact at the articular facets was handled by means of the action of a contact force when such a situation prevailed. It was found that good agreement was obtained in all cases when the predictions of the
present model were compared with the test data, even though the individual element parameters could not be matched precisely to each of the various cases studied. A mechanical criterion of a limit of damage has been suggested based on test results of both physiological replacement materials and excised human soft tissue. Stresses on soft tissue and loads on the vertebrae were computed, and the areas of most likely damage under the loading conditions described were delineated. It was found that several tissues exceeded this artificial
limit substantially and for an extended period of time, pointing to the possibility of serious spinal cord injury.
Acknowledgements - The authors would like to extend their appreciation to Michael Kabo and Professor R. L. Taylor for significant technical assistance. The work was performed under NIH Grant No. 5ROl AM18384 supplemented by Contract No. NOCO14-76-C-0686from the Office of Naval Research.
REFERENCES
Ewing, C. L. and Thomas, D. J. (1972) Human head and neck responses to impact acceleration. Joint Army--Navy Monograph, No. 21. Goldsmith, W. (1979) Some aspects of head and neck injury and protection. To appear in Progress in Biomechanics, Proc. NATO Advanced Study Institute. Noordhoff. Hodgson, V. (1970) In Impact Injury and Crash Protection (Edited by Gurdjian, E. S. et al.), pp. 2755307. Thomas, Springfield. Hopkins, G. R. (1971) Nonlinear lumped-parameter mathematical model of dynamic response of the human body. Symposium on Biodynamic Models and Their Applications, pp. 843-849. Aerospace Medical Research Lab, WPAFB,
Ohio (AMRL-TR-7 l-29). Huston, J. C., Passerella, C. E. and Huston, R. L. (X976) Numerical prediction of head/neck response to shock impact. In Measurement and Prediction of Structuraland Biodynamic Crash-Impact Response (Edited by Saczalski, K. J. and Pilkev. W. D.), DV.1377149. ASME. New York. Huston, R. L., Huston, J.-C. and Harlow, M’. W. (1978) Comprehensive, three dimensional head-neck model for impact and high-acceleration studies. Aoiat.. Space Environ. Med. 49, 205-210. Kaleps, I., von Gierke, H. E. and Weis, E. B. (1971) A fivedegree-of-freedom mathematical model of the body. Symposium on Biodynamic Models and Their Applications, pp. 211-231. Aerospace Medical Research Lab, WPAFB, Ohio (AMROL-TR-71-29). King, A. I. and Chou, C. C. (1976) Mathematical modeling, simulation and experimental testing of biomechanical system crash response. J. Biomech. 9, 301~317. Kulak, R. F., Belytschko, T. B. and Schultz, A. B. (1976) Nonlinear behavior of the human intervertebral disc under axial load. J. Biomech. 9, 301-317. Landkof, B., Goldsmith, W. and Sackman, J. L. (19761Impact on a head-neck structure. J. Biomech. 9, 141-151. Li, T. G., Advani, S. H. and Lee, U. C. (1971) The effect of initial curvature on the dynamic response of the spine to axial accelerations. Symposium on Biodynamic Models and Their Applications, pp. 5555569. Aerospace Medical Research Lab. WPAFB, Ohio (AMRL-TR-71-29). Liu, Y. I. and Murray, J. D. (1966) A theoretical study of the effect of impulse on the human torso. In Biomechanics (Edited by Fung. Y. C.). vv. 167-186. ASME. New York. McElhaney; J., Roberts, VI *L. and Hilyard, J. (1976) Biomechanics of Trauma. Duke University Press. Durham. North Carolina. Mertz, H. J. and Patrick, L. M. (1971) Strength and resoonse ofthe human neck. In Proc. 15th Stapp Ca;Crash Co& pp. 2077255. SAE Paner No. 710955. Muzzy, W. H. III and Lustick, L. (1976) Comparison of kinematic parameters between hybrid II head and neck system with human volunteers for -G, acceleration profiles. In Proc. 20th Stapp Car Crash Co@. pp. 45-74. SAE Paper No. 760801. Ohnishi, T. (1963) Rheology Biorheoloyy 1, 83-90.
of glycerinated
muscle
fibers.
JAN G. REBERand WERNER GOLDSMITH
222
stiffness coefficient length mass displacement in local inertial plane torque strain energy velocity acceleration displacement stiffness force slope location of local inertial origin in global plane displacement in moving reference frame time coordinates of mass center of the segments in the local inertial frame damping force damping torque velocity components of mass segment centers in local inertial frame angular acceleration deflection angular displacement in local inertial frame angular velocity in local inertial frame contact duration angular velocity angular displacement in moving reference frame.
Ome, D. and Liu, Y. K. (1971) A mathematical model of spinal response to impact. J. Biomech. 4, 49-72. Prasad, P. and King, A. 1. (1974) An experimentally validated dynamic model of the spine. J. appl. Mech., Trans. A.S.M.E. 96, (E), 5466550. Prasad, P., King, A. I. and Ewing, C. L. (1974) The role of the articular facets during +G, acceleration. J. appl. Mech., Truns. A.S.M.E. 41 (E), 321-326. Reber, J. G. (1978) Lumped parameter model of whiplash. M.S. Thesis, University of California, Berkeley. Rybicki, E. F. and Hopper, A. T. (1971) A dynamic model of the spine using a porous elastic material. Symposium on Biodvnamic Models und Their Applications, pp. 851-875. Aerospace Medical Research Lab, WPAFB. Ohio (AMRL-TR-71-29). Shirazi, M. (1971) Response of the spine in biodynamic environments. Symposium on Biodynamic Models and Their Applications, pp. 843-849. Aerospace Medical Research Lab, WPAFB, Ohio (AMRL-TR-71-29). Terry, C. T. and Roberts, V. L. (1968) A viscoelastic model of the human spine subjected to +G, accelerations. J. Biomech. 1, 161-168.
Toth, R. (1966) Multiple-degree-of-freedom, nonlinear spinal model. In Proc. 19rh Ann. Co@ Engng Med. Biol. 8, 102. Wickstrom, J., Martinez, J., Johnston, D. and Tappen, N. C. (1965) Acceleration-deceleration injuries of the cervical spine in animals. In Proc. 7th Stapp Car Crash Conference (Edited by Severy, D.), pp. 284-301. Thomas, Springfield. Yamada, H. (1970) Strength of Biological Materials (Edited by Evans, F. G.). Williams & Wilkins, Baltimore.
Subscripts i.nt j.k
EO” co NOMENCLATURE
A C E F I
area damping coefficient Young’s modulus force moment of inertia
mass segment number attachment point number contact effective total.
Superscripts x.)’ in x, $’direction
I I SUP
functional relationship time derivative absolute value supremum.
J Btmnechunics Vol 12. pp. 21 I-222 Q PergamonPress Ltd 1979 Printedm Great Britam
ANALYSIS OF LARGE HEAD-NECK
MOTIONS*
JAN G. RErn:Kt and WERNERGOLDSMITH College
of Engineering,
Berkeley,
University of California, California, U.S.A.
Abstract An eleven-element lumped-parameter two-dimensional model of the human head -neck system was developed to predict the motion of the head and neck under impact loadings and upper torso accelerations that represent “whiplash” conditions. This model incorporales mechanical analogies for the intervertebral disks. muscles, ligaments and articular facets. These simulations are composed ofsprings and dashpots with both linear and nonlinear moduli. The forces developed in these substructures are calculated as functions of the deformation of the system, and employed in the equations of motion. The model was tested by comparing the experimental and analytical results for the head motion determined from force inputs and data reported by other researchers for three different experiments. In general, highly satisfactory agreement was found to exist between predictions and data despite parameter matching problems. Theforces and stresses extant in the individual subsystems were also examined for one of the test runs. Force histories were plotted for elements which were representative of the maximum force levels calculated. A procedure was also developed to approximate the strain in the spinal nerve. Strains greater than the current biological material test ranges were predicted, indicating the need for extended testing of these materials and the acquisition of more accurate element strength parameters. Methods to improve and expand the model are discussed. with particular reference to the extension of the model to three dimensions.
and Roberts (1968), Rybicki and Hopper (l971), Landkof et al. (1976) and Li er al. (1971). A number of related discrete systems subjected primarily to +G, accelerations have been evaluated by Toth (1966). Hopkins (1971), Kaleps et al. (l971), Orne and Liu (1971), Prasad and King (1974), Huston et al. (1976) and Huston et al. (1978). Summaries of investigations involving loads and motions for the cervical region are given in King and Chou (1976), McElhaney CI al. (1976) and Goldsmith (1979). Only Landkof er u[. (1976), Prasad er ul. (1974), Huston er ul. (1976) and Huston et al. (1978) compared their predicted numerical results to experimental data. The current approach differs from previous analyses in the nature of its mathematical considerations which allow an individual evaluation of the loads acting on all elements incorporated. The present technique provides good agreement with experimental results when compared to other predictions and, further, has the advantage of permitting relatively simple extensions to incorporate both a greater number and more complicated representation of the elements. Also. the capability to develop a three-dimensional model using the present techniques can be effectively implemented.
This paper is concerned with the development of a two-dimensional lumped-parameter model of the motion of a head-neck system under impact or impulsive loading and comparison of its predictions with measured results. The basic elements employed in this analysis included a rigid-body representation of the head and the first 10 vertebrae, a viscoelastic model of the muscles, and an elastic scheme for the ligaments. Disks are represented by combinations of linear springs and dampers; a similar model is used for the articular facets together with an appropriate contact mechanism between their respective surfaces. Impact loads to the head as well as sudden accelerations of the torso were utilized as inputs to the system. Previous work in this area has consisted of experimental as well as analytical and/or numerical work. In the first category, measurements were made on inanimate models (Landkof et uI., 1976; Mertz and Patrick, 1971),animals(Hodgson, 1970: Wickstromer u/., 1965), cadavers (Prasad er al.. 1974) and human volunteers (Ewing and Thomas, 1972; Muzzy and Lustick, 1976). Measurements have included kinematic variables as well as forcing inputs including acceleration, impact and impulsive loadings. Models have consisted necessarily of relatively crude continuum representations with closed-form solutions obtained by Liu and Murray (1966) and Shirazi (1971). and finite difference solutions by Terry
* Received f
15 June
SYSTEM
The present model employs vertebra T3 as its reference for the head-neck system. This vertebra and the torso at and below the T3 level are considered to be a rigid base and the combination is denoted as mass segment 1. The remaining IO rigid mass segments are vertebrae Tl, T2, Cl-C7 and the head. The mass and inertia of these segments were derived from the data of
1978.
Current address: P.O. Box IlZA, RD 1, Mohrsville, PA
19541, U.S.A. IS\,
I2
i
/,
211
212
JAN G. REBER and WERNERGOLDSMITH
HEAD
/
DISK STIFFNESS STRUCTURE DISK DAMPING STRUCTURE
MUSCLE ELEMENTS (TYP),
AL 'URE
/
I
I
TO&O
Fig. 1. Vertebral composite with typical soft tissue element models.
Prasad and King (1974) and included the total mass of the neck segment at the appropriate level. The soft tissues which connect these mass segments, consisting of the disks, ligaments, muscles and the synovial joints between the vertebral articular facets were modeled by combinations of springs and dashpots attached to the vertebral bodies at the appropriate points. Figure 1 shows the vertebral bodies employed, together with one representation of each type of soft tissue attachment.
Each type of element is represented by the same model configuration consisting of a combination of springs and dashpots. The individual elements within one element type differ from each other only in their geometric and resistive parameters. Table 1 illustrates the various element type configurations together with their constitutive equations. Figure 2 presents an anatomical section of the skull and vertebrae Cl-C3 together with the equivalent model showing typical ligament, disk and muscle replacements and their attachments to the hard tissue members. A total of 76 muscular elements were utilized; attachment points for the soft tissues onto the head and vertebrae C4 and T3, considered typical, are depicted in Fig. 3. Table 2 lists the number and type of attachment points on each mass segment. The masses and inertial moments for each mass segment are given in Table 3. The stiffnesses of the linear springs employed, which represent the synovial joint and the disk compliances, ranged from 4.5 to 445.0KN/m. The constants of the damping elements that modeled the viscous effects of the disks and the synovial joints ranged from 0.015 to 3 N/m; the ligament areas varied from 5 to 6400 mm’ and the muscle areas ranged from 10 to 640 mm’. A complete listing of the parameters used for the constitutive equations is given in Reber (1978). The above data were obtained from a variety of sources, with particular emphasis on values determined from a physical replica of the neck currently under construction at the University of California, Berkeley. Specifically, this physical model provided the dimensions of the vertebrae, their centers of gravity, and the distances from the center of gravity to the various soft tissue attachment points.
Table 1. Substructure representation with constitutive equations for axial motion Structure Type and Element Description
Modeling Representation
F,d Ligament Nonlinear
NUSCle
Kelvin-Voigt Solid with Nonlinear Spring
t
Constitutive Equations (F in Newtons. Area in mn') F = Area[3x1051n(3.45~+l) + 2.77~lO~(e'.~~'-l)]
F = Area[4477e(e3"'
2
-l)]+VC
F = (Kl+K2+K3)d+(K4+K5)' Disk Stiffness Structure with Linear Elanents
Disk Damping Structure with Linear Elements
Articular Facet Structure with Linear Elements
.~(2L,)2(L2+d)2]'+L;l~Isi"(e)
F = ClV+(C2+C3)* L2+d)sin(e)V [(2Ll)2+(L2 +d214
F = (K, +K2)d +(Cl+C2)V(L4+d)sins [L32+
(L,+d121’
[(2L,)*
Large head-neck
Fig. 2. Anatomical
and equivalent
model of lower head and upper cervical region showing replacement.
1%
--_--
i:s”i
213
motions
HEAD __-
---
The ligaments of the physical model were simulated by strips of natural rubber. A general nonlinear stress-strain equation was developed for this material by experimentation and data curve fitting. The crosssectional areas of the natural rubber strips actually employed in the physical model were then measured and tabulated. The parameters of the disk model were approximated by means ofdata tabulated by Prasad and King (1974), who presented overall linear axial and rotational stiffness factors for all intervertebral disks. This reference also provided the total articular facet stiffness parameters used for the synovial joint model. The nonlinear stress-strain relation employed for the spring element in the Kelvin-Voigt model used for the musculature was derived by curve-fitting a stress-elongation relation reported by Yamada (1970), who tested human abdominal and neck muscles. The dissipative constants used for the dash-
Table 2. Soft tissue attachments per mass segment:
+Yz
8 TORSO
+ M2T
Fig. 3. Attachment points for the head. C4. T3 and torso. I I Multiae attachments at one kndpoint A Articular facet stiffness D Disk stiffness L Ligament element M Muscle element R Pure damping element.
number
and type
+M31
T3
soft tissue element
~____._ Mass no.
Pure stiffness
1 2 3 4 5 6 7 8 9 IO 11
12 12 12 12 12 12 12 I4 I4 6
6
Pure damping 3 6 6 6 6 6 6 6 5 6 4
Muscle element 37 9 9 9 10 11 IO 10 9 11 21
214
JAN G. REBER and WERNER GOLDSMITH
Table 3. Effective
masses and inertial
moments
Mass no.
Mass (kg)
Moment of inertia kg - m’ x IO’
1 2 3 4 5 6 7 8 9 10 II
1.860 1.860 0.772 0.772 0.772 0.772 0.172 0.772 2.189 2.189 3.454
28.81 20.79 7.457 4.519 3.390 3.390 2.825 1.808 1.469 1.469 318.2
pot elements in the muscle model were obtained by deriving the stress equation of the model as a function of strain (using the stress-strain relation of the nonlinear spring) and strain rate. The resulting equation was solved for the strain rate, ri, and the damping constant was varied until the solution approximated the results obtained by Ohnishi (1963) for creep in muscle fibers.
requires that the composite structure can be separated into elements in which the properties of interest can be assumed to be homogeneous. (ii) The treatment of the vertebrae as rigid bodies as implied by (i). This requires that the actual strains induced in the other elements of the system be much greater than those produced in the vertebrae. (iii) The dynamics of the actual system can be approximated by a two-dimensional analysis. This requires mid-sagital plane symmetry and specification of the orientation of the body elements parallel to this plane. (iv) The anisotropic physical properties of complicated body substructures can be approximated by combinations of uniaxial elements. This is accomplished by inclining the individual modeling elements of the substructures at different angles so as to allow mutually perpendicular motions of the total substructure to induce different defqrmations in the individual elements.
STRUCTURAL ANALYSIS AND EQUATIONS OF MOTION
Reference frames ASSUMPTIONS
Certain assumptions have been employed to reduce the complexity of the structure. These are: (i) Lumped-parameter-analysis applicability. This
The formulation of this model requires the use of three sets of two-dimensional reference frames: two sets of inertial frames and one moving frame set (Fig. 4). The first inertial frame is global and locates the origins of the second inertial frame set (referred to as
Y
Sm MOVING FRAME
m
\
I
MASS m ,
-
, GLOBAL
ORIGIN
Fig. 4. Reference frame relationships.
-.
715
Large head -neck motions
the local inertial frames) with respect to a common This point is arbitrarily positioned 52.5 mm anteriorly and 7.5 mm inferiorly to the center of gravity of vertebra T3. The local inertial frames (one for each mass segment) are located with their origins at the centers of gravity of the mass segments when the system is initially at its no-load equilibrium position. These frames serve as references for the motion of the individual mass bodies; the inclinations of their axes are identical to that of the axes of the global frame. This implies that only x and J coordinates are necessary to locate each local frame within the global plane. The third set of frames is attached to the mass segments, and translates and rotates with the segments. The geometry of the vertebral bodies and the attachment points of the spring-dashpot systems are defined within these moving frames delineated by two translational and one rotational coordinate relative to the local inertial frames. This arrangement determines the position of any point in the model by the addition of the vectors which locate the frames described above with respect to the global origin.
where
point.
Equations
of motion
The motion of the system is determined by the application of Newton’s second law. To evaluate the forces developed in the individual elements, the total deflection of each spring and the velocity difference between the endpoints of the damping units must be computed. The components of the displacement, dij, between any two arbitrary points,,j and k, on masses i and m, respectively, are given by
vi; = ctj = rij fli sin(Oi + iJij) - ci - r,& 0, sin(8,
+ $,,&) + &
(4)
and “Tj f $ = rmk 0, cos(U,
+ IJJ,~)
+ C, - [rij Bi cos(0, + tiij)
+ Pi].
The velocity, Vc. which is directed along dij (the line of action of the two damping elements), is the projection of Vii on dij and is given by
dij vt=Vij’m,
or
(5)
The corresponding
force, .Fij, for linear damping
is
.P,j = CijV$
(6)
where Cij is a damping constant. The torques generated by these forces about the centers ofgravity of the individual vertebral bodies are expressed as the vector cross product of the forces and the distance to the points of application. for example Ti, = rij x f,, = rijfTj cos(6, + tiij) - rijffj
sin (Oi + I/I,,).
(7)
The equations of motion are developed by summation of the forces caused by all the elements attached to each mass segment or
tr:, = rrnkcos(Q, + tirnk) + u, + p: - [ri,Cos(O, +
$ij)
+
Ui +
pf].
and
(1) d:, = rmt sin (0, + $I& + L’, + -
[rij
Sin(Bi
+
$ij)
+
I‘, +
p’,
pf],
where r and j/ are the polar coordinates of the pointsj and k in their respective moving reference frames for masses i and m. The coordinates u, L’and 0 locate the origins of the moving frames with respect to those of the local inertial frames and p” and pY are the components of the displacement from the global origin to the respective local inertial origins (Fig. 4). The force fij developed in the stiffness element attached to point ij is then expressed in terms of this displacement as fij = Jj(Kij, dij).
to point j on mass i is
4, = 0;.
(8)
where Mi and Ii are the mass and inertial moment of the ith mass segment, a and a are linear and rotational accelerations, and ‘/,, Y and 0 are the two linear and one rotational state velocities, respectively. There are 11 sets of such relations, one for each mass segment (i = 1-l 1) in the system, providing a set of 66 first-order nonlinear differential equations describing the system motion.
(2)
where xj represents a functional relationship with respect to point j on mass i and Kij represents the stiffness parameters of the element. The magnitude of the velocity, Vij, of point k on mass m with respect
I& = ‘N,, Ci = Y ‘,, and
Substructure
delineation
The individual substructures used in the model are listed in Table 1 along with the constitutive equations employed. The ligaments are represented as nonlinear springs and the muscles as Kelvin-Voigt solids. The parameter values of these elements were determined as described previously. The nonlinear anisotropic be-
216
JAN G. REBERand WERNERGOLDSMITH
havior under strain exhibited by the intervertebral disks was modeled by the structures indicated in Table 1. The two crosslinked members allow sufficient freedom in the choice of stiffness and damping functions to permit good matching of the overall response of the unit to the actual disk force-deflection curves. As a first approximation, all the elements in these models were assumed to be linear. Because of geometric considerations, the complete response of the structure is nonlinear and anisotropic, a necessity for proper disk representation. The synovial joints at the interfaces of the articular facets are represented in the model by two springs and two dampers (Table 1). One disadvantage of this model is the occurrence of a phenomenon termed “snap through.” Because the forces are calculated strictly from the deflection and velocities of the individual elements without regard to the actual mass boundaries, there is a possibility that, in the current analytical system, the articulating surfaces of the facets could “snap through” each other. To circumvent this problem, a contact criterion was developed for the articular facet surfaces. Contact between the two line segments representing the two surfaces of the articular facets was considered to occur only if these intersected within the segment lengths (Fig. 5). This criterion can be expressed mathematically by IG
where R indicates a displacement with respect to the local inertial frame referring to mass segment i, subscripts 1 and 2 denote the two line segment endpoints on lines i and m, and subscript “con” designates the point of intersection. If both these inequalities are satisfied, contact has occurred. For this analysis, the coordinates of R,,, must be determined from the equations for the two lines containing the relevant segments, yielding in terms of slope m,
R&,6,= $R:..
- R:,) + R:,.
(10)
‘I (F sin a)’ dy
+
1
2
A2E2 L2 [F(L2
s0
-
y)cosa12dy E212
Ln(F~~sa)2 .
,
mi
These equations are evaluated in the computer program at every time step for each facet junction and the contact criterion is checked. To calculate the force due to contact, the pedicles supporting the facets are assumed to be elasticallydeformable objects, rigidly attached to the body of the vertebra. The pedicles are represented as two round bars as shown in cross-section in Fig. 6. The principle of strain energy is applied to this representation to deduce an equivalent stiffness for the deflection at the facets. The total strain energy, V, in the upper pedicle is evaluated as
I0
“i
mj+-
and
fJ=)
I&-R:,
1
>( >
mi
-K,
and
Rml + R&m,,, R&, - R& + -
R&. =
dx
+t s0 L8 [(F
AI&
sin a)&
- x) - F(cos a)L,]’
dx
+* s 0
~511 (11)
= LOCAL INERTIAL ORIGIN i
“CONTACT”
where A is cross-sectional area, E is Young’s Modulus and F is a force (dimensions L are found in Fig. 6). The value of V for the lower pedicle is given by changing the limits in the first two integrals to read from -L/2 to zero and inverting the sign of the second term in the last integral. From strain energy theory
(12) and after integration and partial differentiation of equation (1 l), the deflection, 6, is found to be equal to the product of the applied force, F, and a constant. The reciprocal of this constant is interpreted as an effective stiffness, K,,. The contact force generated is then evaluated as F = 6K,,,. “NO
CONTACT”
Fig. 5. Contact
criterion.
(13)
The deflection 6 is interpreted as being half the perpendicular distance from the facet endpoint closest
217
Large head-neck motions
neck as a three-parameter viscoelastic solid calculate the this model to and used head-neck-junction displacement that was then compared to the experimental data. Figure 7 shows the head accelerations predicted by the present model for which no corresponding test data is available. Figure 8 shows the predicted motion of the head-neck junction and the results obtained by Landkof (1976). It may be noted that the present model quite faithfully reproduces the observed motion of the junction. Huston et al. (1976) presented data from experiments with seated cadavers subjected to pendulum
artificial
/
E, .I, . A,
/, / 1
Ll LOWER ARTICULAR
/_
J
LI
/”
PROCESS
rl --
/
---1 E, 91, ,A, _i
I
r2
/
L2
UPPERARTICULARPROCESS
Fig. 6. Model approximating pedicle stiffnesses.
to the intersection point to the other intersecting line segment (Fig. 5).
RESULTS AND DISCUSSION t, TIME, ms
A computer program was written incorporating the above equations in Fortran IV that was solved on a CDC 6400 system employing a fourth-order Runge-Kutta scheme for the integration. The inputs from three sets of experimental runs were used in the computational procedure to predict the kinematic histories of these tests. A comparison of the predicted and experimental results allows an assessment of the validity of the model. The first comparison was based on the data of Landkof et al. (1976) who conducted non-destructive radial impact experiments on a thin water-filled spherical lucite shell attached to a viscoelastic articular artificial neck mounted on a rigid base by a suspended thin shell of aluminum. The input force to the system was measured by a sandwiched quartz crystal yielding a contact variation given by F(t), N =
nt
1272 sin .s ( >
o
= 0
t 2 7,
Acceleration components of the head for Run 1.
(14)
where r is the contact duration of 2.5 ms. Among other measurements, the displacement history of the junction between the lucite sphere and the artificial neck was determined. Landkof (1976) approximated the
1, TIME,ms
Fig. 8. Displacement components of Cl for Run 1.
JANG. REBER and WERNER GOLDSMITH
218
head impact and a corresponding open-chain type system analysis. The representation used for the muscles and ligaments is similar to that of the present work. The reference does not indicate the methodology of representing the attachments and uses a different scheme for denoting the disks and articular facets than employed here. Figure 9 portrays the reported pulse input and Fig. 10 displays the experimental displacement, velocity and acceleration of the head in the x-direction as well as that predicted both by Huston ef al. (1976) and by the present analysis. The
-
PRESENT
---
EXPERIMENT,
I
0
MODEL HUSTON
2
3
t , TIME,
Fig. 9. Input
IO a
et 01
4
5
rns
pulse for Run 2.
-
-
PRESENT MODEL - MODEL, HUSTON et 01
-
-
-
EXPERIMENT.HUSTON
agreement of the present model can only be regarded as fair with respect to magnitude when compared with the prediction reported by Huston er al. (1976), although the shapes of the responses are generally similar. This could be a result of parameter inaccuracies in the present model with respect to Huston’s experiment such as that involving the mass and inertia of the head. However, the experimental acceleration curve reported by Huston seems inconsistent with the reported input force. The reported acceleration curve implies that the maximum positive x component of force occurs at 6 ms, after the end of the 5 ms duration impact. A second input force cpnsisting of a linear rise to 8 KN in 5 ms followed by a sudden drop to zero force level was used in the present analysis to attempt to duplicate this effect. The present model again predicts acceleration reversal after the end of impact as seen in Fig. II. However, it can be noted that the x components of velocity and displacement are now in much better agreement with the experimental data due to the large impulse imparted by the ramp input function. The last comparison will be made with the data of Ewing and Thomas (1972) where human volunteers strapped in a sled on a high-speed track were subjected to accelerations at levels comparable to those causing hyperextensions. The data, obtained by means of accelerometers and high-speed framing camera photographs and reported in terms of the motion of the head relative to that at the top of the torso, are considered to be very reliable. The input is given as a linear rise to an acceleration of 7.4 g at 14.2 ms followed by a linear decay to zero at 340ms. Figures 12, 13 and 14 portray experimental accelerations, velocities and displacements of the head from
et d
q-/
,
,
,
,
ar
‘r_t,-..------‘--;: 0 I
IO
5 I
I
t.
15 I
20 I
25 ,
TIME.ms
Fig. 10. Comparison of predicted and measured histories of the horizontal (x) components of displacement, velocity and acceleration
for Run. 2. Original
input.
t. Fig.
IL Results corresponding
TIME,
fns
to Fig. 10 for
revised
input.
219
Large head- neck motions
__
MuGEL
---
FXPERIMEhl
I
0 0 ,
//
100 / t,
Fig. 12. Comparison of predicted and measured acceleration components of the head for Run 3.
3r
-
MODEL
---
EXPERIMENT
1
,
50
TIME,
150
I-zoo
250
rns
Fig. 14. Comparison of predicted and measured displacement components of the head relative to T3 (base) for Run 3.
T!L pq#-3 L
;_;i_i:‘!-:;,,;:+
,
..__,’ 3 =
-2 1 0 I
50
100 (
150 ,
200
250
I, TIME,ms Fig. 13. Comparison of predicted and measured velocity components of the head relative to T3 (base) for Run 3.
Fig. 15. Deformationofthestructurefor
Ewing to Thomas (1972) compared to the predictions of the present analysis. Agreement is considered to be generally good except for high frequency spikes in the curves which might be due to the contact criterion built into the articular facet analysis. Figure 15 shows the predicted maximum excursion of the spinal column for this case where numerous contacts between the articular facets are occurring. A force analysis was executed for the run corresponding to the data of Ewing and Thomas (1972) for three instants during the
motion, including the time of maximum excursion of the spine. The force histories for the two ligament and muscle elements, that exhibit the largest magnitudes for all such elements at these three instances, respectively, have been plotted in their entirety in Figs. 16 and 17. The stress peaks of 10 and 4 MPa correspond to strains of 3.5 and 1.1, respectively. The strain for the first case is believed to be too high by at least a factor of 2 due to an inaccurate representation of the ligament area in the physical model.
Run 3 at
f = 140ms.
220
JAN G.
REBERand WERNERGOLDSMITH
limit beyond which the soft tissues should not be loaded without serious potential for injury. It was found from the present analysis that the supraspinous ligament as well as the three ligaments attaching Cl to the head were loaded beyond this bound by as much as a factor of 10 for many ms, as indicated in Fig. 16. The region between the head and Cl also exceed the stated limit in the case of muscles, as well as the longus colli muscle. Typical computed load histories for the disks are presented in Figs. 18, 19 and 20; the maximum loads are within the test ranges reported in the literature (Kulak er al., 1976). The sharp shear reversal noted in Fig. 18 may very well represent a potential source of damage.
LIGAMENT ELEMENT IO -
5-
0
I
-I
a” “c
6
::
LIGAMENT NO. 109
ELEMENT
50 I
100 1
nw kl? d
3
0
I 0 I
t,
150 I
200 I
IO ---AXIAL
250 I
I-’ ,’
---SHEAR [“’
TIME.ms
Fig. 16. Representative histories involving maximum stresses
,J
5
/ ‘I \,
‘1
\ ‘_. ‘1 ‘J\\
of two ligament elements.
MUSCLE ELEMENl NO. 104
n
Fig. 18. Axial and shear force histories for the intervertebral disk at level T2-T3.
-
2r MUSCLE ELEMENT NO. 124
z 2 &I d
I-
I
0.
-I-
O I
50 I
100 I t,
150 1
200 1
250
TIME,ms
Fig. 19. Axial and shear force histories for the intervertebral disk at level C7-Tl.
Fig. 17. Representative histories involving maximum stresses of two muscle elements. ---
Replications for the ligaments in the physical model currently under construction consist of natural rubber, and this material as well as excised muscles (King and Chou, 1976) were tested in quasi-static tension to an upper limit of about 0.5 MPa, well into the nonlinear stress-strain regions of these materials. Although not definitive from a physiological damage point of view, this magnitude was felt to be a representative upper
AXIAL SHEAR
Fig. 20. Axial and shear force histories for the intervertebral disk at level C4-C5.
221
Large head-neck motions One of the most frequent complaints of injury in the cervical region involves trauma to nerves, particularly to those within the spinal cord. As an estimate of mechanical damage due to excessive elongation, a criterion might be placed on the permissible extension of the spinal cord. In the present model, this extension is hypothesized to be the same as that of the posterior longitudinal ligament which is also contained within the vertebral foramen and is explicitly represented in the present model. An examination of the angular deformation and strain experienced by the spinal cord on the basis of the above hypothesis points to three regions whose loadings indicate that they are most likely damaged before others. The first, between C6 and C7, exhibits the largest positive strain of I.52 together with a large angular dislocation of 52.4” at a time 100 ms beyond the initiation of loading. The other two regions are between C2 and C3, and the head and C 1, respectively, where the strains are of the order of 0.4 coupled with angular deformations of the order of 30”. The present system can be modified with respect to the number of elements and element parameters by adjustment of the input data to the system and may also be extended to three dimensions. The latter would require the derivation of general displacement and velocity equations in three dimensions and the inclusion of the third linear acceleration and the two rotational acceleration components to the equations of motion. It is estimated that the modifications mentioned above would increase computation time by a factor of eight.
CONCLUSIONS
A lumped parameter model of the head and the cervical spine region encompassing the first 9 vertebrae fixed on a rigid base has been developed and solved by a fourth-order Runge-Kutta integration scheme. The elements were modeled by combinations of rigid masses, nonlinear springs and linear dashpots and comprise the head, the vertebral bodies, the disks, ligaments, muscles and articular facets. The program was evaluated for three situations where experimental results had been obtained by other investigators. Vertebral contact at the articular facets was handled by means of the action of a contact force when such a situation prevailed. It was found that good agreement was obtained in all cases when the predictions of the
present model were compared with the test data, even though the individual element parameters could not be matched precisely to each of the various cases studied. A mechanical criterion of a limit of damage has been suggested based on test results of both physiological replacement materials and excised human soft tissue. Stresses on soft tissue and loads on the vertebrae were computed, and the areas of most likely damage under the loading conditions described were delineated. It was found that several tissues exceeded this artificial
limit substantially and for an extended period of time, pointing to the possibility of serious spinal cord injury.
Acknowledgements - The authors would like to extend their appreciation to Michael Kabo and Professor R. L. Taylor for significant technical assistance. The work was performed under NIH Grant No. 5ROl AM18384 supplemented by Contract No. NOCO14-76-C-0686from the Office of Naval Research.
REFERENCES
Ewing, C. L. and Thomas, D. J. (1972) Human head and neck responses to impact acceleration. Joint Army--Navy Monograph, No. 21. Goldsmith, W. (1979) Some aspects of head and neck injury and protection. To appear in Progress in Biomechanics, Proc. NATO Advanced Study Institute. Noordhoff. Hodgson, V. (1970) In Impact Injury and Crash Protection (Edited by Gurdjian, E. S. et al.), pp. 2755307. Thomas, Springfield. Hopkins, G. R. (1971) Nonlinear lumped-parameter mathematical model of dynamic response of the human body. Symposium on Biodynamic Models and Their Applications, pp. 843-849. Aerospace Medical Research Lab, WPAFB,
Ohio (AMRL-TR-7 l-29). Huston, J. C., Passerella, C. E. and Huston, R. L. (X976) Numerical prediction of head/neck response to shock impact. In Measurement and Prediction of Structuraland Biodynamic Crash-Impact Response (Edited by Saczalski, K. J. and Pilkev. W. D.), DV.1377149. ASME. New York. Huston, R. L., Huston, J.-C. and Harlow, M’. W. (1978) Comprehensive, three dimensional head-neck model for impact and high-acceleration studies. Aoiat.. Space Environ. Med. 49, 205-210. Kaleps, I., von Gierke, H. E. and Weis, E. B. (1971) A fivedegree-of-freedom mathematical model of the body. Symposium on Biodynamic Models and Their Applications, pp. 211-231. Aerospace Medical Research Lab, WPAFB, Ohio (AMROL-TR-71-29). King, A. I. and Chou, C. C. (1976) Mathematical modeling, simulation and experimental testing of biomechanical system crash response. J. Biomech. 9, 301~317. Kulak, R. F., Belytschko, T. B. and Schultz, A. B. (1976) Nonlinear behavior of the human intervertebral disc under axial load. J. Biomech. 9, 301-317. Landkof, B., Goldsmith, W. and Sackman, J. L. (19761Impact on a head-neck structure. J. Biomech. 9, 141-151. Li, T. G., Advani, S. H. and Lee, U. C. (1971) The effect of initial curvature on the dynamic response of the spine to axial accelerations. Symposium on Biodynamic Models and Their Applications, pp. 5555569. Aerospace Medical Research Lab. WPAFB, Ohio (AMRL-TR-71-29). Liu, Y. I. and Murray, J. D. (1966) A theoretical study of the effect of impulse on the human torso. In Biomechanics (Edited by Fung. Y. C.). vv. 167-186. ASME. New York. McElhaney; J., Roberts, VI *L. and Hilyard, J. (1976) Biomechanics of Trauma. Duke University Press. Durham. North Carolina. Mertz, H. J. and Patrick, L. M. (1971) Strength and resoonse ofthe human neck. In Proc. 15th Stapp Ca;Crash Co& pp. 2077255. SAE Paner No. 710955. Muzzy, W. H. III and Lustick, L. (1976) Comparison of kinematic parameters between hybrid II head and neck system with human volunteers for -G, acceleration profiles. In Proc. 20th Stapp Car Crash Co@. pp. 45-74. SAE Paper No. 760801. Ohnishi, T. (1963) Rheology Biorheoloyy 1, 83-90.
of glycerinated
muscle
fibers.
JAN G. REBERand WERNER GOLDSMITH
222
stiffness coefficient length mass displacement in local inertial plane torque strain energy velocity acceleration displacement stiffness force slope location of local inertial origin in global plane displacement in moving reference frame time coordinates of mass center of the segments in the local inertial frame damping force damping torque velocity components of mass segment centers in local inertial frame angular acceleration deflection angular displacement in local inertial frame angular velocity in local inertial frame contact duration angular velocity angular displacement in moving reference frame.
Ome, D. and Liu, Y. K. (1971) A mathematical model of spinal response to impact. J. Biomech. 4, 49-72. Prasad, P. and King, A. 1. (1974) An experimentally validated dynamic model of the spine. J. appl. Mech., Trans. A.S.M.E. 96, (E), 5466550. Prasad, P., King, A. I. and Ewing, C. L. (1974) The role of the articular facets during +G, acceleration. J. appl. Mech., Truns. A.S.M.E. 41 (E), 321-326. Reber, J. G. (1978) Lumped parameter model of whiplash. M.S. Thesis, University of California, Berkeley. Rybicki, E. F. and Hopper, A. T. (1971) A dynamic model of the spine using a porous elastic material. Symposium on Biodvnamic Models und Their Applications, pp. 851-875. Aerospace Medical Research Lab, WPAFB. Ohio (AMRL-TR-71-29). Shirazi, M. (1971) Response of the spine in biodynamic environments. Symposium on Biodynamic Models and Their Applications, pp. 843-849. Aerospace Medical Research Lab, WPAFB, Ohio (AMRL-TR-71-29). Terry, C. T. and Roberts, V. L. (1968) A viscoelastic model of the human spine subjected to +G, accelerations. J. Biomech. 1, 161-168.
Toth, R. (1966) Multiple-degree-of-freedom, nonlinear spinal model. In Proc. 19rh Ann. Co@ Engng Med. Biol. 8, 102. Wickstrom, J., Martinez, J., Johnston, D. and Tappen, N. C. (1965) Acceleration-deceleration injuries of the cervical spine in animals. In Proc. 7th Stapp Car Crash Conference (Edited by Severy, D.), pp. 284-301. Thomas, Springfield. Yamada, H. (1970) Strength of Biological Materials (Edited by Evans, F. G.). Williams & Wilkins, Baltimore.
Subscripts i.nt j.k
EO” co NOMENCLATURE
A C E F I
area damping coefficient Young’s modulus force moment of inertia
mass segment number attachment point number contact effective total.
Superscripts x.)’ in x, $’direction
I I SUP
functional relationship time derivative absolute value supremum.