Free vibrations of the viscoelastic human skull

PII:

Int. J. Engng Sci. Vol. 36, No. 5/6, pp. 565±576, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0020-7225/98 $19.00 + 0.00 S0020-7225(97)00087-6

FREE VIBRATIONS OF THE VISCOELASTIC HUMAN SKULL A. CHARALAMBOPOULOS Institute of Chemical Engineering and High Temperature Chemical Processes, GR 26500 Patras, Greece D. I. FOTIADIS Department of Computer Science, University of Ioannina, GR 45110 Ioannina, Greece C. V. MASSALAS Department of Mathematics, University of Ioannina, GR 45110 Ioannina, Greece (Communicated by

E. S. S° UHUBIÇ )

AbstractÐIn this work we deal with the free vibrations of a viscoelastic skull. The analysis is based on the three-dimensional theory of viscoelasticity and the representation of the displacement ®eld in terms of the Navier eigenvectors. The frequency equation was solved numerically and results are presented for the eigenfrequency and the attenuation spectra. # 1998 Elsevier Science Ltd. All rights reserved

1. INTRODUCTION

Several investigators have studied the frequency spectrum of the human skull both theoretically and experimentally. It is well known that a knowledge of the frequency spectrum can be applied for better understanding of head injuries and in the diagnosis of head diseases. One subject which attracted the attention of many researchers is the free vibration analysis of the human skull. From the experimental point of view, Khalil et al. [1] have presented experimental measurements for the resonance frequencies of the human skull for two kinds of dry skulls and extrapolated their results to living skulls taking into account known and estimated di€erences in mechanical properties. HaÊkansson et al. [2] investigated the free damped natural frequencies of the human skull in vivo and they are the only ones to our knowledge who reported damping coecients. In previous communications we have studied the dynamic characteristics of the human dry skull [3], as well as the e€ect of various geometrical parameters [4, 5]. We restricted our study to the sensitivity of the frequency spectrum from the various parameters assuming that the skull is a linear isotropic elastic material. Several other researchers have paid attention to the dissipative material behaviour of the human head system. Misra [6] studied the steady-state response of a prolate spheroidal shell made of a linear viscoelastic solid containing a viscoelastic ¯uid due to an axisymmetric load varying harmonically with time. Misra and Chakravarty [7] have shown that the dissipative material of the skull as well as the brain have a signi®cant e€ect on the frequency spectrum of the freely vibrating cranial system. The three-dimensional equations of linear viscoelasticity have been used by Hickling and Wenner [8] to describe the behaviour of both the brain and the skull to an axisymmetric impact. However, in the above works no clear statement has been made about the cause of dissipative behaviour of the human head system which matches qualitatively the results of HaÊkansson et al. [2]. In the present work we deal with the role of material behaviour on the frequency spectrum of the human skull. The mathematical modelling is based on the three-dimensional theory of elasticity and the approximation of the skull by a linear viscoelastic material occupying the region bounded by two concentric spheres. We assume steady vibrations with angular frequency o1 and attenuation o2 and, taking advantage of the Fourier transform properties, we obtain the variables which replace LameÁ's constants of the elastic material case. The mathematical analysis 565

566

A. CHARALAMBOPOULOS et al.

is based on the representation of the displacement ®eld in terms of the Navier eigenvectors. The frequency equation is constructed by imposing the satisfaction of the boundary conditions and it is solved numerically. Due to the existence of complex terms this leads to zero-crossing ®nding m of two functions in the complex plane. From the analysis adopted we obtained o m 1 and o 2 , m = 1,2,3, . . . for real and complex  (Poisson ratio) and we also performed a sensitivity analysis of our results for various parameters entering the problem.

2. PROBLEM FORMULATION

We consider two concentric spheres with centre 0, radii r0, r1 and surfaces S0, S1 respectively (Fig. 1). The region V between the two spheres is assumed to be ®lled with a viscoelastic material obeying the following constitutive equation:  …t  @ 1@ 1 T t…r,t† ˆ ‰ru…r,t† ‡ …ru…r,t†† Š ‡ G2 …t ÿ t† ÿ G1 …t ÿ t† …r
u…r,t††I dt G1 …t ÿ t† 2 @t 3 @t ÿ1 …1† where u stands for the displacement vector of the viscoelastic medium, t…r,t† the stress tensor, I the identity tensor and G1(x), G2(x) are independent functions, having zero values for negative argument and de®ning the viscoelastic properties of the medium. The kinematic behaviour of the system is described by the equation ui , i,j ˆ 1,2,3 tij,j ˆ r

