Study of a spherical head model. Analytical and numerical solutions in fluid–structure interaction approach

International Journal of Engineering Science 51 (2012) 1–13

Contents lists available at SciVerse ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

Review

Study of a spherical head model. Analytical and numerical solutions in fluid–structure interaction approach Adil El Baroudi a,⇑, Fulgence Razafimahery b, Lalaonirina Rakotomanana-Ravelonarivo b a b

M3M Laboratory, Belfort-Montbliard University of Technology, 90010 Belfort Cedex, France IRMAR, Mechanics Group, University of Rennes 1, 35042 Rennes Cedex, France

a r t i c l e

i n f o

Article history: Received 19 June 2011 Received in revised form 28 October 2011 Accepted 14 November 2011 Available online 9 December 2011 Keywords: Fluid–structure interaction Modal analysis Confinement study Finite element method

a b s t r a c t This paper deals with the analytical and the numerical computation of elastoacoustic vibration modes. We determine in the present study the eigenfrequencies and the modal shapes of the system of brain, cerebro-spinal fluid (CSF) and skull. Two models are presented in this work: an elastic-acoustic model assuming a rigid skull and an elastic-acoustic-elastic model assuming a deformable skull. The analytical results are compared with the numerical solution obtained using the software Comsol Multiphysics. It is shown that eigenfrequencies and more significantly the modal shapes are strongly influenced by the interaction between solid phases (brain and skull) and the cerebro-spinal fluid. Finally the influence of the CSF compressibility and thickness on the natural frequencies was investigated.  2011 Elsevier Ltd. All rights reserved.

Contents 1. 2. 3.

4.

5.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Basic models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3.1. Fluid domain: dynamic equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3.2. Solid domain: dynamic equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 3.3. Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Elastic-acoustic-elastic model: deformable skull . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 4.1. Analytical solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 4.2. Numerical solution: variational formulation in (u(1), u(2), p) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 4.3. Analytical and numerical results and discussion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 4.3.1. Influence of the CSF compressibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 4.3.2. Influence of the CSF thickness: Effect of confinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4.3.3. Modal shapes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Elastic-acoustic model: rigid skull . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 5.1. Analytical solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 5.2. Numerical solution: variational formulation in (u(1), p) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 5.3. Analytical and numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 5.3.1. Influence of the CSF compressibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5.3.2. Influence of the CSF thickness: Effect of confinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5.3.3. Modal shapes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

⇑ Corresponding author. E-mail address: [email protected] (A. El Baroudi). 0020-7225/$ - see front matter  2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijengsci.2011.11.007

2

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

6.

Conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Appendix A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1. Introduction On the one hand, brain injuries constitute one of the major cause of death in road accidents. To understand how the brain gets injured during an accident, the mechanical response of the contents of the head during an impact has to be known. Several models have been developed to study this phenomenon. Some authors used numerical models using the finite element method (Marjoux, Baumgartner, Deck, & Willinger, 2008; Rotha, Vappoua, Raula, & Willinger, 2009), while others developed experimental devices that simulate a shock (Hault-Dubrulle, 2007; Verschuerena et al., 2007). Most numerical models used did not always show the coupled character of the problem studied (Gong, Lee, & Lu, 2008; Ho & Kleiven, 2009; Zonga, Leeb, & Lu, 2006). The consequences of the skulls vibrations are still poorly understood. It is probable that the low-frequency skull vibrations (below 200 Hz) mainly cause deep cerebral lesions, while higher frequency vibrations have more consequences on the superficial cerebral structures (Will & Berg, 1996). Due to the presence of a thin layer of the cerebro-spinal fluid (CSF) between the brain and the skull, relative motion can occur and may explain many types of brain injury. First, intracerebral hematomas may appear due to bridging veins rupture. Then, contusions or bruising is likely to occur when the brain hits against the inner surface of the skull. Last, coup and contre-coup injuries can happen when the head is suddenly accelerated. The coup injury is caused by the brain hitting the interior of the skull; the contre-coup injury occurs directly opposite the blow due to a process called cavitation. One of the difficulties to apprehend the trauma brain injury (TBI) is the lack of informations about the stress distribution within the brain, the CSF and the skull because of the presence of a fluid phase. Solving fluid–structure mechanical interaction is thus expected to further improving our understanding of the impact of the head during an accident. On the other hand, the knowledge of the flow dynamics and oscillations of the CSF in the intra-cranial space plays an important role for various human pathology in the everyday activities. More recent studies on the brain tend more and more to point out the role of wave within the brain tissue e.g. (Will & Berg, 2007). They experimentally showed that synchronization responses of brainwaves to periodic auditory stimuli had three components corresponding to different spectra: low frequencies [1–5 Hz], intermediate frequencies [3–8 Hz], and higher frequencies [11–44 Hz]. In addition, mechanical vibrations also play important role during measurement as Diffusion Tensor Imaging by means of Magnetic Resonance Imaging. Recent studies showed that the mechanical vibrations cannot be ignored and they should be considered when choosing the sequence parameters for Diffusion Tensor Imaging (Hiltunen, Hari, Jousmaki, & Mueller, 2006). All these observations suggest us to better determine the natural frequencies and modal shapes of the brain, the CSF and the skull. Analytical methods were conducted to determine the frequency spectrum of the head–neck system (Charalambopoulos, Dassios, Fotiadis, & Massalas, 1997). Finite element method has been intensively used to model the head impact (Kleiven, 2002). Most of them neglected the mechanical role of the CSF. A basic question remain: what would be the role of the cerebro-spinal fluid for free and forced vibrations of the brain-CSF-skull system at low frequencies 6200 Hz? Numerical aspects of fluid–structure problems have been studied in the past (Morand & Ohayon, 1979). Most classical methods use the modal shapes of the solid phase for deriving the dynamical behavior of the whole (solid and fluid domains). In a general manner, classical methods may be sufficient if we are only interested in searching for the natural frequencies. However, these methods do not permit to obtain the modal shapes of the whole system which is essential for applying the modal projection when facing dynamical situation (forced vibrations or shocks). The goal of the present work was to search for analytical and numerical modal shapes of the interacting brain-CSF-skull assuming elastic behavior for the brain and the skull and an acoustic wave propagation within the CSF. The model is designed to directly solve the coupled solid–fluid problem, conversely to classical method that only uses the modal shapes of the solid phase. The influence of the CSF compressibility and thickness on the natural frequencies and modal shapes was investigated. Last but not least, the use of simple models was deemed necessary since validation of the fluid-structure problems still remains a challenge. 2. Basic models We adopt the assumption of linear elasticity in three dimensions and assume an inertial coupling, meaning that the longitudinal wave velocity is very large compared with the characteristic velocity of the fluid. The Helmholtz decomposition of the displacement field of a structure (Morse & Feshbach, 1946) reduces the problem into a system of scalar partial differential equations. The two basic models developed in the present study are represented in Fig. 1. The model is elastic-acoustic when the skull assumed to be rigid and elastic-acoustic-elastic when the skull is elastic. During a previous study, an experimental device was developed for a three-dimensional analysis of a head during an impact (Hault-Dubrulle, 2007). Geometry and material properties of the present study were derived from this previous study (cf. Table 1).

