THE U N I V E R S I T Y OF M I C H I G A N COLLEGE OF ENGINEERING Department of Engineering Mechanics Progress Report THE HUMAN BODY UNDER TIME-DEPENDENT BOUNDARY CONDITIONS (Period: September 1, 1966 - March 15, 1968) Y. King Liu ORA Project 01279 supported by: U. S. PUBLIC HEALTH SERVICE DEPARTMENT OF HEALTH, EDUCATION, AND WELFARE NATIONAL INSTITUTES OF HEALTH GRANT NO. UI-00025-02 BETHESDA, MARYLAND administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR March 1968

PREFACE This project was supported by The National Center for Urban and Industrial Health of the Public Health Service, Department of Health, Education, and Welfare under Research Grant UI-00025-02 (formerly AC-00291-01). The author wishes to extend his gratitude to Professor J. D. Murray of the Department of Mathematics, New York University for many helpful discussions while he was a member of the Department of Engineering Mechanics at The University of Michigan. Appreciation is also extended to Professor F. G. Evans of the Department of Anatomy for his consulting on aspects of anatomy. Assistance was ably provided by Mr. David Orne, a Research Assistant in the Department of Engineering Mechanics, on the material in Chapters I, TI, and IV. Mr. Ali E. Engin played a similar role for the material on Head Injury (Chapter III) while employed as a Research Assistant by the Highway Safety Research Institute of The University of Michigan under contract No, 43-67-1136, sponsored by the National Institute for Neurological Diseases and Blindness. Programming and other assistance were also provided by Mr. Chung-Yi Wang, Mr. James E. Tuttle, and Mr. James D. Harris. iii

TABLE OF CONTENTS Page LIST OF FIGURES vi INTRODUCTION 1 Chapter I. THE EFFECTS OF VERTICAL IMPULSE ON THE HUMAN TORSO A. The Problem 3 B. Previous Studies 3 C. Lumped-Parameter Models 3 D. Experimental Studies, Clinical Observations and Biormechanical Data 5 E. Elastic Continuum Models 14 F'. Viscoelastic Continuum Models 23 G. Love-Rayleigh Model 27 H. A Note on the Inversion of the Transformed Solution 28 I. The Tapered Elastic Model 31 J. Nonlinear Elastic Model 32 K. Effects of Nonlongitudinal Loading 34 References to Chapter I 35 II. THE EFFECTS OF HORIZONTAL IMPULSE ON THE HUMAN TORSO 37 A. The Problem 37 B. Previous Studies 37 C. Lumped-Parameter Models 41 D. The Elastic Beam Model 42 References to Chapter II 47 III. A MODEL FOR HEAD INJURY-THE AXISYMMETRIC RESPONSE OF A FLUIDFILLED SPHERICAL SHELL TO A LOCAL RADIAL IMPULSE 48 A. The Problem and Previous Studies 51 B. The In Vacuo Case 53 C. The Fluid-Solid Interaction 61 References to Chapter III 66 iv

TABLE OF CONTENTS (Concluded) Page IV. NUMERICAL AND ANALOG SIMULATION IN BIOMECHANICS 68 A. Introduction and Previous Studies 68 B. Numerical Simulation of the Spine under Vertical Inpulse 69 References to Chapter IV 76

LIST OF FIGURES Figure Page 1-1. The vertical acceleration of the human body. 4 1-2. The lumped-parameter model of vertical acceleration. 4 1-3. Transmissibility of longitudinal vertical vibration. 6 1-4. Load-deflection curve of the intact lumbar spine. 7 1-5. Typical experimental strain and acceleration record. 8 1-6. Mechanical analog of test seat and test subject in a drop-tower test. 10 1-7. Incidence, distribution and level of vertebral fractures in pilot ejection. 10 1-8. Load-deflection results for intervertebral disc. 12 1-9. Failure load of vertebrae. 13 1-10. Typical creep curve for disc. 13 1-11. Creep response of Kelvin model. 15 1-12. Uniform elastic rod model of the spine. 16 1-13. Continuum model of torso under vertical acceleration. 16 1-14. Stress-time history for X =.1 and.2. 20 1-15. Stress-time history for A =.5. 21 1-16. Stress-stime history for X = 1. 22 2-1(a), Lumped-parameter model of the whiplash problem. 38 2-1(b). Typical experimental data. 39 2-1(c). Typical computed results. 40 vi

LISTS OF FIGURES (Concluded) Figure Page 2-2. Lumped-parameter model of the seat-belt syndrome. 43 2-3. Continuum model of the whiplash problem. 44 3-1. Head injury model. 52 3-2. Momentum considerations for the shell. 57 4-1. Sketch of vertebral column curvature. 71 4-2. Lumped-parameter deformation models for the verterbral column (many-degrees-of-freedom). 73 vii

INTRODUCTION Any external movement of the skin constitutes a mechanical input, either transient or steady-state, into the human body. The basic premise used is that the human body will respond to the input according to the laws of mechanics whether or not the human sub ject has preknowledge of the impeding action. Some topical examples are: automobile accidents, ejection from disabled aircraft, manned launchings, and recoveries of space capsules on water or land. The research results reported here are "theoretical?" in nature. The mathematical models developed are such that: (1) most of the essential features in experimental results are present and accounted for, and (2) no interesting and/ or valuable qualitative and quantitative results have been overlooked. In short, the research is aimed at explaining the experimental data as well ae suggest additional experiments to be performed. No real problem in biomechanics can or even should be treated in all its complexity. Usually, the most that can be hoped for is -the identification of the more important variables in the problem and their effects on the results of interest. The process of developing models to describe given physical systems, such as the human body, which meet this requirement is not easy. The development of meaningful experiments must necessarily be guided by some theory, no matter how simple. On the other hand, the construction of a meaningful theory requires the availability of experimental data and observations. If data are not found or provided, there is always the danger that theory will be built up from an abstract philosophical viewpoint based more on how the human body could or should respond rather than how it does. In this initial effort the aim has been to develop those particular problems which have the most topical relevance and are simple analytically. It was felt that the most influential parameters and their variations can be reasonably understood with such a beginning~ The complex real situations in the biodynamic response of the human body have also been treated. Numerical methods, e.g., multi-degree-of-freedom lumped parameter models, finite difference and/or finite element analysis, have been initiated to tackle problems, which are intractable analytically. While these techniques are discussed toward the end of this report (see Chapter IV), the main emphasis has been placed on'"severely"' idealized models, which are justifiable on two counts: (i) they are guides to the most important elements in the'real" problem, and (2) the mathematical models treated are far superior to those which have entered the literature in their ability to correlate with experimental datao 1

A number of biodynamic problems were considered. From among these, three significant results and approaches deserve special mention: (a) A continuum model for the vertical ejection problem was proposed and solved exactly. The model consists of an elastic rod capped by a mass at one end and given an impulse at the other, see Figure 1-13. The elastic rod represents the spinal column and the end mass, either the head or the head, upper limbs and thoracic cage combination. The results were excellent qualitatively in that it was possible to keep track of the progress of the wave front and its effects at all locations and times. Quantitatively, however, the results were on the low side. That this would be the case was anticipated, since the assumption of a linear material (constitutive relationship) was a mathematical convenience and does not correspond to the known load-deformation characteristics of the human spine. However, the locations of high stress and/or acceleration became clear, while effects of such parameters as the time of wave travel, the head to spine mass ratio, etc. became clearly delineated. The effects of viscoelastic materials are found from the elastic solution through a correspondence principle. (b) The problem of head injury received considerable attention. The model consists of a closed elastic spherical shell filled with fluid subjected to a blow, Fig. 3-1. As a first approximation, the shell is considered to be thin and elastic, the fluid is inviscid and irrotational while the blow is manifested by a sudden initial velocity input. Even with these simplifying assumptions, the problem is quite complex and it was decided to break the problem down to its constituent components: (1) a rigid closed spherical shell with irrotational inviscid fluid as the wave guide and (2) the elastic closed spherical shell as a solid wave guide and finally (3) the effect of the fluid-solid interaction of the model is considered. (c) The initiation of "numerical experiments" to ascertain the relative importance of the various parameters in a given biodynamical problem. The essential concept here is to collect as much of the available biomechanical data and clinical observations as possible, then based on anatomical guidance construct a rather detailed mathematical model, which because of its complexity can only be treated numerically. This computer program has merit by itself and also as a means of checking the validity of various "severely" idealized models. 2

CHAPTER I THE EFFECTS OF VERTICAL IMPULSE ON THE HUMAN TORSO A. THE PROBLEM Consider a "pilot" to be seated as shown in Fig. 1-1, At time t = 0, he is propelled from the aircraft by an explosive charge, F(t), acting upon the bottom of the ejection seat. This action can and frequently does lead to, among other effects, damage of the spinal column. We would like to examine the biodynamics of the problem. Bo PREVIOUS STUDIES In 1961, Goldman and Von Gierke1 summarized the data on the effects of shock and vibration on man. In 1964, Von Gierke reviewed the current researches in the biodynamic response of the human bodyo These two reviews, with minor exceptions, are quite complete. Between 1964 and the present date of this report, two papers have appeared, which supplement the previous two survey articles. The first, due to Roberts, Stech, and Terry, reviewed the variety of mathematical models used4in connection with biomechanics up to 1965. The second by the writer and Murray, gave a summary of lumped-parameter models and then examined the use of continuum models in the dynamics of body ballisticso The present chapter will highlight some of the contents of4, which appears as Appendix A, and also report on the progress made since its appearance. Co LUMPED-PARAMETER MODELS It might appear far-fetched to model the given problem by a single-degreeof-freedom lumped-parameter system. Its justification was most succinctly stated by Payne5: "For short duration accelerations, however, the soft, low frequency parts of the body (such as the viscera) do not deflect far enough to load up the spine,.o o The lumped-parameter model is shown in Fig. 1-2. Mathematically, the model can be represented by the ordinary differential equation: mx + f(x) = -mg(t), (1-1)

SPINE SEAT F(t ) Fig. 1-l. The vertical acceleration of the human body. m f(x) g(t) CHAIRFig. 1-2. The lumped-parameter model of vertical acceleration.