…2†

where indices indicate components of the corresponding tensors while indices after the comma indicate di€erentiation with respect to the corresponding Cartesian coordinate. Assuming steady vibrations of the system under consideration, with angular frequency o1 and attenuation (due to viscous properties) o2, we apply Fourier transform analysis to the problem de®ning … ‡1 ^ u…r,t†eÿio t dt …3† u…r,o † ˆ ÿ1

t…r,o † ˆ

… ‡1 ÿ1

t…r,o †eÿio t dt

and

Fig. 1. Geometry.

…4†

Free vibrations of the human skull

gj …o † ˆ

… ‡1 0

Gj …t†eÿio t dt, j ˆ 1,2

567

…5†

with o = o1+io2. Combining equations (1)±(5) and taking advantage of Fourier transform properties we obtain 1 1 io g1 …o †‰r2 ^ u…r,o † ‡ r…r
^ u…r,o ††Š ‡ io ‰g2 …o † ÿ g1 …o †Šr…r
^u…r,o †† ‡ ro 2 ^u…r,o † ˆ 0 2 3 or u ‡ ‰G * …o † ‡ l* …o †Šr…r
^u† ‡ ro 2 ^u ˆ 0 G * …o †r2 ^

…6†

1 G * …o † ˆ io g1 …o † 2

…7†

1 l* …o † ˆ io ‰g2 …o † ÿ g1 …o †Š 3

…8†

where

are the frequency functions replacing LameÁ's constants m, l, respectively, in the case of elastic materials. It is convenient to express G*(o), l*(o) in terms of the generalised Young's modulus E…o † ˆ io g1 …o †

3g2 …o † g1 …o † ‡ 2g2 …o †

and the Poisson ratio …o † ˆ ‰g2 …o † ÿ g1 …o †Š=‰ g1 …o † ‡ 2g2 …o †Š through the relations G * …o † ˆ ‰E…o †Š=‰2…1 ‡ …o ††Š and l* …o † ˆ ‰E…o †…o †Š=‰…1 ‡ …o ††…1 ÿ 2…o ††Š. So, in the frequency domain, the viscoelastic properties of the medium under consideration are completely described by the functions E(o), (o). Inserting dimensionless variables, the equation of elasticity, equation (6), takes the form u…r 0 † ‡ ‰G * …o † ‡ l* …o †Šr 0 …r 0
^u…r 0 †† ‡ O2 ^u…r 0 † ˆ 0 G~ * …o †r 02 ^

…9†

where r 0 ˆ r=a …a ˆ r1 †, r 0 ˆ ar, G~ * …o † ˆ G * …o †=G * …o † ˆ 1 p l~ * …o † ˆ l* …o †=G * …o † ˆ 2…o †=…1 ÿ 2…o ††, O ˆ o a 2r…1 ‡ …o ††=E…o †: In equation (9) the indication of the dependence of u^…r,o † with respect to o has, for simplicity, been suppressed. The physical characteristics of the system enter the mathematical formulation of the problem through the boundary conditions that are described on surfaces S0, S1. More precisely, these two surfaces are stress free, that is Tu…r 0 † ˆ 0, r 0 2 S0 , S1

…10†

T ˆ 2G~ * …o †^r 0
r 0 ‡ l~ * …o †^r…r 0
† ‡ G~ * …o †^r  r 0 

…11†

where

stands for the dimensionless surface stress operator in the medium V, ^r is the unit normal vector on surfaces S0, S1 and m0 ˆ

G * …o † ˆ1 G * …o †

…12†

568

A. CHARALAMBOPOULOS et al.

l0 ˆ

l* …o † 2…o † : ˆ G * …o † 1 ÿ 2…o †

…13†

We note that the system of equations (6)±(10) constitutes a well-posed ``corresponding'' problem. With this term is meant the identical problem except that the body concerned is viscoelastic instead of elastic. The method adopted for the solution of this problem is based on the representation of the displacement ®eld u(r) in terms of the Navier eigenvectors [9]. As is well known, Navier eigenvectors are created through Helmholtz decomposition and constitute a complete set of vector functions in the space of solutions of the time-independent Navier equation. The Navier eigenvectors have the following form: 0 _ ln …Or 0 †Pm r† ‡ Lm,l n …r † ˆ g n …^ 0 Mm,l n …r † ˆ 0 Nm,l n …r † ˆ n…n ‡ 1†

