Calculation of Guided Wave Dispersion Characteristics Using a Three-Transducer Measurement System

Guided ultrasonic waves are of significant interest in the health monitoring of thin structures, and dispersion curves are important tools in the deployment of any guided wave application. Most methods of determining dispersion curves require accurate knowledge of the material properties and thickness of the structure to be inspected, or extensive experimental tests. This paper presents an experimental technique that allows rapid generation of dispersion curves for guided wave applications when knowledge of the material properties and thickness of the structure to be inspected are unknown. The technique uses a single source and measurements at two points, making it experimentally simple. A formulation is presented that allows calculation of phase and group velocities if the wavepacket propagation time and relative phase shift can be measured. The methodology for determining the wavepacket propagation time and relative phase shift from the acquired signals is described. The technique is validated using synthesized signals, finite element model-generated signals and experimental signals from a 3 mm-thick aluminium plate. Accuracies to within 1% are achieved in the experimental measurements.


Introduction
Guided waves have been extensively investigated in the recent years as a nondestructive testing (NDT) technique because of their advantages compared to ultrasonic testing (UT) inspection. The key benefit of this technology is the ability to interrogate the entire thickness of thin-walled structures over large areas from a single location.
During guided wave propagation, waves experience dispersion that distorts the wave shape as the wave propagates, due to a dependence of velocity on frequency [1,2]. Dispersion curves show the relation of phase and group velocities against frequency, for a particular geometry and material. Determination of dispersion curves is important for any guided wave application; accurate dispersion curves enable wave modes of received wavepackets to be distinguished, specific wave modes to be cancelled to clean the acquired signal, the propagation of wave modes in a particular direction using phased arrays [3] and the spatial location of damage in the structure to be determined based on time-of-flight measurements [4]. Guided waves are being used in commercial products, for instance to detect wall thickness losses due to corrosion for pipeline inspection in the oil and gas industry. These commercial products generate a unique wave mode which propagates along the structure avoiding the creation of multimodal excitations, as this would increase the complexity of the analysis of the signals [5,6]. Problems with the application of guided waves occur when the material properties or thickness of the structure to be inspected are unknown due to imprecise records, commercial confidentiality or manufacturing and maintenance uncertainty. This is a particular problem when evaluating composite structures, like wind turbine blades, where elastic constants are not often available. To inspect this kind of structure in an industrial environment, it becomes impracticable as currently there are no means of generating the dispersion curves for such situations. Therefore, there is a need to be able to calculate in situ the phase and group velocities at frequencies of interest in a quick and reliable way without requiring prior knowledge of material properties or thickness.
There has been much work published on the determination of dispersion curves describing many analytical, numerical and experimental methods. For relatively simple structures like plates or pipes, dispersion curves can be predicted using the commercially available software Disperse ® , which uses the global matrix method [7]. Other methods, such as the transfer matrix method [8,9] or pseudospectral collocation method [10], have been utilized to determine dispersion curves. Finite element (FE) [11] and semi-analytical finite element [12,13] methods have been also used. All these numerical and analytical methods require knowledge of material properties and thickness of the structure to determine the dispersion curves.
Experimental methods do not suffer this limitation, since guided wave data is directly acquired from the inspected structure. The most widely used experimental technique to measure dispersion curves is the 2D fast Fourier transform (FFT) [14]. This signal processing technique requires the acquisition of many signals along the wave propagation direction in order to carry out a double FFT in time and space. The result is a wavenumber-frequency matrix, which can be rearranged to give the phase velocity dispersion map. The acquisition of the signals from all required locations is manually prohibitive; therefore, a laser scanning vibrometer (LSV) is commonly used to automatically obtain the signals from preselected points. This device is highly sensitive and bulky, being restricted to controlled areas like laboratories [15] and difficult to use in situ.
Other experimental techniques have been proposed recently. Harb and Yuan [16,17] presented a noncontact technique using an air-coupled transducer (ACT) to generate the wave mode on the plate and an LSV to acquire the propagating mode. The technique requires precise control of the incidence angle of the ultrasonic pressure from the ACT upon the surface. By relating the frequency to the incidence angle at which the wave amplitude is maximized, it is possible to calculate the phase velocity using Snell's law. However, the technique has limited applicability outside of the laboratory. In works by Mažeika et al. [18,19], a zero-crossing technique is used to calculate the phase velocity based on the measurement of the position of constant phase points of the pulse corresponding to zero amplitude as a function of time. A similar approach is taken in the phase velocity method in [20], tracking the peaks of the pulse rather than the zero crossings. To produce convincing results, both techniques require the acquisition of a large number of signals at different distances, making these techniques time-consuming. For the sake of reducing time and sampling points in space, Adams et al. [21] presented a method using an array probe from which a time-space matrix can be created. A 2D filter is applied to the matrix to extract the phase velocities. The technique was demonstrated to be valid in simulation using FE analysis. However, for realistic probe sizes, the technique has limited experimental applications. In [22], sparse wavenumber analysis was used to experimentally recover the dispersion curves from an aluminum plate; however, seventeen sensors were needed to deploy the technique. An alternative approach to getting dispersion curves from experimental signals is to use time-frequency representations [23], where group velocity can be directly determined. However, the calculation of phase velocity involves the integration of group velocity requiring precise values of ω and k at the lower limit, which are not easy to obtain experimentally. For medical purposes, two techniques, SVD and LRT, were evaluated for extracting the dispersion curves of cortical bones in [24]. Both techniques require a large number of receivers and transmitters.
Experimental calculation of bulk wave velocities in dispersive materials has been studied [25]. Sachse and Pao used the method of phase spectral analysis, where the phase delay is calculated using the real and imaginary parts of the Fourier transform of the pulse, and then obtaining the group delay through the differentiation of the phase characteristics of the pulse with respect to frequency. Pialucha et al. presented the amplitude spectrum technique overcoming the requirement of the phase spectrum technique for the separation of successive pulses in the time domain for bulk waves [26]. However, the phase spectral technique provided better results when pulses can be separated.
In this paper, an experimental technique for straightforward calculations of phase and group velocities of guided waves at frequencies of interest is presented. The technique requires just the acquisition of two signals spaced a few cm using conventional transducers, enabling its application on in-service structures located in poorly controlled environments. There is no requirement for prior knowledge of the material properties or thickness of the inspected structure. In Section 2, the formulation on which the experimental technique relies is presented. Section 3 describes in detail the methodology to extract the group delay and phase shift from experimental signals. Successive sections present validation tests using synthesized signals and simulated signals from FE analysis. Finally, experimental signals from a 3 mm-thick aluminium plate are used to validate the technique as a solution to the problem presented in the introduction.