3

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

Fig. 1. Elastic-acoustic-elastic and elastic acoustic models: brain (X1), cerebrospinal fluid (Xf) and skull (X1). Top view of a human brain-CSF-skull system.

Table 1 Geometry and material properties.

m

q½kg=m3 

R½m

c0 ½m=s

6,77.10

0.48 0.33

R1 = 0.0858 R2 = 0.0969 R0 = 0.1000

1450

7.1010

1150 1000 2750

E½Pa Brain Fluid Skull

(X1) (Xf) (X2)

5

These values take into account the fluid rotational movement and then differ from those of Charalambopoulos et al. (1997) which considered the brain tissue as an inviscid irrotational fluid. 3. Governing equations 3.1. Fluid domain: dynamic equation The fluid (CSF) is assumed non-viscous and isotropic which satisfy the acoustic wave equation. The equation of motion of the fluid can be written by Morse and Ingard (1968) as

r2 p 

1 @2p ¼0 c20 @t2

ð1Þ

where p is the acoustic pressure and c0 is the compressional wave velocity in the fluid and

r2 ¼

    1 @ @ 1 @ @ 1 @2 r2 þ 2 sin h þ 2 2 r @r @r r sin h @h @h r2 sin h @ u2

is the laplacian in spherical coordinates (r, h, u) Fig. 2.

ð2Þ

4

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

Fig. 2. Spherical polar coordinates system.

In order to compute the natural vibration mode of the coupled system, we consider harmonic solutions for the displacements and for the pressure p; that is, we seek solutions of Eq. (1) of the forme

pðr; h; u; tÞ ¼ Z 1 ðrÞZ 2 ðhÞZ 3 ðuÞejxt ;

2

j ¼ 1

ð3Þ

with x is the coupled angular frequency. We subsitute these expression in Eq. (1) we obtain the solution (see for example, Reference Abramowitz & Stegun, 1970)

p¼

rffiffiffiffiffiffiffiffih

p

2dr

i GJnþ1 ðdr Þ þ HY nþ1 ðdr Þ Pm n ðcos hÞ sinðmuÞ; 2

2

d¼

x

ð4Þ

c0

where n and m are the separation constants. The functions J nþ1 ðdrÞ and Y nþ1 ðdrÞ are spherical Bessel functions of the first and 2 2 second kind, respectively. Pm n ðxÞ is the associated Legendre function. G and H are arbitrary constants that may be determined from boundary conditions. 3.2. Solid domain: dynamic equation The spheres are assumed to be made of an elastic solids with Young’s modulus Ei, Poisson ratio mi, mass density qi and i = 1, 2 indicates the brain (i = 1) and the skull (i = 2), respectively. For motion in homogeneous isotropic elastic solid, the displacement vector field u(i), satisfies Navier’s equation

@ 2 uðiÞ ¼ r  rðuðiÞ Þ @t 2

qi

ð5Þ

The stress tensor in the solids r(u(i)) is given by Hooke’s law as

rðuðiÞ Þ ¼ ki ðr  uðiÞ ÞI þ 2li eðuðiÞ Þ

ð6Þ

(i)

Here, I is a unit tensor, e(u ) is the deformation tensor, li,ki are Lamé’s elastic constants, qi is the mass density. By using the Helmholtz decomposition (Morse & Feshbach, 1946)

uðiÞ ¼ r/ðiÞ þ r  ðwðiÞ ez Þ þ r  r  ðvðiÞ ez Þ (i)