p gl …Or 0 † m n…n ‡ 1† n 0 Bn …^r† Or

p l 0 0 m n…n ‡ 1†gn …k s r †Cn …^r†

p gln …k 0s r 0 † m gl …k 0 r 0 † Pn …^r† ‡ n…n ‡ 1†g_ ln …k 0s r 0 † n 0 s 0 Bm r† n …^ 0 0 k sr k sr

…14† …15† …16†

n ˆ 0,1,2, . . . ; m ˆ ÿn, ÿ n ‡ 1, . . . ,n ÿ 1,n; l ˆ 1,2 where k 0s ˆ gln …x†

 ˆ

Jn …x† Yn …x†

O G~ *1=2

l ˆ 1 …spherical Bessel function of order n† l ˆ 2 …spherical Neumann function of order n†

…17†

…18†

where g_ ln …x† stands for the derivative of gln …x† with respect to its argument. Finally, the functions r†, Bm r†, Cm r† de®ned on the unit sphere are the vector spherical harmonics, introduced by Pm n …^ n …^ n …^ Hansen [9] and given by the relations r† ˆ ^rY m r† Pm n …^ n …^

…19†

r† Bm n …^

  1 @ j^ @ ^ ˆ p W ‡ r† Ym n …^ n…n ‡ 1† @W sin W @j

…20†

r† Cm n …^

 ^  W @ 1 @ ÿ j^ Ym ˆ p r† n …^ @W n…n ‡ 1† sin W @j

…21†

imj r† ˆ P m are the spherical harmonics and P m where Y m n …^ n …cos W†e n …cos W† the well-known Legendre functions. Consequently, the displacement ®eld u(r) has the representation

u…r† ˆ

‡1 X n X 2  X m,l m,l m,l m,l m,l am,l n Ln …r† ‡ bn Mn …r† ‡ gn Nn …r† :

…22†

nˆ0 mˆÿn lˆ1

Now, the problem of the determination of u(r) is transferred to the determination of the coef®cients of the expansion (22) in Navier eigenfunctions. The expression (22) satis®es Navier equation. What remains to do is to force this expansion to satisfy the boundary conditions (10). Consequently, applying the stress operator T on the expansion (22), expressing the functions m m,l m,l m m TLm,l n , TMn , TNn in terms of the orthogonal set of vector harmonics Pn , Bn , Cn and taking

Free vibrations of the human skull

569

advantage of the independence of the last functions, we conclude that for every speci®c pair of integers (n,m), with vmv Rn, the six coecients involved in expansion (22) satisfy a linear homogeneous system of six equations of the form m Dm n xn ˆ 0

…23†

where m,1 m,2 m,1 m,2 m,1 m,2 T xm n ˆ ‰an ,an ,bn ,bn ,gn ,gn Š

and

2

A1n …r 01 † 6 0 6 1 0 6 B …r † 6 1n 01 ˆ Dm n 6 A …r † 6 n 0 4 0 B 1n …r 00 † where

A2n …r 01 † 0 0 0 C 1n …r 01 † C 2n …r 01 † B 2n …r 01 † 0 0 A2n …r 00 † 0 0 0 C 1n …r 00 † C 2n …r 00 † 2 0 B n …r 0 † 0 0

3 D1n …r 01 † D2n …r 01 † 0 0 7 7 1 0 2 0 7 E n …r 1 † E n …r 1 † 7 D1n …r 00 † D2n …r 00 † 7 7 0 0 5 E 1n …r 00 † E 2n …r 00 †

    4 n…n ‡ 1† l 0 0 l 0 …Or † ‡ l Og …Or † g Aln …r 0 † ˆ ÿ 0 g_ ln …Or 0 † ‡ 2O 1 ÿ n n r O2 r 02  p 1 gl …Or 0 † B ln …r 0 † ˆ 2 n…n ‡ 1† g_ ln …Or 0 † 0 ÿ n 2 02 r O r C ln …r 0 † ˆ

  p 1 n…n ‡ 1† …2 ‡ l 0 †1=2 Og_ ln ……2 ‡ l 0 †1=2 Or 0 † ÿ 0 gln ……2 ‡ l 0 †1=2 Or 0 † r