Theoretical Basis
Assuming that the propagating pulse is sufficiently narrow-band that dispersion does not significantly distort the wavepacket over the propagation distance, the signal can be approximated as [27], where x is propagation distance, t is time, ω = 2π f is the angular frequency of the harmonic wave at frequency f and k = 2π/λ is the wavenumber. G(t) is a function defining the temporal pulse shape. G(t) is defined such that the pulse has a peak at t = 0. v g is the group velocity at which the packet propagates, and so if the trajectory of the wavepacket peak is defined as the point whereby The wavepacket is shown schematically in Figure 1. The harmonic part of the signal has a phase of and propagates at the phase velocity [6,27].
The phase of the harmonic signal at the centre of the wavepacket x 0 , t 0 is Rearranging and substituting for v p and v g gives For nondispersive waves, v p = v g and Equation (7) gives the expected result that the phase at the wavepacket center is φ 0 = 0 for all propagation times. Where the system is dispersive there is a finite phase difference between the harmonic part of the signal and the wavepacket, φ 0 , which increases linearly with propagation time and is positive if v g > v p , negative if v g < v p . Assuming that this phase difference can be measured experimentally, Equations (4) and (5) further yield two useful equations for evaluating the phase velocity: In Figure 2, an illustrative example of a S 0 -mode wavepacket at 150 kHz is shown. The solid line is the S 0 mode after 1.2-m propagation including the dispersion. The dashed line represents the S 0 mode after the propagation without dispersion, namely, the wavepacket was calculated using the same phase velocity for all frequencies. In this example, the phase shift is negative (-1.26 radians), which means the phase velocity is higher than the group velocity. Using (9), the calculation of the phase velocity would be 1.2 m divided by the sum of 223 µs, which is the time of the wavepacket to travel 1.2 m at group velocity, plus -1.34 µs which is the time that the signal at an angular frequency, ω = 2π·150000, takes to cover -1.26 radians. Group velocity and phase velocity are calculated using (2) and (9), respectively: v g = 5381 ms −1 and v p = 5414 ms −1 .