where / , w

(i)

(i)

and v

ð7Þ

are scalar potential functions. Here the total motion may be decomposed into three parts:

 r/(i) represents the dilatational (longitudinal) motion propagating with speed cLi and wave number ai ¼ cxL . i  r  (w(i)ez) and r  r  (v(i)ez) are the rotational (transverse) parts, moving with speed cTi and wave number bi ¼ cxT . i

(i)

Therefore, the scalar components of the displacement vector u ðiÞ

@/ @ ðr v Þ  r r2 vðiÞ þ @r 2 @r  ðiÞ   @ rv 1 @ 1 @wðiÞ þ /ðiÞ þ ¼ r @h sin h @ u @r  ðiÞ   @ r v 1 @ @wðiÞ ¼ /ðiÞ þ  @r r sin h @ u @h

urðiÞ ¼ ðiÞ

uh

uðiÞ u

2

in Spherical coordinates can be expressed by

ðiÞ

ð8Þ

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

5

and the some components of stress tensor are related to the displacement potentials as

(

!)   2 @ 2 /ðiÞ ki 2 ðiÞ @ @ r vðiÞ 2 ðiÞ þ r / þ 2  r r v @r @r 2 li @r 2 h

i @ 1 @ ðiÞ rrhðiÞ ¼ li UðiÞ þ  ðiÞ þ W @h r sin h @ u

i 1 @ 1 @ h ðiÞ ðiÞ rr u ¼ li U þ  ðiÞ  WðiÞ sin h @ u r @h

rrrðiÞ ¼ li 2

ð9Þ

where

" # 2 @/ðiÞ ðiÞ r  / r2 @r   @ 1 @ðr vðiÞ Þ  r2 vðiÞ ¼2 @r r @r

UðiÞ ¼  ðiÞ

WðiÞ ¼ r

ð10Þ

@wðiÞ  wðiÞ @r

These equations should be completed by boundary conditions and by equations of motions of the structure and continuity conditions at the fluid–structure interface. By expanding Eq. (7) and substituting it into Navier’s Equation (5), three partial differential equations are obtained in terms of the three potential functions

r2 /ðiÞ 

1 @ 2 /ðiÞ ¼0 cL2 @t 2 i

r2 wðiÞ 

1 @ 2 wðiÞ ¼0 cT 2 @t2 i

r2 vðiÞ 

1 @ 2 vðiÞ ¼0 cT 2 @t 2 i

with cLi ¼

qffiffiffiffiffiffiffiffiffiffi ki þ2li

qi

and cT i ¼

ð11Þ

qffiffiffiffi li qi

are the compression (longitudinal) and shear (transverse) wave velocities in the spheres,

respectively. Lamé’s coefficients may also be written as

ki ¼

Ei m i ; ð1 þ mi Þð1  2mi Þ

li ¼

Ei 2ð1 þ mi Þ

ð12Þ

Similarly, the solutions for Eq. (11) which belong to the longitudinal and shear waves inside the brain (i = 1) and the skull (i = 2), respectively, are represented in spherical coordinates by:

rffiffiffiffiffiffiffiffiffih

i p Ai J nþ1 ðai rÞ þ Bi Y nþ1 ðai rÞ Pm n ðcos hÞ sinðmuÞ 2 2 2ai r rffiffiffiffiffiffiffiffiffih i p wðiÞ ¼ C i J nþ1 ðbi rÞ þ Di Y nþ1 ðbi rÞ Pm n ðcos hÞ cosðmuÞ 2 2

/ðiÞ ¼

vðiÞ ¼

2bi r rffiffiffiffiffiffiffiffiffih p 2bi r

ð13Þ

i K i J nþ1 ðbi rÞ þ F i Y nþ1 ðbi rÞ Pm n ðcos hÞ sinðmuÞ 2

2

with Ai, Bi, Ci, Di, Ki, Fi are arbitrary constants that may be determined from boundary conditions. 3.3. Boundary conditions In expressions (4) and (13), the terms Ai, Bi, . . . , H correspond to the unknown coefficients to be determined from the following boundary conditions: 1. On the interface C1 between the fluid (CSF) and the brain (i = 1), the displacements (velocities) and normal stresses must be continuous, leading to

(

8 @p 2 ð1Þ > > > @r ¼ q0 x ur > > ð1Þ < r ¼ p rp  n ¼ q0 x2 uð1Þ  n rr ) ð1Þ ð1Þ > rðu Þn ¼ pn > rrh ¼ 0 > > > : rð1Þ ¼ 0 ru

ð14Þ

6

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

2. On the interface C2 between the fluid (CSF) and the skull (i = 2), the displacements (velocities) and normal stresses must be continuous, leading to

8 @p 2 ð2Þ > > > @r ¼ q0 x ur > ( > ð2Þ < rp  n ¼ q0 x2 uð2Þ  n rrr ¼ p ) ð2Þ > rðuð2Þ Þn ¼ pn > rrh ¼ 0 > > > : rð2Þ ¼ 0

ð15Þ

ru

3. At the outside boundary C0 (the outside interface of the skull), the normal stresses must be zero, leading to:

8 ð2Þ > < rrr ¼ 0 rðuð2Þ Þn ¼ 0 ) rð2Þ rh ¼ 0 > : ð2Þ rr u ¼ 0