Dln …r 0 †



g_ l ……2 ‡ l 0 †1=2 Or 0 † gln ……2 ‡ l 0 †1=2 Or 0 † ˆ 2n…n ‡ 1† n ÿ r0 …2 ‡ l 0 †1=2 Or 0



 p g_ l ……2 ‡ l 0 †1=2 Or 0 † ÿ …2 ‡ l 0 †1=2 Ogln ……2 ‡ l 0 †1=2 Or 0 † E ln …r 0 † ˆ n…n ‡ 1† ÿ 2 n r0  n…n ‡ 1† ÿ 1 l 0 1=2 0 g ……2 ‡ l † Or † : ‡2 n …2 ‡ l 0 †1=2 Or 02

Fig. 2. Solution of equations (27a)±(b) in the complex plane (O1,O2) for r=0.25, r1=0.082 m, m r0=0.076 m and n = 2. (1: Re‰det…Dm n …O1 ,O2 ††Š ˆ 0; 2: Im‰det…Dn …O1 ,O2 ††Š ˆ 0).

Fig. 3. Solution of equations (27a)±(b) for  = 1+i2 (1=constant) in the complex plane (O1,O2) for r1=0.082 m, r0=0.076 m and n = 2. (A: 2=0.001; B: 2=0.01; C: 2=0.1; D: 2=0.2).

570 A. CHARALAMBOPOULOS et al.

Free vibrations of the human skull

571

Table 1. Frequency spectrum for  = 0.25 and r1=0.082 m, r0=0.076 m O [3]

Elastic skull

0.7083 0.8522 0.9486 1.063 1.1971 1.2184 1.4210 1.5479 1.6695 1.8924 2.5324 2.6253

o 2709.282 3259.638 2628.364 4065.941 4578.870 4660.341 5435.280 5920.669 6385.785 7238.370 9686.350 10041.689

o1

Viscoelastic skull

2709.291 3259.711 3628.446 4066.033 4578.972 4660.446 5435.402 5920.802 6385.928 7238.533 9686.567 10041.915

o2 18.159 21.849 24.320 27.253 30.691 31.237 36.431 39.685 42.802 48.517 64.925 67.307

In order for the system (23) to have a non-trivial solution, the following condition must be satis®ed: det…Dm n …O†† ˆ 0:

…24†

This condition provides the characteristic (frequency) equation from the roots of which we m obtain the eigenfrequencies om 1 and the attenuation coecients o 2 , m = 1,2,3, . . ., of the system under discussion.

3. SOLUTION ± NUMERICAL RESULTS m The frequency equation (24) is treated numerically in order for its solution Om ˆ Om 1 ‡ iO2 to be obtained. The spherical Bessel functions and their derivatives for complex arguments are computed using backward recurrence. The spherical Neumann functions and their derivatives are computed using a mixed recurrence because its value decreases ®rst and increases as its order increases [10]. A complex LU decomposition routine is used along with a determinant computation routine which results in the following: m Re‰det…Dm n …O††Š ‡ i Im‰det…Dn …O††Š ˆ 0

…25†

m Re‰det…Dm n …O1 ,O2 ††Š ‡ i Im‰det…Dn …O1 ,O2 ††Š ˆ 0

…26†

or

which is equivalent to Re‰det…Dm n …O1 ,O2 ††Š ˆ 0

…27a†

Im‰det…Dm n …O1 ,O2 ††Š ˆ 0:

…27b†

The solution of the system of equations (27) is obtained by a bisection grid method. A rectangular grid is established in the (O1,O2) plane and bisection searches for zero crossings are done along each of the equally spaced grid lines. The spacing chosen for the present solution is 10ÿ3, the same in both axes, and each grid line involves a separate one-dimensional bisection search with accuracy 10ÿ6. Such a solution, obtained for  = 0.25, is graphically shown in Fig. 2 and the crossings of zero-crossing curves are shown to be on the real axis (O2=0). Similar zerocrossing curves for complex Poisson ratios (0.01R 2R0.2) are given in Fig. 3. The crossing of curves gives the complex solution (O = O1+iO2) of the system of equations (27). To ®nd the solution a re®ned grid search method along with re®nement of search in those squares where there are zero crossings of the real and imaginary part of the determinant on their boundaries is used with computation accuracy 10ÿ8 [11].

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

