Modelling the Generation and Propagation of Ultrasonic Signals in Cylindrical Waveguides

Elongated cylindrical structures like rods, pipes, cable strands or fibers, support the propagation of mechanical waves at ultrasonic frequencies along their axes. This waveguide behaviour is used in a number of scientific and engineering applications: the Non Destructive Evaluation (NDE) of the structural health of civil engineering elements for safety purposes (Rose, 2000), in linear displacement sensors (Seco et al., 2009) for high accuracy absolute linear position estimation, in the evaluation of material properties of metal wires, optical fibers or composites (Nayfeh & Nagy, 1996), and as fluid sensors in pipes transporting liquids (Ma et al., 2007). These applications demand exact quantitative models of the processes of wave generation, propagation and reception of the ultrasonic signals in the waveguides.


Introduction
Elongated cylindrical structures like rods, pipes, cable strands or fibers, support the propagation of mechanical waves at ultrasonic frequencies along their axes. This waveguide behaviour is used in a number of scientific and engineering applications: the Non Destructive Evaluation (NDE) of the structural health of civil engineering elements for safety purposes (Rose, 2000), in linear displacement sensors (Seco et al., 2009) for high accuracy absolute linear position estimation, in the evaluation of material properties of metal wires, optical fibers or composites , and as fluid sensors in pipes transporting liquids (Ma et al., 2007). These applications demand exact quantitative models of the processes of wave generation, propagation and reception of the ultrasonic signals in the waveguides.
The mathematical treatment of mechanical wave propagation in cylindrical structures was provided by J. Pochhammer and C. Chree at the end of the XIX century, but its complexity prevented researchers from obtaining quantitative results until the advent of computers. D. Gazis (Gazis, 1959) reported the first exact solutions of the Pochhammer-Chree frequency equation, as well as a complete description of propagation modes and displacement and stress distributions for an isotropic elastic tube, found with an IBM 704 computer. Since then, the literature on the topic has grown steadily, and references are too numerous for this book chapter. We will only mention a few landmark developments: the study of multilayered waveguides beginning with a composite (two-layer) cylinder by H. D. McNiven in 1963; the extension of Gazis' results to anisotropic waveguides, initiated by I. Mirsky in 1965; the consideration of fluids and media with losses surrounding, or contained in the waveguides, beginning with V. A. Del Grosso in 1968; and finally, the demonstration of ultrasonic guided waves generated with electromagnetic transducers by W. Mohr and P. Holler in 1976, and piezoelectrically by M. Silk and K. Bainton in 1979, for the nondestructive testing of pipes.

Modelling the response of the waveguide to external excitation
Of particular importance for transducer design is the determination of the mechanical response of a waveguide when subjected to an external excitation. Several approaches exist to consider this problem.
Integral transform methods (Graff, 1991) convert the differential equations that physically model the excitation forces and the behaviour of the waveguide into a set of algebraic equations, which are more easily solvable. However, in order to find the actual distribution of the elastic field excited in the waveguide, inverse contour integration in the complex plane has to be performed, which is usually complicated. Due to the complexity of the Pochhammer-Chree equations, this procedure is only practical with simplified versions of the wave equation, which in general are not accurate enough for ultrasonic frequencies. See for example, Folk's solution for the transient response of a semi-infinite rod to a step pressure applied to its end (Folk et al., 1958).
The Semi-Analytical Finite Element (SAFE) method is a modification of Finite Element Methods (FEM) in which the elastic field is expanded as a superposition of harmonic waves in the azimuthal-axial (θ-z) plane, while discretized mechanical equations are used in the radial (r) direction of the waveguide. This reduction of the number of dimensions permits a much higher efficiency in the computation of the elastic fields (Hayashi et al., 2003). Waveguides surrounded by infinite media (like a pipe submerged in soil) can be handled by SAFE techniques with proper discretized elements (Jia et al., 2011), as well as waveguides with arbitrary profiles: for example, a railroad rail in (Damljanovic & Weaver, 2004). Although finite element methods are powerful and flexible, they have the shortcoming of great requirements on computer memory and processing time when large structures or high frequencies of operation are considered, and the difficulty encountered in the parameterization of transducer designs (for example, the determination of the transfer function of the transducer-waveguide coupling).
Spectral methods are another numerical technique which approximate the differential elastic equations of the waveguide (Helmholtz equations) by differentiation operators, turning the problem of finding the wavenumber-frequency roots into a matrix eigenvalues determination (Doyle, 1997). This numerical method, which is computationally simple and reportedly does not suffer from the problems associated with large diameter waveguides at high frequencies, has been recently applied to model multi-layered cylindrical waveguides (Karpfinger et al., 2008).
Modal analysis is an analytical method based on the expansion of the forcing terms acting in the waveguide into the set of its proper modes (Auld, 1973). In (Ditri & Rose, 1992), modal analysis is employed to model the loading of a waveguide by a transducer array. This treatment is extended to more general transducers and antisymmetric modes by (Li & Rose, 2001). Modal analysis is a mathematically exact technique that leads to a closed form integral equation for the elastic fields in the waveguide, and which incorporates in a natural way the issue of mode selectivity, offering insight on the physics of waveguide behaviour. For these reasons, modal analysis will be the approach used in this work.