Methodology
In this section, the methodology for determining the phase shift and hence the phase velocity is presented. This general methodology is complemented by a description of a specific experiment setup described in Section 6.1. The technique is based on the analysis of two signals acquired at different distances in the direction of wave propagation. One of them will be treated as the baseline and the other signal will be modified in time and phase in order to match with the baseline. The methodology is composed of different signal processing techniques which are presented in a block diagram in Figure 3.
Methodology's block diagram to extract the optimum phase shift and time delay. 1 2 as can be seen in Figure 6. Afterwards, the Hilbert Transform (H(…)) [2 The experimental setup is composed of one transmitter and two receivers along the propagation direction of the wave to be measured, as shown in Figure 4. For isotropic systems the wave velocity is independent of direction. The technique analyses the effect of propagation between the first receiver and the second receiver. Spacing between receivers of a few centimeters is selected in order to avoid large phase shifts and waveform deformations due to dispersion.
Methodology's block diagram to extract the optimum phase shift and time delay. 1 2 as can be seen in Figure 6. Afterwards, the Hilbert Transform (H(…)) [2 The first step in analyzing the data is the isolation of the desired wavepacket by removing the rest of the signal. The boundaries of each wavepacket are determined by studying the slope of its envelope. The boundaries are set at the closest zero slopes before and after each peak, see example of real signals in Figure 5. Once the region containing the wavepacket is identified, the remainder of the signal is discarded. Both signals are truncated with the same boundaries in order to maintain the time difference between the wavepackets, w 1 and w 2 . The boundaries used are the left limit of the signal of the closest receiver to the transmitter and the second limit of the signal of the further receiver, as can be seen in Figure 6. Afterwards, the Hilbert Transform (H( . . . )) [28] is applied to both wavepackets and the magnitude extracted in order to determine the envelope of each wavepacket: These two envelopes are then cross-correlated to determine the time delay between them:

Transmitter Receivers
where τ is the applied time delay and ∆t the value that maximizes the correlation. Once this time delay is known, the phase delay to the harmonic part of w 2 (t) that maximizes the cross-correlation between the delayed phase-shifted wavepacket w * 2 (t − ∆t) and w 1 (t) is calculated. This process consists of the modification in phase of the signal from the furthest receiver by a phase angle, φ, and finding the phase shift that maximizes the cross correlation. The phase shift, φ, is applied to w * 2 ; then, the cross correlation between w 1 and w * 2 is calculated. The optimum phase shift is the phase angle when the correlation coefficient with time delay, ∆t, is highest.
This methodology enables the calculation of the lag time and phase shift between the two signals acquired at different distances. By knowing these three variables (distance, time and phase) and also the excitation frequency, (2) and (9) can be used to determine the group and phase velocity, respectively.
The proposed methodology in this paper is based on the acquisition of two signals spaced a few centimeters apart, since for such short distances the distortion for any wave mode is relatively small, enabling a straightforward determination of the phase shift.