ð16Þ

4. Elastic-acoustic-elastic model: deformable skull A more realistic problem would let the boundary C0 of the skull free. The free vibrations of the brain-CSF-skull system are governed by the Eq. (5) with boundary conditions Eqs. (14)–(16). 4.1. Analytical solutions Introducing the solutions (cf. Eqs. (13) and (4)) into the formulae for the displacements and the stresses and taking into account Eqs. (14)–(16), we obtain for each mode number n a linear system of equations for the unknown coefficients which can be written as

R1 J 0nþ1 ðb1 R1 Þ  J nþ1 ðb1 R1 Þ ¼ 0 2 2 h ih i R2 J 0nþ1 ðb2 R2 Þ  J nþ1 ðb2 R2 Þ R0 Y 0nþ1 ðb2 R0 Þ  Y nþ1 ðb2 R0 Þ 2 2 2 2 h ih i  R0 J 0nþ1 ðb2 R0 Þ  J nþ1 ðb2 R0 Þ R2 Y 0nþ1 ðb2 R2 Þ  Y nþ1 ðb2 R2 Þ ¼ 0

ð17Þ

Mv ¼ 0

ð19Þ

2

2

2

2

ð18Þ

where v is the vector of the unknown coefficients and M is a 8  8 matrix of the linear system. The matrix M is written as

0

a11

a12

a13

a14

0

0

0

B a21 B B B a31 B B 0 B M¼B B 0 B B 0 B B @ 0

a22

a23

a24

0

0

0

0

1

a32

0

0

0

0

0

0

a43

a44

a45

a46

a47

0

a53

a54

a55

a56

a57

0 0

0 0

0 0

a65 a75

a66 a76

a67 a77

0 C C C 0 C C a48 C C C a58 C C a68 C C C a78 A

0

0

0

a85

a86

a87

a88

0

ð20Þ

where the elements of V are (A1, K1, G, H, A2, B2, K2, F2). The coefficients of the previous quantities are given in Appendix A. For non-trivial solution, the determinant of the matrix M must be equal to zero

det MðxÞ ¼ 0

ð21Þ

Thus, the eigenfrequency equation can be obtained. The Eqs. (17)–(21) indicates a relationship of the celerity of sound waves in fluid and the elastic constants. The roots of Eq. (21) give the frequencies of the spheroidal oscillations and those of Eqs. (17) and (18) the torsional oscillations. This equation can be solved numerically to obtain the eigenfrequencies f ¼ 2xp. 4.2. Numerical solution: variational formulation in (u(1), u(2), p) Finite element method is used to extract the eigenfrequencies and modal shapes. To compute the natural vibration modes of a fluid alone, the fluid is typically described either by pressure or by displacement potential variables. When the fluid is coupled with a solid, standard methods to solve Eqs. (1)–(5) consist of eliminating either the pressure or the displacement potential. However, in both cases non-symmetric eigenvalue problems are obtained (see, for instance, Zienkiewich & Taylor, 1989). To avoid this drawback, introduce in Morand and Ohayon (1979) an alternative procedure which consists of using

7

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

pressure and displacement potential simultaneously. In this section we summarize their approach; further details and discussions can be found in their book (Morand & Ohayon, 1995). The variational formulation of Eqs. (1)–(5) holds

Z

Z

rðuð1Þ Þ : eðv 1 Þdx  x2

X1

Z

Z

rðuð2Þ Þ : eðv 2 Þdx  x2

X2

Z

rp  rqdx  x2

X2

Z

Xf

X1

q1 uð1Þ  v 1 dx  q2 uð2Þ  v 2 dx 

pq dx þ q0 c20

Xf

Z C1

Z Z

pv 1  ndC ¼ 0 C1

pv 2  ndC ¼ 0 C2

uð1Þ  nqdC þ q0

Z

ð22Þ !

uð2Þ  nqdC

¼0

C2

"(v, q) 2 V  Q in which

8 n < V ¼ ðv 1 ; v 2 Þ 2 H1 ðX1 Þ  H1 ðX2 Þ; :

Q ¼ H1 ðXf Þ;

o

v 2 ¼ 0ðC0 Þ

ð23Þ

ðv ; qÞ are a test functions ðiÞ

By using again Lagrange elements, where uh 2 P2  P2 et ph 2 P1, discretization of the variational formulation induces a non symmetrical system

2

K1

6 4O O

O K2 O

3 2 M1 U1 76 7 26 B2 54 U2 5 ¼ x 4 O Ma1 Kp P B1

32

O M2 Ma2

3 U1 76 7 O 54 U2 5 Mp P O

32

ð24Þ

where U1, U2 and P are the vectors of nodal values for u(1), u(2) and p, respectively, and the submatrices of Eq. (24) are defined by

Z Xi

Z

Ci

Z

rðuðiÞ Þ : eðv i Þdx ¼ VTi Ki Ui ; pv i  ndC ¼ VTi Bi P;

Z

rp  rqdx ¼ Q T Mp P;

Xf

Ci

Z Xi

qi uðiÞ  v i dx ¼ VTi Mi Ui

q0 ui  nqdC ¼ Q T Mai Ui Z Xf

ð25Þ

pq dx ¼ Q T Kp P c20

