- Email: [email protected]
Measurement of the dispersion relation of guided non-axisymmetric waves in filament-wound cylindrical structures
Ultrasonics 38 (2000) 517–521 www.elsevier.nl/locate/ultras
Measurement of the dispersion relation of guided non-axisymmetric waves in filament-wound cylindrical structures D. Gsell*, D. Profunser, J. Dual Institute of Mechanics, ETH Zurich, Swiss Federal Institute of Technology, CH-8092 Zurich, Switzerland
Abstract To determine the dispersion relation, guided waves are excited in specimens over a broad frequency range. The surface displacements are measured over time and space. The recorded data are analysed using a quasi-three-dimensional spectrum estimation algorithm. In the time domain a fast Fourier transform is used to extract the frequencies. To obtain the wave numbers, in space a two-dimensional matrix-pencil approach is applied to the data set. Using a suitable constitutive model (transversely isotropic or orthotropic) dispersion curves are calculated. Good agreement was found between the experimental and the numerically calculated dispersion relations after adjusting the material parameters. Since the dispersion relation of a structure depends on the mechanical material properties frequency-dependent material parameters can be extracted from the above-mentioned relation between frequency and wave number. © 2000 Elsevier Science B.V. All rights reserved. Keywords: Cylindrical structures; Dispersion relation; Filament wound; Matrix-pencil; Waves
1. Introduction Frequency-dependent material properties are used by design engineers to predict the dynamic behaviour of structures. One possible approach to obtain the material parameters can be outlined as follows: 1. a wave propagation phenomena is chosen; 2. measurements are performed; 3. acquired data are processed in order to obtain characteristic features of the structure–wave interaction such as dispersion relations; 4. dispersion relations are calculated by a constitutive simulation model; 5. inversion: the parameters of the simulation model are fitted to the experimental data. In this paper a method will be outlined, which enables the measurement of the dispersion relation from nonaxisymmetric wave propagation in cylindrical structures.
2. Wave propagation in orthotropic cylindrical shells In this section a brief description of the calculation of the dispersion relation of non-axisymmetric harmonic * Corresponding author. Tel.: +41-1-632-77-50; fax: +41-1-632-11-45. E-mail address: [email protected] (D. Gsell )
waves in orthotropic cylindrical shells will be made. In order to determine all nine elastic constants, non-axisymmetric waves need to be analysed. To calculate the frequency–wave number relation, one considers a differential element of a shell (Fig. 1). The equations are formulated with respect to a cylindrical coordinate system (r, Q, x), which coincides with the principal orthotropy directions of the material. The constitutive law for an orthotropic material: s=C · e
(1)
(where C is the stiffness matrix with nine independent material parameters), the equilibrium conditions and the kinematic relations in cylindrical coordinates result in a set of three coupled partial differential equations for the displacements u , u and u . According to Shul’ga r Q x and Ramsakaya [1], this system is solved by functions for the displacements of the form: u (r, Q, x)=f (r) exp[i(k x+k RQ−vt)], j j x Q
j=r, Q, x, (2)
where k is the axial wave number, k the circumferential x Q wave number and v is the angular frequency. Typically one would expect periodicity in the direction. In addition, the six boundary conditions at the free surfaces
0041-624X/00/$ - see front matter © 2000 Elsevier Science B.V. All rights reserved. PII: S0 0 4 1 -6 2 4 X ( 9 9 ) 0 0 06 5 - 7
518
D. Gsell et al. / Ultrasonics 38 (2000) 517–521
frequency range from 5 to 500 kHz. These discrete time signals are stored in a computer to perform further data analysis.
4. Signal processing Fig. 3 illustrates the digital data processing, leading from the displacement measurement to the dispersion relation. In the time domain a fast Fourier transform is applied to each measurement point in order to obtain the complex frequency spectrum. From these spectra, measured for various values of x and Q, the parameters k and k need to be determined, according to: xm Qm Fig. 1. Differential shell element and stresses involved in non-axisymmetric wave propagation.
M u(k, p)= ∑ bˆ exp[i(k x+k RQ)]. m xm Qm m=1
(5)
must be fulfilled: s (r=R±H, Q, x, t)=0, j=r, Q, x, (3) rj where R+H is the outer radius and R−H the inner radius of the cylindrical shell. From these equations the dispersion relation v=f(k, p) can be calculated. Therefore the measured surface displacements can be mathematically described by the sum of harmonic functions: M u(k , k , v)= ∑ b exp[i(k x+k RQ−vt)], x Q m xm Qm m=1 where M is the unknown number of wave modes.
(4)
3. Experimental set-up Waves are excited with a piezoelectric transducer (Fig. 2), acting radially, over a broad frequency range. At several points, arranged circumferentially and axially, the surface displacements versus time are measured with a heterodyne laser interferometer over a wide
4.1. Estimating one-dimensional wave numbers In the past few years several one- or two-dimensional algorithms have been developed to estimate parameters of superimposed exponential signals, such as ‘linear prediction’, ‘maximum likelihood’, … and a ‘matrixpencil’ approach. The matrix-pencil algorithm developed by Hua and Sarkar [2] is chosen because it is a very fast and precise tool for estimating the unknown signal parameters. A brief outline is given below for the onedimensional case. A data vector x containing the complex amplitudes (received by the FFT ) at one particular frequency in the space domain is built: x=[x , x , ,, x ]. 0 1 N−1
The complex amplitude of the signal at one particular location can be described as: M M x = ∑ bˆ zn +w = ∑ bˆ exp(ik Dx n)+w , (7) n m m n m xm n m=1 m=1 where bˆ are the complex amplitudes, z the signal m m poles, k the wave numbers and Dx the increment xm between two measured points in the x-direction. w n denotes noise at the chosen frequency. From the data vector x two signal matrices are built:
X = 0 Fig. 2. Experimental set-up.
(6)
C
x
L−1 x L e
x
L−2
,
x
x
L−1 e
,
x
x N−1
x
N−2
,
0
1 e
x
N−L−1
D
,
(8)
519
D. Gsell et al. / Ultrasonics 38 (2000) 517–521
Fig. 3. Visualization of the major steps in signal processing.
and
part (damping factor) can be calculated from:
C
X = 1
x x
L
L−1 e
x N
x L−1 x L−2 e
,
x
,
x
x
,
N−1
1
2 e
x N−L
D
Re(k )= m ,
(9)
where L is called the pencil parameter. L is typically chosen as L=N/2. With the property of the underlying exponential functions, exp(ik Dx n)=exp[ik Dx (n−1)] exp(ik Dx), (10) m m m one can show, that for each of {z=z ; m=1, ,, M} m the rank of the matrix-pencil X −z · X is reduced by 1 0 one. In other words, the signal poles z are the solutions m of the generalized eigenvalue problem: (X −zX ) · q=0, (11) 1 0 or after multiplying Eq. (11) by the pseudo-inverse X# , of the eigenvalue problem: 0 (X# · X −zI ) · q=0, (12) 0 1 where the superscript # denotes the pseudo-inverse and I denotes the identity matrix. In the noiseless case, there exist exactly M non-zero eigenvalues corresponding to the M wave modes, because the matrix-pencil has rank M. For a noisy data set x the matrix-pencil has full rank. In a first step one has to estimate the number M of superimposed signals in the data matrices. This can be done by some statistically based model selection criterion, such as the AIC criterion of Akaike or the MDL criterion of Schwartz [3]. Afterwards in Eq. (12) the pseudo-inverse X# is 0 replaced (or approximated) by the rank-M-truncated ˜ pseudo-inverse X# . This low-rank matrix approximation 0 can be obtained by singular value decomposition. The real part of the wave number and the imaginary
angle(z ) m , Dx
Im(k )=− m
ln(|z |) m . Dx
(13)
An example with experimental data is given in (Fig. 4). A carbon-fibre-reinforced shell with a Epoxy matrix is used as specimen, with inner radius r =14.01 mm and i outer radius r =15.60 mm.The angle of winding was a alternately ±67.5° with respect to the axis of the tube. It consists of 15 layers. The shell is examined in a frequency range 50–400 kHz. The full curves are the theoretical ones, as given by the theory of Shul’ga, using material parameters according to Dual [4]. The gray dots are estimated from the experimental data. Very good agreement between theory and experiment was found. The measurements outside the modes are attributed to noise or other wave modes.
Fig. 4. Dispersion relations of a composite tube: experiment versus theory. The axisymmetric wave modes are shown in the figure.
520
D. Gsell et al. / Ultrasonics 38 (2000) 517–521
4.2. Quasi two-dimensional signal processing As mentioned above, there exist several truly twodimensional algorithms which estimate the parameters of a two-dimensional wave field. In our case, the extraction of wave numbers in a tube, the amount of CPU time becomes very difficult to handle, since at one particular frequency low and high wave numbers exist simultaneously. High wave numbers require a high spatial sampling rate (Nyquist: k ∏p/Dx) and low max ones need sampling over a large distance, especially when noise is present. Therefore many measurement points are necessary, yielding large matrices in Eq. (11). This problem is avoided by the following new approach. The two-dimensional problem is analysed by a combination of four one-dimensional measurements. In a first step measurements are performed along the tube axis x and in the circumferential direction Q ( Fig. 5). Based on this the wave numbers k and k are estimated ix iQ independently. If more than one wave number is found in each direction, one needs to pair the extracted wave numbers correctly. The following pairing algorithm is proposed. $ Perform further measurements along two coordinate axes j and f. With respect to the coordinate system (Q, x), the system (j, f) is rotated by an angle a ( Fig. 6). $ Estimate the wave numbers k and k independently. j f $ Choose a wave vector {k , k } and project it onto Q x the axes (j, f). The resulting wave numbers will be called kˆ and kˆ : j f kˆ cos a sin a k j = · Q . (14) −sin a cos a kˆ k f x If the chosen pair {k , k } is a really an existing pair, Q x the projections kˆ and kˆ must exist as estimated wave j f numbers in the rotated coordinate system (j, f). $ Repeat the steps above, until every combination of the pairs (k , k ) has been tested. Q x
CDC
DC D
Fig. 5. Measurement axes of the quasi-two-dimensional algorithm.
Fig. 6. Example with two wave vectors k and k . The wave vectors 1 2 are projected onto the four axes, yielding the one-dimensional components.
A numerical experiment is given below ( Fig. 7). For one particular frequency a two-dimensional wave field is generated, containing eleven wavemodes. The chosen wave numbers (marked as circles) simulate a composite shell. The field is sampled at 240 points along each of the four axes ( Fig. 6), with increments R DQ=0.79 mm, Dx=1.5 mm, Dj=0.87 mm and Df=1.16 mm, where the outer radius of the shell is set to R=15.0 mm. The rotation angle a between the two coordinate systems is set to 30°. Gaussian noise is added to the simulated field. The signal to noise ratio amounts to 5 dB.
Fig. 7. Numerical experiment with eleven wavemodes, analysed by the quasi two-dimensional matrix-pencil algorithm.
D. Gsell et al. / Ultrasonics 38 (2000) 517–521
The algorithm proposed above estimates wave number pairs (black crosses), with a precision. Compared with a truly two-dimensional two major advantages shall be mentioned: $ fewer data are necessary to obtain a good $ the amount of CPU time is reduced by magnitude.
all eleven very high
521
dimensional case (waves in plates, earthquake waves, …). The proposed algorithm can be easily extended to the three-dimensional DOA problem.
algorithm result; orders of
5. Conclusions A new algorithm to estimate signal parameters in a two-dimensional wave field has been developed. The algorithm is specially designed to estimate signal parameters from the travelling wavemodes in circular cylindrical shells. The proposed method can also be applied to the estimation of the direction of arrival (DOA) in the two-
References [1] N.A. Shul’ga, E.I. Ramsakaya, Propagation of nonaxisymmetrical elastic waves in an orthotropic hollow cylinder, Sov. Appl. Mech. 19 (9) (1983) 748–752. [2] Y. Hua, T.K. Sarkar, Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise, IEEE Trans. Signal Process. 38 (5) (1990) 814–824. [3] V.U. Reddy, L.S. Biradar, SVD-based information theoretic criteria for detection of the number of damped/undamped sinusoids and their performance analyses, IEEE Trans. Signal Process. 41 (9) (1993) 2872–2881. [4] J. Dual, Quantitative nondestructive evaluation using guided waves, in: T. Kishi, T. Saito, C. Ruud, R. Green ( Eds.), Nondestructive Characterization of Materials, Proceedings 5th International Symposium on Nondestructive Characterization of Materials, Gordon and Breach, Montreux, 1991, p. 1061.
Measurement of the dispersion relation of guided non-axisymmetric waves in filament-wound cylindrical structures D. Gsell*, D. Profunser, J. Dual Institute of Mechanics, ETH Zurich, Swiss Federal Institute of Technology, CH-8092 Zurich, Switzerland
Abstract To determine the dispersion relation, guided waves are excited in specimens over a broad frequency range. The surface displacements are measured over time and space. The recorded data are analysed using a quasi-three-dimensional spectrum estimation algorithm. In the time domain a fast Fourier transform is used to extract the frequencies. To obtain the wave numbers, in space a two-dimensional matrix-pencil approach is applied to the data set. Using a suitable constitutive model (transversely isotropic or orthotropic) dispersion curves are calculated. Good agreement was found between the experimental and the numerically calculated dispersion relations after adjusting the material parameters. Since the dispersion relation of a structure depends on the mechanical material properties frequency-dependent material parameters can be extracted from the above-mentioned relation between frequency and wave number. © 2000 Elsevier Science B.V. All rights reserved. Keywords: Cylindrical structures; Dispersion relation; Filament wound; Matrix-pencil; Waves
1. Introduction Frequency-dependent material properties are used by design engineers to predict the dynamic behaviour of structures. One possible approach to obtain the material parameters can be outlined as follows: 1. a wave propagation phenomena is chosen; 2. measurements are performed; 3. acquired data are processed in order to obtain characteristic features of the structure–wave interaction such as dispersion relations; 4. dispersion relations are calculated by a constitutive simulation model; 5. inversion: the parameters of the simulation model are fitted to the experimental data. In this paper a method will be outlined, which enables the measurement of the dispersion relation from nonaxisymmetric wave propagation in cylindrical structures.
2. Wave propagation in orthotropic cylindrical shells In this section a brief description of the calculation of the dispersion relation of non-axisymmetric harmonic * Corresponding author. Tel.: +41-1-632-77-50; fax: +41-1-632-11-45. E-mail address: [email protected] (D. Gsell )
waves in orthotropic cylindrical shells will be made. In order to determine all nine elastic constants, non-axisymmetric waves need to be analysed. To calculate the frequency–wave number relation, one considers a differential element of a shell (Fig. 1). The equations are formulated with respect to a cylindrical coordinate system (r, Q, x), which coincides with the principal orthotropy directions of the material. The constitutive law for an orthotropic material: s=C · e
(1)
(where C is the stiffness matrix with nine independent material parameters), the equilibrium conditions and the kinematic relations in cylindrical coordinates result in a set of three coupled partial differential equations for the displacements u , u and u . According to Shul’ga r Q x and Ramsakaya [1], this system is solved by functions for the displacements of the form: u (r, Q, x)=f (r) exp[i(k x+k RQ−vt)], j j x Q
j=r, Q, x, (2)
where k is the axial wave number, k the circumferential x Q wave number and v is the angular frequency. Typically one would expect periodicity in the direction. In addition, the six boundary conditions at the free surfaces
0041-624X/00/$ - see front matter © 2000 Elsevier Science B.V. All rights reserved. PII: S0 0 4 1 -6 2 4 X ( 9 9 ) 0 0 06 5 - 7
518
D. Gsell et al. / Ultrasonics 38 (2000) 517–521
frequency range from 5 to 500 kHz. These discrete time signals are stored in a computer to perform further data analysis.
4. Signal processing Fig. 3 illustrates the digital data processing, leading from the displacement measurement to the dispersion relation. In the time domain a fast Fourier transform is applied to each measurement point in order to obtain the complex frequency spectrum. From these spectra, measured for various values of x and Q, the parameters k and k need to be determined, according to: xm Qm Fig. 1. Differential shell element and stresses involved in non-axisymmetric wave propagation.
M u(k, p)= ∑ bˆ exp[i(k x+k RQ)]. m xm Qm m=1
(5)
must be fulfilled: s (r=R±H, Q, x, t)=0, j=r, Q, x, (3) rj where R+H is the outer radius and R−H the inner radius of the cylindrical shell. From these equations the dispersion relation v=f(k, p) can be calculated. Therefore the measured surface displacements can be mathematically described by the sum of harmonic functions: M u(k , k , v)= ∑ b exp[i(k x+k RQ−vt)], x Q m xm Qm m=1 where M is the unknown number of wave modes.
(4)
3. Experimental set-up Waves are excited with a piezoelectric transducer (Fig. 2), acting radially, over a broad frequency range. At several points, arranged circumferentially and axially, the surface displacements versus time are measured with a heterodyne laser interferometer over a wide
4.1. Estimating one-dimensional wave numbers In the past few years several one- or two-dimensional algorithms have been developed to estimate parameters of superimposed exponential signals, such as ‘linear prediction’, ‘maximum likelihood’, … and a ‘matrixpencil’ approach. The matrix-pencil algorithm developed by Hua and Sarkar [2] is chosen because it is a very fast and precise tool for estimating the unknown signal parameters. A brief outline is given below for the onedimensional case. A data vector x containing the complex amplitudes (received by the FFT ) at one particular frequency in the space domain is built: x=[x , x , ,, x ]. 0 1 N−1
The complex amplitude of the signal at one particular location can be described as: M M x = ∑ bˆ zn +w = ∑ bˆ exp(ik Dx n)+w , (7) n m m n m xm n m=1 m=1 where bˆ are the complex amplitudes, z the signal m m poles, k the wave numbers and Dx the increment xm between two measured points in the x-direction. w n denotes noise at the chosen frequency. From the data vector x two signal matrices are built:
X = 0 Fig. 2. Experimental set-up.
(6)
C
x
L−1 x L e
x
L−2
,
x
x
L−1 e
,
x
x N−1
x
N−2
,
0
1 e
x
N−L−1
D
,
(8)
519
D. Gsell et al. / Ultrasonics 38 (2000) 517–521
Fig. 3. Visualization of the major steps in signal processing.
and
part (damping factor) can be calculated from:
C
X = 1
x x
L
L−1 e
x N
x L−1 x L−2 e
,
x
,
x
x
,
N−1
1
2 e
x N−L
D
Re(k )= m ,
(9)
where L is called the pencil parameter. L is typically chosen as L=N/2. With the property of the underlying exponential functions, exp(ik Dx n)=exp[ik Dx (n−1)] exp(ik Dx), (10) m m m one can show, that for each of {z=z ; m=1, ,, M} m the rank of the matrix-pencil X −z · X is reduced by 1 0 one. In other words, the signal poles z are the solutions m of the generalized eigenvalue problem: (X −zX ) · q=0, (11) 1 0 or after multiplying Eq. (11) by the pseudo-inverse X# , of the eigenvalue problem: 0 (X# · X −zI ) · q=0, (12) 0 1 where the superscript # denotes the pseudo-inverse and I denotes the identity matrix. In the noiseless case, there exist exactly M non-zero eigenvalues corresponding to the M wave modes, because the matrix-pencil has rank M. For a noisy data set x the matrix-pencil has full rank. In a first step one has to estimate the number M of superimposed signals in the data matrices. This can be done by some statistically based model selection criterion, such as the AIC criterion of Akaike or the MDL criterion of Schwartz [3]. Afterwards in Eq. (12) the pseudo-inverse X# is 0 replaced (or approximated) by the rank-M-truncated ˜ pseudo-inverse X# . This low-rank matrix approximation 0 can be obtained by singular value decomposition. The real part of the wave number and the imaginary
angle(z ) m , Dx
Im(k )=− m
ln(|z |) m . Dx
(13)
An example with experimental data is given in (Fig. 4). A carbon-fibre-reinforced shell with a Epoxy matrix is used as specimen, with inner radius r =14.01 mm and i outer radius r =15.60 mm.The angle of winding was a alternately ±67.5° with respect to the axis of the tube. It consists of 15 layers. The shell is examined in a frequency range 50–400 kHz. The full curves are the theoretical ones, as given by the theory of Shul’ga, using material parameters according to Dual [4]. The gray dots are estimated from the experimental data. Very good agreement between theory and experiment was found. The measurements outside the modes are attributed to noise or other wave modes.
Fig. 4. Dispersion relations of a composite tube: experiment versus theory. The axisymmetric wave modes are shown in the figure.
520
D. Gsell et al. / Ultrasonics 38 (2000) 517–521
4.2. Quasi two-dimensional signal processing As mentioned above, there exist several truly twodimensional algorithms which estimate the parameters of a two-dimensional wave field. In our case, the extraction of wave numbers in a tube, the amount of CPU time becomes very difficult to handle, since at one particular frequency low and high wave numbers exist simultaneously. High wave numbers require a high spatial sampling rate (Nyquist: k ∏p/Dx) and low max ones need sampling over a large distance, especially when noise is present. Therefore many measurement points are necessary, yielding large matrices in Eq. (11). This problem is avoided by the following new approach. The two-dimensional problem is analysed by a combination of four one-dimensional measurements. In a first step measurements are performed along the tube axis x and in the circumferential direction Q ( Fig. 5). Based on this the wave numbers k and k are estimated ix iQ independently. If more than one wave number is found in each direction, one needs to pair the extracted wave numbers correctly. The following pairing algorithm is proposed. $ Perform further measurements along two coordinate axes j and f. With respect to the coordinate system (Q, x), the system (j, f) is rotated by an angle a ( Fig. 6). $ Estimate the wave numbers k and k independently. j f $ Choose a wave vector {k , k } and project it onto Q x the axes (j, f). The resulting wave numbers will be called kˆ and kˆ : j f kˆ cos a sin a k j = · Q . (14) −sin a cos a kˆ k f x If the chosen pair {k , k } is a really an existing pair, Q x the projections kˆ and kˆ must exist as estimated wave j f numbers in the rotated coordinate system (j, f). $ Repeat the steps above, until every combination of the pairs (k , k ) has been tested. Q x
CDC
DC D
Fig. 5. Measurement axes of the quasi-two-dimensional algorithm.
Fig. 6. Example with two wave vectors k and k . The wave vectors 1 2 are projected onto the four axes, yielding the one-dimensional components.
A numerical experiment is given below ( Fig. 7). For one particular frequency a two-dimensional wave field is generated, containing eleven wavemodes. The chosen wave numbers (marked as circles) simulate a composite shell. The field is sampled at 240 points along each of the four axes ( Fig. 6), with increments R DQ=0.79 mm, Dx=1.5 mm, Dj=0.87 mm and Df=1.16 mm, where the outer radius of the shell is set to R=15.0 mm. The rotation angle a between the two coordinate systems is set to 30°. Gaussian noise is added to the simulated field. The signal to noise ratio amounts to 5 dB.
Fig. 7. Numerical experiment with eleven wavemodes, analysed by the quasi two-dimensional matrix-pencil algorithm.
D. Gsell et al. / Ultrasonics 38 (2000) 517–521
The algorithm proposed above estimates wave number pairs (black crosses), with a precision. Compared with a truly two-dimensional two major advantages shall be mentioned: $ fewer data are necessary to obtain a good $ the amount of CPU time is reduced by magnitude.
all eleven very high
521
dimensional case (waves in plates, earthquake waves, …). The proposed algorithm can be easily extended to the three-dimensional DOA problem.
algorithm result; orders of
5. Conclusions A new algorithm to estimate signal parameters in a two-dimensional wave field has been developed. The algorithm is specially designed to estimate signal parameters from the travelling wavemodes in circular cylindrical shells. The proposed method can also be applied to the estimation of the direction of arrival (DOA) in the two-
References [1] N.A. Shul’ga, E.I. Ramsakaya, Propagation of nonaxisymmetrical elastic waves in an orthotropic hollow cylinder, Sov. Appl. Mech. 19 (9) (1983) 748–752. [2] Y. Hua, T.K. Sarkar, Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise, IEEE Trans. Signal Process. 38 (5) (1990) 814–824. [3] V.U. Reddy, L.S. Biradar, SVD-based information theoretic criteria for detection of the number of damped/undamped sinusoids and their performance analyses, IEEE Trans. Signal Process. 41 (9) (1993) 2872–2881. [4] J. Dual, Quantitative nondestructive evaluation using guided waves, in: T. Kishi, T. Saito, C. Ruud, R. Green ( Eds.), Nondestructive Characterization of Materials, Proceedings 5th International Symposium on Nondestructive Characterization of Materials, Gordon and Breach, Montreux, 1991, p. 1061.