Intention and scope of the research
With this book chapter we contribute a numerical simulation treatment of the ultrasonic behaviour of cylindrical waveguides, based on the Pochhammer-Chree (PC) theory, and covering the aspects of assembly of the description matrix of the waveguide, tracing of the frequency-wavenumber curves, computation of mode shapes, use of modal analysis to determine the response of the waveguide to external excitations, and the dispersive propagation of signals.
The work described here has resulted in a software package, named PCDISP, written in the Matlab environment (Matlab, 2004), and freely available 1 to be adapted to particular circumstances. The main features of the PCDISP software will be introduced in this chapter alongside with the theoretical concepts upon which it is based.
The purpose of PCDISP is freeing the researchers from the numerically delicate, time consuming issues arising in the solution of the PC equations, such as the creation of the waveguide matrix, the numerical instabilities encountered when the thickness of the waveguide or the operating frequency are high, the determination of proper modes and the tracing of the dispersive wavenumber-frequency curves. In this way, the researcher can concentrate in the study of the waveguide/transducer interaction as such.
As far as we are aware of, only two other software suites specifically designed for modelling elastic wave propagation in cylindrical waveguides exist. Disperse (Pavlakovic & Lowe, 1999) is a commercial package, based on matrix techniques, capable of analyzing cylindrical or plate waveguides made of perfectly elastic or damped solids, as well as fluids. GUIGUW (Bocchini et al., 2011) is a Matlab-based software which utilizes a SAFE-based approach to model ultrasonic propagation in cylindrical, plate, and arbitrary cross section waveguides. However, none of these computer solutions permit to model the waveguide response to external excitations.
The organization of this chapter is detailed next. Section 2 briefly reminds the mathematical background of the PC theory. Section 3 properly describes the main features of our methodology and how it is implemented in the PCDISP package. Two common transducer setups for the generation of ultrasonic waves are studied in section 4 with the help of PCDISP. Finally, we will offer some conclusions and point to lines in which this research could be further extended.

Background and nomenclature
In this section we present a summarized theoretical background on wave propagation in cylindrical waveguides, treating such aspects as relevant for our purposes; standard references can be consulted for further information (Graff, 1991;Meeker & Meitzler, 1972;Rose, 1999).
A waveguide is a physical structure which supports the propagation of mechanical waves along its elongated direction z, and modifies the behaviour of such waves with respect to free propagation in the bulk material. There are two fundamental characteristics of waveguide propagation. The first is the discretization of waves into propagating modes,o fw h i c ho n l y a finite number are permitted for a given frequency, and whose properties are determined by the shape of the cross section and boundary conditions of the waveguide. The second is the existence of dispersion, which is the nonlinear relationship between wavenumber and frequency. As a consequence, signals with a significant bandwidth are distorted as they travel along the waveguide, because their spectral components propagate at different phase speeds.
The solutions of the wave equation in a cylindrical material are readily found by the use of potentials and the technique of separation of variables, arriving at the following general form for the displacement vector ( u) and stress tensor ( σ): u(r, θ, z)= u(r, θ)e jkz = u(r)e jnθ e jkz σ(r, θ, z)= σ(r, θ)e jkz = σ(r)e jnθ e jkz ,( 1 ) where the cylindrical system is used (with coordinates (r, θ, z),a n du n i tv e c t o r s(e r , e θ , e z )), harmonic time variation e −jωt is assumed, and ω is the angular frequency, k the wavenumber, and integer n is a separation constant called the circumferential order, which determines the symmetry of the solutions in the azimuthal direction.
The radial dependent part of the displacement vector and stress tensor is expressed in matrix form as (Gazis, 1959): , shear vertical (SV), and horizontal (SH) deformation components, and the + and − terms stand for perturbations propagating in the direction of increasing and decreasing radius, respectively.
The coefficients of matrices D u and D σ are of the general form D ij (r; n, k, ω, c),withc being the elastic compliance tensor of the solid. These matrices are given explicitly in tables 1 and 2, for the case of a mechanically isotropic material. They have been obtained from the equations of motion with a symbolic computation program (Maple, 2007), and match those found in (Gazis, 1959), except for some typographical errors in the paper, also propagated to later works as (Graff, 1991). Matrices D u and D σ are implemented in the PCDISP package in function pcmat.
In tables 1 and 2, α 2 = ω 2 /c 2 vol − k 2 and β 2 = ω 2 /c 2 rot − k 2 ,wherec vol and c rot are respectively the volumetric and rotational speeds of the solid (Rose, 1999). Functions Z n (x) and W n (x) are two independent solutions of Bessel's differential equation, with, in general, complex arguments x = αr, βr. Of the possible choices for Z n (x) and W n (x), the numerical stability of the frequency equation determinant is increased when Bessel's ordinary functions J n (x) and Y n (x) are employed for real arguments, and the modified Bessel functions I n (x) and K n (x) for purely imaginary arguments. With this choice, the contributions to the elastic field from the standing waves propagating towards increasing (+) and decreasing (-) radius are separated (see section 3.1.4 for more on the stability of the frequency equation determinant).
To cope with the fact that the recurrence relationships between Bessel's ordinary functions are different from those of the modified functions (Abramowitz & Stegun, 1964), signs λ 1 and λ 2 are introduced, following the scheme of table 3. For complex wavenumbers, PCDISP uses the ordinary Bessel functions J n (x), Y n (x) with complex values, and λ 1 = λ 2 = 1.
The solutions of the wave equation are classified in family modes according to their symmetry properties, which depend on the circumferential index n of equation 1. Modes for which n = 0 have no dependence on the azimuthal coordinate θ and are labelled axisymmetric. They are further divided into torsional modes T(0, m) (which only involve the azimuthal component, and can be thought of as waves which twist the waveguide), and longitudinal modes L(0, m) (with both radial and axial components). Antisymmetric modes (n ≥ 1) are labelled flexural F(n, m), and involve all three components of the displacement vector. In general, multiple propagating modes exist for a given circumferential order and frequency, so a second index m is used to order them. Table 4 summarizes this information.

