MODELING OF LOW SHAFT SPEED BEARING FAULTS FOR CONDITION MONITORING

MODELING OF LOW SHAFT SPEED BEARING FAULTS FOR CONDITION MONITORING

Mechanical Systems and Signal Processing (1998) 12(3), 415–426 MODELING OF LOW SHAFT SPEED BEARING FAULTS FOR CONDITION MONITORING Y.-F. W*  P...

286KB Sizes 4 Downloads 48 Views

Mechanical Systems and Signal Processing (1998) 12(3), 415–426

MODELING OF LOW SHAFT SPEED BEARING FAULTS FOR CONDITION MONITORING Y.-F. W*  P. J. K CRASys, Department of Systems Engineering, ANU, Canberra, ACT 0200, Australia (Received February 1997, accepted after revisions November 1997) A general model of faulty rolling element bearing vibration signals is established. The efficacy of the envelope-autocorrelation technique for condition monitoring of such bearings leads us to derive the envelope-autocorrelation function of the model in the case of very low shaft speed. For application to the bearing condition monitoring, a simplified model is given. Finally the comparison of simulated data using the model with in situ data confirms the satisfactory performance of the model and the effectiveness of envelope-autocorrelation as a fault detection technique for low shaft speed bearings. 7 1998 Academic Press Limited

1. INTRODUCTION

The condition monitoring of rotating machinery is important in terms of system maintenance and process automation, especially for the condition monitoring of rolling element bearings which are the most common components in industrial rotating machinery. They are used in industries from agriculture to aerospace, in equipment as diverse as descaler pinch rolls to the space shuttle main engine turbopumps. For this reason, a variety of bearing fault detection techniques have been proposed; but as far as we know very few articles have reported the theoretical modeling of bearing faults with mathematical descriptions. Braun [1] gave a theoretical model of the rolling element bearing with multipoint bearing defects. In his model, the repetitious impulse responses induced by bearing defects and the modulation by the bearing structural vibration were introduced. McFadden and Smith [2] improved the model by taking into account the impulse series, the modulation of the periodic signal caused by non-uniform load distribution and of the vibration transmission of rolling element bearing, as well as the exponential decay of vibration. Wang and Harrap [3] improved the model further by considering the impulse series, the modulation of the periodic signal from non-uniform load distribution and of the first bearing vibration mode with mathematical descriptions. They introduce the envelope autocorrelation analysis for both simulated signals and experimental data. Both [2] and [3] consider the modulation of cage frequency caused by a non-uniform load distribution. If the fault occurs in a rolling element, such as a ball or roller, then the cage frequency is what needs to be considered. If the fault occurs elsewhere, then due to the non-uniform load distribution, the period of the signal will be different [4, 5]. * On leave from the Department of Mechanical Engineering, Zhejiang University of Technology, Hangzhou, 310014, P. R. China. Present address: Department of Aerospace and Mechanical Engineering, University College, UNSW, Australian Defence Force Academy, Canberra, ACT 2600, Australia. 0888–3270/98/030415 + 12 $25.00/0/pg970149

7 1998 Academic Press Limited

416

.-.   . . 

This paper is motivated by the poor fit of previous signal models to an industrial case examined by the authors: condition monitoring of descaler pinch rolls in a steel mill [6]. Several attempts at standard (and not so standard) analysis of accelerometer signals had been unable to predict failure in the rollers—failure which occurred rapidly and, in some cases, catastrophically. The remainder of this paper is organised as follows. The next section presents our modified model for rolling element bearing fault signals. Section 3 derives the envelope-autocorrelation of this signal model, and also give a model simplification for application to the bearing condition monitoring. Section 4 compares signals generated using our proposed model with signals obtained in an industrial condition monitoring environment. It is concluded that the signal model presented accurately mimics the industrial example. 2. A MODEL FOR ROLLING ELEMENT BEARING FAULTS

Suppose x(t) is an original signal from a rolling element bearing with a single point fault, it can be expressed as x(t) = xf (t) · xq (t) · xbs (t) + xs (t) + n(t),

(1)