0.0278

1.5946

ÿO2

0

n=0

0

O1

n=1

1.9495

0

O1

0.1252

0

ÿO2

2.7003

0.0904

0.1530

1.2274

ÿO2 0.0792

n=2

0.7243

O1

n=3

3.6207

1.9407

0.8776

O1

0.1607

0.2418

0.0853

ÿO2

n=4

4.5881

2.6037

0.9770

O1

0.2289

0.3245

0.0879

ÿO2

n=5

5.5710

3.2474

1.0946

O1

0.2950

0.4047

0.0922

ÿO2

n=6

6.5603

3.8813

1.2539

O1

0.3603

0.4837

0.9995

ÿO2

…n† Table 2. Eigenfrequency coecients O…n† ˆ O…n† 1 ‡ iO2 , for  = 0.25 + i0.1, r1=0.082 m, r0=0.076 m and n = 0,1,...,8

n=7

7.5523

4.5097

1.4619

O1

0.4256

0.5620

0.1010

ÿO2

n=8

8.5455

5.1344

1.7175

O1

0.4913

0.6398

0.1252

ÿO2

572 A. CHARALAMBOPOULOS et al.

Free vibrations of the human skull

573

Table 3. Eigenfrequency and attenuation spectra for  = 0.25 + i0.1, r1=0.082 m and r0=0.076 m Elastic skull [3]

O 0.7083 0.8522 0.9486 1.0631 1.1971 1.2184 1.4210 1.5479 1.6695 1.8924 1.8924 2.5324 2.6253

(n = 2) (n = 3) (n = 4) (n = 5) (n = 2) (n = 6) (n = 7) (n = 0) (n = 8) (n = 1) (n = 3) (n = 4) (n = 2)

o

O1

2784.727 3355.030 3735.108 4184.679 4692.286 4793.594 5588.646 6096.023 6565.997 7452.738 7419.119 9953.767 10323.095

o1

0.7243 0.8776 0.9770 1.0946 1.2274 1.2539 1.4619 1.5946 1.7175 1.9495 1.9407 2.6037 2.7003

Viscoelastic skull

2711.926 3262.434 3628.634 4062.345 4579.019 4650.483 5419.097 5866.431 6364.619 7172.260 7240.025 9713.487 9955.172

O2

o2

ÿ0.0792 ÿ0.0853 ÿ0.0879 ÿ0.0922 ÿ0.1530 ÿ0.0995 ÿ0.0110 ÿ0.0278 ÿ0.1252 ÿ0.0345 ÿ0.2418 ÿ0.3245 ÿ0.0904

60.662 110.256 148.535 189.571 30.675 239.799 300.132 667.104 369.188 813.968 48.526 65.107 970.888

Numerical results were obtained for material properties analogous to those proposed in the literature [8] for the viscoelastic human skull, that is: E = 6.895  1010+i0.092  103 N/m2;  = 0.25 and  = 0.25 + i2, 2 variable; r = 2.132  103 kg/m3; r1=0.082 m, r0$[0.040,0.076 m]. The solution of equations (27a)±(b) for real  corresponds to that of the elastic case which has been presented in ref. [3] and the results are cited in Table 1. In the case of complex ,  = 0.25 + i0.1 and di€erent n (where n is the order of Bessel functions), the roots of equations (27a)±(b) are given in Table 2. From the results obtained we observe that the ordering of the roots is similar to that of the elastic case and the degeneracy for n = 1 disappears. In Table 3 the results for  = 0.25 and  = 0.25 + i0.1 are presented. It is observed that the increase in the imaginary part (o2) of the ®rst o for each n follows the decrement of the real part (o1). The dependence of o1 and o2 on 2 (1 is constant) is shown in Table 4 and graphically, for the ®rst four eigenfrequencies, in Fig. 4. The variations of o1 are small but those of o2 are substantially large. m The variations of o m 1 , o 2 , m = 1,2,3,4 with h/r1 (h = r1ÿr0) are shown in Fig. 5. It is m m observed that o 1 , o 2 increase with the increment of the skull thickness (for 0 Rh/r1R0.040), while for higher values of h/r1 a decrement is observed for o m 2.