Methodology for the simulation of the waveguide behaviour
In this section we describe the methodology used to study waveguide generation and propagation of ultrasonic signals, discuss the numerical issues encountered, and the approach used in the PCDISP package. The full process consists in four stages: assembly of the waveguide description matrix, tracing of the dispersion curves, modal analysis of the excited modes, and modelling of the signal propagation along the waveguide.
While dealing with these topics, we will introduce the relevant PCDISP routines that should be used for the computations. Table 5 shows the components of the PCDISP software, arranged by their functionality. Throughout this chapter, we will use monospace fonts (like pcmat) to refer to programs of the package.

Wavenumber Frequency range Coefficients Bessel functions
Real  Table 4. Notation and non-null components of the amplitude coefficients, displacement vector, and stress tensor, for the three family modes of a cylindrical waveguide.

Single layer waveguides
Consider an isotropic tube of inner radius r int and outer radius r ext in vacuum or air. The boundary conditions specify that the traction part of the stress tensor is null in both surfaces of the tube (Gazis, 1959), so: which leads to the following matrix determinant equation: Equation 4 is called the frequency or characteristic equation of the waveguide, and its roots (ω, k) determine the proper modes supported by it. Once these roots are known, the vector of amplitude coefficients A is determined (up to a multiplicative constant) by solving the following homogeneous system of equations: Since matrix D(ω, k) is singular at the mode's frequency-wavenumber roots, a robust method for computing the amplitude, like the singular value decomposition (SVD), is recommended (Press et al., 1992). With A determined, the distribution of u(r) and σ(r) is computed by routine pcmatdet.
The original Pochhammer-Chree formulation was developed for the simple case of a one-layer isotropic waveguide in vacuum; this, however, represents just a fraction of the waveguides of practical importance. Waveguides may be constituted by several layers, might be built with  anisotropic materials, be embedded in the ground, or transport (or be surrounded by, or both) fluids. We will consider next the extensions of the PC theory which permit to model these situations.

Multilayered waveguides
Some waveguides are formed by several layers: for example, a metallic rod with external insulation, or a tube embedded in rock. The Pochhammer-Chree approach was first used to analyze a two-layer waveguide in (McNiven et al., 1963;Whittier & Jones, 1967), and later extended to laminated waveguides (formed by an arbitrary number of layers) in (Nelson et al., 1971). The modern technique to simulate multilayered waveguides is called the Multiple Layer Matrix (MLM) (Lowe, 1995), and is adapted from the transfer matrix and global matrix techniques developed by W.T. Thomson and L. Knopoff in the period 1950-1964 to study wave propagation in stratified media in seismology.
Following the MLM approach, we assemble a system of linear equations for the complete waveguide, which includes the equations of the elastic waves for each individual layer (obtained with pcmat), the equations which match the displacement and traction stresses at the interface between adjacent layers, and the boundary conditions. Consider the example of a multilayered waveguide shown in figure 1, where a solid cylinder of radius r 1 is enclosed by a tube of inner radius r 1 and outer radius r 2 , in turn surrounded by an infinite medium.
The corresponding system of equations is: Note that the radiation conditions are used to simplify the system matrix, leading to discard the + terms in region 1, since no waves can emanate from r = 0; similarly, the -terms in region 3 are not considered, as no energy comes from r = ∞; however, both outgoing and incoming terms are allowed in the middle region. In each case the unneeded columns of matrix D σt are removed. Next, this equation is solved by the SVD method to find the mode's amplitude vector, in the same way as the single layer waveguide of section 3.1.1.