Synthetic Signal Analysis
The first validation was carried out using synthesized signals. The signal synthesis is based on the adjustment in frequency and wavenumber of the input signal in the frequency domain. Knowing the waveform at one point, in this case the input signal at the transmitting point, and the dispersion characteristics of each wave mode, the wave modes can be reconstructed in the time domain at any distance [29]. If u(x, t) is the reconstructed wave mode at a distance x from the transmitting point, t is the time, F(ω) is the Fourier transform of the input signal and k(ω) is the wavenumber as a function of the angular frequency, then The relation between wavenumber and angular frequency is extracted directly from the dispersion curves provided by Disperse ® and introduced in the integral.
As an initial test of the signal processing method, developed signals are synthesized that are highly pure, without noise and overlapping between modes. The higher the sampling frequency is set, the better accuracy the technique provides. The sampling frequency was set at 10 MHz, since it has been observed that it provides a good balance between computational time and accuracy.
A set of signals for symmetric (S 0 ), antisymmetric (A 0 ) and shear horizontal (SH 0 ) wave modes were created by evaluating equation (14) in MATLAB for two propagation distances, 30 cm and 35 cm: the analysed propagated distance (∆d) in this case is 5 cm. The analysed frequencies are from 60 to 370 kHz with steps of 10 kHz.
The results after applying the signal processing algorithm described in Section 3 to the synthesized signals are presented in Figure 7. The results (black circles) accurately reconstruct the theoretical values from Disperse ® (grey lines). In Figure 8, the extracted phase shift has been represented against the frequency for the three fundamental modes. For the S 0 mode, the phase values have negative values decaying from -2 degrees at 60 kHz to -52 degrees at 370 kHz. For the A 0 mode, conversely, the phase is positive, increasing along the frequency, since the phase velocity is lower than the group velocity; from 377 degrees at 60 kHz to 673 degrees at 370 kHz. In this case, the phase values are much higher than the S 0 and SH 0 ones, because the A 0 mode is highly dispersive at these frequencies. For the SH 0 mode, the extracted phase values are 0 degrees at 60 kHz and 1 degree at 370 kHz. The SH 0 is a nondispersive mode so the phase value remains at near 0 degrees. The sampling frequency is an important factor, since it determines the degree of resolution of the group velocity. As these examples are performed for short propagation distances, the time that the wave modes take to cover that distance is small. Thus, low sampling frequencies are not able to determine the group velocity accurately. Faster modes require higher sampling frequencies to get the same resolution as slower modes. In this case of 5 cm spacing, the propagation time for the S 0 mode at 200 kHz is 9.3 µs; a sampling frequency of 1 MHz will have a time resolution of 1 µs which yields a velocity resolution for this particular example of 556 m/s; with a sampling frequency of 10 MHz, the velocity resolution will be 58.4 m/s; and with a sampling frequency of 100 MHz, the velocity resolution will be 5.8 m/s. This resolution issue can be observed in Figure 7; at lower velocities the curves are smoother, and at higher velocities (S 0 mode) the curves have poorer velocity resolution.