for i = 1, 2; where V1, V2 and Q are the vectors of nodal values for v1, v2 and q, respectively. Ma is the added mass matrix (symmetric and positive definite) (Morand & Ohayon, 1995). 4.3. Analytical and numerical results and discussion The eigenfrequency Eqs. (17)–(21) and (24) were solved. Maple and Comsol Multiphysics softwares are used to extract the first 12 eigenfrequencies fnl ¼ fi ¼ x2pnl [Hz] for different values of n. l is eigenfrequency’s order. According to the results of Table 2, we notice that the numerical and analytical results are in good agreement. This shows that the reliability of the implemented algorithm in Comsol Multiphysics software for numerical computation. 4.3.1. Influence of the CSF compressibility The physical properties of CSF may differ from one person to another, we investigated the influence of the CSF compressibility by varying the sound velocity c0 into the eigenfrequencies Eqs. (21) and (24). Analytical and numerical solutions are reported in Table 3. Eigenfrequencies do not significantly vary even within a quite large band c0 2 [300, 1450] [m  s1],

Table 2 Eigenfrequencies for various mode shapes of deformable skull. fi

n

Order (l)

Exact

FEM

Mode shape

f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 f12

2 2 3 4 1 3 5 4 2 6 1 5

1 1 1 1 1 1 1 1 2 1 1 1

39.01 65.45 65.95 92.27 100.52 101.13 118.22 133.32 139.62 143.76 150.82 163.97

39.01 65.45 65.95 92.27 100.52 101.13 118.22 133.32 139.62 143.76 150.82 163.97

Spheroidal Torsional Spheroidal Spheroidal Spheroidal Torsional Spheroidal Torsional Spheroidal Spheroidal Torsional Torsional

8

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

Table 3 Natural frequencies vs. CSF compressibility (deformable skull). fi

f1 f2 f3 f4 f5 f6 f7 f8 f9 f10

c0 = 1000

c0 = 700

c0 = 500

Exact

FEM

Exact

FEM

Exact

FEM

39.01 65.45 65.95 92.27 100.52 101.13 118.22 133.32 139.62 143.76

39.01 65.45 65.95 92.27 100.52 101.13 118.22 133.32 139.62 143.76

39.01 65.45 65.95 92.27 100.52 101.13 118.22 133.32 139.62 143.76

39.01 65.45 65.95 92.27 100.52 101.13 118.22 133.32 139.62 143.76

39.01 65.45 65.95 92.27 100.52 101.13 118.22 133.32 139.62 143.76

39.01 65.45 65.95 92.27 100.52 101.13 118.22 133.32 139.62 143.76

Table 4 Natural frequencies vs. CSF thickness (deformable skull). fi

f1 f2 f3 f4 f5 f6 f7 f8 f9 f10

e = 0.0088

e = 0.0066

e = 0.0043

Exact

FEM

Exact

FEM

Exact

FEM

34.77 59.54 63.74 84.39 97.91 98.52 109.42 129.42 135.83 137.62

34.77 59.54 63.74 84.39 97.91 98.52 109.42 129.42 135.83 137.62

30.18 52.35 62.19 75.10 95.52 96.11 98.42 122.19 130.33 132.36

30.18 52.35 62.19 75.10 95.52 96.11 98.42 122.19 130.33 132.36

24.44 43.01 60.64 62.56 83.06 93.15 93.72 104.31 128.89 123.10

24.44 43.01 60.64 62.56 83.06 93.15 93.72 104.31 128.89 123.10

Fig. 3. Some brain modal shapes (free skull).

showing that in the fluid–structure interaction problem we can neglect the term of CSF compressibility x2 =c20 and the pressure Eq. (1) becomes r2p = 0. 4.3.2. Influence of the CSF thickness: Effect of confinement Among other factors, cerebrospinal fluid (CSF) thickness has been studied. The thickness of the annular space between the two solids are variable under normal circumstances. We investigated the influence of the CSF thickness on eigenfrequencies. To this purpose, we take for the CSF radius R2 ¼ 0:0969 [m] and we vary e = R2  R1. Analytical and numerical results are reported in Table 4 showing that eigenfrequencies decrease with the CSF thickness.

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

9

Fig. 4. Some skull modal shapes.

4.3.3. Modal shapes Six modal shapes are reported in Figs. 3 and 4. The Von Mises stress field is represented since it plays a key role for understanding either the deep brain injury at low frequency loading. Modal shapes illustrate that the highest strain occur in the central region when accounting for the fluid-interaction. This is not the case when the skull is assumed rigid. It seems to suggest that traumatic brain injury at low frequency is localized at the central region and results into a diffuse axonal injury. 5. Elastic-acoustic model: rigid skull In this section, the skull is assumed to be rigid. The problem is governed by (5) with boundary conditions (14) and on the interface C2 between the fluid (CSF) and the skull we have

rp  n ¼ 0 )

@p ¼0 @r

ð26Þ

5.1. Analytical solution In this case, the Eqs. (17)–(19) are reduced to

R1 J 0nþ1 ðb1 R1 Þ  J nþ1 ðb1 R1 Þ ¼ 0

ð27Þ

Mv ¼ 0

ð28Þ

2

2

The matrix M is reduced to four rows and columns and is written as

0

a11 Ba B 21 M¼B @ a31

a12 a22

a13 a23

a32

0

a43

a44

a45

1 a14 a24 C C C 0 A