Anisotropic and fluid-loaded waveguides
Materials with mechanical anisotropy are routinely employed to build waveguides, and their behaviour is modelled by taking into account the adequate compliance tensor for the material (Pollard, 1977). The important case of hexagonal symmetry is found in waveguides with transverse isotropy (with respect to the propagation axis z), like beryllium, or in fiber reinforced composite cylinders. Although in this case there are five elastic constants (up from two for an isotropic material), the mechanical fields are still separable with the treatment of Gazis for isotropic waveguides, with different coefficients for the D u and D σ matrices, obviously (Mirsky, 1965). In the case of purely orthotropic symmetry (three orthogonal planes of symmetry), the solution of the wave equation is not separable; still, closed solutions can be achieved in the form of a Frobenius power series (obtained in (Mirsky, 1964) for the axisymmetric case, and extended later in (Markus & Mead, 1995) to the asymmetric problem). A recent state of the art in the theory of mechanical waves in anisotropic cylindrical waveguides is found in (Grigorenko, 2005).
Materials with elastic losses (damping) are treated by the Kelvin-Voigt viscoelastic model (Lowe, 1995), which replaces the elements of the compliance tensor c of the material by operators which contain time derivatives, therefore modifying Hooke's law: where σ and ǫ are the stress and strain tensors, and component c ′′ of the compliance tensor models the viscoelastic losses. The solutions of the frequency equation for a waveguide with viscoelastic losses, are, by default, complex wavenumbers.
A case of particular practical importance is that of waveguides including a fluid layer (for example, a pipe carrying a fluid, submerged in a fluid, or both). The theoretical treatment depends on the viscosity of the fluid.
The influence of inviscid fluids (which do not support shear stresses) on cylindrical wave propagation is treated theoretically and experimentally in (Sinha et al., 1992). If the waveguide is submerged in a liquid, propagating modes with complex wavenumber appear, which radiate (leak) energy into the surrounding fluid. That does not happen for waveguides containing fluids and surrounded by vacuum, although the propagating modes themselves are modified from the unloaded situation.
A treatment for Newtonian viscous fluids which is compatible with the PC based formulation of wave propagation in cylinders is introduced in . In a form similar to the Kelvin-Voigt model, the viscous liquid is modelled as an isotropic solid whose compliance tensor includes complex elements: where λ is the compressibility of the fluid, and η its viscosity. This simple model has shown a good accuracy in predicting propagation in waveguides with viscous fluids (Aristegui et al., 2001). Indeed, the changes in the propagation of ultrasonic waves in a pipe caused by the presence of a fluid in its interior can be used to measure the longitudinal wave speed and the viscosity of the fluid (Ma et al., 2007).

The case of large frequency × thickness product
Solutions of the frequency equation become numerically unstable when the product frequency times thickness of the waveguide is high. Physically, this phenomenon arises because the standing waves established in the radial direction of the waveguide are formed by a combination of terms which increase exponentially with the radius r and others that decrease exponentially with it. Since it is the sum of both terms which must match the boundary conditions at r = r int and r = r ext , the dynamic range of positive and negative exponentials in the frequency equation will eventually overflow the numerical capacity of the machine if the radius or frequency increase. With the 64 bits double precision arithmetic of the IEEE 754 standard, we have found that a direct implementation of the frequency equation determinant fails when the product f · (r ext − r int ) is higher than approximately 30 MHz·mm. This threshold is easily reached in the NDE of piping with ultrasonic waves, where frequencies of a few megahertz are common (Rose, 2000).
From table 3, we can see that the problem arises when a mode's phase velocity falls below the volumetric (c vol ) or rotational (c rot ) speeds of the solid, making parameters α or β, respectively, become imaginary. This changes the radial dependence of the mode amplitude from the bounded Bessel functions J and Y, to the exponentially varying I and K (you can use pcviewmatdet to see the individual entries of the waveguide description matrix). The behaviour is different depending on the frequency range: •F o r c ph < c rot ,b ot hα and β are imaginary, and there exist two solutions of the frequency equation, corresponding to two Rayleigh modes propagating close to the outer and inner surfaces of the waveguide, and with amplitudes decaying exponentially from them. As the frequency increases, these two modes become decoupled, the opposite waveguide surface can ultimately be ignored, and the solutions are numerically stable. •F o r c rot < c ph < c vol , α is imaginary and β real. While the SV ± and SH ± terms remain bounded, the terms L + and L − decrease and increase, respectively, with the radius, as described above. The condition number of the waveguide matrix grows with the frequency, and its solution will eventually become unstable. •F o r c ph > c vol ,b o t hα and β are real, and the solution is stable (this is also the case with purely imaginary wavenumber).
Since the detection of the large fd problem in 1965, several solutions have been proposed to increase the numerical stability in plane waveguides (Lowe, 1995). We have not been able to locate similar studies for cylindrical waveguides, so, for the PCDISP software, we have developed an algorithm adapted from the transfer matrix and global matrix approaches and discussed here. We consider this method as a new contribution to the literature.
The cross section of the pipe is divided into L layers of equal thickness, where the l-th layer is given by r l−1 < r < r l ,andr 0 = r int and r L = r ext are the inner and outer radii of the pipe (see figure 2). In the l-th layer, the displacement vector and traction part of the stress tensor are given by: where the vector of amplitude coefficients We use the shorthand notation: and scale this matrix by columns for each layer as: where G l is a diagonal matrix whose entries are taken as 1/ max col|D l |.
The elastic field [u σ] T is propagated from the inner to the outer part of a layer by the following equation: where P l is the propagator matrix of layer l, and is given by: Applying equation 10 to all the layers of the waveguide, we can assemble a global matrix: where I 6 is the identity matrix of size 6 × 6.
The boundary conditions to match the wave fields to the surrounding medium (or to other layers in multilayered waveguides) are introduced at radii r int and r ext . In the case of stress free boundaries, the terms σ 0 and σ L are zero, and their corresponding columns are simply removed from the global matrix.
Once equation 12 has been solved, and we have determined the displacements and traction stresses at each layer boundary [u l , σ l ], the vector of amplitude coefficients for each layer is found by solving: for l = 1,...L. The unscaled amplitude vector is simply: A l = G l · A s l . The algorithm described increases the fdstability limit by a factor proportional to the number of layers L, at the expense of larger matrices and longer computational time. In PCIDSP, the algorithm described is written into routine pcmatdet, and is activated automatically if the waveguide defined in pcwaveguide consists of N L > 1 layers of the same material.

