On the dynamic characteristics of the human skull

Pergamon

ON THE

Vol. 34, No. 12, pp. 1339-1348, 1996 Copyright C) 1996 Elsevier Science lad Printed in Great Britain. All rights reserved 002o-7225/96 $15.00+ 0.0O

Int. Z Engng Sci.

Pll: S0020-72Z~96)00047-X

DYNAMIC CHARACTERISTICS HUMAN SKULL

OF THE

A. CHARALAMBOPOULOS Institute of Chemical Engineering and High Temperature Chemical Processes, GR 265 00, Patras, Greece

G. DASSIOS Department of Chemical Engineering, University of Patras and Institute of Chemical Engineering and High Temperature Chemical Processes, GR 265 00, Patras, Greece

D.I. FOTIADIS Department of Computer Science, University of Ioannina, GR 451 10, loannina, Greece

V. KOSTOPOULOS Department of Mechanical Engineering, University of Patras and Institute of Chemical Engineering and High Temperature Chemical Processes, GR 265 00, Patras, Greece

C.V. MASSALAS Department of Mathematics, University of Ioannina, GR 451 10, Ioannina, Greece AIWlr~l--ln this work an attempt is made to study the dynamic characteristics of the human dry skull. The analysis is based on the three-dimensional theory of elasticity and the representation of the displacement field in terms of the Navier eigenvectors. The frequency equation was solved numerically and the results obtained are fairly good, in comparison to the experimental ones. Copyright ~) 1996 Elsevier Science Ltd

1. INTRODUCTION Cranial biomechanics has drawn the attention of many researchers as a consequence of the observation that fatalities resulting from accidents often involve injury to the head. Information about the physical processes and material response involved is of enormous importance to neurosurgeons and also in the design and construction of protective environments. For the analysis of human body processes there are two types of models in use: physical models based on experimental analysis and mathematical models based on theoretical analysis. The mathematical modelling of the cranial system has received considerable attention in recent years and this is evidenced by the number of publications [1-9]. The human skull is a complex structure made up of several bones, each with its own unique internal and external geometry. In analytical investigations, geometrical approximations and mainly the geometry of spherical and prolate spheroidal shells have been used [2-4]. A parametric study of head models by utilising an analytically based numerical technique was presented in [10]. The mechanical properties of cranial bone were studied in [11]. Obviously, for the collision and head injury analyses we need accurate dynamic response characteristics of the cranial system. In [12] an experimental investigation was presented to identify the dynamic characteristics of freely vibrating human skulls. In the present investigation, adopting the approximation of the human skull by a hollow sphere, an attempt is made to present an analysis for the vibrational characteristics of the human skull. The mathematical analysis is based on the three-dimensional theory of elasticity and the representation of the displacement field in terms of the Navier eigenvectors [13]. It is noted that an exact study of the dynamic characteristics of a closed spherical shell for arbitrary values of its thickness is desirable not only from the theoretical point of view, as well as for the purpose of the human skull dynamic characteristics, but also to estimate the range of applicability of the various shell theories [14-16]. In the framework of 1339

1340

A. C H A R A L A M B O P O U L O S et al.

the present analysis the frequency equation was solved numerically and results for the eigenfrequencies and mode shapes of the system considered are presented in tables and graphs. It is important to note that the results obtained are in agreement with the analogous ones by using shell theories [17], in the range of their validity, as well as with experimental results cited in Ref. [12]. 2. P R O B L E M

FORMULATION

For motions in homogeneous isotropic, elastic solids, the displacement vector field u satisfies the equation:

rff. V i~V2u(r,t) + (A + p,)VV.u(r,t) = p 02u(r't), Ot2

(1)

where A, ~t are Lam~'s constants, p is the mass density, V is the usual del operator, and t is time. In what follows we shall discuss the free vibrations problem of a hollow sphere in the framework of three-dimensional theory of elasticity (Fig. 1). For the problem under discussion we assume that:

u(r,t) = Re[uo(r)e i'°']

(2)

where to denotes the angular frequency measured in radians/sec and i = x / - 1. Replacing (2) into the equation of motion we obtain the well-known Navier frequency equation: /~V2uo(r) + (A + #)VV'uo(r) + ptO2uo(r) : O, r ~. V~. (3) Introducing the dimensionless variables ¥, = r rm

tog I Cp

into (3) we lead to

(4)

c'2W2uo(r ') + (1 - c'2)V'W.Uo(r ') + ~2uo(r' ) = O, where V' = rlV, c"

cp

The interaction of the system considered with its own environment enters the mathematical formulation through the boundary conditions on the surfaces So and $1 (Fig. 1).

Y

r

X

Fig. 1. Problem geometry.

Dynamic characteristics of the human skull

1341

We assume that the surfaces S, i = 0,1 are stress free, that is Tuo(r') = O, r' = r~ T u o ( r ' ) = O, r' = r~

(5)

T = 2#'~.V' + A'~V'. +/~'~ × V'

(6)

where

stands for the dimensionless surface traction operator in V, P is the unit outward normal vector on S,, i = 0,1 and (A',/~') = (A//~,I). We note that the problem described by the equation of motion, (4), and the boundary conditions (5) is a well-posed mathematical problem.

3. FREQUENCY EQUATION For the solution of the problem under discussion we adopted a method based on the representation of the displacement field Uo(r) in terms of the Navier eigenvectors. It is well known that Navier eigenvectors are a result of the Helmholtz decomposition and constitute a complete set of vector functions in the space of solutions of time-independent Navier equation. The Navier eigenvectors have the following form: o ' ( o ~ lr '

L'~'t(r ') = g~(f~r')P~. (~) + ~

~

B~'~(~)

~r

l~'t(r ') = Vn(n + 1)gt~(k'r')C~.(~) t , . 'r" ) r. ,. (. r. ) + ~ { g t l~,"(r') = n(n + 1 ) ~gntxs r, n = 0 , 1 , 2 .....

1=1,2,

I tr, (k;r,)-~ g'(k'rk, r, )~B~. " ~ ,' l ~ )

m= -n,-n+l

.... , n - l ,

n, r e V

(7)

where

~; = - , c"

d g'~(z) = ~ (g'~)

(8)

and g~(z) and g~(z) represent the spherical Bessel functions of the first kind, j.(z), and second kind, y.(z), respectively. The functions/'~ (~),/tr~ (P) and C~'~(~), defined on the unit sphere, are the vector spherical harmonics introduced by Hansen in spherical polar coordinates (r,O,~p), and are given as follows: (~) = ~ ( ~ ) (9)

~(~) - ~

1

{ '~ ~a + # s m - ~ } ~ ( ~ )

1

C~'~(P)=~I

sinO0~

~O~} I'~"(P)

(10)

(11)

where ~ and ~ are the unit vectors in 0 and ~p-directions respectively, Y~'~(~) = P~'~(cosO)e ~ are the spherical harmonics and P~'~(cosO) are the well known Legendre functions. We assume that the displacement vector field uo(r) has the representation 2 , {a,m , .Ln (r) + 13"~'tlt~'t(r) + ?/n"~'q~'t(r)). n--O m - - - - n / - - I

(12)

1342

A. CHARALAMBOPOULOS et al.

The spherical polar coordinates of the real part of uo(r) are given as Uo. = ~~o , ~==-

oc

~

,=, {ot~.(k'p,r')cos(m~p)P~(cosO)+ fa:-g~

~

....

/k

cos(m¢)P~ (cosO)} 3'~(n(n + 1 ) )g~(k'r') ~

cos(m~p)~/~ (COSO)- , t g t . ( k : r ' ) ~ t p )

~l

t .t ,, g.(k. O +,~(g.(k.r)÷ ,_-777_,))cos(m~p)-~P~(cosO) Ksr

\

/

/]

Off

P~(cosO)- fltng:(k'r')cos(mtp)O-fiz.Zooso>l

- a : g:(-~_~,' ) - ~

t 'r' )'~ ms in (m~o) . . . . t [ .t "k'r'" Jr g,(k,

....

] (12a)

The problem now is reduced to the determination of the coefficients ot,m'~,/3~''~ and ~,t. Since the expression (12) satisfies the equation of motion, (4), it remains to ask the boundary conditions (5) to be satisfied. We note that the presence of the surface traction operator T in the boundary conditions (5) requires us to know how this operator acts on the Navier eigenvectors. After tedious and extended manipulations we obtained the following relations:

TL'.'t(~) = - [ ~ , ' g~(i'~r') + 2/z',(1

n(n+l)\t-~2~i )g.(.r') + A',g~(f~r')]/~. (,)

g.t, -t'flr''-~r,2) j/ ~ ( e )

g"'' ( Or' )

+ 2p.'~[

(13)

TM'~'~"(') = I ~ ' ~ [ k ' ~ . ( k ' r ' ) - l g : ( k : r ' ) ] C ~ ' ~ ( ' )

Tl~'t(~)= 21x'n(n + 1)[

gl(k'r') ~

[

+a ' ~

-

(14)

/ g-t -"k'r"" ~ s ) ]P~'.(~)

2g.(ks"'r')

r'

, t

,,

- k~g.(k.r ) +

2n(n + 1 ) - lg~(k;r,)]B,~(~). (15) k,r, 2

Inserting (12) into the expression for the boundary conditions (5). and taking into account the relations (13)-(15) as well as the advantage of the orthogonality arguments for the vector spherical harmonic functions we conclude that for every specific pair of integers (n,m) (with Iml ~ n) the six coefficients involved in expression (12) satisfy a linear algebraic homogeneous system with six equations. This system in matrix form is given as -

I A.(r 0t 0 I B.(r 0t I

r

I

t

A.(ro) 0 B.(r0)

2 A.(r 0p 0

0 1 , C.(rO

2 B.(r 0r 2 P A.(r0) 0

0 0 1 , C.(ro)

2

t

B.(r0)

0

0 2 , C,(r0 0 0

1 w D.(rO 0

2 r" D.(rO

I t E.(rl) I t O.(r0)

2 r E.(rO 2 ; D.(ro) 0 E.(ro)

C2(ro) 0

0 1

t E.(ro)

2

a.

I0

m,2

0

e,~

t

't

iO 0

(16)

0 0 0

or in compact form,

O__x=0

(17)

Dynamic characteristics of the human skull

1343

where Al(r,) = - [ 4 g t (DJ") + 2~(1

Bt,(r ') = 2 ~

C~(r') = ~

n(n 1)~27,/~ + gl(Da.,) + A,f~gl (fb,,)]

[ 1 ~,t(f~r, ) _

[

- 2

f~,2 j

[ k'~,(k" r') - 1 gt (k~r,) ]

[ 'r') O~(r') = 2n(n + D IL g"`--r' Et,(r ') = ~

gt,(f~r')]

g',(k'r' ' ' )] -Ksr , ,2

/J

g,(ksr ) , t , , 2n(n + 1 ) - l g t ( k , r , ) ] " _'______~' r' k,g,(ksr ) + k,r, 2

In order for system (17) to have a nontrivial solution the following condition must be satisfied: det(D) = 0.

(18)

This condition provides the characteristic (frequency) equation, the roots of which are the eigenfrequency coefficients f~ of the system under discussion.

4. N U M E R I C A L

RESULTS--DISCUSSION

The frequency equation (18) has been solved numerically and for this purpose a matrix determinant computation routine was used for different frequency coefficients t2, along with a bisection method to refine steps close to its roots. The root finding algorithm is followed by an LU decomposition and back-substitution routine to determine the eigenvector x_ whose elements are used for the computation of the corresponding displacement components. Given that for every particular n we get independent systems, we deduce that finally the eigenvectors are obtained by (12a) without the external summation over n. The computations were made for material properties analogous to those of human skull [11]: E = 1.379 x

10 9 (N/m2),

v = 0.20 - 0.35, p = 2.0 × 103 (kg/m3),

rl = 0.040 m--0.082 m, ro = 0.036 m---0.08199 m. In Table 1 are cited the first three roots of (18) for different n (n is the order of the Bessel function). From these results we observe that we have degeneracy, that is the frequency aJ = 1.8924cflrl, corresponds to two different mode shapes which are shown in Fig. 1 for n = 1 and n = 3, respectively. It is obvious that the dynamic characteristics of the human skull are influenced by its geometry as well as from its own material properties. The variation of f =
f = 8,,,/rl where the first nine values of 8,. are cited in Table 2. The variations of mint2n,n = 1,2,3 with h/r~, h = r~ - ro, and Poisson's ratio v are shown in Figs 3 and 4 respectively. From the literature [15] it is known that for h/r~-<0.10 the improved theory for spherical

1344

A. CHARALAMBOPOULOS et al. Table 1. Eigenfrequency coefficients ~ for rn = 0.082 m, ro = 0.076 m and n = 0,1,...,8

1

n=0

n=l

0

0

n=2

2 3 4 5 6 7

n=3

n=4

n=5

n=6

n=8

0.7083 0.8522 0.9486 1.063 1.1971 1.2184 1.4210

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

n=7

1.5479 1.6695 1.8924

1.8924 2.5324 2.6253 3.1672 3.5229 3.7855 4.3983 4.4651 5.0077 5.4218 6.3843 7.3490 8.3143

Table 2. Values of 8mfor E = 1.379 × 109N/m 2, v = 0.25,0 = Z0 × lO~Kg/m3,(rl- ~)/rl = 0.0732 m

1

2

8m

102.56

3

121.46

4

134.80

153.92

5

6

173.32

7

173.88

204.88

8

9

224.08

240.48

shells [16] is q u i t e s a t i s f a c t o r y for thin shells. B e y o n d this t h e c o u p l i n g effects with thickness s t r e t c h m o d e b e c o m e i m p o r t a n t a n d a n y shell t h e o r y s h o u l d i n c l u d e this effect for b e t t e r results. F o r t h e s a k e o f c o m p a r i s o n w e p r e s e n t in T a b l e 3 t h e results o f o u r analysis which testify t h e p r e v i o u s conclusions. T h e a u t h o r s o f [12] i n t r o d u c e d an e x p e r i m e n t a l m e t h o d to i d e n t i f y t h e d y n a m i c c h a r a c t e r i s tics ( r e s o n a n t f r e q u e n c i e s a n d m o d e s h a p e s ) o f t h e h u m a n d r y skull. F o r t h e p u r p o s e o f c o m p a r i s o n with t h e results o f t h e p r e s e n t analysis w e s e l e c t e d two cases o f [12] which a r e r e f e r r e d to as skull I a n d skull II. Skull I a n d its c h a r a c t e r i s t i c s c o r r e s p o n d to t h o s e o f a 50th p e r c e n t i l e a d u l t m a l e , while t h o s e o f skull II c o r r e s p o n d to a 5th p e r c e n t i l e a d u l t f e m a l e .

n:3

n=l . ."" .. ...

........

.... .o. "...

Up"" u9

?"

:.. )

.,

i -.

..........

..'"

;

Fig. 2. u, and u~, for f2 = 1.8924, O = x/2 and n = 1,3.

"'"-. .......... ..-""

Dynamic characteristics of the human skull I

I

I

1345

I

I

O

1.5"

0

#3~

D

[3

O

1

0 D ×

8 x

0.5 0.00

n

N

o

#2

D

×

X X

#I

X

i

i

i

i

i

0.10

0.20

0.30

0.40

0.50

0.60

h/r, Fig. 3. Variation of t3 with

h/rl.

Table 3. Eigenfrequency coefficients f l , for n = 0,1,...,8, n

Present analysis

0 1 2 3 4 5 6 7 8

h/r,

= 0.10, v = 0.30

Shell theory Improved [17]

Classical theory

0 0 0.707 0.870 1.004 1.174 1.399 1.682 2.014

0 0 0.708 0.873 1.014 1.202 1.462 1.798 2.208

0 0 0.7072 0.8683 0.9992 1.1687 1.3953 1.6791 2.0134

T h e equivalence of g e o m e t r y used in o u r c o m p u t a t i o n s is as follows: Skull I:

ro = 0.0691 m, r~ = 0.07735 m

Skull II:

ro = 0.0602 m, r~ = 0.06845 m

T h e results o b t a i n e d as well as the analogous experimental ones of Ref. [12] are presented in T a b l e 4 and graphically in Fig. 5. Table 4. Comparison with experimental results Skull I Experiment [12] 1 2 3 4 5 6 7 8 9 10 11 12 13 14

1385 1786 1903 2449 2857 3386 3523 3845 4069 4245 4636

Skull II Present analysis 1331 1649 1911 2248 2772 2922 3371 3577 4063 4812 4962 5911 6655 7170

(n (n (n (n (n (n (n (n (n (n (n (n (n (n

= 2) = 3) = 4) = 2,5) = 6) -- 0) = 6) = 1,3) = 8) = 4) = 2) = 5) = 3) = 6)

Experiment [12] 1641 2344 2969 3477 4453 5000

Present analysis 1525 1906 2265 2561 2752 3323 3386 4063 4157 5015 5459 5628 6812

(n (n (n (n (n (n (n (n (n (n (n (n (n

= 2) = 3) = 4) = 2) = 5) = 0) = 6) = 1,3) = 7) = 8) = 4) = 2) = 5)

1346

A. CHARALAMBOPOULOS et al. 1.2

I

I

I

I

I

I

I

I

0.20

0.25

0,30

0.35

1.1 1

0.9 0.8 0.7 0.6 0.5 0.15

0.40

Y

Fig. 4. Variation of f~ with v. For the case of skull I the results obtained are in good agreement with the experimental ones. In the first eigenfrequency predicted by our analysis, )~ = 1.331 Hz, a shift of about 4% with respect to the experimental one (~exp -- 1.385 Hz) appears, while in f~ = 4962 Hz a shift of about 6.5% in the opposite direction from the experimental one, J~exp--4636 Hz. In the case of skull II the measured frequencies are predicted by the present analysis. However, in the interval [1525, 5015] more eigenfrequencies exist than are presented in Ref. [12]. This fact needs more study from experimental and theoretical point of view. In Fig. 6 mode shapes ( u , , u , , u o ) are presented corresponding to the first five eigenfrequencies given in Table 1 (1: n = 2, ~ = 0.7083; 2: n - - 3 , Q = 0.8522; 3: n = 4, fl = 0.9486; 4: n = 5, ~ = 1.063; 5: n = 2, ~ = 1.1971).

5. C O N C L U D I N G R E M A R K S In the analysis presented we dealt with the dynamic characteristics of the human dry skull. The results obtained are favourable compared with the experimental ones and this fact leads us I

8.0 10 3

I

I

I

I

I

I

I

I

I

1

I

7.0 10 3-

I

D

6.0 10~ []

5.0 10 3" W (Hz)

\

i

4.0 10 3-

c]

3.0 10 3" 2.0 10 3"

I 1.0 10 3

1

0

I

2

I

I

4

I

I

6

I

I

8

I

I

10

I

I

12

Mode #

Fig. 5. Skull l: comparison with experimental results.

i

1

14

Dynamic characteristics of the human skull NO.

1347

tl

U¢

Ur

UlI,

(O = ~r/ 2)

(O = Jr/ 2)

(¢ = x)

(9 TM ~r)

I

\ /

.-

......

i

\,

/

\

Fig. 6. Mode shapes (u,, u¢, uo).

to the conclusion that the geometry considered as well as the analysis adopted (threedimensional theory of elasticity) could be applied to the study of the cranial system. The results of the dynamic analysis of the human skull-brain system, as well as those of head-neck system are in preparation and will be presented in the future. Acknowledgement--The present work forms part of the project 'New Systems for Early Medical Diagnosis and Bioteclmological Applications' which is supported by the Greek General Secretariat for Research and Technology through the EU funded R&D Program EPET IL REFERENCES [1] R. HICKLING and M.L. WENNER, J. Biomechanics 6, 115 (1973). [2] J.C. MISRA, Ingenieur-Archiv 4"1, 11 (1978).

1348 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]

A. CHARALAMBOPOULOS et al.

J.C. MISRA, C. HARTUNG and O. MAHRENHOLTZ, Ingenieur-Archiv 47, 329 (1978). J.C. MISRA and S. CHAKRAVARTY, Acta Mechanica 44, 159 (1982). J.C. MISRA and S. CHAKRAVARTY, J. Biomechanics 15, 635 (1982). CHUNG-SHENG CHU, MISH-SHYAN LIN, HAW-MING HUANG and MAW-CHANG LEE, J. Biomechanics 27, 187 (1994). H.J. WOLTRING, K. LONG, P.J. OSTERBAUER and A.W. FUHR, J. Biomechanics 27, 1415 (1994). S.P. NANDA, Bull. Cal. Math. Soc. 83, 337 (1991). A.E. ENGIN, J. Biomechanics 2, 325 (1969). T.B. KHALIL and R.P. HUBBARD, J. Biomechanics 10, 119 (1977). J.H. McELHANEY, J.L. FOGLE, J.W. MELVIN, R.R. HAYNES, V.L. ROBERTS and N.M. ALEM, J. Biomechanics 3, 495 (1970). T.B. KHAL1L, D.C. VIANO and D.L. SMITH, J. Sound and Vibration 63, 351 (1979). W.W. HANSEN, Phys. Rev., 47, 139 (1935). A.H. SHAH, C.V. RAMKRISHNAN and S.K. DATrA, Trans. ASME, J. Appl. Mech. 36, 431 (1969). A.H. SHAH, C.V. RAMKRISHNAN and S.K. DATI'A, Trans. ASME, J. Appl. Mech. 36, 440 (1969). F. NIORDSON, Int. J. Solids Structures 20, 667 (1984). E. MAGRAB, Vibrations o f Elastic Structural Members. Sirthoff and Noordhoff, Netherlands (1979). (Received I March 1996; accepted 12 March 1996)