ð29Þ

a46

The elements of v are reduced to [A1 K1 G H]. 5.2. Numerical solution: variational formulation in (u(1), p) The variational formulation (22) is reduced to

Z

rðuð1Þ Þ : eðv 1 Þdx  x2

X1

Z

Xf

rp  rqdx  x2

Z Xf

Z X1

q1 uð1Þ  v 1 dx 

pq dx þ q0 c20

Z C1

Z C1

pv 1  ndC ¼ 0 !

uð1Þ  nqdC

¼0

ð30Þ

10

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

Table 5 Eigenfrequencies: analytical vs. numerical (rigid skull). fi

n

Order (l)

Exact

FEM

Mode shape

f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 f12

2 2 3 4 1 3 5 4 2 6 1 5

1 1 1 1 1 1 1 1 2 1 1 1

39.01 65.45 65.93 92.27 100.46 101.13 118.22 133.32 139.62 143.76 150.82 163.97

39.01 65.45 65.95 92.27 100.46 101.13 118.22 133.32 139.62 143.76 150.82 163.97

Spheroidal Torsional Spheroidal Spheroidal Spheroidal Torsional Spheroidal Torsional Spheroidal Spheroidal Torsional Torsional

Table 6 Natural frequencies vs. CSF compressibility (rigid skull). fi

c0 = 1000

f1 f2 f3 f4 f5 f6 f7 f8 f9 f10

c0 = 700

c0 = 500

Exact

FEM

Exact

FEM

Exact

FEM

39.01 65.45 65.95 92.27 100.46 101.13 118.22 133.32 139.62 143.76

39.01 65.45 65.93 92.27 100.46 101.13 118.22 133.32 139.62 143.76

39.01 65.45 65.95 92.27 100.46 101.13 118.22 133.32 139.62 143.76

39.01 65.45 65.95 92.27 100.46 101.13 118.22 133.32 139.62 143.76

39.01 65.45 65.95 92.27 100.46 101.13 118.22 133.32 139.62 143.76

39.01 65.45 65.95 92.27 100.46 101.13 118.22 133.32 139.62 143.76

ð1Þ

"(v1, q) 2 V  Q in which V = {v1 2 H1(X1)} and Q = H1(Xf). By using again Lagrange elements, where uh 2 P2  P2 et ph 2 P1, discretization of the variational formulation induces a non symmetrical system



K1

B1

O

Kp



U1



P

 x2



M1

O

Ma1

Mp



U1 P

 ¼

  0

ð31Þ

0

where U1 and P are the vectors of nodal values for u1 and p, respectively, and the submatrices of (31) are defined in (25) for i = 1. 5.3. Analytical and numerical results An appropriate descending analysis technique was used in the Comsol software for avoiding numerical spurious pressure modes. Analytical solutions was obtained by Maple software. The first 12 eigenfrequencies are reported in Table 5, fi ¼ 2xpi ½Hz. When the skull was assumed to be elastic, we found the same eigenfrequencies and modal shapes as for the rigid skull. Table 7 Natural frequencies vs. CSF thickness (rigid skull). fi

f1 f2 f3 f4 f5 f6 f7 f8 f9 f10

e = 0.0088

e = 0.0066

e = 0.0043

Exact

FEM

Exact

FEM

Exact

FEM

34.78 59.54 65.95 63.74 97.84 98.49 109.11 129.42 135.83 137.62

34.78 59.54 63.74 84.29 97.84 98.49 109.11 129.42 135.83 137.62

30.18 52.34 62.19 75.04 95.45 96.09 98.22 122.19 130.33 132.36

30.18 52.34 62.19 75.04 95.45 96.09 98.22 122.19 130.33 132.36

24.44 43.02 60.64 62.56 83.02 93.08 93.71 104.31 128.89 123.10

24.44 43.02 60.64 62.56 83.02 93.08 93.71 104.31 128.89 123.10

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

11

Fig. 5. Some modal shapes of the brain.

5.3.1. Influence of the CSF compressibility The influence of the sound velocity c0 on the eigenfrequencies is obtained in Table 6. The same simplification as previously holds for the acoustical problem. 5.3.2. Influence of the CSF thickness: Effect of confinement Let us now consider the influence of the CSF thickness e on the natural frequencies. To that purpose, we took the inner radius R2 of X2 to 0:0969 [m] and we varied the radius R1 of X1. According to the Table 7, the eigenfrequencies are strongly dependent on the CSF thickness. 5.3.3. Modal shapes Six modal shapes are reported in Fig. 5. We can see that the modal shapes and eigenfrequencies are similar to the previous model. 6. Conclusions The knowledge of the natural frequencies may be very important either in understanding Trauma Brain Injury during impacted head (Will & Berg, 1996); or neuroimaging (Hiltunen et al., 2006; Will & Berg, 2007). Very simplified models have been developed in the present study for analyzing the influence of the cerebrospinal fluid phase on the frequency spectrum on the brain-CSF-skull system. Some concluding remarks could be drawn:  It seems that the incompressibility condition is justified by simulating the sound velocity c0 of the fluid. CSF compressibility might be neglected during the modal analysis. Indeed, within the interval c0 2 [300, 1450] [m  s1] of the sound velocity, an incompressibility assumption can be used for fluid–structure in the acoustical approximation. However, compressibility could not be neglected in a dynamical simulation as Trauma Brain Injury.  Modal shapes of the rigid skull are the same as those of the elastic skull case. This is probably due to the high stiffness of the skull compared to that of the brain and due to the concentric shape of the domain, which automatically diagonalizes the mass matrix. In such a case, it is therefore possible to replace the ‘‘stiff skull’’ by a infinitely rigid wall. But in all cases, it seems not a very realistic model for analyzing the head impact.  It is shown that the CSF thickness has a strong influence on the natural frequencies. Nevertheless, results seem to be in contradiction with the intuition. Indeed, increasing the fluid thickness tends to decrease the frequencies since the added mass increases accordingly. This apparent contradiction can be explained with the greater contribution of the skull stiffness compared to the added mass of the CSF.  It is clearly stated that analytical and numerical results are in good agreement in this study. By the way, this validates the use of Comsol Multiphysics software for further dynamical analysis of the brain-CSF-skull system. Indeed, care should be taken when using finite element models for fluid–structure interaction simulation since non-symmetrical system may induce spurious modes (Bermudez, Duran, Muschieti, Rodrigez, & Solomin, 1995). These spurious modes are systematically eliminated in the algorithm used in the software.