Computation of the dispersion curves
The roots of the waveguide's frequency equation represent the mechanical modes which satisfy the boundary conditions. The procedure for computing such solutions is described in this section.

Nature and ordering of the solutions of the frequency equation
For a given frequency f , the roots of the frequency equation are in general complex wavenumbers k = k r + jk i . Wavenumbers with k i = 0 correspond to the propagating or proper modes of the waveguide; those with k i = 0r e p r e s e n tevanescent modes which attenuate along the axial distance z. Due to the symmetry of the coefficients of the frequency equation, purely real or imaginary solutions appear in pairs (±k r and ±jk i ), while complex s o l u t i o n sd os oi nq u a r t e t s( ±k r ± jk i ). The principle of the conservation of energy dictates which solutions are valid in a waveguide problem. For example, waves which propagate towards the z+ axis of an infinite waveguide, will necessarily have Im{k}≥0. Although signal propagation does not occur for imaginary or complex wavenumbers, those solutions are needed to fulfill the condition of completeness which will be stated in section 3.3. Furthermore, stationary waves formed by combination of the two wavenumbers k r + jk i and −k r + jk i (provided that k i > 0), can exist locally, storing but not dissipating energy (Meeker & Meitzler, 1972). In ultrasonic applications, these waves are important in the region of generation of ultrasonic waves, at waveguide discontinuities (like defects), and at the waveguide ends.
Solutions of the frequency equation can be traced like continuous curves in the k-f space, where k = k r + jk i . As an example, we show the dispersion curves of the longitudinal modes L(0,m) of a sample waveguide with the data shown in table 6 (this waveguide will be further used in the examples of section 4). The complete spectrum, up to 3 MHz, has been computed with the pckfcurves routine of PCDISP, and is shown in figure 3. As it can be seen, the wavenumber curve of a given mode changes from a real value (shown in blue colour) to purely imaginary (coloured red, and projected into the negative wavenumber axis), or complex (plotted in green, with the real part on the positive k axis, and the imaginary part on the negative k axis), in a complicated fashion.
After finding all possible roots of the frequency equation, branches L(0,m) have to be ordered in such a way that each mode is assigned a unique, continuous wavenumber between the maximum frequency f max and zero frequency. Part (b) of figure 3 shows a reduced frequency range of the spectrum in part (a), and illustrates the convention for labelling modes. PCDISP uses the following rules about the behaviour of dispersion curves (Meeker & Meitzler, 1972): • Only the first torsional T(0,1), longitudinal L(0,1), and first flexural F(1,1), modes propagate down to zero frequency with real wavenumber. • Higher order modes switch from real to imaginary wavenumber at the cutoff frequencies when the wavenumber becomes null (k = 0, f = f cutoff ), and the phase speed infinite. • Miminum points (dω/dk = 0, d 2 ω/dk 2 > 0) of either real or imaginary wavenumber branches are also cutoff frequencies (with finite phase speed), from which complex wavenumber branches start. • Complex branches terminate either at zero frequency or at the maximum points (dω/dk = 0, d 2 ω/dk 2 < 0) of imaginary branches. • The branches are ordered such that k i is positive for modes propagating in the z+ direction, and that the sign of the group speed does not change along the curve (although it becomes null at the cutoff frequencies, and at the purely imaginary branches).

Algorithm for curve tracing in PCDISP
In order to generate continuous curves pckfcurves first traces the real wavenumber branches, which start at solutions found with pcsolvebisection at the f = f max and k = k max axes, and finish at the k = 0a x i s . I f one is only interested in propagating modes, this finishes the procedure, and the pcdisp routine can be used to plot the phase and group speeds of propagating modes in the specified frequency range. Otherwise, the next step consists in tracing branches with purely imaginary wavenumber, which start at points (0, f cutoff ) on the vertical axis of zero wavenumber, and finish when they reach the f = f max , k = jk max axes, or at another cutoff frequency in the k = 0axis. Parts (c) and (d) of figure 3 show the reason why a robust algorithm for curve tracing of the dispersion branches is needed, in order to avoid the apparent crossings between branches, such as modes L(0,7) and L(0,8), with real wavenumbers, and modes L(0,9) and L(0,10), with imaginary wavenumbers. The curve tracing method used in pckfcurves is shown in figure 4. The dispersion curve being traced is extrapolated from the three last computed points {i − 2, i − 1, i} to define an angular interval of width dθ in a circle of radius dr(i) centered in the last correctly determined point (i)( dr is the step size of the algorithm, in normalized coordinates of the k-f space). If a sign change is found in this interval, the algorithm proceeds with a bisection method to accurately estimate the position of point i + 1 (this is the normal situation). Otherwise, the dispersion curve might have undergone a sudden change of curvature, or another mode might have come very close to the one being traced, provoking multiple sign changes. In this case, the step dr(i) is decreased, or, if needed, points i − 1, i − 2, etc, are recomputed with a smaller step dr. Summarizing, the tracing algorithm of PCDISP keeps track of the curvature of the branch and the proximity of neighbouring branches, adjusting the interval step between consecutive points accordingly.
It must be pointed out that the frequency equation has spurious solutions at the lines with slope equal to the volumetric and rotational speeds of the solid (ω/k = c vol and ω/k = c rot ), which have to be removed by the root finding algorithm. In the case of multilayered waveguides, the same phenomenon happens for the speeds of each layer.
The dispersion curves are completed by tracing the complex wavenumber branches k = k r + jk i , starting from the extrema points of the real/imaginary wavenumber branches. Tracing the complex wavenumber branches needs more computational effort, since it is required to solve simultaneously for the real and imaginary parts of the wavenumber; we have obtained satisfactory results with Muller's method (Press et al., 1992) using also an adaptive step. the minimum point of the imaginary branch. At that point, it is joined by the branch corresponding to mode L(0,3), and both go down to zero frequency with negative complex conjugate wavenumbers k r + jk i and −k r + jk i .
Similarly, mode L(0,4) is cut off at the minimum point of the real wavenumber branch (k = 345 m −1 , f = 621.5 kHz). Mode L(0,5) is cut off at (k = 0, f = 698.6 kHz), changes to imaginary wavenumber until it reaches point (k = 0, f = 669.4 kHz), and then again to real, but negative,wavenumberdownto(k = 345 m −1 , f = 621.5 kHz), where it joins mode L(0,4). Mode L(0,5) has a negative wavenumber (and consequently, negative phase speed ω/k), in order to maintain a positive group velocity dω/dk, since this is a propagating mode in the z+ direction of the waveguide, according to the last rule of section 3.2.1. Below 621.5 kHz, the dispersion curves of modes L(0,4) and L(0,5) descend to zero frequency with negative complex conjugate wavenumbers, in the same way as modes L(0,2) and L(0,3).