where: xf (t) is the basic impulse series produced by the fault which impacts repetitiously with another surface in the bearing (the repetition period is Tf = 1/ff and ff is the characteristic frequency of the bearing fault); xq (t) is the modulation effect caused by non-uniform load distribution of the bearings and the cyclic variation of transmission path between the fault impact site and the transducer; and xbs (t) is the bearing induced vibration. This modulating effect is caused partly by the bearing geometry, but more significantly by the position of the vibration sensor relative to the fault location. It is this element which, if not taken into account, significantly effects the efficacy of the signal model in real-life situations where optimal sensor position is not possible or perhaps undesirable. xs (t) is the machinery induced vibration. The signal xs (t) is an unwanted, structured, and predictable damped harmonic signal which is also considered to be a source of noise. n(t) is a Gaussian white noise sequence with variance sn2 . This is unpredictable measurement noise present in any practical measurement system. 2.1.     The basic impulse series due to the bearing fault is assumed to be [1–3] a

Nf

j=1

l=1

xf (t) = d0 s d(t − jTf ) 1 A0 s cos (2plff t + ffl ),

(2a)

where A0 = d0 ff is the amplitude, ffl is the initial phase of lth harmonic, and Nf is the number of harmonics induced by impulse series. The characteristic frequency of a bearing fault can be the characteristic frequencies of an inner race fault (fbpfi ), an outer race fault (fbpfo ), and a rolling element fault (fbsf ), depending on where the fault occurs. 2.2. -  The modulation effect due to non-uniform loading is assumed to be a cosine pulse series [2, 3]; therefore, one has Nq

xq (t) = s Aqk cos (2pkfq t + fqk ), k=1

(2b)

      

417

T 1 Periodic characteristics of bearings with faults under various loadings and transmission path effects Causes of periodicity Stationary loading Shaft unbalance Rolling element diameter variations Transmission path variations

Outer race faults

Inner race faults

Rolling element faults

No effect fr fc

fr No effect fr − fc

fc fr − fc No effect

No effect

fr

fc fbsf /2

fr , Shaft rotational frequency; fc , cage frequency; fbsf , rolling element fault of frequency.

which is a periodic function with the period of Tq (fq = 1/Tq ). The amplitude and initial phase of the kth harmonic are Aqk and fqk , respectively, and Nq is the number of harmonics. The periodicity fq depends on the elements of bearing fault as shown in Table 1 [4, 5]. 2.3. -  If the vibration induced by the bearing is assumed to be the superposition of all vibrational modes, then Nbs

xbs (t) = s Absj e−abs t cos (2pfbsj t + fbsj ), j

(2c)

j=1

where fbsj , Absj , fbsj , and absj are the resonance frequency, amplitude, initial phase, and damping factor of jth vibrational mode of the bearing, respectively, and Nbs is the modal order of bearing vibration signal. 2.4. -  If the vibration due to machinery other than the bearing in question has the same form as equation (2c), then Ns

xs (t) = s Asn e−as t cos (2pfsn t + fsn ), n

(2d)

n=1

where fsn , Asn , fsn , and asn are the resonance frequency, amplitude, initial phase, and damping factor of nth vibrational mode of the machinery system, respectively, and Ns is the modal order of machinery. Theoretically, Nq , Nf , Nbs , and Ns are infinite. The band-limited nature of most measurements means that it is reasonable to assume finite values for these in practice. 3. THE ENVELOPE-AUTOCORRELATION OF THE PROPOSED MODEL

In the very low shaft speed case, the following assumptions are made. The modulation frequencies due to non-uniform load distribution and fault frequency are much lower than the resonance frequencies from the bearing and machinery: min {fbsj , fsn }max {ff , fq },

(j = 1, 2, . . . , Nbs ; n = 1, 2, . . . , Ns ,)

(3)

.-.   . . 

418

and the basic impulse series are dominant in the vibrational signal from the faulty bearing: A0max {Absj , Aqk , Asn }. (j = 1, 2, . . . , Nbs ; k = 1, 2, . . . , Nq ; n = 1, 2, . . . , Ns )

(4)

Because of inequality (3), xs (t) in equation (1) can be filtered out easily. Furthermore it is assumed that damping and initial phase may be neglected in the derivation of the envelope-autocorrelation, because they only affect the magnitudes and do not affect the frequency components. The concepts of envelope and autocorrelation are defined as follows. The signal envelope of a real-valued signal x(t) is =z(t)=, where z(t),x(t) + ix˜ (t),

(5)

and z(t) is the analytic signal [7] associated with x(t), i is the complex operator, x˜ (t) is the Hilbert transform [8] of x(t) 1 p