where f(x) is the spring force, g(t) is the acceleration pulse applied to the base, and x is the relative displacement between the base and the articulated mass, m. The problem is linear or nonlinear according as f(x) is linear or nonlinear. For the case of a rectangular pulse, 4 showed that in the linear case the dynamic load factor is two while it can be larger than two for the nonlinear case. The dynamic load factor is defined as the ratio of the maximum dynamic spring force to the static spring force. Lumped-parameter models have also been used for describing the response of the human body to steady-state vibrations. Coermann et al., 6 have used a seven-mass system for predicting the response of the human body in the sitting or standing position on a platform subject to vertical vibrations at low frequencies. DieckmannT has obtained the transmissibility of longitudinal vibrations from a shake table to various parts of the body of a seated human subject as a function of frequency, as indicated in Fig. 1-3 Observe the resonances of the head/shoulder and the shoulder/table motions at fairly low frequencies. ItJ was stated in. Refo 6 that with respect to steady-state vibrations the body can be described as a lumped-parameter system at frequencies less than 100 cops but at higher frequencies up to 100,000 cps the body behaves more as a complex distributed parameter system, the type of wave propagation being strongly influenced by boundaries and geometrical configurations. Do EXPERIMENTAL STUDIES, CLINICAL OBSERVATIONS AND BIOMECHANICAL DATA 8 Lissner and Evans measured the response of the human cadaver to static loading and to caudo-cephalad accelerations. Typical of their results are: (1) The load-deflection curve of the intact lumbar spine as shown in Fig. 1-4. Its "'hardening" characteristics is quite evident, i.e., the load is a monotonically increasing function of the deformation of the spine. (2) A typical experimental record of the dynamic response of a cadaver, when strain-gage and/or accelerometer measurements are made at selected locations of the skeleton, is shown in Fig. 1-5. The experimental results when compared to the analytic response of a simple linear spring-mass system show good correlation for motion measured at the crest of the ilium and at the sternum. The experimental response of the head, however, indicated much higher accelerations than those predicted by the simple model. Also, strain-gage measurements made on selected vertebrae show that the locations of maximum strain depend upon the mode of bending of the vertebral column, bending being more severe in the cervical and thoracic regions than in the lumbar region. It was concluded that while gross motions of the body may'be approximated by sauch a simple lumped-parameter model, this model is inadequate for predicting injury. Another one degree-of-freedom, lurnped-parameter model was used by Wit:Lman9 in an attempt to determine the influence of the human body on its supporting

3.5 I\ 3.0 I HEAD/ SHOULDER I 2.5 o 2.0 o20 2.0SHOULDER/ / H/ \ /HEAD/ oIt~ ~/ _i\/ TABLE 0 1 2 3 4 6 8 10 20 30 40 FREQUENCY, CPS Fig. 1-3. Transmissibility of longitudinal vertical vibration

LOAD DEFLECTION TESTS RESULTS FOR INTACT LUMBAR SPINE CADAVER 182 300 250 200 cot lli g~~~ 150 100 00 -J 0 10 20 50 40 50 60 70 80 90 DEFLECTION - INXXIOFig. 1-4. Load-deflection curve of the intact lumbar spine.

LUMBAR 5 STRAIN MILLISEC. 0 40 80 120 160 200 0 0 0 0 1.63 CT c~~~J o 40 so 120 MILLISEC. U)o~~~~~ ~CHAIR ACC. 0 40 80 120 200 MILLISEC. Fig. 1-5. Typical experimental strain and acceleration record.

structure. From data collected in drop-tower tests, the parameter values for a spring-mass-damper system representing the human body were progressively refined (by trial and error) to develop an analytical model capable of predicting the dynamic response of a drop-tower carriage to which human subjects were strapped. Figure 1-6 shows a sketch of the carriage and the simple spring-massdamper model. For the conditions of the drop tests it was found that analytic predictions of seat accelerations corresponded well with experimentally measured seat accelerations if the model of the human body had a natural frequency of 10 cps, a critical damping ratio of 0.30, and a mass equal to the weight of the test subject. Maximum accelerations varied from 12 to 14 gYs, when short- and long-pulse durations of 55 and 120 msec were used. The above study has shown that a lumped-parameter model of the human body can be used to determine the accelerations of a fairly rigid supporting seat assembly subjected to impulsive loadings. In the mechanical modeling of the human system one must have access to the relevant biomechanical data before the modeling process can be meaningful. Frequently, however, such data is simply not available, or at least not in a form suitable for direct use by the engineer or physical scientist interested in developing a mathematical model. With regard to the spine, more attention appears to have been given to the behavior of the cervical and lumbar regions than to the thoracic region. As will be pointed out later, vertebrae in the thoracic region are as likely to fracture during pilot ejection as, say, are the lumbar vertebrae. Although much of the motivation for performing biomechanical tests on the vertebrae and the intervertebral discs has been the need to explain such medical conditions as sciatica and lower back pain, much of the experimental data thus obtained can be interpreted by the engineer for his own applications. In some cases, however, only qualitative properties such as the presence of strain-rate effects can be gleaned from the experimental evidence presented. The remainder of this section presents the results of some of the typical biomechanical tests performed on parts of the spinal column, with representative data being given. 10 Hirsch and Nachemson have summarized the condition of the spine of 55 Swedish aviators who were forced to escape from their aircraft by catapult ejection during the period from 1957 to 1960. Thirteen, or about one fourth, of the subjects were found to have incurred vertebral fractures. The injuries were invariably compression fractures of the vertebral bodies. In five cases, only one vertebra was damaged, whereas eight cases showed multiple fractures. The incidence, distribution, and level of vertebral fractures that were detected are summarized in Fig. 1-7. The middle part of the thoracic spine was subjected to the highest frequency of fracture. 11 Brown and others conducted experiments on the lumbosacral spine (five lowermost vertebrae with attached sacrum) with particular emphasis on the 9

FI2 F2 Y K4,'X 3 F3 Fig. 1-6. Mechanical analog of test seat and test subject in a drop-tower test. Cn m 64 L | [After Hirsch] a 4 3- 3- - Th 2 4 6 8 10 12 3 LEVEL OF INJURED VERTEBRAE Fig. 1-7. Incidence, distribution and level of vertebral fractures in pilot ejection. 10

intevertebral discs. They found that the ultimate static compressive load for the lumbar discs ranges from 1000 to 1300 lb. Figure 1-8 summarizes the findings of axial compressive tests conducted on five specimens taken fresh from human cadavers. Stiffness values for these discs range from about 500 to 8000 lb/in. initially to 12,000 to 20,000 lb/in. after the applied load reached 200 to 400 lb. A typical load-deflection curve is shown in Fig. 1-8. Failure under axial load took place invariably in the vertebral end plates even when well-developed ruptures of the annulus of the discs were present. Failure ranged from imperceptible cracks to more or less complete collapse of the end plate, depending on the condition of the bone and the magnitude of the applied loado For values of the breaking strength of individual vertebrae of a much larger number of individuals, see Fig. 1-9, which is taken from Ruff. The units of failure loads are given in kilograms. It can be seen that age is not a discernible factor in the strength of the vertebrae and that the data reflects surprisingly little spread in the failure load for individual vertebrae. Evans and Lissner described the results of axial load and bending tests on the intact lumbar spine and pelvis. Specimens from fourteen embalmed and six unembalmed adult human male cadavers were used. Under axial load, the embalmed specimens' average maximum load was 882 lb (610-1350), but in the embalmed specimens it was only 544 lb (290-690). However, the embalmed specimens absorbed more energy, 594-in,-lb, compared to the unembalmed specimens which absorbed on the average only 567- in.-lb. A typical load-deflection curve for the lumbar spine is shown in Fig. 1-4. A characteristic bilinearity exists which is quite similar to the bilinearity of curves obtained for intervertebral discs taken alone. In the bending tests conducted, the intact lumbar spine was supported as a simply supported beam with a concentrated load applied (through a steel cylinder) to the midpoint of the "beam." Five specimens were tested in anterior bending and four specimens were tested in lateral bending. In anterior bending little difference was observed between the bending moment at failure for the embalmed and unembalmed specimens, the average being about 700 in.-lb. The bending tests of the lumbosacral spine reveal the "integrated" behavior of all of the constituent elements comprising that part of the spine, thus providing some over-all phenomenological check on the plausibility of any specific mechanical models one might propose. For example, we see that the typical load-deflection curve for the bending tests reflect again the same bilinear characteristic as do the load-deflection curves for individual intervertebral discs, see Fig. 1-4. Most of the mechanical properties of intervertebral discs were obtained under fairly slow rates of strain and thus represent essentially the static properties of the discs. Hirsch,13 however, has tested discs under constant load and measured the resulting deformations (creep response) as a function of time. A typical creep curve for a disc subjected to a 100 kg load is shown in Fig. 1-10. Hirsch made the following interesting observations: 11

Load Vertical 1200 Deflection G-2 G. I000 Anterior G-3 G-4 800 80 IG-2K C G-I / Major Slope = O- 600 L LaG~er R / 18,000 lbs/in 2600 G-3 Lateral- / Deflection /:180 lbs 400 -— J.01 Initial Slope= 2500 lbs/in 200 *1 ~ 50 lbs 0.02. 0 0 0.020 0.060 0.100 Lateral Deflection - Inches Fig. 1-8. Load-deflection results for intervertebral disc.

1000 978 lbs. 6 800 600 400 200.3.4.5.t5 -I INCHES Fig. 1-9. Failure load of vertebrae., i. I _- - Load=loo'.,;' I, Q 0.9,, X, I. e o / 3 4. 5 7 rrime, minn. D/sc Creep Test Fig. 1-10. Typical creep curve for disc. 13

(a) the curves are essentially level after 5 to 10 min of load application, (b) most of the deformation occurs in the first 3 0 seconds. (c) the nature of the curves is independent of the level of load, and (d) the discs recovered completely after the load was removed. Unfortunately, the detailed behavior at early times (such as the presence of an instantaneous elastic response) was not given. Nevertheless, the over-all behavior of the discs can probably be adequately described by a one-dimensional, linear viscoelastic model. For example, the creep response of a Kelvin model, shown in Fig. 1-11, has a striking similarity to the creep response of the intervertebral discs. E. ELASTIC CONTINUUM MODELS A continuum model capable of predicting the motion of the head when the human torso is subjected to longitudinal accelerations has been used by Hess 14 and Lombard. They proposed the uniform linear elastic rod model shown in Fig. 1-12. By varying the wave speed c of the rod they were able to obtain a "close" fit to the measured response of the head of test subjects. The time of wave arrival thus determined averaged 25 msec within a range of from 20 to 30 msec. Thus, a very simple continuum model with an adjustable parameter, time of wave travel, was developed for determining the motion of the head when the vertebral column is subjected to known longitudinal accelerations at its lower end. However, since "matching" of experimental data and analytic results was obtained at only one point on the rod, namely at the head, such a continuum model cannot be justified as a model for predicting injury to the vertebral column, not at least by virture of the above study alone. Carrying the continuum approach further, Liu and Murray proposed and studied a model consisting of a uniform rod capped at one end by a rigid mass, Fig. 1-13. The rigid end-mass represents either the head or the head and part of the upper torso, and the rod represents the spinal column. The longitudinal motion of the rod is governed by the equilibrium and kinematic equations: =-,(1-2) i' = u', (1-3) Xt where subscripts denote partial differentiation, x' is the longitudinal coordinate, ut is the displacement, c' is the strain, a' is the stress, p is the mass density, and t' is the time. The above equations are, of course, independent of material properties. For the case of the elastic rod, the constitutive equation is 14

6: _):;J Creep Compltince- ke/v/n MAfer4/J C/?eep Comp/andce - 3Pararneb'r A/oie' J'(~>J =7~ro n respo n~ +~~30 3/-e Fig. 1-11. Creep response of Kelvin model. 15

4'6/ostrc Rod t pe/vls tIvSod Fig. 1-12. Uniform elastic rod model of the spine. A,Ele, L Fig. 1-13. Continuum model of torso under vertical acceleration. 16