Modal analysis
Modal analysis is a mathematical technique which permits to compute the dynamic response of a waveguide subject to arbitrary external forces, as an expansion of the excited wave over the set of normal modes of the waveguide, as defined in section 2 (Auld, 1973). Modal analysis is based upon two properties of normal modes: orthogonality, the existence of a scalar product which is null for any two different modes; and completeness, the capacity of the set of normal modes to span arbitrary waveforms in the waveguide.
For two different modes of the waveguide (1) and (2) of the form given in equation 1, the orthogonality relationship (Auld, 1973) establishes that: which is applicable to linear elastic materials and also to piezoelectric or magnetostrictive linear materials, assuming no elastic or dielectric losses. Later on this result will be generalized to include external forces and stresses.
Discarding the common factor e j(k 1 −k * 2 )z , integrating over the cross section D of the waveguide, and applying the divergence theorem, we find that: In equation 16, ∂D = ∂D int ∪ ∂D ext represents the inner and outer surfaces of the waveguide, and the normal unit vector e n is taken on each surface pointing out of the waveguide's interior, as shown in figure 5.
Because for proper modes the surface traction stress is null ( σ t = σ · e n = 0in∂D), the first integral of equation 16 is zero. Then a suitable scalar product of modes (1) and (2) is: In the right part of equation 17 we have assumed that the circumferential order n of modes (1) and (2) is the same; otherwise, P 12 is zero automatically due to the integration over the θ coordinate. The factor −jω/4 is introduced so that the quantity P 11 equals to the integral of the acoustic Poynting vector in the cross section of the waveguide, i.e., the power transported by the mode. For nonpropagating modes with k i = 0, P 11 is zero.
With this notation, equation 16 reduces to: which implies that P 12 = 0u n l e s sk 1 = k * 2 . In PCDISP, mode orthogonality can be verified with routine pcorthogonalcheck.
The second condition for modal analysis is completeness, which is based on the premise that an arbitrary perturbation in the waveguide can be expanded in the set of normal modes: where the modes are indexed by p, and no distinction has been made between propagating and evanescent modes. The problem lies in computing the set of coefficients a p (z), when the waveguide is under an arbitrary excitation composed of: 1. A vector force field f e (r, θ, z) acting on the bulk material of the tube (region D). In PCDISP, this vector force field is defined in pcextvolumforce. 2. A traction stress σ e (r, θ, z) applied to the tube surfaces (region ∂D).InPCDISP ,thissurface stress is defined in pcextsurfacestress.
In the case of existence of external fields, the orthogonality relationship, equation 14, must be generalized to (Auld, 1973): where f 1,2 (r, θ, z) represent the forcing terms. and equation 25 is changed into an ordinary differential equation, solvable by a standard change of variables, resulting in: where the integration takes place in the region R g where the generating terms f s and f v are not null, and z is the point where the ultrasonic signal is observed, in the direction of increasing z from region R g .
If p is a non-propagating mode, our computation method is changed slightly, since P p = 0. However, P p,p * = 0, and we can set q = p * , k q = k * p , and modify equations 26-28 accordingly.
As a summary, we have established the equations that permit to find the amplitude of the proper modes excited in the waveguide by an arbitrary set of external driving forces. These equations are used by routine pcmodalanalysis of PCDISP .