Finite Element Analysis
A three-dimensional FE analysis was performed in Abaqus to simulate signals. While the previous analysis required the input of the dispersion curves to create the propagated signals, this FE analysis does not require prior knowledge of the dispersion characteristics.
FE models were created simulating guided wave propagation at five different frequencies in an aluminium plate (1000 × 1000 × 3 mm). The transmitter and receivers were modelled as ideal point transducers. The transmitting point was placed at the centre of the plate, where the coordinate system was established. A force was applied in the x-direction to simulate a shear transducer. This generated S 0 and A 0 modes along the x-axis and the SH 0 mode along the y-axis. The input signal was a 5-cycle sine with a Hanning window at central frequencies of 60, 80, 100, 120 and 140 kHz. Multiple receiver points were located on the positive x-axis and positive y-axis from the centre of the plate to the edge every 1 cm, resulting in 50 receivers at each axis. The time and spatial discretization for the finite element model was established at 10 -2 µs and 1 mm, respectively, in order to ensure convergence and to obey the Nyquist sample-rate criterion. The mesh was formed by C3D8 elements.
The FE analysis software provides the results decoupled for each axis direction, so it is straightforward to evaluate each wave mode separately. In Figure 9, the results have been computed in cylindrical coordinates to show clearly the three wave modes of propagation at each axis. In Figure 9a, showing radial displacement, S 0 and A 0 modes are depicted, S 0 being the faster mode. In Figure 9b showing tangential displacement, the SH 0 mode propagates along the y-axis, and in Figure 9c, showing out-of-plane displacement, the A 0 mode is clearly present. Since the plate is relatively small, overlaps between wave modes and echoes from edges are produced, as can be seen in Figure 10. This is a particular problem at lower frequencies (60 and 80 kHz) due to the length of the pulses. Therefore, an additional optimisation step has been added to minimize the error introduced by the overlapping. The FE model has 50 receiver points, every 1 cm along the axis so the five least-overlapped signals are selected for analysis. The selection is performed by analyzing the amplitude at the beginning and end of the envelope of the wavepacket of interest; the signals with the lowest amplitude at those regions are selected. Then, the signal processing described in Section 3 is applied for all the pair combinations between the five selected receivers. Once the phase and group velocities are calculated for each combination, the velocities are averaged to give the final velocity estimates. Any outliers that are introduced during selection of the five signals are eliminated. As shown in Figure 11, the phase velocity matches very well with the theoretical values from Disperse ® . In the group velocity graph, the results correlate well with the theoretical velocities, having a noticeable variation in the S 0 value at 60 kHz. This variation occurs at 60 kHz due to the overlap between the S 0 mode and the A 0 mode, which is more pronounced for the longer pulses at this frequency. This error does not appear for the A 0 mode because the signals used for the analysis were out-of-plane displacements, where the S 0 mode is practically imperceptible. Note that the results improve as the frequency increases, since the wavelength decreases and mode separation increases. Overlapping of pulses is an issue for this technique, as it changes the phase and waveform of the analyzed mode, highly distorting the results and impeding its application.  The results for the phase velocity are better than for the group velocity, more notably when overlapping occurs. This is because of the phase-shift term in the denominator in (9). This contribution in the phase velocity equation minimizes the erroneous calculated value of the propagation time and also provides more resolution on the calculation of the phase velocity, removing the stepped shape seen in the group velocity.

Test Setup
The specimen is a 4 m by 2 m aluminium plate of 3 mm thickness. Two piezoelectric transducers were used, one transmitter and one receiver. The transmitter was fixed at the centre of the plate using a tool that provides force over the transducer; and the receiver was manually attached at various distances from the transmitter. For the evaluation of the S 0 mode, the receiver was placed at 45, 50 and 55 cm from the transmitter. For the case of the A 0 and SH 0 modes, the receiver was set at 25, 30 and 35 cm from the transmitter. Those distances were selected to minimize the wavepacket distortion due to overlaps.
A pulser-receiver (Teletest ® unit) was used to excite the transducer and acquire the signals. Two independent channels were used for the excitation and acquisition. A PC with the Teletest ® software was used to control the pulser-receiver and configure the experimental parameters. Input signals of five-cycle Hanning-windowed bursts with centre frequencies from 40 kHz to 140 kHz in steps of 5 kHz were used to excite the transmitter. This frequency range was chosen since it is the operating frequency range of the transducers used in the experiment. Signals were recorded with a 1 MHz sampling frequency; 32 averages were taken, and the record length of the signals was 1 ms. The experimental setup is presented in Figure 12. Two types of piezoelectric transducers (PI Ceramic GmbH) were used in the experiment, a shear piezoelectric element (d 15 type) and a compressional piezoelectric element (d 33 type), both with an active contact area of 13 mm by 3 mm. Figure 13 shows the direction of poling and corresponding shear deformation for the shear transducer. The shear transducer was used to evaluate the S 0 and SH 0 wave modes, and the A 0 wave mode was analyzed using the compressional transducer, which mainly produces an out-of-plane vibration. Using the shear transducer, S 0 waves are generated in the direction of the poling axis, and SH 0 in the perpendicular axis. A small amount of A 0 is also generated on the poling axis. In the case of the compressional transducer, S 0 and A 0 are created omnidirectionally. The averaged amplitudes of the S 0 and SH 0 modes using the shear transducer at the propagation distances specified above are 13 mV and 33 mV, respectively; and the amplitude of the A 0 mode is 63 mV using the compressional transducer. Experimental signals were acquired at a sampling frequency of 1 MHz providing a poor time resolution for the wave velocities of interest. Therefore, the sampling rate of the signals was increased by a factor of 10 using cubic spline interpolation to give a time accuracy of 0.1 µm (10 MHz). The interpolation is a necessary step when the propagation distance is a few cm, as the propagation time of the faster modes, such as S 0 , is of the same order of magnitude as the original sampling period. By reducing the sampling period by means of interpolation, the accuracy of the group velocity is increased. However, a small error is introduced in the signal due to the interpolation. Overall, the improvements gained for the velocity calculation are more significant than the errors introduced in the signals. The signal-to-noise ratios of the S 0 , A 0 and SH 0 signals were around 30 dB, 37 dB and 35 dB, respectively.