Table 4. Eigenfrequency and attenuation spectra for 2 $ [0.0,0.2] and r1=0.082 m, r0=0.076 m  = 0.25 + i0.01 o1 o2 2709.384 3259.797 3628.370 4066.343 4579.010 4660.153 5435.448 5919.867 6388.110 7237.320 7240.027 9713.626 10040.479

 = 0.25 + i0.05 o1 o2

22.430 2728.532 30.720 3284.036 36.791 3655.601 43.579 4096.624 30.693 4607.786 52.241 4694.414 63.018 5474.888 103.033 5906.034 75.736 6434.002 125.798 7293.849 48.528 7295.524 57.460 9774.518 158.562 10114.546

39.479 66.158 86.613 108.755 30.692 136.072 169.073 355.641 207.038 433.976 48.529 65.099 522.437

 = 0.25 + i0.1 o1 o2 2711.926 3262.434 3628.634 4062.345 4579.019 4650.483 5419.097 5866.431 6364.619 7172.260 7240.025 9713.487 9955.172

60.662 110.256 148.535 189.571 30.675 239.799 300.132 667.104 369.188 813.968 48.526 65.107 970.888

 = 0.25 + i0.15 o1 o2 2715.282 3266.196 3629.611 4058.324 4579.010 4639.610 5400.230 5803.093 6337.188 7095.605 7240.025 9713.469 9854.097

81.527 153.956 209.747 269.092 30.691 341.427 428.193 969.904 527.415 1183.607 48.528 65.105 1406.711

 = 0.25 + i0.2 o1 o2 2719.995 3271.605 3631.354 4053.349 4579.011 4625.411 5375.302 5726.791 6300.660 7001.674 7240.027 9713.488 9718.613

101.955 197.043 269.933 346.764 30.695 440.309 552.007 1244.055 680.074 1499.185 48.527 65.110 1824.672

Fig. 4. Variation of o n1 , o n2 , n = 1,2,3,4 with 2 for r1=0.082 m, r0=0.076 m.

Free vibrations of the human skull

575

Fig. 5. Variation of o n1 and o n2 with h/r1 for  = 0.25 + i0.1 and r1=0.082 m.

4. CONCLUSIONS

In the previous analysis we presented the in¯uence of the viscoelastic nature of the human skull material on its frequency spectrum. From the results obtained we were led to the conclusion that for  $ R and E = (
) + i(
) the frequency spectrum of the skull remains almost that of the elastic case and the attenuation is immaterial. For  and E complex the ordering of the eigenfrequencies remains but there is a decrement of o1 and a signi®cant increment of o2. We also observed that o m 1 increases with the increment of the thickness of the skull while there is a critical h beyond which the attenuation o m 2 decreases. AcknowledgementsÐThe present work forms part of the project ``New Systems for Early Medical Diagnosis and Biotechnological Applications'' which is supported by the Greek General Secretariat for Research and Technology through the EU funded R&D Program EPET II.

REFERENCES 1. Khalil, T. B., Viano, D. C. and Smith, D. L., J. Sound and Vibration, 1979, 63, 351. 2. HaÊkansson, B., Brandt, A. and Carlsson, P., J. Acoust. Soc. Am., 1994, 95(3), 1474. 3. Charalambopoulos, A., Dassios, G., Fotiadis, D. I., Kostopoulos, V. and Massalas, C. V., Int. J. Engng. Sci., 1996, 34, 1339. 4. Charalambopoulos, A., Fotiadis, D. I. and Massalas, C. V., Int. J. Engng. Sci., in press. 5. Charalambopoulos, A., Fotiadis, D. I. and Massalas, C. V., Acta Mechanica, in press. 6. Misra, J. C., Ingenieur ± Archiv, 1978, 47, 11.

576 7. 8. 9. 10. 11.

A. CHARALAMBOPOULOS et al. Misra, J. C. and Chakravarty, S., J. Biomechanics, 1982, 15, 635. Hickling, R. and Wenner, M. L., J. Biomechanics, 1973, 6, 115. Hansen, W. W., Phys. Rev., 1935, 47, 139. Zhang, S. and Jin, J., Computation of Special Functions. Wiley-Interscience, New York, 1996. Hostetter, G. H., Santina, M. S. and Caprio-Montalvo, P. D', Analytical, Numerical, and Computational Methods for Science and Engineering. Prentice Hall, Englewood Cli€s, NJ, 1991. (Received 19 June 1997)