Propagation of waveforms in the waveguide
The modal analysis equations discussed in section 3.3 permit to obtain the frequency response of a transducer exciting the waveguide. In many applications we want to predict what ultrasonic waveforms will be obtained at a certain distance z from the excitation source, when the transducer is excited by a finite length time signal, i.e., to model the transient behaviour of the system.
The method to study the propagation of signal waveforms is relatively straightforward (Doyle, 1997). Let u(0, t) be the input signal in the transducer (placed in region R g ), and U(0, ω)=F [u(0, t)] its Fourier transform. When this signal excites the waveguide, we determine the corresponding volumetric forces and surface stresses, and evaluate terms a p (z) from equation 28 for each significant frequency component of U a n da l ln o r m a lm o d e s of the waveguide. Note that the terms a p (z) incorporate both the frequency response of the transducer itself (inside the integral term) and the effect caused by signal propagation (exp(jk p (ω)z)). Thus, the frequency components at a distance z from the generating region are given by: where the values of the Fourier transform for negative frequencies are taken as complex conjugate of the positive ones, in order to obtain a real time signal. The waveform at z is recovered by the inverse Fourier transform: The imaginary and complex wavenumber parts of the spectrum are required in equation 29 if the exciting signal has significant frequency content below cutoff of the propagating mode, and the measurement point is not far away from the transducer, as was described in section 3.2.1. The dispersive effect characteristic of waveguide propagation is frequently undesired in applications, since it implies a distortion of the original signal. One practical way of minimizing dispersion is to employ signals with narrow spectral content (typically by windowing a sine pulse train) with a central frequency in the region where the curve c ph (ω)= ω/k(ω) is relatively flat (Lowe et al., 1998). Where this is not feasible, compensation methods based on inverting the nonlinear k = k(ω) dependence have been developed (Wilcox, 2003).
The propagation of ultrasonic signals in the waveguide is simulated in PCDISP with routine pcsignalpropagation. When customising this routine, the user must be careful to zero pad the excitation signal u(0, t) at the end such that the resulting time window has enough duration to allow propagation of the signal for the distance between transducer and receiving point. Likewise, a sampling frequency high enough to cover all the spectral content of the signal should be used; in practical applications, oversampling the signal above the Nyquist rate is advantageous since it enhances the signal to noise ratio.

Demonstration of the methodology
In this section we will illustrate the use of the methodology described in this chapter and the developed software, to model the performance of a given transducer. The complete procedure is summarized in table 7, along with the needed routines of the PCDISP package.
Two transducer setups commonly found in ultrasonic guided wave applications will be analyzed (see figure 6). The first one is an electromagnetic acoustic transducer (EMAT) used to generate ultrasound in metallic waveguides without physical contact between the transducer and the sample; the second, an array of piezoelectric rings which generates ultrasound by mechanically loading the external tube surface. We begin by considering the mechanical behaviour of the sample waveguide.

Dispersive curves of longitudinal modes
We will continue to use the waveguide described in table 6. The complete signal spectrum of the longitudinal modes was already shown in figure 3; in figure 7 we plot the phase and group speeds in the range of 0 to 800 kHz. An important requisite for guided waves applications of ultrasound is the selection and exploitation of a single propagating mode, in a region where dispersive effects are minimum, since, as a general principle, an external force will excite all propagating modes existing within its bandwidth (Lowe et al., 1998). We will consider two possibilities: excitation of mode L(0,2) at frequency f 1 = 250 kHz, and use of mode L(0,3) at frequency f 2 = 565 kHz (see figure 7 b). The dispersion curves of these modes are relatively flat at these frequencies, and their group speeds are higher than those of other coexisting modes.
Computation of mode amplitude a p (z) (equation 28) Determine gain for frequency ω: U p (z, ω)=U(0, ω) · a p (z) End loop over frequency Set negative frequencies of the FFT: U p (z, −ω)=U * p (z, ω) Inverse Fourier transform: u p (z, t)=IFFT [U p (z, ω)] (eq. 30) End loop over propagating modes Sum over propagating modes u(z, t)=∑ p u p (z, t) Table 7. General method used for modal analysis computations, and required PCDISP routines.  The displacement and stress profiles of the selected modes are shown in figure 8. Their determination is important in NDE applications since the sensitivity of a propagating mode to a waveguide defect depends on the matching between the defect's shape and the mode profile (Ditri, 1994).
In order to minimize the influence of dispersion in the propagation of signals, we use as excitation signal a tone burst consisting of n cyc = 16 cycles of a central frequency modulated by a raised cosine window. This waveform does a good job in exciting a single frequency of the waveguide with a finite length signal and minimum sidelobes (Oppenheim et al., 1999).