Experimental Results
The results presented in Figures 14 and 15 have been extracted from signals at two different spacings: 5 cm and 10 cm. The velocities extracted from the experimental signals correlate quite well with the theoretical velocities. However, there are slight variations around the theoretical velocities. These variations can be caused by a number of factors, such as the aforementioned interpolation and overlapping, but also by the manipulation of the transducers, which in this case were handled manually. Small variations in the correct spacing distance between transducers lead to errors in the calculation of the velocities. The velocity-distance relationship is linear so the relative error in the calculated velocity due to misplacement is the same as the relative distance error. An error of 1 mm in an intended 5 cm spacing distance will cause a 2% error in the velocity or an error of 110 ms −1 in the S 0 mode (≈ 5400 m/s), for example. On the other hand, with a 1 mm error transducer placement for a 10 cm spacing, the error will be of 1% of the velocity. Therefore, large errors can occur where small propagation distances are used. A 0 and SH 0 exhibit less absolute error variation than S 0 due to their slower velocities, the generated errors being proportional to the velocity.
In this paper, two propagation distances have been chosen: 5 and 10 cm. From the S 0 mode in Figure 15, it can be seen that the calculated velocities at 10 cm are more accurate, and the variance is also smaller. This is for two reasons: the first is the reduction in the error caused by an incorrect transducer placement, and the second is the implicit increase of the propagation time, which reduces the error in determining that time. Other factors that can introduce errors in the velocity calculation are the variations of the boundary conditions of the transducer, such as the applied pressure over the transducer at each receiving location, as well as the difference of transfer function of each transducer. This problem was overcome utilizing the same transducer to acquire the signals at all locations. Similarly, a change in temperature between measurements will cause velocity changes leading to erroneous values. This can be avoided by acquiring results in quick succession to minimize the impact of the temperature in the results. Using three transducers, one transmitter and two receivers would avoid the issues relating to moving the transducers.

Conclusions
This paper has presented a new methodology for creating dispersion curves based on measuring the phase difference and time lag between two pulses acquired at different distances. Firstly, the theoretical basis of the method has been presented where the phase velocity has been related to the phase at the pulse centre for systems of limited dispersion; then, a signal processing methodology to extract the phase shift and time lag between two signals and subsequently calculate the phase and group velocities has been outlined. Three different tests were performed to validate the formulation and methodology; using synthesized signals, signals from an FE model and experimental signals. In the three cases the performance achieved agrees well with the theoretical values.
There are some limitations to the proposed method for the extraction of the phase shift to calculate the dispersion curves. Highly dispersive modes, like high-order modes, are more challenging to evaluate; to mitigate this issue, much shorter spacing distances between transducers should be established to minimize the mode-shape distortion. Overlapping between wavepackets limits the method's applicability, since it distorts the apparent phase of the mode being analyzed. Hence, relatively large samples are required to avoid overlapping between wavepackets and edge echoes. The size of the sample can be smaller at higher frequencies. The velocity-frequency spectrum that can be evaluated using this method will be determined by the excitation mechanism used to generate the guided waves in the solid media. With appropriate experimental design, errors of less than 1% are obtained.
For future work, the proposed technique will be applied to more complex structures, such as composite plates where the elastic constants are not easy to obtain. In this case, by studying the phase and group velocities as a function of propagation direction, the full angular dependency of the dispersion curves will be obtained.