x˜ (t) =

g

a

x(t) dt. t −t −a

(6)

The envelope-autocorrelation of an envelope signal =z(t)= is defined by Rzz (t),E[=z(t)= · =z(t + t)=],

(7a)

if =z(t)= can be considered ergodic, then the envelope-autocorrelation can be estimated as Rzz (t) =

1 T

g

T

=z(t)= · =z(t + t)= dt,

(7b)

0

where T is some finite time interval. For estimating convenience, here another function Fzz (t) which is related to Rzz (t) is introduced: 1 T

Fzz (t),

g

T

=z(t)=2 · =z(t + t)=2 dt.

(8)

0

Under the assumptions stated above and using the above definitions, one has x(t) 1 xbs (t)xq (t)xf (t) + n(t) =

A0 4

Nbs ,Nq ,Nf

s j,k,l = 1

4

Absj Aqk s Bp (t) + n(t),

(9a)

p=1

where B1 (t) = B1j,k,l (t) = cos [2p(fbsj + kfq + lff )t], B2 (t) = B2j,k,l (t) = cos [2p(fbsj − kfq − lff )t], B3 (t) = B3j,k,l (t) = cos [2p(fbsj + kfq − lff )t], B4 (t) = B4j,k,l (t) = cos [2p(fbsj − kfq + lff )t].

(9b)

      

419

Therefore, the function Fzz (t) can be expressed as (see Appendix A) Fzz (t) = C +

1 2

Nbs ,Nq ,Nf

(Dj,k )2{2 cos [2p(2lff )t] + 2 cos [2p(2kfq )t]

s j,k,l = 1

+ cos [2p(2lff + 2kfq )t] + cos [2p(2lff − 2kfq )t]} Nq

+ s

Nf

7

s (Ek,k 1 )2 s Hm 1 (t),

(except k = k1 + l = l1 )

(10)

m1 = 0

k,k 1 = 1 l,l 1 = 1

where H0,1,...,7 (t) = cos {2p[(l 2 l1 )ff 2 (k 2 k1 )fq ]t},

(11a)

and 2 C,C12 + 2snb ,

(A0 )2 4

C1,

Nbs ,Nq ,Nf

s

(Absj )2(Aqk )2,

j,k,l = 1

(A0 )2(Absj )2(Aqk )2 , 8

Dj,k,

N

bs (A )2 Ek,k 1, 0 Aqk Aqk1 s (Absj )2, 8 j=1

(11b)

where C, C1 , Dj,k , Ek,k 1 are constants. In the rolling element bearings, usually it is true for fq Q ff . By comparing Rzz (t) with Fzz (t), the main frequencies in the envelope-autocorrelation function Rzz (t) are the fault characteristic frequency ff and its harmonics lff (l = 2, 3, . . . , Nf ). kfq (k = 1, 2, . . . , Nq ) are the side frequencies of main frequencies, but their amplitudes are smaller than those of main frequencies. There are also some components of small amplitude with frequencies of 12 (l 2 l1 )ff 2 12 (k + k1 )fq (here (l 2 l1 ) is odd, l, l1 = 1, 2, . . . , Nf ; k, k1 = 1, 2, . . . , Nq , except l = l1 + k = k1 ) in Rzz (t). 3.1.         If the bearings are running under the condition of a light load, and the clearances of bearings are not sufficient, the periodicity of non-uniform load distribution imposes hardly any modulation on the impulse series caused by the bearing fault. In this case, it is supposed xq (t) = 1, (12) then the theoretical model of this special case is x'(t) 1 xbs (t)xf (t) + n(t),

(13)

the function F'zz (t) which is related to envelope-autocorrelation function will be as follows, N ,N

F'zz (t) = C' +