12

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

In this paper the Helmholtz potential method is proposed to determine the natural coupled frequencies with a good accuracy. Those naturals coupled frequencies allows one to use the modal projection method for dynamic impact problem associated. Indeed, the knowledge of the first two frequencies are most important in determining of the Rayleigh damping factor (Yang, 2005). Appendix A The elements of the M matrix are written as

a11 ¼ q0 x2 J 0nþ1 ða1 R1 Þ; a13 ¼ J 0nþ1 ðdR1 Þ; a14 ¼ Y 0nþ1 ðdR1 Þ 2 2 2 h i a12 ¼ q0 x2 R1 J 00nþ1 ðb1 R1 Þ þ 2J 0nþ1 ðb1 R1 Þ þ R1 b21 J nþ1 ðb1 R1 Þ 2

2

2

a21 ¼ 2l1 J 00nþ1 ða1 R1 Þ  k1 a21 J nþ1 ða1 R1 Þ 2 2 h i 00 2 0 2 ðb R Þ þ 3J a22 ¼ 2l1 R1 J 000 1 1 1 nþ nþ1 ðb1 R1 Þ þ R1 b1 J nþ1 ðb1 R1 Þ þ b1 J nþ1 ðb1 R1 Þ 2

2

a23 ¼ J nþ1 ðdR1 Þ;

2

2

a24 ¼ Y nþ1 ðdR1 Þ

2

2

a31 ¼ R1 J 0nþ1 ða1 R1 Þ  J nþ1 ða1 R1 Þ 2

2

a32 ¼ R21 J 00nþ1 ðb1 R1 Þ þ R1 J 0nþ1 ðb1 R1 Þ þ 2

2

! b2 R21  1 J nþ1 ðb1 R1 Þ 2 2

a43 ¼ q0 x2 J 0nþ1 ða2 R2 Þ; a44 ¼ q0 x2 Y 0nþ1 ða2 R2 Þ 2 2 h i a45 ¼ q0 x2 R2 J 00nþ1 ðb2 R2 Þ þ 2J 0nþ1 ðb2 R2 Þ þ R2 b22 J nþ1 ðb2 R2 Þ 2 2 2 h i a46 ¼ q0 x2 R2 Y 00nþ1 ðb2 R2 Þ þ 2Y 0nþ1 ðb2 R2 Þ þ R2 b22 Y nþ1 ðb2 R2 Þ 2

a47 ¼

J 0nþ1 ðdR2 Þ; 2

a53 ¼ 2l

a48 ¼

00 2 J nþ12 ð 2 R2 Þ

a

2

2

Y 0nþ1 ðdR2 Þ 2

 k2 a22 J nþ1 ða2 R2 Þ 2

a54 ¼ 2l2 Y 00nþ1 ða2 R2 Þ  k2 a22 Y nþ1 ða2 R2 Þ 2 2 h i 00 2 0 2 ðb R Þ þ 3J ðb a55 ¼ 2l2 R2 J 000 2 2 2 R2 Þ þ R2 b2 J nþ12 ðb2 R2 Þ þ b2 J nþ12 ðb2 R2 Þ nþ12 nþ12 h i 00 2 0 2 a56 ¼ 2l2 R2 Y 000 nþ1 ðb2 R2 Þ þ 3Y nþ1 ðb2 R2 Þ þ R2 b2 Y nþ1 ðb2 R2 Þ þ b2 Y nþ1 ðb2 R2 Þ 2

a57 ¼ J nþ1 ðdR2 Þ;

2

2

a58 ¼ Y nþ1 ðdR2 Þ

2

a65 ¼

2

2

R2 J 0nþ1 ð 2 R2 Þ 2

a

 J nþ1 ða2 R2 Þ 2

a66 ¼ R2 Y 0nþ1 ða2 R2 Þ  Y nþ1 ða2 R2 Þ 2

2

! b2 R22  1 J nþ1 ðb2 R2 Þ 2 2 ! 2 2 b R2 ¼ R22 Y 00nþ1 ðb2 R2 Þ þ R2 Y 0nþ1 ðb2 R2 Þ þ  1 Y nþ1 ðb2 R2 Þ 2 2 2 2