Generation of ultrasonic waves with an electromagnetic acoustic transducer
Electromagnetic acoustic transducers (EMATs) are used to excite ultrasonic waves in metallic waveguides by non contact means (Cawley et al., 2004). Basically, an EMAT consists of a generating coil, which creates a dynamic field H(t) at ultrasonic frequencies, and a bias magnet providing a constant magnetic field H 0 , as shown in figure 6 (a). The physical phenomenon that couples the magnetic field with the elastic field in our non-ferromagnetic aluminium waveguide is the Lorentz force resulting from the interaction between the eddy currents J(t) induced in the tube by the dynamic field and the bias field H 0 , which create a volumetric force given by:  of its weakened side (Zhu, 2001). A TDPRA, shown in figure 6 (b), consists in a number of piezoelectric rings capable of exerting a pressure loading on the outer surface of the tube. The rings are organized into N p identical periods of N r rings each, with the length of a period matched to the wavelength λ of the mode to be excited. The rings of a period are connected to the same excitation source, but with a relative delay between them proportional to their position within the period of the TDPRA; this scheme is repeated throughout the TDPRA. This reinforces the wavelength matching of a single period and also creates constructive interference in the enhanced direction (the "head" of the array) and destructive in the other (the "tail").
To model numerically the TDPRA, we assume that each ring has width z w , with a separation z s between adjacent rings and vibrates in its thickness mode, exerting a pressure loading, σ rr = −p over the outer surface of the waveguide, constant for all frequencies. This leads to an axisymmetric loading (no θ dependence). In this case, the term corresponding to the surface loading (equation 26) is: 24 Ultrasonic Waves www.intechopen.com while the volumetric term is null. The total pressure p(z) is a sum over the N p periods of N r rings each: where z ci =( z w + z s )i + z w /2 is the position of the center of each ring, and (x)=1f o r |x| < 1/2, (x)=0for|x| > 1/2, is the rectangular function.
The parameters of the TDPRA must be tuned to the frequency to be excited. The load line of the TDPRA, given by c ph = N r (z w + z s ) f , is shown in figure 7 (a), along with the phase speed curves of the aluminum tube given in section 4.1, for two different designs: the first one with N p = 4, N r = 8, z w = 2.2 mm and z s = 0.4 mm, intended to excite mode L(0,2) at 250 kHz, and the second with N p = 5, N r = 6, z w = 1.5 mm, z s = 0.4 mm for excitation of mode L(0,3) at 565 kHz. The intersection points of these lines with the phase speed curves correspond to the frequencies for which the TDPRA achieves maximum efficiency in mode coupling.
The transducer gain of the first TDPRA is shown in figure 10 (a). For 250 kHz, the mode L(0,2) is effectively excited, and the excitation frequency can be fine tuned to make it coincide with a minimum of the amplitude of the L(0,1) mode, improving the dynamic range between the two modes. The simulation of the propagated wave (at a distance z = 1.5 m from the TDPRA) gives the expected results (part (b)). The results obtained with the second TDPRA design are shown in parts (c) and (d) of figure 10. Mode L(0,3) dominates in this case at 565 kHz, with a higher relative amplitude over the L(0,1) and L(0,2) modes. In this case, however, as the group velocities are similar, the received signals appear closer in the time view of part (d). In all cases, ultrasonic generation in the "head" side has higher amplitude than in the "tail" side.
As a conclusion, the TDPRA is an efficient transducer for generating ultrasonic signals in cylindrical waveguides, and its feature of phase and wavelength matching permits to excite modes selectively, fine tune the system gain to a desired frequency, and direct the generated signal in only one direction. Although in this communication we have concerned ourselves only with axisymmetric transducers, PCDISP can also be used to study excitation of nonsymmetric modes by piezoelectric arrays (see reference (Li & Rose, 2001) for an example).

Conclusions
In this chapter we have presented a methodology to model the dynamic response of a waveguide of cylindrical symmetry when subject to an arbitrary set of external forces acting at ultrasonic frequencies, based on the combination of the mechanical Pochhammer-Chree equations and modal analysis techniques. Furthermore, a software package (named PCDISP), created in the Matlab environment, is offered freely with the intention of saving other researchers from the time needed for implementation of the PC theory equations, permitting them to focus on their particular problems.
Throughout this communication, we have paid special attention to the numerical issues of stability of the matrix determinant for large frequency thickness products, provided algorithms for robust root solving and tracing of the dispersion curves, and modelled the dispersive effect of the waveguide on signal propagation. The methods described in this chapter are valid for waveguides formed by any number of layers as long as they have cylindrical symmetry. The PCDISP software can be further extended to consider materials with anisotropy (transversely isotropic and orthotropic), as well as materials with elastic damping and waveguides surrounded by, or containing, fluids. Guidelines for such extensions are given in the text.
The performance of modal analysis is illustrated by studying two common transducers employed in guided wave ultrasonic applications: an electromagnetic-acoustic transducer and a time-delay piezoelectric ring array. We believe that transducer analysis with quantitative results is achieved comparatively easier and faster than with other competing techniques like spectral or finite element methods, obtaining significant time savings in the design stage of ultrasonic transducers. Ultrasonic waves are well-known for their broad range of applications. They can be employed in various fields of knowledge such as medicine, engineering, physics, biology, materials etc. A characteristic presented in all applications is the simplicity of the instrumentation involved, even knowing that the methods are mostly very complex, sometimes requiring analytical and numerical developments. This book presents a number of stateof-the-art applications of ultrasonic waves, developed by the main researchers in their scientific fields from all around the world. Phased array modelling, ultrasonic thrusters, positioning systems, tomography, projection, gas hydrate bearing sediments and Doppler Velocimetry are some of the topics discussed, which, together with materials characterization, mining, corrosion, and gas removal by ultrasonic techniques, form an exciting set of updated knowledge. Theoretical advances on ultrasonic waves analysis are presented in every chapter, especially in those about modelling the generation and propagation of waves, and the influence of Goldberg's number on approximation for finite amplitude acoustic waves. Readers will find this book ta valuable source of information where authors describe their works in a clear way, basing them on relevant bibliographic references and actual challenges of their field of study.