OF VISCOELASTIC ANGULAR ACCELERATION IN A RIGID (KELVIN) MATERIAL ROTATIONAL SPHERICAL SHELL-A HEAD INJURY MODEL* Y. KISG
Biomrchanics
LIL-t,
K. B. CHASDRAS:
and D. U. vos ROSESBERG~
Laboratory, Tulanr University Schools of Medicine New Orleans. Louisiana 701 I?. U.S..-\.
and Engineering.
Abstract-.4 finite difference technique \vas employed to obtain the response of viscoelastic matcrial contltined in rigid spherical shells when subjected to a step angular acceleration about a diametrIca xws. Keeping the input acceleration impulse as unity. the temporal and spatial responses of Ihe \~scoclastic material were compared for Various combinations of rectangular pulse magnitude Lnd dururion. The results for 3 single sinusoidal pulse input were also compared with the infinite series solution obtained by Bbcroft (1973). The finite-difference results for the special case of an elastic material agreed ver) aell with the esnct, closrtl-fbr~t~ solution previously obtained by Liu and Chandran (1973). The primary importance of damping in the brain material is elucidated. A parametric stud) showed that slight variations of the assumed values of the nondimensional viscoelastic parameter shit&d the to!:rante curves considerably. Because of this sensitivity. due caution should be exercised in placing too much contidence in arbitrarily constructed tolerance curves pending the accurate determination of the viscoelnstic properties of the brain o\er the frequency range of interest to trauma.
ISTRODUCTIOS
Mathematical and experimental modelling of the head impact has gained special interest in biomechanics research because craniocerebral trauma accounts for about 75 per cent of all fatalities from all accidents. The Pmredings of he Cot$aw7cr 011 Htmi III~WJ (Edited by Caveness and Walker. 1966). summaiized the state-of-the-art in this area up to 1966. More recent expository and review articles by Liu (1970) and Goldsmith (1970) dealt with the progress made since then and, for the sake of brevity. the review of these papers will not be repeated here. The postulations for the head injury mechanism put forth by the investigators in this area can be broadly classified into two major categories. The first is that the rotational acceleration induced by a given impact causes high shear strains in the brain material, thus rupturing the tethering cerebral blood vessels. neo- and subcortical tissue. This mechanism was first suggested by Holbourn (1943). The second postulation is that, due to direct impact, there exist regions within the brain material where the pressure is considerably reduced. The transcapiliary pressure is suddenly increased at these locations. resulting in the bursting of the capillary. When the pressure is further reduced to near the vapor pressure of the brain substance, cavitation takes place. The subsequent catastrophic collapse of the bubbles thus formed is another possible cause of the brain damage. In this paper, we are concerned onI> with rotationally caused shear
* Rrceircd 9 October 1974. + Professor and Director. Biomechanics Laboratory. : ksistant Professor. Biomechnnics Laborator!. < Professor of Chemxal Engineering,
hence, papers dealing \vith the cavitastrains. tion school will not be revielved in the interest of brevity. Holbourn (1943) has pointed out that the bulk modulus of the brain material is about XOC00 psi. which is about the same as that of water. whereas the shear modulus is relatively small. Hence. the brain material is more susceptible to shear stresses. Ommaya et al. (1968) studied ths effect of uhiplash injury on rhesus monkeys and showed that if the head was subjected to a rotational acceleration above a threshold value, there was subdural and subarachnoid injury to the brain. Unterharnscheidt and Higgins (1969) reported similar findings from experimental rotational acceleration induced in squirrel monkeys. based on thorough and systematic neuropathology. Experiments have also been performed to determine the dynamic shear modulus and the loss tangent for the human brain. Fallenstein er al. I 1969) subjected rectangular samples of human and rhesus monkey brain to a sinusoidally varying shear stress to determine the complex shear modulus. Their measurements. however, were restricted to a resonant frequency of 10Hz and they reported the real part of the complex shear modulus in ths range 0%X-0.16 psi, while the loss tangent (the ratio of imaginary to the real part) varied from 0.40 to 0.55. ;llcElhaney et al. (1971) also performed these tests on human brain at IO Hz, but reported a value of the storage modulus varying from 0.062 to 0.03s psi and a loss modulus in the range O-051-0.057 psi with an average loss tangent value of 0.72. Shuck and Ad\-ani (1972) performed torsional tests on cylinders of brain matter over the frequency I-3jOHz. Their results show that the loss tangent varied from 0.1 to 2.5 in this frequency range and at IO Hz. their value for the shear
356
Y. li;lsG
LIV. K. B. CHASDRAN
and
D. U. vos ROSENBERG
modulus differs from the value reported by Fallenstein ef 111.(1969) by an order of magnitude. Hickling and Wsnner (I 973). in discussing a mathematical model for head impact problem. point out that the values of the complex shear modulus of interest lies in frequency range of 0-23OOHz. It is evident therefore. that further experimental determination of the complex shear modulus is needed for clarification. Noncontact predominantly rotational impacts of the head. such as whiplash, have been simulated mathematically as a rigid spherical shell containing elastic or viscoelastic material subjected to rotational acceleration about its diametrical axis by previous investigators. The dynamic response of a viscoelastic material contained in a rigid sphere subjected to harmonic translational and rotational acceleration has Fig. I. The coordinate system. been investigated by Christensen and Gottenberg (1964). The response is given in the form of complex mechanical impedance in order to determine the for the viscoelastic material as in the case for the material properties from these vibrational tests. Lee elastic material due to the complexity of the Laplace and Advani (1971) and more recently, Bycroft (1973). transform inversion. Hence, a finite difference techhave considered elastic and viscoelastjc spheres subnique was employed. jected to rotational acceleration as a model for head injury. Both these works considered the elastic probPROBLEM FORMULATIOS AND SOLUTIOS lem initially and then extended the problem to viscoelastic material through the correspondence principle. In the case of an elastic material contained in a Lee and Advani (1971) formulated the problem in the rigid spherica shell subject to an angular acceleration rotating coordinate system with the spherical material pulse about its diametrical axis, the equation of being moved with respect to the rigid spherical shell. motion in the axisymmetric spherical coordinate sysThe solution was obtained by the method of modetem (Fig. 1) can be written as superposition and the results were in the form of an infinite series. Bycroft (1973). houever. formulated the fP a2 fq-=p P2ii’- r“ (1) problem in inertial coordinates with the rigid shell 1 being subjected to a rotational acceleration and the no-slip condition at the interface in turn induced the where p is the density, p is the shear modulus, I?’ is the angular displacement and T is the time. The subsequent motion of the material contained within. initial and boundary conditions are specified as The solution was obtained. via Laplace transformation. with the inversion obtained using the standard contour integration procedure. The poles of the integrand were obtained and the residues computed with the results gisen in the form of the classical infinite series. It is well known. however. that for small values of time, the infinite series solution is, in many cases, very slowly convergent, if at all. and is therefore unwhere a is the radius of the sphere and i(T) is the suitable for numerical work. The nature of this diffiinput angular acceleration. The exact solution to the culty was dealt with in some detail by Liu (1970). ,problem specified above has been obtained by Liu Liu and Chandran (1971) pointed out that the infinite and Chandran (1973). For a linear viscoelastic series solution presented by Bycroft (1973) contained (Kelvin) material, equation (1) is modified as this defect. Furthermore. Liu and Chandran (1973) obtained the e.YRcc.closed form wave propagation a a2*6’iv- i’ solution for the case of purely elastic material con/I”2 F-G% 1 tained in a rigid spherical shell subjected to a rotational acceleration which is a general function of time. =p-,For a given time and radius. only a firlife number (3) atof terms were necessary for the numerical computation and thus the problem of convergence menwhere tioned above was completely eliminated. A rigid spherical shell containing linear viscoelastic (Kelvin) material subjected to rotational acceleration is considered as a model for head injury in this paper. and p’ is the viscosity coefficient. It is in general difficult to obtain the exact solution
r+
[
C
w -1+p[
iv
Angular accslsrstion of viscoelastic material Equation (3) can be rewritten as
0.15 1
/
_
-/
,,OJ
,’ ’ 0.4
I’ I:,.p~‘0.6
where c = (p/p)’ ’ is the elastic shear wave velocity- of the medium respectively. Equation (1) is evidently the wave equation which is hyperbolic, but for the Kelvin material, equation (1) becomes parabolic (Fliigge. 1967). Hence. there exists no characteristic velocity for this problem. The equation is nondimensionalized using the relations r = i/a. r = c;n. Lt. = &‘a, where a is the radius of the spherical shell. Following the exact solution to the elastic problem. it is assumed that Pt’(r. 19.t) = w(r. t)
sin 0.
.‘Fw!us , 08
0
= 2
(r. 6.0)
2D
3.0
4.0
Gw
Fig. 3. Isometric stress distribution as a function of radius and time for step function acceleration Input. RESLLTS
ASD
DISCC’SSIOU
Assuming a suitable value of the viscoelastic parameter E. the solution is obtained through the finite difference technique on the IBM 704-l computer. The details of the finite difference technique employed are discussed in the Appendix. A three-dimensional plot of the stress distribution as a function of radius and time was obtained with the aid of computer graphics for E = 0.122 and is shown in Fig.?. A comparison between the stress distribution for the viscoslastic and elastic case is shown in Figs. 3 and 1 for r I= I.0 and 0.2 respectively. As expected. the stress values near the center are considerably attenuated in the viscoelastic case as compared to the elastic case. T/W
(5)
The nondimensional governing equation for the viscoelastic material, along with the boundary and initial conditions are given below:
w(r. 6,O)
LO
= 0
higher
stresses
are
alrva~~
obrairled
at
r/:z srrtji~c~
the sphere.
q/
This was found to be generalI> true when computations were made Lvith various values of E z(t), ranging from 0.001 to 10.0. It can be seen from Fig.4 that the effect of damping is rather dramatic near where E = $/pea is the viscoelastic parameter. the center of the sphere. In the case of the elastic A finite difference scheme is employed to obtain material. where there is no damping of the energy the solution for the governing equation and the imparted to the sphere. the fluctuation of the strsss details of the scheme are included in the Appendix. is considerable at nondimensional time near I and The result is checked as follows. In equation (6) if 3. These fluctuations are the result of the superposiwe let E = 0, the equation reduces to the second order tion of waves travelling from the surface to the center wave equation. i.e. the case of purely elastic medium and vice versa. Thus. one can easily visualize the contained in rigid spherical shells. for which the tx~~cr, peaks in the response to be where the standing waves cIosvd fbrm solution has been obtained by Liu and add and valleys to be where they subtract. However. Chandran (1973). The finite difference solution for this since the brain substance has a low shear modulus. tax was in excellent agreement with the enact solurelatively low values of angular accelsrarion could tion. result in high and, perhaps traumatic. shear stress,
$(l.e.t) =
04
,“\ ,r---. 02;
,,’
Elastic / ‘\
‘\
g 5
mcterd
/
r=1.0
I
‘.
--___’
/
*-
‘\ \
/’
\
‘\
.=_____ - ’
, /
/’
--\ \
1
\
\ 4! ‘i /
v6cc.?klste materidE=o.i22
02
0.0
, I.0
I 2.0
3.0
4.0
Time
Fig. 3. Comparison of maximum shear stress between rlastic and viscorlastic material. I’ = 1.1).
Y.
KISG LIU. K. B.
CHASDRAN and
ROSESBERG
D. L‘. vos
/
04‘
r=0.2 0.2 c
q G=l -0.2
L
-0.4 0.0
c El& I\-
i \ i ‘, i 1 I \
v
htenol
, / I ,
_
II / I / : I / /
---, -- ____---.-/ ~scoebstc Mokrlol c =0.122 ; : :
1.0
2D
30
CO
Tim Fig. 4. Comparison of maximum shear stress between elastic and viscoelastic material. r = 0.2. i 0.4-
0.3 -
rz1.0
GO.2 &
-_0.1 1 .Ol
0.1
E
IO
10.0
Fig. 5. Maximum shear stress as a function of E for r = I.0 and r = 0.2. except for damping. The maximum shear stress for various values ofe is plotted in Fig. 5 and it illustrates the decreasing values of the maximum shear stress at I’= 0.2 and r = 1.0 for increasing values of E. It is well known that not only is the magnitude of the input acceleration. but the duration of acceleration is also an important consideration in studying the head injury levels. A meaningful comparison can be made only if we subject the spherical she11 to a unit impulse. This result can be achieved by replacing the step function accelerationflt) by ACf(t) -fit - T)], where .4 and r are the ampIitude and duration of the rectangular acceleration pulse respectively. Stated another way, the area under the input impulse curve is always maintained at unity by appropriately reducing the amplitude when the duration of acceleration is correspondingly increased. The unit impuIse corresponds to an actual value of about 0.02 rad/sec for the values of the parameters assumed in this case. The radial distribution of the maximum value of stress with the ratio r/r as a parameter is plotted in Fig. 6. A value of t equal to 1 represents theoreticalIy the time required for the elastic stress wave to travel from the surface of the sphere to the center. Hence. ti’r represents the ratio of the input pulse duration to the transit time for the stress wave. A similar comparison was made for the case of the elastic material by Liu and Chandran (1973). where it was shown that for small values of r/t. the masimum shear stress was obtained in the central region of the sphere and that as the I r ratio was increased, the maximum value
shifted towards the surface of the sphere. It was further shown that such a shift occurred between values of r/s of 0.1 and 1.0. In contrast, in the present case of the viscoelastic material, even for small values
0.1
02
a4
0.6
OE
I.0
mdius
Fig. 6. Radial distribution of maximum shear stress for unit impulse with c/r as parameter. Numbers in parenthesis represent the dimensionless time at which the maximum stress is obtained.
;\ngular acceleration of viscoclastic material of r :. the matlmum shear stress always occurs at the surface of the sphere. This is to be expected since in ths \iscoslastic case. the stress is greatly attenuated before It reaches the core region. In Fig. 6. the number in parenthesis represents the nondimensional time at ivhich the maximum value of the shear stress is obtained. Again. in comparison with the elastic case. the maximum value is obtained considerably earlier (at a nondimensional time of about 0.1 compared to 2.0 for the elastic case). Figure 7 gives a distribution of maximum shear stress as a function of t’T with the nondimensional radius as the parameter. Bycroft (1973) used a single sinusoidal pulse acceleration input in his analysis for the viscoelastic model and constructed theoretical tolerance curves for concussion. By corrslating with experimental whiplash studies. he chose an arbitrary value of shear strain equal to 0.05 at a nondimensional radius of 0.3 as a threshold value for concussion. He compared the theoretically obtained tolerance curves for concussion with the experimentally obtained threshold of concussion and later on extended the tolerance curse to the human. using a value of E equal to 0.177. based on the experiments by Shuck and Advani (1972). From this tolerance curve for humans. he predicts that angular acceleration of 3500 rad,‘sec’ or greater. will be required to produce concussion in man. Using experimental whiplash on primates and with the aid of Holbourn’s (1943) scaling “law”, Ommaya and Hirsch (1971) predicted a threshold level for concussion of lSO0 rad see’. Both the modulus and loss tangent values which were ass~n& for the theoretical computation are functions of frequencies which can only be determined by experiments. and as indicated earlier. more accurate measurements are needed to determine these values. Thus. one must exercise due
0.15
B
0.10
005
0 0
02
0.4
0.6
a.3
1.0
Fig. S. Comparison of the numerical results of the present finite difference and infinite szricr solution of B>croft (1973) forE = 0.122.
caution in placing too much confidence in theoretically constructed tolerance curves for head injury because even slight variations in the assumed values for the shear modulus and the loss tangent will shift the tolerance curve considerably. In this work. the finite difference solution for a sinusoidal pulse acceleration input similar to that posed by Bycroft (1973) was obtained. In Fig. S. his infinite series solution is compared with the finite diiference solution and it is seen that the finite difference solution gives higher values of shear strain at the surJace of the sphere. It is well known that the classical infinite series is very slowly convergent. if at all, for small values of time, which are the quantities of interest in this problem. As was pointed out by Liu and Chandran 1197-l). this &as indeed the cast’ for B!croft’s solution. The present finite difference solution. on the other hand, agrees very well with the exact solution and hence more confidence can be placed upon it in as much as it has satisfied the WC~SSCII-I conditiorl for accuiacy. The theoretical tolerance curves obtained by Bycroft (1973) superposed on the experimental threshold values for concussion on primates are compared with the ones obtained by the finite difference technique in Figs. 9 and 10. it IS evident_ especially from Fig. 10. that the curve obtained
‘i
Fig. 7. >laGmum shear stress as a function of duration of acceleration for unit impulse with radius as parameter. Numbers in parenthesis represent the dimensionless time at which the maximum stress is obtained.
!
/I/l,I,I,I,l,/,I
2
Fig. 9. Threshold
4
6 8 IO I2 pulseanotion-msec
of concussion:
rhesus
14
I6
monks!. E = O.j06.
290
Lrti. K. B. CHASDRU
Y. KING
\
s
ocmusslvc 0 Non ealmsm”C
\
\ ‘0
k!!!!b ‘L__
--__ -_ ---_-
___
0
I
I
0
Fi g.
-Flntc
‘\
j2 -t
Series blutla, oiffererce S&Awl
-----Infink
‘0
$3 G 2
I
I
2
4
10. Threshold
of
0
I
I
I
6 6 IO Puke Cumiiin-msaz
concussion: E = 0525.
squirrel
I2
monkey.
from the present work yields a more conservative value for the threshold for concussion in primates: A tolerance curve is also constructed for the viscoelastic material with E = 0121 (value used by Bycroft for humans) in Fig. II. with the arbitrary criterion that a shear strain of 005 at a nondimensional radius of 0.3 will result in concussion. This tolerance curve is compared with the one drawn for the elastic material, E = 0. to illustrate the effect of damping in increasing the threshold value of angular rotation which can be tolerated by the brain material. In Fig. 12, the theoretical tolerance curves constructed for various values of E are compared to indicate the importance of accurate determination of the viscoelastic properties of the IC
and D. U. VOXROSENBERG
1’
I
I
r
I
I
0
Fig. I?.
Tolerance curves with E as a parameter.
human brain over the appropriate frequency range (IS3 200 Hz) before a confident prediction of the threshold level for rotational damage is possible.
(1) The damping of the brain material is of critical importance in delineating the shear stresses induced by a rotational acceleration input. (2) The model of a closed spherical shell with a Kelvin viscoelastic core predicts the possibility of subdural and subarachnoid injuries without subcortical involvement but the reverse is impossible. This agrees with the neuropathological study of Unterharnscheidt and Higgins (1969). (3) To consistently explain the phenomena of a subcortical injury without subdural and subarachnoid involvement would probably require the inclusion of the brain stem as well as the head-neck dynamics in the model. (4) Slight variations in the assumed values of the shear modulus and loss tansent of the brain material shifts the so-called tolerance curves considerably. (5) Due caution should be exercised in placing too much confidence in arbitrarily constructed tolerance curves pending the accurate determination of the viscoelastic properties of the brain over the appropriate frequency range. Ackrrowtedgemellts-This work was completed while rhe first author was an NIH Research Career Development Awardee (Grant No. GM -!0723-02). The support of the research reported above by the NIH (Grant No. GM 19107-01) and NSF (Grant So. GK 32M7) is gratefully acknowledged. REFEREYCES
Bycrok G. N. (1973) Mathematical model of a head subjected to an angular acceleration. J. Bionwchanics 6. 487-495.
Fig. I I. Comparison
of tolerance curve between elastic and viscoelastic material.
Caveness. W. F. and Walker. A. E. (Editors) (1966) Proc. Head Injury Co$ (University of Chicago). Lippincott. Philadelphia. PA. Christensen. R. M. and Gorrenberg. W. G. (1964) The dynamic response of a solid. viscoelastic sphere to translational and rotational excitation. J. appl. .\~ec/z. 31. 27% 275.
.Angular acceleration Fallenstein. Dynamrc
G. T.. Hulce. V. D. and M&in. J. W. (1969) mechanical properties of human brain tissue. J Biornrchnics 2. Z 17- 2’6. Fltiggc. W. ( 1967) I’iscodusticity. Blaisdell. M.altham. ?.lass. Goldsmith. W. (1970)Biomechanics of head injury. Biomrchanics: Its foundations ad Ohjecriws (Edited by Fung. Y C.. Perrone. S. and -\nliker. M.). pp. SYC6Z-l. Prentics-Hall. New Jersey. R. and Wenner. M. L. (1973) !+lathematical Hickling. model of a head subjected to an avisymmerric impact. J. Biornrchics 6, I I C 132. Hoibourn. A. H. S. (19-I:) Mechanics of head injuries. f.nrrctt 2. -t3,Y-&tl. Lee. Y. C. and Advani, S. H. (1970) Transient response of a sphere to symmetrical torsional loading-a head injury model. Jfnrh. Biosci. 6, 472-486. Liu. Y. King (1970) The biomechanics of spinal and head impact: problems of mathematical simulation. hoc. Swrp. Biorl.w. .Motld/hy ad Its .?.pplicarions. pp. 7017%. Dayton. Ohio. Liu. Y. King and Chandran. K. B. (1973) The exact soluuon to the rotational acceleration of rigid spherical shells filled with a brain-like elastic material. Proc. IOrh .&mii-erw-y .\/ecci~~g qf‘ rhr Sot. o/ Enqq Sci. Raleigh. S.C. Liu. Y. King and Chandran. K. B. (197-L) Comment on paper: Mathematical model of a head subjected to an angular acceleration. bv G. N. Bycroft. J. Biomechanics 7. 319-321. McElhaney. J. H.. bvlelvin. J. W.. Roberts. 1’. L. and Portnay. H. D. (1972) Dynamic characteristics oi the tissue of the head. .S.t~rrp. Persprcrires in Biomrti. Err~gmg (The Univ. of Strathclyde). Macmillan. London. Ommaja. A. K.. Fax. F. and Yarnell. P. (196s) LVhiplash injury and brain damage. J. .1m. Med. .-lssoc. 204. 235. Ommaya. X. K. and Hirsch. A. E. (1971) Tolerance for cerebral concussion from head impact and whiplash in primates. J. Biomrchanics 1. 13-Z. Shuck. L. F. and .Advani. S. H. (1972) Rheological response of human brain tissue. Presented at ASXlE Annual IVinter Meeting. New York. J. Busic Enguq ASME T,aus. 94, Series D. Paper No. 72 WA BHF-2. Lnterharnscheidt. F. and Higgins. L. S. (1969) Traumatic lesions of brain and spinal cord due to nondeforming angular acceleration of the head. LSuir. Terns Reports BIO/. .Lfrri. 27( 1). van Rosenberg. D. U. (1969) .Lfethotis /br rite .~~rmerical Sohrtio~ 01 Partial Diferetxial Eqwuions. Elsevier. New E’ork.
of viscoelastic
291
maternal viscoetastic parameter shear modulus. lb i in’ viscosity coefficient. Ibm set-in densitv.. Ibm’in’ , shear stress frequency of sine pulse acceleration
E( = p’ pad p p’ P :0
‘at the top indicates dimensional variables. variables have been nondimensionalized.
input
IVithout
: the
APPESDIT Xumrricai
solution
~iriscoelastfc
The differential
=z
s
equation
)
i
h
r-x.
Two new variables
material
ecpiutior~s
can be put in the form
(A-1)
i
are defined as 9
s=rzG” Sr and
The original
equation
then becomes
and the two new variables
are related
sv ss -_=--f-r.
2
Sr
r
?t
by (A-5)
Equations (A-3-A-5) comprise a set of equations in which the highest order derivative in either independent variable is of the first order. The initial conditions for the three variables are ru=v=s=O The boundary
forallratr=O.
conditions
w = r = 0
for all
at r = 0 are
t. and at r = I
SC - = 1 for a unit step function.
c’t
ll
radius of the spherical shell. in. amplitude ofexternal acceleration shear wave velocity. in:sec auisymmetric spherical coordinates time variable displacement in B direction external acceleration input
.-I C
i-. 0
t. : II:&. l(f)
This set of equations is solved by the centered-difference method described by van Rosenberg (1969). A grid of points is set up in the radial direction. and values of the dependent variables are computed at these points at set intervals of time. A set of points. with the indexing nomenclature. is shown in Fig. 13. For conv-enience in using the solution algorithm, the space index is given one set of values for sand another set for r and w as shown.
(Fig. I)
S
I
2
L-1
L‘.H
1
?
i
L-l
i
i+
1
Fig. I?. Grid of points.
L
L Lfl
292
Y.
k&c
LK. K. B.
CHANDRAS
The finite difference analog to equation (A-31 IS centered half-way bctwscn two consecutive time levels and is
This equation is second order correct in C. The finite difference analogs to the terms in equations (A-4 and A-S) are centered in time and space. The derivati\es of c are replaced by
and D. Li. VONROSESBERG When these finite difference analogs are substituted into the original partial di&rential equations. (A-4 and A-j). two sets of finite difference equations in the unknowns si,“_ / and c,,,_ , result. These equations are
[T]s
I_,. n-i + [I +(l
;,I”“;;;;):)]%~
1
At - 7~ = [s,,, - xi_ ,..I + [ - IAr bc,., + wi_ ,,.I [ Ar
and
1
i. cc -4 ! C,-I.“-1 - CL-l.” c ~i.n-l - Ci.” i Sr 1i++.“-* 1 At A’ (A-7)
+
1 _ At[(Ar - 2~) - x(At + ZE)] 2(i - Z)‘(Ar)’
[
At[(At + 2~) - r(Ar - ZE)]
The analogs for s are similar. and all are second order correct in space and time. The analog to the last term in equation (A-j) is
Z(i - l)‘(Ar)’
1CL” -I.n 1c,
(A-ila)
and
si- I.“_ + [
1
&
1.(A-S)
[(I - 4L.“-l A (’ + x’r:.J
s,1.” _, + si-
+[&(_ Ar
The parameter. x. can be varied to form any of three analogs; namely. the four-point analog for x = 0. the forward diagonal analog for x = I, and the backward diagonal analog when r = -I. For the boundary conditions of this problem. the backward diagonal (z = - I) is used Tar t 2 I, after which time the four-point analog (2 = 01 is used. The analogs to the terms Zc,r’ and ZW in equation (A-5) are handled in a similar manner. The values of w at the r._, level are eliminated from the final equation by equation (A-3a). All these analogs are also second order correct in space and in time. An analog, which is also second order correct in space and time. for the cross-derivative of s is given by 2
I +
,.”
-
Si.” -
(Ar)(Ar)
si-
I.“-
1
(A-9)
+[-$(I
[$(I+;)]b+, +
Se..-!
, , I, + “)lC,_, I-I
+~)],..+~(I
“i, = s
I.”
I s,_,n
-SJ]t++,.“. (A-5a)
These equations are written for the L-l values of i from i = Z to i = L. Thus. there are 2(L-I) equations in the set. All values at time f, are known either from the initial conditions or from the computations of the previous time step. Thus, there are ZL unknowns in the equations. However, two of these values are given by the boundary conditions, and the equations can be solved for the dependent variables s and c. An efficient algorithm for this solution has been given by van Rosenberg (1969).
Biomrchanics
LIL-t,
K. B. CHASDRAS:
and D. U. vos ROSESBERG~
Laboratory, Tulanr University Schools of Medicine New Orleans. Louisiana 701 I?. U.S..-\.
and Engineering.
Abstract-.4 finite difference technique \vas employed to obtain the response of viscoelastic matcrial contltined in rigid spherical shells when subjected to a step angular acceleration about a diametrIca xws. Keeping the input acceleration impulse as unity. the temporal and spatial responses of Ihe \~scoclastic material were compared for Various combinations of rectangular pulse magnitude Lnd dururion. The results for 3 single sinusoidal pulse input were also compared with the infinite series solution obtained by Bbcroft (1973). The finite-difference results for the special case of an elastic material agreed ver) aell with the esnct, closrtl-fbr~t~ solution previously obtained by Liu and Chandran (1973). The primary importance of damping in the brain material is elucidated. A parametric stud) showed that slight variations of the assumed values of the nondimensional viscoelastic parameter shit&d the to!:rante curves considerably. Because of this sensitivity. due caution should be exercised in placing too much contidence in arbitrarily constructed tolerance curves pending the accurate determination of the viscoelnstic properties of the brain o\er the frequency range of interest to trauma.
ISTRODUCTIOS
Mathematical and experimental modelling of the head impact has gained special interest in biomechanics research because craniocerebral trauma accounts for about 75 per cent of all fatalities from all accidents. The Pmredings of he Cot$aw7cr 011 Htmi III~WJ (Edited by Caveness and Walker. 1966). summaiized the state-of-the-art in this area up to 1966. More recent expository and review articles by Liu (1970) and Goldsmith (1970) dealt with the progress made since then and, for the sake of brevity. the review of these papers will not be repeated here. The postulations for the head injury mechanism put forth by the investigators in this area can be broadly classified into two major categories. The first is that the rotational acceleration induced by a given impact causes high shear strains in the brain material, thus rupturing the tethering cerebral blood vessels. neo- and subcortical tissue. This mechanism was first suggested by Holbourn (1943). The second postulation is that, due to direct impact, there exist regions within the brain material where the pressure is considerably reduced. The transcapiliary pressure is suddenly increased at these locations. resulting in the bursting of the capillary. When the pressure is further reduced to near the vapor pressure of the brain substance, cavitation takes place. The subsequent catastrophic collapse of the bubbles thus formed is another possible cause of the brain damage. In this paper, we are concerned onI> with rotationally caused shear
* Rrceircd 9 October 1974. + Professor and Director. Biomechanics Laboratory. : ksistant Professor. Biomechnnics Laborator!. < Professor of Chemxal Engineering,
hence, papers dealing \vith the cavitastrains. tion school will not be revielved in the interest of brevity. Holbourn (1943) has pointed out that the bulk modulus of the brain material is about XOC00 psi. which is about the same as that of water. whereas the shear modulus is relatively small. Hence. the brain material is more susceptible to shear stresses. Ommaya et al. (1968) studied ths effect of uhiplash injury on rhesus monkeys and showed that if the head was subjected to a rotational acceleration above a threshold value, there was subdural and subarachnoid injury to the brain. Unterharnscheidt and Higgins (1969) reported similar findings from experimental rotational acceleration induced in squirrel monkeys. based on thorough and systematic neuropathology. Experiments have also been performed to determine the dynamic shear modulus and the loss tangent for the human brain. Fallenstein er al. I 1969) subjected rectangular samples of human and rhesus monkey brain to a sinusoidally varying shear stress to determine the complex shear modulus. Their measurements. however, were restricted to a resonant frequency of 10Hz and they reported the real part of the complex shear modulus in ths range 0%X-0.16 psi, while the loss tangent (the ratio of imaginary to the real part) varied from 0.40 to 0.55. ;llcElhaney et al. (1971) also performed these tests on human brain at IO Hz, but reported a value of the storage modulus varying from 0.062 to 0.03s psi and a loss modulus in the range O-051-0.057 psi with an average loss tangent value of 0.72. Shuck and Ad\-ani (1972) performed torsional tests on cylinders of brain matter over the frequency I-3jOHz. Their results show that the loss tangent varied from 0.1 to 2.5 in this frequency range and at IO Hz. their value for the shear
356
Y. li;lsG
LIV. K. B. CHASDRAN
and
D. U. vos ROSENBERG
modulus differs from the value reported by Fallenstein ef 111.(1969) by an order of magnitude. Hickling and Wsnner (I 973). in discussing a mathematical model for head impact problem. point out that the values of the complex shear modulus of interest lies in frequency range of 0-23OOHz. It is evident therefore. that further experimental determination of the complex shear modulus is needed for clarification. Noncontact predominantly rotational impacts of the head. such as whiplash, have been simulated mathematically as a rigid spherical shell containing elastic or viscoelastic material subjected to rotational acceleration about its diametrical axis by previous investigators. The dynamic response of a viscoelastic material contained in a rigid sphere subjected to harmonic translational and rotational acceleration has Fig. I. The coordinate system. been investigated by Christensen and Gottenberg (1964). The response is given in the form of complex mechanical impedance in order to determine the for the viscoelastic material as in the case for the material properties from these vibrational tests. Lee elastic material due to the complexity of the Laplace and Advani (1971) and more recently, Bycroft (1973). transform inversion. Hence, a finite difference techhave considered elastic and viscoelastjc spheres subnique was employed. jected to rotational acceleration as a model for head injury. Both these works considered the elastic probPROBLEM FORMULATIOS AND SOLUTIOS lem initially and then extended the problem to viscoelastic material through the correspondence principle. In the case of an elastic material contained in a Lee and Advani (1971) formulated the problem in the rigid spherica shell subject to an angular acceleration rotating coordinate system with the spherical material pulse about its diametrical axis, the equation of being moved with respect to the rigid spherical shell. motion in the axisymmetric spherical coordinate sysThe solution was obtained by the method of modetem (Fig. 1) can be written as superposition and the results were in the form of an infinite series. Bycroft (1973). houever. formulated the fP a2 fq-=p P2ii’- r“ (1) problem in inertial coordinates with the rigid shell 1 being subjected to a rotational acceleration and the no-slip condition at the interface in turn induced the where p is the density, p is the shear modulus, I?’ is the angular displacement and T is the time. The subsequent motion of the material contained within. initial and boundary conditions are specified as The solution was obtained. via Laplace transformation. with the inversion obtained using the standard contour integration procedure. The poles of the integrand were obtained and the residues computed with the results gisen in the form of the classical infinite series. It is well known. however. that for small values of time, the infinite series solution is, in many cases, very slowly convergent, if at all. and is therefore unwhere a is the radius of the sphere and i(T) is the suitable for numerical work. The nature of this diffiinput angular acceleration. The exact solution to the culty was dealt with in some detail by Liu (1970). ,problem specified above has been obtained by Liu Liu and Chandran (1971) pointed out that the infinite and Chandran (1973). For a linear viscoelastic series solution presented by Bycroft (1973) contained (Kelvin) material, equation (1) is modified as this defect. Furthermore. Liu and Chandran (1973) obtained the e.YRcc.closed form wave propagation a a2*6’iv- i’ solution for the case of purely elastic material con/I”2 F-G% 1 tained in a rigid spherical shell subjected to a rotational acceleration which is a general function of time. =p-,For a given time and radius. only a firlife number (3) atof terms were necessary for the numerical computation and thus the problem of convergence menwhere tioned above was completely eliminated. A rigid spherical shell containing linear viscoelastic (Kelvin) material subjected to rotational acceleration is considered as a model for head injury in this paper. and p’ is the viscosity coefficient. It is in general difficult to obtain the exact solution
r+
[
C
w -1+p[
iv
Angular accslsrstion of viscoelastic material Equation (3) can be rewritten as
0.15 1
/
_
-/
,,OJ
,’ ’ 0.4
I’ I:,.p~‘0.6
where c = (p/p)’ ’ is the elastic shear wave velocity- of the medium respectively. Equation (1) is evidently the wave equation which is hyperbolic, but for the Kelvin material, equation (1) becomes parabolic (Fliigge. 1967). Hence. there exists no characteristic velocity for this problem. The equation is nondimensionalized using the relations r = i/a. r = c;n. Lt. = &‘a, where a is the radius of the spherical shell. Following the exact solution to the elastic problem. it is assumed that Pt’(r. 19.t) = w(r. t)
sin 0.
.‘Fw!us , 08
0
= 2
(r. 6.0)
2D
3.0
4.0
Gw
Fig. 3. Isometric stress distribution as a function of radius and time for step function acceleration Input. RESLLTS
ASD
DISCC’SSIOU
Assuming a suitable value of the viscoelastic parameter E. the solution is obtained through the finite difference technique on the IBM 704-l computer. The details of the finite difference technique employed are discussed in the Appendix. A three-dimensional plot of the stress distribution as a function of radius and time was obtained with the aid of computer graphics for E = 0.122 and is shown in Fig.?. A comparison between the stress distribution for the viscoslastic and elastic case is shown in Figs. 3 and 1 for r I= I.0 and 0.2 respectively. As expected. the stress values near the center are considerably attenuated in the viscoelastic case as compared to the elastic case. T/W
(5)
The nondimensional governing equation for the viscoelastic material, along with the boundary and initial conditions are given below:
w(r. 6,O)
LO
= 0
higher
stresses
are
alrva~~
obrairled
at
r/:z srrtji~c~
the sphere.
q/
This was found to be generalI> true when computations were made Lvith various values of E z(t), ranging from 0.001 to 10.0. It can be seen from Fig.4 that the effect of damping is rather dramatic near where E = $/pea is the viscoelastic parameter. the center of the sphere. In the case of the elastic A finite difference scheme is employed to obtain material. where there is no damping of the energy the solution for the governing equation and the imparted to the sphere. the fluctuation of the strsss details of the scheme are included in the Appendix. is considerable at nondimensional time near I and The result is checked as follows. In equation (6) if 3. These fluctuations are the result of the superposiwe let E = 0, the equation reduces to the second order tion of waves travelling from the surface to the center wave equation. i.e. the case of purely elastic medium and vice versa. Thus. one can easily visualize the contained in rigid spherical shells. for which the tx~~cr, peaks in the response to be where the standing waves cIosvd fbrm solution has been obtained by Liu and add and valleys to be where they subtract. However. Chandran (1973). The finite difference solution for this since the brain substance has a low shear modulus. tax was in excellent agreement with the enact solurelatively low values of angular accelsrarion could tion. result in high and, perhaps traumatic. shear stress,
$(l.e.t) =
04
,“\ ,r---. 02;
,,’
Elastic / ‘\
‘\
g 5
mcterd
/
r=1.0
I
‘.
--___’
/
*-
‘\ \
/’
\
‘\
.=_____ - ’
, /
/’
--\ \
1
\
\ 4! ‘i /
v6cc.?klste materidE=o.i22
02
0.0
, I.0
I 2.0
3.0
4.0
Time
Fig. 3. Comparison of maximum shear stress between rlastic and viscorlastic material. I’ = 1.1).
Y.
KISG LIU. K. B.
CHASDRAN and
ROSESBERG
D. L‘. vos
/
04‘
r=0.2 0.2 c
q G=l -0.2
L
-0.4 0.0
c El& I\-
i \ i ‘, i 1 I \
v
htenol
, / I ,
_
II / I / : I / /
---, -- ____---.-/ ~scoebstc Mokrlol c =0.122 ; : :
1.0
2D
30
CO
Tim Fig. 4. Comparison of maximum shear stress between elastic and viscoelastic material. r = 0.2. i 0.4-
0.3 -
rz1.0
GO.2 &
-_0.1 1 .Ol
0.1
E
IO
10.0
Fig. 5. Maximum shear stress as a function of E for r = I.0 and r = 0.2. except for damping. The maximum shear stress for various values ofe is plotted in Fig. 5 and it illustrates the decreasing values of the maximum shear stress at I’= 0.2 and r = 1.0 for increasing values of E. It is well known that not only is the magnitude of the input acceleration. but the duration of acceleration is also an important consideration in studying the head injury levels. A meaningful comparison can be made only if we subject the spherical she11 to a unit impulse. This result can be achieved by replacing the step function accelerationflt) by ACf(t) -fit - T)], where .4 and r are the ampIitude and duration of the rectangular acceleration pulse respectively. Stated another way, the area under the input impulse curve is always maintained at unity by appropriately reducing the amplitude when the duration of acceleration is correspondingly increased. The unit impuIse corresponds to an actual value of about 0.02 rad/sec for the values of the parameters assumed in this case. The radial distribution of the maximum value of stress with the ratio r/r as a parameter is plotted in Fig. 6. A value of t equal to 1 represents theoreticalIy the time required for the elastic stress wave to travel from the surface of the sphere to the center. Hence. ti’r represents the ratio of the input pulse duration to the transit time for the stress wave. A similar comparison was made for the case of the elastic material by Liu and Chandran (1973). where it was shown that for small values of r/t. the masimum shear stress was obtained in the central region of the sphere and that as the I r ratio was increased, the maximum value
shifted towards the surface of the sphere. It was further shown that such a shift occurred between values of r/s of 0.1 and 1.0. In contrast, in the present case of the viscoelastic material, even for small values
0.1
02
a4
0.6
OE
I.0
mdius
Fig. 6. Radial distribution of maximum shear stress for unit impulse with c/r as parameter. Numbers in parenthesis represent the dimensionless time at which the maximum stress is obtained.
;\ngular acceleration of viscoclastic material of r :. the matlmum shear stress always occurs at the surface of the sphere. This is to be expected since in ths \iscoslastic case. the stress is greatly attenuated before It reaches the core region. In Fig. 6. the number in parenthesis represents the nondimensional time at ivhich the maximum value of the shear stress is obtained. Again. in comparison with the elastic case. the maximum value is obtained considerably earlier (at a nondimensional time of about 0.1 compared to 2.0 for the elastic case). Figure 7 gives a distribution of maximum shear stress as a function of t’T with the nondimensional radius as the parameter. Bycroft (1973) used a single sinusoidal pulse acceleration input in his analysis for the viscoelastic model and constructed theoretical tolerance curves for concussion. By corrslating with experimental whiplash studies. he chose an arbitrary value of shear strain equal to 0.05 at a nondimensional radius of 0.3 as a threshold value for concussion. He compared the theoretically obtained tolerance curves for concussion with the experimentally obtained threshold of concussion and later on extended the tolerance curse to the human. using a value of E equal to 0.177. based on the experiments by Shuck and Advani (1972). From this tolerance curve for humans. he predicts that angular acceleration of 3500 rad,‘sec’ or greater. will be required to produce concussion in man. Using experimental whiplash on primates and with the aid of Holbourn’s (1943) scaling “law”, Ommaya and Hirsch (1971) predicted a threshold level for concussion of lSO0 rad see’. Both the modulus and loss tangent values which were ass~n& for the theoretical computation are functions of frequencies which can only be determined by experiments. and as indicated earlier. more accurate measurements are needed to determine these values. Thus. one must exercise due
0.15
B
0.10
005
0 0
02
0.4
0.6
a.3
1.0
Fig. S. Comparison of the numerical results of the present finite difference and infinite szricr solution of B>croft (1973) forE = 0.122.
caution in placing too much confidence in theoretically constructed tolerance curves for head injury because even slight variations in the assumed values for the shear modulus and the loss tangent will shift the tolerance curve considerably. In this work. the finite difference solution for a sinusoidal pulse acceleration input similar to that posed by Bycroft (1973) was obtained. In Fig. S. his infinite series solution is compared with the finite diiference solution and it is seen that the finite difference solution gives higher values of shear strain at the surJace of the sphere. It is well known that the classical infinite series is very slowly convergent. if at all, for small values of time, which are the quantities of interest in this problem. As was pointed out by Liu and Chandran 1197-l). this &as indeed the cast’ for B!croft’s solution. The present finite difference solution. on the other hand, agrees very well with the exact solution and hence more confidence can be placed upon it in as much as it has satisfied the WC~SSCII-I conditiorl for accuiacy. The theoretical tolerance curves obtained by Bycroft (1973) superposed on the experimental threshold values for concussion on primates are compared with the ones obtained by the finite difference technique in Figs. 9 and 10. it IS evident_ especially from Fig. 10. that the curve obtained
‘i
Fig. 7. >laGmum shear stress as a function of duration of acceleration for unit impulse with radius as parameter. Numbers in parenthesis represent the dimensionless time at which the maximum stress is obtained.
!
/I/l,I,I,I,l,/,I
2
Fig. 9. Threshold
4
6 8 IO I2 pulseanotion-msec
of concussion:
rhesus
14
I6
monks!. E = O.j06.
290
Lrti. K. B. CHASDRU
Y. KING
\
s
ocmusslvc 0 Non ealmsm”C
\
\ ‘0
k!!!!b ‘L__
--__ -_ ---_-
___
0
I
I
0
Fi g.
-Flntc
‘\
j2 -t
Series blutla, oiffererce S&Awl
-----Infink
‘0
$3 G 2
I
I
2
4
10. Threshold
of
0
I
I
I
6 6 IO Puke Cumiiin-msaz
concussion: E = 0525.
squirrel
I2
monkey.
from the present work yields a more conservative value for the threshold for concussion in primates: A tolerance curve is also constructed for the viscoelastic material with E = 0121 (value used by Bycroft for humans) in Fig. II. with the arbitrary criterion that a shear strain of 005 at a nondimensional radius of 0.3 will result in concussion. This tolerance curve is compared with the one drawn for the elastic material, E = 0. to illustrate the effect of damping in increasing the threshold value of angular rotation which can be tolerated by the brain material. In Fig. 12, the theoretical tolerance curves constructed for various values of E are compared to indicate the importance of accurate determination of the viscoelastic properties of the IC
and D. U. VOXROSENBERG
1’
I
I
r
I
I
0
Fig. I?.
Tolerance curves with E as a parameter.
human brain over the appropriate frequency range (IS3 200 Hz) before a confident prediction of the threshold level for rotational damage is possible.
(1) The damping of the brain material is of critical importance in delineating the shear stresses induced by a rotational acceleration input. (2) The model of a closed spherical shell with a Kelvin viscoelastic core predicts the possibility of subdural and subarachnoid injuries without subcortical involvement but the reverse is impossible. This agrees with the neuropathological study of Unterharnscheidt and Higgins (1969). (3) To consistently explain the phenomena of a subcortical injury without subdural and subarachnoid involvement would probably require the inclusion of the brain stem as well as the head-neck dynamics in the model. (4) Slight variations in the assumed values of the shear modulus and loss tansent of the brain material shifts the so-called tolerance curves considerably. (5) Due caution should be exercised in placing too much confidence in arbitrarily constructed tolerance curves pending the accurate determination of the viscoelastic properties of the brain over the appropriate frequency range. Ackrrowtedgemellts-This work was completed while rhe first author was an NIH Research Career Development Awardee (Grant No. GM -!0723-02). The support of the research reported above by the NIH (Grant No. GM 19107-01) and NSF (Grant So. GK 32M7) is gratefully acknowledged. REFEREYCES
Bycrok G. N. (1973) Mathematical model of a head subjected to an angular acceleration. J. Bionwchanics 6. 487-495.
Fig. I I. Comparison
of tolerance curve between elastic and viscoelastic material.
Caveness. W. F. and Walker. A. E. (Editors) (1966) Proc. Head Injury Co$ (University of Chicago). Lippincott. Philadelphia. PA. Christensen. R. M. and Gorrenberg. W. G. (1964) The dynamic response of a solid. viscoelastic sphere to translational and rotational excitation. J. appl. .\~ec/z. 31. 27% 275.
.Angular acceleration Fallenstein. Dynamrc
G. T.. Hulce. V. D. and M&in. J. W. (1969) mechanical properties of human brain tissue. J Biornrchnics 2. Z 17- 2’6. Fltiggc. W. ( 1967) I’iscodusticity. Blaisdell. M.altham. ?.lass. Goldsmith. W. (1970)Biomechanics of head injury. Biomrchanics: Its foundations ad Ohjecriws (Edited by Fung. Y C.. Perrone. S. and -\nliker. M.). pp. SYC6Z-l. Prentics-Hall. New Jersey. R. and Wenner. M. L. (1973) !+lathematical Hickling. model of a head subjected to an avisymmerric impact. J. Biornrchics 6, I I C 132. Hoibourn. A. H. S. (19-I:) Mechanics of head injuries. f.nrrctt 2. -t3,Y-&tl. Lee. Y. C. and Advani, S. H. (1970) Transient response of a sphere to symmetrical torsional loading-a head injury model. Jfnrh. Biosci. 6, 472-486. Liu. Y. King (1970) The biomechanics of spinal and head impact: problems of mathematical simulation. hoc. Swrp. Biorl.w. .Motld/hy ad Its .?.pplicarions. pp. 7017%. Dayton. Ohio. Liu. Y. King and Chandran. K. B. (1973) The exact soluuon to the rotational acceleration of rigid spherical shells filled with a brain-like elastic material. Proc. IOrh .&mii-erw-y .\/ecci~~g qf‘ rhr Sot. o/ Enqq Sci. Raleigh. S.C. Liu. Y. King and Chandran. K. B. (197-L) Comment on paper: Mathematical model of a head subjected to an angular acceleration. bv G. N. Bycroft. J. Biomechanics 7. 319-321. McElhaney. J. H.. bvlelvin. J. W.. Roberts. 1’. L. and Portnay. H. D. (1972) Dynamic characteristics oi the tissue of the head. .S.t~rrp. Persprcrires in Biomrti. Err~gmg (The Univ. of Strathclyde). Macmillan. London. Ommaja. A. K.. Fax. F. and Yarnell. P. (196s) LVhiplash injury and brain damage. J. .1m. Med. .-lssoc. 204. 235. Ommaya. X. K. and Hirsch. A. E. (1971) Tolerance for cerebral concussion from head impact and whiplash in primates. J. Biomrchanics 1. 13-Z. Shuck. L. F. and .Advani. S. H. (1972) Rheological response of human brain tissue. Presented at ASXlE Annual IVinter Meeting. New York. J. Busic Enguq ASME T,aus. 94, Series D. Paper No. 72 WA BHF-2. Lnterharnscheidt. F. and Higgins. L. S. (1969) Traumatic lesions of brain and spinal cord due to nondeforming angular acceleration of the head. LSuir. Terns Reports BIO/. .Lfrri. 27( 1). van Rosenberg. D. U. (1969) .Lfethotis /br rite .~~rmerical Sohrtio~ 01 Partial Diferetxial Eqwuions. Elsevier. New E’ork.
of viscoelastic
291
maternal viscoetastic parameter shear modulus. lb i in’ viscosity coefficient. Ibm set-in densitv.. Ibm’in’ , shear stress frequency of sine pulse acceleration
E( = p’ pad p p’ P :0
‘at the top indicates dimensional variables. variables have been nondimensionalized.
input
IVithout
: the
APPESDIT Xumrricai
solution
~iriscoelastfc
The differential
=z
s
equation
)
i
h
r-x.
Two new variables
material
ecpiutior~s
can be put in the form
(A-1)
i
are defined as 9
s=rzG” Sr and
The original
equation
then becomes
and the two new variables
are related
sv ss -_=--f-r.
2
Sr
r
?t
by (A-5)
Equations (A-3-A-5) comprise a set of equations in which the highest order derivative in either independent variable is of the first order. The initial conditions for the three variables are ru=v=s=O The boundary
forallratr=O.
conditions
w = r = 0
for all
at r = 0 are
t. and at r = I
SC - = 1 for a unit step function.
c’t
ll
radius of the spherical shell. in. amplitude ofexternal acceleration shear wave velocity. in:sec auisymmetric spherical coordinates time variable displacement in B direction external acceleration input
.-I C
i-. 0
t. : II:&. l(f)
This set of equations is solved by the centered-difference method described by van Rosenberg (1969). A grid of points is set up in the radial direction. and values of the dependent variables are computed at these points at set intervals of time. A set of points. with the indexing nomenclature. is shown in Fig. 13. For conv-enience in using the solution algorithm, the space index is given one set of values for sand another set for r and w as shown.
(Fig. I)
S
I
2
L-1
L‘.H
1
?
i
L-l
i
i+
1
Fig. I?. Grid of points.
L
L Lfl
292
Y.
k&c
LK. K. B.
CHANDRAS
The finite difference analog to equation (A-31 IS centered half-way bctwscn two consecutive time levels and is
This equation is second order correct in C. The finite difference analogs to the terms in equations (A-4 and A-S) are centered in time and space. The derivati\es of c are replaced by
and D. Li. VONROSESBERG When these finite difference analogs are substituted into the original partial di&rential equations. (A-4 and A-j). two sets of finite difference equations in the unknowns si,“_ / and c,,,_ , result. These equations are
[T]s
I_,. n-i + [I +(l
;,I”“;;;;):)]%~
1
At - 7~ = [s,,, - xi_ ,..I + [ - IAr bc,., + wi_ ,,.I [ Ar
and
1
i. cc -4 ! C,-I.“-1 - CL-l.” c ~i.n-l - Ci.” i Sr 1i++.“-* 1 At A’ (A-7)
+
1 _ At[(Ar - 2~) - x(At + ZE)] 2(i - Z)‘(Ar)’
[
At[(At + 2~) - r(Ar - ZE)]
The analogs for s are similar. and all are second order correct in space and time. The analog to the last term in equation (A-j) is
Z(i - l)‘(Ar)’
1CL” -I.n 1c,
(A-ila)
and
si- I.“_ + [
1
&
1.(A-S)
[(I - 4L.“-l A (’ + x’r:.J
s,1.” _, + si-
+[&(_ Ar
The parameter. x. can be varied to form any of three analogs; namely. the four-point analog for x = 0. the forward diagonal analog for x = I, and the backward diagonal analog when r = -I. For the boundary conditions of this problem. the backward diagonal (z = - I) is used Tar t 2 I, after which time the four-point analog (2 = 01 is used. The analogs to the terms Zc,r’ and ZW in equation (A-5) are handled in a similar manner. The values of w at the r._, level are eliminated from the final equation by equation (A-3a). All these analogs are also second order correct in space and in time. An analog, which is also second order correct in space and time. for the cross-derivative of s is given by 2
I +
,.”
-
Si.” -
(Ar)(Ar)
si-
I.“-
1
(A-9)
+[-$(I
[$(I+;)]b+, +
Se..-!
, , I, + “)lC,_, I-I
+~)],..+~(I
“i, = s
I.”
I s,_,n
-SJ]t++,.“. (A-5a)
These equations are written for the L-l values of i from i = Z to i = L. Thus. there are 2(L-I) equations in the set. All values at time f, are known either from the initial conditions or from the computations of the previous time step. Thus, there are ZL unknowns in the equations. However, two of these values are given by the boundary conditions, and the equations can be solved for the dependent variables s and c. An efficient algorithm for this solution has been given by van Rosenberg (1969).