a67 ¼ a68

R22 J 00nþ1 ðb2 R2 Þ 2

þ

R2 J 0nþ1 ðb2 R2 Þ 2

þ

a75 ¼ 2l2 J 00nþ1 ða2 R0 Þ  k2 a22 J nþ1 ða2 R0 Þ 2

2

a76 ¼ 2l2 Y 00nþ1 ða2 R0 Þ  k2 a22 Y nþ1 ða2 R0 Þ 2 2 h i 00 2 0 2 ðb R Þ þ 3J ðb a77 ¼ 2l2 R2 J 000 1 1 2 0 2 R0 Þ þ R0 b2 J nþ12 ðb2 R0 Þ þ b2 J nþ12 ðb2 R0 Þ nþ2 nþ2 h i 00 2 0 2 a78 ¼ 2l2 R2 Y 000 nþ1 ðb2 R0 Þ þ 3Y nþ1 ðb2 R0 Þ þ R0 b2 Y nþ1 ðb2 R0 Þ þ b2 Y nþ1 ðb2 R0 Þ 2

2

2

a85 ¼ R0 J 0nþ1 ða2 R0 Þ  J nþ1 ða2 R0 Þ 2

2

a86 ¼ R0 Y 0nþ1 ða2 R0 Þ  Y nþ1 ða2 R0 Þ 2

2

! b2 R20  1 J nþ1 ðb2 R0 Þ 2 2 2 2 ! b2 R20 ¼ R20 Y 00nþ1 ðb2 R0 Þ þ R0 Y 0nþ1 ðb2 R0 Þ þ  1 Y nþ1 ðb2 R0 Þ 2 2 2 2

a87 ¼ R20 J 00nþ1 ðb2 R0 Þ þ R0 J 0nþ1 ðb2 R0 Þ þ a88

2

A. El Baroudi et al. / International Journal of Engineering Science 51 (2012) 1–13

13

References Abramowitz, M., & Stegun, I. A. (1970). Handbook of Mathematical Functions. New York: Dover. Bermudez, A., Duran, R., Muschieti, M. A., Rodrigez, R., & Solomin, J. (1995). Finite element vibration analysis of fluid–solid systems without spurious modes. SIAM Journal of Numerical Analysis, 32, 1280–1295. Charalambopoulos, A., Dassios, G., Fotiadis, D., & Massalas, C. (1997). Frequency spectrum of the human head–neck system. International Journal of Engineering Science, 35, 753–768. Gong, S., Lee, H., & Lu, C. (2008). Computational simulation of the human head response to non-contact impact. Computers and Structures, 86, 758–770. Hault-Dubrulle, A. (2007). Towards a better understanding of knowledge about fluid/structure interaction phenomena inside the head subjected to a dynamical loading: Experimental and numerical studies, Ph.D. thesis, Université de Valenciennes et du Hainaut Cambrésis, 6 juillet. Hiltunen, J., Hari, R., Jousmaki, V., & Mueller, K. (2006). Quantification of mechanical vibration during diffusion tensor imaging at 3t. NeuroImage, 32, 93–103. Ho, J., & Kleiven, S. (2009). Can sulci protect the brain from traumatic injury? Journal of Biomechanics, 42, 2074–2080. Kleiven, S. (2002). Finite element modeling of the human head, Ph.D. thesis, Royal Institute of Technology. Marjoux, D., Baumgartner, D., Deck, C., & Willinger, R. (2008). Head injury prediction capability of the hic, hip, simon and ulp criteria. Accident Analysis and Prevention, 40, 1135–1148. Morand, H. J. P., & Ohayon, R. (1979). Substructure variational analysis of the vibrations of coupled fluid-structure systems. Finite element results. International Journal of Numerical Methods and Engineering, 14, 741–755. Morand, H. J. P., & Ohayon, R. (1995). Fluid Structure Interaction. John Wiley & Sons. Morse, M., & Feshbach, H. (1946). Methods of Theoretical Physics. New York: McGraw-Hill. Morse, P. M., & Ingard, K. U. (1968). Theoretical Acoustics. New York: McGraw-Hill Book Company. Roth, S., Vappou, J., Raul, J.-S., & Willinger, R. (2009). Child head injury criteria investigation through numerical simulation of real world trauma. Computer Methods and Programs in Biomedicine, 93, 32–45. Verschuerena, P., Delyeb, H., Depreitere, B., Lierde, C. V., Haex, B., Berckmans, D., et al (2007). A new test set-up for skull fracture characterisation. Journal of Biomechanics, 40, 3389–3396. Will, U., & Berg, E. (1996). Modal and temporal analysis of head mathematical models, Traumatic brain injury: Bioscience and mechanics. Mary Ann Liebert Inc Larchmont, 265–276. Will, U., & Berg, E. (2007). Brain wave synchronization and entrainment to periodic acoustic stimuli. Neuroscience Letters, 424, 55–60. Yang, B. (2005). Stress, Strain, and Structural Dynamics. Elsevier Academic Press. Zienkiewich, O. C., & Taylor, R. L. (1989). The Finite Element Method (Vol. 2). McGraw-Hill. Zonga, Z., Leeb, H., & Lu, C. (2006). A three-dimensional human head finite element model and power flow in a human head subject to impact loading. Journal of Biomechanics, 39, 284–292.