- Email: [email protected]
Three-dimensional mathematical model of the human spine structure
J. Riomcxhonics,
1973, Vol. 6. pp. 67 I-680.
Pergamon
Press.
Printed
in Great Britain
THREE-DIMENSIONAL OF THE HUMAN
MATHEMATICAL MODEL SPINE STRUCTURE”?
MANOHAR
M. PANJABI
Section of Orthopaedic Surgery, Yale University School of Medicine, New Haven. Connecticut 065 10. U.S.A.
Abstract-A method is presented for constructing
discrete-parameter type three-dimensional mathematical models and governing equations of motion of the spine structure. The anatomic structure is represented by any combination of rigid bodies, springs, and dashpots in space. These are positioned, orientated, and connected in a manner to represent the true mechanical function of the structure. The rigid bodies are of any shape and have 6 degrees-of-freedom, allowing three-dimensional motion. The springs and dashpots may have up to twenty-one stiffness and damping coefficients respectively to precisely represent the three-dimensional coupled behavior. The method is straightforward and simple to apply. The governing equations are in the matrix form and are easily generated and solved by computer techniques.
INTRODUCTION
The continuum models have been developed by Hess and Lombard (1958); Liu and Murray (1966) and Terry and Roberts (1968). The latest model of this type considers the spine in the form of a viscoelastic rod of Maxwell type material with no end-masses. The discrete-parameter models have been developed by many people and only some are mentioned here. Latham (1957) was the first one, followed by Payne (1962); Toth (1967). and finally by Ome and Liu (1969). The latest discrete-parameter model considers the spine as consisting of a chain of rigid bodies connected in a series by springs and dashpots. Each vertebrae is given 3 degrees-of-freedom, i.e. horizontal and vertical translations and rotation in the sagittal plane. Springs and dashpots are given three stiffness and damping coefficients respectively. The external loads can be applied only on the two end-bodies. Although the model of Ome and Liu has come a long way from Latham’s one degreeof-freedom model, it is still quite restricted in its representation of the real spine structure. It is a series model and thus cannot
models of the spine are of two kinds: the continuum and the discreteparameter type. The former considers the spine as a rod having an infinite number of degrees-of-freedom. It is, however, too different from the real spine structure because of its limitations in taking into account such important anatomic features as the difference in properties of vertebrae versus disc, the coupling of axial rotation of vertebrae with lateral bending, and the effects of ribs and muscles. On the other hand, a discreteparameter model considers the spine as a structure formed by various anatomic elements like vertebrae, ligaments, muscles, and articulating surfaces. Vertebrae are represented by rigid bodies and ligamentous structures, muscles, and articulating facets by massless springs and dashpots. A high degree of sophistication can be accomplished in representing various anatomic elements mathematically and in applying loads to the bodies. This is the type of model which is a good deal closer to the real spine structure. MATHEMATICAL
“Keceicrd27
Nocrmber
1972.
+Work Supported by: Fluid Grant PHS-RR-05358-10. The Commonwealth of New York.
The Crippled Children’s Aid Society. The Yale Program and 671
672
MANOHAR
represent muscles and ribs. Loads can be applied only at the end-bodies. Rigid bodies are constrained to move in the sagittal plane; therefore, no three-dimensional studies of the spine can be carried out. These and other restrictions make it impossible to study with the presently available models some very basic biomechanical problems in orthopaedics. A few examples of such problems are the three-dimensional evaluation of the classical problems: the pilot ejection and the whiplash, the complicated automobile accidents, athletic injuries, the effects of clinical instability, vertebral fusion, description of scoliosis, and the evaluation of the relative efficiency of the various surgical and nonsurgical correcting techniques of the spine. The general mathematical model presented here has the capabilities of studying these and similar problems. THE MODEL
General description
The basic procedure for constructing any mathematical model of a musculoskeletal system is to formulate a set of governing equations of motion representing the mechanical behavior of the model to external loads and displacements. A general method is presented here which makes it rather simple to write the set of governing equations of motion for a discrete-parameter type mathematical model of a musculoskeletal structure. Each rigid body has all the 6 degrees-of-freedom (allowing general three-dimensional motion) and each connecting element (spring and dashpot), all the twenty-one coefficients of stiffness/damping. The model is basically restricted to linear geometry and material properties, but these restrictions can be relaxed to a great extent by special computer techniques. This work is based upon the theories developed by Panjabi (197 1). The spine structure to be studied is idealized to a discrete-parameter type mathematical model. All the anatomic structures are idealized by the three mathematical ele-
M. PANJABI
ments: rigid bodies, massless springs, or massless viscous dashpots. The features of the model and its elements are described below. 1. The model may consist of any number of bodies, springs, and dashpots connected to each other in any manner required to form an anatomically true model. 2. The bodies are three-dimensional rigid bodies and each may have all of the possible 6 degrees-of-freedom of motion. They may be of any size and shape and positioned and orientated in space in any manner to simulate truly the anatomic features of the bony elements they represent. No symmetry of any kind is required. 3. The springs are massless and each may have all of the possible twenty-one (six main and fifteen coupling) spring stiffness coefficients. However, it may suffice to work with a much smaller number of coefficients due to symmetry of some kind. The end conditions may be any combination of fixed, hinged, and ball joints. The springs may be positioned and orientated in any manner to simulate true anatomic features. 4. The dashpots, similar to springs, are massless and each may have all of the twentyone damping coefficients. Again, due to symmetry, it may suffice to work with a much smaller number of coefficients. End conditions may be combinations of fixed, hinged, and ball joints. Position and orientation can be tailored to simulate true anatomic features. 5. The dashpots and springs may be connected in any combination, for example: in parallel to form Kelvin viscoelastic solids or in parallel cum series combination to form the 3-parameter solids. 6. The system may have any initial configuration. 7. External loads, static and dynamic, consisting of force and torque vectors may be applied singularly or simultaneously to any number of the bodies. 8. Matrix methods are used throughout which make the method versatile and suitable for the computer.
MATHEMATICAL
MODEL
OF THE HUMAN
Every method has limitations, and it is important to know these so that one may not only restrict its applications to suitable problems, but also may find methods to relax these limitations. The underlying assumptions of the method are 1. Displacements are small. It can, nevertheless, be relaxed to a great extent, whereever required, by dividing a large displacement into small displacement-steps. The resulting configuration of the system, after a small step, becomes the initial configuration of the system for the next step. 2. Bodies are rigid. 3. Springs and dashpots are massless and linear. Again, this assumption of linear characteristics can be relaxed by adopting small displacement-steps technique. After every step, not only the geometry is altered, but also the stiffness and the damping properties are changed to suit the nonlinear characteristics. Defining
the structure
SPINE
STRUCTURE
613
I
Fig. 1. This shows the method of representing a given musculoskeletal structure by a discrete-parameter type mathematical model.* Rigid bodies are numbered in sequence from 1 upwards and each has local body coordinate system (~yz)~, (~yz)~,. . Connecting elements (springs and dashpots) are also sequentially numbered and have their local coordinate systems (abc),, c&c),... Orientations of the connecting elements are such that the u-axes always point towards the rigid bodies with the higher sequence numbers. Finally there is the global coordinate system (Imn) which defines all the other coordinate systems.
A general model consisting of bodies, springs, and dashpots is shown in Fig. 1. Three types of Cartesian right-hand coordinate systems are used to define the structure: global system, local body systems, and local spring/dashpot systems. The global system, called (fmn), is conveniently orientated and direction towards the body with the higher fixed to the ground. It defines the centers of number. The other two axes, for the sake of gravity (c.g.) of all the bodies and the orientaminimizing the number of coupling coeffictions of the other coordinate systems. The ients, may be orientated in the directions of bodies are numbered from one upwards, with the principal elastic/damping axes of the ground being zero. Local body coordinate spring/dashpot respectively. system (x.v.& is fixed in space at the c.g. of Figure 2 shows the ith body isolated from the ith body and orientated along its principal the rest of the system. Center of gravity of the axes of inertia. (This is required only for the ith body is defined by position vector oi in the study of the dynamic behavior of the system.) global coordinate system (Imn). The points of Springs and dashpots are also sequentially attachment of springs/dashpots are defined numbered. The kth springldashpot has the in local body coordinate system (XyZ)i. A /ith local spring/dashpot coordinate system position vector cljlcdefines the position of the called (crb& It is orientated so that the CZ~- end of the kth spring/dashpot attached to the axis is parallel to the line joining the two ends ith body. The external loads applied to the of the kth spring/dashpot and with its positive body are also defined by similar position vec..*Bony elements are represented by rigid bodies while connecting elements like ligaments. muscles, articulating facets and joints are represented by springs and dashpots.
674
MANOHAR
M. PANJABI
load, di and qj, for each body and in the respective local inertia systems. Thus d=
[d, j dz j-j
di i -/
d,Jr,and
4=[91~q2j-~9i~--i9nlT. The subvectors, components
in turn, are defined by their
[rxr~z+&hl iT, 9i = LfxfiP,P~P~liT di =
where r, 4, f, and p are respectively translation, rotation, force, and moment vectors. The system matrices M, K, and C and the vector q are now treated. Fig. 2. This shows an arbitrary body of the musculoskeletal system, called the ith body, along with one connecting element (spring/dashpot) and one external load, isolated from the system. Local body coordinate system (xy& is fixed in space at the center of gravity i and orientated along the directions of principal inertia axes of the body. Positions of the connecting elements and the external loads are defined by radius vectors Uik and e,K respectively where k is the sequence number of the connecting element.
1. System inertia matrix
The local inertia coordinate systems are orientated along the principal directions of inertia of the respective bodies. This leads to a diagonal system inertia matrix. Thus
(2) tars. eik is the position vector for the point of application of the kth external load vector on the ith body. where System governing
submatrices
for each of the bodies
equations of motion
The procedure described here is given without the underlying arguments. For details, one may consult the original work Panjabi (1971). The system governing equations of motion for any spine structure in the set of local inertia coordinate systems (XYZ)~to (xyz), in the matrix form is tid+Cd’+Kd=q
(1)
where d and q are 6n X 1 system displacement and load vectors respectively. M, C, and K are 6n X 6n system inertia, damping and stiffness matrices respectively. Letter n stands for the number of bodies in the system. The system vectors d and q are defined by a set of subvectors of displacement and
Here, mi stands for the mass and I,, Z,, and stand for the principal inertias of the ith body. Zzi
2. System stifness matrix System stiffness matrix matrix equation
is given by the
K =
(4)
MATHEMATICAL
where submatrices Kij
MODEL
OF THE HUMAN
SPINE
STRUCTURE
675
Kij (6 X 1) are given by =
x k
(5)
K,k.
Each individual submatrix Kz belongs to bodies i and j and the spring k. The rule for summation in equation (5) is If i # j; summation is taken over all the springs connecting the body i with j, and (6)
d*
If i = j; summation is taken over all the springs connected to the body i Further, submatrices Kzare defined in terms of the basic matrices (6 X 1) by matrix equation K”, = ,rJ.“?‘@.r~krS(E.~‘QJJ.~ I 1 (7) 1.1 2.1 ., .r . Here some explanation of the function of individual matrices on the right hand side of equation (7) is justified. We will apply, for illustration’s sake, this equation to a single body attached to the kth single spring, as shown in Fig. 3. Its application in multibody systems is similar but more complex. Thus, equation (7) becomes K” =
~k’l’@T~kT~k~k@
uk_
(8)
If d,,,,, is the displacement produced at the the center of gravity of the body in the xyz coordinate system by the load qtruz) applied to the body at the center of gravity in the same coordinate system, then q(.rwz,= Kkd,,w,,.
(9)
Thus, four basic matrices Uk, Cp,q’” and Sk, of equation (8) represent transformation of vector dClyz)to vector qtS.+ Their individual functions are described below. Uk is the position matrix. It transforms vector 4,,,, to &.vz). The latter is the body displacement at the point of attachment. 4r is the rotation matrix. It transforms vecThe latter is the attachment tor &,,, to dfElmnr point displacement in the global system (Imn). qk is another rotation matrix. It transforms
Fig. 3. Displacement of an isolated single body attached to an arbitrary kth spring is analyzed. Translation and rotation of the center of gravity represented by displacement vector d,,,, produces displacement dk of the spring attachment point. Reaction of the spring is to apply force and moment vectors represented by load vector q* to the body at the attachment point. qczuz,is the load vector at the center of gravity which is equivalent to q” at the spring attachment point.
vector dfT,,, to d&bc,k.The latter is the attachment point displacement in the local spring coordinate system (abc)k. Sk is the spring stiffness matrix. It transforms the vector c&,+ to c&~)~. The latter is the load vector acting at the attachment point on the Spring in the ( UbC) k SyStem. UkTaTqkT is the transpose of ‘Pk@Uk. It transforms q::lhcjk to qtsuxjr the vector on the left hand side of equation (9). Returning back to multibody systems, contents of matrices U”, 0, and qk for both indices i and j are, omitting for convenience indices i and j:
MANOHAR M. PANJABI
676
a=
cosxl cos yl cos xm cos ym cos xn - cos yn _--_
1
cos 21 I cos zm 1 cos zn 1
ZERO
Tk=
cos al
cos
(cosxm
cos cl cos cm ------------_ ZERO
-1
k
ZERO
(15)
S$=Sk
(16) 3. System
+1
ZERO
ZERO
-1
-1
+1
and I 7
k
(18)
I cos an cos bn cos cn1 .
damping matrix
The method of setting up system damping matrix C is the same as that for the system stiffness matrix K. Therefore, equations (419) and Table 1 give the system damping matrix, if word “spring” is exchanged for “displacement” for “velocity” “dashpot”, and “stiffness’ by “damping”. 4. System load vector Let q(lmnl be the system external vector defined in the global system. Then 4(Z?n?l, = c41 i -i
-s11s12 s22
s13
s14
s15
823
St4
825
333
334
s35
s44
s45
SYMMETRIC
s55
4i i - i s;ll~,,,
load
(20)
where qi is the load vector acting on the bodies 1 to n respectively. There may be more than one load vector on a single body. Thus,
The local spring stiffness matrix Sk is Sk =
(12)
The elements sll, s12, - - - x66are the stiffness coefficients along the 6 degrees-of-freedom and their coupled motions. Table 1 gives the information regarding the matrix elements of the basic matrices of equations (10, 11, 12, 17, 18 and 19). Knowing these elements, we can construct the system stiffness matrix K as defined by equation (4).
Sk51 = [‘sq’ 13
where -1
, and
cos 02’L___----__-_--_i cos al cosam I cos bl cos bm I cos cl cos cm
But matrices Si for the kth spring connecting bodies i and j, are different for different combinations of iand j. Here i < j. These are:
I’=
cos zm
-7,
cos an4 cos bnj
I
1
cos ym
, cos x12
cos am cos bm
bl
(11)
---t-------, cos xl
k 9i
(19)
The
=
C k
4ik-
(21)
point of application of qi”, the kth load vector on the ith body, is defined in the local body coordinate system (XyZ)i by the position vector eik.
MATHEMATICAL
MODEL
OF THE HUMAN
SPINE
STRUCTURE
611
Table 1. This table gives information about the elements forming the basic matrices which in turn form the system stiffness matrix K. Equation
How to obtain
Matrix elements
( 10)
lli-, uy, 4
(11)
cos xl, cos yl
(12)
Cos al. cos am,
(16)
s,,,s,~,.. . . . . . . . . . . . . SIX
(19)
1
.. . ..
These are the components of the position vector from the center of gravity to the attachment point in the xyz coordinate system. These are the cosines of the angles between axes x and 1,y and 1,. . . , i and n respectively. These are the cosines of angles between axes a and 1, a and m, . . . . . , c and n respectively. These are the 21 stiffness coefficients representing the stiffness properties of a spring at its j-end in the spring coordinate system. Length of the spring
, cos zn , cos CII
Therefore, transformation of 9(1mn,to the local body coordinate systems xyz and to the respective centers of gravity are given by matrix equation
system behavior under a given set of loading conditions. For a static load vector 9, the solution will be the elastostatic displacement of the bodies of the system, if the effect of damping matrix C is neglected. On the other (22) hand, if transient behavior of the system is of where interest, damping matrix C is retained and a constant load 9 is suddenly applied. In a @= certain study 9 may be a function of time; then the solution will be the dynamic behavior L _ -)- L g?iJ’r-, (23) of the system. In all such cases, solutions of the equation (1) will be in the form of system ZERO L--+_ displacement vector d. Once the system displacement is obtained, the system velocity and and acceleration can be easily computed. r Knowing the displacement and velocity vectors, loads or stresses at the attachment (24) - - - -$y_: __ points of springs and dashpots can be comZERO ;__‘iL_ _ puted. We shall now obtain these relations. Displace the c.g.‘s of the two bodies i and j as shown in Fig. 4 by displacement vectors The submatrices 4ii are given by equation (11) d (syz). and 4zj,z,T Let 9kho and 9tcnbr) be the and submatrices Eii are given by internal loads developed ok the kth sgring at the attachment points. Then we have Eii = 1
_*__
1
_q=o
L
0 0 0
0 - ezi
e,,
0
- egi
e
L
1
.ri
9h’= Q”d
1 ZERO eyi 1 -esi
O
l
0
0
0
(261
where 1
35) I
System behavior
Solution of the system governing equations of motion, i.e. equation (l), represents the
(29)
MANOHAR
Fig. 4. An arbitrary spring k connecting two arbitrary bodies i and j is shown. Unequal displacement of the two attachment points produces forces and moments, represented by load vectors qik and qjk, on the two ends of the spring.
The submatrices
Qi to Qi are given by
Qt = St’Pk@jUjk.
(30)
The basic matrices of equation (30) have been discussed earlier and are given by equations (10-19). Same method is applied to determine the loads at the dpshpot attachment points by replacing d by d and Sk by Dk in the equations (26-30). THE SPINE STRUCTURE
The spine is part of the musculoskeletal system. In a true dynamic analysis, the spini: cannot be isolated and studied separately from the rest of the skeletal system. A good dynamic analysis would at least include the spine structure. Figure 5 shows the spine structure and its various coordinate systems.
M. PANJABI
Fig. 5. The spine structure and three types of coordinate systems required to completely detine its position and motion. Global coordinate system imn is tixed to the ground. Local body coordinate system (XyZ)iis fixed at the center of gravity of rigid body i and orientated in the directions of the principal inertia axes. Local connecting element coordinate system (ab~)~ is orientated along the principal elastic/damping axes of the kth connectingelement. The ak axis points towards the rigid body i+ 1.
For the sake of clarity muscles are not shown, but they may be included in the analysis. Bony elements of the structure (pelvis, vertebrae, head, and sternum) are represented by rigid bodies. Connecting elements joining the bony elements consist of ligaments, articulating joints, cartilage, ribs, and muscles. In general, connecting tissues have viscoelastic properties and are most suitably represented by either a Kelvin type solid or the three-parameter solid. In the latter case, junction between the Kelvin solid and the series spring may be considered as a body of zero dimensions. For the details see Appendix. Single ligaments may be represented by uniaxial three-parameter solids which resist only when pulled. On the other hand, complex behavior of the vertebral motion segments may be represented by three-dimensional solids, Articulating joints “three-parameter” may best be approximated by fictitious springs with zero stiffness in the plane of
MATHEMATICAL
MODEL
OF THE HUMAN
contact and large stiffness in all other directions. As a first approximation passive muscles may be represented by Kelvin or threeparameter uniaxial solids combined with external (muscle) forces acting on the bodies at the points of attachment. The effect of the rest of the musculoskeletal system and the outside constraints (safety harness, seat, etc.) may be neglected or taken into account by including these within the system of analysis. DISCUSSION
No illustrative examples have been given here because a good example will need too much space, while a simple example will not do justice to the generality of this method. The method, however, has been applied in two cases in the field of mechanical engineering. These are the vibration study of a gramophone, an elastodynamic system, panjabi (197 I), and the study of dynamic stresses in the bali bearings of a washing machine, a viscoelastic dynamic system, Eriksson and Jansson (1972). I am confident that the same general method could be used successfully to study the mechanical behavior of complex musculoskeletal systems. Because of its complex three-dimensional geometry and its limited motion, the spine structure is a good example of such a system. The model presented here is general. It can be used with any data: calculated strengthof-material formulae, incomplete experimental data available at present, or the three-dimensional, including coupling effects, experimental data which we hope will be available soon. Of course, the results obtained from the model will reflect the accuracy of the data. We are working in our laboratory to obtain more complete data of the spine segment behavior. We can simulate complex behavior where
not only the structure is three-dimensional but also the external loads. Force and torque may
be applied
to any point
on any body
in
SPINE
STRUCTURE
679
the system, and we can study the static. transient, and dynamic behavior. Work is in progress on some of the clinical projects, mentioned earlier in the introduction, concerning the behavior of the spine structure. Results will be reported in future publications. REFERENCES Eriksson, S. and Jansson, L. (1972) Forces in freely suspended washing machine. (In Swedish) MS. Thesis. Division of Machine Elements, Chalmers University of Technology, Goteborg, Sweden. Fliigge, W. (1967) Viscoelustirity.BIasdell, Waltham, Mass. Frisen M., M%gi M., Sonnerup, L. and Viidik, A. ( 1969) Rheological analysis of soft-collagenous tissue. J. Biomechunics 2. 13. Hess, J. L. and Lombard, C. F. (1958) Theoretical investigations of dynamic response of man to high vertical accelerations. J. Aoiation Med. 29,66-75. Latham. F. (19.57) A study in body ballistics: Seat ejection. Proc. R. SW. B-147. Liu, Y. K. and Murray. J. D. (1966) A theoretical study of the effect of impulse on the human torso. A.S.M.E. Biomechanics Symposium. pp. 167-186. Ome, D. and Liu, Y. K. ( 1969) A mathematical model of spinal response to impact. Thesis 02370, Department of Engineering Mechanics, Univ. of Michigan. Paniabi, M. M. (197 1) Theoretical treatment of vibrations in”single and multiple body suspension systems based on matrix methods. Thesis, Division of Machine Ele-
ments, Chalmers University of Technology, Goteborg. Sweden. Payne, P. R. (1962) The dynamics of human restraint systems. In Impact Acceleration Stress. National Academy of Sciences, National Research Council. Ptiblication 977, pp. 195-258. Terry, C. T. and Roberts, V. L. (1968) A viscoelastic model of the human spine subjected to i a: accelerations. J. Biomechanics 1, 16 1. Toth. R. (I 967) Multiple degree-of-freedom nonlinear spinal model. Presented 19th Annual Conference on Engineering in Medicine rmd Biolo,qy, San Francisco. California. NOMENCLATURE MutriceJ C damping matrix n local dashpot matrix E position matrix for external load I’ 180”rotation matrix K stiffness matrix L force transformation matrix M inertia matrix Q spring-/dashpot-end load matrix s local spring stiffness matrix u displacement transformation matrix cp rotation matrix: xvi to lmir \Ir rotation matrix: lmrr to trhc
680
MANOHAR
M. PANJABI rz
Vectors d
f2
displacement vector position vector for external load ; force vector 4 load vector r translation vector IA position vector moment vector rotation vector f
2
Scalars c damping coefficients I principal moment of inertia of the body 1 initial length of spring m mass of the body n number of bodies in the system s stiffness coefficients i
Indices a, b, c Lj k
axes of the spring coordinate system two arbitrary bodies, i and j arbitrary spring, arbitrary load 1, m, n axes of the global coordinate system x, Y, z axes of the local body coordinate system
Special
(b) (a) Solid; uniaxial response. Fig. 6. (a) Three-parameter Solid as incorporated in a 6(b) \-, “Three-parameter” degree-of-freedom System. A fictitious body h of zero dimensions is introduced to represent the junction of the Kelvin Solid and the Series Spring.
Symbols
- derivative withrespect to time T transpose. APPENDIX Three-parameter
ligament
model
The ligaments may be represented either by Kelvin Solid or by a more complex model, the three-parameter solid. The latter model simulates closely the viscoelastic behaviour of ligaments in the stationery phase i.e. when all deformation is recoverable as demonstrated by Frisen et al. (1969). We can incorporate this ligament model by considering the junction between the Kelvin solid and the series spring as a fictitious rigid body of no dimensions. In this manner it is possible to incorporate into the mathematical model a true three-dimensional “3-paramesolid for representing ligamenteous ter” viscoelastic structure. It is known that the behavior of the uniaxial 3-parameter solid as shown in Fig. 6a is given by the following differential equation c3ji+
(k, + k2)fz = k2cJiZ + k,kpr2.
(31)
Symbols are as in Fig. 6a and formulation is according to Fliigge (1967). To incorporate this model in a threedimensional multi-body system we can build a model as in Fig. 6b. Bodies i and j are connected by the threedimensional “3-parameter” viscoelastic solid consisting of springs 1 and 2, dashpot 3 and the fictitious body h. Here i < h
as body number 1 and the free end of the 3-parameter solid as number 2. The other end of the 3-parameter solid is fixed to ground numbered 0. Then equation (1) becomes Cd+Kd=q
(32)
where
d=
j dJT
[d,
and
q=
j q2]T
[q,
With the use of equations (4- 19) for the simple uniaxial case of Fig. 6a, contents of equation (32) reduce to C=[; d=
;I.“=[“;& [r,
rllT
and
-21, [O f2]‘.
Which when expanded gives (k,+k2)r,-kpre+csr, - k2r, + kzr2
=0
(33)
=h.
(34)
It can be easily seen that the two equations (33 and 34) are the same as the equation (3 1) if r, is eliminated. Thus by introducing an additional displacement for the junction of the Kelvin solid (rl in the uniaxial case) and the series spring we can always formulate system goveming equations in the standard form of equation (1). The price paid for this convenience is the increment in size of the system stiffness and damping matrices.
1973, Vol. 6. pp. 67 I-680.
Pergamon
Press.
Printed
in Great Britain
THREE-DIMENSIONAL OF THE HUMAN
MATHEMATICAL MODEL SPINE STRUCTURE”?
MANOHAR
M. PANJABI
Section of Orthopaedic Surgery, Yale University School of Medicine, New Haven. Connecticut 065 10. U.S.A.
Abstract-A method is presented for constructing
discrete-parameter type three-dimensional mathematical models and governing equations of motion of the spine structure. The anatomic structure is represented by any combination of rigid bodies, springs, and dashpots in space. These are positioned, orientated, and connected in a manner to represent the true mechanical function of the structure. The rigid bodies are of any shape and have 6 degrees-of-freedom, allowing three-dimensional motion. The springs and dashpots may have up to twenty-one stiffness and damping coefficients respectively to precisely represent the three-dimensional coupled behavior. The method is straightforward and simple to apply. The governing equations are in the matrix form and are easily generated and solved by computer techniques.
INTRODUCTION
The continuum models have been developed by Hess and Lombard (1958); Liu and Murray (1966) and Terry and Roberts (1968). The latest model of this type considers the spine in the form of a viscoelastic rod of Maxwell type material with no end-masses. The discrete-parameter models have been developed by many people and only some are mentioned here. Latham (1957) was the first one, followed by Payne (1962); Toth (1967). and finally by Ome and Liu (1969). The latest discrete-parameter model considers the spine as consisting of a chain of rigid bodies connected in a series by springs and dashpots. Each vertebrae is given 3 degrees-of-freedom, i.e. horizontal and vertical translations and rotation in the sagittal plane. Springs and dashpots are given three stiffness and damping coefficients respectively. The external loads can be applied only on the two end-bodies. Although the model of Ome and Liu has come a long way from Latham’s one degreeof-freedom model, it is still quite restricted in its representation of the real spine structure. It is a series model and thus cannot
models of the spine are of two kinds: the continuum and the discreteparameter type. The former considers the spine as a rod having an infinite number of degrees-of-freedom. It is, however, too different from the real spine structure because of its limitations in taking into account such important anatomic features as the difference in properties of vertebrae versus disc, the coupling of axial rotation of vertebrae with lateral bending, and the effects of ribs and muscles. On the other hand, a discreteparameter model considers the spine as a structure formed by various anatomic elements like vertebrae, ligaments, muscles, and articulating surfaces. Vertebrae are represented by rigid bodies and ligamentous structures, muscles, and articulating facets by massless springs and dashpots. A high degree of sophistication can be accomplished in representing various anatomic elements mathematically and in applying loads to the bodies. This is the type of model which is a good deal closer to the real spine structure. MATHEMATICAL
“Keceicrd27
Nocrmber
1972.
+Work Supported by: Fluid Grant PHS-RR-05358-10. The Commonwealth of New York.
The Crippled Children’s Aid Society. The Yale Program and 671
672
MANOHAR
represent muscles and ribs. Loads can be applied only at the end-bodies. Rigid bodies are constrained to move in the sagittal plane; therefore, no three-dimensional studies of the spine can be carried out. These and other restrictions make it impossible to study with the presently available models some very basic biomechanical problems in orthopaedics. A few examples of such problems are the three-dimensional evaluation of the classical problems: the pilot ejection and the whiplash, the complicated automobile accidents, athletic injuries, the effects of clinical instability, vertebral fusion, description of scoliosis, and the evaluation of the relative efficiency of the various surgical and nonsurgical correcting techniques of the spine. The general mathematical model presented here has the capabilities of studying these and similar problems. THE MODEL
General description
The basic procedure for constructing any mathematical model of a musculoskeletal system is to formulate a set of governing equations of motion representing the mechanical behavior of the model to external loads and displacements. A general method is presented here which makes it rather simple to write the set of governing equations of motion for a discrete-parameter type mathematical model of a musculoskeletal structure. Each rigid body has all the 6 degrees-of-freedom (allowing general three-dimensional motion) and each connecting element (spring and dashpot), all the twenty-one coefficients of stiffness/damping. The model is basically restricted to linear geometry and material properties, but these restrictions can be relaxed to a great extent by special computer techniques. This work is based upon the theories developed by Panjabi (197 1). The spine structure to be studied is idealized to a discrete-parameter type mathematical model. All the anatomic structures are idealized by the three mathematical ele-
M. PANJABI
ments: rigid bodies, massless springs, or massless viscous dashpots. The features of the model and its elements are described below. 1. The model may consist of any number of bodies, springs, and dashpots connected to each other in any manner required to form an anatomically true model. 2. The bodies are three-dimensional rigid bodies and each may have all of the possible 6 degrees-of-freedom of motion. They may be of any size and shape and positioned and orientated in space in any manner to simulate truly the anatomic features of the bony elements they represent. No symmetry of any kind is required. 3. The springs are massless and each may have all of the possible twenty-one (six main and fifteen coupling) spring stiffness coefficients. However, it may suffice to work with a much smaller number of coefficients due to symmetry of some kind. The end conditions may be any combination of fixed, hinged, and ball joints. The springs may be positioned and orientated in any manner to simulate true anatomic features. 4. The dashpots, similar to springs, are massless and each may have all of the twentyone damping coefficients. Again, due to symmetry, it may suffice to work with a much smaller number of coefficients. End conditions may be combinations of fixed, hinged, and ball joints. Position and orientation can be tailored to simulate true anatomic features. 5. The dashpots and springs may be connected in any combination, for example: in parallel to form Kelvin viscoelastic solids or in parallel cum series combination to form the 3-parameter solids. 6. The system may have any initial configuration. 7. External loads, static and dynamic, consisting of force and torque vectors may be applied singularly or simultaneously to any number of the bodies. 8. Matrix methods are used throughout which make the method versatile and suitable for the computer.
MATHEMATICAL
MODEL
OF THE HUMAN
Every method has limitations, and it is important to know these so that one may not only restrict its applications to suitable problems, but also may find methods to relax these limitations. The underlying assumptions of the method are 1. Displacements are small. It can, nevertheless, be relaxed to a great extent, whereever required, by dividing a large displacement into small displacement-steps. The resulting configuration of the system, after a small step, becomes the initial configuration of the system for the next step. 2. Bodies are rigid. 3. Springs and dashpots are massless and linear. Again, this assumption of linear characteristics can be relaxed by adopting small displacement-steps technique. After every step, not only the geometry is altered, but also the stiffness and the damping properties are changed to suit the nonlinear characteristics. Defining
the structure
SPINE
STRUCTURE
613
I
Fig. 1. This shows the method of representing a given musculoskeletal structure by a discrete-parameter type mathematical model.* Rigid bodies are numbered in sequence from 1 upwards and each has local body coordinate system (~yz)~, (~yz)~,. . Connecting elements (springs and dashpots) are also sequentially numbered and have their local coordinate systems (abc),, c&c),... Orientations of the connecting elements are such that the u-axes always point towards the rigid bodies with the higher sequence numbers. Finally there is the global coordinate system (Imn) which defines all the other coordinate systems.
A general model consisting of bodies, springs, and dashpots is shown in Fig. 1. Three types of Cartesian right-hand coordinate systems are used to define the structure: global system, local body systems, and local spring/dashpot systems. The global system, called (fmn), is conveniently orientated and direction towards the body with the higher fixed to the ground. It defines the centers of number. The other two axes, for the sake of gravity (c.g.) of all the bodies and the orientaminimizing the number of coupling coeffictions of the other coordinate systems. The ients, may be orientated in the directions of bodies are numbered from one upwards, with the principal elastic/damping axes of the ground being zero. Local body coordinate spring/dashpot respectively. system (x.v.& is fixed in space at the c.g. of Figure 2 shows the ith body isolated from the ith body and orientated along its principal the rest of the system. Center of gravity of the axes of inertia. (This is required only for the ith body is defined by position vector oi in the study of the dynamic behavior of the system.) global coordinate system (Imn). The points of Springs and dashpots are also sequentially attachment of springs/dashpots are defined numbered. The kth springldashpot has the in local body coordinate system (XyZ)i. A /ith local spring/dashpot coordinate system position vector cljlcdefines the position of the called (crb& It is orientated so that the CZ~- end of the kth spring/dashpot attached to the axis is parallel to the line joining the two ends ith body. The external loads applied to the of the kth spring/dashpot and with its positive body are also defined by similar position vec..*Bony elements are represented by rigid bodies while connecting elements like ligaments. muscles, articulating facets and joints are represented by springs and dashpots.
674
MANOHAR
M. PANJABI
load, di and qj, for each body and in the respective local inertia systems. Thus d=
[d, j dz j-j
di i -/
d,Jr,and
4=[91~q2j-~9i~--i9nlT. The subvectors, components
in turn, are defined by their
[rxr~z+&hl iT, 9i = LfxfiP,P~P~liT di =
where r, 4, f, and p are respectively translation, rotation, force, and moment vectors. The system matrices M, K, and C and the vector q are now treated. Fig. 2. This shows an arbitrary body of the musculoskeletal system, called the ith body, along with one connecting element (spring/dashpot) and one external load, isolated from the system. Local body coordinate system (xy& is fixed in space at the center of gravity i and orientated along the directions of principal inertia axes of the body. Positions of the connecting elements and the external loads are defined by radius vectors Uik and e,K respectively where k is the sequence number of the connecting element.
1. System inertia matrix
The local inertia coordinate systems are orientated along the principal directions of inertia of the respective bodies. This leads to a diagonal system inertia matrix. Thus
(2) tars. eik is the position vector for the point of application of the kth external load vector on the ith body. where System governing
submatrices
for each of the bodies
equations of motion
The procedure described here is given without the underlying arguments. For details, one may consult the original work Panjabi (1971). The system governing equations of motion for any spine structure in the set of local inertia coordinate systems (XYZ)~to (xyz), in the matrix form is tid+Cd’+Kd=q
(1)
where d and q are 6n X 1 system displacement and load vectors respectively. M, C, and K are 6n X 6n system inertia, damping and stiffness matrices respectively. Letter n stands for the number of bodies in the system. The system vectors d and q are defined by a set of subvectors of displacement and
Here, mi stands for the mass and I,, Z,, and stand for the principal inertias of the ith body. Zzi
2. System stifness matrix System stiffness matrix matrix equation
is given by the
K =
(4)
MATHEMATICAL
where submatrices Kij
MODEL
OF THE HUMAN
SPINE
STRUCTURE
675
Kij (6 X 1) are given by =
x k
(5)
K,k.
Each individual submatrix Kz belongs to bodies i and j and the spring k. The rule for summation in equation (5) is If i # j; summation is taken over all the springs connecting the body i with j, and (6)
d*
If i = j; summation is taken over all the springs connected to the body i Further, submatrices Kzare defined in terms of the basic matrices (6 X 1) by matrix equation K”, = ,rJ.“?‘@.r~krS(E.~‘QJJ.~ I 1 (7) 1.1 2.1 ., .r . Here some explanation of the function of individual matrices on the right hand side of equation (7) is justified. We will apply, for illustration’s sake, this equation to a single body attached to the kth single spring, as shown in Fig. 3. Its application in multibody systems is similar but more complex. Thus, equation (7) becomes K” =
~k’l’@T~kT~k~k@
uk_
(8)
If d,,,,, is the displacement produced at the the center of gravity of the body in the xyz coordinate system by the load qtruz) applied to the body at the center of gravity in the same coordinate system, then q(.rwz,= Kkd,,w,,.
(9)
Thus, four basic matrices Uk, Cp,q’” and Sk, of equation (8) represent transformation of vector dClyz)to vector qtS.+ Their individual functions are described below. Uk is the position matrix. It transforms vector 4,,,, to &.vz). The latter is the body displacement at the point of attachment. 4r is the rotation matrix. It transforms vecThe latter is the attachment tor &,,, to dfElmnr point displacement in the global system (Imn). qk is another rotation matrix. It transforms
Fig. 3. Displacement of an isolated single body attached to an arbitrary kth spring is analyzed. Translation and rotation of the center of gravity represented by displacement vector d,,,, produces displacement dk of the spring attachment point. Reaction of the spring is to apply force and moment vectors represented by load vector q* to the body at the attachment point. qczuz,is the load vector at the center of gravity which is equivalent to q” at the spring attachment point.
vector dfT,,, to d&bc,k.The latter is the attachment point displacement in the local spring coordinate system (abc)k. Sk is the spring stiffness matrix. It transforms the vector c&,+ to c&~)~. The latter is the load vector acting at the attachment point on the Spring in the ( UbC) k SyStem. UkTaTqkT is the transpose of ‘Pk@Uk. It transforms q::lhcjk to qtsuxjr the vector on the left hand side of equation (9). Returning back to multibody systems, contents of matrices U”, 0, and qk for both indices i and j are, omitting for convenience indices i and j:
MANOHAR M. PANJABI
676
a=
cosxl cos yl cos xm cos ym cos xn - cos yn _--_
1
cos 21 I cos zm 1 cos zn 1
ZERO
Tk=
cos al
cos
(cosxm
cos cl cos cm ------------_ ZERO
-1
k
ZERO
(15)
S$=Sk
(16) 3. System
+1
ZERO
ZERO
-1
-1
+1
and I 7
k
(18)
I cos an cos bn cos cn1 .
damping matrix
The method of setting up system damping matrix C is the same as that for the system stiffness matrix K. Therefore, equations (419) and Table 1 give the system damping matrix, if word “spring” is exchanged for “displacement” for “velocity” “dashpot”, and “stiffness’ by “damping”. 4. System load vector Let q(lmnl be the system external vector defined in the global system. Then 4(Z?n?l, = c41 i -i
-s11s12 s22
s13
s14
s15
823
St4
825
333
334
s35
s44
s45
SYMMETRIC
s55
4i i - i s;ll~,,,
load
(20)
where qi is the load vector acting on the bodies 1 to n respectively. There may be more than one load vector on a single body. Thus,
The local spring stiffness matrix Sk is Sk =
(12)
The elements sll, s12, - - - x66are the stiffness coefficients along the 6 degrees-of-freedom and their coupled motions. Table 1 gives the information regarding the matrix elements of the basic matrices of equations (10, 11, 12, 17, 18 and 19). Knowing these elements, we can construct the system stiffness matrix K as defined by equation (4).
Sk51 = [‘sq’ 13
where -1
, and
cos 02’L___----__-_--_i cos al cosam I cos bl cos bm I cos cl cos cm
But matrices Si for the kth spring connecting bodies i and j, are different for different combinations of iand j. Here i < j. These are:
I’=
cos zm
-7,
cos an4 cos bnj
I
1
cos ym
, cos x12
cos am cos bm
bl
(11)
---t-------, cos xl
k 9i
(19)
The
=
C k
4ik-
(21)
point of application of qi”, the kth load vector on the ith body, is defined in the local body coordinate system (XyZ)i by the position vector eik.
MATHEMATICAL
MODEL
OF THE HUMAN
SPINE
STRUCTURE
611
Table 1. This table gives information about the elements forming the basic matrices which in turn form the system stiffness matrix K. Equation
How to obtain
Matrix elements
( 10)
lli-, uy, 4
(11)
cos xl, cos yl
(12)
Cos al. cos am,
(16)
s,,,s,~,.. . . . . . . . . . . . . SIX
(19)
1
.. . ..
These are the components of the position vector from the center of gravity to the attachment point in the xyz coordinate system. These are the cosines of the angles between axes x and 1,y and 1,. . . , i and n respectively. These are the cosines of angles between axes a and 1, a and m, . . . . . , c and n respectively. These are the 21 stiffness coefficients representing the stiffness properties of a spring at its j-end in the spring coordinate system. Length of the spring
, cos zn , cos CII
Therefore, transformation of 9(1mn,to the local body coordinate systems xyz and to the respective centers of gravity are given by matrix equation
system behavior under a given set of loading conditions. For a static load vector 9, the solution will be the elastostatic displacement of the bodies of the system, if the effect of damping matrix C is neglected. On the other (22) hand, if transient behavior of the system is of where interest, damping matrix C is retained and a constant load 9 is suddenly applied. In a @= certain study 9 may be a function of time; then the solution will be the dynamic behavior L _ -)- L g?iJ’r-, (23) of the system. In all such cases, solutions of the equation (1) will be in the form of system ZERO L--+_ displacement vector d. Once the system displacement is obtained, the system velocity and and acceleration can be easily computed. r Knowing the displacement and velocity vectors, loads or stresses at the attachment (24) - - - -$y_: __ points of springs and dashpots can be comZERO ;__‘iL_ _ puted. We shall now obtain these relations. Displace the c.g.‘s of the two bodies i and j as shown in Fig. 4 by displacement vectors The submatrices 4ii are given by equation (11) d (syz). and 4zj,z,T Let 9kho and 9tcnbr) be the and submatrices Eii are given by internal loads developed ok the kth sgring at the attachment points. Then we have Eii = 1
_*__
1
_q=o
L
0 0 0
0 - ezi
e,,
0
- egi
e
L
1
.ri
9h’= Q”d
1 ZERO eyi 1 -esi
O
l
0
0
0
(261
where 1
35) I
System behavior
Solution of the system governing equations of motion, i.e. equation (l), represents the
(29)
MANOHAR
Fig. 4. An arbitrary spring k connecting two arbitrary bodies i and j is shown. Unequal displacement of the two attachment points produces forces and moments, represented by load vectors qik and qjk, on the two ends of the spring.
The submatrices
Qi to Qi are given by
Qt = St’Pk@jUjk.
(30)
The basic matrices of equation (30) have been discussed earlier and are given by equations (10-19). Same method is applied to determine the loads at the dpshpot attachment points by replacing d by d and Sk by Dk in the equations (26-30). THE SPINE STRUCTURE
The spine is part of the musculoskeletal system. In a true dynamic analysis, the spini: cannot be isolated and studied separately from the rest of the skeletal system. A good dynamic analysis would at least include the spine structure. Figure 5 shows the spine structure and its various coordinate systems.
M. PANJABI
Fig. 5. The spine structure and three types of coordinate systems required to completely detine its position and motion. Global coordinate system imn is tixed to the ground. Local body coordinate system (XyZ)iis fixed at the center of gravity of rigid body i and orientated in the directions of the principal inertia axes. Local connecting element coordinate system (ab~)~ is orientated along the principal elastic/damping axes of the kth connectingelement. The ak axis points towards the rigid body i+ 1.
For the sake of clarity muscles are not shown, but they may be included in the analysis. Bony elements of the structure (pelvis, vertebrae, head, and sternum) are represented by rigid bodies. Connecting elements joining the bony elements consist of ligaments, articulating joints, cartilage, ribs, and muscles. In general, connecting tissues have viscoelastic properties and are most suitably represented by either a Kelvin type solid or the three-parameter solid. In the latter case, junction between the Kelvin solid and the series spring may be considered as a body of zero dimensions. For the details see Appendix. Single ligaments may be represented by uniaxial three-parameter solids which resist only when pulled. On the other hand, complex behavior of the vertebral motion segments may be represented by three-dimensional solids, Articulating joints “three-parameter” may best be approximated by fictitious springs with zero stiffness in the plane of
MATHEMATICAL
MODEL
OF THE HUMAN
contact and large stiffness in all other directions. As a first approximation passive muscles may be represented by Kelvin or threeparameter uniaxial solids combined with external (muscle) forces acting on the bodies at the points of attachment. The effect of the rest of the musculoskeletal system and the outside constraints (safety harness, seat, etc.) may be neglected or taken into account by including these within the system of analysis. DISCUSSION
No illustrative examples have been given here because a good example will need too much space, while a simple example will not do justice to the generality of this method. The method, however, has been applied in two cases in the field of mechanical engineering. These are the vibration study of a gramophone, an elastodynamic system, panjabi (197 I), and the study of dynamic stresses in the bali bearings of a washing machine, a viscoelastic dynamic system, Eriksson and Jansson (1972). I am confident that the same general method could be used successfully to study the mechanical behavior of complex musculoskeletal systems. Because of its complex three-dimensional geometry and its limited motion, the spine structure is a good example of such a system. The model presented here is general. It can be used with any data: calculated strengthof-material formulae, incomplete experimental data available at present, or the three-dimensional, including coupling effects, experimental data which we hope will be available soon. Of course, the results obtained from the model will reflect the accuracy of the data. We are working in our laboratory to obtain more complete data of the spine segment behavior. We can simulate complex behavior where
not only the structure is three-dimensional but also the external loads. Force and torque may
be applied
to any point
on any body
in
SPINE
STRUCTURE
679
the system, and we can study the static. transient, and dynamic behavior. Work is in progress on some of the clinical projects, mentioned earlier in the introduction, concerning the behavior of the spine structure. Results will be reported in future publications. REFERENCES Eriksson, S. and Jansson, L. (1972) Forces in freely suspended washing machine. (In Swedish) MS. Thesis. Division of Machine Elements, Chalmers University of Technology, Goteborg, Sweden. Fliigge, W. (1967) Viscoelustirity.BIasdell, Waltham, Mass. Frisen M., M%gi M., Sonnerup, L. and Viidik, A. ( 1969) Rheological analysis of soft-collagenous tissue. J. Biomechunics 2. 13. Hess, J. L. and Lombard, C. F. (1958) Theoretical investigations of dynamic response of man to high vertical accelerations. J. Aoiation Med. 29,66-75. Latham. F. (19.57) A study in body ballistics: Seat ejection. Proc. R. SW. B-147. Liu, Y. K. and Murray. J. D. (1966) A theoretical study of the effect of impulse on the human torso. A.S.M.E. Biomechanics Symposium. pp. 167-186. Ome, D. and Liu, Y. K. ( 1969) A mathematical model of spinal response to impact. Thesis 02370, Department of Engineering Mechanics, Univ. of Michigan. Paniabi, M. M. (197 1) Theoretical treatment of vibrations in”single and multiple body suspension systems based on matrix methods. Thesis, Division of Machine Ele-
ments, Chalmers University of Technology, Goteborg. Sweden. Payne, P. R. (1962) The dynamics of human restraint systems. In Impact Acceleration Stress. National Academy of Sciences, National Research Council. Ptiblication 977, pp. 195-258. Terry, C. T. and Roberts, V. L. (1968) A viscoelastic model of the human spine subjected to i a: accelerations. J. Biomechanics 1, 16 1. Toth. R. (I 967) Multiple degree-of-freedom nonlinear spinal model. Presented 19th Annual Conference on Engineering in Medicine rmd Biolo,qy, San Francisco. California. NOMENCLATURE MutriceJ C damping matrix n local dashpot matrix E position matrix for external load I’ 180”rotation matrix K stiffness matrix L force transformation matrix M inertia matrix Q spring-/dashpot-end load matrix s local spring stiffness matrix u displacement transformation matrix cp rotation matrix: xvi to lmir \Ir rotation matrix: lmrr to trhc
680
MANOHAR
M. PANJABI rz
Vectors d
f2
displacement vector position vector for external load ; force vector 4 load vector r translation vector IA position vector moment vector rotation vector f
2
Scalars c damping coefficients I principal moment of inertia of the body 1 initial length of spring m mass of the body n number of bodies in the system s stiffness coefficients i
Indices a, b, c Lj k
axes of the spring coordinate system two arbitrary bodies, i and j arbitrary spring, arbitrary load 1, m, n axes of the global coordinate system x, Y, z axes of the local body coordinate system
Special
(b) (a) Solid; uniaxial response. Fig. 6. (a) Three-parameter Solid as incorporated in a 6(b) \-, “Three-parameter” degree-of-freedom System. A fictitious body h of zero dimensions is introduced to represent the junction of the Kelvin Solid and the Series Spring.
Symbols
- derivative withrespect to time T transpose. APPENDIX Three-parameter
ligament
model
The ligaments may be represented either by Kelvin Solid or by a more complex model, the three-parameter solid. The latter model simulates closely the viscoelastic behaviour of ligaments in the stationery phase i.e. when all deformation is recoverable as demonstrated by Frisen et al. (1969). We can incorporate this ligament model by considering the junction between the Kelvin solid and the series spring as a fictitious rigid body of no dimensions. In this manner it is possible to incorporate into the mathematical model a true three-dimensional “3-paramesolid for representing ligamenteous ter” viscoelastic structure. It is known that the behavior of the uniaxial 3-parameter solid as shown in Fig. 6a is given by the following differential equation c3ji+
(k, + k2)fz = k2cJiZ + k,kpr2.
(31)
Symbols are as in Fig. 6a and formulation is according to Fliigge (1967). To incorporate this model in a threedimensional multi-body system we can build a model as in Fig. 6b. Bodies i and j are connected by the threedimensional “3-parameter” viscoelastic solid consisting of springs 1 and 2, dashpot 3 and the fictitious body h. Here i < h
as body number 1 and the free end of the 3-parameter solid as number 2. The other end of the 3-parameter solid is fixed to ground numbered 0. Then equation (1) becomes Cd+Kd=q
(32)
where
d=
j dJT
[d,
and
q=
j q2]T
[q,
With the use of equations (4- 19) for the simple uniaxial case of Fig. 6a, contents of equation (32) reduce to C=[; d=
;I.“=[“;& [r,
rllT
and
-21, [O f2]‘.
Which when expanded gives (k,+k2)r,-kpre+csr, - k2r, + kzr2
=0
(33)
=h.
(34)
It can be easily seen that the two equations (33 and 34) are the same as the equation (3 1) if r, is eliminated. Thus by introducing an additional displacement for the junction of the Kelvin solid (rl in the uniaxial case) and the series spring we can always formulate system goveming equations in the standard form of equation (1). The price paid for this convenience is the increment in size of the system stiffness and damping matrices.