C' = EE', (1-4) where E is the modulus of elasticity of the spine. The combination of (1-2), (1-3), and (1-4) yields the usual one-dimensional wave equation: u ttt = c u' 1 (1- 5) where c is the velocity of sound in the elastic medium and is (E/p)1/2 The boundary and initial conditions of the problem are: AE u' Ix(,t) = -M ut,tt(e,t') (1-6) u'(O, t') = f(t') u'(x',O) = ut, (x',O) = 0, (1-7) where A is the cross-sectional area of the spine. Eqs. (1-6) and (1-7), respectively, state that the "head" times its acceleration is equal to the force at x' = ~ and the displacement of the seat is prescribed by f(t'). Introducing dimensionless variables: u = u'/, x = x, t = ct:!/~ and X = M/pA, (1-8) where X is the "head"-to-stem mass ratio. Substitution of (1-8) into (1-5) to (1-7) gives the dimensionless form of the equation as: u = u, (1-9) tt xx with boundary and initial conditions: x(l,t) -X tt(l,t) x - utt u(O,t) = f(t) 17

u(x,O) = t (X, ) = 0. (1-10) Taking the Laplace transform of (1-9) and (1-10), we get = 2-p XX u(x,p) = u e pu(x,t)dt (1-11) 0 Ux(l.p) = -Xp2 (T(lp) and u(O,p) = i(p) (1-12) Combining (1-11) and (1-12) we get u(x,p) = f(p)[[e- PX+E(p)-p(2-x) ][l+E(p)e 2]1/2 (1-13) where E(p) = (l-Xp)/(l+\p). Application of the binomial, shifting and Faltung theorems yield the exact solution, which is given in Ref. 4. We shall only reproduce the main relevant results here. The nondimensional displacement solution is 00 u(x,t) = f(t-x) + 7 f(t-x-2n) - f(t+x-2n) n=l + ZE 2 (-l)v 2 n) St [f(-2n-x) - f(~-2n+x)](g-t)v le d( t), n=l v=l X (v-l) V0 (1-14) where f(t) = 0 for t < O. The dimensionless stress field at any time is o(x,t) = a'(x,t)/E = -ux(x,t), (1-15) 18

where u(x,t) is given by (1-14). The acceleration is given by utt, which at the "head" is X -l(l,t), see (1-12). Typical stress-time histories with x as a parameter are given in Figs. 1-14 to 1-16. A continuum model answers several questions which heretofore have been unexplained~ It is clear from the experimental and biomechanical data that stress is the cause of injury and not the acceleration directly. These two quantities are related differently for different boundary value problems. The use of the maximum value of acceleration as the criterion of tolerance is, by itself, invalid. Furthermore, the locations of high stress and/or acceleration become clear and the effects of such system parameters as the time of wave travel, "headt"-to-stem mass ratio, etc., became delineated. The results are best exemplified by the dimensionless stress-time plot for a "head"-to-stem ratio of X = 0.5, Fig. 1-15. One unit on the time scale is physically the time for the impulse to travel the length of the rod, Fig. 1-12. Notice that the stress at the mass end is zero prior to the arrival of the wavefront. At the loaded end, the dynamic stress is twice Young's modulus of elasticity just when the reflected wave has returned to its origin. If there were no head mass, i.e., = O0, the waves would be bounced back and forth following the dotted lines. The addition of the "head" mass introduces dispersion into the wave propagation with interesting consequences. Notice that the stress is not necessarily negligible. Furthermore, there exists a tensile stress at the head end near t = 6. Such a phenomenon is difficult to realize in a singledegree-of-freedom, lumped-parameter system. It occurs quite naturally in a continuum model. As far as the author is aware there is little, if any, data on the tensile strength of the spine for either the static or dynamic case. Anatomically, the spine is basically a compression member. It would not be too surprising if its tensile strength is rather low. Only additional biomechanical data, from experiments can confirm or deny such an assertation. 10 Hirsch and Nachemson gave injury statistics of Swedish pilots who were forced to eject from their aircraft, Fig. 1-7. If the results in Fig. 1-15 are representative of the mechanics of pilot ejection, then one should expect predominantly compression failure near the lumbar region: Thl2, L1, L2 and L3. However, the severe and multiple injuries occurred in the thoracic region near Th5 and 6. This can be partially attributed to the inertial effects of the "head" mass consisting of the thoracic cage, upper limbs and the head, The multiplicity of the fracture is probably caused by initial compressive failures and then the subsequent arrival of the tensile stress. The simple continuum model given above do yield some guidelines to such basic questions as where, when and how damage might occur. In some important cases, it is far easier to derive and solve the partial differential equations governing the more realistic distributed-parameter model than to create the equivalent lumped-parameter system and then solve the coupled system of ordinary differential equations. 19

c'(x,t) 3 2 I(Ot) =1:0. 2 \ /\I I ~~~/ F - Sss frX.cr1o, X~~~~~T /_.)tX-0.1 I ["(~ O~~~~~~~~ / 0 2 3 4 5 6 7 8 9 10 Fig. 1-14. Stress-time history for k =.1 and.2.

a-(Xtt) O'(x,t) I _ a\ V //+\/ / XX 0=0.O 0 I 2 3 4 5 7 8 9 1 0 Fig. 1-15. Stress-time history for X =.5.

o~~~~~~~ rooo / k ___ _ ___ z~~~7 - A — -,........._-_.. C:3D -J.. I -\- -_ \ C\ I X-oI A0.2 / LL\ / COO C) — CD I~~~.00 1.00 2,00 3.00 4.00 5.00 6.00 7.00 8.00 TIMEoT THEOFETIOFL STHESS-TIME HISTOHY Fig. 1-16. Stress-time history for. = 1.

Admittedly, many of the system parameters such as damping, initial crookedness of the spine, nonuniform mass distribution and section properties, nonlinear constitutive equations for material behavior, etc., have not been includedo We shall like to consider some of these effects singularly and in groups. F. VISCOELASTIC CONTINUUM MODELS It was pointed out in the previous section that the equilibrium and kinematic equations, i.e., (1-2) and (1-3), are independent of material properties. For a linear viscoelastic rod, the constitutive equation is generally N dnc M dmY n-O n dt n m m d m (1-16) dtvn dt' where the coefficient p0 is usually taken equal to unity, The transformed solution of (1-16), assuming quiescent initial conditions for both a' and c', is: N n M n O Pns ) qmsm) where (x, s) = re a(x, t')dtv 0 Let M N E*(s) - Z m/n Pnsn (1-17) m=O then = = E*(s) E? (1-18) Laplace transform equations (1-2), (1-3), and combining with (1-18), we get the transformed equation for viscoelastic rods [E"()] 9 9 s2a. =(1-19) 23

The solution of (1-19) depends on the loading, initial, and boundary conditions. If these conditions are identical for both the elastic and viscoelastic problems then their transforms will also be the same. Hence, if E in the transformed solution of the elastic problem is replaced by E*(s) then the result will be the transform of the viscoelastic problem. This is the so called "correspondence principle" 15 as applied to the present one-dimensional wave-propagation problem. In the present case, we have solved the elastic problem in the previous section, hence we can apply the correspondence principle to get the transformed solution to the viscoelastic problem. The greatest difficulties arise in attempting to invert the transformed solution into the real time domain. Of course, one can always formulate and solve the viscoelastic problem directly. Two examples will be given below to illustrate the power of the technique. Example 1: Kelvin Material, Finite Rod A viscoelastic rod made of a Kelvin material of length ~ is subjected to a displacement impulse at one end with the other end free. Use the correspondence principle to obtain its transformed solutions. The constitutive equation for a Kelvin material is: at = Ec' + Ms' (1-20a) (t= (E+ Ais) i' = E*(s) ET (1-20b) The transformed solution of the corresponding elastic problem, ioe., (1-13) with X = O, is: u(x,p) = f(p)[[e Px + e P(2-x) ][1 + ep]-] (1-21) where the transformed variable p is with respect to dimensionless time, to From the definition of the dimensionless parameters, (1-8), the constitutive equation, (1-20) and the scaling theorem, i.e., L[g(at')] = g () (122) a a we get, after a little manipulation with a = ~/c, 24

U(x,S) = I(s(-)5[e- /C)X + e-(s/c)(2-x) ][l-e (-2/c)] }) (1-23) But Is/c = s(p/E)1/2 (1-24) Replacing E by E*(s), (1-24) becomes Qs/c = Qs[p/(E + lS)] (1-25) Substitution of (1-25) into (1-23) gives the transformed solution in real time tKo To get the transformed solution in nondimensional time, we apply the scaling theorem once more but with a = c/Q in (1-22), we get ) f(p)[e -X(p)x +-X(p)(2-x) -2[1 + e (p) ]_} (1-26) E(xp) =f(p)[e + e I ][ + e where X(p) = p(l + Ep)' 1/2 and = It/pI c Equation (1N26) is a special case (X = O) of the viscoelastic problem considered in Ref 4o The inversion of (1-26) is, in general, not easy, but if c is considered small, ioeo, c < 1, it was shown in Refo 4 that the original partial differentiation equation leading up to (1-26), ioeo, cu + u u = 0 (1-27) xxt xx tt is of the singular perturbation type, since ~ = O reduces the order of the equat;ion~ An order of magnitude analysis was presented in Refo 4 which showed that if the assumnption that c << 1 is justifiable then the increase in stress at x = O is also of the same order of magnitude 25

Example 2: Three-parameter viscoelastic rod capped by a rigid mass, otherwise it is the same as Example 1. The figure pertaining to the example is the same as Fig. 1-12 except that the constitutive equation of the material is G' + Po~' = qoc? + q1 C, (1-28) where Pl, qo, and ql are material constants. Its Laplace transform with respect to real time t' is E*(s) - (ql/pl)[(qo/ql) + s][(l/pl) + s] - (1-29) Again applying the scaling theorem, (1-22) to (1-28) yields, I.e s e-(is/c)x is~ -(Qs/ )(2-x) Ui(x,s) = f-) (1-30) C 1+ I2s -(21s/c) But Is/c = Is/(E/p)1/2 and replacing E by E*(s) yields s/ c = (ppi/qi)s[(l/pl)+s]/ [ (qo/ql )+s-1"2 (1-31) To get the transformed solution to the viscoelastic problem in terms of dimensionless time, let tv = cvt'// where cv = (ql/ppl)1/2 (1-32) Substituting (1-31) into (1-30) and applying the scaling theorem (1-22) once again, we get e e = -A(P)x -A(p)(2-x)'(x, p) = e P -2A(p) l+e(p)e where 26

A(p) = P(Ca+P)l/2 (Ca2+p)1/2 and = /P1v; 2 = qo2/qlcv c(p) = [el-(p)][l+A(p)][l+A(p)] The physical significance of the viscoelastic constitutive equations is to think of the continuum as consisting of a series of infinitesimal masses, each one connected to its neighbors by a "stress-strain" unit consisting of infinitesimal springs and dashpoto By virtue of physical arguments we must preclude those models which do not return to their original size upon removal of the stresso A comparison of the propagation characteristics of different semi-finite viscoelastic rods for longitudinal waves was made by Morrison and Lee16 Their conclusions were that the Kelvin material gives a diffusion process that is not a wave motion at all while the three-parameter solid gives damped waves, which traveled at the elastic wave velocity or a higher velocity depending on the arrangement of the "stress-strain" units. The solutions given in Ref. 16 are valid for the duration before any reflection takes place, i.e., t < lo G. LOVE-RAYLEIGH MODEL In both the elastic and viscoelastic models, the lateral effects were ignored. If it is assumed that plane sections remain plane and that the radial strains are related to axial strains through Poisson's ratio, then the differential equations for the rod was shown by Love17 to be: u ~tat - (vk)2 ut'tXtVt = C2 u',,, (1- 4) where v = Poisson's ratio, k = radius of gyration of the cross-sectional area. It was shown by Plass et al.,18 that the stress at any section is given by: a(x',t') = pA (vk)2 u,'t't' (x',t') + AE u'X, (x',t') (1-35) For the problem of the rod subjected to a displacement pulse at x' = O and with a rigid mass attached at x' = Q, the boundary conditions at x = ~ and x = O, are: 27

pA (vk)2 u' + AE u', -M u' (1-36a) u'(O,t') = f(t'). (1-36b) The initial conditions are assumed quiescent. The nondimensional version of (1-34) and (1-36) becomes 52U + U = 0 (1-37a) xxtt xx tt subject to the boundary conditions 2 ut (i,t) + u (i,t) = 4 Ut (1, t); (1-37b) where X = M/pA~ and k = vk/~. The solution of (1-37) is, in the usual way u(x,p) = f(p)([e (P)X+c(p)e- (P(2-X) ] [1+e (p)e-2(P) where C(p) = (1 - Xp)/(l + Xp) yI(p) = p(l + 52 p2) 1/2 (1-38) Note that as 5 + O, q(p) -+ p, we obtain the simple elastic model, as we must. Again, the inversion of (1-38) appears to be a difficult problem. 19 A theory due to Mindlin and Herrmann accounts for both radial inertia and radial shear stress and leads to a complicated coupled system of partial differential equations. No attempt was made to consider this case. H. A NOTE ON THE INVERSION OF THE TRANSFORMED SOLUTION In the partial differential equations considered so far, the transforms of the solutions all have quotients of exponential functions, i e., 28

'X(P)x -X(p)(2-x) u(x,p) = f(p) e +a(px)e (1-39) l+c(p,X)e In the Table shown below we note the differences in X(p) and a(X,p) between the various models: Model X(p) c(pk) I. Elastic p [1-XX(p)][l+XX(p)]II. Kelvin p(l+cp) -1/ -do(oj+p)1/2 IIIo 3-parameter solid P -— /' -do(b2+p) 2 IV. Love-Rayleigh p(l+t2p2 )1/2 [1-xpj[l+Xp]-' where c, k, cl, o 2 are system constants. The quotients (1-39) can be expanded in a series of negative exponentials as was done in Ref. 4 and the solutions obtained from the series term by term. This procedure has the further advantage of a simple physical interpretation since each term of the series corresponds to the solution of a related problem for a semi-infinite region, and thus the solution for the finite region can be interpreted in terms of successively reflected waves.20 As a typical. example of such a procedure, it was shown in Ref. 4 that the Kelvin model with a rigid mass at one end and a displacement input at the other can be written as: -X(p)x oo n n - X(p)(x+2n) -X(p)(x-2n) u(x,p) = f(p)+ (e X(x+ -l) ncn(p,X) [e -e) (x-]) (1-)o) The solution u(x,p) is such that ep occurs naturally and so if we change the variable from p to Ep and write x = x/e, t = t/E, E = Cp, (1-41) the solution for u(x,p), (1-40) may be written as u(x,a ) = f(/c)Es (e ) + nl (-l1) na(h). X()(x+(2n/)) eX(n)(x(2n/)) 29

where x and q (and hence t5) occur naturally and X(O) = \(q+1) a(nlX) [= 1-X(:n) ] [ l+Xx(n)01 1 (1-42b) In view of (1-42) the appropriate stretched variables are x, t and so the equation of motion becomes u — + u__- u_ O (1-4) xxt xx tt with the boundary conditions u(O,t) - f(t); uM(.,0) = u(R,0O) 0 (1-44a) and 1. 1~ 1 u(E) + U,( Et) -(X/E) Ut; t (1.44b) We note that for c << 1, (1-27) is certainly of singular perturbation type (since c = 0 reduces the order of the equation). The first-order singular boundary conditions are (1-44a) and U (, T) + +U (, ) = o, xt x which implies that we need consider only u(X,t) + 0, as x + oo (1-45) The latter condition (1-45) merely implies that the first-order singular solution is the solution of (1-43) in the case of a semi-infinite rod and so is given by the first term of (1-42), namely, U(.,l) = f(rl/E)e-Xrl[ +l1]1/2 (1-46) 3o

It is in the forms of (1-43) and (1-46) that Morrison obtained the exact solution to the semi-infinite rod problem, i.e., the satisfaction of the boundary condition (1-45): u(xT) + 0 as x + o. The details are lengthy and are therefore omitted here. The solution for several transit times become increasingly difficult since the coefficients of the successive exponentials in the series for the transforms are progressively more and more complicated functions of p, so that only the first few terms of the series can be used, see (1-40). Fortunately, these are precisely the times of most importance in the present problem since it is within these durations that the maximum stress is attained. A systematic program is in progress to perform these inversions. It is to be noted that the remarks above are similarly applicable to all the models mentioned up to this point. I. THE TAPERED ELASTIC MODEL Physically this model is identical to the elastic one except now the crosssectional area is no longer constant but varies according to some distribution, say A(x') = Ao(x'/k)m, (1-47) where A is the area at x' = kQ, m is an arbitrary real number and k is an 0 arbitrary constant added for convenience in prescribing the total amount of variation of the area in a bar of given length. The bar occupies the region k~ ~ x' < (kl + Q)O The differential equation of the rod is given by Bu + mu = utt (1-48a) subject to the boundary conditions: u (k,t) = - Xutt(kt) and u(l,t) = f(t) (1 -48b) 31

where X M/pAolkk and k = (l+k)/k (1-49a) and = x'/k; u - u'/ki; t = ct'/k; c2 = E/p (149b) The Laplace transform of the equations in the boundary-value problem is given by t E(, p) + rmu (ep) _p2u(Wp) = O u(l,p) = f(p) k (k,p) = -p(k p) (1-50) After (1-50) has been solved for u(W,p), then the required response is obtained by the inversion of u(W,p), which is not an easy task. The level of the difficulty is an order higher than the previous models due to the fact that our equation of motion, (1-48a), is a linear partial differential equation with variable coefficients restrained by time-dependent boundary conditionso The "correspondence principle" of viscoelasticity still applies because of the linearity of the governing differential equation. JO NONLINEAR ELASTIC MODEL Suppose that in the otherwise simple rod of Section E, we have a nonlinear constitutive relationship for the material properties of the spine, ioe., functionally: = f(c') (1-51) In view of (1-2) and (1-3) and letting 32

do'/de' = Eg(c') = df/dc', (1-52) we get the differential equation of motion: C2g(u?' )u'II ='(1-53) x' x'x' t't' In the dimensionless variables defined by (1-8), we get g(x )uxx = tt (1-54) This is a nonlinear partial differential equation because g(u ) is considered nonlinear. It was pointed out in Section C, that the spine has a hardening characteristic, i.e., g(E) is a monotonically increasing function of E. 8 Many interesting results can be obtained by a qualitative examination of (1-54) The immediate observation is that a particular form of (1-54) i.e., [(1 +u)]u = ut' t > (1-55) x xx 22 was shown by Rosen to be formally equivalent to one-dimensional flow of a polytropic gaso As is well-known, normal shocks develop under these circumstanceso A shock is defined here as the accumulation of small effects which produce steep gradients in the waveguide. The pressure jump at the normal shock in gas dynamics would correspond to a jump in the compressive stress in our problem. It is not too far-fetched to conjecture that the shock would occur near the mass end since time and distance are needed for the accumulation of small effects to materialize~ This was confirmed analytically by Murray23 recently. Of course, it is also possible to use the method of characteristics to compute its approximate location and magnitude, see either Fox24 or Fong25 The multiple and severe fractures of the thoracic spine (near Th 5 and 6) as reported by Hirsch and Nachemsonl0 now can be attributed theoretically to two causes: (a) The inertial effects of the "head" mass consisting of the thoracic cage, upper limbs and the head producing high stress levels in the thoracic region. 33

(b) The stress "jump" (due to the hardening nonlinearity) producing initial compression failure and then multiple overextension or separation upon the arrival of the tensile wave. There appears to be x-ray evidence of the same in Hirsch and Nachemson.1O Although the formal analytical closed-form solution of (1-54) is not expected, it was possible to gain much qualitative information about the behavior of the model. Work is in progress to improve these approximate resultsO Ko EFFECTS OF NONLONGITUDINAL LOADING In any given ejection problem, either due to practical considerations or the randomness of the system parameters, the displacement or any other related input would have a horizontal component. In fact, in military aircraft ejection, an eccentricity of 153 either side of the "vertical," is a standard assumption. Of course, one can formulate this dynamic beam-column problem but it is instructive to begin with just the horizontal input, both because it is an easier problem and is also meaningful in its own right. Some examples are: the whiplash problem in rear-end automobile collisions and the seat-belt (lap type) syndrome. These problems, however, are material for the next chapter. 34

References to Chapter I 1. Do Goldman and H. E. Von Gierke, "The Effects of Shock and Vibration on Man," Shock and Vibration Handbook, Chapter 44, McGraw-Hill, 19610 2. H. E. Von Gierke, "Biodynamics Response of the Human Body," Applied Mechanics Review, 17, 1964, pp. 961-958. 31 V. L. Roberts, E. L. Stech, and C. T. Terry, "Review of Mathematical Models which Describe Human Response to Acceleration," ASME Paper 66WA/BHF-13, 1966. 4. Y. King Liu and J. D. Murray, "A Theoretical Study of the Effect of Impulse on the Human Torso," Biomechanics, Symposium, ASME Publication, 1966, ppo 167-186. 5. PO R. Payne, "The Dynamics of Human Restraint Systems," National Academy of Science Symp. on "Impact Acceleration Stress," Brooks AFB, San Antonio, Texas, Nov. 1961, 6. R. R. Coermann, et al,, Aerospace Medicine, 31, 1960, p. 443. 7. D. Dieckmann, Intern. Z. Angew. Physio. Einschl. Arbeitsphysio. 16, 1957, p. 519. 8. Ho R. Lissner and F. G. Evans, "Effects of Acceleration on the Human Skeleton," Progress Report to Public Health Service for Grant No. AC00054-06, Oct. 1963. 9. T. Jo Wittman, "An Analytical Model to Duplicate Human Dynamic Force Response to Impact," ASME Paper 66-WA/BHF-12. 10. Co Hirsch and A. Hachemson, "Clinical Observationson the Spine in Ejected Pilots," Acta Orth. Scand,, 31, Fas. 2, 1961, pp. 135-145. 11. T. Brown, et alo "Some Mechanical Test on the Lumbosacral Spine with Particular Reference to the Intervertebral Discs," Jouro of Bone and Joint Surgery,, No. 5. Oct. 1957, pp. 1135-1164. 120 S. Ruff, "Brief Accelerations: Less than One Second," German Aviation Medicine, World War II, 1:584-597, Chapter VI-C, Dept. of Air Force (U.S. Government Printing Office, 1950). 155

135 C. Hirsch, "The Reaction of Intervertebral Discs to Compression Forces," Jour. of Bone and Joint Surgery, 39A, ppo 1135-1164. 14. J. Lo Hess and C. F. Lombard, "Theoretical Investigations of Dynamic Response of Man to High Vertical Accelerations," Jo of Aviation Medicine, 29, 1958, pp. 66-74. 15. W. Fligge, "Viscoelasticity," Blaidell Book Co., 1967. 16. J. A. Morrison and E. H. Lee, "A Comparison of the Propagation of Longitudinal Waves in Rods of Viscoelastic Materials," Jouro of Polymer Science, XIX, 1956, pp. 93-110. 17. A.EH. Love, "Mathematical Theory of Elasticity," Dover, p. 428. 18. H. Jo Plass, Jr., et alo, "Studies on Longitudinal and Bending Waves in Long Elastic Rods," Defense Research Lab., University of Texas Report DRL-376, CM-860 (1956). 19, R. D. Mindlin and G. Herrmann, "A One Dimensional Theory of Compression Waves in an Elastic Rod," Proc. of First U, S. Natl.o Congr. of Applied Mech,, Chicago, 1951. 200 H. S. Carslaw and J, C. Jaeger, "Operational Methods in Applied Mathematics," Dover, 1963-republication of 2nd edition, by Oxford University Press, 1948. 21. Jo A. Morrison, "Wave Propagation in Rods of Voight Material and Viscoelastic Materials with Three-Parameter Models,"tt Quart. of App Math., 14, 1956, ppo 153-169. 22. G. Rosen, t'Formal Equivalence of the Nonlinear String and One-Dimensional Fluid Flow," Quart. of Applo Math.,_ 2, No. 2, 1965, pp. 286-287e 23. J. D. Murray, private communication. 24. Po Fox, "Perturbation Theory of Wave Propagation Based on the Method of Characteristics," tJ. Math. and Physics, 34, 1955, P. 133. 25. J. T. Fong, "Elastic Plastic Wave in a Half-Space of a Linearly WorkHardening Material for Coupled Shear Loadings,t" Division of Engineering Mechanics, Stanford University, Report No. 161, May 1966.

CHAPTER II THE EFFECTS OF HORIZONTAL IMPULSE ON THE HUMAN TORSO A. THE PROBLEM Consider a passenger in an automobile constrained by a lap-type seat belt involved in either a "head-on" or "rear-end" collision. In the former instance, the action leads to the so-called "seat-belt syndrome" while the latter, to the well-known whiplash injury of the head-neck region. B, PREVIOUS STUDIES The lumped-parameter model shown in Fig. 2-la was used by Martinez et al.,' for describing the whiplash phenomenon. The mass of the head is represented by ml and the neck is represented by the cantilevered spring (of spring constant), kl, The resistance to rotation of the head relative to the neck is represented by the torsional spring, k. The mass of the seat and the body of the passenger are represented by m2 and the spring constant of the seat back is given by k2. The vehicle was subjected to a half-sine-wave acceleration pulse, somewhat indicative of what takes place during a rear-end collision, and the rotational and translational response of the head, and the translational response of the shoulders (mass m2) were computed. Figure 2-lb gives typical experimental data and Fig. 2-lc gives typical computed results. The model has many limitations. All springs are assumed to be linear and the head rotation to be small, thus allowing for linearization of the equations of motion by fiat. It was indicated in Refo 1 that these restrictions would be removed later. Like most lumped-parameters models, the results are applicable to gross motions onlyo If there is anything significant about the whiplash problem, it is that the head rotations are not small. Either an analog or digital computer solution of the full set of nonlinear equations, (2-4), would have been more meaningful. The seat-belt syndrome arises from a front-end collision. The deceleration pitches the head and upper torso forward pivoting about the ball-joint connecting the pelvis and the femur, producing high tensile and flexural stress as on the lumbar spine. As far as the writer is aware no significant analysis had been attempted for this problem except those due to Aquino, 2 which is presently in progress at the Highway Safety Research Institute of The University of Michigan, Ann Arbor, Michigan. 37

X, X0 XI xo 4 -6~ r k 0 A~~~~k c, m2~~~~~~~~~~k c (t)=Asincot CAR Lro Fig. 2-1a. Lumped-parameter model of the whiplash problem.

+3 CAR z H a-HEAD n + 2 w + UI SHOULDER -~o W O L -_ _ _ _ __ _ 5, a _ I F \ f/~~~ |Head-Human Front Car Shoulder-Human Front Car -2 Car Body- Rear Car | _' ~,,.~- I__Car Body- Front Car -3,!.. - I, i.... I, r I' 50 100 150 200 250 300 350 TIME (MILLESECONDS) Fig. 2-lb. Typical experimental data.

2 4 z. HEAD (3100 wL 2 U2 SHOULDER z cE -., ROTATION z U U <.-2 TIME IN MILLISECONDS -100 K1 =15 K2= 30 KT =285 -200 Fig. 2-1c. Typical computed results.

C. LUMPED- PARAMETER MODELS Using the notation shown in Fig. 2-la one can show easily that the equations of motion for finite head rotation are: mlxl + kl(xl-x2) - mlQj2 sinG + mlj cosO 0 = 0 m2x2 + kl(x2-xl) + k2(X2-Xg) = 0 JOG + k G - mlgi sing + mlxl1 cos = 0 (2-1) where the springs are all considered to be linear and JO is moment of inertia about point 0. Let Y1 = x1 - x and Y2 = x2 - x (2-2) then, we get mlyl + kl(yl-y2) - mmlg2 sing + mll cos0 = -= mdg m2y2 + kl(Y2-Yl) + k2Y2 = -m2xg II II Jo0 + ktQ - mlg~ sinG + mlyl~ cosg = -ml~ cosQ xg (2-3) One obvious generalization would be for the spring to be nonlinear, i.e., mlyl + f(yl-Y2) - mjl~2 sing + mil cosQ = -mlxg m2Y2 - f(Yl-Y2) + g(Y2) = m2xg Jio + h(G) - mlg~ sing + mlyjlcosQ -ml cosQg, (2-4) 0 9~~~~~~~~~c

where f(yl - Y2) and g(y2) and h(Z-`) are all nonlinear functions of their arguments. There is not much hope of solving (2-4) analytically. Either a RungaKutta procedure can be used to obtain a numerical solution or an analog solution would also be satisfactory. In its simplest form, the seat-belt syndrome problem, Fig. 2-2, can be shown to yield the following differential equations of motion: G + (21g/~) + ( i cos0/ ) - (g sinG/i) = 0 (2-5) Q + Q sing - ~2 + [k(e-2o)/m] + g cos - O 0 where ~o is the undeformed length of the spring, (spine) which is rather difficult to determine. Let x = 2 - ~o + (mg/k); ~e = O - (mg/k) and 2 k/m; 2 - g/e (2-6) we rewrite (2-5) as (1+ -) + + -n cosQ - 2 sinG = 0 Xe +L @2 X + \ sing - (Xe - g7l-cos@) = (2-7) I + x + rj sin0- (x+~e)42 - g(l-cosG) = 0 Again, a digital or analog solution seems most appropriate, although the chances for an approximate solution is decidedly better than in the whiplash problem, e.g., (2-4). Of course, one can always complicate things by requiring the spring to be nonlinear. D. THE ELASTIC BEAM MODEL The model is the same as in Fig. 1-1 except that the input is a horizontal displacement or acceleration, see Fig. 2-3. The elastic rod represents the 42

... — Impact ae Fig. 2-2. Lumped-parameter model of the seat-belt syndrome. 43

I)(et, V() i(Xt IL I A _- h (t) Fig. 2-3. Continuum model of the whiplash problem.

cervical spine. The obvious deficiency of the model is that the "real" problem is a so-called mixed boundary-value problem, i.e., as the load is being applied more and more of the torso comes into contact with the seat surface and thus redistributes the forces and moments acting on the rod both in space and time. The justification for our model is that the details of the motion from the time of impact until the shoulders are in contact with the seat surface is short. As far as the whiplash injury is concerned, we might as well consider everything below the cervical spine as a lumped mass which moves with the same displacement (or acceleration) as the seat. Let the end x = O be considered fixed to the seat which is given an acceleration h (t) The formulation of a similar problem was given by Karnopp. From tt Fig. 2-35, we note that the transverse motion of the rod is given by y(x,t). It is convenient to define a relative coordinate, 5(x,t), where -(x,t) = y(x,t) - h(t) (2-8) The governing differential equation, followingt is: EI xxxx +pA tt - pAhtt (2-9) The boundary conditions are: EI xx (~Q t) = M[ tt (n,t) + qxtt ( t)] + ht (t) (2-10) I x(x t) q= [Mqtt ( t)+ (Jc+Mq )i xtt(t) - Mqhtt (t)], (2-11) and the initial conditions are (O, t) = ~t(o,t) =- (2-12) where I is the moment of inertia of the beam cross-section, Jc is the polar moment of inertia of the head about its center of mass and q is the distance from the atlas (C1) to the center of mass of the head. Eqs. (2-10) and (2-11) represent the shear force and moment balance at the interface between the rigid mass and the elastic rod.

For any given acceleration pulse, h (t), waves propagate up the rod. For impact-type loadings, however, the Eulertternoulli assumptions leading up to (2-9) to (2-12) was shown to be inadequate by Lamb4 as early as 1917. He showed that these assumptions gave the physical untenable result that the effect of a suddenly applied load is propagated infinitely rapidly. Fluigge and Zajac5 showed that the inclusion of the effect of rotary inertia and shear deformation, i.e., making Timoshenko beam assumptions, would free us from the objections of the elementary theory. The revised problem, although quite complex, nevertheless, still possesses an additional objection. In the whiplash problem we are dealing with large deformations and rotations. Our present formulation is valid only for small deformations even with the Timoshenko beam assumptions. It is desirable to formulate the problem on the basis of finite displacements but small strains. but that will- be another story all together. The above remarks are similarly applicable to the seat-belt syndrome problem. Because of the enormous difficulties encountered just in the problem formulation, it was decided to initially model the problem as a lumped-parameter system of many-degrees-of-freedom, capable of both translation and rotation. The approach is discussed and amplified in Chapter IV. Only after a better understanding has been achieved will continuum models be constructed.

References to Chapter II 1, J. L. Martinez, J. K. Wickstrom, and B. T. Barcelo, "The Whiplash Injury-A Study of Head-Neck Action and Injuries in Animals," Biomechanics Monograph, ASME, 1967, pp. 171-183. 2. C. F. Aquino, Dept. of Mechanical Engineering, The University of Michigan, private communication. 5. Dean Karnopp, "Random Vibrations," Vol. II, Ed. S. Crandall, Chapter 1, MIT Press, 19635. 4. Lamb, H., "On Waves in an Elastic Plate," Proc. Roy. Soc., A93, 1917, pp. 114-128. 5. FlUgge, W. and Zajac, E. E., "Bending Impact Waves in Beams," IngenieurArchiv., 28, 1959, pp. 59-70. 47

CHAPTER III A MODEL FOR HEAD INJURY-THE AXISYMMETRIC RESPONSE OF A FLUID-FILLED SPHERICAL SHELL TO A LOCAL RADIAL IMPULSE Nomenclature A,A A Arbitrary constants (Constants of Integration) o nm nc E Young's Modulus Fe External force distribution on the shell Mz Moment along polar axis Pn(cos 5) Legendre Polynomials of the first kind. P (cos ) Associated Legendre Polynomials of the first kind and first order S Midsurface of the shell T Kinetic Energy V Potential Energy VFo, Fs Volumes of fluid and shell respectively U Strain energy density a Radius of spherical shell ao, an Coefficients of Legendre polynomial expansion of G bn Coefficients of associated Legendre polynomial expansion of 1 c Compressional wave speed in the fluid Cs Wave velocity, [E/ps(l1v2)]1/2 f Shell-fluid parameter h Shell thickness in Spherical Bessel function k Wave number 48

pa Fluid pressure on the surface of the shell r,G, ~ Spherical coordinates of the deformed midsurface s Speed ratio, c/cs t Time u Meridional displacement with respect to center of mass of the shell w Radial displacement, a-r w Radial displacement with respect to center of mass of the shell z Distance from the midsurface a2 Thickness parameter, h2/12a2 ao,'Cnc,Cnm Phase angles (Constants of Integration) B Angle between midsurface normal and radial ray y0 Midsurface shear strain &nc, 5nm Amplitude ratios OEO- cor Midsurface normal strains E,C Z-surface normal strains Nondimensional radial displacement, l/a(a-r) Nondimensional angular meridional displacement $- a, q, Spherical coordinates of the undeformed midsurface ko _, koq Midsurface curvature,%n n(n+l) v Poisson' s ratio po ps Mass density of fluid and shell respectively oa -- Normal stresses T Dimensionless time, ct/a 49

nc nm Angular frequencies of breathing, composite, and membrane modes respectively. 0 Nondimensional frequency, (i4/c 37 Velocity potential for the fluid

A. THE PROBLEM AND PREVIOUS STUDIES We shall attempt to model the situation when the skull is subjected to a blow. The idealized model consists of a closed elastic spherical shell filled with fluid, Figure 3-1. As a first approximation, the shell is considered to be thin and elastic while the fluid is inviscid and irrotational and the blow is manifested by a sudden initial velocity input. Even with these simplifying assumptions, the problem is quite complex and it was decided to break the problem down to its constituent components: (1) A rigid closed spherical shell with irrotational, inviscid fluid as the wave guide. (2) The elastic closed spherical shell as a solid wave guide. (3) The fluid-solid interaction of the complete problem. 1. 2 Part (1) has had some attention by Anzelius and Guttinger. They concluded that the initial velocity input produced a compression wave at the point of impact; however, because the shell was assumed rigid, the effect was instantaneously transmitted to the counterpole, whereupon a tension wave is simultaneously emitted. The collision of the two waves at the center producing caviationtype phenomena, which was considered the cause of damage. The obvious defect in the model led Goldsmith5 to suggest the analytical or numerical solution of a fluid-filled elastic shell. Gold-smith's paper5 is both a tutorial and a survey. An exhaustive review is given of theoretical and experimental approaches previously (up to 1965) employed to delineate the mechanisms and to measure mechanical quantities believed to be important in the production of brain trauma. As a problem in theoretical and applied mec anics, the proposed models have had a long history. As early as 1882, Lamb used an extensional formulation in his study of the radial motion of closed spherical shells. Love's5 formulation includes both flexural and extensional effects and this has become classical bending theory of shells now known as Love's approximation. Based on the extensional theory of Love, Silbiger obtained certain analytical and Baker7 some experimental results. Naghdi and Kalnins8 obtained a solution for the torsionless axisymmetric motion using the classical bending theory and also made a study of asymmetric motion based on extensional theory. Kalnins9, using classical theory, explained certain paradoxes in Love's frequency analysis in terms of the effects of bending. Recently, McIvor and SonstegardlO studied the axisymmetric response of a closed spherical shell to a nearly uniform radial impulse and the associated parametric stability problem. The problem of a fluid-filled shell submerged in another fluid has attracted many researchers, especially in the discipline of acoustics. Jungerll calculated both the reflection from an air-filled shell submerged in a fluid and the 4t-ransmission through the interior of an incidental plane wave. Goodman and Sternl2, using elasticity theory, studied a plane wave impinging on a fluidfilled spherical shell, which is itself immersed in a fluid. Both11 andl2 51

HEAD IMTUR'N MODEL L0CITu. AL SECto LO ITU SECTSKULL OF SKULL Fig. 3-1. Head injury model. 52

are based on the numerical integration of the system of ordinary differential equations, which were obtained from the wave equations governing the propagation of disturbances in the 3 different media occupying the eleastic space. The results were valid for steady state motion since the time-dependent term eiot. had been assumed for each field variable. Hickling13 extended the results to a pressure pulse emanating from a point source. The transient response to such an impulsive pressure can be found by integrating, over a suitable range of frequencies, the product of the steady-state response and the Four'er transform of the applied impulsive pressure. Recently, Rand and DiMaggio obtained frequency equations and mode shapes for the axisymmetric, extensional, torsionless vibrations of fluid-filled elastic spherical shells and rigid prolate spheriodal shells. B. THE IN VACUO CASE The equations of motion for free vibration of a closed spherical shell may be derived using an energy formulation as was done inlO. In this reference, Lagrangaian representation of spherical shell deformation was used. If the undeformed and deformed configurations of the midsurface are represented by spherical coordinates: a,r, q and r,@, respectively, then in the Lagrangian sense, the deformed state coordinates are taken as functions of initial configuration and time, i.e., r = r(S,rl,t). 0 = @(S,,t) and A = 1(,1,~,t). Defining the following nondimensional variables: = ~ - 5 (angular meridional displacement) = (a-r)/a (radial displacement) T = [E/ps(1-v)/2 if c = [E/Ps(l-v)]/2 then T cst/a a2 = h2/12a2 (a thickness parameter) For torsionless axisymmetric motion, r = 1(,:T) and T = ~(i,T) i.e., 0/. = O. Denoting 3/)5 by a dot, we get two linear partial differential equations if both the displacement and its derivative are assumed small: + cot, 5 - t(v+ cot2S) - i(l+v) + 62[*+j* cot 5 - *V(v+ cot2 b) +' +' cot ( - (v + cot2 3)] = a2/aT2 (3-1) 53

(l+v)(llr + 4 cot t -25) - C2[t + 2;i cot 5 - t(l+v+ cot2 5) + x cot 5 (2-v+ cot2S) +' "+2' cot ('+(l++ cot2)+ cot V(2-v+ cot2 5) = 2 (3-2) Assume the displacements ~ and r are represented by; 00 = Z an(T)Pn(cos 5) (3-3) n=O 00 E bn(T) PA(cos ~) (-4) n=l where Pn(cos 5) are Legendre Polynomials of the first kind, (second kind is singular at the poles) and Pn(cos ~) are the associated Legendre Polynomials of the first kind and first order. One can show that10 the square of the frequencies are: = [2(1+v)] (3- 5a) 2n2 [1 + v - 2( v) + n(l+va2) + a2 ~ ([1 + 3v - c2(l-v) + k (l+va2) + a2X2]2 4[-2(1-v 2)(l+c2) + n(l v2) + a2 _n(V ) 4a2 X + ~ X3] // + (3- 5b) where positive and negative signs give frequencies Thm and wnc respectively and kn = n(n+l). The amplitude ratios are: Bnm P1 1 + v - 2[1-v-n(n+l) ql ~m (3-6a' Anm A ql2m m + (1+ )[ln- n(n+l)] Bnc = P1 _ 1 + v - ac2[l-v-n(n+l)] nci 1v + (+on2)[v-n(n+l)] -7) A,, ql-cu~54

The solutions for ~ and' are: 00 5(5,T) = AoSin(W0T+o 0) + Z [AnmSin(n T+nm) + Ac Sin(ncT+anc)]Pn(cos ) (3-7a) 00 = ) n=l [EnmAnmSin(LhmT+gnm) + ncAncSin((nlcT+cnc)]P' (cos ) (3-7b) For this problem, the displacement initial conditions are: ~(,0) 0= (3-8a) v(, o) = 0 j (3-8b) The two velocity initial conditions are not obvious. Its determination involves several steps: (1) Express the local velocity input by a series expansion If w = a-r, let t = vnPn (cos ) <:t n=0 nn (5o < ~ < ~ (3-9) Multiply both sides of (3-9) by Pm(cos 5) s'in n and integrate, keeping in mind the orthog(cna_'ity of the Legendre polynomials in the interval [-1,1], we get v ~ Pn(cos ) sin ~ d$ Vn = o (3-10) Pn p2 (cos 5) sin 5 d( 2 The norm of Pn; N(Pn) = +1 Pn(x)d 2 n = 0,1,2... -1 2n+l or N[Pm(cs ~)1 = O~-P2 (cos ) sin = d 2 (3-11) From (3-10) and (3-11), 55

(2n+l) v jo Pn (cos ~) sin ~ d~ (3-12) However, using the transformation x = cos 5 in (3-12), it is easy to show that = 2 (1 - cos o0) Vn = [Pn-1 (Cs o) - Pn+l (cos o)] (3-13) Hence, vn is known tor every n. (2) The linear momentum along the polar axis PP', Mz, is the component of the radial momentum due to initial radial velocity v for 0 < 5 < +~o, see Fig. 3-2. Let vz be the component of v along the polar axis PP', then the linear momentum along the polar axis is: M vz pdV f cos phdA, Z p.s.v. p.s.v. at where dA = 2ra2 sin 5 dS p.s.v. = partial shellvolume. Or in view of (3-9) 00 M 2ira2ph O v P (cos n) sin ( cos g d5. (3-14) z 0 n=O nn The mass of the shell consistent with thin shell theory is: 4 h3 mtotl = 4- pr[3ha2 + -] 4Tpa2h Hence, the velocity imparted to the center of mass of the shell (in the negative polar axis direction) is 0i Z 2inZO VnPn (cos 0) sin 5 cos e d- (3-15) Using the transformation x = cos i, one can show after a little manipulation that 56

p z Fig. 3-2. Momentum considerations for the shell. 57

1 Vn Vn-1 Vn+l vz vo(1l- cos 250) + Zvn cos ] 8Z = 7 vo8 - COS + n=l (2n+l)v [V COS -o + (2n-l) (2n+3) 16) where use was made of (3-13). (3) To find the relative distribution of the deformation velocity, divide the close spherical shell into two portions: (A) 0 < 5 < +t0 (Region where the velocity input is applied) (B) Outside of region (A) + 5o < 5 < T For region (A), call the deformation velocity in radial direction aw/6t- and in the meridional direction au/ t, then, ~a ~ -— v Cos v = z (cos w at at - az z (cos) (3-17) - - vz sin = -vz PI' (cos 5) Writing (3-17) more explicitly (using (3-9)), we get'"- co c'~w 00 at =vo + (Vz-vz) P1(cos ) + n vnPn (cos 5) (318a) at - VzPi'(co s ) (-18b) wherel vn is determined from (3-13) and vz from (3-16). Initially, for region (r) the shell suffers only the velocity which is imparted to the center of mass of the shell in the negative polar-axis direction. This initial velocity is v. Hence, the two velocity initial conditions in terms of the variables of the problem are: aT cs at for O < " < + o Cs st.0 aT s atCst for 0 < <+ t This completes the determination of the initial conditions. 58

Application of the initial conditions (3.8a) and (3.8b) into solution (3-7a,b) lead to o = Cm = nc = 0~ Hence, from (3-7a,b),''____A_ [A T A (3-19a) aT = Aozo + n [Anmhnm + Anc(nc] Pn(cos ) c at __~(t ) = [- A w + o A ]P'(cos ) = -- (3-19b) T - n=l nm nm m n nc n n Cs at From (3-18a) and (3-19a), and similarly (3-18b) and (3-19b), one can conclude that; vo Ao = for n = 0 Almclm + Alcwzc = CsV1 (3-20a) for n = 1 1 51m Alm wlm + 51c Alc c C v (3-20b) Anm %nm + Anc "ncc = - Vn (3-21a) From (o-20a) and (3-20b), we get Aic = 5+Y(l- 6m) and Aim= 5v+vz(lc5 csauim(3im-Flc) cscnIM(3ic-Sim) Similarly, from (3-2la) and (3-21b), for n > 2 Anm n nc Vn Am- and A;nc "nc cs(Snm-Snc) nm - = m cs(3nc-3nm) Since all the Anc- Anm- bnc- 6bnmy C'nc and oUnm are known, the solutions (3-Ya) and (3-7b) are completely determinate, i.e., explicitly:

Vo 5ic vl + vz(l-51c) t(~,T) = C no sin c (c m) sin (WImT) PI(cos. ) 00 6rim, 6 + n-2 sis [cnc(nm-$nc) sin (nc) + sin ( T)] Pn(cos ~) (% nc- nm)wnm (3.22a) 51c vl + Vz(1-31c) 1~(~, T) = -5lm Cs Vjm ( icm) sin (jmT) sin 5 oov F nnmsin (uncT) sin(w nmT) n cos P(cos S)-nP l(Cos C] s.n.m co.... -.. n=2 Cs L nc n nc) + Wnm(mnc.nm L sin. (3.22b) where v'z, vn, v nm and 6nc, enm and %nc are compu'ted from Eqs. (3-16), (3-13), (3-6) and (3-5) respectively. It was shown in Ref. [10] that for a/h > 20 and n 2 2, the coefficient in the second term under the suammation in (3-22) can be neglected. The linear midsurface strain and rotation quantities, in terms of (, T) and *(~,T), are: c t(gT) - - (3-23a) midsurface normal strains Cor(T) = - +'f cot (3-23b) k5(j,7T) = ll + 1 + _2_ (3-24a)? midsurface ao( ) 1 sin 1 curvatures (3-24b) k T(,a) sin ) (1+ cot ),-24b) (gT) ==the angle between t.he surface normal and -the radial ray. (3-29) For isotropic and homogeneous continuum, Hooke's Law yields the biaxial stress-strain relations E CY(,T) (E + (V-26a) (t, T) = 12 (c + y 2T) (3-26b) 60

where Ee and cr are the z-surface strain quantities. The z-surface strain quantities are: 1 z E(S, T) = Eot + z(ko-a-)(l+oS)(1-a) (3-27a) i z E (ST) = + z[kon sin (4-))-a](l+EC)(l- ( ) (3-27b) where ) = T(5,) + Two sets of graphs will be made from (3-26): (a) Ca(ST) and ao(j,T) will be plotted versus 5 with a progression in time T, i.e., we wish to show the progress of the initial velocity input with respect to the stresses it creates in the shell. (b) at (,T) and O (n,T) will be plotted versus time, To at a given angular meridional displacement, i.e., how the stresses fluctuate at different points in the shell. The computations and automatic plotting are presently in progress. C. THE FLUID-SOLID INTERACTION (1) Formulation of Equations of Motion In order to use Hamilton's principle to obtain the governing differential equations of a fluid-filled spherical shell, it is necessary to calculate the kinetic and potential energies of the region occupied by an ideal fluid and the thin elastic shell surrounding it. The potential energy of the region is -VPo a )f2 v = f UdS - fS p wdS + fS FewdS + 2c2 ( 3-28) where 61

U = strain energy density per unit area Pa - fluid pressure on. the inner surface of'the shell w radial displacement of the shell Fe = external force distribution on the shell p mass density of the fluid c = compressional wave speed in the fluid = velocity potential for the fluid S = midsurface of the shell V -= volume of the fluid The first term in (3-28) is the strain energy of the shell, the second and third terms are the work done due to the internal fluid pressure and the external pressure pulse respectively while the last term is compressive strain energy stored in the fluid. The strain energy, per unit area of the thin shell is well-known, see Novozhilov5, i.e Eh Eh3 U 2(lV2)[( c) _2(1-v)(e -7g )1 ein+ Eh( v2)[(k1+k0) - 2(1- v)(k~ k-T2)]. (3-29) The midsurface strain-displacement relations in the case of axisymmetric torsionless motion of the shell are: lEu 1 62w 6u E q = a( - w) k a 2 ( +Z EQ = ~(u cot;A-w) k 2 + u) s o = O (3-30) The kinetic energy of the region is: 1 LU 2 6w 21 2 T =- p [(-) +a-)]d + Sd P cV2 db0 (3-31) 5 0 where Ps and po are the density of the shell and fluid respectively. Applying Hamilton' s principle f J (T-V)dt = 0 together wiith eqs. (3-28) to (5-31), results in the equations of motion of t;he region along with some natural bouadary conditions: 62

(l+o12)[-I- -u4cot 4) + (v+ cot2 4)u] -o e - oCW cot 4 + [Io(cot2 4+v)+l+v]w + 12 sa Utt = 0 (3-32)'u,<, + 2'uo, cot 4 - [(l+v)(l+cX2) + o cot2 q]!u~ + [Y2 cot3 4 + 3a Cot; - (l+v)(l+cx) cot 0]u + oE[w~O, + 2 cot 4 w... - (l+v+ cot2 4)w', 2 + (2 cot 4 + cot3 4 - v cot ))w0] + 2(l+v)w + NE psa2Wtt + 1-v2 a[Fe(O,t) - Po t(a')t)] (3-33) 1 2 1 2[r2 (r rr + r2 sin ) (sin -))] - t = (3-4) where, as before, the subscript notation is used for partial differentiation. In (3-33), Pa = -PoOt(a,4,t) is the dynamic fluid pressure acting on shell surface. (2) The Frequency Equation From (3-34) the form of velocity potential 0 can be expanded as 00 (r,,t) = n cn(t)Pn(cos $)jn(kr) (3-35) The boundary condition between fluid and shell can be stated as the continuity of normal velocitiesi.e., wt($,t) = ~r(a,4,t) (3-36) Eqs. (3-32) and (3-33) are nondimensionalized by using t = u/a, f = -w/a, t = cst/a, where cs is the wave velocity in the shell. They become ~+~ cot V - 4(v+cot2 9) - (l+v)~ + d[V+fcot 5 - i(v+cot2 ) + + * cot 5 - (v+cot2 )1 ] = a2/aT2 (3-37)

(l+v)(Q.+t cot - 25) - c-['*;+2;r cot,- (l+v+ cot2,) + 4 coto, (2-v+ cot2,) + ~ + 2' cot ~ - S(i+v+ cot2 ) + 5 cot 5 (2-v + cot2 5)] po; a 2 2) + a' (a, =, T)/' Fa (3-38) Ps TEh F() We note that eqs. (3-37) and (3-38) are idenrtical to (3-1) and (3-2) except for the two terms involving:/6T and Fe in (3-38), In the absence of the external pressure pulse, we would like to obtain the frequency equation including the effect of the fluid. Analogous to Section B, we assume the nondimensional displacements, ~ and., are expanded as in (3'-3) and (3-4). The system of differential equations for the coefficients an(T) and bn(T) are: fon2a,(T) + 2(l+v) for n _= 0 dT2 1 __ ao(T) = 0 (3-39) K n ( ~) d2an(T) F L for n > 1, [1+ I + l+v n(1v) b( + ~2(l+v) + C[n- n(1- v)]2an(T) 0 (3 d2b(T)n() 0, d 2 + s(l+V) - [(l-V) Xn] ( an(T) - [(lv) - A1l](1+ bn(T) = 0 (3-41) where f = poa/psh - nondimensional fluid-shell parameter and kn = (n+l)n. We assume the solution of (3-39), (3-40) and (3-41) to be of the form an(T) = An e s bn(T) = Bn eST, (3-42) where s = c/cs ratio of the compressive wave speed in the fluid to the wave speed in bending. Substit'ution of (3-42) int;o (3-39) to (3-41) yields: 64

for n - 0 s22 2l+v) O (-4 [ 1+fa j — j]l (l+v) J 1-l-v- l ] (1+027) 1+n+ n<2(l-v) ] (2-n n On nJ folr n 1 p + (Q ) {jl+V) (2[l-v-X%~.](l [~~)+cY2)+X +v l-V-hnfQj+C2i*~2-hi2I- (2344) It i iinterest;ing to note that the limiting cases of eqs. (>43) and (3-44) agree with the results obtained previously, as it mcust. Case 1 f = O corresponds to the absence of fluid, i.e., we get the iin-vacuo case given by eq. (-355) in Section B. Case 2 s i O correspords to a rigid shell and the equation degeneJ.ates to which is easily shown to be equal to Equat;ion. (3-45b) wa.s obtained by GUttirnger.2 Case 35 2 = O yields the frequency equation correspc:l dri g it-. tJhe membraune (exgtensional) theory. frequ2enc~y' spectra and mode shapes.'These are, of course, a prelude to the.,ransie:nt response probleme 65

References to Chapter III 1. Anzelius, A., "The Effect of an Impact on a Spherical Liquid Mass," Acta Path. et Microbio. Scand., Supplement 48, 1943. 2. Gluttinger, W., "Der Stosseffekt auf eine Flussigkeitskugel also Grundlage einer physikalischen Theorie der Entstehung von Gehirnverletzungen," Zeit. F. Naturforschung, Teil A, Vol. 5, 1950, pp. 622-628. 3. Goldsmith, W., "The Physical Processes Producing Head Injury," Proceedings of the Head Injury Conference, J. B. Lippincott Co., 1966, pp. 350-382. 4. Lamb, H., "On the Vibrations of a Spherical Shell" Proceedings of the London Mathematical Society" Vol. XIV page 50. 5. Love, A.E.H., "On Free and Forced Vibration of an Elastic Spherical Shell Containing a Given Mass of Liquid" Proceedings of the London Mathematical Society Vol. XIX pages 170-207. 6. Silbiger, A., "Free and Forced Vibrations of a Spherical Shell" ONR Report U-106-48, Dec. 1960. 7. Baker, W. E.,'tAxisymmetric Modes of Vibration of Thin Spherical Shell" The Journal of the Acou-stical Society of America 33 (1961) page 1749. 8. Naghdi, P. M. and Kalninrs, A., "On Vibrations of Elastic Spherical Shells" Journal of Applied Mechanics 29 (1962) page 65. 9. Kalnins, A., "Effect of Bending on Vibrations of Spherical Shells" The Journal of the Acoustical Society of America' 36 (1964) page 74. 10. McIvor, I. K. and Sonstegard, D. A., "Axisymmetric Response of a closed Spherical Shell to a nearly Uniform Radial Impulse" Jour. of the Acoustical Soc. of America, Vol. 40, No. 6, Dec 1966., pp. 1540-1547. 11. Junger, M. C., "Vibration of Elastic Shells in a Fluid Medium and the Associated Radiation of Sound" Jour. of Applied Mech., Vol. 74, 1952, pp. 439-445. 12. Goodman, R. R. and Stern, K., "Reflection and Transmission of Sound by Elastic Spherical Shells" Jour. of the Acoustical Soc. of America., Vol. 34, No. 3, March, 1962, pp. 338-344. 66

13. Hickling, R., "Analysis of Echoes from a Solid Elastic Sphere in Water" Jour. of the Acoustical Soc. of America, Vol. 34, No. 10, Oct., 1962 pp. 1582-1592. 14. Rand, R. and DiMaggio, F., "Vibrations of Fluid-filled Spherical and Spherical Shells" Jour. of the Acoustical Soc. of America, Vol. 42, NM. 6, Dec. 1967. 15. Novozhilov, V. V., "Thin Shell Theory" 2nd Ed., Noordhoff Lid., GronInger, Holland, 1964. 67

CHAPTER IV NUMERICAL AND ANALOG SIMULATION IN BIOMECHANICS A. INTRODUCTION AND PREVIOUS STUDIES The concept of numerical and analog simulation of biomechanics problems occurred quite naturally because of the apparent complexity of the system. The type of analyses carried out in Chapters I and through III, despite its steady improvement, is capable of solving only severely idealized situations. For each additional constraint introduced into the boundary value problem, the intractability of the analysis increases, sometimes by several orders of magnitude. It may turn out that this additional constraint might make only a minor difference in the quantities of interest. Is there anyway of knowing this a priori? The answer is a qualified yes. Before a numerical or analog, computer-based solution of a complex biomechanics problem can be achieved, it is necessary to limit their infinite degrees-of-freedom to a finite, if large, number of unknowns. One possible discretization is the concept of the multiple-degree-of-freedom lumped-parameter system analyses. In this type of approach, one would have available a very large number of parameters, whether or not the system is linear or nonlinear. With apparently inexhaustiable combinations and/or permutations of the system parameters, one can almost always fit any set of data. To avoid such curvefitting, which is of dubious value, one can resort to either of the following techniques: (1) A careful and systematic compilation or collection of biomechanri.cal data of the components of the system. We can then "claim" knowledge of the system dynamics. Any given input into this system will produce a computable output. This is the so-called forward analysis or output; problem. (2) Having available both the input and output records one is to determine the intervening system dynamics. This is the so-called identification problem. This technique demands very carefully planned experliments and fairly good ideas of the models for the systems dynamics, i.e., is the system linear or nonlinear? Most inputs can be classified as either transient, periodic or stochastic. Typical of the forward analysis problems are simulations due to Coermann et;. al.,1 and Naab2. In Ref. [1], a seven-mass spring dashpot system was used to represent the low frequency response of a sitting or standing human body. Experimentally obtained impedance versus frequency curves were given with position and muscle tone as parameter. No attempt was made to correlate the experimental results with the quantitative (numerical) values;to be assigned to the 68

parameter elements in the system. In Ref. [2], an 11-degree-of-freedom, lumpedparameter system was used to simulate the automobile crash victim. Unfortunately, through an oversight, no information was given concerning the range of numerical values assigned to the system parameters, hence the output is meaningless to the reader. Any parametric study must be based on some biomechanical data, no matter how crude, otherwise it is of little value. The problem connected with system identification theory has been in the forefront of recent research activity in system dynamics. Typical of such researches in linear, lumped-parameter models is the paper due to Ho and Kalman3. In it, they gave a new algorithm for the effective construction of a linear, finitedimensional dynamical system from its input-output description. 4 Shinbrot has derived a method for determining constant parameters in ordinary differential equations. Perdreauville and Goodson5 have extended Shinbrot's concepts -to partial differential equations. Using either the normal operating records or experimental data, the arbitrary parameters of an assumed (partial differential equation) model of the physical system can be determined. The method is good for linear and/or nonlinear equations with variable coefficients. The accuracy of the results depends on the "exactness" of the model, the amount of data used, the error in numerical integration and the noise in the data. If the system is well defined in terms of differential equations, then discretizations are possible. The well-known method of finite differences comes to mind immediately. There exists a calculus of finite differences, which allows for the approximate solution of the exact differential equa'tions of motion. On the other hand, the finite element methodb allows for the "exact" solution of an approximate system of differential equations. Either one of the above techniques could be used, for example, on possible improvements of the head-injury problem in Chapter III.:B. NUMERICAL SIMULATION OF THE SPINE UNDER VERTICAL IMPULSE (1) Anatomy The supporting structure of the body is a joined framework of bones called the skeleton, which assists in body movement, supports the surrounding tissue, protects the vital organs and other soft tissues, manufactures blood cells and provides a storage area for mineral salts to supply the body needs. The vertebral column is the main load-carrying part of the skeleton under the type of loading associated with the ejection and whiplash problems. A detailed study of the vertebral column reveals it to be an extremely complicated structure consisting of fairly rigid bone segments, the vertebrae, connected together by means of ligaments and intervertebral discs. It is an efficienrt, 69

pliable structure, at least under static loads associated with an erect posture. The spine is tapered,being much larger at the base (where it. is supported by the sacrum-pelvis) than at its upper end (where it supports the head). Note in Fig. 4-1 that the spine is far from straight;, the curvature in the thoracic (chest) region being opposite to the curvature in the cervical (neck) and lumbar (lower back) regions. In the seated position the human frame is supported at the ischial tuberosities of the pelvis. In this position, the pelvis is subjected to some rotation since the ischial tuberosities in general do not lie on a plumb line through the lumbo-sacral joint. However, the presence of a restraint system, such as lap seat belts, usually restricts rotation of the pelvis so that longitudinal impact to the spine through the pelvic region is transferred to the spineproper with little rotational imput. The head is supported by the atlas (first cervical vertebra) at essentially a hinge. The eccentric mass of the head, Fig. 1-1, is balanced by the muscle tension in the neck muscle. The intervertebral discs are frequently thought- of as the "shock absorbers' of the spinal column since they have a great; capacity -to absorb energy. They comprise about 25 percent of the total length of the spine between the head and the pelvis. The vertebrae themselves are quite rigid in comparison to the discs. (2) Factors Not Previously Considered The continuum models of Chapter I, even if they can be solved explicitly can not possibly account for all the factors brought out by even a cursory examination of the anatomy. There is some question as to the relative importance of some of these factors, either singly or in combination. A few examples of these factors are: (a) nonuniform distribution of mass and section properties (b) initial curvature of the spine (c) viscoelastic constitutive equations of the intervertebral discs (d) large deformations (e) effects of radial inertia and shear deformation (f) eccentricity of the head and thoracic cage (g) acceleration pulse shape and rise time The answer to this and many other questions requires a thorough parametric study of a model with enough detail such that predictions of where, when and how injuries occur would achieve reliability. The model must include much of the currently available biomechanical data. Obviously, such a model has merits all its own in terms of the actual solution of a given biomechanics problem. 70

TI (:wni ~T T~7 1 In ~~a m Per-~

(3) The multi-degree-of-freedom lumped-parameter model This is an articulated system which consists of vertebrae, intervertebral discs, the head, and perhaps, the pelvis. The vertebrae are considered to be rigid bodies in which the effective local masses are lumped, and the intervertebral discs are the visocelastic restoring elements. The most crucial aspect of this model is the nature of the kinematical constraints which will be discussed in some detail below. As in the continuum models, the lumped-parameter system is confined to motion in the sagittal plane, i.e., the plane which passes through the center of the spinal column and divides the body into two more-or-less symmetrical halves. The spinal column is initially curved prior to application of external load to the body. The configuration of the spine at any time can be given by the coordinates of two points on each of the vertebrae, namely the intersection of the vertebrae axes with the vertebrae end plates, as shown in Fig. 4-2. These coordinates thus define the size of both the vertebrae and the intervening discs. Because the articular processes of the vertebrae prevent translation of one vertebra relative to the vertebra above or below, the vertebrae are constrained to move in such a way that the'terminal points on the axes of the vertebrae must move along lines coincident with the axes of the adjacent discs, as is illustrated in Fig. 4-2. Thus, rotation and compression (or tension) of the discs occur, but shear deformation in the discs is neglected, somewhat as in the Euler-Bernoulli beam theory. The horizontal, vertical, and rotational motions of the rigid masses are all coupled because of the nature of the above constraints. This can be seen by reference to Fig. 4-2 in which only local coordinates are used. For a typical vertebra, point P is constrained to move along a line passing through PP'. Similarly Q moves along a line QQ'. In the figure, the line 2 represents the length of the vertebra, @ represents the angle made by the intersection of the axes of two successive discs, and a is the angle made by the intersection of the axis of the lower disc and the axis of the vertebra. The quantities a, b, 0, and 2 are fixed and can be computed from the instantaneous configuration of the system at any time. Infinitesimal changes in a, b, and a are coupled as follows: 2 sin a = b sin 0 so that 72

/ doI J2 / ~~~~~~~Disc Axis / p, —__p'_/~ M/.... da Fig. 4-2. Lumped-parameter deformation models for the vertebral column ( many-degrees-of-freedom). 75

da ~ Q cos c = sin Q db or sin Q da = ( ) db (4-1) C 0os a Also, ~2 = a2 + b2 - 2ab cos G or, 0 = 2ada + 2bdb cos @ (2adb + 2bda) (4-2) from which da -= -(b - a cos )/(a b cos @)3db Let Cab = sin Q/f cos a Cab = -(b - a cos Q)/(a - b cos 0) Then, da = Cab - db da = Cab - db (4-3) Thus, for a typical vertebra, a small displacement bb at the upper disc is accompanied by a small displacement 5a at the lower disc, and a small rigid body motion ba of the vertebra. The rotation as well as the "axial" displacements of the vertebra are resisted by the discs according to some force-deformati;on law which can be derived from the constil.tutive equat~ion of the discs and the geometry of the disc cross-section. 74

Thus, the equations of motion for the individual masses can be written and solved numerically for each small time increment during which the displacements are small. Then, a new configuration- must be computed at the end of each time increment to provide new initial conditions for the subsequent time increment. Such a computer program is being written as of this writing. 75

References to Chapter IV 1. Coermann, R. R., et al., Aerospace Medicine, Vol. 31, p. 443, 1960. 2. Naab, K. N., "A Computer Simulation of the Automobile Crash Victim, " Report No. CAL VJ-1823-R18, Cornell Aeronautical Lab., Buffalo, New York. 3. Ho, B. L. and Kalman, R. E., "Effective Construction of Linear StateVariable Models from Input/Output Functions, TRegelungstechnik, Vol. 14, 12, Dec. 1966, pp. 545-548. 4. Shinbrot, M., "On the Analysis of Linear and Nonlinear Dynamical Systems from Transient-Response Data," NACA Technical Notes, TN 3288, 1954. 5. Perdreauville, F. J. and Goodson, R. E., "Identification of Systems Described by Partial Differential Equations, " Journ. of Basic Engineering, Vol. 88, Series D, No. 2, June 1966, pp. 463-468. 6. Zienkiewicz, O. C., The Finite-Element Method, McGraw-Hill, 1967. 76

UNIVERSITY OF MICHIGAN 311 9III11111 11111111115111111111 1111 3 9015 03465 8651