1 bs f s (D' )2 cos [2p(2lff )t] 2 j,l = 1 j Nbs ,Nf ,Nf

+ s j,l,l 1 = 1

(D'j )2{cos [2p(l + l1 )ff t] + cos [2p(l − l1 )ff t]},

(l $ l1 )

(14)

420

.-.   . . 

where C' = (C'1 )2 + 2(s'nb )2, N ,N

C'1 =

(A0 )2 bs f j 2 s (A ) , 2 j,l = 1 bs

D'j =

(A0 )2(Absj )2 , 2

(15)

where C', C'1 , and D'j are constants. In this special case, the main frequencies in the envelope-autocorrelation function are the fault characteristic frequency ff and its harmonics lff (l = 2, 3, . . . , Nf ). They are dominant in the spectrum. There are still some small-amplitude components with the frequencies of 12 (l 2 l1 )ff (here l 2 l1 is odd, l, l1 = 1, 2, . . . , Nf , l $ l1 ). Only if the bearings run under the heavy load and their clearances are big enough, will the general model be effective.

Figure 1. Simulated signal of a bearing with an inner race fault (fsamp = 5000 Hz): (a) simulated signal generated by model; (b) power spectrum of (a); (c) envelope of (a); (d) envelope-autocorrelation of (a); (e) envelope-autocorrelation power spectrum of (a).

      

421

Figure 2. Simulated signal of a bearing in good condition (fsamp = 5000 Hz): (a) simulated signal generated by model; (b) power spectrum of (a); (c) envelope of (a); (d) envelope-autocorrelation of (a); (e) envelope-autocorrelation power spectrum of (a).

4. COMPARISON OF MODEL WITH IN SITU DATA

4.1.      The in situ vibration data are the signals from a bearing housing which supports the descaler pinch rolls. The mean rotational speed of the pinch rolls is 36.06 rpm (0.601 Hz), and the shaft speed varies randomly in a small range about the mean. In our case, the clearance of bearing is small. The pinch rolls are used to move a 25-mm diameter bar, and the bearings run under a light load. A simulated signal from bearing with an inner race fault corresponding to the in situ condition is generated with the model [shown in Fig. 1(a)] with damping neglected and the sampling frequency fsamp = 5000 Hz. The modal orders of vibrational signals from the bearing, machinery system, and impulse series are Nbs = 5, Ns = 5 and Nf = 7, respectively, with a Gaussian white noise sequence (standard deviation sn = 0.0122). The amplitude of the basic impulse series is dominant. The fault frequency and its harmonics are submerged completely in the white noise in power spectrum [shown in Fig. 1(b)]. The envelope contains the white noise. However, the period Tf and the fault characteristic frequency ff are obvious in the envelope-autocorrelation and its power spectrum [shown in Fig. 1(d) and (e)].

422

.-.   . . 

4.2.      In a simulation of a bearing in good condition, no impulse series are present in the signal Fig. 2. Without other excitations, the vibrational signals are only of small amplitude. In this case, the white noise sequence is large in comparison to the vibrational signal. Here the simulated signal is generated with xf (t) = 1, with other parameters as those in Section 4.1. 4.3.   :   The in situ vibration data are from a bearing housing of the descaler pinch rolls. The rolls have a very low shaft speed with a small range of random variation about the mean. The bearings are SKF23226, double-row, spherical roller bearings. The signal in Fig. 3 is from a bearing with inner race fault. The mean shaft speed of rolls is 36.06 RPM (0.601 Hz). In the power spectrum the fault frequency and its harmonics are completely submerged in the background noise [Fig. 3(b)]. The peak frequencies f1 = 8.5 Hz, f2 = 50 Hz and their harmonics are from the background noise. These components are also

Figure 3. In situ data from a bearing with an inner race fault (fsamp = 5000 Hz): (a) in situ data from bar 1365; (b) power spectrum of (a); (c) envelope of (a); (d) envelope-autocorrelation of (a); (e) envelope-autocorrelation power spectrum of (a).

      

423

Figure 4. In situ data from a bearing in good condition (fsamp = 5000 Hz): (a) in situ data from bar 1370; (b) power spectrum of (a); (c) envelope of (a); (d) envelope-autocorrelation of (a); (e) envelope-autocorrelation power spectrum of (a).

present in the signals from healthy bearings [Fig. 4(b)]. From Fig. 3(d), there is an obvious periodicity in the envelope-autocorrelation: its period is Tf = 0.154 s, hence ff = 1/ Tf = 6.5 Hz, which coincides with the characteristic frequency from Fig. 3(e). Comparing simulated signal and in situ data from the bearing with an inner race fault, they have similar frequency components in the envelope-autocorrelation power spectra [Figs 1(e) and 3(e)], and their envelope-autocorrelations [Figs 1(d) and 3(d)] have the same basic periodicity. 4.4.   : -  Comparing the simulated signal shown in Fig. 2 and in situ data shown in Fig. 4 from the good condition bearing, each of the original signals [Figs 2(a) and 4(a)], envelope-autocorrelation [Figs 2(d) and 4(d)] and envelope-autocorrelation power spectra [Figs 2(e) and 4(e)], is very close, which proves the hypothesis. For the healthy bearing, the other excitations are small, and the white noise sequence is relatively dominant. In the power spectra of in situ data, there are some peak frequency components from background noise [Figs 3(b) and 4(b)].

424

.-.   . .  5. CONCLUSIONS

In this paper, a theoretical model of bearing fault is established. The function Fzz (t), which is related to the envelope-autocorrelation of signal in the case of very low shaft speed, is derived. According to the in situ condition, a special case of the model is given for condition monitoring. By comparing the simulated signals and in situ data, the following conclusions can be obtained. , The simulated data from the model confirms the in situ data, and the model has a satisfactory performance. , It is proved theoretically that the main frequencies in the envelope-autocorrelation of signals from faulty bearings running under a very low shaft speed are the fault characteristic frequency ff and its harmonics lff (l = 1, 2, . . . . , Nf ) [equation (14)]. , Envelope-autocorrelation [Figs 1(c) and 3(c)] from faulty bearing in the very low shaft speed case is periodical, and the period is Tf = 1/ff . However, for healthy bearings, there is no periodicity in the envelope-autocorrelation [Figs 2(d) and 4(d)]. , From the envelope-autocorrelation power spectra, the fault frequency and its harmonics are dominant, but in the power spectra, they are completely submerged in the background noise, and there are only some peak frequency components from the background noise [Figs 1(b) and 4(b)]. ACKNOWLEDGEMENTS

The authors acknowledge the funding of the activities of the Cooperative Research Centre for Robust and Adaptive Systems by the Australian Commonwealth Government under the Cooperative Research Centres Programme. REFERENCES 1. S. G. B 1980 IEEE Transactions on Sonics and Ultrasonics SU-27, 317–328. The signature analysis of sonic bearing vibrations. 2. P. D. MF and J. D. S 1984 Journal of Sound and Vibration 96, 69–82. Model for the vibration produced by a single point defect in a rolling element bearing. 3. W. Y. W and M. J. H 1996 Machine Vibration 5, 34–44. Condition monitoring of ball bearings using envelope autocorrelation technique. 4. Y.-T. S, S.-J. L 1992 Journal of Sound and Vibration 155, 75–84. On initial fault detection on a tapered roller bearing: frequency domain analysis. 5. I. M. H 1994 Noise and Vibration Conference, Perth. A review of factors affecting bearing vibration monitoring. 6. J. M. S, P. J. S, A. J. C, F. B and L. B. W 1995 Proceedings of Control ’95, Melbourne, Australia, 117–121. Spectral estimation of variable speed rotating machinery at a hot strip mill. 7. J. V 1948 Cables et Transmissions 2, 61–74. Theorie et applications de la notion de signal analytique. 8. A. V. O and R. W. S 1989 Discrete-Time Signal Processing. Englewood Cliffs, NJ: Prentice-Hall.

APPENDIX A: DERIVATION OF ENVELOPE-AUTOCORRELATION FOR THE PROPOSED MODEL

For a given signal shown in equation (9), with equation (6), one has x˜ (t) 1

A0 4

Nbs ,Nq ,Nf

s j,k,l = 1

4

Absj Aqk s B p (t) + n˜ (t), p=1

(A1a)

      

425

where B 1 (t) = B 1j,k,l (t) = sin [2p(fbsj + kfq + lff )t], B 2 (t) = B 2j,k,l (t) = sin [2p(fbsj − kfq − lff )t], B 3 (t) = B 3j,k,l (t) = sin [2p(fbsj + kfq − lff )t], B 4 (t) = B 4j,k,l (t) = sin [2p(fbsj − kfq + lff )t],

(A1b)

and n˜ (t) is the Hilbert transform of n(t). Then the envelope =z(t)= of the original signal x(t) can be calculated as follows =z(t)=2 = x 2(t) + x˜ 2(t)

=

(A0 )2 16

6

Nbs ,Nq ,Nf

s j,k,l = 1

3

4

(Absj )2(Aqk )2 s [Bp2 (t) + B p2 (t)] p=1

7

4

+2 s

s

[Bp (t)Bp 1 (t) + B p (t)B p 1 (t)]

p = 1 p1 = p + 1

+

(A0 )2 8

6

Nbs ,Nq ,Nf Nbs ,Nq ,Nf

s

s

Absj Aqk Absj1 Aqk1

j,k,l = 1 j 1 ,k 1 ,l 1 = 1

4

4

4

4

p=1

p=1

p=1

p=1

× s Bp (t) s B'p (t) + s B p (t) s B p' (t)

+

A0 n(t) 2

Nbs ,Nq ,Nf

s j,k,l = 1

+ n 2(t) + n˜ 2(t),

4

Absj Aqk s Bp (t) + p=1

A0 n˜ (t) 2

7

Nbs ,Nq ,Nf

4

s

Absj Aqk s B p (t)

j,k,l = 1

p=1

(except j = j1 + k = k1 + l = l1 )

(A2)

where B'p (t) = Bpj1 ,k1 ,l1 (t),

(p = 1, 2, 3, 4)

(A3a)

B p' (t) = B pj1 ,k1 ,l1 (t).

(p = 1, 2, 3, 4)

(A3b)

and

.-.   . . 

426

In the presence of a bearing fault, n(t) is very small compared with the signal xbs (t)xq (t)xf (t), and n 2(t) + n˜ 2(t) is rather smaller and can be ignored. Expanding and simplifying equation (A2), one has =z(t)=2 =

(A0 )2 8

Nbs ,Nq ,Nf

(Absj )2(Aqk )2{2 + 2 cos [2p(2lff )t] + 2 cos [2p(2kfq )t]

s j,k,l = 1

+ cos [2p(2lff + 2kfq )t] + cos [2p(2lff − 2kfq )t]} +

+

(A0 )2 8

Nbs ,Nq ,Nf Nbs ,Nq ,Nf

s

s

15

Absj Aqk Absj1 Aqk1 s H'm 2 (t)

j,k,l = 1 j 1 ,k 1 ,l 1 = 1

A0 n(t) 2

m2 = 0

Nbs ,Nq ,Nf

4

s

Absj Aqk s Bp (t) +

j,k,l = 1

p=1

A0 n˜ (t) 2

Nbs ,Nq ,Nf

4

s

Absj Aqk s B p (t),

j,k,l = 1

p=1

(except j = j1 + k = k1 + l = l1 )

(A4)

where H'0,1,...,15 (t) = cos {2p[(fbsj − fbsj1 ) 2 (l 2 l1 )ff 2 (k 2 k1 )fq ]t}.

(A5)

If j $ j1 , with inequality (3), the frequencies [(fbsj − fbsj1 ) 2 (l 2 l1 )ff 2 (k 2 k1 )fq ](l, l1 = 1, 2, . . . , Nf ; k, k1 = 1, 2, . . . , Nq ) are in high-frequency, so these components can be ignored in the very low shaft speed case. In the second term of equation (A4), only the components with j = j1 are remained, and one denotes A0 n(t) 2

nb(t),

A 0 nb (t),

Nbs ,Nq ,Nf

4

s

Absj Aqk s Bp (t),

j,k,l = 1

p=1

Nbs ,Nq ,Nf 0

2

n˜ (t)

4

s

Absj Aqk s B p (t).

j,k,l = 1

p=1

(A6)

If the autocorrelation is calculated for the envelope signal as shown in equation (8), because of 2 E[nb(t)nb(t + t)] = snb , (A7) 2 E[0 nb (t)0 nb (t + t)] = snb ,

(A8)

is the variance of the signals nb(t) and 0 nb (t), with equation (11b), the function

where s Fzz (t) can be expressed as 2 nb

Fzz (t) = C +

1 2

Nbs ,Nq ,Nf

(Dj,k )2{2 cos [2p(2lff )t] + 2 cos [2p(2kfq )t]

s j,k,l = 1

+ cos [2p(2lff + 2kfq )t] + cos [2p(2lff − 2kfq )t]} Nq

+ s

Nf

7

s (Ek,k 1 )2 s Hm 1 (t),

k,k 1 = 1 l,l 1 = 1

(except k = k1 + l = l1 )

(A9)

m1 = 0

where H0,1,...,7 (t) = cos {2p[(l 2 l1 )ff 2 (k 2 k1 )fq ]t}, and Fzz (t) is a function which is related to Rzz (t).

(A10)